PRO xrt_deconvolve,idx,da,n_idx,n_da,niter=niter,verbose=verbose,quiet=quiet,psf1keV=psf1keV ; ========================================================================= ;+ ; PROJECT: ; Solar-B / XRT ; ; NAME: ; ; XRT_DECONVOLVE ; ; CATEGORY: ; ; Image Processing ; ; PURPOSE: ; ; Process XRT data with the XRT mirror model point spread function (PSF) ; ; CALLING SEQUENCE: ; ; XRT_DECONVOLVE, idx,da,out_idx,out_da, ; [niter=niter] [,/verbose] ; [psf1keV=psf1keV] ; ; INPUTS: ; ; idx - [Mandatory] (structure array) ; The index structure for the input data ; Can be Level-0 or Level-1, as long as it matches the data. ; da - [Mandatory] (2D number array, [Nx,Ny]) ; This is the 2D data array. If da is not float or double then routine changes da to float. ; Can be Level-0 or Level-1, as long as it matches idx. ; niter - [Optional] (integer > 0) ; Set this value to the number of iterations to perform. Use ; this only when you do not want the default N=5. If N=0, the ; routine performs the default number of iterations. ; ; OUTPUTS: ; ; N_IDX - [Mandatory] (structure array) ; Index for the output data. It will ; have the same content as IDX (after it has been ; through ), except that ; there will be an update to the HISTORY field. ; N_DA - [Mandatory] (2D float, [Nx,Ny]) ; This is the 2D data array for the processed image. ; ; ; KEYWORDS: ; /QUIET [Optional] (Boolean) Set for no messages or errors. Default is ; some messages and all errors. ; ; /VERBOSE - [Optional] (Boolean) If set, print out extra messages, some relevant, others not. ; ; /psf1kev - [Optional] (Boolean) If set, use the 1.0 keV PSF instead of the default ; 560 eV PSF. ; ; ; ; EXAMPLES: ; ; Basic usage: ; IDL> xrt_deconvolve, idx, da, n_idx, n_da ; ; IDL> xrt_deconvolve, idx, da, n_idx, n_da,niter=1,/psf1kev ; ; COMMON BLOCKS: ; ; none ; ; NOTES: ; ; 1) When handling binned data, the PSFs will be binned and re-normalized prior to processing. ; ; 2) When handling small fields of view. The routine will place the input array at the center ; of a float array of the same size as the binned PSF. Users who don't want this option should ; use the PSF with aia_deconvolution_richardsonlucy.pro or with their preferred method. The XRT PSF ; works equally well with the AIA deconvolution routine. ; ; 3) This routine alters the data using fast Fourier transforms. Inspect your data after deconvolution. ; ; 4) There could be edge effects from plopping the a sub FOV image into a zero array. Users who want to ; avoid this can use aia_deconvolution_richardsonlucy.pro with their favorite sub-FOV method. ; ; ; CONTACT: ; ; Comments, feedback, and bug reports regarding this routine may be ; directed to this email address: ; xrt_manager ~at~ head.cfa.harvard.edu ; ; MODIFICATION HISTORY: ; progver = 'v2017-April-14' ;--- (P.Jibben) Written. Most of this code was copied and then modified from xrt_prep.pro, ; mk_xrt_composite.pro and perhaps other XRT routines. ; ; ;- ; ========================================================================= ;=== Initial setup ==================================================== ;===== Initialize program constants. prognam = 'XRT_DECONVOLVE' path0=getenv('SSW_XRT')+'/idl/util/' psf0560='XRT20170324_151721.0.PSF560.fits' psf1000='XRT20170324_161721.0.PSF1000.fits' ;===== Addressing Keywords ;===== Keyword "verbose" overrules "quiet". q_vb = keyword_set(verbose) q_qt = keyword_set(quiet) and (not q_vb) if keyword_set(niter) then niter=niter else niter=5 if keyword_set(psf1kev) then use_psf=psf1000 else use_psf=psf0560 ;===== Announce version of . if (not q_qt) then box_message, 'XRT_DECONVOLVE: Running ' + prognam + ' ' $ + progver + '.' ;=== Check inputs ===================================================== ;===== This block checks to see if the parameters are an index + data array. If no, then abort. if q_vb then print, 'XRT_DECONVOLVE: Checking input parameters...' if (n_params() lt 4) then begin if (not q_qt) then box_message,'XRT_DECONVOLVE: Requires >= 4 parameters.' $ + ' Aborting...' qabort = 1B return ;=== If n_params() ge 4. endif else begin idx_type = datatype(idx) da_type = datatype(da) ;=== Is var=idx an index, and var=da a float? if ( (idx_type eq 'STC') and $ ((da_type eq 'FLO') or (da_type eq 'DOU')) ) then begin if q_vb then print, 'XRT_DECONVOLVE: Recognize inputs as ' $ + 'index + data array.' ;=== Ensure da is of a number type. num_check = total(is_number(da)) eq n_elements(da) ;=== Ensure da has 2 or 3 dimensions. dim_check = (size(da,/n_dim) eq 2) or (size(da,/n_dim) eq 3) if ((idx_type eq 'STC') and num_check and dim_check) then begin q_incase1 = 0B ;=== Does idx and da have the same number of images? Nimg = n_elements(idx) if (size(da,/n_dim) eq 3) $ then Nimg2 = n_elements(da[0,0,*]) $ else Nimg2 = 1 if (Nimg ne Nimg2) then begin if (not q_qt) then box_message, ['XRT_DECONVOLVE: Input ' $ + 'parameters do not match on the number of images.', $ 'Aborting...'] qabort = 1B return endif ;=== Check whether the data array is at least 8x8 in size. Nx = n_elements(da[*,0,0]) Ny = n_elements(da[0,*,0]) if ((Nx lt 8) or (Ny lt 8)) then begin if (not q_qt) then box_message, 'XRT_DECONVOLVE: Input data ' $ + 'images must be at least 8x8 in size. Aborting...' qabort = 1B return endif endif ;=== Not an index, data array. endif else begin if (not q_qt) then box_message, 'XRT_DECONVOLVE: Input parameters ' $ + 'must be index structure and float data array.' $ + ' Aborting...' qabort = 1B return endelse endelse if q_vb then print, ' ...OK' ;===== Verify that idx has necessary tags. if q_vb then print, 'XRT_DECONVOLVE: Checking whether this data has ' $ + 'necessary XRT header tags...' if not required_tags(idx, 'HISTORY,NAXIS1,NAXIS2,' + $ 'CHIP_SUM', missing=missing) then begin if (not q_qt) then box_message, ['XRT_DECONVOLVE: This input header ' $ + 'does not have necessary XRT tags:', $ ' ' + missing, $ ' Aborting...'] qabort = 1B return endif if q_vb then print, ' ' $ + ' ...OK' ;===== This block checks input data size and chip sum. if ( (min(idx.chip_sum) ne max(idx.chip_sum) ) or (min(idx.naxis1) ne max(idx.naxis1) ) $ or ( min(idx.naxis2) ne max(idx.naxis2) ) ) then begin if (not q_qt) then box_message,['XRT_DECONVOLVE: Input data must all be same size ' $ +'and chip sum. Aborting...'] qabort=1B return endif ;===== This block checks whether the data has already been processed by . if q_vb then print, 'XRT_DECONVOLVE: Checking whether this data ' $ + 'has been deconvolved...' hist0=get_history(idx,'XRT_DECONVOLVE',found=found) if found then begin if (not q_qt) then box_message,['XRT_DECONVOLVE: Some of this data has' $ +'already been run through XRT_DECONVOLVE.pro ' $ +'check header tag "history". Aborting...'] qabort = 1B return endif if q_vb then print, ' ' $ + ' ...OK' ;===== This block gets the PSFs for deconvolution. if q_vb then print, 'XRT_DECONVOLVE: Getting PSF... ' ;=== Where in the world is the PSF? psfexist = file_exist(path0+use_psf) if (psfexist ne 1) then begin if (not q_qt) then box_message, $ ['XRT_DECONVOLVE: Cannot find PSF in '+path0+'. Aborting...'] qabort = 1B return endif read_xrt,path0+use_psf,pin,psf,/quiet if q_vb then print,'XRT_DECONVOLVE: Using '+use_psf+' PSF' ;===== place PSF & data on same chip sum. if q_vb then print,'XRT_DECONVOLVE: Checking PSF and input data chip sums.' if (pin.chip_sum ne idx[0].chip_sum) then begin if (not q_qt) then box_message, $ ['XRT_DECONVOLVE: Input data and PSF have different chip sums. Binning PSF...'] psf1=shift_img(psf,[.5,.5]) psf=psf1/total(psf1) tmp=congrid(psf,pin.naxis1/fix(idx[0].chip_sum),pin.naxis2/fix(idx[0].chip_sum),/center,cubic=-.5) psf=tmp/total(tmp) szpsf=size(psf,/dimension) pin.chip_sum=idx[0].chip_sum pin.naxis1=szpsf[0] & pin.naxis2=szpsf[1] endif if q_vb then print,' ...OK' ;===== Data must be nonnegative for science but can be negative for artistic effect. if (min(da) le 0.0) then begin if q_vb then print, 'XRT_DECONVOLVE: Data must be non-negative...' $ +' Setting input data non-negative.' da=da>0.0 endif ;==== This block deconvolves the data with the PSF. n_idx=idx tagval0='(XRT_DECONVOLVE) Deconvolved with '+use_psf+', bin='+strcompress(idx[0].chip_sum,/remove_all) $ +', niter='+strcompress(niter,/remove_all)+'.' for kk=0, Nimg-1 do begin data=da[*,*,kk] ;===== This block checks input data size and drops in center of zero array if not same size as psf if q_vb then print,'XRT_DECONVOLVE: Checking data array size.' extractda=0 if ( (idx[0].naxis1*idx[0].chip_sum ne pin.naxis1*pin.chip_sum) or $ (idx[0].naxis2*idx[0].chip_sum ne pin.naxis1*pin.chip_sum) ) then begin extractda=1 if (not q_qt) then box_message,['Input data not same size as PSF. ' $ +'Dropping image in zero array.'] tmp_da=fltarr(pin.naxis1,pin.naxis2) ddx=.5*idx[0].naxis1 ddy=.5*idx[0].naxis2 xcen=pin.naxis1/2-1 ycen=pin.naxis2/2-1 tmp_da[xcen-ddx:xcen+ddx-1,ycen-ddy:ycen+ddy-1]=data endif else tmp_da=data if q_vb then verbose=1 else verbose=0 tmp=xrt_deconvolve_richardsonlucy(tmp_da,psf,niter=niter,verbose=verbose) if (extractda eq 1) then tmp=tmp[xcen-ddx:xcen+ddx-1,ycen-ddy:ycen+ddy-1] tmp = (tmp < 2500.0) if kk eq 0 then n_da=tmp else n_da=[ [[n_da]],[[tmp]] ] hist_prefix = prognam + ' ' + progver + ': (' + systim() + ') t_hist = hist_prefix + tagval0 if kk eq 0 then hist_entry=t_hist else hist_entry=[hist_entry,t_hist] endfor ;==== Update the history for the image. if q_vb then print, 'XRT_DECONVOLVE: Updating HISTORY...' update_history,n_idx,hist_entry,/noroutine,/mode ;==== This block says goodbye. mathfacts=['You are wrong if you think Mathematics is not fun.','One can cut a pie into 8 pieces with three movements.' , $ 'If you write pi to 2 decimal places backwards it spells pie.','A pizza that has radius z'+string(44B)+' and height a' $ +string(44B)+' has volume Pi x z x z x a.','In a room with 23 people there is a 50% chance that two have the same birthday.', $ 'Mathematics is an anagram of '+string(96b)+'me asthmatic'+string(96b)+'.','Did you know that 111,111,111 x 111,111,111 = ' $ +'12,345,678,987,654,321','You can remember the value of pi=3.1415926 by counting the letters in: May I have a large ' $ +'container of coffee?','Math skills may be something you are born with (or not)'+string(44B)+' study find.', $ 'An icosagon is a shape with 20 sides','A jiffy is 1/1000th of a second.','The Hairy Ball Theorem states that that given a hairy ball, ' $ +'it is impossible to comb the hairs continuously and have them all lay flat. Some hair must stick up.', $ 'An opinion without 3.14159 is just an onion.', 'You can cut a sphere up into pieces, and then reassemble those pieces to get ' $ +'two spheres which are exactly the same as the first sphere.','If you fold a standard piece of paper in half 103 times, the ' $ +'thickness of it will be greater than the size of the observable universe.'] nfacts=n_elements(mathfacts) fact_no=fix(nfacts*randomu(seed,1)) if q_vb then box_message,mathfacts[fact_no] END ;======================================================================