/* @(#)mo_bifit.c 17.1.1.1 (ES0-DMD) 01/25/02 17:49:06 */ /*=========================================================================== Copyright (C) 1995 European Southern Observatory (ESO) This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, MA 02139, USA. Corresponding concerning ESO-MIDAS should be addressed as follows: Internet e-mail: midas@eso.org Postal address: European Southern Observatory Data Management Division Karl-Schwarzschild-Strasse 2 D 85748 Garching bei Muenchen GERMANY ===========================================================================*/ /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ .COPYRIGHT (c) 1995 European Southern Observatory .Copyright (c) 1986 Association of Universities for Research in Astronomy Inc. .IDENTIFIER mo_bifit.c .AUTHOR R.H. Warmels IPG-ESO Garching .KEYWORDS fitting of 2dim data .LANGUAGE C .PURPOSE Routines for fitting the mosaics .ENVIRONment MIDAS #include Symbols used by the ccd package .VERSION 1.0 11-Jul-1995 creation .NOTE ------------------------------------------------------------*/ /* * Define _POSIX_SOURCE to indicate * that this is a POSIX program */ #define _POSIX_SOURCE 1 /* * definition of the used functions in this module */ #include #include #include #include #include /* II_BINEAREST -- Procedure to evaluate the nearest neighbour interpolant. The real array coeff contains the coefficients of the 2D interpolant. The procedure assumes that 1 <= x <= nxpix and 1 <= y <= nypix and that coeff[1+first_point] = datain[1,1]. */ II_BINEAREST(coeff, first_point, len_coeff, x, y, zfit, npts) float *coeff; /* 1D coefficient array */ int first_point; /* offset of first data point */ int len_coeff; /* row length of coeff */ float *x; /* array of x values */ float *y; /* array of y values */ float *zfit; /* array of interpolated values */ int npts; /* number of points to be evaluated */ { int nx, ny; int index; int i; for (i = 0; i < npts; i++) { nx = x[i] + 0.5; ny = y[i] + 0.5; /* Define pointer to data[nx,ny] */ index = first_point + (ny - 1) * len_coeff + nx; zfit[i] = coeff[index]; } } /* II_BILINEAR -- Procedure to evaluate the bilinear interpolant. The real array coeff contains the coefficients of the 2D interpolant. The procedure assumes that 1 <= x <= nxpix and 1 <= y <= nypix and that coeff[1+first_point] = datain[1,1]. */ II_BILINEAR(coeff, first_point, len_coeff, x, y, zfit, npts) float *coeff; /* 1D coefficient array */ int first_point; /* offset of first data point */ int len_coeff; /* row length of coeff */ float *x; /* array of x values */ float *y; /* array of y values */ float *zfit; /* array of interpolated values */ int npts; /* number of points to be evaluated */ { int nx, ny; int index; int i; float sx, sy, tx, ty; for (i = 0; i < npts; i++) { nx = x[i]; ny = y[i]; sx = x[i] - nx; tx = 1. - sx; sy = y[i] - ny; ty = 1. - sy; /* define pointer to data[nx,ny] */ index = first_point + (ny - 1) * len_coeff + nx; zfit[i] = tx * ty * coeff[index] + sx * ty * coeff[index + 1] + sy * tx * coeff[index+len_coeff] + sx * sy * coeff[index+len_coeff+1]; } } /* II_BIPOLY3 -- Procedure to evaluate the bicubic polynomial interpolant. The real array coeff contains the coefficients of the 2D interpolant. The procedure assumes that 1 <= x <= nxpix and 1 <= y <= nypix and that coeff[1+first_point] = datain[1,1]. The interpolant is evaluated using Everett's central difference formula. */ II_BIPOLY3(coeff, first_point, len_coeff, x, y, zfit, npts) float *coeff; /* 1D coefficient array */ int first_point; /* offset of first data point */ int len_coeff; /* row length of coeff */ float *x; /* array of x values */ float *y; /* array of y values */ float *zfit; /* array of interpolated values */ int npts; /* number of points to be evaluated */ { int nxold, nyold, nx, ny; int first_row, index; int i, j; float sx, tx, sx2m1, tx2m1, sy, ty; float cd20[4], cd21[4], ztemp[4]; float cd20y, cd21y; nxold = -1; nyold = -1; for (i = 0; i < npts; i++) { nx = x[i]; sx = x[i] - nx; tx = 1. - sx; sx2m1 = sx * sx - 1.; tx2m1 = tx * tx - 1.; ny = y[i]; sy = y[i] - ny; ty = 1. - sy; /* Calculate value of pointer to data[nx,ny-2] */ first_row = first_point + (ny - 2) * len_coeff + nx; /* Calculate the central differences in x at each value of y */ index = first_row; if (nx != nxold || ny != nyold) { for (j = 0; j < 4; j++) { cd20[j] = 1./6. * (coeff[index+1] - 2. * coeff[index] + coeff[index-1]); cd21[j] = 1./6. * (coeff[index+2] - 2. * coeff[index+1] + coeff[index]); index = index + len_coeff; } } /* Interpolate in x at each value of y */ index = first_row; for (j = 0; j < 4; j++) { ztemp[j] = sx * (coeff[index+1] + sx2m1 * cd21[j]) + tx * (coeff[index] + tx2m1 * cd20[j]); index = index + len_coeff; } /* Central differences in y */ cd20y = 1./6. * (ztemp[2] - 2. * ztemp[1] + ztemp[0]); cd21y = 1./6. * (ztemp[3] - 2. * ztemp[2] + ztemp[1]); /* Interpolate in y */ zfit[i] = sy * (ztemp[2] + (sy * sy - 1.) * cd21y) + ty * (ztemp[1] + (ty * ty - 1.) * cd20y); nxold = nx; nyold = ny; } } /* II_BIPOLY5 -- Procedure to evaluate a biquintic polynomial. The real array coeff contains the coefficents of the 2D interpolant. The routine assumes that 1 <= x <= nxpix and 1 <= y <= nypix and that coeff[1+first_point] = datain[1,1]. The interpolant is evaluated using Everett's central difference formula. */ II_BIPOLY5(coeff, first_point, len_coeff, x, y, zfit, npts) float *coeff; /* 1D coefficient array */ int first_point; /* offset of first data point */ int len_coeff; /* row length of coeff */ float *x; /* array of x values */ float *y; /* array of y values */ float *zfit; /* array of interpolated values */ int npts; /* number of points to be evaluated */ { int nxold, nyold, nx, ny; int first_row, index; int i, j; float sx, sx2, sx2m1, sx2m4, tx, tx2, tx2m1, tx2m4, sy, sy2, ty, ty2; float cd20[6], cd21[6], cd40[6], cd41[6], ztemp[6]; float cd20y, cd21y, cd40y, cd41y; nxold = -1; nyold = -1; for (i = 0; i < npts; i++) { nx = x[i]; sx = x[i] - nx; sx2 = sx * sx; sx2m1 = sx2 - 1.; sx2m4 = sx2 - 4.; tx = 1. - sx; tx2 = tx * tx; tx2m1 = tx2 - 1.; tx2m4 = tx2 - 4.; ny = y[i]; sy = y[i] - ny; sy2 = sy * sy; ty = 1. - sy; ty2 = ty * ty; /* Calculate value of pointer to data[nx,ny-2] */ first_row = first_point + (ny - 3) * len_coeff + nx; /* Calculate the central differences in x at each value of y */ index = first_row; if (nx != nxold || ny != nyold) { for (j = 0; j < 6; j++) { cd20[j] = 1./6. * (coeff[index+1] - 2. * coeff[index] + coeff[index-1]); cd21[j] = 1./6. * (coeff[index+2] - 2. * coeff[index+1] + coeff[index]); cd40[j] = 1./120. * (coeff[index-2] - 4. * coeff[index-1] + 6. * coeff[index] - 4. * coeff[index+1] + coeff[index+2]); cd41[j] = 1./120. * (coeff[index-1] - 4. * coeff[index] + 6. * coeff[index+1] - 4. * coeff[index+2] + coeff[index+3]); index = index + len_coeff; } } /* Interpolate in x at each value of y */ index = first_row; for (j = 0; j < 6; j++) { ztemp[j] = sx * (coeff[index+1] + sx2m1 * (cd21[j] + sx2m4 * cd41[j])) + tx * (coeff[index] + tx2m1 * (cd20[j] + tx2m4 * cd40[j])); index = index + len_coeff; } /* Central differences in y */ cd20y = 1./6. * (ztemp[3] - 2. * ztemp[2] + ztemp[1]); cd21y = 1./6. * (ztemp[4] - 2. * ztemp[3] + ztemp[2]); cd40y = 1./120. * (ztemp[0] - 4. * ztemp[1] + 6. * ztemp[2] - 4. * ztemp[3] + ztemp[4]); cd41y = 1./120. * (ztemp[1] - 4. * ztemp[2] + 6. * ztemp[3] - 4. * ztemp[4] + ztemp[5]); /* Interpolate in y */ zfit[i] = sy * (ztemp[3] + (sy2 - 1.) * (cd21y + (sy2 - 4.) * cd41y)) + ty * (ztemp[2] + (ty2 - 1.) * (cd20y + (ty2 - 4.) * cd40y)); nxold = nx; nyold = ny; } } /* II_BISPLINE3 -- Procedure to evaluate a bicubic spline. The real array coeff contains the B-spline coefficients. The procedure assumes that 1 <= x <= nxpix and 1 <= y <= nypix and that coeff[1+first_point] = B-spline[2]. */ II_BISPLINE3(coeff, first_point, len_coeff, x, y, zfit, npts) float *coeff; /* 1D coefficient array */ int first_point; /* offset of first data point */ int len_coeff; /* row length of coeff */ float *x; /* array of x values */ float *y; /* array of y values */ float *zfit; /* array of interpolated values */ int npts; /* number of points to be evaluated */ { int nx, ny; int first_row, index; int i, j; float sx, tx, sy, ty; float bx[4], by[4], accum, sum; for (i = 0; i < npts; i++) { nx = x[i]; sx = x[i] - nx; tx = 1. - sx; ny = y[i]; sy = y[i] - ny; ty = 1. - sy; /* calculate the x B-splines */ bx[0] = tx * tx * tx; bx[1] = 1. + tx * (3. + tx * (3. - 3. * tx)); bx[2] = 1. + sx * (3. + sx * (3. - 3. * sx)); bx[3] = sx * sx * sx; /* Calculate the y B-splines */ by[0] = ty * ty * ty; by[1] = 1. + ty * (3. + ty * (3. - 3. * ty)); by[2] = 1. + sy * (3. + sy * (3. - 3. * sy)); by[3] = sy * sy * sy; /* Calculate the pointer to data[nx,ny-1] */ first_row = first_point + (ny - 2) * len_coeff + nx; /* evaluate spline */ accum = 0.; index = first_row; for (j = 0; j < 4; j++) { sum = coeff[index-1] * bx[0] + coeff[index] * bx[1] + coeff[index+1] * bx[2] + coeff[index+2] * bx[3]; accum = accum + sum * by[j]; index = index + len_coeff; } zfit[i] = accum; } } /* II_SPLINE2D -- This procedure calculates the univariate B-spline coefficients for each row of data. The data are assumed to be uniformly spaced with a spacing of 1. The first element of each row of data is assumed to contain the second derivative of the data at x = 1. The nxpix + 2-th element of each row is assumed to contain the second derivative of the function at x = nxpix. Therfore if each row of data contains nxpix points, nxpix+2 B-spline coefficients will be calculated. The univariate B-spline coefficients for the i-th row of data are output to the i-th column of coeff. Therefore two calls to II_SPLINE2D are required to calculate the 2D B-spline coefficients. */ II_SPLINE2D(data, coeff, nxpix, nvectors, len_data, len_coeff) float *data; /* input data array */ float *coeff; /* output array of univariate coefficients in x */ int nxpix; /* number of x data points */ int nvectors; /* number of univariate splines to calculate */ int len_data; /* row dimension of data */ int len_coeff; /* row dimension of coeff */ { int i, ic, ic_1; int j, jd; int k; float diag[MAXPIX+2]; /* calculate off-diagonal elements by Gaussian elimination */ diag[0] = -2.; diag[1] = 0.; for (i = 2; i < nxpix+2; i++) diag[i] = 1. / (4. - diag[i-1]); /* loop over the nvectors rows of input data */ jd = -1; for (j = 0; j < nvectors; j++) { /* Copy the j-th row of data to the j-th column of coeff */ for (i = 0; i < nxpix + 2; i++) *(coeff+j+i*len_coeff) = *(data+i+jd+1); /* Forward substitution */ *(coeff+j) = *(coeff+j) / 6.; *(coeff+j+len_coeff) = ( *(coeff+j+len_coeff) - *(coeff+j) ) / 6.; for (i = 2; i < nxpix+2; i++) { ic = i*len_coeff; ic_1 = (i-1)*len_coeff; *(coeff+j+ic) = diag[i] * (*(coeff+j+ic) - *(coeff+j+ic_1)); } /* Back subsitution */ *(coeff+j+(nxpix+1)*len_coeff) = ( (diag[nxpix-1] + 2.) * *(coeff+j+nxpix*len_coeff) - *(coeff+j+(nxpix-1)*len_coeff) + *(coeff+j+(nxpix+1)*len_coeff) / 6.) / (1. + diag[nxpix] * ( diag[nxpix-1] + 2.)); for (k = 2; k < nxpix+2; k++) { i = nxpix + 2 - k; *(coeff+j+i*len_coeff) = *(coeff+j+i*len_coeff) - diag[i] * *(coeff+j+(i+1)*len_coeff); } *(coeff+j) = *(coeff+j) + 2. * *(coeff+j+len_coeff) - *(coeff+j+2*len_coeff); jd += len_data; } }