/* @(#)mo_msi.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 .IDENTIFIER mo_msi.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 ------------------------------------------------------------*/ /* * 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 #define FNROWS 5 /* maximum number or rows involved in boundary extension low side */ #define LNROWS 7 /* maximum number of rows involved in high side boundary extension */ /* MO_MSIINIT -- Procedure to initialize the sewquential 2D image interpolation package. MSIINIT checks that the interpolant is one of the permitted types and allocates space for the interpolant descriptor structure. MO_MSIINIT returns the pointer to the interpolant descriptor structure. */ MO_MSIINIT(interp_type, nxpix, nypix) int interp_type; /* interpolant type */ int nxpix; int nypix; { if (interp_type < 1 || interp_type > MO_NINTERP) SCETER(66,"*** FATAL: Illegal interpolant"); else { MO_MSI_TYPE = interp_type; MO_MSI_COEFF = (float *) osmmget((MAXPIX+2) * (MAXPIX+2) * sizeof( float )); *MO_MSI_COEFF = MO_NULL; } } /* MSIFIT -- MSIFIT calculates the coefficients of the interpolant. With the exception of the bicubic spline interpolant the coefficients are stored as the data points. The 2D B-spline coefficients are calculated using the routines II_SPLINE2D. MSIFIT checks that the dimensions of the data array are appropriate for the interpolant selected and allocates space for the coefficient array. Boundary extension is performed using boundary projection. */ MO_MSIFIT(datain, nxpix, nypix, len_datain) float *datain; /* data array */ int nxpix; /* number of points in the x dimension */ int nypix; /* number of points in the y dimension */ int len_datain; /* row length of datain */ { int i, j, k; int fptr, rptr, nptr; int rptrf[FNROWS]; int rptrl[LNROWS]; float *tmp; if (len_datain < nxpix) SCETER(66,"*** FATAL: Row length of datain too small."); /* Check that the number of data points in x and y is appropriate for the interpolant type selected and allocate space for the coefficient array allowing sufficient storage for boundary extension */ switch (MO_MSI_TYPE) { case MO_BINEAREST: if (nxpix < 1 || nypix < 1) SCETER(66,"*** FATAL: Too few data points binearest fit."); else { MO_MSI_NXCOEFF = nxpix; MO_MSI_NYCOEFF = nypix; MO_MSI_FSTPNT = 0; if (*MO_MSI_COEFF != MO_NULL) { osmmfree((char *) MO_MSI_COEFF); MO_MSI_COEFF = (float *) osmmget((MAXPIX+2) * (MAXPIX+2) * sizeof( float )); } } break; case MO_BILINEAR: if (nxpix < 2 || nypix < 2) SCETER(66,"*** FATAL: Too few data points for bilinear fit."); else { MO_MSI_NXCOEFF = nxpix + 1; MO_MSI_NYCOEFF = nypix + 1; MO_MSI_FSTPNT = 0; if (*MO_MSI_COEFF != MO_NULL) { osmmfree((char *) MO_MSI_COEFF); MO_MSI_COEFF = (float *) osmmget((MAXPIX+2) * (MAXPIX+2) * sizeof( float )); } } break; case MO_BIPOLY3: if (nxpix < 4 || nypix < 4) SCETER(66,"*** FATAL: Too few data points for bipolynomial_3 fit."); else { MO_MSI_NXCOEFF = nxpix + 3; MO_MSI_NYCOEFF = nypix + 3; MO_MSI_FSTPNT = MO_MSI_NXCOEFF + 1; if (*MO_MSI_COEFF != MO_NULL) { osmmfree((char *) MO_MSI_COEFF); MO_MSI_COEFF = (float *) osmmget((MAXPIX+2) * (MAXPIX+2) * sizeof( float )); } } break; case MO_BIPOLY5: if (nxpix < 6 || nypix < 6) SCETER(66,"*** FATAL: Too few data points for bipolynomial_5 fit."); else { MO_MSI_NXCOEFF = nxpix + 5; MO_MSI_NYCOEFF = nypix + 5; MO_MSI_FSTPNT = 2 * MO_MSI_NXCOEFF + 2; if (*MO_MSI_COEFF != MO_NULL) { osmmfree((char *) MO_MSI_COEFF); MO_MSI_COEFF = (float *) osmmget((MAXPIX+2) * (MAXPIX+2) * sizeof( float )); } } break; case MO_BISPLINE3: if (nxpix < 4 || nypix < 4) SCETER(66,"*** FATAL: Too few data points for bispline_3 fit."); else { MO_MSI_NXCOEFF = nxpix + 3; MO_MSI_NYCOEFF = nypix + 3; MO_MSI_FSTPNT = MO_MSI_NXCOEFF + 1; if (*MO_MSI_COEFF != MO_NULL) { osmmfree((char *) MO_MSI_COEFF); MO_MSI_COEFF = (float *) osmmget((MAXPIX+2) * (MAXPIX+2) * sizeof( float )); } } } /* Index the coefficient pointer so that MO_COEFF(fptr+1) points to the first data point in the coefficient array */ fptr = MO_MSI_FSTPNT - 1; /* index first element */ /* Load data into coefficient array */ rptr = fptr; for (j = 0; j < nypix; j++) { for (k = 0; k < nxpix; k++) *(MO_MSI_COEFF+rptr+1+k) = *(datain+j*len_datain+k); rptr += MO_MSI_NXCOEFF; } /* Calculate the coefficients of the interpolant boundary extension is performed using boundary projection */ switch (MO_MSI_TYPE) { case MO_BINEAREST: /* No end conditions necessary, coefficients stored as data */ break; case MO_BILINEAR: /* Extend the rows */ rptr = fptr + nxpix; for (j = 0; j < nypix; j++) { *(MO_MSI_COEFF+rptr+1) = 2. * *(MO_MSI_COEFF+rptr) - *(MO_MSI_COEFF+rptr-1); rptr += MO_MSI_NXCOEFF; } /* Define the pointers to the last, 2nd last and third last rows */ rptrl[0] = (MO_MSI_NYCOEFF - 1) * MO_MSI_NXCOEFF; for (i = 1; i < 3; i++) rptrl[i] = rptrl[i-1] - MO_MSI_NXCOEFF; /* Define the last row by extending the columns */ for (k = 0; k < MO_MSI_NXCOEFF; k++) *(MO_MSI_COEFF+rptrl[0]+k) = 2. * *(MO_MSI_COEFF+rptrl[1]+k) - *(MO_MSI_COEFF+rptrl[2]+k); break; case MO_BIPOLY3: /* Extend the rows */ rptr = fptr; nptr = fptr + nxpix; for (j = 0; j < nypix; j++) { *(MO_MSI_COEFF+rptr) = 2. * *(MO_MSI_COEFF+rptr+1) - *(MO_MSI_COEFF+rptr+2); *(MO_MSI_COEFF+nptr+1) = 2. * *(MO_MSI_COEFF+nptr) - *(MO_MSI_COEFF+nptr-1); *(MO_MSI_COEFF+nptr+2) = 2. * *(MO_MSI_COEFF+nptr) - *(MO_MSI_COEFF+nptr-2); rptr = rptr + MO_MSI_NXCOEFF; nptr = nptr + MO_MSI_NXCOEFF; } /* Define pointers to first, second and third row */ rptrf[0] = 0; for (i = 1; i < 3; i++) rptrf[i] = rptrf[i-1] + MO_MSI_NXCOEFF; /* Extend the columns, define first row */ for (k = 0; k < MO_MSI_NXCOEFF; k++) *(MO_MSI_COEFF+rptrf[0]+k) = 2 * *(MO_MSI_COEFF+rptrf[1]+k) - *(MO_MSI_COEFF+rptrf[2]+k); /* Define pointers to last to the fifth last row */ rptrl[0] = (MO_MSI_NYCOEFF - 1) * MO_MSI_NXCOEFF; for (i = 1; i < 5; i++) rptrl[i] = rptrl[i-1] - MO_MSI_NXCOEFF; /* Extend the columns, define 2nd last row */ for (k = 0; k < MO_MSI_NXCOEFF; k++) *(MO_MSI_COEFF+rptrl[1]+k) = 2 * *(MO_MSI_COEFF+rptrl[2]+k) - *(MO_MSI_COEFF+rptrl[3]+k); /* Extend the columns, define last row */ for (k = 0; k < MO_MSI_NXCOEFF; k++) *(MO_MSI_COEFF+rptrl[0]+k) = 2 * *(MO_MSI_COEFF+rptrl[2]+k) - *(MO_MSI_COEFF+rptrl[4]+k); break; case MO_BIPOLY5: /* Extend the rows */ rptr = fptr; nptr = fptr + nxpix; for (j = 0; j < nypix; j++) { *(MO_MSI_COEFF+rptr-1) = 2. * *(MO_MSI_COEFF+rptr+1) - *(MO_MSI_COEFF+rptr+3); *(MO_MSI_COEFF+rptr) = 2. * *(MO_MSI_COEFF+rptr+1) - *(MO_MSI_COEFF+rptr+2); *(MO_MSI_COEFF+nptr+1) = 2. * *(MO_MSI_COEFF+nptr) - *(MO_MSI_COEFF+nptr-1); *(MO_MSI_COEFF+nptr+2) = 2. * *(MO_MSI_COEFF+nptr) - *(MO_MSI_COEFF+nptr-2); *(MO_MSI_COEFF+nptr+3) = 2. * *(MO_MSI_COEFF+nptr) - *(MO_MSI_COEFF+nptr-3); rptr = rptr + MO_MSI_NXCOEFF; nptr = nptr + MO_MSI_NXCOEFF; } /* Define pointers to first five rows */ rptrf[0] = 0; for (i = 1; i < 5; i++) rptrf[i] = rptrf[i-1] + MO_MSI_NXCOEFF; /* Extend the columns, define first row */ for (k = 0; k < MO_MSI_NXCOEFF; k++) *(MO_MSI_COEFF+rptrf[0]+k) = 2 * *(MO_MSI_COEFF+rptrf[2]+k) - *(MO_MSI_COEFF+rptrf[4]+k); /* Extend the columns, define second row */ for (k = 0; k < MO_MSI_NXCOEFF; k++) *(MO_MSI_COEFF+rptrf[1]+k) = 2 * *(MO_MSI_COEFF+rptrf[2]+k) - *(MO_MSI_COEFF+rptrf[3]+k); /* Define pointers last seven rows */ rptrl[0] = (MO_MSI_NYCOEFF - 1) * MO_MSI_NXCOEFF; for (i = 1; i < 7; i++) rptrl[i] = rptrl[i-1] + MO_MSI_NXCOEFF; /* Extend the columns, last row */ for (k = 0; k < MO_MSI_NXCOEFF; k++) *(MO_MSI_COEFF+rptrl[0]+k) = 2 * *(MO_MSI_COEFF+rptrl[3]+k) - *(MO_MSI_COEFF+rptrl[6]+k); /* Extend the columns, 2nd last row */ for (k = 0; k < MO_MSI_NXCOEFF; k++) *(MO_MSI_COEFF+rptrl[1]+k) = 2 * *(MO_MSI_COEFF+rptrl[3]+k) - *(MO_MSI_COEFF+rptrl[5]+k); /* Extend the columns, 3rd last row */ for (k = 0; k < MO_MSI_NXCOEFF; k++) *(MO_MSI_COEFF+rptrl[2]+k) = 2 * *(MO_MSI_COEFF+rptrl[3]+k) - *(MO_MSI_COEFF+rptrl[4]+k); break; case MO_BISPLINE3: /* allocate space for a temporary work arrays */ tmp = (float *) osmmget( MO_MSI_NXCOEFF * MO_MSI_NYCOEFF * sizeof(float)); /* The B-spline coefficients are calculated using the natural end conditions, end coefficents are set to zero */ /* Calculate the univariate B_spline coefficients in x */ II_SPLINE2D(MO_MSI_COEFF, tmp, nxpix, MO_MSI_NYCOEFF, MO_MSI_NXCOEFF, MO_MSI_NYCOEFF); /* Calculate the univariate B-spline coefficients in y to results of x interpolation */ II_SPLINE2D(tmp, MO_MSI_COEFF, nypix, MO_MSI_NXCOEFF, MO_MSI_NYCOEFF, MO_MSI_NXCOEFF); /* free the memory tmp */ osmmfree((char *) tmp); } } /* MSIVECTOR -- Procedure to evaluate the interpolant at an array of arbitrarily spaced points. The routines assume that 1 <= x <= nxpix and 1 <= y <= nypix. Checking for out of bounds pixels is the responsibility of the calling program. */ MO_MSIVECTOR(x, y, zfit, npts) 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 */ { switch (MO_MSI_TYPE) { case MO_BINEAREST: II_BINEAREST(MO_MSI_COEFF, MO_MSI_FSTPNT, MO_MSI_NXCOEFF, x, y, zfit, npts); break; case MO_BILINEAR: II_BILINEAR(MO_MSI_COEFF, MO_MSI_FSTPNT, MO_MSI_NXCOEFF, x, y, zfit, npts); break; case MO_BIPOLY3: II_BIPOLY3(MO_MSI_COEFF, MO_MSI_FSTPNT, MO_MSI_NXCOEFF, x, y, zfit, npts); break; case MO_BIPOLY5: II_BIPOLY5(MO_MSI_COEFF, MO_MSI_FSTPNT, MO_MSI_NXCOEFF, x, y, zfit, npts); break; case MO_BISPLINE3: II_BISPLINE3(MO_MSI_COEFF, MO_MSI_FSTPNT, MO_MSI_NXCOEFF, x, y, zfit, npts); } }