; docformat = 'rst' ; ; NAME: ; multiAperture2D ; PURPOSE: ; Fit a Gaussian or Moffat 2D profile to a source, for different aperture sizes. ;+ ; :Description: ; Fit a Gaussian or Moffat 2D profile to a source, for different aperture sizes. ; ; :Categories: ; Images ; ; :Returns: A structure with all the fitted parameters for each ; aperture. 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 ; ; :Params: ; xy: in, required, type=fltarr(2) ; Coordinates of the source to analyse. ; ; :Keywords: ; rMin: in, optional, type=float, default=3 ; Minimum radius for aperture. ; rMax: in, optional, type=float, default=40 ; Maximum radius for aperture ; gaussian: in, optional, type=boolean, default=1 ; Use Gaussian 2D peak ; moffat: in, optional, type=boolean, default=0 ; Use Moffat 2D peak ; ; :Examples: ; ; :Uses: ; ; :Author: ; B.Carry (IMCCE) ; ; :History: ; Change History:: ; Written in April 2015, B. Carry (IMCCE) ; 2015-Apr. : B.Carry (IMCCE) - Keyword sepXY added ;- function multiAperture2D, im, xy, rMin=rMin, rMax=rMax, gaussian=gaussian, moffat=moffat, sepXY=sepXY ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; COMPILE_OPT hidden ;--I.1-- Check Input Parameters ------------------------------------------------------- if n_params() lt 2 then begin message, 'Syntax: params = multiAperture2d( image, xy [,min=, max=, /gaussian, /moffat])' return, -1 endif ;--I.2-- Image Dimension and Search Radii --------------------------------------------- dim = size(im) if not keyword_set(rMin) then rMin = 3 if not keyword_set(rMax) then rMax = min([ xy(0), xy(1), dim(1)-xy(0), dim(2)-xy(1) ]) $ else rMax = min([ xy(0), xy(1), dim(1)-xy(0), dim(2)-xy(1), rMax ]) nbR = rMax-rMin rad = rMin + indgen(nbR) ;--I.3-- Peak Model: Gauss or Moffat -------------------------------------------------- ;--I.3.1-- At least One, but Not Both if not keyword_set(gaussian) and not keyword_set(moffat) then gaussian=1 if keyword_set(gaussian) and keyword_set(moffat) then begin message, 'Cannot fit Gaussian AND Moffat Peak, please choose' endif ;--I.3.2-- Number of Parameters if keyword_set(gaussian) then begin gaussian = 1 & moffat = 0 nbMPF = 7 endif else begin gaussian = 0 & moffat = 1 nbMPF = 8 endelse ;--I.4-- Parameter Indices ------------------------------------------------------------ nbParam = nbMPF + 3 indRad = 0 indPar = 1+indgen(nbMPF) indEE = 1+nbMPF indRMS = 1+nbMPF+1 ;--I.4-- Output Structure ------------------------------------------------------------- ;--I.4.1-- Structure Definition eParam = {id: '', name:'', unit:'', format: '(E8.2)', $ sym: 'filled circle', color:'black', $ val:fltarr(nbR), mean:0., dev:0., box:fltarr(5)} param = replicate( eParam, nbParam ) ;--I.4.2-- Parameter Labels ;--I.4.2/A-- In-house Parameters param(indRad).id = 'radius' param(indRad).name = 'Aperture radius' param(indRad).unit = 'pix' param(indRad).format= '(I)' param(indEE).id = 'EE' param(indEE).name = 'Encircled Energy' param(indEE).unit = 'ADU' param(indRMS).id = 'RMS' param(indRMS).name = 'RMS residuals' param(indRMS).unit = 'ADU' ;--I.4.2/B-- Inherited from MPFIT2DPEAK mpFitID = ['base','peak','sigx','sigy','xc','yc','tilt','index'] mpFitName = ['Constant baseline level', 'Peak value', $ 'Peak half-width (x)', 'Peak half-width (y)', $ 'Peak centroid (x)', 'Peak centroid (y)', $ 'Rotation angle', 'Moffat power law index'] mpFitUnit = [ 'ADU','ADU','pix','pix','pix','pix','radian',''] mpFitFormat= ['(E8.2)','(E8.2)','(F5.2)','(F5.2)','(F7.3)','(F7.3)','(F6.1)','(E8.2)'] param(indPar).id = mpFitID( 0:nbMPF-1) param(indPar).name = mpFitName( 0:nbMPF-1) param(indPar).unit = mpFitUnit( 0:nbMPF-1) param(indPar).format= mpFitFormat(0:nbMPF-1) param(indPar([2,4])).color = 'Red' param(indPar([3,5])).color = 'Cornflower blue' ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Loop over Aperture and Fit 2D Peak -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Aperture Loop ---------------------------------------------------------------- for kR=0, nbR-1 do begin ;--II.2-- Crop Image to Aperture ------------------------------------------------------- zoom = cropFrame( im, xy, 2*rad(kR) ) nbZoom = n_elements(zoom) param(indRad).val(kR) = rad(kR) ;--II.3-- Fit 2-D Peak ----------------------------------------------------------------- model2d = mpfit2dpeak( zoom, vals, chisq=rms, gaussian=gaussian, moffat=moffat, /tilt ) ;--II.4-- Store Parameters in Output Structure ----------------------------------------- ;--II.4.1-- Blind Storage param(indPar).val(kR) = vals ;--II.4.2-- Centroid Offset to Account for Zoom param(indPar(4:5)).val(kR) += 1.0*(xy - rad(kR)) ;--II.4.3-- Smooth Rotation if kR ge 1 then begin offset = abs(param(indPar(6)).val(kR) - param(indPar(6)).val(kR-1)) if abs( offset-!PI/2 ) le 0.01 then begin swap = param(indPar(2)).val(kr) param(indPar(2)).val(kr) = param(indPar(3)).val(kr) param(indPar(3)).val(kr) = swap param(indPar(6)).val(kr) -= !PI/2 endif endif ;--II.5-- Encircled Energy ------------------------------------------------------------- param(indEE).val(kR) = total( zoom-vals(0) ) ;--II.6-- RMS Residuals ---------------------------------------------------------------- param(indRMS).val(kR) = sqrt( rms/nbZoom ) endfor ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Quick Statistics on Parameters -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- MPFit Parameters ------------------------------------------------------------- for kP=0, nbMPF-1 do begin param(indPar(kP)).box = createBoxPlotData( param(indPar(kP)).val ) param(indPar(kP)).mean= mean(param(indPar(kP)).val) param(indPar(kP)).dev = stddev(param(indPar(kP)).val) endfor ;--III.2-- Encircle Energy -------------------------------------------------------------- param(indEE).box = createboxplotdata( param(indEE).val ) param(indEE).mean= mean(param(indEE).val) param(indEE).dev = stddev(param(indEE).val) ;--III.3-- Return Parameter Structure --------------------------------------------------- return, param end