Initial commit

1 files changed, 320 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/daofind/apconvolve.x b/noao/digiphot/apphot/daofind/apconvolve.x
new file mode 100644
index 00000000..0c6f58aa
--- /dev/null
+++ b/noao/digiphot/apphot/daofind/apconvolve.x
@@ -0,0 +1,320 @@
+include <imhdr.h>
+include <imset.h>
+include "../lib/find.h"
+
+# AP_FCONVOLVE --  Solve for the density enhancement image and optionally
+# the sky enhancement image in the case where datamin and datamax are not
+# defined.
+
+procedure ap_fconvolve (im, den, sky, kernel1, kernel2, skip, nxk, nyk, const2)
+
+pointer	im			# pointer to the input image
+pointer	den			# pointer to the output density image
+pointer	sky			# pointer to the output sky image
+real	kernel1[nxk,nyk]	# the first convolution kernel
+real	kernel2[nxk,nyk]	# the second convolution kernel
+int	skip[nxk,nyk]		# the skip array
+int	nxk, nyk		# dimensions of the kernel
+real	const2			# subtraction constant for the skyimage
+
+int	i, ncols, nlines, col1, col2, inline, outline
+pointer	sp, lineptrs, outbuf1, outbuf2
+pointer	imgs2r(), impl2r()
+errchk	imgs2r, impl2r, imflush
+
+begin
+	# Set up an array of linepointers.
+	call smark (sp)
+	call salloc (lineptrs, nyk, TY_POINTER)
+
+	# Set the number of image buffers.
+	call imseti (im, IM_NBUFS, nyk)
+
+	ncols = IM_LEN(den,1)
+	nlines = IM_LEN(den,2)
+
+	# Set input image column limits.
+	col1 = 1 - nxk / 2
+	col2 = IM_LEN(im,1) + nxk / 2
+
+	# Initialise the line buffers.
+	inline = 1 - nyk / 2
+	do i = 1 , nyk - 1 {
+	    Memi[lineptrs+i] = imgs2r (im, col1, col2, inline, inline)
+	    inline = inline + 1
+	}
+
+	# Generate the output image line by line.
+	do outline = 1, nlines {
+
+	    # Scroll the input buffers.
+	    do i = 1, nyk - 1
+		Memi[lineptrs+i-1] = Memi[lineptrs+i]
+
+	    # Read in new image line.
+	    Memi[lineptrs+nyk-1] = imgs2r (im, col1, col2, inline,
+	        inline)
+
+	    # Get first output image line.
+	    outbuf1 = impl2r (den, outline)
+	    if (outbuf1 == EOF)
+		call error (0, "Error writing first output image.")
+
+	    # Generate first output image line.
+	    call aclrr (Memr[outbuf1], ncols)
+	    do i = 1, nyk
+		call ap_skcnvr (Memr[Memi[lineptrs+i-1]], Memr[outbuf1],
+		    ncols, kernel1[1,i], skip[1,i], nxk)
+
+	    if (sky != NULL) {
+
+	        # Get second output image line.
+	        outbuf2 = impl2r (sky, outline)
+	        if (outbuf2 == EOF)
+		    call error (0, "Error writing second output image.")
+
+	        # Generate second output image line.
+	        call aclrr (Memr[outbuf2], ncols)
+	        do i = 1, nyk
+		    call ap_skcnvr (Memr[Memi[lineptrs+i-1]], Memr[outbuf2],
+		        ncols, kernel2[1,i], skip[1,i], nxk)
+		call ap_w1sur (Memr[outbuf2], Memr[outbuf1], Memr[outbuf2],
+		    ncols, -const2)
+	    }
+
+	    inline = inline + 1
+	}
+
+	# Flush the output image(s).
+	call imflush (den)
+	if (sky != NULL)
+	    call imflush (sky)
+
+	# Free the image buffer pointers.
+	call sfree (sp)
+end
+
+
+# AP_GCONVOLVE --  Solve for the density enhancement image and optionally
+# the sky enhancement image in the case where datamin and datamax are defined.
+
+procedure ap_gconvolve (im, den, sky, kernel1, skip, nxk, nyk, gsums,
+	datamin, datamax)
+
+pointer	im			# pointer to the input image
+pointer	den			# pointer to the output density image
+pointer	sky			# pointer to the output sky image
+real	kernel1[nxk,nyk]	# the first convolution kernel
+int	skip[nxk,nyk]		# the sky array
+int	nxk, nyk		# dimensions of the kernel
+real	gsums[ARB]		# array of kernel sums
+real	datamin, datamax	# the good data minimum and maximum
+
+int	i, ncols, nlines, col1, col2, inline, outline
+pointer	sp, lineptrs, sd, sgd, sg, sgsq, p, outbuf2
+pointer	imgs2r(), impl2r()
+errchk	imgs2r, impl2r, imflush
+
+begin
+	# Set up an array of linepointers.
+	call smark (sp)
+	call salloc (lineptrs, nyk, TY_POINTER)
+
+	# Set the number of image buffers.
+	call imseti (im, IM_NBUFS, nyk)
+
+	ncols = IM_LEN(den,1)
+	nlines = IM_LEN(den,2)
+
+	# Allocate some working space.
+	call salloc (sd, ncols, TY_REAL)
+	call salloc (sgsq, ncols, TY_REAL)
+	call salloc (sg, ncols, TY_REAL)
+	call salloc (p, ncols, TY_REAL)
+
+	# Set input image column limits.
+	col1 = 1 - nxk / 2
+	col2 = IM_LEN(im,1) + nxk / 2
+
+	# Initialise the line buffers.
+	inline = 1 - nyk / 2
+	do i = 1 , nyk - 1 {
+	    Memi[lineptrs+i] = imgs2r (im, col1, col2, inline, inline)
+	    inline = inline + 1
+	}
+
+	# Generate the output image line by line.
+	do outline = 1, nlines {
+
+	    # Scroll the input buffers.
+	    do i = 1, nyk - 1
+		Memi[lineptrs+i-1] = Memi[lineptrs+i]
+
+	    # Read in new image line.
+	    Memi[lineptrs+nyk-1] = imgs2r (im, col1, col2, inline,
+	        inline)
+
+	    # Get first output image line.
+	    sgd = impl2r (den, outline)
+	    if (sgd == EOF)
+		call error (0, "Error writing first output image.")
+
+	    # Generate first output image line.
+	    call aclrr (Memr[sgd], ncols)
+	    call aclrr (Memr[sd], ncols)
+	    call amovkr (gsums[GAUSS_SUMG], Memr[sg], ncols)
+	    call amovkr (gsums[GAUSS_SUMGSQ], Memr[sgsq], ncols)
+	    call amovkr (gsums[GAUSS_PIXELS], Memr[p], ncols)
+	    do i = 1, nyk
+		call ap_gdsum (Memr[Memi[lineptrs+i-1]], Memr[sgd], Memr[sd],
+		    Memr[sg], Memr[sgsq], Memr[p], ncols, kernel1[1,i],
+		    skip[1,i], nxk, datamin, datamax)
+	    call ap_gdavg (Memr[sgd], Memr[sd], Memr[sg], Memr[sgsq],
+	        Memr[p], ncols, gsums[GAUSS_PIXELS], gsums[GAUSS_DENOM],
+		gsums[GAUSS_SGOP])
+
+	    if (sky != NULL) {
+
+	        # Get second output image line.
+	        outbuf2 = impl2r (sky, outline)
+	        if (outbuf2 == EOF)
+		    call error (0, "Error writing second output image.")
+
+	        # Generate second output image line.
+		call ap_davg (Memr[sd], Memr[sgd], Memr[sg], Memr[p],
+		    Memr[outbuf2], ncols)
+	    }
+
+	    inline = inline + 1
+	}
+
+	# Flush the output image(s).
+	call imflush (den)
+	if (sky != NULL)
+	    call imflush (sky)
+
+	# Free the image buffer pointers.
+	call sfree (sp)
+end
+
+
+# AP_SKCNVR -- Compute the convolution kernel using a skip array.
+
+procedure ap_skcnvr (in, out, npix, kernel, skip, nk)
+
+real	in[npix+nk-1]		# the input vector
+real	out[npix]		# the output vector
+int	npix			# the size of the vector
+real	kernel[ARB]		# the convolution kernel
+int	skip[ARB]		# the skip array
+int	nk			# the size of the convolution kernel
+
+int	i, j
+real	sum
+
+begin
+	do i = 1, npix {
+	    sum = out[i]
+	    do j = 1, nk {
+		if (skip[j] == YES)
+		    next
+		sum = sum + in[i+j-1] * kernel[j]
+	    }
+	    out[i] = sum
+	}
+end
+
+
+# AP_GDSUM -- Compute the vector sums required to do the convolution.
+
+procedure  ap_gdsum (in, sgd, sd, sg, sgsq, p, npix, kernel, skip, nk,
+	datamin, datamax)
+
+real	in[npix+nk-1]		# the input vector
+real	sgd[ARB]		# the computed input/output convolution vector
+real	sd[ARB]			# the computed input/output sum vector
+real	sg[ARB]			# the input/ouput first normalization factor
+real	sgsq[ARB]		# the input/ouput second normalization factor
+real	p[ARB]			# the number of points vector
+int	npix			# the size of the vector
+real	kernel[ARB]		# the convolution kernel
+int	skip[ARB]		# the skip array
+int	nk			# the size of the convolution kernel
+real	datamin, datamax	# the good data limits.
+
+int	i, j
+real	data
+
+begin
+	do i = 1, npix {
+	    do j = 1, nk {
+		if (skip[j] == YES)
+		    next
+		data = in[i+j-1]
+		if (data < datamin || data > datamax) {
+		    sgsq[i] = sgsq[i] - kernel[j] ** 2
+		    sg[i] = sg[i] - kernel[j]
+		    p[i] = p[i] - 1.0
+		} else {
+		    sgd[i] = sgd[i] + kernel[j] * data 
+		    sd[i] = sd[i] + data
+		}
+	    }
+	}
+end
+
+
+# AP_GDAVG -- Compute the vector averages required to do the convolution.
+
+procedure  ap_gdavg (sgd, sd, sg, sgsq, p, npix, pixels, denom, sgop)
+
+real	sgd[ARB]		# the computed input/output convolution vector
+real	sd[ARB]			# the computed input/output sum vector
+real	sg[ARB]			# the input/ouput first normalization factor
+real	sgsq[ARB]		# the input/ouput second normalization factor
+real	p[ARB]			# the number of points vector
+int	npix			# the size of the vector
+real	pixels			# number of pixels
+real	denom			# kernel normalization factor
+real	sgop			# kernel normalization factor
+
+int	i
+
+begin
+	do i = 1, npix {
+	    if (p[i] > 1.5) {
+		if (p[i] < pixels) {
+		    sgsq[i] = sgsq[i] - (sg[i] ** 2) / p[i]
+		    if (sgsq[i] != 0.0)
+			sgd[i] = (sgd[i] - sg[i] * sd[i] / p[i]) / sgsq[i]
+		    else
+			sgd[i] = 0.0
+		} else
+		    sgd[i] = (sgd[i] - sgop * sd[i]) / denom
+	    } else
+		sgd[i] = 0.0
+	}
+end
+
+
+# AP_DAVG -- Generate the results the optional sky output image.
+
+procedure ap_davg (sd, sgd, sg, p, out, npix)
+
+real	sd[ARB]			# the computed input/output sum vector
+real	sgd[ARB]		# the computed input/output convolution vector
+real	sg[ARB]			# the input/ouput first normalization factor
+real	p[ARB]			# the number of points vector
+real	out[ARB]		# the output array
+int	npix			# the size of the vector
+
+int	i
+
+begin
+	do i = 1, npix {
+	    if (p[i] > 0.0)
+	        out[i] = (sd[i] - sgd[i] * sg[i]) / p[i]
+	    else
+		out[i] = 0.0
+	}
+end