aboutsummaryrefslogtreecommitdiff
path: root/math/iminterp/mrieval.x
diff options
context:
space:
mode:
Diffstat (limited to 'math/iminterp/mrieval.x')
-rw-r--r--math/iminterp/mrieval.x303
1 files changed, 303 insertions, 0 deletions
diff --git a/math/iminterp/mrieval.x b/math/iminterp/mrieval.x
new file mode 100644
index 00000000..6d89c456
--- /dev/null
+++ b/math/iminterp/mrieval.x
@@ -0,0 +1,303 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "im2interpdef.h"
+include <math/iminterp.h>
+
+# MRIEVAL -- Procedure to evaluate the 2D interpolant at a given value
+# of x and y. MRIEVAL allows the interpolation of a few interpolated
+# points without the computing time and storage required for the
+# sequential version. The routine assumes that 1 <= x <= nxpix and
+# 1 <= y <= nypix.
+
+real procedure mrieval (x, y, datain, nxpix, nypix, len_datain, 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
+int interp_type # interpolant type
+
+int nx, ny, nterms, row_length
+int xindex, yindex, first_row, last_row
+int kx, ky
+int i, j
+pointer tmp
+real coeff[SPLPTS+3,SPLPTS+3]
+real hold21, hold12, hold22
+real sx, sy, tx, ty
+real xval, yval, value
+errchk malloc, calloc, mfree
+
+begin
+ switch (interp_type) {
+
+ case II_BINEAREST:
+ return (datain[int (x[1]+0.5), int (y[1]+0.5)])
+
+ 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)
+ hold22 = 2. * hold21 - datain[nx+1,ny-1]
+ else
+ hold22 = datain[nx+1,ny+1]
+
+ # evaluate the interpolant
+ value = tx * ty * datain[nx,ny] + sx * ty * hold21 +
+ sy * tx * hold12 + sx * sy * hold22
+
+ return (value)
+
+ case II_BIDRIZZLE:
+ call ii_bidriz1 (datain, 0, len_datain, x, y, value, 1, BADVAL)
+
+ return (value)
+
+ case II_BIPOLY3:
+ row_length = SPLPTS + 3
+ nterms = 4
+ nx = x[1]
+ ny = y[1]
+
+ # major problem is that near the edge the interior polynomial
+ # must be defined
+
+ # 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 == (ny + 2)) {
+
+ # extend the rows
+ xindex = 1
+ for (i = nx - 1; i <= nx + 2; i = i + 1) {
+ if (i < 1)
+ coeff[xindex,yindex] = 2. * datain[1,nypix-2] -
+ datain[2-i,nypix-2]
+ else if (i > nxpix)
+ coeff[xindex,yindex] = 2. * datain[nxpix,nypix-2] -
+ datain[2*nxpix-i,nypix-2]
+ else
+ coeff[xindex,yindex] = 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], nterms, 2., -1.)
+ }
+
+ last_row = min (nterms, nypix - ny + 2)
+ if (last_row < nterms) {
+ for (j = last_row + 1; j <= nterms - 1; j = j + 1)
+ call awsur (coeff[1,last_row], coeff[1,2*last_row-j],
+ coeff[1,j], nterms, 2., -1.)
+ if (last_row == 2)
+ call awsur (coeff[1,last_row], coeff[1,4], coeff[1,4],
+ nterms, 2., -1.)
+ else
+ call awsur (coeff[1,last_row], coeff[1,2*last_row-4],
+ coeff[1,4], nterms, 2., -1.)
+ }
+
+
+ # center the x value and call evaluation routine
+ xval = 2 + (x[1] - nx)
+ yval = 2 + (y[1] - ny)
+ call ii_bipoly3 (coeff, 0, row_length, xval, yval, value, 1)
+
+ return (value)
+
+ case II_BIPOLY5:
+ row_length = SPLPTS + 3
+ nterms = 6
+ nx = x[1]
+ ny = y[1]
+
+ # major problem is to define interior polynomial near the edge
+
+ # loop over the 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], nterms, 2., -1.)
+ }
+
+ last_row = min (nterms, nypix - ny + 3)
+ if (last_row < nterms) {
+ for (j = last_row + 1; j <= nterms - 1; j = j + 1)
+ call awsur (coeff[1,last_row], coeff[1,2*last_row-j],
+ coeff[1,j], nterms, 2., -1.)
+ if (last_row == 3)
+ call awsur (coeff[1,last_row], coeff[1,6], coeff[1,6],
+ nterms, 2., -1.)
+ else
+ call awsur (coeff[1,last_row], coeff[1,2*last_row-6],
+ coeff[1,6], nterms, 2., -1.)
+ }
+
+ # call evaluation routine
+ xval = 3 + (x[1] - nx)
+ yval = 3 + (y[1] - ny)
+ call ii_bipoly5 (coeff, 0, row_length, xval, yval, value, 1)
+
+ return (value)
+
+ case II_BISPLINE3:
+ row_length = SPLPTS + 3
+ nx = x[1]
+ ny = y[1]
+
+ # 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)
+
+ # evaluate spline
+ xval = xindex + 1 + (x[1] - nx)
+ yval = yindex + 1 + (y[1] - ny)
+ call ii_bispline3 (coeff, 0, row_length, xval, yval, value, 1)
+
+ # free space
+ call mfree (tmp, TY_REAL)
+
+ return (value)
+
+ case II_BISINC, II_BILSINC:
+ call ii_bisinc (datain, 0, len_datain, nypix, x, y, value, 1,
+ NSINC, DX, DY)
+
+ return (value)
+ }
+end
generated by cgit (git 2.34.1) at 2025-09-20 07:11:21 -0400