PRO choose_psf_stars, xc, yc, flux, peak, corr, range, imxsize, imysize, bpm, det_thresh, corr_thresh, $
; Description: This module chooses a set of suitable PSF stars from an input star list containing
; coordinates "xc" and "yc", flux estimates "flux", local peak pixel values "peak" and
; correlation coefficients with the image PSF "corr". Suitable PSF stars are chosen using
; the following selection criteria:
;
; (i) A PSF star must not be within a distance of "range" along each coordinate axis from
; the image edge.
; (ii) A PSF star must have no bad pixels within a distance of "range" along each
; coordinate axis from its centre.
; (iii) A PSF star must have a local peak pixel value that is greater than or equal to
; the threshold "det_thresh".
; (iv) A PSF star must have a correlation coefficient with the image PSF that is greater
; than or equal to the threshold "corr_thresh".
; (v) A PSF star must not have a brighter star within a distance of "range" from its
; centre.
; (vi) A PSF star must not have a companion star with a flux greater than or equal to
; "psf_comp_flux" times its own flux that is within a distance of
; "psf_comp_dist*range" from its centre.
;
; The chosen PSF stars are returned to the user via a set of subscripts "psf_subs"
; corresponding to the star list.
; The input star list should contain all non-saturated stars down to below the brightness
; threshold for the PSF stars "det_thresh", in order to properly reject candidate PSF stars
; based on the companion star selection criteria (v) and (vi).
;
; Input Parameters:
;
; xc - FLOAT/DOUBLE VECTOR - A one-dimensional vector of x coordinates for the input star list (pix).
; yc - FLOAT/DOUBLE VECTOR - A one-dimensional vector of y coordinates for the input star list (pix).
; flux - FLOAT/DOUBLE VECTOR - A one-dimensional vector of flux estimates for the input star list (ADU).
; peak - FLOAT/DOUBLE VECTOR - A one-dimensional vector of local peak image pixel values for the input
; star list (ADU). If this vector does not have the same number of
; elements as the parameter "xc", then the module assumes that all stars
; from the input star list have a local peak image pixel value that is
; equal to "det_thresh", which means that all stars automatically pass
; the selection criterion (iii) described above.
; corr - FLOAT/DOUBLE VECTOR - A one-dimensional vector of correlation coefficients between each star
; from the input star list and the image PSF. If this vector does not have
; the same number of elements as the parameter "xc", then the module
; assumes that all stars from the input star list have a correlation
; coefficient of "1", which means that all stars automatically pass the
; selection criterion (iv) described above.
; range - FLOAT/DOUBLE - The distance to be used to test for the presence of bad pixels, proximity to
; the image edge, brighter stars and companion stars. If this parameter is set
; to zero, or a negative value, then the module will not apply the selection
; criteria (i), (ii), (v) and (vi) described above.
; imxsize - INTEGER/LONG - The image size along the x-axis (pix). If this parameter is set to zero, or
; a negative value, then the module will not apply the selection criteria
; (i) and (ii) described above.
; imysize - INTEGER/LONG - The image size along the y-axis (pix). If this parameter is set to zero, or
; a negative value, then the module will not apply the selection criteria
; (i) and (ii) described above.
; bpm - BYTE/INTEGER/LONG ARRAY - A two-dimensional bad pixel mask array of size "imxsize" by "imysize"
; pixels, which flags bad pixels with a value of "0" and good pixels
; with a value of "1". If this array does not have the correct
; dimensions and/or size, then the module will consider all image
; pixels to be good, which means that all stars automatically pass the
; selection criterion (ii) described above.
; det_thresh - FLOAT/DOUBLE - The brightness threshold used to select PSF stars as described in the
; selection criterion (iii) above (Default value = 0.0). Negative values
; will be reset to the default value.
; corr_thresh - FLOAT/DOUBLE - The correlation coefficient threshold used to select PSF stars as
; described in the selection criterion (iv) above (Default value = 0.0).
; Negative values, or values of "1" or greater than "1", will be reset to
; the default value.
; psf_comp_dist - FLOAT/DOUBLE - The distance, as a fraction of the parameter "range", within which a
; star is considered to be a companion star for use in the PSF selection
; criterion (vi) described above (Default value = 0.0). Negative values
; will be reset to the default value.
; psf_comp_flux - FLOAT/DOUBLE - The fractional flux above which a companion star will discount a
; candidate PSF star as described in the selection criterion (vi) above
; (Default value = 0.0). Negative values will be reset to the default
; value.
;
; Output Parameters:
;
; psf_subs - LONG VECTOR - A one-dimensional vector of subscripts from the input star list corresponding
; to the chosen PSF stars.
; npsf_subs - LONG - The number of chosen PSF stars, which happens to also be the number of elements in
; the output parameter "psf_subs".
;
; Keywords:
;
; If the keyword SORT_BY_FLUX is set (as "/SORT_BY_FLUX"), then the program will sort the subscripts
; "psf_subs" such that the chosen PSF stars are ordered from the brightest to the faintest.
;
; Author: Dan Bramich (dan.bramich@hotmail.co.uk)
;
; History:
;
; 24/04/2009 - Further module optimisation performed (dmb).
; 27/08/2008 - Module created (dmb).
;Check that "xc" is a vector containing number data
;Check that "yc" is a vector containing number data with the same number of elements as "xc"
;Check that "flux" is a vector containing number data with the same number of elements as "xc"
;Check the values of the parameters "det_thresh", "corr_thresh", "psf_comp_dist" and "psf_comp_flux", and
;reset them to the default values if necessary
;Check that "peak" is a vector containing number data with the same number of elements as "xc". If
;this is not the case, then prepare a vector of peak pixel values that are all set to "det_thresh".
;Check that "corr" is a vector containing number data with the same number of elements as "xc". If
;this is not the case, then prepare a vector of correlation coefficients that are all set to "1.0"
;Check that "range", "imxsize" and "imysize" are all numbers
;Check that "bpm" is a two-dimensional byte/integer/long array of size "imxsize" by "imysize" pixels. If this is
;not the case then ignore the bad pixel mask from now on.
;Create an array of subscripts for the input star list in order to index candidate PSF stars
;Remove candidate PSF stars that do not have a local peak pixel value that is greater than or equal to
;the threshold "det_thresh" or that do not have a correlation coefficient with the image PSF that is
;greater than or equal to the threshold "corr_thresh"
;If required, then remove candidate PSF stars that are within a distance of "range" along each coordinate
;axis from the image edge
;If required, then remove candidate PSF stars that have at least one bad pixel within a distance "range"
;along each coordinate axis from the star centre
;If required, then remove candidate PSF stars that have either:
;(a) A brighter companion star within a distance of "range" from the star centre, or,
;(b) A companion star with a flux greater than or equal to "psf_comp_flux" times the star flux that is within
; a distance of "psf_comp_dist*range" from the star centre.
;Make a copy of the input star list, and sort the list such that the x coordinate is in ascending order
;For each candidate PSF star
;Use a fast binary search to find the subset of stars from the input star list that have x coordinates
;within a distance D along the x axis, where D is the larger of "range" and "psf_comp_dist*range"
;Calculate the square of the distance between the current candidate PSF star and each of the stars in
;the subset of stars from the input star list
;Determine the number of stars in the subset of stars from the input star list that are brighter than
;the current candidate PSF star and that are within a distance of "range" from the centre of the
;current candidate PSF star
;Determine the number of stars in the subset of stars from the input star list that are brighter than
;a factor of "psf_comp_flux" times the star flux of the current candidate PSF star and that are within
;a distance of "psf_comp_dist*range" from the centre of the current candidate PSF star
;Save the subscript and flux of the current candidate PSF star
;If there are no candidate PSF stars left, then return with an empty list of PSF stars
;If necessary, update the subscripts and fluxes of the list of candidate PSF stars
;If the the keyword SORT_BY_FLUX is set, then sort the final list of suitable PSF stars such that they are
;ordered from the brightest to the faintest
;Save the subscripts of the final list of suitable PSF stars