#929: g.transform: dump coeffs and transform sparse points --------------------------+------------------------------------------------- Reporter: hamish | Owner: [email protected] 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>
_______________________________________________ grass-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-dev
