This is splin2.c in view mode; [Download] [Up]
/* 2-D bicubic spline */
#include <stdio.h>
#include <stdlib.h>
#include "splin2.h"
#define APPNAME "ContourView"
/* 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.
*/
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,n);
yytmp=fvector(1,n);
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,n);
free_fvector(ytmp,1,n);
}
/* 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);
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);
}
/* memory allocation routines */
/* ###### 1-d float */
float *fvector(int nl, int nh)
{
float *v;
v=(float *)malloc((unsigned)(nh-nl+1)*sizeof(float));
if(!v) fprintf(stderr, "allocation failure in fvector()\n");
return(v-nl);
}
void free_fvector(float *v, int nl, int nh)
{
free((void *)(v+nl));
}
/* ##### 2-d float */
float **fmatrix(int nrl, int nrh, int ncl, int nch)
{
int i=2,j;
int error=0;
float **m;
m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float *));
if (!m) fprintf(stderr, "allocation failure 1 in fmatrix()\n");
else{
m -= nrl;
i=nrl;
do{
m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float));
if (!m[i]) {
fprintf(stderr, "allocation failure 2 in fmatrix()\n");
error=1;
}else
m[i] -= ncl;
i++;
}while( (i<=nrh) && (!error));
}
if(error){
fprintf(stderr, "Freeing allocated submatrices...\n");
for(j=i-2;j>=nrl;j--) /* free those allocated */
free((void *)(m[j]+ncl));
free((void *)(m+nrl));
m=NULL;
}
return m;
}
void free_fmatrix(float **m, int nrl, int nrh, int ncl, int nch)
{
int i;
for(i=nrh;i>=nrl;i--)free((void *)(m[i]+ncl));
free((void *) (m+nrl));
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.