Affine Transformation Parameters

[Bill Thoen]

Does anyone here know how to calculate the affine transformation
coefficients (A,B,C,D,E,F as used in the "CoordSys Earth Affine
Units..." clause) given a set of coordinates in the original system and
a matching set in the transformed system? I don't want a utility to do
this; I want to know how it's done.

TIA,
- Bill Thoen

[Dragos]

Try this:
http://earth-info.nga.mil/GandG/coordsys/datums/helmert.html
Dragos

[Bill Thoen]

That's interesting but not really what I'm looking for. What I'm trying
to learn is how to determine the *factors* shown in those matrices, not
the matrices themselves. Bill Wemple pointed me to a KB article
(http://testdrive.mapinfo.com/techsupp/miprod.nsf/kbase_by_product/0B57EA4620B5EBF285256A010082A16A)
on the MapInfo site that explains how to do this if you have three,
non-collinear points. However, the extra problem for me is that I have a
collection of points --not just three-- so I need some sort of a best
fit technique applied to determining the 6 factors in MapInfo's affine
transform matrix.

[Spencer Simpson]

Bill: 

 

I've sent this in HTML because of the mathematical formulae.  I also typed this in out of memory so there may be some technical errors.

 

I have some Pascal code to do the least squares fit, which involves some basic linear algebra (i.e. matrix manipulation) to do a two-dimensional linear regression.   I may have tried to port this to MapBasic once, but I've never had a MapBasic project that actually required it, so any MapBasic code I have for this is of dubious value.

 

The basic theory:

 

Given N sets of control points, where each point has coordinates in a "source" coordinate system and coordinates in a "destination" coordinate system, you structure the coordinates into two matrices:

 

A is an Nx3 matrix of coordinates in the "source" coordinate system.   Each row represents one control point.  The first column contains the X coordinates of each control point, the second column contains the Y coordinates, and the third column always contains the value 1 (We're turning Cartesian coordinates in 2-D Euclidean space into homogenous coordinates for "ordinary" points, a finite distance from the origin, in real projective space).

 

B is an Nx3 matrix of coordinates in the "destination" coordinate system, similarly structured.  

 

Then the matrix to transform points from the "source" coordinate system to the "destination" coordinate system is going to be

 

(A^T*A)^-1*(A^T*B)

 

You need at least three points, otherwise the matrix (A^T*A) is singular (i.e. doesn't have an inverse), and the result is undefined.  Some interesting stuff goes on if there are more than three points (which I'll go into later).   The result is a 3x3 matrix with the least-squares homogenous transformation. 

 

[ a  b  c ]

[ d  e  f ]

[ u  v  w ]

 

If you're creating an affine transformation for a MapInfo projection, the "source" coordinate system is the raw un-transformed map projection, and the "destination" is the projection with the transformation applied.  So, use a, b, c, d, e and f in your Coordsys clause (I may have this transposed).

 

Don't worry about u, v, and w, except as an indicator of how much error you can expect.

 

Why?  This method doesn't create an affine transformation so much as it creates a perspectivity transformation between two real projective spaces. u, v, and w represent the homogenous coordinates of a line in the "source" space that becomes the "line at infinity" (with coordinates [0, 0, 1]) in the "destination" space.

 

However, the set of affine transformations is a subset of the set of perspectivities, and since all of the points are "ordinary" points (i.e. the third homogenous coordinate is always 1), you theoretically wind up with an affine transformation, where the line at infinity doesn't change, so in an ideal situation with no rounding error and no measurement error in any of the coordinates,  u=0, v=0, and w=1.   This should be the case (or nearly so due to rounding error) if there are exactly three points.   With measurement error, u, v, and w will deviate from their "ideal" values.  In particular, w is the bivariate correlation coefficient.  If w is significantly different from 1, you probably have a bad point, which you'll have to exclude.  But multiplying a, b, c, d, e and f by 1/w may give better results.

 

I do need to warn you that you can run into significant problems with rounding error if you're not careful (i.e. expecting high precision for widely dispersed points).  Rounding error accumulates quickly, especially when coordinate values are high, so you may need to translate to the centroid of your control points in the "source" coordinate system, and translate back from the centroid of the control points in the "destination" system. Of course, the matrices should contain with values with the appropriate centroid coordinates subtracted.  Perform the above calculation, then calculate the final matrix with

 

[ 1  0  -cx_source ] [ a  b  c ] [ 1  0  cx_dest ]

[ 0  1  -cy_source ] [ d  e  f ] [ 0  1  cy_dest ]

[ 0  0    1     ] [ 0  0  1 ] [ 0  0  1 ]

 

As a side note, you could theoretically use points an infinite distance from the origin as control points (e.g. you could use the North Pole in a Mercator projection).  The third column for that point would have a 0 instead of a 1.  In actual practice, however, this wouldn't give you results you could use in MapInfo.

 

HTH

Spencer

[Bill Thoen]

Very interesting Spencer... But I just want to make sure I understand your notation here:

 

(A^T*A)^-1*(A^T*B)

Does this mean to transpose A and multiply it by A, then take the inverse of that result and multiply that by the transpose of A multiplied by B?

Thanks,
- Bill Thoen

[Spencer Simpson]

Yes, that's exactly how to interpret that formula.  In theory. If you were writing actual software, there are simple ways to optimize the collection of (A^T*A) and (A^T*B).

 

Consider the matrix A for a single control point

 

A = [ x y 1 ]

 

(A^T*A) is

 

[ x² xy x ]

[ yx y² y ]

[ x   y  1 ]

 

Because of the way matrix multiplication works, (A^T*A) for all N points is just

 

[ Σx²  Σxy Σx ]

[ Σyx Σy²  Σy  ]

[ Σx   Σy   N  ]

 

So, to accumulate that matrix, you run through all of the points, adding x², xy, x, y², y, and 1 for each point in the appropriate places. (A^T*B) is accumulated in an analogous fashion.   Take the inverse of the accumulated (A^T*A)  (using DETERMINANTS, since this is only a 3x3 matrix), multiply by the accumulated (A^T*B), and Bob's your uncle.

 

S

 

 

---

[Bill Thoen]

Amazing! This works perfectly, except that I have to transpose the result. Much easier than what I had worked out yesterday. This is just what I was looking for. Thanks!

-- 
- Bill Thoen
  GISnet - www.gisnet.com
  303-786-9961