Consider changing sense of transformation

[timholy]

In #24 we began discussing a possible breaking change: inverting the "sense" in which a transformation is used/applied. Currently pixel values are computed as img[inv(tfm)(x)], but perhaps we should consider switching to img[tfm(x)]. For a rotation, the former can be thought of as rotating the image, whereas the latter is rotating the coordinate axes. Let me acknowledge that I was the one who established the current convention; if we decide that the current system is a mess (I fear it may be), you know who to blame...

At the outset, let's acknowledge that from a standpoint of pure mathematics, either choice is viable. This is because (in my opinion) we should focus on diffeomorphisms and all diffeomorphisms are invertible. However, from a practical standpoint some are easier to invert than others. For example, affine transformations are easy to invert, but diffeomorphisms like this one are less straightforward. The latter might arise from a deformable image registration algorithm, so I think we have to take such seriously.

So, let me make the case for switching. Note here that I'm assuming that the goal of registration is to find a diffeomorphism ϕ such that fixed[x] ≈ moving[ϕ(x)]; if someone has a preferred framework for thinking about registration, now is a really, really good time to mention it.

For the function warp in general you need both the forward and inverse, because (given that framework above)

warped[x] = moving[ϕ(x)]

and consequently warped should be defined for all x such that ϕ(x) is within the domain of moving, i.e., the bounding box of ϕ⁻¹(y) for all y within the domain of moving. Since we need both the forward and the reverse, again it doesn't make a great deal of difference which we choose, with one exception: given that ϕ is a diffeomorphism, we can compute that bounding box using only the edges of moving, and so consequently the inversion only needs to be computed for a small fraction of points. If the inversion is computationally expensive (e.g., invokes a nonlinear solver), then it would be a win to supply ϕ rather than ϕ⁻¹

However, for image registration we are likely to only care about warped over the domain of fixed, in which case we can allocate an empty warped and then call warp!(warped, moving, ϕ). In this case, the behavior should be to iterate over all entries in warped, in which case we don't even need ϕ⁻¹. In this case, imposing the requirement that we can (and do) compute ϕ⁻¹ seems like an unnecessary burden.

[timholy]

Let's also CC @tejus-gupta and @animesh2809, since they're not watching this repo but it may be relevant for them too.

[timholy]

Oops, try again: @annimesh2809

[Evizero]

I agree on all points.

Furthermore I think we should make this breaking change sooner rather than later. (Would probably also need a fix in your blog post?)

I prototyped the changes to warp already in #24. I specifically addressed the problem with fixed and warped by allowing the indices of fixed to be pre-specified (see #24 (comment)). So the mutating version of warp! isn't the only way to avoid ϕ⁻¹.

I think I can propose a complete switch by the end of today (its ~17:00 right now here). This way we may have some code to inform us

[c42f]

What you've written here about ϕ vs ϕ⁻¹ makes a lot of sense Tim.

I've not had a huge amount of experience with image registration itself, but I've seen papers in the area where people compute ϕ as a thin plate spline, for example. In practical situations I expect this can be close to non-invertible. Especially if you're registering an image of a real 3D object in a complex environment - complete with 3D rotation and occlusion - the true warp will fail to be any of invertible, 1:1 or differentiable!

@Evizero is right that this need not stop you, as long as the user specifies the desired output resampling by giving a bounding box in the domain of ϕ rather than the range of ϕ.

Personally, the only thing I find unfortunate about changing this is that I prefer thinking about transformations as moving an object in space and therefore imagining ϕ as the forward warp; warping an image like a bedsheet. Maybe this is just a matter of documentation: if you wanted to keep discussing transformations in the active sense, you could discuss the action of warp in terms of ϕ, but just document that the input of warp must be ϕ⁻¹?

Ah, discussing this stuff brings back some fun memories of old projects :-)

[andyferris]

Requiring the inverse transformation to "warp" an image also seems backwards to me. If one wishes to be artistic, say, then non-invertible transformations might be desirable. I'm thinking something like mirroring an image about it's centre - to achieve this unambiguously, clearly the transformation must be a map from the new coordinate system to the old one.

Here's one from Google:

I always loved that image @c42f !

[timholy]

Since there hasn't yet been disagreement, perhaps we should discuss deprecation strategy. This package may not have a huge number of users yet, but it's been sufficiently advertised by now that there are surely some, and I think that changing the meaning of an argument needs to be done in a way that won't trigger catastrophe.

So, some options:

• introduce a new function name for the new meaning of warp, and deprecate warp itself.
• introduce warp_old and warp_new, don't export either one of them, and export const warp = warp_old. Users can override by defining const warp = ImageTransformations.warp_new before the first usage. See below for more detail.
• have warp(img, tfm) and warp!(dest, img, tfm) issue a warning that can be circumvented by warp(img, Forward(tfm)) (where Forward would be a new dummy wrapper introduced in the next release). Without the Forward, it would still take the inverse of the user-supplied transform.
• keep warp, but decide whether to live in the "new world" or "old world" depending on a file in deps/. Old world issues a deprecation warning.

I'm not very fond of the first option because I think warp is the best name for this function. The second and third will trigger two rounds of depwarns, once for the absence of Forward and later (in 6months-1year) when we deprecate Forward to make it no longer necessary. The last option would be the least intrusive, but runs a substantial risk if there multiple "users" (including packages) that are in different stages of migration.

Here's a demo of how the 2nd option might work:

tim@diva:~$ julia -q
julia> module W
       export w
       wold() = (warn("switch to wnew"); 1)
       wnew() = 2
       const w = wold
       end
W

julia> using W

julia> w()
WARNING: switch to wnew
1

julia> 
tim@diva:~$ julia -q
julia> module W
       export w
       wold() = (warn("switch to wnew"); 1)
       wnew() = 2
       const w = wold
       end
W

julia> using W

julia> const w = W.wnew
wnew (generic function with 1 method)

julia> w()
2

On balance I suspect the Forward wrapper is the safest choice; it ensures that each caller is updated. But it does require the most lines of user code changed.

I should say that I'm traveling and will need to focus on the conference and my own (belated) talk preparations. @Evizero, if you're in a hurry, know that I trust your judgment as much as my own. If the blog post needs a rapid update, I should be able to put some time in by Wednesday, or others are authorized to make changes.

[Evizero]

I am not in a hurry. Just trying to be productive when my time allows it.

How about we temporarily deprecate

@deprecate warp(img, tform) warp(Backward(), img, inv(tform))

warp(::Backward, img, tform) = ...
warp(::Forward, img, tform)  = ...

So for a while we keep the default to backward warping while specifying the forward transform (it is not really "forward warping" because in forward warping we would iterate over the source image, while right now we still iterate over the output image using the inverse of the forward transform)

And at some point in the future we make a switch to default to "proper" backward warping or forward warping. Whatever we decide then

[andyferris]

I think I slightly prefer using Backward() and Forward() as a singleton/trait rather than as a wrapper for a transformation.

That is, unless we can use this to solve the passive vs. active transformation situation within CoordinateTransformations itself - e.g. letting users describe whether a LinearMap is to be interpreted as active or passive. Like I said earlier, I think you are having this issue here because of a design flaw in CoordinateTransformations. Another thing that has come to my mind is to perhaps label transformations with a trait describing whether they are active/passive (forward/backward).

For 2D images warp_canvas would probably describe the new behaviour pretty well (but unfortunately it's not a fantastic description for scientific work on higher dimensional images).

[andyferris]

FYI I created an issue in CoordinateTransformations.

[andyferris]

I've thought about this some more - this package is really young, and my opinion is that it seems simplest to just make a breaking change.

We would need to document quite strongly that warp takes the transformation in the inverse sense (from the new indexing coordinates to the old), because that's what the implementation of warp needs to work with. The Forward and Backward thing could turn out to be a distraction, when it might be more productive to educate the user about the existence of inv(::Transformation) (i.e. warp(img, inv(tfm)) is clearer to me than warp(img, Forward(tfm)) or whatever, and learning that inv works on transformations could be useful...).

[timholy]

As much as I dislike the prospect of two rounds of deprecation warnings, I am not sure it's fair to just change it to its opposite without any mechanism for warning users.

[Evizero]

closed with #24