The method implemented is from Goshtasby (1986), “Piecewise linear mapping functions for image registration”. In Section 2.4 (Summary), they outline three steps of which skimage implements only the first two. The third step is:

Determine functions for mapping of points outside the convex hulls by extending the boundary triangle planes. Intersections of neighboring planes limit the extent of the boundary planes that should be used in the extrapolation.

It is not immediately obvious to me how to do this.

I played around with two approaches: a 4th-degree polynomial transform, and thin plate splines. The first worked OK, the second was great. The implementation is by Zachary Pincus, who is an emeritus skimage developer (he said that he’d consider releasing the code under BSD if someone found it useful).

I’d like to see thin plate splines in scikit-image. It has been requested as early as 2017 already.

```
# https://github.com/zpincus/celltool/blob/9f2da1a905413fedcca89ad71ea94c51b1f2da18/celltool/numerics/image_warp.py
# https://github.com/scikit-image/scikit-image/issues/2429#issuecomment-270493362
from image_warp import warp_images
import matplotlib.pyplot as plt
import numpy as np
# Load data
grid_picture = np.load('grid_picture.npy')
original_points = np.load('original_points.npy')
transformed_points = np.load('transformed_points.npy')
grid_picture_unwarped = warp_images(
original_points[:, ::-1],
transformed_points[:, ::-1],
[grid_picture],
[0, 0, grid_picture.shape[1],
grid_picture.shape[0]]
)[0]
plt.imshow(grid_picture_unwarped)
plt.scatter(transformed_points[:, 0], transformed_points[:, 1], s=2)
plt.show()
```

P.S. Please note the order of coordinates. Our `transform`

module relies on coordinates being in XY instead of row-column format. This was an unfortunate decision on my part, which we are working on remedying.