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

[David Fogel]

As I have indicated previously, a change in method beyond a more
robust linear solver is required.  I have taken the time to fool with
Armadillo.  At the very least I can pass along a trick to improve the
numerical accuracy.. some.

It would help me if you would share the data because I could actually
generate code and an example with real data.  At the same time, I have
never used the TPS option with gdalwarp so it would speed the learning
curve for me (though it looks simple enough).

If you are uncomfortable sharing all the data, the considering sending
the points.  Even as an attachment.  I can do a check to see if you
have a chance with Armadillo w/o additional work.

Please understand that I am not trying to be difficult, but rather
efficient.  I will share comment with Even R. and Frank W. later
Monday or Tuesday with enough information for some triage.  A more
robust solution will follow.  I have the code for this problem that
merely needs some additional notes so that we can all move ahead in a
productive manner.

Best
= David

On Sat, Sep 17, 2011 at 4:10 AM, Jan Hartmann <j.l.h.hartmann at uva.nl> wrote:
>
>
> On 09/09/2011 08:11 AM, Even Rouault wrote:
>
> 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.
>
> I tried to compile the SVN version a few times the last couple of days, but
> it always broke down: Can someone do a "make clean" in their sandbox and
> repair the error?
>
>  g++ -fPIC -L/home/a907hart/local/lib -L/usr/lib  gdaldem.o commonutils.o
> -L/home/a907hart/src/gdal -lgdal  -larmadillo -L/home/a907hart/local/lib
> -lgeos_c -L/home/a907hart/local/pgsql/lib -lpq -lpthread -lm -lrt -ldl
> -L/home/a907hart/local/lib   -lcurl            -o gdaldem
> gdaldem.o:(.data.rel.ro._ZTV21GDALGeneric3x3Dataset[vtable for
> GDALGeneric3x3Dataset]+0x30): undefined reference to
> `GDALDataset::CloseDependentDatasets()'
> gdaldem.o:(.data.rel.ro._ZTV22GDALColorReliefDataset[vtable for
> GDALColorReliefDataset]+0x30): undefined reference to
> `GDALDataset::CloseDependentDatasets()'
>
>
> Jan
>
> 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
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>