aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/imutil/src/t_imsum.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/imutil/src/t_imsum.x')
-rw-r--r--pkg/images/imutil/src/t_imsum.x320
1 files changed, 320 insertions, 0 deletions
diff --git a/pkg/images/imutil/src/t_imsum.x b/pkg/images/imutil/src/t_imsum.x
new file mode 100644
index 00000000..6e4d0c61
--- /dev/null
+++ b/pkg/images/imutil/src/t_imsum.x
@@ -0,0 +1,320 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+
+# IMSUM -- Sum or average images with optional high and low pixel rejection.
+
+procedure t_imsum ()
+
+int list # Input image list
+pointer image # Output image
+pointer hparams # Header parameter list
+pointer option # Output option
+int pixtype # Output pixel datatype
+int calctype # Internal calculation type
+real low_reject # Number or frac of low pix to reject
+real high_reject # Number or frac of high pix to reject
+
+int i, nimages, nlow, nhigh
+pointer sp, str, im_in, im_out
+
+bool clgetb(), streq()
+real clgetr()
+int imtopenp(), imtlen(), imtgetim(), clgwrd()
+pointer immap()
+
+errchk imsum_set, immap, imunmap
+
+begin
+ # Get the input image list. Check that there is at least 1 image.
+ list = imtopenp ("input")
+ nimages = imtlen (list)
+ if (nimages < 1) {
+ call imtclose (list)
+ call error (0, "No input images in list")
+ }
+
+ # Allocate strings and get the parameters.
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (hparams, SZ_LINE, TY_CHAR)
+ call salloc (option, SZ_LINE, TY_CHAR)
+
+ i = clgwrd ("option", Memc[option], SZ_LINE, "|sum|average|median|")
+ if (streq (Memc[option], "median")) {
+ nlow = nimages / 2
+ nhigh = nimages - nlow - 1
+ } else {
+ # If the rejection value is less than 1 then it is a fraction of the
+ # input images otherwise it is the number of pixels to be rejected.
+ low_reject = clgetr ("low_reject")
+ high_reject = clgetr ("high_reject")
+
+ if (low_reject < 1.)
+ nlow = low_reject * nimages
+ else
+ nlow = low_reject
+
+ if (high_reject < 1.)
+ nhigh = high_reject * nimages
+ else
+ nhigh = high_reject
+
+ if (nlow + nhigh >= nimages) {
+ call sfree (sp)
+ call imtclose (list)
+ call error (0, "Number of pixels rejected >= number of images")
+ }
+ }
+ call clgstr ("hparams", Memc[hparams], SZ_LINE)
+
+ # Map the output image and set the title and pixel type.
+ # Check all images have the same number and length of dimensions.
+
+ call imsum_set (list, pixtype, calctype)
+
+ i = imtgetim (list, Memc[image], SZ_FNAME)
+ im_in = immap (Memc[image], READ_ONLY, 0)
+ call clgstr ("output", Memc[image], SZ_FNAME)
+ im_out = immap (Memc[image], NEW_COPY, im_in)
+ call new_title ("title", im_out)
+ IM_PIXTYPE (im_out) = pixtype
+
+ call imtrew (list)
+
+ # Print verbose info.
+ if (clgetb ("verbose")) {
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call printf ("IMSUM:\n")
+ call printf (" Input images:\n")
+ while (imtgetim (list, Memc[str], SZ_LINE) != EOF) {
+ call printf (" %s\n")
+ call pargstr (Memc[str])
+ }
+ call imtrew (list)
+ call printf (" Output image: %s\n")
+ call pargstr (Memc[image])
+ call printf (" Header parameters: %s\n")
+ call pargstr (Memc[hparams])
+ call printf (" Output pixel datatype: %s\n")
+ call dtstring (pixtype, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ call printf (" Calculation type: %s\n")
+ call dtstring (calctype, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ call printf (" Option: %s\n")
+ call pargstr (Memc[option])
+ call printf (" Low rejection: %d\n High rejection: %d\n")
+ call pargi (nlow)
+ call pargi (nhigh)
+ call flush (STDOUT)
+ }
+
+ # Do the image average. Switch on the calculation type.
+ switch (calctype) {
+ case TY_SHORT:
+ call imsums (list, Memc[image], im_out, nlow, nhigh, Memc[option])
+ case TY_INT:
+ call imsumi (list, Memc[image], im_out, nlow, nhigh, Memc[option])
+ case TY_LONG:
+ call imsuml (list, Memc[image], im_out, nlow, nhigh, Memc[option])
+ case TY_REAL:
+ call imsumr (list, Memc[image], im_out, nlow, nhigh, Memc[option])
+ case TY_DOUBLE:
+ call imsumd (list, Memc[image], im_out, nlow, nhigh, Memc[option])
+ default:
+ call imsumr (list, Memc[image], im_out, nlow, nhigh, Memc[option])
+ }
+ call imunmap (im_out)
+ call imunmap (im_in)
+
+ # Set the header parameters.
+ call imtrew (list)
+ call imsum_hparam (list, Memc[image], Memc[hparams], Memc[option])
+
+ call imtclose (list)
+ call sfree (sp)
+end
+
+# IMSUM_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.
+
+define NTYPES 5
+
+procedure imsum_set (list, pixtype, calctype)
+
+int list # List of input images
+int pixtype # Pixel datatype of output image
+int calctype # Pixel datatype for calculations
+
+int i, j, nimages, max_type
+pointer sp, str, im1, im2
+
+int imtgetim(), imtlen()
+bool xt_imleneq()
+pointer immap()
+errchk immap, imunmap
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Determine maximum precedence datatype.
+ # Also check that the images are the same dimension and size.
+
+ nimages = imtlen (list)
+ j = imtgetim (list, Memc[str], SZ_LINE)
+ im1 = immap (Memc[str], READ_ONLY, 0)
+ max_type = IM_PIXTYPE (im1)
+
+ do i = 2, nimages {
+ j = imtgetim (list, Memc[str], SZ_LINE)
+ im2 = immap (Memc[str], READ_ONLY, 0)
+
+ if ((IM_NDIM(im1) != IM_NDIM(im2)) || !xt_imleneq (im1, im2)) {
+ call imunmap (im1)
+ call imunmap (im2)
+ call error (0, "Images have different dimensions or sizes")
+ }
+
+ switch (IM_PIXTYPE (im2)) {
+ case TY_SHORT:
+ if (max_type == TY_USHORT)
+ max_type = TY_INT
+ case TY_USHORT:
+ if (max_type == TY_SHORT)
+ max_type = TY_INT
+ case TY_INT:
+ if (max_type == TY_USHORT || max_type == TY_SHORT)
+ max_type = IM_PIXTYPE (im2)
+ case TY_LONG:
+ if (max_type == TY_USHORT || max_type == TY_SHORT ||
+ max_type == TY_INT)
+ max_type = IM_PIXTYPE (im2)
+ case TY_REAL:
+ if (max_type != TY_DOUBLE)
+ max_type = IM_PIXTYPE (im2)
+ case TY_DOUBLE:
+ max_type = IM_PIXTYPE (im2)
+ default:
+ }
+ call imunmap (im2)
+ }
+
+ call imunmap (im1)
+ call imtrew (list)
+
+ # Set calculation datatype.
+ call clgstr ("calctype", Memc[str], SZ_LINE)
+ switch (Memc[str]) {
+ case EOS:
+ calctype = max_type
+ case 's':
+ calctype = TY_SHORT
+ case 'i':
+ calctype = TY_INT
+ case 'l':
+ calctype = TY_LONG
+ case 'r':
+ calctype = TY_REAL
+ case 'd':
+ calctype = TY_DOUBLE
+ default:
+ call error (0, "Unrecognized datatype")
+ }
+
+ # Set output pixel datatype.
+ call clgstr ("pixtype", Memc[str], SZ_LINE)
+ switch (Memc[str]) {
+ case EOS:
+ pixtype = calctype
+ case 'u':
+ pixtype = TY_USHORT
+ case 's':
+ pixtype = TY_SHORT
+ case 'i':
+ pixtype = TY_INT
+ case 'l':
+ pixtype = TY_LONG
+ case 'r':
+ pixtype = TY_REAL
+ case 'd':
+ pixtype = TY_DOUBLE
+ default:
+ call error (0, "Unrecognized datatype")
+ }
+
+ call sfree (sp)
+end
+
+# IMSUM_HPARM -- Arithmetic on image header parameters.
+#
+# This program is limited by a lack of a rewind procedure for the image
+# header fields list. Thus, a static array of field names is used
+# to require only one pass through the list and the images.
+
+define NFIELDS 10 # Maximum number of fields allowed.
+
+procedure imsum_hparam (list, output, hparams, option)
+
+int list # List of input images.
+char output[ARB] # Output image
+char hparams[ARB] # List of header parameters
+char option[ARB] # Sum option
+
+int i, nfields, flist
+pointer sp, field, dvals, image, in, out
+
+int imofnlu(), imgnfn(), imtgetim(), imtlen()
+bool strne(), streq()
+double imgetd()
+pointer immap()
+
+errchk immap, imofnlu, imgetd, imputd, imunmap
+
+begin
+ # Return if median.
+ if (strne (option, "average") && strne (option, "sum"))
+ return
+
+ # Allocate memory.
+ call smark (sp)
+ call salloc (field, NFIELDS*SZ_FNAME, TY_CHAR)
+ call salloc (dvals, NFIELDS, TY_DOUBLE)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+
+ # Map the fields.
+ out = immap (output, READ_WRITE, 0)
+ flist = imofnlu (out, hparams)
+ i = 0
+ while ((i < NFIELDS) &&
+ (imgnfn (flist, Memc[field+i*SZ_FNAME], SZ_FNAME) != EOF))
+ i = i + 1
+ call imcfnl (flist)
+
+ # Accumulate values from each image.
+
+ nfields = i
+ call aclrd (Memd[dvals], nfields)
+
+ while (imtgetim (list, Memc[image], SZ_FNAME) != EOF) {
+ in = immap (Memc[image], READ_ONLY, 0)
+ do i = 1, nfields
+ Memd[dvals+i-1] = Memd[dvals+i-1] +
+ imgetd (in, Memc[field+(i-1)*SZ_FNAME])
+ call imunmap (in)
+ }
+
+ # Output the sums or average.
+ if (streq (option, "average")) {
+ i = imtlen (list)
+ call adivkd (Memd[dvals], double (i), Memd[dvals], nfields)
+ }
+
+ do i = 1, nfields
+ call imputd (out, Memc[field+(i-1)*SZ_FNAME], Memd[dvals+i-1])
+
+ call imunmap (out)
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 10:48:21 -0400