This, in fact, is the problem of restoration in homography. What you are doing is given by the coordinates in one image and the corresponding coordinates in another image, you are trying to restore the combined translation and rotation matrix, which was used to deform points from one image to another.
You can essentially combine rotation and translation into a single matrix by multiplying two matrices together. Multiplication simply combines two operations. You would get:
H = [cos(theta) -sin(theta) tx] [sin(theta) cos(theta) ty] [ 0 0 1]
The idea is to find the parameters by minimizing the error through the least squares between each pair of points.
Basically, you want to find the following relationship:
xi_after = H*xi_before
H is the combined rotation and translation matrix needed to compare coordinates from one image to another. H also a 3 × 3 matrix and, knowing that the bottom right entry (row 3, column 3) is 1, this simplifies the task. In addition, assuming that your points are in an extended coordinate system, we really want to find this ratio for each coordinate pair from the first image (x_i, y_i) to the other (x_i', y_i') :
[p_i*x_i'] [h11 h12 h13] [x_i] [p_i*y_i'] = [h21 h22 h23] * [y_i] [ p_i ] [h31 h32 1 ] [ 1 ]
The p_i scale is to take into account the scaling and disappearance of points in terms of geography. Let matrix vector multiplication of this equation be performed. We can ignore the 3rd element, since this is not useful to us (for now):
p_i*x_i' = h11*x_i + h12*y_i + h13 p_i*y_i' = h21*x_i + h22*y_i + h23
Now let's look at the third element. We know that p_i = h31*x_i + h32*y_i + 1 . Thus, substituting p_i in each of the equations and rearranging for the solutions for x_i' and y_i' , we obtain:
x_i' = h11*x_i + h12*y_i + h13 - h31*x_i*x_i' - h32*y_i*x_i' y_i' = h21*x_i + h22*y_i + h23 - h31*x_i*y_i' - h32*y_i*y_i'
You now have two equations for each unique pair of points. Now we can build an overdetermined system of equations. Take each pair and build two equations from it. Then you put it in matrix form, i.e.:
Ah = b
A would be a matrix of coefficients constructed from each system of equations using the coordinates from the first image, b would be each pair of points for the second image, and H would be the parameters that you decide. Ultimately, you finally solve this linear system of equations reformulated in matrix form:
You would decide for the vector H , which can be done using least squares. In MATLAB, you can do this via:
h = A \ b;
Sits for you: if moving between images is really just a rotation and translation, then h31 and h32 will be zero after we enable the parameters. However, I always like to be thorough, and so I will still allow for h31 and h32.
NB: This method will work only if you have at least 4 unique pairs of points. Since there are 8 parameters to solve, and there are 2 equations for each point, A must have at least rank 8 in order for the system to be consistent (if you want to use the terminology of linear algebra in a loop). You cannot solve this problem if you have less than 4 points.
If you need MATLAB code, suppose your points are stored in sourcePoints and targetPoints . sourcePoints from the first image, and targetPoints from the second image. Obviously, there should be the same number of dots between the two images. It is assumed that both sourcePoints and targetPoints are stored as M x 2 matrices. The first columns contain your x coordinates, while the second columns contain your y coordinates.
numPoints = size(sourcePoints, 1); %// Cast data to double to be sure sourcePoints = double(sourcePoints); targetPoints = double(targetPoints); %//Extract relevant data xSource = sourcePoints(:,1); ySource = sourcePoints(:,2); xTarget = targetPoints(:,1); yTarget = targetPoints(:,2); %//Create helper vectors vec0 = zeros(numPoints, 1); vec1 = ones(numPoints, 1); xSourcexTarget = -xSource.*xTarget; ySourcexTarget = -ySource.*xTarget; xSourceyTarget = -xSource.*yTarget; ySourceyTarget = -ySource.*yTarget; %//Build matrix A = [xSource ySource vec1 vec0 vec0 vec0 xSourcexTarget ySourcexTarget; ... vec0 vec0 vec0 xSource ySource vec1 xSourceyTarget ySourceyTarget]; %//Build RHS vector b = [xTarget; yTarget]; %//Solve homography by least squares h = A \ b; %// Reshape to a 3 x 3 matrix (optional) %// Must transpose as reshape is performed %// in column major format h(9) = 1; %// Add in that h33 is 1 before we reshape hmatrix = reshape(h, 3, 3)';
Once you're done, you have a combined rotation and translation matrix. If you need x and y translations, just select column 3, rows 1 and 2 in hmatrix . However, we can also work with the vector H , so h13 will be element 3, and h23 will be element 6. If you want the rotation angle, just take the corresponding inverse trigonometric function in rows 1, 2 and columns 1, 2. For vector H this there will be elements 1, 2, 4, and 5. Depending on which elements you select, there will be minimal inconsistency, as this was allowed by least squares. One way to get a good overall angle would probably be to find the angles of all four elements, and then make some kind of middle one. In any case, this is a good starting point.
References
I learned about homography some time ago at the Leow Wee Kheng Computer Vision course. What I told you is based on his slides: http://www.comp.nus.edu.sg/~cs4243/lecture/camera.pdf . Take a look at slides 30-32 if you want to know where I got this stuff from. However, the MATLAB code I wrote myself :)