aboutsummaryrefslogtreecommitdiff
path: root/pkg/proto/maskexpr
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/proto/maskexpr')
-rw-r--r--pkg/proto/maskexpr/gettok.h22
-rw-r--r--pkg/proto/maskexpr/gettok.x922
-rw-r--r--pkg/proto/maskexpr/megeom.x72
-rw-r--r--pkg/proto/maskexpr/megsym.x31
-rw-r--r--pkg/proto/maskexpr/memkmask.x839
-rw-r--r--pkg/proto/maskexpr/meregfuncs.x1449
-rw-r--r--pkg/proto/maskexpr/meregmask.x753
-rw-r--r--pkg/proto/maskexpr/mesetexpr.x36
-rw-r--r--pkg/proto/maskexpr/mesetreg.x292
-rw-r--r--pkg/proto/maskexpr/mkpkg26
-rw-r--r--pkg/proto/maskexpr/mskexpand.x261
-rw-r--r--pkg/proto/maskexpr/peregfuncs.h131
-rw-r--r--pkg/proto/maskexpr/peregfuncs.x877
-rw-r--r--pkg/proto/maskexpr/peregufcn.x808
-rw-r--r--pkg/proto/maskexpr/t_mskexpr.x286
-rw-r--r--pkg/proto/maskexpr/t_mskregions.x264
16 files changed, 7069 insertions, 0 deletions
diff --git a/pkg/proto/maskexpr/gettok.h b/pkg/proto/maskexpr/gettok.h
new file mode 100644
index 00000000..90980fa1
--- /dev/null
+++ b/pkg/proto/maskexpr/gettok.h
@@ -0,0 +1,22 @@
+# GETTOK.H -- External definitions for gettok.h
+
+define GT_IDENT (-99)
+define GT_NUMBER (-98)
+define GT_STRING (-97)
+define GT_COMMAND (-96)
+define GT_PLUSEQ (-95)
+define GT_COLONEQ (-94)
+define GT_EXPON (-93)
+define GT_CONCAT (-92)
+define GT_SE (-91)
+define GT_LE (-90)
+define GT_GE (-89)
+define GT_EQ (-88)
+define GT_NE (-87)
+define GT_LAND (-86)
+define GT_LOR (-85)
+
+# Optional flags.
+define GT_NOSPECIAL 0003
+define GT_NOFILE 0001
+define GT_NOCOMMAND 0002
diff --git a/pkg/proto/maskexpr/gettok.x b/pkg/proto/maskexpr/gettok.x
new file mode 100644
index 00000000..a0975300
--- /dev/null
+++ b/pkg/proto/maskexpr/gettok.x
@@ -0,0 +1,922 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <syserr.h>
+include <error.h>
+include <ctype.h>
+include <fset.h>
+include "gettok.h"
+
+.help gettok
+.nf --------------------------------------------------------------------------
+GETTOK -- Lexical input routines. Used to return tokens from input text,
+performing macro expansion and file expansion. The input text may be either
+an open file descriptor or a text string.
+
+ nchars = gt_expandtext (text, obuf, len_obuf, gsym, gsym_data)
+
+ gt = gt_open (fd, gsym, gsym_data, pbblen, flags)
+ gt = gt_opentext (text, gsym, gsym_data, pbblen, flags)
+ gt_close (gt)
+
+ nchars = gt_expand (gt, obuf, len_obuf)
+ token = gt_gettok (gt, tokbuf, maxch)
+ gt_ungettok (gt, tokbuf)
+ token = gt_rawtok (gt, tokbuf, maxch)
+ token = gt_nexttok (gt)
+
+The client get-symbol routine has the following calling sequence, where
+"nargs" is an output argument which should be set to the number of macro
+arguments, if any. Normally this routine will call SYMTAB to do the
+symbol lookup, but this is not required. GSYM may be set to NULL if no
+macro replacement is desired.
+
+ textp = gsym (gsym_data, symbol, &nargs)
+
+PBBLEN is the size of the pushback buffer used for macro expansion, and
+determines the size of the largest macro replacement string that can be
+pushed back. FLAGS may be used to disable certain types of pushback.
+Both PBBLEN and FLAGS may be given as zero if the client is happy with the
+builtin defaults.
+
+Access to the package is gained by opening a text string with GT_OPENTEXT.
+This returns a descriptor which is passed to GT_GETTOK to read successive
+tokens, which may come from the input text string or from any macros,
+include files, etc., referenced in the text or in any substituted text.
+GT_UNGETTOK pushes a token back into the GT_GETTOK input stream, to be
+returned in the next GT_GETTOK call (following macro expansion). GT_EXPAND
+will process the entire input text string, expanding any macro references
+therein, returning the fully resolved text in the output buffer. A more
+macroscopic version of this is GT_EXPANDTEXT, which does the opentext,
+expand, and close operations internally, using the builtin defaults.
+
+GT_RAWTOK returns the next physical token from an input stream (without
+macro expansion), and GT_NEXTTOK returns the type of the next *physical*
+token (no macro expansion) without actually fetching it (for look ahead
+decision making).
+
+The tokens that can be returned are as follows:
+
+ GT_IDENT [a-zA-Z][a-zA-Z0-9_]*
+ GT_NUMBER [0-9][0-9a-zA-Z.]*(e|E)?(+|-)?[0-9]*
+ GT_STRING if "abc" or 'abc', the abc
+ 'c' other characters, e.g., =+-*/,;:()[] etc
+ EOF at end of input
+
+Macro replacement syntax:
+
+ macro push macro with null arglist
+ macro(arg,arg,...) push macro with argument substitution
+ @file push contents of file
+ @file(arg,arg,...) push file with argument substitution
+ `cmd` substitute output of CL command "cmd"
+
+where
+ macro is an identifier, the name of a global macro
+ or a datafile local macro (parameter)
+
+In all cases, occurences of $N in the replacement text are replaced by the
+macro arguments if any, and macros are recursively expanded. Whitespace,
+including newline, equates to a single space, as does EOF (hence always
+delimits tokens). Comments (# to end of line) are ignored. All identifiers
+in scanned text are checked to see if they are references to predefined
+macros, using the client supplied symbol lookup routine.
+.endhelp ---------------------------------------------------------------------
+
+# General definitions.
+define MAX_LEVELS 20 # max include file nesting
+define MAX_ARGS 9 # max arguments to a macro
+define SZ_CMD 80 # `cmd`
+define SZ_IBUF 8192 # buffer for macro replacement
+define SZ_OBUF 8192 # buffer for macro replacement
+define SZ_ARGBUF 256 # argument list to a macro
+define SZ_TOKBUF 1024 # token buffer
+define DEF_MAXPUSHBACK 16384 # max pushback, macro replacement
+define INC_TOKBUF 4096 # increment if expanded text fills
+
+# The gettok descriptor.
+define LEN_GTDES 50
+define GT_FD Memi[$1] # current input stream
+define GT_UFD Memi[$1+1] # user (client) input file
+define GT_FLAGS Memi[$1+2] # option flags
+define GT_PBBLEN Memi[$1+3] # pushback buffer length
+define GT_DEBUG Memi[$1+4] # for debug messages
+define GT_GSYM Memi[$1+5] # get symbol routine
+define GT_GSYMDATA Memi[$1+6] # client data for above
+define GT_NEXTCH Memi[$1+7] # lookahead character
+define GT_FTEMP Memi[$1+8] # file on stream is a temp file
+define GT_LEVEL Memi[$1+9] # current nesting level
+define GT_SVFD Memi[$1+10+$2-1]# stacked file descriptors
+define GT_SVFTEMP Memi[$1+30+$2-1]# stacked ftemp flags
+
+# Set to YES to enable debug messages.
+define DEBUG NO
+
+
+# GT_EXPANDTEXT -- Perform macro expansion on a text string returning the
+# fully resolved text in the client's output buffer. The number of chars
+# in the output string is returned as the function value.
+
+int procedure gt_expandtext (text, obuf, len_obuf, gsym, gsym_data)
+
+char text[ARB] #I input text to be expanded
+pointer obuf #U output buffer
+int len_obuf #U size of output buffer
+int gsym #I epa of client get-symbol routine
+int gsym_data #I client data for above
+
+pointer gt
+int nchars
+int gt_expand()
+pointer gt_opentext()
+errchk gt_opentext
+
+begin
+ gt = gt_opentext (text, gsym, gsym_data, 0, 0)
+ nchars = gt_expand (gt, obuf, len_obuf)
+ call gt_close (gt)
+
+ return (nchars)
+end
+
+
+# GT_EXPAND -- Perform macro expansion on a GT text stream returning the
+# fully resolved text in the client's output buffer. The number of chars
+# in the output string is returned as the function value.
+
+int procedure gt_expand (gt, obuf, len_obuf)
+
+pointer gt #I gettok descriptor
+pointer obuf #U output buffer
+int len_obuf #U size of output buffer
+
+int token, nchars
+pointer sp, tokbuf, op, otop
+int gt_gettok(), strlen(), gstrcpy()
+errchk realloc
+
+begin
+ call smark (sp)
+ call salloc (tokbuf, SZ_TOKBUF, TY_CHAR)
+
+ # Open input text for macro expanded token input.
+ otop = obuf + len_obuf
+ op = obuf
+
+ # Copy tokens to the output, inserting a space after every token.
+ repeat {
+ token = gt_gettok (gt, Memc[tokbuf], SZ_TOKBUF)
+ if (token != EOF) {
+ if (op + strlen(Memc[tokbuf]) + 3 > otop) {
+ nchars = op - obuf
+ len_obuf = len_obuf + INC_TOKBUF
+ call realloc (obuf, len_obuf, TY_CHAR)
+ otop = obuf + len_obuf
+ op = obuf + nchars
+ }
+
+ if (token == GT_STRING) {
+ Memc[op] = '"'
+ op = op + 1
+ }
+ op = op + gstrcpy (Memc[tokbuf], Memc[op], otop-op)
+ if (token == GT_STRING) {
+ Memc[op] = '"'
+ op = op + 1
+ }
+ Memc[op] = ' '
+ op = op + 1
+ }
+ } until (token == EOF)
+
+ # Cancel the trailing blank and add the EOS.
+ if (op > 1 && op < otop)
+ op = op - 1
+ Memc[op] = EOS
+
+ call sfree (sp)
+ return (op - 1)
+end
+
+
+# GT_OPEN -- Open the GETTOK descriptor on a file descriptor.
+
+pointer procedure gt_open (fd, gsym, gsym_data, pbblen, flags)
+
+int fd #I input file
+int gsym #I epa of client get-symbol routine
+int gsym_data #I client data for above
+int pbblen #I pushback buffer length
+int flags #I option flags
+
+pointer gt
+int sz_pbbuf
+errchk calloc
+
+begin
+ call calloc (gt, LEN_GTDES, TY_STRUCT)
+
+ GT_GSYM(gt) = gsym
+ GT_GSYMDATA(gt) = gsym_data
+ GT_FLAGS(gt) = flags
+ GT_DEBUG(gt) = DEBUG
+
+ GT_FD(gt) = fd
+ GT_UFD(gt) = fd
+
+ if (pbblen <= 0)
+ sz_pbbuf = DEF_MAXPUSHBACK
+ else
+ sz_pbbuf = pbblen
+ call fseti (GT_FD(gt), F_PBBSIZE, sz_pbbuf)
+ GT_PBBLEN(gt) = sz_pbbuf
+
+ return (gt)
+end
+
+
+# GT_OPENTEXT -- Open the GT_GETTOK descriptor. The descriptor is initially
+# opened on the user supplied string buffer (which is opened as a file and
+# which must remain intact while token input is in progress), but include file
+# processing etc. may cause arbitrary nesting of file descriptors.
+
+pointer procedure gt_opentext (text, gsym, gsym_data, pbblen, flags)
+
+char text[ARB] #I input text to be scanned
+int gsym #I epa of client get-symbol routine
+int gsym_data #I client data for above
+int pbblen #I pushback buffer length
+int flags #I option flags
+
+pointer gt
+int sz_pbbuf
+int stropen(), strlen()
+errchk stropen, calloc
+
+begin
+ call calloc (gt, LEN_GTDES, TY_STRUCT)
+
+ GT_GSYM(gt) = gsym
+ GT_GSYMDATA(gt) = gsym_data
+ GT_FLAGS(gt) = flags
+ GT_DEBUG(gt) = DEBUG
+
+ GT_FD(gt) = stropen (text, strlen(text), READ_ONLY)
+ GT_UFD(gt) = 0
+
+ if (pbblen <= 0)
+ sz_pbbuf = DEF_MAXPUSHBACK
+ else
+ sz_pbbuf = pbblen
+ call fseti (GT_FD(gt), F_PBBSIZE, sz_pbbuf)
+ GT_PBBLEN(gt) = sz_pbbuf
+
+ return (gt)
+end
+
+
+# GT_GETTOK -- Return the next token from the input stream. The token ID
+# (a predefined integer code or the character value) is returned as the
+# function value. The text of the token is returned as an output argument.
+# Any macro references, file includes, etc., are performed in the process
+# of scanning the input stream, hence only fully resolved tokens are output.
+
+int procedure gt_gettok (gt, tokbuf, maxch)
+
+pointer gt #I gettok descriptor
+char tokbuf[maxch] #O receives the text of the token
+int maxch #I max chars out
+
+pointer sp, bp, cmd, ibuf, obuf, argbuf, fname, textp
+int fd, token, level, margs, nargs, nchars, i_fd, o_fd, ftemp
+
+int strmac(), open(), stropen()
+int gt_rawtok(), gt_nexttok(), gt_arglist(), zfunc3()
+errchk gt_rawtok, close, ungetci, ungetline, gt_arglist,
+errchk clcmdw, stropen, syserr, zfunc3
+define pushfile_ 91
+
+
+begin
+ call smark (sp)
+
+ # Allocate some buffer space.
+ nchars = SZ_CMD + SZ_IBUF + SZ_OBUF + SZ_ARGBUF + SZ_FNAME + 5
+ call salloc (bp, nchars, TY_CHAR)
+
+ cmd = bp
+ ibuf = cmd + SZ_CMD + 1
+ obuf = ibuf + SZ_IBUF + 1
+ argbuf = obuf + SZ_OBUF + 1
+ fname = argbuf + SZ_ARGBUF + 1
+
+ # Read raw tokens and push back macro or include file text until we
+ # get a fully resolved token.
+
+ repeat {
+ fd = GT_FD(gt)
+
+ # Get a raw token.
+ token = gt_rawtok (gt, tokbuf, maxch)
+
+ # Process special tokens.
+ switch (token) {
+ case EOF:
+ # EOF has been reached on the current stream.
+ level = GT_LEVEL(gt)
+ if (GT_FTEMP(gt) == YES) {
+ call fstats (fd, F_FILENAME, Memc[fname], SZ_FNAME)
+ if (level > 0)
+ call close (fd)
+ iferr (call delete (Memc[fname]))
+ call erract (EA_WARN)
+ } else if (level > 0)
+ call close (fd)
+
+ if (level > 0) {
+ # Restore previous stream.
+ GT_FD(gt) = GT_SVFD(gt,level)
+ GT_FTEMP(gt) = GT_SVFTEMP(gt,level)
+ GT_LEVEL(gt) = level - 1
+ GT_NEXTCH(gt) = NULL
+ } else {
+ # Return EOF token to caller.
+ call strcpy ("EOF", tokbuf, maxch)
+ break
+ }
+
+ case GT_IDENT:
+ # Lookup the identifier in the symbol table.
+ textp = NULL
+ if (GT_GSYM(gt) != NULL)
+ textp = zfunc3 (GT_GSYM(gt), GT_GSYMDATA(gt), tokbuf, margs)
+
+ # Process a defined macro.
+ if (textp != NULL) {
+ # If macro does not have any arguments, merely push back
+ # the replacement text.
+
+ if (margs == 0) {
+ if (GT_NEXTCH(gt) > 0) {
+ call ungetci (fd, GT_NEXTCH(gt))
+ GT_NEXTCH(gt) = 0
+ }
+ call ungetline (fd, Memc[textp])
+ next
+ }
+
+ # Extract argument list, if any, perform argument
+ # substitution on the macro, and push back the edited
+ # text to be rescanned.
+
+ if (gt_nexttok(gt) == '(') {
+ nargs = gt_arglist (gt, Memc[argbuf], SZ_ARGBUF)
+ if (nargs != margs) {
+ call eprintf ("macro `%s' called with ")
+ call pargstr (tokbuf)
+ call eprintf ("wrong number of arguments\n")
+ }
+
+ # Pushback the text of a macro with arg substitution.
+ nchars = strmac (Memc[textp], Memc[argbuf],
+ Memc[obuf], SZ_OBUF)
+ if (GT_NEXTCH(gt) > 0) {
+ call ungetci (fd, GT_NEXTCH(gt))
+ GT_NEXTCH(gt) = 0
+ }
+ call ungetline (fd, Memc[obuf])
+ next
+
+ } else {
+ call eprintf ("macro `%s' called with no arguments\n")
+ call pargstr (tokbuf)
+ }
+ }
+
+ # Return a regular identifier.
+ break
+
+ case GT_COMMAND:
+ # Send a command to the CL and push back the output.
+ if (and (GT_FLAGS(gt), GT_NOCOMMAND) != 0)
+ break
+
+ # Execute the command, spooling the output in a temp file.
+ call mktemp ("tmp$co", Memc[fname], SZ_FNAME)
+ call sprintf (Memc[cmd], SZ_LINE, "%s > %s")
+ call pargstr (tokbuf)
+ call pargstr (Memc[fname])
+ call clcmdw (Memc[cmd])
+
+ # Open the output file as input text.
+ call strcpy (Memc[fname], tokbuf, maxch)
+ nargs = 0
+ ftemp = YES
+ goto pushfile_
+
+ case '@':
+ # Pushback the contents of a file.
+ if (and (GT_FLAGS(gt), GT_NOFILE) != 0)
+ break
+
+ token = gt_rawtok (gt, tokbuf, maxch)
+ if (token != GT_IDENT && token != GT_STRING) {
+ call eprintf ("expected a filename after the `@'\n")
+ next
+ } else {
+ nargs = 0
+ if (gt_nexttok(gt) == '(') # )
+ nargs = gt_arglist (gt, Memc[argbuf], SZ_ARGBUF)
+ ftemp = NO
+ }
+pushfile_
+ # Attempt to open the file.
+ iferr (i_fd = open (tokbuf, READ_ONLY, TEXT_FILE)) {
+ call eprintf ("cannot open `%s'\n")
+ call pargstr (tokbuf)
+ next
+ }
+
+ call fseti (i_fd, F_PBBSIZE, GT_PBBLEN(gt))
+
+ # Cancel lookahead.
+ if (GT_NEXTCH(gt) > 0) {
+ call ungetci (fd, GT_NEXTCH(gt))
+ GT_NEXTCH(gt) = 0
+ }
+
+ # If the macro was called with a nonnull argument list,
+ # attempt to perform argument substitution on the file
+ # contents. Otherwise merely push the fd.
+
+ if (nargs > 0) {
+ # Pushback file contents with argument substitution.
+ o_fd = stropen (Memc[ibuf], SZ_IBUF, NEW_FILE)
+
+ call fcopyo (i_fd, o_fd)
+ nchars = strmac (Memc[ibuf],Memc[argbuf],Memc[obuf],SZ_OBUF)
+ call ungetline (fd, Memc[obuf])
+
+ call close (o_fd)
+ call close (i_fd)
+
+ } else {
+ # Push a new input stream.
+ level = GT_LEVEL(gt) + 1
+ if (level > MAX_LEVELS)
+ call syserr (SYS_FPBOVFL)
+
+ GT_SVFD(gt,level) = GT_FD(gt)
+ GT_SVFTEMP(gt,level) = GT_FTEMP(gt)
+ GT_LEVEL(gt) = level
+
+ fd = i_fd
+ GT_FD(gt) = fd
+ GT_FTEMP(gt) = ftemp
+ }
+
+ default:
+ break
+ }
+ }
+
+ if (GT_DEBUG(gt) > 0) {
+ call eprintf ("token=%d(%o), `%s'\n")
+ call pargi (token)
+ call pargi (max(0,token))
+ if (IS_PRINT(tokbuf[1]))
+ call pargstr (tokbuf)
+ else
+ call pargstr ("")
+ }
+
+ call sfree (sp)
+ return (token)
+end
+
+
+# GT_UNGETTOK -- Push a token back into the GT_GETTOK input stream, to be
+# returned as the next token by GT_GETTOK.
+
+procedure gt_ungettok (gt, tokbuf)
+
+pointer gt #I gettok descriptor
+char tokbuf[ARB] #I text of token
+
+int fd
+errchk ungetci
+
+begin
+ fd = GT_FD(gt)
+
+ if (GT_DEBUG(gt) > 0) {
+ call eprintf ("unget token `%s'\n")
+ call pargstr (tokbuf)
+ }
+
+ # Cancel lookahead.
+ if (GT_NEXTCH(gt) > 0) {
+ call ungetci (fd, GT_NEXTCH(gt))
+ GT_NEXTCH(gt) = 0
+ }
+
+ # First push back a space to ensure that the token is recognized
+ # when the input is rescanned.
+
+ call ungetci (fd, ' ')
+
+ # Now push the token text.
+ call ungetline (fd, tokbuf)
+end
+
+
+# GT_RAWTOK -- Get a raw token from the input stream, without performing any
+# macro expansion or file inclusion. The text of the token in returned in
+# tokbuf, and the token type is returened as the function value.
+
+int procedure gt_rawtok (gt, outstr, maxch)
+
+pointer gt #I gettok descriptor
+char outstr[maxch] #O receives text of token.
+int maxch #I max chars out
+
+int token, delim, fd, ch, last_ch, op
+define again_ 91
+int getci()
+
+begin
+ fd = GT_FD(gt)
+again_
+ # Get lookahead char if we don't already have one.
+ ch = GT_NEXTCH(gt)
+ GT_NEXTCH(gt) = NULL
+ if (ch <= 0 || IS_WHITE(ch) || ch == '\n') {
+ while (getci (fd, ch) != EOF)
+ if (!(IS_WHITE(ch) || ch == '\n'))
+ break
+ }
+
+ # Output the first character.
+ op = 1
+ if (ch != EOF && ch != '"' && ch != '\'' && ch != '`') {
+ outstr[op] = ch
+ op = op + 1
+ }
+
+ # Accumulate token. Some of the token recognition logic used here
+ # (especially for numbers) is crude, but it is not clear that rigour
+ # is justified for this application.
+
+ if (ch == EOF) {
+ call strcpy ("EOF", outstr, maxch)
+ token = EOF
+
+ } else if (ch == '#') {
+ # Ignore a comment.
+ while (getci (fd, ch) != '\n')
+ if (ch == EOF)
+ break
+ goto again_
+
+ } else if (IS_ALPHA(ch) || ch == '_' || ch == '$' || ch == '.') {
+ # Identifier.
+ token = GT_IDENT
+ while (getci (fd, ch) != EOF)
+ if (IS_ALNUM(ch) || ch == '_' || ch == '$' || ch == '.') {
+ outstr[op] = ch
+ op = min (maxch, op+1)
+ } else
+ break
+
+ } else if (IS_DIGIT(ch)) {
+ # Number.
+ token = GT_NUMBER
+
+ # Get number.
+ while (getci (fd, ch) != EOF)
+ if (IS_ALNUM(ch) || ch == '.') {
+ outstr[op] = ch
+ last_ch = ch
+ op = min (maxch, op+1)
+ } else
+ break
+
+ # Get exponent if any.
+ if (last_ch == 'E' || last_ch == 'e') {
+ outstr[op] = ch
+ op = min (maxch, op+1)
+ while (getci (fd, ch) != EOF)
+ if (IS_DIGIT(ch) || ch == '+' || ch == '-') {
+ outstr[op] = ch
+ op = min (maxch, op+1)
+ } else
+ break
+ }
+
+ } else if (ch == '"' || ch == '\'' || ch == '`') {
+ # Quoted string or command.
+
+ if (ch == '`')
+ token = GT_COMMAND
+ else
+ token = GT_STRING
+
+ delim = ch
+ while (getci (fd, ch) != EOF)
+ if (ch==delim && (op>1 && outstr[op-1] != '\\') || ch == '\n')
+ break
+ else {
+ outstr[op] = ch
+ op = min (maxch, op+1)
+ }
+ ch = getci (fd, ch)
+
+ } else if (ch == '+') {
+ # May be the += operator.
+ if (getci (fd, ch) != EOF)
+ if (ch == '=') {
+ token = GT_PLUSEQ
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = '+'
+
+ } else if (ch == ':') {
+ # May be the := operator.
+ if (getci (fd, ch) != EOF)
+ if (ch == '=') {
+ token = GT_COLONEQ
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = ':'
+
+ } else if (ch == '*') {
+ if (getci (fd, ch) != EOF)
+ if (ch == '*') {
+ token = GT_EXPON
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = '*'
+
+ } else if (ch == '/') {
+ if (getci (fd, ch) != EOF)
+ if (ch == '/') {
+ token = GT_CONCAT
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = '/'
+
+ } else if (ch == '?') {
+ if (getci (fd, ch) != EOF)
+ if (ch == '=') {
+ token = GT_SE
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = '?'
+
+ } else if (ch == '<') {
+ if (getci (fd, ch) != EOF)
+ if (ch == '=') {
+ token = GT_LE
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = '<'
+
+ } else if (ch == '>') {
+ if (getci (fd, ch) != EOF)
+ if (ch == '=') {
+ token = GT_GE
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = '>'
+
+ } else if (ch == '=') {
+ if (getci (fd, ch) != EOF)
+ if (ch == '=') {
+ token = GT_EQ
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = '='
+
+ } else if (ch == '!') {
+ if (getci (fd, ch) != EOF)
+ if (ch == '=') {
+ token = GT_NE
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = '!'
+
+ } else if (ch == '&') {
+ if (getci (fd, ch) != EOF)
+ if (ch == '&') {
+ token = GT_LAND
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = '&'
+
+ } else if (ch == '|') {
+ if (getci (fd, ch) != EOF)
+ if (ch == '|') {
+ token = GT_LOR
+ outstr[op] = ch
+ op = op + 1
+ ch = getci (fd, ch)
+ } else
+ token = '|'
+
+ } else {
+ # Other characters.
+ token = ch
+ ch = getci (fd, ch)
+ }
+
+ # Process the lookahead character.
+ if (IS_WHITE(ch) || ch == '\n') {
+ repeat {
+ ch = getci (fd, ch)
+ } until (!(IS_WHITE(ch) || ch == '\n'))
+ }
+
+ if (ch != EOF)
+ GT_NEXTCH(gt) = ch
+
+ outstr[op] = EOS
+ return (token)
+end
+
+
+# GT_NEXTTOK -- Determine the type of the next raw token in the input stream,
+# without actually fetching the token. Operators such as GT_EQ etc. are not
+# recognized at this level. Note that this is at the same level as
+# GT_RAWTOK, i.e., no macro expansion is performed, and the lookahead token
+# is that which would be returned by the next gt_rawtok, which is not
+# necessarily what gt_gettok would return after macro replacement.
+
+int procedure gt_nexttok (gt)
+
+pointer gt #I gettok descriptor
+
+int token, fd, ch
+int getci()
+
+begin
+ fd = GT_FD(gt)
+
+ # Get lookahead char if we don't already have one.
+ ch = GT_NEXTCH(gt)
+ if (ch <= 0 || IS_WHITE(ch) || ch == '\n')
+ while (getci (fd, ch) != EOF)
+ if (!(IS_WHITE(ch) || ch == '\n'))
+ break
+
+ if (ch == EOF)
+ token = EOF
+ else if (IS_ALPHA(ch) || ch == '_' || ch == '$' || ch == '.')
+ token = GT_IDENT
+ else if (IS_DIGIT(ch))
+ token = GT_NUMBER
+ else if (ch == '"' || ch == '\'')
+ token = GT_STRING
+ else if (ch == '`')
+ token = GT_COMMAND
+ else
+ token = ch
+
+ if (GT_DEBUG(gt) > 0) {
+ call eprintf ("nexttok=%d(%o) `%c'\n")
+ call pargi (token)
+ call pargi (max(0,token))
+ if (IS_PRINT(ch))
+ call pargi (ch)
+ else
+ call pargi (0)
+ }
+
+ return (token)
+end
+
+
+# GT_CLOSE -- Close the gettok descriptor and any files opened thereon.
+
+procedure gt_close (gt)
+
+pointer gt #I gettok descriptor
+
+int level, fd
+pointer sp, fname
+
+begin
+ call smark (sp)
+ call salloc (fname, SZ_FNAME, TY_CHAR)
+
+ for (level=GT_LEVEL(gt); level >= 0; level=level-1) {
+ fd = GT_FD(gt)
+ if (GT_FTEMP(gt) == YES) {
+ call fstats (fd, F_FILENAME, Memc[fname], SZ_FNAME)
+ call close (fd)
+ iferr (call delete (Memc[fname]))
+ call erract (EA_WARN)
+ } else if (fd != GT_UFD(gt))
+ call close (fd)
+
+ if (level > 0) {
+ GT_FD(gt) = GT_SVFD(gt,level)
+ GT_FTEMP(gt) = GT_SVFTEMP(gt,level)
+ }
+ }
+
+ call mfree (gt, TY_STRUCT)
+ call sfree (sp)
+end
+
+
+# GT_ARGLIST -- Extract a paren and comma delimited argument list to be used
+# for substitution into a macro replacement string. Since the result will be
+# pushed back and rescanned, we do not have to perform macro substitution on
+# the argument list at this level.
+
+int procedure gt_arglist (gt, argbuf, maxch)
+
+pointer gt #I gettok descriptor
+char argbuf[maxch] #O receives parsed arguments
+int maxch #I max chars out
+
+int level, quote, nargs, op, ch, fd
+int getci()
+
+begin
+ fd = GT_FD(gt)
+
+ # Get lookahead char if we don't already have one.
+ ch = GT_NEXTCH(gt)
+ if (ch <= 0 || IS_WHITE(ch) || ch == '\n')
+ while (getci (fd, ch) != EOF)
+ if (!(IS_WHITE(ch) || ch == '\n'))
+ break
+
+ quote = 0
+ level = 1
+ nargs = 0
+ op = 1
+
+ if (ch == '(') {
+ while (getci (fd, ch) != EOF) {
+ if (ch == '"' || ch == '\'') {
+ if (quote == 0)
+ quote = ch
+ else if (quote == ch)
+ quote = 0
+
+ } else if (ch == '(' && quote == 0) {
+ level = level + 1
+ } else if (ch == ')' && quote == 0) {
+ level = level - 1
+ if (level <= 0) {
+ if (op > 1 && argbuf[op-1] != EOS)
+ nargs = nargs + 1
+ break
+ }
+
+ } else if (ch == ',' && level == 1 && quote == 0) {
+ ch = EOS
+ nargs = nargs + 1
+ } else if (ch == '\n') {
+ ch = ' '
+ } else if (ch == '\\' && quote == 0) {
+ ch = getci (fd, ch)
+ next
+ } else if (ch == '#' && quote == 0) {
+ while (getci (fd, ch) != EOF)
+ if (ch == '\n')
+ break
+ next
+ }
+
+ argbuf[op] = ch
+ op = min (maxch, op + 1)
+ }
+
+ GT_NEXTCH(gt) = NULL
+ }
+
+ argbuf[op] = EOS
+ return (nargs)
+end
diff --git a/pkg/proto/maskexpr/megeom.x b/pkg/proto/maskexpr/megeom.x
new file mode 100644
index 00000000..602493f8
--- /dev/null
+++ b/pkg/proto/maskexpr/megeom.x
@@ -0,0 +1,72 @@
+include <math.h>
+
+# ME_ELLGEOM -- Given the semi-major axis, ratio of semi-minor to semi-major
+# axes, and position angle, compute the parameters of the equation of the
+# ellipse, where the ellipse is defined as A * X ** 2 + B * x * y +
+# C * Y ** 2 - F = 0.
+
+procedure me_ellgeom (a, ratio, theta, aa, bb, cc, ff)
+
+real a #I the semi-major axis
+real ratio #I the ratio of semi-minor to semi-major axes
+real theta #I the position angle of the major axis
+real aa #O the coefficient of x ** 2
+real bb #O the coefficient of x * y
+real cc #O the coefficient of y ** 2
+real ff #O the constant term
+
+real cost, sint, costsq, sintsq
+real asq, bsq
+
+begin
+ # Get the angles.
+ cost = cos (DEGTORAD(theta))
+ sint = sin (DEGTORAD(theta))
+ costsq = cost ** 2
+ sintsq = sint ** 2
+
+ # Compute the parameters of the outer ellipse.
+ asq = a ** 2
+ bsq = (ratio * a) ** 2
+ aa = bsq * costsq + asq * sintsq
+ bb = 2.0 * (bsq - asq) * cost * sint
+ cc = asq * costsq + bsq * sintsq
+ ff = asq * bsq
+end
+
+
+# ME_RECTGEOM -- Construct a polygon representation of a rotated rectangle
+# givev the half-width of the long axis, the ratio of the half-width of the
+# short axis to the long axis, and the rotation angle.
+
+procedure me_rectgeom (hwidth, ratio, theta, xout, yout)
+
+real hwidth #I the half-width of the long axis of the rectangle
+real ratio #I the ratio of short to long axes of the rectangle
+real theta #I the rotation angle
+real xout[ARB] #O the x coordinates of the output vertices
+real yout[ARB] #O the y coordinates of the output vertices
+
+real cost, sint, x, y
+
+begin
+ cost = cos (DEGTORAD(theta))
+ sint = sin (DEGTORAD(theta))
+ x = hwidth
+ y = ratio * x
+ xout[1] = x * cost - y * sint
+ yout[1] = x * sint + y * cost
+ x = -x
+ y = y
+ xout[2] = x * cost - y * sint
+ yout[2] = x * sint + y * cost
+ x = x
+ y = -y
+ xout[3] = x * cost - y * sint
+ yout[3] = x * sint + y * cost
+ x = -x
+ y = y
+ xout[4] = x * cost - y * sint
+ yout[4] = x * sint + y * cost
+end
+
diff --git a/pkg/proto/maskexpr/megsym.x b/pkg/proto/maskexpr/megsym.x
new file mode 100644
index 00000000..44a810bc
--- /dev/null
+++ b/pkg/proto/maskexpr/megsym.x
@@ -0,0 +1,31 @@
+include <ctotok.h>
+include <ctype.h>
+include "gettok.h"
+
+
+# Expression database symbol.
+define LEN_SYM 2
+define SYM_TEXT Memi[$1]
+define SYM_NARGS Memi[$1+1]
+
+
+
+# ME_GSYM -- Get symbol routine for the gettok package.
+
+pointer procedure me_gsym (st, symname, nargs)
+
+pointer st #I symbol table
+char symname[ARB] #I symbol to be looked up
+int nargs #O number of macro arguments
+
+pointer sym
+pointer strefsbuf(), stfind()
+
+begin
+ sym = stfind (st, symname)
+ if (sym == NULL)
+ return (NULL)
+
+ nargs = SYM_NARGS(sym)
+ return (strefsbuf (st, SYM_TEXT(sym)))
+end
diff --git a/pkg/proto/maskexpr/memkmask.x b/pkg/proto/maskexpr/memkmask.x
new file mode 100644
index 00000000..f8af553c
--- /dev/null
+++ b/pkg/proto/maskexpr/memkmask.x
@@ -0,0 +1,839 @@
+include <mach.h>
+include <ctype.h>
+include <imhdr.h>
+include <imset.h>
+include <pmset.h>
+include <evvexpr.h>
+
+define DEF_LINELEN 8192
+
+define LEN_MSKEXPR 42
+define ME_PMIM Memi[$1] # the output mask image
+define ME_REFIM Memi[$1+1] # the reference image
+define ME_REFMSK Memi[$1+2] # the reference mask image
+define ME_REFDAT Memi[$1+3] # current reference image line
+define ME_REFTYPE Memi[$1+4] # the input pixel type
+define ME_REFPMDAT Memi[$1+5] # current mask image line
+define ME_PMV Meml[P2L($1+6+($2)-1)] # position in mask image
+define ME_REFV Meml[P2L($1+13+($2)-1)] # position in reference image
+define ME_REFPMV Meml[P2L($1+20+($2)-1)] # position in reference mask
+
+
+# ME_MKMASK -- Given an expression, a reference image descriptor, a reference
+# mask descriptor, the number of dimensions, size of each dimension, and depth
+# in bits create a mask image and return an imio pointer to the mask.
+
+pointer procedure me_mkmask (expr, mskname, refim, refmsk, ndim, axlen, depth)
+
+char expr[ARB] #I the input expression
+char mskname[ARB] #I the optional input mask name
+pointer refim #I the imio pointer to the reference image
+pointer refmsk #I the imio pointer to the reference mask
+int ndim #I the number of output mask dimensions
+long axlen[ARB] #I the size of the output mask
+int depth #I the pixel depth of the output mask
+
+pointer sp, tmpname, pm, pmim, me, obuf, oexpr
+pointer pm_create(), im_pmmap(), evvexpr(), immap()
+int i, npix, nlines, pmaxval
+int imstati()
+int imgnli(), imgnll(), imgnlr(), imgnld()
+int impnli(), impnls(), impnll()
+int locpr()
+extern me_getop(), me_fcn()
+
+begin
+ # Open the output mask and map it as a virtual image or a disk
+ # image depending on whether or not you wish to save the mask.
+ if (mskname[1] == EOS) {
+ call smark (sp)
+ call salloc (tmpname, SZ_FNAME, TY_CHAR)
+ call mktemp ("tmpmsk", Memc[tmpname], SZ_FNAME)
+ if (refim != NULL) {
+ pmim = im_pmmap (Memc[tmpname], NEW_COPY, refim)
+ } else if (refmsk != NULL) {
+ pmim = im_pmmap (Memc[tmpname], NEW_COPY, refmsk)
+ } else {
+ pmim = im_pmmap (Memc[tmpname], NEW_IMAGE, NULL)
+ IM_NDIM(pmim) = ndim
+ call amovl (axlen, IM_LEN(pmim,1), ndim)
+ }
+ call sfree (sp)
+ } else {
+ if (refim != NULL) {
+ pmim = immap (mskname, NEW_COPY, refim)
+ } else if (refmsk != NULL) {
+ pmim = immap (mskname, NEW_COPY, refmsk)
+ } else {
+ pmim = immap (mskname, NEW_IMAGE, 0)
+ IM_NDIM(pmim) = ndim
+ call amovl (axlen, IM_LEN(pmim,1), ndim)
+ }
+ }
+ IM_PIXTYPE(pmim) = TY_INT
+
+ # Initialize the mask.
+ pm = imstati (pmim, IM_PLDES)
+ call pl_close (pm)
+ pm = pm_create (IM_NDIM(pmim), IM_LEN(pmim,1), depth)
+ call imseti (pmim, IM_PLDES, pm)
+
+ # Determine the mask depth.
+ if (depth > 0) {
+ pmaxval = min (depth, PL_MAXDEPTH)
+ pmaxval = 2 ** pmaxval - 1
+ } else {
+ pmaxval = 2 ** PL_MAXDEPTH - 1
+ }
+
+ # Allocate space for the mask expression structure.
+ call calloc (me, LEN_MSKEXPR, TY_STRUCT)
+ ME_PMIM(me) = pmim
+ ME_REFIM(me) = refim
+ ME_REFMSK(me) = refmsk
+
+ # Determine the input image type.
+ if (refim != NULL) {
+ switch (IM_PIXTYPE(refim)) {
+ case TY_BOOL, TY_SHORT, TY_INT:
+ ME_REFTYPE(me) = TY_INT
+ case TY_LONG:
+ ME_REFTYPE(me) = TY_LONG
+ case TY_REAL:
+ ME_REFTYPE(me) = TY_REAL
+ case TY_DOUBLE:
+ ME_REFTYPE(me) = TY_DOUBLE
+ case TY_COMPLEX:
+ ME_REFTYPE(me) = TY_REAL
+ }
+ }
+
+ # Initalize the i/o pointers.
+ call amovkl (long(1), ME_PMV(me,1), IM_MAXDIM)
+ call amovkl (long(1), ME_REFV(me,1), IM_MAXDIM)
+ call amovkl (long(1), ME_REFPMV(me,1), IM_MAXDIM)
+
+ # Compute the total number of output image lines.
+ npix = IM_LEN(pmim,1)
+ nlines = 1
+ do i = 2, IM_NDIM(pmim)
+ nlines = nlines * IM_LEN(pmim, i)
+
+ # Loop over the mask output image lines which are by default always
+ # integer.
+ do i = 1, nlines {
+
+ # Get the correct reference image line.
+ if (refim != NULL) {
+ switch (ME_REFTYPE(me)) {
+ case TY_INT:
+ if (imgnli (refim, ME_REFDAT(me), ME_REFV(me,1)) == EOF)
+ call error (1, "Error reading reference image data")
+ case TY_LONG:
+ if (imgnll (refim, ME_REFDAT(me), ME_REFV(me,1)) == EOF)
+ call error (1, "Error reading reference image data")
+ case TY_REAL:
+ if (imgnlr (refim, ME_REFDAT(me), ME_REFV(me,1)) == EOF)
+ call error (1, "Error reading reference image data")
+ case TY_DOUBLE:
+ if (imgnld (refim, ME_REFDAT(me), ME_REFV(me,1)) == EOF)
+ call error (1, "Error reading reference image data")
+ case TY_COMPLEX:
+ if (imgnlr (refim, ME_REFDAT(me), ME_REFV(me,1)) == EOF)
+ call error (1, "Error reading reference image data")
+ }
+ }
+
+ # Get the correct reference mask line.
+ if (refmsk != NULL) {
+ if (imgnli (refmsk, ME_REFPMDAT(me), ME_REFPMV(me,1)) == EOF)
+ call error (1, "Error reading reference mask data")
+ }
+
+ # Evalute the expression.
+ oexpr = evvexpr (expr, locpr(me_getop), me, locpr(me_fcn), me, 0)
+ if (O_TYPE(oexpr) == ERR) {
+ call eprintf ("Error evaluting expression\n")
+ break
+ }
+
+ # Copy the evaluated expression to the image.
+ if (O_LEN(oexpr) == 0) {
+ switch (O_TYPE(oexpr)) {
+ case TY_BOOL:
+ if (impnli (pmim, obuf, ME_PMV(me,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropi (NULL, 1, MAX_INT, Memi[obuf], 1, pmaxval,
+ npix, PIX_CLR + PIX_VALUE(O_VALI(oexpr)))
+ case TY_SHORT:
+ if (impnls (pmim, obuf, ME_PMV(me,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixrops (NULL, 1, MAX_SHORT, Mems[obuf], 1,
+ pmaxval, npix, PIX_CLR + PIX_VALUE(O_VALS(oexpr)))
+ case TY_INT:
+ if (impnli (pmim, obuf, ME_PMV(me,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropi (NULL, 1, MAX_INT, Memi[obuf], 1,
+ pmaxval, npix, PIX_CLR + PIX_VALUE(O_VALI(oexpr)))
+ case TY_LONG:
+ if (impnll (pmim, obuf, ME_PMV(me,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropl (NULL, 1, MAX_LONG, Meml[obuf], 1,
+ pmaxval, npix, PIX_CLR + PIX_VALUE(O_VALL(oexpr)))
+ case TY_REAL:
+ call error (3, "Type real expressions are not supported")
+ case TY_DOUBLE:
+ call error (3, "Type double expressions are not supported")
+ default:
+ call error (3, "Unknown expression value type")
+ }
+
+ } else {
+ switch (O_TYPE(oexpr)) {
+ case TY_BOOL:
+ if (impnli (pmim, obuf, ME_PMV(me,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropi (Memi[O_VALP(oexpr)], 1, MAX_INT,
+ Memi[obuf], 1, pmaxval, npix, PIX_SRC)
+ case TY_SHORT:
+ if (impnls (pmim, obuf, ME_PMV(me,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixrops (Mems[O_VALP(oexpr)], 1, MAX_SHORT,
+ Mems[obuf], 1, pmaxval, npix, PIX_SRC)
+ case TY_INT:
+ if (impnli (pmim, obuf, ME_PMV(me,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropi (Memi[O_VALP(oexpr)], 1, MAX_INT,
+ Memi[obuf], 1, pmaxval, npix, PIX_SRC)
+ case TY_LONG:
+ if (impnll (pmim, obuf, ME_PMV(me,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropl (Meml[O_VALP(oexpr)], 1, MAX_LONG,
+ Meml[obuf], 1, pmaxval, npix, PIX_SRC)
+ case TY_REAL:
+ call error (3, "Type real expressions are not supported")
+ case TY_DOUBLE:
+ call error (3, "Type double expressions are not supported")
+ default:
+ call error (3, "Unknown expression value type")
+ }
+ }
+
+ call evvfree (oexpr)
+ }
+
+ # Cleanup.
+ call mfree (me, TY_STRUCT)
+
+ return (pmim)
+end
+
+
+# ME_GETOP -- Called by evvexpr to fetch an input image operand.
+
+procedure me_getop (me, opname, o)
+
+pointer me #I mskexpr descriptor
+char opname[ARB] #I operand name
+pointer o #I output operand to be filled in
+
+pointer sp, param, data, im
+int i, axis
+int imgftype(), btoi()
+double imgetd()
+int imgeti()
+bool imgetb()
+errchk malloc
+define err_ 91
+
+begin
+ call smark (sp)
+
+ # Reference image operand.
+ if ((opname[1] == 'i') && (opname[2] == EOS)) {
+
+ if (ME_REFIM(me) == NULL)
+ goto err_
+
+ O_TYPE(o) = ME_REFTYPE(me)
+ O_LEN(o) = IM_LEN(ME_REFIM(me), 1)
+ O_FLAGS(o) = 0
+ O_VALP(o) = ME_REFDAT(me)
+
+ call sfree (sp)
+ return
+
+ # Reference mask operand.
+ } else if ((opname[1] == 'm') && (opname[2] == EOS)) {
+
+ if (ME_REFMSK(me) == NULL)
+ goto err_
+
+ O_TYPE(o) = TY_INT
+ O_LEN(o) = IM_LEN(ME_REFMSK(me), 1)
+ O_FLAGS(o) = 0
+ O_VALP(o) = ME_REFPMDAT(me)
+
+ call sfree (sp)
+ return
+
+ # Reference image header parameter operand.
+ } else if ((opname[1] == 'i' || opname[1] == 'm') &&
+ (opname[2] == '.')) {
+
+ if (opname[1] == 'i')
+ im = ME_REFIM(me)
+ else
+ im = ME_REFMSK(me)
+ if (im == NULL)
+ goto err_
+
+ # Get the parameter value and set up operand struct.
+ call salloc (param, SZ_FNAME, TY_CHAR)
+ call strcpy (opname[3], Memc[param], SZ_FNAME)
+ iferr (O_TYPE(o) = imgftype (im, Memc[param]))
+ goto err_
+
+ switch (O_TYPE(o)) {
+ case TY_BOOL:
+ O_LEN(o) = 0
+ iferr (O_VALI(o) = btoi (imgetb (im, Memc[param])))
+ goto err_
+
+ case TY_CHAR:
+ O_LEN(o) = SZ_LINE
+ O_FLAGS(o) = O_FREEVAL
+ iferr {
+ call malloc (O_VALP(o), SZ_LINE, TY_CHAR)
+ call imgstr (im, Memc[param], O_VALC(o), SZ_LINE)
+ } then
+ goto err_
+
+ case TY_SHORT, TY_INT, TY_LONG:
+ iferr (O_VALI(o) = imgeti (im, Memc[param]))
+ goto err_
+
+ case TY_REAL, TY_DOUBLE:
+ O_TYPE(o) = TY_DOUBLE
+ iferr (O_VALD(o) = imgetd (im, Memc[param]))
+ goto err_
+
+ default:
+ goto err_
+ }
+
+ call sfree (sp)
+ return
+
+ # The current pixel coordinate [I,J,K,...]. The line coordinate
+ # is a special case since the image is computed a line at a time.
+ # If "I" is requested return a vector where v[i] = i. For J, K,
+ # etc. just return the scalar index value.
+
+ } else if (IS_UPPER(opname[1]) && opname[2] == EOS) {
+
+ axis = opname[1] - 'I' + 1
+ if (axis == 1) {
+ O_TYPE(o) = TY_INT
+ if (IM_LEN(ME_PMIM(me), 1) > 0)
+ O_LEN(o) = IM_LEN(ME_PMIM(me), 1)
+ else
+ O_LEN(o) = DEF_LINELEN
+ call malloc (data, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[data+i-1] = i
+ O_VALP(o) = data
+ O_FLAGS(o) = O_FREEVAL
+ } else {
+ O_TYPE(o) = TY_INT
+ if (IM_LEN(ME_PMIM(me), 1) > 0)
+ O_LEN(o) = IM_LEN(ME_PMIM(me), 1)
+ else
+ O_LEN(o) = DEF_LINELEN
+ call malloc (data, O_LEN(o), TY_INT)
+ if (axis < 1 || axis > IM_MAXDIM)
+ call amovki (1, Memi[data], O_LEN(o))
+ else
+ call amovki (ME_PMV(me,axis), Memi[data], O_LEN(o))
+ O_VALP(o) = data
+ O_FLAGS(o) = O_FREEVAL
+ }
+
+ call sfree (sp)
+ return
+ }
+
+err_
+ O_TYPE(o) = ERR
+ call sfree (sp)
+end
+
+
+# define the builtin functions
+
+define ME_FUNCS "|circle|ellipse|box|rectangle|polygon|cols|lines|\
+vector|pie|cannulus|eannulus|rannulus|pannulus|point|"
+
+define ME_CIRCLE 1
+define ME_ELLIPSE 2
+define ME_BOX 3
+define ME_RECTANGLE 4
+define ME_POLYGON 5
+define ME_COLS 6
+define ME_LINES 7
+define ME_VECTOR 8
+define ME_PIE 9
+define ME_CANNULUS 10
+define ME_EANNULUS 11
+define ME_RANNULUS 12
+define ME_PANNULUS 13
+define ME_POINT 14
+
+
+# ME_FCN -- Called by evvexpr to execute a mskexpr special function.
+
+procedure me_fcn (me, fcn, args, nargs, o)
+
+pointer me #I imexpr descriptor
+char fcn[ARB] #I function name
+pointer args[ARB] #I input arguments
+int nargs #I number of input arguments
+pointer o #I output operand to be filled in
+
+real width
+pointer sp, ufunc, rval1, rval2, orval1, orval2, ix, iy
+int i, ip, func, v_nargs, nver
+int strdic(), ctor()
+bool strne()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (ufunc, SZ_LINE, TY_CHAR)
+
+ # Get the function.
+ func = strdic (fcn, Memc[ufunc], SZ_LINE, ME_FUNCS)
+ if (func > 0 && strne (fcn, Memc[ufunc]))
+ func = 0
+
+ # Test the function.
+ if (func <= 0) {
+ O_TYPE(o) = ERR
+ call sfree (sp)
+ return
+ }
+
+ # Determine number of arguments. This is a separate case statement.
+ # in case we need to deal with a variable number of arguments
+ # function at a later point.
+ switch (func) {
+ case ME_POINT, ME_CIRCLE, ME_ELLIPSE, ME_BOX, ME_RECTANGLE, ME_POLYGON:
+ v_nargs = -1
+ case ME_CANNULUS, ME_EANNULUS, ME_RANNULUS, ME_PANNULUS:
+ v_nargs = -1
+ case ME_COLS, ME_LINES:
+ v_nargs = -1
+ case ME_VECTOR, ME_PIE:
+ v_nargs = -1
+ default:
+ v_nargs = 0
+ }
+
+ # Check the number of arguments.
+ if (v_nargs > 0 && nargs != v_nargs) {
+ O_TYPE(o) = ERR
+ call sfree (sp)
+ return
+ }
+ if (v_nargs < 0 && nargs < abs (v_nargs)) {
+ O_TYPE(o) = ERR
+ call sfree (sp)
+ return
+ }
+
+ if (func == ME_POLYGON && nargs < 6) {
+ O_TYPE(o) = ERR
+ call sfree (sp)
+ return
+ }
+
+ # Type convert the arguments appropriately. At the moment this is
+ # simple if we assume that all the required arguments are real.
+ call salloc (rval1, nargs, TY_REAL)
+ call salloc (rval2, nargs, TY_REAL)
+ do i = 1, nargs {
+ switch (O_TYPE(args[i])) {
+ case TY_CHAR:
+ ip = 1
+ if (ctor (O_VALC(args[i]), ip, Memr[rval1+i-1]) == 0)
+ Memr[rval1+i-1] = 0.
+ case TY_INT:
+ Memr[rval1+i-1] = O_VALI(args[i])
+ case TY_REAL:
+ Memr[rval1+i-1] = O_VALR(args[i])
+ case TY_DOUBLE:
+ Memr[rval1+i-1] = O_VALD(args[i])
+ }
+ }
+
+ # Evaluate the function. Worry about some duplication of code later.
+ switch (func) {
+
+ case ME_CIRCLE:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 5) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_circle (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4])
+ } else if (nargs == 3) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_circle (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_ELLIPSE:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 7) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_ellipse (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], Memr[rval1+6])
+ } else if (nargs == 5) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_ellipse (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_BOX:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 6) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_box (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5])
+ } else if (nargs == 4) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_box (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_RECTANGLE:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 7) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_rectangle (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], Memr[rval1+6])
+ } else if (nargs == 5) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_rectangle (Memi[ix], Memi[iy], Memi[O_VALP(o)],
+ O_LEN(o), Memr[rval1], Memr[rval1+1], Memr[rval1+2],
+ Memr[rval1+3], Memr[rval1+4])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_POLYGON:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs < 6) {
+ O_TYPE(o) = ERR
+ } else if (O_LEN(args[1]) > 0 && O_LEN(args[2]) > 0) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ nver = (nargs - 2) / 2
+ do i = 1, nver
+ #Memr[rval2+i-1] = Memr[rval1+2*i]
+ Memr[rval2+i-1] = Memr[rval1+2*i+1]
+ do i = 1, nver
+ #Memr[rval1+i-1] = Memr[rval1+2*i+1]
+ Memr[rval1+i-1] = Memr[rval1+2*i]
+ call me_polygon (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1], Memr[rval2], nver)
+ } else {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ nver = nargs / 2
+ do i = 1, nver
+ Memr[rval2+i-1] = Memr[rval1+2*i-1]
+ do i = 1, nver
+ Memr[rval1+i-1] = Memr[rval1+2*i-2]
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_polygon (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval2], nver)
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ }
+
+ case ME_COLS:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 2) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_cols (Memi[O_VALP(args[1])], Memi[O_VALP(o)], O_LEN(o),
+ O_VALC(args[2]))
+ } else if (nargs == 1) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_cols (Memi[ix], Memi[O_VALP(o)], O_LEN(o),
+ O_VALC(args[1]))
+ call mfree (ix, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_LINES:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 2) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_lines (Memi[O_VALP(args[1])], Memi[O_VALP(o)], O_LEN(o),
+ O_VALC(args[2]))
+ } else if (nargs == 1) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call amovki (ME_PMV(me,2), Memi[ix], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_lines (Memi[ix], Memi[O_VALP(o)], O_LEN(o),
+ O_VALC(args[1]))
+ call mfree (ix, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_VECTOR:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 7) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_vector (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], Memr[rval1+6])
+ } else if (nargs == 5) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_vector (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_PIE:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 6) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_pie (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], IM_LEN(ME_PMIM(me),1),
+ IM_LEN(ME_PMIM(me),2))
+ } else if (nargs == 4) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_pie (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3],
+ IM_LEN(ME_PMIM(me),1), IM_LEN(ME_PMIM(me),2))
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_CANNULUS:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 6) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_cannulus (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5])
+ } else if (nargs == 4) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_cannulus (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_EANNULUS:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 8) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_eannulus (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], Memr[rval1+6], Memr[rval1+7])
+ } else if (nargs == 6) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_eannulus (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_RANNULUS:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 8) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_rannulus (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], Memr[rval1+6], Memr[rval1+7])
+ } else if (nargs == 6) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_rannulus (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case ME_PANNULUS:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs < 7) {
+ O_TYPE(o) = ERR
+ } else if (O_LEN(args[1]) > 0 && O_LEN(args[2]) > 0) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ width = Memr[rval1+2]
+ nver = (nargs - 3) / 2
+ do i = 1, nver
+ #Memr[rval2+i-1] = Memr[rval1+2*i+1]
+ Memr[rval2+i-1] = Memr[rval1+2*i+2]
+ do i = 1, nver
+ #Memr[rval1+i-1] = Memr[rval1+2*i+2]
+ Memr[rval1+i-1] = Memr[rval1+2*i+1]
+ call salloc (orval1, nver, TY_REAL)
+ call salloc (orval2, nver, TY_REAL)
+ call me_pyexpand (Memr[rval1], Memr[rval2], Memr[orval1],
+ Memr[orval2], nver, width)
+ call me_apolygon (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1], Memr[rval2],
+ Memr[orval1], Memr[orval2], nver)
+ } else {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ width = Memr[rval1]
+ nver = (nargs - 1) / 2
+ do i = 1, nver
+ Memr[rval2+i-1] = Memr[rval1+2*i]
+ do i = 1, nver
+ Memr[rval1+i-1] = Memr[rval1+2*i-1]
+ call salloc (orval1, nver, TY_REAL)
+ call salloc (orval2, nver, TY_REAL)
+ call me_pyexpand (Memr[rval1], Memr[rval2], Memr[orval1],
+ Memr[orval2], nver, width)
+ call me_apolygon (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval2], Memr[orval1], Memr[orval2], nver)
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ }
+
+ case ME_POINT:
+ O_LEN(o) = IM_LEN(ME_PMIM(me),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 4) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_point (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3])
+ } else if (nargs == 2) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (ME_PMV(me,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_point (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+ default:
+ O_TYPE(o) = ERR
+ }
+
+ call sfree (sp)
+end
+
diff --git a/pkg/proto/maskexpr/meregfuncs.x b/pkg/proto/maskexpr/meregfuncs.x
new file mode 100644
index 00000000..467bd1d0
--- /dev/null
+++ b/pkg/proto/maskexpr/meregfuncs.x
@@ -0,0 +1,1449 @@
+include <mach.h>
+include <ctype.h>
+include <math.h>
+
+
+# ME_POINT -- Compute which pixels are equal to a point.
+
+procedure me_point (ix, iy, stat, npts, x1, y1)
+
+int ix[ARB] #I the integer x coordinates
+int iy[ARB] #I the integer y coordinates
+int stat[ARB] #O the integer status array containing YES or NO
+int npts #I the number of points
+real x1, y1 #I the coordinates of the point
+
+int i
+
+begin
+ do i = 1, npts {
+ if (ix[i] == nint(x1) && iy[i] == nint(y1))
+ stat[i] = YES
+ else
+ stat[i] = NO
+ }
+end
+
+
+# ME_CIRCLE -- Compute which pixels are within or on a circle.
+
+procedure me_circle (ix, iy, stat, npts, xc, yc, r)
+
+int ix[ARB] #I the integer x coordinates
+int iy[ARB] #I the integer y coordinates
+int stat[ARB] #O the integer status array containing YES or NO
+int npts #I the number of points
+real xc, yc #I the center of the circle
+real r #I the radius of the circle
+
+real r2, rdist
+int i
+
+begin
+ r2 = r * r
+ do i = 1, npts {
+ rdist = (ix[i] - xc) ** 2 + (iy[i] - yc) ** 2
+ if (rdist <= r2)
+ stat[i] = YES
+ else
+ stat[i] = NO
+ }
+end
+
+
+# ME_CANNULUS -- Compute which pixels are within or on a circular annulus
+# boundary.
+
+procedure me_cannulus (ix, iy, stat, npts, xc, yc, r1, r2)
+
+int ix[ARB] #I the integer x coordinates
+int iy[ARB] #I the integer y coordinates
+int stat[ARB] #O the integer status array containing YES or NO
+int npts #I the number of points
+real xc, yc #I the center of the circle
+real r1, r2 #I the radius of the circular annulus
+
+real r12, r22, rdist
+int i
+
+begin
+ r12 = r1 * r1
+ r22 = r2 * r2
+ do i = 1, npts {
+ rdist = (ix[i] - xc) ** 2 + (iy[i] - yc) ** 2
+ if (rdist >= r12 && rdist <= r22)
+ stat[i] = YES
+ else
+ stat[i] = NO
+ }
+end
+
+
+# ME_ELLIPSE -- Compute which pixels lie within or on an ellipse.
+
+procedure me_ellipse (ix, iy, stat, npts, xc, yc, a, ratio, theta)
+
+int ix[ARB] #I the integer x coordinates
+int iy[ARB] #I the integer y coordinates
+int stat[ARB] #O the integer status array (YES/NO)
+int npts #I the number of points
+real xc, yc #I the center of the ellipse
+real a #I the semi-major axis of the ellipse
+real ratio #I the semi-minor / semi-minor axis
+real theta #I the position angle of the ellipse
+
+real asq, bsq, cost, sint, costsq, sintsq, rdist
+real dx, dy, aa, bb, cc, rr
+int i
+
+begin
+ asq = a * a
+ bsq = (ratio * a) * (ratio * a)
+ cost = cos (DEGTORAD(theta))
+ sint = sin (DEGTORAD(theta))
+ costsq = cost * cost
+ sintsq = sint * sint
+ aa = bsq * costsq + asq * sintsq
+ bb = 2.0 * (bsq - asq) * cost * sint
+ cc = asq * costsq + bsq * sintsq
+ rr = asq * bsq
+
+ do i = 1, npts {
+ dx = (ix[i] - xc)
+ dy = (iy[i] - yc)
+ rdist = aa * dx * dx + bb * dx * dy + cc * dy * dy
+ if (rdist <= rr)
+ stat[i] = YES
+ else
+ stat[i] = NO
+ }
+end
+
+
+# ME_EANNULUS -- Compute which pixels lie within or on an elliptical annular
+# boundary.
+
+procedure me_eannulus (ix, iy, stat, npts, xc, yc, a1, a2, ratio, theta)
+
+int ix[ARB] #I the integer x coordinates
+int iy[ARB] #I the integer y coordinates
+int stat[ARB] #O the integer status array (YES/NO)
+int npts #I the number of points
+real xc, yc #I the center of the ellipse
+real a1, a2 #I the semi-major axis of the i/o ellipse
+real ratio #I the semi-minor / semi-major axis of ellipse
+real theta #I the position angle of the ellipse
+
+real a1sq, b1sq, aa1, bb1, cc1, rr1, rdist1
+real a2sq, b2sq, aa2, bb2, cc2, rr2, rdist2
+real dx, dy, cost, sint, costsq, sintsq
+int i
+
+begin
+ # First ellipse.
+ a1sq = a1 * a1
+ b1sq = (ratio * a1) ** 2
+ cost = cos (DEGTORAD(theta))
+ sint = sin (DEGTORAD(theta))
+ costsq = cost * cost
+ sintsq = sint * sint
+ aa1 = b1sq * costsq + a1sq * sintsq
+ bb1 = 2.0 * (b1sq - a1sq) * cost * sint
+ cc1 = a1sq * costsq + b1sq * sintsq
+ rr1 = a1sq * b1sq
+
+ # Second ellipse.
+ a2sq = a2 * a2
+ b2sq = (ratio * a2) ** 2
+ aa2 = b2sq * costsq + a2sq * sintsq
+ bb2 = 2.0 * (b2sq - a2sq) * cost * sint
+ cc2 = a2sq * costsq + b2sq * sintsq
+ rr2 = a2sq * b2sq
+
+ # Elliptical annulus.
+ do i = 1, npts {
+ dx = (ix[i] - xc)
+ dy = (iy[i] - yc)
+ rdist1 = aa1 * dx * dx + bb1 * dx * dy + cc1 * dy * dy
+ rdist2 = aa2 * dx * dx + bb2 * dx * dy + cc2 * dy * dy
+ if (rdist1 >= rr1 && rdist2 <= rr2)
+ stat[i] = YES
+ else
+ stat[i] = NO
+ }
+end
+
+
+# ME_RECTANGLE -- Compute which pixels lie within or on a rectangle.
+
+procedure me_rectangle (ix, iy, stat, npts, xc, yc, a, ratio, theta)
+
+int ix[ARB] #I the integer x coordinates
+int iy[ARB] #I the integer y coordinates
+int stat[ARB] #O the integer status array (YES/NO)
+int npts #I the number of points
+real xc, yc #I the center of the rectangle
+real a #I the semi-major axis width of the rectangle
+real ratio #I the semi-minor axis / semi-major axis
+real theta #I the position angle of the rectangle
+
+real cost, sint, x, y
+real xver[4], yver[4]
+
+begin
+ # Compute the corners of the equivalent polygon.
+ cost = cos (DEGTORAD(theta))
+ sint = sin (DEGTORAD(theta))
+ x = a
+ y = ratio * a
+ xver[1] = xc + x * cost - y * sint
+ yver[1] = yc + x * sint + y * cost
+ x = -x
+ y = y
+ xver[2] = xc + x * cost - y * sint
+ yver[2] = yc + x * sint + y * cost
+ x = x
+ y = -y
+ xver[3] = xc + x * cost - y * sint
+ yver[3] = yc + x * sint + y * cost
+ x = -x
+ y = y
+ xver[4] = xc + x * cost - y * sint
+ yver[4] = yc + x * sint + y * cost
+
+ # Call the polygon routine.
+ call me_polygon (ix, iy, stat, npts, xver, yver, 4)
+end
+
+
+# ME_RANNULUS -- Compute which pixels lie within or on a rectangular annulus.
+
+procedure me_rannulus (ix, iy, stat, npts, xc, yc, r1, r2, ratio, theta)
+
+int ix[ARB] #I the integer x coordinates
+int iy[ARB] #I the integer y coordinates
+int stat[ARB] #O the integer status array (YES/NO)
+int npts #I the number of points
+real xc, yc #I the center of the rectangle
+real r1, r2 #I the semi-major axis width of the rectangle
+real ratio #I the semi-minor / semi-major axis ratio
+real theta #I the position angle of the rectangle
+
+real cost, sint, x, y, xver1[4], yver1[4], xver2[4], yver2[4]
+
+begin
+ # Compute the corners of the equivalent polygon inner and outer
+ # polygons.
+ cost = cos (DEGTORAD(theta))
+ sint = sin (DEGTORAD(theta))
+
+ # The corners of the inner polygon.
+ x = r1
+ y = ratio * r1
+ xver1[1] = xc + x * cost - y * sint
+ yver1[1] = yc + x * sint + y * cost
+ x = -x
+ y = y
+ xver1[2] = xc + x * cost - y * sint
+ yver1[2] = yc + x * sint + y * cost
+ x = x
+ y = -y
+ xver1[3] = xc + x * cost - y * sint
+ yver1[3] = yc + x * sint + y * cost
+ x = -x
+ y = y
+ xver1[4] = xc + x * cost - y * sint
+ yver1[4] = yc + x * sint + y * cost
+
+ # The corners of the outer polygon.
+ x = r2
+ y = ratio * r2
+ xver2[1] = xc + x * cost - y * sint
+ yver2[1] = yc + x * sint + y * cost
+ x = -x
+ y = y
+ xver2[2] = xc + x * cost - y * sint
+ yver2[2] = yc + x * sint + y * cost
+ x = x
+ y = -y
+ xver2[3] = xc + x * cost - y * sint
+ yver2[3] = yc + x * sint + y * cost
+ x = -x
+ y = y
+ xver2[4] = xc + x * cost - y * sint
+ yver2[4] = yc + x * sint + y * cost
+
+ # Write a routine to determine which pixels are inside the polygon
+ # defined by 2 sets of vertices.
+ call me_apolygon (ix, iy, stat, npts, xver1, yver1, xver2, yver2, 4)
+end
+
+
+# ME_BOX -- Compute which pixels lie within or on a box.
+
+procedure me_box (ix, iy, stat, npts, x1, y1, x2, y2)
+
+int ix[ARB] #I the integer x coordinates
+int iy[ARB] #I the integer y coordinates
+int stat[ARB] #O the integer status array (YES/NO)
+int npts #I the number of points
+real x1, y1 #I first box corner
+real x2, y2 #I first box corner
+
+real xmin, xmax, ymin, ymax
+int i
+
+begin
+ xmin = min (x1, x2)
+ xmax = max (x1, x2)
+ ymin = min (y1, y2)
+ ymax = max (y1, y2)
+
+ do i = 1, npts {
+ if (ix[i] >= xmin && ix[i] <= xmax &&
+ iy[i] >= ymin && iy[i] <= ymax)
+ stat[i] = YES
+ else
+ stat[i] = NO
+ }
+end
+
+
+# ME_POLYGON -- Determine which points lie in or on a specified polygon.
+
+procedure me_polygon (ix, iy, stat, npts, xver, yver, nver)
+
+int ix[ARB] #I the x image pixel coordinates
+int iy[ARB] #I the y image pixel coordinates
+int stat[ARB] #O the output status array
+int npts #I the number of image pixel coordinates
+real xver[ARB] #I the x polygon vertices coordinates
+real yver[ARB] #I the y polygon vertices coordinates
+int nver #I the number of polygon coordinates
+
+real lx, ld
+pointer sp, txver, tyver, work1, work2, xintr
+int i, j, ixmin, ixmax, nintr
+int me_pyclip()
+
+begin
+ call smark (sp)
+ call salloc (txver, nver + 1, TY_REAL)
+ call salloc (tyver, nver + 1, TY_REAL)
+ call salloc (work1, nver + 1, TY_REAL)
+ call salloc (work2, nver + 1, TY_REAL)
+ call salloc (xintr, nver + 1, TY_REAL)
+
+ # Close the polygon.
+ call amovr (xver, Memr[txver], nver)
+ call amovr (yver, Memr[tyver], nver)
+ Memr[txver+nver] = xver[1]
+ Memr[tyver+nver] = yver[1]
+
+ # Loop over the points.
+ call alimi (ix, npts, ixmin, ixmax)
+ lx = ixmax - ixmin + 1
+ do i = 1, npts {
+
+ # Compute the intersection points of the line segment which
+ # spans an image line with the polygon. Sort the line segments.
+ ld = iy[i]
+ if (i == 1) {
+ nintr = me_pyclip (Memr[txver], Memr[tyver], Memr[work1],
+ Memr[work2], Memr[xintr], nver + 1, lx, ld)
+ call asrtr (Memr[xintr], Memr[xintr], nintr)
+ } else if (iy[i] != iy[i-1]) {
+ nintr = me_pyclip (Memr[txver], Memr[tyver], Memr[work1],
+ Memr[work2], Memr[xintr], nver + 1, lx, ld)
+ call asrtr (Memr[xintr], Memr[xintr], nintr)
+ }
+
+ # Are the intersection points in range ?
+ if (nintr <= 0)
+ stat[i] = NO
+ else {
+ stat[i] = NO
+ do j = 1, nintr, 2 {
+ if (ix[i] >= Memr[xintr+j-1] && ix[i] <= Memr[xintr+j])
+ stat[i] = YES
+ }
+ }
+
+ }
+
+ call sfree (sp)
+end
+
+
+# ME_APOLYGON -- Determine which points lie in or on a specified polygonal
+# annulus.
+
+procedure me_apolygon (ix, iy, stat, npts, ixver, iyver, oxver, oyver, nver)
+
+int ix[ARB] #I the x image pixel coordinates
+int iy[ARB] #I the y image pixel coordinates
+int stat[ARB] #O the output status array
+int npts #I the number of image pixel coordinates
+real ixver[ARB] #I the x polygon vertices coordinates
+real iyver[ARB] #I the y polygon vertices coordinates
+real oxver[ARB] #I the x polygon vertices coordinates
+real oyver[ARB] #I the y polygon vertices coordinates
+int nver #I the number of polygon coordinates
+
+real lx, ld
+pointer sp, tixver, tiyver, toxver, toyver, work1, work2, ixintr, oxintr
+int i, j, jj, ixmin, ixmax, inintr, onintr, ibegin, iend
+int me_pyclip()
+
+begin
+ call smark (sp)
+ call salloc (tixver, nver + 1, TY_REAL)
+ call salloc (tiyver, nver + 1, TY_REAL)
+ call salloc (toxver, nver + 1, TY_REAL)
+ call salloc (toyver, nver + 1, TY_REAL)
+ call salloc (work1, nver + 1, TY_REAL)
+ call salloc (work2, nver + 1, TY_REAL)
+ call salloc (ixintr, nver + 1, TY_REAL)
+ call salloc (oxintr, nver + 1, TY_REAL)
+
+ # Close the polygons.
+ call amovr (ixver, Memr[tixver], nver)
+ call amovr (iyver, Memr[tiyver], nver)
+ Memr[tixver+nver] = ixver[1]
+ Memr[tiyver+nver] = iyver[1]
+ call amovr (oxver, Memr[toxver], nver)
+ call amovr (oyver, Memr[toyver], nver)
+ Memr[toxver+nver] = oxver[1]
+ Memr[toyver+nver] = oyver[1]
+
+ # Loop over the points.
+ call alimi (ix, npts, ixmin, ixmax)
+ lx = ixmax - ixmin + 1
+ do i = 1, npts {
+
+ stat[i] = NO
+
+ # Compute the intersection points of the line segment with the
+ # outer polygon.
+ ld = iy[i]
+ if (i == 1) {
+ onintr = me_pyclip (Memr[toxver], Memr[toyver], Memr[work1],
+ Memr[work2], Memr[oxintr], nver + 1, lx, ld)
+ call asrtr (Memr[oxintr], Memr[oxintr], onintr)
+ } else if (iy[i] != iy[i-1]) {
+ onintr = me_pyclip (Memr[toxver], Memr[toyver], Memr[work1],
+ Memr[work2], Memr[oxintr], nver + 1, lx, ld)
+ call asrtr (Memr[oxintr], Memr[oxintr], onintr)
+ }
+ if (onintr <= 0)
+ next
+
+ # Compute the intersection points of the line segment with the
+ # inner polygon.
+ if (i == 1) {
+ inintr = me_pyclip (Memr[tixver], Memr[tiyver], Memr[work1],
+ Memr[work2], Memr[ixintr], nver + 1, lx, ld)
+ call asrtr (Memr[ixintr], Memr[ixintr], inintr)
+ } else if (iy[i] != iy[i-1]) {
+ inintr = me_pyclip (Memr[tixver], Memr[tiyver], Memr[work1],
+ Memr[work2], Memr[ixintr], nver + 1, lx, ld)
+ call asrtr (Memr[ixintr], Memr[ixintr], inintr)
+ }
+
+ # Are the intersection points in range ?
+ if (inintr <= 0) {
+ do j = 1, onintr, 2 {
+ if (ix[i] >= Memr[oxintr+j-1] && ix[i] <= Memr[oxintr+j]) {
+ stat[i] = YES
+ break
+ }
+ }
+ } else {
+ do j = 1, onintr, 2 {
+ do jj = 1, inintr, 2 {
+ if ((Memr[ixintr+jj-1] > Memr[oxintr+j-1]) &&
+ (Memr[ixintr+jj-1] < Memr[oxintr+j])) {
+ ibegin = jj
+ break
+ }
+ }
+ do jj = inintr, 1, -2 {
+ if ((Memr[ixintr+jj-1] > Memr[oxintr+j-1]) &&
+ (Memr[ixintr+jj-1] < Memr[oxintr+j])) {
+ iend = jj
+ break
+ }
+ }
+ if ((ix[i] >= Memr[oxintr+j-1]) &&
+ (ix[i] <= Memr[ixintr+ibegin-1])) {
+ stat[i] = YES
+ } else if ((ix[i] >= Memr[ixintr+iend-1]) &&
+ (ix[i] <= Memr[oxintr+j])) {
+ stat[i] = YES
+ } else {
+ do jj = ibegin + 1, iend - 1, 2 {
+ if ((ix[i] >= Memr[ixintr+jj-1]) &&
+ (ix[i] <= Memr[ixintr+jj])) {
+ stat[i] = YES
+ break
+ }
+ }
+ }
+
+ }
+ }
+
+ }
+
+ call sfree (sp)
+end
+
+
+define MAX_NRANGES 100
+
+# ME_COLS -- Determine which pixels are in the specified column ranges.
+
+procedure me_cols (ix, stat, npts, rangstr)
+
+int ix[ARB] #I the x image pixel coordinates
+int stat[ARB] #O the output status array
+int npts #I the number of image pixel coordinates
+char rangstr[ARB] #I the input range specification string
+
+pointer sp, ranges
+int index, nvals
+int me_decode_ranges(), me_next_number()
+
+begin
+ # Allocate space for storing the ranges.
+ call smark (sp)
+ call salloc (ranges, 3 * MAX_NRANGES + 1, TY_INT)
+
+ # Decode the ranges string. If there was an error set up the ranges
+ # so as to include everything.
+ if (me_decode_ranges (rangstr, Memi[ranges], MAX_NRANGES,
+ nvals) == ERR) {
+ if (me_decode_ranges ("-", Memi[ranges], MAX_NRANGES, nvals) != ERR)
+ ;
+ }
+
+ # Set the status array.
+ call amovki (NO, stat, npts)
+ index = 0
+ while (me_next_number (Memi[ranges], index) != EOF)
+ stat[index] = YES
+
+ call sfree (sp)
+end
+
+
+# ME_LINES -- Determine which pixels are in the specified column ranges.
+
+procedure me_lines (ix, stat, npts, rangstr)
+
+int ix[ARB] #I the x image pixel coordinates
+int stat[ARB] #O the output status array
+int npts #I the number of image pixel coordinates
+char rangstr[ARB] #I the input range specification string
+
+pointer sp, ranges
+int i, lastix, nvals
+int me_decode_ranges()
+bool me_is_in_range()
+
+begin
+ # Allocate space for storing the ranges.
+ call smark (sp)
+ call salloc (ranges, 3 * MAX_NRANGES + 1, TY_INT)
+
+ # Decode the ranges string. If there was an error set up the ranges
+ # so as to include everything.
+ if (me_decode_ranges (rangstr, Memi[ranges], MAX_NRANGES,
+ nvals) == ERR) {
+ if (me_decode_ranges ("-", Memi[ranges], MAX_NRANGES, nvals) != ERR)
+ ;
+ }
+
+ # Set the line numbers.
+ call amovki (NO, stat, npts)
+ lastix = 0
+ do i = 1, npts {
+ if (ix[i] == lastix) {
+ stat[i] = YES
+ } else if (me_is_in_range (Memi[ranges], ix[i])) {
+ lastix = ix[i]
+ stat[i] = YES
+ }
+ }
+
+ call sfree (sp)
+end
+
+
+# ME_VECTOR -- Determine which pixels are on the specified line.
+
+procedure me_vector (ix, iy, stat, npts, x1, y1, x2, y2, width)
+
+int ix[ARB] #I the x image pixel coordinates
+int iy[ARB] #I the y image pixel coordinates
+int stat[ARB] #O the output status array
+int npts #I the number of image pixel coordinates
+real x1, y1 #I coordinates of the first point
+real x2, y2 #I coordinates of the first point
+real width #I the vector width
+
+real x, y, xc, yc, theta, cost, sint
+real xver[4], yver[4]
+
+begin
+ # Compute the corners of the equivalent polygon.
+ xc = (x2 + x1) / 2.0
+ yc = (y2 + y1) / 2.0
+ x = sqrt ((x2 - x1) ** 2 + (y2 - y1) ** 2) / 2.0
+ y = width / 2.0
+ theta = atan2 (y2 - y1, x2 - x1)
+ cost = cos (theta)
+ sint = sin (theta)
+ xver[1] = xc + x * cost - y * sint
+ yver[1] = yc + x * sint + y * cost
+ x = -x
+ y = y
+ xver[2] = xc + x * cost - y * sint
+ yver[2] = yc + x * sint + y * cost
+ x = x
+ y = -y
+ xver[3] = xc + x * cost - y * sint
+ yver[3] = yc + x * sint + y * cost
+ x = -x
+ y = y
+ xver[4] = xc + x * cost - y * sint
+ yver[4] = yc + x * sint + y * cost
+
+ # Call the polygon routine.
+ call me_polygon (ix, iy, stat, npts, xver, yver, 4)
+end
+
+
+define SMALL_NUMBER 1.0e-24
+
+# ME_PIE -- Determine which pixels are inside a pie shaped wedge that
+# intersects the image boundaries.
+
+procedure me_pie (ix, iy, stat, npts, xc, yc, angle1, angle2, width, height)
+
+int ix[ARB] #I the x pixel coordinates
+int iy[ARB] #I the y pixel coordinates
+int stat[ARB] #O the output status array
+int npts #I the number of data points
+real xc, yc #I the center of the wedge
+real angle1, angle2 #I the wedge angles
+int width, height #I the image mask width and height
+
+real sweep, x2, y2, vx[7], vy[7]
+int count, intrcpt1, intrcpt2
+int me_pie_intercept(), me_corner_vertex()
+
+begin
+ # Set the first vertex
+ vx[1] = xc
+ vy[1] = yc
+ sweep = angle2 - angle1
+
+ # If the sweep is too small to be noticed don't bother.
+ if (abs (sweep) < SMALL_NUMBER) {
+ call amovki (NO, stat, npts)
+ return
+ }
+ if (sweep < 0.0)
+ sweep = sweep + 360.0
+
+ # Get the second vertext by computing the intersection of the
+ # first ray with the image boundaries.
+ intrcpt1 = me_pie_intercept (width, height, xc, yc, angle1,
+ vx[2], vy[2])
+
+ # Compute the second intercept.
+ intrcpt2 = me_pie_intercept (width, height, xc, yc, angle2, x2, y2)
+
+ # If angles intercept same side and slice is between them, no corners
+ # else, mark corners until reaching side with second angle intercept.
+ count = 3
+ if ((intrcpt1 != intrcpt2) || (sweep > 180.0)) {
+ repeat {
+ intrcpt1 = me_corner_vertex (intrcpt1, width, height, vx[count],
+ vy[count])
+ count = count + 1
+ } until (intrcpt1 == intrcpt2)
+ }
+
+ # Set last vertex.
+ vx[count] = x2
+ vy[count] = y2
+
+ # Fill in the polygon
+ call me_polygon (ix, iy, stat, npts, vx, vy, count)
+end
+
+
+# ME_PIE_INTERCEPT -- Determine which side is intercepted by a vertex (given
+# center and angle) and set edge intercept point and return index of side.
+
+int procedure me_pie_intercept (width, height, xcen, ycen, angle, xcept, ycept)
+
+int width, height #I the dimensions of the image field
+real xcen, ycen #I the base pivot point of the ray
+real angle #I the angle of ray
+real xcept, ycept #I coordinates of intercept with edge of image
+
+real angl, slope
+
+begin
+ # Put angles in normal range.
+ angl = angle
+ while (angl < 0.0)
+ angl = angl + 360.0
+ while (angl >= 360.0)
+ angl = angl - 360.0
+
+ # Check for a horizontal angle.
+ if (abs (angl) < SMALL_NUMBER) {
+ #xcept = 0
+ xcept = width + 1
+ ycept = ycen
+ #return (2)
+ return (4)
+ }
+ if (abs (angl - 180.0) < SMALL_NUMBER) {
+ #xcept = width + 1
+ xcept = 0
+ ycept = ycen
+ #return (4)
+ return (2)
+ }
+
+# # Convert to a Cartesian angle
+# angl = angl + 90.0
+# if (angl >= 360.0)
+# angl = angl - 360.0
+
+ # Check for vertical angle.
+ if (angl < 180.0) {
+ ycept = height + 1
+ # rule out vertical line
+ if (abs(angl - 90.0) < SMALL_NUMBER) {
+ x_cept = xcen
+ return (1)
+ }
+ } else {
+ ycept = 0.0
+ # rule out vertical line
+ if (abs(angl - 270.0) < SMALL_NUMBER) {
+ xcept = xcen
+ return (3)
+ }
+ }
+
+ # Convert to radians.
+ angl = (angl / 180.0) * PI
+
+ # Calculate slope.
+ slope = tan (angl)
+
+ # Calculate intercept with designated y edge.
+ xcept = xcen + ((ycept - ycen) / slope)
+ if (xcept < 0) {
+ ycept = (ycen - (xcen * slope))
+ xcept = 0.0
+ return (2)
+ } else if (xcept > (width + 1)) {
+ ycept = (ycen + ((width + 1 - xcen) * slope))
+ xcept = width + 1
+ return (4)
+ } else {
+ if (ycept < height)
+ return (3)
+ else
+ return (1)
+ }
+end
+
+
+# ME_CORNER_VERTEX -- Set points just beyond corner to mark the corner in a
+# polygon. Note: 1=top, 2=left, 3=bottom, 4=right, corner is between current
+# and next advance index to next side and also return its value.
+
+int procedure me_corner_vertex (index, width, height, x, y)
+
+int index #I code of side before corner
+int width, height #I dimensions of image field
+real x, y #O coords of corner
+
+begin
+ # Set the corner coordinates.
+ switch (index) {
+ case 1:
+ x = 0.0
+ y = height + 1
+ case 2:
+ x = 0.0
+ y = 0.0
+ case 3:
+ x = width + 1
+ y = 0.0
+ case 4:
+ x = width + 1
+ y = height + 1
+ default:
+ ; #call error (1, "index error in mark_corner")
+ }
+
+ # Set the corner index.
+ index = index + 1
+ if (index > 4)
+ index = 1
+
+ return (index)
+end
+
+
+# ME_PYEXPAND -- Expand a polygon given a list of vertices and an expansion
+# factor in pixels.
+
+procedure me_pyexpand (xin, yin, xout, yout, nver, width)
+
+real xin[ARB] #I the x coordinates of the input vertices
+real yin[ARB] #I the y coordinates of the input vertices
+real xout[ARB] #O the x coordinates of the output vertices
+real yout[ARB] #O the y coordinates of the output vertices
+int nver #I the number of vertices
+real width #I the width of the expansion region
+
+real xcen, ycen, m1, b1, m2, b2, xp1, yp1, xp2, yp2
+int i
+real asumr()
+
+begin
+ # Find the center of gravity of the polygon.
+ xcen = asumr (xin, nver) / nver
+ ycen = asumr (yin, nver) / nver
+
+ do i = 1, nver {
+
+ # Compute the equations of the line segments parallel to the
+ # line seqments composing a single vertex.
+ if (i == 1) {
+ call me_psegment (xcen, ycen, xin[nver], yin[nver], xin[1],
+ yin[1], width, m1, b1, xp1, yp1)
+ call me_psegment (xcen, ycen, xin[1], yin[1], xin[2], yin[2],
+ width, m2, b2, xp2, yp2)
+ } else if (i == nver) {
+ call me_psegment (xcen, ycen, xin[nver-1], yin[nver-1],
+ xin[nver], yin[nver], width, m1, b1, xp1, yp1)
+ call me_psegment (xcen, ycen, xin[nver], yin[nver], xin[1],
+ yin[1], width, m2, b2, xp2, yp2)
+ } else {
+ call me_psegment (xcen, ycen, xin[i-1], yin[i-1], xin[i],
+ yin[i], width, m1, b1, xp1, yp1)
+ call me_psegment (xcen, ycen, xin[i], yin[i], xin[i+1],
+ yin[i+1], width, m2, b2, xp2, yp2)
+ }
+
+ # The new vertex is the intersection of the two new line
+ # segments.
+ if (m1 == m2) {
+ xout[i] = xp2
+ yout[i] = yp2
+ } else if (IS_INDEFR(m1)) {
+ xout[i] = xp1
+ yout[i] = m2 * xp1 + b2
+ } else if (IS_INDEFR(m2)) {
+ xout[i] = xp2
+ yout[i] = m1 * xp2 + b1
+ } else {
+ xout[i] = (b2 - b1) / (m1 - m2)
+ yout[i] = (m2 * b1 - m1 * b2) / (m2 - m1)
+ }
+ }
+end
+
+
+# ME_PSEGMENT -- Construct a line segment parallel to an existing line segment
+# but a specified distance from it in a direction away from a fixed reference
+# point.
+
+procedure me_psegment (xcen, ycen, xb, yb, xe, ye, width, m, b, xp, yp)
+
+real xcen, ycen #I the position of the reference point
+real xb, yb #I the starting coordinates of the line segment
+real xe, ye #I the ending coordinates of the line segment
+real width #I the distance of new line segment from old
+real m #O the slope of the new line segment
+real b #O the intercept of the new line segment
+real xp, yp #O the coordinates of a points on new line
+
+real x1, y1, x2, y2, d1, d2
+
+begin
+ # Compute the slope of the line segment.
+ m = (xe - xb)
+ if (m == 0.0)
+ m = INDEFR
+ else
+ m = (ye - yb) / m
+
+ # Construct the perpendicular to the line segement and locate two
+ # points which are equidistant from the line seqment
+ if (IS_INDEFR(m)) {
+ x1 = xb - width
+ y1 = yb
+ x2 = xb + width
+ y2 = yb
+ } else if (m == 0.0) {
+ x1 = xb
+ y1 = yb - width
+ x2 = xb
+ y2 = yb + width
+ } else {
+ x1 = xb - sqrt ((m * width) ** 2 / (m ** 2 + 1))
+ y1 = yb - (x1 - xb) / m
+ x2 = xb + sqrt ((m * width) ** 2 / (m ** 2 + 1))
+ y2 = yb - (x2 - xb) / m
+ }
+
+ # Choose the point farthest away from the reference point.
+ d1 = (x1 - xcen) ** 2 + (y1 - ycen) ** 2
+ d2 = (x2 - xcen) ** 2 + (y2 - ycen) ** 2
+ if (d1 <= d2) {
+ xp = x2
+ yp = y2
+ } else {
+ xp = x1
+ yp = y1
+ }
+
+ # Compute the intercept.
+ if (IS_INDEFR(m))
+ b = INDEFR
+ else
+ b = yp - m * xp
+end
+
+
+# ME_PYCLIP -- Compute the intersection of an image line with a polygon defined
+# by a list of vertices. The output is a list of ranges stored in the array
+# xranges. Two additional work arrays xintr and slope are required for
+# the computation.
+
+int procedure me_pyclip (xver, yver, xintr, slope, xranges, nver, lx, ld)
+
+real xver[ARB] #I the x vertex coords
+real yver[ARB] #I the y vertex coords
+real xintr[ARB] #O the array of x intersection points
+real slope[ARB] #O the array of y slopes at intersection points
+real xranges[ARB] #O the x line segments
+int nver #I the number of vertices
+real lx, ld #I the equation of the image line
+
+real u1, u2, u1u2, dx, dy, dd, xa, wa
+int i, j, nintr, nplus, nzero, nneg, imin, imax, nadd
+bool collinear
+
+begin
+ # Initialize.
+ collinear = false
+ nplus = 0
+ nzero = 0
+ nneg = 0
+ nintr = 0
+
+ # Compute the intersection points of the image line and the polygon.
+ u1 = lx * (- yver[1] + ld)
+ do i = 2, nver {
+
+ u2 = lx * (- yver[i] + ld)
+ u1u2 = u1 * u2
+
+ # Does the polygon side intersect the image line ?
+ if (u1u2 <= 0.0) {
+
+
+ # Compute the x intersection coordinate if the point of
+ # intersection is not a vertex.
+
+ 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 = lx * (dx * ld - dd)
+ wa = dy * lx
+ nintr = nintr + 1
+ xranges[nintr] = xa / wa
+ slope[nintr] = -dy
+ if (slope[nintr] < 0.0)
+ nneg = nneg + 1
+ else if (slope[nintr] > 0.0)
+ nplus = nplus + 1
+ else
+ nzero = nzero + 1
+ collinear = false
+
+ # For each collinear line segment add two intersection
+ # points. Remove interior collinear intersection points.
+
+ } else if (u1 == 0.0 && u2 == 0.0) {
+
+ if (! collinear) {
+ nintr = nintr + 1
+ xranges[nintr] = xver[i-1]
+ if (i == 2)
+ slope[nintr] = yver[1] - yver[nver-1]
+ else
+ slope[nintr] = yver[i-1] - yver[i-2]
+ if (slope[nintr] < 0.0)
+ nneg = nneg + 1
+ else if (slope[nintr] > 0.0)
+ nplus = nplus + 1
+ else
+ nzero = nzero + 1
+ nintr = nintr + 1
+ xranges[nintr] = xver[i]
+ slope[nintr] = 0.0
+ nzero = nzero + 1
+ } else {
+ xranges[nintr] = xver[i]
+ slope[nintr] = 0.0
+ nzero = nzero + 1
+ }
+ collinear = true
+
+ # If the intersection point is a vertex add it to the
+ # list if it is not collinear with the next point. Add
+ # another point to the list if the vertex is at the
+ # apex of an acute angle.
+
+ } else if (u1 != 0.0) {
+
+ if (i == nver) {
+ dx = (xver[2] - xver[nver])
+ dy = (yver[2] - yver[nver])
+ dd = dy * (yver[nver-1] - yver[nver])
+ } else {
+ dx = (xver[i+1] - xver[i])
+ dy = (yver[i+1] - yver[i])
+ dd = dy * (yver[i-1] - yver[i])
+ }
+
+ # Test whether the point is collinear with the point
+ # ahead. If it is not include the intersection point.
+
+ if (dy != 0.0) {
+ nintr = nintr + 1
+ xranges[nintr] = xver[i]
+ slope[nintr] = yver[i] - yver[i-1]
+ if (slope[nintr] < 0.0)
+ nneg = nneg + 1
+ else if (slope[nintr] > 0.0)
+ nplus = nplus + 1
+ else
+ nzero = nzero + 1
+ }
+
+ # If the intersection point is an isolated vertex add
+ # another point to the list.
+
+ if (dd > 0.0) {
+ nintr = nintr + 1
+ xranges[nintr] = xver[i]
+ slope[nintr] = dy
+ if (slope[nintr] < 0.0)
+ nneg = nneg + 1
+ else if (slope[nintr] > 0.0)
+ nplus = nplus + 1
+ else
+ nzero = nzero + 1
+ }
+
+ collinear = false
+
+ } else
+ collinear = false
+ } else
+ collinear = false
+
+ u1 = u2
+ }
+
+ # Join up any split collinear line segments.
+ if (collinear && (slope[1] == 0.0)) {
+ xranges[1] = xranges[nintr-1]
+ slope[1] = slope[nintr-1]
+ nintr = nintr - 2
+ nzero = nzero - 2
+ }
+
+ # Return the number of intersection points if there are no interior
+ # collinear line segments.
+ if (nzero == 0 || nplus == 0 || nneg == 0)
+ return (nintr)
+
+ # Find the minimum and maximum intersection points.
+ call me_alimr (xranges, nintr, u1, u2, imin, imax)
+
+ # Check for vertices at the ends of the ranges.
+
+ u1 = xranges[min(imin,imax)] - xranges[1]
+ u2 = xranges[nintr] - xranges[max(imin,imax)]
+
+ # Vertices were traversed in order of increasing x.
+ if ((u1 >= 0.0 && u2 > 0.0) || (u1 > 0.0 && u2 >= 0.0) ||
+ (u1 == u2 && imax > imin)) {
+ do i = imax + 1, nintr {
+ if (xranges[i] != xranges[i-1])
+ break
+ imax = i
+ }
+ do i = imin - 1, 1, -1 {
+ if (xranges[i] != xranges[i+1])
+ break
+ imin = i
+ }
+ }
+
+ # Vertices were traversed in order of decreasing x.
+ if ((u1 <= 0.0 && u2 < 0.0) || (u1 < 0.0 && u2 <= 0.0) ||
+ (u1 == u2 && imax < imin)) {
+ do i = imin + 1, nintr {
+ if (xranges[i] != xranges[i-1])
+ break
+ imin = i
+ }
+ do i = imax - 1, 1, -1 {
+ if (xranges[i] != xranges[i+1])
+ break
+ imax = i
+ }
+ }
+
+ # Reorder the x ranges and slopes if necessary.
+ if ((imax < imin) && ! (imin == nintr && imax == 1)) {
+ call amovr (xranges, xintr, nintr)
+ do i = 1, imax
+ xranges[nintr-imax+i] = xintr[i]
+ do i = imin, nintr
+ xranges[i-imax] = xintr[i]
+ call amovr (slope, xintr, nintr)
+ do i = 1, imax
+ slope[nintr-imax+i] = xintr[i]
+ do i = imin, nintr
+ slope[i-imax] = xintr[i]
+ } else if ((imin < imax) && ! (imin == 1 && imax == nintr)) {
+ call amovr (xranges, xintr, nintr)
+ do i = 1, imin
+ xranges[nintr-imin+i] = xintr[i]
+ do i = imax, nintr
+ xranges[i-imin] = xintr[i]
+ call amovr (slope, xintr, nintr)
+ do i = 1, imin
+ slope[nintr-imin+i] = xintr[i]
+ do i = imax, nintr
+ slope[i-imin] = xintr[i]
+ }
+
+ # Add any extra intersection points that are required to deal with
+ # the collinear line segments.
+
+ nadd = 0
+ for (i = 1; i <= nintr-2; ) {
+ if (slope[i] * slope[i+2] > 0.0) {
+ i = i + 2
+ } else {
+ nadd = nadd + 1
+ xranges[nintr+nadd] = xranges[i+1]
+ for (j = i + 3; j <= nintr; j = j + 1) {
+ if (slope[i] * slope[j] > 0)
+ break
+ nadd = nadd + 1
+ xranges[nintr+nadd] = xranges[j-1]
+ }
+ i = j
+ }
+ }
+
+ return (nintr + nadd)
+end
+
+
+# ME_ALIMR -- Compute the maximum and minimum data values and indices of a
+# 1D array.
+
+procedure me_alimr (data, npts, mindat, maxdat, imin, imax)
+
+real data[npts] #I the input data array
+int npts #I the number of points
+real mindat, maxdat #O the minimum and maximum data values
+int imin, imax #O the indices of the minimum and maximum data values
+
+int i
+
+begin
+ imin = 1
+ imax = 1
+ mindat = data[1]
+ maxdat = data[1]
+
+ do i = 2, npts {
+ if (data[i] > maxdat) {
+ imax = i
+ maxdat = data[i]
+ }
+ if (data[i] < mindat) {
+ imin = i
+ mindat = data[i]
+ }
+ }
+end
+
+
+define FIRST 1 # Default starting range
+define LAST MAX_INT # Default ending range
+define STEP 1 # Default step
+define EOLIST 0 # End of list
+
+# ME_DECODE_RANGES -- Parse a string containing a list of integer numbers or
+# ranges, delimited by either spaces or commas. Return as output a list
+# of ranges defining a list of numbers, and the count of list numbers.
+# Range limits must be positive nonnegative integers. ERR is returned as
+# the function value if a conversion error occurs. The list of ranges is
+# delimited by EOLIST.
+
+int procedure me_decode_ranges (range_string, ranges, max_ranges, nvalues)
+
+char range_string[ARB] #I range string to be decoded
+int ranges[3, max_ranges] #O output range array
+int max_ranges #I maximum number of ranges
+int nvalues #O the number of values in the ranges
+
+int ip, nrange, first, last, step, ctoi()
+
+begin
+ ip = 1
+ nvalues = 0
+
+ do nrange = 1, max_ranges - 1 {
+ # Defaults to all nonnegative integers
+ first = FIRST
+ last = LAST
+ step = STEP
+
+ # Skip delimiters
+ while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+ ip = ip + 1
+
+ # Get first limit.
+ # Must be a number, '-', 'x', or EOS. If not return ERR.
+ if (range_string[ip] == EOS) { # end of list
+ if (nrange == 1) {
+ # Null string defaults
+ ranges[1, 1] = first
+ ranges[2, 1] = last
+ ranges[3, 1] = step
+ ranges[1, 2] = EOLIST
+ nvalues = MAX_INT
+ return (OK)
+ } else {
+ ranges[1, nrange] = EOLIST
+ return (OK)
+ }
+ } else if (range_string[ip] == '-')
+ ;
+ else if (range_string[ip] == 'x')
+ ;
+ else if (IS_DIGIT(range_string[ip])) { # ,n..
+ if (ctoi (range_string, ip, first) == 0)
+ return (ERR)
+ } else
+ return (ERR)
+
+ # Skip delimiters
+ while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+ ip = ip + 1
+
+ # Get last limit
+ # Must be '-', or 'x' otherwise last = first.
+ if (range_string[ip] == 'x')
+ ;
+ else if (range_string[ip] == '-') {
+ ip = ip + 1
+ while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+ ip = ip + 1
+ if (range_string[ip] == EOS)
+ ;
+ else if (IS_DIGIT(range_string[ip])) {
+ if (ctoi (range_string, ip, last) == 0)
+ return (ERR)
+ } else if (range_string[ip] == 'x')
+ ;
+ else
+ return (ERR)
+ } else
+ last = first
+
+ # Skip delimiters
+ while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+ ip = ip + 1
+
+ # Get step.
+ # Must be 'x' or assume default step.
+ if (range_string[ip] == 'x') {
+ ip = ip + 1
+ while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+ ip = ip + 1
+ if (range_string[ip] == EOS)
+ ;
+ else if (IS_DIGIT(range_string[ip])) {
+ if (ctoi (range_string, ip, step) == 0)
+ ;
+ if (step == 0)
+ return (ERR)
+ } else if (range_string[ip] == '-')
+ ;
+ else
+ return (ERR)
+ }
+
+ # Output the range triple.
+ ranges[1, nrange] = first
+ ranges[2, nrange] = last
+ ranges[3, nrange] = step
+ nvalues = nvalues + abs (last-first) / step + 1
+ }
+
+ return (ERR) # ran out of space
+end
+
+
+# ME_NEXT_NUMBER -- Given a list of ranges and the current file number,
+# find and return the next file number. Selection is done in such a way
+# that list numbers are always returned in monotonically increasing order,
+# regardless of the order in which the ranges are given. Duplicate entries
+# are ignored. EOF is returned at the end of the list.
+
+int procedure me_next_number (ranges, number)
+
+int ranges[ARB] #I the range array
+int number #U both input and output parameter
+
+int ip, first, last, step, next_number, remainder
+
+begin
+ # If number+1 is anywhere in the list, that is the next number,
+ # otherwise the next number is the smallest number in the list which
+ # is greater than number+1.
+
+ number = number + 1
+ next_number = MAX_INT
+
+ for (ip=1; ranges[ip] != EOLIST; ip=ip+3) {
+ first = min (ranges[ip], ranges[ip+1])
+ last = max (ranges[ip], ranges[ip+1])
+ step = ranges[ip+2]
+ if (step == 0)
+ call error (1, "Step size of zero in range list")
+ if (number >= first && number <= last) {
+ remainder = mod (number - first, step)
+ if (remainder == 0)
+ return (number)
+ if (number - remainder + step <= last)
+ next_number = number - remainder + step
+ } else if (first > number)
+ next_number = min (next_number, first)
+ }
+
+ if (next_number == MAX_INT)
+ return (EOF)
+ else {
+ number = next_number
+ return (number)
+ }
+end
+
+
+# ME_PREVIOUS_NUMBER -- Given a list of ranges and the current file number,
+# find and return the previous file number. Selection is done in such a way
+# that list numbers are always returned in monotonically decreasing order,
+# regardless of the order in which the ranges are given. Duplicate entries
+# are ignored. EOF is returned at the end of the list.
+
+int procedure me_previous_number (ranges, number)
+
+int ranges[ARB] # Range array
+int number # Both input and output parameter
+
+int ip, first, last, step, next_number, remainder
+
+begin
+ # If number-1 is anywhere in the list, that is the previous number,
+ # otherwise the previous number is the largest number in the list which
+ # is less than number-1.
+
+ number = number - 1
+ next_number = 0
+
+ for (ip=1; ranges[ip] != EOLIST; ip=ip+3) {
+ first = min (ranges[ip], ranges[ip+1])
+ last = max (ranges[ip], ranges[ip+1])
+ step = ranges[ip+2]
+ if (step == 0)
+ call error (1, "Step size of zero in range list")
+ if (number >= first && number <= last) {
+ remainder = mod (number - first, step)
+ if (remainder == 0)
+ return (number)
+ if (number - remainder >= first)
+ next_number = number - remainder
+ } else if (last < number) {
+ remainder = mod (last - first, step)
+ if (remainder == 0)
+ next_number = max (next_number, last)
+ else if (last - remainder >= first)
+ next_number = max (next_number, last - remainder)
+ }
+ }
+
+ if (next_number == 0)
+ return (EOF)
+ else {
+ number = next_number
+ return (number)
+ }
+end
+
+
+# ME_IS_IN_RANGE -- Test number to see if it is in range. If the number is
+# INDEFI then it is mapped to the maximum integer.
+
+bool procedure me_is_in_range (ranges, number)
+
+int ranges[ARB] # range array
+int number # number to be tested against ranges
+
+int ip, first, last, step, num
+
+begin
+ if (IS_INDEFI (number))
+ num = MAX_INT
+ else
+ num = number
+
+ for (ip = 1; ranges[ip] != EOLIST; ip = ip+3) {
+ first = min (ranges[ip], ranges[ip+1])
+ last = max (ranges[ip], ranges[ip+1])
+ step = ranges[ip+2]
+ if (num >= first && num <= last)
+ if (mod (num - first, step) == 0)
+ return (true)
+ }
+
+ return (false)
+end
diff --git a/pkg/proto/maskexpr/meregmask.x b/pkg/proto/maskexpr/meregmask.x
new file mode 100644
index 00000000..45db9079
--- /dev/null
+++ b/pkg/proto/maskexpr/meregmask.x
@@ -0,0 +1,753 @@
+include <mach.h>
+include <ctype.h>
+include <imhdr.h>
+include <imset.h>
+include <pmset.h>
+include <evvexpr.h>
+
+define DEF_LINELEN 8192
+
+define LEN_RGEXPR 25
+define RG_PMIM Memi[$1] # the mask image
+define RG_PMIBUF Memi[$1+1] # the mask input data
+define RG_IPMV Meml[P2L($1+2+($2)-1)] # input position in mask image
+define RG_OPMV Meml[P2L($1+9+($2)-1)] # output position in mask image
+
+
+# ME_RGMASK -- Given a region expression, a condition equals true expression,
+# a condition equals false expression, and an existing pixel mask imio
+# descriptor of dimensions, size of each dimension, and depth in bits create
+# a mask image and return an imio pointer to the mask.
+
+int procedure me_rgmask (rexpr, texpr, fexpr, pmim)
+
+char rexpr[ARB] #I the boolean region expression
+char texpr[ARB] #I the condition equals true expression
+char fexpr[ARB] #I the condition equals true expression
+pointer pmim #I the pixel mask imio descriptor
+
+pointer sp, rg, oexpr, expr, obuf
+int i, npix, nlines, depth, pmaxval, stat
+
+pointer evvexpr()
+int imstati(), locpr(), pm_stati()
+int imgnli(), impnli(), impnls(), impnll()
+extern rg_getop(), rg_fcn()
+
+begin
+ # Allocate some work space.
+ call smark (sp)
+ call salloc (expr, 3 * SZ_LINE, TY_CHAR)
+
+ # Allocate space for the mask expression structure.
+ call calloc (rg, LEN_RGEXPR, TY_STRUCT)
+ RG_PMIM(rg) = pmim
+
+ # Initalize the i/o pointers.
+ call amovkl (long(1), RG_OPMV(rg,1), IM_MAXDIM)
+ call amovkl (long(1), RG_IPMV(rg,1), IM_MAXDIM)
+
+ # Create the conditional expression to be evaluated.
+ call sprintf (Memc[expr], 3 * SZ_LINE, "(%s) ? %s : %s")
+ call pargstr (rexpr)
+ call pargstr (texpr)
+ call pargstr (fexpr)
+
+ # Compute the total number of output image lines.
+ npix = IM_LEN(pmim,1)
+ nlines = 1
+ do i = 2, IM_NDIM(pmim)
+ nlines = nlines * IM_LEN(pmim, i)
+ depth = INDEFI
+
+ # Loop over the mask output image lines which are by default always
+ # integer.
+ stat = OK
+ do i = 1, nlines {
+
+ # Get the input mask lines.
+ if (imgnli (pmim, RG_PMIBUF(rg), RG_IPMV(rg,1)) == EOF)
+ call error (2, "Error reading input mask data")
+
+ # Determine the depth of the mask.
+ if (IS_INDEFI(depth)) {
+ depth = pm_stati (imstati (pmim, IM_PLDES), P_DEPTH)
+ if (depth > 0) {
+ pmaxval = min (depth, PL_MAXDEPTH)
+ pmaxval = 2 ** depth - 1
+ } else
+ pmaxval = 2 ** PL_MAXDEPTH - 1
+ }
+
+ # Evalute the expression.
+ oexpr = evvexpr (Memc[expr], locpr(rg_getop), rg, locpr(rg_fcn),
+ rg, 0)
+ if (O_TYPE(oexpr) == ERR) {
+ call eprintf ("Error evaluting expression\n")
+ stat = ERR
+ break
+ }
+
+ # Copy the evaluated expression to the image.
+ if (O_LEN(oexpr) == 0) {
+ switch (O_TYPE(oexpr)) {
+ case TY_BOOL:
+ if (impnli (pmim, obuf, RG_OPMV(rg,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropi (NULL, 1, MAX_INT, Memi[obuf], 1, pmaxval,
+ npix, PIX_CLR + PIX_VALUE(O_VALI(oexpr)))
+ case TY_SHORT:
+ if (impnls (pmim, obuf, RG_OPMV(rg,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixrops (NULL, 1, MAX_SHORT, Mems[obuf], 1,
+ pmaxval, npix, PIX_CLR + PIX_VALUE(O_VALS(oexpr)))
+ case TY_INT:
+ if (impnli (pmim, obuf, RG_OPMV(rg,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropi (NULL, 1, MAX_INT, Memi[obuf], 1,
+ pmaxval, npix, PIX_CLR + PIX_VALUE(O_VALI(oexpr)))
+ case TY_LONG:
+ if (impnll (pmim, obuf, RG_OPMV(rg,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropl (NULL, 1, MAX_LONG, Meml[obuf], 1,
+ pmaxval, npix, PIX_CLR + PIX_VALUE(O_VALL(oexpr)))
+ case TY_REAL:
+ call error (3, "Type real expressions are not supported")
+ case TY_DOUBLE:
+ call error (3, "Type double expressions are not supported")
+ default:
+ call error (3, "Unknown expression value type")
+ }
+
+ } else {
+ switch (O_TYPE(oexpr)) {
+ case TY_BOOL:
+ if (impnli (pmim, obuf, RG_OPMV(rg,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropi (Memi[O_VALP(oexpr)], 1, MAX_INT,
+ Memi[obuf], 1, pmaxval, npix, PIX_SRC)
+ case TY_SHORT:
+ if (impnls (pmim, obuf, RG_OPMV(rg,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixrops (Mems[O_VALP(oexpr)], 1, MAX_SHORT,
+ Mems[obuf], 1, pmaxval, npix, PIX_SRC)
+ case TY_INT:
+ if (impnli (pmim, obuf, RG_OPMV(rg,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropi (Memi[O_VALP(oexpr)], 1, MAX_INT,
+ Memi[obuf], 1, pmaxval, npix, PIX_SRC)
+ case TY_LONG:
+ if (impnll (pmim, obuf, RG_OPMV(rg,1)) == EOF)
+ call error (2, "Error writing output mask data")
+ call pl_pixropl (Meml[O_VALP(oexpr)], 1, MAX_LONG,
+ Meml[obuf], 1, pmaxval, npix, PIX_SRC)
+ case TY_REAL:
+ call error (3, "Type real expressions are not supported")
+ case TY_DOUBLE:
+ call error (3, "Type double expressions are not supported")
+ default:
+ call error (3, "Unknown expression value type")
+ }
+ }
+
+ call evvfree (oexpr)
+ }
+
+ # Cleanup.
+ call mfree (rg, TY_STRUCT)
+
+ call sfree (sp)
+
+ return (stat)
+end
+
+
+# RG_GETOP -- Called by evvexpr to fetch an input image operand.
+
+procedure rg_getop (rg, opname, o)
+
+pointer rg #I mskexpr descriptor
+char opname[ARB] #I operand name
+pointer o #I output operand to be filled in
+
+pointer sp, param, data, im
+int i, axis
+int imgftype(), btoi()
+double imgetd()
+int imgeti()
+bool imgetb()
+errchk malloc
+define err_ 91
+
+begin
+ call smark (sp)
+
+ # Pixel image operand.
+ if ((opname[1] == 'p') && (opname[2] == EOS)) {
+
+ if (RG_PMIM(rg) == NULL)
+ goto err_
+
+ O_TYPE(o) = TY_INT
+ O_LEN(o) = IM_LEN(RG_PMIM(rg), 1)
+ O_FLAGS(o) = 0
+ O_VALP(o) = RG_PMIBUF(rg)
+
+ call sfree (sp)
+ return
+
+ # Reference image header parameter operand.
+ } else if ((opname[1] == 'p') && (opname[2] == '.')) {
+
+ im = RG_PMIM(rg)
+ if (im == NULL)
+ goto err_
+
+ # Get the parameter value and set up operand struct.
+ call salloc (param, SZ_FNAME, TY_CHAR)
+ call strcpy (opname[3], Memc[param], SZ_FNAME)
+ iferr (O_TYPE(o) = imgftype (im, Memc[param]))
+ goto err_
+
+ switch (O_TYPE(o)) {
+
+ case TY_BOOL:
+ O_LEN(o) = 0
+ iferr (O_VALI(o) = btoi (imgetb (im, Memc[param])))
+ goto err_
+
+ case TY_CHAR:
+ O_LEN(o) = SZ_LINE
+ O_FLAGS(o) = O_FREEVAL
+ iferr {
+ call malloc (O_VALP(o), SZ_LINE, TY_CHAR)
+ call imgstr (im, Memc[param], O_VALC(o), SZ_LINE)
+ } then
+ goto err_
+
+ case TY_SHORT, TY_INT, TY_LONG:
+ iferr (O_VALI(o) = imgeti (im, Memc[param]))
+ goto err_
+
+ case TY_REAL, TY_DOUBLE:
+ O_TYPE(o) = TY_DOUBLE
+ iferr (O_VALD(o) = imgetd (im, Memc[param]))
+ goto err_
+
+ default:
+ goto err_
+ }
+
+ call sfree (sp)
+ return
+
+ # The current pixel coordinate [I,J,K,...]. The line coordinate
+ # is a special case since the image is computed a line at a time.
+ # If "I" is requested return a vector where v[i] = i. For J, K,
+ # etc. just return the scalar index value.
+
+ } else if (IS_UPPER(opname[1]) && opname[2] == EOS) {
+
+ axis = opname[1] - 'I' + 1
+ if (axis == 1) {
+ O_TYPE(o) = TY_INT
+ if (IM_LEN(RG_PMIM(rg), 1) > 0)
+ O_LEN(o) = IM_LEN(RG_PMIM(rg), 1)
+ else
+ O_LEN(o) = DEF_LINELEN
+ call malloc (data, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[data+i-1] = i
+ O_VALP(o) = data
+ O_FLAGS(o) = O_FREEVAL
+ } else {
+ O_TYPE(o) = TY_INT
+ if (IM_LEN(RG_PMIM(rg), 1) > 0)
+ O_LEN(o) = IM_LEN(RG_PMIM(rg), 1)
+ else
+ O_LEN(o) = DEF_LINELEN
+ call malloc (data, O_LEN(o), TY_INT)
+ if (axis < 1 || axis > IM_MAXDIM)
+ call amovki (1, Memi[data], O_LEN(o))
+ else
+ call amovki (RG_OPMV(rg,axis), Memi[data], O_LEN(o))
+ O_VALP(o) = data
+ O_FLAGS(o) = O_FREEVAL
+ }
+
+ call sfree (sp)
+ return
+ }
+
+err_
+ O_TYPE(o) = ERR
+ call sfree (sp)
+end
+
+
+# define the builtin functions
+
+define RG_FUNCS "|circle|ellipse|box|rectangle|polygon|cols|lines|\
+vector|pie|cannulus|eannulus|rannulus|pannulus|point|"
+
+define RG_CIRCLE 1
+define RG_ELLIPSE 2
+define RG_BOX 3
+define RG_RECTANGLE 4
+define RG_POLYGON 5
+define RG_COLS 6
+define RG_LINES 7
+define RG_VECTOR 8
+define RG_PIE 9
+define RG_CANNULUS 10
+define RG_EANNULUS 11
+define RG_RANNULUS 12
+define RG_PANNULUS 13
+define RG_POINT 14
+
+
+# RG_FCN -- Called by evvexpr to execute a mskexpr special function.
+
+procedure rg_fcn (rg, fcn, args, nargs, o)
+
+pointer rg #I imexpr descriptor
+char fcn[ARB] #I function name
+pointer args[ARB] #I input arguments
+int nargs #I number of input arguments
+pointer o #I output operand to be filled in
+
+real width
+pointer sp, ufunc, rval1, rval2, orval1, orval2, ix, iy
+int i, ip, func, v_nargs, nver
+int strdic(), ctor()
+bool strne()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (ufunc, SZ_LINE, TY_CHAR)
+
+ # Get the function.
+ func = strdic (fcn, Memc[ufunc], SZ_LINE, RG_FUNCS)
+ if (func > 0 && strne (fcn, Memc[ufunc]))
+ func = 0
+
+ # Test the function.
+ if (func <= 0) {
+ O_TYPE(o) = ERR
+ call sfree (sp)
+ return
+ }
+
+ # Determine number of arguments. This is a separate case statement.
+ # in case we need to deal with a variable number of arguments
+ # function at a later point.
+ switch (func) {
+ case RG_POINT, RG_CIRCLE, RG_ELLIPSE, RG_BOX, RG_RECTANGLE, RG_POLYGON:
+ v_nargs = -1
+ case RG_CANNULUS, RG_EANNULUS, RG_RANNULUS, RG_PANNULUS:
+ v_nargs = -1
+ case RG_COLS, RG_LINES:
+ v_nargs = -1
+ case RG_VECTOR, RG_PIE:
+ v_nargs = -1
+ default:
+ v_nargs = 0
+ }
+
+ # Check the number of arguments.
+ if (v_nargs > 0 && nargs != v_nargs) {
+ O_TYPE(o) = ERR
+ call sfree (sp)
+ return
+ }
+ if (v_nargs < 0 && nargs < abs (v_nargs)) {
+ O_TYPE(o) = ERR
+ call sfree (sp)
+ return
+ }
+
+ if (func == RG_POLYGON && nargs < 6) {
+ O_TYPE(o) = ERR
+ call sfree (sp)
+ return
+ }
+
+ # Type convert the arguments appropriately. At the moment this is
+ # simple if we assume that all the required arguments are real.
+ call salloc (rval1, nargs, TY_REAL)
+ call salloc (rval2, nargs, TY_REAL)
+ do i = 1, nargs {
+ switch (O_TYPE(args[i])) {
+ case TY_CHAR:
+ ip = 1
+ if (ctor (O_VALC(args[i]), ip, Memr[rval1+i-1]) == 0)
+ Memr[rval1+i-1] = 0.
+ case TY_INT:
+ Memr[rval1+i-1] = O_VALI(args[i])
+ case TY_REAL:
+ Memr[rval1+i-1] = O_VALR(args[i])
+ case TY_DOUBLE:
+ Memr[rval1+i-1] = O_VALD(args[i])
+ }
+ }
+
+ # Evaluate the function. Worry about some duplication of code later.
+ switch (func) {
+
+ case RG_CIRCLE:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 5) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_circle (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4])
+ } else if (nargs == 3) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_circle (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_ELLIPSE:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 7) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_ellipse (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], Memr[rval1+6])
+ } else if (nargs == 5) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_ellipse (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_BOX:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 6) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_box (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5])
+ } else if (nargs == 4) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_box (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_RECTANGLE:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 7) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_rectangle (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], Memr[rval1+6])
+ } else if (nargs == 5) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_rectangle (Memi[ix], Memi[iy], Memi[O_VALP(o)],
+ O_LEN(o), Memr[rval1], Memr[rval1+1], Memr[rval1+2],
+ Memr[rval1+3], Memr[rval1+4])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_POLYGON:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs < 6) {
+ O_TYPE(o) = ERR
+ } else if (O_LEN(args[1]) > 0 && O_LEN(args[2]) > 0) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ nver = (nargs - 2) / 2
+ do i = 1, nver
+ Memr[rval2+i-1] = Memr[rval1+2*i+1]
+ do i = 1, nver
+ Memr[rval1+i-1] = Memr[rval1+2*i]
+ call me_polygon (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1], Memr[rval2], nver)
+ } else {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ nver = nargs / 2
+ do i = 1, nver
+ Memr[rval2+i-1] = Memr[rval1+2*i-1]
+ do i = 1, nver
+ Memr[rval1+i-1] = Memr[rval1+2*i-2]
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_polygon (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval2], nver)
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ }
+
+ case RG_COLS:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 2) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_cols (Memi[O_VALP(args[1])], Memi[O_VALP(o)], O_LEN(o),
+ O_VALC(args[2]))
+ } else if (nargs == 1) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_cols (Memi[ix], Memi[O_VALP(o)], O_LEN(o),
+ O_VALC(args[1]))
+ call mfree (ix, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_LINES:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 2) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_lines (Memi[O_VALP(args[1])], Memi[O_VALP(o)], O_LEN(o),
+ O_VALC(args[2]))
+ } else if (nargs == 1) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call amovki (RG_OPMV(rg,2), Memi[ix], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_lines (Memi[ix], Memi[O_VALP(o)], O_LEN(o),
+ O_VALC(args[1]))
+ call mfree (ix, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_VECTOR:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 7) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_vector (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], Memr[rval1+6])
+ } else if (nargs == 5) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_vector (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_PIE:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 6) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_pie (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], IM_LEN(RG_PMIM(rg),1),
+ IM_LEN(RG_PMIM(rg),2))
+ } else if (nargs == 4) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_pie (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3],
+ IM_LEN(RG_PMIM(rg),1), IM_LEN(RG_PMIM(rg),2))
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_CANNULUS:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 6) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_cannulus (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5])
+ } else if (nargs == 4) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_cannulus (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_EANNULUS:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 8) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_eannulus (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], Memr[rval1+6], Memr[rval1+7])
+ } else if (nargs == 6) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_eannulus (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_RANNULUS:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 8) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_rannulus (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5], Memr[rval1+6], Memr[rval1+7])
+ } else if (nargs == 6) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_rannulus (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1], Memr[rval1+2], Memr[rval1+3],
+ Memr[rval1+4], Memr[rval1+5])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+
+ case RG_PANNULUS:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs < 7) {
+ O_TYPE(o) = ERR
+ } else if (O_LEN(args[1]) > 0 && O_LEN(args[2]) > 0) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ width = Memr[rval1+2]
+ nver = (nargs - 3) / 2
+ do i = 1, nver
+ #Memr[rval2+i-1] = Memr[rval1+2*i+1]
+ Memr[rval2+i-1] = Memr[rval1+2*i+2]
+ do i = 1, nver
+ #Memr[rval1+i-1] = Memr[rval1+2*i+2]
+ Memr[rval1+i-1] = Memr[rval1+2*i+1]
+ call salloc (orval1, nver, TY_REAL)
+ call salloc (orval2, nver, TY_REAL)
+ call me_pyexpand (Memr[rval1], Memr[rval2], Memr[orval1],
+ Memr[orval2], nver, width)
+ call me_apolygon (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1], Memr[rval2],
+ Memr[orval1], Memr[orval2], nver)
+ } else {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ width = Memr[rval1]
+ nver = (nargs - 1) / 2
+ do i = 1, nver
+ Memr[rval2+i-1] = Memr[rval1+2*i]
+ do i = 1, nver
+ Memr[rval1+i-1] = Memr[rval1+2*i-1]
+ call salloc (orval1, nver, TY_REAL)
+ call salloc (orval2, nver, TY_REAL)
+ call me_pyexpand (Memr[rval1], Memr[rval2], Memr[orval1],
+ Memr[orval2], nver, width)
+ call me_apolygon (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval2], Memr[orval1], Memr[orval2], nver)
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ }
+
+ case RG_POINT:
+ O_LEN(o) = IM_LEN(RG_PMIM(rg),1)
+ O_TYPE(o) = TY_BOOL
+ if (nargs == 4) {
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_point (Memi[O_VALP(args[1])], Memi[O_VALP(args[2])],
+ Memi[O_VALP(o)], O_LEN(o), Memr[rval1+2], Memr[rval1+3])
+ } else if (nargs == 2) {
+ call malloc (ix, O_LEN(o), TY_INT)
+ call malloc (iy, O_LEN(o), TY_INT)
+ do i = 1, O_LEN(o)
+ Memi[ix+i-1] = i
+ call amovki (RG_OPMV(rg,2), Memi[iy], O_LEN(o))
+ call malloc (O_VALP(o), O_LEN(o), TY_INT)
+ call me_point (Memi[ix], Memi[iy], Memi[O_VALP(o)], O_LEN(o),
+ Memr[rval1], Memr[rval1+1])
+ call mfree (ix, TY_INT)
+ call mfree (iy, TY_INT)
+ } else {
+ O_TYPE(o) = ERR
+ }
+ default:
+ O_TYPE(o) = ERR
+ }
+
+ call sfree (sp)
+end
+
diff --git a/pkg/proto/maskexpr/mesetexpr.x b/pkg/proto/maskexpr/mesetexpr.x
new file mode 100644
index 00000000..40c2495f
--- /dev/null
+++ b/pkg/proto/maskexpr/mesetexpr.x
@@ -0,0 +1,36 @@
+# ME_SETEXPR -- Set the pixel mask region to the appropriate number.
+
+procedure me_setexpr (expr, pmim, pregno, pregval, verbose)
+
+char expr[ARB] #I the region expression
+pointer pmim #I the pixelmask image descriptor
+int pregno #I the current region number
+int pregval #I the current region value
+bool verbose #I print status messages ?
+
+pointer sp, chregval
+int nchars, stat
+int itoc(), me_rgmask()
+
+begin
+ call smark (sp)
+ call salloc (chregval, SZ_FNAME, TY_CHAR)
+ nchars = itoc (pregval, Memc[chregval], SZ_FNAME)
+ if (nchars <= 0) {
+ if (verbose) {
+ call printf (" Region value %d cannot be encoded\n")
+ call pargi (pregval)
+ }
+ } else {
+ stat = me_rgmask (expr, Memc[chregval], "p", pmim)
+ if (stat == ERR) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ }
+ }
+
+ call sfree (sp)
+end
+
diff --git a/pkg/proto/maskexpr/mesetreg.x b/pkg/proto/maskexpr/mesetreg.x
new file mode 100644
index 00000000..3fbe3f7b
--- /dev/null
+++ b/pkg/proto/maskexpr/mesetreg.x
@@ -0,0 +1,292 @@
+include <imset.h>
+include <plset.h>
+
+define RG_REGIONS "|circle|ellipse|box|rectangle|polygon|vector|columns|\
+lines|pie|cannulus|eannulus|rannulus|pannulus|point|"
+
+define RG_CIRCLE 1
+define RG_ELLIPSE 2
+define RG_BOX 3
+define RG_RECTANGLE 4
+define RG_POLYGON 5
+define RG_VECTOR 6
+define RG_COLUMNS 7
+define RG_LINES 8
+define RG_PIE 9
+define RG_CANNULUS 10
+define RG_EANNULUS 11
+define RG_RANNULUS 12
+define RG_PANNULUS 13
+define RG_POINT 14
+
+define MAX_NVERTICES 100
+
+# RG_SETREG -- Set the pixel mask region to the appropriate number.
+
+procedure me_setreg (region, pmim, pregno, pregval, verbose)
+
+char region[ARB] #I the region description
+pointer pmim #I the pixelmask image descriptor
+int pregno #I the current region number
+int pregval #I the current region value
+bool verbose #I print status messages ?
+
+real xc, yc, a, b, ratio, theta
+real x1, y1, x2, y2, width
+pointer sp, function, ufunction, pl, xver, yver, rangestr
+int nfuncs, nver, nold
+int strdic(), imstati(), nscan()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (function, SZ_FNAME, TY_CHAR)
+ call salloc (ufunction, SZ_FNAME, TY_CHAR)
+ call salloc (xver, MAX_NVERTICES, TY_REAL)
+ call salloc (yver, MAX_NVERTICES, TY_REAL)
+ call salloc (rangestr, SZ_FNAME, TY_CHAR)
+
+ # Determine the type of region.
+ call sscan (region)
+ call gargwrd (Memc[function], SZ_FNAME)
+ nfuncs = strdic (Memc[function], Memc[ufunction], SZ_FNAME, RG_REGIONS)
+ if (nfuncs <= 0) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ call sfree (sp)
+ return
+ }
+
+ pl = imstati (pmim, IM_PLDES)
+
+ switch (nfuncs) {
+
+ case RG_CIRCLE:
+ call gargr (xc)
+ call gargr (yc)
+ call gargr (a)
+ if (nscan() < 4) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_circle (pl, xc, yc, a, PIX_SRC+PIX_VALUE(pregval))
+ }
+
+ case RG_ELLIPSE:
+ call gargr (xc)
+ call gargr (yc)
+ call gargr (a)
+ call gargr (ratio)
+ call gargr (theta)
+ if (nscan() < 6) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_ellipse (pl, xc, yc, a, ratio, theta,
+ PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_BOX:
+ call gargr (x1)
+ call gargr (y1)
+ call gargr (x2)
+ call gargr (y2)
+ if (nscan() < 5) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_box (pl, x1, y1, x2, y2, PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_RECTANGLE:
+ call gargr (xc)
+ call gargr (yc)
+ call gargr (a)
+ call gargr (ratio)
+ call gargr (theta)
+ if (nscan() < 6) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_rectangle (pl, xc, yc, a, ratio, theta,
+ PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_POLYGON:
+ nver = 0
+ repeat {
+ nold = nscan()
+ call gargr (Memr[xver+nver])
+ call gargr (Memr[yver+nver])
+ if ((nscan() - nold) == 2)
+ nver = nver + 1
+ else
+ break
+ } until ((nscan() - nold) < 2)
+ if (nver <3 ) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_polygon (pl, Memr[xver], Memr[yver], nver,
+ PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_VECTOR:
+ call gargr (x1)
+ call gargr (y1)
+ call gargr (x2)
+ call gargr (y2)
+ call gargr (width)
+ if (nscan() < 6) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_vector (pl, x1, y1, x2, y2, width,
+ PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_COLUMNS:
+ call gargwrd (Memc[rangestr], SZ_FNAME)
+ if (nscan() < 2) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_cols (pl, Memc[rangestr],
+ PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_LINES:
+ call gargwrd (Memc[rangestr], SZ_FNAME)
+ if (nscan() < 2) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_lines (pl, Memc[rangestr],
+ PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_PIE:
+ call gargr (xc)
+ call gargr (yc)
+ call gargr (a)
+ call gargr (b)
+ if (nscan() < 5) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_pie (pl, xc, yc, a, b, PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_CANNULUS:
+ call gargr (xc)
+ call gargr (yc)
+ call gargr (a)
+ call gargr (b)
+ if (nscan() < 5) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_cannulus (pl, xc, yc, a, b,
+ PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_EANNULUS:
+ call gargr (xc)
+ call gargr (yc)
+ call gargr (a)
+ call gargr (b)
+ call gargr (ratio)
+ call gargr (theta)
+ if (nscan() < 7) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_eannulus (pl, xc, yc, a, b, ratio, theta,
+ PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_RANNULUS:
+ call gargr (xc)
+ call gargr (yc)
+ call gargr (a)
+ call gargr (b)
+ call gargr (ratio)
+ call gargr (theta)
+ if (nscan() < 7) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_rannulus (pl, xc, yc, a, b, ratio, theta,
+ PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_PANNULUS:
+ call gargr (b)
+ if (nscan () < 2) {
+ nver = 0
+ } else {
+ nver = 0
+ repeat {
+ nold = nscan()
+ call gargr (Memr[xver+nver])
+ call gargr (Memr[yver+nver])
+ if ((nscan() - nold) == 2)
+ nver = nver + 1
+ else
+ break
+ } until ((nscan() - nold) < 2)
+ }
+ if (nver < 3) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_apolygon (pl, b, Memr[xver], Memr[yver], nver,
+ PIX_SRC + PIX_VALUE(pregval))
+ }
+
+ case RG_POINT:
+ call gargr (xc)
+ call gargr (yc)
+ if (nscan() < 3) {
+ if (verbose) {
+ call printf (" Region %d cannot be decoded\n")
+ call pargi (pregno)
+ }
+ } else {
+ call pe_point (pl, xc, yc, PIX_SRC+PIX_VALUE(pregval))
+ }
+
+ default:
+ ;
+ }
+
+ call sfree (sp)
+end
diff --git a/pkg/proto/maskexpr/mkpkg b/pkg/proto/maskexpr/mkpkg
new file mode 100644
index 00000000..ee3e86db
--- /dev/null
+++ b/pkg/proto/maskexpr/mkpkg
@@ -0,0 +1,26 @@
+# Make the MSKEXPR and MSKREGIONS tasks
+
+$checkout libpkg.a ../
+$update libpkg.a
+$checkin libpkg.a ../
+$exit
+
+libpkg.a:
+ t_mskexpr.x <fset.h> <ctype.h> <imhdr.h>
+ memkmask.x <mach.h> <ctype.h> <imhdr.h> <imset.h> <pmset.h> \
+ <evvexpr.h>
+
+ t_mskregions.x <fset.h> <ctype.h> <imhdr.h> <imset.h> <pmset.h>
+ mesetreg.x <plset.h> <imset.h>
+ mesetexpr.x
+ meregmask.x <mach.h> <ctype.h> <imhdr.h> <imset.h> <pmset.h> \
+ <evvexpr.h>
+ peregfuncs.x <plset.h> <plio.h> <math.h> peregfuncs.h
+ peregufcn.x <plset.h> <plio.h> <math.h> peregfuncs.h
+ megeom.x <math.h>
+
+ meregfuncs.x <mach.h> <ctype.h> <math.h>
+ mskexpand.x <ctotok.h> <ctype.h> gettok.h
+ megsym.x <ctotok.h> <ctype.h> gettok.h
+ gettok.x <syserr.h> <error.h> <fset.h> <ctype.h> gettok.h
+ ;
diff --git a/pkg/proto/maskexpr/mskexpand.x b/pkg/proto/maskexpr/mskexpand.x
new file mode 100644
index 00000000..5fb6cc9d
--- /dev/null
+++ b/pkg/proto/maskexpr/mskexpand.x
@@ -0,0 +1,261 @@
+include <ctotok.h>
+include <ctype.h>
+include "gettok.h"
+
+# Some definitions.
+
+# Default symbol table size limits.
+define DEF_LENINDEX 97
+define DEF_LENSTAB 1024
+define DEF_LENSBUF 8192
+
+# Expression database symbol.
+define LEN_SYM 2
+define SYM_TEXT Memi[$1]
+define SYM_NARGS Memi[$1+1]
+
+# Argument list symbol
+define LEN_ARGSYM 1
+define ARGNO Memi[$1]
+
+
+# ME_GETEXPRDB -- Read the expression database into a symbol table. The
+# input file has the following structure:
+#
+# <symbol>['(' arg-list ')'][':'|'='] replacement-text
+#
+# Symbols must be at the beginning of a line. The expression text is
+# terminated by a nonempty, noncomment line with no leading whitespace.
+
+pointer procedure me_getexprdb (fname)
+
+char fname[ARB] #I file to be read
+
+pointer sym, sp, lbuf, st, a_st, ip, symname, tokbuf, text
+int tok, fd, line, nargs, op, token, buflen, offset, stpos, n
+pointer stopen(), stenter()
+int open(), getlline(), ctotok(), stpstr()
+errchk open, getlline, stopen, stenter, me_puttok
+
+define skip_ 91
+
+begin
+ call smark (sp)
+ call salloc (lbuf, SZ_COMMAND, TY_CHAR)
+ call salloc (text, SZ_COMMAND, TY_CHAR)
+ call salloc (tokbuf, SZ_COMMAND, TY_CHAR)
+ call salloc (symname, SZ_FNAME, TY_CHAR)
+
+ fd = open (fname, READ_ONLY, TEXT_FILE)
+ st = stopen ("imexpr", DEF_LENINDEX, DEF_LENSTAB, DEF_LENSBUF)
+ a_st = stopen ("args", DEF_LENINDEX, DEF_LENSTAB, DEF_LENSBUF)
+ line = 0
+
+ while (getlline (fd, Memc[lbuf], SZ_COMMAND) != EOF) {
+ line = line + 1
+ ip = lbuf
+
+ # Skip comments and blank lines.
+ while (IS_WHITE(Memc[ip]))
+ ip = ip + 1
+ if (Memc[ip] == '\n' || Memc[ip] == '#')
+ next
+
+ # Get symbol name.
+ if (ctotok (Memc,ip,Memc[symname],SZ_FNAME) != TOK_IDENTIFIER) {
+ call eprintf ("exprdb: expected identifier at line %d\n")
+ call pargi (line)
+skip_ while (getlline (fd, Memc[lbuf], SZ_COMMAND) != EOF) {
+ line = line + 1
+ if (Memc[lbuf] == '\n')
+ break
+ }
+ }
+
+ call stmark (a_st, stpos)
+
+ # Check for the optional argument-symbol list. Allow only a
+ # single space between the symbol name and its argument list,
+ # otherwise we can't tell the difference between an argument
+ # list and the parenthesized expression which follows.
+
+ if (Memc[ip] == ' ')
+ ip = ip + 1
+
+ if (Memc[ip] == '(') {
+ ip = ip + 1
+ n = 0
+ repeat {
+ tok = ctotok (Memc, ip, Memc[tokbuf], SZ_FNAME)
+ if (tok == TOK_IDENTIFIER) {
+ sym = stenter (a_st, Memc[tokbuf], LEN_ARGSYM)
+ n = n + 1
+ ARGNO(sym) = n
+ } else if (Memc[tokbuf] == ',') {
+ ;
+ } else if (Memc[tokbuf] != ')') {
+ call eprintf ("exprdb: bad arglist at line %d\n")
+ call pargi (line)
+ call stfree (a_st, stpos)
+ goto skip_
+ }
+ } until (Memc[tokbuf] == ')')
+ }
+
+ # Check for the optional ":" or "=".
+ while (IS_WHITE(Memc[ip]))
+ ip = ip + 1
+ if (Memc[ip] == ':' || Memc[ip] == '=')
+ ip = ip + 1
+
+ # Accumulate the expression text.
+ buflen = SZ_COMMAND
+ op = 1
+
+ repeat {
+ repeat {
+ token = ctotok (Memc, ip, Memc[tokbuf], SZ_COMMAND)
+ if (Memc[tokbuf] == '#')
+ break
+ else if (token != TOK_EOS && token != TOK_NEWLINE)
+ call me_puttok (a_st, text, op, buflen, Memc[tokbuf])
+ } until (token == TOK_EOS)
+
+ if (getlline (fd, Memc[lbuf], SZ_COMMAND) == EOF)
+ break
+ else
+ line = line + 1
+
+ for (ip=lbuf; IS_WHITE(Memc[ip]); ip=ip+1)
+ ;
+ if (ip == lbuf) {
+ call ungetline (fd, Memc[lbuf])
+ line = line - 1
+ break
+ }
+ }
+
+ # Free any argument list symbols.
+ call stfree (a_st, stpos)
+
+ # Scan the expression text and count the number of $N arguments.
+ nargs = 0
+ for (ip=text; Memc[ip] != EOS; ip=ip+1)
+ if (Memc[ip] == '$' && IS_DIGIT(Memc[ip+1])) {
+ nargs = max (nargs, TO_INTEG(Memc[ip+1]))
+ ip = ip + 1
+ }
+
+ # Enter symbol in table.
+ sym = stenter (st, Memc[symname], LEN_SYM)
+ offset = stpstr (st, Memc[text], 0)
+ SYM_TEXT(sym) = offset
+ SYM_NARGS(sym) = nargs
+ }
+
+ call stclose (a_st)
+ call sfree (sp)
+
+ return (st)
+end
+
+
+# ME_PUTTOK -- Append a token string to a text buffer.
+
+procedure me_puttok (a_st, text, op, buflen, token)
+
+pointer a_st #I argument-symbol table
+pointer text #U text buffer
+int op #U output pointer
+int buflen #U buffer length, chars
+char token[ARB] #I token string
+
+pointer sym
+int ip, ch1, ch2
+pointer stfind()
+errchk realloc
+
+begin
+ # Replace any symbolic arguments by "$N".
+ if (a_st != NULL && IS_ALPHA(token[1])) {
+ sym = stfind (a_st, token)
+ if (sym != NULL) {
+ token[1] = '$'
+ token[2] = TO_DIGIT(ARGNO(sym))
+ token[3] = EOS
+ }
+ }
+
+ # Append the token string to the text buffer.
+ for (ip=1; token[ip] != EOS; ip=ip+1) {
+ if (op + 1 > buflen) {
+ buflen = buflen + SZ_COMMAND
+ call realloc (text, buflen, TY_CHAR)
+ }
+
+ # The following is necessary because ctotok parses tokens such as
+ # "$N", "==", "!=", etc. as two tokens. We need to rejoin these
+ # characters to make one token.
+
+ if (op > 1 && token[ip+1] == EOS) {
+ ch1 = Memc[text+op-3]
+ ch2 = token[ip]
+
+ if (ch1 == '$' && IS_DIGIT(ch2))
+ op = op - 1
+ else if (ch1 == '*' && ch2 == '*')
+ op = op - 1
+ else if (ch1 == '/' && ch2 == '/')
+ op = op - 1
+ else if (ch1 == '<' && ch2 == '=')
+ op = op - 1
+ else if (ch1 == '>' && ch2 == '=')
+ op = op - 1
+ else if (ch1 == '=' && ch2 == '=')
+ op = op - 1
+ else if (ch1 == '!' && ch2 == '=')
+ op = op - 1
+ else if (ch1 == '?' && ch2 == '=')
+ op = op - 1
+ else if (ch1 == '&' && ch2 == '&')
+ op = op - 1
+ else if (ch1 == '|' && ch2 == '|')
+ op = op - 1
+ }
+
+ Memc[text+op-1] = token[ip]
+ op = op + 1
+ }
+
+ # Append a space to ensure that tokens are delimited.
+ Memc[text+op-1] = ' '
+ op = op + 1
+
+ Memc[text+op-1] = EOS
+end
+
+
+# ME_EXPANDTEXT -- Scan an expression, performing macro substitution on the
+# contents and returning a fully expanded string.
+
+pointer procedure me_expandtext (st, expr)
+
+pointer st #I symbol table (macros)
+char expr[ARB] #I input expression
+
+pointer buf, gt
+int buflen, nchars
+int locpr(), gt_expand()
+pointer gt_opentext()
+extern me_gsym()
+
+begin
+ buflen = SZ_COMMAND
+ call malloc (buf, buflen, TY_CHAR)
+
+ gt = gt_opentext (expr, locpr(me_gsym), st, 0, GT_NOFILE)
+ nchars = gt_expand (gt, buf, buflen)
+ call gt_close (gt)
+
+ return (buf)
+end
diff --git a/pkg/proto/maskexpr/peregfuncs.h b/pkg/proto/maskexpr/peregfuncs.h
new file mode 100644
index 00000000..cc777a9a
--- /dev/null
+++ b/pkg/proto/maskexpr/peregfuncs.h
@@ -0,0 +1,131 @@
+# PEREGFUNCS.H -- Structure definitions.
+
+# Circle definitions.
+
+define LEN_CIRCLEDES 5
+define C_PL Memi[$1] # reference mask
+define C_XCEN Memr[P2R($1+1)] # X center of circle
+define C_YCEN Memr[P2R($1+2)] # Y center of circle
+define C_RADIUS Memr[P2R($1+3)] # radius of circle
+define C_PV Memi[$1+4] # pixel value
+
+
+# Ellipse definitions.
+
+define LEN_ELLDES 10
+define E_PL Memi[$1] # reference mask
+define E_XCEN Memr[P2R($1+1)] # X center of ellipse
+define E_YCEN Memr[P2R($1+2)] # Y center of ellipse
+define E_AA Memr[P2R($1+3)] # aa parameter
+define E_BB Memr[P2R($1+4)] # bb parameter
+define E_CC Memr[P2R($1+5)] # cc parameter
+define E_FF Memr[P2R($1+6)] # ff paramater
+define E_DXMAX Memr[P2R($1+7)] # the maximum x offset
+define E_DYMAX Memr[P2R($1+8)] # the maximum x offset
+define E_PV Memi[$1+9] # pixel value
+
+
+# Box definitions.
+
+define LEN_BOXDES 6
+define B_PL Memi[$1] # reference mask
+define B_X1 Memr[P2R($1+1)] # X1 lower left corner of box
+define B_Y1 Memr[P2R($1+2)] # Y1 lower left corner of box
+define B_X2 Memr[P2R($1+3)] # X2 upper right corner of box
+define B_Y2 Memr[P2R($1+4)] # Y2 upper right corner of box
+define B_PV Memi[$1+5] # pixel value
+
+
+# Polygon definitions.
+
+define TOL 0.0001 # pixel units
+define swapi {tempi=$2;$2=$1;$1=tempi}
+define swapr {tempr=$2;$2=$1;$1=tempr}
+define equal (abs($1-$2)<TOL)
+
+define LEN_PGONDES 7
+define P_PL Memi[$1] # pointer to X vector
+define P_XP Memi[$1+1] # pointer to X vector
+define P_YP Memi[$1+2] # pointer to Y vector
+define P_OO Memi[$1+3] # pointer to previous range list
+define P_OY Memi[$1+4] # y value of previous range list
+define P_NS Memi[$1+5] # number of line segments
+define P_PV Memi[$1+6] # pixel value
+
+
+# Circular annulus definitions.
+
+define LEN_CANNDES 6
+define CA_PL Memi[$1] # reference mask
+define CA_XCEN Memr[P2R($1+1)] # x center of circle
+define CA_YCEN Memr[P2R($1+2)] # y center of circle
+define CA_RADIUS1 Memr[P2R($1+3)] # inner radius of annulus
+define CA_RADIUS2 Memr[P2R($1+4)] # outer radius of annulus
+define CA_PV Memi[$1+5] # pixel value
+
+
+# Elliptical annulus defintiions.
+
+define LEN_EANNDES 16
+define EA_PL Memi[$1] # reference mask
+define EA_XCEN Memr[P2R($1+1)] # x center of ellipse
+define EA_YCEN Memr[P2R($1+2)] # y center of ellipse
+define EA_AA1 Memr[P2R($1+3)] # aa parameter for inner ellipse
+define EA_BB1 Memr[P2R($1+4)] # bb parameter for inner ellipse
+define EA_CC1 Memr[P2R($1+5)] # cc parameter for inner ellipse
+define EA_FF1 Memr[P2R($1+6)] # ff parameter for inner ellipse
+define EA_DXMAX1 Memr[P2R($1+7)] # max dx value for inner ellipse
+define EA_DYMAX1 Memr[P2R($1+8)] # max dy value for inner ellipse
+define EA_AA2 Memr[P2R($1+9)] # aa parameter for outer ellipse
+define EA_BB2 Memr[P2R($1+10)] # bb parameter for outer ellipse
+define EA_CC2 Memr[P2R($1+11)] # cc parameter for outer ellipse
+define EA_FF2 Memr[P2R($1+12)] # ff parameter for outer ellipse
+define EA_DXMAX2 Memr[P2R($1+13)] # max dx value for outer ellipse
+define EA_DYMAX2 Memr[P2R($1+14)] # max dy value for outer ellipse
+define EA_PV Memi[$1+15] # pixel value
+
+
+# Rasterop annulus definitions.
+
+define LEN_RANNDES 7
+define RA_PL Memi[$1] # the mask descriptor
+define RA_IXP Memi[$1+1] # pointer to inner polygon X vector
+define RA_IYP Memi[$1+2] # pointer to inner Y polygon vector
+define RA_OXP Memi[$1+3] # pointer to outer X polygon vector
+define RA_OYP Memi[$1+4] # pointer to outer Y polygon vector
+define RA_NVER Memi[$1+5] # number of vertices
+define RA_PV Memi[$1+6] # mask pixel value
+
+
+# Polygon annulus definitions.
+
+define LEN_PAGONDES 7
+define PA_PL Memi[$1] # the mask descriptor
+define PA_IXP Memi[$1+1] # pointer to inner polygon X vector
+define PA_IYP Memi[$1+2] # pointer to inner Y polygon vector
+define PA_OXP Memi[$1+3] # pointer to outer X polygon vector
+define PA_OYP Memi[$1+4] # pointer to outer Y polygon vector
+define PA_NVER Memi[$1+5] # number of vertices
+define PA_PV Memi[$1+6] # mask pixel value
+
+
+# Column definitions.
+
+define LEN_COLSDES 4
+define L_PL Memi[$1] # reference mask
+define L_RANGES Memi[$1+1] # pointer to the ranges
+define L_NRANGES Memi[$1+2] # the number of ranges
+define L_XS Memi[$1+3] # the starting x coordinate value
+define L_NPIX Memi[$1+4] # the number of pixels value
+define L_PV Memi[$1+5] # pixel value
+
+
+# Line definitions.
+
+define LEN_LINESDES 3
+define L_PL Memi[$1] # reference mask
+define L_RANGES Memi[$1+1] # pointer to the ranges
+define L_PV Memi[$1+2] # pixel value
+
+define MAX_NRANGES 100
+define SMALL_NUMBER 1.0e-24
diff --git a/pkg/proto/maskexpr/peregfuncs.x b/pkg/proto/maskexpr/peregfuncs.x
new file mode 100644
index 00000000..9e79d422
--- /dev/null
+++ b/pkg/proto/maskexpr/peregfuncs.x
@@ -0,0 +1,877 @@
+include <math.h>
+include <plset.h>
+include <plio.h>
+include "peregfuncs.h"
+
+
+# PE_POINT -- Rasterop between a point region as source and an existing
+# mas as destination.
+
+procedure pe_point (pl, x, y, rop)
+
+pointer pl #I mask descriptor
+real x,y #I center coords of circle
+int rop #I rasterop
+
+begin
+ call pl_point (pl, nint(x), nint(y), rop)
+end
+
+
+# PE_CIRCLE -- Rasterop between a circular region as source and an existing
+# mask as destination. It is not necessary for the center of the circle to
+# be inside the mask; if it is outside, the boundary of the circle will be
+# clipped to the boundary of the mask. This is a 2-dim operator. If the
+# image dimensionality is greater than two the pl_setplane procedure should
+# be called first to specify the plane to be modified. These routines are
+# a modification of the ones in plio$plcircle. The main difference is
+# that the x, y, radius parameters are initially set to real numbers not
+# integers.
+
+procedure pe_circle (pl, x, y, radius, rop)
+
+pointer pl #I mask descriptor
+real x,y #I center coords of circle
+real radius #I radius of circle
+int rop #I rasterop
+
+real y1r, y2r, x1r, x2r
+int y1, y2
+pointer sp, ufd
+extern pe_ucircle()
+
+begin
+ # Not sure why we need to call this routine here.
+ #call plvalid (pl)
+
+ # Test the line and column limits.
+ y1r = y - radius
+ y2r = y + radius
+ x1r = x - radius
+ x2r = x + radius
+ if ((y2r < 0.5) || (y1r > PL_AXLEN(pl,2) + 0.5))
+ return
+ if ((x2r < 0.5) || (x1r > PL_AXLEN(pl,1) + 0.5))
+ return
+
+ call smark (sp)
+ call salloc (ufd, LEN_CIRCLEDES, TY_STRUCT)
+
+ y1 = max ( 1, min (PL_AXLEN(pl,2), int(y1r)))
+ y2 = max (y1, min (PL_AXLEN(pl,2), int(y2r)))
+
+ C_PL(ufd) = pl
+ C_XCEN(ufd) = x
+ C_YCEN(ufd) = y
+ C_RADIUS(ufd) = radius
+ C_PV(ufd) = 1
+
+ call pl_regionrop (pl, pe_ucircle, ufd, y1, y2, rop)
+
+ call sfree (sp)
+end
+
+
+# PE_ELLIPSE -- Rasterop between an elliptical region as source and an existing
+# mask as destination. It is not necessary for the center of the ellipse to
+# be inside the mask; if it is outside, the boundary of the ellipse will be
+# clipped to the boundary of the mask. This is a 2-dim operator. If the
+# image dimensionality is greater than two the pl_setplane procedure should
+# be called first to specify the plane to be modified. These routines are
+# a modification of the ones in plio$plcircle. The main difference is
+# that the x, y, radius parameters are initially set to real numbers not
+# integers.
+
+procedure pe_ellipse (pl, x, y, radius, ratio, theta, rop)
+
+pointer pl #I mask descriptor
+real x,y #I center coords of ellipse
+real radius #I semi-major axis of ellipse
+real ratio #I the ratio semi-minor / semi-major axes
+real theta #I position angle in degrees
+int rop #I rasterop
+
+real aa, bb, cc, ff, dx, dy
+real y1r, y2r, x1r, x2r, r2
+int y1, y2
+pointer sp, ufd
+extern pe_uellipse()
+
+begin
+ # Not sure why we need to call this routine here.
+ #call plvalid (pl)
+
+ # Get ellipse parameters.
+ call me_ellgeom (radius, ratio, theta, aa, bb, cc, ff)
+ r2 = radius * radius
+ dx = ff / (aa - bb * bb / 4.0 / cc)
+ if (dx > 0.0)
+ dx = sqrt (dx)
+ else
+ dx = 0.0
+ dy = ff / (cc - bb * bb / 4.0 / aa)
+ if (dy > 0.0)
+ dy = sqrt (dy)
+ else
+ dy = 0.0
+
+ # Test the line and column limits.
+ y1r = y - dy
+ y2r = y + dy
+ x1r = x - dx
+ x2r = x + dx
+ if ((y2r < 0.5) || (y1r > PL_AXLEN(pl,2) + 0.5))
+ return
+ if ((x2r < 0.5) || (x1r > PL_AXLEN(pl,1) + 0.5))
+ return
+
+ call smark (sp)
+ call salloc (ufd, LEN_ELLDES, TY_STRUCT)
+ y1 = max ( 1, min (PL_AXLEN(pl,2), int(y1r)))
+ y2 = max (y1, min (PL_AXLEN(pl,2), int(y2r)))
+
+ E_PL(ufd) = pl
+ E_XCEN(ufd) = x
+ E_YCEN(ufd) = y
+ E_DXMAX(ufd) = dx
+ E_DYMAX(ufd) = dy
+ E_AA(ufd) = aa / r2
+ E_BB(ufd) = bb / r2
+ E_CC(ufd) = cc / r2
+ E_FF(ufd) = ff / r2
+ E_PV(ufd) = 1
+
+ call pl_regionrop (pl, pe_uellipse, ufd, y1, y2, rop)
+
+ call sfree (sp)
+end
+
+
+# PE_BOX -- Rasterop between a rectangular region as source and an existing
+# mask as destination. It is not necessary for the corners of the box to
+# be inside the mask; if they are outside, the boundary of the box will be
+# clipped to the boundary of the mask. This is a 2-dim operator. If the
+# image dimensionality is greater than two the pl_setplane procedure should
+# be called first to specify the plane to be modified. These routines are
+# a modification of the ones in plio$plbox. The main difference is
+# that the x, y, radius parameters are initially set to real numbers not
+# integers.
+
+procedure pe_box (pl, x1, y1, x2, y2, rop)
+
+pointer pl #I mask descriptor
+real x1, y1 #I lower left corner of box
+real x2, y2 #I upper right corner of box
+int rop #I rasterop
+
+pointer sp, ufd
+extern pe_ubox()
+
+begin
+ # Not sure why we need to call this routine here.
+ #call plvalid (pl)
+
+ # Test the line and column limits.
+ if ((y2 < 0.5) || (y1 > PL_AXLEN(pl,2) + 0.5))
+ return
+ if ((x2 < 0.5) || (x1 > PL_AXLEN(pl,1) + 0.5))
+ return
+
+ call smark (sp)
+ call salloc (ufd, LEN_BOXDES, TY_STRUCT)
+
+ B_PL(ufd) = pl
+ B_X1(ufd) = max (1, min (PL_AXLEN(pl,1), nint(x1)))
+ B_Y1(ufd) = max (1, min (PL_AXLEN(pl,2), nint(y1)))
+ B_X2(ufd) = max (1, min (PL_AXLEN(pl,1), nint(x2)))
+ B_Y2(ufd) = max (1, min (PL_AXLEN(pl,2), nint(y2)))
+ B_PV(ufd) = 1
+
+ call pl_regionrop (pl, pe_ubox, ufd, int(B_Y1(ufd)), int(B_Y2(ufd)),
+ rop)
+
+ call sfree (sp)
+end
+
+
+# PE_RECTANGLE -- Rasterop between a rectangular region as source and an
+# existing mask as destination. It is not necessary for the center of the
+# rectangle to be inside the mask; if it is outside, the boundary of the
+# rectangle will be clipped to the boundary of the mask. This is a 2-dim
+# operator. If the image dimensionality is greater than two the pl_setplane
+# procedure should be called first to specify the plane to be modified.
+# These routines are a modification of the ones in plio$plcircle. The main
+# difference is that the x, y, radius parameters are initially set to real
+# numbers not integers.
+
+procedure pe_rectangle (pl, x, y, radius, ratio, theta, rop)
+
+pointer pl #I mask descriptor
+real x,y #I center coords of rectangle
+real radius #I semi-major axis of rectangle
+real ratio #I the ratio semi-minor / semi-major axes
+real theta #I position angle in degrees
+int rop #I rasterop
+
+real xr[4], yr[4]
+int i
+
+begin
+ # Get rectangle vertices.
+ call me_rectgeom (radius, ratio, theta, xr, yr)
+ do i = 1, 4 {
+ xr[i] = x + xr[i]
+ yr[i] = y + yr[i]
+ }
+
+ # Mark the polygon.
+ call pe_polygon (pl, xr, yr, 4, rop)
+end
+
+
+# PE_VECTOR -- Rasterop between a rectangular vector region as source and an
+# existing mask as destination. It is not necessary for the center of the
+# rectangle to be inside the mask; if it is outside, the boundary of the
+# rectangle will be clipped to the boundary of the mask. This is a 2-dim
+# operator. If the image dimensionality is greater than two the pl_setplane
+# procedure should be called first to specify the plane to be modified.
+# These routines are a modification of the ones in plio$plcircle. The main
+# difference is that the x, y, radius parameters are initially set to real
+# numbers not integers.
+
+procedure pe_vector (pl, x1, y1, x2, y2, width, rop)
+
+pointer pl #I mask descriptor
+real x1, y1 #I beginning point of vector
+real x2, y2 #I ending point of vector
+real width #I width of vector
+int rop #I rasterop
+
+real xr[4], yr[4]
+real xc, yc, radius, ratio, theta
+int i
+
+begin
+ # Compute the center of the rectangle.
+ xc = (x1 + x2) / 2.0
+ yc = (y1 + y2) / 2.0
+
+ # Compute the semi-major axis, axis ratio, and position angle.
+ radius = sqrt ((x2 - x1) ** 2 + (y2 - y1) ** 2) / 2.0
+ if (radius <= 0.0)
+ return
+ ratio = width / radius
+ theta = RADTODEG (atan2 (y2 - y1, x2 - x1))
+
+ # Get rectangle vertices.
+ call me_rectgeom (radius, ratio, theta, xr, yr)
+
+ # Add back in the center coordinates.
+ do i = 1, 4 {
+ xr[i] = xc + xr[i]
+ yr[i] = yc + yr[i]
+ }
+
+ # Mark the polygon.
+ call pe_polygon (pl, xr, yr, 4, rop)
+end
+
+
+# PE_POLYGON -- Perform a rasterop operation on the area enclosed by a polygon
+# drawn in a 2-dimensional plane of a mask. If the dimensionality of the mask
+# exceeds 2, the pl_setplane() procedure should be called first to define the
+# plane of the mask to be modified.
+
+procedure pe_polygon (pl, x, y, npts, rop)
+
+pointer pl #I mask descriptor
+real x[npts] #I polygon x-vertices
+real y[npts] #I polygon y-vertices
+int npts #I number of points in polygon
+int rop #I rasterop defining operation
+
+real line_1r, line_2r
+pointer sp, ufd, xp, yp, oo
+int line_1, line_2, i
+extern pe_upolygon()
+errchk plvalid
+
+begin
+ # Note sure why this is called.
+ #call plvalid (pl)
+ if (npts < 3)
+ return
+
+ call smark (sp)
+ call salloc (ufd, LEN_PGONDES, TY_STRUCT)
+ call salloc (oo, RL_FIRST + (npts+1)*3, TY_INT)
+ call salloc (xp, npts + 1, TY_REAL)
+ call salloc (yp, npts + 1, TY_REAL)
+
+ # Initialize the region descriptor.
+ P_PL(ufd) = pl
+ P_XP(ufd) = xp
+ P_YP(ufd) = yp
+ P_PV(ufd) = 1
+ P_OO(ufd) = oo
+ P_OY(ufd) = -1
+ P_NS(ufd) = npts - 1
+ RLI_LEN(oo) = 0
+
+ # Copy the user supplied polygon vertices into the descriptor,
+ # normalizing the polygon in the process.
+
+ do i = 1, npts {
+ Memr[xp+i-1] = x[i]
+ Memr[yp+i-1] = y[i]
+ }
+
+ if (abs(x[1]-x[npts]) > TOL || abs(y[1]-y[npts]) > TOL) {
+ Memr[xp+npts] = x[1]
+ Memr[yp+npts] = y[1]
+ P_NS(ufd) = npts
+ }
+
+ # Compute the range in Y in which the polygon should be drawn.
+ call alimr (y, npts, line_1r, line_2r)
+ line_1 = max (1, min (PL_AXLEN(pl,2), int (line_1r)))
+ line_2 = max (line_1, min (PL_AXLEN(pl,2), int (line_2r)))
+
+ call pl_regionrop (pl, pe_upolygon, ufd, line_1, line_2, rop)
+
+ call sfree (sp)
+end
+
+
+# PE_CANNULUS -- Rasterop between a circular annular region as source and an
+# existing mask as destination. It is not necessary for the center of the
+# annulus to be inside the mask; if it is outside, the boundary of the
+# annulus will be clipped to the boundary of the mask. This is a 2-dim
+# operator. If the image dimensionality is greater than two the pl_setplane
+# procedure should be called first to specify the plane to be modified. These
+# routines are a modification of the ones in plio$plcircle. The main difference
+# is that the x, y, radius1, radius2, parameters are initially set to real
+# numbers not integers.
+
+procedure pe_cannulus (pl, x, y, radius1, radius2, rop)
+
+pointer pl #I mask descriptor
+real x,y #I center coords of circular annulus
+real radius1 #I inner radius of circular annulus
+real radius2 #I outer radius of circular annulus
+int rop #I rasterop
+
+real y1r, y2r, x1r, x2r
+int y1, y2
+pointer sp, ufd
+extern pe_ucannulus()
+
+begin
+ # Not sure why we need to call this routine here
+ #call plvalid (pl)
+
+ # The outer annulus must be greater than or equal to the inner annulus
+ if (radius2 < radius1)
+ return
+
+ # Test image limits.
+ y1r = y - radius2
+ y2r = y + radius2
+ if ((y2r < 0.5) || (y1r > (PL_AXLEN(pl,2) + 0.5)))
+ return
+ x1r = x - radius2
+ x2r = x + radius2
+ if ((x2r < 0.5) || (x1r > (PL_AXLEN(pl,1) + 0.5)))
+ return
+
+ call smark (sp)
+ call salloc (ufd, LEN_CANNDES, TY_STRUCT)
+
+ y1 = max ( 1, min (PL_AXLEN(pl,2), int(y1r)))
+ y2 = max (y1, min (PL_AXLEN(pl,2), int(y2r)))
+
+ CA_PL(ufd) = pl
+ CA_XCEN(ufd) = x
+ CA_YCEN(ufd) = y
+ CA_RADIUS1(ufd) = radius1
+ CA_RADIUS2(ufd) = radius2
+ CA_PV(ufd) = 1
+
+ call pl_regionrop (pl, pe_ucannulus, ufd, y1, y2, rop)
+
+ call sfree (sp)
+end
+
+
+# PE_EANNULUS -- Rasterop between an elliptical annular region as source and an
+# existing mask as destination. It is not necessary for the center of the
+# annulus to be inside the mask; if it is outside, the boundary of the
+# annulus will be clipped to the boundary of the mask. This is a 2-dim
+# operator. If the image dimensionality is greater than two the pl_setplane
+# procedure should be called first to specify the plane to be modified. These
+# routines are a modification of the ones in plio$plcircle. The main difference
+# is that the x, y, radius1, radius2, parameters are initially set to real
+# numbers not integers.
+
+procedure pe_eannulus (pl, x, y, radius1, radius2, ratio, theta, rop)
+
+pointer pl #I mask descriptor
+real x,y #I center coords of circular annulus
+real radius1 #I inner radius of circular annulus
+real radius2 #I outer radius of circular annulus
+real ratio #I the semi-minor / semi-major axis ratio
+real theta #I the position angle in degrees
+int rop #I rasterop
+
+real aa, bb, cc, ff, r2, dx, dy
+real y1r, y2r, x1r, x2r
+int y1, y2
+pointer sp, ufd
+extern pe_ueannulus()
+
+begin
+ # Not sure why we need to call this routine here
+ #call plvalid (pl)
+
+ # The outer annulus must be greater than or equal to the inner annulus
+ if (radius2 < radius1)
+ return
+
+ # Get the outer ellipse parameters.
+ call me_ellgeom (radius2, ratio, theta, aa, bb, cc, ff)
+ r2 = radius2 * radius2
+ dx = ff / (aa - bb * bb / 4.0 / cc)
+ if (dx > 0.0)
+ dx = sqrt (dx)
+ else
+ dx = 0.0
+ dy = ff / (cc - bb * bb / 4.0 / aa)
+ if (dy > 0.0)
+ dy = sqrt (dy)
+ else
+ dy = 0.0
+
+ # Test image limits.
+ y1r = y - dy
+ y2r = y + dy
+ if ((y2r < 0.5) || (y1r > (PL_AXLEN(pl,2) + 0.5)))
+ return
+ x1r = x - dx
+ x2r = x + dx
+ if ((x2r < 0.5) || (x1r > (PL_AXLEN(pl,1) + 0.5)))
+ return
+
+ call smark (sp)
+ call salloc (ufd, LEN_EANNDES, TY_STRUCT)
+
+ EA_PL(ufd) = pl
+ EA_XCEN(ufd) = x
+ EA_YCEN(ufd) = y
+ EA_AA2(ufd) = aa / r2
+ EA_BB2(ufd) = bb / r2
+ EA_CC2(ufd) = cc / r2
+ EA_FF2(ufd) = ff / r2
+ EA_DXMAX2(ufd) = dx
+ EA_DYMAX2(ufd) = dy
+ EA_PV(ufd) = 1
+
+ # Get the inner ellipse parameters.
+ call me_ellgeom (radius1, ratio, theta, aa, bb, cc, ff)
+ r2 = radius1 * radius1
+ dx = ff / (aa - bb * bb / 4.0 / cc)
+ if (dx > 0.0)
+ dx = sqrt (dx)
+ else
+ dx = 0.0
+ dy = ff / (cc - bb * bb / 4.0 / aa)
+ if (dy > 0.0)
+ dy = sqrt (dy)
+ else
+ dy = 0.0
+
+ EA_AA1(ufd) = aa / r2
+ EA_BB1(ufd) = bb / r2
+ EA_CC1(ufd) = cc / r2
+ EA_FF1(ufd) = ff / r2
+ EA_DXMAX1(ufd) = dx
+ EA_DYMAX1(ufd) = dy
+
+ y1 = max ( 1, min (PL_AXLEN(pl,2), int(y1r)))
+ y2 = max (y1, min (PL_AXLEN(pl,2), int(y2r)))
+ call pl_regionrop (pl, pe_ueannulus, ufd, y1, y2, rop)
+
+ call sfree (sp)
+end
+
+
+# PE_RANNULUS -- Perform a rasterop operation on the area enclosed by a
+# rectangular annulus drawn in a 2-dimensional plane of a mask. If the
+# dimensionality of the mask exceeds 2, the pl_setplane() procedure should be
+# called first to define the plane of the mask to be modified.
+
+procedure pe_rannulus (pl, x, y, radius1, radius2, ratio, theta, rop)
+
+pointer pl #I mask descriptor
+real x, y #I the center of the rectangular annulus
+real radius1, radius2 #I inner and outer semi-major axes
+real ratio #I ratio of the semi-minor / semi-major axes
+real theta #I position angle
+int rop #I rasterop defining operation
+
+real line_1r, line_2r
+pointer sp, ufd, ixp, iyp, oxp, oyp
+int line_1, line_2, i
+extern pe_uarect()
+errchk plvalid
+
+begin
+ # Note sure why this is called.
+ #call plvalid (pl)
+
+ # Initialize the
+ call smark (sp)
+ call salloc (ufd, LEN_RANNDES, TY_STRUCT)
+ call salloc (ixp, 5, TY_REAL)
+ call salloc (iyp, 5, TY_REAL)
+ call salloc (oxp, 5, TY_REAL)
+ call salloc (oyp, 5, TY_REAL)
+
+ # Copy and close the inner polygon.
+ call me_rectgeom (radius1, ratio, theta, Memr[ixp], Memr[iyp])
+ do i = 1, 4 {
+ Memr[ixp+i-1] = Memr[ixp+i-1] + x
+ Memr[iyp+i-1] = Memr[iyp+i-1] + y
+ }
+ Memr[ixp+4] = Memr[ixp]
+ Memr[iyp+4] = Memr[iyp]
+
+ # Create and close the outer polygon.
+ call me_rectgeom (radius2, ratio, theta, Memr[oxp], Memr[oyp])
+ do i = 1, 4 {
+ Memr[oxp+i-1] = Memr[oxp+i-1] + x
+ Memr[oyp+i-1] = Memr[oyp+i-1] + y
+ }
+ Memr[oxp+4] = Memr[oxp]
+ Memr[oyp+4] = Memr[oyp]
+
+ # Compute the range in X in which the polygon should be drawn
+ # and reject polygons that are off the image.
+ call alimr (Memr[oxp], 4, line_1r, line_2r)
+ line_1 = max (1, min (PL_AXLEN(pl,1), int (line_1r)))
+ line_2 = max (line_1, min (PL_AXLEN(pl,1), int (line_2r)))
+ if (line_2 < 1 || line_1 > PL_AXLEN(pl,1)) {
+ call sfree (sp)
+ return
+ }
+
+ # Compute the range in Y in which the polygon should be drawn
+ # and reject polygons that are off the image.
+ call alimr (Memr[oyp], 4, line_1r, line_2r)
+ line_1 = max (1, min (PL_AXLEN(pl,2), int (line_1r)))
+ line_2 = max (line_1, min (PL_AXLEN(pl,2), int (line_2r)))
+ if (line_2 < 1 || line_1 > PL_AXLEN(pl,2)) {
+ call sfree (sp)
+ return
+ }
+
+ # Initialize the region descriptor.
+ RA_PL(ufd) = pl
+ RA_IXP(ufd) = ixp
+ RA_IYP(ufd) = iyp
+ RA_OXP(ufd) = oxp
+ RA_OYP(ufd) = oyp
+ RA_NVER(ufd) = 4
+ RA_PV(ufd) = 1
+
+ call pl_regionrop (pl, pe_uarect, ufd, line_1, line_2, rop)
+
+ call sfree (sp)
+end
+
+
+# PE_APOLYGON -- Perform a rasterop operation on the area enclosed by a
+# polygonal annulus drawn in a 2-dimensional plane of a mask. If the
+# dimensionality of the mask exceeds 2, the pl_setplane() procedure should be
+# called first to define the plane of the mask to be modified.
+
+procedure pe_apolygon (pl, width, x, y, npts, rop)
+
+pointer pl #I mask descriptor
+real width #I width of the polygonal annulus
+real x[npts] #I the inner polygon x-vertices
+real y[npts] #I outer polygon y-vertices
+int npts #I number of points in polygon
+int rop #I rasterop defining operation
+
+real line_1r, line_2r
+pointer sp, ufd, ixp, iyp, oxp, oyp
+int line_1, line_2, i
+extern pe_uapolygon()
+errchk plvalid
+
+begin
+ # Note sure why this is called.
+ #call plvalid (pl)
+ if (npts < 3)
+ return
+
+ # Initialize the
+ call smark (sp)
+ call salloc (ufd, LEN_PAGONDES, TY_STRUCT)
+ call salloc (ixp, npts + 1, TY_REAL)
+ call salloc (iyp, npts + 1, TY_REAL)
+ call salloc (oxp, npts + 1, TY_REAL)
+ call salloc (oyp, npts + 1, TY_REAL)
+
+ # Copy and close the inner polygon.
+ do i = 1, npts {
+ Memr[ixp+i-1] = x[i]
+ Memr[iyp+i-1] = y[i]
+ }
+ Memr[ixp+npts] = x[1]
+ Memr[iyp+npts] = y[1]
+
+ # Create and close the outer polygon.
+ call me_pyexpand (Memr[ixp], Memr[iyp], Memr[oxp], Memr[oyp],
+ npts, width)
+ Memr[oxp+npts] = Memr[oxp]
+ Memr[oyp+npts] = Memr[oyp]
+
+ # Compute the range in X in which the polygon should be drawn
+ # and reject polygons that are off the image.
+ call alimr (Memr[oxp], npts, line_1r, line_2r)
+ line_1 = max (1, min (PL_AXLEN(pl,1), int (line_1r)))
+ line_2 = max (line_1, min (PL_AXLEN(pl,1), int (line_2r)))
+ if (line_2 < 1 || line_1 > PL_AXLEN(pl,1)) {
+ call sfree (sp)
+ return
+ }
+
+ # Compute the range in Y in which the polygon should be drawn
+ # and reject polygons that are off the image.
+ call alimr (Memr[oyp], npts, line_1r, line_2r)
+ line_1 = max (1, min (PL_AXLEN(pl,2), int (line_1r)))
+ line_2 = max (line_1, min (PL_AXLEN(pl,2), int (line_2r)))
+ if (line_2 < 1 || line_1 > PL_AXLEN(pl,2)) {
+ call sfree (sp)
+ return
+ }
+
+ # Initialize the region descriptor.
+ PA_PL(ufd) = pl
+ PA_IXP(ufd) = ixp
+ PA_IYP(ufd) = iyp
+ PA_OXP(ufd) = oxp
+ PA_OYP(ufd) = oyp
+ PA_NVER(ufd) = npts
+ PA_PV(ufd) = 1
+
+ call pl_regionrop (pl, pe_uapolygon, ufd, line_1, line_2, rop)
+
+ call sfree (sp)
+end
+
+
+# PE_COLS -- Rasterop between a set of column ranges as source and an existing
+# mask as destination. It is not necessary for the ranges to be inside the
+# mask; if they are outside, the boundary of the line ranges will be
+# clipped to the boundary of the mask. This is a 2-dim operator. If the
+# image dimensionality is greater than two the pl_setplane procedure should
+# be called first to specify the plane to be modified. These routines are
+# a modification of the ones in plio$plcircle. The main difference is
+# that the x, y, radius parameters are initially set to real numbers not
+# integers.
+
+procedure pe_cols (pl, rangestr, rop)
+
+pointer pl #I mask descriptor
+char rangestr[ARB] #I the input ranges string
+int rop #I rasterop
+
+int npts, nvalues, colno, x1, x2, nregions
+pointer sp, ufd, rgptr, lineptr
+int me_decode_ranges(), me_next_number(), me_previous_number(), pl_p2ri()
+extern pe_ucols()
+
+begin
+ # Not sure why we need to call this routine here.
+ #call plvalid (pl)
+ npts = PL_AXLEN(pl,1)
+
+ call smark (sp)
+ call salloc (ufd, LEN_COLSDES, TY_STRUCT)
+ call salloc (rgptr, 3 * MAX_NRANGES + 1, TY_INT)
+ call salloc (lineptr, npts, TY_INT)
+
+ # Decode the ranges string
+ if (me_decode_ranges (rangestr, Memi[rgptr], MAX_NRANGES,
+ nvalues) == ERR) {
+ call sfree (sp)
+ return
+ }
+
+ # Get the column limits.
+ x1 = INDEFI
+ x2 = INDEFI
+ colno = 0
+ if (me_next_number (Memi[rgptr], colno) != EOF)
+ x1 = colno
+ colno = npts + 1
+ if (me_previous_number (Memi[rgptr], colno) != EOF)
+ x2 = colno
+ if (IS_INDEFI(x1) || IS_INDEFI(x2)) {
+ call sfree (sp)
+ return
+ }
+
+ # Set the pixel values.
+ call aclri (Memi[lineptr], npts)
+ colno = 0
+ while (me_next_number (Memi[rgptr], colno) != EOF) {
+ if (colno < 1 || colno > npts)
+ next
+ Memi[lineptr+colno-1] = 1
+ }
+
+ # Convert the pixel list to a ranges list.
+ nregions = pl_p2ri (Memi[lineptr], 1, Memi[rgptr], npts)
+
+ L_PL(ufd) = pl
+ L_RANGES(ufd) = rgptr
+ L_NRANGES(ufd) = nregions
+ L_XS(ufd) = 1
+ L_NPIX(ufd) = npts
+ L_PV(ufd) = 1
+
+ # Call the regions operator.
+ call pl_regionrop (pl, pe_ucols, ufd, 1, PL_AXLEN(pl,2), rop)
+
+ call sfree (sp)
+end
+
+
+# PE_LINES -- Rasterop between a set of line ranges as source and an existing
+# mask as destination. It is not necessary for the ranges to be inside the
+# mask; if they are outside, the boundary of the line ranges will be
+# clipped to the boundary of the mask. This is a 2-dim operator. If the
+# image dimensionality is greater than two the pl_setplane procedure should
+# be called first to specify the plane to be modified. These routines are
+# a modification of the ones in plio$plcircle. The main difference is
+# that the x, y, radius parameters are initially set to real numbers not
+# integers.
+
+procedure pe_lines (pl, rangestr, rop)
+
+pointer pl #I mask descriptor
+char rangestr[ARB] #I the input ranges string
+int rop #I rasterop
+
+int i, y1, y2, nvalues
+pointer sp, rgptr, ufd
+int me_decode_ranges()
+bool me_is_in_range()
+extern pe_ulines()
+
+begin
+ # Not sure why we need to call this routine here.
+ #call plvalid (pl)
+
+ call smark (sp)
+ call salloc (ufd, LEN_LINESDES, TY_STRUCT)
+ call salloc (rgptr, 3 * MAX_NRANGES + 1, TY_INT)
+
+ # Decode the ranges string
+ if (me_decode_ranges (rangestr, Memi[rgptr], MAX_NRANGES,
+ nvalues) == ERR) {
+ call sfree (sp)
+ return
+ }
+
+ # Find the line limits.
+ y1 = INDEFI
+ y2 = INDEFI
+ do i = 1, PL_AXLEN(pl,2) {
+ if (me_is_in_range (Memi[rgptr], i)) {
+ y1 = i
+ break
+ }
+ }
+ if (IS_INDEFI(y1)) {
+ call sfree (sp)
+ return
+ }
+ do i = PL_AXLEN(pl,2), 1, -1 {
+ if (me_is_in_range (Memi[rgptr], i)) {
+ y2 = i
+ break
+ }
+ }
+ if (IS_INDEFI(y2)) {
+ call sfree (sp)
+ return
+ }
+
+ L_PL(ufd) = pl
+ L_RANGES(ufd) = rgptr
+ L_PV(ufd) = 1
+
+ call pl_regionrop (pl, pe_ulines, ufd, y1, y2, rop)
+
+ call sfree (sp)
+end
+
+
+# PE_PIE -- Determine which pixels are inside a pie shaped wedge that
+# intersects the image boundaries.
+
+procedure pe_pie (pl, xc, yc, angle1, angle2, rop)
+
+pointer pl #I the pixel mask descriptor
+real xc, yc #I the center of the wedge
+real angle1, angle2 #I the wedge angles
+int rop #I the mask raster op
+
+real sweep, x2, y2, vx[7], vy[7]
+int count, intrcpt1, intrcpt2
+int me_pie_intercept(), me_corner_vertex()
+
+begin
+ # Set the first vertex
+ vx[1] = xc
+ vy[1] = yc
+ sweep = angle2 - angle1
+
+ # If the sweep is too small to be noticed don't bother.
+ if (abs (sweep) < SMALL_NUMBER) {
+ return
+ }
+ if (sweep < 0.0)
+ sweep = sweep + 360.0
+
+ # Get the second vertext by computing the intersection of the
+ # first ray with the image boundaries.
+ intrcpt1 = me_pie_intercept (PL_AXLEN(pl,1), PL_AXLEN(pl,2), xc, yc,
+ angle1, vx[2], vy[2])
+
+ # Compute the second intercept.
+ intrcpt2 = me_pie_intercept (PL_AXLEN(pl,1), PL_AXLEN(pl,2), xc, yc,
+ angle2, x2, y2)
+
+ # If angles intercept same side and slice is between them, no corners
+ # else, mark corners until reaching side with second angle intercept.
+ count = 3
+ if ((intrcpt1 != intrcpt2) || (sweep > 180.0)) {
+ repeat {
+ intrcpt1 = me_corner_vertex (intrcpt1, PL_AXLEN(pl,1),
+ PL_AXLEN(pl,2), vx[count], vy[count])
+ count = count + 1
+ } until (intrcpt1 == intrcpt2)
+ }
+
+ # Set last vertex.
+ vx[count] = x2
+ vy[count] = y2
+
+ # Fill in the polygon
+ call pe_polygon (pl, vx, vy, count, rop)
+end
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
diff --git a/pkg/proto/maskexpr/t_mskexpr.x b/pkg/proto/maskexpr/t_mskexpr.x
new file mode 100644
index 00000000..9a1aa912
--- /dev/null
+++ b/pkg/proto/maskexpr/t_mskexpr.x
@@ -0,0 +1,286 @@
+include <fset.h>
+include <ctype.h>
+include <imhdr.h>
+
+# T_MSKEXPR -- Create a list of pixel masks using an expression and a list of
+# reference images.
+#
+# The input expression may be any legal EVVEXPR expression which can be
+# converted to a valid integer mask pixel value. The input operands must be one
+# of, i for the reference image, i.keyword for a reference image header
+# keyword, m for the input mask image, m.keyword for the input mask image
+# header keyword a numeric constant, a builtin function, or the pixel operands
+# I, J, K, etc. May be desirable to replace the reference image operand
+# i with $I. This is a detail however.
+#
+# This task uses the get tokens library in images to expand the macros.
+# This library should probably be removed from images and put in xtools
+# for the applications programmers or in the core system as a useful
+# library maybe in fmtio like imexpr. There is a similar if not identical(?)
+# library in qpoe.
+
+procedure t_mskexpr()
+
+pointer expr, st, xexpr, refim, pmim, refmsk
+pointer sp, exprdb, dims, uaxlen, mskname, imname, refname
+int i, ip, op, msklist, imlist, rmsklist, len_exprbuf, fd, nchars, ch
+int undim, npix, depth
+bool verbose
+
+pointer me_getexprdb(), me_expandtext(), immap(), yt_mappm(), me_mkmask()
+long fstatl()
+int imtopenp(), imtlen(), open(), getci(), imtgetim(), ctoi()
+int clgeti(), strmatch(), imaccess()
+bool clgetb(), strne()
+errchk immap(), yt_pmmap()
+
+begin
+ # Get the expression parameter.
+ call malloc (expr, SZ_COMMAND, TY_CHAR)
+ call clgstr ("expr", Memc[expr], SZ_COMMAND)
+
+ # Get the output mask list.
+ msklist = imtopenp ("masks")
+ if (imtlen (msklist) <= 0) {
+ call eprintf ("The output mask list is empty\n")
+ call imtclose (msklist)
+ call mfree (expr, TY_CHAR)
+ return
+ }
+
+ # Get the input reference image list.
+ imlist = imtopenp ("refimages")
+ if (imtlen (imlist) > 0 && imtlen (imlist) != imtlen (msklist)) {
+ call eprintf (
+ "The reference image and output mask lists are not the same size\n")
+ call imtclose (imlist)
+ call imtclose (msklist)
+ call mfree (expr, TY_CHAR)
+ return
+ }
+
+ # Get the input reference mask list.
+ rmsklist = imtopenp ("refmasks")
+ if (imtlen (rmsklist) > 0 && imtlen (rmsklist) != imtlen (msklist)) {
+ call eprintf (
+ "The reference image and output mask lists are not the same size\n")
+ call imtclose (rmsklist)
+ call imtclose (msklist)
+ call imtclose (imlist)
+ call mfree (expr, TY_CHAR)
+ return
+ }
+
+ # Get some working space.
+ call smark (sp)
+ call salloc (exprdb, SZ_FNAME, TY_CHAR)
+ call salloc (dims, SZ_FNAME, TY_CHAR)
+ call salloc (uaxlen, IM_MAXDIM, TY_LONG)
+ call salloc (mskname, SZ_FNAME, TY_CHAR)
+ call salloc (imname, SZ_FNAME, TY_CHAR)
+ call salloc (refname, SZ_FNAME, TY_CHAR)
+
+ # Get remaining parameters,
+ call clgstr ("exprdb", Memc[exprdb], SZ_PATHNAME)
+ call clgstr ("dims", Memc[dims], SZ_FNAME)
+ depth = clgeti ("depth")
+ verbose = clgetb ("verbose")
+
+ # Load the expression database if any.
+ if (strne (Memc[exprdb], "none"))
+ st = me_getexprdb (Memc[exprdb])
+ else
+ st = NULL
+
+ # Get the expression to be evaluated and expand any file inclusions
+ # or macro references.
+ len_exprbuf = SZ_COMMAND
+ if (Memc[expr] == '@') {
+ fd = open (Memc[expr+1], READ_ONLY, TEXT_FILE)
+ nchars = fstatl (fd, F_FILESIZE)
+ if (nchars > len_exprbuf) {
+ len_exprbuf = nchars
+ call realloc (expr, len_exprbuf, TY_CHAR)
+ }
+ for (op = expr; getci(fd, ch) != EOF; op = op + 1) {
+ if (ch == '\n')
+ Memc[op] = ' '
+ else
+ Memc[op] = ch
+ }
+ Memc[op] = EOS
+ call close (fd)
+ }
+ if (st != NULL) {
+ xexpr = me_expandtext (st, Memc[expr])
+ call mfree (expr, TY_CHAR)
+ expr = xexpr
+ }
+ if (verbose) {
+ call printf ("Expr: %s\n")
+ call pargstr (Memc[expr])
+ call flush (STDOUT)
+ }
+
+ # Determine the default dimension and size of the output image. If the
+ # reference image is defined then the dimensions of the reference
+ # image determine the dimensions of the output mask. Otherwise the
+ # default dimensions are used.
+
+ undim = 0
+ call aclrl (Meml[uaxlen], IM_MAXDIM)
+ for (ip = 1; ctoi (Memc[dims], ip, npix) > 0; ) {
+ Meml[uaxlen+undim] = npix
+ undim = undim + 1
+ for (ch = Memc[dims+ip-1]; IS_WHITE(ch) || ch == ',';
+ ch = Memc[dims+ip-1])
+ ip = ip + 1
+ }
+
+ # Loop over the output mask names.
+ while (imtgetim (msklist, Memc[mskname], SZ_FNAME) != EOF) {
+
+ # Add .pl to output mask name.
+ if (strmatch (Memc[mskname], ".pl$") == 0)
+ call strcat (".pl", Memc[mskname], SZ_FNAME)
+
+ # Check whether the output mask already exists.
+ if (imaccess (Memc[mskname], 0) == YES) {
+ if (verbose) {
+ call printf ("Mask %s already exists\n")
+ call pargstr (Memc[mskname])
+ }
+ next
+ }
+
+ # Open the reference image.
+ if (imtlen (imlist) > 0) {
+ if (imtgetim (imlist, Memc[imname], SZ_FNAME) != EOF) {
+ iferr (refim = immap (Memc[imname], READ_ONLY, 0)) {
+ refim = NULL
+ call printf (
+ "Cannot open reference image %s for mask %s\n")
+ call pargstr (Memc[imname])
+ call pargstr (Memc[mskname])
+ next
+ }
+ } else {
+ refim = NULL
+ call printf ("Cannot open reference image for mask %s\n")
+ call pargstr (Memc[mskname])
+ next
+ }
+ } else
+ refim = NULL
+
+ # Open the reference mask.
+ if (imtlen (rmsklist) > 0) {
+ if (imtgetim (rmsklist, Memc[refname], SZ_FNAME) != EOF) {
+ if (refim != NULL) {
+ iferr (refmsk = yt_mappm (Memc[refname], refim,
+ "logical", Memc[refname], SZ_FNAME))
+ refmsk = NULL
+ } else {
+ iferr (refmsk = immap (Memc[refname], READ_ONLY, 0))
+ refmsk = NULL
+ }
+ if (refmsk == NULL) {
+ call printf (
+ "Cannot open reference mask %s for mask %s\n")
+ call pargstr (Memc[refname])
+ call pargstr (Memc[mskname])
+ if (refim != NULL)
+ call imunmap (refim)
+ next
+ } else if (refim != NULL) {
+ if (IM_NDIM(refim) != IM_NDIM(refmsk)) {
+ call printf (
+ "Reference image and mask for %s don't match\n")
+ call pargstr (Memc[mskname])
+ call imunmap (refmsk)
+ if (refim != NULL)
+ call imunmap (refim)
+ next
+ } else {
+ do i = 1, IM_NDIM(refim) {
+ if (IM_LEN(refim,i) == IM_LEN(refmsk,i))
+ next
+ else
+ break
+ }
+ if (i <= IM_NDIM(refim)) {
+ call printf (
+ "Reference image and mask for %s don't match\n")
+ call pargstr (Memc[mskname])
+ call imunmap (refmsk)
+ if (refim != NULL)
+ call imunmap (refim)
+ next
+ }
+ }
+ }
+ } else {
+ refmsk = NULL
+ call printf ("Cannot open reference mask for mask %s\n")
+ call pargstr (Memc[refname])
+ if (refim != NULL)
+ call imunmap (refim)
+ next
+ }
+ } else
+ refmsk = NULL
+
+ if (verbose) {
+ if (refim != NULL && refmsk != NULL) {
+ call printf ("Creating mask %s\n")
+ call pargstr (Memc[mskname])
+ call printf (" Using reference image %s and mask %s\n")
+ call pargstr (Memc[imname])
+ call pargstr (Memc[refname])
+ } else if (refim != NULL) {
+ call printf ("Creating mask %s using reference image %s\n")
+ call pargstr (Memc[mskname])
+ call pargstr (Memc[imname])
+ } else if (refmsk != NULL) {
+ call printf ("Creating mask %s using reference image %s\n")
+ call pargstr (Memc[mskname])
+ call pargstr (Memc[refname])
+ } else {
+ call printf ("Creating mask %s\n")
+ call pargstr (Memc[mskname])
+ }
+ }
+
+ # Evalute the expression return a mask image pointer.
+ if (refim != NULL)
+ pmim = me_mkmask (Memc[expr], Memc[mskname], refim, refmsk,
+ IM_NDIM(refim), IM_LEN(refim,1), depth)
+ else if (refmsk != NULL)
+ pmim = me_mkmask (Memc[expr], Memc[mskname], refim, refmsk,
+ IM_NDIM(refmsk), IM_LEN(refmsk,1), depth)
+ else
+ pmim = me_mkmask (Memc[expr], Memc[mskname], refim, refmsk,
+ undim, Meml[uaxlen], depth)
+
+ # Save the mask.
+ call imunmap (pmim)
+
+ # Close the reference image.
+ if (refim != NULL)
+ call imunmap (refim)
+
+ # Close the reference mask.
+ if (refmsk != NULL)
+ call imunmap (refmsk)
+ }
+
+ # Cleanup.
+ call mfree (expr, TY_CHAR)
+ if (st != NULL)
+ call stclose (st)
+ call imtclose (rmsklist)
+ call imtclose (msklist)
+ call imtclose (imlist)
+ call sfree (sp)
+end
+
diff --git a/pkg/proto/maskexpr/t_mskregions.x b/pkg/proto/maskexpr/t_mskregions.x
new file mode 100644
index 00000000..0313055d
--- /dev/null
+++ b/pkg/proto/maskexpr/t_mskregions.x
@@ -0,0 +1,264 @@
+include <fset.h>
+include <ctype.h>
+include <imhdr.h>
+include <imset.h>
+include <pmset.h>
+
+define RG_NUMOPTIONS "|constant|number|"
+define RG_CONSTANT 1
+define RG_NUMBER 2
+
+# T_MSKREGIONS -- Create or edit a list of pixel masks using regions
+# descriptors and a list of reference images.
+#
+# The regions descriptor may define a single region or a region expression.
+# For example a circle may be defined as a single region, e.g.
+#
+# circle xc yc radius
+#
+# whereas the overlap of two circular regions may be defined as a region
+# expression
+#
+# circle (xc1, yc1, r1) && circle (xc2, yc2, r2)
+#
+# note that brackets are necessary in one case and not the other and can
+# be used to decide whether or not to send the regions descriptor off to
+# the parser as opposed to sending it off to a simple interpreter.
+#
+# The regions input operands must be one of the builtin region functions.
+#
+
+procedure t_mskregions()
+
+pointer sp, exprdb, dims, regnumber, uaxlen, mskname, imname, regfname
+pointer st, refim, pmim, expr, xexpr
+int reglist, msklist, imlist, undim, regval, depth, regfd, pregval
+int ip, npix, ch, regno, pregno
+char lbrackett
+bool verbose, append
+
+pointer pl
+
+pointer me_getexprdb(), immap(), me_expandtext(), pl_create()
+int clpopnu(), imtopenp(), clplen(), imtlen(), clgeti(), ctoi(), clgfil()
+int imtgetim(), imaccess(), strmatch(), imstati(), fscan(), open()
+int strdic(), stridx()
+bool clgetb(), strne()
+data lbrackett /'('/
+errchk immap()
+
+begin
+ # Get the regions file list.
+ reglist = clpopnu ("regions")
+ if (clplen (reglist) <= 0) {
+ call eprintf ("The regions file list is empty\n")
+ call clpcls (reglist)
+ return
+ }
+
+ # Get the output mask list.
+ msklist = imtopenp ("masks")
+ if (imtlen (msklist) <= 0) {
+ call eprintf ("The output mask list is empty\n")
+ call imtclose (msklist)
+ call clpcls (reglist)
+ return
+ } else if (clplen (reglist) > 1 && clplen (reglist) !=
+ imtlen (msklist)) {
+ call eprintf ("The regions and mask list have different sizes\n")
+ call imtclose (msklist)
+ call clpcls (reglist)
+ return
+ }
+
+ # Get the output image list.
+ imlist = imtopenp ("refimages")
+ if (imtlen (imlist) > 0 && imtlen (imlist) != imtlen (msklist)) {
+ call eprintf (
+ "The reference image and mask lists are not the same size\n")
+ call imtclose (imlist)
+ call imtclose (msklist)
+ call clpcls (reglist)
+ return
+ }
+
+ # Get some working space.
+ call smark (sp)
+ call salloc (exprdb, SZ_FNAME, TY_CHAR)
+ call salloc (dims, SZ_FNAME, TY_CHAR)
+ call salloc (regnumber, SZ_FNAME, TY_CHAR)
+ call salloc (uaxlen, IM_MAXDIM, TY_LONG)
+ call salloc (mskname, SZ_FNAME, TY_CHAR)
+ call salloc (imname, SZ_FNAME, TY_CHAR)
+ call salloc (regfname, SZ_FNAME, TY_CHAR)
+
+ # Get remaining parameters,
+ call clgstr ("dims", Memc[dims], SZ_FNAME)
+ call clgstr ("regnumber", Memc[regnumber], SZ_FNAME)
+ regno = strdic (Memc[regnumber], Memc[regnumber], SZ_FNAME,
+ RG_NUMOPTIONS)
+ regval = clgeti ("regval")
+ depth = clgeti ("depth")
+ call clgstr ("exprdb", Memc[exprdb], SZ_PATHNAME)
+ append = clgetb ("append")
+ verbose = clgetb ("verbose")
+
+ # Load the expression database if any.
+ if (strne (Memc[exprdb], "none"))
+ st = me_getexprdb (Memc[exprdb])
+ else
+ st = NULL
+
+ # Determine the default dimension and size of the output image. If the
+ # reference image is defined then the dimensions of the reference
+ # image determine the dimensions of the output mask. Otherwise the
+ # default dimensions are used.
+
+ undim = 0
+ call aclrl (Meml[uaxlen], IM_MAXDIM)
+ for (ip = 1; ctoi (Memc[dims], ip, npix) > 0; ) {
+ Meml[uaxlen+undim] = npix
+ undim = undim + 1
+ for (ch = Memc[dims+ip-1]; IS_WHITE(ch) || ch == ',';
+ ch = Memc[dims+ip-1])
+ ip = ip + 1
+ }
+
+ # Loop over the output mask names.
+ regfd = NULL
+ while (imtgetim (msklist, Memc[mskname], SZ_FNAME) != EOF) {
+
+ # Add .pl to output mask name.
+ if (strmatch (Memc[mskname], ".pl$") == 0)
+ call strcat (".pl", Memc[mskname], SZ_FNAME)
+
+ # Check whether the output mask already exists.
+ if (imaccess (Memc[mskname], 0) == YES) {
+ if (! append) {
+ if (verbose) {
+ call printf ("Mask %s already exists\n")
+ call pargstr (Memc[mskname])
+ }
+ next
+ }
+ }
+
+ # Open the reference image.
+ if (imtlen (imlist) > 0) {
+ if (imtgetim (imlist, Memc[imname], SZ_FNAME) != EOF) {
+ iferr (refim = immap (Memc[imname], READ_ONLY, 0)) {
+ refim = NULL
+ call printf (
+ "Cannot open reference image %s for mask %s\n")
+ call pargstr (Memc[imname])
+ call pargstr (Memc[mskname])
+ next
+ }
+ } else {
+ refim = NULL
+ call printf ("Cannot open reference image for mask %s\n")
+ call pargstr (Memc[mskname])
+ next
+ }
+ } else
+ refim = NULL
+
+ # Open the output mask.
+ if (imaccess (Memc[mskname], 0) == YES) {
+ pmim = immap (Memc[mskname], READ_WRITE, 0)
+ } else {
+ if (refim != NULL) {
+ pmim = immap (Memc[mskname], NEW_COPY, refim)
+ } else {
+ pmim = immap (Memc[mskname], NEW_IMAGE, 0)
+ IM_NDIM(pmim) = undim
+ call amovl (Meml[uaxlen], IM_LEN(pmim,1), undim)
+ }
+ IM_PIXTYPE(pmim) = TY_INT
+ pl = imstati (pmim, IM_PLDES)
+ call pl_close (pl)
+ #pl = pl_create (undim, Meml[uaxlen], depth)
+ pl = pl_create (IM_NDIM(pmim), IM_LEN(pmim,1), depth)
+ call imseti (pmim, IM_PLDES, pl)
+ call imunmap (pmim)
+ pmim = immap (Memc[mskname], READ_WRITE, 0)
+ }
+
+ # Open the regions list.
+ if (clgfil (reglist, Memc[regfname], SZ_FNAME) != EOF) {
+ if (regfd != NULL)
+ call close (regfd)
+ regfd = open (Memc[regfname], READ_ONLY, TEXT_FILE)
+ } else if (regfd != NULL)
+ call seek (regfd, BOF)
+
+ # Print a header banner.
+ if (verbose) {
+ if (refim == NULL) {
+ call printf ("Creating mask %s\n")
+ call pargstr (Memc[mskname])
+ } else {
+ call printf ("Creating mask %s using reference image %s\n")
+ call pargstr (Memc[mskname])
+ call pargstr (Memc[imname])
+ }
+ call printf (" Using regions file %s\n")
+ call pargstr (Memc[regfname])
+ }
+
+ # Loop over the regions file.
+ pregval = regval
+ pregno = 1
+ while (fscan (regfd) != EOF) {
+
+ # Get the expression.
+ call malloc (expr, SZ_LINE, TY_CHAR)
+ call gargstr (Memc[expr], SZ_LINE)
+
+ # Determine whether or not the region specificationis an
+ # expression or a region description. If the string is
+ # an expression expand it as necessary.
+ if (stridx (lbrackett, Memc[expr]) > 0) {
+ if (st != NULL) {
+ xexpr = me_expandtext (st, Memc[expr])
+ call mfree (expr, TY_CHAR)
+ expr = xexpr
+ }
+ call me_setexpr (Memc[expr], pmim, pregno, pregval, verbose)
+ } else {
+ call me_setreg (Memc[expr], pmim, pregno, pregval, verbose)
+ }
+
+ # Increment the region number if appropriate.
+ pregno = pregno + 1
+ if (regno == RG_NUMBER)
+ pregval = pregval + 1
+
+ call mfree (expr, TY_CHAR)
+ }
+
+ # Save the output mask.
+ call imunmap (pmim)
+
+ # Close the reference image.
+ if (refim != NULL)
+ call imunmap (refim)
+
+ }
+
+ # Close the last regions file.
+ if (regfd != NULL)
+ call close (regfd)
+
+ # Close the expression database symbol table.
+ if (st != NULL)
+ call stclose (st)
+
+ # Close the various image and file lists.
+ call imtclose (imlist)
+ call imtclose (msklist)
+ call clpcls (reglist)
+
+ call sfree (sp)
+end
+
generated by cgit (git 2.34.1) at 2025-09-20 10:21:59 -0400