Streamline API for warp and WarpedView #24
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,49 @@ | ||
| # FIXME: upstream https://github.com/JuliaGraphics/ColorVectorSpace.jl/issues/75 | ||
| @inline _nan(::Type{HSV{Float16}}) = HSV{Float16}(NaN16,NaN16,NaN16) | ||
| @inline _nan(::Type{HSV{Float32}}) = HSV{Float32}(NaN32,NaN32,NaN32) | ||
| @inline _nan(::Type{HSV{Float64}}) = HSV{Float64}(NaN,NaN,NaN) | ||
| @inline _nan{T}(::Type{T}) = nan(T) | ||
|
|
||
| # 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 | ||
|
|
||
| function box_extrapolation(itp::AbstractInterpolation, degree::Union{Linear,Constant}, args...) | ||
| throw(ArgumentError("Boxing an interpolation in another interpolation is discouraged. Did you specify the parameter \"$degree\" on purpose?")) | ||
|
|
||
| end | ||
|
|
||
| function box_extrapolation(itp::AbstractExtrapolation, fill::FillType) | ||
| throw(ArgumentError("Boxing an extrapolation in another extrapolation is discouraged. Did you specify the parameter \"$fill\" on purpose?")) | ||
| end | ||
|
|
||
| # This is type-piracy, but necessary if we want Interpolations to be | ||
| # independent of OffsetArrays. | ||
| function AxisAlgorithms.A_ldiv_B_md!(dest::OffsetArray, F, src::OffsetArray, dim::Integer, b::AbstractVector) | ||
| indsdim = indices(parent(src), dim) | ||
| indsF = indices(F)[2] | ||
| if indsF == indsdim | ||
| AxisAlgorithms.A_ldiv_B_md!(parent(dest), F, parent(src), dim, b) | ||
| return dest | ||
| end | ||
| throw(DimensionMismatch("indices $(indices(parent(src))) do not match $(indices(F))")) | ||
| end | ||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,74 @@ | ||
| immutable InvWarpedView{T,N,A<:AbstractArray,F1<:Transformation,I,F2<:Transformation,E<:AbstractExtrapolation} <: AbstractArray{T,N} | ||
|
There was a problem hiding this comment. At some point we might want to define a There was a problem hiding this comment. good idea. I will think on that |
||
| 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 | ||
|
There was a problem hiding this comment. Obviously, this will be important to get right. |
||
| """ | ||
| @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 | ||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,51 +1,35 @@ | ||
| """ | ||
| warp(img, tform, [indices], [degree = Linear()], [fill = NaN]) -> imgw | ||
|
|
||
| Transform the coordinates of `img`, returning a new `imgw` satisfying | ||
| `imgw[I] = img[tform(I)]`. This approach is known as backward | ||
| mode warping. The transformation `tform` should be defined using | ||
| [CoordinateTransformations.jl](https://github.com/FugroRoames/CoordinateTransformations.jl). | ||
|
There was a problem hiding this comment. Not strictly required that it come from CoordinateTransformations, is it? Anything that can take an There was a problem hiding this comment. I'll rewrite it to be more neutral |
||
|
|
||
| # Interpolation scheme | ||
|
|
||
| At off-grid points, `imgw` is calculated by interpolation. The | ||
|
There was a problem hiding this comment. I suggest abandoning the idea that "on grid points" [edit: or, conversely off-grid points] are particularly meaningful. This is only an artifact of certain interpolation schemes. From a camera, a pixel value is not a gridded point sample of some underlying continuous image. (Neither is it a little square "pixel box" we've all grown so used to looking at - pixel boxes are just the zeroth order b-spline interpolant.) In reality, the image is the incoming light field convolved with the point spread function of the optics, and integrated over the physical CCD/cmos sensor pixel. The pixels are then a discrete low dimensional representation of the underlying light field. All of which is to say... I'm not sure what! Perhaps just that we shouldn't confuse the standard (ie, cheap) interpolation schemes with the interpretation of the pixels values. There was a problem hiding this comment. Alternatives:
Other ideas? There was a problem hiding this comment. It sounds like I might be missing something about the way images are treated in this package, or about the way that interpolations works in conjunction with this. Normally (ie, in all other software packages I know about), cubic spline image resampling technically gives an approximating spline rather than an interpolating one. Is that not the case here? There was a problem hiding this comment. For quadratic and cubic, Interpolations solves a tridiagonal/Woodbury system for the coefficients to ensure that they are interpolating. There was a problem hiding this comment. To address @c42f's point more completely, yes, agreed 100% with your physical model of the process. The issue is that without more information there isn't, to my knowledge, a great way to do a lot better than interpolation. So personally I'm fine with talking about a "grid," though I'd also be fine with a slightly more verbose statement "if There was a problem hiding this comment. Thanks for clearing that up. I'm also not sure what more should be done, other than some kind of interpolating or approximating spline. (Or at least, that anything else which could be done is far too fancy to be done behind the user's back.) This makes more sense of @Evizero's worry about using higher-than-linear orders inside I must admit it's not quite clear to me that an interpolating spline is the right choice in general over the much simpler and - at least somewhat - more efficient approximating splines for image resampling. In the alternative kernel-based view of reconstructing a continuous image from discrete pixel values, there's a convenient family of bicubic kernels parameterizing the apparent sharpness of the resulting reconstruction. This formulation makes the reconstruction entirely spatially local, and there's no need to run any kind of solver before sampling in an arbitrary location. Perhaps all I'm quibbling about is that the docs give the idea that the computation is different on vs off the integer lattice. It shouldn't be, other than as an optimization in the case that the reconstruction is done with an actual interpolating spline. There was a problem hiding this comment.
@c42f, a reference would be awesome.
That's absolutely right. Do you have a good way to phrase this? I wondered if that's what you were worried about, but I couldn't quickly come up with a concise way of expressing that. There was a problem hiding this comment. It's mentioned in passing on the wikipedia page about bicubic interpolation, but it's not great. I finally found a better reference in some old code of mine:
The code itself was just an implementation of the piecewise formulae from the paper. There was a problem hiding this comment.
Not exactly :-) How about Obviously I'm trying to avoid the word "interpolation", as there's quite reasonable non-interpolating reconstructions. There was a problem hiding this comment. Musing a bit more about reconstructions and the desiderata from the Mitchell paper -
|
||
| degree of the b-spline can be specified with the optional | ||
| parameter `degree`, which can take the values `Linear()` or | ||
| `Constant()`. | ||
|
There was a problem hiding this comment. This raises an interesting issue: since There was a problem hiding this comment. I was aiming for for consistency with the view. Regardless of consistency, I do think that we shouldn't hide such allocating actions from a user. It should be a choice. It isn't that much more effort to pass The options presented are special in that they are quite conservative in the memory allocation regard. There was a problem hiding this comment. I agree with Tim here, I think it would be very natural to support other interpolation types. (Caveat: I don't think I currently have a good overview of what you're building in this package as a whole.) There was a problem hiding this comment. Having had the long discussion about "off-grid points" I think I now understand the context for your reservations about higher orders, and the view pledge which Tim is referring to. But I also think we could do a cubic spline view just fine, as long as we take a convolution kernel-style viewpoint of the reconstruction (interpolation/approximation) process: in this variant, there's no matrix solver to run after modifying the underlying image, you only need the pixel values from a local 4x4 image patch around the sample point, and the fractional coordinates of the sample point away from the integer lattice. There was a problem hiding this comment. Sure. But we don't have that implemented yet, so let's finish this so that it works for what we have. I can go either way (restrict to There was a problem hiding this comment. Sure, there's lots of great stuff here, I don't mean to hold up getting it merged. |
||
|
|
||
| The b-spline interpolation is used when `img` is a plain array | ||
| and `img[tform(I)]` is inbound. In the case `tform(I)` maps | ||
| to indices outside the original `img`, the value or extrapolation | ||
|
There was a problem hiding this comment. might be clearer to say "those locations are set to a value |
||
| scheme denoted by the optional parameter `fill` (which defaults | ||
| to `NaN` if the element type supports it, and `0` otherwise) is | ||
| used to indicate locations for which `tform(I)` was outside the | ||
| bounds of the input `img`. | ||
|
|
||
| For more control over the interpolation scheme --- and how | ||
| beyond-the-edge points are handled --- pass it in as an | ||
|
|
||
| `AbstractInterpolation` or `AbstractExtrapolation` from | ||
| [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl). | ||
|
There was a problem hiding this comment. It's interesting how these packages fit together. It might be worth noting that For severe warping and large downsampling things can get quite interesting, as you need full anisotropic filtering, usually coupled with image pyramiding to make it efficient. Arguably, There was a problem hiding this comment. (Obviously, this is a side comment - nothing to do with any changes you should be making here in this PR.) There was a problem hiding this comment. We currently have the same kind of problem with see #7 (comment) |
||
|
|
||
| # The meaning of the coordinates | ||
|
|
||
| The output array `imgw` has indices that would result from applying | ||
| `inv(tform)` to the indices of `img`. This can be very handy for keeping | ||
| track of how pixels in `imgw` line up with pixels in `img`. | ||
|
|
||
| If you just want a plain array, you can "strip" the custom indices | ||
|
|
@@ -65,7 +49,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, inv(tfm)); | ||
|
There was a problem hiding this comment. I think it would be better to pass |
||
|
|
||
| julia> indices(imgw) | ||
| (-196:709,-68:837) | ||
|
|
@@ -76,7 +60,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, inv(rot)); | ||
|
|
||
| julia> indices(imgw) | ||
| (-293:612,-293:611) | ||
|
|
@@ -87,28 +71,117 @@ julia> indices(imgr) | |
| (Base.OneTo(906),Base.OneTo(905)) | ||
| ``` | ||
| """ | ||
| function warp_new{T}(img::AbstractExtrapolation{T}, tform, inds::Tuple = autorange(img, inv(tform))) | ||
|
There was a problem hiding this comment. Someday this will have to become There was a problem hiding this comment. Presumably invertibility is more of a trait? I can't see it fitting neatly into the type hierarchy. |
||
| out = OffsetArray(Array{T}(map(length, inds)), inds) | ||
| warp!(out, img, tform) | ||
| end | ||
|
|
||
| # this function was never exported, so no need to deprecate | ||
|
|
||
| function warp!(out, img::AbstractExtrapolation, tform) | ||
| @inbounds for I in CartesianRange(indices(out)) | ||
| out[I] = _getindex(img, tform(SVector(I.I))) | ||
| end | ||
| out | ||
| end | ||
|
|
||
| function warp_new(img::AbstractArray, tform, inds::Tuple, args...) | ||
| etp = box_extrapolation(img, args...) | ||
| warp_new(etp, tform, inds) | ||
| end | ||
|
|
||
| function warp_new(img::AbstractArray, tform, args...) | ||
| etp = box_extrapolation(img, args...) | ||
| warp_new(etp, tform) | ||
| end | ||
|
|
||
| # # after deprecation period: | ||
| # @deprecate warp_new(img::AbstractArray, tform, args...) warp(img, tform, args...) | ||
| # @deprecate warp_old(img::AbstractArray, tform, args...) warp(img, tform, args...) | ||
|
|
||
| """ | ||
| warp(img, tform, [indices], [degree = Linear()], [fill = NaN]) -> imgw | ||
|
|
||
| !!! warning "Deprecation Warning!" | ||
|
|
||
| This method with the signature `warp(img, tform, args...)` is | ||
|
There was a problem hiding this comment. I'm not sure we need to document the Or we can specify how one switches in the docstring (e.g., to help people who run into this at a time when they don't happen to have internet access, e.g., while flying). |
||
| deprecated in favour of the cleaner interpretation `warp(img, | ||
| inv(tform), args...)`. Set `const warp = | ||
| ImageTransformations.warp_new` right after package import to | ||
| change to the new behaviour right away. | ||
|
|
||
| Transform the coordinates of `img`, returning a new `imgw` | ||
| satisfying `imgw[I] = img[inv(tform(I))]`. This approach is known | ||
| as backward mode warping. The transformation `tform` should be | ||
| defined using | ||
| [CoordinateTransformations.jl](https://github.com/FugroRoames/CoordinateTransformations.jl). | ||
|
|
||
| # Interpolation scheme | ||
|
|
||
| At off-grid points, `imgw` is calculated by interpolation. The | ||
| degree of the b-spline can be specified with the optional | ||
| parameter `degree`, which can take the values `Linear()` or | ||
| `Constant()`. | ||
|
|
||
| The b-spline interpolation is used when `img` is a plain array | ||
| and `img[inv(tform(I))]` is inbound. In the case `inv(tform(I))` | ||
| maps to indices outside the original `img`, the value or | ||
| extrapolation scheme denoted by the optional parameter `fill` | ||
| (which defaults to `NaN` if the element type supports it, and `0` | ||
| otherwise) is used to indicate locations for which | ||
| `inv(tform(I))` was outside the bounds of the input `img`. | ||
|
|
||
| For more control over the interpolation scheme --- and how | ||
| beyond-the-edge points are handled --- pass it in as an | ||
| `AbstractInterpolation` or `AbstractExtrapolation` from | ||
| [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl). | ||
|
|
||
| # The meaning of the coordinates | ||
|
|
||
| The output array `imgw` has indices that would result from applying | ||
| `tform` to the indices of `img`. This can be very handy for keeping | ||
| track of how pixels in `imgw` line up with pixels in `img`. | ||
|
|
||
| If you just want a plain array, you can "strip" the custom indices | ||
| with `parent(imgw)`. | ||
|
|
||
| # Examples: a 2d rotation (see JuliaImages documentation for pictures) | ||
|
|
||
| ```jldoctest | ||
| julia> using Images, CoordinateTransformations, TestImages, OffsetArrays | ||
|
|
||
| julia> img = testimage("lighthouse"); | ||
|
|
||
| julia> indices(img) | ||
| (Base.OneTo(512),Base.OneTo(768)) | ||
|
|
||
| # Rotate around the center of `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> indices(imgw) | ||
| (-196:709,-68:837) | ||
|
|
||
| # Alternatively, specify the origin in the image itself | ||
| 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> indices(imgw) | ||
| (-293:612,-293:611) | ||
|
|
||
| julia> imgr = parent(imgw); | ||
|
|
||
| julia> indices(imgr) | ||
| (Base.OneTo(906),Base.OneTo(905)) | ||
| ``` | ||
| """ | ||
| @generated function warp_old(img::AbstractArray, tform, args...) | ||
|
There was a problem hiding this comment. Does this need to be a There was a problem hiding this comment. The reason for |
||
| warn("'warp(img, tform)' is deprecated in favour of the cleaner interpretation 'warp(img, inv(tform))'. Set 'const warp = ImageTransformations.warp_new' right after package import to change to the new behaviour right away.") | ||
|
There was a problem hiding this comment. Shall we use the The web site can provide more justification, e.g., like my post in #25 (comment) except probably using embedded pictures. There was a problem hiding this comment. Not sure what you are referring to with the comment on embedded pictures. |
||
| :(warp_new(img, inv(tform), args...)) | ||
| end | ||
| const warp = warp_old | ||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What does the
box_signify?There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just that the image is "boxed" or "decorated" by an extrapolation. It's supposed to signify that this will be a view.
Can you think of a better name?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I was tempted by
wrap_extrapolation, but I think that could be confusing: people might think it's a typo ofwarp! So let's stay withbox. (EDIT: originally I had a typo in the word typo. Amusing!)