# include # include # include # include # include # include # include "shepard_interp_nd.h" # include "r8lib.h" /******************************************************************************/ double *shepard_interp_nd ( int m, int nd, double xd[], double zd[], double p, int ni, double xi[] ) /******************************************************************************/ /* Purpose: SHEPARD_INTERP_ND evaluates a multidimensional Shepard interpolant. Licensing: This code is distributed under the MIT license. Modified: 02 October 2012 Author: John Burkardt Reference: Donald Shepard, A two-dimensional interpolation function for irregularly spaced data, ACM '68: Proceedings of the 1968 23rd ACM National Conference, ACM, pages 517-524, 1969. Parameters: Input, int M, the spatial dimension. Input, int ND, the number of data points. Input, double XD[M*ND], the data points. Input, double ZD[ND], the data values. Input, double P, the power. Input, int NI, the number of interpolation points. Input, double XI[M*NI], the interpolation points. Output, double SHEPARD_INTERP_ND[NI], the interpolated values. */ { int i; int i2; int j; double s; double t; double *w; int z; double *zi; w = ( double * ) malloc ( nd * sizeof ( double ) ); zi = ( double * ) malloc ( ni * sizeof ( double ) ); for ( i = 0; i < ni; i++ ) { if ( p == 0.0 ) { for ( j = 0; j < nd; j++ ) { w[j] = 1.0 / ( double ) ( nd ); } } else { z = -1; for ( j = 0; j < nd; j++ ) { t = 0.0; for ( i2 = 0; i2 < m; i2++ ) { t = t + pow ( xi[i2+i*m] - xd[i2+j*m], 2 ); } w[j] = sqrt ( t ); if ( w[j] == 0.0 ) { z = j; break; } } if ( z != -1 ) { for ( j = 0; j < nd; j++ ) { w[j] = 0.0; } w[z] = 1.0; } else { for ( j = 0; j < nd; j++ ) { w[j] = 1.0 / pow ( w[j], p ); } s = 0.0; for ( j = 0; j < nd; j++ ) { s = s + w[j]; } for ( j = 0; j < nd; j++ ) { w[j] = w[j] / s; } } } zi[i] = r8vec_dot_product ( nd, w, zd ); } free ( w ); return zi; }