Обратное или обратимое отображение ndimage.map_coordinates на основе планарной гомографии

Я провел рефакторинг вашего кода, чтобы посмотреть, что происходит.

Во-первых, map_coordinates() выполняет «вытягивание», то есть втягивает исходные пиксели в результирующую сетку, используя предоставленные вами индексы. Эти индексы должны быть сгенерированы с использованием обычной сетки (для результата) и обратного преобразования (для исходного кадра). Вот почему кажется, что ваш квадрат расширяется, а не сжимается.

Затем... удаление Z имеет значение, и где вы отбрасываете "Z" (входы/выходы для преобразования 4x4), особенно при инвертировании преобразования.

Учитывая вращение вне плоскости, скажем вокруг Y, вы получите что-то вроде этого:

[[ 0.70711  0.       0.70711  0.     ]
 [ 0.       1.       0.       0.     ]
 [-0.70711  0.       0.70711  0.     ]
 [ 0.       0.       0.       1.     ]]

Обратное этому:

[[ 0.70711  0.      -0.70711  0.     ]
 [ 0.       1.       0.       0.     ]
 [ 0.70711  0.       0.70711  0.     ]
 [ 0.       0.       0.       1.     ]]

Если вы поместите Z в оба (и примените это к 2D-данным, которые у вас есть), вы теперь получите пару преобразований, которые оба сжимают изображение:

[[0.70711 0.      0.     ]
 [0.      1.      0.     ]
 [0.      0.      1.     ]]

[[0.70711 0.      0.     ]
 [0.      1.      0.     ]
 [0.      0.      1.     ]]

(В вашем случае это каждый раз вызывает расширение из-за map_coordinates() и его операции «вытягивания»)

Сжатие — это появление вращения точек в плоскости, но это не вращение. Удаление Z не поддерживает inv(M) @ M == I.

Повернутые точки, вывернутые из своей плоскости, имеют ненулевое Z, что важно, когда вы хотите повернуть их дальше (например, повернуть назад). Удаление Z означает, что у вас больше нет этой информации. Вы должны принять их положение в пространстве, а трансформация в 2D вместо этого должна сжиматься или растягиваться, в зависимости от того, из какой плоскости вы предполагаете, что они пришли и куда им нужно идти.

Вы должны сначала перенести Z в M (4x4), а затем инвертировать полученную матрицу 3x3. Теперь у вас есть правильная инверсия, которая расширяет изображение, что приводит к тождественному преобразованию.

[[0.70711 0.      0.     ]
 [0.      1.      0.     ]
 [0.      0.      1.     ]]

[[1.41421 0.      0.     ]
 [0.      1.      0.     ]
 [0.      0.      1.     ]]

Теперь немного кода:

def translate4(tx=0, ty=0, tz=0):
    T = np.eye(4)
    T[0:3, 3] = (tx, ty, tz)
    return T

def rotate4(rx=0, ry=0, rz=0):
    R = np.eye(4)
    R[0:3, 0:3] = make_rotation_matrix((rx, ry, rz))
    return R

def dropZ(T4):
    "assumes that XYZW inputs have Z=0 and that the result's Z will be ignored"
    tmp = T4[[0,1,3], :]
    tmp = tmp[:, [0,1,3]]
    return tmp

# input data. don't mind the use of OpenCV. I wasn't in the mood to come up with random() calls to give the square some texture.
im_source = cv.imread(cv.samples.findFile("lena.jpg"), cv.IMREAD_GRAYSCALE)
height, width = im_source.shape[:2]

# transformation: rotate around center
cx, cy = (width-1)/2, (height-1)/2
Tin = translate4(-cx, -cy)
Tout = translate4(+cx, +cy)
R = rotate4(ry=45, rz=30) # with a little in-plane rotation
M = Tout @ R @ Tin

M3 = dropZ(M)
Mi3 = inv(M3)
#print(M3)
#print(Mi3)

# coordinates grid
xin = np.arange(width)
yin = np.arange(height)
xin, yin = np.meshgrid(xin, yin)
zin = np.zeros_like(xin)
win = np.ones_like(xin)

points4 = np.vstack((xin.flatten(), yin.flatten(), zin.flatten(), win.flatten()))
print(points4)

points3 = np.vstack((xin.flatten(), yin.flatten(), win.flatten()))
print(points3)

# always: transform inverted because map_coords is backwards/pull
# can't invert right at the map_coords() call because we've already warped the grid by then

points_warped = inv(M3) @ points3 # apply M3 to identity grid, for input image
print("warped:")
print(points_warped)

points_identity = M3 @ points_warped # apply inv(M3) to warped grid, giving identity grid
# it's equal to M3 @ inv(M3) @ points3
# which is I (identity) @ points3
print("unwarped: (identity grid)")
print(points_identity)

points_unwarping = M3 @ points3 # apply inv(M3) to identity grid, suitable for unwarping *warped* image

# map_coordinates() wants indices, so Y,X or I,J
coords_warped = points_warped.reshape((3, height, width))[[1,0]]
coords_identity = points_identity.reshape((3, height, width))[[1,0]]
coords_unwarping = points_unwarping.reshape((3, height, width))[[1,0]]

im_warped = map_coordinates(im_source, coords_warped)
im_identity = map_coordinates(im_source, coords_identity)
im_unwarped = map_coordinates(im_warped, coords_unwarping)

neighbors = np.hstack((im_source, im_warped, im_identity, im_unwarped))
#neighbors = np.hstack((im1, im2, im3))
plt.figure(figsize=(20,20))
plt.imshow(neighbors, cmap='gray')
plt.show()

К счастью, это все линейно (не нелинейно) и inv(M) @ M == I == M @ inv(M).