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.hartm...@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.hartm...@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@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev > > _______________________________________________ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev > > _______________________________________________ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev > _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev