PRO findstars, imdata, bpmdata, sky, thresh, fwhm, sharplim, roundlim, xc, yc, flux, sharp, round, nstars, QUIET = quiet
; Description: This module is based on the DAOFind routine by Peter B. Stetson (Stetson P.B.,
; 1987, PASP, 99, 191). The program detects stars, assumed to be represented
; by circularly symmetric PSFs, by convolving the input image "imdata" with a
; lowered Gaussian of FWHM "fwhm" pixels, and finding all peaks corresponding to
; positive brightness features in the convolved detection image that have values
; larger than "thresh". This is equivalent to finding all peaks corresponding to
; positive brightness features in the input image that have values larger than
; "sky + thresh". Further filtering of the detected objects to protect against
; spurious non-stellar detections is facilitated by the calculation of the "sharp"
; and "round" statistics, which are filtered according to the limits specified by
; "sharplim" and "roundlim". The centroids "xc" and "yc" of all detected stars are
; calculated using the method of thresholded moments, and they are returned along
; with crude flux estimates "flux", the "sharp" and "round" statistics, and the
; number of detected stars "nstars". Note that bad pixels are ignored using the
; bad pixel mask "bpmdata".
; This code was developed directly from the IDL code for the "find.pro" routine
; in the IDL Astronomy Library by Wayne Landsman. The improvements in this code lie
; in the correct treatment of the bad pixel mask, the removal of objects too close
; to the image edge and/or bad pixels (which was causing skewed centroid
; measurements for such objects), and the use of thresholded moments for the
; centroid calculations.
;
; Input Parameters:
;
; imdata - FLOAT/DOUBLE ARRAY - An image array within which stars are to be detected.
; bpmdata - BYTE/INTEGER/LONG ARRAY - A bad pixel mask array of the same size as "imdata"
; which flags bad pixels with a value of "0" and good
; pixels with a value of "1". If this array does not have
; the same dimensions as "imdata", then all image pixels
; will be considered good, and "bpmdata" will be set to an
; array of the same dimensions as "imdata" with all values
; equal to "1".
; sky - FLOAT/DOUBLE - A reliable estimate of the sky background level in the image (ADU).
; thresh - FLOAT/DOUBLE - All local peaks that have values larger than "thresh" (ADU) in the
; convolved detection image are considered as candidate star
; detections. This is equivalent to considering all local peaks that
; have values larger than "sky + thresh" (ADU) in the input image.
; fwhm - FLOAT/DOUBLE - A reliable estimate of the PSF FWHM (pix).
; sharplim - FLOAT/DOUBLE VECTOR - A two-element vector of numbers specifying the lower and
; upper limits (in this order) for the acceptable values of
; the "sharp" statistic calculated for candidate star
; detections (Default value = [0.2, 1.0]). The sharp
; statistic can take any non-negative value.
; roundlim - FLOAT/DOUBLE VECTOR - A two-element vector of numbers specifying the lower and
; upper limits (in this order) for the acceptable values of
; the "round" statistic calculated for candidate star
; detections (Default value = [-1.0, 1.0]). The round
; statistic can take values in the range -2.0 to 2.0
; inclusive.
;
; Output Parameters:
;
; xc - DOUBLE VECTOR - A vector of DOUBLES of length "nstars" specifying the x coordinate
; of each detected star (pix).
; yc - DOUBLE VECTOR - A vector of DOUBLES of length "nstars" specifying the y coordinate
; of each detected star (pix).
; flux - DOUBLE VECTOR - A vector of DOUBLES of length "nstars" specifying the flux estimate
; for each detected star (ADU).
; sharp - FLOAT VECTOR - A vector of FLOATS of length "nstars" specifying the "sharp"
; statistic for each detected star.
; round - FLOAT VECTOR - A vector of FLOATS of length "nstars" specifying the "round"
; statistic for each detected star.
; nstars - LONG - The number of star detections.
;
; Keywords:
;
; If the keyword QUIET is set (as "/QUIET"), then the program will not print any progress
; messages on screen. This functionality is useful if the module is to be called many times
; within an image processing routine for instance.
;
; Author: Dan Bramich (dan.bramich@hotmail.co.uk)
;
; History:
;
; 24/04/2009 - Further module optimisation performed (dmb).
; 17/07/2008 - Module created (dmb).
;
; This module does the following, treating bad pixels correctly at all stages:
;
; 1) Checks that the input image array "imdata" is an image of the correct number type, and
; determines its dimensions. If "imdata" is not an image of the correct number type, then
; the module returns without doing anything.
;
; 2) Checks that "bpmdata" has the same dimensions as "imdata" and is of the correct number type.
; If this is not the case, then the module creates a bad pixel mask that flags every pixel in
; "imdata" as good, using this new bad pixel mask to replace the values stored in "bpmdata".
;
; 3) Checks the values of some of the input parameters, resetting them to default values where
; necessary, and returning without doing anything when parameter values are out of range.
;
; 4) Prepares the convolution kernel for creating the detection image as a lowered Gaussian of
; FWHM given by the parameter "fwhm". The convolution kernel is a pixel grid in the shape of
; a circle with a radius given approximately by "(0.64*fwhm) + 0.5" and limited to a
; minimum value of 2.5 and a maximum value of 12.5.
;
; 5) Convolves the input image with the detection convolution kernel, growing the bad pixel mask
; by the size and shape of the convolution kernel, and setting the bad pixel mask to zero
; within a range of the kernel size from any image edge.
;
; 6) Finds all peaks in the detection image with peak values above the threshold "thresh",
; and that are local maxima within the shape of the convolution kernel. If there are no peaks
; satisfying these constraints, then the module exits with zero star detections. If there
; is at least one peak satisfying these constraints, then the module considers these peaks
; as star candidates and moves on to step (7).
;
; 7) Removes those star candidates with peaks in the detection image that are within a distance
; of the kernel size from the image edge, or that have at least one bad pixel within the kernel
; area around the peak in the detection image. If there are no star candidates left, then the
; module exits with zero star detections. If there is at least one star candidate left, then
; the module considers these star candidates further by moving on to step (8).
;
; 8) For the remaining star candidates, the module calculates the "sharp" statistic, the "round"
; statistic, and attempts to calculate the thresholded centroid. Star candidates are removed
; if either of the "sharp" or "round" statistics do not fall within the specified lower and
; upper limits in the parameters "sharplim" and "roundlim", or if the module was unable to
; calculate a thresholded centroid.
;
; 9) If there are no star candidates left, then the module exits with zero star detections. If
; there is at least one star candidate left, then the output vectors "xc", "yc", "flux",
; "sharp" and "round", each of length "nstars", are returned.
;If the keyword QUIET is set, then do not print progress messages
;Check that "imdata" is an image of the correct number type
;If "bpmdata" does not have the same dimensions as "imdata", then assume that all image pixels
;are good
;Check the parameters "sky", "thresh" and "fwhm"
;Check the parameters "sharplim" and "roundlim", reverting to the default values if the
;input values do not make sense
;Prepare the convolution kernel for creating the detection image
;Convolve the input image with the detection convolution kernel
;Find all peaks in the detection image with peak values above the threshold "thresh",
;and that are local maxima within the shape of the convolution kernel
;Remove detected peaks in the detection image that are within a distance of the kernel size
;from the edge, or that have at least one bad pixel within the kernel area around the peak
;Remove detected peaks in the detection image that do not satisfy both the "sharp"
;and "round" criteria, and for which a thresholded centroid cannot be calculated
;Prepare the convolution kernel for calculating the "round" statistic
;For each detected peak
;Calculate the "sharp" statistic
;Continue to calculate the "round" statistic if the value of the "sharp" statistic is acceptable
;Extract one-dimensional pixel distributions in the x and y directions
;Calculate the heights of the fitted one-dimensional Gaussian profile to the one-dimensional
;pixel distributions using the previously prepared convolution kernel
;Continue to calculate the "round" statistic if both fitted heights are positive
;Calculate the "round" statistic
;Continue to determine the peak centroid if the value of the "round" statistic is acceptable
;Calculate the x and y centroids of the current peak. The threshold chosen when calculating
;the thresholded centroid is 0.3 times the object detection threshold because the image
;cutout used was chosen to have a half side of 1.5 sigma.
;If the centroid calculation was successful, then save the current peak