aboutsummaryrefslogtreecommitdiff
path: root/pkg/proto/maskexpr/peregufcn.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/proto/maskexpr/peregufcn.x')
-rw-r--r--pkg/proto/maskexpr/peregufcn.x808
1 files changed, 808 insertions, 0 deletions
diff --git a/pkg/proto/maskexpr/peregufcn.x b/pkg/proto/maskexpr/peregufcn.x
new file mode 100644
index 00000000..8d4a64e1
--- /dev/null
+++ b/pkg/proto/maskexpr/peregufcn.x
@@ -0,0 +1,808 @@
+include <math.h>
+include <plset.h>
+include <plio.h>
+include "peregfuncs.h"
+
+
+# PE_UCIRCLE -- Regionrop ufcn for a circle (circular region), clipped at
+# the borders of the mask.
+
+bool procedure pe_ucircle (ufd, y, rl_reg, xs, npix)
+
+pointer ufd #I user function descriptor
+int y #I mask line number
+int rl_reg[3,ARB] #O output range list for line Y
+int xs #O first pixel to be edited
+int npix #O number of pixels affected
+
+real radius, dx, dy
+pointer pl
+int rn, axlen, x1, x1_clipped, x2, x2_clipped
+
+begin
+ pl = C_PL(ufd)
+ rn = RL_FIRST
+ axlen = PL_AXLEN(pl,1)
+ radius = C_RADIUS(ufd)
+
+ dy = abs (C_YCEN(ufd) - y)
+ if (dy < radius) {
+ dx = radius * radius - dy * dy
+ if (dx > 0.0)
+ dx = sqrt (dx)
+ else
+ dx = 0.0
+ x1 = int(C_XCEN(ufd) - dx)
+ x2 = int(C_XCEN(ufd) + dx)
+ x1_clipped = max(1, min (axlen, x1))
+ x2_clipped = max(x1, min (axlen, x2))
+ xs = x1_clipped
+ npix = x2_clipped - x1_clipped + 1
+ RL_X(rl_reg,rn) = 1
+ RL_N(rl_reg,rn) = npix
+ RL_V(rl_reg,rn) = C_PV(ufd)
+ rn = rn + 1
+ } else {
+ npix = 0
+ xs = 1
+ }
+
+ RL_LEN(rl_reg) = rn - 1
+ RL_AXLEN(rl_reg) = npix
+
+ return (true)
+end
+
+
+# PE_UELLIPSE -- Regionrop ufcn for an ellipse (elliptical region), clipped at
+# the borders of the mask.
+
+bool procedure pe_uellipse (ufd, y, rl_reg, xs, npix)
+
+pointer ufd #I user function descriptor
+int y #I mask line number
+int rl_reg[3,ARB] #O output range list for line Y
+int xs #O first pixel to be edited
+int npix #O number of pixels affected
+
+real dy, dy2, ady, bb, cc, discr, dx1, dx2
+pointer pl
+int rn, axlen, x1, x1_clipped, x2, x2_clipped
+
+begin
+ pl = E_PL(ufd)
+ rn = RL_FIRST
+ axlen = PL_AXLEN(pl,1)
+
+ dy = y - E_YCEN(ufd)
+ dy2 = dy * dy
+ ady = abs (dy)
+ bb = E_BB(ufd) * dy
+ cc = E_CC(ufd) * dy2
+ if (ady < E_DYMAX(ufd)) {
+ discr = bb * bb - 4.0 * E_AA(ufd) * (cc - E_FF(ufd))
+ if (discr > 0.0)
+ discr = sqrt (discr)
+ else
+ discr = 0.0
+ dx1 = (-bb - discr) / 2.0 / E_AA(ufd)
+ dx2 = (-bb + discr) / 2.0 / E_AA(ufd)
+ x1 = int(E_XCEN(ufd) + min (dx1, dx2))
+ x2 = int(E_XCEN(ufd) + max (dx1, dx2))
+ x1_clipped = max(1, min (axlen, x1))
+ x2_clipped = max(x1, min (axlen, x2))
+ xs = x1_clipped
+ npix = x2_clipped - x1_clipped + 1
+ RL_X(rl_reg,rn) = 1
+ RL_N(rl_reg,rn) = npix
+ RL_V(rl_reg,rn) = E_PV(ufd)
+ rn = rn + 1
+ } else {
+ npix = 0
+ xs = 1
+ }
+
+ RL_LEN(rl_reg) = rn - 1
+ RL_AXLEN(rl_reg) = npix
+
+ return (true)
+end
+
+
+# PE_UBOX -- Regionrop ufcn for an unrotated rectangle (rectangular region),
+# clipped at the borders of the mask.
+
+bool procedure pe_ubox (ufd, y, rl_reg, xs, npix)
+
+pointer ufd #I user function descriptor
+int y #I mask line number
+int rl_reg[3,ARB] #O output range list for line Y
+int xs #O first pixel to be edited
+int npix #O number of pixels affected
+
+int rn
+bool rl_new
+
+begin
+ rl_new = true
+ rn = RL_FIRST
+
+ if (y >= B_Y1(ufd) && y <= B_Y2(ufd)) {
+ xs = B_X1(ufd)
+ npix = B_X2(ufd) - B_X1(ufd) + 1
+ RL_X(rl_reg,rn) = 1
+ RL_N(rl_reg,rn) = npix
+ RL_V(rl_reg,rn) = B_PV(ufd)
+ rn = rn + 1
+ } else {
+ npix = 0
+ xs = 1
+ }
+
+ RL_LEN(rl_reg) = rn - 1
+ RL_AXLEN(rl_reg) = npix
+
+ return (true)
+end
+
+
+# PE_UPOLYGON -- Regionrop ufcn for a general closed polygonal region.
+# This a copy of pl_upolgon which contains the following editorial comment.
+# Surely there must be a simpler way to code this ... I have a polygon
+# routines of my own which I use in the photometry code which may be
+# a bit simpler. Might replace this at some point.
+
+bool procedure pe_upolygon (ufd, line, rl_reg, xs, npix)
+
+pointer ufd #I user function descriptor
+int line #I mask line number
+int rl_reg[3,ARB] #O output range list for line Y
+int xs #O start of edit region in dst mask
+int npix #O number of pixels affected
+
+pointer xp, yp, pl
+bool rl_new, cross
+int nseg, np, low, rn, i1, i2, ii, i, j
+int tempi, axlen, rl_len, p_prev, p_next
+real tempr, y, y1, y2, x1, x2, p1, p2, p_y, n_y
+
+int btoi()
+bool plr_equali()
+define done_ 91
+
+begin
+ pl = P_PL(ufd)
+ axlen = PL_AXLEN(pl,1)
+ rn = RL_FIRST
+ npix = 0
+ xs = 1
+
+ nseg = P_NS(ufd)
+ xp = P_XP(ufd)
+ yp = P_YP(ufd)
+ y = real(line)
+
+ # Find the point(s) of intersection of the current mask line with
+ # the line segments forming the polygon. Note that the line must
+ # cross a segment to go from inside to outside or vice versa; if a
+ # segment (or vertex) is merely touched it should be drawn, but it
+ # is not a point of crossing.
+
+ do i = 1, nseg {
+ # Locate next and previous line segments.
+ if (i == 1)
+ p_prev = nseg
+ else
+ p_prev = i - 1
+ if (i == nseg)
+ p_next = 1
+ else
+ p_next = i + 2
+
+ # Get endpoints of current segment.
+ x1 = Memr[xp+i-1]; x2 = Memr[xp+i]
+ y1 = Memr[yp+i-1]; y2 = Memr[yp+i]
+ if (y1 > y2) {
+ swapr (x1, x2)
+ swapr (y1, y2)
+ swapi (p_next, p_prev)
+ }
+
+ # Does current line intersect the polygon line segment?
+ if (y > y1-TOL && y < y2+TOL) {
+ p_y = Memr[yp+p_prev-1]
+ n_y = Memr[yp+p_next-1]
+
+ if (y2 - y1 > TOL) {
+ # Single point of intersection.
+ p1 = x1 + ((x2 - x1) / (y2 - y1)) * (y - y1)
+ p2 = p1
+
+ if (equal (p1, x1) && equal (y, y1))
+ cross = ((p_y - y1) < 0)
+ else if (equal (p1, x2) && equal (y, y2))
+ cross = ((n_y - y2) > 0)
+ else
+ cross = true
+
+ } else {
+ # Intersection is entire line segment.
+ p1 = x1; p2 = x2
+ cross = (((p_y - y) * (n_y - y)) < 0)
+ }
+
+ i1 = max(1, min(axlen, nint(p1)))
+ i2 = max(1, min(axlen, nint(p2)))
+ if (i1 > i2)
+ swapi (i1, i2)
+
+ np = i2 - i1 + 1
+ if (np > 0) {
+ RL_X(rl_reg,rn) = i1
+ RL_N(rl_reg,rn) = np
+ RL_V(rl_reg,rn) = btoi(cross)
+ rn = rn + 1
+ }
+ }
+ }
+
+ rl_len = rn - 1
+ if (rl_len <= RL_FIRST)
+ goto done_
+
+ # Sort the line intersection-segments in order of increasing X.
+ do j = RL_FIRST, rl_len {
+ # Get low X value of initial segment.
+ i1 = RL_X(rl_reg,j)
+ np = RL_N(rl_reg,j)
+ i1 = min (i1, i1 + np - 1)
+ low = j
+
+ # Find lowest valued segment in remainder of array.
+ do i = j+1, rl_len {
+ i2 = RL_X(rl_reg,i)
+ np = RL_N(rl_reg,i)
+ i2 = min (i2, i2 + np - 1)
+ if (i2 < i1) {
+ i1 = i2
+ low = i
+ }
+ }
+
+ # Interchange the initial segment and the low segment.
+ if (low != j) {
+ swapi (RL_X(rl_reg,j), RL_X(rl_reg,low))
+ swapi (RL_N(rl_reg,j), RL_N(rl_reg,low))
+ swapi (RL_V(rl_reg,j), RL_V(rl_reg,low))
+ }
+ }
+
+ # Combine any segments which overlap.
+ rn = RL_FIRST
+ do i = RL_FIRST + 1, rl_len {
+ i1 = RL_X(rl_reg,rn)
+ i2 = RL_N(rl_reg,rn) + i1 - 1
+ ii = RL_X(rl_reg,i)
+ if (ii >= i1 && ii <= i2) {
+ i2 = ii + RL_N(rl_reg,i) - 1
+ RL_N(rl_reg,rn) = max (RL_N(rl_reg,rn), i2 - i1 + 1)
+ RL_V(rl_reg,rn) = max (RL_V(rl_reg,rn), RL_V(rl_reg,i))
+ } else {
+ rn = rn + 1
+ RL_X(rl_reg,rn) = RL_X(rl_reg,i)
+ RL_N(rl_reg,rn) = RL_N(rl_reg,i)
+ RL_V(rl_reg,rn) = RL_V(rl_reg,i)
+ }
+ }
+ rl_len = rn
+
+ # Now combine successive pairs of intersections to produce the line
+ # segments to be drawn. If all points are crossing points (where the
+ # image line crosses the polygon boundary) then we draw a line between
+ # the first two points, then the second two points, and so on. Points
+ # where the image line touches the polygon boundary but does not cross
+ # it are plotted, but are not joined with other points to make line
+ # segments.
+
+ rn = RL_FIRST
+ ii = RL_FIRST
+
+ do j = RL_FIRST, rl_len {
+ if (j <= ii && j < rl_len) {
+ next
+
+ } else if (RL_V(rl_reg,ii) == YES) {
+ # Skip a vertext that touches but does not cross.
+ if (RL_V(rl_reg,j) == NO && j < rl_len)
+ next
+
+ # Draw a line between the two crossing points.
+ RL_X(rl_reg,rn) = RL_X(rl_reg,ii)
+ RL_N(rl_reg,rn) = max (RL_N(rl_reg,ii),
+ RL_X(rl_reg,j) + RL_N(rl_reg,j) - RL_X(rl_reg,ii))
+ RL_V(rl_reg,rn) = P_PV(ufd)
+ rn = rn + 1
+ ii = j + 1
+
+ } else {
+ # Plot only the first point.
+ RL_X(rl_reg,rn) = RL_X(rl_reg,ii)
+ RL_N(rl_reg,rn) = RL_N(rl_reg,ii)
+ RL_V(rl_reg,rn) = P_PV(ufd)
+ rn = rn + 1
+
+ if (j >= rl_len && j != ii) {
+ # Plot the second point, if and end of list.
+ RL_X(rl_reg,rn) = RL_X(rl_reg,j)
+ RL_N(rl_reg,rn) = RL_N(rl_reg,j)
+ RL_V(rl_reg,rn) = P_PV(ufd)
+ rn = rn + 1
+ } else
+ ii = j
+ }
+ }
+
+done_
+ # Convert the X values in the range list to be relative to the start
+ # of the list. Compute NPIX, the range in pixels spanned by the range
+ # list.
+
+ rl_len = rn - 1
+ xs = RL_X(rl_reg,RL_FIRST)
+ npix = RL_X(rl_reg,rl_len) + RL_N(rl_reg,rl_len) - xs
+
+ do i = RL_FIRST, rl_len
+ RL_X(rl_reg,i) = RL_X(rl_reg,i) - xs + 1
+
+ RL_LEN(rl_reg) = rl_len
+ RL_AXLEN(rl_reg) = npix
+
+ rl_new = true
+ if (P_OY(ufd) == line - 1)
+ rl_new = !plr_equali (rl_reg, Memi[P_OO(ufd)])
+ call amovi (rl_reg, Memi[P_OO(ufd)], rn - 1)
+ P_OY(ufd) = line
+
+ return (rl_new)
+end
+
+
+# PE_UCANNULUS -- Regionrop ufcn for a circular annulus clipped at the borders
+# of the mask.
+
+bool procedure pe_ucannulus (ufd, y, rl_reg, xs, npix)
+
+pointer ufd #I user function descriptor
+int y #I mask line number
+int rl_reg[3,ARB] #O output range list for line Y
+int xs #O first pixel to be edited
+int npix #O number of pixels affected
+
+real radius1, radius2, dx, dy
+pointer pl
+int rn, axlen, x1o, x1o_clipped, x2o, x2o_clipped, x1i, x1i_clipped
+int x2i, x2i_clipped
+
+begin
+ pl = CA_PL(ufd)
+ rn = RL_FIRST
+ axlen = PL_AXLEN(pl,1)
+ radius1 = CA_RADIUS1(ufd)
+ radius2 = CA_RADIUS2(ufd)
+
+ dy = abs (CA_YCEN(ufd) - y)
+ if (dy < radius2) {
+ dx = radius2 * radius2 - dy * dy
+ if (dx > 0.0)
+ dx = sqrt (dx)
+ else
+ dx = 0.0
+ x1o = int (CA_XCEN(ufd) - dx)
+ x2o = int (CA_XCEN(ufd) + dx)
+ x1o_clipped = max(1, min(axlen, x1o))
+ x2o_clipped = max(1, min(axlen, x2o))
+ xs = x1o_clipped
+ if (dy < radius1) {
+ dx = radius1 * radius1 - dy * dy
+ if (dx > 0.0)
+ dx = sqrt (dx)
+ else
+ dx = 0.0
+ x1i = int (CA_XCEN(ufd) - dx)
+ x2i = int (CA_XCEN(ufd) + dx)
+ x1i_clipped = max(1, min (axlen, x1i))
+ x2i_clipped = max(1, min (axlen, x2i))
+ RL_X(rl_reg,rn) = 1
+ RL_N(rl_reg,rn) = x1i_clipped - x1o_clipped + 1
+ RL_V(rl_reg,rn) = CA_PV(ufd)
+ rn = rn + 1
+ RL_X(rl_reg,rn) = x2i_clipped - x1o_clipped + 1
+ RL_N(rl_reg,rn) = x2o_clipped - x2i_clipped + 1
+ RL_V(rl_reg,rn) = CA_PV(ufd)
+ rn = rn + 1
+ npix = x2o_clipped - x1o_clipped + 1
+ } else {
+ RL_X(rl_reg,rn) = 1
+ RL_N(rl_reg,rn) = x2o_clipped - x1o_clipped + 1
+ RL_V(rl_reg,rn) = CA_PV(ufd)
+ npix = RL_N(rl_reg,rn)
+ rn = rn + 1
+ }
+ } else {
+ npix = 0
+ xs = 1
+ }
+
+ RL_LEN(rl_reg) = rn - 1
+ RL_AXLEN(rl_reg) = npix
+
+ return (true)
+end
+
+
+# PE_UEANNULUS -- Regionrop ufcn for a circular annulus clipped at the borders
+# of the mask.
+
+bool procedure pe_ueannulus (ufd, y, rl_reg, xs, npix)
+
+pointer ufd #I user function descriptor
+int y #I mask line number
+int rl_reg[3,ARB] #O output range list for line Y
+int xs #O first pixel to be edited
+int npix #O number of pixels affected
+
+real dy, dy2, ady, bb2, cc2, bb1, cc1, discr, dx1, dx2
+pointer pl
+int rn, axlen, x1o, x1o_clipped, x2o, x2o_clipped, x1i, x1i_clipped
+int x2i, x2i_clipped
+
+begin
+ pl = EA_PL(ufd)
+ rn = RL_FIRST
+ axlen = PL_AXLEN(pl,1)
+
+ dy = y - EA_YCEN(ufd)
+ dy2 = dy * dy
+ ady = abs (dy)
+ bb2 = EA_BB2(ufd) * dy
+ cc2 = EA_CC2(ufd) * dy2
+ bb1 = EA_BB1(ufd) * dy
+ cc1 = EA_CC1(ufd) * dy2
+
+ if (ady < EA_DYMAX2(ufd)) {
+ discr = bb2 * bb2 - 4.0 * EA_AA2(ufd) * (cc2 - EA_FF2(ufd))
+ if (discr > 0.0)
+ discr = sqrt (discr)
+ else
+ discr = 0.0
+ dx1 = (-bb2 - discr) / 2.0 / EA_AA2(ufd)
+ dx2 = (-bb2 + discr) / 2.0 / EA_AA2(ufd)
+ x1o = EA_XCEN(ufd) + min (dx1, dx2)
+ x2o = EA_XCEN(ufd) + max (dx1, dx2)
+ x1o_clipped = max(1, min(axlen, x1o))
+ x2o_clipped = max(1, min(axlen, x2o))
+ xs = x1o_clipped
+ if (ady < EA_DYMAX1(ufd)) {
+ discr = bb1 * bb1 - 4.0 * EA_AA1(ufd) * (cc1 - EA_FF1(ufd))
+ if (discr > 0.0)
+ discr = sqrt (discr)
+ else
+ discr = 0.0
+ dx1 = (-bb1 - discr) / 2.0 / EA_AA1(ufd)
+ dx2 = (-bb1 + discr) / 2.0 / EA_AA1(ufd)
+ x1i = EA_XCEN(ufd) + min (dx1, dx2)
+ x2i = EA_XCEN(ufd) + max (dx1, dx2)
+ x1i_clipped = max(1, min (axlen, x1i))
+ x2i_clipped = max(1, min (axlen, x2i))
+ RL_X(rl_reg,rn) = 1
+ RL_N(rl_reg,rn) = x1i_clipped - x1o_clipped + 1
+ RL_V(rl_reg,rn) = EA_PV(ufd)
+ rn = rn + 1
+ RL_X(rl_reg,rn) = x2i_clipped - x1o_clipped + 1
+ RL_N(rl_reg,rn) = x2o_clipped - x2i_clipped + 1
+ RL_V(rl_reg,rn) = EA_PV(ufd)
+ rn = rn + 1
+ npix = x2o_clipped - x1o_clipped + 1
+ } else {
+ RL_X(rl_reg,rn) = 1
+ RL_N(rl_reg,rn) = x2o_clipped - x1o_clipped + 1
+ RL_V(rl_reg,rn) = EA_PV(ufd)
+ npix = RL_N(rl_reg,rn)
+ rn = rn + 1
+ }
+ } else {
+ npix = 0
+ xs = 1
+ }
+
+ RL_LEN(rl_reg) = rn - 1
+ RL_AXLEN(rl_reg) = npix
+
+ return (true)
+end
+
+
+# PE_UARECT -- Compute the intersection of an image line and a rectangular
+# polygonal annulus and define the region to be masked.
+
+bool procedure pe_uarect (ufd, y, rl_reg, xs, npix)
+
+pointer ufd #I the region descriptor structure
+int y #I the current line
+int rl_reg[3,ARB] #O the output regions list
+int xs #O the starting x value
+int npix #O the number of pixels affected
+
+real lx, ld
+pointer sp, work1, work2, oxintr, ixintr, pl
+int j, jj, rn, onintr, inintr, ix1, ix2, ox1, ox2, ibegin, iend, jx1, jx2
+int me_pyclip()
+
+begin
+ # Allocate working memory.
+ call smark (sp)
+ call salloc (work1, RA_NVER(ufd) + 1, TY_REAL)
+ call salloc (work2, RA_NVER(ufd) + 1, TY_REAL)
+ call salloc (oxintr, RA_NVER(ufd) + 1, TY_REAL)
+ call salloc (ixintr, RA_NVER(ufd) + 1, TY_REAL)
+
+ # Initialize.
+ pl = RA_PL(ufd)
+ rn = RL_FIRST
+ lx = PL_AXLEN(pl,1)
+ ld = y
+
+ # Find the intersection of the outer polygon with the image line.
+ onintr = me_pyclip (Memr[RA_OXP(ufd)], Memr[RA_OYP(ufd)], Memr[work1],
+ Memr[work2], Memr[oxintr], RA_NVER(ufd) + 1, lx, ld)
+ call asrtr (Memr[oxintr], Memr[oxintr], onintr)
+
+ if (onintr > 0) {
+
+ # Find the intersection of the inner polygon with the image line.
+ inintr = me_pyclip (Memr[RA_IXP(ufd)], Memr[RA_IYP(ufd)],
+ Memr[work1], Memr[work2], Memr[ixintr], RA_NVER(ufd) + 1,
+ lx, ld)
+ call asrtr (Memr[ixintr], Memr[ixintr], inintr)
+
+ # Create the region list.
+ xs = max (1, min (int(Memr[oxintr]), PL_AXLEN(pl,1)))
+ if (inintr <= 0) {
+ do j = 1, onintr, 2 {
+ ox1 = max (1, min (int(Memr[oxintr+j-1]), PL_AXLEN(pl,1)))
+ ox2 = max (ox1, min (int(Memr[oxintr+j]), PL_AXLEN(pl,1)))
+ RL_X(rl_reg,rn) = ox1 - xs + 1
+ RL_N(rl_reg,rn) = ox2 - ox1 + 1
+ RL_V(rl_reg,rn) = RA_PV(ufd)
+ rn = rn + 1
+ }
+ npix = RL_X(rl_reg, rn-1) + RL_N(rl_reg,rn-1) - 1
+ } else {
+ do j = 1, onintr, 2 {
+ ox1 = max (1, min (int(Memr[oxintr+j-1]), PL_AXLEN(pl,1)))
+ ox2 = max (ox1, min (int(Memr[oxintr+j]), PL_AXLEN(pl,1)))
+ do jj = 1, inintr, 2 {
+ ix1 = max (1, min (int(Memr[ixintr+jj-1]),
+ PL_AXLEN(pl,1)))
+ if (ix1 > ox1 && ix1 < ox2) {
+ ibegin = jj
+ break
+ }
+
+ }
+ do jj = inintr, 1, -2 {
+ ix2 = max (1, min (int(Memr[ixintr+jj-1]),
+ PL_AXLEN(pl,1)))
+ if (ix2 > ox1 && ix2 < ox2) {
+ iend = jj
+ break
+ }
+ }
+ RL_X(rl_reg,rn) = ox1 - xs + 1
+ RL_N(rl_reg,rn) = ix1 - ox1 + 1
+ RL_V(rl_reg,rn) = RA_PV(ufd)
+ rn = rn + 1
+ do jj = ibegin + 1, iend - 1, 2 {
+ jx1 = max (1, min (int(Memr[ixintr+jj-1]),
+ PL_AXLEN(pl,1)))
+ jx2 = max (jx1, min (int(Memr[ixintr+jj]),
+ PL_AXLEN(pl,1)))
+ RL_X(rl_reg,rn) = jx1 - xs + 1
+ RL_N(rl_reg,rn) = jx2 - jx1 + 1
+ RL_V(rl_reg,rn) = RA_PV(ufd)
+ rn = rn + 1
+ }
+ RL_X(rl_reg,rn) = ix2 - xs + 1
+ RL_N(rl_reg,rn) = ox2 - ix2 + 1
+ RL_V(rl_reg,rn) = RA_PV(ufd)
+ rn = rn + 1
+
+ }
+ npix = RL_X(rl_reg, rn-1) + RL_N(rl_reg,rn-1) - 1
+ }
+
+ } else {
+ xs = 1
+ npix = 0
+ }
+
+ RL_LEN(rl_reg) = rn - 1
+ RL_AXLEN(rl_reg) = npix
+
+ call sfree (sp)
+
+ return (true)
+end
+
+
+# PE_UAPOLYGON -- Compute the intersection of an image line and the polygonal
+# annulus and define the region to be masked.
+
+bool procedure pe_uapolygon (ufd, y, rl_reg, xs, npix)
+
+pointer ufd #I the region descriptor structure
+int y #I the current line
+int rl_reg[3,ARB] #O the output regions list
+int xs #O the starting x value
+int npix #O the number of pixels affected
+
+real lx, ld
+pointer sp, work1, work2, oxintr, ixintr, pl
+int j, jj, rn, onintr, inintr, ix1, ix2, ox1, ox2, ibegin, iend, jx1, jx2
+int me_pyclip()
+
+begin
+ # Allocate working memory.
+ call smark (sp)
+ call salloc (work1, PA_NVER(ufd) + 1, TY_REAL)
+ call salloc (work2, PA_NVER(ufd) + 1, TY_REAL)
+ call salloc (oxintr, PA_NVER(ufd) + 1, TY_REAL)
+ call salloc (ixintr, PA_NVER(ufd) + 1, TY_REAL)
+
+ # Initialize.
+ pl = PA_PL(ufd)
+ rn = RL_FIRST
+ lx = PL_AXLEN(pl,1)
+ ld = y
+
+ # Find the intersection of the outer polygon with the image line.
+ onintr = me_pyclip (Memr[PA_OXP(ufd)], Memr[PA_OYP(ufd)], Memr[work1],
+ Memr[work2], Memr[oxintr], PA_NVER(ufd) + 1, lx, ld)
+ call asrtr (Memr[oxintr], Memr[oxintr], onintr)
+
+ if (onintr > 0) {
+
+ # Find the intersection of the inner polygon with the image line.
+ inintr = me_pyclip (Memr[PA_IXP(ufd)], Memr[PA_IYP(ufd)],
+ Memr[work1], Memr[work2], Memr[ixintr], PA_NVER(ufd) + 1,
+ lx, ld)
+ call asrtr (Memr[ixintr], Memr[ixintr], inintr)
+
+ # Create the region list.
+ xs = max (1, min (int(Memr[oxintr]), PL_AXLEN(pl,1)))
+ if (inintr <= 0) {
+ do j = 1, onintr, 2 {
+ ox1 = max (1, min (int(Memr[oxintr+j-1]), PL_AXLEN(pl,1)))
+ ox2 = max (ox1, min (int(Memr[oxintr+j]), PL_AXLEN(pl,1)))
+ RL_X(rl_reg,rn) = ox1 - xs + 1
+ RL_N(rl_reg,rn) = ox2 - ox1 + 1
+ RL_V(rl_reg,rn) = PA_PV(ufd)
+ rn = rn + 1
+ }
+ npix = RL_X(rl_reg, rn-1) + RL_N(rl_reg,rn-1) - 1
+ } else {
+ do j = 1, onintr, 2 {
+ ox1 = max (1, min (int(Memr[oxintr+j-1]), PL_AXLEN(pl,1)))
+ ox2 = max (ox1, min (int(Memr[oxintr+j]), PL_AXLEN(pl,1)))
+ do jj = 1, inintr, 2 {
+ ix1 = max (1, min (int(Memr[ixintr+jj-1]),
+ PL_AXLEN(pl,1)))
+ if (ix1 > ox1 && ix1 < ox2) {
+ ibegin = jj
+ break
+ }
+
+ }
+ do jj = inintr, 1, -2 {
+ ix2 = max (1, min (int(Memr[ixintr+jj-1]),
+ PL_AXLEN(pl,1)))
+ if (ix2 > ox1 && ix2 < ox2) {
+ iend = jj
+ break
+ }
+ }
+ RL_X(rl_reg,rn) = ox1 - xs + 1
+ RL_N(rl_reg,rn) = ix1 - ox1 + 1
+ RL_V(rl_reg,rn) = PA_PV(ufd)
+ rn = rn + 1
+ do jj = ibegin + 1, iend - 1, 2 {
+ jx1 = max (1, min (int(Memr[ixintr+jj-1]),
+ PL_AXLEN(pl,1)))
+ jx2 = max (jx1, min (int(Memr[ixintr+jj]),
+ PL_AXLEN(pl,1)))
+ RL_X(rl_reg,rn) = jx1 - xs + 1
+ RL_N(rl_reg,rn) = jx2 - jx1 + 1
+ RL_V(rl_reg,rn) = PA_PV(ufd)
+ rn = rn + 1
+ }
+ RL_X(rl_reg,rn) = ix2 - xs + 1
+ RL_N(rl_reg,rn) = ox2 - ix2 + 1
+ RL_V(rl_reg,rn) = PA_PV(ufd)
+ rn = rn + 1
+
+ }
+ npix = RL_X(rl_reg, rn-1) + RL_N(rl_reg,rn-1) - 1
+ }
+
+ } else {
+ xs = 1
+ npix = 0
+ }
+
+ RL_LEN(rl_reg) = rn - 1
+ RL_AXLEN(rl_reg) = npix
+
+ call sfree (sp)
+
+ return (true)
+end
+
+
+# PE_UCOLS -- Regionrop ufcn for a set of column ranges (column regions),
+# clipped at the borders of the mask.
+
+bool procedure pe_ucols (ufd, y, rl_reg, xs, npix)
+
+pointer ufd #I user function descriptor
+int y #I mask line number
+int rl_reg[3,ARB] #O output range list for line Y
+int xs #O first pixel to be edited
+int npix #O number of pixels affected
+
+begin
+ # Copy the ranges.
+ call amovi (Memi[L_RANGES(ufd)], rl_reg, L_NRANGES(ufd) * 3)
+ xs = L_XS(ufd)
+ npix = L_NPIX(ufd)
+
+ return (true)
+end
+
+
+# PE_ULINES -- Regionrop ufcn for a set of lines ranges (line regions),
+# clipped at the borders of the mask.
+
+bool procedure pe_ulines (ufd, y, rl_reg, xs, npix)
+
+pointer ufd #I user function descriptor
+int y #I mask line number
+int rl_reg[3,ARB] #O output range list for line Y
+int xs #O first pixel to be edited
+int npix #O number of pixels affected
+
+pointer pl
+int rn, axlen
+bool me_is_in_range()
+
+begin
+ pl = L_PL(ufd)
+ rn = RL_FIRST
+ axlen = PL_AXLEN(pl,1)
+
+ if (me_is_in_range (Memi[L_RANGES(ufd)], y)) {
+ xs = 1
+ npix = axlen
+ RL_X(rl_reg,rn) = 1
+ RL_N(rl_reg,rn) = axlen
+ RL_V(rl_reg,rn) = L_PV(ufd)
+ rn = rn + 1
+ } else {
+ xs = 1
+ npix = 0
+ }
+
+ RL_LEN(rl_reg) = rn - 1
+ RL_AXLEN(rl_reg) = npix
+
+ return (true)
+end