Initial commit

1 files changed, 489 insertions, 0 deletions

diff --git a/pkg/images/imutil/src/t_imarith.x b/pkg/images/imutil/src/t_imarith.x
new file mode 100644
index 00000000..6d5f6105
--- /dev/null
+++ b/pkg/images/imutil/src/t_imarith.x
@@ -0,0 +1,489 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<error.h>
+include	<lexnum.h>
+
+define	ADD	1			# Opcodes.
+define	SUB	2
+define	MUL	3
+define	DIV	4
+define	MIN	5
+define	MAX	6
+
+# T_IMARITH -- Simple image arithmetic.
+#
+# For each pixel in each image compute:
+#
+#	operand1 op operand2 = result
+#
+# Do the operations as efficiently as possible.  Allow operand1 or operand2
+# to be a constant.  Allow resultant image to have the same name as an
+# operand image.  Allow lists for the operands and the results.
+# Allow one of the operands to have extra dimensions but require that the
+# common dimensions are of the same length.
+
+procedure t_imarith ()
+
+int	list1				# Operand1 list
+int	list2				# Operand2 list
+int	list3				# Result list
+int	op				# Operator
+bool	verbose				# Verbose option
+bool	noact				# Noact option
+double	c1				# Constant for operand1
+double	c2				# Constant for operand2
+double	divzero				# Zero divide replacement
+int	pixtype				# Output pixel datatype
+int	calctype			# Datatype for calculations
+
+int	i, j, pixtype1, pixtype2
+short	sc1, sc2, sdz
+int	hlist
+double	dval1, dval2
+pointer	im1, im2, im3
+pointer	sp, operand1, operand2, result, imtemp
+pointer	opstr, dtstr, field, title, hparams
+
+int	imtopenp(), imtgetim(), imtlen(), imofnlu(), imgnfn()
+double	clgetd(), imgetd()
+bool	clgetb(), streq()
+int	clgwrd()
+int	gctod(), lexnum()
+pointer	immap()
+errchk	immap, imgetd, imputd
+
+begin
+	# Allocate memory for strings.
+	call smark (sp)
+	call salloc (operand1, SZ_FNAME, TY_CHAR)
+	call salloc (operand2, SZ_FNAME, TY_CHAR)
+	call salloc (result, SZ_FNAME, TY_CHAR)
+	call salloc (imtemp, SZ_FNAME, TY_CHAR)
+	call salloc (opstr, SZ_FNAME, TY_CHAR)
+	call salloc (dtstr, SZ_FNAME, TY_CHAR)
+	call salloc (field, SZ_FNAME, TY_CHAR)
+	call salloc (title, SZ_IMTITLE, TY_CHAR)
+	call salloc (hparams, SZ_LINE, TY_CHAR)
+
+	# Get the operands and the operator.
+	list1 = imtopenp ("operand1")
+	op = clgwrd ("op", Memc[opstr], SZ_FNAME, ",+,-,*,/,min,max,")
+	list2 = imtopenp ("operand2")
+	list3 = imtopenp ("result")
+
+	# Get the rest of the options.
+	call clgstr ("hparams", Memc[hparams], SZ_LINE)
+	verbose = clgetb ("verbose")
+	noact = clgetb ("noact")
+	if (op == DIV)
+	    divzero = clgetd ("divzero")
+
+	# Check the number of elements.
+	if (((imtlen (list1) != 1) && (imtlen (list1) != imtlen (list3))) ||
+	    ((imtlen (list2) != 1) && (imtlen (list2) != imtlen (list3)))) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call imtclose (list3)
+	    call error (1, "Wrong number of elements in the operand lists")
+	}
+
+	# Do each operation.
+	while (imtgetim (list3, Memc[result], SZ_FNAME) != EOF) {
+	    if (imtgetim (list1, Memc[imtemp], SZ_FNAME) != EOF)
+		call strcpy (Memc[imtemp], Memc[operand1], SZ_FNAME)
+	    if (imtgetim (list2, Memc[imtemp], SZ_FNAME) != EOF)
+		call strcpy (Memc[imtemp], Memc[operand2], SZ_FNAME)
+
+	    # Image sections in the output are not allowed.
+	    call imgsection (Memc[result], Memc[field], SZ_FNAME)
+	    if (Memc[field] != EOS) {
+		call eprintf (
+	        "imarith: image sections in the output are not allowed (%s)\n")
+		    call pargstr (Memc[result])
+		next
+	    }
+
+	    # To allow purely numeric file names first test if the operand
+	    # is a file.  If it is not then attempt to interpret the operand
+	    # as a numerical constant.  Otherwise it is an error.
+	    iferr {
+		im1 = immap (Memc[operand1], READ_ONLY, 0)
+		pixtype1 = IM_PIXTYPE(im1)
+	    } then {
+		i = 1
+	        j = gctod (Memc[operand1], i, c1)
+		if ((Memc[operand1+i-1]!=EOS) && (Memc[operand1+i-1]!=' ')) {
+		    call eprintf ("%s is not an image or a number\n")
+			call pargstr (Memc[operand1])
+		    next
+		}
+
+		i = 1
+		pixtype1 = lexnum (Memc[operand1], i, j)
+		switch (pixtype1) {
+		case LEX_REAL:
+		    pixtype1 = TY_REAL
+		default:
+		    pixtype1 = TY_SHORT
+		}
+		im1 = NULL
+	    }
+
+	    iferr {
+		im2 = immap (Memc[operand2], READ_ONLY, 0)
+		pixtype2 = IM_PIXTYPE(im2)
+	    } then {
+		i = 1
+	        j = gctod (Memc[operand2], i, c2)
+		if ((Memc[operand2+i-1]!=EOS) && (Memc[operand2+i-1]!=' ')) {
+		    call eprintf ("%s is not an image or a number\n")
+			call pargstr (Memc[operand2])
+		    if (im1 != NULL)
+			call imunmap (im1)
+		    next
+		}
+
+		i = 1
+		pixtype2 = lexnum (Memc[operand2], i, j)
+		switch (pixtype2) {
+		case LEX_REAL:
+		    pixtype2 = TY_REAL
+		default:
+		    pixtype2 = TY_SHORT
+		}
+		im2 = NULL
+	    }
+
+	    # Determine the output pixel datatype and calculation datatype.
+	    call ima_set (pixtype1, pixtype2, op, pixtype, calctype)
+
+	    # If verbose or noact print the operation.
+	    if (verbose || noact) {
+	        call printf ("IMARITH:\n  Operation = %s\n")
+		    call pargstr (Memc[opstr])
+	        call printf ("  Operand1 = %s\n  Operand2 = %s\n")
+		    call pargstr (Memc[operand1])
+		    call pargstr (Memc[operand2])
+	        call printf ("  Result = %s\n  Result pixel type = %s\n")
+		    call pargstr (Memc[result])
+		    call dtstring (pixtype, Memc[dtstr], SZ_FNAME)
+		    call pargstr (Memc[dtstr])
+	        call printf ("  Calculation type = %s\n")
+		    call dtstring (calctype, Memc[dtstr], SZ_FNAME)
+		    call pargstr (Memc[dtstr])
+		if (op == DIV) {
+		    call printf (
+			"  Replacement value for division by zero  = %g\n")
+			call pargd (divzero)
+		}
+	    }
+
+	    # Do the operation if the no act switch is not set.
+	    if (!noact) {
+	        # Check the two operands have the same dimension lengths
+		# over the same dimensions.
+	        if ((im1 != NULL) && (im2 != NULL)) {
+		    j = OK
+		    do i = 1, min (IM_NDIM (im1), IM_NDIM (im2))
+	    	        if (IM_LEN (im1, i) != IM_LEN (im2, i))
+			    j = ERR
+		    if (j == ERR) {
+			call imunmap (im1)
+			call imunmap (im2)
+	                call eprintf (
+			    "Input images have different dimensions\n")
+			next
+		    }
+		}
+
+	        # Create a temporary output image as a copy of one of the
+		# operand images (the one with the highest dimension).
+		# This allows the resultant image to have
+		# the same name as one of the operand images.
+	        if ((im1 != NULL) && (im2 != NULL)) {
+	            call xt_mkimtemp (Memc[operand1], Memc[result],
+			Memc[imtemp], SZ_FNAME)
+		    if (streq (Memc[result], Memc[imtemp]))
+	                call xt_mkimtemp (Memc[operand2], Memc[result],
+			    Memc[imtemp], SZ_FNAME)
+		    if (IM_NDIM(im1) >= IM_NDIM(im2))
+	                im3 = immap (Memc[result], NEW_COPY, im1)
+		    else
+	                im3 = immap (Memc[result], NEW_COPY, im2)
+	        } else if (im1 != NULL) {
+	            call xt_mkimtemp (Memc[operand1], Memc[result],
+			Memc[imtemp], SZ_FNAME)
+	            im3 = immap (Memc[result], NEW_COPY, im1)
+	        } else if (im2 != NULL) {
+	            call xt_mkimtemp (Memc[operand2], Memc[result],
+			Memc[imtemp], SZ_FNAME)
+	            im3 = immap (Memc[result], NEW_COPY, im2)
+	        } else
+	            call error (0, "No operand images")
+
+		# Set the result image title and pixel datatype.
+		call clgstr ("title", Memc[title], SZ_IMTITLE)
+		if (Memc[title] != EOS)
+		    call strcpy (Memc[title], IM_TITLE (im3), SZ_IMTITLE)
+		IM_PIXTYPE (im3) = pixtype
+
+		# Call the appropriate procedure to do the arithmetic
+		# efficiently.
+	        switch (calctype) {
+	        case TY_SHORT:
+		    sc1 = c1
+		    sc2 = c2
+	            switch (op) {
+	            case ADD:
+	                call ima_adds (im1, im2, im3, sc1, sc2)
+	            case SUB:
+	                call ima_subs (im1, im2, im3, sc1, sc2)
+	            case MUL:
+	                call ima_muls (im1, im2, im3, sc1, sc2)
+	            case DIV:
+			sdz = divzero
+	                call ima_divs (im1, im2, im3, sc1, sc2, sdz)
+	            case MIN:
+	                call ima_mins (im1, im2, im3, sc1, sc2)
+	            case MAX:
+	                call ima_maxs (im1, im2, im3, sc1, sc2)
+	            }
+	        case TY_INT:
+	            switch (op) {
+	            case ADD:
+	                call ima_addi (im1, im2, im3, int (c1), int (c2))
+	            case SUB:
+	                call ima_subi (im1, im2, im3, int (c1), int (c2))
+	            case MUL:
+	                call ima_muli (im1, im2, im3, int (c1), int (c2))
+	            case DIV:
+	                call ima_divi (im1, im2, im3, int (c1), int (c2),
+			    int (divzero))
+	            case MIN:
+	                call ima_mini (im1, im2, im3, int (c1), int (c2))
+	            case MAX:
+	                call ima_maxi (im1, im2, im3, int (c1), int (c2))
+	            }
+	        case TY_LONG:
+	            switch (op) {
+	            case ADD:
+	                call ima_addl (im1, im2, im3, long (c1), long (c2))
+	            case SUB:
+	                call ima_subl (im1, im2, im3, long (c1), long (c2))
+	            case MUL:
+	                call ima_mull (im1, im2, im3, long (c1), long (c2))
+	            case DIV:
+	                call ima_divl (im1, im2, im3, long (c1), long (c2),
+			    long (divzero))
+	            case MIN:
+	                call ima_minl (im1, im2, im3, long (c1), long (c2))
+	            case MAX:
+	                call ima_maxl (im1, im2, im3, long (c1), long (c2))
+	            }
+	        case TY_REAL:
+	            switch (op) {
+	            case ADD:
+	                call ima_addr (im1, im2, im3, real (c1), real (c2))
+	            case SUB:
+	                call ima_subr (im1, im2, im3, real (c1), real (c2))
+	            case MUL:
+	                call ima_mulr (im1, im2, im3, real (c1), real (c2))
+	            case DIV:
+	                call ima_divr (im1, im2, im3, real (c1), real (c2),
+			    real (divzero))
+	            case MIN:
+	                call ima_minr (im1, im2, im3, real (c1), real (c2))
+	            case MAX:
+	                call ima_maxr (im1, im2, im3, real (c1), real (c2))
+	            }
+	        case TY_DOUBLE:
+	            switch (op) {
+	            case ADD:
+	                call ima_addd (im1, im2, im3, double(c1), double(c2))
+	            case SUB:
+	                call ima_subd (im1, im2, im3, double(c1), double(c2))
+	            case MUL:
+	                call ima_muld (im1, im2, im3, double(c1), double(c2))
+	            case DIV:
+	                call ima_divd (im1, im2, im3, double(c1), double(c2),
+			    double(divzero))
+	            case MIN:
+	                call ima_mind (im1, im2, im3, double(c1), double(c2))
+	            case MAX:
+	                call ima_maxd (im1, im2, im3, double(c1), double(c2))
+	            }
+		}
+
+	        # Do the header parameters.
+	        iferr {
+		    ifnoerr (dval1 = imgetd (im3, "CCDMEAN"))
+		        call imdelf (im3, "CCDMEAN")
+
+	            hlist = imofnlu (im3, Memc[hparams])
+		    while (imgnfn (hlist, Memc[field], SZ_FNAME) != EOF) {
+	                if (im1 != NULL)
+	                    dval1 = imgetd (im1, Memc[field])
+		        else
+			    dval1 = c1
+	                if (im2 != NULL)
+	                    dval2 = imgetd (im2, Memc[field])
+		        else
+			    dval2 = c2
+
+		        switch (op) {
+	                case ADD:
+			    call imputd (im3, Memc[field], dval1 + dval2)
+	                case SUB:
+			    call imputd (im3, Memc[field], dval1 - dval2)
+	                case MUL:
+			    call imputd (im3, Memc[field], dval1 * dval2)
+	                case DIV:
+			    if (dval2 == 0.) {
+				call eprintf (
+		    "WARNING: Division by zero in header keyword (%s)\n")
+				    call pargstr (Memc[field])
+			    } else
+				call imputd (im3, Memc[field], dval1 / dval2)
+	                case MIN:
+			    call imputd (im3, Memc[field], min (dval1, dval2))
+	                case MAX:
+			    call imputd (im3, Memc[field], max (dval1, dval2))
+		        }
+		    }
+		    call imcfnl (hlist)
+	        } then
+		    call erract (EA_WARN)
+	    }
+
+	    # Unmap images and release the temporary output image.
+	    if (im1 != NULL)
+	        call imunmap (im1)
+	    if (im2 != NULL)
+	        call imunmap (im2)
+	    if (!noact) {
+	        call imunmap (im3)
+	        call xt_delimtemp (Memc[result], Memc[imtemp])
+	    }
+	}
+
+	call imtclose (list1)
+	call imtclose (list2)
+	call imtclose (list3)
+	call sfree (sp)
+end
+
+
+# IMA_SET -- Determine the output image pixel type and the calculation
+# datatype.  The default pixel types are based on the highest arithmetic
+# precendence of the input images or constants.  Division requires
+# a minimum of real.
+
+procedure ima_set (pixtype1, pixtype2, op, pixtype, calctype)
+
+int	pixtype1			# Pixel datatype of operand 1
+int	pixtype2			# Pixel datatype of operand 2
+int	pixtype				# Pixel datatype of resultant image
+int	op				# Operation
+int	calctype			# Pixel datatype for calculations
+
+char	line[1]
+int	max_type
+
+begin
+	# Determine maximum precedence datatype.
+	switch (pixtype1) {
+	case TY_SHORT:
+	    if (op == DIV)
+		max_type = TY_REAL
+	    else if (pixtype2 == TY_USHORT)
+		max_type = TY_LONG
+	    else
+	        max_type = pixtype2
+	case TY_USHORT:
+	    if (op == DIV)
+		max_type = TY_REAL
+	    else if ((pixtype2 == TY_SHORT) || (pixtype2 == TY_USHORT))
+		max_type = TY_LONG
+	    else
+		max_type = pixtype2
+	case TY_INT:
+	    if (op == DIV)
+		max_type = TY_REAL
+	    else if ((pixtype2 == TY_SHORT) || (pixtype2 == TY_USHORT))
+		max_type = pixtype1
+	    else
+		max_type = pixtype2
+	case TY_LONG:
+	    if (op == DIV)
+		max_type = TY_REAL
+	    else if ((pixtype2 == TY_SHORT) || (pixtype2 == TY_USHORT) ||
+	        (pixtype2 == TY_INT))
+		max_type = pixtype1
+	    else
+		max_type = pixtype2
+	case TY_REAL:
+	    if (pixtype2 == TY_DOUBLE)
+		max_type = pixtype2
+	    else
+		max_type = pixtype1
+	case TY_DOUBLE:
+	    max_type = pixtype1
+	}
+
+	# Set calculation datatype.
+	call clgstr ("calctype", line, 1)
+	switch (line[1]) {
+	case '1':
+	    if (pixtype1 == TY_USHORT)
+		calctype = TY_LONG
+	    else
+	        calctype = pixtype1
+	case '2':
+	    if (pixtype2 == TY_USHORT)
+		calctype = TY_LONG
+	    else
+	        calctype = pixtype2
+	case EOS:
+	    calctype = max_type
+	case 's':
+	    calctype = TY_SHORT
+	case 'u':
+	    calctype = TY_LONG
+	case 'i':
+	    calctype = TY_INT
+	case 'l':
+	    calctype = TY_LONG
+	case 'r':
+	    calctype = TY_REAL
+	case 'd':
+	    calctype = TY_DOUBLE
+	default:
+	    call error (6, "Unrecognized datatype")
+	}
+
+	# Set output pixel datatype.
+	call clgstr ("pixtype", line, 1)
+	switch (line[1]) {
+	case '1':
+	    pixtype = pixtype1
+	case '2':
+	    pixtype = pixtype2
+	case EOS:
+	    pixtype = calctype
+	case 's':
+	    pixtype = TY_SHORT
+	case 'u':
+	    pixtype = TY_USHORT
+	case 'i':
+	    pixtype = TY_INT
+	case 'l':
+	    pixtype = TY_LONG
+	case 'r':
+	    pixtype = TY_REAL
+	case 'd':
+	    pixtype = TY_DOUBLE
+	default:
+	    call error (6, "Unrecognized dataype")
+	}
+end