aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/imutil/src/imsum.gx
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/imutil/src/imsum.gx')
-rw-r--r--pkg/images/imutil/src/imsum.gx398
1 files changed, 398 insertions, 0 deletions
diff --git a/pkg/images/imutil/src/imsum.gx b/pkg/images/imutil/src/imsum.gx
new file mode 100644
index 00000000..31afc420
--- /dev/null
+++ b/pkg/images/imutil/src/imsum.gx
@@ -0,0 +1,398 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include "../imsum.h"
+
+define TMINSW 1.00 # Relative timings for nvecs = 5
+define TMXMNSW 1.46
+define TMED3 0.18
+define TMED5 0.55
+
+# IMSUM -- Sum or average images with optional high and low pixel rejection.
+#
+# This procedure has to be clever in not exceeding the maximum number of images
+# which can be mapped at one time. If no pixels are being rejected then the
+# images can be summed (or averaged) in blocks using the output image to hold
+# intermediate results. If pixels are being rejected then lines from all
+# images must be obtained. If the number of images exceeds the maximum
+# then only a subset of the images are kept mapped and the remainder are
+# mapped and unmapped for each line. This, of course, is inefficient but
+# there is no other way.
+
+$for(silrd)
+procedure imsum$t (list, output, im_out, nlow, nhigh, option)
+
+int list # List of input images
+char output[ARB] # Output image
+pointer im_out # Output image pointer
+int nlow # Number of low pixels to reject
+int nhigh # Number of high pixels to reject
+char option[ARB] # Output option
+
+int i, n, nimages, naccept, npix, ndone, pass
+PIXEL const
+pointer sp, input, v1, v2, im, buf, buf1, buf_in, buf_out
+
+bool streq()
+int imtlen(), imtgetim(), imtrgetim()
+pointer immap(), imgnl$t(), impnl$t()
+errchk immap, imunmap, imgnl$t, impnl$t
+
+begin
+ # Initialize.
+ nimages = imtlen (list)
+ naccept = nimages - nlow - nhigh
+ const = naccept
+ npix = IM_LEN(im_out, 1)
+ if (naccept < 1)
+ call error (0, "Number of rejected pixels is too large")
+
+ # Allocate memory.
+ call smark (sp)
+ call salloc (input, SZ_FNAME, TY_CHAR)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (v2, IM_MAXDIM, TY_LONG)
+ call salloc (im, nimages, TY_INT)
+
+ # If there are no pixels to be rejected avoid calls to reject pixels
+ # and do the operation in blocks so that the number of images mapped
+ # does not exceed the maximum. The output image is used to
+ # store intermediate results.
+
+ if ((nlow == 0) && (nhigh == 0)) {
+ pass = 0
+ ndone = 0
+ repeat {
+ n = 0
+ while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+ Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+ n = n + 1
+ if (n == IMS_MAX)
+ break
+ }
+ ndone = ndone + n
+
+ pass = pass + 1
+ if (pass > 1) {
+ call imunmap (im_out)
+ im_out = immap (output, READ_WRITE, 0)
+ }
+
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+ # For each input line compute an output line.
+ while (impnl$t (im_out, buf_out, Meml[v2]) != EOF) {
+
+ # Clear the output buffer during the first pass and
+ # read in the partial sum from the output image during
+ # subsequent passes.
+
+ if (pass == 1)
+ call aclr$t (Mem$t[buf_out], npix)
+ else {
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+ if (imgnl$t (im_out, buf_in, Meml[v2]) == EOF)
+ call error (0, "Error reading input image")
+ call amov$t (Mem$t[buf_in], Mem$t[buf_out], npix)
+ }
+
+ # Accumulate lines from each input image.
+ do i = 1, n {
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+ if (imgnl$t (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+ call error (0, "Error reading input image")
+ call aadd$t (Mem$t[buf_in], Mem$t[buf_out],
+ Mem$t[buf_out], npix)
+ }
+
+ # If all images have been accumulated and averaging then
+ # divide by the number of images.
+ if ((ndone == nimages) && streq (option, "average"))
+ call adivk$t (Mem$t[buf_out], const, Mem$t[buf_out],
+ npix)
+
+ call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+ }
+
+ do i = 1, n
+ call imunmap (Memi[im+i-1])
+ } until (ndone == nimages)
+
+ # Finish up.
+ call sfree (sp)
+ return
+ }
+
+
+ # Map the input images up to the maximum allowed. The remainder
+ # will be mapped during each line.
+ n = 0
+ while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+ Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+ n = n + 1
+ if (n == IMS_MAX - 1)
+ break
+ }
+
+ # Allocate additional buffer space.
+ call salloc (buf, nimages, TY_INT)
+ if (nimages - n > 0)
+ call salloc (buf1, (nimages-n)*npix, TY_PIXEL)
+
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+ # Compute output lines for each input line.
+ while (impnl$t (im_out, buf_out, Meml[v2]) != EOF) {
+
+ # Read lines from the images which remain open.
+ for (i = 1; i <= n; i = i + 1) {
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+ if (imgnl$t (Memi[im+i-1], Memi[buf+i-1], Meml[v2]) == EOF)
+ call error (0, "Error reading input image")
+ }
+
+ # For all additional images map the image, read a line, copy the
+ # data to a buffer since the image buffer is reused, and unmap
+ # the image.
+ for (; i <= nimages; i = i + 1) {
+ if (imtrgetim (list, i, Memc[input], SZ_FNAME) == EOF)
+ break
+ Memi[im+i-1] = immap (Memc[input], READ_ONLY, 0)
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+ if (imgnl$t (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+ call error (0, "Error reading input image")
+ Memi[buf+i-1] = buf1 + (i - n - 1) * npix
+ call amov$t (Mem$t[buf_in], Mem$t[Memi[buf+i-1]], npix)
+ call imunmap (Memi[im+i-1])
+ }
+
+ # Reject pixels.
+ call imrej$t (Memi[buf], nimages, Mem$t[buf_out], npix, nlow, nhigh)
+
+ # If averaging divide the sum by the number of images averaged.
+ if ((naccept > 1) && streq (option, "average")) {
+ const = naccept
+ call adivk$t (Mem$t[buf_out], const, Mem$t[buf_out], npix)
+ }
+
+ call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+ }
+
+ # Finish up.
+ do i = 1, n
+ call imunmap (Memi[im+i-1])
+ call sfree (sp)
+end
+
+
+# IMREJ -- Reject the number of high and low points and sum the rest.
+
+procedure imrej$t (a, nvecs, b, npts, nlow, nhigh)
+
+pointer a[nvecs] # Pointers to set of vectors
+int nvecs # Number of vectors
+PIXEL b[npts] # Output vector
+int npts # Number of points in the vectors
+int nlow # Number of low points to be rejected
+int nhigh # Number of high points to be rejected
+
+int i, j
+int naccept, minrej, npairs, nlow1, nhigh1
+real tmedian, time1, time2
+
+begin
+ naccept = nvecs - nlow - nhigh
+
+ # If no points are rejected return the sum.
+
+ if (naccept == nvecs) {
+ call amov$t (Mem$t[a[1]], b, npts)
+ for (j = 2; j <= naccept; j = j + 1)
+ call aadd$t (Mem$t[a[j]], b, b, npts)
+ return
+ }
+
+ minrej = min (nlow, nhigh)
+ npairs = minrej
+ nlow1 = nlow - npairs
+ nhigh1 = nhigh - npairs
+
+ if ((naccept == 1) && (npairs > 0)) {
+ if (npairs == 1) {
+ tmedian = TMED3
+ npairs = npairs - 1
+ } else {
+ tmedian = TMED5
+ npairs = npairs - 2
+ }
+ } else
+ tmedian = 0
+
+ # Compare the time required to reject the minimum number
+ # of low or high points and extract the number of points to accept
+ # with the time to reject pairs and the excess number of low or
+ # high points to either reach a median of 3 or 5 points or isolate
+ # the acceptable points.
+
+ time1 = TMINSW * (minrej + naccept)
+ time2 = tmedian + TMXMNSW * npairs + TMINSW * (nlow1 + nhigh1)
+
+ i = nvecs
+ if (time1 < time2) {
+
+ # Sort the nlow and naccept points
+ if (nlow < nhigh) {
+ for (j = 1; j <= nlow + naccept; j = j + 1) {
+ call minsw$t (a, i, npts)
+ i = i - 1
+ }
+ call amov$t (Mem$t[a[nhigh+1]], b, npts)
+ for (j = nhigh+2; j <= nhigh+naccept; j = j + 1)
+ call aadd$t (Mem$t[a[j]], b, b, npts)
+
+ # Sort the nhigh and naccept points
+ } else {
+ for (j = 1; j <= nhigh + naccept; j = j + 1) {
+ call maxsw$t (a, i, npts)
+ i = i - 1
+ }
+ call amov$t (Mem$t[a[nlow+1]], b, npts)
+ for (j = nlow+2; j <= nlow+naccept; j = j + 1)
+ call aadd$t (Mem$t[a[j]], b, b, npts)
+ }
+
+ } else {
+ # Reject the npairs low and high points.
+ for (j = 1; j <= npairs; j = j + 1) {
+ call mxmnsw$t (a, i, npts)
+ i = i - 2
+ }
+ # Reject the excess low points.
+ for (j = 1; j <= nlow1; j = j + 1) {
+ call minsw$t (a, i, npts)
+ i = i - 1
+ }
+ # Reject the excess high points.
+ for (j = 1; j <= nhigh1; j = j + 1) {
+ call maxsw$t (a, i, npts)
+ i = i - 1
+ }
+
+ # Check if the remaining points constitute a 3 or 5 point median
+ # or the set of desired points.
+ if (tmedian == 0.) {
+ call amov$t (Mem$t[a[1]], b, npts)
+ for (j = 2; j <= naccept; j = j + 1)
+ call aadd$t (Mem$t[a[j]], b, b, npts)
+ } else if (tmedian == TMED3) {
+ call amed3$t (Mem$t[a[1]], Mem$t[a[2]], Mem$t[a[3]], b, npts)
+ } else {
+ call amed5$t (Mem$t[a[1]], Mem$t[a[2]], Mem$t[a[3]],
+ Mem$t[a[4]], Mem$t[a[5]], b, npts)
+ }
+ }
+end
+
+
+# MINSW -- Given an array of vector pointers for each element in the vectors
+# swap the minimum element with that of the last vector.
+
+procedure minsw$t (a, nvecs, npts)
+
+pointer a[nvecs] # Array of vector pointers
+int nvecs # Number of vectors
+int npts # Number of points in the vectors
+
+int i, j, k, kmin
+PIXEL temp
+
+begin
+ do i = 0, npts - 1 {
+ kmin = a[1] + i
+ do j = 2, nvecs {
+ k = a[j] + i
+ if (Mem$t[k] < Mem$t[kmin])
+ kmin = k
+ }
+ if (k != kmin) {
+ temp = Mem$t[k]
+ Mem$t[k] = Mem$t[kmin]
+ Mem$t[kmin] = temp
+ }
+ }
+end
+
+
+# MAXSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector.
+
+procedure maxsw$t (a, nvecs, npts)
+
+pointer a[nvecs] # Array of vector pointers
+int nvecs # Number of vectors
+int npts # Number of points in the vectors
+
+int i, j, k, kmax
+PIXEL temp
+
+begin
+ do i = 0, npts - 1 {
+ kmax = a[1] + i
+ do j = 2, nvecs {
+ k = a[j] + i
+ if (Mem$t[k] > Mem$t[kmax])
+ kmax = k
+ }
+ if (k != kmax) {
+ temp = Mem$t[k]
+ Mem$t[k] = Mem$t[kmax]
+ Mem$t[kmax] = temp
+ }
+ }
+end
+
+
+# MXMNSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector and the minimum element
+# with that of the next to last vector. The number of vectors must be greater
+# than 1.
+
+procedure mxmnsw$t (a, nvecs, npts)
+
+pointer a[nvecs] # Array of vector pointers
+int nvecs # Number of vectors
+int npts # Number of points in the vectors
+
+int i, j, k, kmax, kmin
+PIXEL temp
+
+begin
+ do i = 0, npts - 1 {
+ kmax = a[1] + i
+ kmin = kmax
+ do j = 2, nvecs {
+ k = a[j] + i
+ if (Mem$t[k] > Mem$t[kmax])
+ kmax = k
+ else if (Mem$t[k] < Mem$t[kmin])
+ kmin = k
+ }
+ temp = Mem$t[k]
+ Mem$t[k] = Mem$t[kmax]
+ Mem$t[kmax] = temp
+ if (kmin == k) {
+ j = a[nvecs - 1] + i
+ temp = Mem$t[j]
+ Mem$t[j] = Mem$t[kmax]
+ Mem$t[kmax] = temp
+ } else {
+ j = a[nvecs - 1] + i
+ temp = Mem$t[j]
+ Mem$t[j] = Mem$t[kmin]
+ Mem$t[kmin] = temp
+ }
+ }
+end
+$endfor
generated by cgit (git 2.34.1) at 2025-09-20 08:55:06 -0400