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

[Even Rouault]

Le vendredi 09 septembre 2011 03:23:29, Big Bear a écrit :
> The linear solver in the TPS routine is naive for any number of
> reasons.1,2,3,...  At the same time, one is going to suffer
> considerably when the number of control points is in the thousands.
> Slow to evaluate so many coefficients for each point.
> 
> It is not such a task to improve the routine.  Best to work with
> someone who is cognizant per the implementation.  Likely to be
> straightforward.  ( Been there, done that. )
> 
> Important to know how the implementation is integrated into gdalwarp
> before getting too enthusiastic :-/
> 
> Thoughts?
> 
> 1 Looks to be Gauss-Jordan w/o pivoting
> 2 Never invert a matrix unless required
> 3 No scaling of the coordinates

You might try latest GDAL trunk, in particular 
http://trac.osgeo.org/gdal/changeset/22876 .

It offers the option to use libarmadillo to speed-up matrix inversion in 
thinplatespline.cpp (speed-up when TPS are in the thousands), but perhaps as a 
side effect, you'd get also more numerical stability.

> 
> On Thu, Sep 8, 2011 at 5:23 AM, Jan Hartmann <j.l.h.hartmann at uva.nl> wrote:
> > 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
> > 
> > 
> > 
> > _______________________________________________
> > gdal-dev mailing list
> > gdal-dev at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/gdal-dev
> 
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev