; $Id: starfinder.pro, v 1.8.2b_dev01 Sep 2024 e.d. $
;
;+
; NAME:
; STARFINDER
;
; PURPOSE:
; Given a stellar field and an approximation of the Point Spread Function
; (PSF), detect the stellar sources and estimate their position and flux.
; The image may be contaminated by a smooth non-uniform background emission.
; If the imaged field is not isoplanatic, it is possible to partition it
; into approximately isoplanatic sub-regions, provided a local estimate of
; the PSF is available. In this case, the field is analyzed as a whole as
; in the constant PSF case, with the only difference that the local PSF is
; used to analyze each star.
; A spatially variable PSF may also be handled by providing a parametric
; PSF model variable across the imaged field.
;
; CATEGORY:
; Signal processing. Stellar fields astrometry and photometry.
;
; CALLING SEQUENCE:
; STARFINDER, Image, Psf, Threshold, Min_correlation, $
; X, Y, Fluxes, Sigma_x, Sigma_y, Sigma_f, Correlation
;
; INPUTS:
; Image: 2D image of the stellar field
;
; Psf: 2D array, representing the Point Spread Function of the stellar
; field. If the PSF is space-variant, Psf is a 3D stack of local
; PSF measurements, relative to a partition of the imaged field into
; sub-regions arranged in a regular grid. In this case it is necessary
; to supply the bounds of the partition (see keyword SV_PAR).
;
; Threshold: Vector of lower detection levels (above the local background,
; which is temporarly removed before detection). These levels are
; assumed to be "absolute" levels by default. If the keyword
; REL_THRESHOLD (see below) is set and the noise standard deviation is
; defined (see the keyword NOISE_STD) below, then the threshold levels
; are assumed to be "relative", i.e. expressed in units of the noise
; standard deviation.
; Example 1: let Threshold be = [10., 3.]
; If the keyword REL_THRESHOLD is not defined or the keyword
; parameter NOISE_STD is undefined, the threshold levels 10 and 3
; are absolute intensity levels, i.e. all the peaks brighter than 10
; and 3 counts are considered as presumed stars. If REL_THRESHOLD is
; set and NOISE_STD is defined, the intensity levels for detection
; are respectively (10 * NOISE_STD) and (3 * NOISE_STD).
; The number of elements of the vector Threshold establishes the number
; of iterations of the "basic step" (see PROCEDURE below).
; The components of the vector Threshold may be all equal or decreasing.
;
; Min_correlation: Minimum value of correlation between an acceptable
; stellar image and the Psf. It must be > 0 and < 1.
;
; KEYWORD PARAMETERS:
; SV_PAR: Set this keyword to a structure defining the bounds of the
; imaged field partition when the PSF is not isoplanatic.
; Let the imaged field be partitioned into Nx*Ny sub-regions, such
; that the (j, i)-th region is bounded by
; Lx[j] < x < Ux[j], Ly[i] < y < Uy[i], where
; j = 0, ..., Nx - 1; i = 0, ..., Ny - 1.
; The structure must be defined as follows:
; SV_PAR = {Lx: Lx, Ux: Ux, Ly: Ly, Uy: Uy}.
; The (j, i)-th sub-region must correspond to the k-th psf in the 3D
; stack, where k = i*Nx + j.
;
; PSF_TYPE: Set this keyword to a string indicating the name of a PSF
; model. Available PSF models are defined in the documentation of the
; routine IMAGE_MODEL (see documentation of the input parameter Psf in
; the file IMAGE_MODEL.PRO). This keyword has to be used together with
; the keyword PSF_DATA.
; The keywords PSF_TYPE and PSF_DATA are useful to define a parametric
; PSF model, that can be made variable across the field of view,
; depending on how it is defined in IMAGE_MODEL.PRO.
; When PSF_TYPE and PSF_DATA are used, the input parameter Psf and the
; keyword parameter SV_PAR described above have to be supplied anyway,
; as they are used in the preliminary phases of the analysis.
;
; PSF_DATA: Use this keyword to provide a pointer to the auxiliary data
; required by the PSF model defined by the keyword PSF_TYPE, when it is
; used. Auxiliary data are defined according to the documentation of the
; input parameter Data in the file IMAGE_MODEL.PRO.
;
; REL_THRESHOLD: Set this keyword to specify that the detection levels
; contained in the input parameter Threshold have to be considered as
; "relative" levels, i.e. they must be multiplied by the noise standard
; deviation. This keyword is effective only if the NOISE_STD is defined
; (see below).
;
; X_BAD, Y_BAD: Coordinates of bad pixels to be excluded.
; It is important to mask bad pixels, especially in the following phases
; of the analysis:
; 1) search for presumed stars
; 2) correlation check
; 3) local fitting.
;
; BACKGROUND: 2D array, containing an initial guess of the image background
; emission. If the keyword is undefined, the background is estimated
; automatically, provided the box size BACK_BOX (see below) is defined.
; The background estimate is subtracted before searching for presumed
; stars in the image and before computing the correlation coefficient.
;
; BACK_BOX: Set this keyword to a scalar value specifying the box size
; to estimate the background emission. For more details see the routine
; ESTIMATE_BACKGROUND in the file "estimate_background.pro".
; The value of this keyword may be conveniently defined as a multiple
; of the PSF Full Width at Half Maximum (e.g. 5 * FWHM).
; If if is not supplied, the program uses the argument of the keyword
; BACKGROUND (if defined of course!) as a background estimate. However
; this will prevent any improvements in the background knowledge in the
; course of iterations. If the stellar field is supposed to have no
; background emission, the keywords BACKGROUND and BACK_BOX should be
; left undefined.
;
; For other keywords concerning background estimation see the routines
; IMAGE_BACKGROUND and ESTIMATE_BACKGROUND.
;
; FOUR: Set this keyword to identify relative maxima referring only to
; the 4-neighbors of each pixel. The default is to use 8-neighbors.
;
; PRE_SMOOTH: Set this keyword to smooth the image with a standard low-pass
; filter before searching for local maxima. This reduces the probability
; to pick up noise spikes, which would increase the computation time for
; further checks and the probability of false detections.
; For more details see the procedure SEARCH_OBJECTS, in the file
; "search_objects.pro".
;
; MINIF: If PRE_SMOOTH is set, the keyword MINIF may be used to enter the
; integer minification factor used to down-sample the image before
; searching for local maxima, as an alternative smoothing strategy to
; the standard low-pass filter used by default when PRE_SMOOTH is set.
;
; CORREL_MAG: Set this keyword to an integer scalar specifying the sub-pixel
; accuracy with which the template PSF should be superposed to a given
; object when computing the correlation coefficient.
; If CORREL_MAG = 1, the PSF is superposed to the
; object just taking the maximum intensity pixel as a reference.
; If CORREL_MAG is > 1, all the possible fractional shifts of the PSF,
; with step size of 1/CORREL_MAG, are considered.
; The default value is CORREL_MAG = 2 (half pixel positioning).
;
; DEBLEND: Set this keyword to a nonzero value to perform de-blending
; of the objects detected iterating the basic step (see PROCEDURE).
; This de-blending strategy allows the algorithm to detect secondary
; components of crowded groups, whose principal component has previously
; been detected as a single star. This strategy may be referred to as
; 'partial de-blending'.
;
; DEBLOST: Set this keyword to a nonzero value to recover all the
; significant intensity enhancement lost in the course of the analysis,
; in order to search for possible blends. This strategy allows the
; algorithm to re-analyze the blended objects whose principal component
; has not been previously found as a single star (see DEBLEND above).
; This strategy, applied after the 'partial de-blending' described with
; the keyword DEBLEND above, may be referred to as 'full de-blending'.
; More details can be found in the PROCEDURE description below.
;
; CUT_THRESHOLD: Set this keyword to a scalar value in the range ]0, 1[
; to specify the relative threshold to use when "cutting" a given
; object to measure its area. The default is CUT_THRESHOLD = 0.5,
; i.e. the object is cut at 50% of the peak height. Notice that the
; background is subtracted beforehands. Notice also that if the
; absolute value of the threshold for binary detection is comparable
; to the gaussian noise level (see also keywords GAUSSIAN_NOISE and
; BIN_N_STDEV), the estimated area of the object is not reliable: thus
; deblending is not applied. This check on the threshold value is
; performed only if the gaussian noise standard deviation is defined
; (see the keyword GAUSSIAN_NOISE).
; For more details on the determination of the area, see the routine
; PEAK_AREA in the file "peak_area.pro".
;
; BIN_N_STDEV: Set this keyword to a positive scalar to specify the
; minimum cutting threshold for binary detection, in units of the
; gaussian noise standard deviation. The default value is
; BIN_N_STDEV = 3.
;
; BIN_THRESHOLD: When a crowded group (in the simplest case a close binary
; star) is analyzed, the deblending algorithm searches for residuals
; representing secondary components by subtracting the known sources:
; The brightest local maximum in the residual sub-image is taken as a
; candidate secondary star. In principle the detection threshold used
; for isolated objects, set by the parameter Threshold described above,
; should not be applied in this case, because the residual is affected
; by the poor knowledge about the subtracted components.
; BIN_THRESHOLD specifies the relative threshold, in units of the
; parameter Threshold, to use for the detection of candidate secondary
; sources in crowded groups. The default is BIN_THRESHOLD = 0, i.e. all
; the positive residuals (remember that the local background level has
; been removed!) are considered. Their intensity is compared to the
; detection threshold after fitting: if the fitted intensity happens to;
; be below the value set by Threshold, the new secondary component is
; rejected.
;
; BIN_TOLERANCE: An object is classified as "blend" if the relative error;
; between its area and the area of the PSF is greater than the threshold
; fixed by this keyword. The default is BIN_TOLERANCE = 0.2, i.e. 20%.
;
; NOISE_STD: Set this keyword to a 2D array, with the same size as Image,
; representing the noise standard deviation in each pixel of the input
; data. This array is used to compute the formal errors on astrometry
; and photometry.
;
; NO_SLANT: Set this keyword to avoid fitting the local background with
; a slanting plane. In this case the background to use for the local
; fit of the stars is derived from the 2D background array and it is
; kept fixed in the fitting. This option should be used only if the
; (input) estimate of the image background is very accurate.
;
; MIN_DISTANCE: Minimum distance, expressed in units of the PSF FWHM, to
; consider two sources as separated. This parameter is used to check the
; outcome of the fitting of a group of neighboring sources: if two sources
; are closer than this limit, the fit is considered unacceptable and the
; last detected source is rejected as a false detection.
; The default value is 1 (which means 1 PSF FWHM).
;
; N_ITER: At the end of the analysis, the detected stars are fitted again
; a number of times equal to the value specified with this keyword.
; The default is N_ITER = 1. To avoid re-fitting, set N_ITER = 0.
; Notice that the final number of re-fitting iterations includes also
; the extra re-fitting which is normally performed after each basic
; step (see also the keyword NO_INTERMEDIATE_ITER).
;
; NO_INTERMEDIATE_ITER: The currently detected stars are re-fitted after
; each basic step. Set NO_INTERMEDIATE_ITER to avoid this re-fitting.
;
; ASTROMETRIC_TOL, PHOTOMETRIC_TOL: When the last re-fitting is iterated
; a lot of times, it is possible to stop the iterations when convergence
; on stellar positions and fluxes is achieved. The default values for
; convergence are ASTROMETRIC_TOL = 0.01 pixels and
; PHOTOMETRIC_TOL = 0.01 (i.e. 1%).
; In practice setting N_ITER = 1 or 2 is generally enough to obtain
; reliable astrometry and photometry, without having to iterate a lot.
;
; X_INPUT, Y_INPUT, F_INPUT: vectors with positions and fluxes of the stars,
; when these are already known. When these parameters are defined, no
; search is performed: the known sources are just fitted a pre-fixed number
; of times (see N_ITER) or until convergence (see ASTROMETRIC_TOL and
; PHOTOMETRIC_TOL). The sources are sorted by decreasing flux.
; Usual checks applied in fitting procedure are applied (minimum flux
; above threshold, minimum distance between any two stars).
; When this option is used, single-component fitting is performed, i.e.
; stars closer than the fitting box width are fit one at a time.
; Input sources for which the fit fails are rejected.
; *** Recommendation. When this option is used, the following keyword
; should be set:
; MIN_DISTANCE = 0.0
; In order to avoid optimizing the input positions set the keyword /FLUXOPT
; (see below). In this case, the positions are kept fixed to the input values,
; only fluxes are optimized.
;
; FLUXOPT: Set this binary keyword to keep the positions of the sources
; fixed in the optimization process when an input list of sources is given
; (see keywords X_INPUT, Y_INPUT, F_INPUT above).
;
; C_INPUT: when an input list of sources is given (see keywords X_INPUT,
; Y_INPUT, F_INPUT above), the keyword C_INPUT may be used to provide the
; correlation coefficients of the sources (detemined in a previous detection
; run). These correlation coefficients are just retained on output, associated
; to the sources for which the fit process is successful.
;
; SVDINV: The fitting process is performed by a Newton-Gauss method.
; This method requires a matrix inversion at every iteration.
; Set this keyword to perform this matrix inversion by SVD (Singular
; Value Decomposition). The default is matrix inversion by Moore-Penrose
; Generalized Inverse (see routine 'ginv.pro').
;
; SILENT: Set this keyword to avoid printing messaged to the output window.
; If not set, STARFINDER prints messages about the current analysis step and
; about the candidate and detected stars. "Candidate stars" are the peaks
; found in the image; "Detected stars" are the peaks which are confirmed to be
; stars in the analysis. The number of "Candidate stars" does not account for
; the additional candidate peaks found with the deblending procedure.
;
; OUTPUTS:
; X, Y: Coordinates of detected stars in pixel units.
; The origin is at (0, 0).
;
; Fluxes: Vector of stellar fluxes, referred to the normalization of the Psf.
;
; Sigma_x, Sigma_y, Sigma_f: Vector of formal errors on positions and
; fluxes, estimated by the fitting procedure.
; These vectors have all components equal to zero if no information is
; supplied about noise: in this case it has no meaning to estimate the
; error propagation in the fit.
;
; Correlation: Vector of correlation coefficients for accepted stars.
; Stars belonging to crowded groups detected with the de-blending
; procedure have a correlation coefficient of -1.
;
; OPTIONAL OUTPUTS:
; BACKGROUND: Last estimate of the image background, computed by the
; program after the iteration corresponding to the lowest
; detection level.
;
; STARS: 2D array, containing one shifted scaled replica of the PSF for
; each detected star.
;
; SIDE EFFECTS:
; If the keyword NOISE_STD is set to a scalar value it is transformed into
; a 2D constant array with the same size as the input image.
;
; RESTRICTIONS:
; 1) The algorithm is based on the assumption that the input Image may
; be considered as a superpositions of stellar sources and a smooth
; background.
; 2) Saturated stars should have been approximately repaired beforehands.
; This is especially important for the reliable detection and analysis of
; fainter sources in their surrounding.
; 3) High accuracy astrometry and photometry are based on sub-pixel
; positioning of the input PSF. In this version of the algorithm this
; relies on interpolation techniques, which are not suited to undersampled
; data. Thus the algorithm, at least in its present form, should not be
; used with poorly sampled data.
; 4) The keyword /FLUXOPT should be used only when a list of input sources
; is given and the input position estimates are very accurate.
;
; PROCEDURE:
; The standard analysis is accomplished as a sequence of "basic steps".
; One of these steps consists of 3 phases:
; 1) detection of presumed stars above a given threshold
; 2) check and analysis of detected objects, sorted by decreasing
; intensity
; 3) re_fitting
; Pixel (x,y) is the approximate center of a presumed star if
; a) Image(x,y) is a local maximum
; b) Image(x,y) > Background(x,y) + Stars(x,y) + Threshold,
; where Image is the stellar field, Background is an approximation of
; the background emission, Stars is a sum of shifted scaled replicas
; of the PSF (one for each detected star) and Threshold is the
; minimum central intensity of a detectable star. In practice the user
; provides a set of Threshold levels (generally decreasing!): the
; number of these levels fixes the number of times the basic step is
; repeated.
; At the first iteration of the basic step, the "image model" called
; Stars is identically zero; in later iterations it contains the detected
; sources. Subtraction of these stars simplifies the detection of fainter
; objects and allows a better estimation of the image Background.
; It should be noticed that the model of detected stars (Stars) is
; subtracted only in order to detecte new presumed stars and refine the
; background, whereas the subsequent astrometric and photometric analysis
; is performed on the original Image.
; Let us describe in detail how the presumed stars are analyzed (Phase
; 2 of each basic step).
; First of all each object in the newly formed list is "re-identified",
; i.e. searched again after subtraction of other brighter sources: this
; allows to reject many spurious detections associated to PSF features
; of bright stars. The object is compared with the PSF by a correlation
; check and accepted only if a sufficient similarity is found. Then it is
; fitted to estimate its position and flux. The local fitting is performed
; by FITSTARS (see the file "fitstars.pro" for more details), which takes
; into consideration all the following contributions:
; a) other known stars in the fitting box, which are fitted along with
; the object under examination (multiple fitting)
; b) known stars lying outside the fitting box: this term is extracted
; from the image model Stars and considered as a fixed contribution
; c) the local background, which may be approximated either with a
; slanting plane (whose coefficients are optimized as well) or with a
; sub-image extracted from the image Background (kept fixed).
; If the fit is successful, the parameters of the stars are saved and
; the image model Stars is updated.
; The basic step (phases 1 + 2 + 3) may be iterated: a new list of
; objects is formed by searching in the image after subtraction of the
; previously detected stars. Then the analysis proceeds on the original
; frame. This iteration is very useful to detect hidden objects, e.g.
; close binaries down to separations comparable to 1 PSF FWHM.
; An optional deblending strategy is available, consisting of two
; separate phases. The first step (see keyword DEBLEND) consists of a search
; for secondary sources around the objects detected in the earlier phases of
; the analysis; in the second step (see keyword DEBLOST) all the significant
; intensity enhancement previously discarded are recovered and a check is
; performed in order to find among these some possible blends. Blends are
; recognized on the basis of their larger extension as compared to the PSF.
; The area of a given object is measured by a thresholding technique, applied
; at a specified threshold below the peak. When an object is recognized as a
; blend, deblending is attempted iteratively by searching for a new residual
; and subsequent fitting. The iteration is stopped when no more residual is
; found or the fitting of the last residual is not successful. The deblending
; procedure is applied at the end of the last iteration of the basic step
; described above, i.e. when all the "resolved" objects are supposed to have
; been detected.
; It should be stressed that step 3 (re-fitting) has the goal to improve
; the astrometric and photometric measurement of the detected stars. If the
; fit of a given source fails, however, the source (which has already been
; identified as a true one) is not rejected: its parameters are simply
; set to their previous value, without any update.
;
; MODIFICATION HISTORY:
; Written by: Emiliano Diolaiti, August-September 1999.
; Updates:
; 1) New partial deblending strategy (Emiliano Diolaiti, October 1999).
; 2) Full deblending strategy (Emiliano Diolaiti, February 2000).
; 3) Added PSF stack option; this version might not handle properly
; very close groups of objects lying on the boundary of two adjacent
; sub-regions (Emiliano Diolaiti, February 2000).
; 4) Modified keyword parameters in auxiliary routines to avoid trouble
; under IDL 5.0 (Emiliano Diolaiti, May 2000).
; 5) Modified STARFINDER_FIT module. Now faster.
; (Emiliano Diolaiti, December 2000).
; 6) Moved auxiliary routines to handle the list of stars in the
; program module STARLIST. (Emiliano Diolaiti, June 2001).
; 7) Modified STARFINDER_DEBLEND module. Now faster.
; (Emiliano Diolaiti, June 2001).
; 8) PSF stack option (see modification no. 3): fixed problem with stars
; lying close to boundary between adjacent sub-regions
; (E.D., December 2004). DEBLEND option not checked with PSF stack option.
; 9) Increased minimum distance between two stars to 1 FWHM.
; 10) Increased box width for distance check in STARFINDER_CHECK module
; (after suggestion by Jessica Lu, May 2005).
; 11) Temporarily removed modification no. 10 and restored previous version
; for trouble at image edge.
; 12) Added keywords X_INPUT, Y_INPUT, F_INPUT for fitting a set of known
; stars. Not yet tested with space-variant PSF (E. D., February 2006).
; 13) Added keyword MIN_DISTANCE, to set minimum acceptable distance between
; two stars (E. D., August 2006).
; 14) Fixed parameter value: fitting box size was not properly set with
; keywords X_INPUT, Y_INPUT, F_INPUT (E. D., March 2012).
; 15) When an input list of sources is given (keywords X_INPUT, Y_INPUT,
; F_INPUT), the sources for which the fit fails are rejected
; (E. D., March 2012).
; 16) Keyword "re_fitting" turned into a positional parameters in modules
; STARFINDER_FIT, STARFINDER_DEBLEND, STARFINDER_ANALYZE (E. D., March 2012).
; 17) Fixed bug about keyword MIN_DIST: 0 value was not accepted before
; (E. D., March 2012).
; 18) Added keyword SVDINV (E. D., March 2012).
; 19) Changed some output messages in VERBOSE mode (E.D., March 2012).
; 20) Created separate source files for all auxiliary routines described
; at point 6) (E. D., March 2012).
; 21) Changed default value of keyword CUT_THRESHOLD = 0.5 (E.D., March 2012).
; 22) Added keyword MAG=3 in call to PEAK_WIDTH in module
; STARFINDER_DEB_PAR (E.D., March 2012).
; 23) Added keyword FLUXOPT (E.D., April 2012).
; 24) Removed check MAX_DIST in module STARFINDER_DEBLEND (E.D., April 2012).
; 25) Set default value of keyword CORREL_MAG = 2 (E.D., April 2012).
; 26) Removed bug in STARFINDER_FIT module: sometimes very close sources
; were swapped by call to COMPARE_LISTS (E.D., April 2012).
; 27) Updated documentation (E. D., May 2014).
; 28) Modified data type of variable "list_of_stars":
; from array of structures to structure of arrays (E.D., September 2024).
;-
;;; Auxiliary procedures/functions to handle list of stars data structure.
;FUNCTION star
;
;PRO update_list, list, SUBSCRIPTS = s, x, y, f, c, $
; sigma_x, sigma_y, sigma_f, IS_STAR = is_star
;
;FUNCTION create_element, x, y, f
;
;FUNCTION merge_list, l1, l2
;
;FUNCTION add_subscript, subscripts, s
;
;FUNCTION delete_element, list
;
;PRO star_param, list, SUBSCRIPTS = s, $
; n, x, y, f, c, sigma_x, sigma_y, sigma_f
;
;FUNCTION where_stars, list, LX = lx, UX = ux, LY = ly, UY = uy, n
;
;FUNCTION extract_stars, list, n
;
;FUNCTION sort_list, list, SUBSCRIPTS = s
;
;FUNCTION reverse_class, list, SUBSCRIPTS = s
;
;FUNCTION starlist, n, x, y, f, c
;
;;; Auxiliary procedures/functions to define the program's parameters.
; STARFINDER_FWHM: compute PSF FWHM.
FUNCTION starfinder_fwhm, psf, n_psf
on_error, 2
fw = fltarr(n_psf)
for n = 0L, n_psf - 1 do fw[n] = fwhm(psf[*,*,n], /CUBIC, MAG = 3)
return, fw
end
; DEFINE_CORR_BOX: define correlation box size.
FUNCTION define_corr_box, psf_fwhm
on_error, 2
box = round(1.5 * psf_fwhm)
box = box + 1 - box mod 2
return, box
end
; DEFINE_FIT_BOX: define fitting box size.
FUNCTION define_fit_box, psf_fwhm, WIDER = wider
on_error, 2
if keyword_set(wider) then b = 2. else b = 1.
box = round(define_corr_box(psf_fwhm) > 5 + b * psf_fwhm)
box = box + 1 - box mod 2
return, box
end
; STARFINDER_ID_PAR: define parameters to re-identify presumed stars.
FUNCTION starfinder_id_par, psf_fwhm, FOUR = four
on_error, 2
box = define_fit_box(psf_fwhm)
neigh4 = keyword_set(four) and 1B
return, {box: box, neigh4: neigh4}
end
; STARFINDER_CORR_PAR: define parameters for correlation.
FUNCTION starfinder_corr_par, psf, psf_fwhm, n_psf, CORREL_MAG = correl_mag, $
_EXTRA = extra
on_error, 2
correlation_box = define_corr_box(psf_fwhm)
search_box = correlation_box / 2
search_box = search_box + 1 - search_box mod 2
if n_elements(correl_mag) eq 0 then correl_mag = 2
if correl_mag gt 1 then edge = 2 else edge = 0
temp_siz = correlation_box + 2 * edge
template = ptrarr(n_psf, /ALLOCATE)
if correl_mag gt 1 then templates = ptrarr(n_psf, /ALLOCATE) $
else begin
templates = 0 & dx = 0 & dy = 0
endelse
for n = 0L, n_psf - 1 do begin
temp = sub_array(psf[*,*,n], temp_siz[n])
if correl_mag gt 1 then begin
shifted_templates, temp, correl_mag, _EXTRA = extra, EDGE = edge, $
templs, dx, dy
temp = sub_array(temp, correlation_box[n])
*templates[n] = templs
endif
*template[n] = temp
endfor
return, {correl_mag: correl_mag, $
correlation_box: correlation_box, $
search_box: search_box, $
template: template, $
templates: templates, dx: dx, dy: dy}
end
; STARFINDER_FIT_PAR: define parameters for fitting.
FUNCTION starfinder_fit_par, psf, psf_fwhm, n_psf, MIN_DISTANCE = min_dist, _EXTRA = extra
on_error, 2
fitting_box = define_fit_box(psf_fwhm, _EXTRA = extra)
psf_size = 2 * max(fitting_box) < (max((size52(psf, /DIM))[0:1]) - 2)
; edge = round((fitting_box - psf_fwhm) / 2.0) > 1
edge = 0.5
if n_elements(min_dist) eq 0 then min_dist = 1.0
min_distance = min_dist * psf_fwhm
fitting_psf = fltarr(psf_size, psf_size, n_psf)
psf_max = fltarr(n_psf)
for n = 0L, n_psf - 1 do begin
fitting_psf[*,*,n] = sub_array(psf[*,*,n], psf_size)
psf_max[n] = max(psf[*,*,n])
endfor
return, {fitting_box: fitting_box, edge: edge, $
min_distance: min_distance, $
fitting_psf: fitting_psf, psf: psf, $
psf_max: psf_max, psf_fwhm: psf_fwhm}
end
; STARFINDER_DEB_PAR: define parameters for deblending.
FUNCTION starfinder_deb_par, psf, psf_fwhm, n_psf, $
BIN_THRESHOLD = bin_threshold, $
CUT_THRESHOLD = cut, $
BIN_N_STDEV = n_std, BIN_TOLERANCE = tol
on_error, 2
if n_elements(bin_threshold) eq 0 then bin_threshold = 0.
if n_elements(cut) eq 0 then cut = 0.5 & cut = cut > 0 < 1 ;cut = cut > 0.2 < 1
if n_elements(n_std) eq 0 then n_std = 3
if n_elements(tol) eq 0 then tol = 0.2 & tol = tol > 0 < 1
box = lonarr(n_psf) & area = fltarr(n_psf)
for n = 0L, n_psf - 1 do begin
box[n] = round(3 * peak_width(psf[*,*,n], REL = cut, MAG = 3))
area[n] = peak_area(psf[*,*,n], REL = cut, MAG = 3)
endfor
return, {bin_threshold: bin_threshold, n_std: n_std, $
box: box, cut: cut, tol: tol, psf_area: area}
end
; STARFINDER_DEALLOC: de-allocate heap variables.
PRO starfinder_dealloc, corr_par, fit_data, model_data
on_error, 2
if n_elements(corr_par) eq 0 then return
if (ptr_valid(corr_par.template))[0] then ptr_free, corr_par.template
if (ptr_valid(corr_par.templates))[0] then ptr_free, corr_par.templates
if ptr_valid(fit_data) then ptr_free, fit_data
if ptr_valid(model_data) then ptr_free, model_data
return
end
;;; Other auxiliary procedures/functions.
; STARFINDER_BAD: find bad pixels in the sub-image [lx:ux,ly:uy].
PRO starfinder_bad, x_bad, y_bad, lx, ux, ly, uy, x_bad_here, y_bad_here
on_error, 2
w = where(x_bad ge lx and x_bad le ux and $
y_bad ge ly and y_bad le uy, n)
if n ne 0 then begin
x_bad_here = x_bad[w] - lx & y_bad_here = y_bad[w] - ly
endif
return
end
; STARFINDER_BOXES: extract boxes from image, background and image model.
PRO starfinder_boxes, image, background, stars, x, y, boxsize, $
lx, ux, ly, uy, i_box, b_box, s_box, $
x_bad, y_bad, x_bad_here, y_bad_here
on_error, 2
i_box = sub_array(image, boxsize, REFERENCE = [x, y], $
LX = lx, UX = ux, LY = ly, UY = uy)
s_box = stars[lx:ux,ly:uy] & b_box = background[lx:ux,ly:uy]
if n_elements(x_bad) ne 0 and n_elements(y_bad) ne 0 then $
starfinder_bad, x_bad, y_bad, lx, ux, ly, uy, x_bad_here, y_bad_here
return
end
; STARFINDER_ID: re-identification of a presumed star.
PRO starfinder_id, image, background, stars, sv_par, id_par, $
min_intensity, x, y, found
on_error, 2
if n_elements(sv_par) ne 0 then $
r = pick_region(sv_par.lx, sv_par.ux, sv_par.ly, sv_par.uy, x, y) $
else r = 0
; Extract boxes
starfinder_boxes, image, background, stars, x, y, id_par.box[r], $
lx, ux, ly, uy, i_box, b_box, s_box
; Search
max_search, i_box - b_box - s_box, min_intensity[lx:ux,ly:uy], $
X0 = x - lx, Y0 = y - ly, n, x, y, $
/NEAREST, /MAXIMUM, FOUR = id_par.neigh4
found = n ne 0
if found then begin
x = x[0] + lx & y = y[0] + ly
endif
return
end
; STARFINDER_CORRELATE: correlation check.
PRO starfinder_correlate, image, background, stars, sv_par, x_bad, y_bad, $
corr_par, min_correlation, x, y, correl, accepted
on_error, 2
if n_elements(sv_par) ne 0 then $
r = pick_region(sv_par.lx, sv_par.ux, sv_par.ly, sv_par.uy, x, y) $
else r = 0
; Extract boxes
boxsize = corr_par.correlation_box[r] + corr_par.search_box[r] + 1
starfinder_boxes, image, background, stars, x, y, boxsize, $
lx, ux, ly, uy, i_box, b_box, s_box, x_bad, y_bad, xb, yb
; Compute correlation
if corr_par.correl_mag gt 1 then begin
templates = *corr_par.templates[r] & dx = corr_par.dx & dy = corr_par.dy
endif
correlate_max, i_box - b_box - s_box, *corr_par.template[r], x - lx, y - ly, $
corr_par.search_box[r], X_BAD = xb, Y_BAD = yb, $
TEMPLATES = templates, DX = dx, DY = dy, correl, x, y
x = x + lx & y = y + ly
accepted = correl ge min_correlation and $
x ge lx and x le ux and y ge ly and y le uy
return
end
; STARFINDER_CHECK: check outcome of local fitting.
FUNCTION starfinder_check, fit_error, x_fit, y_fit, f_fit, $
lx, ux, ly, uy, list, min_intensity, fit_par, r
on_error, 2
; Check convergence of fit, minimum acceptable value of fluxes,
; range of positions
; padd = 2.0 ; added after suggestion by Jessica Lu
padd = 0.0
good = fit_error ge 0 and $
min(x_fit) ge lx-padd and max(x_fit) le ux+padd and $
min(y_fit) ge ly-padd and max(y_fit) le uy+padd
if good then begin
threshold = min_intensity[round(x_fit),round(y_fit)] / fit_par.psf_max[r]
good = min(f_fit - threshold) ge 0
good = min(f_fit - threshold) ge 0
endif
if good then begin
; Consider the subset of stars in a neighborhood of [lx:ux,ly:uy]
; and check their reciprocal distances
s = where_stars(LX = lx - fit_par.min_distance[r], $
UX = ux + fit_par.min_distance[r], $
LY = ly - fit_par.min_distance[r], $
UY = uy + fit_par.min_distance[r], list, n)
star_param, list, SUBSCRIPTS = s, n, x, y, f
if n gt 1 then $
good = min(reciprocal_distance(x, y)) ge fit_par.min_distance[r]
endif
return, good
end
; STARFINDER_FIT: local fitting.
PRO starfinder_fit, list, this_max, image, siz, background, stars, noise_std, $
sv_par, x_bad, y_bad, fit_par, fit_data, model_data, $
min_intensity, NO_SLANT = no_slant, $
re_fitting, _EXTRA = extra, star_here, $
PSF_TYPE = psf_type, PSF_DATA = psf_data
on_error, 2
star_param, list, SUBSCRIPTS = this_max, n, x, y
if n_elements(sv_par) ne 0 then $
r = pick_region(sv_par.lx, sv_par.ux, sv_par.ly, sv_par.uy, x, y) $
else r = 0
; Extract boxes
starfinder_boxes, image, background, stars, x, y, fit_par.fitting_box[r], $
lx, ux, ly, uy, i_box, b_box, s_box, x_bad, y_bad, xb, yb
if n_elements(xb) ne 0 and n_elements(yb) ne 0 then $
w_bad = coord_to_subs(xb, yb, (size52(i_box, /DIM))[0])
; Select known stars into i_box and subtract them from the image model
s = where_stars(LX = lx + fit_par.edge, UX = ux - fit_par.edge, $
LY = ly + fit_par.edge, UY = uy - fit_par.edge, list, n)
if re_fitting ne 0 then begin
; if n eq 0 then s = this_max ; edge star
; Modified March 2012
; if n eq 0 then begin
if n eq 0 or re_fitting eq 2 then begin
s = this_max ; edge star or single-component fitting
n = 1
endif
s_and_this = s
endif else s_and_this = add_subscript(s, this_max)
star_param, list, SUBSCRIPTS = s, n, x_add, y_add, f_add, c
if n ne 0 then $
if keyword_set(psf_type) then $
stars = image_model(x_add, y_add, -f_add, siz[0], siz[1], psf_type, $
psf_data, _EXTRA = extra, MODEL = stars) else $
if n_elements(sv_par) ne 0 then $
stars = image_model(x_add, y_add, -f_add, siz[0], siz[1], fit_par.psf, $
LX = sv_par.lx, UX = sv_par.ux, LY = sv_par.ly, UY = sv_par.uy, $
model_data, _EXTRA = extra, MODEL = stars) else $
stars = image_model(x_add, y_add, -f_add, siz[0], siz[1], fit_par.psf, $
model_data, _EXTRA = extra, MODEL = stars)
; Fixed contribution for local fitting
contrib = stars[lx:ux,ly:uy] ; stars outside fitting box
if keyword_set(no_slant) then begin
contrib = contrib + b_box & b_box = 0
endif
; Local fitting
if re_fitting ne 0 then begin
x = x_add & y = y_add & f0 = f_add
endif else $
; Modified March 2012: n is the number of known sources in the box
star_param, list, SUBSCRIPTS = s_and_this, n_and_this, x, y, f, c
; star_param, list, SUBSCRIPTS = s_and_this, n, x, y, f, c
x0 = x - lx & y0 = y - ly
if n_elements(noise_std) ne 0 then noise = noise_std[lx:ux,ly:uy]
if n_elements(sv_par) ne 0 then begin
sv_lx = sv_par.lx - lx
sv_ux = sv_par.ux - lx
sv_ly = sv_par.ly - ly
sv_uy = sv_par.uy - ly
endif
if keyword_set(psf_type) then begin
xori = (*psf_data).xori
yori = (*psf_data).yori
(*psf_data).xori = lx
(*psf_data).yori = ly
fitstars, i_box, FIXED = contrib, psf_type, $
PSF_DATA = psf_data, x0, y0, F0 = f0, BACKGROUND = b_box, $
NO_SLANT = no_slant, NOISE_STD = noise, _EXTRA = extra, $
BAD_DATA = w_bad, x, y, f, b, fit_error, sigma_x, sigma_y, sigma_f
(*psf_data).xori = xori
(*psf_data).yori = yori
endif else $
fitstars, i_box, FIXED = contrib, fit_par.fitting_psf, $
PSF_DATA = fit_data, x0, y0, F0 = f0, BACKGROUND = b_box, $
NO_SLANT = no_slant, NOISE_STD = noise, _EXTRA = extra, $
BAD_DATA = w_bad, LX = sv_lx, UX = sv_ux, LY = sv_ly, UY = sv_uy, $
x, y, f, b, fit_error, sigma_x, sigma_y, sigma_f
; Modified April 2012
; x_out = x & y_out = y
; compare_lists, x0, y0, x_out, y_out, x0_out, y0_out, x, y, SUBSCRIPTS_2 = w
compare_lists, x0, y0, x, y, SUBSCRIPTS_2 = w
x = x[w] + lx & y = y[w] + ly & f = f[w]
if n_elements(sigma_f) ne 0 then begin
sigma_x = sigma_x[w]
sigma_y = sigma_y[w]
sigma_f = sigma_f[w]
endif
; Is the examined max a star? Assume it is, then perform checks
temp_list = list
update_list, temp_list, SUB = s_and_this, x, y, f, c, $
sigma_x, sigma_y, sigma_f, /IS_STAR
star_here = starfinder_check(fit_error, x, y, f, lx, ux, ly, uy, $
temp_list, min_intensity, fit_par, r)
; Update list of stars and image model
if star_here then begin
list = temp_list
x_add = x & y_add = y & f_add = f
endif else if re_fitting eq 2 then begin
list = reverse_class(list, SUBSCRIPT = this_max)
n = n - 1
if n gt 0 then begin
w = where(s_and_this ne this_max[0])
star_param, list, SUBSCRIPTS = s_and_this[w], n, x_add, y_add, f_add, c
endif
endif
if n gt 0 or re_fitting eq 0 then $ ; should be equivalent to "n gt 0 or re_fitting lt 2"
if keyword_set(psf_type) then $
stars = image_model(x_add, y_add, f_add, siz[0], siz[1], psf_type, $
psf_data, _EXTRA = extra, MODEL = stars) else $
if n_elements(sv_par) ne 0 then $
stars = image_model(x_add, y_add, f_add, siz[0], siz[1], fit_par.psf, $
LX = sv_par.lx, UX = sv_par.ux, LY = sv_par.ly, UY = sv_par.uy, $
model_data, _EXTRA = extra, MODEL = stars) else $
stars = image_model(x_add, y_add, f_add, siz[0], siz[1], fit_par.psf, $
model_data, _EXTRA = extra, MODEL = stars)
return
end
; STARFINDER_BLEND: check the area of a presumed star to identify blends.
FUNCTION starfinder_blend, image, background, stars, noise_std, $
deb_par, x, y, r
on_error, 2
starfinder_boxes, image, background, stars, x, y, deb_par.box[r], $
lx, ux, ly, uy, i_box, b_box, s_box
ima = i_box - b_box - s_box ; after this, it should be positive!
x0 = x - lx & y0 = y - ly
threshold = ima[x0,y0] * deb_par.cut
check = 1B
if n_elements(noise_std) ne 0 then $
check = threshold gt deb_par.n_std * noise_std[x,y]
if check then begin
area = peak_area(ima, X = x0, Y = y0, MAG = 3, ABS_THRESHOLD = threshold)
check = relative_error(deb_par.psf_area[r], area) gt deb_par.tol
endif
return, check
end
; STARFINDER_DEBLEND: de-blend a crowded group of stars
; (partial or full de-blending).
PRO starfinder_deblend, list, this_max, image, siz, background, stars, $
noise_std, sv_par, x_bad, y_bad, id_par, fit_par, $
deb_par, fit_data, model_data, min_intensity, $
re_fitting, AROUND = around, _EXTRA = extra, $
PSF_TYPE = psf_type, PSF_DATA = psf_data
on_error, 2
star_param, list, SUB = this_max, n, x, y, f
x0 = round(x) & y0 = round(y)
if n_elements(sv_par) ne 0 then $
r = pick_region(sv_par.lx, sv_par.ux, sv_par.ly, sv_par.uy, x, y) $
else r = 0
; Is the present object a blend?
if keyword_set(around) then $
if keyword_set(psf_type) then $
stars = image_model(x, y, -f, siz[0], siz[1], psf_type, $
psf_data, _EXTRA = extra, MODEL = stars) else $
if n_elements(sv_par) ne 0 then $
stars = image_model(x, y, -f, siz[0], siz[1], fit_par.psf, $
LX = sv_par.lx, UX = sv_par.ux, LY = sv_par.ly, UY = sv_par.uy, $
model_data, _EXTRA = extra, MODEL = stars) else $
stars = image_model(x, y, -f, siz[0], siz[1], fit_par.psf, $
model_data, _EXTRA = extra, MODEL = stars)
; stars = image_model(x, y, -f, siz[0], siz[1], fit_par.psf[*,*,r], $
; model_data, _EXTRA = extra, MODEL = stars)
blended = starfinder_blend(image, background, stars, noise_std, $
deb_par, x0, y0, r)
if keyword_set(around) then $
if keyword_set(psf_type) then $
stars = image_model(x, y, f, siz[0], siz[1], psf_type, $
psf_data, _EXTRA = extra, MODEL = stars) else $
if n_elements(sv_par) ne 0 then $
stars = image_model(x, y, f, siz[0], siz[1], fit_par.psf, $
LX = sv_par.lx, UX = sv_par.ux, LY = sv_par.ly, UY = sv_par.uy, $
model_data, _EXTRA = extra, MODEL = stars) else $
stars = image_model(x, y, f, siz[0], siz[1], fit_par.psf, $
model_data, _EXTRA = extra, MODEL = stars)
; stars = image_model(x, y, f, siz[0], siz[1], fit_par.psf[*,*,r], $
; model_data, _EXTRA = extra, MODEL = stars)
if not blended then return ; not a blend
ncomp = 0L & found = ncomp eq 0 & star_here = found
if keyword_set(around) then ncomp = ncomp + 1
; max_dist = 1.25 * fit_par.psf_fwhm[r]
while star_here do begin
x = x0 & y = y0
if ncomp ne 0 then begin
; Identify another component of the blend
starfinder_boxes, image, background, stars, x, y, deb_par.box[r], $
lx, ux, ly, uy, i_box, b_box, s_box
max_search, i_box - b_box - s_box, $
deb_par.bin_threshold * min_intensity[lx:ux,ly:uy], $
X0 = x - lx, Y0 = y - ly, CHECK_DIST = max_dist, $
/MAXIMUM, FOUR = id_par.neigh4, n, x, y
found = n ne 0 & star_here = found
if found then begin
x = x[0] + lx & y = y[0] + ly
endif
endif
if found then begin
; Fit the new detected component
if ncomp ne 0 then begin
list = merge_list(list, create_element(x, y, 0))
index = n_elements(list) - 1
endif else index = this_max
starfinder_fit, list, index, image, siz, background, stars, $
noise_std, sv_par, x_bad, y_bad, fit_par, fit_data, $
model_data, min_intensity, re_fitting, _EXTRA = extra, $
star_here, PSF_TYPE = psf_type, PSF_DATA = psf_data
endif
if star_here then ncomp = ncomp + 1
endwhile
if ncomp eq 1 and not keyword_set(around) then begin
; The object was mis-recognized as a blend: it is single
star_param, list, SUB = this_max, n, x, y, f
if keyword_set(psf_type) then $
stars = image_model(x, y, -f, siz[0], siz[1], psf_type, $
psf_data, _EXTRA = extra, MODEL = stars) else $
if n_elements(sv_par) ne 0 then $
stars = image_model(x, y, -f, siz[0], siz[1], fit_par.psf, $
LX = sv_par.lx, UX = sv_par.ux, LY = sv_par.ly, UY = sv_par.uy, $
model_data, _EXTRA = extra, MODEL = stars) else $
stars = image_model(x, y, -f, siz[0], siz[1], fit_par.psf, $
model_data, _EXTRA = extra, MODEL = stars)
; stars = image_model(x, y, -f, siz[0], siz[1], fit_par.psf[*,*,r], $
; model_data, _EXTRA = extra)
list = reverse_class(list, SUBSCRIPT = this_max)
endif
return
end
; STARFINDER_ANALYZE: analyze the object list[this_max].
PRO starfinder_analyze, list, this_max, image, siz, background, stars, $
noise_std, sv_par, x_bad, y_bad, id_par, corr_par, $
fit_par, fit_data, model_data, $
min_intensity, min_correlation, re_fitting, _EXTRA = extra
on_error, 2
star_param, list, SUB = this_max, n, x, y
; Is the maximum list[this_max] a feature of an already detected star?
starfinder_id, image, background, stars, sv_par, id_par, $
min_intensity, x, y, check
if not check then return ; yes, it is
; Correlation check
starfinder_correlate, image, background, stars, sv_par, x_bad, y_bad, $
corr_par, min_correlation, x, y, c, check
if not check then return ; too low correlation
; Fit
update_list, list, SUB = this_max, x, y, 0, c
starfinder_fit, list, this_max, image, siz, background, stars, noise_std, $
sv_par, x_bad, y_bad, fit_par, fit_data, model_data, $
min_intensity, re_fitting, _EXTRA = extra
return
end
; STARFINDER_CONVERG: check convergence of positions and fluxes for Phase 3.
FUNCTION starfinder_converg, list0, list, $
ASTROMETRIC_TOL = a_tol, PHOTOMETRIC_TOL = ph_tol
on_error, 2
if n_elements(a_tol) eq 0 then a_tol = 1e-2;1e-4;0.01
if n_elements(ph_tol) eq 0 then ph_tol = 1e-2;1e-4;0.01
s = where_stars(list, n)
if n eq 0 then return, n eq 0
star_param, list0, SUB = s, n, x0, y0, f0
star_param, list, SUB = s, n, x, y, f
check = convergence(x0, x, a_tol, /ABSOLUTE) and $
convergence(y0, y, a_tol, /ABSOLUTE) and $
convergence(f0, f, ph_tol)
return, check
end
;;; The main routine.
PRO starfinder, $
image, psf, SV_PAR = sv_par, $
X_BAD = x_bad, Y_BAD = y_bad, $
BACKGROUND = background, BACK_BOX = back_box, $
threshold, REL_THRESHOLD = rel_threshold, NOISE_STD = noise_std, $
min_correlation, DEBLEND = deblend, DEBLOST = deblost, _EXTRA = extra, $
N_ITER = n_iter, NO_INTERMEDIATE_ITER = no_intermediate, SILENT = silent, $
x, y, fluxes, sigma_x, sigma_y, sigma_f, correlation, STARS = stars, $
X_INPUT = x_in, Y_INPUT = y_in, F_INPUT = f_in, C_INPUT = c_in, $
PSF_TYPE = psf_type, PSF_DATA = psf_data
catch, error
if error ne 0 then begin
starfinder_dealloc, corr_par, fit_data, model_data
return
endif
; Define background and image model
siz = size52(image, /DIM)
if n_elements(background) eq 0 then $
if n_elements(back_box) eq 0 then $
background = fltarr(siz[0], siz[1]) else $
background = estimate_background(image, back_box, _EXTRA = extra)
if n_elements(noise_std) ne 0 then $
if n_elements(noise_std) ne n_elements(image) then $
noise_std = replicate(noise_std[0], siz[0], siz[1])
stars = fltarr(siz[0], siz[1])
; Define some program parameters and default values
re_fitting = 0 ; re_fitting = 0 --> normal mode
if n_elements(x_in) ne 0 and n_elements(y_in) ne 0 and n_elements(f_in) ne 0 then $
re_fitting = 2 ; re_fitting = 2 --> fit of input list of sources
if n_elements(n_iter) eq 0 then n_iter = 1
if size52(psf, /N_DIM) eq 3 then $
n_psf = (size52(psf, /DIM))[2] else n_psf = 1
psf_fwhm = starfinder_fwhm(psf, n_psf)
id_par = starfinder_id_par(psf_fwhm, _EXTRA = extra)
corr_par = starfinder_corr_par(psf, psf_fwhm, n_psf, _EXTRA = extra)
; modified March 2012
if re_fitting eq 2 then $
fit_par = starfinder_fit_par(psf, psf_fwhm, n_psf, /WIDER, _EXTRA = extra) else $
fit_par = starfinder_fit_par(psf, psf_fwhm, n_psf, _EXTRA = extra)
; fit_par = starfinder_fit_par(psf, psf_fwhm, n_psf, _EXTRA = extra)
if keyword_set(deblend) or keyword_set(deblost) then $
deb_par = starfinder_deb_par(psf, psf_fwhm, n_psf, _EXTRA = extra)
if n_elements(sv_par) eq 0 then begin
fit_data = ptr_new(/ALLOCATE) & model_data = ptr_new(/ALLOCATE)
endif
; Iterative search for suspected stars and subsequent analysis
n_levels = n_elements(threshold)
n_stars = 0L & n_presumed = 0L & n_max = 0L
if re_fitting eq 2 then begin
n_levels = 1
n_stars = n_elements(f_in)
list_of_stars = sort_list(reverse_class(starlist(n_stars, x_in, y_in, f_in)))
if keyword_set(psf_type) then $
stars = image_model(x_in, y_in, f_in, siz[0], siz[1], psf_type, $
psf_data, _EXTRA = extra) else $
if n_elements(sv_par) ne 0 then $
stars = image_model(x_in, y_in, f_in, siz[0], siz[1], psf, $
LX = sv_par.lx, UX = sv_par.ux, LY = sv_par.ly, UY = sv_par.uy, $
_EXTRA = extra) else $
stars = image_model(x_in, y_in, f_in, siz[0], siz[1], psf, _EXTRA = extra)
endif
for n_lev = 0L, n_levels - 1 do begin
; Update detection threshold
threshold_n = float(threshold[n_lev])
if keyword_set(rel_threshold) and n_elements(noise_std) ne 0 then $
threshold_n = threshold_n * noise_std
if n_elements(threshold_n) eq 1 then $
threshold_n = replicate(threshold_n[0], siz[0], siz[1])
; Phase 1: find local maxima having
; central intensity > known stars + background + threshold
if re_fitting eq 0 then begin
if not keyword_set(silent) then $
print, "STARFINDER: search for candidate stars"
search_objects, image - stars, LOW_SURFACE = background, threshold_n, $
_EXTRA = extra, n_max, x0, y0, i0
n_presumed = n_presumed + n_max
endif
; Phase 2: analyze presumed stars with correlation and fitting,
; estimating positions and fluxes.
if n_max ne 0 and re_fitting eq 0 then begin
if not keyword_set(silent) then $
print, "STARFINDER: analysis of candidate stars"
list_of_max = starlist(n_max, x0, y0, i0)
list_of_stars = merge_list(list_of_stars, sort_list(list_of_max))
for n = n_stars, n_stars - 1 + n_max do $
starfinder_analyze, list_of_stars, n, image, siz, background, stars, $
noise_std, sv_par, x_bad, y_bad, id_par, $
corr_par, fit_par, fit_data, model_data, $
threshold_n, min_correlation, re_fitting, _EXTRA = extra, $
PSF_TYPE = psf_type, PSF_DATA = psf_data
list_of_stars = sort_list(extract_stars(list_of_stars, n_stars))
; Update background estimate
if n_elements(back_box) ne 0 then $
background = estimate_background(image - stars, back_box, _EXTRA = extra)
; Modify fitting parameters for subsequent operations
fit_par = starfinder_fit_par(psf, psf_fwhm, n_psf, /WIDER, _EXTRA = extra)
if ptr_valid(fit_data) then ptr_free, fit_data
endif ; n_max ne 0
; Phase 2b: de-blend detected objects.
if keyword_set(deblend) and n_lev eq n_levels - 1 and n_stars ne 0 and re_fitting eq 0 then begin
if not keyword_set(silent) then $
print, "STARFINDER: deblending of detected stars"
; Check each detected star: is it a blend?
for n = 0L, n_stars - 1 do $
starfinder_deblend, list_of_stars, n, image, siz, background, stars, $
noise_std, sv_par, x_bad, y_bad, id_par, $
fit_par, deb_par, fit_data, model_data, $
threshold_n, re_fitting, _EXTRA = extra, /AROUND, $
PSF_TYPE = psf_type, PSF_DATA = psf_data
list_of_stars = sort_list(extract_stars(list_of_stars, n_stars))
; Update background estimate
if n_elements(back_box) ne 0 then $
background = estimate_background(image - stars, back_box, _EXTRA = extra)
endif
; Phase 2c: recover lost objects to search for blends.
; if keyword_set(deblost) and n_lev eq n_levels - 1 and n_stars ne 0 and re_fitting eq 0 then begin
if keyword_set(deblost) and n_lev eq n_levels - 1 and re_fitting eq 0 then begin
if not keyword_set(silent) then $
print, "STARFINDER: search for blends among lost objects"
search_objects, image - stars, LOW_SURFACE = background, threshold_n, $
_EXTRA = extra, n_max, x0, y0, i0
if n_max ne 0 then begin
list_of_max = starlist(n_max, x0, y0, i0)
list_of_stars = merge_list(list_of_stars, sort_list(list_of_max))
for n = n_stars, n_stars - 1 + n_max do $
starfinder_deblend, list_of_stars, n, image, siz, background, stars, $
noise_std, sv_par, x_bad, y_bad, id_par, fit_par, $
deb_par, fit_data, model_data, threshold_n, $
re_fitting, _EXTRA = extra, $
PSF_TYPE = psf_type, PSF_DATA = psf_data
list_of_stars = sort_list(extract_stars(list_of_stars, n_stars))
; Update background estimate
if n_elements(back_box) ne 0 then $
background = estimate_background(image - stars, back_box, _EXTRA = extra)
endif
endif
; Phase 3: re-determination of positions and fluxes by iterative
; fitting of detected stars.
if n_lev eq n_levels - 1 then maxit = n_iter else $
if keyword_set(no_intermediate) then maxit = 0 else maxit = 1
if n_stars ne 0 then list0 = list_of_stars
if re_fitting eq 0 then re_fitting = 1
iter = 0L & converging = n_stars eq 0
while iter lt maxit and not converging do begin
if not keyword_set(silent) then $
if n_lev lt n_levels - 1 then $
print, "STARFINDER: re-fitting" else $
print, "STARFINDER: re-fitting: iteration", iter + 1
; Modified March 2012
; if n_lev lt n_levels - 1 then $
; print, "STARFINDER: intermediate re-fitting" else $
; print, "STARFINDER: final re-fitting: iteration", iter + 1
for n = 0L, n_stars - 1 do $
; modified March 2012
; if list_of_stars[n].is_a_star or re_fitting eq 1 then $ ;;; if list_of_stars[n].is_a_star or re_fitting eq 1 or iter eq 0 then $
; modified September 2024
if is_star(list_of_stars, n) or re_fitting eq 1 then $ ;;; if list_of_stars[n].is_a_star or re_fitting eq 1 or iter eq 0 then $
starfinder_fit, list_of_stars, n, image, siz, background, stars, $
noise_std, sv_par, x_bad, y_bad, fit_par, fit_data, $
model_data, threshold_n, re_fitting, _EXTRA = extra, $
PSF_TYPE = psf_type, PSF_DATA = psf_data
if maxit gt 1 then $
converging = starfinder_converg(list0, list_of_stars, _EXTRA = extra)
iter = iter + 1 & list0 = list_of_stars
endwhile
if re_fitting eq 1 then re_fitting = 0
fit_par = starfinder_fit_par(psf, psf_fwhm, n_psf, _EXTRA = extra)
if ptr_valid(fit_data) then ptr_free, fit_data
endfor ; end of cycle on threshold levels
; De-allocate pointer heap variables
starfinder_dealloc, corr_par, fit_data, model_data
; Save results
if n_stars ne 0 then begin
list_of_stars = sort_list(extract_stars(list_of_stars, n_stars))
star_param, list_of_stars, n_stars, $
x, y, fluxes, correlation, sigma_x, sigma_y, sigma_f
endif
if not keyword_set(silent) then begin
if re_fitting ne 2 then $
print, "STARFINDER", n_presumed, " candidate stars"
print, "STARFINDER", n_stars, " detected stars"
endif
return
end