;+
function trace_undiffract, index, data, satlevel=satlevelin, $
pedestal=pedestalin, guess=guessin, $
guard=guard
;NAME:
; TRACE_UNDIFFRACT
;PURPOSE:
; Remove the diffraction pattern from TRACE images
;CATEGORY:
;CALLING SEQUENCE:
; newimage = trace_undiffract(index,data)
;INPUTS:
; index = data index
; data = image(s). To get enough orders in the diffraction pattern the
; image should be as large as possible given memory/time
; constraints. fltarr(nx,ny,ni)
;OPTIONAL INPUT PARAMETERS:
;KEYWORD PARAMETERS
; satlevel = Saturation level. Any values above this are considered
; saturated. Default = 4050. This is a reasonable value for
; compressed images, but uncompressed images should have
; satlevel set to 4095.
; pedestal = Pedestal level. Subtracted before the deconvolution and
; re-added at the end. If the pedestal has already been
; subtracted, be sure to set this to zero! Default = 86.
; guess = The guess is only important when there are saturated pixels in
; the image. There are 3 possibilities for the initial guess:
; 1. fltarr(nx,ny,ni). User defined initial guess for deconvolved
; image(es). Must be the same dimension as data.
; 2. If guess is a nonzero scalar (i.e. /guess) then the initial
; guess is set to the deconvolution:
; fimage = fft(image,+1)
; fpsf=fft(shift(psf,-floor(nx/2),-floor(ny/2)),+1)
; guess = float(fft(fimage/fpsf,-1))
; /guess is generally the best choice for speed but does not
; work well if there are saturated pixels.
; 3. Otherwise (default), guess is set to the raw image.
; The default is slow, but generally more accurate when there
; are saturated pixels in the image.
; guard = Puts a guard ring of zeroes around the image as padding for the
; FFT. For example, if guard=64, then a 64 pixel ring is put
; around the image and the size of the image is increased by
; 64*2=128 in each dimension. The guard ring is removed before
; the image is returned. Since the TRACE diffraction pattern
; is quite broad, it is *very important* to include a sufficiently
; large guard band. For speed, the final dimension, including the
; guard band, should be a power of 2, since FFT's are used
; throughout.
;OUTPUTS:
; newimage = image with diffraction pattern removed
;COMMON BLOCKS:
;SIDE EFFECTS:
; Extremely slow.
;RESTRICTIONS:
; Any pixels GE satlevel are considered saturated and treated differently
; that the rest of the image. When there are saturated pixels, the
; intensity at these pixels is estimated from the diffraction pattern,
; but it is not known how good (or bad) this estimate is.
;
; The deconvolution is done with the FFT so the dimension of the
; data, including the guard band, should be a power of two for speed.
;
; When there are saturated pixels, the intensity at these pixels is
; estimated from the diffraction pattern, but it is at best a crude
; estimate. Also, the algorithm can't necessarily tell whether the
; diffraction pattern is in fact a real solar feature. Hence the "other"
; solution is for the saturated areas to have a *low* intensity. Watch
; out for this.
;
; This code is still experimental! In particular, the statistics of
; convergence are not consistent from one image to the next. This needs
; to be understood.
;PROCEDURE:
;MODIFICATION HISTORY:
; T. Metcalf 2000-March-24
; T. Metcalf 2002-Oct-25 Changed the way sigma is calculated
; for saturated pixels. Increase sigm
; of saturated pixels and all affected
; pixels by sqrt(sigma).
; T. Metcalf 2003-Jan-30 Fixed bug which did not allow satlevel
; to be passed in.
;-
; experimental
print
print,' ____ __ __ ____ ____ ____ ___ _ _ ____ _ _ _____ __ _'
print,'| ___)\ \/ / | __ \ | ___)| __ \( )| \ / | | ___)| \| |(_ _)/ \ | |'
print,'| _)_ ) ( | __/ | _)_ | / | | | \\// | | _)_ | \\ | | | | () | | |__'
print,'|____)/_/\_\ |_| |____)|_|\_\(___)|_|\/|_| |____)|_|\_| |_| |_||_| |____)'
print
szdata = size(data)
if szdata[0] NE 2 AND szdata[0] NE 3 then $
message,'data must be 2 or 3 dimensional'
if szdata[0] EQ 2 then nimages = 1 $
else nimages = szdata[3]
nx = szdata[1]
ny = szdata[2]
compressed = 1 ; Flag for compressed data: SHOULD BE TAKEN FROM INDEX!!!
new = float(data)
if n_elements(guard) LE 0 then guard = 0L else guard = long(guard>0L)
gnx = nx+guard*2
gny = ny+guard*2
if guard LE 0 then $
message,/info,'WARNING: You set the guard band to 0. This is not really a good idea.'
if n_elements(satlevelin) LE 0 then begin
if compressed then satlevel = 4050. else satlevel = 4095.
endif else satlevel = float(satlevelin)
message,/info,'Saturation level is set to '+string(satlevel)
if n_elements(pedestalin) LE 0 then pedestal=86. $
else pedestal = float(pedestalin)
message,/info,'Pedestal level is set to '+string(pedestal)
for i=0L,nimages-1 do begin
image = fltarr(gnx,gny)
image[guard,guard] = float(data[*,*,i])
; Could mask the saturated areas, but the algorithm can't tell whether the
; diffraction pattern is in fact a real solar feature. Hence the "other"
; solution is for the saturated areas to have a *low* intensity. To prevent
; that from being a valid solution, the guess from trace_invert_diffpatt
; must be preserved and masking the saturated area will fail.
;;mask = where(image LE 0.0 OR image GE satlevel) ; Mask sat + guard band
mask = where(image LE 0.0) ; mask guard band
saturated = where(image GE satlevel,nsaturated)
unsaturated = where(image LT satlevel AND image GT 0.0)
image[unsaturated] = (image[unsaturated]-pedestal)>1.
if abs(float(index[i].wave_len)-171.) LE 5. then begin
message,/info,'Using FeIX and FeX'
psfFeIX = trace_diffraction_pattern(gnx,gny,171.073, $
index[i].cdelt1)
psfFeX = trace_diffraction_pattern(gnx,gny,174.534, $
index[i].cdelt1)
psf = 0.538*psfFeIX + 0.462*psfFeX ; Ratios from Chianti, could be off
endif else begin
psf = trace_diffraction_pattern(gnx,gny,index[i].wave_len, $
index[i].cdelt1)
endelse
psf = psf/total(psf)
; estimate intensity in saturated areas
if nsaturated GT 0 then $
image = trace_invert_diffpatt(image,psf,saturated,sigma=satsigma,satpattern=satpattern)
; Now add in a real PSF by convolving the diffraction pattern with it.
; This is commented out until I really know what the PSF is. Actually
; this is not really tested yet, but the idea is there.
;psf2 = pxn_psf(nddist2([gnx,gny],[floor(gnx/2),floor(gny/2)]),1.0,1.0)
;psf2 = psf2/total(psf2)
;fpsf2 = fft(shift(psf2,-floor(gnx/2),-floor(gny/2)),+1,/double)
;psf3 = fftconvol(fft(psf,+1,/double),fpsf2)>0.
;psf3 = float(psf3)/total(psf3)
;stop
;psf=psf3
if n_elements(guessin) EQ n_elements(data) then begin
guess = fltarr(gnx,gny)
guess[guard,guard] = float(guessin[*,*,i])
message,/info,'Using supplied initial guess'
endif else if n_elements(guessin) EQ 1 then begin
if guessin[0] NE 0 then begin
fimage = fft(image,+1)
fpsf=fft(shift(psf,-floor(gnx/2),-floor(gny/2)),+1)
guess = float(fft(fimage/fpsf,-1))
message,/info,'Using deconvolution as initial guess'
endif else begin
guess = image
message,/info,'Using raw data as initial guess'
endelse
endif else begin
guess = image
message,/info,'Using raw data as initial guess'
endelse
; There are (gain/qe) EUV photons per DN, but after the lumagen, there
; is one photon per electron. So, the error has to
; be adjusted by the sqrt of the gain.
if strpos(strupcase(index[i].amp),'B') GE 0 then gain=42.0
if strpos(strupcase(index[i].amp),'A') GE 0 then gain=12.0
;if abs(float(index[i].wave_len)-171.) LE 5. then qe=0.09
;if abs(float(index[i].wave_len)-195.) LE 5. then qe=0.09 ; Is this right???
;if abs(float(index[i].wave_len)-284.) LE 5. then qe=0.08
;if abs(float(index[i].wave_len)-304.) LE 5. then qe=0.08
;if abs(float(index[i].wave_len)-1215.) LE 5. then qe=0.09
;if abs(float(index[i].wave_len)-1605.) LE 5. then qe=0.15
if n_elements(gain) LE 0 then $
message,"Can't figure out which amplifier was used: " + $
string(index[i].amp)
;ccdbin = sqrt(float(index[i].sum_ccdx)*float(index[i].sum_ccdy)) * $
; float(index[i].bin_ccd)
; The 20^2 is the 20 electron read noise
;sigma = sqrt(image*gain + 20.0^2)/gain/ccdbin
; The following is Ted Tarbell's suggestion. The 33.0 accounts for read
; noise and for dark current subtraction
sigma = sqrt(image*gain + 33.0^2)/gain/float(index[i].bin_ccd)
if nsaturated GT 0 then begin
sigma[satpattern] = sigma[satpattern]*sqrt(sigma[satpattern])
endif
uselog = 1
usepsfcorr = (nsaturated GT 0)
if nsaturated GT 250 then chi2aim = 3.0 $
else if nsaturated GT 1 then chi2aim = 1.0 $
else begin
if compressed then chi2aim = 0.1 else chi2aim = 1.0
endelse
gnew = image_deconvolve(image,psf,sigma, $
maxit=200L, $
uselog=uselog, $
guess=guess, $
mask=mask, $
/float, $
chi2aim=chi2aim, $
usepsfcorr=usepsfcorr, $
chi2limit=0.0, $
fixmemweight=0)
gnew[unsaturated] = gnew[unsaturated]+pedestal
new[*,*,i] = gnew[guard:guard+nx-1,guard:guard+ny-1]
endfor
return,new
end