Restrict rounding to the computed corners (#143)

Rounding the transformation makes assumptions about the nature of
these transformations. However, in principle any kind of function
mapping points to new points should be supported.  It seems safer to
perform rounding on the *result* of the transformation.

Author: timholy

Project.toml

@@ -19,18 +19,19 @@ ColorVectorSpace = "0.6, 0.7, 0.8, 0.9"
 CoordinateTransformations = "0.5, 0.6"
 ImageBase = "0.1.1"
 ImageCore = "0.8.1, 0.9"
-Interpolations = "0.9, 0.10, 0.11, 0.12, 0.13"
+Interpolations = "0.13.4"
 OffsetArrays = "0.10, 0.11, 1.0.1"
 Rotations = "0.12, 0.13, 1.0"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 julia = "1"
 
 [extras]
+EndpointRanges = "340492b5-2a47-5f55-813d-aca7ddf97656"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 
 [targets]
-test = ["ImageMagick", "LinearAlgebra", "ReferenceTests", "Test", "TestImages"]
+test = ["EndpointRanges", "ImageMagick", "LinearAlgebra", "ReferenceTests", "Test", "TestImages"]

src/autorange.jl

@@ -34,7 +34,6 @@ autorange(A, tform)
 """
 autorange(A::AbstractArray, tform) = autorange(CartesianIndices(A), tform)
 function autorange(R::CartesianIndices, tform)
-    tform = _round(tform)
     mn = mx = tform(SVector(first(R).I))
     for I in CornerIterator(R)
         x = tform(SVector(I.I))
@@ -46,7 +45,7 @@ function autorange(R::CartesianIndices, tform)
     _autorange(Tuple(mn), Tuple(mx))
 end
 
-@noinline _autorange(mn,mx) = map((a,b)->floor(Int,a):ceil(Int,b), mn, mx)
+@noinline _autorange(mn,mx) = map((a,b)->floor(Int,_round(a)):ceil(Int,_round(b)), mn, mx)
 
 ## Iterate over the corner-indices of a rectangular region
 struct CornerIterator{I<:CartesianIndex}
@@ -106,23 +105,13 @@ try_static(tfm::Translation{<:SVector}, img::AbstractArray{T,N}) where {T,N} = t
 try_static(tfm::Translation{<:AbstractVector}, img::AbstractArray{T,N}) where {T,N} =
     Translation(SVector{N}(tfm.translation))
 
-# Slightly round/discretize the transformation so that the warpped image size isn't affected by
+# Slightly round/discretize the output coordinates so that the warped image size isn't affected by
 # numerical stability
 # https://github.com/JuliaImages/ImageTransformations.jl/issues/104
 _default_digits(::Type{T}) where T<:Number = _default_digits(floattype(T))
 # these constants come from eps() digits
-_default_digits(::Type{<:AbstractFloat}) = 15
-_default_digits(::Type{Float64}) = 15
-_default_digits(::Type{Float32}) = 7
+_default_digits(::Type{T}) where T<:AbstractFloat = ceil(Int, -log10(sqrt(eps(T))))
+_default_digits(::Type{Float64}) = 8
+_default_digits(::Type{Float32}) = 4
 
-function _round(tform::T; kwargs...) where T<:CoordinateTransformations.Transformation
-    rounded_fields = map(Base.OneTo(fieldcount(T))) do i
-        __round(getfield(tform, i); kwargs...)
-    end
-    T(rounded_fields...)
-end
-_round(tform; kwargs...) = tform
-
-__round(x; kwargs...) = x
-__round(x::AbstractArray; digits=_default_digits(eltype(x)), kwargs...) = round.(x; digits=digits, kwargs...)
-__round(x::T; digits=_default_digits(T), kwargs...) where T<:Number = round(x; digits=digits, kwargs...)
+_round(x::T; digits=_default_digits(T), kwargs...) where T<:Number = round(x; digits=digits, kwargs...)

src/compat.jl

@@ -1,12 +1,6 @@
 # This file includes two kinds of codes
 #   - Codes for backward compatibility
-#   - Glue codes that might nolonger be necessary in the future
-
-# patch for issue #110
-if isdefined(Base, :ComposedFunction) # Julia >= 1.6.0-DEV.85
-    # https://github.com/JuliaLang/julia/pull/37517
-    _round(tform::ComposedFunction; kwargs...) = _round(tform.outer; kwargs...) ∘ _round(tform.inner; kwargs...)
-end
+#   - Glue codes that might no longer be necessary in the future
 
 @static if !isdefined(Base, :IdentityUnitRange)
     const IdentityUnitRange = Base.Slice

src/invwarpedview.jl

@@ -22,7 +22,7 @@ struct InvWarpedView{T,N,A,F,I,FI<:Transformation,E} <: AbstractArray{T,N}
 end
 
 function InvWarpedView(inner::WarpedView{T,N,TA,F,I,E}) where {T,N,TA,F,I,E}
-    tinv = _round(inv(inner.transform))
+    tinv = inv(inner.transform)
     InvWarpedView{T,N,TA,F,I,typeof(tinv),E}(inner, tinv)
 end
 

src/warp.jl

@@ -54,7 +54,7 @@ This approach is known as backward mode warping. It is called "backward" because
 the internal coordinate transformation is actually an inverse map from `axes(imgr)` to `axes(img)`.
 
 You can manually specify interpolation behavior by constructing `AbstractExtrapolation` object
-and passing it to `warp` as `img`. However, this is usually cumbersome. For this reason, there 
+and passing it to `warp` as `img`. However, this is usually cumbersome. For this reason, there
 are two keywords `method` and `fillvalue` to conveniently construct an `AbstractExtrapolation`
 object during `warp`.
 
@@ -163,7 +163,6 @@ function warp(img::AbstractExtrapolation{T}, tform, inds::Tuple = autorange(img,
 end
 
 function warp!(out, img::AbstractExtrapolation, tform)
-    tform = _round(tform)
     @inbounds for I in CartesianIndices(axes(out))
         # Backward mode:
         #   1. get the target index `I` of `out`

src/warpedview.jl

@@ -19,7 +19,6 @@ function WarpedView(
         tform::Transformation,
         inds=autorange(A, inv(tform)); kwargs...) where {T,N,}
     etp = box_extrapolation(A; kwargs...)
-    tform = _round(tform)
     WarpedView{T,N,typeof(A),typeof(tform),typeof(inds),typeof(etp)}(A, tform, inds, etp)
 end
 

test/warp.jl

@@ -1,4 +1,5 @@
 using CoordinateTransformations, Rotations, TestImages, ImageCore, StaticArrays, OffsetArrays, Interpolations, LinearAlgebra
+using EndpointRanges
 using Test, ReferenceTests
 
 include("twoints.jl")
@@ -181,7 +182,8 @@ img_camera = testimage("camera")
         @test axes(wv2) == axes(img_camera)
         @test eltype(wv2) === eltype(img_camera)
         @test parent(wv2) === img_camera
-        @test wv2 ≈ img_camera
+        @test_skip wv2 ≈ img_camera      # see discussion in #143
+        @test wv2[ibegin+1:iend-1,ibegin+1:iend-1] ≈ img_camera[ibegin+1:iend-1,ibegin+1:iend-1]  # TODO: change to begin/end, drop EndpointRanges
 
         imgr = @inferred(invwarpedview(img_camera, tfm))
         @test imgr == @inferred(InvWarpedView(img_camera, tfm))
@@ -312,6 +314,14 @@ img_camera = testimage("camera")
         @test any(isnan, v)
         @test parent(v) isa InvWarpedView
     end
+
+    @testset "3d warps" begin
+        img = testimage("mri")
+        θ = π/8
+        tfm = AffineMap(RotZ(θ), (I - RotZ(θ))*center(img))
+        imgr = warpedview(img, tfm, axes(img))
+        @test axes(imgr) == axes(img)
+    end
 end
 
 img_pyramid = Gray{Float64}[