Add thin-plate splines warping

[decorouz]

Description

Tracks progress for the TPS Outreachy project.

cc @ana42742, @mkcor

Checklist

• A descriptive but concise pull request title
• Docstrings for all functions
• Unit tests
• A gallery example in ./doc/examples for new features
• Contribution guide is followed

[lagru]

Great progress and work! I've added a few comments to keep you busy. ;)

[mkcor]

I wonder how fitting this title is... When I read "Interpolate images (...)," I expect that images with missing or hidden regions would be reconstructed with a given interpolation method, along the lines of the cornea image restoration example. Instead, at first glance, this example looks more like "Deform an image with thin-plate splines" to me. What do you think? /cc @lagru

[lagru]

This might get into semantics but I think interpolation just means something like finding new values based on existing often neighboring data. That's not necessarily reconstruction only but also deformation / remapping, / relocating. I don't find "deform" quite fitting if TPS can also be used to shift or rotate an image. So I'm okay with "interpolate" being the generalized term. Whether this gallery example demonstrates all use cases of "interpolation" is another thing...

[decorouz]

The example demonstrate interpolating a 2D image over a set src and dst points.
@mkcor I update this section based on #7040 (comment)

[stefanv]

You cannot technically "interpolate images", so how about something like "Image warping using thin-plate splines"?

[mkcor]

Following the convention of writing titles at the imperative mood, I'd suggest either "Warp images using thin-plate splines" or "Use thin-plate splines to warp images" (or "Use thin-plate splines for image warping?").

[mkcor]

I'm happy to improve the writing of the gallery example in subsequent PRs.

As suggested by @lagru in yesterday's community call, we could move this new transform into the future submodule, so we don't commit to specific API choices for now, but users can start using it (and possibly provide feedback). @stefanv what do you think?

[lagru]

@decorouz do you mind if I move this to the skimage.future module or do you want to give it a try? :)

[decorouz]

@decorouz do you mind if I move this to the skimage.future module or do you want to give it a try? :)

Please go ahead and move it.

[lagru]

The CI failure is due to a warning while building the docs on Windows. Not sure why that's showing up, but it seems completely unrelated to this.

[lagru]

@decorouz, @mkcor I made another pass at the gallery example and docstrings. I've put it through https://languagetool.org and also tried to improve the clearness in a few places, especially the descriptions of tps_warp's interpolation_order parameter. I think the previous version was incorrect.

[mkcor]

I've left only minor suggestions! 🚀

[imagesc-bot]

This pull request has been mentioned on Image.sc Forum. There might be relevant details there:

https://forum.image.sc/t/equivalent-for-matlabs-piecewiselineartransformation2d/51035/12

[stefanv]

Thanks for all the hard work @decorouz, and @lagru and @mkcor for guidance.

[stefanv]

You cannot technically "interpolate images", so how about something like "Image warping using thin-plate splines"?

[stefanv]

The phrase interpolation is rather confusing here, since we do two different types of interpolation during image warping. (1) Interpolation of the warp model and (2) interpolation of pixel values. The first is often not even called interpolation, so I'd avoid the term in that context.

[stefanv]

Why can we not fit TPS within the existing transformation framework?
Estimate model + warp should be sufficient, and better not to have multiple ways of doing the same thing.

[lagru]

Passing TpsTransform() to warp already seems to work. Though, I notice that one needs to switch src and dst to receive the same results:

import numpy as np
import skimage as ski

image = np.array(
    [[9, 9, 9, 9, 9, 9, 9],
     [9, 0, 0, 0, 0, 0, 9],
     [9, 0, 1, 1, 1, 0, 9],
     [9, 0, 1, 2, 1, 0, 9],
     [9, 0, 1, 1, 1, 0, 9],
     [9, 0, 0, 0, 0, 0, 9],
     [9, 9, 9, 9, 9, 9, 9]],
    dtype=float,
)
# Zoom in
src = np.array([[2, 2], [2, 4], [4, 4], [4, 2]])
dst = np.array([[1, 1], [1, 5], [5, 5], [5, 1]])

tform = ski.future.TpsTransform()
tform.estimate(src=dst, dst=src)  # switched src & dst!
using_warp = ski.transform.warp(image, tform)

using_tps_warp = ski.future.tps_warp(image, src=src, dst=dst)
np.testing.assert_array_equal(using_warp, using_tps_warp)

However, tps_warp provides grid_scaling as an optional performance optimization at the cost of accuracy. Would you recommend integrating that into ski.transform.warps? It's an optimization we took from the reference implementation...

[decorouz]

Also, the need to switch src and dst we took from the reference implementation
https://github.com/zpincus/celltool/blob/045920089f88a445fa2483500ef867d267ecef13/celltool/numerics/image_warp.py#L41-L43

# make the reverse transform warping from the to_points to the from_points, because we
# do image interpolation in this reverse fashion
transform = _make_warp(to_points, from_points, x, y)

[stefanv]

We can make scaling a property of the ThinPlateSpline transform.

[stefanv]

Can we get rid of this method? warp should already provide its functionality.

[lagru]

See #7040 (comment) above on the parameter grid_scaling.

[lagru]

From the same discussion: let's go with warp instead of adding tpw_warp. We can look into whether we can integrate the grid_scaling optimization in warp itself but it is not that essential to memory optimization (@stefanv can explain that better than me).

I'm not sure when I will get time to look into implementing that right now. @decorouz, feel welcome give it a try if you have time and incentive. :D

[stefanv]

For each pixel in the output image, we need to find out where in the input image to "fetch" it from. We can do that by computing a coordinate transformation for every pixel. That can be done per pixel (very low memory consumption), but in Python we typically vectorize it, i.e. use a bunch of memory to store the coordinate map (output coord -> input coord).

I haven't looked at this carefully, but I suspect the optimization used here is to generate a smaller coordinate map, and to then trick map_coordinate into interpolating the coordinates as needed. Less accurated, but smaller memory consumption, and perhaps, due to caching, a little bit faster.

[stefanv]

I think this information is too scant.

I'd start with something along the lines of:

To warp an image, we start with a set of source and target coordinates.
The goal is to deform the image such that the source points move to the target locations.
Typically, we only know the target positions for a few, select source points.
To calculate the target positions for all other pixel positions, we need a model.
Various such models exist, such as affine or projective transformations.

At this point, we can refer the reader to https://scikit-image.org/docs/stable/auto_examples/transform/plot_transform_types.html

Most transformations are linear (i.e., they preserve straight lines), but sometimes
we need more flexibility. One model that represents a non-linear transformation,
i.e. one where lines can bend, is thin-plate splines.

Thin-plate splines draw on the analogy of a metal sheet, which has inherent rigidity.
Consider our source points: each has to move a certain distance, in both the x and y directions, to land in their corresponding target positions.
First, examine only the x coordinates.
Imagine placing a thin metal plate on top of the image.
Now bend it, such that at each source point, the plate's z offset is the distance,
positive or negative, that that source point has to travel in the x direction
in order to land in its target position.
The plate resists bending, and therefore remains smooth.
We can read offsets for coordinates other than source points from the position of the plate.
The same procedure can be repeated for the y coordinates.

This gives us our thin-plate spline model that maps any (x, y) coordinate to a target position.

We should also refer the reader to an explanation of how that warping is done (output to input image, using e.g. bilinear interpolation). @jni I recall we wrote up a description of that somewhere; Elegant SciPy? skimage tutorials?

[stefanv]

Explicitly inverting matrices is seldom a good idea. Use linalg.solve?

[stefanv]

Also, perhaps better not to reinvent the wheel and use SciPy's RBF to do the work?
Could reduce the code here by quite a bit, if it works.

[lagru]

I remember looking at it in the past, but I don't remember why we didn't follow that approach. Totally possible, that we didn't yet have a good grasp on concepts. @decorouz, do you want to give it a shot?

[lagru]

To add some context from a discussion @stefanv and me just had: linalg.solve is usually more precise in calculating the inverse, while linalg.inv is prone to small numeric errors that accumulate during the calculation. (At least that's what I understood 😃)

[stefanv]

And faster :)

The short mathematical answer.
The longer mathematical unpacking.

[stefanv]

What are "norm distance" and "interlandmark"? Consider more common terminology here.