## Finding center of rotation from two sets of points

Asked by Hans Dohse

### Hans Dohse (view profile)

on 9 Nov 2012
Latest activity Edited by Luis Isaac

### Luis Isaac (view profile)

on 27 Dec 2017
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?

on 9 Nov 2012
Edited by Jan

### Jan (view profile)

on 9 Nov 2012

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 by Matt J

on 9 Nov 2012
Edited by Matt J

### Matt J (view profile)

on 9 Nov 2012

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

### Hans Dohse (view profile)

on 13 Nov 2012
Can this solution be adapted to use the result of one of the cp2tform() transforms
Matt J

### Matt J (view profile)

on 13 Nov 2012
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 by Luis Isaac

### Luis Isaac (view profile)

on 13 Dec 2017
Edited by Luis Isaac

### Luis Isaac (view profile)

on 27 Dec 2017

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