From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- math/iminterp/msifit.x | 275 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 275 insertions(+) create mode 100644 math/iminterp/msifit.x (limited to 'math/iminterp/msifit.x') diff --git a/math/iminterp/msifit.x b/math/iminterp/msifit.x new file mode 100644 index 00000000..27d861ff --- /dev/null +++ b/math/iminterp/msifit.x @@ -0,0 +1,275 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "im2interpdef.h" +include + +# 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. + +procedure msifit (msi, datain, nxpix, nypix, len_datain) + +pointer msi # pointer to interpolant descriptor structure +real datain[len_datain,ARB] # 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 +pointer fptr, nptr, rptr +pointer tmp +pointer rptrf[FNROWS] +pointer rptrl[LNROWS] + +errchk calloc, mfree + +begin + # check the row length of datain + if (len_datain < nxpix) + call error (0, "MSIFIT: 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 (MSI_TYPE(msi)) { + + case II_BINEAREST: + + if (nxpix < 1 || nypix < 1) { + call error (0, "MSIFIT: Too few data points.") + return + } else { + MSI_NXCOEFF(msi) = nxpix + MSI_NYCOEFF(msi) = nypix + MSI_FSTPNT(msi) = 0 + if (MSI_COEFF(msi) != NULL) + call mfree (MSI_COEFF(msi), TY_REAL) + call malloc (MSI_COEFF(msi), nxpix * nypix, TY_REAL) + } + + case II_BILINEAR, II_BIDRIZZLE: + + if (nxpix < 2 || nypix < 2) { + call error (0, "MSIFIT: Too few data points.") + return + } else { + MSI_NXCOEFF(msi) = nxpix + 1 + MSI_NYCOEFF(msi) = nypix + 1 + MSI_FSTPNT(msi) = 0 + if (MSI_COEFF(msi) != NULL) + call mfree (MSI_COEFF(msi), TY_REAL) + call malloc (MSI_COEFF(msi), + MSI_NXCOEFF(msi) * MSI_NYCOEFF(msi), TY_REAL) + } + + case II_BIPOLY3: + + if (nxpix < 4 || nypix < 4) { + call error (0, "MSIFIT: Too few data points.") + return + } else { + MSI_NXCOEFF(msi) = nxpix + 3 + MSI_NYCOEFF(msi) = nypix + 3 + MSI_FSTPNT(msi) = MSI_NXCOEFF(msi) + 1 + if (MSI_COEFF(msi) != NULL) + call mfree (MSI_COEFF(msi), TY_REAL) + call malloc (MSI_COEFF(msi), + MSI_NXCOEFF(msi) * MSI_NYCOEFF(msi), TY_REAL) + } + + case II_BIPOLY5: + + if (nxpix < 6 || nypix < 6) { + call error (0, "MSIFIT: Too few data points.") + return + } else { + MSI_NXCOEFF(msi) = nxpix + 5 + MSI_NYCOEFF(msi) = nypix + 5 + MSI_FSTPNT(msi) = 2 * MSI_NXCOEFF(msi) + 2 + if (MSI_COEFF(msi) != NULL) + call mfree (MSI_COEFF(msi), TY_REAL) + call malloc (MSI_COEFF(msi), + MSI_NXCOEFF(msi) * MSI_NYCOEFF(msi), TY_REAL) + } + + case II_BISPLINE3: + + if (nxpix < 4 || nypix < 4) { + call error (0, "MSIFIT: Too few data points.") + return + } else { + MSI_NXCOEFF(msi) = nxpix + 3 + MSI_NYCOEFF(msi) = nypix + 3 + MSI_FSTPNT(msi) = MSI_NXCOEFF(msi) + 1 + if (MSI_COEFF(msi) != NULL) + call mfree (MSI_COEFF(msi), TY_REAL) + call calloc (MSI_COEFF(msi), + MSI_NXCOEFF(msi) * MSI_NYCOEFF(msi), TY_REAL) + } + + case II_BISINC, II_BILSINC: + + if (nxpix < 1 || nypix < 1) { + call error (0, "MSIFIT: Too few data points.") + return + } else { + MSI_NXCOEFF(msi) = nxpix + MSI_NYCOEFF(msi) = nypix + MSI_FSTPNT(msi) = 0 + if (MSI_COEFF(msi) != NULL) + call mfree (MSI_COEFF(msi), TY_REAL) + call calloc (MSI_COEFF(msi), nxpix * nypix, TY_REAL) + } + + } + + # index the coefficient pointer so that COEFF(fptr+1) points to the + # first data point in the coefficient array + fptr = MSI_COEFF(msi) - 1 + MSI_FSTPNT(msi) + + # load data into coefficient array + rptr = fptr + do j = 1, nypix { + call amovr (datain[1,j], COEFF(rptr+1), nxpix) + rptr = rptr + MSI_NXCOEFF(msi) + } + + # calculate the coefficients of the interpolant + # boundary extension is performed using boundary projection + + switch (MSI_TYPE(msi)) { + + case II_BINEAREST, II_BISINC, II_BILSINC: + + # no end conditions necessary, coefficients stored as data + + case II_BILINEAR, II_BIDRIZZLE: + + # extend the rows + rptr = fptr + nxpix + do j = 1, nypix { + COEFF(rptr+1) = 2. * COEFF(rptr) - COEFF(rptr-1) + rptr = rptr + MSI_NXCOEFF(msi) + } + + # define the pointers to the last, 2nd last and third last rows + rptrl[1] = MSI_COEFF(msi) + (MSI_NYCOEFF(msi) - 1) * + MSI_NXCOEFF(msi) + do i = 2, 3 + rptrl[i] = rptrl[i-1] - MSI_NXCOEFF(msi) + + # define the last row by extending the columns + call awsur (COEFF(rptrl[2]), COEFF(rptrl[3]), COEFF(rptrl[1]), + MSI_NXCOEFF(msi), 2., -1.) + + case II_BIPOLY3: + + # extend the rows + rptr = fptr + nptr = fptr + nxpix + do j = 1, nypix { + COEFF(rptr) = 2. * COEFF(rptr+1) - COEFF(rptr+2) + COEFF(nptr+1) = 2. * COEFF(nptr) - COEFF(nptr-1) + COEFF(nptr+2) = 2. * COEFF(nptr) - COEFF(nptr-2) + rptr = rptr + MSI_NXCOEFF(msi) + nptr = nptr + MSI_NXCOEFF(msi) + } + + # define pointers to first, second and third rows + rptrf[1] = MSI_COEFF(msi) + do i = 2, 3 + rptrf[i] = rptrf[i-1] + MSI_NXCOEFF(msi) + + # extend the columns, first row + call awsur (COEFF(rptrf[2]), COEFF(rptrf[3]), COEFF(rptrf[1]), + MSI_NXCOEFF(msi), 2., -1.) + + # define the pointers to the last to fifth last rows + rptrl[1] = MSI_COEFF(msi) + (MSI_NYCOEFF(msi) - 1) * + MSI_NXCOEFF(msi) + do i = 2, 5 + rptrl[i] = rptrl[i-1] - MSI_NXCOEFF(msi) + + # extend the columns, define 2nd last row + call awsur (COEFF(rptrl[3]), COEFF(rptrl[4]), COEFF(rptrl[2]), + MSI_NXCOEFF(msi), 2., -1.) + + # extend the columns, define last row + call awsur (COEFF(rptrl[3]), COEFF(rptrl[5]), COEFF(rptrl[1]), + MSI_NXCOEFF(msi), 2., -1.) + + case II_BIPOLY5: + + # extend the rows + rptr = fptr + nptr = fptr + nxpix + do j = 1, nypix { + COEFF(rptr-1) = 2. * COEFF(rptr+1) - COEFF(rptr+3) + COEFF(rptr) = 2. * COEFF(rptr+1) - COEFF(rptr+2) + COEFF(nptr+1) = 2. * COEFF(nptr) - COEFF(nptr-1) + COEFF(nptr+2) = 2. * COEFF(nptr) - COEFF(nptr-2) + COEFF(nptr+3) = 2. * COEFF(nptr) - COEFF(nptr-3) + rptr = rptr + MSI_NXCOEFF(msi) + nptr = nptr + MSI_NXCOEFF(msi) + } + + # define pointers to first five rows + rptrf[1] = MSI_COEFF(msi) + do i = 2, 5 + rptrf[i] = rptrf[i-1] + MSI_NXCOEFF(msi) + + # extend the columns, define first row + call awsur (COEFF(rptrf[3]), COEFF(rptrf[5]), COEFF(rptrf[1]), + MSI_NXCOEFF(msi), 2., -1.) + + # extend the columns, define second row + call awsur (COEFF(rptrf[3]), COEFF(rptrf[4]), COEFF(rptrf[2]), + MSI_NXCOEFF(msi), 2., -1.) + + # define pointers last seven rows + rptrl[1] = MSI_COEFF(msi) + (MSI_NYCOEFF(msi) - 1) * + MSI_NXCOEFF(msi) + do i = 2, 7 + rptrl[i] = rptrl[i-1] - MSI_NXCOEFF(msi) + + # extend the columns, last row + call awsur (COEFF(rptrl[4]), COEFF(rptrl[7]), COEFF(rptrl[1]), + MSI_NXCOEFF(msi), 2., -1.) + + # extend the columns, 2nd last row + call awsur (COEFF(rptrl[4]), COEFF(rptrl[6]), COEFF(rptrl[2]), + MSI_NXCOEFF(msi), 2., -1.) + + # extend the columns, 3rd last row + call awsur (COEFF(rptrl[4]), COEFF(rptrl[5]), COEFF(rptrl[3]), + MSI_NXCOEFF(msi), 2., -1.) + + case II_BISPLINE3: + + # allocate space for a temporary work arrays + call calloc (tmp, MSI_NXCOEFF(msi) * MSI_NYCOEFF(msi), TY_REAL) + + # 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 + call ii_spline2d (COEFF(MSI_COEFF(msi)), TEMP(tmp), + nxpix, MSI_NYCOEFF(msi), MSI_NXCOEFF(msi), MSI_NYCOEFF(msi)) + + + # calculate the univariate B-spline coefficients in y to + # results of x interpolation + call ii_spline2d (TEMP(tmp), COEFF(MSI_COEFF(msi)), + nypix, MSI_NXCOEFF(msi), MSI_NYCOEFF(msi), MSI_NXCOEFF(msi)) + + # deallocate storage for temporary arrays + call mfree (tmp, TY_REAL) + } +end -- cgit