Initial commit

1 files changed, 275 insertions, 0 deletions

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 <math/iminterp.h>
+
+# 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