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

[Big Bear]

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

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
>