/************************************************************************** * * * Copyright (c) 1996 Michael Richmond and Richard Treffers * * * * This software may be copied and distributed for * * educational, research and not for profit services * * provided that this copyright and statement are * * included in all such copies. * * * **************************************************************************/ /* *** 2D IMAGE INTERPOLATION ROUTINES *** This file contains a set of functions to be used for interpolating between image pixels, to find the value of the image at fractional pixel locations. Functions: BINSET Set up the interpolation protocol for the other functions SINCBIN Sets up 2 1D sinc arrays for repeated row and column interpolation BINCOL Just interpolate along image columns BINROW Just interpolate along rows bilinear Use 2D bilinear interpolation XBIN Use high accuracy 2D sinc interpolation OUTBIN Use lower (but faster) sinc function Author: Tod R. Lauer 5/5/83 Lick Observatory based on routines developed by Steve Kent at MIT 7/7/94 - fill changed to int. 3/3/96 - fill changed to int16 - finished translating sinc-interpolation and lagrangian interpolation routines from FORTRAN to C - removed "outbin" program - cause "rowbin", "colbin", "sinc" to call "initsinc" */ #include "pcvista.h" #include "interp.h" #include #include #include /* * Compute interpolated value in image array at point (x,y). * This entry point uses a simple bilinear interpolation scheme, * linear interpolaton between the 4 pixels surrounding the X, Y * position. * * Note: x:columns and y:rows * * Return the value "fill" if we can't interpolate properly. */ int16 bilinear(x, y, data, nrow, ncol, fill) double x, y; int16 **data, fill; int nrow, ncol; { double rval; int ix, iy; ix = x; iy = y; if ((ix < 0) || (ix >=ncol) || (iy < 0 ) || ( iy >= nrow)) { return fill; } if (ix == (ncol - 1)) { ix--; } if (iy == (nrow - 1)) { iy--; } rval = data[iy][ix] * (iy + 1.0 - y) * (ix + 1.0 - x) + data[iy+1][ix] * (y - iy) * (ix + 1.0 - x) + data[iy+1][ix+1] * (y - iy) * (x - ix) + data[iy][ix+1] * (iy + 1.0 - y) * (x - ix); return((int16) (rval + 0.5)); } #undef UNFINISHED /* !!!!! This removes the rest of the work */ #ifdef UNFINISHED void initsinc(nrow, ncol) int nrow, ncol; { /*!!! stub */ } int16 bincol(x,y,data,nrow,ncol,fill) double x,y; int16 **data; int nrow,ncol,fill; { fprintf(stderr,"bincol is not finished\n"); exit(1); return 1; } int16 binrow(x,y,data,nrow,ncol,fill) double x,y; int16 **data; int nrow,ncol,fill; { fprintf(stderr,"binrow is not finished\n"); exit(1); return 1; } int16 sinc(x,y,data,nrow,ncol,fill) double x,y; int16 **data; int nrow,ncol,fill; { fprintf(stderr,"sinc is not finished\n"); exit(1); return 1; } int16 lagrange(x,y,data,nrow,ncol,fill) double x,y; int16 **data; int nrow,ncol,fill; { fprintf(stderr,"lagrange is not finished\n"); exit(1); return 1; } #else #define N2 51 #define N (N2/2) #define N02 25 #define N0 (N02/2) #define BETA .25 #ifndef PI #define PI 3.141593 #endif static double ar[N2], ac[N2]; static int initsinc_flag = 0; /* * set up 2 1D sinc arrays for later use * by "bincol", "binrow", and "sinc". * * Set the global variable "initsinc_flag" to 1 to signal * that we have calculate the values for the array. * * Input: DX Fractional shift in columns * DY Fractional shift in rows */ void initsinc(dx, dy) double dx, dy; { int i, ix, iy, dxn, dyn; double sumr, sumc; double sindx, sindy; double ax, ay; double px, py; double beta_1; /* * should probably allow the user to set the value * of "beta" (the amount of tapering in sinc interpolation) * from the command line. */ beta_1 = BETA; ix = (int) dx; iy = (int) dy; dx = dx - ix; dy = dy - iy; /* * make sure that both dx and dy fall in the range * [-0.5, +0.5], adjusting the integral part * if necessary. */ if (dx > 0.5) { dx -= 1.0; ix += 1; } if (dy > 0.5) { dx -= 1.0; iy += 1; } if (N % 2 == 1) { sindx = -1.0; sindy = -1.0; } else { sindx = 1.0; sindy = 1.0; } if (dx < 0.0) { sindx = -sindx; } if (dy < 0.0) { sindy = -sindy; } dxn = 0.0 - N - dx; dyn = 0.0 - N - dy; /* Generate the sinc arrays with exponential tapering */ for(i = 0; i < N2; i++) { ax = i + dxn; ay = i + dyn; if (ax == 0.0) { px = 1.0; } else if (dx == 0.0) { px = 0.0; } else { px = (sindx/ax)*(exp(-(ax*beta_1)*(ax*beta_1))); } if (ay == 0.0) { py = 1.0; } else if (dy == 0.0) { py = 0.0; } else { py = (sindy/ay)*(exp(-(ay*beta_1)*(ay*beta_1))); } ar[i] = py; ac[i] = px; sindx = -sindx; sindy = -sindy; } /* Normalize the arrays to unity */ sumr = 0.0; sumc = 0.0; for (i = 0; i < N2; i++) { sumr += ar[i]; sumc += ac[i]; } sumr /= (double) N2; sumc /= (double) N2; for (i = 0; i < N2; i++) { ar[i] /= sumr; ac[i] /= sumc; } initsinc_flag = 1; } /* * Get image value at new column location * * This function assumes that the user has already called * the function "initsinc" to create the sinc arrays "ac" and "ar". * It is used for repeated interpolation of the same fractional amount. * * Input: X Column location of new pixel * Y Row location of new pixel * * Return the value "fill" if we can't interpolate properly. */ int16 bincol(x, y, data, nrow, ncol, fill) double x, y; int nrow, ncol; int16 **data, fill; { double retval; double dx, dy; int ix, iy; int i, ixn, nx; ix = (int) x; iy = (int) y; dx = x - ix; dy = y - iy; /* * Make SURE that we have set up the sinc-interpolation arrays */ if (initsinc_flag == 0) { initsinc(dx, dy); } if (iy < 0) { iy = 0; } if (iy >= nrow) { iy = nrow - 1; } if (dx > 0.5) { ix++; } ixn = ix - N - 1; /* * Convolve the image row with the 1D sinc array, and sum to get the * new pixel value. If at the image edge, just repeat the edge pixel. */ retval = 0.0; for (i = 0; i < N2; i++) { nx = ixn + 1; if (nx < 0) { nx = 0; } if (nx >= ncol) { nx = ncol - 1; } retval += ac[i]*data[nx][iy]; } /* is it necessary to divide by some factor here? */ /* make sure that the return value lies within the valid range */ retval += 0.5; if (retval < MIN_DATA_VAL) { return(MIN_DATA_VAL); } else if (retval > MAX_DATA_VAL) { return(MAX_DATA_VAL); } else { return((int16) retval); } } /* * Get image value at new row location * * This function assumes that the user has already called * the function "initsinc" to create the sinc arrays "ac" and "ar". * It is used for repeated interpolation of the same fractional amount. * * Input: X Column location of new pixel * Y Row location of new pixel * * Return the value "fill" if we can't interpolate properly. */ int16 binrow(x, y, data, nrow, ncol, fill) double x, y; int nrow, ncol; int16 **data, fill; { int i, iyn, ny; int ix, iy; double dx, dy; double retval; ix = (int) x; iy = (int) y; dx = x - ix; dy = y - iy; /* * Make SURE that we have set up the sinc-interpolation arrays */ if (initsinc_flag == 0) { initsinc(dx, dy); } if (ix < 0) { ix = 0; } if (ix >= ncol) { ix = ncol - 1; } if (dy > 0.5) { iy++; } iyn = iy - N - 1; /* * Convolve the image row with the 1D sinc array, and sum to get the * new pixel value. If at the image edge, just repeat the edge pixel. */ retval = 0.0; for (i = 0; i < N2; i++) { ny = iyn + i; if (ny < 0) { ny = 0; } if (ny >= nrow) { ny = nrow - 1; } retval += ar[i]*data[ix][ny]; } /* is it necessary to divide by some factor here? */ /* make sure that the return value lies within the valid range */ retval += 0.5; if (retval < MIN_DATA_VAL) { return(MIN_DATA_VAL); } else if (retval > MAX_DATA_VAL) { return(MAX_DATA_VAL); } else { return((int16) retval); } } /* * High accuracy 2D interpolation routine. * * Assumes the user has already called "initsinc" to set up * the interpolation arrays "ac" and "ar". * * This function does not check for zones or negative pixels. * Use SIN(X)/X*EXP(-X**2) form. * * Return the value "fill" if we can't interpolate properly. */ int16 sinc(x, y, data, nrow, ncol, fill) double x, y; int16 **data, fill; int nrow, ncol; { int i, j, calcflag; int ix, iy, ixn, iyn, nx, ny; double dx, dy; double sum; double retval; ix = (int) x; iy = (int) y; /* * Make SURE that we have set up the sinc-interpolation arrays */ if (initsinc_flag == 0) { dx = x - ix; dy = y - iy; initsinc(dx, dy); } /* * Convolve the image array with the sinc array and sum to generate the * new pixel value at X and Y. */ retval = 0.0; calcflag = 0; ixn = ix - N - 1; iyn = iy - N - 1; for (i = 0; i < N2; i++) { sum = 0.0; for (j = 0; j < N2; j++) { nx = ixn + j; if (nx < 0) { nx = 0; } if (nx >= ncol) { nx = ncol - 1; } ny = iyn + i; if (ny < 0) { ny = 0; } if (ny >= nrow) { ny = nrow - 1; } sum += ac[j]*data[nx][ny]; calcflag = 1; } retval += ar[i]*sum; } /* * check to see if we didn't actually perform any calculations; * if so, we set the return value equal to "fill" argument. */ if (calcflag == 0) { retval = (double) fill; } /* make sure that the return value lies within the valid range */ retval += 0.5; if (retval < MIN_DATA_VAL) { return(MIN_DATA_VAL); } else if (retval > MAX_DATA_VAL) { return(MAX_DATA_VAL); } else { return((int16) retval); } } /* * Compute the interpolated value in image array at point (x,y). * This entry point uses a 4-pt Lagrangian interpolation scheme, * cubic interpolation between the 16 pixels surrounding the (x,y) * position. (Abramowitz and Stegun, p879, 25.2.13) * * Adapted from a computer program by Robert Braun * R.A.M. Walterbos - 01mar88 * MWR - 03Mar96 * * return "fill" if we can't interpolate properly. */ int16 lagrange(x, y, data, nrow, ncol, fill) double x, y; int16 **data, fill; int nrow, ncol; { int ix, iy; double px, py; double cx0, cx1, cx2, cx3, cy0, cy1, cy2, cy3; double fy0, fy1, fy2, fy3; double retval; ix = (int) x; iy = (int) y; /* test if in the array */ if ((ix < 0) || (ix >= ncol)) { return(fill); } if ((iy < 0) || (iy >= ncol)) { return(fill); } /* * edge values will get same value as adjoining pixel */ if (ix == 0) { ix = 1; } if (ix >= ncol - 2) { ix = ncol - 3; } if (iy == 0) { iy = 1; } if (iy >= nrow - 2) { iy = nrow - 3; } /* * Start interpolation; first define coefficients */ px = x - ix; py = y - iy; cx0 = -px*(px-1.0)*(px-2.0)*0.16666667; cx1 = (px*px-1.0)*(px-2.0)*0.500000; cx2 = -px*(px+1.0)*(px-2.0)*0.500000; cx3 = px*(px*px-1.0)*0.1666667; cy0 = -py*(py-1.0)*(py-2.0)*0.16666667; cy1 = (py*py-1.0)*(py-2.0)*0.500000; cy2 = -py*(py+1.0)*(py-2.0)*0.500000; cy3 = py*(py*py-1.0)*0.1666667; /* * Calculate the 4 interpolated values along rows */ fy0 = cx0*data[ix-1][iy-1] + cx1*data[ix][iy-1] + cx2*data[ix+1][iy-1] + cx3*data[ix+2][iy-1]; fy1 = cx0*data[ix-1][iy] + cx1*data[ix][iy] + cx2*data[ix+1][iy] + cx3*data[ix+2][iy]; fy2 = cx0*data[ix-1][iy+1] + cx1*data[ix][iy+1] + cx2*data[ix+1][iy+1] + cx3*data[ix+2][iy+1]; fy3 = cx0*data[ix-1][iy+2] + cx1*data[ix][iy+2] + cx2*data[ix+1][iy+2] + cx3*data[ix+2][iy+2]; /* * Finally calculate the interpolated value along the column */ retval = cy0*fy0 + cy1*fy1 + cy2*fy2 + cy3*fy3; /* make sure that the return value lies within the valid range */ retval += 0.5; if (retval < MIN_DATA_VAL) { return(MIN_DATA_VAL); } else if (retval > MAX_DATA_VAL) { return(MAX_DATA_VAL); } else { return((int16) retval); } } #endif /* FINISHED */