Image warping with scikit-image and transform.PolynomialTransform

Question

I attach a zip archive with all the files needed to illustrate and reproduce the problem.

(I don't have permissions to upload images yet...)

I have an image (test2.png in the zip archive ) with curved lines.

I try to warp it so the lines are straight. I thought of using scikit-image transform, and in particular transform.PolynomialTransform because the transformation involves high order distortions.

So first I measure the precise position of each line at regular intervals in x to define the input interest points (in the file source_test2.csv). Then I compute the corresponding desired positions, located along a straight line (in the file destination_test2.csv).

The figure correspondence.png shows how it looks like.

Next, I simply call transform.PolynomialTransform() using a polynomial of order 3. It finds a solution, but when I apply it using transform.warp(), the result is crazy, as illustrated in the file Crazy_Warped.png

Anybody can tell what I am doing wrong? I tried polynomial of order 2 without luck... I managed to get a good transformation for a sub-image (the first 400 columns only). Is transform.PolynomialTransform() completely unstable in a case like mine?

Here is the entire code:

import numpy as np
import matplotlib.pyplot as plt
import asciitable
import matplotlib.pylab as pylab
from skimage import io, transform

# read image 
orig=io.imread("test2.png",as_grey=True)
# read tables with reference points and their desired transformed positions
source=asciitable.read("source_test2.csv")
destination=asciitable.read("destination_test2.csv")

# format as numpy.arrays as required by scikit-image
# (need to add 1 because I started to count positions from 0...)
source=np.column_stack((source["x"]+1,source["y"]+1))
destination=np.column_stack((destination["x"]+1,destination["y"]+1))
# Plot
plt.imshow(orig, cmap='gray', interpolation='nearest')
plt.plot(source[:,0],source[:,1],'+r')
plt.plot(destination[:,0],destination[:,1],'+b')
plt.xlim(0,orig.shape[1])
plt.ylim(0,orig.shape[0])

# Compute the transformation
t = transform.PolynomialTransform()
t.estimate(destination,source,3)

# Warping the image
img_warped = transform.warp(orig, t, order=2, mode='constant',cval=float('nan'))

# Show the result
plt.imshow(img_warped, cmap='gray', interpolation='nearest')
plt.plot(source[:,0],source[:,1],'+r')
plt.plot(destination[:,0],destination[:,1],'+b')
plt.xlim(0,img_warped.shape[1])
plt.ylim(0,img_warped.shape[0])
# Save as a file
io.imsave("warped.png",img_warped)

Thanks in advance!

Answer 1

There are a couple of things wrong here, mainly they have to do with coordinate conventions. For example, if we examine the code where you plot the original image, and then put the clicked point on top of it:

plt.imshow(orig, cmap='gray', interpolation='nearest')
plt.plot(source[:,0],source[:,1],'+r')
plt.xlim(0,orig.shape[1])
plt.ylim(0,orig.shape[0])

(I've taken out the destination points to make it cleaner) then we get the following image:

As you can see, the y-axis is flipped, if we invert the y-axis with:

source[:,1] = orig.shape[0] - source[:,1]

before plotting, then we get the following:

So that is the first problem (don't forget to invert the destination points as well), the second has to do with the transform itself:

t.estimate(destination,source,3)

From the documentation we see that the call takes the source points first, then the destination points. So the order of those arguments should be flipped.

Lastly, the clicked points are of the form (x,y), but the image is stored as (y,x), so we have to transpose the image before applying the transform and then transpose back again:

img_warped = transform.warp(orig.transpose(), t, order=2, mode='constant',cval=float('nan'))
img_warped = img_warped.transpose()

When you make these changes, you get the following warped image:

These lines aren't perfectly flat but it makes much more sense.

Answer 2

Thank you very much for the detailed answer! I cannot believe I did not see the axis inversion problem... Thanks for catching it! But I am afraid your final solution does not solve my problem... The image you get is still crazy. It should be continuous, no have such big holes and weird distortions... (see final solution below)

I found I could get a reasonable solution using RANSAC:

from skimage.measure import ransac
t, inliers = ransac((destination,source), transform.PolynomialTransform, min_samples=20,residual_threshold=1.0, max_trials=1000)
outliers = inliers == False

I then get the following result

Note that I think I was right using (destination,source) in that order! I think it has to do with the fact that transform.warp requires the inverse_map as input for the transformation object, not the forward map. But maybe I am wrong? The good result I am getting suggest it's correct.

I guess that Polynomial transforms are too unstable, and using RANSAC allows to get a reasonable solution. My problem was then to find a way to change the polynomial order in the RANSAC call... transform.PolynomialTransform() does not take any parameters, and uses by default a 2nd order polynomial, but from the result I can see I would need a 3rd or 4th order polynomial.

So I opened a new question, and got a solution from Stefan van der Walt. Follow the link to see how to do it.

Thanks again for your help!