[RESOLVED] Register by rotation: a tough one (thawing out though)

[krtxmrtz]

I want to register 2 polygons P and P' of arbitrary shapes with vertices of coordinates
(xk, yk) and (x'k, y'k), k=1,...,N
such that
x'k = xk + a
y'k = yk + b
for all k, where a and b are unknown, i.e. P' has the same shape as P but it is displaced.

However, because the coordinates x' and y' for polygon P' come from experimental measurements they do not exactly match the theoretically calculated (x,y) but they are very close. Therefore, what I do is calculate the center of mass for both polygons to find average values for a and b as the difference of the CoM coordinates.

So far so good. Now the difficult part is that sometimes the P' polygon is slightly rotated with respect to P. So the question is, what parameter could I check to guess the approximate angle of rotation (in a similar manner to the center of mass calculation to account for the shift)?

[zaza]

Hi,

You could get the vector from the centre-of-mass coords to the furthest point (closest point / any point as long as you can be sure that you can get the corresponding point on the other polygon) on each polygon - use the dot product to get the angle between them and that is your rotation angle.
Using the furthest point would be of use if you can't positively identify one point with another straight away but need to "check" - on the other hand it is not so good if your points lie very close to a circle. If that is the case, you might not want to use the centre of mass, but rather a point on the polygon instead.

HTH

zaza

[krtxmrtz]

Thanks zaza for your answer, your suggestions seem sound enough, though the problem is precisily what you have rightly outlined: there are a large number of points (e.g. 1000 or more) and the shape, though not that close to a circle, is quite convex and this makes it very difficult to match any pair of points from each contour. Moreover, the first point (i.e. k=1) could lie on a vertex quite different from the first point on the other contour. This means I can't rely on the farthest point because due to the rotation the farthest point on one contour may not correspond to the farthest point on the other contour.

[zaza]

Really, it depends on how great the spread is and how it compares to the accuracy of your statistics.
If I draw an equilateral triangle in a circle and then move one of the points outside, then copy the figure and rotate it, I know that I can match the two furthest points regardless of how they are numbered, otherwise the figure has been distorted as well as displaced. If the scatter on your experimental points is of the same order as the variation in the points, then obviously this won't work. But, to be honest, if the scatter in the points is larger than the variation in value of the points, then they are the same to all intents and purposes, unless you can identify them using another variable.

For example, if the points form a near-perfect arc, then trying to identify them by their deviation from the arc may be extremely difficult, if not impossible particularly if the deviations are smaller than the scatter. But identifying them in terms of the start and end point of the arc may be a lot easier if you know that one is "the same" as the other.
If you're checking 1000+ points for distribution matching, then...good luck. You might be ok to do a least-squares fit and set the tolerance on the points smaller than the scatter, but it'll be chuggy.

zaza

[krtxmrtz]

...otherwise the figure has been distorted as well as displaced.

As I mentioned above, "... they do not exactly match the theoretically calculated...", i.e. I have some distortion. The thing is the 2 polygons are easily matched by the naked eye, in other words, if I draw them on a sheet of paper and cut out the shape then registration is quite straightforward.

I'll try to make a sample image to show you what the polygons look like. Stay tuned...

[krtxmrtz]

Here is a sample: the 2 displayed curves (blue = theoretical, red = experimental) are actually polygons with a large number of points.

[zaza]

Hi,

Are the angles always likely to be so small? You could calculate a chi-squared value (goodness of fit - basically least-squares) on the basis of the vertical offset of the blue polygon to the red polygon. Although the points won't match up individually, you could calculate this parameter as you vary the angle and minimise it. That way, you're finding the best approximation of one set of points to the other. DO you know least-squares fitting or chi-squared?

HTH

zaza

[krtxmrtz]

Hi,

Are the angles always likely to be so small?
zaza
Yes, I'd say less than 5 deg.

DO you know least-squares fitting or chi-squared?

Nope, I often hear about it but never learned it in the past.

[zaza]

Hi,

Here are some links to least-squares fitting:

Lots of information:
http://mathworld.wolfram.com/LeastSquaresFitting.html

Less information:
http://www.efunda.com/math/leastsquares/leastsquares.cfm

Effectively, you would just "assume" that one of your polygons is the "real" polygon, and calculate a parameter which tells you how different the other one is from that, using the sum of squares of the vertical offset of each point. Find where it is a minimum, and that is your best fit.

Click the links. Read them. You'll like them. :)


HTH

zaza

[krtxmrtz]

Actually I know a good deal about least-squares fitting, it's chi-square that I never got to learn.
To use least-squares with my polygons, don't I have to know the corresponding point in polygon 2 to a specific poit in polygon 1? Because this is something that I don't know either.

[zaza]

Hi,

In a way, it depends on whether what you showed me is actually what you're trying to achieve or not. If that was just an illustration and you're dealing with it mathematically, then I agree that really you need to be comparing point for point. But in the illustration you gave, what you're really after is minimising the area not common to both shapes. The rotation of the shape graphically will not give you the same polygon because you're approximating with the pixels. Hence the best you're ever likely to do is basically a pixel count of the area between the blue curve and the red curve and try to minimise this - that's why I say that the least-squares of the vertical offset would work (the squares is necessary so that you don't subtract areas when the lines cross).

HTH

zaza

[krtxmrtz]

Well, I've made a simple test in an Excel spreadsheet, calculating the vertical scatter with a least squares formula and it appears to work. I'll try to implement it in my VB project. I'll let you know what comes out...

[krtxmrtz]

On second thoughts I don't think least squares fitting can be applied, for the 2 contours have a different number of points.
Maybe I could resample them adding points to either contour at the position where the other contour has a point (using linear interpolation if necessary) but that is not likely to work either, because the lateral extent of the polygons is usually different.

[twanvl]

Just some idea:
If you write the coordinates in polar form (rk,θk) around the center of mass, then the equation becomes:
rk = r'k'
θk = θ'k' + x + 2 n π (for n in Z)
For some x

If you sum r*θ for all points (to keep things shorter I write Σ[r] instead of Σ[rk]:
Σ[r θ] = Σ[r' (θ' + x + 2 nk' π)]
Since n is a whole number, Σ[n] = n' in Z
Σ[r θ] = Σ[r' θ'] + Σ[r' x] + 2 n' π Σ[r']
Factoring out x
Σ[r θ] = Σ[r' θ'] + x Σ[r] + 2 n' π Σ[r']
x Σ[r'] = Σ[r θ] - Σ[r' θ'] - 2 n' π Σ[r']
x = (Σ[r θ] - Σ[r' θ']) / Σ[r'] - 2 n' π
We don't care about a 2π factor, so that leaves us with:
x = (Σ[r θ] - Σ[r' θ']) / Σ[r']

EDIT: This is not correct, because Σ[r' n] ≠ Σ[r']Σ[n]

--------

In case the greek doesn't show up:

Just some idea:
If you write the coordinates in polar form (rk,thetak) around the center of mass, then the equation becomes:
rk = r'k'
thetak = theta'k' + dtheta + 2 n pi (for n in Z)

If you sum r*theta for all points (to keep things shorter I write SUM[r] instead of SUM[rk]:
SUM[r theta] = SUM[r' (theta' + dtheta + 2 nk' pi)]
Since n is a whole number, SUM[n] = n' in Z
SUM[r theta] = SUM[r' theta'] + SUM[r' dtheta] + 2 n' pi SUM[r']
Factoring out dtheta
SUM[r theta] = SUM[r' theta'] + dtheta SUM[r] + 2 n' pi SUM[r']
dtheta SUM[r'] = SUM[r theta] - SUM[r' theta'] - 2 n' pi SUM[r']
dtheta = (SUM[r theta] - SUM[r' theta']) / SUM[r'] - 2 n' pi
We don't care about the 2pi factor, so that leaves us with:
dtheta = (SUM[r theta] - SUM[r' theta']) / SUM[r']

[krtxmrtz]

Just some idea:
If you write the coordinates in polar form ...
Hi! Writing the coordinates in polar form seems a good idea, offhand I think it could provide the means to the resampling I mentioned in my previous post.
On the other hand I can't follow what you're trying to do with your math, I don't understand all the symbols.

[twanvl]

What I was attempting to do was finding the "center of angles", i.e. the center of mass, but each point weighted by its angle. You can already calculate the center of mass, where you weigh each point by its coordinate. Now you do the same, but in polar coordinates instead of cartasian.

The only problem is that method with the area of the trapezium isn't going to work in this case, because finding the center of a trapezium is quite hard. I suggest something using triangles instead (note that triangles would also work for the center of mass). The simple equation is that for the area of a triangle abc = (0,0) (x1,y1) (x2,y2):
A = (x1 * y2 - x2 * y1) / 2
(this is because it is equal to half the cross product of ab and bc). Note that this area is negative if the triangle is counter clockwise, that is if it would fall outside the polygon, this doesn't matter, it will all add up, just like the trapeziums do.

The difficult question is finding the 'center of angle', which is:
∫∫abc theta dx dy
Where ∫∫ is a surface integral, basicly the sum of some function, theta in this case, for all points (x,y) on the triangle abc.
Solving this integral requires some advanced mathematics. I'll give it a shot though.

[twanvl]

The 'center of angle' C is:
C = ∫∫(a,b,c) theta dx dy
First we need a theta, so we rewrite the equation using a Jacobian:
C = ∫theta1theta2∫0R(theta) r theta dr dtheta
Where (r1,theta1) and (r2,theta2) are the two corners of the triangle in polar coordinates. The function R gives us the distance between (0,0) and a point on the line between the two points.

The line between the points is given by:
L(t) = p1 * (1-t) + p2 * t
This applies to both the length, and the angle
Lr(t) = r1 * (1-t) + r2 * t
Ltheta(t) = theta1 * (1-t) + theta2 * t
To find out which t corresponds to theta we solve:
Ltheta(t) = theta
theta1 * (1-t) + theta2 * t = theta
theta1 + (theta2 - theta1) * t = theta
(theta2 - theta1) * t = theta - theta1
t = (theta - theta1) / (theta2 - theta1)
We can use this to contruct our function R(theta)
R(theta) = Lr((theta - theta1) / (theta2 - theta1))
R(theta) = r1 + (r2 - r1) * (theta - theta1) / (theta2 - theta1)

Now we continue with integration, first the inner integral:
∫0R(theta) r theta dr
= theta * [r2/2]0R(t)
= theta * R(t)2
= theta * ( r1 + (r2 - r1) * (theta - theta1) / (theta2 - theta1) )

Now the total integral becomes:
C = ∫theta1theta2 theta * ( r1 + (r2 - r1) * (theta - theta1) / (theta2 - theta1) ) dtheta
This is too much for me, fortunatly the only integrator by wolfram knows the answer (http://integrals.wolfram.com/webMathematica/Integrate.jsp?expr=x%20*%20%28%20r1%20%2B%20%28r2%20-%20r1%29%20*%20%28x%20-%20theta1%29%20/%20%28theta2%20-%20theta1%29%20%29&format=StandardForm&fontsize=Medium) (substitute x for theta), which gives:
C = [ theta^2 (3 r2 theta1 - 3 r1 theta2 + 2 r1 theta - 2 r2 theta) / ( 6 (theta1 - theta2) ) ]theta1theta2
C = [ theta2^2 (3 r2 theta1 - 3 r1 theta2 + 2 r1 theta2 - 2 r2 theta2) / ( 6 (theta1 - theta2) ) ] - [ theta1^2 (3 r2 theta1 - 3 r1 theta2 + 2 r1 theta1 - 2 r1 theta1) / ( 6 (theta1 - theta2) ) ]
C = [ theta2^2 (3 r2 theta1 - 3 r1 theta2 + 2 r1 theta2 - 2 r2 theta2) + theta1^2 (r2 theta1 - 3 r1 theta2 + 2 r1 theta1) ] / ( 6 (theta1 - theta2) )

Good luck implementing that :p, and I probably made a mistake too.

[twanvl]

Since I also have access to mathematica I am also going to use that for calculating the integrals, since it makes less mistakes then I do.

First, defining L(t) to find R(theta)

In[1]:= Lr[t] = r1 * (1-t) + r2 * t

Out[1]= r1 (1 - t) + r2 t

In[2]:= Ltheta[t] = theta1 * (1-t) + theta2 * t

Out[2]= (1 - t) theta1 + t theta2

In[4]:= Solve[Ltheta[t]==theta,t]

-theta + theta1
Out[4]= {{t -> ---------------}}
theta1 - theta2

In[5]:= Lr[t]/.%

-theta + theta1 r2 (-theta + theta1)
Out[5]= {r1 (1 - ---------------) + --------------------}
theta1 - theta2 theta1 - theta2

So far so good, this is the same as my R(theta), had I not rewritten Lr.
Now on to the integral.

In[12]:= R[theta] = r1 + (r2 - r1) * (theta - theta1) / (theta2 - theta1)

(-r1 + r2) (theta - theta1)
Out[12]= r1 + ---------------------------
-theta1 + theta2

In[4]:= Integrate[r theta,{r,0,R[theta]}]

(r1 - r2) (theta - theta1) 2
theta (r1 + --------------------------)
theta1 - theta2
Out[4]= ----------------------------------------
2


In[3]:= Integrate[%,{theta,theta1,theta2}]

Out[3]= -((theta1 - theta2) (2 r1 r2 (theta1 + theta2) +

2 2
> r1 (3 theta1 + theta2) + r2 (theta1 + 3 theta2))) / 24


So there is your answer again.

[twanvl]

Now that you have the 'center of angle' you should be able to rotate the polygon by the difference to align them. You may encounter problems though if the polygon is symetrical, in that case the 'center of angle' will equal 0. To resolve that just choose another point to use as (0,0), it really shouldn't matter. But if you don't use the center of mass, the polygon will move again after rotation.

By the way, see the attachment for an overview of the idea.

[krtxmrtz]

Now I've just realized you're trying a similar approach to what I tried already in the beginning. What I did was to figure out -for each one of the 2 polygons- a sum of all the angles (divided by the number of sides). The difference between these 2 sums was what I would try to minimize by applying slight rotations to one of the contours.
However I soon became aware of a practical difficulty: there is a discontinuity at the origin, as the angle goes from 2*Pi back to 0 (or, if you're using negative angles, at the opposed side where it goes from Pi to -Pi). Therefore, you could have for one vertex of the first contour an angle of 359 deg. but, for the corresponding vertex of the other contour an angle of 0.5 deg. where you would expect to have a very small difference.

[krtxmrtz]

...there is a discontinuity at the origin...
I think I know how to circumvent this. The recipe would be as follows:

1. Set up a new coordinate system at he center of mass (and transform the coordinates accordingly) for each contour.
2. Shift one of the contours so that the CMs coincide.
3. Add up the angles (e.g. in degrees) between the x (horizontal) axis and the line from the center of mass to each vertex. Call this sum S.
4. Calculate S modulus 360 (i.e. one full turn)
5. Divide the result by the number of vertices and call the result Alfa1.
6. Repeat from 3 through 5 for the second contour and call the result Alfa2.
7. Alfa2 - Alfa1 will be the predicted angle that contour #2 is rotated with respect to contour #1.

I have manually tested this for a number of very simple polygons and it seems to work quite well.

[krtxmrtz]

I think I know how to circumvent this. The recipe would be as follows...
...I have manually tested this for a number of very simple polygons and it seems to work quite well.
In actual practice there is yet another difficulty. :(

The number of points in each contour is different and the contour shapes themselves may be (actually are) slightly different, as though the second contour were contructed from the first by slight deformations of the vertices. Also, there may be groups of 3 consecutive points forming a straight line, effectively further increasing the number of points.

Because of this, I think the way to go is to minimize the quantity Alfa2 - Alfa1 described in my previous post, i.e. trying various small angle rotations according to a search algorithm.

[krtxmrtz]

Well, I had almost given up, but after all these months I recently had a new idea in my sleep.

An approach that seems to work in all cases I've tested it is by means of the (planar) moment of inertia, MI.

If an object is rotated about an axis its MI will change unless the contour is symmetric. For example, if a polygon on the x-y plane is rotated about the z axis the MIs about the x axis and about the y axis will vary.
So what I've done is, calculate the MI about the y axis for the first (reference) polygon and for the second (test) contour and compute the difference. Thereafter I use an algorithm to determine an angle such that if the test contour is rotated by it, this difference is a minimum.

Can I rate myself up? :thumb:

[krtxmrtz]

...An approach that seems to work in all cases I've tested it is by means of the (planar) moment of inertia, MI.
...
I had only tested small angles. For large angles it doesn't work because the moment of inertia as a function of the rotation angle has a number of peaks (typically about 4) in the range 0 to 360 deg. i.e. one full turn. Then my minimization algorithm gets trapped in some local function minimum, not necessarily the "good" minimum.
Although I always stick to small (not large than 5 deg.) angles, I'd like to know if some other trick could work for larger angles.