Initial commit

1 files changed, 139 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/daofind/apbfdfind.x b/noao/digiphot/apphot/daofind/apbfdfind.x
new file mode 100644
index 00000000..acf34ae4
--- /dev/null
+++ b/noao/digiphot/apphot/daofind/apbfdfind.x
@@ -0,0 +1,139 @@
+include <imhdr.h>
+include <imio.h>
+include <mach.h>
+include "../lib/apphot.h"
+include "../lib/noise.h"
+include "../lib/find.h"
+
+# AP_BFDFIND -- Find stars in an image using a pattern matching
+# technique.
+
+procedure ap_bfdfind (im, cnv, sky, out, ap, boundary, constant, verbose)
+
+pointer	im		# pointer to the input image
+pointer	cnv		# pointer to the convolved image
+pointer	sky		# pointer to the sky image
+int	out		# the output file descriptor
+pointer	ap		# pointer to the apphot structure
+int	boundary	# type of boundary extension
+real	constant	# constant for constant boundary extension
+int	verbose		# verbose switch
+
+int	norm, nxk, nyk, nstars, stid
+pointer	sp, gker2d, ngker2d, dker2d, skip
+real	a, b, c, f, gsums[LEN_GAUSS], skysigma, skymode, threshold, relerr
+real	dmax, dmin, xsigsq, ysigsq
+int	apstati(), ap_find()
+real	ap_egkernel(), apstatr()
+
+begin
+	# Compute the parameters of the Gaussian kernel.
+	call ap_egparams (FWHM_TO_SIGMA * apstatr (ap, FWHMPSF) *
+	    apstatr (ap, SCALE), apstatr (ap, RATIO), apstatr (ap, THETA),
+	    apstatr (ap, NSIGMA), a, b, c, f, nxk, nyk)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (gker2d, nxk * nyk, TY_REAL)
+	call salloc (ngker2d, nxk * nyk, TY_REAL)
+	call salloc (dker2d, nxk * nyk, TY_REAL)
+	call salloc (skip, nxk * nyk, TY_INT)
+
+	# Compute the 1 and 2 D kernels.
+	if (IS_INDEFR(apstatr(ap, DATAMIN)) && IS_INDEFR(apstatr(ap,
+	    DATAMAX))) {
+	    norm = YES
+	    dmin = -MAX_REAL
+	    dmax = MAX_REAL
+	} else {
+	    norm = NO
+	    if (IS_INDEFR(apstatr (ap, DATAMIN)))
+		dmin = -MAX_REAL
+	    else
+		dmin = apstatr (ap, DATAMIN)
+	    if (IS_INDEFR(apstatr (ap, DATAMAX)))
+		dmax = MAX_REAL
+	    else
+		dmax = apstatr (ap, DATAMAX)
+	}
+	relerr = ap_egkernel (Memr[gker2d], Memr[ngker2d], Memr[dker2d],
+	    Memi[skip], nxk, nyk, gsums, a, b, c, f)
+
+	# Set up the image boundary extension characteristics.
+	call ap_imset (im, boundary, max (1 + nxk / 2, 1 + nyk / 2),
+	    constant)
+	call ap_imset (cnv, boundary, max (1 + nxk / 2, 1 + nyk / 2),
+	    constant)
+
+	# Convolve the input image with the Gaussian kernel. The resultant
+	# picture constains in each pixel the height of the Gaussian
+	# function centered in the subarray which best represents the data
+	# within a circle of nsigma * sigma of the Gaussian.
+
+	if (IM_ACMODE(cnv) != READ_ONLY) {
+	    if (norm == YES)
+	        call ap_fconvolve (im, cnv, sky, Memr[ngker2d], Memr[dker2d],
+		    Memi[skip], nxk, nyk, gsums[GAUSS_SGOP])
+	    else
+	        call ap_gconvolve (im, cnv, sky, Memr[gker2d], Memi[skip],
+		    nxk, nyk, gsums, dmin, dmax)
+	}
+
+	# Save the task parameters in the database file if the savepars
+	# switch is enabled, otherwise a simple list of detected objects
+	# is written to the data base file.
+
+	call ap_wfdparam (out, ap)
+
+	# Find all the objects in the input image with the specified image
+	# characteristics.
+
+	if (verbose == YES) {
+	    call printf ("\nImage: %s  ")
+		call pargstr (IM_HDRFILE(im))
+	    call printf ("fwhmpsf: %g  ratio: %g  theta: %g  nsigma: %g\n\n")
+		call pargr (apstatr (ap, FWHMPSF))
+		call pargr (apstatr (ap, RATIO))
+		call pargr (apstatr (ap, THETA))
+		call pargr (apstatr (ap, NSIGMA))
+	}
+
+	# Get the skymode and threshold.
+	skysigma = apstatr (ap, SKYSIGMA)
+	if (IS_INDEFR(skysigma)) {
+	    skymode = 0.0
+	    threshold = 0.0
+	} else {
+	    skymode = apstatr (ap, EPADU) * (skysigma ** 2 -
+	        (apstatr (ap, READNOISE) / apstatr (ap, EPADU)) ** 2)
+	    skymode = max (0.0, skymode)
+	    threshold = apstatr (ap, THRESHOLD) * skysigma
+	}
+
+	# Compute the x and y sigma.
+	xsigsq = (apstatr (ap, SCALE) * apstatr (ap, FWHMPSF) / 2.35482) ** 2
+	ysigsq = (apstatr (ap, SCALE) * apstatr (ap, RATIO) *
+	    apstatr (ap, FWHMPSF) / 2.35482) ** 2
+
+	# Find the stars.
+	stid = 1
+	nstars = ap_find (ap, im, cnv, out, NULL, Memr[gker2d],
+	    Memi[skip], nxk, nyk, skymode, threshold, relerr,
+	    apstati (ap,POSITIVE), xsigsq, ysigsq, dmin, dmax,
+	    apstatr (ap, SHARPLO), apstatr (ap, SHARPHI), apstatr (ap,
+	    ROUNDLO), apstatr (ap, ROUNDHI), verbose, stid, NO)
+
+	if (verbose == YES) {
+	    call printf (
+	        "\nthreshold: %g relerr: %5.3f  %g <= sharp <= %g  ")
+		call pargr (threshold)
+		call pargr (relerr)
+		call pargr (apstatr (ap, SHARPLO))
+		call pargr (apstatr (ap, SHARPHI))
+	    call printf ("%g <= round <= %g \n\n")
+		call pargr (apstatr (ap, ROUNDLO))
+		call pargr (apstatr (ap, ROUNDHI))
+	}
+
+	call sfree (sp)
+end