; docformat = 'rst'
;
; NAME: 
;    selectAperture2D
; PURPOSE:
;    Select aperture from multiAperture2D for source characterization.
;+
; :Description:
;    Select aperture from multiAperture2D for source characterization.
;
; :Categories:
;    Images
;
; :Returns: The indices of selected aperture, according to input parameter statistics and threshold.
;
; :Params:
;    param: in, required, type=structure
;      The peak parameters, from multiAperture2D. Fields are:
;      .id:   Parameter simple identifyer
;      .name: Name of the parameter
;      .unit: Unit of the parameter
;      .val:  An array of N values, one per aperture
;      .mean: Average value for N apertures
;      .dev:  Standard deviation for N apertures
;    xy: in, optional, type=fltarr(2)
;      Estimated centroid position used to check against source confusion.
;
; :Keywords:
;    cutXY: in, optional, type=float, default=0.0025
;      Threshold for centroid X/Y variation selection.
;    cutEE: in, optional, type=float, default=0.035
;      Threshold for encircled energy variation selection.
;    sepXY: in, optional, type=float, default=5
;      Threshold for centroid X/Y separation from guessed position.
;
; :Examples:
;
; :Uses:
;
; :Author:
;    B.Carry (IMCCE)
;
; :History:
;     Change History::
;        Written in April 2015, B. Carry (IMCCE)
;-
function selectAperture2D, param, xy, cutXY=cutXY, cutEE=cutEE, sepXY=sepXY
  COMPILE_OPT hidden
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--I.1-- Check Input Parameters  -------------------------------------------------------
  if not keyword_set(param) then begin
    message, 'Syntax: index = selectAperture2d( param)'
    return, -1
  endif

 ;--I.2-- Peak Model: Gauss or Moffat  --------------------------------------------------
  nbParam = n_elements( param )
  nbMPF = nbParam - 3

 ;--I.3-- Parameter Indices  ------------------------------------------------------------
  indRad = 0
  indPar = 1+indgen(nbMPF)
  indEE  = 1+nbMPF
  indRMS = 1+nbMPF+1

 ;--I.4-- Selection Thresholds  ---------------------------------------------------------
  if not keyword_set(cutXY) then cutXY = 0.0025
  if not keyword_set(cutEE) then cutEE = 0.035
  if not keyword_set(sepXY) then sepXY = 5


;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Select Appropriate Radius Indices                  -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;--II.1-- Variation Functions  ----------------------------------------------------------
  varXc = shift(param[indPar[4]].val,1)/param[indPar[4]].val - 1.
  varYc = shift(param[indPar[5]].val,1)/param[indPar[5]].val - 1.
  varEE = shift(param[  indEE  ].val,1)/param[  indEE  ].val - 1.

 ;--II.2-- Select  -----------------------------------------------------------------------
 ;--II.2.1-- Original Centroid Provided
  if keyword_set(xy) then begin
    sel = where( param[indPar[4]].val-xy(0) le sepXY and $
                 param[indPar[5]].val-xy(1) le sepXY and $
                 abs(varXc) le cutXY and $
                 abs(varYc) le cutXY and $
                 abs(varEE) le cutEE, nbSel )

 ;--II.2.2-- Original Centroid Unknown
  endif else begin
    sel = where( abs(varXc) le cutXY and $
                 abs(varYc) le cutXY and $
                 abs(varEE) le cutEE, nbSel )
  endelse
  if nbSel eq 0 then return, -1

 ;--II.3-- Parameter Statistics  ---------------------------------------------------------
  for kP=0, nbMPF-1 do begin
;    param(indPar(kP)).box = createBoxPlotData( param(indPar(kP)).val(sel) )
    param[indPar[kP]].mean=   mean(param[indPar[kP]].val[sel]) 
    param[indPar[kP]].dev = stddev(param[indPar[kP]].val[sel])
  endfor



 ;--II.4-- Sigma-Clipping  ---------------------------------------------------------------
  cutSig=3.
  in = where( abs(param[indPar[4]].val-param[indPar[4]].box[2]) le param[indPar[4]].dev*cutSig and $
              abs(param[indPar[5]].val-param[indPar[5]].box[2]) le param[indPar[5]].dev*cutSig and $
              abs(param[  indEE  ].val-param[  indEE  ].box[2]) le param[  indEE  ].dev*cutSig, nbIn )
  if nbIn ne nbSel then begin
    sel=in
    for kP=0, nbMPF-1 do begin
      param[indPar[kP]].box = createBoxPlotData( param[indPar[kP]].val[sel] )
      param[indPar[kP]].mean=   mean(param[indPar[kP]].val[sel]) 
      param[indPar[kP]].dev = stddev(param[indPar[kP]].val[sel])
    endfor
    
  endif




 ;--III-- Return Selection  -------------------------------------------------------------
  return, sel

end