Hi, motivated by the simple transformation to improve the quality of the TPS transformation I had a look at the polynomial counterpart. It turned out that the quality of the approximation not always improves with increasing order of the polynomial as it should be case. So I tried the same trick with the polynomial warping functions. As a benchmark I implemented a solution which uses armadillo and solves the problem ||Ax - b|| ->min directly without using the normal equations ATAx=ATb, which is more stable.
To benchmark the results the mean quadratic deviation of the warped GCPs is computed. This are the results using 1681 GCPs with Gauss Krueger coordinates (UTM coordinates should show almost the same results). Order 1 2 3 4 5 ------------------------------------------- original 17.19 21.94 25.8 - - precond. 17.19 13.98 10.22 - - armadillo 17.19 13.98 10.22 9.24 2729195.82 Using spherical coordinates yields Order 1 2 3 4 5 ------------------------------------------- original 17.18 15.4 15.19 - - precond. 17.18 13.97 10.22 - - armadillo 17.18 13.97 10.22 9.24 7.26 This shows, that a scaling of the coordinates can improve the stability even more. The effect even shows up using less GCPs. With 20 GCPs and Gauss Kruger coordinates the result is as follows: Order 1 2 3 4 5 ------------------------------------------- Original 0.3674 0.3672 0.3685 - - precond. 0.3674 0.3359 0.2572 - - armadillo 0.3674 0.3359 0.2572 0.1216 5.793 I attach a patch file with the modifications and the script to compute the mean deviations. Analyzing he code I presumably found two small bugs which I have corrected. Martin ----Ursprüngliche Nachricht----- Von: gdal-dev [mailto:[email protected]] Im Auftrag von Franzke, Martin Gesendet: Freitag, 18. Mai 2018 17:20 An: [email protected] Betreff: [gdal-dev] Interesting issue with ogr2ogr and -tps > >On mardi 27 juin 2017 07:34:59 CEST Rahkonen Jukka (MML) wrote: >> Hi, >> >> There is a well written question in gis.stackexchange about how the >> -tps option in ogr2ogr behaves >> https://gis.stackexchange.com/questions/245391/ogr2org-tps-spline-met >> hod-cr >> eates-progressive-error >> >> Does anybody have a clue about what happens? > >TPS is supposed to be exact at GCP, and interpolate between them. That >said, it looks from this report like there might be a numerical >stability issue. There are 2 solvers for the TPS >transformer: a built-in that does matrix inversion "at hand" using >Gauss- Jordan elimination method (the one that is taught at school if I >remember well -:) ), and another one that uses the armadillo library >that is faster and relates underneath on well proven numerical >computation libraries (blas/lapack). I don't think osgeo4w builds use >armadillo (actually I don't see anything in the nmake.opt related to >armadillo). >That said, I see that the Gauss-Jordan implementation usd uses the "partial >pivonting" >mentionned in >https://en.wikipedia.org/wiki/Pivot_element#Partial_and_complete_pivoti >ng so this should normally account for most numerical stabilities issues. >So either the user is rather unlikely and has encountered a case where >numerical instability occur, either his test/GCPs are wrong. >And I wouldn't have anticipated that numerical stabilities issues have >such a north-to-south pattern, but perhaps... > >Even > > >-- >Spatialys - Geospatial professional services http://www.spatialys.com > > Hi, I had a similar problem when not using armadillo, but a simple transformation does the trick. The problem is that usually the distance of different gcps compared their absolute value is very small. This leads to a bad conditioned linear equation system (see https://en.wikipedia.org/wiki/Condition_number): The first column consists of 1s and the second and third of the x, y values of the gcps respectively. Now if the difference of the gcps is relatively small the second and third row are nearly colinear to the first one, hence the matrix is nearly singular. The solution is to shift the coordinate system by the mean value of all gcps and solving the modified linear equations. Since the basis function of the TPS only depends on the difference of coordinates the resulting transformation is the same. For the implementation details see the attached patch file. Martin
gdal_crs.diff
Description: gdal_crs.diff
calc_gauss_mean_dev.py
Description: calc_gauss_mean_dev.py
_______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
