Thin Plate Spline (TPS) for large number of tiepoints

[pacesm]

Reporter: martin.paces Date: 2013/04/05 16:01

Although the TPS interpolation method (see (1) and (2)) is superior for a relatively small set of unordered (scattered) tie-points it is not very efficient for a large number of tie-points since:

• a large system of linear equations must be solved to for initialise the the interpolation weight (system matrix: (3+N)x(3+N) full, symmetrical and generally indefinite);
• distances between x_' and each tie-point must be calculated to evaluate the interpolated function at the point '_x . The complexity of interpolation of M point for N tie-points is therefore O(MxN) as opposed to O(M) of the piece-wise polynomial methods.

The TPS method fails, e.g., in case of strip-line MERIS/RR products with dense regular rectangular grid of tie-points (~10,000 of TPs!). Please note that the alternative global polynomial interpolation methods are not an alternative as these are not able to capture the geometry with an acceptable precision.

This fact is further amplified by the naive TPS implementation in the GDAL library. Instead of using memory efficient method for solving of the TPS linear system (provided e.g. by LAPACK) GDAL rather employs its own solver based on full inversion of the system matrix. This way the initialisation requires roughly 4x more memory than necessary. The inefficient memory handling leads to memory allocation failures for a smaller number of tie-points. Fixing of the solver would allow using of more tie-poins but generally does not solve the fundamental issue of the TPS method.

The problem could be solved by implementation of a piecewise polynomial interpolation exploiting the tie-point regular grid ordering or triangulation of general scattered data.

(1) http://en.wikipedia.org/wiki/Polyharmonic_spline (2) http://en.wikipedia.org/wiki/Thin_plate_spline

[pacesm]

Date: 2013/04/22 11:21

I made a small research on the issue and I also run a couple of tests. These are the findings:

1) The real issue is the TPS initialization which require solution of the dense linear system. This requires temporarily O(N^2^) storage and the solution itself is O(N^3^) algorithm.

2) After cleaning the messy excessive memory allocation and employing faster solver (ATLAS library) I found that:

• with nowadays typical 4GiB of RAM we can safely go up to ~20,000 tie-points
• the computing time starts to be prohibitive much earlier before the memory is exceeded. For more than >1,000 tie-points the computing time of the linear solver (on my laptop with decent i7 CPU) escalates quickly above the acceptable levels:

N time 1000 0.13s 2000 1.0s 3000 3.4s 4000 8.4s 5000 19s 6000 33s 7000 63s 8000 107s 9000 160s 10000 217s

Please note that GDAL actually internally performs 2 initialization one for forward (img. to geo. coordinates) and the second for reverse (geo. to img. coordinates) interpolations. Thus the real GDAL TPS initialization will take twice of the figures above.

This leads to conclusion that registration and subsetting of a single MERIS RR product with TPS interpolation even after the GDAL fixing/optimisation would require several minutes to be calculated which is unacceptably high.

[pacesm]

Problem resolved by creation of custom GDAL version with an additional TPS based transformer.

This new transformer is a Thin Plate Spline interpolation on using now, significantly smaller, set of nodes and Least Square (LSQ) Fit of the original tie-point set.

The EOxServer has been adapted to take benefit of the new interpolation method in following 0.3 branch commits:

f31e2b3db65659978a61ce1e197460559c10c93e, 3bb5045f445830d2b61e4ced3f8f78d23bb017dd, 9f4e268b05653c8ac4378cbe2ac42798826024db, 5a0617bda30c3b73ac458e056d3a3cce0611c673, d13fd81467f2c5005cc28b5c23e9b70305b9efea, 33a0a41d3874b33bb58f8a15ad74169d68d9df23

These changes still need to be merged with the master.