Repatch (from linux) of OSX IRAF

1 files changed, 971 insertions, 0 deletions

diff --git a/noao/nproto/ace/convolve.x b/noao/nproto/ace/convolve.x
new file mode 100644
index 00000000..af9734ef
--- /dev/null
+++ b/noao/nproto/ace/convolve.x
@@ -0,0 +1,971 @@
+include	<ctype.h>
+include	<imhdr.h>
+
+
+# ODCNV -- Get a line of data possibly convolved.  Also get the unconvolved
+# data, the sky data, and the sky sigma data.
+#
+# This routine must be called sequentially starting with the first line.
+# It is initialized when the first line.  Memory is freed by using a final
+# call with a line of zero.
+
+procedure convolve (in, bpm, sky, sig, exp, offset, scale, line, cnv,
+	indata, bp, cnvdata, skydata, sigdata, expdata, cnvwt, logfd)
+
+pointer	in[2]		#I Image pointers
+pointer	bpm[2]		#I BPM pointer
+pointer	sky[2]		#I Sky map
+pointer	sig[2]		#I Sigma map
+pointer	exp[2]		#I Exposure map
+int	offset[2]	#I Offsets
+real	scale[2]	#I Image scales
+int	line		#I Line
+char	cnv[ARB]	#I Convolution string
+pointer	indata[2]	#O Pointers to unconvolved image data
+pointer	bp		#O Bad pixel data
+pointer	cnvdata		#O Pointer to convolved image data
+pointer	skydata[2]	#O Pointer to sky data
+pointer	sigdata[2]	#O Pointer to sigma data corrected by exposure map
+pointer	expdata[2]	#O Pointer to exposure map data
+real	cnvwt		#O Weight for convolved sigma
+int	logfd		#I Logfile
+
+int	i, j, k,  nx, ny, nx2, ny2, nc, nl, mode, off
+real	wts, wts1
+pointer	bpm2, kptr, ptr, symptr, symwptr
+bool	dobpm, overlap, fp_equalr()
+
+pointer	kernel, sym, symbuf, symwts, buf, buf2, buf3, bpbuf, bpwts, wtsl, scales
+data	kernel/NULL/, sym/NULL/, symbuf/NULL/, symwts/NULL/
+data	buf/NULL/, buf2/NULL/, buf3/NULL/, bpbuf/NULL/, bpwts/NULL/
+data	wtsl/NULL/, scales/NULL/
+
+errchk	cnvparse, cnvgline2
+
+begin
+	# If no convolution.
+	if (cnv[1] == EOS) {
+	    if (line == 0)
+		return
+
+	    call cnvgline1 (line, offset, in, bpm, indata, bp)
+	    call cnvgline2 (line, offset, in, sky, sig, exp, skydata,
+		sigdata, expdata)
+	    cnvwt = 1
+	    if (in[2] == NULL)
+		cnvdata = indata[1]
+	    else
+		call asubr_scale (Memr[indata[1]], scale[1],
+		    Memr[indata[2]], scale[2], Memr[cnvdata], IM_LEN(in[1],1))
+	    return
+	}
+
+	# Free memory.
+	if (line == 0) {
+	    if (symbuf != NULL) {
+		do i = 0, ARB {
+		    ptr = Memi[symbuf+i]
+		    if (ptr == -1)
+			break
+		    call mfree (ptr, TY_REAL)
+		}
+	    }
+	    if (symwts != NULL) {
+		do i = 0, ARB {
+		    ptr = Memi[symwts+i]
+		    if (ptr == -1)
+			break
+		    call mfree (ptr, TY_REAL)
+		}
+	    }
+	    call mfree (scales, TY_REAL)
+	    call mfree (wtsl, TY_REAL)
+	    call mfree (kernel, TY_REAL)
+	    call mfree (scales, TY_REAL)
+	    call mfree (sym, TY_INT)
+	    call mfree (symbuf, TY_POINTER)
+	    call mfree (symwts, TY_POINTER)
+	    call mfree (buf, TY_REAL)
+	    call mfree (buf2, TY_REAL)
+	    call mfree (buf3, TY_REAL)
+	    call mfree (bpbuf, TY_INT)
+	    call mfree (bpwts, TY_REAL)
+
+	    return
+	}
+
+	# Initialize by getting the kernel coefficients, setting the
+	# image I/O buffers using a scrolling array, and allocate memory.
+
+	if (line == 1 || buf == NULL) {
+	    if (buf != NULL) {
+		if (symbuf != NULL) {
+		    do i = 0, ARB {
+			ptr = Memi[symbuf+i]
+			if (ptr == -1)
+			    break
+			call mfree (ptr, TY_REAL)
+		    }
+		}
+		if (symwts != NULL) {
+		    do i = 0, ARB {
+			ptr = Memi[symwts+i]
+			if (ptr == -1)
+			    break
+			call mfree (ptr, TY_REAL)
+		    }
+		}
+		call mfree (scales, TY_REAL)
+		call mfree (wtsl, TY_REAL)
+		call mfree (kernel, TY_REAL)
+		call mfree (scales, TY_REAL)
+		call mfree (sym, TY_INT)
+		call mfree (symbuf, TY_POINTER)
+		call mfree (symwts, TY_POINTER)
+		call mfree (buf, TY_REAL)
+		call mfree (buf2, TY_REAL)
+		call mfree (buf3, TY_REAL)
+		call mfree (bpbuf, TY_INT)
+		call mfree (bpwts, TY_REAL)
+	    }
+
+	    nc = IM_LEN(in[1],1)
+	    nl = IM_LEN(in[1],2)
+
+	    call cnvparse (cnv, kernel, nx, ny, logfd)
+	    nx2 = nx / 2
+	    ny2 = ny / 2
+	    call malloc (scales, ny, TY_REAL)
+	    call calloc (wtsl, ny, TY_REAL)
+	    call amovkr (1., Memr[scales], ny)
+
+	    # Check for lines which are simple scalings of the first line.
+	    do i = 2, ny {
+		kptr = kernel + (i - 1) * nx
+		wts = 0.
+		do k = 0, nx-1 {
+		    if (Memr[kptr+k] == 0. || Memr[kernel+k] == 0.) {
+			wts = 0.
+			break
+		    }
+		    if (wts == 0.)
+			wts = Memr[kptr+k] / Memr[kernel+k]
+		    else {
+			wts1 = Memr[kptr+k] / Memr[kernel+k]
+			if (!fp_equalr (wts, wts1))
+			    break
+		    }
+		}
+		if (wts != 0. && fp_equalr (wts, wts1)) {
+		    Memr[scales+i-1] = wts
+		    call amovr (Memr[kernel], Memr[kptr], nx)
+		}
+	    }
+
+	    wts = 0
+	    do i = 1, ny {
+		kptr = kernel + (i - 1) * nx
+		wts1 = 0.
+		do j = 1, nx {
+		    wts1 = wts1 + Memr[kptr]
+		    kptr = kptr + 1
+		}
+		Memr[wtsl+i-1] = wts1
+		wts = wts + wts1
+	    }
+	    if (wts != 0.) {
+		call adivkr (Memr[wtsl], wts, Memr[wtsl], ny)
+		call adivkr (Memr[kernel], wts, Memr[kernel], nx*ny)
+	    }
+	    cnvwt = sqrt (wts)
+
+	    if (in[2] == NULL)
+		bpm2 = NULL
+	    else
+		bpm2 = bpm[2] 
+	    if (bpm[1] == NULL && bpm2 == NULL)
+		dobpm = false
+	    else
+		dobpm = true
+	    if (dobpm) {
+		call malloc (bpbuf, nc*ny, TY_INT)
+		call malloc (bpwts, nc, TY_REAL)
+		call calloc (symwts, ny*ny+1, TY_POINTER)
+		Memi[symwts+ny*ny] = -1
+	    }
+
+	    # Check for any line symmetries in the kernel.
+	    call malloc (sym, ny, TY_INT)
+	    call calloc (symbuf, ny*ny+1, TY_POINTER)
+	    Memi[symbuf+ny*ny] = -1
+	    do i = ny, 1, -1 {
+		kptr = kernel + (i - 1) * nx
+		do j = ny, 1, -1 {
+		    ptr = kernel + (j - 1) * nx
+		    do k = 0, nx-1 {
+			if (Memr[kptr+k] != Memr[ptr+k])
+			    break
+		    }
+		    if (k == nx) {
+			Memi[sym+i-1] = j
+			break
+		    }
+		}
+	    }
+	    do i = ny, 1, -1 {
+		k = 0
+		do j = ny, 1, -1
+		    if (Memi[sym+j-1] == i)
+			k = k + 1
+		if (k == 1)
+		    Memi[sym+i-1] = 0
+	    }
+
+	    call malloc (buf, nc*ny, TY_REAL)
+	    if (in[2] != NULL) {
+		call malloc (buf2, nc*ny, TY_REAL)
+		call malloc (buf3, nc*ny, TY_REAL)
+	    }
+
+	    if (in[2] != NULL) {
+		overlap = true
+		if (1-offset[1] < 1 || nc-offset[1] > IM_LEN(in[2],1))
+		    overlap = false
+		if (1-offset[2] < 1 || nl-offset[2] > IM_LEN(in[2],2))
+		    overlap = false
+	    }
+	    do i = 1, ny {
+		call cnvgline1 (i, offset, in, bpm, indata, bp)
+		off = mod (i, ny) * nc
+		call amovr (Memr[indata[1]], Memr[buf+off], nc)
+		if (in[2] != NULL) {
+		    call amovr (Memr[indata[2]], Memr[buf2+off], nc)
+		    call asubr_scale (Memr[buf+off], scale[1],
+			Memr[buf2+off], scale[2], Memr[buf3+off], nc)
+		}
+		if (dobpm)
+		    call amovi (Memi[bp], Memi[bpbuf+off], nc)
+	    }
+	}
+
+	# Get new line.
+	j = line +  ny2
+	if (j > ny && j <= nl) {
+	    call cnvgline1 (j, offset, in, bpm, indata, bp)
+	    off = mod (j, ny) * nc
+	    call amovr (Memr[indata[1]], Memr[buf+off], nc)
+	    if (in[2] != NULL) {
+		call amovr (Memr[indata[2]], Memr[buf2+off], nc)
+		call asubr_scale (Memr[buf+off], scale[1],
+		    Memr[buf2+off], scale[2], Memr[buf3+off], nc)
+	    }
+	    if (dobpm) {
+		ptr = bpbuf + off
+		call amovi (Memi[bp], Memi[ptr], nc)
+	    }
+	}
+
+	# Compute the convolution vector with boundary reflection.
+	# Save and reuse lines with the same kernel weights apart
+	# from a scale factor.
+
+	kptr = kernel
+	call aclrr (Memr[cnvdata], nc)
+	if (dobpm)
+	    call aclrr (Memr[bpwts], nc)
+	do i = 1, ny {
+	    j = line + i - ny2 - 1
+	    if (j < 1)
+		j = 2 - j
+	    else if (j > nl)
+		j = 2 * nl - j
+	    off = mod (j, ny) * nc
+	    if (in[2] == NULL)
+		ptr = buf
+	    else
+		ptr = buf3
+	    k = Memi[sym+i-1]
+	    if (k == 0) {
+		mode = 1
+		symptr = ptr
+		symwptr = bpwts
+	    } else {
+		if (k == i)
+		    mode = 2
+		else
+		    mode = 3
+		symptr = Memi[symbuf+(k-1)*ny+mod(j,ny)]
+		if (symptr == NULL) {
+		    call malloc (symptr, nc, TY_REAL)
+		    Memi[symbuf+(k-1)*ny+mod(j,ny)] = symptr
+		    mode = 2
+		}
+		if (dobpm) {
+		    symwptr = Memi[symwts+(k-1)*ny+mod(j,ny)]
+		    if (symwptr == NULL) {
+			call malloc (symwptr, nc, TY_REAL)
+			Memi[symwts+(k-1)*ny+mod(j,ny)] = symwptr
+		    }
+		}
+	    }
+	    if (dobpm)
+		call convolve2 (Memr[ptr+off], Memr[cnvdata], Memr[symptr],
+		    nc, Memr[kptr], Memr[scales+i-1], nx, Memi[bpbuf+off],
+		    Memr[wtsl+i-1], Memr[bpwts], Memr[symwptr], mode)
+	    else
+		call convolve1 (Memr[ptr+off], Memr[cnvdata], Memr[symptr],
+		    nc, Memr[kptr], Memr[scales+i-1], nx, mode)
+	    kptr = kptr + nx
+	}
+	if (dobpm) {
+	    do i = 0, nc-1
+		if (Memr[bpwts+i] != 0.)
+		    Memr[cnvdata+i] = Memr[cnvdata+i] / Memr[bpwts+i]
+	}
+
+	# Set the output vectors.
+	off = mod (line, ny) * nc
+	indata[1] = buf + off 
+	if (dobpm) {
+	    if (bpm2 == NULL)
+		bp = bpbuf + off
+	    else
+		call amovi (Memi[bpbuf+off], Memi[bp], nc)
+	}
+	if (in[2] != NULL) {
+	    if (overlap)
+		indata[2] = buf2 + off
+	    else
+		call amovr (Memr[buf2+off], Memr[indata[2]], nc)
+	}
+	call cnvgline2 (line, offset, in, sky, sig, exp, skydata, sigdata,
+	   expdata)
+end
+	
+
+
+# ODCNV1 --  One dimensional convolution with boundary reflection.
+#
+# The convolution is added to the output so that it might be used
+# as part of a 2D convolution.
+
+procedure convolve1 (in, out, save, nc, xkernel, scale, nx, mode)
+
+real	in[nc]			#I Input data to be convolved
+real	out[nc]			#O Output convolved data
+real	save[nc]		#U Output saved data
+int	nc			#I Number of data points
+real	xkernel[nx]		#I Convolution weights
+real	scale			#I Scale for saved vector
+int	nx			#I Number of convolution points (must be odd)
+int	mode			#I Mode (1=no save, 2=save, 3=use save)
+
+int	i, j, k, nx2
+real	val
+bool	fp_equalr()
+
+begin
+	if (mode == 1) {
+	    nx2 = nx / 2
+	    do i = 1, nx2 {
+		val = 0
+		do j = 1, nx {
+		    k = i + j - nx2 - 1
+		    if (k < 1)
+			k = 2 - k
+		    val = val + in[k] * xkernel[j]
+		}
+		out[i] = out[i] + val
+	    }
+	    do i = nx2+1, nc-nx2 {
+		k = i - nx2
+		val = 0
+		do j = 1, nx {
+		    val = val + in[k] * xkernel[j]
+		    k = k + 1
+		}
+		out[i] = out[i] + val
+	    }
+	    do i = nc-nx2+1, nc {
+		val = 0
+		do j = 1, nx {
+		    k = i + j - nx2 - 1
+		    if (k > nc)
+			k = 2 * nc - k
+		    val = val + in[k] * xkernel[j]
+		}
+		out[i] = out[i] + val
+	    }
+	} else if (mode == 2) {
+	    nx2 = nx / 2
+	    do i = 1, nx2 {
+		val = 0
+		do j = 1, nx {
+		    k = i + j - nx2 - 1
+		    if (k < 1)
+			k = 2 - k
+		    val = val + in[k] * xkernel[j]
+		}
+		out[i] = out[i] + val
+		save[i] = val
+	    }
+	    do i = nx2+1, nc-nx2 {
+		k = i - nx2
+		val = 0
+		do j = 1, nx {
+		    val = val + in[k] * xkernel[j]
+		    k = k + 1
+		}
+		out[i] = out[i] + val
+		save[i] = val
+	    }
+	    do i = nc-nx2+1, nc {
+		val = 0
+		do j = 1, nx {
+		    k = i + j - nx2 - 1
+		    if (k > nc)
+			k = 2 * nc - k
+		    val = val + in[k] * xkernel[j]
+		}
+		out[i] = out[i] + val
+		save[i] = val
+	    }
+	} else {
+	    if (fp_equalr (1., scale)) {
+		do i = 1, nc
+		    out[i] = out[i] + save[i]
+	    } else {
+		do i = 1, nc
+		    out[i] = out[i] + scale * save[i]
+	    }
+	}
+end
+
+
+# ODCNV2 --  One dimensional convolution with boundary reflection and masking.
+#
+# The convolution is added to the output so that it might be used
+# as part of a 2D convolution.
+
+procedure convolve2 (in, out, save, nc, xkernel, scale, nx, bp,
+	wtssum, wts, wtsave, mode)
+
+real	in[nc]			#I Input data to be convolved
+real	out[nc]			#O Output convolved data
+real	save[nc]		#U Output saved data
+int	nc			#I Number of data points
+real	xkernel[nx]		#I Convolution weights
+real	scale			#I Scale for saved vector
+int	nx			#I Number of convolution points (must be odd)
+int	bp[nc]			#I Bad pixel data
+real	wtssum			#I Sum of weights
+real	wts[nc]			#I Weights
+real	wtsave[nc]		#U Output saved weight data
+int	mode			#I Mode (1=no save, 2=save, 3=use save)
+
+int	i, j, k, nx2
+real	val, wt
+bool	fp_equalr()
+
+begin
+	if (mode == 1) {
+	    nx2 = nx / 2
+	    do i = 1, nx2 {
+		val = 0
+		wt = wtssum
+		do j = 1, nx {
+		    k = i + j - nx2 - 1
+		    if (k < 1)
+			k = 2 - k
+		    if (bp[k] == 0)
+			val = val + in[k] * xkernel[j]
+		    else
+			wt = wt - xkernel[j]
+		}
+		out[i] = out[i] + val
+		wts[i] = wts[i] + wt
+	    }
+	    do i = nx2+1, nc-nx2 {
+		k = i - nx2
+		val = 0
+		wt = wtssum
+		do j = 1, nx {
+		    if (bp[k] == 0)
+			val = val + in[k] * xkernel[j]
+		    else
+			wt = wt - xkernel[j]
+		    k = k + 1
+		}
+		out[i] = out[i] + val
+		wts[i] = wts[i] + wt
+	    }
+	    do i = nc-nx2+1, nc {
+		val = 0
+		wt = wtssum
+		do j = 1, nx {
+		    k = i + j - nx2 - 1
+		    if (k > nc)
+			k = 2 * nc - k
+		    if (bp[k] == 0)
+			val = val + in[k] * xkernel[j]
+		    else
+			wt = wt - xkernel[j]
+		}
+		out[i] = out[i] + val
+		wts[i] = wts[i] + wt
+	    }
+	} else if (mode == 2) {
+	    nx2 = nx / 2
+	    do i = 1, nx2 {
+		val = 0
+		wt = wtssum
+		do j = 1, nx {
+		    k = i + j - nx2 - 1
+		    if (k < 1)
+			k = 2 - k
+		    if (bp[k] == 0)
+			val = val + in[k] * xkernel[j]
+		    else
+			wt = wt - xkernel[j]
+		}
+		out[i] = out[i] + val
+		wts[i] = wts[i] + wt
+		save[i] = val
+		wtsave[i] = wt
+	    }
+	    do i = nx2+1, nc-nx2 {
+		k = i - nx2
+		val = 0
+		wt = wtssum
+		do j = 1, nx {
+		    if (bp[k] == 0)
+			val = val + in[k] * xkernel[j]
+		    else
+			wt = wt - xkernel[j]
+		    k = k + 1
+		}
+		out[i] = out[i] + val
+		wts[i] = wts[i] + wt
+		save[i] = val
+		wtsave[i] = wt
+	    }
+	    do i = nc-nx2+1, nc {
+		val = 0
+		wt = wtssum
+		do j = 1, nx {
+		    k = i + j - nx2 - 1
+		    if (k > nc)
+			k = 2 * nc - k
+		    if (bp[k] == 0)
+			val = val + in[k] * xkernel[j]
+		    else
+			wt = wt - xkernel[j]
+		}
+		out[i] = out[i] + val
+		wts[i] = wts[i] + wt
+		save[i] = val
+		wtsave[i] = wt
+	    }
+	} else {
+	    if (fp_equalr (1., scale)) {
+		do i = 1, nc {
+		    out[i] = out[i] + save[i]
+		    wts[i] = wts[i] + wtsave[i]
+		}
+	    } else {
+		do i = 1, nc {
+		    out[i] = out[i] + scale * save[i]
+		    wts[i] = wts[i] + scale * wtsave[i]
+		}
+	    }
+	}
+end
+
+
+# ASUBR_SCALE -- out = in1 * scale1 - in2 * scale2
+
+procedure asubr_scale (in1, scale1, in2, scale2, out, n)
+
+real	in1[n]			#I Input vector
+real	scale1			#I Scale
+real	in2[n]			#I Input vector
+real	scale2			#I Scale
+real	out[n]			#O Output vector
+int	n			#I Number of points
+
+int	i
+
+begin
+	if (scale1 == 1. && scale2 == 1.)
+	    call asubr (in1, in2, out, n)
+	else if (scale1 == 1.) {
+	    do i = 1, n
+		out[i] = in1[i] - in2[i] * scale2
+	} else if (scale2 == 1.) {
+	    do i = 1, n
+		out[i] = in1[i] * scale1 - in2[i]
+	} else {
+	    do i = 1, n
+		out[i] = in1[i] * scale1 - in2[i] * scale2
+	}
+end
+
+
+procedure cnvgline1 (line, offset, im, bpm, imdata, bp)
+
+int	line			#I Line to be read
+int	offset[2]		#I Offsets
+pointer	im[2]			#I Image pointers
+pointer	bpm[2]			#I Bad pixel mask pointers
+pointer	imdata[2]		#U Image data
+pointer	bp			#U Bad pixel data
+
+bool	overlap
+int	nl1, nl2, loff, l2
+int	nc1, nc2, nc3, off1, off2, off3, c1, c2
+pointer	imgl2r(), imgl2i()
+
+
+begin
+	# Get data for first image.  Use IMIO buffers except the
+	# bad pixel buffer is not used if there is a second image.
+
+	imdata[1] = imgl2r (im[1], line)
+	if (bpm[1] != NULL) {
+	    if (im[2] == NULL)
+		bp = imgl2i (bpm[1], line)
+	    else
+		call amovi (Memi[imgl2i(bpm[1],line)], Memi[bp],
+		    IM_LEN(bpm[1],1))
+	}
+	if (im[2] == NULL)
+	    return
+
+	# Initialize.
+	if (line == 1) {
+	    nc1 = IM_LEN(im[1],1)
+	    nc2 = IM_LEN(im[2],1)
+	    nl1 = IM_LEN(im[1],2)
+	    nl2 = IM_LEN(im[2],2)
+
+	    overlap = true
+	    if (1-offset[1] < 1 || nc1-offset[1] > nc2)
+		overlap = false
+	    if (1-offset[2] < 1 || nl1-offset[2] > nl2)
+		overlap = false
+
+	    off2 = -offset[1]
+	    c1 = max (1, 1+off2)
+	    c2 = min (nc2, nc1+off2)
+	    nc2 = c2 - c1 + 1
+	    off1 = c1 - off2 - 1
+	    off3 = c2 - off2
+	    off2 = max (0, off2)
+	    nc3 = nc1 - off3
+	    if (off1 > 0) {
+		call aclrr (Memr[imdata[2]], off1)
+		if (bpm[1] == NULL)
+		    call amovki (1, Memi[bp], off1)
+	    }
+	    if (nc3 > 0) {
+		call aclrr (Memr[imdata[2]+off3], nc3)
+		if (bpm[1] == NULL)
+		    call amovki (1, Memi[bp+off3], nc3)
+	    }
+
+	    loff = -offset[2]
+	    if (loff < 0)
+		call aclrr (Memr[imdata[2]], nc1)
+	}
+
+	l2 = line + loff
+	if (l2 < 1 || l2 > nl2) {
+	    call amovki (1, Memi[bp], nc1)
+	    return
+	}
+
+	if (overlap) {
+	    imdata[2] = imgl2r (im[2], l2) + off2
+	    if (bpm[1] != NULL && bpm[2] != NULL)
+		call amaxi (Memi[imgl2i(bpm[2],l2)+off2], Memi[bp], Memi[bp],
+		    nc1)
+	    else if (bpm[2] != NULL)
+		call amovi (Memi[imgl2i(bpm[2],l2)+off2], Memi[bp], nc1)
+	} else {
+	    # Copy the overlapping parts of the second image to the output
+	    # buffers which must be allocated externally.  Use the bad pixel
+	    # mask to flag regions where there is no overlap.
+
+	    call amovr (Memr[imgl2r(im[2],l2)+off2], Memr[imdata[2]+off1], nc2)
+	    if (bpm[1] != NULL && bpm[2] != NULL) {
+		call amaxi (Memi[imgl2i(bpm[2],l2)+off2], Memi[bp+off1],
+		    Memi[bp+off1], nc2)
+		if (off1 > 0)
+		    call amovki (1, Memi[bp], off1)
+		if (nc3 > 0)
+		    call amovki (1, Memi[bp+off3], nc3)
+	    } else if (bpm[2] != NULL)
+		call amovi (Memi[imgl2i(bpm[2],l2)+off2], Memi[bp+off1], nc2)
+	}
+end
+
+
+procedure cnvgline2 (line, offset, im, skymap, sigmap, expmap,
+	skydata, sigdata, expdata)
+
+int	line			#I Line to be read
+int	offset[2]		#I Offsets
+pointer	im[2]			#I Image pointers
+pointer	skymap[2]		#I Sky map
+pointer	sigmap[2]		#I Sky sigma map
+pointer	expmap[2]		#I Exposure map
+pointer	skydata[2]		#U Sky data
+pointer	sigdata[2]		#U Sky sigma data
+pointer	expdata[2]		#U Exposure map data
+
+bool	overlap
+int	nl1, nl2, loff, l2
+int	nc1, nc2, nc3, off1, off2, off3, c1, c2
+pointer	ptr
+
+pointer	map_glr()
+errchk	map_glr
+
+begin
+	# Get data for first image.
+
+	skydata[1] = map_glr (skymap[1], line, READ_ONLY)
+	if (expmap[1] == NULL)
+	    sigdata[1] = map_glr (sigmap[1], line, READ_ONLY)
+	else {
+	    sigdata[1] = map_glr (sigmap[1], line, READ_WRITE)
+	    expdata[1] = map_glr (expmap[1], line, READ_ONLY)
+	    call expsigma (Memr[sigdata[1]], Memr[expdata[1]],
+		IM_LEN(im[1],1), 0)
+	}
+	if (im[2] == NULL)
+	    return
+
+	# Initialize.
+	if (line == 1) {
+	    nc1 = IM_LEN(im[1],1)
+	    nc2 = IM_LEN(im[2],1)
+	    nl1 = IM_LEN(im[1],2)
+	    nl2 = IM_LEN(im[2],2)
+
+	    overlap = true
+	    if (1-offset[1] < 1 || nc1-offset[1] > nc2)
+		overlap = false
+	    if (1-offset[2] < 1 || nl1-offset[2] > nl2)
+		overlap = false
+
+	    off2 = -offset[1]
+	    c1 = max (1, 1+off2)
+	    c2 = min (nc2, nc1+off2)
+	    nc2 = c2 - c1 + 1
+	    off1 = c1 - off2 - 1
+	    off3 = c2 - off2
+	    nc3 = nc1 - off3
+	    if (off1 > 0) {
+		call aclrr (Memr[skydata[2]], off1)
+		call aclrr (Memr[sigdata[2]], off1)
+		if (expmap[2] != NULL)
+		    call aclrr (Memr[expdata[2]], off1)
+	    }
+	    if (nc3 > 0) {
+		call aclrr (Memr[skydata[2]+off3], nc3)
+		call aclrr (Memr[sigdata[2]+off3], nc3)
+		if (expmap[2] != NULL)
+		    call aclrr (Memr[expdata[2]+off3], nc3)
+	    }
+
+	    loff = -offset[2]
+	    if (loff < 0) {
+		call aclrr (Memr[skydata[2]], nc1)
+		call aclrr (Memr[sigdata[2]], nc1)
+		if (expmap[2] != NULL)
+		    call aclrr (Memr[expdata[2]], nc1)
+	    }
+	}
+
+	l2 = line + loff
+	if (l2 < 1 || l2 > nl2)
+	    return
+
+	if (overlap) {
+	    skydata[2] = map_glr (skymap[2], l2, READ_ONLY) + off2
+	    if (expmap[2] == NULL)
+		sigdata[2] = map_glr (sigmap[2], l2, READ_ONLY) + off2
+	    else {
+		sigdata[2] = map_glr (sigmap[2], l2, READ_WRITE) + off2
+		expdata[2] = map_glr (expmap[2], l2, READ_ONLY) + off2
+		call expsigma (Memr[sigdata[2]], Memr[expdata[2]], nc2, 0)
+	    }
+	} else {
+	    # Copy the overlapping parts of the second image to the output
+	    # buffers which must be allocated externally.
+
+	    ptr = map_glr(skymap[2],l2,READ_ONLY)
+	    call amovr (Memr[ptr+off2], Memr[skydata[2]+off1], nc2)
+	    ptr = map_glr(sigmap[2],l2,READ_ONLY)
+	    call amovr (Memr[ptr+off2], Memr[sigdata[2]+off1], nc2)
+	    if (expmap[2] != NULL) {
+		ptr = map_glr(expmap[2],l2,READ_ONLY)
+		call amovr (Memr[ptr+off2], Memr[expdata[2]+off1], nc2)
+		call expsigma (Memr[sigdata[2]], Memr[expdata[2]], nc2, 0)
+	    }
+	}
+end
+
+
+# CNVPARSE -- Parse convolution string.
+
+procedure cnvparse (cnvstr, kernel, nx, ny, logfd)
+
+char	cnvstr[ARB]		#I Convolution string
+pointer	kernel			#O Pointer to convolution kernel elements
+int	nx, ny			#O Convolution size
+int	logfd			#I Log file descriptor
+
+int	i, j, nx2, ny2
+int	ip, fd, open(), fscan(), nscan(), ctor(), ctoi(), strncmp()
+real	val, sx, sy
+pointer	ptr
+errchk	open
+
+define	unknown_	10
+
+begin
+	kernel = NULL
+
+	for (ip=1; IS_WHITE(cnvstr[ip]); ip=ip+1)
+	    ;
+
+	if (cnvstr[ip] == EOS) {
+	    nx = 1
+	    ny = 1
+	    call malloc (kernel, 1, TY_REAL)
+	    Memr[kernel] = 1
+	} else if (cnvstr[ip] == '@') {
+	    fd = open (cnvstr[ip+1], READ_ONLY, TEXT_FILE)
+	    call malloc (kernel, 100, TY_REAL)
+	    i = 0
+	    nx = 0
+	    ny = 0
+	    while (fscan (fd) != EOF) {
+		do j = 1, ARB {
+		    call gargr (val)
+		    if (nscan() < j)
+			break
+		    Memr[kernel+i] = val
+		    i = i + 1
+		    if (mod (i, 100) == 0)
+			call realloc (kernel, i+100, TY_REAL)
+		}
+		j = j - 1
+		if (nx == 0)
+		    nx = j
+		else if (j != nx) {
+		    call close (fd)
+		    call error (1,
+			"Number of convolution elements inconsistent")
+		}
+		ny = ny + 1
+	    }
+	    call close (fd)
+	} else if (IS_ALPHA(cnvstr[ip])) {
+	    if (strncmp ("block", cnvstr[ip], 5) == 0) {
+		i = 6
+		if (ctoi (cnvstr[ip], i, nx) == 0 ||
+		    ctoi (cnvstr[ip], i, ny) == 0)
+		    goto unknown_
+		call malloc (kernel, nx*ny, TY_REAL)
+		call amovkr (1., Memr[kernel], nx*ny)
+	    } else if (strncmp ("bilinear", cnvstr[ip], 8) == 0) {
+		i = 9
+		if (ctoi (cnvstr[ip], i, nx) == 0 ||
+		    ctoi (cnvstr[ip], i, ny) == 0)
+		    goto unknown_
+		call malloc (kernel, nx*ny, TY_REAL)
+
+		nx2 = nx / 2
+		ny2 = ny / 2
+		ptr = kernel
+		do j = 0, ny-1 {
+		    do i = 0, nx-1 {
+			Memr[ptr] = (nx2-abs(nx2-i)+1) * (ny2-abs(ny2-j)+1)
+			ptr = ptr + 1
+		    }
+		}
+	    } else if (strncmp ("gauss", cnvstr[ip], 5) == 0) {
+		i = 6
+		if (ctoi (cnvstr[ip], i, nx) == 0 ||
+		    ctoi (cnvstr[ip], i, ny) == 0)
+		    goto unknown_
+		if (ctor (cnvstr[ip], i, sx) == 0 ||
+		    ctor (cnvstr[ip], i, sy) == 0)
+		    goto unknown_
+		call malloc (kernel, nx*ny, TY_REAL)
+
+		nx2 = nx / 2
+		ny2 = ny / 2
+		val = 2 * sx * sy
+		ptr = kernel
+		do j = 0, ny-1 {
+		    do i = 0, nx-1 {
+			Memr[ptr] = exp (-((i-nx2)**2+(j-ny2)**2) / val)
+			ptr = ptr + 1
+		    }
+		}
+	    }
+	} else {
+	    call malloc (kernel, 100, TY_REAL)
+	    i = 0
+	    nx = 0
+	    ny = 0
+	    while (cnvstr[ip] != EOS) {
+		do j = 1, ARB {
+		    if (ctor (cnvstr, ip, val) == 0)
+			break
+		    Memr[kernel+i] = val
+		    i = i + 1
+		    if (mod (i, 100) == 0)
+			call realloc (kernel, i+100, TY_REAL)
+		}
+		j = j - 1
+		if (nx == 0)
+		    nx = j
+		else if (j != nx)
+		    call error (1,
+			"Number of convolution elements inconsistent")
+		ny = ny + 1
+		if (cnvstr[ip] != EOS)
+		    ip = ip + 1
+		for (; IS_WHITE(cnvstr[ip]); ip=ip+1)
+		    ;
+	    }
+	}
+
+	if (kernel == NULL)
+unknown_    call error (1, "Unrecognized convolution")
+
+	if (mod (nx, 2) != 1 || mod (ny, 2) != 1) {
+	    call mfree (kernel, TY_REAL)
+	    call error (1, "Convolution size must be odd")
+	}
+
+	if (logfd != NULL) {
+	    ptr = kernel
+	    call eprintf ("    Convolution:\n")
+	    do j = 1, ny {
+		call eprintf ("     ")
+		do i = 1, nx {
+		    call eprintf (" %7.3g")
+			call pargr (Memr[ptr])
+		    ptr = ptr + 1
+		}
+		call eprintf ("\n")
+	    }
+	}
+
+end