aboutsummaryrefslogtreecommitdiff
path: root/math/iminterp/msigrl.x
diff options
context:
space:
mode:
Diffstat (limited to 'math/iminterp/msigrl.x')
-rw-r--r--math/iminterp/msigrl.x238
1 files changed, 238 insertions, 0 deletions
diff --git a/math/iminterp/msigrl.x b/math/iminterp/msigrl.x
new file mode 100644
index 00000000..7baeac34
--- /dev/null
+++ b/math/iminterp/msigrl.x
@@ -0,0 +1,238 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <mach.h>
+include "im2interpdef.h"
+include <math/iminterp.h>
+
+# MSIGRL -- Procedure to integrate the 2D interpolant over a specified area.
+# The x and y arrays are assumed to describe a polygon which is the domain over
+# which the integration is to be performed. The x and y must describe a closed
+# curve and npts must be >= 4 with the last vertex equal to the first vertex.
+# The routine uses the technique of separation of variables. The restriction on
+# the polygon is that horizontal lines have at most one segment in common with
+# the domain of integration. Polygons which do not fit this restriction can be
+# split into one or more polygons before calling msigrl and the results can
+# then be summed.
+
+real procedure msigrl (msi, x, y, npts)
+
+pointer msi # pointer to the interpolant descriptor structure
+real x[npts] # array of x values
+real y[npts] # array of y values
+int npts # number of points which describe the boundary
+
+int i, interp_type, nylmin, nylmax, offset
+pointer x1lim, x2lim, xintegrl, ptr
+real xmin, xmax, ymin, ymax, accum
+real ii_1dinteg()
+
+begin
+ # set up 1D interpolant type
+ switch (MSI_TYPE(msi)) {
+ case II_BINEAREST:
+ interp_type = II_NEAREST
+ case II_BILINEAR:
+ interp_type = II_LINEAR
+ case II_BIDRIZZLE:
+ interp_type = II_DRIZZLE
+ case II_BIPOLY3:
+ interp_type = II_POLY3
+ case II_BIPOLY5:
+ interp_type = II_POLY5
+ case II_BISPLINE3:
+ interp_type = II_SPLINE3
+ case II_BISINC:
+ interp_type = II_SINC
+ case II_BILSINC:
+ interp_type = II_LSINC
+ }
+
+ # set up temporary storage for x limits and the x integrals
+ call calloc (x1lim, MSI_NYCOEFF(msi), TY_REAL)
+ call calloc (x2lim, MSI_NYCOEFF(msi), TY_REAL)
+ call calloc (xintegrl, MSI_NYCOEFF(msi), TY_REAL)
+
+ # offset of first data point from edge of coefficient array
+ offset = mod (MSI_FSTPNT(msi), MSI_NXCOEFF(msi))
+
+ # convert the (x,y) points which describe the polygon into
+ # two arrays of x limits x1lim and x2lim and two y limits ymin and ymax
+ call ii_find_limits (x, y, npts, 0, 0, MSI_NYCOEFF(msi),
+ Memr[x1lim+offset], Memr[x2lim+offset], ymin, ymax, nylmin, nylmax)
+ nylmin = nylmin + offset
+ nylmax = nylmax + offset
+
+ # integrate in x
+ ptr = MSI_COEFF(msi) + offset + (nylmin - 1) * MSI_NXCOEFF(msi)
+ do i = nylmin, nylmax {
+ xmin = min (Memr[x1lim+i-1], Memr[x2lim+i-1])
+ xmax = max (Memr[x1lim+i-1], Memr[x2lim+i-1])
+ Memr[xintegrl+i-1] = ii_1dinteg (COEFF(ptr), MSI_NXCOEFF(msi),
+ xmin, xmax, interp_type, MSI_NSINC(msi), DX, MSI_XPIXFRAC(msi))
+ ptr = ptr + MSI_NXCOEFF(msi)
+ }
+
+ # integrate in y
+ if (interp_type == II_SPLINE3) {
+ call amulkr (Memr[xintegrl], 6.0, Memr[xintegrl], MSI_NYCOEFF(msi))
+ accum = ii_1dinteg (Memr[xintegrl+offset], MSI_NYCOEFF(msi), ymin,
+ ymax, II_NEAREST, MSI_NSINC(msi), DY, MSI_YPIXFRAC(msi))
+ } else {
+ accum = ii_1dinteg (Memr[xintegrl+offset], MSI_NYCOEFF(msi), ymin,
+ ymax, II_NEAREST, MSI_NSINC(msi), DY, MSI_YPIXFRAC(msi))
+ }
+
+ # free space
+ call mfree (xintegrl, TY_REAL)
+ call mfree (x1lim, TY_REAL)
+ call mfree (x2lim, TY_REAL)
+
+ return (accum)
+end
+
+
+# II_FIND_LIMITS -- Procedure to transform a set of (x,y)'s describing a
+# polygon into a set of limits.
+
+procedure ii_find_limits (x, y, npts, xboff, xeoff, max_nylines, x1lim, x2lim,
+ymin, ymax, nylmin, nylmax)
+
+real x[npts] # array of x values
+real y[npts] # array of y values
+int npts # number of data points
+int xboff, xeoff # boundary extension limits
+int max_nylines # max number of lines to integrate
+real x1lim[ARB] # array of x1 limits
+real x2lim[ARB] # array of x2 limits
+real ymin # minimum y value for integration
+real ymax # maximum y value for integration
+int nylmin # minimum line number for x integration
+int nylmax # maximum line number for x integration
+
+int i, ninter
+pointer sp, xintr, yintr
+real xmin, xmax, lx, ld
+int ii_pyclip()
+
+begin
+ call smark (sp)
+ call salloc (xintr, npts, TY_REAL)
+ call salloc (yintr, npts, TY_REAL)
+
+ # find x and y limits and their indicess
+ call alimr (x, npts, xmin, xmax)
+ call alimr (y, npts, ymin, ymax)
+
+ # calculate the line limits for integration
+ nylmin = max (1, min (int (ymin + 0.5) - xboff, max_nylines))
+ nylmax = min (max_nylines, max (1, int (ymax + 0.5) + xeoff))
+
+ # initialize
+ lx = xmax - xmin
+
+ # calculate the limits
+ for (i = nylmin; i <= nylmax; i = i + 1) {
+
+ if (ymin > i)
+ ld = min (i + 0.5, ymax) * lx
+ else if (ymax < i)
+ ld = max (i - 0.5, ymin) * lx
+ else
+ ld = i * lx
+ ninter = ii_pyclip (x, y, Memr[xintr], Memr[yintr], npts, lx, ld)
+ if (ninter <= 0) {
+ x1lim[i] = xmin
+ x2lim[i] = xmin
+ } else {
+ x1lim[i] = min (Memr[xintr], Memr[xintr+1])
+ x2lim[i] = max (Memr[xintr], Memr[xintr+1])
+ }
+ }
+
+ call sfree (sp)
+end
+
+
+# II_YCLIP -- Procedure to determine the intersection points of a
+# horizontal image line with an arbitrary polygon.
+
+int procedure ii_pyclip (xver, yver, xintr, yintr, nver, lx, ld)
+
+real xver[ARB] # x vertex coords
+real yver[ARB] # y vertex coords
+real xintr[ARB] # x intersection coords
+real yintr[ARB] # y intersection coords
+int nver # number of vertices
+real lx, ld # equation of image line
+
+int i, nintr
+real u1, u2, u1u2, dx, dy, dd, xa, ya, wa
+
+begin
+ nintr = 0
+ u1 = - lx * yver[1] + ld
+ do i = 2, nver {
+
+ u2 = - lx * yver[i] + ld
+ u1u2 = u1 * u2
+
+ # Test whether polygon line segment intersects image line or not.
+ if (u1u2 <= 0.0) {
+
+
+ # Compute the intersection coords.
+ if (u1 != 0.0 && u2 != 0.0) {
+
+ dy = yver[i-1] - yver[i]
+ dx = xver[i-1] - xver[i]
+ dd = xver[i-1] * yver[i] - yver[i-1] * xver[i]
+ xa = (dx * ld - lx * dd)
+ ya = dy * ld
+ wa = dy * lx
+ nintr = nintr + 1
+ xintr[nintr] = xa / wa
+ yintr[nintr] = ya / wa
+
+ # Test for collinearity.
+ } else if (u1 == 0.0 && u2 == 0.0) {
+
+ nintr = nintr + 1
+ xintr[nintr] = xver[i-1]
+ yintr[nintr] = yver[i-1]
+ nintr = nintr + 1
+ xintr[nintr] = xver[i]
+ yintr[nintr] = yver[i]
+
+ } else if (u1 != 0.0) {
+
+ if (i == 1) {
+ dy = (yver[2] - yver[1])
+ dd = (yver[nver-1] - yver[1])
+ } else if (i == nver) {
+ dy = (yver[2] - yver[nver])
+ dd = dy * (yver[nver-1] - yver[nver])
+ } else {
+ dy = (yver[i+1] - yver[i])
+ dd = dy * (yver[i-1] - yver[i])
+ }
+
+ if (dy != 0.0) {
+ nintr = nintr + 1
+ xintr[nintr] = xver[i]
+ yintr[nintr] = yver[i]
+ }
+
+ if (dd > 0.0) {
+ nintr = nintr + 1
+ xintr[nintr] = xver[i]
+ yintr[nintr] = yver[i]
+ }
+
+ }
+ }
+
+ u1 = u2
+ }
+
+ return (nintr)
+end
generated by cgit (git 2.34.1) at 2025-09-20 07:10:49 -0400