path: root/noao/onedspec/odcombine/srcwt/icomb.gx

# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	<imhdr.h>
include	<imset.h>
include	<pmset.h>
include	<error.h>
include	<syserr.h>
include	<mach.h>
include	"../icombine.h"

# The following is for compiling under V2.11.
define	IM_BUFFRAC	IM_BUFSIZE
include	<imset.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 (sird)
procedure icombine$t (in, out, scales, zeros, wts, wtp, offsets, nimages,
	bufsize)

pointer	in[nimages]		# Input images
pointer	out[ARB]		# Output images
real	scales[nimages]		# Scales
real	zeros[nimages]		# Zeros
real	wts[nimages]		# Weights
pointer	wtp[nimages]		# Weight image pointers
int	offsets[nimages,ARB]	# Input image offsets
int	nimages			# Number of input images
int	bufsize			# IMIO buffer size

char	str[1]
int	i, j, k, npts, fd, stropen(), xt_imgnl$t()
pointer	sp, d, w, id, n, m, lflag, v, dbuf, wbuf
pointer	im, buf, xt_opix(), impl1i()
errchk	stropen, xt_cpix, xt_opix, xt_imgnl$t, impl1i, ic_combine$t
$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 (dbuf, nimages, TY_POINTER)
	call salloc (d, nimages, TY_POINTER)
	call salloc (wbuf, nimages, TY_POINTER)
	call salloc (w, 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 (v, IM_MAXDIM, TY_LONG)
	call amovki (NULL, Memi[dbuf], nimages)
	call amovki (NULL, Memi[d], nimages)
	call amovki (NULL, Memi[wbuf], nimages)
	call amovki (NULL, Memi[w], nimages)
	call amovki (D_ALL, Memi[lflag], nimages)
	call amovkl (1, Meml[v], IM_MAXDIM)

	# If not aligned or growing create data buffers of output length
	# otherwise use the IMIO buffers.

	if (!aligned || grow >= 1.) {
	    do i = 1, nimages
		call salloc (Memi[dbuf+i-1], npts, TY_PIXEL)
	} else {
	    do i = 1, nimages {
		im = xt_opix (in[i], i, 0)
		if (im != in[i])
		    call salloc (Memi[dbuf+i-1], npts, TY_PIXEL)
	    }
	}

	if (project) {
	    call imseti (in[1], IM_NBUFS, nimages)
	    call imseti (in[1], IM_BUFFRAC, 0)
	    call imseti (in[1], IM_BUFSIZE, bufsize)
	    do i = 1, 6 {
		if (out[i] != NULL) {
		    call imseti (out[i], IM_BUFFRAC, 0)
		    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, 6 {
		if (out[i] != NULL) {
		    call imseti (out[i], IM_BUFFRAC, 0)
		    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)
	    }
	    if (out[4] != NULL) {
		buf = impl1i (out[4])
		call aclri (Memi[buf], npts)
	    }
	    if (out[5] != NULL) {
		buf = impl1i (out[5])
		call aclri (Memi[buf], npts)
	    }
	    if (out[6] != NULL) {
		buf = impl1i (out[6])
		call aclri (Memi[buf], npts)
	    }

	    # Do I/O for first input image line.
	    if (!project) {
		do i = 1, nimages {
		    call xt_imseti (i, "bufsize", bufsize)
		    j = max (0, offsets[i,1])
		    k = min (npts, IM_LEN(in[i],1) + offsets[i,1])
		    if (k - j < 1)
			call xt_cpix (i)
		    j = 1 - offsets[i,2]
		    if (j < 1 || j > IM_LEN(in[i],2))
			call xt_cpix (i)
		}

		do i = 1, nimages {
		    j = max (0, offsets[i,1])
		    k = min (npts, IM_LEN(in[i],1) + offsets[i,1])
		    if (k - j < 1)
			next
		    j = 1 - offsets[i,2]
		    if (j < 1 || j > IM_LEN(in[i],2))
			next
		    iferr {
			Meml[v+1] = j
			j = xt_imgnl$t (in[i], i, buf, Meml[v], 1)
		    } then {
			call imseti (im, IM_PIXFD, NULL)
			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, scales, zeros, wts,
	    wtp, Memi[wbuf], Memi[w], nimages, npts)

	call sfree (sp)
end


# IC_COMBINE -- Combine images.

procedure ic_combine$t (in, out, dbuf, d, id, n, m, lflag, offsets,
	scales, zeros, wts, wtp, wbuf, w, 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
pointer	wtp[nimages]		# Combining weight image pointers
pointer	wbuf[nimages]		# Weight buffers for nonaligned images
pointer	w[nimages]		# Weight pointers
int	nimages			# Number of input images
int	npts			# Number of points per output line

int	i, ext, ctor(), errcode()
real	r, imgetr()
pointer	sp, fname, imname, v1, v2, v3, work
pointer	outdata, buf, nm, pms
pointer	immap(), impnli(), xt_opix()
$if (datatype == sil)
pointer	impnlr(), imgnlr()
$else
pointer	impnl$t(), imgnl$t
$endif
errchk	immap, ic_scale, xt_opix, imgetr, ic_grow, ic_rmasks
errchk	ic_grow$t, ic_gdata$t

include	"../icombine.com"
data	ext/0/

begin
	call smark (sp)
	call salloc (fname, SZ_FNAME, TY_CHAR)
	call salloc (imname, SZ_FNAME, TY_CHAR)
	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, wtp, nimages, npts)

	# Allocate weight buffers if needed.
	if (wtype == S_WTMAP || wtype == S_SIGMAP) {
	    if (!aligned) {
		do i = 1, nimages
		    call salloc (wbuf[i], npts, TY_REAL)
	    } else {
		do i = 1, nimages {
		    if (wtp[i] != xt_opix (wtp[i], nimages+i, 0))
			call salloc (wbuf[i], npts, TY_REAL)
		}
	    }
	}

	# 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)
		    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
	case PCLIP:
	    mclip = true
	case AVSIGCLIP, SIGCLIP:
	    if (doscale1)
		keepids = true
	case NONE:
	    mclip = false
	}

	if (out[4] != NULL)
	    keepids = true

	if (out[6] != NULL) {
	    keepids = true
	    call ic_einit (in, nimages, Memc[expkeyword], 1., 2**27-1)
	}

	if (grow >= 1.) {
	    keepids = true
	    call salloc (work, npts * nimages, TY_INT)
	}
	pms = NULL

	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, wtp, wbuf, w, 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 (pms == NULL || nkeep > 0) {
		if (docombine) {
		    switch (combine) {
		    case AVERAGE:
			call ic_average$t (d, id, n, wts, w, npts, YES, YES,
			    Memr[outdata])
		    case MEDIAN:
			call ic_median$t (d, n, npts, YES, Memr[outdata])
		    case SUM:
			call ic_average$t (d, id, n, wts, w, npts, YES, NO,
			    Memr[outdata])
		    }
		}
	    }

	    if (grow >= 1.)
		call ic_grow (out, Meml[v2], id, n, Memi[work], nimages, npts,
		    pms)

	    if (pms == NULL) {
		if (out[2] != NULL) {
		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
		    i = impnli (out[2], buf, Meml[v1])
		    do i = 1, npts {
			if (n[i] == 0)
			    Memi[buf] = 1
			else
			    Memi[buf] = 0
		    }
		}

		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, w, npts, Memr[outdata],
			Memr[buf])
		}

		if (out[4] != NULL)
		    call ic_rmasks (out[4], Meml[v2], id, nimages, n, npts)

		if (out[5] != NULL) {
		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
		    i = impnli (out[5], buf, Meml[v1])
		    call amovki (nimages, Memi[buf], npts)
		    call asubi (Memi[buf], n, Memi[buf], npts)
		}

		if (out[6] != NULL)
		    call ic_emask (out[6], Meml[v2], id, nimages, n, wts, w,
			npts)
	    }

	    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, wtp, wbuf, w, 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 (pms == NULL || nkeep > 0) {
		if (docombine) {
		    switch (combine) {
		    case AVERAGE:
			call ic_average$t (d, id, n, wts, w, npts, YES, YES,
			    Mem$t[outdata])
		    case MEDIAN:
			call ic_median$t (d, n, npts, YES, Mem$t[outdata])
		    case SUM:
			call ic_average$t (d, id, n, wts, w, npts, YES, NO,
			    Mem$t[outdata])
		    }
		}
	    }

	    if (grow >= 1.)
		call ic_grow (out, Meml[v2], id, n, Memi[work], nimages, npts,
		    pms)

	    if (pms == NULL) {
		if (out[2] != NULL) {
		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
		    i = impnli (out[2], buf, Meml[v1])
		    do i = 1, npts {
			if (n[i] == 0)
			    Memi[buf] = 1
			else
			    Memi[buf] = 0
			buf = buf + 1
		    }
		}

		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, w, npts, Mem$t[outdata],
			Mem$t[buf])
		}

		if (out[4] != NULL)
		    call ic_rmasks (out[4], Meml[v2], id, nimages, n, npts)

		if (out[5] != NULL) {
		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
		    i = impnli (out[5], buf, Meml[v1])
		    call amovki (nimages, Memi[buf], npts)
		    call asubi (Memi[buf], n, Memi[buf], npts)
		}

		if (out[6] != NULL)
		    call ic_emask (out[6], Meml[v2], id, nimages, n, wts, w,
			npts)
	    }

	    call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
	}
	$endif

	if (pms != NULL) {
	    if (nkeep > 0) {
		call imstats (out[1], IM_IMAGENAME, Memc[fname], SZ_FNAME)
		call imunmap (out[1])
		iferr (buf = immap (Memc[fname], READ_WRITE, 0)) {
		    switch (errcode()) {
		    case SYS_FXFOPNOEXTNV:
			call imgcluster (Memc[fname], Memc[fname], SZ_FNAME)
			ext = ext + 1
			call sprintf (Memc[imname], SZ_FNAME, "%s[%d]")
			    call pargstr (Memc[fname])
			    call pargi (ext)
			iferr (buf = immap (Memc[imname], READ_WRITE, 0)) {
			    buf = NULL
			    ext = 0
			}
			repeat {
			    call sprintf (Memc[imname], SZ_FNAME, "%s[%d]")
				call pargstr (Memc[fname])
				call pargi (ext+1)
			    iferr (outdata = immap (Memc[imname],READ_WRITE,0)) 
				break
			    if (buf != NULL)
				call imunmap (buf)
			    buf = outdata
			    ext = ext + 1
			}
		    default:
			call erract (EA_ERROR)
		    }
		}
		out[1] = buf
	    }

	    call amovkl (long(1), Meml[v1], IM_MAXDIM)
	    call amovkl (long(1), Meml[v2], IM_MAXDIM)
	    call amovkl (long(1), Meml[v3], IM_MAXDIM)
	    $if (datatype == sil)
	    while (impnlr (out[1], outdata, Meml[v1]) != EOF) {
		call ic_gdata$t (in, out, dbuf, d, wtp, wbuf, w, id, n, m,
		    lflag, offsets, scales, zeros, nimages, npts,
		    Meml[v2], Meml[v3])

		call ic_grow$t (Meml[v2], d, id, n, Memi[work], nimages, npts,
		    pms)

		if (nkeep > 0) {
		    do i = 1, npts {
			if (n[i] < nkeep) {
			    Meml[v1+1] = Meml[v1+1] - 1
			    if (imgnlr (out[1], buf, Meml[v1]) == EOF)
				;
			    call amovr (Memr[buf], Memr[outdata], npts)
			    break
			}
		    }
		}

		switch (combine) {
		case AVERAGE:
		    call ic_average$t (d, id, n, wts, w, npts, NO, YES,
			Memr[outdata])
		case MEDIAN:
		    call ic_median$t (d, n, npts, NO, Memr[outdata])
		case SUM:
		    call ic_average$t (d, id, n, wts, w, npts, NO, NO,
			Memr[outdata])
		}

		if (out[2] != NULL) {
		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
		    i = impnli (out[2], buf, Meml[v1])
		    do i = 1, npts {
			if (n[i] == 0)
			    Memi[buf] = 1
			else
			    Memi[buf] = 0
		    }
		}
			
		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, w, npts, Memr[outdata],
			Memr[buf])
		}

		if (out[4] != NULL)
		    call ic_rmasks (out[4], Meml[v2], id, nimages, n, npts)

		if (out[5] != NULL) {
		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
		    i = impnli (out[5], buf, Meml[v1])
		    call amovki (nimages, Memi[buf], npts)
		    call asubi (Memi[buf], n, Memi[buf], npts)
		}

		if (out[6] != NULL)
		    call ic_emask (out[6], Meml[v2], id, nimages, n, wts, w,
			npts)

		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, wtp, wbuf, w, id, n, m,
		    lflag, offsets, scales, zeros, nimages, npts,
		    Meml[v2], Meml[v3])

		call ic_grow$t (Meml[v2], d, id, n, Memi[work], nimages, npts,
		    pms)

		if (nkeep > 0) {
		    do i = 1, npts {
			if (n[i] < nkeep) {
			    Meml[v1+1] = Meml[v1+1] - 1
			    if (imgnl$t (out[1], buf, Meml[v1]) == EOF)
				;
			    call amov$t (Mem$t[buf], Mem$t[outdata], npts)
			    break
			}
		    }
		}

		switch (combine) {
		case AVERAGE:
		    call ic_average$t (d, id, n, wts, w, npts, NO, YES,
			Mem$t[outdata])
		case MEDIAN:
		    call ic_median$t (d, n, npts, NO, Mem$t[outdata])
		case SUM:
		    call ic_average$t (d, id, n, wts, w, npts, NO, NO,
			Mem$t[outdata])
		}

		if (out[2] != NULL) {
		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
		    i = impnli (out[2], buf, Meml[v1])
		    do i = 1, npts {
			if (n[i] == 0)
			    Memi[buf] = 1
			else
			    Memi[buf] = 0
		    }
		}
			
		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, w, npts, Mem$t[outdata],
			Mem$t[buf])
		}

		if (out[4] != NULL)
		    call ic_rmasks (out[4], Meml[v2], id, nimages, n, npts)

		if (out[5] != NULL) {
		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
		    i = impnli (out[5], buf, Meml[v1])
		    call amovki (nimages, Memi[buf], npts)
		    call asubi (Memi[buf], n, Memi[buf], npts)
		}

		if (out[6] != NULL)
		    call ic_emask (out[6], Meml[v2], id, nimages, n, wts, w,
			npts)

		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
	    }
	    $endif

	    do i = 1, nimages
		call pm_close (Memi[pms+i-1])
	    call mfree (pms, TY_POINTER)
	}

	if (wtype == S_WTMAP || wtype == S_SIGMAP) {
	    do i = 1, nimages
		call xt_imunmap (wtp[i], nimages+i)
	}
	call sfree (sp)
end
$endfor