Preusser

The Preusser algorithm is part of the Estimate Scattered Data option in the Non-Parametric menu. It is also offered in the Interpolate Uniform Grid option.

Author
Albrecht Preusser
Fritz-Haber-Inst. der MPG
Faradayweg 4-6
D-14195 Berlin
Germany
E-mail: preusser@fhi-berlin.mpg.de
web: http://www.fhi-berlin.mpg.de/grz/pub/preusser.html

References
"Efficient Formulation of a Bivariate Nonic C2-Hermite Polynomial on Triangles", ACM Transactions on Mathematical Software, Vol. 16, No. 3, Sept. 1990, p246-252.
"Algorithm 684: C1 and C2-Interpolation on Triangles with Quintic and Nonic Bivariate Polynomials", ACM Transactions on Mathematical Software, Vol. 16, No. 3, Sept. 1990, p253-257.

Description
For the Preusser algorithm, Albrecht Preusser's ACM 684 procedure is implemented. This algorithm's innovation is a fast C2 (twice continuously differentiable) interpolant, and is the only procedure among the scattered data algorithms to offer a triangle-based interpolant with this property. The algorithm also offers a C1(once continuously differentiable) quintic interpolant identical to that used in the Akima I and Akima II procedures. The C1 interpolant is a 21 coefficient Hermite bivariate quintic polynomial while the C2 interpolant is a 55 coefficient Hermite bivariate nonic. The algorithm has been supplemented with a Renka I extrapolation procedure.

Revisions to Algorithm
The ACM 624 triangulation used in this algorithm was replaced with the ACM 751 TRIPACK routines. The gradient estimation procedure was replaced with a modified version of the ACM 752 SRFPACK global gradient estimation routine. As a result of these revisions, the first two gradients at each node used in this procedure will exactly match those used in the Renka I option. Extrapolation was added by using a modified form of the C1 cubic extrapolant in ACM 752. The data are also intrinsically scaled to a unit square.

Interpolation
This algorithm is based upon Hermite polynomials on triangles. There is one polynomial for each triangle, whose coefficients are determined by the partial derivatives at the nodes of the triangle mesh. For the nonic case, 14 partial derivatives are computed at each node within the convex hull boundary of the data. This interpolant consists of 55 coefficients for each triangle, 45 values at the nodes with the additional conditions assuring continuity of second partial derivatives at triangle boundaries. With 14 derivatives per node, 42 partial derivatives are used per triangle. For the quintic case, 5 partial derivatives are computed at each node. The quintic interpolant consists of 21 coefficients for every triangle, 18 values at the nodes with 3 additional conditions that assure continuity of first partial derivatives across triangle edges. Here 15 total partial derivatives are used per triangle. The gradients are computed using Robert Renka's global gradient procedure.

Extrapolation
The algorithm's added extrapolation uses the first two partials at the applicable boundary nodes to maintain a C1 boundary at the transition to the extrapolation regions.

Estimated Partial Derivatives
In general, the Preusser C2 procedure will produce the most accurate first order partial derivatives of the triangulation algorithms since these consist of analytic derivatives of the C2 nonic interpolant. The first partial derivatives from the C2 algorithm will be smooth within the data region and continuous at the interpolation-extrapolation boundary and beyond. The second partials are also analytic derivatives that exhibit a continuous surface within the data region, although smoothness and accuracy will be lacking. Estimated partial derivatives for the Preusser C1 procedure are done using analytic derivatives of the C1 quintic interpolant. The first order partials of the C1 interpolant will be continuous but not smooth; the second order partials will lack continuity. The estimated derivatives are used for graphing the partial derivative surfaces and for partial derivatives computed in the Evaluation procedure. The five partial derivatives offered by TableCurve Studio will consist exclusively of analytic derivatives of the interpolant or extrapolant.

Estimated Volumes
The primary method is to interpolate a uniform grid of 10,000 points and to evaluate the double integral using a bicubic B-spline integration procedure. The error reported in the Evaluation procedure will be the fractional difference with a similar 2,500 node integration. A repeat evaluation with the same limits will use an actual double integration procedure using the actual interpolant. A very fast double Gaussian quadrature procedure is first attempted to 1E-5 precision. If this is unsuccessful, a double adaptive quadrature procedure is then used.

Considerations
For scattered data, the Preusser C2 is the only triangulation algorithm that directly produces smooth first derivatives directly from the interpolant. The Preusser C2 algorithm also offers the greatest smoothness of the scattered data triangulation-based interpolants. While the algorithms are still subject to thin triangle drawbacks and the non-unique nature of Delaunay triangulations, the global gradient computations aid in preserving overall shape properties. Note also that the fundamental difference in how gradients are computed with this procedure will result in significantly different surfaces when the Preusser C1 quintic is compared with either the Akima I or Akima II C1 quintic procedures. On the other hand, due to identical first order gradient computation procedures, as well as a common extrapolant, the differences between the C1 and C2 Preusser algorithms and the Renka I procedure will be subtle. For the basic surfaces, as opposed to partial derivatives, the Renka I, Preusser C1, and Preusser C2 procedures can be regarded as a single algorithm with a variable order interpolant.

Algorithm Adjustments
This algorithm's only adjustment is the selection of the C1 quintic or the C2 nonic interpolant.

Choosing C1 or C2 Both choices pass exactly through every data point, so the choice concerns smoothness, not the quality of the fit. The C1 interpolant is continuous in its value and its first derivatives; the C2 interpolant is continuous through the second derivatives as well. To view, contour or export a surface, C1 is sufficient: a discontinuity in curvature cannot be seen on a shaded surface, and C1 requires less computation and memory.

The difference appears when you plot partial derivatives. The first partial derivatives of a C1 surface are continuous but not smooth, and its second partial derivatives are not continuous across triangle edges, so a second-derivative plot can show the triangulation as a faint network of seams. These seams are artifacts of the mesh, not features of the data. Choose C2 when you intend to plot or evaluate second partial derivatives, or when you need smooth first partial derivatives.

The price of C2 is a larger interpolant: 55 coefficients per triangle for the C2 nonic against 21 for the C1 quintic, and 14 partial derivatives at each node against 5. For very large scattered data sets, C1 is the lighter choice.

Because both procedures are interpolants, noise in the data is reproduced in the surface. For noisy data, smooth the data first with the Loess procedure. Values beyond the convex hull of the data come from the extrapolant described above and should be treated with caution.