Rotating a 2D image using change of coordinates and scipy interpolation-CodePudding

The

After Inclination by an angle of 30°: OK

After Inclination by an angle of 30° AND Rotation by an angle of 40°: PROBLEM

CodePudding user response：

There are several issues here. In the example, I've reduced the resolution and transformation angles a little, to make the issues easier to spot:

champ = 80 # mas diam = 40 # mas nb_points = 20 rho = np.linspace(0,champ/2,nb_points) inclin = 50 angle = 30

Issue 1: mesh vs. coordinate-aligned grid

xv in your code is a 1D array, and the data locations are at any combination of each element of xv with any other element -- a regular, rectangular grid, aligned with the coordinate axes. During inclination transforming, x_inc and y_inc are also 1D arrays, and data is moved to every combination of any element of x_inc with any element of y_inc. This works correctly because the transformation only affects the y axis of your image. this means that for every value of y, the x coordinate is still constant, and for each value of x, the y coordinate is constant.

Rotation, however, breaks that rule. This requires you to use a full-factorial grid of data.

This means instead of handing xv, xv to the transformation function, you need to give it a full factorial mesh, created e.g. by np.meshgrid. The following shows what your code is actually doing, and how it interprets the data (reduced resolution to 20 to make data points discernible):

fig1, ax1 = plt.subplots() ax1.set_aspect('equal') ax1.scatter(x_inc, y_inc, marker='.', c='b', s=2) x_inc_full, y_inc_full = np.meshgrid(x_inc, y_inc) ax1.scatter(x_inc_full, y_inc_full, marker='.', c='b', s=0.1) ax1.scatter(x_rot, y_rot, marker='.', c='r', s=2) x_rot_full, y_rot_full = np.meshgrid(x_rot, y_rot) ax1.scatter(x_rot_full, y_rot_full, marker='.', c='r', s=0.1)

Result:

You can see the code takes the diagonal input line and stretches/rotates it, then re-interprets it as a mesh.

If we instead perform the transformations on a full mesh, things start looking more sensible:

# again, but doing the transformations on the mesh fig2, ax2 = plt.subplots() ax2.set_aspect('equal') x_orig, y_orig = np.meshgrid(xv, xv) x_incf, y_incf, xtot_inc = coord_inclination(x_orig, y_orig, inclin) ax2.scatter(x_incf, y_incf, marker='.', c='b', s=0.2) x_rotf, y_rotf, xtot_rot = coord_rotation(x_orig, y_orig, angle) ax2.scatter(x_rotf, y_rotf, marker='.', c='r', s=0.2)

Result:

Now, that's better!

Except...

Issue 2: the transformations don't do what you think they do

If you look at the previous figure, you notice that the x-axis has been compressed, not the y axis, yet in your inclined figure the circle is wider than tall.

That's because you are interpolating the figure onto the new grid, but then plotting the interpolated figure, on a regular grid.

This means what you're seeing is actually the inverse of the transformation you probably meant to apply. The inclination transform condenses the x axis, but the resulting image is stretched along the x axis, by the same amount. That's visible in your original pictures, but becomes obvious with the updated parameters I'm using:

The solution would be to apply the transformation inversely: interpret the image as being on the transformed grid (which is what you want it to look like), then interpolate it onto the regular grid for display:

interp_incf = scipy.interpolate.interp2d(x_incf, y_incf, image) image_incf = interp_incf(xv, xv)

except...

Issue 3: scipy doesn't like interpolations onto or from grids that are not axis-parallel

The only way I've found to do interpolation involving arbitrary meshes is by running (ugh...) loops around all the point coordinates that we want to interpolate to:

# inclined interp_orig = scipy.interpolate.interp2d(xv, xv, image) image_incf = np.array([[interp_orig(x, y) for x, y in zip(xline, yline)]\ for xline, yline in zip(x_incf, y_incf)]\ ).reshape(x_incf.shape)

However, the above still suffers from Issue number 2 (the transformation works the wrong way round). However, we can't interpolate from an arbitrarily-deformed mesh. This means that the transformation itself needs to be inverted:

def inv_coord_rotation(qu, qv, ang_tmp): ang = ang_tmp*np.pi/180 inv_qu = qu*cos(-ang) qv*sin(-ang) inv_qv = -qu*sin(-ang) qv*cos(-ang) new_q = np.sqrt((inv_qu)**2 (inv_qv)**2) return inv_qu, inv_qv, new_q def inv_coord_inclination(qu, qv, inclination_tmp): inclination = inclination_tmp*np.pi/180 inv_qu = qu/cos(inclination) inv_qv = qv inv_q = np.sqrt((inv_qu)**2 (inv_qv)**2) return inv_qu, inv_qv, inv_q

...and now we can finally do what we came here to do:

# and now for the pictures, with interpolation applied in correct direction: fig3, axes = plt.subplots(ncols=3) # interpolator for the original image: interp_orig = scipy.interpolate.interp2d(xv, xv, image) # inclined interp_orig = scipy.interpolate.interp2d(xv, xv, image) image_incf = np.array([[interp_orig(x, y) for x, y in zip(xline, yline)]\ for xline, yline in zip(x_incfi, y_incfi)]\ ).reshape(x_incfi.shape) axes[0].imshow(image_incf, extent=[min(xv),max(xv),min(xv),max(xv)]) # should be identical to the original inclined version, image_inc # rotated version of the inclined image interp_incf = scipy.interpolate.interp2d(xv, xv, image_incf) x_rotfi, y_rotfi, xtot_rotfi = inv_coord_rotation(x_orig, y_orig, angle) image_rotf = np.array([[interp_incf(x, y) for x, y in zip(xline, yline)]\ for xline, yline in zip(x_rotfi, y_rotfi)]\ ).reshape(x_rotfi.shape) axes[1].imshow(image_rotf, extent=[min(xv),max(xv),min(xv),max(xv)])

however

... I just tried combining the two transformations, to avoid interpolating an already-interpolated image, and it's not trivial, because it would require you to apply the "inclination" transform in the rotated space, and then things get weird.

It's probably possible to use implement a function which takes the non-inverted transformed pixel coordinates and works out the inverse transformation from there, but really, at this point we're just working against scipy.interpolate, similar to my own odyssey with it. I think the general recommendation if you're dealing with images is to use Pillow, but that would be outside of my own experience.