;+ ; ; Deconvolve PSF from an SXT image to get a cleaner image. Default ; PSF is delta-function core plus scattering wings, but it's ; possible to use a Moffat core with/without wings too. ; PRO sxt_decon, in_index, in_data, out_index, out_data, xc=xc, yc=yc, $ spower=spower, log10T=log10T, moffat=moffat, no_scat=no_scat, $ scat_fraction=scat_fraction, no_shadow=no_shadow, $ noverbose=noverbose, qstop=qstop ; ; INPUT PARAMETERS: ; in_index = input index structure ; in_data = input image ; ; OPTIONAL INPUT PARAMETERS: ; XC,YC, the center of the FOV in CCD full-res coords (optional) ; LOG10T, estimate of temperature in the region. In SXT_PSF, ; this will override the default scattering amplitude. ; NO_SCAT, sets PSF to be Moffat core only, no scattering. ; The default PSF for SXT_DECON is a delta function core ; plus scattering wings. (Note, the default PSF for ; the program SXT_PSF is a Moffat core plus wings.) ; (Also note that /no_scat will override /moffat) ; MOFFAT, Uses a PSF comprising a Moffat core PLUS SCATTERING. ; If you don't want scattering, then use /no_scat. ; Overrides the default mode, which is delta core plus wings. ; NO_SHADOW, The default mode is to use a scattering wing model ; that includes the "ray" shadows of the filter supports. ; This is currently our best fit to the scattering wing, ; with a slope of -2 in the radial dependence. ; Selecting NO_SHADOW will defeat this, and assume a ; cylindrically symmetric, smooth scattering wing with ; no "rays". ; SPOWER = Slope of the scattering wing, with or without "rays", ; is set by SPOWER. Default SPOWER is -2. ; SCAT_FRACTION = fraction of total signal to be attributed to the ; scattering wings. The default values are based on ; McKenzie's analysis of May 2002; but with this optional ; parameter you can force X percent of the signal into ; the wings; e.g., scat_fraction=0.073 for 7.3% scattering ; ; CALLING SEQUENCE ; sxt_decon,in_index,in_img,out_index,out_img ; sxt_decon,in_ind,indata,out_ind,outdata,xc=144,yc=213,log10t=6.12 ; sxt_decon,in_index,in_image,out_index,out_image,scat_fraction=.044 ; ; OUTPUT PARAMETERS: ; out_index = history-revised index structure (NOT YET!!) ; out_data = deconvolved image ; ; RESTRICTIONS: ; Calls SXT_PSF for point spread function, and then uses Fourier ; transforms for convolution. ; USES BIG ARRAYS, TAKES TIME, USES MEMORY. ; As of May 2002, now allows data cubes (but actually treats ; only one image at a time). ; ; PROGRAMS CALLED: ; GT_FILTB, GT_RES, UNIQ, SORT, N_ELEMENTS, GT_SHAPE, ; WHERE, SXT_PSF, FFT ; ; HISTORY: ; Written June 17, 1994 Barry LaBonte (as SXT_SCATTER) ; Cleaned up the code June 21, 1994 BJL ; Default values are Acton's result for 6-Sep-92 image June 23,1994 BJL ; 30-Oct-96 [LWA] Generalized method of determining pixel size. ; (called LWA_SCATTER) ; 08-Mar-99 DMcK Revised deconvolution from LWA_SCATTER, ; incorporated SXT_PSF. Renamed DMCK_DECON ; 31-Mar-99 DMcK Removed some unused code, renamed SXT_DECON ; 01-Apr-99 DMcK Added DOUBLE=1 to FFT calls ; 07-Nov-99 DMcK Added log10T optional input, for use in SXT_PSF. ; 08-Feb-01 DMcK Included /DELTA, SLEVEL, SPOWER, & SCAT_FRACTION ; as optional inputs for SXT_PSF. Included /WIENER ; for optional Wiener (optimal) filtering. The latter ; calls Tom Metcalf's WIENER.PRO. ; 09-APR-01 DMcK Added switch NO_SHADOW, made SXT_SCATWINGS the ; default scattering-wing model. ;V2.00 15-MAY-02 DMcK Removed defunct SLEVEL and WIENER options; made ; /DELTA the default PSF for deconvolution, but ; allows deconvolution with Moffat core + wings; ; now allows data cubes, and gets image size from ; the input index. Accelerated by calling SXT_PSF ; only once per res/filt/size combination. Force ; total signal to be conserved. Promoted to ; Version 2.00. (Does not yet update history.) ;V2.01 11-JUL-02 DMcK Don't force total signal to be conserved. ;- progverno = 2.010*1000 sttime=systime(1) out_index=in_index ; Setup output arrays out_data=in_data ; ii_N=n_elements(in_index) ss_tote=0 ; For bookkeeping ; Setup PSF core+wing choices if keyword_set(no_scat) then begin ; PSF will be Moffat core only qdelta=0 & qmoff_only=1 endif else if keyword_set(moffat) then begin ; PSF will be Moffat + wings qdelta=0 & qmoff_only=0 endif else begin ; PSF will be delta + wings (default) qdelta=1 & qmoff_only=0 endelse gres=gt_res(in_index) & gfilt=gt_filtb(in_index) resz=gres(uniq(gres,sort(gres))) & n_resz=n_elements(resz) filtz=gfilt(uniq(gfilt,sort(gfilt))) & n_filtz=n_elements(filtz) szz=gt_shape(in_index) & gs1=szz(0,*) s1z=gs1(uniq(gs1,sort(gs1))) & n_s1z=n_elements(s1z) for iires = resz(0),resz(n_resz-1) do begin ; For each res in the set for iifilt = filtz(0),filtz(n_filtz-1) do begin ; For each filtb for iisz = s1z(0),s1z(n_s1z-1) do begin ; For each size ss=where(gres eq iires AND gfilt eq iifilt AND gs1 eq iisz, ss_N) if (ss_N gt 0) then begin ; Call SXT_PSF once per res/filt/size if not(keyword_set(xc)) or not(keyword_set(yc)) then begin psf=sxt_psf(in_index(ss(0)),/four,dim=2*iisz,$ spower=spower,delta=qdelta,moffat_only=qmoff_only,$ log10t=log10t,noverbose=noverbose,$ scat_fraction=scat_fraction,no_shadow=no_shadow) endif else begin psf=sxt_psf(iifilt,xc,yc,res=iires,/four,dim=2*iisz,$ spower=spower,delta=qdelta,moffat_only=qmoff_only,$ log10t=log10t,noverbose=noverbose,$ scat_fraction=scat_fraction,no_shadow=no_shadow) endelse ;print,' Just made PSF for (res,filt,size) = (',$ ;strcompress(iires,/rem),', ',strcompress(iifilt,/rem),', ',$ ;strcompress(iisz,/rem),')' ; Output during debugging ; Take the Fourier Transform psf = FFT( TEMPORARY(psf), -1 ,double=1) for ii = 0,ss_N-1 do begin ; NOW THE IMAGES index=in_index(ss(ii)) in_img=in_data(*,*,ss(ii)) ; renorm_target=total(in_img,/double) ; Prepare to force the ; flux conservation. ; How big is the input? s1 = iisz s2 = szz(1,ss(ii)) IF( s1 NE s2 ) THEN $ PRINT, 'WARNING: Image may not include whole Sun.' ; Insert into double sized array to avoid Fourier wraparound big = FLTARR( 2*s1, 2*s2 ) big(s1/2:s1-1+s1/2, s2/2:s2-1+s2/2) = in_img(0:s1-1,0:s2-1) ; Transform, convolve, and back transform ;if keyword_set(wiener) then begin ; if not(keyword_set(noverbose)) then $ ; box_message,'Implementing Wiener filter' ; w = wiener(big) ; big = FFT( ( FFT(TEMPORARY(big),-1,double=1) * w / psf ),1,double=1) ;endif else begin big = FFT( ( FFT(TEMPORARY(big),-1,double=1) / psf ),1,double=1) ;endelse ; Commented out the Wiener invocation, in case that filter eventually ; is widely distributed. ; Extract the valid region out_img = FLOAT( TEMPORARY( big(s1/2:s1-1+s1/2, s2/2:s2-1+s2/2) ) ) if keyword_set(qstop) then stop ; renorm_begin=total(out_img,/double) ; Forced flux conservation ; renorm_amount=(renorm_begin-renorm_target)/(float(s1)*float(s2)) ; out_img=out_img-renorm_amount ; out_data(0:s1-1,0:s2-1,ss(ii)) = out_img ; Create output variables ;his_index,out_index,ss(ii),'sxt_decon',progverno ss_tote = ss_tote + 1 endfor ; ii loop endif ; ss_N loop endfor ; iisz loop endfor ; iifilt loop endfor ; iires loop if (ss_tote ne ii_N) then begin tbeep,3,0.25 box_message,' Something may be wrong:' print,' Number of input images = ',strcompress(ii_N,/rem) print,' Number of images treated = ',strcompress(ss_tote,/rem) endif if not(keyword_set(noverbose)) then print,' SXT_DECON took ',$ strtrim((systime(1)-sttime)/60.,1),' minutes to complete your request.' END