How to properly shear a line in a 2D image using interpolation in Python?

Question

I'm trying to apply a shearing transformation to a simple 2D line filter (represented as a binary image) using interpolation methods in Python. However, the resulting image after applying the shear transformation looks almost identical to the input filter, with no visible shearing effect.

When I apply the same shearing transformation to a Gaussian function (another 2D image), the shearing works as expected.

Here’s the code I’m using to apply the shear:

[![import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate, ndimage

# Define the original 2D function f(x, y)
def f(x, y, sigma, muu):
    dst = np.sqrt(x ** 2 + y ** 2)
    normal = 1 / (2.0 * np.pi * sigma ** 2)
    gauss = np.exp(-((dst - muu) ** 2 / (2.0 * sigma ** 2))) * normal
    return gauss

# Dimensions
x_dim, y_dim = 5, 19
x_values = np.linspace(-(x_dim // 2), (x_dim // 2), x_dim)
y_values = np.linspace(-(y_dim // 2), (y_dim // 2), y_dim)
xv, yv = np.meshgrid(x_values, y_values, indexing='ij')

# Shearing parameters
a = 2
b = 1
interp_method = "linear"

# Generate image A and sheared version A_ab
A = f(xv, yv, 2, 0)
A_ab = f(xv * a + yv * b, yv, 2, 0)

# Create a line filter
filter = np.zeros((x_dim, y_dim))
filter\[:, y_dim // 2\] = 1

# Interpolation of A and the filter
interp_A = interpolate.RegularGridInterpolator(
    (x_values, y_values), A, bounds_error=False, fill_value=None, method=interp_method)
interp_filter = interpolate.RegularGridInterpolator(
    (x_values, y_values), filter, bounds_error=False, fill_value=None, method=interp_method)

# Shearing transformation
out_xv = xv * a + yv * b
out_yv = yv
sheared_coords_EPI = np.transpose(np.array((out_xv, out_yv)), axes=(1, 2, 0))
sheared_A = interp_A(sheared_coords_EPI)
sheared_filter = interp_filter(sheared_coords_EPI)

# Shift coordinates to positive space for warping
x_shift = x_dim // 2
y_shift = y_dim // 2
shifted_out_xv = out_xv + x_shift
shifted_out_yv = out_yv + y_shift

# Warp the filter using the new coordinates
sheared_filter2 = ndimage.map_coordinates(filter, \[shifted_out_xv, shifted_out_yv\], order=1, mode='constant', cval=0.0)

# Flatten and interpolate the filter
points = np.vstack((xv.flatten(), yv.flatten())).T
values = filter.flatten()
new_points = np.vstack((out_xv.flatten(), out_yv.flatten())).T
interpolated_filter = interpolate.griddata(points, values, new_points, method='linear')
interpolated_filter = interpolated_filter.reshape((x_dim, y_dim))

# Plotting
plt.figure(), plt.imshow(A), plt.title("A")
plt.figure(), plt.imshow(A_ab), plt.title("A_ab")
plt.figure(), plt.imshow(sheared_A), plt.title("sheared_A")
plt.figure(), plt.imshow(filter), plt.title("filter")
plt.figure(), plt.imshow(sheared_filter), plt.title("sheared_filter")
plt.figure(), plt.imshow(sheared_filter2), plt.title("sheared_filter2")
plt.figure(), plt.imshow(interpolated_filter), plt.title("interpolated_filter")
plt.show()]

Answer 1

Your shear works in a direction other than what you expect.

Try to put a line in the filter array as filter[x_dim // 2, :] = 1 and you will see that everything works fine.

But the vertical line at 0 will stay untouched, and 0 is exactly in the center of your picture.

Gaussian has circular symmetry, that's why shear transformation will modify it the same in horizontal and vertical directions.

Probably, the main problem in your code is that you use x in all places where it should be y and vice versa.

I.e. when you create the filter, you call np.zeros((x_dim, y_dim)), where x_dim=5 and y_dim=19 but if you look at the picture or at the array itself, you will see that its height is 5 and its width is 19 - axes are mixed up.

In numpy (and in the field of linear algebra in general) y is the first coordinate, and x is the second, because matrix is represented as a array of rows (and we choose row by first index, and than position in the row by second index).