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
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.
I would give the Iterative closest point algorithm a try.
Here you find an implementation with scaling. (SICP)
Another useful link
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.
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.
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.