aboutsummaryrefslogtreecommitdiff
path: root/noao/imred/ccdred/src/generic
diff options
context:
space:
mode:
Diffstat (limited to 'noao/imred/ccdred/src/generic')
-rw-r--r--noao/imred/ccdred/src/generic/ccdred.h150
-rw-r--r--noao/imred/ccdred/src/generic/cor.x694
-rw-r--r--noao/imred/ccdred/src/generic/icaclip.x1102
-rw-r--r--noao/imred/ccdred/src/generic/icaverage.x163
-rw-r--r--noao/imred/ccdred/src/generic/iccclip.x898
-rw-r--r--noao/imred/ccdred/src/generic/icgdata.x459
-rw-r--r--noao/imred/ccdred/src/generic/icgrow.x148
-rw-r--r--noao/imred/ccdred/src/generic/icmedian.x343
-rw-r--r--noao/imred/ccdred/src/generic/icmm.x300
-rw-r--r--noao/imred/ccdred/src/generic/icombine.x607
-rw-r--r--noao/imred/ccdred/src/generic/icpclip.x442
-rw-r--r--noao/imred/ccdred/src/generic/icsclip.x964
-rw-r--r--noao/imred/ccdred/src/generic/icsigma.x205
-rw-r--r--noao/imred/ccdred/src/generic/icsort.x550
-rw-r--r--noao/imred/ccdred/src/generic/icstat.x444
-rw-r--r--noao/imred/ccdred/src/generic/mkpkg11
-rw-r--r--noao/imred/ccdred/src/generic/proc.x735
17 files changed, 8215 insertions, 0 deletions
diff --git a/noao/imred/ccdred/src/generic/ccdred.h b/noao/imred/ccdred/src/generic/ccdred.h
new file mode 100644
index 00000000..2d370d86
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/ccdred.h
@@ -0,0 +1,150 @@
+# CCDRED Data Structures and Definitions
+
+# The CCD structure: This structure is used to communicate processing
+# parameters between the package procedures. It contains pointers to
+# data, calibration image IMIO pointers, scaling parameters, and the
+# correction flags. The corrections flags indicate which processing
+# operations are to be performed. The subsection parameters do not
+# include a step size. A step size is assumed. If arbitrary subsampling
+# is desired this would be the next generalization.
+
+define LEN_CCD 131 # Length of CCD structure
+
+# CCD data coordinates
+define CCD_C1 Memi[$1] # CCD starting column
+define CCD_C2 Memi[$1+1] # CCD ending column
+define CCD_L1 Memi[$1+2] # CCD starting line
+define CCD_L2 Memi[$1+3] # CCD ending line
+
+# Input data
+define IN_IM Memi[$1+10] # Input image pointer
+define IN_C1 Memi[$1+11] # Input data starting column
+define IN_C2 Memi[$1+12] # Input data ending column
+define IN_L1 Memi[$1+13] # Input data starting line
+define IN_L2 Memi[$1+14] # Input data ending line
+
+# Output data
+define OUT_IM Memi[$1+20] # Output image pointer
+define OUT_C1 Memi[$1+21] # Output data starting column
+define OUT_C2 Memi[$1+22] # Output data ending column
+define OUT_L1 Memi[$1+23] # Output data starting line
+define OUT_L2 Memi[$1+24] # Output data ending line
+
+# Mask data
+define MASK_IM Memi[$1+30] # Mask image pointer
+define MASK_C1 Memi[$1+31] # Mask data starting column
+define MASK_C2 Memi[$1+32] # Mask data ending column
+define MASK_L1 Memi[$1+33] # Mask data starting line
+define MASK_L2 Memi[$1+34] # Mask data ending line
+define MASK_PM Memi[$1+35] # Mask pointer
+define MASK_FP Memi[$1+36] # Mask fixpix data
+
+# Zero level data
+define ZERO_IM Memi[$1+40] # Zero level image pointer
+define ZERO_C1 Memi[$1+41] # Zero level data starting column
+define ZERO_C2 Memi[$1+42] # Zero level data ending column
+define ZERO_L1 Memi[$1+43] # Zero level data starting line
+define ZERO_L2 Memi[$1+44] # Zero level data ending line
+
+# Dark count data
+define DARK_IM Memi[$1+50] # Dark count image pointer
+define DARK_C1 Memi[$1+51] # Dark count data starting column
+define DARK_C2 Memi[$1+52] # Dark count data ending column
+define DARK_L1 Memi[$1+53] # Dark count data starting line
+define DARK_L2 Memi[$1+54] # Dark count data ending line
+
+# Flat field data
+define FLAT_IM Memi[$1+60] # Flat field image pointer
+define FLAT_C1 Memi[$1+61] # Flat field data starting column
+define FLAT_C2 Memi[$1+62] # Flat field data ending column
+define FLAT_L1 Memi[$1+63] # Flat field data starting line
+define FLAT_L2 Memi[$1+64] # Flat field data ending line
+
+# Illumination data
+define ILLUM_IM Memi[$1+70] # Illumination image pointer
+define ILLUM_C1 Memi[$1+71] # Illumination data starting column
+define ILLUM_C2 Memi[$1+72] # Illumination data ending column
+define ILLUM_L1 Memi[$1+73] # Illumination data starting line
+define ILLUM_L2 Memi[$1+74] # Illumination data ending line
+
+# Fringe data
+define FRINGE_IM Memi[$1+80] # Fringe image pointer
+define FRINGE_C1 Memi[$1+81] # Fringe data starting column
+define FRINGE_C2 Memi[$1+82] # Fringe data ending column
+define FRINGE_L1 Memi[$1+83] # Fringe data starting line
+define FRINGE_L2 Memi[$1+84] # Fringe data ending line
+
+# Trim section
+define TRIM_C1 Memi[$1+90] # Trim starting column
+define TRIM_C2 Memi[$1+91] # Trim ending column
+define TRIM_L1 Memi[$1+92] # Trim starting line
+define TRIM_L2 Memi[$1+93] # Trim ending line
+
+# Bias section
+define BIAS_C1 Memi[$1+100] # Bias starting column
+define BIAS_C2 Memi[$1+101] # Bias ending column
+define BIAS_L1 Memi[$1+102] # Bias starting line
+define BIAS_L2 Memi[$1+103] # Bias ending line
+
+define READAXIS Memi[$1+110] # Read out axis (1=cols, 2=lines)
+define CALCTYPE Memi[$1+111] # Calculation data type
+define OVERSCAN_TYPE Memi[$1+112] # Overscan type
+define OVERSCAN_VEC Memi[$1+113] # Pointer to overscan vector
+define DARKSCALE Memr[P2R($1+114)] # Dark count scale factor
+define FRINGESCALE Memr[P2R($1+115)] # Fringe scale factor
+define FLATSCALE Memr[P2R($1+116)] # Flat field scale factor
+define ILLUMSCALE Memr[P2R($1+117)] # Illumination scale factor
+define MINREPLACE Memr[P2R($1+118)] # Minimum replacement value
+define MEAN Memr[P2R($1+119)] # Mean of output image
+define COR Memi[$1+120] # Overall correction flag
+define CORS Memi[$1+121+($2-1)] # Individual correction flags
+
+# The correction array contains the following elements with array indices
+# given by the macro definitions.
+
+define NCORS 10 # Number of corrections
+
+define FIXPIX 1 # Fix bad pixels
+define TRIM 2 # Trim image
+define OVERSCAN 3 # Apply overscan correction
+define ZEROCOR 4 # Apply zero level correction
+define DARKCOR 5 # Apply dark count correction
+define FLATCOR 6 # Apply flat field correction
+define ILLUMCOR 7 # Apply illumination correction
+define FRINGECOR 8 # Apply fringe correction
+define FINDMEAN 9 # Find the mean of the output image
+define MINREP 10 # Check and replace minimum value
+
+# The following definitions identify the correction values in the correction
+# array. They are defined in terms of bit fields so that it is possible to
+# add corrections to form unique combination corrections. Some of
+# these combinations are implemented as compound operations for efficiency.
+
+define O 001B # overscan
+define Z 002B # zero level
+define D 004B # dark count
+define F 010B # flat field
+define I 020B # Illumination
+define Q 040B # Fringe
+
+# The following correction combinations are recognized.
+
+define ZO 003B # zero level + overscan
+define DO 005B # dark count + overscan
+define DZ 006B # dark count + zero level
+define DZO 007B # dark count + zero level + overscan
+define FO 011B # flat field + overscan
+define FZ 012B # flat field + zero level
+define FZO 013B # flat field + zero level + overscan
+define FD 014B # flat field + dark count
+define FDO 015B # flat field + dark count + overscan
+define FDZ 016B # flat field + dark count + zero level
+define FDZO 017B # flat field + dark count + zero level + overscan
+define QI 060B # fringe + illumination
+
+# The following overscan functions are recognized.
+define OVERSCAN_TYPES "|mean|median|minmax|chebyshev|legendre|spline3|spline1|"
+define OVERSCAN_MEAN 1 # Mean of overscan
+define OVERSCAN_MEDIAN 2 # Median of overscan
+define OVERSCAN_MINMAX 3 # Minmax of overscan
+define OVERSCAN_FIT 4 # Following codes are function fits
diff --git a/noao/imred/ccdred/src/generic/cor.x b/noao/imred/ccdred/src/generic/cor.x
new file mode 100644
index 00000000..fd2a8d6b
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/cor.x
@@ -0,0 +1,694 @@
+include "ccdred.h"
+
+
+.help cor Feb87 noao.imred.ccdred
+.nf ----------------------------------------------------------------------------
+cor -- Process CCD image lines
+
+These procedures are the heart of the CCD processing. They do the desired
+set of processing operations on the image line data as efficiently as
+possible. They are called by the PROC procedures. There are four procedures
+one for each readout axis and one for short and real image data.
+Some sets of operations are coded as single compound operations for efficiency.
+To keep the number of combinations managable only the most common
+combinations are coded as compound operations. The combinations
+consist of any set of line overscan, column overscan, zero level, dark
+count, and flat field and any set of illumination and fringe
+correction. The corrections are applied in place to the output vector.
+
+The column readout procedure is more complicated in order to handle
+zero level and flat field corrections specified as one dimensional
+readout corrections instead of two dimensional calibration images.
+Column readout format is probably extremely rare and the 1D readout
+corrections are used only for special types of data.
+.ih
+SEE ALSO
+proc, ccdred.h
+.endhelp -----------------------------------------------------------------------
+
+
+# COR1 -- Correct image lines with readout axis 1 (lines).
+
+procedure cor1s (cors, out, overscan, zero, dark, flat, illum,
+ fringe, n, darkscale, flatscale, illumscale, frgscale)
+
+int cors[ARB] # Correction flags
+short out[n] # Output data
+real overscan # Overscan value
+short zero[n] # Zero level correction
+short dark[n] # Dark count correction
+short flat[n] # Flat field correction
+short illum[n] # Illumination correction
+short fringe[n] # Fringe correction
+int n # Number of pixels
+real darkscale # Dark count scale factor
+real flatscale # Flat field scale factor
+real illumscale # Illumination scale factor
+real frgscale # Fringe scale factor
+
+int i, op
+
+begin
+ op = cors[OVERSCAN] + cors[ZEROCOR] + cors[DARKCOR] + cors[FLATCOR]
+ switch (op) {
+ case O: # overscan
+ do i = 1, n
+ out[i] = out[i] - overscan
+ case Z: # zero level
+ do i = 1, n
+ out[i] = out[i] - zero[i]
+
+ case ZO: # zero level + overscan
+ do i = 1, n
+ out[i] = out[i] - overscan - zero[i]
+
+ case D: # dark count
+ do i = 1, n
+ out[i] = out[i] - darkscale * dark[i]
+ case DO: # dark count + overscan
+ do i = 1, n
+ out[i] = out[i] - overscan - darkscale * dark[i]
+ case DZ: # dark count + zero level
+ do i = 1, n
+ out[i] = out[i] - zero[i] - darkscale * dark[i]
+ case DZO: # dark count + zero level + overscan
+ do i = 1, n
+ out[i] = out[i] - overscan - zero[i] - darkscale * dark[i]
+
+ case F: # flat field
+ do i = 1, n
+ out[i] = out[i] * flatscale / flat[i]
+ case FO: # flat field + overscan
+ do i = 1, n
+ out[i] = (out[i] - overscan) * flatscale / flat[i]
+ case FZ: # flat field + zero level
+ do i = 1, n
+ out[i] = (out[i] - zero[i]) * flatscale / flat[i]
+ case FZO: # flat field + zero level + overscan
+ do i = 1, n
+ out[i] = (out[i] - overscan - zero[i]) * flatscale /
+ flat[i]
+ case FD: # flat field + dark count
+ do i = 1, n
+ out[i] = (out[i] - darkscale * dark[i]) * flatscale / flat[i]
+ case FDO: # flat field + dark count + overscan
+ do i = 1, n
+ out[i] = (out[i] - overscan - darkscale * dark[i]) *
+ flatscale / flat[i]
+ case FDZ: # flat field + dark count + zero level
+ do i = 1, n
+ out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+ flatscale / flat[i]
+ case FDZO: # flat field + dark count + zero level + overscan
+ do i = 1, n
+ out[i] = (out[i] - overscan - zero[i] -
+ darkscale * dark[i]) * flatscale / flat[i]
+ }
+
+ # Often these operations will not be performed so test for no
+ # correction rather than go through the switch.
+
+ op = cors[ILLUMCOR] + cors[FRINGECOR]
+ if (op != 0) {
+ switch (op) {
+ case I: # illumination
+ do i = 1, n
+ out[i] = out[i] * illumscale / illum[i]
+ case Q: # fringe
+ do i = 1, n
+ out[i] = out[i] - frgscale * fringe[i]
+ case QI: # fringe + illumination
+ do i = 1, n
+ out[i] = out[i]*illumscale/illum[i] - frgscale*fringe[i]
+ }
+ }
+end
+
+
+# COR2 -- Correct lines for readout axis 2 (columns). This procedure is
+# more complex than when the readout is along the image lines because the
+# zero level and/or flat field corrections may be single readout column
+# vectors.
+
+procedure cor2s (line, cors, out, overscan, zero, dark, flat, illum,
+ fringe, n, zeroim, flatim, darkscale, flatscale, illumscale, frgscale)
+
+int line # Line to be corrected
+int cors[ARB] # Correction flags
+short out[n] # Output data
+real overscan[n] # Overscan value
+short zero[n] # Zero level correction
+short dark[n] # Dark count correction
+short flat[n] # Flat field correction
+short illum[n] # Illumination correction
+short fringe[n] # Fringe correction
+int n # Number of pixels
+pointer zeroim # Zero level IMIO pointer (NULL if 1D vector)
+pointer flatim # Flat field IMIO pointer (NULL if 1D vector)
+real darkscale # Dark count scale factor
+real flatscale # Flat field scale factor
+real illumscale # Illumination scale factor
+real frgscale # Fringe scale factor
+
+short zeroval
+real flatval
+int i, op
+
+begin
+ op = cors[OVERSCAN] + cors[ZEROCOR] + cors[DARKCOR] + cors[FLATCOR]
+ switch (op) {
+ case O: # overscan
+ do i = 1, n
+ out[i] = out[i] - overscan[i]
+ case Z: # zero level
+ if (zeroim != NULL)
+ do i = 1, n
+ out[i] = out[i] - zero[i]
+ else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = out[i] - zeroval
+ }
+
+ case ZO: # zero level + overscan
+ if (zeroim != NULL)
+ do i = 1, n
+ out[i] = out[i] - overscan[i] - zero[i]
+ else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = out[i] - overscan[i] - zeroval
+ }
+
+ case D: # dark count
+ do i = 1, n
+ out[i] = out[i] - darkscale * dark[i]
+ case DO: # dark count + overscan
+ do i = 1, n
+ out[i] = out[i] - overscan[i] - darkscale * dark[i]
+ case DZ: # dark count + zero level
+ if (zeroim != NULL)
+ do i = 1, n
+ out[i] = out[i] - zero[i] - darkscale * dark[i]
+ else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = out[i] - zeroval - darkscale * dark[i]
+ }
+ case DZO: # dark count + zero level + overscan
+ if (zeroim != NULL)
+ do i = 1, n
+ out[i] = out[i] - overscan[i] - zero[i] -
+ darkscale * dark[i]
+ else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = out[i] - overscan[i] - zeroval -
+ darkscale * dark[i]
+ }
+
+ case F: # flat field
+ if (flatim != NULL) {
+ do i = 1, n
+ out[i] = out[i] * flatscale / flat[i]
+ } else {
+ flatval = flatscale / flat[line]
+ do i = 1, n
+ out[i] = out[i] * flatval
+ }
+ case FO: # flat field + overscan
+ if (flatim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i]) * flatscale / flat[i]
+ } else {
+ flatval = flatscale / flat[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i]) * flatval
+ }
+ case FZ: # flat field + zero level
+ if (flatim != NULL) {
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - zero[i]) * flatscale / flat[i]
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - zeroval) * flatscale / flat[i]
+ }
+ } else {
+ flatval = flatscale / flat[line]
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - zero[i]) * flatval
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - zeroval) * flatval
+ }
+ }
+ case FZO: # flat field + zero level + overscan
+ if (flatim != NULL) {
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zero[i]) *
+ flatscale / flat[i]
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zeroval) *
+ flatscale / flat[i]
+ }
+ } else {
+ flatval = flatscale / flat[line]
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zero[i]) * flatval
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zeroval) * flatval
+ }
+ }
+ case FD: # flat field + dark count
+ if (flatim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - darkscale * dark[i]) * flatscale/flat[i]
+ } else {
+ flatval = flatscale / flat[line]
+ do i = 1, n
+ out[i] = (out[i] - darkscale * dark[i]) * flatval
+ }
+ case FDO: # flat field + dark count + overscan
+ if (flatim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - darkscale * dark[i]) *
+ flatscale / flat[i]
+ } else {
+ flatval = flatscale / flat[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - darkscale * dark[i]) *
+ flatval
+ }
+ case FDZ: # flat field + dark count + zero level
+ if (flatim != NULL) {
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+ flatscale / flat[i]
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - zeroval - darkscale * dark[i]) *
+ flatscale / flat[i]
+ }
+ } else {
+ flatval = flatscale / flat[line]
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+ flatval
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - zeroval - darkscale * dark[i]) *
+ flatval
+ }
+ }
+ case FDZO: # flat field + dark count + zero level + overscan
+ if (flatim != NULL) {
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zero[i] -
+ darkscale * dark[i]) * flatscale / flat[i]
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zeroval -
+ darkscale * dark[i]) * flatscale / flat[i]
+ }
+ } else {
+ flatval = flatscale / flat[line]
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zero[i] -
+ darkscale * dark[i]) * flatval
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zeroval -
+ darkscale * dark[i]) * flatval
+ }
+ }
+ }
+
+ # Often these operations will not be performed so test for no
+ # correction rather than go through the switch.
+
+ op = cors[ILLUMCOR] + cors[FRINGECOR]
+ if (op != 0) {
+ switch (op) {
+ case I: # illumination
+ do i = 1, n
+ out[i] = out[i] * illumscale / illum[i]
+ case Q: # fringe
+ do i = 1, n
+ out[i] = out[i] - frgscale * fringe[i]
+ case QI: # fringe + illumination
+ do i = 1, n
+ out[i] = out[i]*illumscale/illum[i] - frgscale*fringe[i]
+ }
+ }
+end
+
+# COR1 -- Correct image lines with readout axis 1 (lines).
+
+procedure cor1r (cors, out, overscan, zero, dark, flat, illum,
+ fringe, n, darkscale, flatscale, illumscale, frgscale)
+
+int cors[ARB] # Correction flags
+real out[n] # Output data
+real overscan # Overscan value
+real zero[n] # Zero level correction
+real dark[n] # Dark count correction
+real flat[n] # Flat field correction
+real illum[n] # Illumination correction
+real fringe[n] # Fringe correction
+int n # Number of pixels
+real darkscale # Dark count scale factor
+real flatscale # Flat field scale factor
+real illumscale # Illumination scale factor
+real frgscale # Fringe scale factor
+
+int i, op
+
+begin
+ op = cors[OVERSCAN] + cors[ZEROCOR] + cors[DARKCOR] + cors[FLATCOR]
+ switch (op) {
+ case O: # overscan
+ do i = 1, n
+ out[i] = out[i] - overscan
+ case Z: # zero level
+ do i = 1, n
+ out[i] = out[i] - zero[i]
+
+ case ZO: # zero level + overscan
+ do i = 1, n
+ out[i] = out[i] - overscan - zero[i]
+
+ case D: # dark count
+ do i = 1, n
+ out[i] = out[i] - darkscale * dark[i]
+ case DO: # dark count + overscan
+ do i = 1, n
+ out[i] = out[i] - overscan - darkscale * dark[i]
+ case DZ: # dark count + zero level
+ do i = 1, n
+ out[i] = out[i] - zero[i] - darkscale * dark[i]
+ case DZO: # dark count + zero level + overscan
+ do i = 1, n
+ out[i] = out[i] - overscan - zero[i] - darkscale * dark[i]
+
+ case F: # flat field
+ do i = 1, n
+ out[i] = out[i] * flatscale / flat[i]
+ case FO: # flat field + overscan
+ do i = 1, n
+ out[i] = (out[i] - overscan) * flatscale / flat[i]
+ case FZ: # flat field + zero level
+ do i = 1, n
+ out[i] = (out[i] - zero[i]) * flatscale / flat[i]
+ case FZO: # flat field + zero level + overscan
+ do i = 1, n
+ out[i] = (out[i] - overscan - zero[i]) * flatscale /
+ flat[i]
+ case FD: # flat field + dark count
+ do i = 1, n
+ out[i] = (out[i] - darkscale * dark[i]) * flatscale / flat[i]
+ case FDO: # flat field + dark count + overscan
+ do i = 1, n
+ out[i] = (out[i] - overscan - darkscale * dark[i]) *
+ flatscale / flat[i]
+ case FDZ: # flat field + dark count + zero level
+ do i = 1, n
+ out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+ flatscale / flat[i]
+ case FDZO: # flat field + dark count + zero level + overscan
+ do i = 1, n
+ out[i] = (out[i] - overscan - zero[i] -
+ darkscale * dark[i]) * flatscale / flat[i]
+ }
+
+ # Often these operations will not be performed so test for no
+ # correction rather than go through the switch.
+
+ op = cors[ILLUMCOR] + cors[FRINGECOR]
+ if (op != 0) {
+ switch (op) {
+ case I: # illumination
+ do i = 1, n
+ out[i] = out[i] * illumscale / illum[i]
+ case Q: # fringe
+ do i = 1, n
+ out[i] = out[i] - frgscale * fringe[i]
+ case QI: # fringe + illumination
+ do i = 1, n
+ out[i] = out[i]*illumscale/illum[i] - frgscale*fringe[i]
+ }
+ }
+end
+
+
+# COR2 -- Correct lines for readout axis 2 (columns). This procedure is
+# more complex than when the readout is along the image lines because the
+# zero level and/or flat field corrections may be single readout column
+# vectors.
+
+procedure cor2r (line, cors, out, overscan, zero, dark, flat, illum,
+ fringe, n, zeroim, flatim, darkscale, flatscale, illumscale, frgscale)
+
+int line # Line to be corrected
+int cors[ARB] # Correction flags
+real out[n] # Output data
+real overscan[n] # Overscan value
+real zero[n] # Zero level correction
+real dark[n] # Dark count correction
+real flat[n] # Flat field correction
+real illum[n] # Illumination correction
+real fringe[n] # Fringe correction
+int n # Number of pixels
+pointer zeroim # Zero level IMIO pointer (NULL if 1D vector)
+pointer flatim # Flat field IMIO pointer (NULL if 1D vector)
+real darkscale # Dark count scale factor
+real flatscale # Flat field scale factor
+real illumscale # Illumination scale factor
+real frgscale # Fringe scale factor
+
+real zeroval
+real flatval
+int i, op
+
+begin
+ op = cors[OVERSCAN] + cors[ZEROCOR] + cors[DARKCOR] + cors[FLATCOR]
+ switch (op) {
+ case O: # overscan
+ do i = 1, n
+ out[i] = out[i] - overscan[i]
+ case Z: # zero level
+ if (zeroim != NULL)
+ do i = 1, n
+ out[i] = out[i] - zero[i]
+ else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = out[i] - zeroval
+ }
+
+ case ZO: # zero level + overscan
+ if (zeroim != NULL)
+ do i = 1, n
+ out[i] = out[i] - overscan[i] - zero[i]
+ else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = out[i] - overscan[i] - zeroval
+ }
+
+ case D: # dark count
+ do i = 1, n
+ out[i] = out[i] - darkscale * dark[i]
+ case DO: # dark count + overscan
+ do i = 1, n
+ out[i] = out[i] - overscan[i] - darkscale * dark[i]
+ case DZ: # dark count + zero level
+ if (zeroim != NULL)
+ do i = 1, n
+ out[i] = out[i] - zero[i] - darkscale * dark[i]
+ else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = out[i] - zeroval - darkscale * dark[i]
+ }
+ case DZO: # dark count + zero level + overscan
+ if (zeroim != NULL)
+ do i = 1, n
+ out[i] = out[i] - overscan[i] - zero[i] -
+ darkscale * dark[i]
+ else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = out[i] - overscan[i] - zeroval -
+ darkscale * dark[i]
+ }
+
+ case F: # flat field
+ if (flatim != NULL) {
+ do i = 1, n
+ out[i] = out[i] * flatscale / flat[i]
+ } else {
+ flatval = flatscale / flat[line]
+ do i = 1, n
+ out[i] = out[i] * flatval
+ }
+ case FO: # flat field + overscan
+ if (flatim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i]) * flatscale / flat[i]
+ } else {
+ flatval = flatscale / flat[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i]) * flatval
+ }
+ case FZ: # flat field + zero level
+ if (flatim != NULL) {
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - zero[i]) * flatscale / flat[i]
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - zeroval) * flatscale / flat[i]
+ }
+ } else {
+ flatval = flatscale / flat[line]
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - zero[i]) * flatval
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - zeroval) * flatval
+ }
+ }
+ case FZO: # flat field + zero level + overscan
+ if (flatim != NULL) {
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zero[i]) *
+ flatscale / flat[i]
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zeroval) *
+ flatscale / flat[i]
+ }
+ } else {
+ flatval = flatscale / flat[line]
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zero[i]) * flatval
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zeroval) * flatval
+ }
+ }
+ case FD: # flat field + dark count
+ if (flatim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - darkscale * dark[i]) * flatscale/flat[i]
+ } else {
+ flatval = flatscale / flat[line]
+ do i = 1, n
+ out[i] = (out[i] - darkscale * dark[i]) * flatval
+ }
+ case FDO: # flat field + dark count + overscan
+ if (flatim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - darkscale * dark[i]) *
+ flatscale / flat[i]
+ } else {
+ flatval = flatscale / flat[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - darkscale * dark[i]) *
+ flatval
+ }
+ case FDZ: # flat field + dark count + zero level
+ if (flatim != NULL) {
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+ flatscale / flat[i]
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - zeroval - darkscale * dark[i]) *
+ flatscale / flat[i]
+ }
+ } else {
+ flatval = flatscale / flat[line]
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+ flatval
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - zeroval - darkscale * dark[i]) *
+ flatval
+ }
+ }
+ case FDZO: # flat field + dark count + zero level + overscan
+ if (flatim != NULL) {
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zero[i] -
+ darkscale * dark[i]) * flatscale / flat[i]
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zeroval -
+ darkscale * dark[i]) * flatscale / flat[i]
+ }
+ } else {
+ flatval = flatscale / flat[line]
+ if (zeroim != NULL) {
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zero[i] -
+ darkscale * dark[i]) * flatval
+ } else {
+ zeroval = zero[line]
+ do i = 1, n
+ out[i] = (out[i] - overscan[i] - zeroval -
+ darkscale * dark[i]) * flatval
+ }
+ }
+ }
+
+ # Often these operations will not be performed so test for no
+ # correction rather than go through the switch.
+
+ op = cors[ILLUMCOR] + cors[FRINGECOR]
+ if (op != 0) {
+ switch (op) {
+ case I: # illumination
+ do i = 1, n
+ out[i] = out[i] * illumscale / illum[i]
+ case Q: # fringe
+ do i = 1, n
+ out[i] = out[i] - frgscale * fringe[i]
+ case QI: # fringe + illumination
+ do i = 1, n
+ out[i] = out[i]*illumscale/illum[i] - frgscale*fringe[i]
+ }
+ }
+end
diff --git a/noao/imred/ccdred/src/generic/icaclip.x b/noao/imred/ccdred/src/generic/icaclip.x
new file mode 100644
index 00000000..1530145c
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icaclip.x
@@ -0,0 +1,1102 @@
+# 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_aavsigclips (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 = Mems[d[1]+k]
+ high = Mems[d[2]+k]
+ if (low > high) {
+ d1 = low
+ low = high
+ high = d1
+ }
+ sum = 0.
+ do j = 3, n1 {
+ d1 = Mems[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 = Mems[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 + (Mems[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 = Mems[d[1]+k]
+ do j = 2, n1
+ sum = sum + Mems[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 = Mems[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
+ Mems[dp1] = Mems[dp2]
+ Mems[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 = Mems[dp1]
+ r = (d1 - a) / s1
+ if (r < -lsigma || r > hsigma) {
+ Memr[resid+n1] = abs(r)
+ if (j < n1) {
+ dp2 = d[n1] + k
+ Mems[dp1] = Mems[dp2]
+ Mems[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 = Mems[dp1]
+ Mems[dp1] = Mems[dp2]
+ Mems[dp2] = d1
+ mp2 = m[l] + k
+ s = Memi[mp1]
+ Memi[mp1] = Memi[mp2]
+ Memi[mp2] = s
+ }
+ }
+ }
+ sum = sum + Mems[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 = Mems[dp1]
+ Mems[dp1] = Mems[dp2]
+ Mems[dp2] = d1
+ if (keepids) {
+ mp1 = m[j] + k
+ mp2 = m[l] + k
+ s = Memi[mp1]
+ Memi[mp1] = Memi[mp2]
+ Memi[mp2] = s
+ }
+ }
+ }
+ }
+ sum = sum + Mems[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_mavsigclips (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] = Mems[d[1]+k]
+ else {
+ low = Mems[d[1]+k]
+ high = Mems[d[2]+k]
+ median[i] = (low + high) / 2.
+ }
+ next
+ }
+
+ # Median
+ n3 = 1 + n1 / 2
+ if (mod (n1, 2) == 0) {
+ low = Mems[d[n3-1]+k]
+ high = Mems[d[n3]+k]
+ med = (low + high) / 2.
+ } else
+ med = Mems[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 + (Mems[d[j]+k] - med) ** 2 / s1
+ }
+ } else {
+ s1 = max (one, med)
+ do j = 1, n1
+ s = s + (Mems[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 - Mems[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 = (Mems[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 - Mems[d[nl]+k]) / s1
+ if (r <= lsigma)
+ break
+ Memr[resid+nl] = r
+ n1 = n1 - 1
+ }
+ for (; nh >= nl; nh = nh - 1) {
+ r = (Mems[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 = Mems[d[n3-1]+k]
+ high = Mems[d[n3]+k]
+ med = (low + high) / 2.
+ } else
+ med = Mems[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 = Mems[d[n3-1]+k]
+ high = Mems[d[n3]+k]
+ med = (low + high) / 2.
+ } else
+ med = Mems[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) {
+ Mems[d[l]+k] = Mems[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) {
+ Mems[d[l]+k] = Mems[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
+
+# 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/imred/ccdred/src/generic/icaverage.x b/noao/imred/ccdred/src/generic/icaverage.x
new file mode 100644
index 00000000..3646b725
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icaverage.x
@@ -0,0 +1,163 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include "../icombine.h"
+
+
+# IC_AVERAGE -- Compute the average image line.
+# Options include a weight average.
+
+procedure ic_averages (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 = Mems[d[1]+k] * wt
+ do j = 2, n[i] {
+ wt = wts[Memi[m[j]+k]]
+ sum = sum + Mems[d[j]+k] * wt
+ }
+ average[i] = sum
+ }
+ } else {
+ do i = 1, npts {
+ k = i - 1
+ sum = Mems[d[1]+k]
+ do j = 2, n[i]
+ sum = sum + Mems[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 = Mems[d[1]+k] * wt
+ sumwt = wt
+ do j = 2, n[i] {
+ wt = wts[Memi[m[j]+k]]
+ sum = sum + Mems[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 = Mems[d[1]+k]
+ do j = 2, n[i]
+ sum = sum + Mems[d[j]+k]
+ average[i] = sum / n[i]
+ } else
+ average[i] = blank
+ }
+ }
+ }
+end
+
+# 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/imred/ccdred/src/generic/iccclip.x b/noao/imred/ccdred/src/generic/iccclip.x
new file mode 100644
index 00000000..57709064
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/iccclip.x
@@ -0,0 +1,898 @@
+# 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_accdclips (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 = Mems[d[1]+k]
+ do j = 2, n1
+ sum = sum + Mems[d[j]+k]
+ average[i] = sum / n1
+ }
+ }
+ next
+ }
+
+ repeat {
+ if (n1 == 2) {
+ sum = Mems[d[1]+k]
+ sum = sum + Mems[d[2]+k]
+ a = sum / 2
+ } else {
+ low = Mems[d[1]+k]
+ high = Mems[d[2]+k]
+ if (low > high) {
+ d1 = low
+ low = high
+ high = d1
+ }
+ sum = 0.
+ do j = 3, n1 {
+ d1 = Mems[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 = Mems[dp1]
+ r = (d1 - a) / s
+ if (r < -lsigma || r > hsigma) {
+ Memr[resid+n1] = abs(r)
+ if (j < n1) {
+ dp2 = d[n1] + k
+ Mems[dp1] = Mems[dp2]
+ Mems[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 = Mems[dp1]
+ r = (d1 - a) / s
+ if (r < -lsigma || r > hsigma) {
+ Memr[resid+n1] = abs(r)
+ if (j < n1) {
+ dp2 = d[n1] + k
+ Mems[dp1] = Mems[dp2]
+ Mems[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 = Mems[dp1]
+ Mems[dp1] = Mems[dp2]
+ Mems[dp2] = d1
+ mp2 = m[l] + k
+ s = Memi[mp1]
+ Memi[mp1] = Memi[mp2]
+ Memi[mp2] = s
+ }
+ }
+ }
+ sum = sum + Mems[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 = Mems[dp1]
+ Mems[dp1] = Mems[dp2]
+ Mems[dp2] = d1
+ if (keepids) {
+ mp1 = m[j] + k
+ mp2 = m[l] + k
+ s = Memi[mp1]
+ Memi[mp1] = Memi[mp2]
+ Memi[mp2] = s
+ }
+ }
+ }
+ }
+ sum = sum + Mems[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_mccdclips (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 = Mems[d[n3-1]+k]
+ med = (med + Mems[d[n3]+k]) / 2.
+ } else
+ med = Mems[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 - Mems[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 = (Mems[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 - Mems[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 = (Mems[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) {
+ Mems[d[l]+k] = Mems[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) {
+ Mems[d[l]+k] = Mems[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
+
+# 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/imred/ccdred/src/generic/icgdata.x b/noao/imred/ccdred/src/generic/icgdata.x
new file mode 100644
index 00000000..5c6ac18c
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icgdata.x
@@ -0,0 +1,459 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include <mach.h>
+include "../icombine.h"
+
+
+# IC_GDATA -- Get line of image and mask data and apply threshold and scaling.
+# Entirely empty lines are excluded. The data are compacted within the
+# input data buffers. If it is required, the connection to the original
+# image index is keeped in the returned m data pointers.
+
+procedure ic_gdatas (in, out, dbuf, d, id, n, m, lflag, offsets, scales,
+ zeros, nimages, npts, v1, v2)
+
+pointer in[nimages] # Input images
+pointer out[ARB] # Output images
+pointer dbuf[nimages] # Data buffers for nonaligned images
+pointer d[nimages] # Data pointers
+pointer id[nimages] # ID pointers
+int n[npts] # Number of good pixels
+pointer m[nimages] # Mask pointers
+int lflag[nimages] # Empty mask flags
+int offsets[nimages,ARB] # Image offsets
+real scales[nimages] # Scale factors
+real zeros[nimages] # Zero offset factors
+int nimages # Number of input images
+int npts # NUmber of output points per line
+long v1[ARB], v2[ARB] # Line vectors
+
+int i, j, k, l, ndim, nused
+real a, b
+pointer buf, dp, ip, mp, imgnls()
+
+include "../icombine.com"
+
+begin
+ # Get masks and return if there is no data
+ call ic_mget (in, out, offsets, v1, v2, m, lflag, nimages)
+ if (dflag == D_NONE)
+ return
+
+ # Get data and fill data buffers. Correct for offsets if needed.
+ ndim = IM_NDIM(out[1])
+ do i = 1, nimages {
+ if (lflag[i] == D_NONE)
+ next
+ if (aligned) {
+ call amovl (v1, v2, IM_MAXDIM)
+ if (project)
+ v2[ndim+1] = i
+ j = imgnls (in[i], d[i], v2)
+ } else {
+ v2[1] = v1[1]
+ do j = 2, ndim
+ v2[j] = v1[j] - offsets[i,j]
+ if (project)
+ v2[ndim+1] = i
+ j = imgnls (in[i], buf, v2)
+ call amovs (Mems[buf], Mems[dbuf[i]+offsets[i,1]],
+ IM_LEN(in[i],1))
+ d[i] = dbuf[i]
+ }
+ }
+
+ # Apply threshold if needed
+ if (dothresh) {
+ do i = 1, nimages {
+ dp = d[i]
+ if (lflag[i] == D_ALL) {
+ do j = 1, npts {
+ a = Mems[dp]
+ if (a < lthresh || a > hthresh) {
+ Memi[m[i]+j-1] = 1
+ lflag[i] = D_MIX
+ dflag = D_MIX
+ }
+ dp = dp + 1
+ }
+ } else if (lflag[i] == D_MIX) {
+ mp = m[i]
+ do j = 1, npts {
+ if (Memi[mp] == 0) {
+ a = Mems[dp]
+ if (a < lthresh || a > hthresh) {
+ Memi[m[i]+j-1] = 1
+ dflag = D_MIX
+ }
+ }
+ dp = dp + 1
+ mp = mp + 1
+ }
+ }
+
+ # Check for completely empty lines
+ if (lflag[i] == D_MIX) {
+ lflag[i] = D_NONE
+ mp = m[i]
+ do j = 1, npts {
+ if (Memi[mp] == 0) {
+ lflag[i] = D_MIX
+ break
+ }
+ mp = mp + 1
+ }
+ }
+ }
+ }
+
+ # Apply scaling (avoiding masked pixels which might overflow?)
+ if (doscale) {
+ if (dflag == D_ALL) {
+ do i = 1, nimages {
+ dp = d[i]
+ a = scales[i]
+ b = -zeros[i]
+ do j = 1, npts {
+ Mems[dp] = Mems[dp] / a + b
+ dp = dp + 1
+ }
+ }
+ } else if (dflag == D_MIX) {
+ do i = 1, nimages {
+ dp = d[i]
+ a = scales[i]
+ b = -zeros[i]
+ if (lflag[i] == D_ALL) {
+ do j = 1, npts {
+ Mems[dp] = Mems[dp] / a + b
+ dp = dp + 1
+ }
+ } else if (lflag[i] == D_MIX) {
+ mp = m[i]
+ do j = 1, npts {
+ if (Memi[mp] == 0)
+ Mems[dp] = Mems[dp] / a + b
+ dp = dp + 1
+ mp = mp + 1
+ }
+ }
+ }
+ }
+ }
+
+ # Sort pointers to exclude unused images.
+ # Use the lflag array to keep track of the image index.
+
+ if (dflag == D_ALL)
+ nused = nimages
+ else {
+ nused = 0
+ do i = 1, nimages
+ if (lflag[i] != D_NONE) {
+ nused = nused + 1
+ d[nused] = d[i]
+ m[nused] = m[i]
+ lflag[nused] = i
+ }
+ if (nused == 0)
+ dflag = D_NONE
+ }
+
+ # Compact data to remove bad pixels
+ # Keep track of the image indices if needed
+ # If growing mark the end of the included image indices with zero
+
+ if (dflag == D_ALL) {
+ call amovki (nused, n, npts)
+ if (keepids)
+ do i = 1, nimages
+ call amovki (i, Memi[id[i]], npts)
+ } else if (dflag == D_NONE)
+ call aclri (n, npts)
+ else {
+ call aclri (n, npts)
+ if (keepids) {
+ do i = 1, nused {
+ l = lflag[i]
+ dp = d[i]
+ ip = id[i]
+ mp = m[i]
+ do j = 1, npts {
+ if (Memi[mp] == 0) {
+ n[j] = n[j] + 1
+ k = n[j]
+ if (k < i) {
+ Mems[d[k]+j-1] = Mems[dp]
+ Memi[id[k]+j-1] = l
+ } else
+ Memi[ip] = l
+ }
+ dp = dp + 1
+ ip = ip + 1
+ mp = mp + 1
+ }
+ }
+ if (grow > 0) {
+ do j = 1, npts {
+ do i = n[j]+1, nimages
+ Memi[id[i]+j-1] = 0
+ }
+ }
+ } else {
+ do i = 1, nused {
+ dp = d[i]
+ mp = m[i]
+ do j = 1, npts {
+ if (Memi[mp] == 0) {
+ n[j] = n[j] + 1
+ k = n[j]
+ if (k < i)
+ Mems[d[k]+j-1] = Mems[dp]
+ }
+ dp = dp + 1
+ mp = mp + 1
+ }
+ }
+ }
+ }
+
+ # Sort the pixels and IDs if needed
+ if (mclip) {
+ call malloc (dp, nimages, TY_SHORT)
+ if (keepids) {
+ call malloc (ip, nimages, TY_INT)
+ call ic_2sorts (d, Mems[dp], id, Memi[ip], n, npts)
+ call mfree (ip, TY_INT)
+ } else
+ call ic_sorts (d, Mems[dp], n, npts)
+ call mfree (dp, TY_SHORT)
+ }
+end
+
+# IC_GDATA -- Get line of image and mask data and apply threshold and scaling.
+# Entirely empty lines are excluded. The data are compacted within the
+# input data buffers. If it is required, the connection to the original
+# image index is keeped in the returned m data pointers.
+
+procedure ic_gdatar (in, out, dbuf, d, id, n, m, lflag, offsets, scales,
+ zeros, nimages, npts, v1, v2)
+
+pointer in[nimages] # Input images
+pointer out[ARB] # Output images
+pointer dbuf[nimages] # Data buffers for nonaligned images
+pointer d[nimages] # Data pointers
+pointer id[nimages] # ID pointers
+int n[npts] # Number of good pixels
+pointer m[nimages] # Mask pointers
+int lflag[nimages] # Empty mask flags
+int offsets[nimages,ARB] # Image offsets
+real scales[nimages] # Scale factors
+real zeros[nimages] # Zero offset factors
+int nimages # Number of input images
+int npts # NUmber of output points per line
+long v1[ARB], v2[ARB] # Line vectors
+
+int i, j, k, l, ndim, nused
+real a, b
+pointer buf, dp, ip, mp, imgnlr()
+
+include "../icombine.com"
+
+begin
+ # Get masks and return if there is no data
+ call ic_mget (in, out, offsets, v1, v2, m, lflag, nimages)
+ if (dflag == D_NONE)
+ return
+
+ # Get data and fill data buffers. Correct for offsets if needed.
+ ndim = IM_NDIM(out[1])
+ do i = 1, nimages {
+ if (lflag[i] == D_NONE)
+ next
+ if (aligned) {
+ call amovl (v1, v2, IM_MAXDIM)
+ if (project)
+ v2[ndim+1] = i
+ j = imgnlr (in[i], d[i], v2)
+ } else {
+ v2[1] = v1[1]
+ do j = 2, ndim
+ v2[j] = v1[j] - offsets[i,j]
+ if (project)
+ v2[ndim+1] = i
+ j = imgnlr (in[i], buf, v2)
+ call amovr (Memr[buf], Memr[dbuf[i]+offsets[i,1]],
+ IM_LEN(in[i],1))
+ d[i] = dbuf[i]
+ }
+ }
+
+ # Apply threshold if needed
+ if (dothresh) {
+ do i = 1, nimages {
+ dp = d[i]
+ if (lflag[i] == D_ALL) {
+ do j = 1, npts {
+ a = Memr[dp]
+ if (a < lthresh || a > hthresh) {
+ Memi[m[i]+j-1] = 1
+ lflag[i] = D_MIX
+ dflag = D_MIX
+ }
+ dp = dp + 1
+ }
+ } else if (lflag[i] == D_MIX) {
+ mp = m[i]
+ do j = 1, npts {
+ if (Memi[mp] == 0) {
+ a = Memr[dp]
+ if (a < lthresh || a > hthresh) {
+ Memi[m[i]+j-1] = 1
+ dflag = D_MIX
+ }
+ }
+ dp = dp + 1
+ mp = mp + 1
+ }
+ }
+
+ # Check for completely empty lines
+ if (lflag[i] == D_MIX) {
+ lflag[i] = D_NONE
+ mp = m[i]
+ do j = 1, npts {
+ if (Memi[mp] == 0) {
+ lflag[i] = D_MIX
+ break
+ }
+ mp = mp + 1
+ }
+ }
+ }
+ }
+
+ # Apply scaling (avoiding masked pixels which might overflow?)
+ if (doscale) {
+ if (dflag == D_ALL) {
+ do i = 1, nimages {
+ dp = d[i]
+ a = scales[i]
+ b = -zeros[i]
+ do j = 1, npts {
+ Memr[dp] = Memr[dp] / a + b
+ dp = dp + 1
+ }
+ }
+ } else if (dflag == D_MIX) {
+ do i = 1, nimages {
+ dp = d[i]
+ a = scales[i]
+ b = -zeros[i]
+ if (lflag[i] == D_ALL) {
+ do j = 1, npts {
+ Memr[dp] = Memr[dp] / a + b
+ dp = dp + 1
+ }
+ } else if (lflag[i] == D_MIX) {
+ mp = m[i]
+ do j = 1, npts {
+ if (Memi[mp] == 0)
+ Memr[dp] = Memr[dp] / a + b
+ dp = dp + 1
+ mp = mp + 1
+ }
+ }
+ }
+ }
+ }
+
+ # Sort pointers to exclude unused images.
+ # Use the lflag array to keep track of the image index.
+
+ if (dflag == D_ALL)
+ nused = nimages
+ else {
+ nused = 0
+ do i = 1, nimages
+ if (lflag[i] != D_NONE) {
+ nused = nused + 1
+ d[nused] = d[i]
+ m[nused] = m[i]
+ lflag[nused] = i
+ }
+ if (nused == 0)
+ dflag = D_NONE
+ }
+
+ # Compact data to remove bad pixels
+ # Keep track of the image indices if needed
+ # If growing mark the end of the included image indices with zero
+
+ if (dflag == D_ALL) {
+ call amovki (nused, n, npts)
+ if (keepids)
+ do i = 1, nimages
+ call amovki (i, Memi[id[i]], npts)
+ } else if (dflag == D_NONE)
+ call aclri (n, npts)
+ else {
+ call aclri (n, npts)
+ if (keepids) {
+ do i = 1, nused {
+ l = lflag[i]
+ dp = d[i]
+ ip = id[i]
+ mp = m[i]
+ do j = 1, npts {
+ if (Memi[mp] == 0) {
+ n[j] = n[j] + 1
+ k = n[j]
+ if (k < i) {
+ Memr[d[k]+j-1] = Memr[dp]
+ Memi[id[k]+j-1] = l
+ } else
+ Memi[ip] = l
+ }
+ dp = dp + 1
+ ip = ip + 1
+ mp = mp + 1
+ }
+ }
+ if (grow > 0) {
+ do j = 1, npts {
+ do i = n[j]+1, nimages
+ Memi[id[i]+j-1] = 0
+ }
+ }
+ } else {
+ do i = 1, nused {
+ dp = d[i]
+ mp = m[i]
+ do j = 1, npts {
+ if (Memi[mp] == 0) {
+ n[j] = n[j] + 1
+ k = n[j]
+ if (k < i)
+ Memr[d[k]+j-1] = Memr[dp]
+ }
+ dp = dp + 1
+ mp = mp + 1
+ }
+ }
+ }
+ }
+
+ # Sort the pixels and IDs if needed
+ if (mclip) {
+ call malloc (dp, nimages, TY_REAL)
+ if (keepids) {
+ call malloc (ip, nimages, TY_INT)
+ call ic_2sortr (d, Memr[dp], id, Memi[ip], n, npts)
+ call mfree (ip, TY_INT)
+ } else
+ call ic_sortr (d, Memr[dp], n, npts)
+ call mfree (dp, TY_REAL)
+ }
+end
+
diff --git a/noao/imred/ccdred/src/generic/icgrow.x b/noao/imred/ccdred/src/generic/icgrow.x
new file mode 100644
index 00000000..b94e1cbc
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icgrow.x
@@ -0,0 +1,148 @@
+# 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_grows (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] - Mems[d[j2]+k2]) / (n2-1)
+ mp2 = m[n2] + k2
+ if (j2 < n2) {
+ Mems[d[j2]+k2] = Mems[d[n2]+k2]
+ Memi[mp1] = Memi[mp2]
+ }
+ Memi[mp2] = 0
+ n[i2] = n2 - 1
+ break
+ }
+ }
+ }
+ }
+ }
+end
+
+# 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/imred/ccdred/src/generic/icmedian.x b/noao/imred/ccdred/src/generic/icmedian.x
new file mode 100644
index 00000000..ec0166ba
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icmedian.x
@@ -0,0 +1,343 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "../icombine.h"
+
+
+# IC_MEDIAN -- Median of lines
+
+procedure ic_medians (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, j, k, j1, j2, n1, lo, up, lo1, up1
+bool even
+real val1, val2, val3
+short temp, wtemp
+
+include "../icombine.com"
+
+begin
+ # If no data return after possibly setting blank values.
+ if (dflag == D_NONE) {
+ do i = 1, npts
+ median[i]= blank
+ return
+ }
+
+ # If the data were previously sorted then directly compute the median.
+ 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 = Mems[d[j1]+k]
+ val2 = Mems[d[j2]+k]
+ median[i] = (val1 + val2) / 2.
+ } else
+ median[i] = Mems[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 = Mems[d[j1]+k]
+ val2 = Mems[d[j2]+k]
+ median[i] = (val1 + val2) / 2.
+ } else
+ median[i] = Mems[d[j1]+k]
+ } else
+ median[i] = blank
+ }
+ }
+ return
+ }
+
+ # Compute the median.
+ do i = 1, npts {
+ k = i - 1
+ n1 = n[i]
+
+ # If there are more than 3 points use Wirth algorithm. This
+ # is the same as vops$amed.gx except for an even number of
+ # points it selects the middle two and averages.
+ if (n1 > 3) {
+ lo = 1
+ up = n1
+ j = max (lo, min (up, (up+1)/2))
+
+ while (lo < up) {
+ if (! (lo < up))
+ break
+
+ temp = Mems[d[j]+k]; lo1 = lo; up1 = up
+
+ repeat {
+ while (Mems[d[lo1]+k] < temp)
+ lo1 = lo1 + 1
+ while (temp < Mems[d[up1]+k])
+ up1 = up1 - 1
+ if (lo1 <= up1) {
+ wtemp = Mems[d[lo1]+k]
+ Mems[d[lo1]+k] = Mems[d[up1]+k]
+ Mems[d[up1]+k] = wtemp
+ lo1 = lo1 + 1; up1 = up1 - 1
+ }
+ } until (lo1 > up1)
+
+ if (up1 < j)
+ lo = lo1
+ if (j < lo1)
+ up = up1
+ }
+
+ median[i] = Mems[d[j]+k]
+
+ if (mod (n1,2) == 0) {
+ lo = 1
+ up = n1
+ j = max (lo, min (up, (up+1)/2)+1)
+
+ while (lo < up) {
+ if (! (lo < up))
+ break
+
+ temp = Mems[d[j]+k]; lo1 = lo; up1 = up
+
+ repeat {
+ while (Mems[d[lo1]+k] < temp)
+ lo1 = lo1 + 1
+ while (temp < Mems[d[up1]+k])
+ up1 = up1 - 1
+ if (lo1 <= up1) {
+ wtemp = Mems[d[lo1]+k]
+ Mems[d[lo1]+k] = Mems[d[up1]+k]
+ Mems[d[up1]+k] = wtemp
+ lo1 = lo1 + 1; up1 = up1 - 1
+ }
+ } until (lo1 > up1)
+
+ if (up1 < j)
+ lo = lo1
+ if (j < lo1)
+ up = up1
+ }
+ median[i] = (median[i] + Mems[d[j]+k]) / 2
+ }
+
+ # If 3 points find the median directly.
+ } else if (n1 == 3) {
+ val1 = Mems[d[1]+k]
+ val2 = Mems[d[2]+k]
+ val3 = Mems[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
+ }
+
+ # If 2 points average.
+ } else if (n1 == 2) {
+ val1 = Mems[d[1]+k]
+ val2 = Mems[d[2]+k]
+ median[i] = (val1 + val2) / 2
+
+ # If 1 point return the value.
+ } else if (n1 == 1)
+ median[i] = Mems[d[1]+k]
+
+ # If no points return with a possibly blank value.
+ else
+ median[i] = blank
+ }
+end
+
+# 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, j, k, j1, j2, n1, lo, up, lo1, up1
+bool even
+real val1, val2, val3
+real temp, wtemp
+
+include "../icombine.com"
+
+begin
+ # If no data return after possibly setting blank values.
+ if (dflag == D_NONE) {
+ do i = 1, npts
+ median[i]= blank
+ return
+ }
+
+ # If the data were previously sorted then directly compute the median.
+ 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
+ }
+
+ # Compute the median.
+ do i = 1, npts {
+ k = i - 1
+ n1 = n[i]
+
+ # If there are more than 3 points use Wirth algorithm. This
+ # is the same as vops$amed.gx except for an even number of
+ # points it selects the middle two and averages.
+ if (n1 > 3) {
+ lo = 1
+ up = n1
+ j = max (lo, min (up, (up+1)/2))
+
+ while (lo < up) {
+ if (! (lo < up))
+ break
+
+ temp = Memr[d[j]+k]; lo1 = lo; up1 = up
+
+ repeat {
+ while (Memr[d[lo1]+k] < temp)
+ lo1 = lo1 + 1
+ while (temp < Memr[d[up1]+k])
+ up1 = up1 - 1
+ if (lo1 <= up1) {
+ wtemp = Memr[d[lo1]+k]
+ Memr[d[lo1]+k] = Memr[d[up1]+k]
+ Memr[d[up1]+k] = wtemp
+ lo1 = lo1 + 1; up1 = up1 - 1
+ }
+ } until (lo1 > up1)
+
+ if (up1 < j)
+ lo = lo1
+ if (j < lo1)
+ up = up1
+ }
+
+ median[i] = Memr[d[j]+k]
+
+ if (mod (n1,2) == 0) {
+ lo = 1
+ up = n1
+ j = max (lo, min (up, (up+1)/2)+1)
+
+ while (lo < up) {
+ if (! (lo < up))
+ break
+
+ temp = Memr[d[j]+k]; lo1 = lo; up1 = up
+
+ repeat {
+ while (Memr[d[lo1]+k] < temp)
+ lo1 = lo1 + 1
+ while (temp < Memr[d[up1]+k])
+ up1 = up1 - 1
+ if (lo1 <= up1) {
+ wtemp = Memr[d[lo1]+k]
+ Memr[d[lo1]+k] = Memr[d[up1]+k]
+ Memr[d[up1]+k] = wtemp
+ lo1 = lo1 + 1; up1 = up1 - 1
+ }
+ } until (lo1 > up1)
+
+ if (up1 < j)
+ lo = lo1
+ if (j < lo1)
+ up = up1
+ }
+ median[i] = (median[i] + Memr[d[j]+k]) / 2
+ }
+
+ # If 3 points find the median directly.
+ } else 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
+ }
+
+ # If 2 points average.
+ } else if (n1 == 2) {
+ val1 = Memr[d[1]+k]
+ val2 = Memr[d[2]+k]
+ median[i] = (val1 + val2) / 2
+
+ # If 1 point return the value.
+ } else if (n1 == 1)
+ median[i] = Memr[d[1]+k]
+
+ # If no points return with a possibly blank value.
+ else
+ median[i] = blank
+ }
+end
+
diff --git a/noao/imred/ccdred/src/generic/icmm.x b/noao/imred/ccdred/src/generic/icmm.x
new file mode 100644
index 00000000..259759bd
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icmm.x
@@ -0,0 +1,300 @@
+# 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_mms (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
+short 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 = Mems[k]
+ dmax = d1; dmin = d1; jmax = 1; jmin = 1; kmax = k; kmin = k
+ do j = 2, n1 {
+ d2 = d1
+ k = d[j] + i1
+ d1 = Mems[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) {
+ Mems[kmax] = d2
+ Memi[m[jmax]+i1] = Memi[m[j]+i1]
+ } else {
+ Mems[kmax] = d1
+ Memi[m[jmax]+i1] = Memi[m[n1]+i1]
+ }
+ }
+ if (jmin < j) {
+ if (jmax != n1) {
+ Mems[kmin] = d1
+ Memi[m[jmin]+i1] = Memi[m[n1]+i1]
+ } else {
+ Mems[kmin] = d2
+ Memi[m[jmin]+i1] = Memi[m[j]+i1]
+ }
+ }
+ } else {
+ if (jmax < j) {
+ if (jmin != j)
+ Mems[kmax] = d2
+ else
+ Mems[kmax] = d1
+ }
+ if (jmin < j) {
+ if (jmax != n1)
+ Mems[kmin] = d1
+ else
+ Mems[kmin] = d2
+ }
+ }
+ n1 = n1 - 2
+ }
+
+ # Reject the excess low points.
+ do np = 1, nlow {
+ k = d[1] + i1
+ d1 = Mems[k]
+ dmin = d1; jmin = 1; kmin = k
+ do j = 2, n1 {
+ k = d[j] + i1
+ d1 = Mems[k]
+ if (d1 < dmin) {
+ dmin = d1; jmin = j; kmin = k
+ }
+ }
+ if (keepids) {
+ if (jmin < n1) {
+ Mems[kmin] = d1
+ Memi[m[jmin]+i1] = Memi[m[n1]+i1]
+ }
+ } else {
+ if (jmin < n1)
+ Mems[kmin] = d1
+ }
+ n1 = n1 - 1
+ }
+
+ # Reject the excess high points.
+ do np = 1, nhigh {
+ k = d[1] + i1
+ d1 = Mems[k]
+ dmax = d1; jmax = 1; kmax = k
+ do j = 2, n1 {
+ k = d[j] + i1
+ d1 = Mems[k]
+ if (d1 > dmax) {
+ dmax = d1; jmax = j; kmax = k
+ }
+ }
+ if (keepids) {
+ if (jmax < n1) {
+ Mems[kmax] = d1
+ Memi[m[jmax]+i1] = Memi[m[n1]+i1]
+ }
+ } else {
+ if (jmax < n1)
+ Mems[kmax] = d1
+ }
+ n1 = n1 - 1
+ }
+ n[i] = n1
+ }
+
+ if (dflag == D_ALL && npairs + nlow + nhigh > 0)
+ dflag = D_MIX
+end
+
+# 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/imred/ccdred/src/generic/icombine.x b/noao/imred/ccdred/src/generic/icombine.x
new file mode 100644
index 00000000..b4ff60be
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icombine.x
@@ -0,0 +1,607 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include <imset.h>
+include <error.h>
+include <syserr.h>
+include <mach.h>
+include "../icombine.h"
+
+
+# ICOMBINE -- Combine images
+#
+# The memory and open file descriptor limits are checked and an attempt
+# to recover is made either by setting the image pixel files to be
+# closed after I/O or by notifying the calling program that memory
+# ran out and the IMIO buffer size should be reduced. After the checks
+# a procedure for the selected combine option is called.
+# Because there may be several failure modes when reaching the file
+# limits we first assume an error is due to the file limit, except for
+# out of memory, and close some pixel files. If the error then repeats
+# on accessing the pixels the error is passed back.
+
+
+procedure icombines (in, out, offsets, nimages, bufsize)
+
+pointer in[nimages] # Input images
+pointer out[ARB] # Output images
+int offsets[nimages,ARB] # Input image offsets
+int nimages # Number of input images
+int bufsize # IMIO buffer size
+
+char str[1]
+int i, j, npts, fd, stropen(), errcode(), imstati()
+pointer sp, d, id, n, m, lflag, scales, zeros, wts, dbuf
+pointer buf, imgl1s(), impl1i()
+errchk stropen, imgl1s, impl1i
+pointer impl1r()
+errchk impl1r
+
+include "../icombine.com"
+
+begin
+ npts = IM_LEN(out[1],1)
+
+ # Allocate memory.
+ call smark (sp)
+ call salloc (d, nimages, TY_POINTER)
+ call salloc (id, nimages, TY_POINTER)
+ call salloc (n, npts, TY_INT)
+ call salloc (m, nimages, TY_POINTER)
+ call salloc (lflag, nimages, TY_INT)
+ call salloc (scales, nimages, TY_REAL)
+ call salloc (zeros, nimages, TY_REAL)
+ call salloc (wts, nimages, TY_REAL)
+ call amovki (D_ALL, Memi[lflag], nimages)
+
+ # If aligned use the IMIO buffer otherwise we need vectors of
+ # output length.
+
+ if (!aligned) {
+ call salloc (dbuf, nimages, TY_POINTER)
+ do i = 1, nimages
+ call salloc (Memi[dbuf+i-1], npts, TY_SHORT)
+ }
+
+ if (project) {
+ call imseti (in[1], IM_NBUFS, nimages)
+ call imseti (in[1], IM_BUFSIZE, bufsize)
+ do i = 1, 3 {
+ if (out[i] != NULL)
+ call imseti (out[i], IM_BUFSIZE, bufsize)
+ }
+ } else {
+ # Reserve FD for string operations.
+ fd = stropen (str, 1, NEW_FILE)
+
+ # Do I/O to the images.
+ do i = 1, 3 {
+ if (out[i] != NULL)
+ call imseti (out[i], IM_BUFSIZE, bufsize)
+ }
+ buf = impl1r (out[1])
+ call aclrr (Memr[buf], npts)
+ if (out[3] != NULL) {
+ buf = impl1r (out[3])
+ call aclrr (Memr[buf], npts)
+ }
+ if (out[2] != NULL) {
+ buf = impl1i (out[2])
+ call aclri (Memi[buf], npts)
+ }
+
+ do i = 1, nimages {
+ call imseti (in[i], IM_BUFSIZE, bufsize)
+ iferr (buf = imgl1s (in[i])) {
+ switch (errcode()) {
+ case SYS_MFULL:
+ call sfree (sp)
+ call strclose (fd)
+ call erract (EA_ERROR)
+ case SYS_FTOOMANYFILES, SYS_IKIOPIX:
+ if (imstati (in[i], IM_CLOSEFD) == YES) {
+ call sfree (sp)
+ call strclose (fd)
+ call erract (EA_ERROR)
+ }
+ do j = i-2, nimages
+ call imseti (in[j], IM_CLOSEFD, YES)
+ buf = imgl1s (in[i])
+ default:
+ call sfree (sp)
+ call strclose (fd)
+ call erract (EA_ERROR)
+ }
+ }
+ }
+
+ call strclose (fd)
+ }
+
+ call ic_combines (in, out, Memi[dbuf], Memi[d], Memi[id], Memi[n],
+ Memi[m], Memi[lflag], offsets, Memr[scales], Memr[zeros],
+ Memr[wts], nimages, npts)
+end
+
+
+# IC_COMBINE -- Combine images.
+
+procedure ic_combines (in, out, dbuf, d, id, n, m, lflag, offsets,
+ scales, zeros, wts, nimages, npts)
+
+pointer in[nimages] # Input images
+pointer out[ARB] # Output image
+pointer dbuf[nimages] # Data buffers for nonaligned images
+pointer d[nimages] # Data pointers
+pointer id[nimages] # Image index ID pointers
+int n[npts] # Number of good pixels
+pointer m[nimages] # Mask pointers
+int lflag[nimages] # Line flags
+int offsets[nimages,ARB] # Input image offsets
+real scales[nimages] # Scale factors
+real zeros[nimages] # Zero offset factors
+real wts[nimages] # Combining weights
+int nimages # Number of input images
+int npts # Number of points per output line
+
+int i, ctor()
+real r, imgetr()
+pointer sp, v1, v2, v3, outdata, buf, nm, impnli()
+pointer impnlr()
+errchk ic_scale, imgetr
+
+include "../icombine.com"
+
+begin
+ call smark (sp)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (v2, IM_MAXDIM, TY_LONG)
+ call salloc (v3, IM_MAXDIM, TY_LONG)
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovkl (long(1), Meml[v2], IM_MAXDIM)
+ call amovkl (long(1), Meml[v3], IM_MAXDIM)
+
+ call ic_scale (in, out, offsets, scales, zeros, wts, nimages)
+
+ # Set combine parameters
+ switch (combine) {
+ case AVERAGE:
+ if (dowts)
+ keepids = true
+ else
+ keepids = false
+ case MEDIAN:
+ dowts = false
+ keepids = false
+ }
+ docombine = true
+
+ # Set rejection algorithm specific parameters
+ switch (reject) {
+ case CCDCLIP, CRREJECT:
+ call salloc (nm, 3*nimages, TY_REAL)
+ i = 1
+ if (ctor (Memc[rdnoise], i, r) > 0) {
+ do i = 1, nimages
+ Memr[nm+3*(i-1)] = r
+ } else {
+ do i = 1, nimages
+ Memr[nm+3*(i-1)] = imgetr (in[i], Memc[rdnoise])
+ }
+ i = 1
+ if (ctor (Memc[gain], i, r) > 0) {
+ do i = 1, nimages {
+ Memr[nm+3*(i-1)+1] = r
+ Memr[nm+3*(i-1)] =
+ max ((Memr[nm+3*(i-1)] / r) ** 2, 1e4 / MAX_REAL)
+ }
+ } else {
+ do i = 1, nimages {
+ r = imgetr (in[i], Memc[gain])
+ Memr[nm+3*(i-1)+1] = r
+ Memr[nm+3*(i-1)] =
+ max ((Memr[nm+3*(i-1)] / r) ** 2, 1e4 / MAX_REAL)
+ }
+ }
+ i = 1
+ if (ctor (Memc[snoise], i, r) > 0) {
+ do i = 1, nimages
+ Memr[nm+3*(i-1)+2] = r
+ } else {
+ do i = 1, nimages {
+ r = imgetr (in[i], Memc[snoise])
+ Memr[nm+3*(i-1)+2] = r
+ }
+ }
+ if (!keepids) {
+ if (doscale1 || grow > 0)
+ keepids = true
+ else {
+ do i = 2, nimages {
+ if (Memr[nm+3*(i-1)] != Memr[nm] ||
+ Memr[nm+3*(i-1)+1] != Memr[nm+1] ||
+ Memr[nm+3*(i-1)+2] != Memr[nm+2]) {
+ keepids = true
+ break
+ }
+ }
+ }
+ }
+ if (reject == CRREJECT)
+ lsigma = MAX_REAL
+ case MINMAX:
+ mclip = false
+ if (grow > 0)
+ keepids = true
+ case PCLIP:
+ mclip = true
+ if (grow > 0)
+ keepids = true
+ case AVSIGCLIP, SIGCLIP:
+ if (doscale1 || grow > 0)
+ keepids = true
+ case NONE:
+ mclip = false
+ grow = 0
+ }
+
+ if (keepids) {
+ do i = 1, nimages
+ call salloc (id[i], npts, TY_INT)
+ }
+
+ while (impnlr (out[1], outdata, Meml[v1]) != EOF) {
+ call ic_gdatas (in, out, dbuf, d, id, n, m, lflag, offsets,
+ scales, zeros, nimages, npts, Meml[v2], Meml[v3])
+
+ switch (reject) {
+ case CCDCLIP, CRREJECT:
+ if (mclip)
+ call ic_mccdclips (d, id, n, scales, zeros, Memr[nm],
+ nimages, npts, Memr[outdata])
+ else
+ call ic_accdclips (d, id, n, scales, zeros, Memr[nm],
+ nimages, npts, Memr[outdata])
+ case MINMAX:
+ call ic_mms (d, id, n, npts)
+ case PCLIP:
+ call ic_pclips (d, id, n, nimages, npts, Memr[outdata])
+ case SIGCLIP:
+ if (mclip)
+ call ic_msigclips (d, id, n, scales, zeros, nimages, npts,
+ Memr[outdata])
+ else
+ call ic_asigclips (d, id, n, scales, zeros, nimages, npts,
+ Memr[outdata])
+ case AVSIGCLIP:
+ if (mclip)
+ call ic_mavsigclips (d, id, n, scales, zeros, nimages,
+ npts, Memr[outdata])
+ else
+ call ic_aavsigclips (d, id, n, scales, zeros, nimages,
+ npts, Memr[outdata])
+ }
+
+ if (grow > 0)
+ call ic_grows (d, id, n, nimages, npts, Memr[outdata])
+
+ if (docombine) {
+ switch (combine) {
+ case AVERAGE:
+ call ic_averages (d, id, n, wts, npts, Memr[outdata])
+ case MEDIAN:
+ call ic_medians (d, n, npts, Memr[outdata])
+ }
+ }
+
+ if (out[2] != NULL) {
+ call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+ i = impnli (out[2], buf, Meml[v1])
+ call amovki (nimages, Memi[buf], npts)
+ call asubi (Memi[buf], n, Memi[buf], npts)
+ }
+
+ if (out[3] != NULL) {
+ call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+ i = impnlr (out[3], buf, Meml[v1])
+ call ic_sigmas (d, id, n, wts, npts, Memr[outdata],
+ Memr[buf])
+ }
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+ }
+
+ call sfree (sp)
+end
+
+procedure icombiner (in, out, offsets, nimages, bufsize)
+
+pointer in[nimages] # Input images
+pointer out[ARB] # Output images
+int offsets[nimages,ARB] # Input image offsets
+int nimages # Number of input images
+int bufsize # IMIO buffer size
+
+char str[1]
+int i, j, npts, fd, stropen(), errcode(), imstati()
+pointer sp, d, id, n, m, lflag, scales, zeros, wts, dbuf
+pointer buf, imgl1r(), impl1i()
+errchk stropen, imgl1r, impl1i
+pointer impl1r()
+errchk impl1r
+
+include "../icombine.com"
+
+begin
+ npts = IM_LEN(out[1],1)
+
+ # Allocate memory.
+ call smark (sp)
+ call salloc (d, nimages, TY_POINTER)
+ call salloc (id, nimages, TY_POINTER)
+ call salloc (n, npts, TY_INT)
+ call salloc (m, nimages, TY_POINTER)
+ call salloc (lflag, nimages, TY_INT)
+ call salloc (scales, nimages, TY_REAL)
+ call salloc (zeros, nimages, TY_REAL)
+ call salloc (wts, nimages, TY_REAL)
+ call amovki (D_ALL, Memi[lflag], nimages)
+
+ # If aligned use the IMIO buffer otherwise we need vectors of
+ # output length.
+
+ if (!aligned) {
+ call salloc (dbuf, nimages, TY_POINTER)
+ do i = 1, nimages
+ call salloc (Memi[dbuf+i-1], npts, TY_REAL)
+ }
+
+ if (project) {
+ call imseti (in[1], IM_NBUFS, nimages)
+ call imseti (in[1], IM_BUFSIZE, bufsize)
+ do i = 1, 3 {
+ if (out[i] != NULL)
+ call imseti (out[i], IM_BUFSIZE, bufsize)
+ }
+ } else {
+ # Reserve FD for string operations.
+ fd = stropen (str, 1, NEW_FILE)
+
+ # Do I/O to the images.
+ do i = 1, 3 {
+ if (out[i] != NULL)
+ call imseti (out[i], IM_BUFSIZE, bufsize)
+ }
+ buf = impl1r (out[1])
+ call aclrr (Memr[buf], npts)
+ if (out[3] != NULL) {
+ buf = impl1r (out[3])
+ call aclrr (Memr[buf], npts)
+ }
+ if (out[2] != NULL) {
+ buf = impl1i (out[2])
+ call aclri (Memi[buf], npts)
+ }
+
+ do i = 1, nimages {
+ call imseti (in[i], IM_BUFSIZE, bufsize)
+ iferr (buf = imgl1r (in[i])) {
+ switch (errcode()) {
+ case SYS_MFULL:
+ call sfree (sp)
+ call strclose (fd)
+ call erract (EA_ERROR)
+ case SYS_FTOOMANYFILES, SYS_IKIOPIX:
+ if (imstati (in[i], IM_CLOSEFD) == YES) {
+ call sfree (sp)
+ call strclose (fd)
+ call erract (EA_ERROR)
+ }
+ do j = i-2, nimages
+ call imseti (in[j], IM_CLOSEFD, YES)
+ buf = imgl1r (in[i])
+ default:
+ call sfree (sp)
+ call strclose (fd)
+ call erract (EA_ERROR)
+ }
+ }
+ }
+
+ call strclose (fd)
+ }
+
+ call ic_combiner (in, out, Memi[dbuf], Memi[d], Memi[id], Memi[n],
+ Memi[m], Memi[lflag], offsets, Memr[scales], Memr[zeros],
+ Memr[wts], nimages, npts)
+end
+
+
+# IC_COMBINE -- Combine images.
+
+procedure ic_combiner (in, out, dbuf, d, id, n, m, lflag, offsets,
+ scales, zeros, wts, nimages, npts)
+
+pointer in[nimages] # Input images
+pointer out[ARB] # Output image
+pointer dbuf[nimages] # Data buffers for nonaligned images
+pointer d[nimages] # Data pointers
+pointer id[nimages] # Image index ID pointers
+int n[npts] # Number of good pixels
+pointer m[nimages] # Mask pointers
+int lflag[nimages] # Line flags
+int offsets[nimages,ARB] # Input image offsets
+real scales[nimages] # Scale factors
+real zeros[nimages] # Zero offset factors
+real wts[nimages] # Combining weights
+int nimages # Number of input images
+int npts # Number of points per output line
+
+int i, ctor()
+real r, imgetr()
+pointer sp, v1, v2, v3, outdata, buf, nm, impnli()
+pointer impnlr()
+errchk ic_scale, imgetr
+
+include "../icombine.com"
+
+begin
+ call smark (sp)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (v2, IM_MAXDIM, TY_LONG)
+ call salloc (v3, IM_MAXDIM, TY_LONG)
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovkl (long(1), Meml[v2], IM_MAXDIM)
+ call amovkl (long(1), Meml[v3], IM_MAXDIM)
+
+ call ic_scale (in, out, offsets, scales, zeros, wts, nimages)
+
+ # Set combine parameters
+ switch (combine) {
+ case AVERAGE:
+ if (dowts)
+ keepids = true
+ else
+ keepids = false
+ case MEDIAN:
+ dowts = false
+ keepids = false
+ }
+ docombine = true
+
+ # Set rejection algorithm specific parameters
+ switch (reject) {
+ case CCDCLIP, CRREJECT:
+ call salloc (nm, 3*nimages, TY_REAL)
+ i = 1
+ if (ctor (Memc[rdnoise], i, r) > 0) {
+ do i = 1, nimages
+ Memr[nm+3*(i-1)] = r
+ } else {
+ do i = 1, nimages
+ Memr[nm+3*(i-1)] = imgetr (in[i], Memc[rdnoise])
+ }
+ i = 1
+ if (ctor (Memc[gain], i, r) > 0) {
+ do i = 1, nimages {
+ Memr[nm+3*(i-1)+1] = r
+ Memr[nm+3*(i-1)] =
+ max ((Memr[nm+3*(i-1)] / r) ** 2, 1e4 / MAX_REAL)
+ }
+ } else {
+ do i = 1, nimages {
+ r = imgetr (in[i], Memc[gain])
+ Memr[nm+3*(i-1)+1] = r
+ Memr[nm+3*(i-1)] =
+ max ((Memr[nm+3*(i-1)] / r) ** 2, 1e4 / MAX_REAL)
+ }
+ }
+ i = 1
+ if (ctor (Memc[snoise], i, r) > 0) {
+ do i = 1, nimages
+ Memr[nm+3*(i-1)+2] = r
+ } else {
+ do i = 1, nimages {
+ r = imgetr (in[i], Memc[snoise])
+ Memr[nm+3*(i-1)+2] = r
+ }
+ }
+ if (!keepids) {
+ if (doscale1 || grow > 0)
+ keepids = true
+ else {
+ do i = 2, nimages {
+ if (Memr[nm+3*(i-1)] != Memr[nm] ||
+ Memr[nm+3*(i-1)+1] != Memr[nm+1] ||
+ Memr[nm+3*(i-1)+2] != Memr[nm+2]) {
+ keepids = true
+ break
+ }
+ }
+ }
+ }
+ if (reject == CRREJECT)
+ lsigma = MAX_REAL
+ case MINMAX:
+ mclip = false
+ if (grow > 0)
+ keepids = true
+ case PCLIP:
+ mclip = true
+ if (grow > 0)
+ keepids = true
+ case AVSIGCLIP, SIGCLIP:
+ if (doscale1 || grow > 0)
+ keepids = true
+ case NONE:
+ mclip = false
+ grow = 0
+ }
+
+ if (keepids) {
+ do i = 1, nimages
+ call salloc (id[i], npts, TY_INT)
+ }
+
+ while (impnlr (out[1], outdata, Meml[v1]) != EOF) {
+ call ic_gdatar (in, out, dbuf, d, id, n, m, lflag, offsets,
+ scales, zeros, nimages, npts, Meml[v2], Meml[v3])
+
+ switch (reject) {
+ case CCDCLIP, CRREJECT:
+ if (mclip)
+ call ic_mccdclipr (d, id, n, scales, zeros, Memr[nm],
+ nimages, npts, Memr[outdata])
+ else
+ call ic_accdclipr (d, id, n, scales, zeros, Memr[nm],
+ nimages, npts, Memr[outdata])
+ case MINMAX:
+ call ic_mmr (d, id, n, npts)
+ case PCLIP:
+ call ic_pclipr (d, id, n, nimages, npts, Memr[outdata])
+ case SIGCLIP:
+ if (mclip)
+ call ic_msigclipr (d, id, n, scales, zeros, nimages, npts,
+ Memr[outdata])
+ else
+ call ic_asigclipr (d, id, n, scales, zeros, nimages, npts,
+ Memr[outdata])
+ case AVSIGCLIP:
+ if (mclip)
+ call ic_mavsigclipr (d, id, n, scales, zeros, nimages,
+ npts, Memr[outdata])
+ else
+ call ic_aavsigclipr (d, id, n, scales, zeros, nimages,
+ npts, Memr[outdata])
+ }
+
+ if (grow > 0)
+ call ic_growr (d, id, n, nimages, npts, Memr[outdata])
+
+ if (docombine) {
+ switch (combine) {
+ case AVERAGE:
+ call ic_averager (d, id, n, wts, npts, Memr[outdata])
+ case MEDIAN:
+ call ic_medianr (d, n, npts, Memr[outdata])
+ }
+ }
+
+ if (out[2] != NULL) {
+ call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+ i = impnli (out[2], buf, Meml[v1])
+ call amovki (nimages, Memi[buf], npts)
+ call asubi (Memi[buf], n, Memi[buf], npts)
+ }
+
+ if (out[3] != NULL) {
+ call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+ i = impnlr (out[3], buf, Meml[v1])
+ call ic_sigmar (d, id, n, wts, npts, Memr[outdata],
+ Memr[buf])
+ }
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+ }
+
+ call sfree (sp)
+end
+
diff --git a/noao/imred/ccdred/src/generic/icpclip.x b/noao/imred/ccdred/src/generic/icpclip.x
new file mode 100644
index 00000000..da09bb75
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icpclip.x
@@ -0,0 +1,442 @@
+# 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_pclips (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 = Mems[d[n2-1]+j]
+ med = (med + Mems[d[n2]+j]) / 2.
+ } else
+ med = Mems[d[n2]+j]
+
+ if (n1 < max (MINCLIP, maxkeep+1)) {
+ if (combine == MEDIAN)
+ median[i] = med
+ next
+ }
+
+ # Define sigma for clipping
+ sigma = t * (Mems[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 - Mems[d[nl]+j]) / sigma
+ if (r < lsigma)
+ break
+ Memr[resid+nl] = r
+ }
+ for (nh=n1; nh>=1; nh=nh-1) {
+ r = (Mems[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 = Mems[d[n5-1]+j]
+ med = (med + Mems[d[n5]+j]) / 2.
+ } else
+ med = Mems[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) {
+ Mems[d[l]+j] = Mems[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) {
+ Mems[d[l]+j] = Mems[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
+
+# 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/imred/ccdred/src/generic/icsclip.x b/noao/imred/ccdred/src/generic/icsclip.x
new file mode 100644
index 00000000..d7ccfd84
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icsclip.x
@@ -0,0 +1,964 @@
+# 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_asigclips (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 = Mems[d[1]+k]
+ do j = 2, n1
+ sum = sum + Mems[d[j]+k]
+ average[i] = sum / n1
+ }
+ }
+ next
+ }
+
+ # Compute average with the high and low rejected.
+ low = Mems[d[1]+k]
+ high = Mems[d[2]+k]
+ if (low > high) {
+ d1 = low
+ low = high
+ high = d1
+ }
+ sum = 0.
+ do j = 3, n1 {
+ d1 = Mems[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 = Mems[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 = Mems[dp1]
+ r = (d1 - a) / (s * Memr[wp])
+ if (r < -lsigma || r > hsigma) {
+ Memr[resid+n1] = abs (r)
+ if (j < n1) {
+ dp2 = d[n1] + k
+ Mems[dp1] = Mems[dp2]
+ Mems[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 + (Mems[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 = Mems[dp1]
+ r = (d1 - a) / s
+ if (r < -lsigma || r > hsigma) {
+ Memr[resid+n1] = abs (r)
+ if (j < n1) {
+ dp2 = d[n1] + k
+ Mems[dp1] = Mems[dp2]
+ Mems[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 = Mems[dp1]
+ Mems[dp1] = Mems[dp2]
+ Mems[dp2] = d1
+ mp2 = m[l] + k
+ s = Memi[mp1]
+ Memi[mp1] = Memi[mp2]
+ Memi[mp2] = s
+ }
+ }
+ }
+ sum = sum + Mems[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 = Mems[dp1]
+ Mems[dp1] = Mems[dp2]
+ Mems[dp2] = d1
+ if (keepids) {
+ mp1 = m[j] + k
+ mp2 = m[l] + k
+ s = Memi[mp1]
+ Memi[mp1] = Memi[mp2]
+ Memi[mp2] = s
+ }
+ }
+ }
+ }
+ sum = sum + Mems[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_msigclips (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 = (Mems[d[n3-1]+k] + Mems[d[n3]+k]) / 2.
+ else
+ med = Mems[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 + ((Mems[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 - Mems[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 = (Mems[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 + (Mems[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 - Mems[d[nl]+k]) / s
+ if (r <= lsigma)
+ break
+ Memr[resid+nl] = r
+ n1 = n1 - 1
+ }
+ for (; nh >= nl; nh = nh - 1) {
+ r = (Mems[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) {
+ Mems[d[l]+k] = Mems[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) {
+ Mems[d[l]+k] = Mems[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
+
+# 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/imred/ccdred/src/generic/icsigma.x b/noao/imred/ccdred/src/generic/icsigma.x
new file mode 100644
index 00000000..bc0d9788
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icsigma.x
@@ -0,0 +1,205 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include "../icombine.h"
+
+
+# IC_SIGMA -- Compute the sigma image line.
+# The estimated sigma includes a correction for the finite population.
+# Weights are used if desired.
+
+procedure ic_sigmas (d, m, n, wts, npts, average, sigma)
+
+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
+real sigma[npts] # Sigma line (returned)
+
+int i, j, k, n1
+real wt, sigcor, sumwt
+real a, sum
+
+include "../icombine.com"
+
+begin
+ if (dflag == D_ALL) {
+ n1 = n[1]
+ if (dowts) {
+ if (n1 > 1)
+ sigcor = real (n1) / real (n1 - 1)
+ else
+ sigcor = 1.
+ do i = 1, npts {
+ k = i - 1
+ a = average[i]
+ wt = wts[Memi[m[1]+k]]
+ sum = (Mems[d[1]+k] - a) ** 2 * wt
+ do j = 2, n1 {
+ wt = wts[Memi[m[j]+k]]
+ sum = sum + (Mems[d[j]+k] - a) ** 2 * wt
+ }
+ sigma[i] = sqrt (sum * sigcor)
+ }
+ } else {
+ if (n1 > 1)
+ sigcor = 1. / real (n1 - 1)
+ else
+ sigcor = 1.
+ do i = 1, npts {
+ k = i - 1
+ a = average[i]
+ sum = (Mems[d[1]+k] - a) ** 2
+ do j = 2, n1
+ sum = sum + (Mems[d[j]+k] - a) ** 2
+ sigma[i] = sqrt (sum * sigcor)
+ }
+ }
+ } else if (dflag == D_NONE) {
+ do i = 1, npts
+ sigma[i] = blank
+ } else {
+ if (dowts) {
+ do i = 1, npts {
+ n1 = n[i]
+ if (n1 > 0) {
+ k = i - 1
+ if (n1 > 1)
+ sigcor = real (n1) / real (n1 -1)
+ else
+ sigcor = 1
+ a = average[i]
+ wt = wts[Memi[m[1]+k]]
+ sum = (Mems[d[1]+k] - a) ** 2 * wt
+ sumwt = wt
+ do j = 2, n1 {
+ wt = wts[Memi[m[j]+k]]
+ sum = sum + (Mems[d[j]+k] - a) ** 2 * wt
+ sumwt = sumwt + wt
+ }
+ sigma[i] = sqrt (sum / sumwt * sigcor)
+ } else
+ sigma[i] = blank
+ }
+ } else {
+ do i = 1, npts {
+ n1 = n[i]
+ if (n1 > 0) {
+ k = i - 1
+ if (n1 > 1)
+ sigcor = 1. / real (n1 - 1)
+ else
+ sigcor = 1.
+ a = average[i]
+ sum = (Mems[d[1]+k] - a) ** 2
+ do j = 2, n1
+ sum = sum + (Mems[d[j]+k] - a) ** 2
+ sigma[i] = sqrt (sum * sigcor)
+ } else
+ sigma[i] = blank
+ }
+ }
+ }
+end
+
+# IC_SIGMA -- Compute the sigma image line.
+# The estimated sigma includes a correction for the finite population.
+# Weights are used if desired.
+
+procedure ic_sigmar (d, m, n, wts, npts, average, sigma)
+
+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
+real sigma[npts] # Sigma line (returned)
+
+int i, j, k, n1
+real wt, sigcor, sumwt
+real a, sum
+
+include "../icombine.com"
+
+begin
+ if (dflag == D_ALL) {
+ n1 = n[1]
+ if (dowts) {
+ if (n1 > 1)
+ sigcor = real (n1) / real (n1 - 1)
+ else
+ sigcor = 1.
+ do i = 1, npts {
+ k = i - 1
+ a = average[i]
+ wt = wts[Memi[m[1]+k]]
+ sum = (Memr[d[1]+k] - a) ** 2 * wt
+ do j = 2, n1 {
+ wt = wts[Memi[m[j]+k]]
+ sum = sum + (Memr[d[j]+k] - a) ** 2 * wt
+ }
+ sigma[i] = sqrt (sum * sigcor)
+ }
+ } else {
+ if (n1 > 1)
+ sigcor = 1. / real (n1 - 1)
+ else
+ sigcor = 1.
+ do i = 1, npts {
+ k = i - 1
+ a = average[i]
+ sum = (Memr[d[1]+k] - a) ** 2
+ do j = 2, n1
+ sum = sum + (Memr[d[j]+k] - a) ** 2
+ sigma[i] = sqrt (sum * sigcor)
+ }
+ }
+ } else if (dflag == D_NONE) {
+ do i = 1, npts
+ sigma[i] = blank
+ } else {
+ if (dowts) {
+ do i = 1, npts {
+ n1 = n[i]
+ if (n1 > 0) {
+ k = i - 1
+ if (n1 > 1)
+ sigcor = real (n1) / real (n1 -1)
+ else
+ sigcor = 1
+ a = average[i]
+ wt = wts[Memi[m[1]+k]]
+ sum = (Memr[d[1]+k] - a) ** 2 * wt
+ sumwt = wt
+ do j = 2, n1 {
+ wt = wts[Memi[m[j]+k]]
+ sum = sum + (Memr[d[j]+k] - a) ** 2 * wt
+ sumwt = sumwt + wt
+ }
+ sigma[i] = sqrt (sum / sumwt * sigcor)
+ } else
+ sigma[i] = blank
+ }
+ } else {
+ do i = 1, npts {
+ n1 = n[i]
+ if (n1 > 0) {
+ k = i - 1
+ if (n1 > 1)
+ sigcor = 1. / real (n1 - 1)
+ else
+ sigcor = 1.
+ a = average[i]
+ sum = (Memr[d[1]+k] - a) ** 2
+ do j = 2, n1
+ sum = sum + (Memr[d[j]+k] - a) ** 2
+ sigma[i] = sqrt (sum * sigcor)
+ } else
+ sigma[i] = blank
+ }
+ }
+ }
+end
diff --git a/noao/imred/ccdred/src/generic/icsort.x b/noao/imred/ccdred/src/generic/icsort.x
new file mode 100644
index 00000000..a39b68e2
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icsort.x
@@ -0,0 +1,550 @@
+# 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_sorts (a, b, nvecs, npts)
+
+pointer a[ARB] # pointer to input vectors
+short b[ARB] # work array
+int nvecs[npts] # number of vectors
+int npts # number of points in vectors
+
+short 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] = Mems[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] = Mems[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
+ Mems[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_2sorts (a, b, c, d, nvecs, npts)
+
+pointer a[ARB] # pointer to input vectors
+short 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
+
+short 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] = Mems[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]
+ d[3] = d[1]
+ d[1] = itemp
+ }
+ } 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 {
+ Mems[a[i]+l] = b[i]
+ Memi[c[i]+l] = d[i]
+ }
+ }
+end
+
+# 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]
+ d[3] = d[1]
+ d[1] = itemp
+ }
+ } 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/imred/ccdred/src/generic/icstat.x b/noao/imred/ccdred/src/generic/icstat.x
new file mode 100644
index 00000000..41512ccb
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/icstat.x
@@ -0,0 +1,444 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include "../icombine.h"
+
+define NMAX 10000 # Maximum number of pixels to sample
+
+
+# IC_STAT -- Compute image statistics within specified section.
+# The image section is relative to a reference image which may be
+# different than the input image and may have an offset. Only a
+# subsample of pixels is used. Masked and thresholded pixels are
+# ignored. Only the desired statistics are computed to increase
+# efficiency.
+
+procedure ic_stats (im, imref, section, offsets, image, nimages,
+ domode, domedian, domean, mode, median, mean)
+
+pointer im # Data image
+pointer imref # Reference image for image section
+char section[ARB] # Image section
+int offsets[nimages,ARB] # Image section offset from data to reference
+int image # Image index (for mask I/O)
+int nimages # Number of images in offsets.
+bool domode, domedian, domean # Statistics to compute
+real mode, median, mean # Statistics
+
+int i, j, ndim, n, nv
+real a
+pointer sp, v1, v2, dv, va, vb
+pointer data, mask, dp, lp, mp, imgnls()
+short ic_modes()
+real asums()
+
+
+include "../icombine.com"
+
+begin
+ call smark (sp)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (v2, IM_MAXDIM, TY_LONG)
+ call salloc (dv, IM_MAXDIM, TY_LONG)
+ call salloc (va, IM_MAXDIM, TY_LONG)
+ call salloc (vb, IM_MAXDIM, TY_LONG)
+
+ # Determine the image section parameters. This must be in terms of
+ # the data image pixel coordinates though the section may be specified
+ # in terms of the reference image coordinates. Limit the number of
+ # pixels in each dimension to a maximum.
+
+ ndim = IM_NDIM(im)
+ if (project)
+ ndim = ndim - 1
+ call amovki (1, Memi[v1], IM_MAXDIM)
+ call amovki (1, Memi[va], IM_MAXDIM)
+ call amovki (1, Memi[dv], IM_MAXDIM)
+ call amovi (IM_LEN(imref,1), Memi[vb], ndim)
+ call ic_section (section, Memi[va], Memi[vb], Memi[dv], ndim)
+ if (im != imref)
+ do i = 1, ndim {
+ Memi[va+i-1] = Memi[va+i-1] - offsets[image,i]
+ Memi[vb+i-1] = Memi[vb+i-1] - offsets[image,i]
+ }
+
+ do j = 1, 10 {
+ n = 1
+ do i = 0, ndim-1 {
+ Memi[v1+i] = max (1, min (Memi[va+i], Memi[vb+i]))
+ Memi[v2+i] = min (IM_LEN(im,i+1), max (Memi[va+i], Memi[vb+i]))
+ Memi[dv+i] = j
+ nv = max (1, (Memi[v2+i] - Memi[v1+i]) / Memi[dv+i] + 1)
+ Memi[v2+i] = Memi[v1+i] + (nv - 1) * Memi[dv+i]
+ n = n * nv
+ }
+ if (n < NMAX)
+ break
+ }
+
+ call amovl (Memi[v1], Memi[va], IM_MAXDIM)
+ Memi[va] = 1
+ if (project)
+ Memi[va+ndim] = image
+ call amovl (Memi[va], Memi[vb], IM_MAXDIM)
+
+ # Accumulate the pixel values within the section. Masked pixels and
+ # thresholded pixels are ignored.
+
+ call salloc (data, n, TY_SHORT)
+ dp = data
+ while (imgnls (im, lp, Memi[vb]) != EOF) {
+ call ic_mget1 (im, image, offsets[image,1], Memi[va], mask)
+ lp = lp + Memi[v1] - 1
+ if (dflag == D_ALL) {
+ if (dothresh) {
+ do i = Memi[v1], Memi[v2], Memi[dv] {
+ a = Mems[lp]
+ if (a >= lthresh && a <= hthresh) {
+ Mems[dp] = a
+ dp = dp + 1
+ }
+ lp = lp + Memi[dv]
+ }
+ } else {
+ do i = Memi[v1], Memi[v2], Memi[dv] {
+ Mems[dp] = Mems[lp]
+ dp = dp + 1
+ lp = lp + Memi[dv]
+ }
+ }
+ } else if (dflag == D_MIX) {
+ mp = mask + Memi[v1] - 1
+ if (dothresh) {
+ do i = Memi[v1], Memi[v2], Memi[dv] {
+ if (Memi[mp] == 0) {
+ a = Mems[lp]
+ if (a >= lthresh && a <= hthresh) {
+ Mems[dp] = a
+ dp = dp + 1
+ }
+ }
+ mp = mp + Memi[dv]
+ lp = lp + Memi[dv]
+ }
+ } else {
+ do i = Memi[v1], Memi[v2], Memi[dv] {
+ if (Memi[mp] == 0) {
+ Mems[dp] = Mems[lp]
+ dp = dp + 1
+ }
+ mp = mp + Memi[dv]
+ lp = lp + Memi[dv]
+ }
+ }
+ }
+ for (i=2; i<=ndim; i=i+1) {
+ Memi[va+i-1] = Memi[va+i-1] + Memi[dv+i-1]
+ if (Memi[va+i-1] <= Memi[v2+i-1])
+ break
+ Memi[va+i-1] = Memi[v1+i-1]
+ }
+ if (i > ndim)
+ break
+ call amovl (Memi[va], Memi[vb], IM_MAXDIM)
+ }
+
+ n = dp - data
+ if (n < 1) {
+ call sfree (sp)
+ call error (1, "Image section contains no pixels")
+ }
+
+ # Compute only statistics needed.
+ if (domode || domedian) {
+ call asrts (Mems[data], Mems[data], n)
+ mode = ic_modes (Mems[data], n)
+ median = Mems[data+n/2-1]
+ }
+ if (domean)
+ mean = asums (Mems[data], n) / n
+
+ call sfree (sp)
+end
+
+
+define NMIN 10 # Minimum number of pixels for mode calculation
+define ZRANGE 0.8 # Fraction of pixels about median to use
+define ZSTEP 0.01 # Step size for search for mode
+define ZBIN 0.1 # Bin size for mode.
+
+# IC_MODE -- Compute mode of an array. The mode is found by binning
+# with a bin size based on the data range over a fraction of the
+# pixels about the median and a bin step which may be smaller than the
+# bin size. If there are too few points the median is returned.
+# The input array must be sorted.
+
+short procedure ic_modes (a, n)
+
+short a[n] # Data array
+int n # Number of points
+
+int i, j, k, nmax
+real z1, z2, zstep, zbin
+short mode
+bool fp_equalr()
+
+begin
+ if (n < NMIN)
+ return (a[n/2])
+
+ # Compute the mode. The array must be sorted. Consider a
+ # range of values about the median point. Use a bin size which
+ # is ZBIN of the range. Step the bin limits in ZSTEP fraction of
+ # the bin size.
+
+ i = 1 + n * (1. - ZRANGE) / 2.
+ j = 1 + n * (1. + ZRANGE) / 2.
+ z1 = a[i]
+ z2 = a[j]
+ if (fp_equalr (z1, z2)) {
+ mode = z1
+ return (mode)
+ }
+
+ zstep = ZSTEP * (z2 - z1)
+ zbin = ZBIN * (z2 - z1)
+ zstep = max (1., zstep)
+ zbin = max (1., zbin)
+
+ z1 = z1 - zstep
+ k = i
+ nmax = 0
+ repeat {
+ z1 = z1 + zstep
+ z2 = z1 + zbin
+ for (; i < j && a[i] < z1; i=i+1)
+ ;
+ for (; k < j && a[k] < z2; k=k+1)
+ ;
+ if (k - i > nmax) {
+ nmax = k - i
+ mode = a[(i+k)/2]
+ }
+ } until (k >= j)
+
+ return (mode)
+end
+
+# IC_STAT -- Compute image statistics within specified section.
+# The image section is relative to a reference image which may be
+# different than the input image and may have an offset. Only a
+# subsample of pixels is used. Masked and thresholded pixels are
+# ignored. Only the desired statistics are computed to increase
+# efficiency.
+
+procedure ic_statr (im, imref, section, offsets, image, nimages,
+ domode, domedian, domean, mode, median, mean)
+
+pointer im # Data image
+pointer imref # Reference image for image section
+char section[ARB] # Image section
+int offsets[nimages,ARB] # Image section offset from data to reference
+int image # Image index (for mask I/O)
+int nimages # Number of images in offsets.
+bool domode, domedian, domean # Statistics to compute
+real mode, median, mean # Statistics
+
+int i, j, ndim, n, nv
+real a
+pointer sp, v1, v2, dv, va, vb
+pointer data, mask, dp, lp, mp, imgnlr()
+real ic_moder()
+real asumr()
+
+
+include "../icombine.com"
+
+begin
+ call smark (sp)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (v2, IM_MAXDIM, TY_LONG)
+ call salloc (dv, IM_MAXDIM, TY_LONG)
+ call salloc (va, IM_MAXDIM, TY_LONG)
+ call salloc (vb, IM_MAXDIM, TY_LONG)
+
+ # Determine the image section parameters. This must be in terms of
+ # the data image pixel coordinates though the section may be specified
+ # in terms of the reference image coordinates. Limit the number of
+ # pixels in each dimension to a maximum.
+
+ ndim = IM_NDIM(im)
+ if (project)
+ ndim = ndim - 1
+ call amovki (1, Memi[v1], IM_MAXDIM)
+ call amovki (1, Memi[va], IM_MAXDIM)
+ call amovki (1, Memi[dv], IM_MAXDIM)
+ call amovi (IM_LEN(imref,1), Memi[vb], ndim)
+ call ic_section (section, Memi[va], Memi[vb], Memi[dv], ndim)
+ if (im != imref)
+ do i = 1, ndim {
+ Memi[va+i-1] = Memi[va+i-1] - offsets[image,i]
+ Memi[vb+i-1] = Memi[vb+i-1] - offsets[image,i]
+ }
+
+ do j = 1, 10 {
+ n = 1
+ do i = 0, ndim-1 {
+ Memi[v1+i] = max (1, min (Memi[va+i], Memi[vb+i]))
+ Memi[v2+i] = min (IM_LEN(im,i+1), max (Memi[va+i], Memi[vb+i]))
+ Memi[dv+i] = j
+ nv = max (1, (Memi[v2+i] - Memi[v1+i]) / Memi[dv+i] + 1)
+ Memi[v2+i] = Memi[v1+i] + (nv - 1) * Memi[dv+i]
+ n = n * nv
+ }
+ if (n < NMAX)
+ break
+ }
+
+ call amovl (Memi[v1], Memi[va], IM_MAXDIM)
+ Memi[va] = 1
+ if (project)
+ Memi[va+ndim] = image
+ call amovl (Memi[va], Memi[vb], IM_MAXDIM)
+
+ # Accumulate the pixel values within the section. Masked pixels and
+ # thresholded pixels are ignored.
+
+ call salloc (data, n, TY_REAL)
+ dp = data
+ while (imgnlr (im, lp, Memi[vb]) != EOF) {
+ call ic_mget1 (im, image, offsets[image,1], Memi[va], mask)
+ lp = lp + Memi[v1] - 1
+ if (dflag == D_ALL) {
+ if (dothresh) {
+ do i = Memi[v1], Memi[v2], Memi[dv] {
+ a = Memr[lp]
+ if (a >= lthresh && a <= hthresh) {
+ Memr[dp] = a
+ dp = dp + 1
+ }
+ lp = lp + Memi[dv]
+ }
+ } else {
+ do i = Memi[v1], Memi[v2], Memi[dv] {
+ Memr[dp] = Memr[lp]
+ dp = dp + 1
+ lp = lp + Memi[dv]
+ }
+ }
+ } else if (dflag == D_MIX) {
+ mp = mask + Memi[v1] - 1
+ if (dothresh) {
+ do i = Memi[v1], Memi[v2], Memi[dv] {
+ if (Memi[mp] == 0) {
+ a = Memr[lp]
+ if (a >= lthresh && a <= hthresh) {
+ Memr[dp] = a
+ dp = dp + 1
+ }
+ }
+ mp = mp + Memi[dv]
+ lp = lp + Memi[dv]
+ }
+ } else {
+ do i = Memi[v1], Memi[v2], Memi[dv] {
+ if (Memi[mp] == 0) {
+ Memr[dp] = Memr[lp]
+ dp = dp + 1
+ }
+ mp = mp + Memi[dv]
+ lp = lp + Memi[dv]
+ }
+ }
+ }
+ for (i=2; i<=ndim; i=i+1) {
+ Memi[va+i-1] = Memi[va+i-1] + Memi[dv+i-1]
+ if (Memi[va+i-1] <= Memi[v2+i-1])
+ break
+ Memi[va+i-1] = Memi[v1+i-1]
+ }
+ if (i > ndim)
+ break
+ call amovl (Memi[va], Memi[vb], IM_MAXDIM)
+ }
+
+ n = dp - data
+ if (n < 1) {
+ call sfree (sp)
+ call error (1, "Image section contains no pixels")
+ }
+
+ # Compute only statistics needed.
+ if (domode || domedian) {
+ call asrtr (Memr[data], Memr[data], n)
+ mode = ic_moder (Memr[data], n)
+ median = Memr[data+n/2-1]
+ }
+ if (domean)
+ mean = asumr (Memr[data], n) / n
+
+ call sfree (sp)
+end
+
+
+define NMIN 10 # Minimum number of pixels for mode calculation
+define ZRANGE 0.8 # Fraction of pixels about median to use
+define ZSTEP 0.01 # Step size for search for mode
+define ZBIN 0.1 # Bin size for mode.
+
+# IC_MODE -- Compute mode of an array. The mode is found by binning
+# with a bin size based on the data range over a fraction of the
+# pixels about the median and a bin step which may be smaller than the
+# bin size. If there are too few points the median is returned.
+# The input array must be sorted.
+
+real procedure ic_moder (a, n)
+
+real a[n] # Data array
+int n # Number of points
+
+int i, j, k, nmax
+real z1, z2, zstep, zbin
+real mode
+bool fp_equalr()
+
+begin
+ if (n < NMIN)
+ return (a[n/2])
+
+ # Compute the mode. The array must be sorted. Consider a
+ # range of values about the median point. Use a bin size which
+ # is ZBIN of the range. Step the bin limits in ZSTEP fraction of
+ # the bin size.
+
+ i = 1 + n * (1. - ZRANGE) / 2.
+ j = 1 + n * (1. + ZRANGE) / 2.
+ z1 = a[i]
+ z2 = a[j]
+ if (fp_equalr (z1, z2)) {
+ mode = z1
+ return (mode)
+ }
+
+ zstep = ZSTEP * (z2 - z1)
+ zbin = ZBIN * (z2 - z1)
+
+ z1 = z1 - zstep
+ k = i
+ nmax = 0
+ repeat {
+ z1 = z1 + zstep
+ z2 = z1 + zbin
+ for (; i < j && a[i] < z1; i=i+1)
+ ;
+ for (; k < j && a[k] < z2; k=k+1)
+ ;
+ if (k - i > nmax) {
+ nmax = k - i
+ mode = a[(i+k)/2]
+ }
+ } until (k >= j)
+
+ return (mode)
+end
+
diff --git a/noao/imred/ccdred/src/generic/mkpkg b/noao/imred/ccdred/src/generic/mkpkg
new file mode 100644
index 00000000..3d841680
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/mkpkg
@@ -0,0 +1,11 @@
+# Make CCDRED Package.
+
+$checkout libpkg.a ../..
+$update libpkg.a
+$checkin libpkg.a ../..
+$exit
+
+libpkg.a:
+ cor.x ccdred.h
+ proc.x ccdred.h <imhdr.h>
+ ;
diff --git a/noao/imred/ccdred/src/generic/proc.x b/noao/imred/ccdred/src/generic/proc.x
new file mode 100644
index 00000000..242da9c9
--- /dev/null
+++ b/noao/imred/ccdred/src/generic/proc.x
@@ -0,0 +1,735 @@
+include <imhdr.h>
+include "ccdred.h"
+
+
+.help proc Feb87 noao.imred.ccdred
+.nf ----------------------------------------------------------------------------
+proc -- Process CCD images
+
+These are the main CCD reduction procedures. There is one for each
+readout axis (lines or columns) and one for short and real image data.
+They apply corrections for bad pixels, overscan levels, zero levels,
+dark counts, flat field response, illumination response, and fringe
+effects. The image is also trimmed if it was mapped with an image
+section. The mean value for the output image is computed when the flat
+field or illumination image is processed to form the scale factor for
+these calibrations in order to avoid reading through these image a
+second time.
+
+The processing information and parameters are specified in the CCD
+structure. The processing operations to be performed are specified by
+the correction array CORS in the ccd structure. There is one array
+element for each operation with indices defined symbolically by macro
+definitions (see ccdred.h); i.e. FLATCOR. The value of the array
+element is an integer bit field in which the bit set is the same as the
+array index; i.e element 3 will have the third bit set for an operation
+with array value 2**(3-1)=4. If an operation is not to be performed
+the bit is not set and the array element has the numeric value zero.
+Note that the addition of several correction elements gives a unique
+bit field describing a combination of operations. For efficiency the
+most common combinations are implemented as separate units.
+
+The CCD structure also contains the correction or calibration data
+consisting either pointers to data, IMIO pointers for the calibration
+images, and scale factors.
+
+The processing is performed line-by-line. The procedure CORINPUT is
+called to get an input line. This procedure trims and fixes bad pixels by
+interpolation. The output line and lines from the various calibration
+images are read. The image vectors as well as the overscan vector and
+the scale factors are passed to the procedure COR (which also
+dereferences the pointer data into simple arrays and variables). That
+procedure does the actual corrections apart from bad pixel
+corrections.
+
+The final optional step is to add each corrected output line to form a
+mean. This adds efficiency since the operation is done only if desired
+and the output image data is already in memory so there is no I/O
+penalty.
+
+SEE ALSO
+ ccdred.h, cor, fixpix, setfixpix, setoverscan, settrim,
+ setzero, setdark, setflat, setillum, setfringe
+.endhelp ----------------------------------------------------------------------
+
+
+
+# PROC1 -- Process CCD images with readout axis 1 (lines).
+
+procedure proc1s (ccd)
+
+pointer ccd # CCD structure
+
+int line, ncols, nlines, findmean, rep
+int overscan_type, overscan_c1, noverscan
+real overscan, darkscale, flatscale, illumscale, frgscale, mean
+short minrep
+pointer in, out, zeroim, darkim, flatim, illumim, fringeim, overscan_vec
+pointer inbuf, outbuf, zerobuf, darkbuf, flatbuf, illumbuf, fringebuf
+
+real asums()
+real find_overscans()
+pointer imgl2s(), impl2s(), ccd_gls(), xt_fpss()
+
+begin
+ # Initialize. If the correction image is 1D then just get the
+ # data once.
+
+ in = IN_IM(ccd)
+ out = OUT_IM(ccd)
+ ncols = OUT_C2(ccd) - OUT_C1(ccd) + 1
+ nlines = OUT_L2(ccd) - OUT_L1(ccd) + 1
+
+ findmean = CORS(ccd, FINDMEAN)
+ if (findmean == YES)
+ mean = 0.
+ rep = CORS(ccd, MINREP)
+ if (rep == YES)
+ minrep = MINREPLACE(ccd)
+
+ if (CORS(ccd, OVERSCAN) == 0)
+ overscan_type = 0
+ else {
+ overscan_type = OVERSCAN_TYPE(ccd)
+ overscan_vec = OVERSCAN_VEC(ccd)
+ overscan_c1 = BIAS_C1(ccd) - 1
+ noverscan = BIAS_C2(ccd) - overscan_c1
+ }
+
+ if (CORS(ccd, ZEROCOR) == 0) {
+ zeroim = NULL
+ zerobuf = 1
+ } else if (IM_LEN(ZERO_IM(ccd),2) == 1) {
+ zeroim = NULL
+ zerobuf = ccd_gls (ZERO_IM(ccd), ZERO_C1(ccd), ZERO_C2(ccd), 1)
+ } else
+ zeroim = ZERO_IM(ccd)
+
+ if (CORS(ccd, DARKCOR) == 0) {
+ darkim = NULL
+ darkbuf = 1
+ } else if (IM_LEN(DARK_IM(ccd),2) == 1) {
+ darkim = NULL
+ darkbuf = ccd_gls (DARK_IM(ccd), DARK_C1(ccd), DARK_C2(ccd), 1)
+ darkscale = FLATSCALE(ccd)
+ } else {
+ darkim = DARK_IM(ccd)
+ darkscale = DARKSCALE(ccd)
+ }
+
+ if (CORS(ccd, FLATCOR) == 0) {
+ flatim = NULL
+ flatbuf = 1
+ } else if (IM_LEN(FLAT_IM(ccd),2) == 1) {
+ flatim = NULL
+ flatbuf = ccd_gls (FLAT_IM(ccd), FLAT_C1(ccd), FLAT_C2(ccd), 1)
+ flatscale = FLATSCALE(ccd)
+ } else {
+ flatim = FLAT_IM(ccd)
+ flatscale = FLATSCALE(ccd)
+ }
+
+ if (CORS(ccd, ILLUMCOR) == 0) {
+ illumim = NULL
+ illumbuf = 1
+ } else {
+ illumim = ILLUM_IM(ccd)
+ illumscale = ILLUMSCALE(ccd)
+ }
+
+ if (CORS(ccd, FRINGECOR) == 0) {
+ fringeim = NULL
+ fringebuf = 1
+ } else {
+ fringeim = FRINGE_IM(ccd)
+ frgscale = FRINGESCALE(ccd)
+ }
+
+ # For each line read lines from the input. Procedure XT_FPS replaces
+ # bad pixels by interpolation. The trimmed region is copied to the
+ # output. Get lines from the output image and from the zero level,
+ # dark count, flat field, illumination, and fringe images. Call COR1
+ # to do the actual pixel corrections. Finally, add the output pixels
+ # to a sum for computing the mean. We must copy data outside of the
+ # output data section.
+
+ do line = 2 - OUT_L1(ccd), 0
+ call amovs (
+ Mems[imgl2s(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+ Mems[impl2s(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+ do line = 1, nlines {
+ outbuf = impl2s (out, OUT_L1(ccd)+line-1)
+
+ inbuf = xt_fpss (MASK_FP(ccd), in, IN_L1(ccd)+line-1, IN_C1(ccd),
+ IN_C2(ccd), IN_L1(ccd), IN_L2(ccd), NULL)
+ call amovs (Mems[inbuf+IN_C1(ccd)-OUT_C1(ccd)], Mems[outbuf],
+ IM_LEN(out,1))
+
+ outbuf = outbuf + OUT_C1(ccd) - 1
+ if (overscan_type != 0) {
+ if (overscan_type < OVERSCAN_FIT)
+ overscan = find_overscans (Mems[inbuf+overscan_c1],
+ noverscan, overscan_type)
+ else
+ overscan = Memr[overscan_vec+line-1]
+ }
+ if (zeroim != NULL)
+ zerobuf = ccd_gls (zeroim, ZERO_C1(ccd), ZERO_C2(ccd),
+ ZERO_L1(ccd)+line-1)
+ if (darkim != NULL)
+ darkbuf = ccd_gls (darkim, DARK_C1(ccd), DARK_C2(ccd),
+ DARK_L1(ccd)+line-1)
+ if (flatim != NULL)
+ flatbuf = ccd_gls (flatim, FLAT_C1(ccd), FLAT_C2(ccd),
+ FLAT_L1(ccd)+line-1)
+ if (illumim != NULL)
+ illumbuf = ccd_gls (illumim, ILLUM_C1(ccd), ILLUM_C2(ccd),
+ ILLUM_L1(ccd)+line-1)
+ if (fringeim != NULL)
+ fringebuf = ccd_gls (fringeim, FRINGE_C1(ccd), FRINGE_C2(ccd),
+ FRINGE_L1(ccd)+line-1)
+
+ call cor1s (CORS(ccd,1), Mems[outbuf],
+ overscan, Mems[zerobuf], Mems[darkbuf],
+ Mems[flatbuf], Mems[illumbuf], Mems[fringebuf], ncols,
+ darkscale, flatscale, illumscale, frgscale)
+
+ if (rep == YES)
+ call amaxks (Mems[outbuf], minrep, Mems[outbuf], ncols)
+ if (findmean == YES)
+ mean = mean + asums (Mems[outbuf], ncols)
+ }
+
+ do line = nlines+1, IM_LEN(out,2)-OUT_L1(ccd)+1
+ call amovs (
+ Mems[imgl2s(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+ Mems[impl2s(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+ # Compute the mean from the sum of the output pixels.
+ if (findmean == YES)
+ MEAN(ccd) = mean / ncols / nlines
+end
+
+
+# PROC2 -- Process CCD images with readout axis 2 (columns).
+
+procedure proc2s (ccd)
+
+pointer ccd # CCD structure
+
+int line, ncols, nlines, findmean, rep
+real darkscale, flatscale, illumscale, frgscale, mean
+short minrep
+pointer in, out, zeroim, darkim, flatim, illumim, fringeim, overscan_vec
+pointer inbuf, outbuf, zerobuf, darkbuf, flatbuf, illumbuf, fringebuf
+
+real asums()
+pointer imgl2s(), impl2s(), imgs2s(), ccd_gls(), xt_fpss()
+
+begin
+ # Initialize. If the correction image is 1D then just get the
+ # data once.
+
+ in = IN_IM(ccd)
+ out = OUT_IM(ccd)
+ ncols = OUT_C2(ccd) - OUT_C1(ccd) + 1
+ nlines = OUT_L2(ccd) - OUT_L1(ccd) + 1
+
+ findmean = CORS(ccd, FINDMEAN)
+ if (findmean == YES)
+ mean = 0.
+ rep = CORS(ccd, MINREP)
+ if (rep == YES)
+ minrep = MINREPLACE(ccd)
+
+ overscan_vec = OVERSCAN_VEC(ccd)
+
+ if (CORS(ccd, ZEROCOR) == 0) {
+ zeroim = NULL
+ zerobuf = 1
+ } else if (IM_LEN(ZERO_IM(ccd),1) == 1) {
+ zeroim = NULL
+ zerobuf = imgs2s (ZERO_IM(ccd), 1, 1, ZERO_L1(ccd), ZERO_L2(ccd))
+ } else
+ zeroim = ZERO_IM(ccd)
+
+ if (CORS(ccd, DARKCOR) == 0) {
+ darkim = NULL
+ darkbuf = 1
+ } else if (IM_LEN(DARK_IM(ccd),1) == 1) {
+ darkim = NULL
+ darkbuf = imgs2s (DARK_IM(ccd), 1, 1, DARK_L1(ccd), DARK_L2(ccd))
+ darkscale = DARKSCALE(ccd)
+ } else {
+ darkim = DARK_IM(ccd)
+ darkscale = DARKSCALE(ccd)
+ }
+
+ if (CORS(ccd, FLATCOR) == 0) {
+ flatim = NULL
+ flatbuf = 1
+ } else if (IM_LEN(FLAT_IM(ccd),1) == 1) {
+ flatim = NULL
+ flatbuf = imgs2s (FLAT_IM(ccd), 1, 1, FLAT_L1(ccd), FLAT_L2(ccd))
+ flatscale = FLATSCALE(ccd)
+ } else {
+ flatim = FLAT_IM(ccd)
+ flatscale = FLATSCALE(ccd)
+ }
+
+ if (CORS(ccd, ILLUMCOR) == 0) {
+ illumim = NULL
+ illumbuf = 1
+ } else {
+ illumim = ILLUM_IM(ccd)
+ illumscale = ILLUMSCALE(ccd)
+ }
+
+ if (CORS(ccd, FRINGECOR) == 0) {
+ fringeim = NULL
+ fringebuf = 1
+ } else {
+ fringeim = FRINGE_IM(ccd)
+ frgscale = FRINGESCALE(ccd)
+ }
+
+ # For each line read lines from the input. Procedure CORINPUT
+ # replaces bad pixels by interpolation and applies a trim to the
+ # input. Get lines from the output image and from the zero level,
+ # dark count, flat field, illumination, and fringe images.
+ # Call COR2 to do the actual pixel corrections. Finally, add the
+ # output pixels to a sum for computing the mean.
+ # We must copy data outside of the output data section.
+
+ do line = 2 - OUT_L1(ccd), 0
+ call amovs (
+ Mems[imgl2s(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+ Mems[impl2s(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+ do line = 1, nlines {
+ outbuf = impl2s (out, OUT_L1(ccd)+line-1)
+
+ inbuf = xt_fpss (MASK_FP(ccd), in, IN_L1(ccd)+line-1, IN_C1(ccd),
+ IN_C2(ccd), IN_L1(ccd), IN_L2(ccd), NULL)
+ call amovs (Mems[inbuf+IN_C1(ccd)-OUT_C1(ccd)], Mems[outbuf],
+ IM_LEN(out,1))
+
+ outbuf = outbuf + OUT_C1(ccd) - 1
+ if (zeroim != NULL)
+ zerobuf = ccd_gls (zeroim, ZERO_C1(ccd), ZERO_C2(ccd),
+ ZERO_L1(ccd)+line-1)
+ if (darkim != NULL)
+ darkbuf = ccd_gls (darkim, DARK_C1(ccd), DARK_C2(ccd),
+ DARK_L1(ccd)+line-1)
+ if (flatim != NULL)
+ flatbuf = ccd_gls (flatim, FLAT_C1(ccd), FLAT_C2(ccd),
+ FLAT_L1(ccd)+line-1)
+ if (illumim != NULL)
+ illumbuf = ccd_gls (illumim, ILLUM_C1(ccd), ILLUM_C2(ccd),
+ ILLUM_L1(ccd)+line-1)
+ if (fringeim != NULL)
+ fringebuf = ccd_gls (fringeim, FRINGE_C1(ccd), FRINGE_C2(ccd),
+ FRINGE_L1(ccd)+line-1)
+
+ call cor2s (line, CORS(ccd,1), Mems[outbuf],
+ Memr[overscan_vec], Mems[zerobuf], Mems[darkbuf],
+ Mems[flatbuf], Mems[illumbuf], Mems[fringebuf], ncols,
+ zeroim, flatim, darkscale, flatscale, illumscale, frgscale)
+
+ if (rep == YES)
+ call amaxks (Mems[outbuf], minrep, Mems[outbuf], ncols)
+ if (findmean == YES)
+ mean = mean + asums (Mems[outbuf], ncols)
+ }
+
+ do line = nlines+1, IM_LEN(out,2)-OUT_L1(ccd)+1
+ call amovs (
+ Mems[imgl2s(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+ Mems[impl2s(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+ # Compute the mean from the sum of the output pixels.
+ if (findmean == YES)
+ MEAN(ccd) = mean / ncols / nlines
+end
+
+
+# FIND_OVERSCAN -- Find the overscan value for a line.
+# No check is made on the number of pixels.
+# The median is the (npix+1)/2 element.
+
+real procedure find_overscans (data, npix, type)
+
+short data[npix] #I Overscan data
+int npix #I Number of overscan points
+int type #I Type of overscan calculation
+
+int i
+real overscan, d, dmin, dmax
+short asoks()
+
+begin
+ if (type == OVERSCAN_MINMAX) {
+ overscan = data[1]
+ dmin = data[1]
+ dmax = data[1]
+ do i = 2, npix {
+ d = data[i]
+ overscan = overscan + d
+ if (d < dmin)
+ dmin = d
+ else if (d > dmax)
+ dmax = d
+ }
+ overscan = (overscan - dmin - dmax) / (npix - 2)
+ } else if (type == OVERSCAN_MEDIAN)
+ overscan = asoks (data, npix, (npix + 1) / 2)
+ else {
+ overscan = data[1]
+ do i = 2, npix
+ overscan = overscan + data[i]
+ overscan = overscan / npix
+ }
+
+ return (overscan)
+end
+
+# PROC1 -- Process CCD images with readout axis 1 (lines).
+
+procedure proc1r (ccd)
+
+pointer ccd # CCD structure
+
+int line, ncols, nlines, findmean, rep
+int overscan_type, overscan_c1, noverscan
+real overscan, darkscale, flatscale, illumscale, frgscale, mean
+real minrep
+pointer in, out, zeroim, darkim, flatim, illumim, fringeim, overscan_vec
+pointer inbuf, outbuf, zerobuf, darkbuf, flatbuf, illumbuf, fringebuf
+
+real asumr()
+real find_overscanr()
+pointer imgl2r(), impl2r(), ccd_glr(), xt_fpsr()
+
+begin
+ # Initialize. If the correction image is 1D then just get the
+ # data once.
+
+ in = IN_IM(ccd)
+ out = OUT_IM(ccd)
+ ncols = OUT_C2(ccd) - OUT_C1(ccd) + 1
+ nlines = OUT_L2(ccd) - OUT_L1(ccd) + 1
+
+ findmean = CORS(ccd, FINDMEAN)
+ if (findmean == YES)
+ mean = 0.
+ rep = CORS(ccd, MINREP)
+ if (rep == YES)
+ minrep = MINREPLACE(ccd)
+
+ if (CORS(ccd, OVERSCAN) == 0)
+ overscan_type = 0
+ else {
+ overscan_type = OVERSCAN_TYPE(ccd)
+ overscan_vec = OVERSCAN_VEC(ccd)
+ overscan_c1 = BIAS_C1(ccd) - 1
+ noverscan = BIAS_C2(ccd) - overscan_c1
+ }
+
+ if (CORS(ccd, ZEROCOR) == 0) {
+ zeroim = NULL
+ zerobuf = 1
+ } else if (IM_LEN(ZERO_IM(ccd),2) == 1) {
+ zeroim = NULL
+ zerobuf = ccd_glr (ZERO_IM(ccd), ZERO_C1(ccd), ZERO_C2(ccd), 1)
+ } else
+ zeroim = ZERO_IM(ccd)
+
+ if (CORS(ccd, DARKCOR) == 0) {
+ darkim = NULL
+ darkbuf = 1
+ } else if (IM_LEN(DARK_IM(ccd),2) == 1) {
+ darkim = NULL
+ darkbuf = ccd_glr (DARK_IM(ccd), DARK_C1(ccd), DARK_C2(ccd), 1)
+ darkscale = FLATSCALE(ccd)
+ } else {
+ darkim = DARK_IM(ccd)
+ darkscale = DARKSCALE(ccd)
+ }
+
+ if (CORS(ccd, FLATCOR) == 0) {
+ flatim = NULL
+ flatbuf = 1
+ } else if (IM_LEN(FLAT_IM(ccd),2) == 1) {
+ flatim = NULL
+ flatbuf = ccd_glr (FLAT_IM(ccd), FLAT_C1(ccd), FLAT_C2(ccd), 1)
+ flatscale = FLATSCALE(ccd)
+ } else {
+ flatim = FLAT_IM(ccd)
+ flatscale = FLATSCALE(ccd)
+ }
+
+ if (CORS(ccd, ILLUMCOR) == 0) {
+ illumim = NULL
+ illumbuf = 1
+ } else {
+ illumim = ILLUM_IM(ccd)
+ illumscale = ILLUMSCALE(ccd)
+ }
+
+ if (CORS(ccd, FRINGECOR) == 0) {
+ fringeim = NULL
+ fringebuf = 1
+ } else {
+ fringeim = FRINGE_IM(ccd)
+ frgscale = FRINGESCALE(ccd)
+ }
+
+ # For each line read lines from the input. Procedure XT_FPS replaces
+ # bad pixels by interpolation. The trimmed region is copied to the
+ # output. Get lines from the output image and from the zero level,
+ # dark count, flat field, illumination, and fringe images. Call COR1
+ # to do the actual pixel corrections. Finally, add the output pixels
+ # to a sum for computing the mean. We must copy data outside of the
+ # output data section.
+
+ do line = 2 - OUT_L1(ccd), 0
+ call amovr (
+ Memr[imgl2r(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+ Memr[impl2r(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+ do line = 1, nlines {
+ outbuf = impl2r (out, OUT_L1(ccd)+line-1)
+
+ inbuf = xt_fpsr (MASK_FP(ccd), in, IN_L1(ccd)+line-1, IN_C1(ccd),
+ IN_C2(ccd), IN_L1(ccd), IN_L2(ccd), NULL)
+ call amovr (Memr[inbuf+IN_C1(ccd)-OUT_C1(ccd)], Memr[outbuf],
+ IM_LEN(out,1))
+
+ outbuf = outbuf + OUT_C1(ccd) - 1
+ if (overscan_type != 0) {
+ if (overscan_type < OVERSCAN_FIT)
+ overscan = find_overscanr (Memr[inbuf+overscan_c1],
+ noverscan, overscan_type)
+ else
+ overscan = Memr[overscan_vec+line-1]
+ }
+ if (zeroim != NULL)
+ zerobuf = ccd_glr (zeroim, ZERO_C1(ccd), ZERO_C2(ccd),
+ ZERO_L1(ccd)+line-1)
+ if (darkim != NULL)
+ darkbuf = ccd_glr (darkim, DARK_C1(ccd), DARK_C2(ccd),
+ DARK_L1(ccd)+line-1)
+ if (flatim != NULL)
+ flatbuf = ccd_glr (flatim, FLAT_C1(ccd), FLAT_C2(ccd),
+ FLAT_L1(ccd)+line-1)
+ if (illumim != NULL)
+ illumbuf = ccd_glr (illumim, ILLUM_C1(ccd), ILLUM_C2(ccd),
+ ILLUM_L1(ccd)+line-1)
+ if (fringeim != NULL)
+ fringebuf = ccd_glr (fringeim, FRINGE_C1(ccd), FRINGE_C2(ccd),
+ FRINGE_L1(ccd)+line-1)
+
+ call cor1r (CORS(ccd,1), Memr[outbuf],
+ overscan, Memr[zerobuf], Memr[darkbuf],
+ Memr[flatbuf], Memr[illumbuf], Memr[fringebuf], ncols,
+ darkscale, flatscale, illumscale, frgscale)
+
+ if (rep == YES)
+ call amaxkr (Memr[outbuf], minrep, Memr[outbuf], ncols)
+ if (findmean == YES)
+ mean = mean + asumr (Memr[outbuf], ncols)
+ }
+
+ do line = nlines+1, IM_LEN(out,2)-OUT_L1(ccd)+1
+ call amovr (
+ Memr[imgl2r(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+ Memr[impl2r(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+ # Compute the mean from the sum of the output pixels.
+ if (findmean == YES)
+ MEAN(ccd) = mean / ncols / nlines
+end
+
+
+# PROC2 -- Process CCD images with readout axis 2 (columns).
+
+procedure proc2r (ccd)
+
+pointer ccd # CCD structure
+
+int line, ncols, nlines, findmean, rep
+real darkscale, flatscale, illumscale, frgscale, mean
+real minrep
+pointer in, out, zeroim, darkim, flatim, illumim, fringeim, overscan_vec
+pointer inbuf, outbuf, zerobuf, darkbuf, flatbuf, illumbuf, fringebuf
+
+real asumr()
+pointer imgl2r(), impl2r(), imgs2r(), ccd_glr(), xt_fpsr()
+
+begin
+ # Initialize. If the correction image is 1D then just get the
+ # data once.
+
+ in = IN_IM(ccd)
+ out = OUT_IM(ccd)
+ ncols = OUT_C2(ccd) - OUT_C1(ccd) + 1
+ nlines = OUT_L2(ccd) - OUT_L1(ccd) + 1
+
+ findmean = CORS(ccd, FINDMEAN)
+ if (findmean == YES)
+ mean = 0.
+ rep = CORS(ccd, MINREP)
+ if (rep == YES)
+ minrep = MINREPLACE(ccd)
+
+ overscan_vec = OVERSCAN_VEC(ccd)
+
+ if (CORS(ccd, ZEROCOR) == 0) {
+ zeroim = NULL
+ zerobuf = 1
+ } else if (IM_LEN(ZERO_IM(ccd),1) == 1) {
+ zeroim = NULL
+ zerobuf = imgs2r (ZERO_IM(ccd), 1, 1, ZERO_L1(ccd), ZERO_L2(ccd))
+ } else
+ zeroim = ZERO_IM(ccd)
+
+ if (CORS(ccd, DARKCOR) == 0) {
+ darkim = NULL
+ darkbuf = 1
+ } else if (IM_LEN(DARK_IM(ccd),1) == 1) {
+ darkim = NULL
+ darkbuf = imgs2r (DARK_IM(ccd), 1, 1, DARK_L1(ccd), DARK_L2(ccd))
+ darkscale = DARKSCALE(ccd)
+ } else {
+ darkim = DARK_IM(ccd)
+ darkscale = DARKSCALE(ccd)
+ }
+
+ if (CORS(ccd, FLATCOR) == 0) {
+ flatim = NULL
+ flatbuf = 1
+ } else if (IM_LEN(FLAT_IM(ccd),1) == 1) {
+ flatim = NULL
+ flatbuf = imgs2r (FLAT_IM(ccd), 1, 1, FLAT_L1(ccd), FLAT_L2(ccd))
+ flatscale = FLATSCALE(ccd)
+ } else {
+ flatim = FLAT_IM(ccd)
+ flatscale = FLATSCALE(ccd)
+ }
+
+ if (CORS(ccd, ILLUMCOR) == 0) {
+ illumim = NULL
+ illumbuf = 1
+ } else {
+ illumim = ILLUM_IM(ccd)
+ illumscale = ILLUMSCALE(ccd)
+ }
+
+ if (CORS(ccd, FRINGECOR) == 0) {
+ fringeim = NULL
+ fringebuf = 1
+ } else {
+ fringeim = FRINGE_IM(ccd)
+ frgscale = FRINGESCALE(ccd)
+ }
+
+ # For each line read lines from the input. Procedure CORINPUT
+ # replaces bad pixels by interpolation and applies a trim to the
+ # input. Get lines from the output image and from the zero level,
+ # dark count, flat field, illumination, and fringe images.
+ # Call COR2 to do the actual pixel corrections. Finally, add the
+ # output pixels to a sum for computing the mean.
+ # We must copy data outside of the output data section.
+
+ do line = 2 - OUT_L1(ccd), 0
+ call amovr (
+ Memr[imgl2r(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+ Memr[impl2r(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+ do line = 1, nlines {
+ outbuf = impl2r (out, OUT_L1(ccd)+line-1)
+
+ inbuf = xt_fpsr (MASK_FP(ccd), in, IN_L1(ccd)+line-1, IN_C1(ccd),
+ IN_C2(ccd), IN_L1(ccd), IN_L2(ccd), NULL)
+ call amovr (Memr[inbuf+IN_C1(ccd)-OUT_C1(ccd)], Memr[outbuf],
+ IM_LEN(out,1))
+
+ outbuf = outbuf + OUT_C1(ccd) - 1
+ if (zeroim != NULL)
+ zerobuf = ccd_glr (zeroim, ZERO_C1(ccd), ZERO_C2(ccd),
+ ZERO_L1(ccd)+line-1)
+ if (darkim != NULL)
+ darkbuf = ccd_glr (darkim, DARK_C1(ccd), DARK_C2(ccd),
+ DARK_L1(ccd)+line-1)
+ if (flatim != NULL)
+ flatbuf = ccd_glr (flatim, FLAT_C1(ccd), FLAT_C2(ccd),
+ FLAT_L1(ccd)+line-1)
+ if (illumim != NULL)
+ illumbuf = ccd_glr (illumim, ILLUM_C1(ccd), ILLUM_C2(ccd),
+ ILLUM_L1(ccd)+line-1)
+ if (fringeim != NULL)
+ fringebuf = ccd_glr (fringeim, FRINGE_C1(ccd), FRINGE_C2(ccd),
+ FRINGE_L1(ccd)+line-1)
+
+ call cor2r (line, CORS(ccd,1), Memr[outbuf],
+ Memr[overscan_vec], Memr[zerobuf], Memr[darkbuf],
+ Memr[flatbuf], Memr[illumbuf], Memr[fringebuf], ncols,
+ zeroim, flatim, darkscale, flatscale, illumscale, frgscale)
+
+ if (rep == YES)
+ call amaxkr (Memr[outbuf], minrep, Memr[outbuf], ncols)
+ if (findmean == YES)
+ mean = mean + asumr (Memr[outbuf], ncols)
+ }
+
+ do line = nlines+1, IM_LEN(out,2)-OUT_L1(ccd)+1
+ call amovr (
+ Memr[imgl2r(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+ Memr[impl2r(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+ # Compute the mean from the sum of the output pixels.
+ if (findmean == YES)
+ MEAN(ccd) = mean / ncols / nlines
+end
+
+
+# FIND_OVERSCAN -- Find the overscan value for a line.
+# No check is made on the number of pixels.
+# The median is the (npix+1)/2 element.
+
+real procedure find_overscanr (data, npix, type)
+
+real data[npix] #I Overscan data
+int npix #I Number of overscan points
+int type #I Type of overscan calculation
+
+int i
+real overscan, d, dmin, dmax
+real asokr()
+
+begin
+ if (type == OVERSCAN_MINMAX) {
+ overscan = data[1]
+ dmin = data[1]
+ dmax = data[1]
+ do i = 2, npix {
+ d = data[i]
+ overscan = overscan + d
+ if (d < dmin)
+ dmin = d
+ else if (d > dmax)
+ dmax = d
+ }
+ overscan = (overscan - dmin - dmax) / (npix - 2)
+ } else if (type == OVERSCAN_MEDIAN)
+ overscan = asokr (data, npix, (npix + 1) / 2)
+ else {
+ overscan = data[1]
+ do i = 2, npix
+ overscan = overscan + data[i]
+ overscan = overscan / npix
+ }
+
+ return (overscan)
+end
generated by cgit (git 2.34.1) at 2025-09-20 10:54:16 -0400