How to transform a distorted square, captured via xy camera coordinates, to a perfect square in cm coordinates?

Viewed 29

I am very sorry for asking a question that is probably very easy if you know how to solve it, and where many versions of the same question has been asked before. However, I am creating a new post since I have not found an answer to this specific question.

Basically, I have a 200cm x 200cm square that I am recording with a camera above it. However, the camera distorts the square slightly, see example here.. I am wondering how I go from transforming the x,y coordinates in the camera to real-life x,y coordinates (e.g., between 0-200 cm for each side). I understand that I probably need to apply some kind of transformation matrix, but I do not know which one, nor how to determine the transformation matrix. I haven't done any serious linear-algebra in a long time, so I appreciate any pointers for what to read up on, or how to get it done. I am working in python, so if there is some ready code for doing the transformation that would also be useful to know.

Thanks a lot!

1 Answers

I will show this using and .

import numpy as np

First, you have to understand the projection model

def apply_homography(H, p1):
    p = H @ p1.T
    return (p[:2] / p[2]).T

With some algebraic manipulation you can determine the points at the plane z=1 that produced the given points.

def revert_homography(H, p2):
    Hb = np.linalg.inv(H)
    # 1 figure out which z coordinate should be added to p2 
    # order to get z=1 for p1
    z = 1/(Hb[2,2] + (Hb[2,0] * p2[:,0] + Hb[2,1]*p2[:,1]))
    p2 = np.hstack([p2[:,:2] * z[:,None], z[:, None]])
    return p2 @ Hb.T

The projection is not invertible, but under the complanarity assumption it may be inverted successfully.

Now, let's see how to determine the H matrix from the given points (assuming they are coplanar).

If you have the four corners in order in order you can simply specify the (x,y) coordinates of the cornder, then you can use the projection equations to determine the homography matrix like here, or here.

This requires at least 5 points to be determined as there is 9 coefficients, but we can fix one element of the matrix and make it an inhomogeneous equation.

def find_homography(p1, p2):
    A = np.zeros((8, 2*len(p1)))
    # x2'*(H[2,0]*x1+H[2,1]*x2)
    A[6,0::2] = p1[:,0] * p2[:,0]
    A[7,0::2] = p1[:,1] * p2[:,0]
    # - (H[0,0]*x1+H[0,1]*y1+H[0,2])
    A[0,0::2] = -p1[:,0] 
    A[1,0::2] = -p1[:,1]
    A[2,0::2] = -1
    
    # y2'*(H[2,0]*x1+H[2,1]*x2)
    A[6,1::2] = p1[:,0] * p2[:,1]
    A[7,1::2] = p1[:,1] * p2[:,1]
    # - (H[1,0]*x1+H[1,1]*y1+H[1,2])
    A[3,1::2] = -p1[:,0]
    A[4,1::2] = -p1[:,1]
    A[5,1::2] = -1
    
    # assuming H[2,2] = 1 we can pass its coefficient
    # to the independent term making an inhomogeneous
    # equation
    b = np.zeros(2*len(p2))
    b[0::2] = -p2[:,0]
    b[1::2] = -p2[:,1]
    h = np.ones(9)
    h[:8] = np.linalg.lstsq(A.T, b, rcond=None)[0]
    return h.reshape(3,3)

Here a complete usage example. I pick a random H and transform four random points, this is what you have, I show how to find the transformation matrix H_. Next I create a test set of points, and I show how to find the world coordinates from the image coordinates.

# Pick a random Homography
H = np.random.rand(3,3)
H[2,2] = 1

# Pick a set of random points
p1 = np.random.randn(4, 3);
p1[:,2] = 1;

# The coordinates of the points in the image
p2 = apply_homography(H, p1)


# testing
# Create a set of random points
p_test = np.random.randn(20, 3)
p_test[:,2] = 1;
p_test2 = apply_homography(H, p_test)

# Now using only the corners find the homography
# Find a homography transform
H_ = find_homography(p1, p2)

assert np.allclose(H, H_)

# Predict the plane points for the test points
p_test_predicted = revert_homography(H_, p_test2)

assert np.allclose(p_test_predicted, p_test)
Related