From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- noao/onedspec/scombine/generic/icaclip.x | 555 +++++++++++++++++++++++++++++ noao/onedspec/scombine/generic/icaverage.x | 84 +++++ noao/onedspec/scombine/generic/iccclip.x | 453 +++++++++++++++++++++++ noao/onedspec/scombine/generic/icgrow.x | 76 ++++ noao/onedspec/scombine/generic/icmedian.x | 139 ++++++++ noao/onedspec/scombine/generic/icmm.x | 152 ++++++++ noao/onedspec/scombine/generic/icpclip.x | 224 ++++++++++++ noao/onedspec/scombine/generic/icsclip.x | 486 +++++++++++++++++++++++++ noao/onedspec/scombine/generic/icsort.x | 275 ++++++++++++++ noao/onedspec/scombine/generic/mkpkg | 16 + 10 files changed, 2460 insertions(+) create mode 100644 noao/onedspec/scombine/generic/icaclip.x create mode 100644 noao/onedspec/scombine/generic/icaverage.x create mode 100644 noao/onedspec/scombine/generic/iccclip.x create mode 100644 noao/onedspec/scombine/generic/icgrow.x create mode 100644 noao/onedspec/scombine/generic/icmedian.x create mode 100644 noao/onedspec/scombine/generic/icmm.x create mode 100644 noao/onedspec/scombine/generic/icpclip.x create mode 100644 noao/onedspec/scombine/generic/icsclip.x create mode 100644 noao/onedspec/scombine/generic/icsort.x create mode 100644 noao/onedspec/scombine/generic/mkpkg (limited to 'noao/onedspec/scombine/generic') diff --git a/noao/onedspec/scombine/generic/icaclip.x b/noao/onedspec/scombine/generic/icaclip.x new file mode 100644 index 00000000..41432dd7 --- /dev/null +++ b/noao/onedspec/scombine/generic/icaclip.x @@ -0,0 +1,555 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + +define MINCLIP 3 # Minimum number of images for this algorithm + + +# IC_AAVSIGCLIP -- Reject pixels using an average sigma about the average +# The average sigma is normalized by the expected poisson sigma. + +procedure ic_aavsigclipr (d, m, n, scales, zeros, nimages, npts, average) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +int nimages # Number of images +int npts # Number of output points per line +real average[npts] # Average + +int i, j, k, l, jj, n1, n2, nin, nk, maxkeep +real d1, low, high, sum, a, s, s1, r, one +data one /1.0/ +pointer sp, sums, resid, dp1, dp2, mp1, mp2 + +include "../icombine.com" + +begin + # If there are insufficient pixels go on to the combining. + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + call smark (sp) + call salloc (sums, npts, TY_REAL) + call salloc (resid, nimages+1, TY_REAL) + + # Since the unweighted average is computed here possibly skip combining + if (dowts || combine != AVERAGE) + docombine = true + else + docombine = false + + # Compute the unweighted average with the high and low rejected and + # the poisson scaled average sigma. There must be at least three + # pixels at each point to define the average and contributions to + # the mean sigma. Corrections for differences in the image + # scale factors are selected by the doscale1 flag. + + nin = n[1] + s = 0. + n2 = 0 + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (n1 < 3) + next + + # Unweighted average with the high and low rejected + low = Memr[d[1]+k] + high = Memr[d[2]+k] + if (low > high) { + d1 = low + low = high + high = d1 + } + sum = 0. + do j = 3, n1 { + d1 = Memr[d[j]+k] + if (d1 < low) { + sum = sum + low + low = d1 + } else if (d1 > high) { + sum = sum + high + high = d1 + } else + sum = sum + d1 + } + a = sum / (n1 - 2) + sum = sum + low + high + + # Poisson scaled sigma accumulation + if (doscale1) { + do j = 1, n1 { + dp1 = d[j] + k + mp1 = m[j] + k + + d1 = Memr[dp1] + l = Memi[mp1] + s1 = max (one, (a + zeros[l]) / scales[l]) + s = s + (d1 - a) ** 2 / s1 + } + } else { + s1 = max (one, a) + do j = 1, n1 + s = s + (Memr[d[j]+k] - a) ** 2 / s1 + } + n2 = n2 + n1 + + # Save the average and sum for later. + average[i] = a + Memr[sums+k] = sum + } + + # Here is the final sigma. + if (n2 > 1) + s = sqrt (s / (n2 - 1)) + + # Reject pixels and compute the final average (if needed). + # There must be at least three pixels at each point for rejection. + # Iteratively scale the mean sigma and reject pixels + # Compact the data and keep track of the image IDs if needed. + + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + if (n1 <= max (2, maxkeep)) { + if (!docombine) { + if (n1 == 0) + average[i] = blank + else { + sum = Memr[d[1]+k] + do j = 2, n1 + sum = sum + Memr[d[j]+k] + average[i] = sum / n1 + } + } + next + } + + a = average[i] + sum = Memr[sums+k] + + repeat { + n2 = n1 + if (s > 0.) { + if (doscale1) { + for (j=1; j<=n1; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + + d1 = Memr[dp1] + l = Memi[mp1] + s1 = s * sqrt (max (one, (a+zeros[l]) / scales[l])) + r = (d1 - a) / s1 + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs(r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + mp2 = m[n1] + k + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } else { + s1 = s * sqrt (max (one, a)) + for (j=1; j<=n1; j=j+1) { + dp1 = d[j] + k + d1 = Memr[dp1] + r = (d1 - a) / s1 + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs(r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[n1] + k + l = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + } + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } + } + if (n1 > 1) + a = sum / n1 + } until (n1 == n2 || n1 <= max (2, maxkeep)) + + # If too many are rejected add some back in. + # Pixels with equal residuals are added together. + if (n1 < maxkeep) { + nk = maxkeep + if (doscale1) { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } else { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } + if (n1 > 1) + a = sum / n1 + } + + # Save the average if needed. + n[i] = n1 + if (!docombine) { + if (n1 > 0) + average[i] = a + else + average[i] = blank + } + } + + # Check if the data flag has to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + call sfree (sp) +end + + +# IC_MAVSIGCLIP -- Reject pixels using an average sigma about the median +# The average sigma is normalized by the expected poisson sigma. + +procedure ic_mavsigclipr (d, m, n, scales, zeros, nimages, npts, median) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +int nimages # Number of images +int npts # Number of output points per line +real median[npts] # Median + +int i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep +pointer sp, resid, mp1, mp2 +real med, low, high, r, s, s1, one +data one /1.0/ + +include "../icombine.com" + +begin + # If there are insufficient pixels go on to the combining. + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + + # Compute the poisson scaled average sigma about the median. + # There must be at least three pixels at each point to define + # the mean sigma. Corrections for differences in the image + # scale factors are selected by the doscale1 flag. + + s = 0. + n2 = 0 + nin = n[1] + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (n1 < 3) { + if (n1 == 0) + median[i] = blank + else if (n1 == 1) + median[i] = Memr[d[1]+k] + else { + low = Memr[d[1]+k] + high = Memr[d[2]+k] + median[i] = (low + high) / 2. + } + next + } + + # Median + n3 = 1 + n1 / 2 + if (mod (n1, 2) == 0) { + low = Memr[d[n3-1]+k] + high = Memr[d[n3]+k] + med = (low + high) / 2. + } else + med = Memr[d[n3]+k] + + # Poisson scaled sigma accumulation + if (doscale1) { + do j = 1, n1 { + l = Memi[m[j]+k] + s1 = max (one, (med + zeros[l]) / scales[l]) + s = s + (Memr[d[j]+k] - med) ** 2 / s1 + } + } else { + s1 = max (one, med) + do j = 1, n1 + s = s + (Memr[d[j]+k] - med) ** 2 / s1 + } + n2 = n2 + n1 + + # Save the median for later. + median[i] = med + } + + # Here is the final sigma. + if (n2 > 1) + s = sqrt (s / (n2 - 1)) + else + return + + # Compute individual sigmas and iteratively clip. + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + if (n1 < max (3, maxkeep+1)) + next + nl = 1 + nh = n1 + med = median[i] + + repeat { + n2 = n1 + n3 = nl + n1 / 2 + + if (n1 >= max (MINCLIP, maxkeep+1) && s > 0.) { + if (doscale1) { + for (; nl <= n2; nl = nl + 1) { + l = Memi[m[nl]+k] + s1 = s * sqrt (max (one, (med+zeros[l])/scales[l])) + r = (med - Memr[d[nl]+k]) / s1 + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + l = Memi[m[nh]+k] + s1 = s * sqrt (max (one, (med+zeros[l])/scales[l])) + r = (Memr[d[nh]+k] - med) / s1 + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } else { + s1 = s * sqrt (max (one, med)) + for (; nl <= n2; nl = nl + 1) { + r = (med - Memr[d[nl]+k]) / s1 + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + r = (Memr[d[nh]+k] - med) / s1 + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } + + # Recompute median + if (n1 < n2) { + if (n1 > 0) { + n3 = nl + n1 / 2 + if (mod (n1, 2) == 0) { + low = Memr[d[n3-1]+k] + high = Memr[d[n3]+k] + med = (low + high) / 2. + } else + med = Memr[d[n3]+k] + } else + med = blank + } + } + } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1)) + + # If too many are rejected add some back in. + # Pixels with equal residuals are added together. + while (n1 < maxkeep) { + if (nl == 1) + nh = nh + 1 + else if (nh == n[i]) + nl = nl - 1 + else { + r = Memr[resid+nl-1] + s = Memr[resid+nh+1] + if (r < s) { + nl = nl - 1 + r = r + TOL + if (s <= r) + nh = nh + 1 + if (nl > 1) { + if (Memr[resid+nl-1] <= r) + nl = nl - 1 + } + } else { + nh = nh + 1 + s = s + TOL + if (r <= s) + nl = nl - 1 + if (nh < n2) { + if (Memr[resid+nh+1] <= s) + nh = nh + 1 + } + } + } + n1 = nh - nl + 1 + + # Recompute median + if (n1 < n2) { + if (n1 > 0) { + n3 = nl + n1 / 2 + if (mod (n1, 2) == 0) { + low = Memr[d[n3-1]+k] + high = Memr[d[n3]+k] + med = (low + high) / 2. + } else + med = Memr[d[n3]+k] + } else + med = blank + } + } + + # Only set median and reorder if needed + n[i] = n1 + if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) { + j = max (nl, n1 + 1) + if (keepids) { + do l = 1, min (n1, nl-1) { + Memr[d[l]+k] = Memr[d[j]+k] + if (grow > 0) { + mp1 = m[l] + k + mp2 = m[j] + k + id = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = id + } else + Memi[m[l]+k] = Memi[m[j]+k] + j = j + 1 + } + } else { + do l = 1, min (n1, nl - 1) { + Memr[d[l]+k] = Memr[d[j]+k] + j = j + 1 + } + } + } + + if (combine == MEDIAN) + median[i] = med + } + + # Check if data flag needs to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + # Flag that the median is computed. + if (combine == MEDIAN) + docombine = false + else + docombine = true + + call sfree (sp) +end + diff --git a/noao/onedspec/scombine/generic/icaverage.x b/noao/onedspec/scombine/generic/icaverage.x new file mode 100644 index 00000000..6c5c870b --- /dev/null +++ b/noao/onedspec/scombine/generic/icaverage.x @@ -0,0 +1,84 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include +include "../icombine.h" + + +# IC_AVERAGE -- Compute the average image line. +# Options include a weight average. + +procedure ic_averager (d, m, n, wts, npts, average) + +pointer d[ARB] # Data pointers +pointer m[ARB] # Image ID pointers +int n[npts] # Number of points +real wts[ARB] # Weights +int npts # Number of output points per line +real average[npts] # Average (returned) + +int i, j, k +real sumwt, wt +real sum + +include "../icombine.com" + +begin + # If no data has been excluded do the average without checking the + # number of points and using the fact that the weights are normalized. + # If all the data has been excluded set the average to the blank value. + + if (dflag == D_ALL) { + if (dowts) { + do i = 1, npts { + k = i - 1 + wt = wts[Memi[m[1]+k]] + sum = Memr[d[1]+k] * wt + do j = 2, n[i] { + wt = wts[Memi[m[j]+k]] + sum = sum + Memr[d[j]+k] * wt + } + average[i] = sum + } + } else { + do i = 1, npts { + k = i - 1 + sum = Memr[d[1]+k] + do j = 2, n[i] + sum = sum + Memr[d[j]+k] + average[i] = sum / n[i] + } + } + } else if (dflag == D_NONE) { + do i = 1, npts + average[i] = blank + } else { + if (dowts) { + do i = 1, npts { + if (n[i] > 0) { + k = i - 1 + wt = wts[Memi[m[1]+k]] + sum = Memr[d[1]+k] * wt + sumwt = wt + do j = 2, n[i] { + wt = wts[Memi[m[j]+k]] + sum = sum + Memr[d[j]+k] * wt + sumwt = sumwt + wt + } + average[i] = sum / sumwt + } else + average[i] = blank + } + } else { + do i = 1, npts { + if (n[i] > 0) { + k = i - 1 + sum = Memr[d[1]+k] + do j = 2, n[i] + sum = sum + Memr[d[j]+k] + average[i] = sum / n[i] + } else + average[i] = blank + } + } + } +end diff --git a/noao/onedspec/scombine/generic/iccclip.x b/noao/onedspec/scombine/generic/iccclip.x new file mode 100644 index 00000000..26b17ba2 --- /dev/null +++ b/noao/onedspec/scombine/generic/iccclip.x @@ -0,0 +1,453 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + +define MINCLIP 2 # Mininum number of images for algorithm + + +# IC_ACCDCLIP -- Reject pixels using CCD noise parameters about the average + +procedure ic_accdclipr (d, m, n, scales, zeros, nm, nimages, npts, average) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +real nm[3,nimages] # Noise model parameters +int nimages # Number of images +int npts # Number of output points per line +real average[npts] # Average + +int i, j, k, l, jj, n1, n2, nin, nk, maxkeep +real d1, low, high, sum, a, s, r, zero +data zero /0.0/ +pointer sp, resid, dp1, dp2, mp1, mp2 + +include "../icombine.com" + +begin + # If there are no pixels go on to the combining. Since the unweighted + # average is computed here possibly skip the combining later. + + # There must be at least max (1, nkeep) pixels. + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } else if (dowts || combine != AVERAGE) + docombine = true + else + docombine = false + + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + + # There must be at least two pixels for rejection. The initial + # average is the low/high rejected average except in the case of + # just two pixels. The rejections are iterated and the average + # is recomputed. Corrections for scaling may be performed. + # Depending on other flags the image IDs may also need to be adjusted. + + nin = n[1] + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + if (n1 <= max (MINCLIP-1, maxkeep)) { + if (!docombine) { + if (n1 == 0) + average[i] = blank + else { + sum = Memr[d[1]+k] + do j = 2, n1 + sum = sum + Memr[d[j]+k] + average[i] = sum / n1 + } + } + next + } + + repeat { + if (n1 == 2) { + sum = Memr[d[1]+k] + sum = sum + Memr[d[2]+k] + a = sum / 2 + } else { + low = Memr[d[1]+k] + high = Memr[d[2]+k] + if (low > high) { + d1 = low + low = high + high = d1 + } + sum = 0. + do j = 3, n1 { + d1 = Memr[d[j]+k] + if (d1 < low) { + sum = sum + low + low = d1 + } else if (d1 > high) { + sum = sum + high + high = d1 + } else + sum = sum + d1 + } + a = sum / (n1 - 2) + sum = sum + low + high + } + n2 = n1 + if (doscale1) { + for (j=1; j<=n1; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + + l = Memi[mp1] + s = scales[l] + d1 = max (zero, s * (a + zeros[l])) + s = sqrt (nm[1,l] + d1/nm[2,l] + (d1*nm[3,l])**2) / s + + d1 = Memr[dp1] + r = (d1 - a) / s + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs(r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + mp2 = m[n1] + k + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } else { + if (!keepids) { + s = max (zero, a) + s = sqrt (nm[1,1] + s/nm[2,1] + (s*nm[3,1])**2) + } + for (j=1; j<=n1; j=j+1) { + if (keepids) { + l = Memi[m[j]+k] + s = max (zero, a) + s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2) + } + dp1 = d[j] + k + d1 = Memr[dp1] + r = (d1 - a) / s + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs(r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[n1] + k + l = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + } + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } + } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1)) + + if (n1 < maxkeep) { + nk = maxkeep + if (doscale1) { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } else { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } + } + + n[i] = n1 + if (!docombine) + if (n1 > 0) + average[i] = sum / n1 + else + average[i] = blank + } + + # Check if the data flag has to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + call sfree (sp) +end + + +# IC_CCDCLIP -- Reject pixels using CCD noise parameters about the median + +procedure ic_mccdclipr (d, m, n, scales, zeros, nm, nimages, npts, median) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +real nm[3,nimages] # Noise model +int nimages # Number of images +int npts # Number of output points per line +real median[npts] # Median + +int i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep +real r, s +pointer sp, resid, mp1, mp2 +real med, zero +data zero /0.0/ + +include "../icombine.com" + +begin + # There must be at least max (MINCLIP, nkeep+1) pixels. + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + + # Compute median and sigma and iteratively clip. + nin = n[1] + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + nl = 1 + nh = n1 + + repeat { + n2 = n1 + n3 = nl + n1 / 2 + + if (n1 == 0) + med = blank + else if (mod (n1, 2) == 0) { + med = Memr[d[n3-1]+k] + med = (med + Memr[d[n3]+k]) / 2. + } else + med = Memr[d[n3]+k] + + if (n1 >= max (MINCLIP, maxkeep+1)) { + if (doscale1) { + for (; nl <= n2; nl = nl + 1) { + l = Memi[m[nl]+k] + s = scales[l] + r = max (zero, s * (med + zeros[l])) + s = sqrt (nm[1,l] + r/nm[2,l] + (r*nm[3,l])**2) / s + r = (med - Memr[d[nl]+k]) / s + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + l = Memi[m[nh]+k] + s = scales[l] + r = max (zero, s * (med + zeros[l])) + s = sqrt (nm[1,l] + r/nm[2,l] + (r*nm[3,l])**2) / s + r = (Memr[d[nh]+k] - med) / s + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } else { + if (!keepids) { + s = max (zero, med) + s = sqrt (nm[1,1] + s/nm[2,1] + (s*nm[3,1])**2) + } + for (; nl <= n2; nl = nl + 1) { + if (keepids) { + l = Memi[m[nl]+k] + s = max (zero, med) + s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2) + } + r = (med - Memr[d[nl]+k]) / s + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + if (keepids) { + l = Memi[m[nh]+k] + s = max (zero, med) + s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2) + } + r = (Memr[d[nh]+k] - med) / s + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } + } + } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1)) + + while (n1 < maxkeep) { + if (nl == 1) + nh = nh + 1 + else if (nh == n[i]) + nl = nl - 1 + else { + r = Memr[resid+nl-1] + s = Memr[resid+nh+1] + if (r < s) { + nl = nl - 1 + r = r + TOL + if (s <= r) + nh = nh + 1 + if (nl > 1) { + if (Memr[resid+nl-1] <= r) + nl = nl - 1 + } + } else { + nh = nh + 1 + s = s + TOL + if (r <= s) + nl = nl - 1 + if (nh < n2) { + if (Memr[resid+nh+1] <= s) + nh = nh + 1 + } + } + } + n1 = nh - nl + 1 + } + + # Only set median and reorder if needed + n[i] = n1 + if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) { + j = max (nl, n1 + 1) + if (keepids) { + do l = 1, min (n1, nl-1) { + Memr[d[l]+k] = Memr[d[j]+k] + if (grow > 0) { + mp1 = m[l] + k + mp2 = m[j] + k + id = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = id + } else + Memi[m[l]+k] = Memi[m[j]+k] + j = j + 1 + } + } else { + do l = 1, min (n1, nl - 1) { + Memr[d[l]+k] = Memr[d[j]+k] + j = j + 1 + } + } + } + + if (combine == MEDIAN) + median[i] = med + } + + # Check if data flag needs to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + # Flag that the median is computed. + if (combine == MEDIAN) + docombine = false + else + docombine = true + + call sfree (sp) +end + diff --git a/noao/onedspec/scombine/generic/icgrow.x b/noao/onedspec/scombine/generic/icgrow.x new file mode 100644 index 00000000..074bd8c3 --- /dev/null +++ b/noao/onedspec/scombine/generic/icgrow.x @@ -0,0 +1,76 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + + +# IC_GROW -- Reject neigbors of rejected pixels. +# The rejected pixels are marked by having nonzero ids beyond the number +# of included pixels. The pixels rejected here are given zero ids +# to avoid growing of the pixels rejected here. The unweighted average +# can be updated but any rejected pixels requires the median to be +# recomputed. When the number of pixels at a grow point reaches nkeep +# no further pixels are rejected. Note that the rejection order is not +# based on the magnitude of the residuals and so a grow from a weakly +# rejected image pixel may take precedence over a grow from a strongly +# rejected image pixel. + +procedure ic_growr (d, m, n, nimages, npts, average) + +pointer d[ARB] # Data pointers +pointer m[ARB] # Image id pointers +int n[npts] # Number of good pixels +int nimages # Number of images +int npts # Number of output points per line +real average[npts] # Average + +int i1, i2, j1, j2, k1, k2, l, is, ie, n2, maxkeep +pointer mp1, mp2 + +include "../icombine.com" + +begin + if (dflag == D_NONE) + return + + do i1 = 1, npts { + k1 = i1 - 1 + is = max (1, i1 - grow) + ie = min (npts, i1 + grow) + do j1 = n[i1]+1, nimages { + l = Memi[m[j1]+k1] + if (l == 0) + next + if (combine == MEDIAN) + docombine = true + + do i2 = is, ie { + if (i2 == i1) + next + k2 = i2 - 1 + n2 = n[i2] + if (nkeep < 0) + maxkeep = max (0, n2 + nkeep) + else + maxkeep = min (n2, nkeep) + if (n2 <= maxkeep) + next + do j2 = 1, n2 { + mp1 = m[j2] + k2 + if (Memi[mp1] == l) { + if (!docombine && n2 > 1) + average[i2] = + (n2*average[i2] - Memr[d[j2]+k2]) / (n2-1) + mp2 = m[n2] + k2 + if (j2 < n2) { + Memr[d[j2]+k2] = Memr[d[n2]+k2] + Memi[mp1] = Memi[mp2] + } + Memi[mp2] = 0 + n[i2] = n2 - 1 + break + } + } + } + } + } +end diff --git a/noao/onedspec/scombine/generic/icmedian.x b/noao/onedspec/scombine/generic/icmedian.x new file mode 100644 index 00000000..e7607340 --- /dev/null +++ b/noao/onedspec/scombine/generic/icmedian.x @@ -0,0 +1,139 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + + +# IC_MEDIAN -- Median of lines + +procedure ic_medianr (d, n, npts, median) + +pointer d[ARB] # Input data line pointers +int n[npts] # Number of good pixels +int npts # Number of output points per line +real median[npts] # Median + +int i, j1, j2, j3, k, n1 +bool even +real val1, val2, val3 + +include "../icombine.com" + +begin + if (dflag == D_NONE) { + do i = 1, npts + median[i]= blank + return + } + + # Check for previous sorting + if (mclip) { + if (dflag == D_ALL) { + n1 = n[1] + even = (mod (n1, 2) == 0) + j1 = n1 / 2 + 1 + j2 = n1 / 2 + do i = 1, npts { + k = i - 1 + if (even) { + val1 = Memr[d[j1]+k] + val2 = Memr[d[j2]+k] + median[i] = (val1 + val2) / 2. + } else + median[i] = Memr[d[j1]+k] + } + } else { + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (n1 > 0) { + j1 = n1 / 2 + 1 + if (mod (n1, 2) == 0) { + j2 = n1 / 2 + val1 = Memr[d[j1]+k] + val2 = Memr[d[j2]+k] + median[i] = (val1 + val2) / 2. + } else + median[i] = Memr[d[j1]+k] + } else + median[i] = blank + } + } + return + } + + # Repeatedly exchange the extreme values until there are three + # or fewer pixels. + + do i = 1, npts { + k = i - 1 + n1 = n[i] + while (n1 > 3) { + j1 = 1 + j2 = 1 + val1 = Memr[d[j1]+k] + val2 = val1 + do j3 = 2, n1 { + val3 = Memr[d[j3]+k] + if (val3 > val1) { + j1 = j3 + val1 = val3 + } else if (val3 < val2) { + j2 = j3 + val2 = val3 + } + } + j3 = n1 - 1 + if (j1 < j3 && j2 < j3) { + Memr[d[j1]+k] = val3 + Memr[d[j2]+k] = Memr[d[j3]+k] + Memr[d[j3]+k] = val1 + Memr[d[n1]+k] = val2 + } else if (j1 < j3) { + if (j2 == j3) { + Memr[d[j1]+k] = val3 + Memr[d[n1]+k] = val1 + } else { + Memr[d[j1]+k] = Memr[d[j3]+k] + Memr[d[j3]+k] = val1 + } + } else if (j2 < j3) { + if (j1 == j3) { + Memr[d[j2]+k] = val3 + Memr[d[n1]+k] = val2 + } else { + Memr[d[j2]+k] = Memr[d[j3]+k] + Memr[d[j3]+k] = val2 + } + } + n1 = n1 - 2 + } + + if (n1 == 3) { + val1 = Memr[d[1]+k] + val2 = Memr[d[2]+k] + val3 = Memr[d[3]+k] + if (val1 < val2) { + if (val2 < val3) # abc + median[i] = val2 + else if (val1 < val3) # acb + median[i] = val3 + else # cab + median[i] = val1 + } else { + if (val2 > val3) # cba + median[i] = val2 + else if (val1 > val3) # bca + median[i] = val3 + else # bac + median[i] = val1 + } + } else if (n1 == 2) { + val1 = Memr[d[1]+k] + val2 = Memr[d[2]+k] + median[i] = (val1 + val2) / 2 + } else if (n1 == 1) + median[i] = Memr[d[1]+k] + else + median[i] = blank + } +end diff --git a/noao/onedspec/scombine/generic/icmm.x b/noao/onedspec/scombine/generic/icmm.x new file mode 100644 index 00000000..1c314241 --- /dev/null +++ b/noao/onedspec/scombine/generic/icmm.x @@ -0,0 +1,152 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + + +# IC_MM -- Reject a specified number of high and low pixels + +procedure ic_mmr (d, m, n, npts) + +pointer d[ARB] # Data pointers +pointer m[ARB] # Image ID pointers +int n[npts] # Number of good pixels +int npts # Number of output points per line + +int n1, ncombine, npairs, nlow, nhigh, np +int i, i1, j, jmax, jmin +pointer k, kmax, kmin +real d1, d2, dmin, dmax + +include "../icombine.com" + +begin + if (dflag == D_NONE) + return + + if (dflag == D_ALL) { + n1 = n[1] + nlow = flow * n1 + 0.001 + nhigh = fhigh * n1 + 0.001 + ncombine = n1 - nlow - nhigh + npairs = min (nlow, nhigh) + nlow = nlow - npairs + nhigh = nhigh - npairs + } + + do i = 1, npts { + i1 = i - 1 + n1 = n[i] + if (dflag == D_MIX) { + nlow = flow * n1 + 0.001 + nhigh = fhigh * n1 + 0.001 + ncombine = max (ncombine, n1 - nlow - nhigh) + npairs = min (nlow, nhigh) + nlow = nlow - npairs + nhigh = nhigh - npairs + } + + # Reject the npairs low and high points. + do np = 1, npairs { + k = d[1] + i1 + d1 = Memr[k] + dmax = d1; dmin = d1; jmax = 1; jmin = 1; kmax = k; kmin = k + do j = 2, n1 { + d2 = d1 + k = d[j] + i1 + d1 = Memr[k] + if (d1 > dmax) { + dmax = d1; jmax = j; kmax = k + } else if (d1 < dmin) { + dmin = d1; jmin = j; kmin = k + } + } + j = n1 - 1 + if (keepids) { + if (jmax < j) { + if (jmin != j) { + Memr[kmax] = d2 + Memi[m[jmax]+i1] = Memi[m[j]+i1] + } else { + Memr[kmax] = d1 + Memi[m[jmax]+i1] = Memi[m[n1]+i1] + } + } + if (jmin < j) { + if (jmax != n1) { + Memr[kmin] = d1 + Memi[m[jmin]+i1] = Memi[m[n1]+i1] + } else { + Memr[kmin] = d2 + Memi[m[jmin]+i1] = Memi[m[j]+i1] + } + } + } else { + if (jmax < j) { + if (jmin != j) + Memr[kmax] = d2 + else + Memr[kmax] = d1 + } + if (jmin < j) { + if (jmax != n1) + Memr[kmin] = d1 + else + Memr[kmin] = d2 + } + } + n1 = n1 - 2 + } + + # Reject the excess low points. + do np = 1, nlow { + k = d[1] + i1 + d1 = Memr[k] + dmin = d1; jmin = 1; kmin = k + do j = 2, n1 { + k = d[j] + i1 + d1 = Memr[k] + if (d1 < dmin) { + dmin = d1; jmin = j; kmin = k + } + } + if (keepids) { + if (jmin < n1) { + Memr[kmin] = d1 + Memi[m[jmin]+i1] = Memi[m[n1]+i1] + } + } else { + if (jmin < n1) + Memr[kmin] = d1 + } + n1 = n1 - 1 + } + + # Reject the excess high points. + do np = 1, nhigh { + k = d[1] + i1 + d1 = Memr[k] + dmax = d1; jmax = 1; kmax = k + do j = 2, n1 { + k = d[j] + i1 + d1 = Memr[k] + if (d1 > dmax) { + dmax = d1; jmax = j; kmax = k + } + } + if (keepids) { + if (jmax < n1) { + Memr[kmax] = d1 + Memi[m[jmax]+i1] = Memi[m[n1]+i1] + } + } else { + if (jmax < n1) + Memr[kmax] = d1 + } + n1 = n1 - 1 + } + n[i] = n1 + } + + if (dflag == D_ALL && npairs + nlow + nhigh > 0) + dflag = D_MIX +end diff --git a/noao/onedspec/scombine/generic/icpclip.x b/noao/onedspec/scombine/generic/icpclip.x new file mode 100644 index 00000000..d9028a93 --- /dev/null +++ b/noao/onedspec/scombine/generic/icpclip.x @@ -0,0 +1,224 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + +define MINCLIP 3 # Minimum number for clipping + + +# IC_PCLIP -- Percentile clip +# +# 1) Find the median +# 2) Find the pixel which is the specified order index away +# 3) Use the data value difference as a sigma and apply clipping +# 4) Since the median is known return it so it does not have to be recomputed + +procedure ic_pclipr (d, m, n, nimages, npts, median) + +pointer d[ARB] # Data pointers +pointer m[ARB] # Image id pointers +int n[npts] # Number of good pixels +int nimages # Number of input images +int npts # Number of output points per line +real median[npts] # Median + +int i, j, k, l, id, n1, n2, n3, n4, n5, nl, nh, nin, maxkeep +bool even, fp_equalr() +real sigma, r, s, t +pointer sp, resid, mp1, mp2 +real med + +include "../icombine.com" + +begin + # There must be at least MINCLIP and more than nkeep pixels. + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + + # Set sign of pclip parameter + if (pclip < 0) + t = -1. + else + t = 1. + + # If there are no rejected pixels compute certain parameters once. + if (dflag == D_ALL) { + n1 = n[1] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + n2 = 1 + n1 / 2 + even = (mod (n1, 2) == 0) + if (pclip < 0.) { + if (even) + n3 = max (1, nint (n2 - 1 + pclip)) + else + n3 = max (1, nint (n2 + pclip)) + } else + n3 = min (n1, nint (n2 + pclip)) + nin = n1 + } + + # Now apply clipping. + do i = 1, npts { + # Compute median. + if (dflag == D_MIX) { + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + if (n1 == 0) { + if (combine == MEDIAN) + median[i] = blank + next + } + n2 = 1 + n1 / 2 + even = (mod (n1, 2) == 0) + if (pclip < 0) { + if (even) + n3 = max (1, nint (n2 - 1 + pclip)) + else + n3 = max (1, nint (n2 + pclip)) + } else + n3 = min (n1, nint (n2 + pclip)) + } + + j = i - 1 + if (even) { + med = Memr[d[n2-1]+j] + med = (med + Memr[d[n2]+j]) / 2. + } else + med = Memr[d[n2]+j] + + if (n1 < max (MINCLIP, maxkeep+1)) { + if (combine == MEDIAN) + median[i] = med + next + } + + # Define sigma for clipping + sigma = t * (Memr[d[n3]+j] - med) + if (fp_equalr (sigma, 0.)) { + if (combine == MEDIAN) + median[i] = med + next + } + + # Reject pixels and save residuals. + # Check if any pixels are clipped. + # If so recompute the median and reset the number of good pixels. + # Only reorder if needed. + + for (nl=1; nl<=n1; nl=nl+1) { + r = (med - Memr[d[nl]+j]) / sigma + if (r < lsigma) + break + Memr[resid+nl] = r + } + for (nh=n1; nh>=1; nh=nh-1) { + r = (Memr[d[nh]+j] - med) / sigma + if (r < hsigma) + break + Memr[resid+nh] = r + } + n4 = nh - nl + 1 + + # If too many pixels are rejected add some back in. + # All pixels with the same residual are added. + while (n4 < maxkeep) { + if (nl == 1) + nh = nh + 1 + else if (nh == n[i]) + nl = nl - 1 + else { + r = Memr[resid+nl-1] + s = Memr[resid+nh+1] + if (r < s) { + nl = nl - 1 + r = r + TOL + if (s <= r) + nh = nh + 1 + if (nl > 1) { + if (Memr[resid+nl-1] <= r) + nl = nl - 1 + } + } else { + nh = nh + 1 + s = s + TOL + if (r <= s) + nl = nl - 1 + if (nh < n2) { + if (Memr[resid+nh+1] <= s) + nh = nh + 1 + } + } + } + n4 = nh - nl + 1 + } + + # If any pixels are rejected recompute the median. + if (nl > 1 || nh < n1) { + n5 = nl + n4 / 2 + if (mod (n4, 2) == 0) { + med = Memr[d[n5-1]+j] + med = (med + Memr[d[n5]+j]) / 2. + } else + med = Memr[d[n5]+j] + n[i] = n4 + } + if (combine == MEDIAN) + median[i] = med + + # Reorder if pixels only if necessary. + if (nl > 1 && (combine != MEDIAN || grow > 0)) { + k = max (nl, n4 + 1) + if (keepids) { + do l = 1, min (n1, nl-1) { + Memr[d[l]+j] = Memr[d[k]+j] + if (grow > 0) { + mp1 = m[l] + j + mp2 = m[k] + j + id = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = id + } else + Memi[m[l]+j] = Memi[m[k]+j] + k = k + 1 + } + } else { + do l = 1, min (n1, nl - 1) { + Memr[d[l]+j] = Memr[d[k]+j] + k = k + 1 + } + } + } + } + + # Check if data flag needs to be reset for rejected pixels. + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + # Flag whether the median has been computed. + if (combine == MEDIAN) + docombine = false + else + docombine = true + + call sfree (sp) +end diff --git a/noao/onedspec/scombine/generic/icsclip.x b/noao/onedspec/scombine/generic/icsclip.x new file mode 100644 index 00000000..e38f7935 --- /dev/null +++ b/noao/onedspec/scombine/generic/icsclip.x @@ -0,0 +1,486 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + +define MINCLIP 3 # Mininum number of images for algorithm + + +# IC_ASIGCLIP -- Reject pixels using sigma clipping about the average +# The initial average rejects the high and low pixels. A correction for +# different scalings of the images may be made. Weights are not used. + +procedure ic_asigclipr (d, m, n, scales, zeros, nimages, npts, average) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +int nimages # Number of images +int npts # Number of output points per line +real average[npts] # Average + +int i, j, k, l, jj, n1, n2, nin, nk, maxkeep +real d1, low, high, sum, a, s, r, one +data one /1.0/ +pointer sp, resid, w, wp, dp1, dp2, mp1, mp2 + +include "../icombine.com" + +begin + # If there are insufficient pixels go on to the combining + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + # Flag whether returned average needs to be recomputed. + if (dowts || combine != AVERAGE) + docombine = true + else + docombine = false + + # Save the residuals and the sigma scaling corrections if needed. + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + if (doscale1) + call salloc (w, nimages, TY_REAL) + + # Do sigma clipping. + nin = n[1] + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + + # If there are not enough pixels simply compute the average. + if (n1 < max (3, maxkeep)) { + if (!docombine) { + if (n1 == 0) + average[i] = blank + else { + sum = Memr[d[1]+k] + do j = 2, n1 + sum = sum + Memr[d[j]+k] + average[i] = sum / n1 + } + } + next + } + + # Compute average with the high and low rejected. + low = Memr[d[1]+k] + high = Memr[d[2]+k] + if (low > high) { + d1 = low + low = high + high = d1 + } + sum = 0. + do j = 3, n1 { + d1 = Memr[d[j]+k] + if (d1 < low) { + sum = sum + low + low = d1 + } else if (d1 > high) { + sum = sum + high + high = d1 + } else + sum = sum + d1 + } + a = sum / (n1 - 2) + sum = sum + low + high + + # Iteratively reject pixels and compute the final average if needed. + # Compact the data and keep track of the image IDs if needed. + + repeat { + n2 = n1 + if (doscale1) { + # Compute sigma corrected for scaling. + s = 0. + wp = w - 1 + do j = 1, n1 { + dp1 = d[j] + k + mp1 = m[j] + k + wp = wp + 1 + + d1 = Memr[dp1] + l = Memi[mp1] + r = sqrt (max (one, (a + zeros[l]) / scales[l])) + s = s + ((d1 - a) / r) ** 2 + Memr[wp] = r + } + s = sqrt (s / (n1 - 1)) + + # Reject pixels. Save the residuals and data values. + wp = w - 1 + if (s > 0.) { + for (j=1; j<=n1; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + wp = wp + 1 + + d1 = Memr[dp1] + r = (d1 - a) / (s * Memr[wp]) + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs (r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + Memr[wp] = Memr[w+n1-1] + mp2 = m[n1] + k + l = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } + } else { + # Compute the sigma without scale correction. + s = 0. + do j = 1, n1 + s = s + (Memr[d[j]+k] - a) ** 2 + s = sqrt (s / (n1 - 1)) + + # Reject pixels. Save the residuals and data values. + if (s > 0.) { + for (j=1; j<=n1; j=j+1) { + dp1 = d[j] + k + d1 = Memr[dp1] + r = (d1 - a) / s + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs (r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[n1] + k + l = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + } + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } + } + + # Recompute the average. + if (n1 > 1) + a = sum / n1 + } until (n1 == n2 || n1 <= max (2, maxkeep)) + + # If too many pixels are rejected add some back. + # All pixels with equal residuals are added back. + if (n1 < maxkeep) { + nk = maxkeep + if (doscale1) { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } else { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } + + # Recompute the average. + if (n1 > 1) + a = sum / n1 + } + + # Save the average if needed. + n[i] = n1 + if (!docombine) { + if (n1 > 0) + average[i] = a + else + average[i] = blank + } + } + + # Check if the data flag has to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + call sfree (sp) +end + + +# IC_MSIGCLIP -- Reject pixels using sigma clipping about the median + +procedure ic_msigclipr (d, m, n, scales, zeros, nimages, npts, median) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +int nimages # Number of images +int npts # Number of output points per line +real median[npts] # Median + +int i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep +real r, s +pointer sp, resid, w, mp1, mp2 +real med, one +data one /1.0/ + +include "../icombine.com" + +begin + # If there are insufficient pixels go on to the combining + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + # Save the residuals and sigma scaling corrections if needed. + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + if (doscale1) + call salloc (w, nimages, TY_REAL) + + # Compute median and sigma and iteratively clip. + nin = n[1] + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + nl = 1 + nh = n1 + + repeat { + n2 = n1 + n3 = nl + n1 / 2 + + if (n1 == 0) + med = blank + else if (mod (n1, 2) == 0) + med = (Memr[d[n3-1]+k] + Memr[d[n3]+k]) / 2. + else + med = Memr[d[n3]+k] + + if (n1 >= max (MINCLIP, maxkeep+1)) { + if (doscale1) { + # Compute the sigma with scaling correction. + s = 0. + do j = nl, nh { + l = Memi[m[j]+k] + r = sqrt (max (one, (med + zeros[l]) / scales[l])) + s = s + ((Memr[d[j]+k] - med) / r) ** 2 + Memr[w+j-1] = r + } + s = sqrt (s / (n1 - 1)) + + # Reject pixels and save the residuals. + if (s > 0.) { + for (; nl <= n2; nl = nl + 1) { + r = (med - Memr[d[nl]+k]) / (s * Memr[w+nl-1]) + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + r = (Memr[d[nh]+k] - med) / (s * Memr[w+nh-1]) + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } + } else { + # Compute the sigma without scaling correction. + s = 0. + do j = nl, nh + s = s + (Memr[d[j]+k] - med) ** 2 + s = sqrt (s / (n1 - 1)) + + # Reject pixels and save the residuals. + if (s > 0.) { + for (; nl <= n2; nl = nl + 1) { + r = (med - Memr[d[nl]+k]) / s + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + r = (Memr[d[nh]+k] - med) / s + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } + } + } + } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1)) + + # If too many pixels are rejected add some back. + # All pixels with equal residuals are added back. + while (n1 < maxkeep) { + if (nl == 1) + nh = nh + 1 + else if (nh == n[i]) + nl = nl - 1 + else { + r = Memr[resid+nl-1] + s = Memr[resid+nh+1] + if (r < s) { + nl = nl - 1 + r = r + TOL + if (s <= r) + nh = nh + 1 + if (nl > 1) { + if (Memr[resid+nl-1] <= r) + nl = nl - 1 + } + } else { + nh = nh + 1 + s = s + TOL + if (r <= s) + nl = nl - 1 + if (nh < n2) { + if (Memr[resid+nh+1] <= s) + nh = nh + 1 + } + } + } + n1 = nh - nl + 1 + } + + # Only set median and reorder if needed + n[i] = n1 + if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) { + j = max (nl, n1 + 1) + if (keepids) { + do l = 1, min (n1, nl-1) { + Memr[d[l]+k] = Memr[d[j]+k] + if (grow > 0) { + mp1 = m[l] + k + mp2 = m[j] + k + id = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = id + } else + Memi[m[l]+k] = Memi[m[j]+k] + j = j + 1 + } + } else { + do l = 1, min (n1, nl - 1) { + Memr[d[l]+k] = Memr[d[j]+k] + j = j + 1 + } + } + } + + if (combine == MEDIAN) + median[i] = med + } + + # Check if data flag needs to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + # Flag that the median has been computed. + if (combine == MEDIAN) + docombine = false + else + docombine = true + + call sfree (sp) +end + diff --git a/noao/onedspec/scombine/generic/icsort.x b/noao/onedspec/scombine/generic/icsort.x new file mode 100644 index 00000000..f3d2fb21 --- /dev/null +++ b/noao/onedspec/scombine/generic/icsort.x @@ -0,0 +1,275 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +define LOGPTR 32 # log2(maxpts) (4e9) + + +# IC_SORT -- Quicksort. This is based on the VOPS asrt except that +# the input is an array of pointers to image lines and the sort is done +# across the image lines at each point along the lines. The number of +# valid pixels at each point is allowed to vary. The cases of 1, 2, and 3 +# pixels per point are treated specially. + +procedure ic_sortr (a, b, nvecs, npts) + +pointer a[ARB] # pointer to input vectors +real b[ARB] # work array +int nvecs[npts] # number of vectors +int npts # number of points in vectors + +real pivot, temp, temp3 +int i, j, k, l, p, npix, lv[LOGPTR], uv[LOGPTR] +define swap {temp=$1;$1=$2;$2=temp} +define copy_ 10 + +begin + do l = 0, npts-1 { + npix = nvecs[l+1] + if (npix <= 1) + next + + do i = 1, npix + b[i] = Memr[a[i]+l] + + # Special cases + if (npix <= 3) { + pivot = b[1] + temp = b[2] + if (npix == 2) { + if (temp < pivot) { + b[1] = temp + b[2] = pivot + } else + next + } else { + temp3 = b[3] + if (temp < pivot) { # bac|bca|cba + if (temp < temp3) { # bac|bca + b[1] = temp + if (pivot < temp3) # bac + b[2] = pivot + else { # bca + b[2] = temp3 + b[3] = pivot + } + } else { # cba + b[1] = temp3 + b[3] = pivot + } + } else if (temp3 < temp) { # acb|cab + b[3] = temp + if (pivot < temp3) # acb + b[2] = temp3 + else { # cab + b[1] = temp3 + b[2] = pivot + } + } else + next + } + goto copy_ + } + + # General case + do i = 1, npix + b[i] = Memr[a[i]+l] + + lv[1] = 1 + uv[1] = npix + p = 1 + + while (p > 0) { + if (lv[p] >= uv[p]) # only one elem in this subset + p = p - 1 # pop stack + else { + # Dummy do loop to trigger the Fortran optimizer. + do p = p, ARB { + i = lv[p] - 1 + j = uv[p] + + # Select as the pivot the element at the center of the + # array, to avoid quadratic behavior on an already + # sorted array. + + k = (lv[p] + uv[p]) / 2 + swap (b[j], b[k]) + pivot = b[j] # pivot line + + while (i < j) { + for (i=i+1; b[i] < pivot; i=i+1) + ; + for (j=j-1; j > i; j=j-1) + if (b[j] <= pivot) + break + if (i < j) # out of order pair + swap (b[i], b[j]) # interchange elements + } + + j = uv[p] # move pivot to position i + swap (b[i], b[j]) # interchange elements + + if (i-lv[p] < uv[p] - i) { # stack so shorter done first + lv[p+1] = lv[p] + uv[p+1] = i - 1 + lv[p] = i + 1 + } else { + lv[p+1] = i + 1 + uv[p+1] = uv[p] + uv[p] = i - 1 + } + + break + } + p = p + 1 # push onto stack + } + } + +copy_ + do i = 1, npix + Memr[a[i]+l] = b[i] + } +end + + +# IC_2SORT -- Quicksort. This is based on the VOPS asrt except that +# the input is an array of pointers to image lines and the sort is done +# across the image lines at each point along the lines. The number of +# valid pixels at each point is allowed to vary. The cases of 1, 2, and 3 +# pixels per point are treated specially. A second integer set of +# vectors is sorted. + +procedure ic_2sortr (a, b, c, d, nvecs, npts) + +pointer a[ARB] # pointer to input vectors +real b[ARB] # work array +pointer c[ARB] # pointer to associated integer vectors +int d[ARB] # work array +int nvecs[npts] # number of vectors +int npts # number of points in vectors + +real pivot, temp, temp3 +int i, j, k, l, p, npix, lv[LOGPTR], uv[LOGPTR], itemp +define swap {temp=$1;$1=$2;$2=temp} +define iswap {itemp=$1;$1=$2;$2=itemp} +define copy_ 10 + +begin + do l = 0, npts-1 { + npix = nvecs[l+1] + if (npix <= 1) + next + + do i = 1, npix { + b[i] = Memr[a[i]+l] + d[i] = Memi[c[i]+l] + } + + # Special cases + if (npix <= 3) { + pivot = b[1] + temp = b[2] + if (npix == 2) { + if (temp < pivot) { + b[1] = temp + b[2] = pivot + iswap (d[1], d[2]) + } else + next + } else { + temp3 = b[3] + if (temp < pivot) { # bac|bca|cba + if (temp < temp3) { # bac|bca + b[1] = temp + if (pivot < temp3) { # bac + b[2] = pivot + iswap (d[1], d[2]) + } else { # bca + b[2] = temp3 + b[3] = pivot + itemp = d[2] + d[2] = d[3] + } + } else { # cba + b[1] = temp3 + b[3] = pivot + iswap (d[1], d[3]) + } + } else if (temp3 < temp) { # acb|cab + b[3] = temp + if (pivot < temp3) { # acb + b[2] = temp3 + iswap (d[2], d[3]) + } else { # cab + b[1] = temp3 + b[2] = pivot + itemp = d[2] + d[2] = d[1] + d[1] = d[3] + d[3] = itemp + } + } else + next + } + goto copy_ + } + + # General case + lv[1] = 1 + uv[1] = npix + p = 1 + + while (p > 0) { + if (lv[p] >= uv[p]) # only one elem in this subset + p = p - 1 # pop stack + else { + # Dummy do loop to trigger the Fortran optimizer. + do p = p, ARB { + i = lv[p] - 1 + j = uv[p] + + # Select as the pivot the element at the center of the + # array, to avoid quadratic behavior on an already + # sorted array. + + k = (lv[p] + uv[p]) / 2 + swap (b[j], b[k]); swap (d[j], d[k]) + pivot = b[j] # pivot line + + while (i < j) { + for (i=i+1; b[i] < pivot; i=i+1) + ; + for (j=j-1; j > i; j=j-1) + if (b[j] <= pivot) + break + if (i < j) { # out of order pair + swap (b[i], b[j]) # interchange elements + swap (d[i], d[j]) + } + } + + j = uv[p] # move pivot to position i + swap (b[i], b[j]) # interchange elements + swap (d[i], d[j]) + + if (i-lv[p] < uv[p] - i) { # stack so shorter done first + lv[p+1] = lv[p] + uv[p+1] = i - 1 + lv[p] = i + 1 + } else { + lv[p+1] = i + 1 + uv[p+1] = uv[p] + uv[p] = i - 1 + } + + break + } + p = p + 1 # push onto stack + } + } + +copy_ + do i = 1, npix { + Memr[a[i]+l] = b[i] + Memi[c[i]+l] = d[i] + } + } +end diff --git a/noao/onedspec/scombine/generic/mkpkg b/noao/onedspec/scombine/generic/mkpkg new file mode 100644 index 00000000..4d371363 --- /dev/null +++ b/noao/onedspec/scombine/generic/mkpkg @@ -0,0 +1,16 @@ +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + icaclip.x ../icombine.com ../icombine.h + icaverage.x ../icombine.com ../icombine.h + iccclip.x ../icombine.com ../icombine.h + icgrow.x ../icombine.com ../icombine.h + icmedian.x ../icombine.com ../icombine.h + icmm.x ../icombine.com ../icombine.h + icpclip.x ../icombine.com ../icombine.h + icsclip.x ../icombine.com ../icombine.h + icsort.x + ; -- cgit