Viewing contents of file '../idllib/astron/contrib/varosi/vlibm/allpro/deconv_analyze.pro'
;+
; NAME:
; DeConv_Analyze
; PURPOSE:
; Analyze current deconvolution result by
; computing residuals and its auto-correlation,
; then display result with residuals, and print statistics.
; Called by DeConv_Tool and DeConv_Review.
; CALLING:
; deconv_analyze, image_data, psf_obs, image_deconv, deconv_info, /INIT
; INPUTS:
; image_data = original observed image data.
; psf_obs = point spread function of observation.
; image_deconv = current deconvolution result.
; deconv_info = structure containing image/PSF info and deconv statistics.
; KEYWORDS:
; /INITIALIZE : must be set first call, creates windows and prints Labels.
; /NO_RESIDUALS : skips the computation & display of residual image.
; RE_CONVOL_IMAGE = optional input of image_deconv convolved with PSF,
; usually computed by deconv algorithms for next iteration,
; compared with image_data for goodness of fit. If not supplied
; it will be computed (option is just to save time).
; FT_PSF = optional input of the Fourier transform of PSF (to save time),
; if not supplied it will be computed and passed on for next time.
; OUTPUTS:
; deconv_info = structure containing deconv statistics is updated.
; EXTERNAL CALLS:
; pro tvs
; function histo
; function convolve
; function get_window
; function positivity
; function deconv_struct
; function FullWid_HalfMax
; COMMON BLOCKS:
; common deconv_analyze, windc, winrh, winac
; common deconv_analyz1, magfim, ixsiz, iysiz, magfr, blankim, posxy
; common deconv_analyz2, min_rstd, max_rstd
; common deconv_analyz3, Log_im_fact
; HISTORY:
; Written: Frank Varosi NASA/GSFC 1992.
;-
pro DeConv_Analyze, image_data, psf_obs, image_deconv, deconv_info, INIT=init, $
RE_CONVOL_IMAGE=im_reconv, FT_PSF=psf_FT, NO_FT=noft, $
LAG_CHISQ=Lag, LAG_MIN=Lag_min, NO_RESIDUALS=noresid, $
LOG10_MIN=Log10, TV_MONITOR=tv_mon
common deconv_analyze, windc, winrh, winac
common deconv_analyz1, magfim, ixsiz, iysiz, magfr, blankim, posxy
common deconv_analyz2, min_rstd, max_rstd
common deconv_analyz3, Log_im_fact
if N_struct( deconv_info ) LE 0 then begin
message,"missing: deconv_info (structure), creating blank",/INFO
deconv_info = deconv_struct( 1 )
endif
sim = size( image_data )
Npix = N_elements( image_data )
save_p = !P
!P.font=0
xcs = !D.x_ch_size
ycs = !D.y_ch_size
if keyword_set( init ) then begin
if (deconv_info.noise_model EQ "POISSON") then begin
print," ChiSq Log ChiSq(ratio) "+$
" flux: FWHM"
print,"iter # Resid Likelihood Auto-correl "+$
"Max conserv X Y"
if (Npix GT 9) then $
Log_im_fact = image_data * aLog(image_data>1) $
- image_data
endif else begin
print," Residuals: " + $
" ChiSq(ratio) flux: FWHM"
print,"iter # Mean St.Dev. & ratio" + $
" Auto-correl Max conserv X Y"
endelse
endif
if (N_elements( image_deconv ) LE 9) OR (Npix LE 9) then goto,PRINT
imwork = positivity( image_deconv, EPS=deconv_info.positivity_eps )
imtotdc = total( imwork )
deconv_info.conserv = imtotdc/total( image_data )
deconv_info.fwhm_xy = FullWid_HalfMax( imwork )
if keyword_set( Log10 ) AND keyword_set( tv_mon ) then begin
ims = deconv_info.image_stats
if (ims.sigma_noise GT 0) then $
Log10r = ims.sigma_noise > ims.sky_level $
else Log10r = 10.^( nint( aLog10( ims.max ) ) - 4 )
endif
if keyword_set( init ) AND keyword_set( tv_mon ) then begin
if (max( sim(1:2) ) LT 256) then begin
magfim = nint( 256./max( sim(1:2) ) )
endif else if (max( sim(1:2) ) LT 512) then begin
magfim = nint( 512./max( sim(1:2) ) )
endif else if (max( sim(1:2) ) GT 512) then begin
magfim = 512./max( sim(1:2) )
magfim = 1.0/nint( 1/magfim )
endif else magfim=1
simm = magfim*sim
ixsiz = simm(1)
iysiz = simm(2)
wxsiz = 2*ixsiz
wysiz = iysiz + 1.5*ycs
get_window, windc, TIT="Deconvolution:", XS=wxsiz, XP=0, $
YS=wysiz, YP=!DEVY-wysiz
erase
wshow, windc, ICON=0
tvs, image_data, ixsiz, LOG=Log10r, MAG=magfim, /COLOR_BAR,$
NAME="Data image"
blankim = bytarr( ixsiz, 2 * ycs )
title = "Data: " + deconv_info.data_name
xyouts, 1.5*ixsiz, iysiz + 0.5*ycs, ALIGN=0.5,/DEV, title
if (max( sim(1:2) ) LT 256) then begin
magfr = nint( 512./total( sim(1:2) ) )
endif else if (max( sim(1:2) ) GE 256) then begin
magfr = 256./max( sim(1:2) )
magfr = 1.0/nint( 1/magfr )
endif
simm = magfr*sim
wxsiz = simm(1) + 360
wysiz = simm(2) + 1.5*ycs
get_window, winrh, TIT="Residuals:", XS=wxsiz, XP=!DEVX-wxsiz, $
YS=wysiz, YP=!DEVY-wysiz
erase
wshow,winrh,ICON=0
posxy = [ simm(1)+6*xcs, 2*ycs, wxsiz-2*xcs, wysiz-1.5*ycs ]
wxsiz = simm(1) + ( 256 < simm(2))
get_window,winac,TIT="Auto-Correl:",XS=wxsiz,XP=!DEVX-wxsiz, $
YS=wysiz,YP=!DEVY-2*wysiz-24
erase
wshow,winac,ICON=0
xyouts, xcs, wysiz - ycs, /DEV,"Auto-Correlation of Residuals:"
empty
endif
if keyword_set( tv_mon ) then begin
wset,windc
tvs, imwork, LOG=Log10r, MAG=magfim
tv, blankim, 0, iysiz
title = "DeConv: " + deconv_info.algorithm + $
" - Niter = " + strtrim( deconv_info.Niter, 2 )
xyouts, xcs, !D.y_vsize - ycs, /DEV, title
empty
endif
if keyword_set( noresid ) then goto,PRINT
maximd = max( imwork, MIN=minimd )
deconv_info.mindc = minimd
deconv_info.maxdc = maximd
deconv_info.entropy = $
-total( (aLog( imwork>1 )-aLog( imtotdc ) )* imwork )/imtotdc
if N_elements( im_reconv ) NE N_elements( image_deconv ) then $
im_reconv = convolve( imwork, psf_obs, FT_PSF=psf_FT, NO=noft )
if N_elements( min_rstd ) NE 1 then min_rstd = 3
if N_elements( max_rstd ) NE 1 then max_rstd = 3
if (deconv_info.noise_model EQ "POISSON") then begin
deconv_info.Log_Likelihood = $
total( image_data * aLog(im_reconv>1) -im_reconv -Log_im_fact )
if N_elements( zeps ) NE 1 then zeps = 1.e-29
deconv_info.chisq_fit = $
sqrt( total( (image_data-im_reconv)^2/(im_reconv>zeps) ) )/Npix
;; wp = where( (image_data GT zeps) AND (im_reconv GT zeps), npos)
;; wn = where( (image_data LE zeps) OR (im_reconv LE zeps), nneg)
;; imwork(wp) = image_data(wp)/im_reconv(wp)
;; if (nneg GT 0) then imwork(wn) = 1
;; imwork = (image_data+0.5) * alog(imwork) +im_reconv -image_data
wp = where( im_reconv GT zeps, npos)
wn = where( im_reconv LE zeps, nneg)
imwork(wp) = ( image_data(wp) - im_reconv(wp) )/im_reconv(wp)
if (nneg GT 0) then imwork(wn) = 0
rstdev = stdev( imwork, rmean )
deconv_info.R_Mean = rmean
deconv_info.R_stdev = rstdev
rmin = rmean - min_rstd*rstdev
rmax = rmean + max_rstd*rstdev
sigma = rstdev
imrh = histo( imwork(where( imwork )), imrv, $
NBIN=sqrt( Npix/2 )<40, MIN=rmin, MAX=rmax )
endif else if (deconv_info.noise_model EQ "GAUSSIAN") then begin
imwork = image_data - im_reconv
rstdev = stdev( imwork, rmean )
deconv_info.R_Mean = rmean
deconv_info.R_stdev = rstdev
rmin = rmean - min_rstd*rstdev
rmax = rmean + max_rstd*rstdev
imrh = histo( imwork, imrv, NBIN=sqrt( Npix/2 )<40, $
MIN=rmin, MAX=rmax )
if (deconv_info.image_stats.sigma_noise GT 0) then $
sigma = deconv_info.image_stats.sigma_noise else sigma = rstdev
deconv_info.Log_Likelihood = aLog( 1/sigma/sqrt( 2*!PI ) ) - $
total( imwork^2 )/(2*sigma^2)
deconv_info.chisq_fit = rstdev/sigma
endif else goto,PRINT
if keyword_set( tv_mon ) then begin
wset,winrh
plot, imrv, imrh, PSYM=10, POS=posxy, /DEV, $
XRAN=[rmin,rmax], TIT="Histogram of Residuals"
imrg = exp( -(imrv-rmean)^2/(2*rstdev^2) )
oplot, imrv, imrg*(total( imrh )/total( imrg )), LINE=1
tvs, imwork, MAG=magfr, MIN=rmin, MAX=rmax
xyouts, xcs, !D.y_vsize - ycs, /DEV, title
empty
endif
if N_elements( Lag ) EQ 1 then deconv_info.Lag_chisq = Lag
if N_elements( Lag_min ) EQ 1 then deconv_info.Lag_min = Lag_min
if (deconv_info.Lag_chisq LT 1) then begin
Lag = nint( deconv_info.psf_stats.fwhm>1 ) < (min(sim(1:2))/3)
deconv_info.Lag_chisq = Lag
endif else Lag = deconv_info.Lag_chisq
imwork = convolve( imwork/sigma, /AUTO_CORREL )
s = sim/2 - Lag
L = sim/2 + Lag
imrchi = imwork(s(1):L(1),s(2):L(2))
imrchi(Lag,Lag) = ( imrchi(Lag,Lag) - Npix * sigma^2 )/sqrt(2)
imrchi = imrchi^2 / Npix
if (deconv_info.Lag_min GT 0) then begin
L = deconv_info.Lag_min -1
imrchi( Lag-L:Lag+L, Lag-L:Lag+L ) = 0
endif
deconv_info.R_auto_chi = total( imrchi( Lag:*, Lag ) ) + $
total( imrchi( *, Lag+1:* ) )
if keyword_set( tv_mon ) then begin
wset,winac
tvs, imwork, MAG=magfr
sxy = 256 < fix(!D.y_vsize-1.5*ycs)
tvs, imrchi, !D.x_vsize-sxy, MAG=sxy/(2*Lag+1)
empty
wset,windc
endif
imwork=0
imrchi=0
!P = save_p
PRINT:
Lag = deconv_info.Lag_chisq
mode_chisq = (2*Lag+1)*(Lag+1) -Lag
if (deconv_info.Lag_min GT 0) then begin
L = deconv_info.Lag_min -1
mode_chisq = mode_chisq - ( (2*L+1)*(L+1) -L )
endif
mode_chisq = ( mode_chisq -2 ) > 1
deconv_info.Expect_auto_chi = mode_chisq
if (deconv_info.noise_model EQ "POISSON") then begin
print, deconv_info.Niter, $
deconv_info.ChiSq_Fit, deconv_info.Log_Likelihood, $
deconv_info.R_auto_chi / mode_chisq, $
deconv_info.maxdc, deconv_info.conserv, $
nint( deconv_info.fwhm_xy *10 )/10.,$
FORMAT="( i6, x, 3G10.3, x, G10.3, F8.2, 2F7.1 )"
endif else begin
if N_elements( sigma ) NE 1 then begin
sigma = deconv_info.image_stats.sigma_noise
if (sigma LE 0) then sigma=1
endif
print, deconv_info.Niter, $
deconv_info.R_Mean, deconv_info.R_stdev, $
deconv_info.R_stdev / sigma, $
deconv_info.R_auto_chi / mode_chisq, $
deconv_info.maxdc, deconv_info.conserv, $
nint( deconv_info.fwhm_xy *10 )/10.,$
FORMAT="( i6, x, 2G10.3, F7.2, G12.3, G11.3, F6.2, 2F7.1 )"
endelse
end