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