interp.c

Go to the documentation of this file.

#include "ds_extern.h"

/*
   2d array interpolation routine
   swiped from Darrel, converted to 'C' by Jeff
*/


/*
 version 1.0  mpifr cyber edition  22 MAY 1977.
 see nod 2 manual page m 2.1                                       
 C.G Haslam       mar  1969.                                        
 This routine interpolates the amplitude z at the point y which lies
 between the points b and c. y=0.0 corresponds to b and y=1.0
 corresponds to c. a,b,c,d are the amplitudes at 4 points spaced at
 unit intervals.                                                  
*/


static float yintpn();
static int xinterp();


static float yintpn(a, b, c, d, y)
double a, b, c, d, y;
{
  double e, f, g;

  if(y == 0.0)
    return((float)b);

  e = b-a;      /* exact copy case */
  f = b-c;
  g = e + f + f + d - c;
  g = g*y- e - f - g;
  g = g*y - a + c;

  return( (float)(0.5*g*y + b));
}

/*
  interp interpolates a 2-d array array into parray.
  uses lots of memory; PC gymnastics deleted
*/

int interp(array, nx, ny, parray, npx, npy )
float *array;
int nx, ny;
float *parray;
int npx, npy;
{
  float *one, *two, *three, *four, *save, *done;
  float *xone, *xtwo, *xthree, *xfour;
  int k, i, iy;
  double y, dy, lasty;
  double oversampy, oversampx, oversamp;

  one   = (float *)malloc( (unsigned)(sizeof(float) * npx));
  two   = (float *)malloc( (unsigned)(sizeof(float) * npx));
  three = (float *)malloc( (unsigned)(sizeof(float) * npx));
  four  = (float *)malloc( (unsigned)(sizeof(float) * npx));

  if( !one || !two || !three || !four ) {
    fprintf( stderr, "no memory for interpolation\n");
    return(1);
  }

  xone   = &array[0];
  xtwo   = &array[nx];
  xthree = &array[nx*2];

  oversampy = (npy-1.0)/(ny-1.0);
  oversampx = (npx-1.0)/(nx-1.0);
  oversamp = oversampx < oversampy ? oversampx:oversampy;

/*
  npx = 1 + (int)rint(oversamp * ((float)nx - 1.0));
  npy = 1 + (int)rint(oversamp * ((float)ny - 1.0));
 */

/*  check that we need to interpolate */

  if(fabs(oversamp - 1.0) < 0.001 ) {
     fprintf(stderr,"interpolation not needed - oversamp is %f\n", oversamp );
     return(0);
  }

  xinterp( xone,   nx, two, npx );
  xinterp( xtwo, nx, three, npx );
  xinterp( xthree,  nx, four, npx );
  lasty = 0;

/* create first row by extrapolation */

  for( i=0; i< npx; i++ )
    one[i] = two[i] - three[i] + two[i];

  for( k = 0; k<npy; k++ ) {        /* for each new row, interpolate in y */

    y  = (double)k/oversampy;
    iy = (int)y;

    if( iy != lasty ) { /* check to see if it is time to move the pointers */
      
      save = one;
      one = two;
      two = three;
      three = four;
      four = save;

      if( iy+2 < ny ) {                   /* interpolate along x */
        xfour = &array[(iy+2)*nx];
        xinterp( xfour, nx, four, npx );
      } else {                            /* otherwize extrapolate */
        for( i=0; i< npx; i++ )
          four[i] = three[i] - two[i] + three[i];
      }
      lasty = iy;
    }
    done = &parray[npx*k];
    dy = y - (double)iy;

    for( i=0 ; i<npx; i++ ) 
      done[i] = yintpn( one[i], two[i], three[i], four[i], dy );
  }

  free((char *)one);
  free((char *)two);
  free((char *)three);
  free((char *)four);

  return(0);
}

/*
  subroutine xinterp interpolates the input array x1 onto the output
  array x2, with nx2 points, oversampling the input array by oversamp
*/

static int xinterp(x1, nx1, x2, nx2)
float *x1;
int nx1;
float *x2;
int nx2;
{
  float oversamp, x, a, b, c, d, dx;
  int jx, ix, inx;
  
  oversamp = (float)(nx2-1)/(float)(nx1-1);
  for( jx=0; jx<nx2; jx++ ) {  /* jx is address of interpolated array */
    x = jx/oversamp;
    inx=(int)rint(x);
    if(fabs(x-(float)inx) < 0.0001)
      x2[jx]=x1[inx];
    else {
      ix= (int)x;
      dx=x-ix;
      b=x1[ix];
      c=x1[ix+1];
      if(ix > 0)
        a=x1[ix-1];
      else
        a=b-c+b;         /* linear extrapolation for first point */

      if(ix+2 < nx1 )
        d=x1[ix+2];
      else
        d=c-b+c;         /* linear extrapolation for last point */

      x2[jx]= yintpn(a, b, c, d, dx);
    }
  }

  return(0);
}