;+ ;NAME: iris_sg_deconvolve.pro ; ;PURPOSE: Deconvolve empirically determined spatial point spread ; functions from IRIS level 2 spectra. ; ;METHOD: Fourier implementation of iterative Richardson-Lucy ; algorithm. Typically iterations are truncated early ; using a cost function (see fomod keyword) to prevent ; fitting to noise. Some user experimentation may be ; required to find the optimal number of iterations. ; For 4s FUV (1336 and 1394Å) exposures, we found 9-10 ; iterations appropriate, while the higher SNB ratio 4s ; NUV (2796 and 2814Å) exposures could be pushed to 50 ; iterations. Indvidual results may vary depending on ; the input data. ; ; Note that the Richardson-Lucy routine does not handle negative ; values gracefully; therefore this routine thresholds input images ; so that input_image > 0. ; ;SUBROUTINES: ; fft_conv_1d ; ;CALLING SEQUENCE: dcvim=iris_sg_deconvolve(image) ; ;OUTPUT: returns image array after 1 RL iteartion with 2796Å PSF. ; ;EXAMPLE: ; ; dcvim=iris_sg_deconvolve(image,iterations=50,wave=2796) ; Returns dcvim after 50 RL iterations using the PSF ; derived from the 2796Å spectral window. ; ;CALLS: ; ; iris_sg_psfs.sav ; ;INPUT: ; ; image Image to be deconvolved. ; ;OPTIONAL KEYWORDS: ; ; iterations Number of RL iterations (default is 1). ; ; wave Selects PSF for spectral windows 1336Å, 1394Å ; 2796Å, and 2814Å (default is 2796Å). ; ; fft_div Set to deconvolve by division in fourier space. ; iterations keyword is ignored when this option is ; selected. ; ; double Set to use double precision. ; ; psf Outputs PSF used in deconvolution. ; ; org_fun Original input function. Use to continue iterations ; after stopping. (i.e. to compute cost function). ; ;Restrictions: ; ;HISTORY: ; ; Hans Courrier (HTC) 17 March 2018 ; HTC 26 June 2018: Removed PSFs estimated from Voigt profiles. ; S.L.Freeland - tweak to psf file location for $SSW transportability; use $IRIS_RESPONSE ; ;COMMON BLOCKS: ; ; sg_psf_common Common block that contains PSFs and a flag keyword. ; Prevents multiple calls to RESTORE, e.g. when ; using iris_sg_deconvolve in a for-loop. ; ;TODO: ;- FUNCTION FFT_CONV_1D, func, psf, imsize, DIV=div,REV_PSF=rev_psf,DOUBLE=double ; Short function to do FFT convolution only in image y dimension. ; IDL built in convolution programs won't do convolution in ; 1D for 2D images. ; ;OPTIONAL KEYWORD INPUT: ; ; div Divide in FFT space (deconvolve) ; rev_psf Reverse the (1D) input PSF. ; double Use double precision ;Get size of input psf psize = SIZE(psf) psflen = psize[1] ;Check length of PSF against input image (length of PSF is 420 pixels) pads = imsize[2]-psflen ;Get difference of image size and psf length ;Resize PSF if larger than input image IF pads LE 0 THEN BEGIN rs = ABS(pads)/2 IF ABS(pads) MOD 2 THEN pin = psf[rs+1:psflen-rs-1] $ ;Resize PSF odd ELSE pin = psf[rs:psflen-rs-1] ;Resize PSF even pin /= TOTAL(pin, DOUBLE=double) ;Renormalize PSF ENDIF ;Zero pad PSF if smaller than input image IF pads GT 0 THEN BEGIN rs = pads/2 padl = MAKE_ARRAY(rs,DOUBLE=double) IF pads MOD 2 THEN padr = MAKE_ARRAY(rs+1,DOUBLE=double) $ ;Resize PSF odd ELSE padr = MAKE_ARRAY(rs,DOUBLE=double) ;Resize PSF even pin = [padl,psf,padr] ;Add zero padding by concatinating arrays ENDIF ;Resize and Rebin PSF to input image size pin = REBIN(REFORM(pin,1,imsize[2]),imsize[1],imsize[2]) ;Shift PSF center to zero to center output pin = SHIFT(pin,0,-imsize[2]/2) ;Reverse PSF if keyword set IF KEYWORD_SET(rev_psf) THEN BEGIN pin = REVERSE(pin,2) IF ~(psflen MOD 2) THEN pin = SHIFT(pin,0,1) ENDIF ;FFT inputs fpsf = FFT(pin, 1, DOUBLE=DOUBLE, DIMENSION=2) ffunc = FFT(func, 1, DOUBLE=DOUBLE, DIMENSION=2) ;Multiply (or divide) inputs in f space and FFT back to k space IF ~KEYWORD_SET(div) THEN $ output = REAL_PART(FFT((ffunc*fpsf),-1,DOUBLE=DOUBLE,DIMENSION=2)) $ ELSE $ output = REAL_PART(FFT((ffunc/fpsf),-1,DOUBLE=DOUBLE,DIMENSION=2)) RETURN, output END FUNCTION IRIS_SG_DECONVOLVE,image,ITERATIONS=iterations,WAVE=wave,DOUBLE=double, $ FFT_DIV=fft_div,PSF=psf,ORG_FUN=org_fun ;Do just one iteration of Richardson-Lucy if iterations not set IF ~KEYWORD_SET(iterations) THEN iterations=1 ;Set wave keyword to 2814 if not set IF ~KEYWORD_SET(wave) THEN wave=2814 ;Get image dimensions and size imdim = N_DIMENSIONS(image) ;Image dimensions imsize = SIZE(image) ;Image size psflen = 420 ;Length of PSFs ;Return error if input image is not 2D IF imdim NE 2 THEN BEGIN PRINT, 'Something is wrong, input image is not 2D' RETURN, -1 ENDIF ;Check for PSFs in common block, if not restore and place in ;common block (for multiple iterative calls) COMMON sg_psf_common,sg_mark,psf_1336,psf_1394, $ psf_2796,psf_2814 IF ~EXIST(sg_mark) THEN BEGIN psffile=concat_dir('$IRIS_RESPONSE','iris_sg_psfs') RESTGENX, FILE=psffile,sg_psf_1336,sg_psf_1394, $ sg_psf_2796,sg_psf_2814 ;Place PSFs into common block variables psf_1336 = sg_psf_1336 psf_1394 = sg_psf_1394 psf_2796 = sg_psf_2796 psf_2814 = sg_psf_2814 sg_mark=1 ;Set marker flag ENDIF ;Set deconvolution PSF through WAVE keyword CASE wave OF 1336: psf = psf_1336 '1336': psf = psf_1336 1394: psf = psf_1394 '1394': psf = psf_1394 2796: psf = psf_2796 '2796': psf = psf_2796 2814: psf = psf_2814 '2814': psf = psf_2814 ELSE: psf = psf_2796 ;Set default psf to 2814 ENDCASE ;Now that everything is set up, do the convolution IF KEYWORD_SET(fft_div) THEN dcvim = FFT_CONV_1D(image,psf,imsize,/DIV,DOUBLE=double) $ ;FFT decon. ELSE BEGIN ;Richardson-Lucy deconvolution ;Initialize deconvolution dcvim = image > 0 ;Clip negative values from image ;Allow multiple calls by passing the original function back in through fomod IF ARG_PRESENT(org_fun) THEN datain=org_fun ELSE datain = dcvim FOR i=1,iterations DO BEGIN step1 = datain/(FFT_CONV_1D(dcvim,psf,imsize,DOUBLE=double)) dcvim *= FFT_CONV_1D(step1,psf,imsize,DOUBLE=double,/REV_PSF) ENDFOR ENDELSE RETURN, dcvim END