#929: g.transform: dump coeffs and transform sparse points
--------------------------+-------------------------------------------------
Reporter: hamish | Owner: grass-dev@lists.osgeo.org
Type: enhancement | Status: closed
Priority: normal | Milestone: 7.0.0
Component: Imagery | Version: svn-trunk
Resolution: fixed | Keywords: g.transform
Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by hamish):
I have been pointed to this post by the author:
http://www.cruisersforum.com/forums/f134/charts-31346-36.html#post398493
with GPL3 code linked at the end of the post. (remove the .doc; stupid
webCMS...)
quote:
----
- First for all ref points it converts ref lat/lon in ref E/N using the
exact projection formulas (Mercator and Transverse-Mercator supported)
- Second it tries to approximate the ref X/Y to E/N by polynomial approx.
It tests 23 polynomials to find the best one (minimum error, well
conditioning). It works with at least 2 ref points.
- Third it uses this couple of transformations (lat/lon->E/N->X/Y) to
recalculate 36 new ref points with a disposal that is always well-
conditioned so it can calculate the final third order polys that goes
directly form lat/lon to X/Y and viceversa"
----
/endquote
the tests are:
{{{
/*
PolyEleEval - by Marco Certelli - License: GPL - Version 01beta
This function evaluates the elements of the polynomial that
approximates a set
of ref points placed in the plane xy. This is specific for XY<->EN
transformations
needed in Chart Georeferencing. One approximating poly is defined as:
N = a[0] + a[1]X + a[2]Y + a[3]X^2 + a[4]XY + a[5]Y^2 + a[6]X^3 +
a[7]X^2Y + a[8]XY^2 + a[9]Y^3
E = b[0] + b[1]X + b[2]Y + b[3]X^2 + b[4]XY + b[5]Y^2 + b[6]X^3 +
b[7]X^2Y + b[8]XY^2 + b[9]Y^3
This function takes the int nref, representing number of available ref
points,
2 vectors of double x and y, representing the ref points position
(either
x, y or lat, lon), and returns a vector of int e[10] which i-th value
[1 or 0]
states if the relative polynomial coefficients a[i] is well or ill-
conditioned.
The x and y input values shall be between 0.0 and 1.0 (the caller
shall
scale data). The function itself returns an int 1 if an error
occurred.
*/
#define N_POLY_TEST 23
{
int i,j,n,t,ne,optt;
double *A[2*MAX_REF_POINTS+2], B1[2*MAX_REF_POINTS+2],
B2[MAX_REF_POINTS], C1[20], C2[10];
double eqmN,eqmE,sN,sE,rcond,sC,prjerror,ill,opt,minopt;
// n X Y X X Y X X X Y a
// X Y Y X X Y Y =
// X Y Y Y b
int test[N_POLY_TEST][11] =
{3,1,1,0,0,0,0,0,0,0,1, // 0 - P=C + X + Y
& (Cxx=Cyy, Cxy=-Cyx)
3,1,1,0,0,0,0,0,0,0,0, // 1 - P=C + X + Y
4,1,1,1,0,0,0,0,0,0,1, // 2 - P=C + X + Y + X^2
& (Cxx=Cyy, Cxy=-Cyx)
4,1,1,0,1,0,0,0,0,0,1, // 3 - P=C + X + Y + XY
& (Cxx=Cyy, Cxy=-Cyx)
4,1,1,0,0,1,0,0,0,0,1, // 4 - P=C + X + Y + Y^2
& (Cxx=Cyy, Cxy=-Cyx)
4,1,1,1,0,0,0,0,0,0,0, // 5 - P=C + X + Y + X^2
4,1,1,0,1,0,0,0,0,0,0, // 6 - P=C + X + Y + XY
4,1,1,0,0,1,0,0,0,0,0, // 7 - P=C + X + Y + Y^2
5,1,1,1,1,0,0,0,0,0,0, // 8 - P=C + X + Y + X^2 +
XY
5,1,1,0,1,1,0,0,0,0,0, // 9 - P=C + X + Y + XY +
Y^2
5,1,1,1,0,1,0,0,0,0,0, // 10 - P=C + X + Y + X^2 +
Y^2
6,1,1,1,1,1,0,0,0,0,0, // 11 - P=C + X + Y + X^2 +
XY + Y^2
6,1,1,1,1,0,0,1,0,0,0, // 12 - P=C + X + Y + X^2 +
XY + X^2Y
6,1,1,0,1,1,0,0,1,0,0, // 13 - P=C + X + Y + XY +
Y^2+ Y^2X
7,1,1,1,1,1,1,0,0,0,0, // 14 - P=C + X + Y + X^2 +
XY + Y^2 + X^3
7,1,1,1,1,1,0,1,0,0,0, // 15 - P=C + X + Y + X^2 +
XY + Y^2 + X^2Y
7,1,1,1,1,1,0,0,1,0,0, // 16 - P=C + X + Y + X^2 +
XY + Y^2 + XY^2
7,1,1,1,1,1,0,0,0,1,0, // 17 - P=C + X + Y + X^2 +
XY + Y^2 + Y^3
8,1,1,1,1,1,0,1,1,0,0, // 18 - P=C + X + Y + X^2 +
XY + Y^2 + X^2Y+ XY^2
8,1,1,1,1,1,1,0,0,1,0, // 19 - P=C + X + Y + X^2 +
XY + Y^2 + X^3 + Y^3
9,1,1,1,1,1,1,1,1,0,0, // 20 - P=C + X + Y + X^2 +
XY + Y^2 + X^3 + X^2Y+ XY^2
9,1,1,1,1,1,0,1,1,1,0, // 21 - P=C + X + Y + X^2 +
XY + Y^2 + X^2Y+ XY^2+ Y^3
10,1,1,1,1,1,1,1,1,1,0};// 22 - P=C + X + Y + X^2 +
XY + Y^2 + X^3 + X^2Y+ XY^2 + Y^3
int element[N_POLY_TEST][10];
}}}
I am wondering if this is useful for us but I am reminded that in stats
that the simplest model fit that explains the data is usually the best one
to use. You can use high nth-order polynomials to exactly fit the data,
but what you are really doing is fitting the noise.
So for example I have a good & undistorted scan of a printed map which I
wish to georectify. I can get a smaller RMS by using order=3, but that's
just fitting my measurement error better. It's more appropriate to use
order=2 which will average out my measurement error instead of modelling
the noise -- which is what I ''really'' want.
comments?
Hamish
--
Ticket URL: <https://trac.osgeo.org/grass/ticket/929#comment:4>
GRASS GIS <http://grass.osgeo.org>