diff options
Diffstat (limited to 'noao/digiphot/apphot/daofind')
19 files changed, 2715 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 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 diff --git a/noao/digiphot/apphot/daofind/apegkernel.x b/noao/digiphot/apphot/daofind/apegkernel.x new file mode 100644 index 00000000..3023098d --- /dev/null +++ b/noao/digiphot/apphot/daofind/apegkernel.x @@ -0,0 +1,132 @@ +include <math.h> +include "../lib/find.h" + +# Set up the gaussian fitting structure. + +# AP_EGPARAMS -- Calculate the parameters of the elliptical Gaussian needed +# to compute the kernel. + +procedure ap_egparams (sigma, ratio, theta, nsigma, a, b, c, f, nx, ny) + +real sigma # sigma of Gaussian in x +real ratio # Ratio of half-width in y to x +real theta # position angle of Gaussian +real nsigma # limit of convolution +real a, b, c, f # ellipse parameters +int nx, ny # dimensions of the kernel + +real sx2, sy2, cost, sint, discrim +bool fp_equalr () + +begin + # Define some temporary variables. + sx2 = sigma ** 2 + sy2 = (ratio * sigma) ** 2 + cost = cos (DEGTORAD (theta)) + sint = sin (DEGTORAD (theta)) + + # Compute the ellipse parameters. + if (fp_equalr (ratio, 0.0)) { + if (fp_equalr (theta, 0.0) || fp_equalr (theta, 180.)) { + a = 1. / sx2 + b = 0.0 + c = 0.0 + } else if (fp_equalr (theta, 90.0)) { + a = 0.0 + b = 0.0 + c = 1. / sx2 + } else + call error (0, "AP_EGPARAMS: Cannot make 1D Gaussian.") + f = nsigma ** 2 / 2. + nx = 2 * int (max (sigma * nsigma * abs (cost), RMIN)) + 1 + ny = 2 * int (max (sigma * nsigma * abs (sint), RMIN)) + 1 + } else { + a = cost ** 2 / sx2 + sint ** 2 / sy2 + b = 2. * (1.0 / sx2 - 1.0 / sy2) * cost * sint + c = sint ** 2 / sx2 + cost ** 2 / sy2 + discrim = b ** 2 - 4. * a * c + f = nsigma ** 2 / 2. + nx = 2 * int (max (sqrt (-8. * c * f / discrim), RMIN)) + 1 + ny = 2 * int (max (sqrt (-8. * a * f / discrim), RMIN)) + 1 + } +end + + +# AP_EGKERNEL -- Compute the elliptical Gaussian kernel. + +real procedure ap_egkernel (gkernel, ngkernel, dkernel, skip, nx, ny, gsums, a, + b, c, f) + +real gkernel[nx,ny] # output Gaussian amplitude kernel +real ngkernel[nx,ny] # output normalized Gaussian amplitude kernel +real dkernel[nx,ny] # output Gaussian sky kernel +int skip[nx,ny] # output skip subraster +int nx, ny # input dimensions of the kernel +real gsums[ARB] # output array of gsums +real a, b, c, f # ellipse parameters + +int i, j, x0, y0, x, y +real npts, rjsq, rsq, relerr,ef + +begin + # Initialize. + x0 = nx / 2 + 1 + y0 = ny / 2 + 1 + gsums[GAUSS_SUMG] = 0.0 + gsums[GAUSS_SUMGSQ] = 0.0 + npts = 0.0 + + # Compute the kernel and principal sums. + do j = 1, ny { + y = j - y0 + rjsq = y ** 2 + do i = 1, nx { + x = i - x0 + rsq = sqrt (x ** 2 + rjsq) + ef = 0.5 * (a * x ** 2 + c * y ** 2 + b * x * y) + gkernel[i,j] = exp (-1.0 * ef) + if (ef <= f || rsq <= RMIN) { + #gkernel[i,j] = exp (-ef) + ngkernel[i,j] = gkernel[i,j] + dkernel[i,j] = 1.0 + gsums[GAUSS_SUMG] = gsums[GAUSS_SUMG] + gkernel[i,j] + gsums[GAUSS_SUMGSQ] = gsums[GAUSS_SUMGSQ] + + gkernel[i,j] ** 2 + skip[i,j] = NO + npts = npts + 1.0 + } else { + #gkernel[i,j] = 0.0 + ngkernel[i,j] = 0.0 + dkernel[i,j] = 0.0 + skip[i,j] = YES + } + } + } + + # Store the remaining sums. + gsums[GAUSS_PIXELS] = npts + gsums[GAUSS_DENOM] = gsums[GAUSS_SUMGSQ] - gsums[GAUSS_SUMG] ** 2 / + npts + gsums[GAUSS_SGOP] = gsums[GAUSS_SUMG] / npts + + # Normalize the amplitude kernel. + do j = 1, ny { + do i = 1, nx { + if (skip[i,j] == NO) + ngkernel[i,j] = (gkernel[i,j] - gsums[GAUSS_SGOP]) / + gsums[GAUSS_DENOM] + } + } + + # Normalize the sky kernel + do j = 1, ny { + do i = 1, nx { + if (skip[i,j] == NO) + dkernel[i,j] = dkernel[i,j] / npts + } + } + + relerr = 1.0 / gsums[GAUSS_DENOM] + + return (sqrt (relerr)) +end diff --git a/noao/digiphot/apphot/daofind/apfdcolon.x b/noao/digiphot/apphot/daofind/apfdcolon.x new file mode 100644 index 00000000..cb618bb9 --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfdcolon.x @@ -0,0 +1,282 @@ +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/noise.h" +include "../lib/find.h" + +# AP_FDCOLON -- Process colon commands from the daofind task. + +procedure ap_fdcolon (ap, im, out, stid, cmdstr, newimage, newbuf, newfit) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the iraf image +int out # output file descriptor +int stid # output file sequence number +char cmdstr # command string +int newimage # new mage ? +int newbuf # new center buffer ? +int newfit # new center fit ? + +int cl, junk +pointer sp, incmd, outcmd +int strdic() + +begin + call smark (sp) + call salloc (incmd, SZ_LINE, TY_CHAR) + call salloc (outcmd, SZ_LINE, TY_CHAR) + + # Get the command. + call sscan (cmdstr) + call gargwrd (Memc[incmd], SZ_LINE) + if (Memc[incmd] == EOS) { + call sfree (sp) + return + } + + # Process the command making sure that the pointer to the + # coords file is always NULL. + if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, APCMDS) != 0) { + cl = NULL + call ap_apcolon (ap, im, cl, out, stid, junk, cmdstr, newimage, + junk, junk, junk, junk, newbuf, newfit) + if (cl != NULL) { + call close (cl) + cl = NULL + call apsets (ap, CLNAME, "") + } + } else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, NCMDS) != 0) { + call ap_nscolon (ap, im, out, stid, cmdstr, junk, junk, + junk, junk, newbuf, newfit) + } else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, FCMDS) != 0) { + call ap_fcolon (ap, out, stid, cmdstr, newbuf, newfit) + } else { + call ap_fimcolon (ap, cmdstr) + } + + call sfree (sp) +end + + +# AP_FCOLON -- Process colon commands for setting the find algorithm +# parameters. + +procedure ap_fcolon (ap, out, stid, cmdstr, newbuf, newfit) + +pointer ap # pointer to the apphot structure +int out # output file descriptor +int stid # file number id +char cmdstr[ARB] # command string +int newbuf, newfit # change magnitude parameters + +bool bval +int ncmd +pointer sp, cmd, str +real rval + +bool itob() +int strdic(), nscan(), btoi(), apstati() +real apstatr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Get the command. + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call sfree (sp) + return + } + + # Process the command. + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, FCMDS) + switch (ncmd) { + + case FCMD_NSIGMA: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_NSIGMA) + call pargr (apstatr (ap, NSIGMA)) + call pargstr (UN_FSIGMA) + } else { + call apsetr (ap, NSIGMA, rval) + if (stid > 1) + call ap_rparam (out, KY_NSIGMA, rval, UN_FSIGMA, + "size of kernel in sigma") + newbuf = YES; newfit = YES + } + + case FCMD_RATIO: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_RATIO) + call pargr (apstatr (ap, RATIO)) + call pargstr (UN_FNUMBER) + } else { + call apsetr (ap, RATIO, rval) + if (stid > 1) + call ap_rparam (out, KY_RATIO, rval, UN_FNUMBER, + "sigma y / x of Gaussian kernel") + newbuf = YES; newfit = YES + } + + case FCMD_SHARPLO: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_SHARPLO) + call pargr (apstatr (ap, SHARPLO)) + call pargstr (UN_FNUMBER) + } else { + call apsetr (ap, SHARPLO, rval) + if (stid > 1) + call ap_rparam (out, KY_SHARPLO, rval, UN_FNUMBER, + "lower sharpness bound") + newfit = YES + } + + case FCMD_SHARPHI: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_SHARPHI) + call pargr (apstatr (ap, SHARPHI)) + call pargstr (UN_FNUMBER) + } else { + call apsetr (ap, SHARPHI, rval) + if (stid > 1) + call ap_rparam (out, KY_SHARPHI, rval, UN_FNUMBER, + "upper sharpness bound") + newfit = YES + } + + case FCMD_ROUNDLO: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_ROUNDLO) + call pargr (apstatr (ap, ROUNDLO)) + call pargstr (UN_FNUMBER) + } else { + call apsetr (ap, ROUNDLO, rval) + if (stid > 1) + call ap_rparam (out, KY_ROUNDLO, rval, UN_FNUMBER, + "lower roundness bound") + newfit = YES + } + + case FCMD_ROUNDHI: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_ROUNDHI) + call pargr (apstatr (ap, ROUNDHI)) + call pargstr (UN_FNUMBER) + } else { + call apsetr (ap, ROUNDHI, rval) + if (stid > 1) + call ap_rparam (out, KY_ROUNDHI, rval, UN_FNUMBER, + "upper roundness bound") + newfit = YES + } + + case FCMD_MKDETECTIONS: + call gargb (bval) + if (nscan () == 1) { + call printf ("%s = %b\n") + call pargstr (KY_MKDETECTIONS) + call pargb (itob (apstati (ap, MKDETECTIONS))) + } else + call apseti (ap, MKDETECTIONS, btoi (bval)) + + case FCMD_THETA: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_THETA) + call pargr (apstatr (ap, THETA)) + call pargstr (UN_FDEGREES) + } else { + call apsetr (ap, THETA, rval) + if (stid > 1) + call ap_rparam (out, KY_THETA, rval, UN_FDEGREES, + "position angle") + newbuf = YES; newfit = YES + } + + case FCMD_THRESHOLD: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_THRESHOLD) + call pargr (apstatr (ap, THRESHOLD)) + call pargstr (UN_FSIGMA) + } else { + call apsetr (ap, THRESHOLD, rval) + if (stid > 1) + call ap_rparam (out, KY_THRESHOLD, rval, UN_FSIGMA, + "detection threshold in sigma") + newfit = YES + } + + default: + call printf ("Unknown or ambiguous colon command\7\n") + } + + call sfree (sp) +end + + +# AP_FIMCOLON -- Process colon commands for the daofind task that do +# not affect the data dependent or find parameters. + +procedure ap_fimcolon (ap, cmdstr) + +pointer ap # pointer to the apphot structure +char cmdstr[ARB] # command string + +int ncmd +pointer sp, cmd +int strdic() + +begin + # Get the command. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call sfree (sp) + return + } + + # Process the command. + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, MISC1) + switch (ncmd) { + case ACMD_SHOW: + call gargwrd (Memc[cmd], SZ_LINE) + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, FSHOWARGS) + switch (ncmd) { + case FCMD_DATA: + call printf ("\n") + call ap_nshow (ap) + call printf ("\n") + case FCMD_FIND: + call printf ("\n") + call ap_fshow (ap) + call printf ("\n") + default: + call printf ("\n") + call ap_fdshow (ap) + call printf ("\n") + } + default: + call printf ("Unknown or ambiguous colon command\7\n") + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/daofind/apfdconfirm.x b/noao/digiphot/apphot/daofind/apfdconfirm.x new file mode 100644 index 00000000..ce69cbd9 --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfdconfirm.x @@ -0,0 +1,22 @@ +# AP_FDCONFIRM -- Procedure to confirm the critical daofind parameters. + +procedure ap_fdconfirm (ap) + +pointer ap # pointer to the apphot structure + +real rval +real ap_vfwhmpsf(), ap_vsigma(), ap_vthreshold() +real ap_vdatamin(), ap_vdatamax() + +begin + call printf ("\n") + + # Verify the critical parameters. + rval = ap_vfwhmpsf (ap) + rval = ap_vsigma (ap) + rval = ap_vthreshold (ap) + rval = ap_vdatamin (ap) + rval = ap_vdatamax (ap) + + call printf ("\n") +end diff --git a/noao/digiphot/apphot/daofind/apfdfind.x b/noao/digiphot/apphot/daofind/apfdfind.x new file mode 100644 index 00000000..7778c946 --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfdfind.x @@ -0,0 +1,213 @@ +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/find.h" + +define HELPFILE "apphot$daofind/daofind.key" + +# AP_FDFIND -- Find objects in an image interactively. + +int procedure ap_fdfind (denname, skyname, ap, im, gd, id, out, boundary, + constant, save, skysave, interactive, cache) + +char denname[ARB] # name of density enhancement image +char skyname[ARB] # name of the fitted sky image +pointer ap # pointer to the apphot structure +pointer im # pointer to the IRAF image +pointer gd # pointer to the graphics stream +pointer id # pointer to the image display stream +int out # output file descriptor +int boundary # type of boundary extension +real constant # constatn for constant boundary extension +int save # save convolved image +int skysave # save the sky image +int interactive # interactive mode +int cache # cache the input image pixels + +real wx, wy +pointer sp, cmd, root, den, sky +int wcs, key, newimage, newden, newfit, stid, memstat, req_size, old_size +int buf_size + +real apstatr() +pointer ap_immap() +int clgcur(), apgqverify(), apgtverify(), sizeof(), ap_memstat() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (root, SZ_FNAME, TY_CHAR) + + # Initialize cursor command. + key = ' ' + Memc[cmd] = EOS + call strcpy (" ", Memc[root], SZ_FNAME) + + # Initialize fitting parameters. + den = NULL + sky = NULL + newimage = NO + newden = YES + newfit = YES + memstat = NO + + # Loop over the cursor commands. + stid = 1 + while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) != + EOF) { + + # Store the current cursor coordinates. + call ap_vtol (im, wx, wy, wx, wy, 1) + call apsetr (ap, CWX, wx) + call apsetr (ap, CWY, wy) + + # Process the colon commands. + switch (key) { + + # Quit. + case 'q': + if (interactive == YES && apgqverify ("daofind", + ap, key) == YES) { + call sfree (sp) + if (den != NULL) { + call imunmap (den) + if (save == NO) + call imdelete (denname) + } + if (sky != NULL) + call imunmap (sky) + return (apgtverify (key)) + } else { + if (den != NULL) { + call imunmap (den) + if (save == NO) + call imdelete (denname) + } + if (sky != NULL) + call imunmap (sky) + call sfree (sp) + return (YES) + } + + # Get information on keystroke commands. + case '?': + if ((id != NULL) && (gd == id)) + call gpagefile (id, HELPFILE, "") + else if (interactive == YES) + call pagefile (HELPFILE, "[space=morehelp,q=quit,?=help]") + + # Plot a centered stellar radial profile + case 'd': + if (interactive == YES) + call ap_qrad (ap, im, wx, wy, gd) + + # Interactively set the daofind parameters. + case 'i': + if (interactive == YES) { + call ap_fdradsetup (ap, im, wx, wy, gd, out, stid) + newden = YES + newfit = YES + } + + # Verify the critical daofind parameters. + case 'v': + call ap_fdconfirm (ap) + newden = YES + newfit = YES + + # Save daofind parameters in the pset files. + case 'w': + call ap_fdpars (ap) + + # Process apphot : commands. + case ':': + call ap_fdcolon (ap, im, out, stid, Memc[cmd], newimage, + newden, newfit) + + # Determine the viewport and data window of image display. + if (newimage == YES) { + if ((id != NULL) && (id != gd)) + call ap_gswv (id, Memc[cmd], im, 4) + req_size = MEMFUDGE * (IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + 2 * IM_LEN(im,1) * + IM_LEN(im,2) * sizeof (TY_REAL)) + memstat = ap_memstat (cache, req_size, old_size) + if (memstat == YES) + call ap_pcache (im, INDEFI, buf_size) + } + newimage = NO + + # Find the stars. + case 'f', ' ': + + if (newden == YES) { + + if (den != NULL) { + call imunmap (den) + call imdelete (denname) + } + den = ap_immap (denname, im, ap, save) + if (memstat == YES) + call ap_pcache (den, INDEFI, buf_size) + + if (sky != NULL) + call imunmap (sky) + if (skysave == YES) { + sky = ap_immap (skyname, im, ap, YES) + if (memstat == YES) + call ap_pcache (den, INDEFI, buf_size) + } else + sky = NULL + newden = NO + + if (key == 'f') { + call ap_fdstars (im, ap, den, sky, NULL, id, + boundary, constant, NO, stid) + } else { + call ap_outmap (ap, out, Memc[root]) + call ap_fdstars (im, ap, den, sky, out, id, + boundary, constant, NO, stid) + newfit = NO + } + + } else if (newfit == YES) { + + if (key == 'f') { + call ap_fdstars (im, ap, den, sky, NULL, id, + boundary, constant, YES, stid) + } else { + call ap_outmap (ap, out, Memc[root]) + call ap_fdstars (im, ap, den, sky, out, id, + boundary, constant, YES, stid) + newfit = NO + } + + } else { + call printf ("Detection parameters have not changed\n") + } + + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + + stid = 1 + #newden = NO + #newfit = NO + + default: + # do nothing + call printf ("Unknown or ambiguous keystroke command\n") + } + + # Setup for the next object. + key = ' ' + Memc[cmd] = EOS + call apsetr (ap, WX, apstatr (ap, CWX)) + call apsetr (ap, WY, apstatr (ap, CWY)) + + } + +end diff --git a/noao/digiphot/apphot/daofind/apfdfree.x b/noao/digiphot/apphot/daofind/apfdfree.x new file mode 100644 index 00000000..3f1ea762 --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfdfree.x @@ -0,0 +1,34 @@ +include "../lib/apphotdef.h" + +# AP_FDFREE -- Free the apphot data structure. + +procedure ap_fdfree (ap) + +pointer ap # pointer to the apphot structure + +begin + if (ap == NULL) + return + if (AP_NOISE(ap) != NULL) + call ap_noisecls (ap) + if (AP_PFIND(ap) != NULL) + call ap_fdcls (ap) + if (AP_PDISPLAY(ap) != NULL) + call ap_dispcls (ap) + if (AP_IMBUF(ap) != NULL) + call mfree (AP_IMBUF(ap), TY_REAL) + if (AP_MW(ap) != NULL) + call mw_close (AP_MW(ap)) + call mfree (ap, TY_STRUCT) +end + + +# AP_FDCLS -- Free the find data structure. + +procedure ap_fdcls (ap) + +pointer ap # pointer to the apphot structure + +begin + call mfree (AP_PFIND(ap), TY_STRUCT) +end diff --git a/noao/digiphot/apphot/daofind/apfdgpars.x b/noao/digiphot/apphot/daofind/apfdgpars.x new file mode 100644 index 00000000..95d87d23 --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfdgpars.x @@ -0,0 +1,20 @@ +include "../lib/noise.h" +include "../lib/display.h" + +# AP_FDGPARS -- Open up the apphot data structure and get the daofind input +# parameters. + +procedure ap_fdgpars (ap) + +pointer ap # pointer to the apphot structure + +begin + # Open the apphot structure. + call ap_fdinit (ap, 2.0, AP_NPOISSON) + + # Get the data dependent parameters. + call ap_gdapars (ap) + + # Get the find parameters. + call ap_gfipars (ap) +end diff --git a/noao/digiphot/apphot/daofind/apfdinit.x b/noao/digiphot/apphot/daofind/apfdinit.x new file mode 100644 index 00000000..9c31b494 --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfdinit.x @@ -0,0 +1,59 @@ +include "../lib/apphotdef.h" +include "../lib/finddef.h" + +# AP_FDINIT - Initialize the daofind data structure. + +procedure ap_fdinit (ap, fwhmpsf, noise) + +pointer ap # pointer to the apphot structure +real fwhmpsf # FWHM of the PSF +int noise # noise function + +begin + # Allocate space. + call malloc (ap, LEN_APSTRUCT, TY_STRUCT) + + # Set the default global apphot package parameters. + call ap_defsetup (ap, fwhmpsf) + + # Setup the noise structure. + call ap_noisesetup (ap, noise) + + # Setup the display structure. + call ap_dispsetup (ap) + + # Setup the find structure. + call ap_fdsetup (ap) + + # Unused structures are set to null. + AP_PCENTER(ap) = NULL + AP_PSKY(ap) = NULL + AP_PPSF(ap) = NULL + AP_PPHOT(ap) = NULL + AP_POLY(ap) = NULL + AP_RPROF(ap) = NULL +end + + +# AP_FDSETUP -- Initialize the find structure. + +procedure ap_fdsetup (ap) + +pointer ap # pointer to the apphot strucuture + +pointer fnd + +begin + call malloc (AP_PFIND(ap), LEN_FIND, TY_STRUCT) + fnd = AP_PFIND(ap) + + AP_RATIO(fnd) = DEF_RATIO + AP_THETA(fnd) = DEF_RATIO + AP_NSIGMA(fnd) = DEF_NSIGMA + + AP_THRESHOLD(fnd) = DEF_THRESHOLD + AP_SHARPLO(fnd) = DEF_SHARPLO + AP_SHARPHI(fnd) = DEF_SHARPHI + AP_ROUNDLO(fnd) = DEF_ROUNDLO + AP_ROUNDHI(fnd) = DEF_ROUNDLO +end diff --git a/noao/digiphot/apphot/daofind/apfdpars.x b/noao/digiphot/apphot/daofind/apfdpars.x new file mode 100644 index 00000000..108f205b --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfdpars.x @@ -0,0 +1,16 @@ +include "../lib/display.h" + +# AP_FDPARS -- Write out the current daofind parameters to the current +# parameter files. + +procedure ap_fdpars (ap) + +pointer ap # pointer to apphot structure + +begin + # Write the data dependent parameters. + call ap_dapars (ap) + + # Write the daofind parameters. + call ap_fipars (ap) +end diff --git a/noao/digiphot/apphot/daofind/apfdradsetup.x b/noao/digiphot/apphot/daofind/apfdradsetup.x new file mode 100644 index 00000000..708d2690 --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfdradsetup.x @@ -0,0 +1,82 @@ +define HELPFILE "apphot$daofind/idaofind.key" + +# AP_FDRADSETUP -- Procedure to set up daofind interactively using the radial +# profile plot of a bright star. + +procedure ap_fdradsetup (ap, im, wx, wy, gd, out, stid) + +pointer ap # pointer to apphot structure +pointer im # pointer to the IRAF image +real wx, wy # cursor coordinates +pointer gd # pointer to graphics stream +int out # output file descriptor +int stid # output file sequence number) + +int key, wcs +pointer sp, cmd +real rmin, rmax, imin, imax, xcenter, ycenter, rval +real u1, u2, v1, v2, x1, x2, y1, y2 +int ap_showplot(), clgcur() +real ap_cfwhmpsf(), ap_csigma(), ap_cdatamin(), ap_cdatamax() + +begin + # Check for open graphics stream. + if (gd == NULL) + return + call greactivate (gd, 0) + + # Store the old viewport and window coordinates. + call ggview (gd, u1, u2, v1, v2) + call ggwind (gd, x1, x2, y1, y2) + + # Make the plot. + if (ap_showplot (ap, im, wx, wy, gd, xcenter, ycenter, rmin, rmax, + imin, imax) == ERR) { + call gdeactivate (gd, 0) + return + } + + # Initialize. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + call printf ( + "Waiting for setup menu command (?=help, v=default setup, q=quit):\n") + while (clgcur ("gcommands", xcenter, ycenter, wcs, key, Memc[cmd], + SZ_LINE) != EOF) { + + # Enter the cursor setup loop. + switch (key) { + case 'q': + break + case '?': + call gpagefile (gd, HELPFILE, "") + case 'f': + rval = ap_cfwhmpsf (ap, gd, out, stid, rmin, rmax, imin, imax) + case 's': + rval = ap_csigma (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'l': + rval = ap_cdatamin (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'u': + rval = ap_cdatamax (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'v': + rval = ap_cfwhmpsf (ap, gd, out, stid, rmin, rmax, imin, imax) + rval = ap_csigma (ap, gd, out, stid, rmin, rmax, imin, imax) + default: + call printf ("Unknown or ambiguous keystroke command\007\n") + } + call printf ( + "Waiting for setup menu command (?=help, v=default setup, q=quit):\n") + } + + # Interactive setup is complete. + call printf ( + "Interactive setup is complete. Type w to store parameters.\n") + + # Restore the viewport and window coordinates. + call gsview (gd, u1, u2, v1, v2) + call gswind (gd, x1, x2, y1, y2) + + call gdeactivate (gd, 0) + call sfree (sp) +end diff --git a/noao/digiphot/apphot/daofind/apfdshow.x b/noao/digiphot/apphot/daofind/apfdshow.x new file mode 100644 index 00000000..61be95ca --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfdshow.x @@ -0,0 +1,72 @@ +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/find.h" + +# AP_FDSHOW -- Display the current find parameters. + +procedure ap_fdshow (ap) + +pointer ap # pointer to the apphot strucuture + +begin + call ap_nshow (ap) + call printf ("\n") + call ap_fshow (ap) +end + + +# AP_FSHOW -- Procedure to display the current data parameters. + +procedure ap_fshow (ap) + +pointer ap # pointer to the apphot strucuture + +pointer sp, str +bool itob() +int apstati() +real apstatr() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Set the object charactersitics. + call printf ("Kernel Parameters\n") + call printf (" %s = %g %s %s = %b\n") + call pargstr (KY_FWHMPSF) + call pargr (apstatr (ap, FWHMPSF)) + call pargstr (UN_ASCALEUNIT) + call pargstr (KY_POSITIVE) + call pargb (itob (apstati (ap, POSITIVE))) + + call printf (" %s = %g %s %s = %g %s = %g %s\n") + call pargstr (KY_NSIGMA) + call pargr (apstatr (ap, NSIGMA)) + call pargstr (UN_FSIGMA) + call pargstr (KY_RATIO) + call pargr (apstatr (ap, RATIO)) + call pargstr (KY_THETA) + call pargr (apstatr (ap, THETA)) + call pargstr (UN_FDEGREES) + + # Print the rest of the data dependent parameters. + call printf ("\nDetection Parameters\n") + call printf (" %s = %g %s\n") + call pargstr (KY_THRESHOLD) + call pargr (apstatr (ap, THRESHOLD)) + call pargstr (UN_FSIGMA) + + call printf (" %s = %g %s = %g\n") + call pargstr (KY_SHARPLO) + call pargr (apstatr (ap, SHARPLO)) + call pargstr (KY_SHARPHI) + call pargr (apstatr (ap, SHARPHI)) + + call printf (" %s = %g %s = %g\n") + call pargstr (KY_ROUNDLO) + call pargr (apstatr (ap, ROUNDLO)) + call pargstr (KY_ROUNDHI) + call pargr (apstatr (ap, ROUNDHI)) + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/daofind/apfdstars.x b/noao/digiphot/apphot/daofind/apfdstars.x new file mode 100644 index 00000000..893619f5 --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfdstars.x @@ -0,0 +1,151 @@ +include <imhdr.h> +include <mach.h> +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/display.h" +include "../lib/find.h" + +# AP_FDSTARS -- Find stars in an image using a pattern matching technique. + +procedure ap_fdstars (im, ap, cnv, sky, out, id, boundary, constant, + refit, stid) + +pointer im # pointer to the input image +pointer ap # pointer to the apphot structure +pointer cnv # pointer to the convolved image +pointer sky # pointer to the sky image +int out # the output file descriptor +pointer id # pointer to image display stream +int boundary # type of boundary extension +real constant # constant for constant boundary extension +int refit # detect stars again +int stid # output file sequence number + +int norm, nxk, nyk, nstars +pointer sp, str, gker2d, ngker2d, dker2d, skip +real a, b, c, f, skysigma, skymode, threshold, relerr, gsums[LEN_GAUSS] +real dmin, dmax, xsigsq, ysigsq +int apstati(), ap_find() +real ap_egkernel(), apstatr() +data gker2d/NULL/, ngker2d/NULL/, dker2d/NULL/ skip /NULL/ + +define detect_ 99 + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + if (refit == YES) + goto detect_ + + # 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. + if (gker2d != NULL) + call mfree (gker2d, TY_REAL) + call malloc (gker2d, nxk * nyk, TY_REAL) + if (ngker2d != NULL) + call mfree (ngker2d, TY_REAL) + call malloc (ngker2d, nxk * nyk, TY_REAL) + if (dker2d != NULL) + call mfree (dker2d, TY_REAL) + call malloc (dker2d, nxk * nyk, TY_REAL) + if (skip != NULL) + call mfree (skip, TY_INT) + call malloc (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 (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) + +detect_ + + # Write the output header file. + if (stid <= 1) + call ap_wfdparam (out, ap) + + 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)) + + # Find all the objects in the input image with the specified image + # characteristics. + + 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 + } + xsigsq = (apstatr (ap, SCALE) * apstatr (ap, FWHMPSF) / 2.35482) ** 2 + ysigsq = (apstatr (ap, SCALE) * apstatr (ap, RATIO) * + apstatr (ap, FWHMPSF) / 2.35482) ** 2 + + nstars = ap_find (ap, im, cnv, out, id, 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), YES, stid, apstati (ap, + MKDETECTIONS)) + stid = stid + nstars + + 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 apstats (ap, OUTNAME, Memc[str], SZ_FNAME) + call printf ("Output file: %s\n\n") + call pargstr (Memc[str]) + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/daofind/apfind.x b/noao/digiphot/apphot/daofind/apfind.x new file mode 100644 index 00000000..1bcbdd52 --- /dev/null +++ b/noao/digiphot/apphot/daofind/apfind.x @@ -0,0 +1,626 @@ +include <gset.h> +include <mach.h> +include <imhdr.h> +include "../lib/apphot.h" + +# AP_FIND -- Detect images in the convolved image and then compute image +# characteristics using the original image. + +int procedure ap_find (ap, im, cnv, out, id, ker2d, skip, nxk, nyk, skymode, + threshold, relerr, emission, xsigsq, ysigsq, datamin, datamax, + sharplo, sharphi, roundlo, roundhi, interactive, stid, mkdetections) + +pointer ap # the apphot descriptor +pointer im # pointer to the input image +pointer cnv # pointer to the output image +int out # the output file descriptor +pointer id # pointer to the display stream +real ker2d[nxk,ARB] # 2D Gaussian kernel +int skip[nxk,ARB] # 2D skip kernel +int nxk, nyk # dimensions of the kernel +real skymode # estimate of the sky +real threshold # threshold for image detection +real relerr # the relative error of the convolution kernel +int emission # emission features +real xsigsq, ysigsq # sigma of gaussian in x and y +real datamin, datamax # minimum and maximum good data values +real sharplo, sharphi # sharpness limits +real roundlo,roundhi # roundness parameter limits +int interactive # interactive mode +int stid # sequence number +int mkdetections # mark detections + +int inline, i, j, ncols, col1, col2, line1, line2, index, pos +int xmiddle, ymiddle, nonzero, nobjs, nstars, ntotal +pointer sp, bufptrs, imlbuf, cnvlbuf, imbuf, cnvbuf, cols +pointer satur, sharp, round1, round2, x, y + +int ap_detect(), ap_test(), apstati() +pointer imgs2r() +errchk imgs2r() + +begin + # Set up useful line and column limits. + ncols = IM_LEN(im,1) + nxk - 1 + col1 = 1 - nxk / 2 + col2 = IM_LEN(im,1) + nxk / 2 + line1 = 1 + nyk / 2 + line2 = IM_LEN(im,2) + nyk / 2 + xmiddle = 1 + nxk / 2 + ymiddle = 1 + nyk / 2 + + # Compute find the number of defined elements in the kernel. + nonzero = 0 + #skip[xmiddle,ymiddle] = NO + do j = 1, nyk { + do i = 1, nxk { + if (skip[i,j] == NO) + nonzero = nonzero + 1 + } + } + skip[xmiddle,ymiddle] = YES + nonzero = nonzero - 1 + + # Set up a cylindrical buffers and some working space for + # the detected images. + call smark (sp) + call salloc (bufptrs, nyk, TY_INT) + call salloc (imbuf, nyk * ncols, TY_REAL) + call salloc (cnvbuf, nyk * ncols, TY_REAL) + call salloc (cols, ncols, TY_INT) + call salloc (satur, ncols, TY_INT) + call salloc (sharp, ncols, TY_REAL) + call salloc (round1, ncols, TY_REAL) + call salloc (round2, ncols, TY_REAL) + call salloc (x, ncols, TY_REAL) + call salloc (y, ncols, TY_REAL) + + # Read in the first nyk - 1 lines. + pos = nyk + do inline = 1 - nyk / 2, nyk / 2 { + imlbuf = imgs2r (im, col1, col2, inline, inline) + cnvlbuf = imgs2r (cnv, col1, col2, inline, inline) + if (emission == YES) { + call amovr (Memr[imlbuf], Memr[imbuf+(inline+ymiddle-2)*ncols], + ncols) + call amovr (Memr[cnvlbuf], Memr[cnvbuf+(inline+ymiddle-2)* + ncols], ncols) + } else { + call amulkr (Memr[imlbuf], -1.0, Memr[imbuf+(inline+ymiddle-2)* + ncols], ncols) + call amulkr (Memr[cnvlbuf], -1.0, Memr[cnvbuf+(inline+ymiddle- + 2)* ncols], ncols) + } + Memi[bufptrs+pos-1] = pos - 1 + pos = pos - 1 + } + + # Generate the starlist line by line. + ntotal = 0 + pos = nyk + do inline = line1, line2 { + + # Setup the buffer pointer array. + do j = 2, nyk + Memi[bufptrs+j-2] = Memi[bufptrs+j-1] + Memi[bufptrs+nyk-1] = pos + index = (pos - 1) * ncols + + # Read in new image line. + imlbuf = imgs2r (im, col1, col2, inline, inline) + cnvlbuf = imgs2r (cnv, col1, col2, inline, inline) + + # Copy new lines into cylindrical buffer. + if (emission == YES) { + call amovr (Memr[imlbuf], Memr[imbuf+index], ncols) + call amovr (Memr[cnvlbuf], Memr[cnvbuf+index], ncols) + } else { + call amulkr (Memr[imlbuf], -1.0, Memr[imbuf+index], ncols) + call amulkr (Memr[cnvlbuf], -1.0, Memr[cnvbuf+index], ncols) + } + + # Detect stars in each image line. In order for a given pixel + # to be detected as an image the pixel must be above threshold + # and be greater than any other pixel within nsigma sigma. + + # Increment the cylindrical buffer. + if (mod (pos, nyk) == 0) + pos = 1 + else + pos = pos + 1 + + nobjs = ap_detect (Memr[cnvbuf], Memi[bufptrs], ncols, skip, nxk, + nyk, relerr * threshold, Memi[cols]) + if (nobjs <= 0) + next + + # Compute the sharpness parameter. + call ap_sharp_round (Memr[imbuf], Memr[cnvbuf], Memi[bufptrs], + ncols, skip, nxk, nyk, Memi[cols], Memi[satur], Memr[round1], + Memr[sharp], nobjs, nonzero, skymode, datamin, datamax) + + # Compute the roundness parameters. + call ap_xy_round (Memr[imbuf], Memi[bufptrs], ncols, ker2d, nxk, + nyk, Memi[cols], inline, Memr[round2], Memr[x], + Memr[y], nobjs, skymode, datamin, datamax, xsigsq, ysigsq) + + # Test the image characeteristics of detected objects. + nstars = ap_test (Memi[cols], Memr[x], Memr[y], Memi[satur], + Memr[round1], Memr[round2], Memr[sharp], nobjs, IM_LEN(im,1), + IM_LEN(im,2), sharplo, sharphi, roundlo, roundhi) + + # Mark the stars on the display. + if ((nstars > 0) && (interactive == YES) && (id != NULL) && + (mkdetections == YES)) { + call greactivate (id, 0) + do j = 1, nstars { + #call ap_ltov (im, Memr[x+j-1], Memr[y+j-1], xc, yc, 1) + #call gmark (id, xc, yc, GM_PLUS, 1.0, 1.0) + call gmark (id, Memr[x+j-1], Memr[y+j-1], GM_PLUS, 1.0, 1.0) + } + call gdeactivate (id, 0) + } + + switch (apstati (ap, WCSOUT)) { + case WCS_PHYSICAL: + call ap_ltoo (ap, Memr[x], Memr[y], Memr[x], Memr[y], nstars) + case WCS_TV: + call ap_ltov (im, Memr[x], Memr[y], Memr[x], Memr[y], nstars) + default: + ; + } + + # Print results on the standard output. + if (interactive == YES) + call apstdout (Memr[cnvbuf], Memi[bufptrs], ncols, nyk, + Memi[cols], Memr[x], Memr[y], Memr[sharp], Memr[round1], + Memr[round2], nstars, ntotal, relerr * threshold) + + # Save the results in the file. + call apdtfout (out, Memr[cnvbuf], Memi[bufptrs], ncols, nyk, + Memi[cols], Memr[x], Memr[y], Memr[sharp], Memr[round1], + Memr[round2], nstars, ntotal, relerr * threshold, stid) + + + ntotal = ntotal + nstars + + } + + # Free space + call sfree (sp) + + return (ntotal) +end + + +# AP_DETECT -- Detec stellar objects in an image line. In order to be +# detected as a star the candidate object must be above threshold and have +# a maximum pixel value greater than any pixels within nsigma * sigma. + +int procedure ap_detect (density, ptrs, ncols, skip, nxk, nyk, threshold, cols) + +real density[ncols, ARB] # density array +int ptrs[ARB] # pointer array +int ncols # x dimesnsion of intensity buffer +int skip[nxk,ARB] # skip array +int nxk, nyk # size of convolution kernel +real threshold # density threshold +int cols[ARB] # column numbers of detected stars + +int i, j, k, kk, middle, nhalf, nobjs +define nextpix_ 11 + +begin + middle = 1 + nyk / 2 + nhalf = nxk / 2 + + # Loop over all the columns in an image line. + nobjs = 0 + for (i = 1 + nhalf; i <= ncols - nhalf; ) { + + # Test whether the density enhancement is above threshold. + if (density[i,ptrs[middle]] < threshold) + goto nextpix_ + + # Test whether a given density enhancement is a local maximum. + do j = 1, nyk { + kk = 1 + do k = i - nhalf, i + nhalf { + if (skip[kk,j] == NO) { + if (density[i,ptrs[middle]] < density[k,ptrs[j]]) + goto nextpix_ + } + kk = kk + 1 + } + } + + # Add the detected object to the list. + nobjs = nobjs + 1 + cols[nobjs] = i + + # If a local maximum is detected there can be no need to + # check pixels in this row between i and i + nhalf. + i = i + nhalf +nextpix_ + # Work on the next pixel. + i = i + 1 + } + + return (nobjs) +end + + +# AP_SHARP_ROUND -- Compute an estimate of the roundness and sharpness of the +# detected objects. The roundness parameter is computed by comparing a measure +# of the bilateral symmetry with a measure of the four-fold symmetry. The +# sharpness parameter is defined as the ratio of the difference between the +# height of the central pixel and the mean of the surrounding pixels to the +# density enhancement of the central pixel. + +procedure ap_sharp_round (data, density, ptrs, ncols, skip, nxk, nyk, cols, + satur, round, sharps, nobjs, nonzero, skymode, datamin, datamax) + +real data[ncols,ARB] # image data +real density[ncols,ARB] # density enhancements +int ptrs[ARB] # buffer pointers +int ncols # length of data array +int skip[nxk,ARB] # 2D kernel +int nxk, nyk # size of convolution kernel +int cols[ARB] # array of columns +int satur[ARB] # array of saturated state parameters +real round[ARB] # array of roundness parameters +real sharps[ARB] # array of sharpness parameters +int nobjs # number of objects +int nonzero # number of nonzero kernel elements +real skymode # estimate of the sky mode +real datamin, datamax # minimum and maximum good data values + +int i, j, k, xmiddle, ymiddle, npixels, nhalf +real pixval, midpix, temp, sharp, sum2, sum4 + +begin + # Loop over the detected objects. + nhalf = min (nxk / 2, nyk / 2) + xmiddle = 1 + nxk / 2 + ymiddle = 1 + nyk / 2 + do i = 1, nobjs { + + # Compute the first estimate of roundness. + sum2 = 0.0 + sum4 = 0.0 + do k = 0, nhalf { + do j = 1, nhalf { + sum2 = sum2 + + density[cols[i]-k,ptrs[ymiddle-j]] + + density[cols[i]+k,ptrs[ymiddle+j]] - + density[cols[i]-j,ptrs[ymiddle+k]] - + density[cols[i]+j,ptrs[ymiddle-k]] + sum4 = sum4 + + abs (density[cols[i]-k,ptrs[ymiddle-j]]) + + abs (density[cols[i]+k,ptrs[ymiddle+j]]) + + abs (density[cols[i]-j,ptrs[ymiddle+k]]) + + abs (density[cols[i]+j,ptrs[ymiddle-k]]) + } + } + if (sum2 == 0.0) + round[i] = 0.0 + else if (sum4 <= 0.0) + round[i] = INDEFR + else + round[i] = 2.0 * sum2 / sum4 + + satur[i] = NO + + # Eliminate the sharpness test if the central pixel is bad. + midpix = data[cols[i],ptrs[ymiddle]] + if (midpix > datamax) { + satur[i] = YES + sharps[i] = INDEFR + next + } + if (midpix < datamin) { + sharps[i] = INDEFR + next + } + + # Accumulate the sharpness statistic. + sharp = 0.0 + npixels = nonzero + do j = 1, nyk { + temp = 0.0 + do k = 1, nxk { + if (skip[k,j] == YES) + next + pixval = data[cols[i]-xmiddle+k,ptrs[j]] + if (pixval > datamax) { + satur[i] = YES + npixels = npixels - 1 + } else if (pixval < datamin) { + npixels = npixels - 1 + } else { + temp = temp + (pixval - skymode) + } + } + sharp = sharp + temp + } + + # Compute the sharpness statistic. + if (density[cols[i],ptrs[ymiddle]] <= 0.0 || npixels <= 0) + sharps[i] = INDEFR + else + sharps[i] = (midpix - skymode - sharp / real (npixels)) / + density[cols[i],ptrs[ymiddle]] + + } +end + + +# AP_XY_ROUND -- Estimate the x-y centers and the roundness of the detected +# objects. The height of the equivalent Gaussian function in x and y is fit by +# least squares to the marginal distribution of the image data. If either +# of these of these heights is negative set the roundess characteristic to +# -MAX_REAL, otherwise compute a roundness characteristic. At the same +# time setup the necessary sums for computing the first order corection +# to the centroid of the gaussian profile. + +procedure ap_xy_round (data, ptrs, ncols, ker2d, nxk, nyk, cols, inline, + rounds, x, y, nobjs, skymode, datamin, datamax, xsigsq, ysigsq) + +real data[ncols,ARB] # density enhancements +int ptrs[ARB] # buffer pointers +int ncols # number of columns in cylindrical buffer +real ker2d[nxk,ARB] # the gaussian convolution kernel +int nxk # size of kernel in x +int nyk # size of kernel in y +int cols[ARB] # the input positions +int inline # the input image line +real rounds[ARB] # array of sharpness parameters +real x[ARB] # output x coords +real y[ARB] # output y coords +int nobjs # number of objects +real skymode # estimate of the sky mode +real datamin, datamax # minium and maximum data values +real xsigsq, ysigsq # x-y gaussian sigma squared + +int i, j, k, xmiddle, ymiddle, n +real sumgd, sumgsq, sumg, sumd, sumdx, dgdx, sdgdx, sdgdxsq, sddgdx, sgdgdx +real pixval, p, sg, sd, wt, hx, hy, dx, dy, skylvl, xhalf, yhalf + +begin + xhalf = real (nxk / 2) + 0.5 + yhalf = real (nyk / 2) + 0.5 + xmiddle = 1 + nxk / 2 + ymiddle = 1 + nyk / 2 + + # Loop over the detected objects. + do i = 1, nobjs { + + # Initialize the x fit. + sumgd = 0.0 + sumgsq = 0.0 + sumg = 0.0 + sumd = 0.0 + sumdx = 0.0 + sdgdx = 0.0 + sdgdxsq = 0.0 + sddgdx = 0.0 + sgdgdx = 0.0 + p = 0.0 + n = 0 + + # Compute the sums required for the x fit. + do k = 1, nxk { + + sg = 0.0 + sd = 0.0 + do j = 1, nyk { + wt = real (ymiddle - abs (j - ymiddle)) + pixval = data[cols[i]-xmiddle+k,ptrs[j]] + if (pixval < datamin || pixval > datamax) + next + sd = sd + (pixval - skymode) * wt + sg = sg + ker2d[k,j] * wt + } + + if (sg <= 0.0) + next + wt = real (xmiddle - abs (k - xmiddle)) + sumgd = sumgd + wt * sg * sd + sumgsq = sumgsq + wt * sg ** 2 + sumg = sumg + wt * sg + sumd = sumd + wt * sd + sumdx = sumdx + wt * sd * (xmiddle - k) + p = p + wt + n = n + 1 + dgdx = sg * (xmiddle - k) + sdgdxsq = sdgdxsq + wt * dgdx ** 2 + sdgdx = sdgdx + wt * dgdx + sddgdx = sddgdx + wt * sd * dgdx + sgdgdx = sgdgdx + wt * sg * dgdx + } + + # Need at least three points to estimate the x height, position + # and local sky brightness of the star. + + if (n <= 2 || p <= 0.0) { + x[i] = INDEFR + y[i] = INDEFR + rounds[i] = INDEFR + next + } + + # Solve for the height of the best-fitting gaussian to the + # xmarginal. Reject the star if the height is non-positive. + + hx = sumgsq - (sumg ** 2) / p + if (hx <= 0.0) { + x[i] = INDEFR + y[i] = INDEFR + rounds[i] = INDEFR + next + } + hx = (sumgd - sumg * sumd / p) / hx + if (hx <= 0.0) { + x[i] = INDEFR + y[i] = INDEFR + rounds[i] = INDEFR + next + } + + # Solve for the new x centroid. + skylvl = (sumd - hx * sumg) / p + dx = (sgdgdx - (sddgdx - sdgdx * (hx * sumg + skylvl * p))) / + (hx * sdgdxsq / xsigsq) + if (abs (dx) > xhalf) { + if (sumd == 0.0) + dx = 0.0 + else + dx = sumdx / sumd + if (abs (dx) > xhalf) + dx = 0.0 + } + x[i] = (cols[i] - xmiddle + 1) + dx + + # Initialize y fit. + sumgd = 0.0 + sumgsq = 0.0 + sumg = 0.0 + sumd = 0.0 + sumdx = 0.0 + sdgdx = 0.0 + sdgdxsq = 0.0 + sddgdx = 0.0 + sgdgdx = 0.0 + p = 0.0 + n = 0 + + do j = 1, nyk { + sg = 0.0 + sd = 0.0 + do k = 1, nxk { + wt = real (xmiddle - abs (k - xmiddle)) + pixval = data[cols[i]-xmiddle+k,ptrs[j]] + if (pixval < datamin || pixval > datamax) + next + sd = sd + (pixval - skymode) * wt + sg = sg + ker2d[k,j] * wt + } + if (sg <= 0.0) + next + wt = real (ymiddle - abs (j - ymiddle)) + sumgd = sumgd + wt * sg * sd + sumgsq = sumgsq + wt * sg ** 2 + sumg = sumg + wt * sg + sumd = sumd + wt * sd + sumdx = sumdx + wt * sd * (j - ymiddle) + p = p + wt + n = n + 1 + dgdx = sg * (ymiddle - j) + sdgdx = sdgdx + wt * dgdx + sdgdxsq = sdgdxsq + wt * dgdx ** 2 + sddgdx = sddgdx + wt * sd * dgdx + sgdgdx = sgdgdx + wt * sg * dgdx + } + + # Need at least three points to estimate the y height, position + # and local sky brightness of the star. + + if (n <= 2 || p <= 0.0) { + x[i] = INDEFR + y[i] = INDEFR + rounds[i] = INDEFR + next + } + + # Solve for the height of the best-fitting gaussian to the + # y marginal. Reject the star if the height is non-positive. + + hy = sumgsq - (sumg ** 2) / p + if (hy <= 0.0) { + x[i] = INDEFR + y[i] = INDEFR + rounds[i] = INDEFR + next + } + hy = (sumgd - sumg * sumd / p) / (sumgsq - (sumg ** 2) / p) + if (hy <= 0.0) { + x[i] = INDEFR + y[i] = INDEFR + rounds[i] = INDEFR + next + } + + # Solve for the new x centroid. + skylvl = (sumd - hy * sumg) / p + dy = (sgdgdx - (sddgdx - sdgdx * (hy * sumg + skylvl * p))) / + (hy * sdgdxsq / ysigsq) + if (abs (dy) > yhalf) { + if (sumd == 0.0) + dy = 0.0 + else + dy = sumdx / sumd + if (abs (dy) > yhalf) + dy = 0.0 + } + y[i] = (inline - ymiddle + 1) + dy + + # Compute the roundness. + rounds[i] = 2.0 * (hx - hy) / (hx + hy) + } +end + + +# AP_TEST -- Test the characteristic of the detected images for roundness +# and sharpness. + +int procedure ap_test (cols, x, y, satur, round1, round2, sharps, nobjs, + ncols, nlines, sharplo, sharphi, roundlo, roundhi) + +int cols[ARB] # col IDS of detected images +real x[ARB] # x positions +real y[ARB] # y positions +int satur[ARB] # saturation condition +real round1[ARB] # first roundness parameters +real round2[ARB] # second roundness parameters +real sharps[ARB] # sharpness parameters +int nobjs # number of objects +int ncols, nlines # size of the input image +real sharplo, sharphi # sharpness parameters +real roundlo, roundhi # roundness parameters + +int i, nstars + +begin + # Loop over the detected objects. + nstars = 0 + do i = 1, nobjs { + + # Compute the sharpness statistic + if (! IS_INDEFR(sharps[i]) && (sharps[i] < sharplo || + sharps[i] > sharphi)) + next + if (IS_INDEFR(round1[i]) || round1[i] < roundlo || + round1[i] > roundhi) + next + if (satur[i] == NO) { + if (IS_INDEFR(round2[i]) || round2[i] < roundlo || + round2[i] > roundhi) + next + } + if (IS_INDEFR(x[i]) || x[i] < 0.5 || x[i] > (ncols+0.5)) + next + if (IS_INDEFR(y[i]) || y[i] < 0.5 || y[i] > (nlines+0.5)) + next + + # Add object to the list. + nstars = nstars + 1 + cols[nstars] = cols[i] + x[nstars] = x[i] + y[nstars] = y[i] + sharps[nstars] = sharps[i] + round1[nstars] = round1[i] + round2[nstars] = round2[i] + } + + return (nstars) +end diff --git a/noao/digiphot/apphot/daofind/apimset.x b/noao/digiphot/apphot/daofind/apimset.x new file mode 100644 index 00000000..6f1e8cdc --- /dev/null +++ b/noao/digiphot/apphot/daofind/apimset.x @@ -0,0 +1,17 @@ +include <imset.h> + +# AP_IMSET -- Setup image boundary extension charactersitics. + +procedure ap_imset (im, boundary, npix, constant) + +pointer im # pointer to the image +int boundary # type of boundary condition +int npix # number of pixels of boundary entension +real constant # constant for constant boundary extension + +begin + call imseti (im, IM_TYBNDRY, boundary) + call imseti (im, IM_NBNDRYPIX, npix) + if (boundary == BT_CONSTANT) + call imsetr (im, IM_BNDRYPIXVAL, constant) +end diff --git a/noao/digiphot/apphot/daofind/daofind.key b/noao/digiphot/apphot/daofind/daofind.key new file mode 100644 index 00000000..e68ceb80 --- /dev/null +++ b/noao/digiphot/apphot/daofind/daofind.key @@ -0,0 +1,66 @@ + Interactive Keystroke Commands + +? Print help +: Colon commands +v Verify the critical parameters +w Save the current parameters +d Plot radial profile of star near cursor +i Interactively set parameters using star near cursor +f Find stars in the image +spbar Find stars in the image, output results +q Exit task + + + Colon Commands + +:show [data/center] List the parameters + + Colon Parameter Editing Commands + +# Image and file name parameters + +:image [string] Image name +:output [string] Output file name + +# Data dependent parameters + +:scale [value] Image scale (units per pixel) +:fwhmpsf [value] Full width half maximum of psf (scale units) +:emission [y/n] Emission feature (y), absorption (n) +:sigma [value] Standard deviation of sky (counts) +:datamin [value] Minimum good data value (counts) +:datamax [value] Maximum good data value (counts) + +# Noise description parameters + +:noise [string] Noise model (constant|poisson) +:gain [string] Gain image header keyword +:ccdread [string] Readout noise image header keyword +:epadu [value] Gain (electrons per adu) +:readnoise [value] Readout noise (electrons) + +# Observation parameters + +:exposure [string] Exposure time image header keyword +:airmass [string] Airmass image header keyword +:filter [string] Filter image header keyword +:obstime [string] Time of observation image header keyword +:itime [value] Exposure time (time units) +:xairmass [value] Airmass value (number) +:ifilter [string] Filter id string +:otime [string] Time of observation (time units) + +# Object detection parameters + +:nsigma [value] Size of Gaussian kernel (sigma) +:threshold [value] Detection intensity threshold (counts) +:ratio [value] Sigmay / sigmax of Gaussian kernel +:theta [value] Position angle of Gaussian kernel +:sharplo [value] Lower bound on sharpness +:sharphi [value] Upper bound on sharpness +:roundlo [value] Lower bound on roundness +:roundhi [value] Upper bound on roundness + +# Plotting and marking commands + +:mkdetections [y/n] Mark detections on the image display diff --git a/noao/digiphot/apphot/daofind/idaofind.key b/noao/digiphot/apphot/daofind/idaofind.key new file mode 100644 index 00000000..7da8fe42 --- /dev/null +++ b/noao/digiphot/apphot/daofind/idaofind.key @@ -0,0 +1,8 @@ + Interactive Daofind Setup Menu + + v Mark and verify critical daofind parameters (f,s) + + f Mark and verify the full-width half-maximum of the psf + s Mark and verify the standard deviation of the background + l Mark and verify the minimum good data value + u Mark and verify the maximum good data value diff --git a/noao/digiphot/apphot/daofind/mkpkg b/noao/digiphot/apphot/daofind/mkpkg new file mode 100644 index 00000000..28935c11 --- /dev/null +++ b/noao/digiphot/apphot/daofind/mkpkg @@ -0,0 +1,37 @@ +# DAOFIND task + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + apbfdfind.x <imhdr.h> <mach.h> \ + <imio.h> "../lib/apphot.h" \ + "../lib/noise.h" "../lib/find.h" + apconvolve.x <imhdr.h> <imset.h> \ + "../lib/find.h" + apegkernel.x <math.h> "../lib/find.h" + apfdcolon.x "../lib/apphot.h" "../lib/noise.h" \ + "../lib/find.h" "../lib/display.h" + apfdconfirm.x + apfdfind.x "../lib/apphot.h" "../lib/display.h" \ + "../lib/find.h" <imhdr.h> + apfdfree.x "../lib/apphotdef.h" + apfdgpars.x "../lib/noise.h" "../lib/display.h" + apfdinit.x "../lib/apphotdef.h" "../lib/finddef.h" + apfdpars.x "../lib/display.h" + apfdradsetup.x + apfdshow.x "../lib/apphot.h" "../lib/noise.h" \ + "../lib/find.h" + apfdstars.x <imhdr.h> <mach.h> \ + "../lib/apphot.h" "../lib/noise.h" \ + "../lib/find.h" "../lib/display.h" + apfind.x <imhdr.h> <mach.h> \ + <gset.h> "../lib/apphot.h" + apimset.x <imset.h> + t_daofind.x <fset.h> <imhdr.h> \ + <gset.h> "../lib/apphot.h" \ + "../lib/noise.h" "../lib/find.h" \ + <imhdr.h> + ; diff --git a/noao/digiphot/apphot/daofind/t_daofind.x b/noao/digiphot/apphot/daofind/t_daofind.x new file mode 100644 index 00000000..eaf77fbf --- /dev/null +++ b/noao/digiphot/apphot/daofind/t_daofind.x @@ -0,0 +1,419 @@ +include <imhdr.h> +include <gset.h> +include <fset.h> +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/find.h" + +# T_DAOFIND -- Automatically detect objects in an image given the full +# width half maximum of the image point spread function and a detection +# threshold. + +procedure t_daofind () + +pointer image # pointer to input image +pointer denimage # pointer to density enhancement image +pointer skyimage # pointer to the sky image +int output # the results file descriptor +int boundary # type of boundary extension +real constant # constant for constant boundary extension +int interactive # interactive mode +int verify # verify mode +int update # update critical parameters +int verbose # verbose mode +int cache # cache the image pixels + +pointer im, cnv, sky, sp, outfname, denname, skyname, str +pointer ap, cname, display, graphics, id, gd +int limlist, lolist, densave, skysave, out, root, stat, imlist, olist +int wcs, req_size, old_size, buf_size, memstat + +real clgetr() +pointer gopen(), immap(), ap_immap() +int imtlen(), clplen(), btoi(), clgwrd(), aptmpimage() +int open(), strncmp(), strlen(), fnldir(), ap_fdfind() +int imtopenp(), clpopnu(), clgfil(), imtgetim(), ap_memstat(), sizeof() +bool clgetb(), streq() + +begin + # Flush STDOUT on a new line. + call fseti (STDOUT, F_FLUSHNL, YES) + + # Allocate working space. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (denimage, SZ_FNAME, TY_CHAR) + call salloc (skyimage, SZ_FNAME, TY_CHAR) + call salloc (display, SZ_FNAME, TY_CHAR) + call salloc (graphics, SZ_FNAME, TY_CHAR) + + call salloc (outfname, SZ_FNAME, TY_CHAR) + call salloc (denname, SZ_FNAME, TY_CHAR) + call salloc (skyname, SZ_FNAME, TY_CHAR) + call salloc (cname, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Fetch the image and output file lists. + imlist = imtopenp ("image") + limlist = imtlen (imlist) + olist = clpopnu ("output") + lolist = clplen (olist) + + # Get prefix for density enhancement and sky images. + call clgstr ("starmap", Memc[denimage], SZ_FNAME) + call clgstr ("skymap", Memc[skyimage], SZ_FNAME) + + # Get the image boundary extensions parameters. + boundary = clgwrd ("boundary", Memc[str], SZ_LINE, + ",constant,nearest,reflect,wrap,") + constant = clgetr ("constant") + + call clgstr ("icommands.p_filename", Memc[cname], SZ_FNAME) + if (Memc[cname] != EOS) + interactive = NO + else + interactive = btoi (clgetb ("interactive")) + verbose = btoi (clgetb ("verbose")) + verify = btoi (clgetb ("verify")) + update = btoi (clgetb ("update")) + cache = btoi (clgetb ("cache")) + + # Get the parameters. + call ap_fdgpars (ap) + + # Confirm the algorithm parameters. + if (verify == YES && interactive == NO) { + call ap_fdconfirm (ap) + if (update == YES) + call ap_fdpars (ap) + } + + # Get the wcs info. + wcs = clgwrd ("wcsout", Memc[str], SZ_LINE, WCSOUTSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the output coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call apseti (ap, WCSOUT, wcs) + + # Get the graphics and display devices. + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + call clgstr ("display", Memc[display], SZ_FNAME) + + # Open the graphics and display devices. + if (interactive == YES) { + if (Memc[graphics] == EOS) + gd = NULL + else { + iferr { + gd = gopen (Memc[graphics], APPEND+AW_DEFER, STDGRAPH) + } then { + call eprintf ( + "Warning: Error opening the graphics device.\n") + gd = NULL + } + } + if (Memc[display] == EOS) + id = NULL + else if (streq (Memc[graphics], Memc[display])) + id = gd + else { + iferr { + id = gopen (Memc[display], APPEND, STDIMAGE) + } then { + call eprintf ( + "Warning: Graphics overlay not available for display device.\n") + id = NULL + } + } + } else { + gd = NULL + id = NULL + } + + # Loop over the images. + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open the image and get the required keywords from the header. + + im = immap (Memc[image], READ_ONLY, 0) + call apimkeys (ap, im, Memc[image]) + + # Set the image display viewport. + if ((id != NULL) && (id != gd)) + call ap_gswv (id, Memc[image], im, 4) + + # Cache the input image pixels. + req_size = MEMFUDGE * (IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + 2 * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (TY_REAL)) + memstat = ap_memstat (cache, req_size, old_size) + if (memstat == YES) + call ap_pcache (im, INDEFI, buf_size) + + # Determine the results file name. If output is a null string or + # a directory name then the extension "coo" is added to the + # root image name and the appropriate version number is appended + # in order to construct a default output file name. + + out = NULL + if (lolist == 0) { + call strcpy ("", Memc[outfname], SZ_FNAME) + } else { + stat = clgfil (olist, Memc[output], SZ_FNAME) + root = fnldir (Memc[output], Memc[outfname], SZ_FNAME) + if (strncmp ("default", Memc[output+root], 7) == 0 || root == + strlen (Memc[output])) { + call apoutname (Memc[image], Memc[outfname], "coo", + Memc[outfname], SZ_FNAME) + lolist = limlist + } else if (stat != EOF) { + call strcpy (Memc[output], Memc[outfname], SZ_FNAME) + } else { + call apoutname (Memc[image], Memc[outfname], "coo", + Memc[outfname], SZ_FNAME) + lolist = limlist + } + } + call apsets (ap, OUTNAME, Memc[outfname]) + + # Set up the directory and name for the density enhancement image. + # Note that the default value for denimage is the null string + # indicating the current directory. If denimage is null or + # contains only a directory specification make a temporary image + # name using the prefix "den", otherwise use the user specified + # prefix. + + densave = aptmpimage (Memc[image], Memc[denimage], "den", + Memc[denname], SZ_FNAME) + + # Set up the directory and name for the sky values image. + + skysave = aptmpimage (Memc[image], Memc[skyimage], "sky", + Memc[skyname], SZ_FNAME) + + # Find the stars in an image. + if (interactive == NO) { + + cnv = NULL + if (Memc[cname] != EOS) { + stat = ap_fdfind (Memc[denname], Memc[skyname], ap, im, + NULL, NULL, out, boundary, constant, densave, + skysave, NO, cache) + } else { + cnv = ap_immap (Memc[denname], im, ap, densave) + if (memstat == YES) + call ap_pcache (cnv, INDEFI, buf_size) + if (skysave == YES) { + sky = ap_immap (Memc[skyname], im, ap, skysave) + if (memstat == YES) + call ap_pcache (sky, INDEFI, buf_size) + } else + sky = NULL + if (Memc[outfname] != EOS) + out = open (Memc[outfname], NEW_FILE, TEXT_FILE) + else + out = NULL + call ap_bfdfind (im, cnv, sky, out, ap, boundary, + constant, verbose) + call imunmap (cnv) + if (sky != NULL) + call imunmap (sky) + if (densave == NO) + call imdelete (Memc[denname]) + stat = NO + } + + } else + stat = ap_fdfind (Memc[denname], Memc[skyname], ap, im, gd, + id, out, boundary, constant, densave, skysave, YES, + cache) + + # Clean up files. + call imunmap (im) + call close (out) + + # Uncache memory + call fixmem (old_size) + + if (stat == YES) + break + } + + # Close up the graphics system. + if (id == gd && id != NULL) + call gclose (id) + else { + if (gd != NULL) + call gclose (gd) + if (id != NULL) + call gclose (id) + } + + # Close image and output file lists. + call ap_fdfree (ap) + call imtclose (imlist) + call clpcls (olist) + call sfree (sp) +end + + +# AP_IMMAP -- Map the output image. + +pointer procedure ap_immap (outname, im, ap, save) + +char outname[ARB] # convolved image name +pointer im # pointer to input image +pointer ap # pointer to the apphot structure +int save # save the convolved image + +int newimage +pointer outim +real tmp_scale, tmp_fwhmpsf, tmp_ratio, tmp_theta, tmp_nsigma +real tmp_datamin, tmp_datamax +bool fp_equalr() +int imaccess() +pointer immap() +real apstatr(), imgetr() +errchk imgetr() + +begin + # Check to see if the image already exists. If it does not + # open a new image and write the PSF characteristics into the + # image user area, otherwise fetch the PSF parameters from the image + # header. + + if (save == NO || imaccess (outname, READ_ONLY) == NO) { + + outim = immap (outname, NEW_COPY, im) + IM_PIXTYPE(outim) = TY_REAL + call imaddr (outim, "SCALE", 1.0 / apstatr (ap, SCALE)) + call imaddr (outim, "FWHMPSF", apstatr (ap, FWHMPSF)) + call imaddr (outim, "NSIGMA", apstatr (ap, NSIGMA)) + call imaddr (outim, "RATIO", apstatr (ap, RATIO)) + call imaddr (outim, "THETA", apstatr (ap, THETA)) + call imaddr (outim, "DATAMIN", apstatr (ap, DATAMIN)) + call imaddr (outim, "DATAMAX", apstatr (ap, DATAMIN)) + + } else { + + outim = immap (outname, READ_ONLY, 0) + iferr (tmp_scale = 1.0 / imgetr (outim, "SCALE")) + tmp_scale = INDEFR + iferr (tmp_fwhmpsf = imgetr (outim, "FWHMPSF")) + tmp_fwhmpsf = INDEFR + iferr (tmp_nsigma = imgetr (outim, "NSIGMA")) + tmp_nsigma = INDEFR + iferr (tmp_ratio = imgetr (outim, "RATIO")) + tmp_ratio = INDEFR + iferr (tmp_theta = imgetr (outim, "THETA")) + tmp_theta = INDEFR + iferr (tmp_datamin = imgetr (outim, "DATAMIN")) + tmp_datamin = INDEFR + iferr (tmp_datamax = imgetr (outim, "DATAMAX")) + tmp_datamax = INDEFR + + if (IS_INDEFR(tmp_scale) || ! fp_equalr (tmp_scale, apstatr (ap, + SCALE))) + newimage = YES + else if (IS_INDEFR(tmp_fwhmpsf) || ! fp_equalr (tmp_fwhmpsf, + apstatr (ap, FWHMPSF))) + newimage = YES + else if (IS_INDEFR(tmp_nsigma) || ! fp_equalr (tmp_nsigma, + apstatr (ap, NSIGMA))) + newimage = YES + else if (IS_INDEFR(tmp_ratio) || ! fp_equalr (tmp_ratio, + apstatr (ap, RATIO))) + newimage = YES + else if (IS_INDEFR(tmp_theta) || ! fp_equalr (tmp_theta, + apstatr (ap, THETA))) + newimage = YES + else if (IS_INDEFR(tmp_datamin) || ! fp_equalr (tmp_datamin, + apstatr (ap, DATAMIN))) + newimage = YES + else if (IS_INDEFR(tmp_datamax) || ! fp_equalr (tmp_datamax, + apstatr (ap, DATAMAX))) + newimage = YES + else + newimage = NO + + if (newimage == YES) { + call imunmap (outim) + call imdelete (outname) + outim = immap (outname, NEW_COPY, im) + IM_PIXTYPE(outim) = TY_REAL + call imaddr (outim, "SCALE", 1.0 / apstatr (ap, SCALE)) + call imaddr (outim, "FWHMPSF", apstatr (ap, FWHMPSF)) + call imaddr (outim, "NSIGMA", apstatr (ap, NSIGMA)) + call imaddr (outim, "RATIO", apstatr (ap, RATIO)) + call imaddr (outim, "THETA", apstatr (ap, THETA)) + } + } + + return (outim) +end + + +# AP_OUTMAP -- Manufacture an output coordinate file name. + +procedure ap_outmap (ap, out, root) + +pointer ap # pointer to the apphot structure +int out # the output file descriptor +char root[ARB] # root of the previous output file name, may be + # updated + +int findex, lindex, version +pointer sp, image, outname, newoutname +int strmatch(), gstrmatch(), open(), fnldir(), access() + +begin + if (out != NULL) + call close (out) + + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (outname, SZ_FNAME, TY_CHAR) + call salloc (newoutname, SZ_FNAME, TY_CHAR) + + # Get the old names. + call apstats (ap, IMNAME, Memc[image], SZ_FNAME) + call apstats (ap, OUTNAME, Memc[outname], SZ_FNAME) + + # Manufacture a new name. Search first for existing files with + # the form "outname.coo.*" and create a new version number. + # If the first search fails look for names containing root + # and append a version number to create a new output file + # name. Otherwise simply use the output name. + + if (Memc[outname] == EOS) { + Memc[newoutname] = EOS + } else if (strmatch (Memc[outname], "\.coo\.") > 0) { + findex = fnldir (Memc[outname], Memc[newoutname], SZ_FNAME) + call apoutname (Memc[image], Memc[newoutname], "coo", + Memc[newoutname], SZ_FNAME) + } else if (gstrmatch (Memc[outname], root, findex, lindex) > 0) { + repeat { + version = version + 1 + call strcpy (Memc[outname], Memc[newoutname], SZ_FNAME) + call sprintf (Memc[newoutname+lindex], SZ_FNAME, ".%d") + call pargi (version) + } until (access (Memc[newoutname], 0, 0) == NO) + } else { + version = 1 + call strcpy (Memc[outname], root, SZ_FNAME) + call strcpy (Memc[outname], Memc[newoutname], SZ_FNAME) + } + + # Open the output image. + if (Memc[newoutname] == EOS) + out = NULL + else + out = open (Memc[newoutname], NEW_FILE, TEXT_FILE) + call apsets (ap, OUTNAME, Memc[newoutname]) + + call sfree (sp) +end |