This is splin2.c in view mode; [Download] [Up]
/* 2-D bicubic spline from Numerical Recipes in C by Press. W.H et al. */
/* ###### BUG NOTE ######
NOTE that there is a bug in Numerical Recipes (in both 1-st and 2-nd editions)
in file splin2.c which is used here with the bug fixed.
*/
#include <stdio.h>
#include <stdlib.h>
#include "splin2.h"
#include "matalloc.h"
/* #include "version.h" */
#ifndef APPNAME
#define APPNAME "ContourView"
#endif
/* float *ytmp,*yytmp; for debuggin */
/* Given x1a, x2a, ya (Data), and y2a (computed by ssplie2),
* this routine computes interpolated value y at point x1, x2 by
* bicubic spline. From Press et al., Num Rec for C.
* Original in the book is in error. See the code below.
*/
void splin2(float *x1a, float *x2a, float **ya, float **y2a,
int m, int n, float x1, float x2, float *y)
{
int j;
float *ytmp,*yytmp;
ytmp=fvector(1, m); /* was fvector(1, n) which is incorrect */
yytmp=fvector(1, m); /* was fvector(1, n) which is incorrect */
for (j=1;j<=m;j++)
splint(x2a, ya[j], y2a[j], n, x2, &yytmp[j]);
spline(x1a, yytmp, m, 1.0e30, 1.0e30, ytmp);
splint(x1a, yytmp, ytmp, m, x1, y);
free_fvector(yytmp, 1, m); /* was free_fvector(yytmp, 1, n) which is incorrect */
free_fvector(ytmp, 1, m); /* was free_fvector(ytmp, 1, n) which is incorrect */
}
/* From Press, et al. Num Rec for C.
* Returns second derivatives in array y2a[1..m][1..n]
*/
void splie2(float *x1a, float *x2a, float **ya, int m, int n, float **y2a)
{
int j;
for (j=1;j<=m;j++)
spline(x2a,ya[j],n,1.0e30,1.0e30,y2a[j]);
}
void splint(float *xa, float *ya, float *y2a, int n, float x, float *y)
{
int klo,khi,k;
float h, b, a;
klo=1;
khi=n;
while (khi-klo > 1) {
k=(khi+klo) >> 1;
if (xa[k] > x) khi=k;
else klo=k;
}
h=xa[khi]-xa[klo];
if (h == 0.0) {
fprintf(stderr,"%s: Bad XA input to routine splint()\n", APPNAME);
exit(1);
}
a=(xa[khi]-x)/h;
b=(x-xa[klo])/h;
*y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
}
void spline(float *x, float *y, int n, float yp1, float ypn, float *y2)
{
int i,k;
float p,qn,sig,un,*u;
u=fvector(1,n-1);
if (yp1 > 0.99e30)
y2[1]=u[1]=0.0;
else {
y2[1] = -0.5;
u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
}
for (i=2;i<=n-1;i++) {
sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
p=sig*y2[i-1]+2.0;
y2[i]=(sig-1.0)/p;
u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
}
if (ypn > 0.99e30)
qn=un=0.0;
else {
qn=0.5;
un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
}
y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
for (k=n-1;k>=1;k--)
y2[k]=y2[k]*y2[k+1]+u[k];
free_fvector(u,1,n-1);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.