Transformation between two set of points

Question

I have object let's say on model image. I want to compute transformation (displacement, scale, rotation) between object on model image and object on target image. I want to make assumption that object's can be treated as 2D so only 2D transformations should be computed.

First I want to to it in manually assisted way. The user selects base point on model image and then target point on target image. The number of points should be defined by user (but no less than some minimum 2-3 points). When points gives different information, the transformation should be averaged and for example from this the quality of matching can be computed.

So the questions is rather about computing transformation of two sets of points, but as I want to do it on image I've added image processing tag.

Especially welcomed are references and advices with some pieces of code or pseudocode.

With two points it's very easy issue, only rotation, scale and displacement of line should be taken, but how to do it with more points, and with averaging it and computing some quality factors.

Current solution is:

void transformFnc(std::vector<PointF> basePoints, std::vector<PointF> targetPoints,
                  PointF& offset, double rotation, double scale)
{
    std::vector<Line> basePointsLines;
    std::vector<Line> targetPointsLines;
    assert(basePoints.size() == targetPoints.size());
    int pointsNumber = basePoints.size();

    for(int i = 0; i < pointsNumber; i++)
    {
         for(int j = i + 1; j < pointsNumber; j++)
         {
             basePointsLines.push_back(Line(basePoints[i], basePoints[j]));
             targetPointsLines.push_back(Line(targetPoints[i], targetPoints[j]));
         }
    }
    std::vector<double> scalesVector;
    std::vector<double> rotationsVector;
    double baseCenterX = 0, baseCenterY = 0, targetCenterX = 0, targetCenterY = 0;
    for(std::vector<Line>::iterator it = basePointsLines.begin(), i = targetPointsLines.begin();
        it != basePointsLines.end(), i != targetPointsLines.end(); it++, i++)
    {
        scalesVector.push_back((*i).length()/(*it).length());
        baseCenterX += (*it).pointAt(0.5).x(); 
        baseCenterY += (*it).pointAt(0.5).y();
        targetCenterX += (*i).pointAt(0.5).x();
        targetCenterY += (*i).pointAt(0.5).y();
        double rotation;
        rotation = (*i).angleTo((*it));
        rotationsVector.push_back(rotation);
    }
    baseCenterX = baseCenterX / pointsNumber;
    baseCenterY = baseCenterY / pointsNumber;
    targetCenterX = targetCenterX / pointsNumber;
    targetCenterY = targetCenterY / pointsNumber;

    offset = PointF(targetCenterX - baseCenterX, targetCenterY - baseCenterY);
    scale = sum(scalesVector) / scalesVector.size();
    rotation = sum(rotationsVector) / rotationsVector.size();
}

Only optimization I can find in this code is to eliminate from scales and rotations those values which differs too much from the rest.

I'm looking for codes or pseudocodes of solution propositions. It can also be references to some codes.

So far from answers I know that:

RANSAC algorithm can be used
I need to look for some affine transform computing algorithm in the least square sense

Do you also want other transformations like stretching and shearing? — Generic Human, Jul 27 '12 at 13:33
No. The assumption is that it is solid object. No perspective, stretching, shearing. — krzych, Jul 27 '12 at 13:36
possible duplicate of [Transform a triangle to another triangle](http://stackoverflow.com/questions/1114257/transform-a-triangle-to-another-triangle) — finnw, Jul 27 '12 at 15:03
You need to accept one of the answers here. They seem to answer your question sufficiently. — Void Star, Jan 16 '13 at 01:24
I will accept one when the bounty will end. Unfortunatelly any of them deserves bounty. — krzych, Jan 16 '13 at 08:13

score 29 · Accepted Answer · edited Jul 24 '14 at 09:54

29

First generalize the problem in a simple affine transformation with a 3x3 affine transformation matrix: i.e.

[M11 M12 M13]
[M21 M22 M23]
[M31 M32 M33]

Since we already know that the third row will always be [0 0 1] we can simply disregard it.

Now we can describe the problem as the following matrix equation

[xp0]     [x0 y0 1  0  0  0 ]
[yp0]     [0  0  0  x0 y0 1 ]     [M11]
[xp1]     [x1 y1 1  0  0  0 ]     [M12]
[yp1]  =  [0  0  0  x1 y1 1 ]  *  [M13]
[xp2]     [x2 y2 1  0  0  0 ]     [M21]
[yp2]     [0  0  0  x2 y2 1 ]     [M22]
[xp3]     [x3 y3 1  0  0  0 ]     [M23]
[yp3]     [0  0  0  x3 y3 1 ]

where xp and yp are the projected coordinates and x and y are the original coordinates.

Let's call this

proj = M * trans

We can then calculate a least squares fit for the transformation by

trans = pinv(M) * proj

where pinv is the pseudo inverse.

This gives us an affine transformation that best fits the points given in the least squares sense.

Now obviously this will also give shear, coordinate flips as well as non-uniform scaling which you did not want so we need to limit the affine transformation in some way to avoid shear. This turns out to be quite easy, we can use a single vector to describe the rotation (direction of the vector) and scaling (magnitude of the vector,) the other vector will simply be orthogonal to it. This reduces the degrees of freedom by two.

M21 = -M12
M22 = M11

So reduce to

[xp0]     [x0  y0 1 0]
[yp0]     [y0 -x0 0 1]
[xp1]     [x1  y1 1 0]     [M11]
[yp1]  =  [y1 -x1 0 1]  *  [M12]
[xp2]     [x2  y2 1 0]     [M13]
[yp2]     [y2 -x2 0 1]     [M23]
[xp3]     [x3  y3 1 0]
[yp3]     [y3 -x3 0 1]

and calculate M21 and M22 from M12 and M11 after we have solved the above matrix equation.

edited Jul 24 '14 at 09:54

Adult Ape

3
2

answered Jan 11 '13 at 03:27

wich

16,709
6
47
72

In this method all points are treated equally though it is not robust to noises. Also description of method how I get the points to be matched is very important for this question. – krzych Jan 14 '13 at 21:08
4

In your question you only asked about the case where the matching was already done manually. It is robust to noise as it uses a least squares fit on an overdetermined system, the more points you match the more the noise will be averaged out. If the noise distribution is the same for all points this solution suffices, if you already know that the noise ratios for the various points you can do better by scaling `proj` and `M` according to the known noise ratios. – wich Jan 15 '13 at 02:03
Ok. I see. Noise distribution in this method won't be the same unfortunatelly. Can you compare the errors of method proposed by you to my current solution using lines between points? In my method I can eliminate noises by eliminating points which gives different rotation and scale than others. – krzych Jan 15 '13 at 09:31
I think you are confusing instances with distributions, of course point a will have a different error than point b, but the question is, what is the distribution of errors examined over a large set of problem instances, is the variation in errors for point a over all these instances larger than for point b, or do they have about the same variation? You may want to check http://en.wikipedia.org/wiki/Normal_distribution – wich Jan 18 '13 at 02:24
1

Do note that the least squares method already dampens outliers, if you have 20 points suggesting rotation alpha, and two points suggesting beta, the outcome of the algorithm will be between alpha and beta, but will be much closer to alpha than it is to beta. – wich Jan 18 '13 at 02:25
Maybe I've not precised it clear. But selected points number is user defined, the minimum amount is 3. Case where user wants to select more than 10 points will be rare. Though better elimination of outliers than least squares method is needed. – krzych Jan 18 '13 at 11:09

score 6 · Answer 2 · answered Jan 11 '13 at 07:52

6

I would give the Iterative closest point algorithm a try.

Here you find an implementation with scaling. (SICP)

Another useful link

answered Jan 11 '13 at 07:52

bitWorking

12,485
1
32
38

Second link you provided is resistant to scale. I need to compute scale too. Do you know how to find scale from compute transform code? – krzych Jan 14 '13 at 10:09
This is the canonical, well-understood approach to this problem. It's used in feature based panorama stitching software because it's so effective. – Kaganar Jan 15 '13 at 16:48
This approach relies on mass center of features. If features are not distributed evenly error becomes big. – krzych Jan 15 '13 at 17:37

score 2 · Answer 3 · answered Jul 27 '12 at 13:47

For simplicity, suppose your inputs x1,...,xn and outputs y1,...,yn are complex numbers.

You can compute the average displacement by computing avg(y) - avg(x), and once this is done you can subtract the average to both x and y so that they are centered around 0.
You now want to find a rotation and scale. You can represent both as a single complex number z, so that x*z should be as close as possible to y. But x*z is an R-linear function of the (real) coordinates (zx,zy) of z: so you can use classic linear algebra to solve for z such that x*z is as close as possible to y in the least square sense.

Putting it all together, this gives you the optimal transform in the least square sense.

score 1 · Answer 4 · answered Jul 27 '12 at 13:38

Your transformation is an affine transformation, which can be written in a 3*3 matrix. So your problem is basically to compute the least-mean-square-error affine transformation from one set of points to the others.

This problem is quite simply resolved in common computational geometry literature. A good classic book is this one: http://www.robots.ox.ac.uk/~vgg/hzbook/hzbook1.html (no ads, it's just the reference book for most people). You will find all information about 2D and 3D geometry. A quick google on words like "affine transformation LMSE" will also give you information and maybe codes.

Moreover, you can also use other types of algorithms, robust ones, like RANSAC. Depending on your application, it may be interesting to go further in that direction.

As far as I know affine transformation is something different than what I describet. It takes stretching, shearing etc. — krzych, Jul 27 '12 at 17:53
Right it's a generalization. But the way you get equations from positions is the same, you just have one more constraint (no shear). — Bentoy13, Jul 30 '12 at 11:03

score 0 · Answer 5 · answered Apr 29 '13 at 05:58

0

More simple and clear code in Matlab that can give you transform.

And more complex C++ code(using VXL lib) with python and matlab wrapper included.

Or you can use some modificated ICP(iterative closest point) algorithm that is robust to noise.

answered Apr 29 '13 at 05:58

mrgloom

20,061
36
171
301

Transformation between two set of points

5 Answers5

Linked