Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Clarification of rotations #147

Closed
henry2004y opened this issue Oct 9, 2021 · 7 comments · Fixed by #148
Closed

Clarification of rotations #147

henry2004y opened this issue Oct 9, 2021 · 7 comments · Fixed by #148

Comments

@henry2004y
Copy link

henry2004y commented Oct 9, 2021

Hi,

I have a Python script in which there is a call to scipy.ndimage.rotate for rotating a 2D array, similar to the following

using PyCall
ndimage = pyimport("scipy.ndimage")
# test array
img = Float64[1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
img_90 = ndimage.rotate(img, 90, reshape=false)

I came across this package when searching for the same kind of operation on a plain 2D array in Julia. However, soon I got confused by the rotation methods provided here, especially imrotate.

using ImageTransformations
using Rotations: RotMatrix
using CoordinateTransformations: recenter

# test array
img = Float64[1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]

# define transformation
trfm = recenter(RotMatrix(pi/2), center(img))
imgw = warp(img, trfm) # this is what I expected, same as in scipy

imgr = imrotate(img, pi/2) # this gives a 6x6 OffsetArray with values that I don't quite understand

where imrotate returned

6×6 OffsetArray(::Matrix{Float64}, 0:5, 0:5) with eltype Float64 with indices 0:5×0:5:
 NaN  NaN    NaN        NaN        NaN        NaN
 NaN  NaN    NaN          4.99995  NaN        NaN
 NaN   14.0    9.99998    5.99999  NaN        NaN
 NaN  NaN     11.0        7.00002    3.00003  NaN
 NaN  NaN     12.0      NaN        NaN        NaN
 NaN  NaN    NaN        NaN        NaN        NaN

I tried the example given in the doc for imrotate, and it looked just fine.

Can you kindly point out what I missed here? What is the closest possible translation from scipy.ndimage.rotate?

@johnnychen94
Copy link
Member

johnnychen94 commented Oct 9, 2021

There are two issues here:

  • There are some numerical errors on pi/2 because you can't represent pi arbitrarily precisely. (see also imrotate(img, pi/2) changes size of image #145)
  • The default imrotate output is the minimal array that encloses all pixels in the original image

There are two options to workaround this:

  • use the Base function rotr90(img) (also rotl90 and rot180)
  • clip the output by passing axes: imrotate(img, pi/2, axes(img))

I'm going to close this issue as it's very much duplicated of #145

@henry2004y
Copy link
Author

henry2004y commented Oct 9, 2021

Rotating the 2D array by 90 degree is just for testing purposes: the actual angle is in principle any value in [0, 2pi). What puzzles me is the behavior of imrotate, especially its not NaN structure after rotation. As a first time user to this function, I believe I am not the only one to be confused.

Even with axes:

julia> imrotate(img, pi/2, axes(img))
4×4 Matrix{Float64}:
 NaN    NaN          4.99995  NaN
  14.0    9.99998    5.99999  NaN
 NaN     11.0        7.00002    3.00003
 NaN     12.0      NaN        NaN

I am not sure if using any kind of low order interpolation (e.g. Closest()) will avoid the NaN fillings.

Since my goal is to reproduce the result in scipy, I guess I should go with warp then?

@timholy
Copy link
Member

timholy commented Oct 9, 2021

I'm reopening because the differences are due to the strategy for handling roundoff error in #58 (comment):

julia> using Debugger

julia> @enter imrotate(img, π/2, axes(img))
In #imrotate#9(kwargs, , img, θ, inds) at /home/tim/.julia/packages/ImageTransformations/O9GYH/src/warp.jl:245
 245  function imrotate(img::AbstractArray{T}, θ::Real, inds::Union{Tuple, Nothing} = nothing; kwargs...) where T
 246      # TODO: expose rotation center as a keyword
>247      θ = floor(mod(θ,2pi)*typemax(Int16))/typemax(Int16) # periodic discretezation
 248      tform = recenter(RotMatrix{2}(θ), center(img))
 249      # Use the `nothing` trick here because moving the `autorange` as default value is not type-stable
 250      inds = isnothing(inds) ? autorange(img, inv(tform)) : inds
 251      warp(img, tform, inds; kwargs...)
 252  end

About to run: (*)(2, π)
1|debug> n
In #imrotate#9(kwargs, , img, θ, inds) at /home/tim/.julia/packages/ImageTransformations/O9GYH/src/warp.jl:245
 245  function imrotate(img::AbstractArray{T}, θ::Real, inds::Union{Tuple, Nothing} = nothing; kwargs...) where T
 246      # TODO: expose rotation center as a keyword
 247      θ = floor(mod(θ,2pi)*typemax(Int16))/typemax(Int16) # periodic discretezation
>248      tform = recenter(RotMatrix{2}(θ), center(img))
 249      # Use the `nothing` trick here because moving the `autorange` as default value is not type-stable
 250      inds = isnothing(inds) ? autorange(img, inv(tform)) : inds
 251      warp(img, tform, inds; kwargs...)
 252  end

About to run: (Core.apply_type)(RotMatrix, 2)
1|julia> θ
1.5707876827295755

1|julia> π/2
1.5707963267948966

That's a pretty large difference from what the user entered.

@timholy timholy reopened this Oct 9, 2021
timholy added a commit that referenced this issue Oct 9, 2021
@henry2004y
Copy link
Author

Thanks! Another tiny issue I want to mention is the fillvalue keyword:

julia> imrotate(img, pi/2, axes(img), fillvalue=0.0)
4×4 Matrix{Float64}:
 NaN    NaN          4.99995  NaN
  14.0    9.99998    5.99999  NaN
 NaN     11.0        7.00002    3.00003
 NaN     12.0      NaN        NaN

Seemingly contrary to the method doc

  • fillvalue: the value that used to fill the new region. The default value is NaN if possible, otherwise is 0.

I don't observe any changes to the output in this case.

timholy added a commit that referenced this issue Oct 10, 2021
@timholy
Copy link
Member

timholy commented Oct 10, 2021

I can't replicate that; in a revised version of #148 I added these tests:

@test  any(isnan, imrotate(img, π/3))
@test !any(isnan, imrotate(img, π/3; fillvalue=0.0))
@test !any(isnan, imrotate(img, π/3, fillvalue=0.0))
@test  any(isnan, imrotate(img, π/3, axes(img)))
@test !any(isnan, imrotate(img, π/3, axes(img); fillvalue=0.0))
@test !any(isnan, imrotate(img, π/3, axes(img), fillvalue=0.0))

and they all pass.

@henry2004y
Copy link
Author

Sorry I just found that I was accidentally using an earlier version 0.8.13, probably due to some restrictions of other packages.

@timholy
Copy link
Member

timholy commented Oct 10, 2021

We really need to get Images 0.25 out!

julia> using RegistryCompatTools

help?> RegistryCompatTools
search: RegistryCompatTools

  RegistryCompatTools makes it easier to discover packages that need [compat] updates.

    •  held_back_by("SomePkg") - which packages are holding back SomePkg?

    •  held_back_packages()["SomePkg"] - which packages does SomePkg hold back?

    •  find_julia_packages_github() - discover which packages I can push to

julia> held_back_by("ImageTransformations")
15-element Vector{String}:
 "ADI"
 "Asimov"
 "Augmentor"
 "BasicTextRender"
 "CameraCalibrations"
 "DataAugmentation"
 "EasyML"
 "HCIToolbox"
 "Images"                    # not good...
 "MIRT"
 "MakieGallery"
 "Photometry"
 "UNet"
 "WordCloud"
 "YOLO"

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants