3.6 Interpolation in Two or More Dimensions 127 The implementation ofequation(3.6.6).which performs a bicubic interpolation, gives back the interpolated function value and the two gradient values,and uses the above routine bcucof,is simply: #include "nrutil.h" void bcuint(float y],float yi[],float y2[],float y120],float x1l, float xlu,float x21,float x2u,float x1,float x2,float *ansy, float *ansy1,float *ansy2) Bicubic interpolation within a grid square.Input quantities are y,y1,y2,y12(as described in bcucof):x11 and xiu,the lower and upper coordinates of the grid square in the 1-direction; x21 and x2u likewise for the 2-direction;and x1,x2,the coordinates of the desired point for 8 the interpolation.The interpolated function value is returned as ansy,and the interpolated gradient values as ansy1 and ansy2.This routine calls bcucof. void bcucof(float y[],float y1],float y21],float y12],float d1, float d2,float **c); int i; ICAL float t,u,d1,d2,**c; c=matr1x(1,4,1,4); RECIPES d1=x1u-x11; d2=x2u-x21; bcucof(y,y1,y2,y12,d1,d2,c); Get the c's. if (xlu =x11 II x2u =x21)nrerror("Bad input in routine bcuint"); t=(x1-x11)/d1: Equation (3.6.4). America computer, Press. u=(x2-x21)/d2; 9 *ansy=(*ansy2)=(*ansy1)=0.0; for(1=4;1>=1;1--)[ Equation (3.6.6). 9 *ansy=t*(*ansy)+((c[i][4]*u+c[i][3])*u+c[i][2])*u+c[i][1]; ictly proh Programs *ansy2=t*(*ansy2)+(3.0*c[i][4]*u+2.0*c[i][3])*u+c[i][2]; *ansy1u*(*ansy1)+(3.0*c[4][i]*t+2.0*c[3][1])*t+c[2][1]; SCIENTIFIC *ansy1 /d1; *ansy2 /d2; to dir free_matrix(c,1,4,1,4); 18818920 COMPUTING (ISBN Higher Order for Smoothness:Bicubic Spline 10-:6211 The other common technique for obtaining smoothness in two-dimensional s60二% Numerical Recipes 43108 interpolation is the bicubic spline.Actually.this is equivalent to a special case of bicubic interpolation:The interpolating function is of the same functional form as equation(3.6.6);the values of the derivatives at the grid points are,however, (outside determined "globally"by one-dimensional splines.However,bicubic splines are Software. 首 usually implemented in a form that looks rather different from the above bicubic interpolation routines,instead looking much closer in form to the routine polin2 above:To interpolate one functional value,one performs m one-dimensional splines across the rows of the table,followed by one additional one-dimensional spline down the newly created column.It is a matter of taste (and trade-off between time and memory)as to how much of this process one wants to precompute and store Instead of precomputing and storing all the derivative information(as in bicubic interpolation),spline users typically precompute and store only one auxiliary table. of second derivatives in one direction only.Then one need only do spline evaluations (not constructions)for the m row splines;one must still do a construction and an3.6 Interpolation in Two or More Dimensions 127 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The implementation of equation (3.6.6), which performs a bicubic interpolation, gives back the interpolated function value and the two gradient values, and uses the above routine bcucof, is simply: #include "nrutil.h" void bcuint(float y[], float y1[], float y2[], float y12[], float x1l, float x1u, float x2l, float x2u, float x1, float x2, float *ansy, float *ansy1, float *ansy2) Bicubic interpolation within a grid square. Input quantities are y,y1,y2,y12 (as described in bcucof); x1l and x1u, the lower and upper coordinates of the grid square in the 1-direction; x2l and x2u likewise for the 2-direction; and x1,x2, the coordinates of the desired point for the interpolation. The interpolated function value is returned as ansy, and the interpolated gradient values as ansy1 and ansy2. This routine calls bcucof. { void bcucof(float y[], float y1[], float y2[], float y12[], float d1, float d2, float **c); int i; float t,u,d1,d2,**c; c=matrix(1,4,1,4); d1=x1u-x1l; d2=x2u-x2l; bcucof(y,y1,y2,y12,d1,d2,c); Get the c’s. if (x1u == x1l || x2u == x2l) nrerror("Bad input in routine bcuint"); t=(x1-x1l)/d1; Equation (3.6.4). u=(x2-x2l)/d2; *ansy=(*ansy2)=(*ansy1)=0.0; for (i=4;i>=1;i--) { Equation (3.6.6). *ansy=t*(*ansy)+((c[i][4]*u+c[i][3])*u+c[i][2])*u+c[i][1]; *ansy2=t*(*ansy2)+(3.0*c[i][4]*u+2.0*c[i][3])*u+c[i][2]; *ansy1=u*(*ansy1)+(3.0*c[4][i]*t+2.0*c[3][i])*t+c[2][i]; } *ansy1 /= d1; *ansy2 /= d2; free_matrix(c,1,4,1,4); } Higher Order for Smoothness: Bicubic Spline The other common technique for obtaining smoothness in two-dimensional interpolation is the bicubic spline. Actually, this is equivalent to a special case of bicubic interpolation: The interpolating function is of the same functional form as equation (3.6.6); the values of the derivatives at the grid points are, however, determined “globally” by one-dimensional splines. However, bicubic splines are usually implemented in a form that looks rather different from the above bicubic interpolation routines, instead looking much closer in form to the routine polin2 above: To interpolate one functional value, one performs m one-dimensional splines across the rows of the table, followed by one additional one-dimensional spline down the newly created column. It is a matter of taste (and trade-off between time and memory) as to how much of this process one wants to precompute and store. Instead of precomputing and storing all the derivative information (as in bicubic interpolation), spline users typically precompute and store only one auxiliary table, of second derivatives in one direction only. Then one need only do spline evaluations (not constructions) for the m row splines; one must still do a construction and an