diff options
Diffstat (limited to 'pkg/images/imutil/src/generic/imsum.x')
| -rw-r--r-- | pkg/images/imutil/src/generic/imsum.x | 1902 |
1 files changed, 1902 insertions, 0 deletions
diff --git a/pkg/images/imutil/src/generic/imsum.x b/pkg/images/imutil/src/generic/imsum.x new file mode 100644 index 00000000..fcb43716 --- /dev/null +++ b/pkg/images/imutil/src/generic/imsum.x @@ -0,0 +1,1902 @@ +# 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. + + +procedure imsums (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 +short const +pointer sp, input, v1, v2, im, buf, buf1, buf_in, buf_out + +bool streq() +int imtlen(), imtgetim(), imtrgetim() +pointer immap(), imgnls(), impnls() +errchk immap, imunmap, imgnls, impnls + +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 (impnls (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 aclrs (Mems[buf_out], npix) + else { + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + if (imgnls (im_out, buf_in, Meml[v2]) == EOF) + call error (0, "Error reading input image") + call amovs (Mems[buf_in], Mems[buf_out], npix) + } + + # Accumulate lines from each input image. + do i = 1, n { + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + if (imgnls (Memi[im+i-1], buf_in, Meml[v2]) == EOF) + call error (0, "Error reading input image") + call aadds (Mems[buf_in], Mems[buf_out], + Mems[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 adivks (Mems[buf_out], const, Mems[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_SHORT) + + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + + # Compute output lines for each input line. + while (impnls (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 (imgnls (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 (imgnls (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 amovs (Mems[buf_in], Mems[Memi[buf+i-1]], npix) + call imunmap (Memi[im+i-1]) + } + + # Reject pixels. + call imrejs (Memi[buf], nimages, Mems[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 adivks (Mems[buf_out], const, Mems[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 imrejs (a, nvecs, b, npts, nlow, nhigh) + +pointer a[nvecs] # Pointers to set of vectors +int nvecs # Number of vectors +short 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 amovs (Mems[a[1]], b, npts) + for (j = 2; j <= naccept; j = j + 1) + call aadds (Mems[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 minsws (a, i, npts) + i = i - 1 + } + call amovs (Mems[a[nhigh+1]], b, npts) + for (j = nhigh+2; j <= nhigh+naccept; j = j + 1) + call aadds (Mems[a[j]], b, b, npts) + + # Sort the nhigh and naccept points + } else { + for (j = 1; j <= nhigh + naccept; j = j + 1) { + call maxsws (a, i, npts) + i = i - 1 + } + call amovs (Mems[a[nlow+1]], b, npts) + for (j = nlow+2; j <= nlow+naccept; j = j + 1) + call aadds (Mems[a[j]], b, b, npts) + } + + } else { + # Reject the npairs low and high points. + for (j = 1; j <= npairs; j = j + 1) { + call mxmnsws (a, i, npts) + i = i - 2 + } + # Reject the excess low points. + for (j = 1; j <= nlow1; j = j + 1) { + call minsws (a, i, npts) + i = i - 1 + } + # Reject the excess high points. + for (j = 1; j <= nhigh1; j = j + 1) { + call maxsws (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 amovs (Mems[a[1]], b, npts) + for (j = 2; j <= naccept; j = j + 1) + call aadds (Mems[a[j]], b, b, npts) + } else if (tmedian == TMED3) { + call amed3s (Mems[a[1]], Mems[a[2]], Mems[a[3]], b, npts) + } else { + call amed5s (Mems[a[1]], Mems[a[2]], Mems[a[3]], + Mems[a[4]], Mems[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 minsws (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 +short temp + +begin + do i = 0, npts - 1 { + kmin = a[1] + i + do j = 2, nvecs { + k = a[j] + i + if (Mems[k] < Mems[kmin]) + kmin = k + } + if (k != kmin) { + temp = Mems[k] + Mems[k] = Mems[kmin] + Mems[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 maxsws (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 +short temp + +begin + do i = 0, npts - 1 { + kmax = a[1] + i + do j = 2, nvecs { + k = a[j] + i + if (Mems[k] > Mems[kmax]) + kmax = k + } + if (k != kmax) { + temp = Mems[k] + Mems[k] = Mems[kmax] + Mems[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 mxmnsws (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 +short temp + +begin + do i = 0, npts - 1 { + kmax = a[1] + i + kmin = kmax + do j = 2, nvecs { + k = a[j] + i + if (Mems[k] > Mems[kmax]) + kmax = k + else if (Mems[k] < Mems[kmin]) + kmin = k + } + temp = Mems[k] + Mems[k] = Mems[kmax] + Mems[kmax] = temp + if (kmin == k) { + j = a[nvecs - 1] + i + temp = Mems[j] + Mems[j] = Mems[kmax] + Mems[kmax] = temp + } else { + j = a[nvecs - 1] + i + temp = Mems[j] + Mems[j] = Mems[kmin] + Mems[kmin] = temp + } + } +end + +procedure imsumi (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 +int const +pointer sp, input, v1, v2, im, buf, buf1, buf_in, buf_out + +bool streq() +int imtlen(), imtgetim(), imtrgetim() +pointer immap(), imgnli(), impnli() +errchk immap, imunmap, imgnli, impnli + +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 (impnli (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 aclri (Memi[buf_out], npix) + else { + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + if (imgnli (im_out, buf_in, Meml[v2]) == EOF) + call error (0, "Error reading input image") + call amovi (Memi[buf_in], Memi[buf_out], npix) + } + + # Accumulate lines from each input image. + do i = 1, n { + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + if (imgnli (Memi[im+i-1], buf_in, Meml[v2]) == EOF) + call error (0, "Error reading input image") + call aaddi (Memi[buf_in], Memi[buf_out], + Memi[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 adivki (Memi[buf_out], const, Memi[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_INT) + + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + + # Compute output lines for each input line. + while (impnli (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 (imgnli (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 (imgnli (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 amovi (Memi[buf_in], Memi[Memi[buf+i-1]], npix) + call imunmap (Memi[im+i-1]) + } + + # Reject pixels. + call imreji (Memi[buf], nimages, Memi[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 adivki (Memi[buf_out], const, Memi[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 imreji (a, nvecs, b, npts, nlow, nhigh) + +pointer a[nvecs] # Pointers to set of vectors +int nvecs # Number of vectors +int 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 amovi (Memi[a[1]], b, npts) + for (j = 2; j <= naccept; j = j + 1) + call aaddi (Memi[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 minswi (a, i, npts) + i = i - 1 + } + call amovi (Memi[a[nhigh+1]], b, npts) + for (j = nhigh+2; j <= nhigh+naccept; j = j + 1) + call aaddi (Memi[a[j]], b, b, npts) + + # Sort the nhigh and naccept points + } else { + for (j = 1; j <= nhigh + naccept; j = j + 1) { + call maxswi (a, i, npts) + i = i - 1 + } + call amovi (Memi[a[nlow+1]], b, npts) + for (j = nlow+2; j <= nlow+naccept; j = j + 1) + call aaddi (Memi[a[j]], b, b, npts) + } + + } else { + # Reject the npairs low and high points. + for (j = 1; j <= npairs; j = j + 1) { + call mxmnswi (a, i, npts) + i = i - 2 + } + # Reject the excess low points. + for (j = 1; j <= nlow1; j = j + 1) { + call minswi (a, i, npts) + i = i - 1 + } + # Reject the excess high points. + for (j = 1; j <= nhigh1; j = j + 1) { + call maxswi (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 amovi (Memi[a[1]], b, npts) + for (j = 2; j <= naccept; j = j + 1) + call aaddi (Memi[a[j]], b, b, npts) + } else if (tmedian == TMED3) { + call amed3i (Memi[a[1]], Memi[a[2]], Memi[a[3]], b, npts) + } else { + call amed5i (Memi[a[1]], Memi[a[2]], Memi[a[3]], + Memi[a[4]], Memi[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 minswi (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 +int temp + +begin + do i = 0, npts - 1 { + kmin = a[1] + i + do j = 2, nvecs { + k = a[j] + i + if (Memi[k] < Memi[kmin]) + kmin = k + } + if (k != kmin) { + temp = Memi[k] + Memi[k] = Memi[kmin] + Memi[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 maxswi (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 +int temp + +begin + do i = 0, npts - 1 { + kmax = a[1] + i + do j = 2, nvecs { + k = a[j] + i + if (Memi[k] > Memi[kmax]) + kmax = k + } + if (k != kmax) { + temp = Memi[k] + Memi[k] = Memi[kmax] + Memi[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 mxmnswi (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 +int temp + +begin + do i = 0, npts - 1 { + kmax = a[1] + i + kmin = kmax + do j = 2, nvecs { + k = a[j] + i + if (Memi[k] > Memi[kmax]) + kmax = k + else if (Memi[k] < Memi[kmin]) + kmin = k + } + temp = Memi[k] + Memi[k] = Memi[kmax] + Memi[kmax] = temp + if (kmin == k) { + j = a[nvecs - 1] + i + temp = Memi[j] + Memi[j] = Memi[kmax] + Memi[kmax] = temp + } else { + j = a[nvecs - 1] + i + temp = Memi[j] + Memi[j] = Memi[kmin] + Memi[kmin] = temp + } + } +end + +procedure imsuml (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 +long const +pointer sp, input, v1, v2, im, buf, buf1, buf_in, buf_out + +bool streq() +int imtlen(), imtgetim(), imtrgetim() +pointer immap(), imgnll(), impnll() +errchk immap, imunmap, imgnll, impnll + +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 (impnll (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 aclrl (Meml[buf_out], npix) + else { + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + if (imgnll (im_out, buf_in, Meml[v2]) == EOF) + call error (0, "Error reading input image") + call amovl (Meml[buf_in], Meml[buf_out], npix) + } + + # Accumulate lines from each input image. + do i = 1, n { + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + if (imgnll (Memi[im+i-1], buf_in, Meml[v2]) == EOF) + call error (0, "Error reading input image") + call aaddl (Meml[buf_in], Meml[buf_out], + Meml[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 adivkl (Meml[buf_out], const, Meml[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_LONG) + + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + + # Compute output lines for each input line. + while (impnll (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 (imgnll (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 (imgnll (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 amovl (Meml[buf_in], Meml[Memi[buf+i-1]], npix) + call imunmap (Memi[im+i-1]) + } + + # Reject pixels. + call imrejl (Memi[buf], nimages, Meml[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 adivkl (Meml[buf_out], const, Meml[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 imrejl (a, nvecs, b, npts, nlow, nhigh) + +pointer a[nvecs] # Pointers to set of vectors +int nvecs # Number of vectors +long 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 amovl (Meml[a[1]], b, npts) + for (j = 2; j <= naccept; j = j + 1) + call aaddl (Meml[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 minswl (a, i, npts) + i = i - 1 + } + call amovl (Meml[a[nhigh+1]], b, npts) + for (j = nhigh+2; j <= nhigh+naccept; j = j + 1) + call aaddl (Meml[a[j]], b, b, npts) + + # Sort the nhigh and naccept points + } else { + for (j = 1; j <= nhigh + naccept; j = j + 1) { + call maxswl (a, i, npts) + i = i - 1 + } + call amovl (Meml[a[nlow+1]], b, npts) + for (j = nlow+2; j <= nlow+naccept; j = j + 1) + call aaddl (Meml[a[j]], b, b, npts) + } + + } else { + # Reject the npairs low and high points. + for (j = 1; j <= npairs; j = j + 1) { + call mxmnswl (a, i, npts) + i = i - 2 + } + # Reject the excess low points. + for (j = 1; j <= nlow1; j = j + 1) { + call minswl (a, i, npts) + i = i - 1 + } + # Reject the excess high points. + for (j = 1; j <= nhigh1; j = j + 1) { + call maxswl (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 amovl (Meml[a[1]], b, npts) + for (j = 2; j <= naccept; j = j + 1) + call aaddl (Meml[a[j]], b, b, npts) + } else if (tmedian == TMED3) { + call amed3l (Meml[a[1]], Meml[a[2]], Meml[a[3]], b, npts) + } else { + call amed5l (Meml[a[1]], Meml[a[2]], Meml[a[3]], + Meml[a[4]], Meml[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 minswl (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 +long temp + +begin + do i = 0, npts - 1 { + kmin = a[1] + i + do j = 2, nvecs { + k = a[j] + i + if (Meml[k] < Meml[kmin]) + kmin = k + } + if (k != kmin) { + temp = Meml[k] + Meml[k] = Meml[kmin] + Meml[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 maxswl (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 +long temp + +begin + do i = 0, npts - 1 { + kmax = a[1] + i + do j = 2, nvecs { + k = a[j] + i + if (Meml[k] > Meml[kmax]) + kmax = k + } + if (k != kmax) { + temp = Meml[k] + Meml[k] = Meml[kmax] + Meml[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 mxmnswl (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 +long temp + +begin + do i = 0, npts - 1 { + kmax = a[1] + i + kmin = kmax + do j = 2, nvecs { + k = a[j] + i + if (Meml[k] > Meml[kmax]) + kmax = k + else if (Meml[k] < Meml[kmin]) + kmin = k + } + temp = Meml[k] + Meml[k] = Meml[kmax] + Meml[kmax] = temp + if (kmin == k) { + j = a[nvecs - 1] + i + temp = Meml[j] + Meml[j] = Meml[kmax] + Meml[kmax] = temp + } else { + j = a[nvecs - 1] + i + temp = Meml[j] + Meml[j] = Meml[kmin] + Meml[kmin] = temp + } + } +end + +procedure imsumr (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 +real const +pointer sp, input, v1, v2, im, buf, buf1, buf_in, buf_out + +bool streq() +int imtlen(), imtgetim(), imtrgetim() +pointer immap(), imgnlr(), impnlr() +errchk immap, imunmap, imgnlr, impnlr + +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 (impnlr (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 aclrr (Memr[buf_out], npix) + else { + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + if (imgnlr (im_out, buf_in, Meml[v2]) == EOF) + call error (0, "Error reading input image") + call amovr (Memr[buf_in], Memr[buf_out], npix) + } + + # Accumulate lines from each input image. + do i = 1, n { + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + if (imgnlr (Memi[im+i-1], buf_in, Meml[v2]) == EOF) + call error (0, "Error reading input image") + call aaddr (Memr[buf_in], Memr[buf_out], + Memr[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 adivkr (Memr[buf_out], const, Memr[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_REAL) + + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + + # Compute output lines for each input line. + while (impnlr (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 (imgnlr (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 (imgnlr (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 amovr (Memr[buf_in], Memr[Memi[buf+i-1]], npix) + call imunmap (Memi[im+i-1]) + } + + # Reject pixels. + call imrejr (Memi[buf], nimages, Memr[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 adivkr (Memr[buf_out], const, Memr[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 imrejr (a, nvecs, b, npts, nlow, nhigh) + +pointer a[nvecs] # Pointers to set of vectors +int nvecs # Number of vectors +real 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 amovr (Memr[a[1]], b, npts) + for (j = 2; j <= naccept; j = j + 1) + call aaddr (Memr[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 minswr (a, i, npts) + i = i - 1 + } + call amovr (Memr[a[nhigh+1]], b, npts) + for (j = nhigh+2; j <= nhigh+naccept; j = j + 1) + call aaddr (Memr[a[j]], b, b, npts) + + # Sort the nhigh and naccept points + } else { + for (j = 1; j <= nhigh + naccept; j = j + 1) { + call maxswr (a, i, npts) + i = i - 1 + } + call amovr (Memr[a[nlow+1]], b, npts) + for (j = nlow+2; j <= nlow+naccept; j = j + 1) + call aaddr (Memr[a[j]], b, b, npts) + } + + } else { + # Reject the npairs low and high points. + for (j = 1; j <= npairs; j = j + 1) { + call mxmnswr (a, i, npts) + i = i - 2 + } + # Reject the excess low points. + for (j = 1; j <= nlow1; j = j + 1) { + call minswr (a, i, npts) + i = i - 1 + } + # Reject the excess high points. + for (j = 1; j <= nhigh1; j = j + 1) { + call maxswr (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 amovr (Memr[a[1]], b, npts) + for (j = 2; j <= naccept; j = j + 1) + call aaddr (Memr[a[j]], b, b, npts) + } else if (tmedian == TMED3) { + call amed3r (Memr[a[1]], Memr[a[2]], Memr[a[3]], b, npts) + } else { + call amed5r (Memr[a[1]], Memr[a[2]], Memr[a[3]], + Memr[a[4]], Memr[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 minswr (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 +real temp + +begin + do i = 0, npts - 1 { + kmin = a[1] + i + do j = 2, nvecs { + k = a[j] + i + if (Memr[k] < Memr[kmin]) + kmin = k + } + if (k != kmin) { + temp = Memr[k] + Memr[k] = Memr[kmin] + Memr[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 maxswr (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 +real temp + +begin + do i = 0, npts - 1 { + kmax = a[1] + i + do j = 2, nvecs { + k = a[j] + i + if (Memr[k] > Memr[kmax]) + kmax = k + } + if (k != kmax) { + temp = Memr[k] + Memr[k] = Memr[kmax] + Memr[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 mxmnswr (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 +real temp + +begin + do i = 0, npts - 1 { + kmax = a[1] + i + kmin = kmax + do j = 2, nvecs { + k = a[j] + i + if (Memr[k] > Memr[kmax]) + kmax = k + else if (Memr[k] < Memr[kmin]) + kmin = k + } + temp = Memr[k] + Memr[k] = Memr[kmax] + Memr[kmax] = temp + if (kmin == k) { + j = a[nvecs - 1] + i + temp = Memr[j] + Memr[j] = Memr[kmax] + Memr[kmax] = temp + } else { + j = a[nvecs - 1] + i + temp = Memr[j] + Memr[j] = Memr[kmin] + Memr[kmin] = temp + } + } +end + +procedure imsumd (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 +double const +pointer sp, input, v1, v2, im, buf, buf1, buf_in, buf_out + +bool streq() +int imtlen(), imtgetim(), imtrgetim() +pointer immap(), imgnld(), impnld() +errchk immap, imunmap, imgnld, impnld + +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 (impnld (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 aclrd (Memd[buf_out], npix) + else { + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + if (imgnld (im_out, buf_in, Meml[v2]) == EOF) + call error (0, "Error reading input image") + call amovd (Memd[buf_in], Memd[buf_out], npix) + } + + # Accumulate lines from each input image. + do i = 1, n { + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + if (imgnld (Memi[im+i-1], buf_in, Meml[v2]) == EOF) + call error (0, "Error reading input image") + call aaddd (Memd[buf_in], Memd[buf_out], + Memd[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 adivkd (Memd[buf_out], const, Memd[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_DOUBLE) + + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovl (Meml[v1], Meml[v2], IM_MAXDIM) + + # Compute output lines for each input line. + while (impnld (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 (imgnld (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 (imgnld (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 amovd (Memd[buf_in], Memd[Memi[buf+i-1]], npix) + call imunmap (Memi[im+i-1]) + } + + # Reject pixels. + call imrejd (Memi[buf], nimages, Memd[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 adivkd (Memd[buf_out], const, Memd[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 imrejd (a, nvecs, b, npts, nlow, nhigh) + +pointer a[nvecs] # Pointers to set of vectors +int nvecs # Number of vectors +double 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 amovd (Memd[a[1]], b, npts) + for (j = 2; j <= naccept; j = j + 1) + call aaddd (Memd[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 minswd (a, i, npts) + i = i - 1 + } + call amovd (Memd[a[nhigh+1]], b, npts) + for (j = nhigh+2; j <= nhigh+naccept; j = j + 1) + call aaddd (Memd[a[j]], b, b, npts) + + # Sort the nhigh and naccept points + } else { + for (j = 1; j <= nhigh + naccept; j = j + 1) { + call maxswd (a, i, npts) + i = i - 1 + } + call amovd (Memd[a[nlow+1]], b, npts) + for (j = nlow+2; j <= nlow+naccept; j = j + 1) + call aaddd (Memd[a[j]], b, b, npts) + } + + } else { + # Reject the npairs low and high points. + for (j = 1; j <= npairs; j = j + 1) { + call mxmnswd (a, i, npts) + i = i - 2 + } + # Reject the excess low points. + for (j = 1; j <= nlow1; j = j + 1) { + call minswd (a, i, npts) + i = i - 1 + } + # Reject the excess high points. + for (j = 1; j <= nhigh1; j = j + 1) { + call maxswd (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 amovd (Memd[a[1]], b, npts) + for (j = 2; j <= naccept; j = j + 1) + call aaddd (Memd[a[j]], b, b, npts) + } else if (tmedian == TMED3) { + call amed3d (Memd[a[1]], Memd[a[2]], Memd[a[3]], b, npts) + } else { + call amed5d (Memd[a[1]], Memd[a[2]], Memd[a[3]], + Memd[a[4]], Memd[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 minswd (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 +double temp + +begin + do i = 0, npts - 1 { + kmin = a[1] + i + do j = 2, nvecs { + k = a[j] + i + if (Memd[k] < Memd[kmin]) + kmin = k + } + if (k != kmin) { + temp = Memd[k] + Memd[k] = Memd[kmin] + Memd[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 maxswd (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 +double temp + +begin + do i = 0, npts - 1 { + kmax = a[1] + i + do j = 2, nvecs { + k = a[j] + i + if (Memd[k] > Memd[kmax]) + kmax = k + } + if (k != kmax) { + temp = Memd[k] + Memd[k] = Memd[kmax] + Memd[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 mxmnswd (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 +double temp + +begin + do i = 0, npts - 1 { + kmax = a[1] + i + kmin = kmax + do j = 2, nvecs { + k = a[j] + i + if (Memd[k] > Memd[kmax]) + kmax = k + else if (Memd[k] < Memd[kmin]) + kmin = k + } + temp = Memd[k] + Memd[k] = Memd[kmax] + Memd[kmax] = temp + if (kmin == k) { + j = a[nvecs - 1] + i + temp = Memd[j] + Memd[j] = Memd[kmax] + Memd[kmax] = temp + } else { + j = a[nvecs - 1] + i + temp = Memd[j] + Memd[j] = Memd[kmin] + Memd[kmin] = temp + } + } +end + |