Finding center of rotation from two sets of points

Given a set of 2D points and the same set of 2D points rotated about an unknown point, how can I find the center of rotation and the amount of rotation?

Answer 1

[Jan]

Dr. Yong-San Yoon, Korea Advanced Institute of Science and Technology KAIST has suggested the following method:

Solve the following linear system to get the least square solution:

 m points, n frames
 [Xc, Yc, Zc]: center of rotation
 xji, yji, zji: Cartesian coordinates of the j.th point in the i.th frame.
 P<j>: Residuals of radii: P<j> = R<j>^2 - (Xc^2 + Yc^2 + Zc^2)

 / 2*x11 2*y11 2*z11 1 0 ... \          / x11^2 + y11^2 + z11^2 \
 | 2*x12 2*y12 2*z12 1 0 ... | / Xc \   | x12^2 + y12^2 + z12^2 |
 | 2*x13 2*y13 2*z13 1 0 ... | | Yc |   | x13^2 + y13^2 + z13^2 |
 | ...                       | | Zc |   | ...                   |
 | 2*x1n 2*y1n 2*z1n 1 0 ... | | P1 |   | x1n^2 + y1n^2 + z1n^2 |
 | ...                       | | P2 |   | ...                   |
 | 2*x21 2*y21 2*z21 0 1 ... | | .. | = | x21^2 + y21^2 + z21^2 |
 | 2*x22 2*y22 2*z22 0 1 ... | | .. |   | x22^2 + y22^2 + z22^2 |
 | 2*x23 2*y23 2*z23 0 1 ... | | .. |   | x23^2 + y23^2 + z23^2 |
 | ...                       | | .. |   | ...                   |
 | 2*xm1 2*ym1 2*zm1 ... 0 1 | | .. |   | xm1^2 + ym1^2 + zm1^2 |
 | ...                       | \ Pm /   | ...                   |
 \ 2*xmn 2*ymn 2*zmn ... 0 1 /          \ xmn^2 + ymn^2 + zmn^2 /

For the implementation in Matlab I assume the 3D case - then 3 frames are required, because a rotation between 2 frames defines an axes but not one center of rotation. For 2D two frames are enough and the formula can be adjusted accordingly:

p1 = rand(3, 3);  % [3 frames, 3 coordinates]
p2 = rand(3, 3);
p3 = rand(3, 3);
p4 = rand(3, 3);
nPoint = 4;
nFrame = 3;
MPM = cat(1, p1, p2, p3, p4);      % Points as matrix
A   = cat(2, 2 .* MPM, zeros(nFrame * nPoint, nPoint));
for iN = 1:nPoint
  A(((iN - 1) * nFrame + 1):(iN * nFrame), 3 + iN) = 1;
end

b   = sum(MPM .* MPM, 2);
x   = A \ b;
C   = transpose(x(1:3));
Res = A * x - b;

The more points, the better the results. But noise and translation will be a severe problem: If you mark a bunch of points on a rigid body, rotate it around several axes (or around a point in 2D). If there is a translation in addition, the found "center of rotation" depends on the position of the marked points...

Answer 2

[Matt J]

Using this file you can find the rotation matrix and R and translation vector t that maps the first set of points to the second one. It will also give you the angle of rotation in degrees, assuming we're talking about 2D. Once you have R and t, you need to solve for the point x that doesn't change under the transformation,

    x=R*x+t

which would be

    x=(eye(2)-R)\t

[Hans Dohse]

Can this solution be adapted to use the result of one of the cp2tform() transforms

[Matt J]

If you have a known affine transform

y=A*x+b

you can try to find a fixed point for it by solving

x=A*x+b

Answer 3

[Luis Isaac]

There is a very elegant solution of this problem using complex numbers in 2D

function [S,Theta,t,c]=FindSTR3(c1,c0)
% calculate the finite centre of rotation, the angle of rotation,
% and the scale factor
% Input: E, F - lists of 2D corresponding points. These are n by 2 matrices
% where n is the number of points and the x val is in column 1 y val in column 2.
% Output: S -     Scale
%         c -     Centre of rotation coordinates.
%         theta - Angle of rotation between the two sets about the centre.
%         t -     Traslation vector

% first do some error checking
    [rE, cE] = size(c1);
    [rF, cF] = size(c0);
    if (cE ~= 2 ) || (cF ~= 2 )
        error('E, F should be an nx2matrix')
    end
    if (rE ~= rF)
        error('matrices E and F are of different size')
    end

      A=c1(:,1)+c1(:,2)*1j;       % Complexify the data
      B=c0(:,1)+c0(:,2)*1j;   
      meanA=sum(A)/rE;            % Mean of each of the pointsets
      meanB=sum(B)/rF;
      A=A-meanA;                  % Translate both pointsets to the origin
      B=B-meanB;
      x = pinv(B)*A;              % LS solution to the rotation between A and B

      Theta = angle(x);           % Optimal angle of rotation
      S = abs(x);                 % Optimal magnification
      v=meanA-(x/S)*meanB;        % Optimal translation vector as complex number
      fcr=v/(1-x/S);              % Optimal origin of rotation as complex number

      t=[real(v),imag(v)];
      c=[real(fcr),imag(fcr)];
end