How can I implement OpenCV's perspectiveTransform in Python

Question

I'm mapping 2d points from a source rectangle to a destination rectangle. I'd like to be able to do this without requiring OpenCV. I've already got getPerspectiveTransform implemented but I'm having trouble finding a source for the math required for perspectiveTransform. Everything I've found is either about using cv2.perspectiveTransform or how to implement cv2.getPerspectiveTransform.

So far I have this and it works:

import numpy as np
import cv2

def getPerspectiveTransform(sourcePoints, destinationPoints):
    """
    Calculates the 3x3 matrix to transform the four source points to the four destination points

    Comment copied from OpenCV:
    /* Calculates coefficients of perspective transformation
    * which maps soruce (xi,yi) to destination (ui,vi), (i=1,2,3,4):
    *
    *      c00*xi + c01*yi + c02
    * ui = ---------------------
    *      c20*xi + c21*yi + c22
    *
    *      c10*xi + c11*yi + c12
    * vi = ---------------------
    *      c20*xi + c21*yi + c22
    *
    * Coefficients are calculated by solving linear system:
    *             a                         x    b
    * / x0 y0  1  0  0  0 -x0*u0 -y0*u0 \ /c00\ /u0\
    * | x1 y1  1  0  0  0 -x1*u1 -y1*u1 | |c01| |u1|
    * | x2 y2  1  0  0  0 -x2*u2 -y2*u2 | |c02| |u2|
    * | x3 y3  1  0  0  0 -x3*u3 -y3*u3 |.|c10|=|u3|,
    * |  0  0  0 x0 y0  1 -x0*v0 -y0*v0 | |c11| |v0|
    * |  0  0  0 x1 y1  1 -x1*v1 -y1*v1 | |c12| |v1|
    * |  0  0  0 x2 y2  1 -x2*v2 -y2*v2 | |c20| |v2|
    * \  0  0  0 x3 y3  1 -x3*v3 -y3*v3 / \c21/ \v3/
    *
    * where:
    *   cij - matrix coefficients, c22 = 1
    */

    """
    if sourcePoints.shape != (4,2) or destinationPoints.shape != (4,2):
        raise ValueError("There must be four source points and four destination points")

    a = np.zeros((8, 8))
    b = np.zeros((8))
    for i in range(4):
        a[i][0] = a[i+4][3] = sourcePoints[i][0]
        a[i][1] = a[i+4][4] = sourcePoints[i][1]
        a[i][2] = a[i+4][5] = 1
        a[i][3] = a[i][4] = a[i][5] = 0
        a[i+4][0] = a[i+4][1] = a[i+4][2] = 0
        a[i][6] = -sourcePoints[i][0]*destinationPoints[i][0]
        a[i][7] = -sourcePoints[i][1]*destinationPoints[i][0]
        a[i+4][6] = -sourcePoints[i][0]*destinationPoints[i][1]
        a[i+4][7] = -sourcePoints[i][1]*destinationPoints[i][1]
        b[i] = destinationPoints[i][0]
        b[i+4] = destinationPoints[i][1]

    x = np.linalg.solve(a, b)
    x.resize((9,), refcheck=False)
    x[8] = 1 # Set c22 to 1 as indicated in comment above
    return x.reshape((3,3))

if __name__ == "__main__":
    # Create a transform to change table coordinates in inches to projector coordinates
    sourceCorners = np.array([[0.0, 0.0],[120.0,0.0],[120.0,63.0],[0.0,63.0]])
    destinationCorners = np.array([[4095.0,0],[3071,4095],[1024,4095],[0,0]])
    perspectiveTransform = getPerspectiveTransform(sourceCorners, destinationCorners)

    points = np.array([0,0,120,63,120,0,0,63,120,63,0,0,120,0], dtype=float).reshape(-1,1,2)
    perspectivePoints = cv2.perspectiveTransform(points, perspectiveTransform)

    print(perspectivePoints)

Result:

[[[4095.    0.]]    
 [[1024. 4095.]]    
 [[3071. 4095.]]    
 [[   0.    0.]]    
 [[1024. 4095.]]    
 [[4095.    0.]]    
 [[3071. 4095.]]]

Graphs:

The OpenCV C source for perspectiveTransform has cryptic variable names, no comments, and is pretty unreadable.

Can anyone point me to a good source for how to implement perspectiveTransform?

Answer 1

Basically the perspective transformation is

[                  ]   [source_x]   [target_x * w]  
[perspectivesMatrix] x [source_y] = [target_y * w]
[                  ]   [    1   ]   [      w     ]

where perspectiveMatrix is a 3x3 matrix of the form

[c00 c01 c02]
[c10 c11 c12]
[c20 c21 c22]

Since you already have perspectiveMatrix, all we need to do is replicate the previous formula.

def perspectiveTransform(perspectiveMatrix, sourcePoints):
    '''
    perspectiveMatrix as above
    sourcePoints has shape (n,2)
    '''

    # first we extend source points by a column of 1
    # augment has shape (n,1)
    augment = np.ones((sourcePoints.shape[0],1))
    # projective_corners is a 3xn matrix with last row all 1
    # note that we transpose the concatenation
    projective_corners = np.concatenate( (sourceCorners, augment), axis=1).T

    # projective_points has shape 3xn
    projective_points = perspectiveMatrix.dot(projective_corners)

    # obtain the target_points by dividing the projective_points 
    # by its last row (where it is non-zero)
    # target_points has shape (3,n).
    target_points = np.true_divide(projective_points, projective_points[-1])

    # so we want return points in row form
    return target_points[:2].T


if __name__=='__main__':
    # Create a transform to change table coordinates in inches to projector coordinates
    sourceCorners = np.array([[0.0, 0.0],[120.0,0.0],[120.0,63.0],[0.0,63.0]],dtype=np.float32)
    destinationCorners = np.array([[4095.0,0],[3071,4095],[1024,4095],[0,0]],dtype=np.float32)
    perspectiveMatrix = getPerspectiveTransform(sourceCorners, destinationCorners)

    # test points
    points = np.array([0,0,120,63,120,0,0,63,120,63,0,0,120,0], dtype=float)

    # perspectiveTransform by cv2
    cv_perspectivePoints = cv2.perspectiveTransform(points.reshape(-1,1,2), perspectiveMatrix)

    # our implementation of perspectiveTransform
    perspectivePoints = perspectiveTransform(perspectiveMatrix, points)

    # should yields something close to 0.0
    print(cv_perspectivePoints.reshape(-1,2) - perspectivePoints)