[gdal-dev] Interesting issue with ogr2ogr and -tps

[Martin.Franzke at telekom.de]

>
>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-method-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_pivoting 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: thinplatespline.patch
Type: application/octet-stream
Size: 2136 bytes
Desc: thinplatespline.patch
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20180518/578000c4/attachment.obj>