Implementing from scratch cv2.warpPerspective()

Question

I was making some experimentations with the OpenCV function cv2.warpPerspective when I decided to code it from scratch to better understand its pipeline. Though I followed (hopefully) every theoretical step, it seems I am still missing something and I am struggling a lot to understand what. Could you please help me?

SRC image (left) and True DST Image (right)

Output of the cv2.warpPerspective overlapped on the True DST

# Invert the homography SRC->DST to DST->SRC
hinv = np.linalg.inv(h)
src = gray1
dst = np.zeros(gray2.shape)
h, w = src.shape

# Remap back and check the domain
for ox in range(h):
    for oy in range(w):

        # Backproject from DST to SRC
        xw, yw, w = hinv.dot(np.array([ox, oy, 1]).T)

        # cv2.INTER_NEAREST
        x, y = int(xw/w), int(yw/w)

        # Check if it falls in the src domain
        c1 = x >= 0 and y < h
        c2 = y >= 0 and y < w

        if c1 and c2:
            dst[x, y] = src[ox, oy]

cv2.imshow(dst + gray2//2)

Output of my code

PS: The output images are the overlapping of Estimated DST and the True DST to better highlight differences.

Oudh

welcome. if you haven't already, please take the tour. – Modeling

What's the idea behind cv2.imshow(dst + src//2)? Do you really want to mix warped and not-warped image? Didnt you want to mix dst with true-warped? – Sumrall

I'd say that's blending, purely for visualization, to compare input and result. I would have recommended dividing both by 2 so the expected average brightness stays the same – Modeling

In his comment he says "overlapping of Estimated DST and the True DST to better highlight differences" but he doesn't blend dst and true-dst but dst and src. So probably his test is just wrong. – Sumrall

a minimal reproducible example would have been nice. there are no values for the homography. – Modeling

@Sumrall my bad. It was a typo. Now I edited the right cv2.imshow(). The output does not follow any blending purpose. I summed half of the ground truth simply to visualize it in the background and check how far my estimate is from the 'truth'. – Oudh

I'm seeing the edits... for ox in range(h) that's gonna get you in trouble when your pictures aren't square... because x goes with w, not h, and y goes with h, not w. and now you've mixed up the indices in the assignment too. you want dst[oy, ox] = src[iy, ix] – Modeling

Answer 1

Your issue amounts to a typo. You mixed up the naming of your coordinates. The homography assumes (x,y,1) order, which would correspond to (j,i,1).

Just use (x, y, 1) in the calculation, and (xw, yw, w) in the result of that (then x,y = xw/w, yw/w). the w factor mirrors the math, when formulated properly.

Avoid indexing into .shape. The indices don't "speak". Just do (height, width) = src.shape[:2] and use those.

I'd recommend to fix the naming scheme, or define it up top in a comment. I'd recommend sticking with x,y instead of i,j,u,v, and then extend those with prefixes/suffixes for the space they're in ("src/dst/in/out"). Perhaps something like ox,oy for iterating, just xw,yw,w for the homography result, which turns into x,y via division, and ix,iy (integerized) for sampling in the input? Then you can use dst[oy, ox] = src[iy, ix]

Modeling

Hey Chris, thank you for your suggestions. I think I am going to improve the namespace. I found your observation very useful and now it works! Thanks! – Oudh

Answer 2

This is a working sample code with minor modifications and typo fixes from the OP following the suggestions of @Christoph for anyone interested.

src_img = post_img
dst = np.zeros(src_img.shape)
height, width = src_img.shape

for ox in range(height):
    for oy in range(width):
        xw, yw, w = homography_matrix.dot(np.array([ox, oy, 1]).T)

        tx, ty = int(xw / w), int(yw / w)

        # Check if it falls in the src domain
        c1 = tx >= 0 and tx < height
        c2 = ty >= 0 and ty < width

        if c1 and c2:
            dst[ty, tx] = src_img[oy, ox]

The homography_matrix was calculated using the OpenCV function cv2.findHomography

Medorra