[gdal-dev] approximation quality of gcp transformation

[Even Rouault]

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:gdal-dev-bounces at lists.osgeo.org] Im Auftrag von
> Franzke, Martin Gesendet: Freitag, 18. Mai 2018 17:20
> An: gdal-dev at lists.osgeo.org
> 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