Initial commit

1 files changed, 395 insertions, 0 deletions

diff --git a/noao/imred/ccdred/src/icombine.gx b/noao/imred/ccdred/src/icombine.gx
new file mode 100644
index 00000000..d6e93ef0
--- /dev/null
+++ b/noao/imred/ccdred/src/icombine.gx
@@ -0,0 +1,395 @@
+# 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.
+
+$for (sr)
+procedure icombine$t (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, imgl1$t(), impl1i()
+errchk	stropen, imgl1$t, impl1i
+$if (datatype == sil)
+pointer	impl1r()
+errchk	impl1r
+$else
+pointer	impl1$t()
+errchk	impl1$t
+$endif
+
+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_PIXEL)
+	}
+
+	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)
+	    }
+	    $if (datatype == sil)
+	    buf = impl1r (out[1])
+	    call aclrr (Memr[buf], npts)
+	    if (out[3] != NULL) {
+		buf = impl1r (out[3])
+		call aclrr (Memr[buf], npts)
+	    }
+	    $else
+	    buf = impl1$t (out[1])
+	    call aclr$t (Mem$t[buf], npts)
+	    if (out[3] != NULL) {
+		buf = impl1$t (out[3])
+		call aclr$t (Mem$t[buf], npts)
+	    }
+	    $endif
+	    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 = imgl1$t (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 = imgl1$t (in[i])
+		    default:
+			call sfree (sp)
+			call strclose (fd)
+			call erract (EA_ERROR)
+		    }
+		}
+	    }
+
+	    call strclose (fd)
+	}
+
+	call ic_combine$t (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_combine$t (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()
+$if (datatype == sil)
+pointer	impnlr()
+$else
+pointer	impnl$t()
+$endif
+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)
+	}
+
+	$if (datatype == sil)
+	while (impnlr (out[1], outdata, Meml[v1]) != EOF) {
+	    call ic_gdata$t (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_mccdclip$t (d, id, n, scales, zeros, Memr[nm],
+			nimages, npts, Memr[outdata])
+		else
+		    call ic_accdclip$t (d, id, n, scales, zeros, Memr[nm],
+			nimages, npts, Memr[outdata])
+	    case MINMAX:
+		call ic_mm$t (d, id, n, npts)
+	    case PCLIP:
+		call ic_pclip$t (d, id, n, nimages, npts, Memr[outdata])
+	    case SIGCLIP:
+		if (mclip)
+		    call ic_msigclip$t (d, id, n, scales, zeros, nimages, npts,
+			Memr[outdata])
+		else
+		    call ic_asigclip$t (d, id, n, scales, zeros, nimages, npts,
+			Memr[outdata])
+	    case AVSIGCLIP:
+		if (mclip)
+		    call ic_mavsigclip$t (d, id, n, scales, zeros, nimages,
+			npts, Memr[outdata])
+		else
+		    call ic_aavsigclip$t (d, id, n, scales, zeros, nimages,
+			npts, Memr[outdata])
+	    }
+
+	    if (grow > 0)
+		call ic_grow$t (d, id, n, nimages, npts, Memr[outdata])
+
+	    if (docombine) {
+		switch (combine) {
+		case AVERAGE:
+		    call ic_average$t (d, id, n, wts, npts, Memr[outdata])
+		case MEDIAN:
+		    call ic_median$t (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_sigma$t (d, id, n, wts, npts, Memr[outdata],
+		    Memr[buf])
+	    }
+	    call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+	}
+	$else
+	while (impnl$t (out[1], outdata, Meml[v1]) != EOF) {
+	    call ic_gdata$t (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_mccdclip$t (d, id, n, scales, zeros, Memr[nm],
+			nimages, npts, Mem$t[outdata])
+		else
+		    call ic_accdclip$t (d, id, n, scales, zeros, Memr[nm],
+			nimages, npts, Mem$t[outdata])
+	    case MINMAX:
+		call ic_mm$t (d, id, n, npts)
+	    case PCLIP:
+		call ic_pclip$t (d, id, n, nimages, npts, Mem$t[outdata])
+	    case SIGCLIP:
+		if (mclip)
+		    call ic_msigclip$t (d, id, n, scales, zeros, nimages, npts,
+			Mem$t[outdata])
+		else
+		    call ic_asigclip$t (d, id, n, scales, zeros, nimages, npts,
+			Mem$t[outdata])
+	    case AVSIGCLIP:
+		if (mclip)
+		    call ic_mavsigclip$t (d, id, n, scales, zeros, nimages,
+			npts, Mem$t[outdata])
+		else
+		    call ic_aavsigclip$t (d, id, n, scales, zeros, nimages,
+			npts, Mem$t[outdata])
+	    }
+
+	    if (grow > 0)
+		call ic_grow$t (d, id, n, nimages, npts, Mem$t[outdata])
+
+	    if (docombine) {
+		switch (combine) {
+		case AVERAGE:
+		    call ic_average$t (d, id, n, wts, npts, Mem$t[outdata])
+		case MEDIAN:
+		    call ic_median$t (d, n, npts, Mem$t[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 = impnl$t (out[3], buf, Meml[v1])
+		call ic_sigma$t (d, id, n, wts, npts, Mem$t[outdata],
+		    Mem$t[buf])
+	    }
+	    call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+	}
+	$endif
+
+	call sfree (sp)
+end
+$endfor