OpenCV warpPerspective implementation

I've implemented the least squares method to find the homomorphic image to fix the rotation and projection in an image.

Now I'm trying to implement the OpenCV warpPerspective method in order to "fix" my image, my python implementation is like:

def fix_image(img, t):
new_image_map = {}
minx, miny = img.shape[0], img.shape[1]
maxx, maxy = 0, 0
for i in range(img.shape[0]):
for j in range(img.shape[1]):
xy = np.array([i, j, 1], np.float64)
uv = np.matmul(t, xy)
uv = uv / uv[2]
minx = min(minx, uv[0])
maxx = max(maxx, uv[0])
miny = min(miny, uv[1])
maxy = max(maxy, uv[1])

new_image_map[int(uv[0]), int(uv[1])] = (i, j)

minx, miny = int(minx), int(miny)
maxx, maxy = int(maxx), int(maxy)
final_img = np.zeros((maxx - minx + 1, maxy - miny + 1)) \
if len(img.shape) == 2 else np.zeros((maxx - minx + 1, maxy - miny + 1, img.shape[2]))

for k, v in new_image_map.items():
final_img[k[0] - minx, k[1] - miny] = img[v]

return final_img


I know that I still need to interpolate the empty points but the problem is that the shape is not right, I'm checking the results by comparing with the actual OpenCV implementation.

dst = cv2.warpPerspective(storm_img, tr, (1448, 1456))


As you can see, it is far from the expected.

resp = fix_image(storm_img, transformation)


I do not want to just use the OpenCV method because I want to learn how to implement it. What am I getting wrong?

The problem in your implementation is that it returns the tensor of float, while an image must be a tensor of int. Because of that, your rendering library, which I assume is matplotlib, cannot correctly plot an image.

To fix that you need to specify the type of final_img explicitly. That is, you need to add a parameter dtype=np.int when creating final_img, giving us the following:

final_img = np.zeros((maxx - minx + 1, maxy - miny + 1), dtype=np.int) \
if len(img.shape) == 2 else np.zeros((maxx - minx + 1, maxy - miny + 1, img.shape[2]), dtype=np.int)


You can also refer to this python implementation of cv2.warpPerspective for more details.

Here is a vectorized implementation for the same (don't need to loop through each pixel and apply homography transform):

def warpPerspective(dst, H):
h, w = dst.shape[:2]
y, x = np.indices((h, w))  # meshgrid
dst_hom_pts = np.stack((x.ravel(), y.ravel(), np.ones(y.size)))
src_hom_pts = np.dot(np.linalg.inv(H), dst_hom_pts)
src_hom_pts /= src_hom_pts[2] # convert from homogenous coordinates
src_hom_pts = np.round(src_hom_pts).astype(int)
src = np.zeros((h, w, 3))
xx, yy = np.maximum(np.minimum(src_hom_pts[1], h-1), 0), \
np.maximum(np.minimum(src_hom_pts[0], w-1), 0)
for ch in range(3): # for each color channel
src[xx, yy, ch] = dst[...,ch].ravel()
return src

src = np.array([[0,0],[0,219],[225,219],[225,0]])
dst = np.array([[82,17],[7,165],[194,202],[219,33]])

# estimate homography matrix H for the forward transform
H = np.array([[ 5.46738946e-01, -3.51284991e-01,  8.20000000e+01],
[ 6.17460513e-02,  4.67917068e-01,  1.70000000e+01],
[-2.83789692e-04, -1.25989102e-03,  1.00000000e+00]]))

dst_img = dst_img / dst_img.max() # normalize

src_img = warpPerspective(dst_img, H)

plt.figure(figsize=(12,5))
plt.imshow(np.hstack((dst_img, src_img))), plt.axis('off')
plt.show()


Use the storm input image and the function warpPerspective() to obtain the image with perspective correction, as shown above.

The homography matrix can be estimated from (src, dst) pair of 4 points using the function cv2.findHomography() or from scratch (as it is done here and here).

Note that there are artifacts since the pixels missed out need to be interpolated, but the basic implementation should be as per the expectation.