39 views (last 30 days)

Jan
on 9 Nov 2012

Edited: Jan
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...

Matt J
on 9 Nov 2012

Edited: Matt J
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

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

Luis Isaac
on 13 Dec 2017

Edited: Luis Isaac
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

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.