Given Three Points Compute Affine Transformation
Solution 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!
Solution 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.
Solution 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
Solution 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;
}
Related videos on Youtube
rotsner
Updated on September 23, 2022Comments
-
rotsner over 1 year
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.
-
rotsner about 10 years@chappjc If only that was the same class xD
-
-
rotsner about 10 yearsThank you! Life saver. I was unfortunately at an interview. If only his slides were as concise as your post :)
-
chappjc about 10 yearsGreat answer, and you took the high road; just one thing to add - a convenient way to go from
x
->X
ifx
isN-by-2
is withX = kron(eye(2),[x ones(size(x,1),1)])
. Also, for theX
anda
sizes you're using, you'll needmldivide
to solve the system (a = X \ x_prime
), unless you transpose both (x_prime' / X'
). -
Shai about 10 years@chappjc thanks for these important comments. I'll leave the
kron
part out - to let the poor student do some of his homework himself, but I'll correct themldivide
(I always get confused between them...) -
rayryeng about 7 yearsThanks for the review Shai. I have pairs of matches between two images and I had to implement RANSAC to find the best affine transform to describe the relationship between them...Couldn't quite remember how to set up the system. +1.