9.2.1. class MpShepard2d
Bivariate Shepard Interpolation of Scattered Data

Class definition

To use the class in your program include the class definitions by

  #include "mpshepard2d.h"

The public members of the class are given below

  class MpShepard2d 
  {
    public:
      enum { 
  	NoError = 0,		// success
  	VectorSizeMismatch,	// vectors don't have same dimension (1..n)
  	TooFewNodes,		// 6 <= n
  	FitNodesOutOfRange,	// 5 <= nq <= min(40,n-1)
  	WeightNodesOutOfRange,	// 1 <= nw <= min(40,n-1)
  	AbsentInterpolation,	// you don't have generated the interplation
  	CollinearNodes,		// collinear nodes - interpolation fails
  	DuplicateNodes,		// duplicate nodes - interpolation fails
  	InvalidInput,		// invalid input parameters (can't happen)
  	CellsOutOfRange		// nr > 1  (can't happen)
  	WorkSizeMismatch,	// work array has wrong size  (can't happen)
      };
  
      MpShepard2d (void);
     ~MpShepard2d (void);
      void Remove (void);
  
      int Interpolate (const Vector &X, const Vector &Y, const Vector &F,
  		       int nq, int nw, int nr = 0);
      int GetValue    (double px, double py, double &q);
      int GetGradient (double px, double py, double &q, double &qx, double &qy);
  };

Public member functions

int Interpolate (const Vector &X, const Vector &Y, const Vector &F, int nq, int nw, int nr = 0);

Arguments

const Vector &X

const Vector &Y

Arrays of length n containing the cartesian coordinates of the nodes.

const Vector &F

Array of length n containing the data values in one-to-one correspondence with the nodes.

int nq

Number of data points to be used in the least squares fit for coefficients defining the nodal functions q(k). A highly recommended value is nq = 13. 5 <= nq <= min(40,n-1).

int nw

Number of nodes within (and defining) the radii of influence r(k) which enter into the weights w(k). For n sufficiently large, a recommended value is nw = 19. 1 <= nw <= min(40,n-1).

int nr = 0

Number of rows and columns in the cell grid. A rectangle containing the nodes is partitioned into cells in order to increase search efficiency. nr = sqrt(n/3) is recommended. If a value less or equal than 0 is given then nr = sqrt(n/3) is taken. The minimum value for nr is 1. If a smaller value is given then nr = 1 is taken.

return value

VectorSizeMismatch	If X or Y or F have non-conformant dimensions
TooFewNodes	If X contains less than 6 data elements
FitNodesOutOfRange	If nq < 5 or nq > min(40,n-1), n = dim(X)
WeightNodesOutOfRange	If nw < 1 or nw > min(40,n-1)
CollinearNodes	No unique solution due to collinear nodes
DuplicateNodes	Duplicate nodes were encountered
NoError	On success

Description

This method computes a set of parameters a and rsq defining a smooth (once continuously differentiable) bi- variate function q(x,y) which interpolates data values f at scattered nodes (x,y). The interpolant q may be eval- uated at an arbitrary point by the method GetValue(), and its first derivatives are computed by the method GetGradient(). The interpolation scheme is a modified quadratic Shepard method

   q = (w(1)*q(1)+w(2)*q(2)+..+w(n)*q(n))/(w(1)+w(2)+..+w(n))
  

for bivariate functions w(k) and q(k). The nodal functions are given by

   q(k)(x,y) = a(1,k)*(x-x(k))**2 + a(2,k)*(x-x(k))*(y-y(k))
  	       + a(3,k)*(y-y(k))**2 + a(4,k)*(x-x(k))
  	       + a(5,k)*(y-y(k))	  + f(k)
  

Thus, q(k) is a quadratic function which interpolates the data value at node k. Its coefficients a(,k) are obtained by a weighted least squares fit to the closest nq data points with weights similar to w(k). Note that the radius of influence for the least squares fit is fixed for each k, but varies with k. The weights are taken to be

  
   w(k)(x,y) = ( (r(k)-d(k))+ / r(k)*d(k) )**2
  

where (r(k)-d(k))+ = 0 if r(k) <= d(k) and d(k)(x,y) is the Euclidean distance between (x,y) and (x(k),y(k)). The radius of influence r(k) varies with k and is chosen so that nw nodes are within the radius. Note that w(k) is not defined at node (x(k),y(k)), but q(x,y) has limit f(k) as (x,y) approaches (x(k),y(k)).

int GetValue (double px, double py, double &q);

Arguments

double px

double py

Cartesian coordinates of the point p at which q is to be evaluated.

double &q

function value q(px,py) unless an error occured, in which case 0.0 is returned and the error flag is set.

return value

AbsentInterpolation	You don't have generated the interpolation
NoError	On success

int GetGradient (double px, double py, double &q, double &qx, double &qy);

Arguments

double px

double py

Cartesian coordinates of the point p at which q is to be evaluated.

double &q

function value q(px,py) unless an error occured, in which case 0.0 is returned and the error flag is set.

double &qx

double &qy

qx,qy = first partial derivatives of q at (px,py) unless an error occured.

return value

AbsentInterpolation	You don't have generated the interpolation
NoError	On success

Notes

Examples

References

1. Robert J. Renka, Transactions on Mathematical Software, Vol. 14, No. 2, 149 (1988), Algorithm 660, QSHEP2D: Quadratic Shepard Method for Bivariate, Interpolation of Scattered Data.