;+ ; NAME: ; psffit ; PURPOSE: (one line only) ; Fit a numerical PSF to one or more sources in an image. ; DESCRIPTION: ; The PSF must be provided as an input to this routine. Suitable ; PSFs can be generated with PSFSTACK. For this routine to work, the ; resulting PSF must be amenable to shifting by interpolation. Thus ; an under-sampled PSF will not work with this tool. ; Generally speaking, if you ask to fit more than one source at a time, ; the sources should have overlapping PSFs. Otherwise, you are better ; off with separate calls. This is NOT enforced. ; CATEGORY: ; CCD Data Processing ; CALLING SEQUENCE: ; psffit,image,psf,region,xin,yin,xout,yout,flux ; INPUTS: ; image - Input image to be fitted. ; psf - Numerical PSF to fit to image. This must be normalized to have ; unit volume. ; region - Fitting region for chi-sq calculation. Provide ; [x1,x2,y1,y2] in image index coordinates. ; xin - Starting x location for source(s) ; yin - Starting y location for source(s) ; OPTIONAL INPUT PARAMETERS: ; KEYWORD INPUT PARAMETERS: ; SILENT - Flag, if set will suppress all diagnostic output ; PSFMAX - Two element vector containing the x,y position of the peak of ; the input PSF. The position is in the native array coordinates ; of the psf array. If not provided, this will be computed. ; It will save time for many fits with the same psf to provide ; this value. ; SKY - Background sky value (from entire image). If this is provided ; and calcsky is not set then this value is used for sky. If ; calcsky is set then sky is computed and returned in this variable. ; CALCSKY - Flag, if set forces this routine to compute sky background. ; GAIN - Gain of the image, e-/ADU. This is used to generate an array ; of uncertainties for each pixel based on photon statistics. ; default=1.0 ; RDNOISE - Readout noise of CCD (default=1) [e-] ; NMAX - Maximum number of iterations for the amoeba call. Default=5000 ; If this limit is reached, the source will be tagged as ; not fittable, ie., chisq = -1. ; OUTPUTS: ; xout - Fitted x location for source(s) ; yout - Fitted y location for source(s) ; flux - Integrated flux for source(s) in DN ; KEYWORD OUTPUT PARAMETERS: ; CHISQ - final goodness of fit for the returned values ; COMMON BLOCKS: ; PSFFIT_COM, for internal use only ; SIDE EFFECTS: ; RESTRICTIONS: ; PROCEDURE: ; MODIFICATION HISTORY: ; Written by Marc W. Buie, Southwest Research Institute, 2012/09/28 ; 2012/10/04, MWB, added weighted fitting. ; 2012/12/06, MWB, fixed a problem with the SKY and CALCSKY keywords ; 2013/03/05, MWB, added NMAX keyword and proper trap of bad fit from amoeba ;- function psffit_func,vals common psffit_com,info nvals=n_elements(vals) nsrc = nvals/3 x=vals[0:nsrc-1] y=vals[nsrc:2*nsrc-1] f=abs(vals[2*nsrc:3*nsrc-1]) model=fltarr(info.nx,info.ny) addpsf,model,x,y,f,info.psf,psfmax=[info.xm,info.ym] model += info.sky dimage=info.image[info.x0:info.x1,info.y0:info.y1]- $ model[info.x0:info.x1,info.y0:info.y1] ; chisq = total(dimage^2) / info.npix sigarr=info.simage[info.x0:info.x1,info.y0:info.y1] chisq = total((dimage/sigarr)^2) / info.npix if chisq lt info.chibest then begin if not info.silent then print,info.chicount,chisq,x[0],y[0],f[0] info.chibest = chisq endif info.chicount++ return,chisq end pro psffit,image,psf,region,xin,yin,xout,yout,flux, $ SILENT=silent,PSFMAX=psfmax,SKY=sky,CALCSKY=calcsky,CHISQ=chisq, $ GAIN=gain,RDNOISE=rdnoise,NMAX=nmax common psffit_com,info self='psffit: ' if badpar(image,[2,3,4,5],2,caller=self+'(image) ') then return if badpar(psf,[2,3,4,5],2,caller=self+'(psf) ') then return if badpar(region,[2,3],1,caller=self+'(region) ') then return if badpar(xin,[2,3,4,5],1,caller=self+'(xin) ',npts=npts) then return if badpar(yin,[2,3,4,5],1,caller=self+'(yin) ') then return if badpar(silent,[0,1,2,3],0,caller=self+'(SILENT) ',default=0) then return if badpar(psfmax,[0,2,3],1,caller=self+'(PSFMAX) ',default=[-1,-1]) then return if badpar(calcsky,[0,2,3,4,5],0,caller=self+'(CALCSKY) ',default=0) then return if badpar(gain,[0,2,3,4,5],0,caller=self+'(GAIN) ',default=1.) then return if badpar(rdnoise,[0,2,3,4,5],0,caller=self+'(RDNOISE) ',default=1.) then return if badpar(nmax,[0,2,3],0,caller=self+'(NMAX) ',default=5000) then return sz=size(image,/dimen) if min(psfmax) lt 0 then begin psfdim=size(psf,/dimen) boxm,psf,psfdim[0]/2,psfdim[1]/2,psfdim[0]/2-1,psfdim[1]/2-1,xm,ym findmax,xm,ym,psf,xmax,ymax,fm endif else begin xmax=psfmax[0] ymax=psfmax[1] endelse if calcsky then begin skysclim,image,lowval,hival,meanval,sigma, $ lowclip=0.1,hiclip=0.9,npts=60000 sky=meanval endif else begin if badpar(sky,[2,3,4,5],0,caller=self+'(SKY) ') then return endelse simage=sqrt(image*gain+rdnoise^2)/gain npix = (region[1]-region[0]+1) * (region[3]-region[2]+1) if region[0] ge 0 and region[0] lt sz[0] and $ region[1] ge 0 and region[1] lt sz[0] and $ region[2] ge 0 and region[2] lt sz[1] and $ region[3] ge 0 and region[3] lt sz[1] then begin info= { image: image, $ simage: simage, $ nx: sz[0], $ ny: sz[1], $ psf: psf, $ xm: xmax, $ ym: ymax, $ x0: region[0], $ x1: region[1], $ y0: region[2], $ y1: region[3], $ sky: sky, $ silent: silent, $ chicount: 0L, $ chibest: 1.0e29, $ npix: npix $ } ftol=1.0e-4 flux0=image[round(xin),round(yin)]/max(psf)*0.8 start=[xin,yin,flux0] if not silent then print,'start ',start scale=[replicate(0.1,npts*2),flux0*0.01] if not silent then print,'scale ',scale fitval=amoeba(ftol,FUNCTION_NAME='psffit_func', $ P0=start, SCALE=scale, NMAX=nmax) if n_elements(fitval) gt 1 then begin chisq=psffit_func(fitval) xout=trimrank(fitval[0:npts-1]) yout=trimrank(fitval[npts:2*npts-1]) flux=trimrank(fitval[2*npts:3*npts-1]) endif else begin xout=xin yout=yin flux=fltarr(npts) chisq=-1.0 endelse endif else begin xout=xin yout=yin flux=fltarr(npts) chisq=-1.0 endelse end