aboutsummaryrefslogtreecommitdiff
path: root/noao/imred/crutil/src/t_craverage.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/imred/crutil/src/t_craverage.x')
-rw-r--r--noao/imred/crutil/src/t_craverage.x847
1 files changed, 847 insertions, 0 deletions
diff --git a/noao/imred/crutil/src/t_craverage.x b/noao/imred/crutil/src/t_craverage.x
new file mode 100644
index 00000000..f7b82113
--- /dev/null
+++ b/noao/imred/crutil/src/t_craverage.x
@@ -0,0 +1,847 @@
+include <error.h>
+include <imhdr.h>
+include <mach.h>
+
+define MAXBUF 500000 # Maximum pixel buffer
+
+define PLSIG 15.87 # Low percentile
+define PHSIG 84.13 # High percentile
+
+
+# T_CRAVERAGE -- Detect, fix, and flag cosmic rays. Also detect objects.
+# Deviant pixels relative to a local average with the candidate pixel
+# excluded and sigma are detected and replaced by the average value
+# and/or written to a cosmic ray mask. Average values above a the median
+# of a background annulus are detected as objects and cosmic rays are
+# excluded. The object positions may be output in the mask.
+
+procedure t_craverage ()
+
+int inlist # Input image list
+int outlist # Output image list
+int crlist # Output mask list
+int avglist # Output average list
+int siglist # Output sigma list
+int crval # Output cosmic ray mask value
+int objval # Output object mask value
+int navg # Averaging box size
+int nrej # Number of high pixels to reject from average
+int nbkg # Background width
+int nsig # Sigma box size
+real lobjsig, hobjsig # Object threshold sigmas
+real lcrsig, hcrsig # CR threshold sigmas outside of object
+real var0 # Variance coefficient for DN^0 term
+real var1 # Variance coefficient for DN^1 term
+real var2 # Variance coefficient for DN^2 term
+real crgrw # Cosmic ray grow radius
+real objgrw # Object grow radius
+
+int i, nc, nl, nlstep, nbox, l1, l2, l3, l4, nl1, pmmode
+pointer sp, input, output, crmask, crmask1, extname, average, sigma
+pointer in, out, pm, aim, sim
+pointer inbuf, pinbuf, outbuf, pbuf, abuf, sbuf
+
+real clgetr()
+int clgeti(), imtopenp(), imtgetim()
+pointer immap(), imgs2s(), imgs2r(), imps2r(), imps2s()
+errchk immap, imgs2s, imgs2r, imps2r, imps2s, craverage, crgrow, imgstr
+
+begin
+ call smark (sp)
+ call salloc (input, SZ_FNAME, TY_CHAR)
+ call salloc (output, SZ_FNAME, TY_CHAR)
+ call salloc (crmask, SZ_FNAME, TY_CHAR)
+ call salloc (crmask1, SZ_FNAME, TY_CHAR)
+ call salloc (average, SZ_FNAME, TY_CHAR)
+ call salloc (sigma, SZ_FNAME, TY_CHAR)
+ call salloc (extname, SZ_FNAME, TY_CHAR)
+
+ # Get parameters.
+ inlist = imtopenp ("input")
+ outlist = imtopenp ("output")
+ crlist = imtopenp ("crmask")
+ avglist = imtopenp ("average")
+ siglist = imtopenp ("sigma")
+ crval = clgeti ("crval")
+ objval = clgeti ("objval")
+ navg = max (1, clgeti ("navg") / 2)
+ nrej = min (clgeti ("nrej"), navg-1)
+ nbkg = clgeti ("nbkg")
+ nsig = clgeti ("nsig")
+ lobjsig = clgetr ("lobjsig")
+ hobjsig = clgetr ("hobjsig")
+ lcrsig = clgetr ("lcrsig")
+ hcrsig = clgetr ("hcrsig")
+ nbox = 2 * (navg + nbkg) + 1
+ var0 = clgetr ("var0")
+ var1 = clgetr ("var1")
+ var2 = clgetr ("var2")
+ crgrw = clgetr ("crgrow")
+ objgrw = clgetr ("objgrow")
+
+ # Do the input images.
+ Memc[crmask1] = EOS
+ while (imtgetim (inlist, Memc[input], SZ_FNAME) != EOF) {
+ if (imtgetim (outlist, Memc[output], SZ_FNAME) == EOF)
+ Memc[output] = EOS
+ if (imtgetim (crlist, Memc[crmask], SZ_FNAME) == EOF)
+ call strcpy (Memc[crmask1], Memc[crmask], SZ_FNAME)
+ else if (Memc[crmask] == '!')
+ call strcpy (Memc[crmask], Memc[crmask1], SZ_FNAME)
+ if (imtgetim (avglist, Memc[average], SZ_FNAME) == EOF)
+ Memc[average] = EOS
+ if (imtgetim (siglist, Memc[sigma], SZ_FNAME) == EOF)
+ Memc[sigma] = EOS
+
+ # Map the input and output images.
+ iferr {
+ in = NULL; out = NULL; pm = NULL; aim = NULL; sim = NULL
+ inbuf = NULL; pinbuf = NULL; outbuf = NULL; pbuf = NULL;
+ abuf = NULL; sbuf=NULL
+
+ in = immap (Memc[input], READ_ONLY, 0)
+ if (Memc[output] != EOS)
+ out = immap (Memc[output], NEW_COPY, in)
+ if (Memc[crmask] != EOS) {
+ if (Memc[crmask] == '!')
+ call imgstr (in, Memc[crmask+1], Memc[crmask], SZ_FNAME)
+ pmmode = READ_WRITE
+ iferr (call imgstr (in, "extname", Memc[extname], SZ_FNAME))
+ call strcpy ("pl", Memc[extname], SZ_FNAME)
+ call xt_maskname (Memc[crmask], Memc[extname], pmmode,
+ Memc[crmask], SZ_FNAME)
+ iferr (pm = immap (Memc[crmask], pmmode, 0)) {
+ pmmode = NEW_COPY
+ pm = immap (Memc[crmask], pmmode, in)
+ }
+ }
+ if (Memc[average] != EOS)
+ aim = immap (Memc[average], NEW_COPY, in)
+ if (Memc[sigma] != EOS)
+ sim = immap (Memc[sigma], NEW_COPY, in)
+
+ # Go through the input in large blocks of lines. If the
+ # block is smaller than the whole image overlap the blocks
+ # so the average only has boundaries at the ends of the image.
+ # However, the output is done in non-overlapping blocks with
+ # the pointers are adjusted so that addresses can be in the
+ # space of the input block. CRAVERAGE does not address
+ # outside of the output data block. Set the mask values
+ # based on the distances to the nearest good pixels.
+
+ nc = IM_LEN(in,1)
+ nl = IM_LEN(in,2)
+ nlstep = max (1, MAXBUF / nc - nbox)
+
+ do i = 1, nl, nlstep {
+ l1 = i
+ l2 = min (nl, i + nlstep - 1)
+ l3 = max (1, l1 - nbox / 2)
+ l4 = min (nl, l2 + nbox / 2)
+ nl1 = l4 - l3 + 1
+ inbuf = imgs2r (in, 1, nc, l3, l4)
+ if (out != NULL)
+ outbuf = imps2r (out, 1, nc, l1, l2) - (l1 - l3) * nc
+ if (pm != NULL) {
+ if (pmmode == READ_WRITE) {
+ pinbuf = imgs2s (pm, 1, nc, l3, l4)
+ pbuf = imps2s (pm, 1, nc, l1, l2)
+ call amovs (Mems[pinbuf+(l1-l3)*nc],
+ Mems[pbuf], nc*(l2-l1+1))
+ pbuf = pbuf - (l1 - l3) * nc
+ } else {
+ pinbuf = NULL
+ pbuf = imps2s (pm, 1, nc, l1, l2)
+ call aclrs (Mems[pbuf], nc*(l2-l1+1))
+ pbuf = pbuf - (l1 - l3) * nc
+ }
+ }
+ if (aim != NULL)
+ abuf = imps2r (aim, 1, nc, l1, l2) - (l1 - l3) * nc
+ if (sim != NULL)
+ sbuf = imps2r (sim, 1, nc, l1, l2) - (l1 - l3) * nc
+ if (pinbuf == NULL)
+ call craverage (inbuf, outbuf, pbuf, abuf, sbuf,
+ nc, nl1, l1-l3+1, l2-l3+1, navg, nrej, nbkg,
+ var0, var1, var2, nsig, lcrsig, hcrsig,
+ lobjsig, hobjsig, crval, objval)
+ else
+ call craverage1 (inbuf, pinbuf, outbuf, pbuf, abuf,
+ sbuf, nc, nl1, l1-l3+1, l2-l3+1, navg, nrej, nbkg,
+ var0, var1, var2, nsig, lcrsig, hcrsig,
+ lobjsig, hobjsig, crval, objval)
+ }
+
+ # Grow regions if desired. The routines are nops if the
+ # grow is zero.
+
+ if (pm != NULL) {
+ if (pmmode != READ_WRITE) {
+ call imunmap (pm)
+ iferr (pm = immap (Memc[crmask], READ_WRITE, 0))
+ call error (1, "Can't reopen mask for growing")
+ }
+
+ if (crval == objval)
+ call crgrow (pm, max (crgrw, objgrw), crval, crval)
+ else {
+ call crgrow (pm, crgrw, crval, crval)
+ call crgrow (pm, objgrw, objval, objval)
+ }
+ }
+ } then
+ call erract (EA_WARN)
+
+ if (sim != NULL)
+ call imunmap (sim)
+ if (aim != NULL)
+ call imunmap (aim)
+ if (pm != NULL)
+ call imunmap (pm)
+ if (out != NULL)
+ call imunmap (out)
+ call imunmap (in)
+ }
+
+ call imtclose (inlist)
+ call imtclose (outlist)
+ call imtclose (crlist)
+ call imtclose (avglist)
+ call imtclose (siglist)
+
+ call sfree (sp)
+end
+
+
+# CRAVERAGE -- Detect, replace, and flag cosmic rays.
+# A local background is computed using moving box averages to avoid
+# contaminating bad pixels. If variance model is given then that is
+# used otherwise a local sigma is computed in blocks (it is not a moving box
+# for efficiency) by using a percentile point of the sorted pixel values to
+# estimate the width of the distribution uncontaminated by bad pixels). Once
+# the background and sigma are known deviant pixels are found by using sigma
+# threshold factors.
+
+procedure craverage (in, out, pout, aout, sout, nc, nl, nl1, nl2,
+ navg, nrej, nbkg, var0, var1, var2, nsig, lcrsig, hcrsig,
+ lobjsig, hobjsig crval, objval)
+
+pointer in #I Input data
+pointer out #O Output data
+pointer pout #O Output mask (0=good, 1=bad)
+pointer aout #O Output averages
+pointer sout #O Output sigmas
+int nc, nl #I Number of columns and lines
+int nl1, nl2 #I Lines to compute
+int navg #I Averaging box half-size
+int nrej #I Number of high pixels to reject from average
+int nbkg #I Median background width
+real var0 #I Variance coefficient for DN^0 term
+real var1 #I Variance coefficient for DN^1 term
+real var2 #I Variance coefficient for DN^2 term
+int nsig #I Sigma box size
+real lcrsig, hcrsig #I Threshold sigmas outside of object
+real lobjsig, hobjsig #I Object threshold sigmas
+int crval #I CR mask value
+int objval #I Object mask value
+
+int i, j, c, c1, c2, c3, c4, l, l1, l2, l3, l4, n1, n2
+int navg2, nbkg2, nsig2, plsig, phsig
+real data, avg, bkg, sigma, losig, hosig
+real low, high, cravg(), amedr()
+pointer stack, avgs, bkgs, sigs, work1, work2
+pointer ptr1, ptr2, ip, op, pp, ap, sp
+
+begin
+ navg2 = (2 * navg + 1) ** 2
+ nbkg2 = (2 * (navg + nbkg) + 1) ** 2 - navg2
+ nsig2 = nsig * nsig
+
+ call smark (stack)
+ call salloc (avgs, nc, TY_REAL)
+ call salloc (bkgs, nc, TY_REAL)
+ call salloc (sigs, nc, TY_REAL)
+ call salloc (work1, navg2, TY_REAL)
+ call salloc (work2, max (nsig2, nbkg2), TY_REAL)
+
+ if (var0 != 0. && var1 == 0. && var2 ==0.)
+ call amovkr (sqrt(var0), Memr[sigs], nc)
+
+ avgs = avgs - 1
+ sigs = sigs - 1
+ bkgs = bkgs - 1
+
+ plsig = nint (PLSIG*nsig2/100.-1)
+ phsig = nint (PHSIG*nsig2/100.-1)
+ losig = lobjsig / sqrt (real(navg2-1))
+ hosig = hobjsig / sqrt (real(navg2-1))
+
+ do l = nl1, nl2 {
+ # Compute statistics.
+ l1 = max (1, l-navg-nbkg)
+ l2 = max (1, l-navg)
+ l3 = min (nl, l+navg)
+ l4 = min (nl, l+navg+nbkg)
+ ap = aout + (l - 1) * nc
+ do c = 1, nc {
+ c1 = max (1, c-navg-nbkg)
+ c2 = max (1, c-navg)
+ c3 = min (nc, c+navg)
+ c4 = min (nc, c+navg+nbkg)
+ ptr1 = work1
+ ptr2 = work2
+ n1 = 0
+ n2 = 0
+ do j = l1, l2-1 {
+ ip = in + (j - 1) * nc + c1 - 1
+ do i = c1, c4 {
+ Memr[ptr2] = Memr[ip]
+ n2 = n2 + 1
+ ptr2 = ptr2 + 1
+ ip = ip + 1
+ }
+ }
+ do j = l2, l3 {
+ ip = in + (j - 1) * nc + c1 - 1
+ do i = c1, c2-1 {
+ Memr[ptr2] = Memr[ip]
+ n2 = n2 + 1
+ ptr2 = ptr2 + 1
+ ip = ip + 1
+ }
+ do i = c2, c3 {
+ if (j != l || i != c) {
+ Memr[ptr1] = Memr[ip]
+ n1 = n1 + 1
+ ptr1 = ptr1 + 1
+ }
+ ip = ip + 1
+ }
+ do i = c3+1, c4 {
+ Memr[ptr2] = Memr[ip]
+ n2 = n2 + 1
+ ptr2 = ptr2 + 1
+ ip = ip + 1
+ }
+ }
+ do j = l3+1, l4 {
+ ip = in + (j - 1) * nc + c1 - 1
+ do i = c1, c4 {
+ Memr[ptr2] = Memr[ip]
+ n2 = n2 + 1
+ ptr2 = ptr2 + 1
+ ip = ip + 1
+ }
+ }
+ avg = cravg (Memr[work1], n1, nrej)
+ bkg = amedr (Memr[work2], n2)
+ Memr[bkgs+c] = bkg
+ Memr[avgs+c] = avg
+ if (aout != NULL) {
+ Memr[ap] = avg - bkg
+ ap = ap + 1
+ }
+ }
+
+ # Compute sigmas and output if desired.
+ if (var0 != 0. || var1 != 0. || var2 != 0.) {
+ if (var1 != 0.) {
+ if (var2 != 0.) {
+ do c = 1, nc {
+ data = max (0., Memr[avgs+c])
+ Memr[sigs+c] = sqrt (var0+var1*data+var2*data**2)
+ }
+ } else {
+ do c = 1, nc {
+ data = max (0., Memr[avgs+c])
+ Memr[sigs+c] = sqrt (var0 + var1 * data)
+ }
+ }
+ } else if (var2 != 0.) {
+ do c = 1, nc {
+ data = max (0., Memr[avgs+c])
+ Memr[sigs+c] = sqrt (var0 + var2 * data**2)
+ }
+ }
+ } else {
+ # Compute sigmas from percentiles. This is done in blocks.
+ if (mod (l-nl1, nsig) == 0 && l<nl-nsig+1) {
+ do c = 1, nc-nsig+1, nsig {
+ ptr2 = work2
+ n2 = 0
+ do j = l, l+nsig-1 {
+ ip = in + (j - 1) * nc + c - 1
+ do i = 1, nsig {
+ Memr[ptr2] = Memr[ip]
+ n2 = n2 + 1
+ ptr2 = ptr2 + 1
+ ip = ip + 1
+ }
+ }
+ call asrtr (Memr[work2], Memr[work2], n2)
+ sigma = (Memr[work2+phsig] - Memr[work2+plsig]) / 2.
+ call amovkr (sigma, Memr[sigs+c], nsig)
+ }
+ call amovkr (sigma, Memr[sigs+c], nc-c+1)
+ }
+ }
+ if (sout != NULL) {
+ sp = sout + (l - 1) * nc
+ do c = 1, nc {
+ Memr[sp] = Memr[sigs+c]
+ sp = sp + 1
+ }
+ }
+
+ # Detect, fix, and flag cosmic rays.
+ if (pout == NULL && out == NULL)
+ ;
+ else if (pout == NULL) {
+ ip = in + (l - 1) * nc
+ op = out + (l - 1) * nc
+ do c = 1, nc {
+ data = Memr[ip]
+ avg = Memr[avgs+c]
+ bkg = Memr[bkgs+c]
+ sigma = Memr[sigs+c]
+ low = bkg - losig * sigma
+ high = bkg + hosig * sigma
+ if (avg < low || avg > high) {
+ Memr[op] = data
+ } else {
+ low = avg - lcrsig * sigma
+ high = avg + hcrsig * sigma
+ if (data < low || data > high)
+ Memr[op] = avg
+ else
+ Memr[op] = data
+ }
+ ip = ip + 1
+ op = op + 1
+ }
+ } else if (out == NULL) {
+ ip = in + (l - 1) * nc
+ pp = pout + (l - 1) * nc
+ do c = 1, nc {
+ data = Memr[ip]
+ avg = Memr[avgs+c]
+ bkg = Memr[bkgs+c]
+ sigma = Memr[sigs+c]
+ low = bkg - losig * sigma
+ high = bkg + hosig * sigma
+ if (avg < low || avg > high)
+ Mems[pp] = objval
+ else {
+ low = avg - lcrsig * sigma
+ high = avg + hcrsig * sigma
+ if (data < low || data > high)
+ Mems[pp] = crval
+ }
+ ip = ip + 1
+ pp = pp + 1
+ }
+ } else {
+ ip = in + (l - 1) * nc
+ op = out + (l - 1) * nc
+ pp = pout + (l - 1) * nc
+ do c = 1, nc {
+ data = Memr[ip]
+ avg = Memr[avgs+c]
+ bkg = Memr[bkgs+c]
+ sigma = Memr[sigs+c]
+ low = bkg - losig * sigma
+ high = bkg + hosig * sigma
+ if (avg < low || avg > high) {
+ Memr[op] = data
+ Mems[pp] = objval
+ } else {
+ low = avg - lcrsig * sigma
+ high = avg + hcrsig * sigma
+ if (data < low || data > high) {
+ Memr[op] = avg
+ Mems[pp] = crval
+ } else
+ Memr[op] = data
+ }
+ ip = ip + 1
+ op = op + 1
+ pp = pp + 1
+ }
+ }
+ }
+
+ call sfree (stack)
+end
+
+
+# CRAVERAGE1 -- Detect, replace, and flag cosmic rays checking input mask.
+# A local background is computed using moving box averages to avoid
+# contaminating bad pixels. If variance model is given then that is
+# used otherwise a local sigma is computed in blocks (it is not a moving box
+# for efficiency) by using a percentile point of the sorted pixel values to
+# estimate the width of the distribution uncontaminated by bad pixels). Once
+# the background and sigma are known deviant pixels are found by using sigma
+# threshold factors.
+
+procedure craverage1 (in, pin, out, pout, aout, sout, nc, nl, nl1, nl2,
+ navg, nrej, nbkg, var0, var1, var2, nsig, lcrsig, hcrsig,
+ lobjsig, hobjsig crval, objval)
+
+pointer in #I Input data
+pointer pin #I Pixel mask data
+pointer out #O Output data
+pointer pout #O Output mask (0=good, 1=bad)
+pointer aout #O Output averages
+pointer sout #O Output sigmas
+int nc, nl #I Number of columns and lines
+int nl1, nl2 #I Lines to compute
+int navg #I Averaging box half-size
+int nrej #I Number of high pixels to reject from average
+int nbkg #I Median background width
+real var0 #I Variance coefficient for DN^0 term
+real var1 #I Variance coefficient for DN^1 term
+real var2 #I Variance coefficient for DN^2 term
+int nsig #I Sigma box size
+real lcrsig, hcrsig #I Threshold sigmas outside of object
+real lobjsig, hobjsig #I Object threshold sigmas
+int crval #I CR mask value
+int objval #I Object mask value
+
+int i, j, c, c1, c2, c3, c4, l, l1, l2, l3, l4, n1, n2
+int navg2, nbkg2, nsig2, plsig, phsig
+real data, avg, bkg, sigma, losig, hosig
+real low, high, cravg(), amedr()
+pointer stack, avgs, bkgs, sigs, work1, work2
+pointer ptr1, ptr2, ip, mp, op, pp, ap, sp
+
+begin
+ navg2 = (2 * navg + 1) ** 2
+ nbkg2 = (2 * (navg + nbkg) + 1) ** 2 - navg2
+ nsig2 = nsig * nsig
+
+ call smark (stack)
+ call salloc (avgs, nc, TY_REAL)
+ call salloc (bkgs, nc, TY_REAL)
+ call salloc (sigs, nc, TY_REAL)
+ call salloc (work1, navg2, TY_REAL)
+ call salloc (work2, max (nsig2, nbkg2), TY_REAL)
+
+ if (var0 != 0. && var1 == 0. && var2 ==0.)
+ call amovkr (sqrt(var0), Memr[sigs], nc)
+
+ avgs = avgs - 1
+ sigs = sigs - 1
+ bkgs = bkgs - 1
+
+ losig = lobjsig / sqrt (real(navg2-1))
+ hosig = hobjsig / sqrt (real(navg2-1))
+
+ do l = nl1, nl2 {
+ # Compute statistics.
+ l1 = max (1, l-navg-nbkg)
+ l2 = max (1, l-navg)
+ l3 = min (nl, l+navg)
+ l4 = min (nl, l+navg+nbkg)
+ ap = aout + (l - 1) * nc
+ do c = 1, nc {
+ c1 = max (1, c-navg-nbkg)
+ c2 = max (1, c-navg)
+ c3 = min (nc, c+navg)
+ c4 = min (nc, c+navg+nbkg)
+ ptr1 = work1
+ ptr2 = work2
+ n1 = 0
+ n2 = 0
+ do j = l1, l2-1 {
+ ip = in + (j - 1) * nc + c1 - 1
+ mp = pin + (j - 1) * nc + c1 - 1
+ do i = c1, c4 {
+ if (Mems[mp] == 0) {
+ Memr[ptr2] = Memr[ip]
+ n2 = n2 + 1
+ ptr2 = ptr2 + 1
+ }
+ ip = ip + 1
+ mp = mp + 1
+ }
+ }
+ do j = l2, l3 {
+ ip = in + (j - 1) * nc + c1 - 1
+ mp = pin + (j - 1) * nc + c1 - 1
+ do i = c1, c2-1 {
+ if (Mems[mp] == 0) {
+ Memr[ptr2] = Memr[ip]
+ n2 = n2 + 1
+ ptr2 = ptr2 + 1
+ }
+ ip = ip + 1
+ mp = mp + 1
+ }
+ do i = c2, c3 {
+ if ((j != l || i != c) && Mems[mp] == 0) {
+ Memr[ptr1] = Memr[ip]
+ n1 = n1 + 1
+ ptr1 = ptr1 + 1
+ }
+ ip = ip + 1
+ mp = mp + 1
+ }
+ do i = c3+1, c4 {
+ if (Mems[mp] == 0) {
+ Memr[ptr2] = Memr[ip]
+ n2 = n2 + 1
+ ptr2 = ptr2 + 1
+ }
+ ip = ip + 1
+ mp = mp + 1
+ }
+ }
+ do j = l3+1, l4 {
+ ip = in + (j - 1) * nc + c1 - 1
+ mp = pin + (j - 1) * nc + c1 - 1
+ do i = c1, c4 {
+ if (Mems[mp] == 0) {
+ Memr[ptr2] = Memr[ip]
+ n2 = n2 + 1
+ ptr2 = ptr2 + 1
+ }
+ ip = ip + 1
+ }
+ }
+ if (n1 > 0)
+ avg = cravg (Memr[work1], n1, nrej)
+ else
+ avg = INDEFR
+ if (n2 > 0)
+ bkg = amedr (Memr[work2], n2)
+ else
+ bkg = INDEFR
+ Memr[bkgs+c] = bkg
+ Memr[avgs+c] = avg
+ if (aout != NULL) {
+ if (IS_INDEFR(avg) || IS_INDEFR(bkg))
+ Memr[ap] = 0.
+ else
+ Memr[ap] = avg - bkg
+ ap = ap + 1
+ }
+ }
+
+ # Compute sigmas and output if desired.
+ if (var0 != 0. || var1 != 0. || var2 != 0.) {
+ if (var1 != 0.) {
+ if (var2 != 0.) {
+ do c = 1, nc {
+ data = max (0., Memr[avgs+c])
+ Memr[sigs+c] = sqrt (var0+var1*data+var2*data**2)
+ }
+ } else {
+ do c = 1, nc {
+ data = max (0., Memr[avgs+c])
+ Memr[sigs+c] = sqrt (var0 + var1 * data)
+ }
+ }
+ } else if (var2 != 0.) {
+ do c = 1, nc {
+ data = max (0., Memr[avgs+c])
+ Memr[sigs+c] = sqrt (var0 + var2 * data**2)
+ }
+ }
+ } else {
+ # Compute sigmas from percentiles. This is done in blocks.
+ if (mod (l-nl1, nsig) == 0 && l<nl-nsig+1) {
+ do c = 1, nc-nsig+1, nsig {
+ ptr2 = work2
+ n2 = 0
+ do j = l, l+nsig-1 {
+ ip = in + (j - 1) * nc + c - 1
+ mp = pin + (j - 1) * nc + c - 1
+ do i = 1, nsig {
+ if (Mems[mp] == 0) {
+ Memr[ptr2] = Memr[ip]
+ n2 = n2 + 1
+ ptr2 = ptr2 + 1
+ }
+ ip = ip + 1
+ mp = mp + 1
+ }
+ }
+ if (n2 > 10) {
+ call asrtr (Memr[work2], Memr[work2], n2)
+ plsig = nint (PLSIG*n2/100.-1)
+ phsig = nint (PHSIG*n2/100.-1)
+ sigma = (Memr[work2+phsig]-Memr[work2+plsig])/2.
+ } else
+ sigma = INDEFR
+ call amovkr (sigma, Memr[sigs+c], nsig)
+ }
+ call amovkr (sigma, Memr[sigs+c], nc-c+1)
+ }
+ }
+ if (sout != NULL) {
+ sp = sout + (l - 1) * nc
+ do c = 1, nc {
+ sigma = Memr[sigs+c]
+ if (IS_INDEFR(sigma))
+ Memr[sp] = 0.
+ else
+ Memr[sp] = sigma
+ sp = sp + 1
+ }
+ }
+
+ # Detect, fix, and flag cosmic rays.
+ if (pout == NULL && out == NULL)
+ ;
+ if (pout == NULL) {
+ ip = in + (l - 1) * nc
+ mp = pin + (l - 1) * nc
+ op = out + (l - 1) * nc
+ do c = 1, nc {
+ data = Memr[ip]
+ avg = Memr[avgs+c]
+ bkg = Memr[bkgs+c]
+ sigma = Memr[sigs+c]
+ if (!(Mems[mp] != 0 || IS_INDEFR(avg) ||
+ IS_INDEFR(bkg) || IS_INDEFR(sigma))) {
+ low = bkg - losig * sigma
+ high = bkg + hosig * sigma
+ if (avg < low || avg > high) {
+ Memr[op] = data
+ } else {
+ low = avg - lcrsig * sigma
+ high = avg + hcrsig * sigma
+ if (data < low || data > high)
+ Memr[op] = avg
+ else
+ Memr[op] = data
+ }
+ } else
+ Memr[op] = data
+ ip = ip + 1
+ mp = mp + 1
+ op = op + 1
+ }
+ } else if (out == NULL) {
+ ip = in + (l - 1) * nc
+ mp = pin + (l - 1) * nc
+ pp = pout + (l - 1) * nc
+ do c = 1, nc {
+ data = Memr[ip]
+ avg = Memr[avgs+c]
+ bkg = Memr[bkgs+c]
+ sigma = Memr[sigs+c]
+ if (!(Mems[mp] != 0 || IS_INDEFR(avg) ||
+ IS_INDEFR(bkg) || IS_INDEFR(sigma))) {
+ low = bkg - losig * sigma
+ high = bkg + hosig * sigma
+ if (avg < low || avg > high)
+ Mems[pp] = objval
+ else {
+ low = avg - lcrsig * sigma
+ high = avg + hcrsig * sigma
+ if (data < low || data > high)
+ Mems[pp] = crval
+ }
+ }
+ ip = ip + 1
+ mp = mp + 1
+ pp = pp + 1
+ }
+ } else {
+ ip = in + (l - 1) * nc
+ mp = pin + (l - 1) * nc
+ op = out + (l - 1) * nc
+ pp = pout + (l - 1) * nc
+ do c = 1, nc {
+ data = Memr[ip]
+ avg = Memr[avgs+c]
+ bkg = Memr[bkgs+c]
+ sigma = Memr[sigs+c]
+ if (!(Mems[mp] != 0 || IS_INDEFR(avg) ||
+ IS_INDEFR(bkg) || IS_INDEFR(sigma))) {
+ low = bkg - losig * sigma
+ high = bkg + hosig * sigma
+ if (avg < low || avg > high) {
+ Memr[op] = data
+ Mems[pp] = objval
+ } else {
+ low = avg - lcrsig * sigma
+ high = avg + hcrsig * sigma
+ if (data < low || data > high) {
+ Memr[op] = avg
+ Mems[pp] = crval
+ } else
+ Memr[op] = data
+ }
+ } else
+ Memr[op] = data
+ ip = ip + 1
+ mp = mp + 1
+ op = op + 1
+ pp = pp + 1
+ }
+ }
+ }
+
+ call sfree (stack)
+end
+
+
+# CRAVG -- Compute average with the highest nrej points excluded.
+# When nrej is greater than 2 the data array will be returned sorted.
+
+real procedure cravg (data, npts, nrej)
+
+real data[npts] #I Input data (will be sorted if nrej>2)
+int npts #I Number of data points
+int nrej #I Number of data points to reject
+
+int i
+real sum, max1, max2, val
+
+begin
+ if (npts <= nrej)
+ return (INDEFR)
+
+ switch (nrej) {
+ case 0:
+ sum = 0.
+ do i = 1, npts
+ sum = sum + data[i]
+ case 1:
+ sum = 0.
+ max1 = data[1]
+ do i = 2, npts {
+ val = data[i]
+ if (val > max1) {
+ sum = sum + max1
+ max1 = val
+ } else
+ sum = sum + val
+ }
+ case 2:
+ sum = 0.
+ max1 = min (data[1], data[2])
+ max2 = max (data[1], data[2])
+ do i = 3, npts {
+ val = data[i]
+ if (val > max1) {
+ sum = sum + max1
+ if (val > max2) {
+ max1 = max2
+ max2 = val
+ } else
+ max1 = val
+ } else
+ sum = sum + val
+ }
+ default:
+ call asrtr (data, data, npts)
+ sum = 0.
+ do i = 1, npts-nrej
+ sum = sum + data[i]
+ }
+
+ return (sum / (npts - nrej))
+end
generated by cgit (git 2.34.1) at 2025-09-20 07:12:01 -0400