refactor warp and WarpedView to proper backward mode

src/ImageTransformations.jl

@@ -21,13 +21,13 @@ export
     imresize,
     center,
     warp,
-    WarpedView,
-    warpedview
+    InvWarpedView,
+    invwarpedview
 
 include("autorange.jl")
 include("resizing.jl")
 include("warp.jl")
-include("warpedview.jl")
+include("invwarpedview.jl")
 
 @inline _getindex(A, v::StaticVector) = A[convert(Tuple, v)...]
 

src/invwarpedview.jl

@@ -0,0 +1,74 @@
+immutable InvWarpedView{T,N,A<:AbstractArray,F1<:Transformation,I,F2<:Transformation,E<:AbstractExtrapolation} <: AbstractArray{T,N}
+    parent::A
+    transform::F1
+    indices::I
+    inverse::F2
+    extrapolation::E
+
+    function (::Type{InvWarpedView{T,N,TA,F,I}}){T,N,TA<:AbstractArray,F<:Transformation,I<:Tuple}(
+            parent::TA,
+            tform::F,
+            indices::I)
+        @assert eltype(parent) == T
+        etp = _box_extrapolation(parent)
+        tinv = inv(tform)
+        new{T,N,TA,F,I,typeof(tinv),typeof(etp)}(parent, tform, indices, tinv, etp)
+    end
+end
+
+function InvWarpedView(inner::InvWarpedView, outer_tform::Transformation)
+    tform = compose(outer_tform, inner.transform)
+    A = parent(inner)
+    inds = autorange(A, tform)
+    InvWarpedView(A, tform, inds)
+end
+
+function InvWarpedView{T,N,F<:Transformation,I<:Tuple}(
+        A::AbstractArray{T,N},
+        tform::F,
+        inds::I = autorange(A, tform))
+    InvWarpedView{T,N,typeof(A),F,I}(A, tform, inds)
+end
+
+Base.parent(A::InvWarpedView) = A.parent
+@inline Base.indices(A::InvWarpedView) = A.indices
+
+@compat Compat.IndexStyle{T<:InvWarpedView}(::Type{T}) = IndexCartesian()
+@inline Base.getindex{T,N}(A::InvWarpedView{T,N}, I::Vararg{Int,N}) =
+    _getindex(A.extrapolation, A.inverse(SVector(I)))
+
+Base.size(A::InvWarpedView)    = OffsetArrays.errmsg(A)
+Base.size(A::InvWarpedView, d) = OffsetArrays.errmsg(A)
+
+function ShowItLikeYouBuildIt.showarg(io::IO, A::InvWarpedView)
+    print(io, "InvWarpedView(")
+    showarg(io, parent(A))
+    print(io, ", ")
+    print(io, A.transform)
+    print(io, ')')
+end
+
+Base.summary(A::InvWarpedView) = summary_build(A)
+
+"""
+TODO
+"""
+@inline invwarpedview(A::AbstractArray, tform::Transformation, args...) =
+    InvWarpedView(A, tform, args...)
+
+function invwarpedview{T}(
+        A::AbstractArray{T},
+        tform::Transformation,
+        degree::Union{Linear,Constant},
+        fill::FillType = _default_fill(T),
+        args...)
+    invwarpedview(_box_extrapolation(A, degree, fill), tform, args...)
+end
+
+function invwarpedview(
+        A::AbstractArray,
+        tform::Transformation,
+        fill::FillType,
+        args...)
+    invwarpedview(A, tform, Linear(), fill, args...)
+end

src/warp.jl

@@ -1,5 +1,29 @@
+# The default values used by extrapolation for off-domain points
+@compat const FillType = Union{Number,Colorant,Flat,Periodic,Reflect}
+@compat const FloatLike{T<:AbstractFloat} = Union{T,AbstractGray{T}}
+@compat const FloatColorant{T<:AbstractFloat} = Colorant{T}
+@inline _default_fill{T<:FloatLike}(::Type{T}) = convert(T, NaN)
+@inline _default_fill{T<:FloatColorant}(::Type{T}) = nan(T)
+@inline _default_fill{T}(::Type{T}) = zero(T)
+
+_box_extrapolation(etp::AbstractExtrapolation) = etp
+
+function _box_extrapolation{T}(itp::AbstractInterpolation{T}, fill::FillType = _default_fill(T))
+    etp = extrapolate(itp, fill)
+    _box_extrapolation(etp)
+end
+
+function _box_extrapolation{T,N,D<:Union{Linear,Constant}}(parent::AbstractArray{T,N}, degree::D = Linear(), args...)
+    itp = Interpolations.BSplineInterpolation{T,N,typeof(parent),BSpline{D},OnGrid,0}(parent)
+    _box_extrapolation(itp, args...)
+end
+
+function _box_extrapolation(parent::AbstractArray, fill::FillType)
+    _box_extrapolation(parent, Linear(), fill)
+end
+
 """
-    warp(img, tform, [degree = Linear()], [fill = NaN]) -> imgw
+    warp(img, tform, [indices], [degree = Linear()], [fill = NaN]) -> imgw
 
 Transform the coordinates of `img`, returning a new `imgw` satisfying
 `imgw[x] = img[tform(x)]`. `tform` should be defined using
@@ -48,7 +72,7 @@ julia> indices(img)
 julia> tfm = recenter(RotMatrix(pi/4), center(img))
 AffineMap([0.707107 -0.707107; 0.707107 0.707107], [347.01,-68.7554])
 
-julia> imgw = warp(img, tfm);
+julia> imgw = warp(img, inv(tfm));
 
 julia> indices(imgw)
 (-196:709,-68:837)
@@ -59,7 +83,7 @@ julia> img0 = OffsetArray(img, -30:481, -384:383);  # origin near top of image
 julia> rot = LinearMap(RotMatrix(pi/4))
 LinearMap([0.707107 -0.707107; 0.707107 0.707107])
 
-julia> imgw = warp(img0, rot);
+julia> imgw = warp(img0, inv(rot));
 
 julia> indices(imgw)
 (-293:612,-293:611)
@@ -70,39 +94,26 @@ julia> indices(imgr)
 (Base.OneTo(906),Base.OneTo(905))
 ```
 """
-function warp{T}(img::AbstractExtrapolation{T}, tform)
-    inds = autorange(img, tform)
+function warp{T}(img::AbstractExtrapolation{T}, tform, inds::Tuple = autorange(img, inv(tform)))
     out = OffsetArray(Array{T}(map(length, inds)), inds)
     warp!(out, img, tform)
 end
 
 function warp!(out, img::AbstractExtrapolation, tform)
-    tinv = inv(tform)
     @inbounds for I in CartesianRange(indices(out))
-        out[I] = _getindex(img, tinv(SVector(I.I)))
+        out[I] = _getindex(img, tform(SVector(I.I)))
     end
     out
 end
 
-# The default values used by extrapolation for off-domain points
-@compat const FillType = Union{Number,Colorant,Flat,Periodic,Reflect}
-@compat const FloatLike{T<:AbstractFloat} = Union{T,AbstractGray{T}}
-@compat const FloatColorant{T<:AbstractFloat} = Colorant{T}
-@inline _default_fill{T<:FloatLike}(::Type{T}) = convert(T, NaN)
-@inline _default_fill{T<:FloatColorant}(::Type{T}) = nan(T)
-@inline _default_fill{T}(::Type{T}) = zero(T)
-
-function warp{T,N}(img::AbstractArray{T,N}, tform, degree::Union{Linear,Constant} = Linear(), args...)
-    itp = Interpolations.BSplineInterpolation{T,N,typeof(img),Interpolations.BSpline{typeof(degree)},OnGrid,0}(img)
-    warp(itp, tform, args...)
-end
-
-function warp{T,N}(img::AbstractArray{T,N}, tform, fill::FillType, args...)
-    warp(img, tform, Linear(), fill)
+function warp{T,N}(img::AbstractArray{T,N}, tform, inds::Tuple, args...)
+    etp = _box_extrapolation(img, args...)
+    warp(etp, tform, inds)
 end
 
-function warp{T}(img::AbstractInterpolation{T}, tform, fill::FillType = _default_fill(T))
-    warp(extrapolate(img, fill), tform)
+function warp{T,N}(img::AbstractArray{T,N}, tform, args...)
+    etp = _box_extrapolation(img, args...)
+    warp(etp, tform)
 end
 
 # This is type-piracy, but necessary if we want Interpolations to be

test/warp.jl

@@ -15,36 +15,32 @@ img_camera = testimage("camera")
     ref_inds = (-78:591, -78:591)
 
     @testset "warp" begin
-        imgr = @inferred(warp(img_camera, tfm))
+        imgr = @inferred(warp(img_camera, inv(tfm)))
         @test indices(imgr) == ref_inds
         @test eltype(imgr) == eltype(img_camera)
 
-        imgr = @inferred(warp(img_camera, tfm, 1))
+        imgr = @inferred(warp(img_camera, inv(tfm), 1))
         @test eltype(imgr) == eltype(img_camera)
-        imgr2 = @inferred warp(imgr, inv(tfm))
+        imgr2 = @inferred warp(imgr, tfm)
         @test eltype(imgr2) == eltype(img_camera)
         # look the same but are not similar enough to pass test
         # @test imgr2[indices(img_camera)...] ≈ img_camera
     end
 
-    @testset "WarpedView" begin
-        wv = @inferred(WarpedView(img_camera, tfm))
-        @test summary(wv) == "-78:591×-78:591 WarpedView(::Array{Gray{N0f8},2}, AffineMap([0.92388 0.382683; -0.382683 0.92388], [-78.6334,$(SPACE)117.683])) with element type ColorTypes.Gray{FixedPointNumbers.Normed{UInt8,8}}"
+    @testset "InvWarpedView" begin
+        wv = @inferred(InvWarpedView(img_camera, tfm))
+        @test summary(wv) == "-78:591×-78:591 InvWarpedView(::Array{Gray{N0f8},2}, AffineMap([0.92388 0.382683; -0.382683 0.92388], [-78.6334,$(SPACE)117.683])) with element type ColorTypes.Gray{FixedPointNumbers.Normed{UInt8,8}}"
         @test_throws ErrorException size(wv)
         @test_throws ErrorException size(wv, 1)
         @test indices(wv) == ref_inds
         @test eltype(wv) === eltype(img_camera)
-        @test typeof(parent(wv)) <: Interpolations.AbstractExtrapolation
-        @test typeof(parent(wv).itp) <: Interpolations.AbstractInterpolation
-        @test parent(wv).itp.coefs === img_camera
+        @test parent(wv) === img_camera
 
         # check nested transformation using the inverse
-        wv2 = @inferred(WarpedView(wv, inv(tfm)))
+        wv2 = @inferred(InvWarpedView(wv, inv(tfm)))
         @test indices(wv2) == indices(img_camera)
         @test eltype(wv2) === eltype(img_camera)
-        @test typeof(parent(wv2)) <: Interpolations.AbstractExtrapolation
-        @test typeof(parent(wv2).itp) <: Interpolations.AbstractInterpolation
-        @test parent(wv2).itp.coefs === img_camera
+        @test parent(wv2) === img_camera
         @test wv2 ≈ img_camera
     end
 end
@@ -82,42 +78,43 @@ ref_img_pyramid_quad = Float64[
     tfm2 = LinearMap(RotMatrix(-pi/4))
 
     @testset "warp" begin
-        imgr = warp(img_pyramid, tfm1)
+        imgr = warp(img_pyramid, inv(tfm1))
         @test indices(imgr) == (0:6, 0:6)
         @test eltype(imgr) == eltype(img_pyramid)
         # Use map and === because of the NaNs
         @test nearlysame(round.(Float64.(parent(imgr)),3), round.(ref_img_pyramid,3))
 
         @testset "OffsetArray" begin
-            imgr_cntr = warp(img_pyramid_cntr, tfm2)
+            imgr_cntr = warp(img_pyramid_cntr, inv(tfm2))
             @test indices(imgr_cntr) == (-3:3, -3:3)
             @test nearlysame(parent(imgr_cntr), parent(imgr))
         end
 
         @testset "Quadratic Interpolation" begin
             itp = interpolate(img_pyramid_cntr, BSpline(Quadratic(Flat())), OnCell())
-            imgrq_cntr = warp(itp, tfm2)
+            imgrq_cntr = warp(itp, inv(tfm2))
             @test indices(imgrq_cntr) == (-3:3, -3:3)
             @test nearlysame(round.(Float64.(parent(imgrq_cntr)),3), round.(ref_img_pyramid_quad,3))
         end
     end
 
-    @testset "WarpedView" begin
-        imgr = WarpedView(img_pyramid, tfm1)
+    @testset "InvWarpedView" begin
+        imgr = InvWarpedView(img_pyramid, tfm1)
         @test indices(imgr) == (0:6, 0:6)
         # Use map and === because of the NaNs
         @test nearlysame(round.(Float64.(imgr[0:6, 0:6]),3), round.(ref_img_pyramid,3))
 
         @testset "OffsetArray" begin
-            imgr_cntr = WarpedView(img_pyramid_cntr, tfm2)
+            imgr_cntr = InvWarpedView(img_pyramid_cntr, tfm2)
             @test indices(imgr_cntr) == (-3:3, -3:3)
             @test nearlysame(imgr_cntr[indices(imgr_cntr)...], imgr[indices(imgr)...])
         end
 
         @testset "Quadratic Interpolation" begin
             itp = interpolate(img_pyramid_cntr, BSpline(Quadratic(Flat())), OnCell())
-            imgrq_cntr = WarpedView(itp, tfm2)
-            @test summary(imgrq_cntr) == "-3:3×-3:3 WarpedView(interpolate(::OffsetArray{Gray{Float64},2}, BSpline(Quadratic(Flat())), OnCell()), LinearMap([0.707107 0.707107; -0.707107 0.707107])) with element type ColorTypes.Gray{Float64}"
+            imgrq_cntr = InvWarpedView(itp, tfm2)
+            @test parent(imgrq_cntr) === itp
+            @test summary(imgrq_cntr) == "-3:3×-3:3 InvWarpedView(interpolate(::OffsetArray{Gray{Float64},2}, BSpline(Quadratic(Flat())), OnCell()), LinearMap([0.707107 0.707107; -0.707107 0.707107])) with element type ColorTypes.Gray{Float64}"
             @test indices(imgrq_cntr) == (-3:3, -3:3)
             @test nearlysame(round.(Float64.(imgrq_cntr[indices(imgrq_cntr)...]),3), round.(ref_img_pyramid_quad,3))
         end