Piecewise affine transformation + strain output looks weird

Question

I have an image that I am trying to use using skimage.PiecewiseAffineTransformand skimage.warp. I have a set of breakpoints ( true) mapped to a new set of breakpoints ( ideal), but the warp does not return the expected.

In this example, I have a simple wavelength gradient that I am trying to “straighten” into columns. (You may ask why I find outlines and interpolation, but that is because I actually apply this code to a more complex use case. I just wanted to reproduce all the code for this simple example, which leads to the same weird conclusion.)

Is there any reason my output image has an input image distorted into a square and an insert? I am using Python 2.7.12 and matplotlib 1.5.1. Here is the code.

import matplotlib.pyplot as plt
import numpy as np
from skimage import measure, transform

true = np.array([range(i,i+10) for i in range(20)])
ideal = np.array([range(10)]*20)

# Find contours of ideal and true images and create list of control points for warp
true_control_pts = []
ideal_control_pts = []

for lam in ideal[0]:
    try:
        # Get the isowavelength contour in the true and ideal images
        tc = measure.find_contours(true, lam)[0]
        ic = measure.find_contours(ideal, lam)[0]
        nc = np.ones(ic.shape)

        # Use the y coordinates of the ideal contour
        nc[:, 0] = ic[:, 0]

        # Interpolate true contour onto ideal contour y axis so there are the same number of points
        nc[:, 1] = np.interp(ic[:, 0], tc[:, 0], tc[:, 1])

        # Add the control points to the appropriate list
        true_control_pts.append(nc.tolist())
        ideal_control_pts.append(ic.tolist())

    except (IndexError,AttributeError):
        pass

true_control_pts = np.array(true_control_pts)
ideal_control_pts = np.array(ideal_control_pts)

length = len(true_control_pts.flatten())/2
true_control_pts = true_control_pts.reshape(length,2)
ideal_control_pts = ideal_control_pts.reshape(length,2)

# Plot the original image
image = np.array([range(i,i+10) for i in range(20)]).astype(np.int32)
plt.figure()
plt.imshow(image, origin='lower', interpolation='none')
plt.title('Input image')

# Warp the actual image given the transformation between the true and ideal wavelength maps
tform = transform.PiecewiseAffineTransform()
tform.estimate(true_control_pts, ideal_control_pts)
out = transform.warp(image, tform)

# Plot the warped image!
fig, ax = plt.subplots()
ax.imshow(out, origin='lower', interpolation='none')
plt.title('Should be parallel lines')

The result for this looks like this:

Any help would be greatly appreciated!

Answer 1

, . , . 45- , , , , . , , /* */ ( , Python, :)).

, ! .

:

• interp
• , (, )
• ( 0 - 1e-9).

, , - 3D-, , "" "" . , , , . interp.

, , "before" "after" ideal true:). , - .

import pdb
import matplotlib.pyplot as plt
import numpy as np
from skimage import measure, transform, img_as_float

from mpl_toolkits.mplot3d import Axes3D # /**/

#/**/
# From https://stackoverflow.com/a/14491059/2877364 by
# https://stackoverflow.com/users/1355221/dansalmo
def flatten_list(L):
    for item in L:
        try:
            for i in flatten_list(item): yield i
        except TypeError:
            yield item
#end flatten_list

true_input = np.array([range(i,i+10) for i in range(20)])  # /** != True **/
ideal = np.array([range(10)]*20)

#pdb.set_trace()
# Find contours of ideal and true_input images and create list of control points for warp
true_control_pts = []
ideal_control_pts = []
OLD_true=[]     # /**/ for debugging
OLD_ideal=[]

for lam in [x+0.5 for x in ideal[0]]:   # I tried the 0.5 offset just to see,
    try:                                # but it didn't make much difference

        # Get the isowavelength contour in the true_input and ideal images
        tc = measure.find_contours(true_input, lam)[0]
        ic = measure.find_contours(ideal, lam)[0]
        nc = np.zeros(ic.shape) # /** don't need ones() **/

        # Use the y /** X? **/ coordinates of the ideal contour
        nc[:, 0] = ic[:, 0]

        # Interpolate true contour onto ideal contour y axis so there are the same number of points

        # /** Have to sort first - https://docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html#numpy-interp **/
        tc_sorted = tc[tc[:,0].argsort()]
            # /** Thanks to https://stackoverflow.com/a/2828121/2877364 by
            # https://stackoverflow.com/users/208339/steve-tjoa **/

        nc[:, 1] = np.interp(ic[:, 0], tc_sorted[:, 0], tc_sorted[:, 1],
            left=np.nan, right=np.nan)
            # /** nan: If the interpolation is out of range, we're not getting
            #     useful data.  Therefore, flag it with a nan. **/

        # /** Filter out the NaNs **/
        # Thanks to https://stackoverflow.com/a/11453235/2877364 by
        # https://stackoverflow.com/users/449449/eumiro
        #pdb.set_trace()
        indices = ~np.isnan(nc).any(axis=1)
        nc_nonan = nc[indices]
        ic_nonan = ic[indices]

        # Add the control points to the appropriate list.
        # /** Flattening here since otherwise I wound up with dtype=object
        #     in the numpy arrays. **/
        true_control_pts.append(nc_nonan.flatten().tolist())
        ideal_control_pts.append(ic_nonan.flatten().tolist())

        OLD_true.append(nc)     # /** for debug **/
        OLD_ideal.append(ic)

    except (IndexError,AttributeError):
        pass

#pdb.set_trace()
# /** Make vectors of all the control points. **/
true_flat = list(flatten_list(true_control_pts))
ideal_flat = list(flatten_list(ideal_control_pts))
true_control_pts = np.array(true_flat)
ideal_control_pts = np.array(ideal_flat)

# Make the vectors 2d
length = len(true_control_pts)/2
true_control_pts = true_control_pts.reshape(length,2)
ideal_control_pts = ideal_control_pts.reshape(length,2)

#pdb.set_trace()

# Plot the original image
image = np.array([range(i,i+10) for i in range(20)]) / 30.0 # /**.astype(np.int32)**/
    # /** You don't want int32 images!  See
    #     http://scikit-image.org/docs/dev/user_guide/data_types.html .
    #     Manually rescale the image to [0.0,1.0] by dividing by 30. **/
image_float = img_as_float(image) #/** make sure skimage is happy */ 
fig = plt.figure()
plt.imshow(image_float, origin='lower', interpolation='none')
plt.title('Input image')
fig.show()  # /** I needed this on my test system **/

# Warp the actual image given the transformation between the true and ideal wavelength maps
tform = transform.PiecewiseAffineTransform()
tform.estimate(true_control_pts, ideal_control_pts)
out = transform.warp(image, tform)
    # /** since we started with float, and this is float, too, the two are
    #     comparable. **/

pdb.set_trace()

# Plot the warped image!
fig, ax = plt.subplots()
ax.imshow(out, origin='lower', interpolation='none')    # /**note: float**/
plt.title('Should be parallel lines')
fig.show()

# /** Show the control points.
#     The z=0 plane will be the "true" control points (before), and the
#     z=1 plane will be the "ideal" control points (after). **/
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
fig.show()

for rowidx in range(length):
    ax.plot([true_control_pts[rowidx,0], ideal_control_pts[rowidx,0]],
            [true_control_pts[rowidx,1], ideal_control_pts[rowidx,1]],
            [0,1])

input() # /** because I was running from the command line **/

:

:

, , .

import pdb
import matplotlib.pyplot as plt
import numpy as np
from skimage import measure, transform, img_as_float

from mpl_toolkits.mplot3d import Axes3D # /**/

#/**/
# From https://stackoverflow.com/a/14491059/2877364 by
# https://stackoverflow.com/users/1355221/dansalmo
def flatten_list(L):
    for item in L:
        try:
            for i in flatten_list(item): yield i
        except TypeError:
            yield item
#end flatten_list

#/**/
# Modified from /questions/490589/how-to-generate-equispaced-interpolating-values/2094239#2094239 by
# https://stackoverflow.com/users/2588210/christian-k
def equispace(data, npts):
    x,y = data.T
    xd =np.diff(x)
    yd = np.diff(y)
    dist = np.sqrt(xd**2+yd**2)
    u = np.cumsum(dist)
    u = np.hstack([[0],u])

    t = np.linspace(0,u.max(),npts)
    xn = np.interp(t, u, x)
    yn = np.interp(t, u, y)
    return np.column_stack((xn, yn))

true_input = np.array([range(i,i+10) for i in range(20)])  # /** != True **/
ideal = np.array([range(10)]*20)

#pdb.set_trace()
# Find contours of ideal and true_input images and create list of control points for warp
true_control_pts = []
ideal_control_pts = []
OLD_true=[]     # /**/ for debugging
OLD_ideal=[]

for lam in [x+0.5 for x in ideal[0]]:   # I tried the 0.5 offset just to see,
    try:                                # but it didn't make much difference

        # Get the isowavelength contour in the true_input and ideal images
        tc = measure.find_contours(true_input, lam)[0]
            # /** So this might not have very many numbers in it. **/
        ic = measure.find_contours(ideal, lam)[0]
            # /** CAUTION: this is assuming the contours are going the same
            #       direction.  If not, you'll need to make it so. **/
        #nc = np.zeros(ic.shape) # /** don't need ones() **/

        # /** We just want to find points on _tc_ to match _ic_.  That's
        #       interpolation _within_ a curve. **/
        #pdb.set_trace()
        nc_by_t = equispace(tc,ic.shape[0])
        ic_by_t = equispace(ic,ic.shape[0])


        ### /** Not this **/
        ## Use the y /** X? **/ coordinates of the ideal contour
        #nc[:, 0] = ic[:, 0]
        #
        ## Interpolate true contour onto ideal contour y axis so there are the same number of points
        #
        ## /** Have to sort first - https://docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html#numpy-interp **/
        #tc_sorted = tc[tc[:,0].argsort()]
        #    # /** Thanks to https://stackoverflow.com/a/2828121/2877364 by
        #    # https://stackoverflow.com/users/208339/steve-tjoa **/
        #
        #nc[:, 1] = np.interp(ic[:, 0], tc_sorted[:, 0], tc_sorted[:, 1],
        #    left=np.nan, right=np.nan)
        #    # /** nan: If the interpolation is out of range, we're not getting
        #    #     useful data.  Therefore, flag it with a nan. **/

        # /** Filter out the NaNs **/
        # Thanks to https://stackoverflow.com/a/11453235/2877364 by
        # https://stackoverflow.com/users/449449/eumiro
        #pdb.set_trace()
        #indices = ~np.isnan(nc).any(axis=1)
        #nc_nonan = nc[indices]
        #ic_nonan = ic[indices]
        #

        # Add the control points to the appropriate list.
        ## /** Flattening here since otherwise I wound up with dtype=object
        ##     in the numpy arrays. **/
        #true_control_pts.append(nc_nonan.flatten().tolist())
        #ideal_control_pts.append(ic_nonan.flatten().tolist())

        #OLD_true.append(nc)     # /** for debug **/
        #OLD_ideal.append(ic)

        true_control_pts.append(nc_by_t)
        ideal_control_pts.append(ic_by_t)

    except (IndexError,AttributeError):
        pass

pdb.set_trace()
# /** Make vectors of all the control points. **/
true_flat = list(flatten_list(true_control_pts))
ideal_flat = list(flatten_list(ideal_control_pts))
true_control_pts = np.array(true_flat)
ideal_control_pts = np.array(ideal_flat)

# Make the vectors 2d
length = len(true_control_pts)/2
true_control_pts = true_control_pts.reshape(length,2)
ideal_control_pts = ideal_control_pts.reshape(length,2)

#pdb.set_trace()

# Plot the original image
image = np.array([range(i,i+10) for i in range(20)]) / 30.0 # /**.astype(np.int32)**/
    # /** You don't want int32 images!  See
    #     http://scikit-image.org/docs/dev/user_guide/data_types.html .
    #     Manually rescale the image to [0.0,1.0] by dividing by 30. **/
image_float = img_as_float(image) #/** make sure skimage is happy */ 
fig = plt.figure()
plt.imshow(image_float, origin='lower', interpolation='none')
plt.title('Input image')
fig.show()  # /** I needed this on my test system **/

# Warp the actual image given the transformation between the true and ideal wavelength maps
tform = transform.PiecewiseAffineTransform()
tform.estimate(true_control_pts, ideal_control_pts)
out = transform.warp(image, tform)
    # /** since we started with float, and this is float, too, the two are
    #     comparable. **/

pdb.set_trace()

# Plot the warped image!
fig, ax = plt.subplots()
ax.imshow(out, origin='lower', interpolation='none')    # /**note: float**/
plt.title('Should be parallel lines')
fig.show()

# /** Show the control points.
#     The z=0 plane will be the "true" control points (before), and the
#     z=1 plane will be the "ideal" control points (after). **/
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
fig.show()

for rowidx in range(length):
    ax.plot([true_control_pts[rowidx,0], ideal_control_pts[rowidx,0]],
            [true_control_pts[rowidx,1], ideal_control_pts[rowidx,1]],
            [0,1])

input() # /** because I was running from the command line **/

-: X Y , , ? , , . :

nc_by_t = equispace(tc,ic.shape[0])
ic_by_t = equispace(ic,ic.shape[0])

ic.shape[0] , . equispace . , , , , , , ic. , ; 100 .

, . , . X Y? ? ( )?

,

? !

, .

( - . ===========)

import pdb
import matplotlib.pyplot as plt
import numpy as np
from skimage import measure, transform, img_as_float

from mpl_toolkits.mplot3d import Axes3D # /**/

#/**/
def test_figure(data,title):
    image_float = img_as_float(data) #/** make sure skimage is happy */ 
    fig = plt.figure()
    plt.imshow(image_float, origin='lower', interpolation='none')
    plt.title(title)
    fig.show()

#/**/
# From https://stackoverflow.com/a/14491059/2877364 by
# https://stackoverflow.com/users/1355221/dansalmo
def flatten_list(L):
    for item in L:
        try:
            for i in flatten_list(item): yield i
        except TypeError:
            yield item
#end flatten_list

#/**/
# Modified from /questions/490589/how-to-generate-equispaced-interpolating-values/2094239#2094239 by
# https://stackoverflow.com/users/2588210/christian-k
def equispace(data, npts):
    x,y = data.T
    xd =np.diff(x)
    yd = np.diff(y)
    dist = np.sqrt(xd**2+yd**2)
    u = np.cumsum(dist)
    u = np.hstack([[0],u])

    t = np.linspace(0,u.max(),npts)
    xn = np.interp(t, u, x)
    yn = np.interp(t, u, y)
    return np.column_stack((xn, yn))

# ======================  BIGGER
true_input = np.array([range(i,i+20) for i in range(30)])  # /** != True **/
ideal = np.array([range(20)]*30)
# ======================    
test_figure(true_input / 50.0, 'true_input')
test_figure(ideal / 20.0, 'ideal')

#pdb.set_trace()
# Find contours of ideal and true_input images and create list of control points for warp
true_control_pts = []
ideal_control_pts = []
OLD_true=[]     # /**/ for debugging
OLD_ideal=[]

for lam in [x+0.5 for x in ideal[0]]:   # I tried the 0.5 offset just to see,
    try:                                # but it didn't make much difference

        # Get the isowavelength contour in the true_input and ideal images
        tc = measure.find_contours(true_input, lam)[0]
            # /** So this might not have very many numbers in it. **/
        ic = measure.find_contours(ideal, lam)[0]
            # /** CAUTION: this is assuming the contours are going the same
            #       direction.  If not, you'll need to make it so. **/
        #nc = np.zeros(ic.shape) # /** don't need ones() **/

        # /** We just want to find points on _tc_ to match _ic_.  That's
        #       interpolation _within_ a curve. **/
        #pdb.set_trace()
        nc_by_t = equispace(tc,ic.shape[0])
        ic_by_t = equispace(ic,ic.shape[0])


        ### /** Not this **/
        ## Use the y /** X? **/ coordinates of the ideal contour
        #nc[:, 0] = ic[:, 0]
        #
        ## Interpolate true contour onto ideal contour y axis so there are the same number of points
        #
        ## /** Have to sort first - https://docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html#numpy-interp **/
        #tc_sorted = tc[tc[:,0].argsort()]
        #    # /** Thanks to https://stackoverflow.com/a/2828121/2877364 by
        #    # https://stackoverflow.com/users/208339/steve-tjoa **/
        #
        #nc[:, 1] = np.interp(ic[:, 0], tc_sorted[:, 0], tc_sorted[:, 1],
        #    left=np.nan, right=np.nan)
        #    # /** nan: If the interpolation is out of range, we're not getting
        #    #     useful data.  Therefore, flag it with a nan. **/

        # /** Filter out the NaNs **/
        # Thanks to https://stackoverflow.com/a/11453235/2877364 by
        # https://stackoverflow.com/users/449449/eumiro
        #pdb.set_trace()
        #indices = ~np.isnan(nc).any(axis=1)
        #nc_nonan = nc[indices]
        #ic_nonan = ic[indices]
        #

        # Add the control points to the appropriate list.
        ## /** Flattening here since otherwise I wound up with dtype=object
        ##     in the numpy arrays. **/
        #true_control_pts.append(nc_nonan.flatten().tolist())
        #ideal_control_pts.append(ic_nonan.flatten().tolist())

        #OLD_true.append(nc)     # /** for debug **/
        #OLD_ideal.append(ic)

        true_control_pts.append(nc_by_t)
        ideal_control_pts.append(ic_by_t)

    except (IndexError,AttributeError):
        pass

#pdb.set_trace()
# /** Make vectors of all the control points. **/
true_flat = list(flatten_list(true_control_pts))
ideal_flat = list(flatten_list(ideal_control_pts))
true_control_pts = np.array(true_flat)
ideal_control_pts = np.array(ideal_flat)

# Make the vectors 2d
length = len(true_control_pts)/2
true_control_pts = true_control_pts.reshape(length,2)
ideal_control_pts = ideal_control_pts.reshape(length,2)

#pdb.set_trace()

# Plot the original image
image = np.array([range(i,i+10) for i in range(20)]) / 30.0 # /**.astype(np.int32)**/
    # /** You don't want int32 images!  See
    #     http://scikit-image.org/docs/dev/user_guide/data_types.html .
    #     Manually rescale the image to [0.0,1.0] by dividing by 30. **/

# ======================    
# /** Pad from 10x20 to 20x30 just for grins **/
#pdb.set_trace()
image = np.concatenate( (np.zeros((20,5)),image,np.zeros((20,5))), 1)
    # now it 20x20
image = np.concatenate( (np.zeros((5,20)),image,np.zeros((5,20))), 0)
# ======================    

#Plot it
image_float = img_as_float(image) #/** make sure skimage is happy */ 
fig = plt.figure()
plt.imshow(image_float, origin='lower', interpolation='none')
plt.title('Input image')
fig.show()  # /** I needed this on my test system **/

# Warp the actual image given the transformation between the true and ideal wavelength maps
tform = transform.PiecewiseAffineTransform()
tform.estimate(true_control_pts, ideal_control_pts)
out = transform.warp(image, tform)
    # /** since we started with float, and this is float, too, the two are
    #     comparable. **/

pdb.set_trace()

# Plot the warped image!
fig, ax = plt.subplots()
ax.imshow(out, origin='lower', interpolation='none')    # /**note: float**/
plt.title('Should be parallel lines')
fig.show()

# /** Show the control points.
#     The z=0 plane will be the "true" control points (before), and the
#     z=1 plane will be the "ideal" control points (after). **/
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
fig.show()

for rowidx in range(length):
    ax.plot([true_control_pts[rowidx,0], ideal_control_pts[rowidx,0]],
            [true_control_pts[rowidx,1], ideal_control_pts[rowidx,1]],
            [0,1])

input() # /** because I was running from the command line **/