Viewing contents of file '../idllib/astron/contrib/varosi/vlibm/allpro/deconv_get_psf.pro'
;+
; NAME:
; deconv_get_psf
; PURPOSE:
; Function extracts and models a PSF from the input image.
; CALLING EXAMPLE:
; psf = deconv_get_psf( image, psf_stats, /MENUS, /TV )
; INPUTS:
; image
; KEYWORDS:
; /MENUS : allow interactive processing options via menu queries,
; then other keywords need not be set.
; /SELECT_PSF : interactively select PSF data from supplied image,
; otherwise whole image is used as PSF.
; /MODEL_PSF : replace data with model (user selected) fit PSF function.
; /EXTRAPOLATE_PSF : replace wings of PSF (below noise level) by model.
; /TV_MONITOR : display images in windows.
; SIZE_PSF = optional 2 element array specifying the size to extract.
; NAME = optional string giving a name to the PSF
; NOISE_MODEL = optional string: "GAUSSIAN", "POISSON", "NONE"
; OUTPUTS:
; psf_stats = structure containing info/statistics about PSF.
;
; Function returns the PSF, extracted and modeled from the input image.
;
; EXTERNAL CALLS:
; function extract_region
; function im_stats
; function FullWid_HalfMax
; function psf_gaussian
; function psf_Lorentzian
; function psf_merge
; pro centroid
; pro show_psf
; PROCEDURE:
; Complicated, might as well look at the code.
; HISTORY:
; Written, Frank Varosi NASA/GSFC 1992.
;-
function deconv_get_psf, image, psf_stats, MENUS=menus, TV_MONITOR=tv_mon, $
SIZE_PSF=sxy, NAME=name_psf, NOISE_MODEL=noise_model, $
EXTRAPOLATE_PSF=extrap, MODEL_PSF=model, SELECT_PSF=select
START:
if keyword_set( menus ) then begin
menu = ["PSF processing options:" ,$
" " ,$
"SELECT subregion of PSF image" ,$
" " ,$
"MODEL by fitting profiles" ,$
" " ,$
"EXTRAPOLATE by fitting profiles" ,$
" " ,$
"continue processing PSF" ]
opt = ""
select=0
model=0
extrap=0
while (opt NE "continue") do begin
msel = wmenu( menu, INIT=8, TIT=0 ) > 0
request = menu(msel)
opt = next_word( request )
CASE opt OF
"SELECT": BEGIN
menu(msel) ="cancel subregion selection"
select=1
print," will select subregion of PSF"
END
"MODEL": BEGIN
menu(msel) = "cancel modeling of PSF"
model=1
print," PSF will be modeled"
END
"EXTRAPOLATE": BEGIN
menu(msel)="cancel extrapolation of PSF"
extrap=1
print," PSF will be modeled"
END
"no": BEGIN
CASE next_word( request ) OF
"selection": BEGIN
menu(msel) = "SELECT subregion"
select=0
END
"modeling": BEGIN
menu(msel) = "MODEL by fitting PSF"
model=0
END
"extrapolation": BEGIN
menu(msel) ="EXTRAPOLATE by fitting PSF"
extrap=0
END
else:
ENDCASE
END
else:
ENDCASE
endwhile
endif
SELECT:
if keyword_set( model ) OR $
keyword_set( extrap ) OR $
keyword_set( select ) then begin
minLog = 10.^( nint( aLog10( max( image ) ) ) - 4 )
psf = extract_region( image,/DISPLAY_IMAGE, LOG=minLog ) > 0
spsf = size( psf )
sxy = spsf(1:2)
if min( sxy ) LT 7 then begin
print," "
message,"selected region is too small...",/INFO
goto,SELECT
endif
endif else psf = image > 0
psf_stats = im_stats( psf, /GET_FWHM, /PRINT, $
NAME=name_psf, NOISE=noise_model )
centroid, psf, cx,cy, XY_P=psf_stats.peak_xy
print," centroid at (x,y) :",cx,cy," pixels"
psf_stats.cent_xy = [cx,cy]
cxy = round_off( [cx,cy]*2 )/2.
print," using PSF center pixel = ",cxy
cxy_frac = abs( cxy - fix( cxy ) )
if N_elements( sxy ) NE 2 then begin
spsf = size( psf )
print," current PSF size (default): ",spsf(1:2)
input=""
read," enter desired x & y size of PSF : ",input
if strlen( input ) LE 0 then sxy = spsf(1:2) else begin
sxy = fix( input )
if N_elements( sxy ) EQ 1 then sxy = [sxy,sxy]
endelse
; Center the PSF in the array,
; so the size is an odd # of pixels if centroid is in middle of a pixel,
; or the size is an even # of pixels if centroid is between pixels:
sxy = sxy < ( 2*( cxy < (spsf(1:2)-cxy) ) )
we = where( cxy_frac EQ 0, neven )
if (neven GT 0) then sxy(we) = sxy(we) - (sxy(we) MOD 2)
wo = where( cxy_frac, nodd )
if (nodd GT 0) then sxy(wo) = sxy(wo) - ((sxy(wo)+1) MOD 2)
sxy2 = sxy/2.
b = cxy - sxy2
L = cxy + sxy2 -1
psf = psf(b(0):L(0),b(1):L(1))
endif else begin
; Center the PSF in the array, but do not decrease array size:
sxy = 2*( cxy > (sxy - cxy) )
b = ( sxy/2 - cxy ) > 0
psfn = fltarr( sxy(0), sxy(1) )
psfn(b(0),b(1)) = psf
psf = psfn
show_psf, psf, NAME=psf_stats.name,/PRINT_INFO
if yes_no_menu("reselect PSF region",/BIN,/NO) then goto,SELECT
endelse
if keyword_set( model ) OR keyword_set( extrap ) then begin
if min( psf_stats.fwhm_xy ) LT 1 then begin
message,"cannot determine FWHM for modeling/extrap",/INF
goto,START
endif
if keyword_set( extrap ) then $
menu_tit=" PSF model for extrapolation ?" $
else menu_tit=" PSF model ?"
MODEL: if keyword_set( tv_mon ) OR keyword_set( menus ) then begin
menu = [menu_tit ,$
" " ,$
"L = generalized Lorentzian" ,$
" " ,$
"G = Gaussian" ]
psf_type = $
strmid( menu( wmenu( menu,INIT=2,TIT=0 )>0 ), 0,1 )
endif else begin
psf_type = ""
print, menu_tit
print," G = Gaussian"
print," L = generalized Lorentzian (default)"
read, psf_type
endelse
minfit = psf_stats.sky_Level + psf_stats.sigma_noise
if strupcase( psf_type ) EQ "G" then begin
FWHM = FullWid_HalfMax( psf,/GAUSSIAN, FIT_PAR=parms, $
MINV=minfit )
mpix = ( fix( 5*FWHM - sxy )/2 + 1 ) > 0
npix = 2 * mpix + sxy
print," model PSF size = 5 * FWHM = ",npix
parms(*,0) = psf_stats.max
parms(0,1) = parms(*,1) + mpix
psf_mod = psf_gaussian( parms, NPIX=npix ) > 0
parms(0,2) = parms(*,2) * 2 * sqrt( 2* aLog(2) )
psf_type = "Gaussian"
get=2
endif else begin
FWHM = FullWid_HalfMax( psf, /LORENTZ, FIT_PAR=parms, $
MINV=minfit )
mpix = ( fix( 13*FWHM - sxy )/2 + 1 ) > 0
npix = 2 * mpix + sxy
print," model PSF size = 13 * FWHM = ",npix
parms(*,0) = psf_stats.max
parms(0,1) = parms(*,1) + mpix
psf_mod = psf_Lorentzian( parms, NPIX=npix ) > 0
parms(0,2) = parms(*,2) * 2
psf_type = "Lorentz"
get=3
endelse
if keyword_set( extrap ) then begin
psf_mod = psf_merge( psf, psf_mod, LEVEL_SPLICE=minfit )
psf_name =psf_stats.name+" ("+psf_type+" extrapolation)"
endif else psf_name = psf_stats.name+" ("+psf_type+" model)"
sxym = size( psf_mod )
b = [ mpix, mpix ]
sxym = size( psf_mod )
L = sxym(1:2) -mpix -1
psf_diff = psf - psf_mod(b(0):L(0),b(1):L(1))
print," model PSF error (scale 0 to 1) = ", $
sqrt( total( psf_diff^2 )/total( psf^2 ) )
if keyword_set( tv_mon ) OR keyword_set( menus ) then begin
show_psf, psf_mod, NAME=psf_name, /PRINT_INFO
if yes_no_menu( "reselect PSF model",/BIN,/NO ) then $
goto,MODEL
endif else begin
print," " + psf_type + " parms:"
print," FWHM x y"
print, parms(*,2:*)
text=""
read," PSF model okay [Y]/N ? ",text
if strlen( text ) GT 0 then begin
if strupcase( strmid( text,0,1 ) ) EQ "N" then $
goto,MODEL
endif
endelse
psf = psf_mod
psf_stats = im_stats( psf, GET_FWHM=get, NOISE="none", $
NAME=psf_name )
endif
return, psf/total( psf )
end