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.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
