diff options
Diffstat (limited to 'pkg/images/imutil/src/generic')
| -rw-r--r-- | pkg/images/imutil/src/generic/imaadd.x | 255 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/imadiv.x | 347 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/imamax.x | 212 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/imamin.x | 212 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/imamul.x | 257 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/imanl.x | 159 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/imasub.x | 252 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/imfuncs.x | 1613 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/imjoin.x | 527 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/imrep.x | 1423 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/imsum.x | 1902 | ||||
| -rw-r--r-- | pkg/images/imutil/src/generic/mkpkg | 21 |
12 files changed, 7180 insertions, 0 deletions
diff --git a/pkg/images/imutil/src/generic/imaadd.x b/pkg/images/imutil/src/generic/imaadd.x new file mode 100644 index 00000000..cd492467 --- /dev/null +++ b/pkg/images/imutil/src/generic/imaadd.x @@ -0,0 +1,255 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + + +# IMA_ADD -- Image arithmetic addition. + +procedure ima_adds (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +short a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nls() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do a vector/scalar + # addition to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nls (im, buf, v, 2) != EOF) { + if (a == 0) + call amovs (Mems[buf[2]], Mems[buf[1]], len) + else + call aaddks (Mems[buf[2]], a, Mems[buf[1]], len) + } + + # If imageb is constant then read imagea and do a vector/scalar + # addition to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nls (im, buf, v, 2) != EOF) { + if (b == 0) + call amovs (Mems[buf[2]], Mems[buf[1]], len) + else + call aaddks (Mems[buf[2]], b, Mems[buf[1]], len) + } + + # Read imagea and imageb and do a vector addition into imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nls (im, buf, v, 3) != EOF) + call aadds (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], len) + } +end + +# IMA_ADD -- Image arithmetic addition. + +procedure ima_addi (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +int a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nli() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do a vector/scalar + # addition to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nli (im, buf, v, 2) != EOF) { + if (a == 0) + call amovi (Memi[buf[2]], Memi[buf[1]], len) + else + call aaddki (Memi[buf[2]], a, Memi[buf[1]], len) + } + + # If imageb is constant then read imagea and do a vector/scalar + # addition to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nli (im, buf, v, 2) != EOF) { + if (b == 0) + call amovi (Memi[buf[2]], Memi[buf[1]], len) + else + call aaddki (Memi[buf[2]], b, Memi[buf[1]], len) + } + + # Read imagea and imageb and do a vector addition into imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nli (im, buf, v, 3) != EOF) + call aaddi (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], len) + } +end + +# IMA_ADD -- Image arithmetic addition. + +procedure ima_addl (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +long a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nll() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do a vector/scalar + # addition to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nll (im, buf, v, 2) != EOF) { + if (a == 0) + call amovl (Meml[buf[2]], Meml[buf[1]], len) + else + call aaddkl (Meml[buf[2]], a, Meml[buf[1]], len) + } + + # If imageb is constant then read imagea and do a vector/scalar + # addition to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nll (im, buf, v, 2) != EOF) { + if (b == 0) + call amovl (Meml[buf[2]], Meml[buf[1]], len) + else + call aaddkl (Meml[buf[2]], b, Meml[buf[1]], len) + } + + # Read imagea and imageb and do a vector addition into imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nll (im, buf, v, 3) != EOF) + call aaddl (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], len) + } +end + +# IMA_ADD -- Image arithmetic addition. + +procedure ima_addr (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +real a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nlr() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do a vector/scalar + # addition to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nlr (im, buf, v, 2) != EOF) { + if (a == 0.0) + call amovr (Memr[buf[2]], Memr[buf[1]], len) + else + call aaddkr (Memr[buf[2]], a, Memr[buf[1]], len) + } + + # If imageb is constant then read imagea and do a vector/scalar + # addition to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nlr (im, buf, v, 2) != EOF) { + if (b == 0.0) + call amovr (Memr[buf[2]], Memr[buf[1]], len) + else + call aaddkr (Memr[buf[2]], b, Memr[buf[1]], len) + } + + # Read imagea and imageb and do a vector addition into imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nlr (im, buf, v, 3) != EOF) + call aaddr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], len) + } +end + +# IMA_ADD -- Image arithmetic addition. + +procedure ima_addd (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +double a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nld() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do a vector/scalar + # addition to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nld (im, buf, v, 2) != EOF) { + if (a == 0.0D0) + call amovd (Memd[buf[2]], Memd[buf[1]], len) + else + call aaddkd (Memd[buf[2]], a, Memd[buf[1]], len) + } + + # If imageb is constant then read imagea and do a vector/scalar + # addition to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nld (im, buf, v, 2) != EOF) { + if (b == 0.0D0) + call amovd (Memd[buf[2]], Memd[buf[1]], len) + else + call aaddkd (Memd[buf[2]], b, Memd[buf[1]], len) + } + + # Read imagea and imageb and do a vector addition into imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nld (im, buf, v, 3) != EOF) + call aaddd (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], len) + } +end + diff --git a/pkg/images/imutil/src/generic/imadiv.x b/pkg/images/imutil/src/generic/imadiv.x new file mode 100644 index 00000000..1de8b194 --- /dev/null +++ b/pkg/images/imutil/src/generic/imadiv.x @@ -0,0 +1,347 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +# IMA_DIV -- Image arithmetic division. + + +procedure ima_divs (im_a, im_b, im_c, a, b, c) + +pointer im_a, im_b, im_c +short a, b, c + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nls() +short ima_efncs() +extern ima_efncs + +short divzero +common /imadcoms/ divzero + +begin + # Loop through all of the image lines. + divzero = c + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do a vector + # reciprical to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nls (im, buf, v, 2) != EOF) + call arczs (a, Mems[buf[2]], Mems[buf[1]], len, + ima_efncs) + + # If imageb is constant then read imagea. If the constant + # is 1 do a vector move to imagec otherwise do a vector/scalar + # divide to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nls (im, buf, v, 2) != EOF) { + if (b == 0) + call amovks (divzero, Mems[buf[1]], len) + else if (b == 1) + call amovs (Mems[buf[2]], Mems[buf[1]], len) + else + call adivks (Mems[buf[2]], b, Mems[buf[1]], len) + } + + # Read imagea and imageb and do the vector divide to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nls (im, buf, v, 3) != EOF) + call advzs (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], + len, ima_efncs) + } +end + + +# IMA_EFNC -- Error function for division by zero. + +short procedure ima_efncs (a) + +short a +short divzero +common /imadcoms/ divzero + +begin + return (divzero) +end + +procedure ima_divi (im_a, im_b, im_c, a, b, c) + +pointer im_a, im_b, im_c +int a, b, c + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nli() +int ima_efnci() +extern ima_efnci + +int divzero +common /imadcomi/ divzero + +begin + # Loop through all of the image lines. + divzero = c + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do a vector + # reciprical to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nli (im, buf, v, 2) != EOF) + call arczi (a, Memi[buf[2]], Memi[buf[1]], len, + ima_efnci) + + # If imageb is constant then read imagea. If the constant + # is 1 do a vector move to imagec otherwise do a vector/scalar + # divide to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nli (im, buf, v, 2) != EOF) { + if (b == 0) + call amovki (divzero, Memi[buf[1]], len) + else if (b == 1) + call amovi (Memi[buf[2]], Memi[buf[1]], len) + else + call adivki (Memi[buf[2]], b, Memi[buf[1]], len) + } + + # Read imagea and imageb and do the vector divide to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nli (im, buf, v, 3) != EOF) + call advzi (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], + len, ima_efnci) + } +end + + +# IMA_EFNC -- Error function for division by zero. + +int procedure ima_efnci (a) + +int a +int divzero +common /imadcomi/ divzero + +begin + return (divzero) +end + +procedure ima_divl (im_a, im_b, im_c, a, b, c) + +pointer im_a, im_b, im_c +long a, b, c + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nll() +long ima_efncl() +extern ima_efncl + +long divzero +common /imadcoml/ divzero + +begin + # Loop through all of the image lines. + divzero = c + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do a vector + # reciprical to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nll (im, buf, v, 2) != EOF) + call arczl (a, Meml[buf[2]], Meml[buf[1]], len, + ima_efncl) + + # If imageb is constant then read imagea. If the constant + # is 1 do a vector move to imagec otherwise do a vector/scalar + # divide to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nll (im, buf, v, 2) != EOF) { + if (b == 0) + call amovkl (divzero, Meml[buf[1]], len) + else if (b == 1) + call amovl (Meml[buf[2]], Meml[buf[1]], len) + else + call adivkl (Meml[buf[2]], b, Meml[buf[1]], len) + } + + # Read imagea and imageb and do the vector divide to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nll (im, buf, v, 3) != EOF) + call advzl (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], + len, ima_efncl) + } +end + + +# IMA_EFNC -- Error function for division by zero. + +long procedure ima_efncl (a) + +long a +long divzero +common /imadcoml/ divzero + +begin + return (divzero) +end + +procedure ima_divr (im_a, im_b, im_c, a, b, c) + +pointer im_a, im_b, im_c +real a, b, c + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nlr() +real ima_efncr() +extern ima_efncr + +real divzero +common /imadcomr/ divzero + +begin + # Loop through all of the image lines. + divzero = c + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do a vector + # reciprical to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nlr (im, buf, v, 2) != EOF) + call arczr (a, Memr[buf[2]], Memr[buf[1]], len, + ima_efncr) + + # If imageb is constant then read imagea. If the constant + # is 1 do a vector move to imagec otherwise do a vector/scalar + # divide to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nlr (im, buf, v, 2) != EOF) { + if (b == 0.0) + call amovkr (divzero, Memr[buf[1]], len) + else if (b == 1.0) + call amovr (Memr[buf[2]], Memr[buf[1]], len) + else + call adivkr (Memr[buf[2]], b, Memr[buf[1]], len) + } + + # Read imagea and imageb and do the vector divide to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nlr (im, buf, v, 3) != EOF) + call advzr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], + len, ima_efncr) + } +end + + +# IMA_EFNC -- Error function for division by zero. + +real procedure ima_efncr (a) + +real a +real divzero +common /imadcomr/ divzero + +begin + return (divzero) +end + +procedure ima_divd (im_a, im_b, im_c, a, b, c) + +pointer im_a, im_b, im_c +double a, b, c + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nld() +double ima_efncd() +extern ima_efncd + +double divzero +common /imadcomd/ divzero + +begin + # Loop through all of the image lines. + divzero = c + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do a vector + # reciprical to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nld (im, buf, v, 2) != EOF) + call arczd (a, Memd[buf[2]], Memd[buf[1]], len, + ima_efncd) + + # If imageb is constant then read imagea. If the constant + # is 1 do a vector move to imagec otherwise do a vector/scalar + # divide to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nld (im, buf, v, 2) != EOF) { + if (b == 0.0D0) + call amovkd (divzero, Memd[buf[1]], len) + else if (b == 1.0D0) + call amovd (Memd[buf[2]], Memd[buf[1]], len) + else + call adivkd (Memd[buf[2]], b, Memd[buf[1]], len) + } + + # Read imagea and imageb and do the vector divide to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nld (im, buf, v, 3) != EOF) + call advzd (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], + len, ima_efncd) + } +end + + +# IMA_EFNC -- Error function for division by zero. + +double procedure ima_efncd (a) + +double a +double divzero +common /imadcomd/ divzero + +begin + return (divzero) +end + diff --git a/pkg/images/imutil/src/generic/imamax.x b/pkg/images/imutil/src/generic/imamax.x new file mode 100644 index 00000000..36fec944 --- /dev/null +++ b/pkg/images/imutil/src/generic/imamax.x @@ -0,0 +1,212 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +# IMA_MAX -- Image arithmetic maximum value. + + +procedure ima_maxs (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +short a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nls() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do the vector/scalar + # maximum to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nls (im, buf, v, 2) != EOF) + call amaxks (Mems[buf[2]], a, Mems[buf[1]], len) + + # If imageb is constant then read imagea and do the vector/scalar + # maximum to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nls (im, buf, v, 2) != EOF) + call amaxks (Mems[buf[2]], b, Mems[buf[1]], len) + + # Read imagea and imageb and do a vector-vector maximum + # operation to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nls (im, buf, v, 3) != EOF) + call amaxs (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], len) + } +end + +procedure ima_maxi (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +int a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nli() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do the vector/scalar + # maximum to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nli (im, buf, v, 2) != EOF) + call amaxki (Memi[buf[2]], a, Memi[buf[1]], len) + + # If imageb is constant then read imagea and do the vector/scalar + # maximum to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nli (im, buf, v, 2) != EOF) + call amaxki (Memi[buf[2]], b, Memi[buf[1]], len) + + # Read imagea and imageb and do a vector-vector maximum + # operation to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nli (im, buf, v, 3) != EOF) + call amaxi (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], len) + } +end + +procedure ima_maxl (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +long a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nll() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do the vector/scalar + # maximum to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nll (im, buf, v, 2) != EOF) + call amaxkl (Meml[buf[2]], a, Meml[buf[1]], len) + + # If imageb is constant then read imagea and do the vector/scalar + # maximum to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nll (im, buf, v, 2) != EOF) + call amaxkl (Meml[buf[2]], b, Meml[buf[1]], len) + + # Read imagea and imageb and do a vector-vector maximum + # operation to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nll (im, buf, v, 3) != EOF) + call amaxl (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], len) + } +end + +procedure ima_maxr (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +real a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nlr() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do the vector/scalar + # maximum to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nlr (im, buf, v, 2) != EOF) + call amaxkr (Memr[buf[2]], a, Memr[buf[1]], len) + + # If imageb is constant then read imagea and do the vector/scalar + # maximum to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nlr (im, buf, v, 2) != EOF) + call amaxkr (Memr[buf[2]], b, Memr[buf[1]], len) + + # Read imagea and imageb and do a vector-vector maximum + # operation to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nlr (im, buf, v, 3) != EOF) + call amaxr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], len) + } +end + +procedure ima_maxd (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +double a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nld() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do the vector/scalar + # maximum to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nld (im, buf, v, 2) != EOF) + call amaxkd (Memd[buf[2]], a, Memd[buf[1]], len) + + # If imageb is constant then read imagea and do the vector/scalar + # maximum to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nld (im, buf, v, 2) != EOF) + call amaxkd (Memd[buf[2]], b, Memd[buf[1]], len) + + # Read imagea and imageb and do a vector-vector maximum + # operation to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nld (im, buf, v, 3) != EOF) + call amaxd (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], len) + } +end + diff --git a/pkg/images/imutil/src/generic/imamin.x b/pkg/images/imutil/src/generic/imamin.x new file mode 100644 index 00000000..5124db41 --- /dev/null +++ b/pkg/images/imutil/src/generic/imamin.x @@ -0,0 +1,212 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +# IMA_MIN -- Image arithmetic minimum value. + + +procedure ima_mins (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +short a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nls() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do the vector/scalar + # minimum to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nls (im, buf, v, 2) != EOF) + call aminks (Mems[buf[2]], a, Mems[buf[1]], len) + + # If imageb is constant then read imagea and do the vector/scalar + # minimum to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nls (im, buf, v, 2) != EOF) + call aminks (Mems[buf[2]], b, Mems[buf[1]], len) + + # Read imagea and imageb and do a vector-vector minimum operation + # to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nls (im, buf, v, 3) != EOF) + call amins (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], len) + } +end + +procedure ima_mini (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +int a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nli() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do the vector/scalar + # minimum to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nli (im, buf, v, 2) != EOF) + call aminki (Memi[buf[2]], a, Memi[buf[1]], len) + + # If imageb is constant then read imagea and do the vector/scalar + # minimum to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nli (im, buf, v, 2) != EOF) + call aminki (Memi[buf[2]], b, Memi[buf[1]], len) + + # Read imagea and imageb and do a vector-vector minimum operation + # to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nli (im, buf, v, 3) != EOF) + call amini (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], len) + } +end + +procedure ima_minl (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +long a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nll() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do the vector/scalar + # minimum to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nll (im, buf, v, 2) != EOF) + call aminkl (Meml[buf[2]], a, Meml[buf[1]], len) + + # If imageb is constant then read imagea and do the vector/scalar + # minimum to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nll (im, buf, v, 2) != EOF) + call aminkl (Meml[buf[2]], b, Meml[buf[1]], len) + + # Read imagea and imageb and do a vector-vector minimum operation + # to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nll (im, buf, v, 3) != EOF) + call aminl (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], len) + } +end + +procedure ima_minr (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +real a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nlr() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do the vector/scalar + # minimum to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nlr (im, buf, v, 2) != EOF) + call aminkr (Memr[buf[2]], a, Memr[buf[1]], len) + + # If imageb is constant then read imagea and do the vector/scalar + # minimum to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nlr (im, buf, v, 2) != EOF) + call aminkr (Memr[buf[2]], b, Memr[buf[1]], len) + + # Read imagea and imageb and do a vector-vector minimum operation + # to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nlr (im, buf, v, 3) != EOF) + call aminr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], len) + } +end + +procedure ima_mind (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +double a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nld() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb and do the vector/scalar + # minimum to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nld (im, buf, v, 2) != EOF) + call aminkd (Memd[buf[2]], a, Memd[buf[1]], len) + + # If imageb is constant then read imagea and do the vector/scalar + # minimum to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nld (im, buf, v, 2) != EOF) + call aminkd (Memd[buf[2]], b, Memd[buf[1]], len) + + # Read imagea and imageb and do a vector-vector minimum operation + # to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nld (im, buf, v, 3) != EOF) + call amind (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], len) + } +end + diff --git a/pkg/images/imutil/src/generic/imamul.x b/pkg/images/imutil/src/generic/imamul.x new file mode 100644 index 00000000..05fdf8a4 --- /dev/null +++ b/pkg/images/imutil/src/generic/imamul.x @@ -0,0 +1,257 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +# IMA_MUL -- Image arithmetic multiplication. + + +procedure ima_muls (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +short a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nls() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb. If the constant + # is 1 do a vector move to imagec otherwise do a vector + # multiply to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nls (im, buf, v, 2) != EOF) { + if (a == 1) + call amovs (Mems[buf[2]], Mems[buf[1]], len) + else + call amulks (Mems[buf[2]], a, Mems[buf[1]], len) + } + + # If imageb is constant then read imagea. If the constant + # is 1 do a vector move to imagec otherwise do a vector + # multiply to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nls (im, buf, v, 2) != EOF) { + if (b == 1) + call amovs (Mems[buf[2]], Mems[buf[1]], len) + else + call amulks (Mems[buf[2]], b, Mems[buf[1]], len) + } + + # Read imagea and imageb and do the vector multiply to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nls (im, buf, v, 3) != EOF) + call amuls (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], len) + } +end + +procedure ima_muli (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +int a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nli() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb. If the constant + # is 1 do a vector move to imagec otherwise do a vector + # multiply to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nli (im, buf, v, 2) != EOF) { + if (a == 1) + call amovi (Memi[buf[2]], Memi[buf[1]], len) + else + call amulki (Memi[buf[2]], a, Memi[buf[1]], len) + } + + # If imageb is constant then read imagea. If the constant + # is 1 do a vector move to imagec otherwise do a vector + # multiply to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nli (im, buf, v, 2) != EOF) { + if (b == 1) + call amovi (Memi[buf[2]], Memi[buf[1]], len) + else + call amulki (Memi[buf[2]], b, Memi[buf[1]], len) + } + + # Read imagea and imageb and do the vector multiply to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nli (im, buf, v, 3) != EOF) + call amuli (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], len) + } +end + +procedure ima_mull (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +long a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nll() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb. If the constant + # is 1 do a vector move to imagec otherwise do a vector + # multiply to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nll (im, buf, v, 2) != EOF) { + if (a == 1) + call amovl (Meml[buf[2]], Meml[buf[1]], len) + else + call amulkl (Meml[buf[2]], a, Meml[buf[1]], len) + } + + # If imageb is constant then read imagea. If the constant + # is 1 do a vector move to imagec otherwise do a vector + # multiply to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nll (im, buf, v, 2) != EOF) { + if (b == 1) + call amovl (Meml[buf[2]], Meml[buf[1]], len) + else + call amulkl (Meml[buf[2]], b, Meml[buf[1]], len) + } + + # Read imagea and imageb and do the vector multiply to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nll (im, buf, v, 3) != EOF) + call amull (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], len) + } +end + +procedure ima_mulr (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +real a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nlr() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb. If the constant + # is 1 do a vector move to imagec otherwise do a vector + # multiply to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nlr (im, buf, v, 2) != EOF) { + if (a == 1.0) + call amovr (Memr[buf[2]], Memr[buf[1]], len) + else + call amulkr (Memr[buf[2]], a, Memr[buf[1]], len) + } + + # If imageb is constant then read imagea. If the constant + # is 1 do a vector move to imagec otherwise do a vector + # multiply to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nlr (im, buf, v, 2) != EOF) { + if (b == 1.0) + call amovr (Memr[buf[2]], Memr[buf[1]], len) + else + call amulkr (Memr[buf[2]], b, Memr[buf[1]], len) + } + + # Read imagea and imageb and do the vector multiply to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nlr (im, buf, v, 3) != EOF) + call amulr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], len) + } +end + +procedure ima_muld (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +double a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nld() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb. If the constant + # is 1 do a vector move to imagec otherwise do a vector + # multiply to imagec. + if (im_a == NULL) { + im[2] = im_b + while (ima_nld (im, buf, v, 2) != EOF) { + if (a == 1.0D0) + call amovd (Memd[buf[2]], Memd[buf[1]], len) + else + call amulkd (Memd[buf[2]], a, Memd[buf[1]], len) + } + + # If imageb is constant then read imagea. If the constant + # is 1 do a vector move to imagec otherwise do a vector + # multiply to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nld (im, buf, v, 2) != EOF) { + if (b == 1.0D0) + call amovd (Memd[buf[2]], Memd[buf[1]], len) + else + call amulkd (Memd[buf[2]], b, Memd[buf[1]], len) + } + + # Read imagea and imageb and do the vector multiply to imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nld (im, buf, v, 3) != EOF) + call amuld (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], len) + } +end + diff --git a/pkg/images/imutil/src/generic/imanl.x b/pkg/images/imutil/src/generic/imanl.x new file mode 100644 index 00000000..8ec958c4 --- /dev/null +++ b/pkg/images/imutil/src/generic/imanl.x @@ -0,0 +1,159 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +# IMA_NL -- For each line in the output image lines from the input images +# are returned. The input images are repeated as necessary. EOF is returned +# when the last line of the output image has been reached. One dimensional +# images are read only once and the data pointers are assumed to be unchanged +# from previous calls. The image line vectors must be initialized externally +# and then left untouched. +# +# This procedure is typically used when operations upon lines or pixels +# make sense in mixed dimensioned images. For example to add a one dimensional +# image to all lines of a higher dimensional image or to subtract a +# two dimensional image from all bands of three dimensional image. +# The lengths of the common dimensions should generally be checked +# for equality with xt_imleneq. + + +int procedure ima_nls (im, data, v, nimages) + +pointer im[nimages] # IMIO pointers; the first one is the output +pointer data[nimages] # Returned data pointers +long v[IM_MAXDIM, nimages] # Line vectors +int nimages # Number of images + +int i + +int impnls(), imgnls() + +begin + if (impnls (im[1], data[1], v[1,1]) == EOF) + return (EOF) + + for (i=2; i <= nimages; i=i+1) { + if (imgnls (im[i], data[i], v[1,i]) == EOF) { + if (IM_NDIM(im[i]) > 1) { + call amovkl (long(1), v[1,i], IM_MAXDIM) + if (imgnls (im[i], data[i], v[1,i]) == EOF) + call error (0, "Error reading image line") + } + } + } + + return (OK) +end + +int procedure ima_nli (im, data, v, nimages) + +pointer im[nimages] # IMIO pointers; the first one is the output +pointer data[nimages] # Returned data pointers +long v[IM_MAXDIM, nimages] # Line vectors +int nimages # Number of images + +int i + +int impnli(), imgnli() + +begin + if (impnli (im[1], data[1], v[1,1]) == EOF) + return (EOF) + + for (i=2; i <= nimages; i=i+1) { + if (imgnli (im[i], data[i], v[1,i]) == EOF) { + if (IM_NDIM(im[i]) > 1) { + call amovkl (long(1), v[1,i], IM_MAXDIM) + if (imgnli (im[i], data[i], v[1,i]) == EOF) + call error (0, "Error reading image line") + } + } + } + + return (OK) +end + +int procedure ima_nll (im, data, v, nimages) + +pointer im[nimages] # IMIO pointers; the first one is the output +pointer data[nimages] # Returned data pointers +long v[IM_MAXDIM, nimages] # Line vectors +int nimages # Number of images + +int i + +int impnll(), imgnll() + +begin + if (impnll (im[1], data[1], v[1,1]) == EOF) + return (EOF) + + for (i=2; i <= nimages; i=i+1) { + if (imgnll (im[i], data[i], v[1,i]) == EOF) { + if (IM_NDIM(im[i]) > 1) { + call amovkl (long(1), v[1,i], IM_MAXDIM) + if (imgnll (im[i], data[i], v[1,i]) == EOF) + call error (0, "Error reading image line") + } + } + } + + return (OK) +end + +int procedure ima_nlr (im, data, v, nimages) + +pointer im[nimages] # IMIO pointers; the first one is the output +pointer data[nimages] # Returned data pointers +long v[IM_MAXDIM, nimages] # Line vectors +int nimages # Number of images + +int i + +int impnlr(), imgnlr() + +begin + if (impnlr (im[1], data[1], v[1,1]) == EOF) + return (EOF) + + for (i=2; i <= nimages; i=i+1) { + if (imgnlr (im[i], data[i], v[1,i]) == EOF) { + if (IM_NDIM(im[i]) > 1) { + call amovkl (long(1), v[1,i], IM_MAXDIM) + if (imgnlr (im[i], data[i], v[1,i]) == EOF) + call error (0, "Error reading image line") + } + } + } + + return (OK) +end + +int procedure ima_nld (im, data, v, nimages) + +pointer im[nimages] # IMIO pointers; the first one is the output +pointer data[nimages] # Returned data pointers +long v[IM_MAXDIM, nimages] # Line vectors +int nimages # Number of images + +int i + +int impnld(), imgnld() + +begin + if (impnld (im[1], data[1], v[1,1]) == EOF) + return (EOF) + + for (i=2; i <= nimages; i=i+1) { + if (imgnld (im[i], data[i], v[1,i]) == EOF) { + if (IM_NDIM(im[i]) > 1) { + call amovkl (long(1), v[1,i], IM_MAXDIM) + if (imgnld (im[i], data[i], v[1,i]) == EOF) + call error (0, "Error reading image line") + } + } + } + + return (OK) +end + diff --git a/pkg/images/imutil/src/generic/imasub.x b/pkg/images/imutil/src/generic/imasub.x new file mode 100644 index 00000000..1a0fcb2c --- /dev/null +++ b/pkg/images/imutil/src/generic/imasub.x @@ -0,0 +1,252 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +# IMA_SUB -- Image arithmetic subtraction. + + +procedure ima_subs (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +short a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nls() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb. Do a vector/scalar + # subtraction and then negate the result. + if (im_a == NULL) { + im[2] = im_b + while (ima_nls (im, buf, v, 2) != EOF) { + if (a != 0) { + call asubks (Mems[buf[2]], a, Mems[buf[1]], len) + call anegs (Mems[buf[1]], Mems[buf[1]], len) + } else + call anegs (Mems[buf[2]], Mems[buf[1]], len) + } + + # If imageb is constant then read imagea and do a vector/scalar + # subtraction to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nls (im, buf, v, 2) != EOF) { + if (b == 0) + call amovs (Mems[buf[2]], Mems[buf[1]], len) + else + call asubks (Mems[buf[2]], b, Mems[buf[1]], len) + } + + # Read imagea and imageb and do a vector subtraction into imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nls (im, buf, v, 3) != EOF) + call asubs (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], len) + } +end + +procedure ima_subi (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +int a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nli() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb. Do a vector/scalar + # subtraction and then negate the result. + if (im_a == NULL) { + im[2] = im_b + while (ima_nli (im, buf, v, 2) != EOF) { + if (a != 0) { + call asubki (Memi[buf[2]], a, Memi[buf[1]], len) + call anegi (Memi[buf[1]], Memi[buf[1]], len) + } else + call anegi (Memi[buf[2]], Memi[buf[1]], len) + } + + # If imageb is constant then read imagea and do a vector/scalar + # subtraction to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nli (im, buf, v, 2) != EOF) { + if (b == 0) + call amovi (Memi[buf[2]], Memi[buf[1]], len) + else + call asubki (Memi[buf[2]], b, Memi[buf[1]], len) + } + + # Read imagea and imageb and do a vector subtraction into imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nli (im, buf, v, 3) != EOF) + call asubi (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], len) + } +end + +procedure ima_subl (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +long a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nll() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb. Do a vector/scalar + # subtraction and then negate the result. + if (im_a == NULL) { + im[2] = im_b + while (ima_nll (im, buf, v, 2) != EOF) { + if (a != 0) { + call asubkl (Meml[buf[2]], a, Meml[buf[1]], len) + call anegl (Meml[buf[1]], Meml[buf[1]], len) + } else + call anegl (Meml[buf[2]], Meml[buf[1]], len) + } + + # If imageb is constant then read imagea and do a vector/scalar + # subtraction to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nll (im, buf, v, 2) != EOF) { + if (b == 0) + call amovl (Meml[buf[2]], Meml[buf[1]], len) + else + call asubkl (Meml[buf[2]], b, Meml[buf[1]], len) + } + + # Read imagea and imageb and do a vector subtraction into imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nll (im, buf, v, 3) != EOF) + call asubl (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], len) + } +end + +procedure ima_subr (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +real a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nlr() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb. Do a vector/scalar + # subtraction and then negate the result. + if (im_a == NULL) { + im[2] = im_b + while (ima_nlr (im, buf, v, 2) != EOF) { + if (a != 0.0) { + call asubkr (Memr[buf[2]], a, Memr[buf[1]], len) + call anegr (Memr[buf[1]], Memr[buf[1]], len) + } else + call anegr (Memr[buf[2]], Memr[buf[1]], len) + } + + # If imageb is constant then read imagea and do a vector/scalar + # subtraction to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nlr (im, buf, v, 2) != EOF) { + if (b == 0.0) + call amovr (Memr[buf[2]], Memr[buf[1]], len) + else + call asubkr (Memr[buf[2]], b, Memr[buf[1]], len) + } + + # Read imagea and imageb and do a vector subtraction into imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nlr (im, buf, v, 3) != EOF) + call asubr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], len) + } +end + +procedure ima_subd (im_a, im_b, im_c, a, b) + +pointer im_a, im_b, im_c +double a, b + +int len +pointer im[3], buf[3] +long v[IM_MAXDIM, 3] + +int ima_nld() + +begin + # Loop through all of the image lines. + im[1] = im_c + len = IM_LEN (im[1], 1) + call amovkl (long(1), v, 3 * IM_MAXDIM) + + # If imagea is constant then read imageb. Do a vector/scalar + # subtraction and then negate the result. + if (im_a == NULL) { + im[2] = im_b + while (ima_nld (im, buf, v, 2) != EOF) { + if (a != 0.0D0) { + call asubkd (Memd[buf[2]], a, Memd[buf[1]], len) + call anegd (Memd[buf[1]], Memd[buf[1]], len) + } else + call anegd (Memd[buf[2]], Memd[buf[1]], len) + } + + # If imageb is constant then read imagea and do a vector/scalar + # subtraction to imagec. + } else if (im_b == NULL) { + im[2] = im_a + while (ima_nld (im, buf, v, 2) != EOF) { + if (b == 0.0D0) + call amovd (Memd[buf[2]], Memd[buf[1]], len) + else + call asubkd (Memd[buf[2]], b, Memd[buf[1]], len) + } + + # Read imagea and imageb and do a vector subtraction into imagec. + } else { + im[2] = im_a + im[3] = im_b + while (ima_nld (im, buf, v, 3) != EOF) + call asubd (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], len) + } +end + diff --git a/pkg/images/imutil/src/generic/imfuncs.x b/pkg/images/imutil/src/generic/imfuncs.x new file mode 100644 index 00000000..67bc4ed5 --- /dev/null +++ b/pkg/images/imutil/src/generic/imfuncs.x @@ -0,0 +1,1613 @@ +include <imhdr.h> +include <mach.h> +include <math.h> + + + +# IF_LOG10 -- Compute the base 10 logarithm of image1 and write the results to +# image2. + +procedure if_log10r (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +real if_elogr() +extern if_elogr() +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call alogr (Memr[buf1], Memr[buf2], npix, if_elogr) +end + + +# IF_ELOG -- The error function for log10. Note that MAX_EXPONENT is +# currently an integer so it is converted to the appropriate data type +# before being returned. + +real procedure if_elogr (x) + +real x # the input pixel value + +begin + return (real(-MAX_EXPONENT)) +end + + +# IF_ALOG10 -- Take the power of 10 of image1 and write the results to image2. + +procedure if_alog10r (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_va10r (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VA10 -- Take the antilog (base 10) of a vector. + +procedure if_va10r (a, b, n) + +real a[n] # the input vector +real b[n] # the output vector +int n # the number of points + +int i +real maxexp, maxval + +begin + maxexp = MAX_EXPONENT + maxval = MAX_REAL + + do i = 1, n { + if (a[i] >= maxexp) + b[i] = maxval + else if (a[i] <= (-maxexp)) + b[i] = 0.0 + else + b[i] = 10.0 ** a[i] + } +end + + +# IF_LN -- Take the natural log of the pixels in image1 and write the results +# to image2. + +procedure if_lnr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 + +real if_elnr() +extern if_elnr() +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call allnr (Memr[buf1], Memr[buf2], npix, if_elnr) +end + + +# IF_ELN -- The error function for the natural logarithm. + +real procedure if_elnr (x) + +real x # input value + +begin + return (real (LN_10) * real(-MAX_EXPONENT)) +end + + +# IF_ALN -- Take the natural antilog of the pixels in image1 and write the +# results to image2. + +procedure if_alnr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_valnr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VALN -- Take the natural antilog of a vector. + +procedure if_valnr (a, b, n) + +real a[n] # the input vector +real b[n] # the output vector +int n # the number of pixels + +int i +real maxexp, maxval, eval + +begin + maxexp = log (10.0 ** real (MAX_EXPONENT)) + maxval = MAX_REAL + eval = real (BASE_E) + + do i = 1, n { + if (a[i] >= maxexp) + b[i] = maxval + else if (a[i] <= -maxexp) + b[i] = 0.0 + else + b[i] = eval ** a[i] + } +end + + +# IF_SQR -- Take the square root of pixels in image1 and write the results +# to image2. + +procedure if_sqrr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +real if_esqrr() +extern if_esqrr() +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call asqrr (Memr[buf1], Memr[buf2], npix, if_esqrr) +end + + +# IF_ESQR -- Error function for the square root. + +real procedure if_esqrr (x) + +real x # input value + +begin + return (0.0) +end + + +# IF_SQUARE -- Take the square of the pixels in image1 and write to image2. +procedure if_squarer (im1, im2) + +pointer im1 # the input image pointer +pointer im2 # the output image pointer + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call apowkr (Memr[buf1], 2, Memr[buf2], npix) +end + + +# IF_CBRT -- Take the cube root of the pixels in image1 and write the results +# to image2. + +procedure if_cbrtr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_vcbrtr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VCBRT -- Compute the cube root of a vector. + +procedure if_vcbrtr (a, b, n) + +real a[n] # the input vector +real b[n] # the output vector +int n # the number of pixels + +int i +real onethird + +begin + onethird = 1.0 / 3.0 + do i = 1, n { + if (a[i] >= 0.0) { + b[i] = a[i] ** onethird + } else { + b[i] = -a[i] + b[i] = - (b[i] ** onethird) + } + } +end + + +# IF_CUBE -- Take the cube of the pixels in image1 and write the results to +# image2. + +procedure if_cuber (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call apowkr (Memr[buf1], 3, Memr[buf2], npix) +end + + +# IF_COS -- Take cosine of pixels in image1 and write the results to image2. + +procedure if_cosr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_vcosr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VCOS - Compute the cosine of a vector. + +procedure if_vcosr (a, b, n) + +real a[n] # the input vector +real b[n] # the output vector +int n # the number of pixels + +int i + +begin + do i = 1, n + b[i] = cos(a[i]) +end + + +# IF_SIN -- Take sine of the pixels in image1 and write the results to image2. + +procedure if_sinr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +pointer buf1, buf2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_vsinr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VSIN - Take the sine of a vector. + +procedure if_vsinr (a, b, n) + +real a[n] # the input vector +real b[n] # the output vector +int n # the number of pixels + +int i + +begin + do i = 1, n + b[i] = sin(a[i]) +end + + +# IF_TAN -- Take tangent of pixels in image1 and write the results to image2. + +procedure if_tanr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_vtanr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VTAN - Take the tangent of a vector. + +procedure if_vtanr (a, b, n) + +real a[n] # the input vector +real b[n] # the output vector +int n # the number of pixels + +int i + +begin + do i = 1, n + b[i] = tan(a[i]) +end + + +# IF_ACOS -- Take arccosine of pixels in image1 and write the results to image2. + +procedure if_acosr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +pointer buf1, buf2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_vacosr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VACOS - Take the arccosine of a vector. + +procedure if_vacosr (a, b, n) + +real a[n] # the input vector +real b[n] # the output vector +int n # the number of pixels + +int i + +begin + do i = 1, n { + if (a[i] > 1.0) + b[i] = acos (1.0) + else if (a[i] < -1.0) + b[i] = acos (-1.0) + else + b[i] = acos(a[i]) + } +end + + +# IF_ASIN -- Take arcsine of pixels in image1 and write the results to image2. + +procedure if_asinr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +pointer buf1, buf2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_vasinr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VASIN - Take arcsine of vector + +procedure if_vasinr (a, b, n) + +real a[n] +real b[n] +int n + +int i + +begin + do i = 1, n { + if (a[i] > 1.0) + b[i] = asin (1.0) + else if (a[i] < -1.0) + b[i] = asin (-1.0) + else + b[i] = asin(a[i]) + } +end + + +# IF_ATAN -- Take arctangent of pixels in image1 and write the results to +# image2. + +procedure if_atanr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_vatanr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VATAN - Take the arctangent of a vector. + +procedure if_vatanr (a, b, n) + +real a[n] +real b[n] +int n + +int i + +begin + do i = 1, n + b[i] = atan(a[i]) +end + + +# IF_HCOS -- Take the hyperbolic cosine of pixels in image1 and write the +# results to image2. + +procedure if_hcosr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_vhcosr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VHCOS - Take the hyperbolic cosine of a vector. + +procedure if_vhcosr (a, b, n) + +real a[n] # the input vector +real b[n] # the output vector +int n # the number of pixels + +int i +real maxexp, maxval + +begin + maxexp = log (10.0 ** real(MAX_EXPONENT)) + maxval = MAX_REAL + + do i = 1, n { + if (abs (a[i]) >= maxexp) + b[i] = maxval + else + b[i] = cosh (a[i]) + } +end + + +# IF_HSIN -- Take the hyperbolic sine of pixels in image1 and write the +# results to image2. + +procedure if_hsinr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +pointer buf1, buf2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_vhsinr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VHSIN - Take the hyperbolic sine of a vector. + +procedure if_vhsinr (a, b, n) + +real a[n] # the input vector +real b[n] # the output vector +int n # the number of pixels + +int i +real maxexp, maxval + +begin + maxexp = log (10.0 ** real(MAX_EXPONENT)) + maxval = MAX_REAL + + do i = 1, n { + if (a[i] >= maxexp) + b[i] = maxval + else if (a[i] <= -maxexp) + b[i] = -maxval + else + b[i] = sinh(a[i]) + } +end + + +# IF_HTAN -- Take the hyperbolic tangent of pixels in image1 and write the +# results to image2. + +procedure if_htanr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +pointer buf1, buf2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call if_vhtanr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_VHTAN - Take the hyperbolic tangent of a vector. + +procedure if_vhtanr (a, b, n) + +real a[n] # the input vector +real b[n] # the output vector +int n # the number of pixels + +int i + +begin + do i = 1, n + b[i] = tanh(a[i]) +end + + +# IF_RECIP -- Take the reciprocal of the pixels in image1 and write the +# results to image2. + +procedure if_recipr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +real if_erecipr() +extern if_erecipr() +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call arczr (1.0, Memr[buf1], Memr[buf2], npix, if_erecipr) +end + + +# IF_ERECIP -- Error function for the reciprocal computation. + +real procedure if_erecipr (x) + +real x + +begin + return (0.0) +end + + + +# IF_LOG10 -- Compute the base 10 logarithm of image1 and write the results to +# image2. + +procedure if_log10d (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +double if_elogd() +extern if_elogd() +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call alogd (Memd[buf1], Memd[buf2], npix, if_elogd) +end + + +# IF_ELOG -- The error function for log10. Note that MAX_EXPONENT is +# currently an integer so it is converted to the appropriate data type +# before being returned. + +double procedure if_elogd (x) + +double x # the input pixel value + +begin + return (double(-MAX_EXPONENT)) +end + + +# IF_ALOG10 -- Take the power of 10 of image1 and write the results to image2. + +procedure if_alog10d (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_va10d (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VA10 -- Take the antilog (base 10) of a vector. + +procedure if_va10d (a, b, n) + +double a[n] # the input vector +double b[n] # the output vector +int n # the number of points + +int i +double maxexp, maxval + +begin + maxexp = MAX_EXPONENT + maxval = MAX_REAL + + do i = 1, n { + if (a[i] >= maxexp) + b[i] = maxval + else if (a[i] <= (-maxexp)) + b[i] = 0.0D0 + else + b[i] = 10.0D0 ** a[i] + } +end + + +# IF_LN -- Take the natural log of the pixels in image1 and write the results +# to image2. + +procedure if_lnd (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 + +double if_elnd() +extern if_elnd() +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call allnd (Memd[buf1], Memd[buf2], npix, if_elnd) +end + + +# IF_ELN -- The error function for the natural logarithm. + +double procedure if_elnd (x) + +double x # input value + +begin + return (double (LN_10) * double(-MAX_EXPONENT)) +end + + +# IF_ALN -- Take the natural antilog of the pixels in image1 and write the +# results to image2. + +procedure if_alnd (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_valnd (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VALN -- Take the natural antilog of a vector. + +procedure if_valnd (a, b, n) + +double a[n] # the input vector +double b[n] # the output vector +int n # the number of pixels + +int i +double maxexp, maxval, eval + +begin + maxexp = log (10.0D0 ** double (MAX_EXPONENT)) + maxval = MAX_REAL + eval = double (BASE_E) + + do i = 1, n { + if (a[i] >= maxexp) + b[i] = maxval + else if (a[i] <= -maxexp) + b[i] = 0.0D0 + else + b[i] = eval ** a[i] + } +end + + +# IF_SQR -- Take the square root of pixels in image1 and write the results +# to image2. + +procedure if_sqrd (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +double if_esqrd() +extern if_esqrd() +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call asqrd (Memd[buf1], Memd[buf2], npix, if_esqrd) +end + + +# IF_ESQR -- Error function for the square root. + +double procedure if_esqrd (x) + +double x # input value + +begin + return (0.0D0) +end + + +# IF_SQUARE -- Take the square of the pixels in image1 and write to image2. +procedure if_squared (im1, im2) + +pointer im1 # the input image pointer +pointer im2 # the output image pointer + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call apowkd (Memd[buf1], 2, Memd[buf2], npix) +end + + +# IF_CBRT -- Take the cube root of the pixels in image1 and write the results +# to image2. + +procedure if_cbrtd (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_vcbrtd (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VCBRT -- Compute the cube root of a vector. + +procedure if_vcbrtd (a, b, n) + +double a[n] # the input vector +double b[n] # the output vector +int n # the number of pixels + +int i +double onethird + +begin + onethird = 1.0D0 / 3.0D0 + do i = 1, n { + if (a[i] >= 0.0D0) { + b[i] = a[i] ** onethird + } else { + b[i] = -a[i] + b[i] = - (b[i] ** onethird) + } + } +end + + +# IF_CUBE -- Take the cube of the pixels in image1 and write the results to +# image2. + +procedure if_cubed (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call apowkd (Memd[buf1], 3, Memd[buf2], npix) +end + + +# IF_COS -- Take cosine of pixels in image1 and write the results to image2. + +procedure if_cosd (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_vcosd (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VCOS - Compute the cosine of a vector. + +procedure if_vcosd (a, b, n) + +double a[n] # the input vector +double b[n] # the output vector +int n # the number of pixels + +int i + +begin + do i = 1, n + b[i] = cos(a[i]) +end + + +# IF_SIN -- Take sine of the pixels in image1 and write the results to image2. + +procedure if_sind (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +pointer buf1, buf2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_vsind (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VSIN - Take the sine of a vector. + +procedure if_vsind (a, b, n) + +double a[n] # the input vector +double b[n] # the output vector +int n # the number of pixels + +int i + +begin + do i = 1, n + b[i] = sin(a[i]) +end + + +# IF_TAN -- Take tangent of pixels in image1 and write the results to image2. + +procedure if_tand (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_vtand (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VTAN - Take the tangent of a vector. + +procedure if_vtand (a, b, n) + +double a[n] # the input vector +double b[n] # the output vector +int n # the number of pixels + +int i + +begin + do i = 1, n + b[i] = tan(a[i]) +end + + +# IF_ACOS -- Take arccosine of pixels in image1 and write the results to image2. + +procedure if_acosd (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +pointer buf1, buf2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_vacosd (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VACOS - Take the arccosine of a vector. + +procedure if_vacosd (a, b, n) + +double a[n] # the input vector +double b[n] # the output vector +int n # the number of pixels + +int i + +begin + do i = 1, n { + if (a[i] > 1.0D0) + b[i] = acos (1.0D0) + else if (a[i] < -1.0D0) + b[i] = acos (-1.0D0) + else + b[i] = acos(a[i]) + } +end + + +# IF_ASIN -- Take arcsine of pixels in image1 and write the results to image2. + +procedure if_asind (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +pointer buf1, buf2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_vasind (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VASIN - Take arcsine of vector + +procedure if_vasind (a, b, n) + +double a[n] +double b[n] +int n + +int i + +begin + do i = 1, n { + if (a[i] > 1.0D0) + b[i] = asin (1.0D0) + else if (a[i] < -1.0D0) + b[i] = asin (-1.0D0) + else + b[i] = asin(a[i]) + } +end + + +# IF_ATAN -- Take arctangent of pixels in image1 and write the results to +# image2. + +procedure if_atand (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_vatand (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VATAN - Take the arctangent of a vector. + +procedure if_vatand (a, b, n) + +double a[n] +double b[n] +int n + +int i + +begin + do i = 1, n + b[i] = atan(a[i]) +end + + +# IF_HCOS -- Take the hyperbolic cosine of pixels in image1 and write the +# results to image2. + +procedure if_hcosd (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_vhcosd (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VHCOS - Take the hyperbolic cosine of a vector. + +procedure if_vhcosd (a, b, n) + +double a[n] # the input vector +double b[n] # the output vector +int n # the number of pixels + +int i +double maxexp, maxval + +begin + maxexp = log (10.0D0 ** double(MAX_EXPONENT)) + maxval = MAX_REAL + + do i = 1, n { + if (abs (a[i]) >= maxexp) + b[i] = maxval + else + b[i] = cosh (a[i]) + } +end + + +# IF_HSIN -- Take the hyperbolic sine of pixels in image1 and write the +# results to image2. + +procedure if_hsind (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +pointer buf1, buf2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_vhsind (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VHSIN - Take the hyperbolic sine of a vector. + +procedure if_vhsind (a, b, n) + +double a[n] # the input vector +double b[n] # the output vector +int n # the number of pixels + +int i +double maxexp, maxval + +begin + maxexp = log (10.0D0 ** double(MAX_EXPONENT)) + maxval = MAX_REAL + + do i = 1, n { + if (a[i] >= maxexp) + b[i] = maxval + else if (a[i] <= -maxexp) + b[i] = -maxval + else + b[i] = sinh(a[i]) + } +end + + +# IF_HTAN -- Take the hyperbolic tangent of pixels in image1 and write the +# results to image2. + +procedure if_htand (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +pointer buf1, buf2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call if_vhtand (Memd[buf1], Memd[buf2], npix) +end + + +# IF_VHTAN - Take the hyperbolic tangent of a vector. + +procedure if_vhtand (a, b, n) + +double a[n] # the input vector +double b[n] # the output vector +int n # the number of pixels + +int i + +begin + do i = 1, n + b[i] = tanh(a[i]) +end + + +# IF_RECIP -- Take the reciprocal of the pixels in image1 and write the +# results to image2. + +procedure if_recipd (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +double if_erecipd() +extern if_erecipd() +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call arczd (1.0, Memd[buf1], Memd[buf2], npix, if_erecipd) +end + + +# IF_ERECIP -- Error function for the reciprocal computation. + +double procedure if_erecipd (x) + +double x + +begin + return (0.0D0) +end + + + + + +# IF_ABS -- Take the absolute value of pixels in image1 and write the results +# to image2. + +procedure if_absl (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnll(), impnll() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnll (im1, buf1, v1) != EOF) && + (impnll (im2, buf2, v2) != EOF)) + call aabsl (Meml[buf1], Meml[buf2], npix) +end + + +# IF_NEG -- Take negative of pixels in image1 and write the results to image2. + +procedure if_negl (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnll(), impnll() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnll (im1, buf1, v1) != EOF) && + (impnll (im2, buf2, v2) != EOF)) + call anegl (Meml[buf1], Meml[buf2], npix) +end + + + +# IF_ABS -- Take the absolute value of pixels in image1 and write the results +# to image2. + +procedure if_absr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call aabsr (Memr[buf1], Memr[buf2], npix) +end + + +# IF_NEG -- Take negative of pixels in image1 and write the results to image2. + +procedure if_negr (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnlr(), impnlr() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnlr (im1, buf1, v1) != EOF) && + (impnlr (im2, buf2, v2) != EOF)) + call anegr (Memr[buf1], Memr[buf2], npix) +end + + + +# IF_ABS -- Take the absolute value of pixels in image1 and write the results +# to image2. + +procedure if_absd (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call aabsd (Memd[buf1], Memd[buf2], npix) +end + + +# IF_NEG -- Take negative of pixels in image1 and write the results to image2. + +procedure if_negd (im1, im2) + +pointer im1 # pointer to the input image +pointer im2 # pointer to the output image + +int npix +long v1[IM_MAXDIM], v2[IM_MAXDIM] +pointer buf1, buf2 +int imgnld(), impnld() + +begin + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im1, 1) + while ((imgnld (im1, buf1, v1) != EOF) && + (impnld (im2, buf2, v2) != EOF)) + call anegd (Memd[buf1], Memd[buf2], npix) +end + + diff --git a/pkg/images/imutil/src/generic/imjoin.x b/pkg/images/imutil/src/generic/imjoin.x new file mode 100644 index 00000000..83b02541 --- /dev/null +++ b/pkg/images/imutil/src/generic/imjoin.x @@ -0,0 +1,527 @@ +include <imhdr.h> + +define VPTR Memi[$1+$2-1] # Array of axis vector pointers + + + +# IMJOIN -- Join the set of input images into an output image along the +# specified axis, any dimension. + +procedure imjoins (inptr, nimages, out, joindim, outtype) + +pointer inptr[nimages] #I Input IMIO pointers +int nimages #I Number of input images +pointer out #I Output IMIO pointer +int joindim #I Dimension along which to join images +int outtype #I Output datatype + +int i, image, line, nlines, nbands, stat, cum_len +pointer sp, vin, vout, in, inbuf, outbuf + +pointer imgnls() +pointer impnls() + +begin + # Allocate working space. + call smark (sp) + call salloc (vin, nimages, TY_INT) + call salloc (vout, IM_MAXDIM, TY_LONG) + + # Initialize the v vectors. + call amovkl (long(1), Meml[vout], IM_MAXDIM) + do image = 1, nimages { + call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG) + call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM) + } + + # Join input images along the specified dimension. Joins along + # columns and lines require processing in special order, all others + # in the same order. In the first two cases we process all input + # images in inner loops, so we have to keep all those image + # descriptors open. + + switch (joindim) { + case 1: # join columns + nlines = 1 + do i = 2, IM_NDIM(out) + nlines = nlines * IM_LEN(out,i) + do i = 1, nlines { + stat = impnls (out, outbuf, Meml[vout]) + cum_len = 0 + do image = 1, nimages { + in = inptr[image] + stat = imgnls (in, inbuf, Meml[VPTR(vin,image)]) + call amovs (Mems[inbuf], Mems[outbuf+cum_len], + IM_LEN(in,1)) + cum_len = cum_len + IM_LEN(in,1) + } + } + + case 2: # join lines + nbands = 1 + do i = 3, IM_NDIM(out) + nbands = nbands * IM_LEN(out,i) + do i = 1, nbands { + do image = 1, nimages { + in = inptr[image] + do line = 1, IM_LEN(in,2) { + stat = impnls (out, outbuf, Meml[vout]) + stat = imgnls (in, inbuf, Meml[VPTR(vin,image)]) + call amovs (Mems[inbuf], Mems[outbuf], IM_LEN(in,1)) + } + } + } + + default: # join bands or higher + do image = 1, nimages { + in = inptr[image] + nlines = 1 + do i = 2, IM_NDIM(in) + nlines = nlines * IM_LEN(in,i) + do i = 1, nlines { + stat = impnls (out, outbuf, Meml[vout]) + stat = imgnls (in, inbuf, Meml[VPTR(vin,image)]) + call amovs (Mems[inbuf], Mems[outbuf], IM_LEN(in,1)) + } + } + } + + call sfree (sp) +end + + + +# IMJOIN -- Join the set of input images into an output image along the +# specified axis, any dimension. + +procedure imjoini (inptr, nimages, out, joindim, outtype) + +pointer inptr[nimages] #I Input IMIO pointers +int nimages #I Number of input images +pointer out #I Output IMIO pointer +int joindim #I Dimension along which to join images +int outtype #I Output datatype + +int i, image, line, nlines, nbands, stat, cum_len +pointer sp, vin, vout, in, inbuf, outbuf + +pointer imgnli() +pointer impnli() + +begin + # Allocate working space. + call smark (sp) + call salloc (vin, nimages, TY_INT) + call salloc (vout, IM_MAXDIM, TY_LONG) + + # Initialize the v vectors. + call amovkl (long(1), Meml[vout], IM_MAXDIM) + do image = 1, nimages { + call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG) + call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM) + } + + # Join input images along the specified dimension. Joins along + # columns and lines require processing in special order, all others + # in the same order. In the first two cases we process all input + # images in inner loops, so we have to keep all those image + # descriptors open. + + switch (joindim) { + case 1: # join columns + nlines = 1 + do i = 2, IM_NDIM(out) + nlines = nlines * IM_LEN(out,i) + do i = 1, nlines { + stat = impnli (out, outbuf, Meml[vout]) + cum_len = 0 + do image = 1, nimages { + in = inptr[image] + stat = imgnli (in, inbuf, Meml[VPTR(vin,image)]) + call amovi (Memi[inbuf], Memi[outbuf+cum_len], + IM_LEN(in,1)) + cum_len = cum_len + IM_LEN(in,1) + } + } + + case 2: # join lines + nbands = 1 + do i = 3, IM_NDIM(out) + nbands = nbands * IM_LEN(out,i) + do i = 1, nbands { + do image = 1, nimages { + in = inptr[image] + do line = 1, IM_LEN(in,2) { + stat = impnli (out, outbuf, Meml[vout]) + stat = imgnli (in, inbuf, Meml[VPTR(vin,image)]) + call amovi (Memi[inbuf], Memi[outbuf], IM_LEN(in,1)) + } + } + } + + default: # join bands or higher + do image = 1, nimages { + in = inptr[image] + nlines = 1 + do i = 2, IM_NDIM(in) + nlines = nlines * IM_LEN(in,i) + do i = 1, nlines { + stat = impnli (out, outbuf, Meml[vout]) + stat = imgnli (in, inbuf, Meml[VPTR(vin,image)]) + call amovi (Memi[inbuf], Memi[outbuf], IM_LEN(in,1)) + } + } + } + + call sfree (sp) +end + + + +# IMJOIN -- Join the set of input images into an output image along the +# specified axis, any dimension. + +procedure imjoinl (inptr, nimages, out, joindim, outtype) + +pointer inptr[nimages] #I Input IMIO pointers +int nimages #I Number of input images +pointer out #I Output IMIO pointer +int joindim #I Dimension along which to join images +int outtype #I Output datatype + +int i, image, line, nlines, nbands, stat, cum_len +pointer sp, vin, vout, in, inbuf, outbuf + +pointer imgnll() +pointer impnll() + +begin + # Allocate working space. + call smark (sp) + call salloc (vin, nimages, TY_INT) + call salloc (vout, IM_MAXDIM, TY_LONG) + + # Initialize the v vectors. + call amovkl (long(1), Meml[vout], IM_MAXDIM) + do image = 1, nimages { + call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG) + call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM) + } + + # Join input images along the specified dimension. Joins along + # columns and lines require processing in special order, all others + # in the same order. In the first two cases we process all input + # images in inner loops, so we have to keep all those image + # descriptors open. + + switch (joindim) { + case 1: # join columns + nlines = 1 + do i = 2, IM_NDIM(out) + nlines = nlines * IM_LEN(out,i) + do i = 1, nlines { + stat = impnll (out, outbuf, Meml[vout]) + cum_len = 0 + do image = 1, nimages { + in = inptr[image] + stat = imgnll (in, inbuf, Meml[VPTR(vin,image)]) + call amovl (Meml[inbuf], Meml[outbuf+cum_len], + IM_LEN(in,1)) + cum_len = cum_len + IM_LEN(in,1) + } + } + + case 2: # join lines + nbands = 1 + do i = 3, IM_NDIM(out) + nbands = nbands * IM_LEN(out,i) + do i = 1, nbands { + do image = 1, nimages { + in = inptr[image] + do line = 1, IM_LEN(in,2) { + stat = impnll (out, outbuf, Meml[vout]) + stat = imgnll (in, inbuf, Meml[VPTR(vin,image)]) + call amovl (Meml[inbuf], Meml[outbuf], IM_LEN(in,1)) + } + } + } + + default: # join bands or higher + do image = 1, nimages { + in = inptr[image] + nlines = 1 + do i = 2, IM_NDIM(in) + nlines = nlines * IM_LEN(in,i) + do i = 1, nlines { + stat = impnll (out, outbuf, Meml[vout]) + stat = imgnll (in, inbuf, Meml[VPTR(vin,image)]) + call amovl (Meml[inbuf], Meml[outbuf], IM_LEN(in,1)) + } + } + } + + call sfree (sp) +end + + + +# IMJOIN -- Join the set of input images into an output image along the +# specified axis, any dimension. + +procedure imjoinr (inptr, nimages, out, joindim, outtype) + +pointer inptr[nimages] #I Input IMIO pointers +int nimages #I Number of input images +pointer out #I Output IMIO pointer +int joindim #I Dimension along which to join images +int outtype #I Output datatype + +int i, image, line, nlines, nbands, stat, cum_len +pointer sp, vin, vout, in, inbuf, outbuf + +pointer imgnlr() +pointer impnlr() + +begin + # Allocate working space. + call smark (sp) + call salloc (vin, nimages, TY_INT) + call salloc (vout, IM_MAXDIM, TY_LONG) + + # Initialize the v vectors. + call amovkl (long(1), Meml[vout], IM_MAXDIM) + do image = 1, nimages { + call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG) + call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM) + } + + # Join input images along the specified dimension. Joins along + # columns and lines require processing in special order, all others + # in the same order. In the first two cases we process all input + # images in inner loops, so we have to keep all those image + # descriptors open. + + switch (joindim) { + case 1: # join columns + nlines = 1 + do i = 2, IM_NDIM(out) + nlines = nlines * IM_LEN(out,i) + do i = 1, nlines { + stat = impnlr (out, outbuf, Meml[vout]) + cum_len = 0 + do image = 1, nimages { + in = inptr[image] + stat = imgnlr (in, inbuf, Meml[VPTR(vin,image)]) + call amovr (Memr[inbuf], Memr[outbuf+cum_len], + IM_LEN(in,1)) + cum_len = cum_len + IM_LEN(in,1) + } + } + + case 2: # join lines + nbands = 1 + do i = 3, IM_NDIM(out) + nbands = nbands * IM_LEN(out,i) + do i = 1, nbands { + do image = 1, nimages { + in = inptr[image] + do line = 1, IM_LEN(in,2) { + stat = impnlr (out, outbuf, Meml[vout]) + stat = imgnlr (in, inbuf, Meml[VPTR(vin,image)]) + call amovr (Memr[inbuf], Memr[outbuf], IM_LEN(in,1)) + } + } + } + + default: # join bands or higher + do image = 1, nimages { + in = inptr[image] + nlines = 1 + do i = 2, IM_NDIM(in) + nlines = nlines * IM_LEN(in,i) + do i = 1, nlines { + stat = impnlr (out, outbuf, Meml[vout]) + stat = imgnlr (in, inbuf, Meml[VPTR(vin,image)]) + call amovr (Memr[inbuf], Memr[outbuf], IM_LEN(in,1)) + } + } + } + + call sfree (sp) +end + + + +# IMJOIN -- Join the set of input images into an output image along the +# specified axis, any dimension. + +procedure imjoind (inptr, nimages, out, joindim, outtype) + +pointer inptr[nimages] #I Input IMIO pointers +int nimages #I Number of input images +pointer out #I Output IMIO pointer +int joindim #I Dimension along which to join images +int outtype #I Output datatype + +int i, image, line, nlines, nbands, stat, cum_len +pointer sp, vin, vout, in, inbuf, outbuf + +pointer imgnld() +pointer impnld() + +begin + # Allocate working space. + call smark (sp) + call salloc (vin, nimages, TY_INT) + call salloc (vout, IM_MAXDIM, TY_LONG) + + # Initialize the v vectors. + call amovkl (long(1), Meml[vout], IM_MAXDIM) + do image = 1, nimages { + call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG) + call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM) + } + + # Join input images along the specified dimension. Joins along + # columns and lines require processing in special order, all others + # in the same order. In the first two cases we process all input + # images in inner loops, so we have to keep all those image + # descriptors open. + + switch (joindim) { + case 1: # join columns + nlines = 1 + do i = 2, IM_NDIM(out) + nlines = nlines * IM_LEN(out,i) + do i = 1, nlines { + stat = impnld (out, outbuf, Meml[vout]) + cum_len = 0 + do image = 1, nimages { + in = inptr[image] + stat = imgnld (in, inbuf, Meml[VPTR(vin,image)]) + call amovd (Memd[inbuf], Memd[outbuf+cum_len], + IM_LEN(in,1)) + cum_len = cum_len + IM_LEN(in,1) + } + } + + case 2: # join lines + nbands = 1 + do i = 3, IM_NDIM(out) + nbands = nbands * IM_LEN(out,i) + do i = 1, nbands { + do image = 1, nimages { + in = inptr[image] + do line = 1, IM_LEN(in,2) { + stat = impnld (out, outbuf, Meml[vout]) + stat = imgnld (in, inbuf, Meml[VPTR(vin,image)]) + call amovd (Memd[inbuf], Memd[outbuf], IM_LEN(in,1)) + } + } + } + + default: # join bands or higher + do image = 1, nimages { + in = inptr[image] + nlines = 1 + do i = 2, IM_NDIM(in) + nlines = nlines * IM_LEN(in,i) + do i = 1, nlines { + stat = impnld (out, outbuf, Meml[vout]) + stat = imgnld (in, inbuf, Meml[VPTR(vin,image)]) + call amovd (Memd[inbuf], Memd[outbuf], IM_LEN(in,1)) + } + } + } + + call sfree (sp) +end + + + +# IMJOIN -- Join the set of input images into an output image along the +# specified axis, any dimension. + +procedure imjoinx (inptr, nimages, out, joindim, outtype) + +pointer inptr[nimages] #I Input IMIO pointers +int nimages #I Number of input images +pointer out #I Output IMIO pointer +int joindim #I Dimension along which to join images +int outtype #I Output datatype + +int i, image, line, nlines, nbands, stat, cum_len +pointer sp, vin, vout, in, inbuf, outbuf + +pointer imgnlx() +pointer impnlx() + +begin + # Allocate working space. + call smark (sp) + call salloc (vin, nimages, TY_INT) + call salloc (vout, IM_MAXDIM, TY_LONG) + + # Initialize the v vectors. + call amovkl (long(1), Meml[vout], IM_MAXDIM) + do image = 1, nimages { + call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG) + call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM) + } + + # Join input images along the specified dimension. Joins along + # columns and lines require processing in special order, all others + # in the same order. In the first two cases we process all input + # images in inner loops, so we have to keep all those image + # descriptors open. + + switch (joindim) { + case 1: # join columns + nlines = 1 + do i = 2, IM_NDIM(out) + nlines = nlines * IM_LEN(out,i) + do i = 1, nlines { + stat = impnlx (out, outbuf, Meml[vout]) + cum_len = 0 + do image = 1, nimages { + in = inptr[image] + stat = imgnlx (in, inbuf, Meml[VPTR(vin,image)]) + call amovx (Memx[inbuf], Memx[outbuf+cum_len], + IM_LEN(in,1)) + cum_len = cum_len + IM_LEN(in,1) + } + } + + case 2: # join lines + nbands = 1 + do i = 3, IM_NDIM(out) + nbands = nbands * IM_LEN(out,i) + do i = 1, nbands { + do image = 1, nimages { + in = inptr[image] + do line = 1, IM_LEN(in,2) { + stat = impnlx (out, outbuf, Meml[vout]) + stat = imgnlx (in, inbuf, Meml[VPTR(vin,image)]) + call amovx (Memx[inbuf], Memx[outbuf], IM_LEN(in,1)) + } + } + } + + default: # join bands or higher + do image = 1, nimages { + in = inptr[image] + nlines = 1 + do i = 2, IM_NDIM(in) + nlines = nlines * IM_LEN(in,i) + do i = 1, nlines { + stat = impnlx (out, outbuf, Meml[vout]) + stat = imgnlx (in, inbuf, Meml[VPTR(vin,image)]) + call amovx (Memx[inbuf], Memx[outbuf], IM_LEN(in,1)) + } + } + } + + call sfree (sp) +end + + diff --git a/pkg/images/imutil/src/generic/imrep.x b/pkg/images/imutil/src/generic/imrep.x new file mode 100644 index 00000000..bcc29d0a --- /dev/null +++ b/pkg/images/imutil/src/generic/imrep.x @@ -0,0 +1,1423 @@ +include <imhdr.h> +include <mach.h> + + + +# IMREP -- Replace pixels in an image between lower and upper by value. + +procedure imreps (im, lower, upper, value, img) + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real value # Replacement value +real img # Imaginary value for complex + +pointer buf1, buf2 +int npix, junk +real ilower +short floor, ceil, newval +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnls(), impnls() + +bool fp_equalr() + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im, 1) + newval = double (value) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnls (im, buf2, v2) != EOF) + call amovks (newval, Mems[buf2], npix) + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + ceil = int (upper) + while (imgnls (im, buf1, v1) != EOF) { + junk = impnls (im, buf2, v2) + call amovs (Mems[buf1], Mems[buf2], npix) + call arles (Mems[buf2], npix, ceil, newval) + } + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + while (imgnls (im, buf1, v1) != EOF) { + junk = impnls (im, buf2, v2) + call amovs (Mems[buf1], Mems[buf2], npix) + call arges (Mems[buf2], npix, floor, newval) + } + + # Replace pixels between lower and upper by value. + } else { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + ceil = int (upper) + while (imgnls (im, buf1, v1) != EOF) { + junk = impnls (im, buf2, v2) + call amovs (Mems[buf1], Mems[buf2], npix) + call areps (Mems[buf2], npix, floor, ceil, newval) + } + } +end + + +# IMRREP -- Replace pixels in an image between lower and upper by value +# and a radius around those pixels. + +procedure imrreps (im, lower, upper, radius, value, img) + + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real radius # Radius +real value # Replacement value +real img # Imaginary value for complex + +pointer buf, buf1, buf2, ptr +int i, j, k, l, nc, nl, nradius, nbufs +real ilower +short floor, ceil, newval, val1, val2 +real radius2, y2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] # IMIO vectors +int imgnls(), impnls() +bool fp_equalr() + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + nc = IM_LEN(im, 1) + if (IM_NDIM(im) > 1) + nl = IM_LEN(im,2) + else + nl = 1 + newval = double (value) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnls (im, buf2, v2) != EOF) + call amovks (newval, Mems[buf2], nc) + return + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + floor = -MAX_SHORT + ceil = int (upper) + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + ceil = MAX_SHORT + + # Replace pixels between lower and upper by value. + } else { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + ceil = int (upper) + } + + # Initialize buffering. + radius2 = radius * radius + nradius = int (radius) + nbufs = min (1 + 2 * nradius, nl) + call calloc (buf, nc*nbufs, TY_SHORT) + + while (imgnls (im, buf1, v1) != EOF) { + j = v1[2] - 1 + buf2 = buf + mod (j, nbufs) * nc + do i = 1, nc { + val1 = Mems[buf1] + val2 = Mems[buf2] + if ((val1 >= floor) && (val1 <= ceil)) { + do k = max(1,j-nradius), min (nl,j+nradius) { + ptr = buf + mod (k, nbufs) * nc - 1 + y2 = (k - j) ** 2 + do l = max(1,i-nradius), min (nc,i+nradius) { + if ((l-i)**2 + y2 > radius2) + next + Mems[ptr+l] = INDEFS + } + } + } else { + if (!IS_INDEFS(val2)) + Mems[buf2] = val1 + } + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + + if (j > nradius) { + while (impnls (im, buf2, v2) != EOF) { + k = v2[2] - 1 + buf1 = buf + mod (k, nbufs) * nc + do i = 1, nc { + val1 = Mems[buf1] + if (IS_INDEFS(Mems[buf1])) + Mems[buf2] = newval + else + Mems[buf2] = val1 + Mems[buf1] = 0. + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + if (j != nl) + break + } + } + } + + call mfree (buf, TY_SHORT) +end + + +# AREP -- Replace array values which are between floor and ceil by value. + +procedure areps (a, npts, floor, ceil, newval) + +short a[npts] # Input arrays +int npts # Number of points +short floor, ceil # Replacement limits +short newval # Replacement value + +int i + +begin + + do i = 1, npts { + if ((a[i] >= floor) && (a[i] <= ceil)) + a[i] = newval + } +end + + +# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL. + +procedure arles (a, npts, floor, newval) + +short a[npts] +int npts +short floor, newval + +int i + +begin + + do i = 1, npts + if (a[i] <= floor) + a[i] = newval +end + + +# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL. + +procedure arges (a, npts, ceil, newval) + +short a[npts] +int npts +short ceil, newval + +int i + +begin + + do i = 1, npts + if (a[i] >= ceil) + a[i] = newval +end + + + +# IMREP -- Replace pixels in an image between lower and upper by value. + +procedure imrepi (im, lower, upper, value, img) + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real value # Replacement value +real img # Imaginary value for complex + +pointer buf1, buf2 +int npix, junk +real ilower +int floor, ceil, newval +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnli(), impnli() + +bool fp_equalr() + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im, 1) + newval = double (value) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnli (im, buf2, v2) != EOF) + call amovki (newval, Memi[buf2], npix) + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + ceil = int (upper) + while (imgnli (im, buf1, v1) != EOF) { + junk = impnli (im, buf2, v2) + call amovi (Memi[buf1], Memi[buf2], npix) + call arlei (Memi[buf2], npix, ceil, newval) + } + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + while (imgnli (im, buf1, v1) != EOF) { + junk = impnli (im, buf2, v2) + call amovi (Memi[buf1], Memi[buf2], npix) + call argei (Memi[buf2], npix, floor, newval) + } + + # Replace pixels between lower and upper by value. + } else { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + ceil = int (upper) + while (imgnli (im, buf1, v1) != EOF) { + junk = impnli (im, buf2, v2) + call amovi (Memi[buf1], Memi[buf2], npix) + call arepi (Memi[buf2], npix, floor, ceil, newval) + } + } +end + + +# IMRREP -- Replace pixels in an image between lower and upper by value +# and a radius around those pixels. + +procedure imrrepi (im, lower, upper, radius, value, img) + + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real radius # Radius +real value # Replacement value +real img # Imaginary value for complex + +pointer buf, buf1, buf2, ptr +int i, j, k, l, nc, nl, nradius, nbufs +real ilower +int floor, ceil, newval, val1, val2 +real radius2, y2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] # IMIO vectors +int imgnli(), impnli() +bool fp_equalr() + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + nc = IM_LEN(im, 1) + if (IM_NDIM(im) > 1) + nl = IM_LEN(im,2) + else + nl = 1 + newval = double (value) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnli (im, buf2, v2) != EOF) + call amovki (newval, Memi[buf2], nc) + return + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + floor = -MAX_INT + ceil = int (upper) + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + ceil = MAX_INT + + # Replace pixels between lower and upper by value. + } else { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + ceil = int (upper) + } + + # Initialize buffering. + radius2 = radius * radius + nradius = int (radius) + nbufs = min (1 + 2 * nradius, nl) + call calloc (buf, nc*nbufs, TY_INT) + + while (imgnli (im, buf1, v1) != EOF) { + j = v1[2] - 1 + buf2 = buf + mod (j, nbufs) * nc + do i = 1, nc { + val1 = Memi[buf1] + val2 = Memi[buf2] + if ((val1 >= floor) && (val1 <= ceil)) { + do k = max(1,j-nradius), min (nl,j+nradius) { + ptr = buf + mod (k, nbufs) * nc - 1 + y2 = (k - j) ** 2 + do l = max(1,i-nradius), min (nc,i+nradius) { + if ((l-i)**2 + y2 > radius2) + next + Memi[ptr+l] = INDEFI + } + } + } else { + if (!IS_INDEFI(val2)) + Memi[buf2] = val1 + } + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + + if (j > nradius) { + while (impnli (im, buf2, v2) != EOF) { + k = v2[2] - 1 + buf1 = buf + mod (k, nbufs) * nc + do i = 1, nc { + val1 = Memi[buf1] + if (IS_INDEFI(Memi[buf1])) + Memi[buf2] = newval + else + Memi[buf2] = val1 + Memi[buf1] = 0. + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + if (j != nl) + break + } + } + } + + call mfree (buf, TY_INT) +end + + +# AREP -- Replace array values which are between floor and ceil by value. + +procedure arepi (a, npts, floor, ceil, newval) + +int a[npts] # Input arrays +int npts # Number of points +int floor, ceil # Replacement limits +int newval # Replacement value + +int i + +begin + + do i = 1, npts { + if ((a[i] >= floor) && (a[i] <= ceil)) + a[i] = newval + } +end + + +# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL. + +procedure arlei (a, npts, floor, newval) + +int a[npts] +int npts +int floor, newval + +int i + +begin + + do i = 1, npts + if (a[i] <= floor) + a[i] = newval +end + + +# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL. + +procedure argei (a, npts, ceil, newval) + +int a[npts] +int npts +int ceil, newval + +int i + +begin + + do i = 1, npts + if (a[i] >= ceil) + a[i] = newval +end + + + +# IMREP -- Replace pixels in an image between lower and upper by value. + +procedure imrepl (im, lower, upper, value, img) + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real value # Replacement value +real img # Imaginary value for complex + +pointer buf1, buf2 +int npix, junk +real ilower +long floor, ceil, newval +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnll(), impnll() + +bool fp_equalr() + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im, 1) + newval = double (value) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnll (im, buf2, v2) != EOF) + call amovkl (newval, Meml[buf2], npix) + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + ceil = int (upper) + while (imgnll (im, buf1, v1) != EOF) { + junk = impnll (im, buf2, v2) + call amovl (Meml[buf1], Meml[buf2], npix) + call arlel (Meml[buf2], npix, ceil, newval) + } + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + while (imgnll (im, buf1, v1) != EOF) { + junk = impnll (im, buf2, v2) + call amovl (Meml[buf1], Meml[buf2], npix) + call argel (Meml[buf2], npix, floor, newval) + } + + # Replace pixels between lower and upper by value. + } else { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + ceil = int (upper) + while (imgnll (im, buf1, v1) != EOF) { + junk = impnll (im, buf2, v2) + call amovl (Meml[buf1], Meml[buf2], npix) + call arepl (Meml[buf2], npix, floor, ceil, newval) + } + } +end + + +# IMRREP -- Replace pixels in an image between lower and upper by value +# and a radius around those pixels. + +procedure imrrepl (im, lower, upper, radius, value, img) + + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real radius # Radius +real value # Replacement value +real img # Imaginary value for complex + +pointer buf, buf1, buf2, ptr +int i, j, k, l, nc, nl, nradius, nbufs +real ilower +long floor, ceil, newval, val1, val2 +real radius2, y2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] # IMIO vectors +int imgnll(), impnll() +bool fp_equalr() + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + nc = IM_LEN(im, 1) + if (IM_NDIM(im) > 1) + nl = IM_LEN(im,2) + else + nl = 1 + newval = double (value) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnll (im, buf2, v2) != EOF) + call amovkl (newval, Meml[buf2], nc) + return + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + floor = -MAX_LONG + ceil = int (upper) + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + ceil = MAX_LONG + + # Replace pixels between lower and upper by value. + } else { + ilower = int (lower) + if (fp_equalr(lower,ilower)) + floor = int (lower) + else + floor = int (lower+1.0) + ceil = int (upper) + } + + # Initialize buffering. + radius2 = radius * radius + nradius = int (radius) + nbufs = min (1 + 2 * nradius, nl) + call calloc (buf, nc*nbufs, TY_LONG) + + while (imgnll (im, buf1, v1) != EOF) { + j = v1[2] - 1 + buf2 = buf + mod (j, nbufs) * nc + do i = 1, nc { + val1 = Meml[buf1] + val2 = Meml[buf2] + if ((val1 >= floor) && (val1 <= ceil)) { + do k = max(1,j-nradius), min (nl,j+nradius) { + ptr = buf + mod (k, nbufs) * nc - 1 + y2 = (k - j) ** 2 + do l = max(1,i-nradius), min (nc,i+nradius) { + if ((l-i)**2 + y2 > radius2) + next + Meml[ptr+l] = INDEFL + } + } + } else { + if (!IS_INDEFL(val2)) + Meml[buf2] = val1 + } + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + + if (j > nradius) { + while (impnll (im, buf2, v2) != EOF) { + k = v2[2] - 1 + buf1 = buf + mod (k, nbufs) * nc + do i = 1, nc { + val1 = Meml[buf1] + if (IS_INDEFL(Meml[buf1])) + Meml[buf2] = newval + else + Meml[buf2] = val1 + Meml[buf1] = 0. + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + if (j != nl) + break + } + } + } + + call mfree (buf, TY_LONG) +end + + +# AREP -- Replace array values which are between floor and ceil by value. + +procedure arepl (a, npts, floor, ceil, newval) + +long a[npts] # Input arrays +int npts # Number of points +long floor, ceil # Replacement limits +long newval # Replacement value + +int i + +begin + + do i = 1, npts { + if ((a[i] >= floor) && (a[i] <= ceil)) + a[i] = newval + } +end + + +# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL. + +procedure arlel (a, npts, floor, newval) + +long a[npts] +int npts +long floor, newval + +int i + +begin + + do i = 1, npts + if (a[i] <= floor) + a[i] = newval +end + + +# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL. + +procedure argel (a, npts, ceil, newval) + +long a[npts] +int npts +long ceil, newval + +int i + +begin + + do i = 1, npts + if (a[i] >= ceil) + a[i] = newval +end + + + +# IMREP -- Replace pixels in an image between lower and upper by value. + +procedure imrepr (im, lower, upper, value, img) + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real value # Replacement value +real img # Imaginary value for complex + +pointer buf1, buf2 +int npix, junk +real floor, ceil, newval +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnlr(), impnlr() + + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im, 1) + newval = double (value) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnlr (im, buf2, v2) != EOF) + call amovkr (newval, Memr[buf2], npix) + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + ceil = double (upper) + while (imgnlr (im, buf1, v1) != EOF) { + junk = impnlr (im, buf2, v2) + call amovr (Memr[buf1], Memr[buf2], npix) + call arler (Memr[buf2], npix, ceil, newval) + } + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + floor = double (lower) + while (imgnlr (im, buf1, v1) != EOF) { + junk = impnlr (im, buf2, v2) + call amovr (Memr[buf1], Memr[buf2], npix) + call arger (Memr[buf2], npix, floor, newval) + } + + # Replace pixels between lower and upper by value. + } else { + floor = double (lower) + ceil = double (upper) + while (imgnlr (im, buf1, v1) != EOF) { + junk = impnlr (im, buf2, v2) + call amovr (Memr[buf1], Memr[buf2], npix) + call arepr (Memr[buf2], npix, floor, ceil, newval) + } + } +end + + +# IMRREP -- Replace pixels in an image between lower and upper by value +# and a radius around those pixels. + +procedure imrrepr (im, lower, upper, radius, value, img) + + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real radius # Radius +real value # Replacement value +real img # Imaginary value for complex + +pointer buf, buf1, buf2, ptr +int i, j, k, l, nc, nl, nradius, nbufs +real floor, ceil, newval, val1, val2 +real radius2, y2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] # IMIO vectors +int imgnlr(), impnlr() + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + nc = IM_LEN(im, 1) + if (IM_NDIM(im) > 1) + nl = IM_LEN(im,2) + else + nl = 1 + newval = double (value) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnlr (im, buf2, v2) != EOF) + call amovkr (newval, Memr[buf2], nc) + return + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + floor = -MAX_REAL + ceil = double (upper) + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + floor = double (lower) + ceil = MAX_REAL + + # Replace pixels between lower and upper by value. + } else { + floor = double (lower) + ceil = double (upper) + } + + # Initialize buffering. + radius2 = radius * radius + nradius = int (radius) + nbufs = min (1 + 2 * nradius, nl) + call calloc (buf, nc*nbufs, TY_REAL) + + while (imgnlr (im, buf1, v1) != EOF) { + j = v1[2] - 1 + buf2 = buf + mod (j, nbufs) * nc + do i = 1, nc { + val1 = Memr[buf1] + val2 = Memr[buf2] + if ((val1 >= floor) && (val1 <= ceil)) { + do k = max(1,j-nradius), min (nl,j+nradius) { + ptr = buf + mod (k, nbufs) * nc - 1 + y2 = (k - j) ** 2 + do l = max(1,i-nradius), min (nc,i+nradius) { + if ((l-i)**2 + y2 > radius2) + next + Memr[ptr+l] = INDEFR + } + } + } else { + if (!IS_INDEFR(val2)) + Memr[buf2] = val1 + } + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + + if (j > nradius) { + while (impnlr (im, buf2, v2) != EOF) { + k = v2[2] - 1 + buf1 = buf + mod (k, nbufs) * nc + do i = 1, nc { + val1 = Memr[buf1] + if (IS_INDEFR(Memr[buf1])) + Memr[buf2] = newval + else + Memr[buf2] = val1 + Memr[buf1] = 0. + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + if (j != nl) + break + } + } + } + + call mfree (buf, TY_REAL) +end + + +# AREP -- Replace array values which are between floor and ceil by value. + +procedure arepr (a, npts, floor, ceil, newval) + +real a[npts] # Input arrays +int npts # Number of points +real floor, ceil # Replacement limits +real newval # Replacement value + +int i + +begin + + do i = 1, npts { + if ((a[i] >= floor) && (a[i] <= ceil)) + a[i] = newval + } +end + + +# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL. + +procedure arler (a, npts, floor, newval) + +real a[npts] +int npts +real floor, newval + +int i + +begin + + do i = 1, npts + if (a[i] <= floor) + a[i] = newval +end + + +# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL. + +procedure arger (a, npts, ceil, newval) + +real a[npts] +int npts +real ceil, newval + +int i + +begin + + do i = 1, npts + if (a[i] >= ceil) + a[i] = newval +end + + + +# IMREP -- Replace pixels in an image between lower and upper by value. + +procedure imrepd (im, lower, upper, value, img) + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real value # Replacement value +real img # Imaginary value for complex + +pointer buf1, buf2 +int npix, junk +double floor, ceil, newval +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnld(), impnld() + + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im, 1) + newval = double (value) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnld (im, buf2, v2) != EOF) + call amovkd (newval, Memd[buf2], npix) + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + ceil = double (upper) + while (imgnld (im, buf1, v1) != EOF) { + junk = impnld (im, buf2, v2) + call amovd (Memd[buf1], Memd[buf2], npix) + call arled (Memd[buf2], npix, ceil, newval) + } + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + floor = double (lower) + while (imgnld (im, buf1, v1) != EOF) { + junk = impnld (im, buf2, v2) + call amovd (Memd[buf1], Memd[buf2], npix) + call arged (Memd[buf2], npix, floor, newval) + } + + # Replace pixels between lower and upper by value. + } else { + floor = double (lower) + ceil = double (upper) + while (imgnld (im, buf1, v1) != EOF) { + junk = impnld (im, buf2, v2) + call amovd (Memd[buf1], Memd[buf2], npix) + call arepd (Memd[buf2], npix, floor, ceil, newval) + } + } +end + + +# IMRREP -- Replace pixels in an image between lower and upper by value +# and a radius around those pixels. + +procedure imrrepd (im, lower, upper, radius, value, img) + + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real radius # Radius +real value # Replacement value +real img # Imaginary value for complex + +pointer buf, buf1, buf2, ptr +int i, j, k, l, nc, nl, nradius, nbufs +double floor, ceil, newval, val1, val2 +real radius2, y2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] # IMIO vectors +int imgnld(), impnld() + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + nc = IM_LEN(im, 1) + if (IM_NDIM(im) > 1) + nl = IM_LEN(im,2) + else + nl = 1 + newval = double (value) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnld (im, buf2, v2) != EOF) + call amovkd (newval, Memd[buf2], nc) + return + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + floor = -MAX_DOUBLE + ceil = double (upper) + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + floor = double (lower) + ceil = MAX_DOUBLE + + # Replace pixels between lower and upper by value. + } else { + floor = double (lower) + ceil = double (upper) + } + + # Initialize buffering. + radius2 = radius * radius + nradius = int (radius) + nbufs = min (1 + 2 * nradius, nl) + call calloc (buf, nc*nbufs, TY_DOUBLE) + + while (imgnld (im, buf1, v1) != EOF) { + j = v1[2] - 1 + buf2 = buf + mod (j, nbufs) * nc + do i = 1, nc { + val1 = Memd[buf1] + val2 = Memd[buf2] + if ((val1 >= floor) && (val1 <= ceil)) { + do k = max(1,j-nradius), min (nl,j+nradius) { + ptr = buf + mod (k, nbufs) * nc - 1 + y2 = (k - j) ** 2 + do l = max(1,i-nradius), min (nc,i+nradius) { + if ((l-i)**2 + y2 > radius2) + next + Memd[ptr+l] = INDEFD + } + } + } else { + if (!IS_INDEFD(val2)) + Memd[buf2] = val1 + } + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + + if (j > nradius) { + while (impnld (im, buf2, v2) != EOF) { + k = v2[2] - 1 + buf1 = buf + mod (k, nbufs) * nc + do i = 1, nc { + val1 = Memd[buf1] + if (IS_INDEFD(Memd[buf1])) + Memd[buf2] = newval + else + Memd[buf2] = val1 + Memd[buf1] = 0. + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + if (j != nl) + break + } + } + } + + call mfree (buf, TY_DOUBLE) +end + + +# AREP -- Replace array values which are between floor and ceil by value. + +procedure arepd (a, npts, floor, ceil, newval) + +double a[npts] # Input arrays +int npts # Number of points +double floor, ceil # Replacement limits +double newval # Replacement value + +int i + +begin + + do i = 1, npts { + if ((a[i] >= floor) && (a[i] <= ceil)) + a[i] = newval + } +end + + +# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL. + +procedure arled (a, npts, floor, newval) + +double a[npts] +int npts +double floor, newval + +int i + +begin + + do i = 1, npts + if (a[i] <= floor) + a[i] = newval +end + + +# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL. + +procedure arged (a, npts, ceil, newval) + +double a[npts] +int npts +double ceil, newval + +int i + +begin + + do i = 1, npts + if (a[i] >= ceil) + a[i] = newval +end + + + +# IMREP -- Replace pixels in an image between lower and upper by value. + +procedure imrepx (im, lower, upper, value, img) + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real value # Replacement value +real img # Imaginary value for complex + +pointer buf1, buf2 +int npix, junk +complex floor, ceil, newval +long v1[IM_MAXDIM], v2[IM_MAXDIM] +int imgnlx(), impnlx() + + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + npix = IM_LEN(im, 1) + newval = complex (value, img) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnlx (im, buf2, v2) != EOF) + call amovkx (newval, Memx[buf2], npix) + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + ceil = double (upper) + while (imgnlx (im, buf1, v1) != EOF) { + junk = impnlx (im, buf2, v2) + call amovx (Memx[buf1], Memx[buf2], npix) + call arlex (Memx[buf2], npix, ceil, newval) + } + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + floor = double (lower) + while (imgnlx (im, buf1, v1) != EOF) { + junk = impnlx (im, buf2, v2) + call amovx (Memx[buf1], Memx[buf2], npix) + call argex (Memx[buf2], npix, floor, newval) + } + + # Replace pixels between lower and upper by value. + } else { + floor = double (lower) + ceil = double (upper) + while (imgnlx (im, buf1, v1) != EOF) { + junk = impnlx (im, buf2, v2) + call amovx (Memx[buf1], Memx[buf2], npix) + call arepx (Memx[buf2], npix, floor, ceil, newval) + } + } +end + + +# IMRREP -- Replace pixels in an image between lower and upper by value +# and a radius around those pixels. + +procedure imrrepx (im, lower, upper, radius, value, img) + + +pointer im # Image descriptor +real lower, upper # Range to be replaced +real radius # Radius +real value # Replacement value +real img # Imaginary value for complex + +pointer buf, buf1, buf2, ptr +int i, j, k, l, nc, nl, nradius, nbufs +complex floor, ceil, newval, val1, val2 +real abs_floor, abs_ceil +real radius2, y2 +long v1[IM_MAXDIM], v2[IM_MAXDIM] # IMIO vectors +int imgnlx(), impnlx() + +begin + # Setup start vector for sequential reads and writes. + call amovkl (long(1), v1, IM_MAXDIM) + call amovkl (long(1), v2, IM_MAXDIM) + + nc = IM_LEN(im, 1) + if (IM_NDIM(im) > 1) + nl = IM_LEN(im,2) + else + nl = 1 + newval = complex (value, img) + + # If both lower and upper are INDEF then replace all pixels by value. + if (IS_INDEFR (lower) && IS_INDEFR (upper)) { + while (impnlx (im, buf2, v2) != EOF) + call amovkx (newval, Memx[buf2], nc) + return + + # If lower is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (lower)) { + floor = 0 + ceil = real (upper) + abs_floor = abs (floor) + abs_ceil = abs (ceil) + + # If upper is INDEF then all pixels below upper are replaced by value. + } else if (IS_INDEFR (upper)) { + floor = real (lower) + ceil = MAX_REAL + abs_floor = abs (floor) + abs_ceil = abs (ceil) + + # Replace pixels between lower and upper by value. + } else { + floor = real (lower) + ceil = real (upper) + abs_floor = abs (floor) + abs_ceil = abs (ceil) + } + + # Initialize buffering. + radius2 = radius * radius + nradius = int (radius) + nbufs = min (1 + 2 * nradius, nl) + call calloc (buf, nc*nbufs, TY_COMPLEX) + + while (imgnlx (im, buf1, v1) != EOF) { + j = v1[2] - 1 + buf2 = buf + mod (j, nbufs) * nc + do i = 1, nc { + val1 = Memx[buf1] + val2 = Memx[buf2] + if ((abs (val1) >= abs_floor) && (abs (val1) <= abs_ceil)) { + do k = max(1,j-nradius), min (nl,j+nradius) { + ptr = buf + mod (k, nbufs) * nc - 1 + y2 = (k - j) ** 2 + do l = max(1,i-nradius), min (nc,i+nradius) { + if ((l-i)**2 + y2 > radius2) + next + Memx[ptr+l] = INDEFX + } + } + } else { + if (!IS_INDEFX(val2)) + Memx[buf2] = val1 + } + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + + if (j > nradius) { + while (impnlx (im, buf2, v2) != EOF) { + k = v2[2] - 1 + buf1 = buf + mod (k, nbufs) * nc + do i = 1, nc { + val1 = Memx[buf1] + if (IS_INDEFX(Memx[buf1])) + Memx[buf2] = newval + else + Memx[buf2] = val1 + Memx[buf1] = 0. + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + if (j != nl) + break + } + } + } + + call mfree (buf, TY_COMPLEX) +end + + +# AREP -- Replace array values which are between floor and ceil by value. + +procedure arepx (a, npts, floor, ceil, newval) + +complex a[npts] # Input arrays +int npts # Number of points +complex floor, ceil # Replacement limits +complex newval # Replacement value + +int i +real abs_floor +real abs_ceil + +begin + abs_floor = abs (floor) + abs_ceil = abs (ceil) + + do i = 1, npts { + if ((abs (a[i]) >= abs_floor) && (abs (a[i]) <= abs_ceil)) + a[i] = newval + } +end + + +# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL. + +procedure arlex (a, npts, floor, newval) + +complex a[npts] +int npts +complex floor, newval + +int i +real abs_floor + +begin + abs_floor = abs (floor) + + do i = 1, npts + if (abs (a[i]) <= abs_floor) + a[i] = newval +end + + +# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL. + +procedure argex (a, npts, ceil, newval) + +complex a[npts] +int npts +complex ceil, newval + +int i +real abs_ceil + +begin + abs_ceil = abs (ceil) + + do i = 1, npts + if (abs (a[i]) >= abs_ceil) + a[i] = newval +end + + 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 + diff --git a/pkg/images/imutil/src/generic/mkpkg b/pkg/images/imutil/src/generic/mkpkg new file mode 100644 index 00000000..9878bc7b --- /dev/null +++ b/pkg/images/imutil/src/generic/mkpkg @@ -0,0 +1,21 @@ +# Make IMUTIL. + +$checkout libpkg.a ../../../ +$update libpkg.a +$checkin libpkg.a ../../../ +$exit + +libpkg.a: + imaadd.x <imhdr.h> + imadiv.x <imhdr.h> + imamax.x <imhdr.h> + imamin.x <imhdr.h> + imamul.x <imhdr.h> + imanl.x <imhdr.h> + imasub.x <imhdr.h> + imfuncs.x <imhdr.h> <mach.h> <math.h> + imjoin.x <imhdr.h> + imrep.x <imhdr.h> <mach.h> + imsum.x ../imsum.h <imhdr.h> + ; + |