MORE ON THIN-PLATE SPLINES

A Fortran listing

The following code executes the thin-plate spline interpolation introduced in Section 2.2. The calling program first calls TPSOLV with arguments as follows:

IDIM – number of dimensions of the data: 1, 2, or 3

NPTS (passed in common /FORMS/) – number of landmarks

POINTS (passed in common /FORMS/) – Cartesian coordinates of the landmarks. First subscript is Cartesian coordinate, second landmark number, third the form (1 or 2). Here, dimensioned to a maximum count of 50.

The subroutine returns the coefficients wt and a1. ax,ayin entries TPCOEF (1 … NPTS) and TPCOEF (NPTS + 1 … NPTS + 3), separately for each coordinate of the data. The net integral bending energy is also returned in the value SCORE.

The subroutine TPSOLV is essentially a setup and loop over calls to a linear equation solver. That incorporated here is from the package LINPACK; one may substitute other matrix inversion routines as convenient.

Once the coefficients TPCOEF have been computed, the effect of the spline mapping on any point PT of the starting image is computed by calling TPEVAL with arguments IDIM and TPCOEF copied from TPSOLV and input point PT as shown. The image of the argument point is returned in EVAL.

Differential constraints and relaxation

The material in this section previously appeared in Bookstein (1989a). The landmarks of the scene analyzed in Figure 7.5.1 were extracted from the pair of simulated edge images (actually, doodles on a desk pad) shown in Figure A.1.1. (The width of the X's at each digitized point is presumed to encode the uncertainty of edge location.)