Given Three Points Compute Affine Transformation

Question

I have two images and found three similar 2D points using a sift. I need to compute the affine transformation between the images. Unfortunately, I missed lecture and the information out there is a little dense for me. What would the general method be for computing this 2x3 matrix?

I have the matrix of points in a 2x3 matrix [x1 y1;x2 y2;x3 y3] but I am lost from there. Thanks for any help.

Answer 1

Usually, an affine transormation of 2D points is experssed as

x' = A*x

Where x is a three-vector [x; y; 1] of original 2D location and x' is the transformed point. The affine matrix A is

A = [a11 a12 a13;
     a21 a22 a23;
       0   0   1]

This form is useful when x and A are known and you wish to recover x'.

However, you can express this relation in a different way. Let

X = [xi yi 1  0  0  0;
      0  0 0 xi yi  1 ]

and a is a column vector

a = [a11; a12; a13; a21; a22; a23]

Then

X*a = [xi'; yi']

Holds for all pairs of corresponding points x_i, x_i'.

This alternative form is very useful when you know the correspondence between pairs of points and you wish to recover the paramters of A.
Stacking all your points in a large matrix X (two rows for each point) you'll have 2*n-by-6 matrix X multiplyied by 6-vector of unknowns a equals a 2*n-by-1 column vector of the stacked corresponding points (denoted by x_prime):

X*a = x_prime

Solving for a:

a = X \ x_prime

Recovers the parameters of a in a least-squares sense.

Good luck and stop skipping class!

Answer 2

Sorry for not using Matlab, but I only worked with Python a little bit. I think this code may help you (sorry for bad codestyle -- I'm mathematician, not programmer)

import numpy as np
# input data
ins = [[1, 1], [2, 3], [3, 2]]  # <- points
out = [[0, 2], [1, 2], [-2, -1]] # <- mapped to
# calculations
l = len(ins)
B = np.vstack([np.transpose(ins), np.ones(l)])
D = 1.0 / np.linalg.det(B)
entry = lambda r,d: np.linalg.det(np.delete(np.vstack([r, B]), (d+1), axis=0))
M = [[(-1)**i * D * entry(R, i) for i in range(l)] for R in np.transpose(out)]
A, t = np.hsplit(np.array(M), [l-1])
t = np.transpose(t)[0]
# output
print("Affine transformation matrix:\n", A)
print("Affine transformation translation vector:\n", t)
# unittests
print("TESTING:")
for p, P in zip(np.array(ins), np.array(out)):
  image_p = np.dot(A, p) + t
  result = "[OK]" if np.allclose(image_p, P) else "[ERROR]"
  print(p, " mapped to: ", image_p, " ; expected: ", P, result)

This code demonstrates how to recover affine transformation as matrix and vector and tests that initial points are mapped to where they should. You can test this code with Google colab, so you don't have to install anything. Probably, you can translate it to Matlab.

Regarding theory behind this code: it is based on equation presented in "Beginner's guide to mapping simplexes affinely", matrix recovery is described in section "Recovery of canonical notation". The same authors published "Workbook on mapping simplexes affinely" that contains many practical examples of this kind.

Answer 3

for implementation purposes in OpenCV you can use cv2.getAffineTransform

getAffineTransform() [1/2]

Mat cv::getAffineTransform  (   const Point2f   src[],
    const Point2f   dst[] 
)   

Python:

cv.getAffineTransform(  src, dst    ) ->    retval

Answer 4

This is for anyone who wants to do it in C. The algorithm is based on this MathStackExchange post. The author shows how to use Gauss-Jordan elimination to find the formula for the Affine transformation coefficients,

/*
*Function: Affine Solver
*Role: Finds Affine Transforming mapping (X,Y) to (X',Y')
*Input: double array[A,B,C,D,E,F], 
*   int array[X-Coordinates], int array[Y-Coordinates],
*   int array[X'-Coordinates],int array[Y'-Coordinates]
*Output:void - Fills double array[A,B,C,D,E,F]
*/

void AffineSolver(double *AtoF, int *X, int *Y, int *XP, int *YP)
{
    AtoF[0] = (double)(XP[1]*Y[0] - XP[2]*Y[0] - XP[0]*Y[1] + XP[2]*Y[1] + XP[0]*Y[2] -XP[1]*Y[2]) / 
          (double)(X[1]*Y[0] - X[2]*Y[0] - X[0]*Y[1] + X[2]*Y[1] + X[0]*Y[2] - X[1]*Y[2]);
          
    AtoF[1] = (double)(XP[1]*X[0] - XP[2]*X[0] - XP[0]*X[1] + XP[2]*X[1] + XP[0]*X[2] -XP[1]*X[2]) / 
          (double)(-X[1]*Y[0] + X[2]*Y[0] + X[0]*Y[1] - X[2]*Y[1] - X[0]*Y[2] +X[1]*Y[2]);
          
    AtoF[2] = (double)(YP[1]*Y[0] - YP[2]*Y[0] - YP[0]*Y[1] + YP[2]*Y[1] + YP[0]*Y[2] -YP[1]*Y[2]) / 
          (double)(X[1]*Y[0] - X[2]*Y[0] - X[0]*Y[1] + X[2]*Y[1] + X[0]*Y[2] - X[1]*Y[2]);
          
    AtoF[3] = (double)(YP[1]*X[0] - YP[2]*X[0] - YP[0]*X[1] + YP[2]*X[1] + YP[0]*X[2] -YP[1]*X[2]) / 
          (double)(-X[1]*Y[0] + X[2]*Y[0] + X[0]*Y[1] - X[2]*Y[1] - X[0]*Y[2] +X[1]*Y[2]);
          
    AtoF[4] = (double)(XP[2]*X[1]*Y[0] - XP[1]*X[2]*Y[0]-XP[2]*X[0]*Y[1] + XP[0]*X[2]*Y[1]+
           XP[1]*X[0]*Y[2] - XP[0]*X[1]*Y[2]) / 
          (double)(X[1]*Y[0] - X[2]*Y[0] - X[0]*Y[1] + X[2]*Y[1] + X[0]*Y[2] - X[1]*Y[2]);
          
    AtoF[5] = (double)(YP[2]*X[1]*Y[0] - YP[1]*X[2]*Y[0]-YP[2]*X[0]*Y[1] + YP[0]*X[2]*Y[1] + YP[1]*X[0]*Y[2] -            YP[0]*X[1]*Y[2]) / 
          (double)(X[1]*Y[0] - X[2]*Y[0] - X[0]*Y[1] + X[2]*Y[1] + X[0]*Y[2] - X[1]*Y[2]);
}

/*
*Function: PrintMatrix
*Role: Prints 2*3 matrix as //a b e
                //c d f
*Input: double array[ABCDEF]
*Output: voids
*/

void PrintMatrix(double *AtoF)
{
    printf("a = %f ",AtoF[0]);
    printf("b = %f ",AtoF[1]);
    printf("e = %f\n",AtoF[4]);
    printf("c = %f ",AtoF[2]);
    printf("d = %f ",AtoF[3]);
    printf("f = %f ",AtoF[5]);
}

int main()
{
/*Test*/
/*Find transform mapping (0,10),(0,0),(10,0) to (0,5)(0,0)(5,0)*/
/*Expected Output*/
//a = 0.500000 b = 0.000000 e = -0.000000
//c = -0.000000 d = 0.500000 f = -0.000000
/*Test*/
    double *AtoF = calloc(6, sizeof(double));
    int  X[] = { 0, 0,10};
    int  Y[] = {10, 0, 0};
    int XP[] = { 0, 0, 5};
    int YP[] = { 5, 0, 0};
    AffineSolver(AtoF,X,Y,XP,YP);
    PrintMatrix(AtoF);
    free(AtoF);
    return 0;
}