From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- pkg/proto/maskexpr/gettok.h | 22 + pkg/proto/maskexpr/gettok.x | 922 +++++++++++++++++++++++ pkg/proto/maskexpr/megeom.x | 72 ++ pkg/proto/maskexpr/megsym.x | 31 + pkg/proto/maskexpr/memkmask.x | 839 +++++++++++++++++++++ pkg/proto/maskexpr/meregfuncs.x | 1449 +++++++++++++++++++++++++++++++++++++ pkg/proto/maskexpr/meregmask.x | 753 +++++++++++++++++++ pkg/proto/maskexpr/mesetexpr.x | 36 + pkg/proto/maskexpr/mesetreg.x | 292 ++++++++ pkg/proto/maskexpr/mkpkg | 26 + pkg/proto/maskexpr/mskexpand.x | 261 +++++++ pkg/proto/maskexpr/peregfuncs.h | 131 ++++ pkg/proto/maskexpr/peregfuncs.x | 877 ++++++++++++++++++++++ pkg/proto/maskexpr/peregufcn.x | 808 +++++++++++++++++++++ pkg/proto/maskexpr/t_mskexpr.x | 286 ++++++++ pkg/proto/maskexpr/t_mskregions.x | 264 +++++++ 16 files changed, 7069 insertions(+) create mode 100644 pkg/proto/maskexpr/gettok.h create mode 100644 pkg/proto/maskexpr/gettok.x create mode 100644 pkg/proto/maskexpr/megeom.x create mode 100644 pkg/proto/maskexpr/megsym.x create mode 100644 pkg/proto/maskexpr/memkmask.x create mode 100644 pkg/proto/maskexpr/meregfuncs.x create mode 100644 pkg/proto/maskexpr/meregmask.x create mode 100644 pkg/proto/maskexpr/mesetexpr.x create mode 100644 pkg/proto/maskexpr/mesetreg.x create mode 100644 pkg/proto/maskexpr/mkpkg create mode 100644 pkg/proto/maskexpr/mskexpand.x create mode 100644 pkg/proto/maskexpr/peregfuncs.h create mode 100644 pkg/proto/maskexpr/peregfuncs.x create mode 100644 pkg/proto/maskexpr/peregufcn.x create mode 100644 pkg/proto/maskexpr/t_mskexpr.x create mode 100644 pkg/proto/maskexpr/t_mskregions.x (limited to 'pkg/proto/maskexpr') 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 +include +include +include +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 + +# 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 +include +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 +include +include +include +include +include + +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 +include +include + + +# 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 +include +include +include +include +include + +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 +include + +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 + memkmask.x \ + + + t_mskregions.x + mesetreg.x + mesetexpr.x + meregmask.x \ + + peregfuncs.x peregfuncs.h + peregufcn.x peregfuncs.h + megeom.x + + meregfuncs.x + mskexpand.x gettok.h + megsym.x gettok.h + gettok.x 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 +include +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: +# +# ['(' 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) +include +include +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 +include +include +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 +include +include + +# 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 +include +include +include +include + +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 + -- cgit