[gdal-dev] approximation quality of gcp transformation

[Martin.Franzke at telekom.de]

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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: gdal_crs.diff
Type: application/octet-stream
Size: 10695 bytes
Desc: gdal_crs.diff
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20180614/fd71eb5d/attachment-0002.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: calc_gauss_mean_dev.py
Type: application/octet-stream
Size: 1718 bytes
Desc: calc_gauss_mean_dev.py
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20180614/fd71eb5d/attachment-0003.obj>