[gdal-dev] Numerical instability with thin plate spline transform

[Jan Hartmann]

Not sure whether this can be considered a bug, so I give it for what it 
is worth.

I'm doing thin plate spline transformation from one set of projected 
coordinates to another. Both sets have values between -600000 and 600000 
(meters). A typical set of gcps looks like:

  -gcp 62402 -74383 181915 315371 \
  -gcp -100169 -24213 19685 366661 \
  -gcp 118323 233822 239994 623190 \
...

When transforming the the first two columns of this list with 
gdaltransform -tps, using the same gcp-list, the results should be 
exactly equal to the last two columns:

  (from gdal_tps.cpp:)
* The thin plate spline transformer produces exact transformation
  * at all control points and smoothly varying transformations between
  * control points with greatest influence from local control points.
  * It is suitable for for many applications not well modelled by polynomial
  * transformations.
  *


However, they aren't. With a few control points the differences are only 
in millimeters, but with a few hundred gcps the differences become 
meters, and above thousand points the error can get to fifty meters. The 
problem gets worse when some control points are very near to other ones. 
I was happy to discover that dividing all numbers by 1000 solves the 
problem, but there certainly is a numerical instability in the algorithm.

Of course, thin plate spline transformations normally use scan 
coordinates as input, and these are almost always small numbers. Even 
so, I can imagine problems with very large rasters and the rubber 
sheeting applications I am working with.

Jan


-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20110908/dd55132c/attachment.html