include include include include include include include "starfind.h" # SF_FIND -- Find stars in an image using a pattern matching technique and # a circularly symmetric Gaussian pattern. procedure sf_find (im, out, sf, nxblock, nyblock, wcs, wxformat, wyformat, boundary, constant, verbose) pointer im #I pointer to the input image int out #I the output file descriptor pointer sf #I pointer to the apphot structure int nxblock #I the x dimension blocking factor int nyblock #I the y dimension blocking factor char wcs[ARB] #I the world coordinate system char wxformat[ARB] #I the x axis world coordinate format char wyformat[ARB] #I the y axis world coordinate format int boundary #I type of boundary extension real constant #I constant for constant boundary extension int verbose #I verbose switch int i, j, fwidth, swidth, norm int l1, l2, c1, c2, ncols, nlines, nxb, nyb, nstars, stid pointer sp, gker2d, ngker2d, skip, fmtstr, twxformat, twyformat pointer imbuf, denbuf, str, mw, ct real sigma, nsigma, a, b, c, f, gsums[LEN_GAUSS], relerr, dmin, dmax real maglo, maghi bool streq() int sf_stfind() pointer mw_openim(), mw_sctran() real sf_egkernel() errchk mw_openim(), mw_sctran(), mw_gattrs() begin # Allocate working space. call smark (sp) call salloc (twxformat, SZ_FNAME, TY_CHAR) call salloc (twyformat, SZ_FNAME, TY_CHAR) call salloc (fmtstr, SZ_LINE, TY_CHAR) call salloc (str, SZ_LINE, TY_CHAR) # Compute the parameters of the Gaussian kernel. sigma = HWHM_TO_SIGMA * SF_HWHMPSF(sf) nsigma = SF_FRADIUS(sf) / HWHM_TO_SIGMA call sf_egparams (sigma, 1.0, 0.0, nsigma, a, b, c, f, fwidth, fwidth) # Compute the separation parameter swidth = max (2, int (SF_SEPMIN(sf) * SF_HWHMPSF(sf) + 0.5)) # Compute the minimum and maximum pixel values. if (IS_INDEFR(SF_DATAMIN(sf)) && IS_INDEFR(SF_DATAMAX(sf))) { norm = YES dmin = -MAX_REAL dmax = MAX_REAL } else { norm = NO if (IS_INDEFR(SF_DATAMIN(sf))) dmin = -MAX_REAL else dmin = SF_DATAMIN(sf) if (IS_INDEFR(SF_DATAMAX(sf))) dmax = MAX_REAL else dmax = SF_DATAMAX(sf) } # Compute the magnitude limits if (IS_INDEFR(SF_MAGLO(sf))) maglo = -MAX_REAL else maglo = SF_MAGLO(sf) if (IS_INDEFR(SF_MAGHI(sf))) maghi = MAX_REAL else maghi = SF_MAGHI(sf) # Open the image WCS. if (wcs[1] == EOS) { mw = NULL ct = NULL } else { iferr { mw = mw_openim (im) } then { call erract (EA_WARN) mw = NULL ct = NULL } else { iferr { ct = mw_sctran (mw, "logical", wcs, 03B) } then { call erract (EA_WARN) ct = NULL call mw_close (mw) mw = NULL } } } # Set the WCS formats. if (ct == NULL) call strcpy (wxformat, Memc[twxformat], SZ_FNAME) else if (wxformat[1] == EOS) { if (mw != NULL) { iferr (call mw_gwattrs (mw, 1, "format", Memc[twxformat], SZ_FNAME)) { if (streq (wcs, "world")) call strcpy ("%11.8g", Memc[twxformat], SZ_FNAME) else call strcpy ("%9.3f", Memc[twxformat], SZ_FNAME) } } else call strcpy ("%9.3f", Memc[twxformat], SZ_FNAME) } else call strcpy (wxformat, Memc[twxformat], SZ_FNAME) if (ct == NULL) call strcpy (wyformat, Memc[twyformat], SZ_FNAME) else if (wyformat[1] == EOS) { if (mw != NULL) { iferr (call mw_gwattrs (mw, 2, "format", Memc[twyformat], SZ_FNAME)) { if (streq (wcs, "world")) call strcpy ("%11.8g", Memc[twyformat], SZ_FNAME) else call strcpy ("%9.3f", Memc[twyformat], SZ_FNAME) } } else call strcpy ("%9.3f", Memc[twyformat], SZ_FNAME) } else call strcpy (wyformat, Memc[twyformat], SZ_FNAME) # Create the output format string. call sprintf (Memc[fmtstr], SZ_LINE, " %s %s %s %s %s %s %s %s %s %s %s\n") call pargstr ("%9.3f") call pargstr ("%9.3f") call pargstr (Memc[twxformat]) call pargstr (Memc[twyformat]) call pargstr ("%7.2f") call pargstr ("%6d") call pargstr ("%6.2f") call pargstr ("%6.3f") call pargstr ("%6.1f") call pargstr ("%7.3f") call pargstr ("%6d") # Set up the image boundary extension characteristics. call imseti (im, IM_TYBNDRY, boundary) call imseti (im, IM_NBNDRYPIX, 1 + fwidth / 2 + swidth) if (boundary == BT_CONSTANT) call imsetr (im, IM_BNDRYPIXVAL, constant) # Set up the blocking factor. # Compute the magnitude limits if (IS_INDEFI(nxblock)) nxb = IM_LEN(im,1) else nxb = nxblock if (IS_INDEFI(nyblock)) nyb = IM_LEN(im,2) else nyb = nyblock # Print the detection criteria on the standard output. if (verbose == YES) { call fstats (out, F_FILENAME, Memc[str], SZ_LINE) call printf ("\nImage: %s Output: %s\n") call pargstr (IM_HDRFILE(im)) call pargstr (Memc[str]) call printf ("Detection Parameters\n") call printf ( " Hwhmpsf: %0.3f (pixels) Threshold: %g (ADU) Npixmin: %d\n") call pargr (SF_HWHMPSF(sf)) call pargr (SF_THRESHOLD(sf)) call pargi (SF_NPIXMIN(sf)) call printf (" Datamin: %g (ADU) Datamax: %g (ADU)\n") call pargr (SF_DATAMIN(sf)) call pargr (SF_DATAMAX(sf)) call printf (" Fradius: %0.3f (HWHM) Sepmin: %0.3f (HWHM)\n\n") call pargr (SF_FRADIUS(sf)) call pargr (SF_SEPMIN(sf)) } if (out != NULL) { call fstats (out, F_FILENAME, Memc[str], SZ_LINE) call fprintf (out, "\n# Image: %s Output: %s\n") call pargstr (IM_HDRFILE(im)) call pargstr (Memc[str]) call fprintf (out, "# Detection Parameters\n") call fprintf (out, "# Hwhmpsf: %0.3f (pixels) Threshold: %g (ADU) Npixmin: %d\n") call pargr (SF_HWHMPSF(sf)) call pargr (SF_THRESHOLD(sf)) call pargi (SF_NPIXMIN(sf)) call fprintf (out, "# Datamin: %g (ADU) Datamax: %g (ADU)\n") call pargr (SF_DATAMIN(sf)) call pargr (SF_DATAMAX(sf)) call fprintf (out, "# Fradius: %g (HWHM) Sepmin: %g (HWHM)\n") call pargr (SF_FRADIUS(sf)) call pargr (SF_SEPMIN(sf)) call fprintf (out, "# Selection Parameters\n") call pargi (SF_NPIXMIN(sf)) call fprintf (out, "# Maglo: %0.3f Maghi: %0.3f\n") call pargr (SF_MAGLO(sf)) call pargr (SF_MAGHI(sf)) call fprintf (out, "# Roundlo: %0.3f Roundhi: %0.3f\n") call pargr (SF_ROUNDLO(sf)) call pargr (SF_ROUNDHI(sf)) call fprintf (out, "# Sharplo: %0.3f Sharphi: %0.3f\n") call pargr (SF_SHARPLO(sf)) call pargr (SF_SHARPHI(sf)) call fprintf (out, "# Columns\n") call fprintf (out, "# 1: X 2: Y \n") if (ct == NULL) { call fprintf (out, "# 3: Mag 4: Area\n") call fprintf (out, "# 5: Hwhm 6: Roundness\n") call fprintf (out, "# 7: Pa 8: Sharpness\n\n") } else { call fprintf (out, "# 3: Wx 4: Wy \n") call fprintf (out, "# 5: Mag 6: Area\n") call fprintf (out, "# 7: Hwhm 8: Roundness\n") call fprintf (out, "# 9: Pa 10: Sharpness\n\n") } } # Process the image block by block. stid = 1 nstars = 0 do j = 1, IM_LEN(im,2), nyb { l1 = j l2 = min (IM_LEN(im,2), j + nyb - 1) nlines = l2 - l1 + 1 + 2 * (fwidth / 2 + swidth) do i = 1, IM_LEN(im,1), nxb { # Allocate space for the convolution kernel. call malloc (gker2d, fwidth * fwidth, TY_REAL) call malloc (ngker2d, fwidth * fwidth, TY_REAL) call malloc (skip, fwidth * fwidth, TY_INT) # Allocate space for the data and the convolution. c1 = i c2 = min (IM_LEN(im,1), i + nxb - 1) ncols = c2 - c1 + 1 + 2 * (fwidth / 2 + swidth) call malloc (imbuf, ncols * nlines, TY_REAL) call malloc (denbuf, ncols * nlines, TY_REAL) # Compute the convolution kernels. relerr = sf_egkernel (Memr[gker2d], Memr[ngker2d], Memi[skip], fwidth, fwidth, gsums, a, b, c, f) # Do the convolution. if (norm == YES) call sf_fconvolve (im, c1, c2, l1, l2, swidth, Memr[imbuf], Memr[denbuf], ncols, nlines, Memr[ngker2d], Memi[skip], fwidth, fwidth, gsums[GAUSS_SGOP]) else call sf_gconvolve (im, c1, c2, l1, l2, swidth, Memr[imbuf], Memr[denbuf], ncols, nlines, Memr[gker2d], Memi[skip], fwidth, fwidth, gsums, dmin, dmax) # Find the stars. nstars = sf_stfind (out, Memr[imbuf], Memr[denbuf], ncols, nlines, c1, c2, l1, l2, swidth, Memi[skip], fwidth, fwidth, SF_HWHMPSF(sf), SF_THRESHOLD(sf), dmin, dmax, ct, SF_NPIXMIN(sf), maglo, maghi, SF_ROUNDLO(sf), SF_ROUNDHI(sf), SF_SHARPLO(sf), SF_SHARPHI(sf), Memc[fmtstr], stid, verbose) # Increment the sequence number. stid = stid + nstars # Free the memory. call mfree (imbuf, TY_REAL) call mfree (denbuf, TY_REAL) call mfree (gker2d, TY_REAL) call mfree (ngker2d, TY_REAL) call mfree (skip, TY_INT) } } # Print out the selection parameters. if (verbose == YES) { call printf ("\nSelection Parameters\n") call printf ( " Maglo: %0.3f Maghi: %0.3f\n") call pargr (SF_MAGLO(sf)) call pargr (SF_MAGHI(sf)) call printf ( " Roundlo: %0.3f Roundhi: %0.3f\n") call pargr (SF_ROUNDLO(sf)) call pargr (SF_ROUNDHI(sf)) call printf ( " Sharplo: %0.3f Sharphi: %0.3f\n") call pargr (SF_SHARPLO(sf)) call pargr (SF_SHARPHI(sf)) } if (mw != NULL) { call mw_ctfree (ct) call mw_close (mw) } call sfree (sp) end # SF_STFIND -- Detect images in the convolved image and then compute image # characteristics using the original image. int procedure sf_stfind (out, imbuf, denbuf, ncols, nlines, c1, c2, l1, l2, sepmin, skip, nxk, nyk, hwhmpsf, threshold, datamin, datamax, ct, nmin, maglo, maghi, roundlo, roundhi, sharplo, sharphi, fmtstr, stid, verbose) int out #I the output file descriptor real imbuf[ncols,nlines] #I the input data buffer real denbuf[ncols,nlines] #I the input density enhancements buffer int ncols, nlines #I the dimensions of the input buffers int c1, c2 #I the image columns limits int l1, l2 #I the image lines limits int sepmin #I the minimum object separation int skip[nxk,ARB] #I the pixel fitting array int nxk, nyk #I the dimensions of the fitting array real hwhmpsf #I the HWHM of the PSF in pixels real threshold #I the threshold for object detection real datamin, datamax #I the minimum and maximum good data values pointer ct #I the coordinate transformation pointer int nmin #I the minimum number of good object pixels real maglo,maghi #I the magnitude estimate limits real roundlo,roundhi #I the ellipticity estimate limits real sharplo, sharphi #I the sharpness estimate limits char fmtstr[ARB] #I the format string int stid #U the object sequence number int verbose #I verbose mode int line1, line2, inline, xmiddle, ymiddle, ntotal, nobjs, nstars pointer sp, cols, sharp, x, y, ellip, theta, npix, mag, size int sf_detect(), sf_test() begin # Set up useful line and column limits. line1 = 1 + sepmin + nyk / 2 line2 = nlines - sepmin - nyk / 2 xmiddle = 1 + nxk / 2 ymiddle = 1 + nyk / 2 # Set up a cylindrical buffers and some working space for # the detected images. call smark (sp) call salloc (cols, ncols, TY_INT) call salloc (x, ncols, TY_REAL) call salloc (y, ncols, TY_REAL) call salloc (mag, ncols, TY_REAL) call salloc (npix, ncols, TY_INT) call salloc (size, ncols, TY_REAL) call salloc (ellip, ncols, TY_REAL) call salloc (theta, ncols, TY_REAL) call salloc (sharp, ncols, TY_REAL) # Generate the starlist line by line. ntotal = 0 do inline = line1, line2 { # Detect local maximum in the density enhancement buffer. nobjs = sf_detect (denbuf[1,inline-nyk/2-sepmin], ncols, sepmin, nxk, nyk, threshold, Memi[cols]) if (nobjs <= 0) next # Do not skip the middle pixel in the moments computation. call sf_moments (imbuf[1,inline-nyk/2], denbuf[1,inline-nyk/2], ncols, skip, nxk, nyk, Memi[cols], Memr[x], Memr[y], Memi[npix], Memr[mag], Memr[size], Memr[ellip], Memr[theta], Memr[sharp], nobjs, datamin, datamax, threshold, hwhmpsf, real (-sepmin - nxk / 2 + c1 - 1), real (inline - sepmin - nyk + l1 - 1)) # Test the image characeteristics of detected objects. nstars = sf_test (Memi[cols], Memr[x], Memr[y], Memi[npix], Memr[mag], Memr[size], Memr[ellip], Memr[theta], Memr[sharp], nobjs, real (c1 - 0.5), real (c2 + 0.5), real (l1 - 0.5), real (l2 + 0.5), nmin, maglo, maghi, roundlo, roundhi, sharplo, sharphi) # Print results on the standard output. if (verbose == YES) call sf_write (STDOUT, Memi[cols], Memr[x], Memr[y], Memr[mag], Memi[npix], Memr[size], Memr[ellip], Memr[theta], Memr[sharp], nstars, ct, fmtstr, ntotal + stid) # Save the results in the file. call sf_write (out, Memi[cols], Memr[x], Memr[y], Memr[mag], Memi[npix], Memr[size], Memr[ellip], Memr[theta], Memr[sharp], nstars, ct, fmtstr, ntotal + stid) ntotal = ntotal + nstars } # Free space call sfree (sp) return (ntotal) end # SF_DETECT -- Detect 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 sepmin pixels. int procedure sf_detect (density, ncols, sepmin, nxk, nyk, threshold, cols) real density[ncols, ARB] #I the input density enhancements array int ncols #I the x dimension of the input array int sepmin #I the minimum separation in pixels int nxk, nyk #I size of the fitting area real threshold #I density threshold int cols[ARB] #O column numbers of detected stars int i, j, k, ymiddle, nxhalf, nyhalf, ny, b2, nobjs, rj2, r2 define nextpix_ 11 begin ymiddle = 1 + nyk / 2 + sepmin nxhalf = nxk / 2 nyhalf = nyk / 2 ny = 2 * sepmin + 1 b2 = sepmin ** 2 # Loop over all the columns in an image line. nobjs = 0 for (i = 1 + nxhalf + sepmin; i <= ncols - nxhalf - sepmin; ) { # Test whether the density enhancement is above threshold. if (density[i,ymiddle] < threshold) goto nextpix_ # Test whether a given density enhancement satisfies the # separation criterion. do j = 1, ny { rj2 = (j - sepmin - 1) ** 2 do k = i - sepmin, i + sepmin { r2 = (i - k) ** 2 + rj2 if (r2 <= b2) { if (density[i,ymiddle] < density[k,j+nyhalf]) goto nextpix_ } } } # 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 + sepmin. i = i + sepmin nextpix_ # Work on the next pixel. i = i + 1 } return (nobjs) end # SF_MOMENTS -- Perform a moments analysis on the dectected objects. procedure sf_moments (data, den, ncols, skip, nxk, nyk, cols, x, y, npix, mag, size, ellip, theta, sharp, nobjs, datamin, datamax, threshold, hwhmpsf, xoff, yoff) real data[ncols,ARB] #I the input data array real den[ncols,ARB] #I the input density enhancements array int ncols #I the x dimension of the input buffer int skip[nxk,ARB] #I the input fitting array int nxk, nyk #I the dimensions of the fitting array int cols[ARB] #I the input initial positions real x[ARB] #O the output x coordinates real y[ARB] #O the output y coordinates int npix[ARB] #O the output area in number of pixels real mag[ARB] #O the output magnitude estimates real size[ARB] #O the output size estimates real ellip[ARB] #O the output ellipticity estimates real theta[ARB] #O the output position angle estimates real sharp[ARB] #O the output sharpness estimates int nobjs #I the number of objects real datamin, datamax #I the minium and maximum good data values real threshold #I threshold for moments computation real hwhmpsf #I the HWHM of the PSF real xoff, yoff #I the x and y coordinate offsets int i, j, k, xmiddle, ymiddle, sumn double pixval, sumix, sumiy, sumi, sumixx, sumixy, sumiyy, r2, dx, dy, diff double mean begin # Initialize xmiddle = 1 + nxk / 2 ymiddle = 1 + nyk / 2 # Compute the pixel sum, number of pixels, and the x and y centers. do i = 1, nobjs { # Estimate the background using the input data and the # best fitting Gaussian amplitude sumn = 0 sumi = 0.0 do j = 1, nyk { do k = 1, nxk { if (skip[k,j] == NO) next pixval = data[cols[i]-xmiddle+k,j] if (pixval < datamin || pixval > datamax) next sumi = sumi + pixval sumn = sumn + 1 } } if (sumn <= 0) mean = data[cols[i],ymiddle] - den[cols[i],ymiddle] else mean = sumi / sumn # Compute the first order moments. sumi = 0.0 sumn = 0 sumix = 0.0d0 sumiy = 0.0d0 do j = 1, nyk { do k = 1, nxk { if (skip[k,j] == YES) next pixval = data[cols[i]-xmiddle+k,j] if (pixval < datamin || pixval > datamax) next pixval = pixval - mean if (pixval <= 0.0) next sumi = sumi + pixval sumix = sumix + (cols[i] - xmiddle + k) * pixval sumiy = sumiy + j * pixval sumn = sumn + 1 } } # Use the first order moments to estimate the positions # magnitude, area, and amplitude of the object. if (sumi <= 0.0) { x[i] = cols[i] y[i] = (1.0 + nyk) / 2.0 mag[i] = INDEFR npix[i] = 0 } else { x[i] = sumix / sumi y[i] = sumiy / sumi mag[i] = -2.5 * log10 (sumi) npix[i] = sumn } # Compute the second order central moments using the results of # the first order moment analysis. sumixx = 0.0d0 sumiyy = 0.0d0 sumixy = 0.0d0 do j = 1, nyk { dy = j - y[i] do k = 1, nxk { if (skip[k,j] == YES) next pixval = data[cols[i]-xmiddle+k,j] if (pixval < datamin || pixval > datamax) next pixval = pixval - mean if (pixval <= 0.0) next dx = cols[i] - xmiddle + k - x[i] sumixx = sumixx + pixval * dx ** 2 sumixy = sumixy + pixval * dx * dy sumiyy = sumiyy + pixval * dy ** 2 } } # Use the second order central moments to estimate the size, # ellipticity, position angle, and sharpness of the objects. if (sumi <= 0.0) { size[i] = 0.0 ellip[i] = 0.0 theta[i] = 0.0 sharp[i] = INDEFR } else { sumixx = sumixx / sumi sumixy = sumixy / sumi sumiyy = sumiyy / sumi r2 = sumixx + sumiyy if (r2 <= 0.0) { size[i] = 0.0 ellip[i] = 0.0 theta[i] = 0.0 sharp[i] = INDEFR } else { size[i] = sqrt (LN_2 * r2) sharp[i] = size[i] / hwhmpsf diff = sumixx - sumiyy ellip[i] = sqrt (diff ** 2 + 4.0d0 * sumixy ** 2) / r2 if (diff == 0.0d0 && sumixy == 0.0d0) theta[i] = 0.0 else theta[i] = RADTODEG (0.5d0 * atan2 (2.0d0 * sumixy, diff)) if (theta[i] < 0.0) theta[i] = theta[i] + 180.0 } } # Convert the computed coordinates to the image system. x[i] = x[i] + xoff y[i] = y[i] + yoff } end # SF_TEST -- Check that the detected objects are in the image, contain # enough pixels above background to be measurable objects, and are within # the specified magnitude, roundness and sharpness range. int procedure sf_test (cols, x, y, npix, mag, size, ellip, theta, sharps, nobjs, c1, c2, l1, l2, nmin, maglo, maghi, roundlo, roundhi, sharplo, sharphi) int cols[ARB] #U the column ids of detected object real x[ARB] #U the x position estimates real y[ARB] #U the y positions estimates int npix[ARB] #U the area estimates real mag[ARB] #U the magnitude estimates real size[ARB] #U the size estimates real ellip[ARB] #U the ellipticity estimates real theta[ARB] #U the position angle estimates real sharps[ARB] #U sharpness estimates int nobjs #I the number of detected objects real c1, c2 #I the image column limits real l1, l2 #I the image line limits int nmin #I the minimum area real maglo, maghi #I the magnitude limits real roundlo, roundhi #I the roundness limits real sharplo, sharphi #I the sharpness limits int i, nstars begin # Loop over the detected objects. nstars = 0 do i = 1, nobjs { if (x[i] < c1 || x[i] > c2) next if (y[i] < l1 || y[i] > l2) next if (npix[i] < nmin) next if (mag[i] < maglo || mag[i] > maghi) next if (ellip[i] < roundlo || ellip[i] > roundhi) next if (! IS_INDEFR(sharps[i]) && (sharps[i] < sharplo || sharps[i] > sharphi)) next # Add object to the list. nstars = nstars + 1 cols[nstars] = cols[i] x[nstars] = x[i] y[nstars] = y[i] mag[nstars] = mag[i] npix[nstars] = npix[i] size[nstars] = size[i] ellip[nstars] = ellip[i] theta[nstars] = theta[i] sharps[nstars] = sharps[i] } return (nstars) end # SF_WRITE -- Write the results to the output file. procedure sf_write (fd, cols, x, y, mag, npix, size, ellip, theta, sharp, nstars, ct, fmtstr, stid) int fd #I the output file descriptor int cols[ARB] #I column numbers real x[ARB] #I xcoords real y[ARB] #I y coords real mag[ARB] #I magnitudes int npix[ARB] #I number of pixels real size[ARB] #I object sizes real ellip[ARB] #I ellipticities real theta[ARB] #I position angles real sharp[ARB] #I sharpnesses int nstars #I number of detected stars in the line pointer ct #I coordinate transformation char fmtstr[ARB] #I the output format string int stid #I output file sequence number double lx, ly, wx, wy int i begin if (fd == NULL) return do i = 1, nstars { call fprintf (fd, fmtstr) call pargr (x[i]) call pargr (y[i]) if (ct != NULL) { lx = x[i] ly = y[i] call mw_c2trand (ct, lx, ly, wx, wy) call pargd (wx) call pargd (wy) } call pargr (mag[i]) call pargi (npix[i]) call pargr (size[i]) call pargr (ellip[i]) call pargr (theta[i]) call pargr (sharp[i]) call pargi (stid + i - 1) } end