aboutsummaryrefslogtreecommitdiff
path: root/math/iminterp/mrider.x
diff options
context:
space:
mode:
Diffstat (limited to 'math/iminterp/mrider.x')
-rw-r--r--math/iminterp/mrider.x420
1 files changed, 420 insertions, 0 deletions
diff --git a/math/iminterp/mrider.x b/math/iminterp/mrider.x
new file mode 100644
index 00000000..8413df55
--- /dev/null
+++ b/math/iminterp/mrider.x
@@ -0,0 +1,420 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "im2interpdef.h"
+include <math/iminterp.h>
+
+# MRIDER -- Procedure to evaluate the derivatives of the interpolant
+# without the storage overhead required by the sequential version.
+# The derivatives are stored such that der[1,1] = the value of the
+# interpolant at x and y, der[2,1] = the first derivative in x and
+# der[2,1] = the first derivative in y.
+
+procedure mrider (x, y, datain, nxpix, nypix, len_datain, der, nxder, nyder,
+ len_der, interp_type)
+
+real x[ARB] # x value
+real y[ARB] # y value
+real datain[len_datain,ARB] # data array
+int nxpix # number of x data points
+int nypix # number of y data points
+int len_datain # row length of datain
+real der[len_der, ARB] # array of derivatives
+int nxder # number of derivatives in x
+int nyder # number of derivatives in y
+int len_der # row length of der, len_der >= nxder
+int interp_type # interpolant type
+
+int nx, ny, nxterms, nyterms, row_length
+int index, xindex, yindex, first_row, last_row
+int i, j, ii, jj, kx, ky
+pointer tmp
+real coeff[SPLPTS+3,SPLPTS+3], pcoeff[MAX_NTERMS,MAX_NTERMS]
+real pctemp[MAX_NTERMS,MAX_NTERMS], sum[MAX_NTERMS]
+real hold21, hold12, hold22, accum, deltax, deltay, tmpx[4], tmpy[4]
+real xmin, xmax, ymin, ymax, sx, sy, tx, ty
+
+errchk malloc, calloc, mfree
+
+begin
+ if (nxder < 1 || nyder < 1)
+ return
+
+ # zero the derivatives
+ do j = 1, nyder {
+ do i = 1, nxder
+ der[i,j] = 0.
+ }
+
+ switch (interp_type) {
+
+ case II_BINEAREST:
+
+ der[1,1] = datain[int (x[1]+.5), int (y[1]+.5)]
+
+ return
+
+ case II_BISINC, II_BILSINC:
+
+ call ii_bisincder (x[1], y[1], der, nxder, nyder, len_der, datain,
+ 0, len_datain, nypix, NSINC, DX, DY)
+
+ return
+
+ case II_BILINEAR:
+
+ nx = x[1]
+ sx = x[1] - nx
+ tx = 1. - sx
+
+ ny = y[1]
+ sy = y[1] - ny
+ ty = 1. - sy
+
+ # protect against the case where x = nxpix and/or y = nypix
+ if (nx >= nxpix)
+ hold21 = 2. * datain[nx,ny] - datain[nx-1,ny]
+ else
+ hold21 = datain[nx+1,ny]
+ if (ny >= nypix)
+ hold12 = 2. * datain[nx,ny] - datain[nx,ny-1]
+ else
+ hold12 = datain[nx,ny+1]
+ if (nx >= nxpix && ny >= nypix)
+ hold22 = 2. * hold21 - (2. * datain[nx,ny-1] -
+ datain[nx-1,ny-1])
+ else if (nx >= nxpix)
+ hold22 = 2. * hold12 - datain[nx-1,ny+1]
+ else if (ny >= nypix)
+ hold 22 = 2. * hold21 - datain[nx+1,ny-1]
+ else
+ hold22 = datain[nx+1,ny+1]
+
+ # evaluate the derivatives
+ der[1,1] = tx * ty * datain[nx,ny] + sx * ty * hold21 +
+ sy * tx * hold12 + sx * sy * hold22
+ if (nxder > 1)
+ der[2,1] = - ty * datain[nx,ny] + ty * hold21 -
+ sy * hold12 + sy * hold22
+ if (nyder > 1)
+ der[1,2] = - tx * datain[nx,ny] - sx * hold21 +
+ tx * hold12 + sx * hold22
+ if (nxder > 1 && nyder > 1)
+ der[2,2] = datain[nx,ny] - hold21 - hold12 + hold22
+
+
+ return
+
+ case II_BIDRIZZLE:
+ call ii_bidriz1 (datain, 0, len_datain, x, y, der[1,1], 1, BADVAL)
+ if (nxder > 1) {
+ xmax = max (x[1], x[2], x[3], x[4])
+ xmin = min (x[1], x[2], x[3], x[4])
+ ymax = max (y[1], y[2], y[3], y[4])
+ ymin = min (y[1], y[2], y[3], y[4])
+ deltax = xmax - xmin
+ if (deltax == 0.0)
+ der[2,1] = 0.0
+ else {
+ tmpx[1] = xmin; tmpy[1] = ymin
+ tmpx[2] = (xmax - xmin) / 2.0; tmpy[2] = ymin
+ tmpx[3] = (xmax - xmin) / 2.0; tmpy[3] = ymax
+ tmpx[4] = xmin; tmpy[4] = ymax
+ call ii_bidriz1 (datain, 0, len_datain, tmpx, tmpy,
+ accum, 1, BADVAL)
+ tmpx[1] = (xmax - xmin) / 2.0; tmpy[1] = ymin
+ tmpx[2] = xmax; tmpy[2] = ymin
+ tmpx[3] = xmax; tmpy[3] = ymax
+ tmpx[4] = (xmax - xmin) / 2.0; tmpy[4] = ymax
+ call ii_bidriz1 (datain, 0, len_datain, tmpx, tmpy,
+ der[2,1], 1, BADVAL)
+ der[2,1] = 2.0 * (der[2,1] - accum) / deltax
+ }
+ }
+ if (nyder > 1) {
+ deltay = ymax - ymin
+ if (deltay == 0.0)
+ der[1,2] = 0.0
+ else {
+ tmpx[1] = xmin; tmpy[1] = ymin
+ tmpx[2] = xmax; tmpy[2] = ymin
+ tmpx[3] = xmax; tmpy[3] = (ymax - ymin) / 2.0
+ tmpx[4] = xmin; tmpy[4] = (ymax - ymin) / 2.0
+ call ii_bidriz1 (datain, 0, len_datain, tmpx, tmpy,
+ accum, 1, BADVAL)
+ tmpx[1] = xmin; tmpy[1] = (ymax - ymin) / 2.0
+ tmpx[2] = xmax; tmpy[2] = (ymax - ymin) / 2.0
+ tmpx[3] = xmax; tmpy[3] = ymax
+ tmpx[4] = xmin; tmpy[4] = ymax
+ call ii_bidriz1 (datain, 0, len_datain, tmpx, tmpy,
+ der[1,2], 1, BADVAL)
+ der[1,2] = 2.0 * (der[1,2] - accum) / deltay
+ }
+ }
+
+ return
+
+ case II_BIPOLY3:
+
+ row_length = SPLPTS + 3
+
+ nxterms = 4
+ nyterms = 4
+
+ nx = x[1]
+ ny = y[1]
+
+ sx = x[1] - nx
+ sy = y[1] - ny
+
+ # use boundary projection to extend the data rows
+ yindex = 1
+ for (j = ny - 1; j <= ny + 2; j = j + 1) {
+
+ # check that the data row is defined
+ if (j >= 1 && j <= nypix) {
+
+ # extend the rows
+ xindex = 1
+ for (i = nx - 1; i <= nx + 2; i = i + 1) {
+ if (i < 1)
+ coeff[xindex,yindex] = 2. * datain[1,j] -
+ datain[2-i,j]
+ else if (i > nxpix)
+ coeff[xindex,yindex] = 2. * datain[nxpix,j] -
+ datain[2*nxpix-i,j]
+ else
+ coeff[xindex,yindex] = datain[i,j]
+ xindex = xindex + 1
+ }
+ } else if (j == (nypix + 2)) {
+
+ # allow for the final row
+ xindex = 1
+ for (i = nx - 1; i <= nx + 2; i = i + 1) {
+ if (i < 1)
+ coeff[xindex,nyterms] = 2. * datain[1,nypix-2] -
+ datain[2-i,nypix-2]
+ else if (i > nxpix)
+ coeff[xindex,nyterms] = 2. * datain[nxpix,nypix-2] -
+ datain[2*nxpix-i,nypix-2]
+ else
+ coeff[xindex,nyterms] = datain[i,nypix-2]
+ xindex = xindex + 1
+ }
+
+ }
+
+ yindex = yindex + 1
+ }
+
+
+ # project columns
+ first_row = max (1, 3 - ny)
+ if (first_row > 1) {
+ for (j = 1; j < first_row; j = j + 1)
+ call awsur (coeff[1, first_row], coeff[1, 2*first_row-j],
+ coeff[1,j], nxterms, 2., -1.)
+ }
+
+ last_row = min (nxterms, nypix - ny + 2)
+ if (last_row < nxterms) {
+ for (j = last_row + 1; j <= nxterms - 1; j = j + 1)
+ call awsur (coeff[1,last_row], coeff[1,2*last_row-j],
+ coeff[1,j], nxterms, 2., -1.)
+ if (last_row == 2)
+ call awsur (coeff[1,last_row], coeff[1,4], coeff[1,4],
+ nxterms, 2., -1.)
+ else
+ call awsur (coeff[1,last_row], coeff[1,2*last_row-4],
+ coeff[1,4], nxterms, 2., -1.)
+ }
+
+ # calculate the coefficients of the bicubic polynomial
+ call ii_pcpoly3 (coeff, 2, row_length, pcoeff, MAX_NTERMS)
+
+ case II_BIPOLY5:
+ row_length = SPLPTS + 3
+
+ nxterms = 6
+ nyterms = 6
+
+ nx = x[1]
+ ny = y[1]
+
+ sx = x[1] - nx
+ sy = y[1] - ny
+
+ # extend rows of data
+ yindex = 1
+ for (j = ny - 2; j <= ny + 3; j = j + 1) {
+
+ # select the rows containing data
+ if (j >= 1 && j <= nypix) {
+
+ # extend the rows
+ xindex = 1
+ for (i = nx - 2; i <= nx + 3; i = i + 1) {
+ if (i < 1)
+ coeff[xindex,yindex] = 2. * datain[1,j] -
+ datain[2-i,j]
+ else if (i > nxpix)
+ coeff[xindex,yindex] = 2. * datain[nxpix,j] -
+ datain[2*nxpix-i,j]
+ else
+ coeff[xindex,yindex] = datain[i,j]
+ xindex = xindex + 1
+ }
+
+ } else if (j == (ny + 3)) {
+
+ # extend the rows
+ xindex = 1
+ for (i = nx - 2; i <= nx + 3; i = i + 1) {
+ if (i < 1)
+ coeff[xindex,yindex] = 2. * datain[1,nypix-3] -
+ datain[2-i,nypix-3]
+ else if (i > nxpix)
+ coeff[xindex,yindex] = 2. * datain[nxpix,nypix-3] -
+ datain[2*nxpix-i,nypix-3]
+ else
+ coeff[xindex,yindex] = datain[i,nypix-3]
+ xindex = xindex + 1
+ }
+
+ }
+
+ yindex = yindex + 1
+ }
+
+ # project columns
+
+ first_row = max (1, 4 - ny)
+ if (first_row > 1) {
+ for (j = 1; j < first_row; j = j + 1)
+ call awsur (coeff[1,first_row], coeff[1,2*first_row-j],
+ coeff[1,j], nxterms, 2., -1.)
+ }
+
+ last_row = min (nxterms, nypix - ny + 3)
+ if (last_row < nxterms) {
+ for (j = last_row + 1; j <= nxterms - 1; j = j + 1)
+ call awsur (coeff[1,last_row], coeff[1,2*last_row-j],
+ coeff[1,j], nxterms, 2., -1.)
+ if (last_row == 3)
+ call awsur (coeff[1,last_row], coeff[1,6], coeff[1,6],
+ nxterms, 2., -1.)
+ else
+ call awsur (coeff[1,last_row], coeff[1,2*last_row-6],
+ coeff[1,6], nxterms, 2., -1.)
+ }
+
+ # caculate the polynomial coeffcients
+ call ii_pcpoly5 (coeff, 3, row_length, pcoeff, MAX_NTERMS)
+
+
+ case II_BISPLINE3:
+ row_length = SPLPTS + 3
+
+ nxterms = 4
+ nyterms = 4
+
+ nx = x[1]
+ ny = y[1]
+
+ sx = x[1] - nx
+ sy = y[1] - ny
+
+ # allocate space for temporary array and 0 file
+ call calloc (tmp, row_length * row_length, TY_REAL)
+
+ ky = 0
+ # maximum number of points used in each direction is SPLPTS
+ for (j = ny - SPLPTS/2 + 1; j <= ny + SPLPTS/2; j = j + 1) {
+
+ if (j < 1 || j > nypix)
+ ;
+ else {
+ ky = ky + 1
+ if (ky == 1)
+ yindex = ny - j + 1
+
+ kx = 0
+ for (i = nx - SPLPTS/2 + 1; i <= nx + SPLPTS/2; i = i + 1) {
+ if (i < 1 || i > nxpix)
+ ;
+ else {
+ kx = kx + 1
+ if (kx == 1)
+ xindex = nx - i + 1
+ coeff[kx+1,ky+1] = datain[i,j]
+ }
+ }
+
+ coeff[1,ky+1] = 0.
+ coeff[kx+2,ky+1] = 0.
+ coeff[kx+3,ky+1] = 0.
+
+ }
+ }
+
+ # zero out 1st and last 2 rows
+ call amovkr (0., coeff[1,1], kx+3)
+ call amovkr (0., coeff[1,ky+2], kx+3)
+ call amovkr (0., coeff[1,ky+3],kx+3)
+
+ # calculate the spline coefficients
+ call ii_spline2d (coeff, Memr[tmp], kx, ky+2, row_length,
+ row_length)
+ call ii_spline2d (Memr[tmp], coeff, ky, kx+2, row_length,
+ row_length)
+
+ # calculate the polynomial coefficients
+ index = (yindex - 1) * row_length + xindex + 1
+ call ii_pcspline3 (coeff, index, row_length, pcoeff, MAX_NTERMS)
+
+ # free space
+ call mfree (tmp, TY_REAL)
+ }
+
+ # evaluate the derivatives of the higher order interpolants
+ do j = 1, nyder {
+
+ # set pctemp
+ do jj = nyterms, j, -1 {
+ do ii = 1, nxterms
+ pctemp[ii,jj] = pcoeff[ii,jj]
+ }
+
+ do i = 1, nxder {
+
+ # accumulate the partial sums in x
+ do jj = nyterms, j, -1 {
+ sum[jj] = pctemp[nxterms,jj]
+ do ii = nxterms - 1, i, -1
+ sum[jj] = pctemp[ii,jj] + sum[jj] * sx
+ }
+
+ # accumulate the sum in y
+ accum = sum[nyterms]
+ do jj = nyterms - 1, j, -1
+ accum = sum[jj] + accum * sy
+
+ # evaulate the derivative
+ der[i,j] = accum
+
+ # differentiate in x
+ do jj = nyterms, j, -1 {
+ do ii = nxterms, i + 1, -1
+ pctemp[ii,jj] = (ii - i) * pctemp[ii,jj]
+ }
+
+ }
+
+ # differentiate in y
+ do jj = 1, nxterms {
+ do ii = nyterms, j + 1, -1
+ pcoeff[jj,ii] = (ii - j) * pcoeff[jj,ii]
+ }
+
+ }
+end
generated by cgit (git 2.34.1) at 2025-09-20 10:15:45 -0400