6

I am trying to fit a transformation from one set of coordinates to another.

x' = R + Px + Qy
y' = S - Qx + Py
Where P,Q,R,S are constants, P = scale*cos(rotation). Q=scale*sin(rotation)

There is a well known 'by hand' formula for fitting P,Q,R,S to a set of corresponding points. But I need to have an error estimate on the fit - so I need a least squares solution.

Read 'Numerical Recipes' but I'm having trouble working out how to do this for data sets with x and y in them.

Can anyone point to an example/tutorial/code sample of how to do this ?
Not too bothered about the language.
But - just use built in feature of Matlab/Lapack/numpy/R probably not helpful !

edit: I have a large set of old(x,y) new(x,y) to fit to. The problem is overdetermined (more data points than unknowns) so simple matrix inversion isn't enough - and as I said I really need the error on the fit.

3
  • 3
    Do you have a set of (x_i,y_i,x'_i,y'_i)s or x+/-dx, y+/-dy ... or what? If you have exactly one each of x, y, x', y' you can only do an exact solution, and there is no way to extract an error estimate... Aug 7, 2009 at 17:25
  • None of the minimise type functions LMA/Gauss-Newton give a direct error. I suppose I could calculate the best fit and then just work out the error from each point. I had thought it was much simpler than this (ie. a simple mode of a linerar LSquers) and I was just being dumb Aug 7, 2009 at 18:53
  • Interesting problem! I posted some code that should do the trick, but the fun part was breaking out my rusty old math skills again :)
    – e.James
    Aug 7, 2009 at 23:52

6 Answers 6

4

The following code should do the trick. I used the following formula for the residuals:

residual[i] =   (computed_x[i] - actual_x[i])^2
              + (computed_y[i] - actual_y[i])^2

And then derived the least-squares formulae based on the general procedure described at Wolfram's MathWorld.

I tested out this algorithm in Excel and it performs as expected. I Used a collection of ten random points which were then rotated, translated and scaled by a randomly generated transformation matrix.

With no random noise applied to the output data, this program produces four parameters (P, Q, R, and S) which are identical to the input parameters, and an rSquared value of zero.

As more and more random noise is applied to the output points, the constants start to drift away from the correct values, and the rSquared value increases accordingly.

Here is the code:

// test data
const int N = 1000;
float oldPoints_x[N] = { ... };
float oldPoints_y[N] = { ... };
float newPoints_x[N] = { ... };
float newPoints_y[N] = { ... };

// compute various sums and sums of products
// across the entire set of test data
float Ex =  Sum(oldPoints_x, N);
float Ey =  Sum(oldPoints_y, N);
float Exn = Sum(newPoints_x, N);
float Eyn = Sum(newPoints_y, N);
float Ex2 = SumProduct(oldPoints_x, oldPoints_x, N);
float Ey2 = SumProduct(oldPoints_y, oldPoints_y, N);
float Exxn = SumProduct(oldPoints_x, newPoints_x, N);
float Exyn = SumProduct(oldPoints_x, newPoints_y, N);
float Eyxn = SumProduct(oldPoints_y, newPoints_x, N);
float Eyyn = SumProduct(oldPoints_y, newPoints_y, N);

// compute the transformation constants
// using least-squares regression
float divisor = Ex*Ex + Ey*Ey - N*(Ex2 + Ey2);
float P = (Exn*Ex + Eyn*Ey - N*(Exxn + Eyyn))/divisor;
float Q = (Exn*Ey + Eyn*Ex + N*(Exyn - Eyxn))/divisor;
float R = (Exn - P*Ex - Q*Ey)/N;
float S = (Eyn - P*Ey + Q*Ex)/N;

// compute the rSquared error value
// low values represent a good fit
float rSquared = 0;
float x;
float y;
for (int i = 0; i < N; i++)
{
    x = R + P*oldPoints_x[i] + Q*oldPoints_y[i];
    y = S - Q*oldPoints_x[i] + P*oldPoints_y[i];
    rSquared += (x - newPoints_x[i])^2;
    rSquared += (y - newPoints_y[i])^2;
}
1
  • I would normally use more descriptive variable names, but in this case I think long names like sumOfNewXValues would actually make everything harder to read. Mathematical formulae seem to be a special case.
    – e.James
    Aug 7, 2009 at 23:56
4

To find P, Q, R, and S, then you can use least squares. I think the confusing thing is that that usual description of least squares uses x and y, but they don't match the x and y in your problem. You just need translate your problem carefully into the least squares framework. In your case the independent variables are the untransformed coordinates x and y, the dependent variables are the transformed coordinates x' and y', and the adjustable parameters are P, Q, R, and S. (If this isn't clear enough, let me know and I'll post more detail.)

Once you've found P, Q, R, and S, then scale = sqrt(P^2 + Q^2) and you can then find the rotation from sin rotation = Q / scale and cos rotation = P / scale.

5
  • Yeah, the usual presentation of linear least squares is for fitting the scalar equation y = F(x) = \sum_i c_i f_i(x); this problem is fitting the vector equation r' = F(r) = \sum_i c_i f_i(r).
    – las3rjock
    Aug 7, 2009 at 17:49
  • The problem is that linear LSF is only in one variable + 2 unknowns. I have 2 variables and 4 unknowns (well three if you assume scale is the same) Aug 7, 2009 at 18:07
  • Linear least squares (in the usual presentation) actually works for one variable and n unknowns. See en.wikipedia.org/wiki/Linear_least_squares , where the first graphical example is for fitting a second-order polynomial (n = 3).
    – las3rjock
    Aug 7, 2009 at 18:45
  • Least squares can handle m > 1 independent variables and n > 1 dependent variables. In my edition of Numerical Recipes its covered in section 15.4 General Linear Least Squares or look at the wikipedia article referenced above. Aug 7, 2009 at 19:24
  • Thanks I think that chapter now makes more sense! When I did this before there was always some trick to just turning it into a linear equation. The by-hand solution is just a simple sum of square of differnces process I assumed it was simple and I was missing something. Aug 7, 2009 at 19:58
1

You can use the levmar program to calculate this. Its tested and integrated into multiple products including mine. Its licensed under the GPL, but if this is a non-opensource project, he will change the license for you (for a fee)

4
  • Thanks, there are lots of free LMA implmentations. I hadn't appreciated that it was a necessary approach for what looked like a simple fit. Aug 7, 2009 at 18:50
  • Its the error estimate that is hard. that requires the Jacobian of the solution matrix (IIRC). Keep in mind, that LS errors are not errors in the sense that the world thinks about errors. They're a measure of the numerical stability of the answer. So, something could be the correct solution, but not particularly stable (i.e. changing one value slightly leads to a large change in the objective functions).
    – Steve
    Aug 7, 2009 at 19:22
  • Yes I think the best approach to the error from the user's point of view is to calculate the best position and then the residual on each point and take the standard deviation of these. Aug 7, 2009 at 19:49
  • @mgb - Be aware that sophisticated users will not be expecting that (and may not want it). The LM errors are what are typically reported by mathematical packages
    – Steve
    Aug 7, 2009 at 20:36
1

Define the 3x3 matrix T(P,Q,R,S) such that (x',y',1) = T (x,y,1). Then compute

A = \sum_i |(T (x_i,y_i,1)) - (x'_i,y'_i,1)|^2

and minimize A against (P,Q,R,S).

Coding this yourself is a medium to large sized project unless you can guarntee that the data are well conditioned, especially when you want good error estimates out of the procedure. You're probably best off using an existing minimizer that supports error estimates..

Particle physics type would use minuit either directly from CERNLIB (with the coding most easily done in fortran77), or from ROOT (with the coding in c++, or it should be accessible though the python bindings). But that is a big installation if you don't have one of these tools already.

I'm sure that others can suggest other minimizers.

1

Thanks eJames, thats almost exaclty what I have. I coded it from an old army surveying manual that was based on an earlier "Instructions to Surveyors" note that must be 100years old! (It uses N and E for North and East rather than x/y )

The goodness of fit parameter will be very useful - I can interactively throw out selected points if they make the fit worse.

FindTransformation(vector<Point2D> known,vector<Point2D> unknown) {
{
    // sums
    for (unsigned int ii=0;ii<known.size();ii++) {
       sum_e += unknown[ii].x;
       sum_n += unknown[ii].y;
       sum_E += known[ii].x;
       sum_N += known[ii].y;                            
       ++n;         
    }

    // mean position
    me = sum_e/(double)n;
    mn = sum_n/(double)n;
    mE = sum_E/(double)n;
    mN = sum_N/(double)n;

    // differences
    for (unsigned int ii=0;ii<known.size();ii++) {

       de = unknown[ii].x - me;
       dn = unknown[ii].y - mn;

       // for P
       sum_deE += (de*known[ii].x);
       sum_dnN += (dn*known[ii].y);
       sum_dee += (de*unknown[ii].x);
       sum_dnn += (dn*unknown[ii].y);

       // for Q
       sum_dnE += (dn*known[ii].x);
       sum_deN += (de*known[ii].y);                     
   }

double P = (sum_deE + sum_dnN) / (sum_dee + sum_dnn);
double Q = (sum_dnE - sum_deN) / (sum_dee + sum_dnn);

double R = mE - (P*me) - (Q*mn);
double S = mN + (Q*me) - (P*mn);
}
1
  • Cool. I shudder to think of how they would have calculated it 100 years ago :)
    – e.James
    Aug 8, 2009 at 2:57
0

One issue is that numeric stuff like this is often tricky. Even when the algorithms are straightforward, there's often problems that show up in actual computation.

For that reason, if there is a system you can get easily that has a built-in feature, it might be best to use that.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.