formula for best approximation for center of 2D rotation with small angles

This is not a homework. I am asking to see if problem is classical (trivial) or non-trivial. It looks simple on a surface, and I hope it is truly a simple problem.

1. Have N points (N >= 2) with coordinates Xn, Yn on a surface of 2D solid body.
2. Solid body has some small rotation (below Pi/180) combined with small shifts (below 1% of distanc开发者_开发技巧e between any 2 points of N). Possibly some small deformation too (<<0.001%)
3. Same N points have new coordinates named XXn, YYn
4. Calculate with best approximation the location of center of rotation as point C with coordinates XXX, YYY.

Thank you

If you know correspondence (i.e. you know which points are the same before and after the transformation), and you choose to allow scaling, then the problem is a set of linear equations. If you have 2 or more points then you can find a least-squares solution with little difficulty.

For initial points (xi,yi) and transformed points (xi',yi') you have equations of the form

xi' = a xi + b yi + c
yi' =-b xi + a yi + d

which you can rearrange into a linear system

A x = y 

where

A = | x1  y1 1 0 | 
    | y1 -x1 0 1 |
    | x2  y2 1 0 |
    | y2 -x2 0 1 |
    |    ...     |
    | xn  yn 1 0 |
    | yn -xn 0 1 |

x = | a |
    | b |
    | c |
    | d |

y = | x1' |
    | y1' |
    | x2' |
    | y2' |
    | ... |
    | xn' |
    | yn' |

the standard "least-squares" form of which is

A^T A x = A^T y

and has the solution

x = (A^T A)^-1 A^T y

with A^T as the transpose of A and A^-1 as the inverse of A. Normally you would use an SVD or QR decomposition to compute the solution as they ought to be more stable and less computationally intensive than the inverse.

Once you've found x (and so the four elements of the transformation a, b, c and d) then the various elements of the transformation are given by

scale       = sqrt(a*a+b*b)
rotation    = atan2(b,a)
translation = (c,d)/scale

If you don't include scaling then the system is non-linear, and requires an iterative solution (but isn't too difficult to solve). If you do not know correspondence then the problem is substantially harder, for small transformations something like iterated closest point works, for large transformations it's a lot harder.

Edit: I forgot to include the centre of rotation. A rotation theta about an arbitrary point p is a sequence

translate(p) rotate(theta) translate(-p)

if you expand it all out as an affine transformation (essentially what we have above) then the translation terms come to

dx = px - cos(theta)*px + sin(theta)*py
dy = py - sin(theta)*px - cos(theta)*py

we know theta (rotation), dx (c) and dy (d) from the equations above. With a little bit of fiddling we can solve for px and py

px = 0.5*(dx - sin(theta)*dy/(1-cos(theta)))
py = 0.5*(dy + sin(theta)*dx/(1-cos(theta)))

You'll notice that the equations are undefined if theta is zero, because there is no centre of rotation when no rotation is performed.

I think I have all that correct, but I don't have time to double check it all right now.

Look up the "Kabsch Algorithm". It is a general-purpose algorithm for creating rotation matrices using N known pairs. Kabsch created it to assist denoising stereo photographs. You rotate a feature in picture A to picture B, and if it is not in the target position, the feature is noise. http://en.wikipedia.org/wiki/Kabsch_algorithm

See On calculating the finite centre of rotation for rigid planar motion for a relatively simple solution. I say "relatively simple" because it still uses things like psuedo-inverses and SVD (singular value decomposition). And here's a wikipedia article on Instant centre of rotation. And another paper: ESTIMATION OF THE FINITE CENTER OF ROTATION IN PLANAR MOVEMENTS.

If you can handle stiffer stuff, try Least Squares Estimation of Transformation Parameters Between Two Point Patterns.

First of all, the problem is non-trivial.

A "simple" solition. It works best when the polygon resembles circle, and points are distributed evenly.

• iterate through N
• For both old and new dataset, find the 2 farthest points of the point N. So now you have the triangle before and after the transformation. Use the clockwise direction from the center of each triangle to number its vertices as [0] (=the N-th point in the original dataset), [1], and [2] (the 2 farthest points).
• Calculate center of rotation, and deformation (both x and y) of this triangle. If the deformation is more then your 0.001% - drop the data for this triangle, otherwise save it.
• Calculate the average for the centers of rotation.

The right solution: define the function Err(Point BEFORE[N], Point AFTER[N], double TFORM[3][3]), where BEFORE - constant old data points, AFTER - constant new data points, TFORM[3][3] affine transformation matrix, Err(...) function that returns the scalar error value, 0.0 when the TFORM translated BEFORE to exact AFTER, or some >0.0 error value. Then use any numeric math you want to find the minimum of the Err(TFORM): e.g. gradient search.

Calculate polygon centers O1 and O2. Determine line formulae for O1 with (X0, Y0) and O2 with (XX0, YY0). Find intersection of lines to get C.

If I understand your problem correctly, this could be solved in this way:

• find extremities (furthest points, probably on several axises)
• scale either one to match
• their rotation should now be trivial (?)

Choose any 2 points on the body, P1, P2, before and after rotation. Find vectors between these before and after points. Cross these vectors with a vector normal to the plane of rotation. This results in two new vectors, the intersection of the lines formed by the initial points and these two new vectors is the center of the rotation.

{
if P1after = P1before return P1after
if P2after = P2before return P2after
Vector V1 = P1after - P1before
Vector V2 = P2after - P2before
normal = Vn // can be messy to create for arbitrary 3d orientation but is simple if you know orientation, for instance, normal = (0,0,1) for an object in the x,y plane)

Vector VL1 = V1 x Vn //Vector V1 cross product with Vn
Vector VL2 = V2 x Vn
return intersectLines(P1after,VL1,P2after,VL2) //Center of rotation is intersection of two lines
}

intersectLines(Point P1, Vector V1, Point P2, Vector V2)
{
//intersect two lines using point, direction form of a line
//returns a Point
}