Martin, would you mind submitting your patch as a pull request against https://github.com/OSGeo/gdal so that regression tests and other quality checks are automatically passed ?
Ideall you would do two commits: the first one with the bugfixes, and a second one with your improvements (would be fine as part as the same PR, as I see there would be conflicts otherwise) I'd note that the polynomial order 3 implementation (or maybe it is a fundamental problem) is suspected to have instabilities that make it not fit for practical uses (despite the issue you have identified with the outlier elimination, that was a later addition) Even > 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 -- Spatialys - Geospatial professional services http://www.spatialys.com _______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
