pro desat::get_background_estimation, self_gen, str, info, bg ;;; deconvolution step on un-saturated data ;bg_real = self -> deconv( *str.bg_ind, str, it = str.it, lev = str.lev, dwavelength=dwavelength, /fill) self -> deconv, self_gen, *str.bg_ind, str, it = str.it, lev = str.lev, dwavelength=dwavelength, /fill , bg_real ;;; background estimation routine bg_real_info = (self_gen -> get(/index))[*str.bg_ind] bg = dst_bkg_fourier_interp(bg_real, bg_real_info, info, str) ;return, bg end pro desat::deconv, self_gen, index, str, fill = fill, it = it, lev = lev, dwavelength = dwavelength, flux ;+ ; NAME: ; desat::deconv ; PURPOSE: ; Deconvolution procedure EM based for SDO/AIA images ; EXPLANATION: ; Deconvolution were computed with the whole psf over a reduced FOV for SDO/AIA data ; CALLING SEQUENCE: ; dst_result = obj -> desat( index, fill=fill ) ; INPUTS: ; index = indices of data images in the obj on wich deconvolution has to be performed. ; OUTPUT: ; flux = deconvolved images ; KEYWORDS: ; fill = Deconvolved data were substituted into the obj in index position ; it = Maximum number of iterations ; lev = Stopping rule parameter for EM ; CALLS: ; CALLED BY: ; DESATURATION ; CALLS: ; dst_em_fft ;- default, it , 1000 default, lev, 0.1 info = (self_gen-> get(/index))[0] imag = (self_gen-> get(/data))[*,*,index] flux = imag * 0 dim = [[info.naxis1],[info.naxis2]] ;;; DECONVOLUTION FOR BACKGROUND ESTIMATION for i=0, n_elements(index)-1 do begin print, 'Deconvolution: ' + strtrim(i+1,1), ' of ', strtrim(size(/n_e, index),1) in = { s: ptr_new(findgen(product(dim))),$ g: ptr_new(findgen(product(dim))),$ c_exp: ptr_new(/allocate_heap),$ x: ptr_new(/allocate_heap),$ y: ptr_new(/allocate_heap),$ im: ptr_new(imag[*,*,i]),$ bg_bkup: ptr_new(flux[*,*,i]),$ bg: ptr_new(flux[*,*,i]),$ ns: product(dim),$ npix: str.npix,$ c_stat: 0,$ lev: lev,$ it: it} if size((*self.psf).psf,/n_d) eq 3 then begin message, 'Not Enabled' ; dst_em_fft_mw, in, (*self.psf).psf, it=in.it, level=level ; flux[*,*,i] = reform( total(*in.x, 2), dim ) endif else begin self->dst_em_fft, in, (*self.psf).psf, it=in.it, level=level flux[*,*,i] = reform( *in.x, dim ) endelse imag[*,*,i] = convolve(flux[*,*,i], ((*self.psf).opsf)[*,*,0])>0 endfor if keyword_set(fill) eq 1 then begin dec_info = self_gen-> get(/index) dec = self_gen-> get(/data) dec[*,*,index] = imag index2map, dec_info, dec, map self_gen-> set, index=dec_info, map=map endif ;return, flux end pro desat::identify_saturation_regions, info, str, loud ;+ ; NAME: ; desat::identify_saturation_regions ; PURPOSE: ; By means of the correlation product computed with dst_em_fft routine return the ; position for primary saturated pixels, bloomed pixels, diffraction fringes. ; EXPLANATION: ; Selecting two regions of the image (I_1 with intensities > 15000 and the I_2 with ; intensities < 10000) the correlation into I_1 is computed using I_2 as data. Note that ; in our formulation the correlation product correspond to the first iteration of ; dst_em_fft procedure. Convolving the result of the cross-correlation with the dispersion ; core of the PSF we recognize the primary saturation site S (the one with intensity > saturation ; level) and as a consequence the I_2 pixels that are not in S are the bloomed ones B. ; Moreover, considering the support of the convolution product between mask_s, defined as ; mask_s[s] = 1. and 0 elsewhere, and the psf diffraction component is possible to find out ; the diffraction fringes directly generated by the primary saturation site. ; CALLING SEQUENCE: ; desat::identify_saturation_regions, str, psf_str ; INPUTS: ; str = dst_str global structure for the desaturation routine. ; psf_str = structure containing all the psf components, diffraction, dispersion and ; the complete one. This structure is returned as output from dst_psf_gen.pro ; routine. ; OUTPUT: ; str.s = primary saturated pixels position ; str.z = flaring pixels (not saturated pixels, with high intensity avoided by the method) ; str.g = diffraction fringes pixels position ; str.b = bloomed pixels positions ; str.ns = number of saturated pixel ; str.sat_lev = saturation intensity level ; CALLS: ; CALLED BY: ; DESATURATION ;- cpsf = (*self.psf).cpsf opsf = ((*self.psf).opsf)[*,*,0] im = *str.im bkg = *str.bg bg_bkup = *str.bg_bkup Q_2 = im lt 10000 Q_1 = im gt 15000 I_1 = where( Q_1 ) I_2 = where( Q_2 ) ; data selection mask for correlation mask = im * 0. & mask[I_1] = 1 if size(cpsf, /n_d) eq 3 then begin corr_fring = im * 0 for iw = 0, (size(cpsf, /dim))[2] - 1 do corr_fring += convolve( mask, cpsf[*,*,iw] ) corr_fring = corr_fring ge (1.e-3 * max(cpsf)) endif else begin corr_fring = convolve( mask, cpsf ) ge (1.e-3 * max(cpsf)) endelse corr_g = where(corr_fring and Q_2) in = {s : ptr_new(I_1) ,$ g : ptr_new(corr_g) ,$ im : ptr_new(im) ,$ bg : ptr_new(bkg) ,$ bg_bkup : ptr_new(bg_bkup) ,$ x : ptr_new(/allocate_heap) ,$ y : ptr_new(/allocate_heap) ,$ c_exp : ptr_new(/allocate_heap) ,$ ns : total(Q_1) ,$ npix : str.npix ,$ c_stat : 0 ,$ it : 1 ,$ lev : 1 } if size(cpsf,/n_d) eq 3 then begin message, 'Not Enabled' ; dst_em_fft_mw, in, cpsf, it=it, level=level ; c = im *0. & c[I_1] = total(*in.x, 2) endif else begin self->dst_em_fft, in, cpsf, it=it, level=level c = im *0. & c[I_1] = *in.x endelse ;;; SATURATED PIXEL POSITION c1 = sigma_filter( c, 5, n_sigma = 1, /iterate) ;;;Necessary to avoid artifacts in the cc_bg = bkg * 0 & cc_bg[I_1] = bkg[I_1] cc = (convolve(c1, opsf)) + cc_bg ;;;correlation map print, 'Maximum correlation map:', max(cc) q_s = max(cc) gt str.sat_lev ? cc ge 14000 : cc ge 0.1*max(cc) s = where(q_s) ;;; BLOOMED PIXELS b = where( Q_1 and ~q_s, nbloom) mask_s = im * 0. & mask_s[s] = 1. if size(cpsf, /n_d) eq 3 then begin tmp = im * 0 for iw = 0, (size(cpsf, /dim))[2] - 1 do tmp += convolve( mask_s, cpsf[*,*,iw] ) qfwd_ind = tmp ge (1.e-3 * max(cpsf)) endif else begin qfwd_ind = convolve( mask_s, cpsf ) ge 1.e-3 * max(cpsf) endelse fwd_ind = where( qfwd_ind ) mask_z = im * 0 + 1 & mask_z[ I_2 ] = 0. q_z = convolve(mask_z, opsf) gt min(opsf[where(opsf gt 0)]) z= where(q_z) ;;; DATA (Fringe) PIXELS g = where(Q_2 and qfwd_ind and ~q_z) mg = byte( im * 0 ) & mg[z] = 1b & gg = where( Q_2 and Qfwd_ind and ~mg ) if ~same_data2( gg, g) then stop mask_g = im * 0. & mask_g[g] = 1. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;PLOT;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if loud eq 1 then begin tt = im*0 & tt[I_2] = 1 index2map, info ,((im-bkg)*tt)^0.3, map plot_map, map, /positive, /square ,title = 'Core Masked Bkg-subtracted data', thick=1.5 index2map, info ,cc^0.3, map plot_map, map, /positive, /square ,title = 'Correlation', thick=1.5 endif ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; str.s = ptr_new(s) str.z = ptr_new(z) str.g = ptr_new(g) str.b = ptr_new(b) str.ns = total(q_s) end pro desat::dst_em_fft , str , cpsf , it=it , level=level ;+ ; NAME: ; dst_em_fft ; PURPOSE: ; Deconvolution/desaturation routine based on EM method implemented by using FFT convolution. ; EXPLANATION: ; This routine has a multiple usage in the desaturation method: ; 1) deconvolution method routine for the not saturated image; ; 2) compute the correlation product for determining primary saturated pixels; ; 3) desaturation routine in the case of large saturated areas ( > 1000 ); ; CALLING SEQUENCE: ; dst_em_fft , str , cpsf , it=it , level=level , pow=pow ; INPUTS: ; str = dst_str global structure for the desaturation routine. ; cpsf = diffraction component of the PSF ; OPTIONAL: ; it = maximum iteration number for EM method (default 300) ; level = EM stopping rule parameter (default 0.01) ; OUTPUT: ; str.x = retrieved intensities for saturated pixels ; str.y = array corresponding to the diffraction fringes intensities ; str.c_exp = expectation computed by the method (cpsf*x + bg at the last iteration) ; str.c_stat = c_statistic computed between y and c_exp ; CALLS: ; CALLED BY: ; DST_X_Y ; DECONV ; DESATURATION ; CALLS: ; CONVOLVE, C_SATATISTIC, AVG ;- default , it , str.it default , level , str.lev default , pow , 1 s = *str.s g = *str.g psf = cpsf ;;;ZEROPADDING THE DATA ARRAY TO AN N*N ARRAY WHERE N IS A POWER OF 2 (INCREASES FFT SPEED) dim = size(/dim,*str.im) check = alog(dim)/alog(2) if round(check[0]) ne check[0] or round(check[1]) ne check[1] then dim = [2^(floor(check[0])+1),2^(floor(check[1])+1)] im = fltarr(dim) bkg = fltarr(dim) mask_g = fltarr(dim) x = fltarr(dim) ;;; emp_kkt = fltarr(it+1) exp_kkt = fltarr(it+1) mask = im * 0 & mask[dim[0]/2,dim[1]/2] = 1 im[0,0] = *str.im bkg[0,0] = *str.bg psf = convolve(mask , psf)>0 psf2 = psf * psf x_tmp = *str.im*0. x_tmp[s] = total(*str.bg_bkup) gt 0 ? (*str.bg_bkup)[s] : 1 x[0,0] = x_tmp sat_pos = where(f_div(x , x) gt 0) mask_g_4_zeropad = *str.im * 0. & mask_g_4_zeropad[g] = 1. mask_g[0,0] = mask_g_4_zeropad y = ( im * mask_g )>0. bkg = ( bkg * mask_g )>0. ;;; GT modification V = convolve(mask_g , psf , /corr )>0. for I = 0, it-1 do begin y_i = (convolve( x , psf ) + bkg)>0 Z = (f_div(y , y_i)) * mask_g U = convolve( z , psf , /corr )>0. x[sat_pos] *= ( f_div( U , V ) )[sat_pos] emp_kkt[i] = total( ( x * (V - U))^2.) exp_kkt[i] = total( x^2. * (convolve( f_div(1. , y_i * mask_g)>0 , psf2 , /corr )>0.) ) if emp_kkt[i] le level * exp_kkt[i] then break print, i ,' --- ', emp_kkt[i]/exp_kkt[i] , level endfor dim = size(/dim , *str.im) x = x[0:dim[0]-1 , 0:dim[1]-1] y = y[0:dim[0]-1 , 0:dim[1]-1] y_i = y_i[0:dim[0]-1 , 0:dim[1]-1] C_stat = c_statistic( y , y_i ) print , i , ' R_flx:' , total( y_i[g]) , ' Exp_flx:' , total( y[g] ) , ' C_stat:' , c_stat , ' ' str.x = ptr_new(x[s]) str.y = ptr_new(y) str.c_exp = ptr_new(y_i[g]) str.c_stat = c_stat end pro desat::dst_mtrx , str , cpsf ;+ ; NAME: ; dst_mtrx ; PURPOSE: ; It computes the mtrx matrix for dst_em routine ; EXPLANATION: ; select the components of the convolution matrix that links the saturated pixels str.s ; with the diffraction fringes in str.g ; CALLING SEQUENCE: ; dst_mtrx , str , cpsf ; INPUTS: ; str = dst_str global structure for the desaturation routine. \ ; cps = diffraction component of the PSF ; OUTPUT: ; str.mtrx = mtrx matrix for dst_em routine ; str.y = vector containing the diffraction fringes intensities ; str.bg = vector containing the estimated background intensities for the same pixels ; considered in y ; CALLS: ; CALLED BY: ; DESATURATION ; History: ; RAS, 5-oct-2017, changed name to dst_mtrx from dst_sqm as it's a matrix but not a square matrix ;- s = *str.s g = *str.g n_sat = str.ns dim = size(/dim , *str.im) dim_psf = size(/dim , cpsf) n_g = size(/n_e , g) mtrx = fltarr(n_g , n_sat) ij_sat = array_indices(*str.im , s) im1 = *str.im * 0.0 & im1[dim[0]/2, dim[1]/2] = 1 & c0 = convolve(/corr,im1, cpsf) > 0.0 for i=0, n_sat-1 do mtrx[*,i] = (shift(c0 , ij_sat[0,i] - dim[0]/2 , ij_sat[1,i] - dim[1]/2 ))[g] str.mtrx = ptr_new(mtrx) str.y = ptr_new((*str.im)[g] > 0.) str.bg = ptr_new((*str.bg)[g] > 0.) end pro desat::dst_em, str ,level=level, it=it ;+ ; NAME: ; dst_em ; PURPOSE: ; Computes the intensities of saturated pixels by using an EM algorithm. ; EXPLANATION: ; It computes the intensities for the primary saturated pixels by means EM algorithm ; implemented using the convolution matrix (mtrx). This routine is used for small saturated ; region (x < 1000). ; CALLING SEQUENCE: ; dst_em, str ,level=level, it=it ; INPUTS: ; str = dst_str global structure for the desaturation routine. ; cpsf = diffraction component of the psf ; OPTIONAL: ; it = maximum iteration number for EM method (default 300) ; level = EM stopping rule parameter (default 0.01) ; OUTPUT: ; str.x = retrieved intensities for saturated pixels ; str.y = array corresponding to the diffraction fringes intensities ; str.c_exp = expectation computed by the method (cpsf*x + bg at the last iteration) ; str.c_stat = c_statistic computed between y and c_exp ; CALLS: ; CALLED BY: ; DESATURATION ;- y = *str.y mtrx = *str.mtrx bg = *str.bg default , it , str.it default , level , str.lev dim_mtrx = size(/dim , mtrx) emp_kkt = fltarr(it) exp_kkt = fltarr(it) mtrx2 = mtrx * mtrx ;x = fltarr(str.ns) + 1. ;x = (*str.bg_bkup)[*str.s] x = total((*str.bg_bkup)[*str.s]) gt 0 ? (*str.bg_bkup)[*str.s] : fltarr(str.ns) + 1. ;;; GT modification V = ((y * 0.) + 1.) # mtrx for i = 0 , it-1 do begin y_i = isarray(x) eq 1 ? mtrx # x + bg : mtrx * x + bg z = f_div( y , y_i ) U = reform( z # mtrx ) x = x * f_div( U , V ) emp_kkt[i] = total( ( x * (V - U) )^2. ) exp_kkt[i] = total( x^2. * ( f_div(1.,y_i) # mtrx2 ) ) if emp_kkt[i] le level * exp_kkt[i] then break print, i ,' --- ', emp_kkt[i]/exp_kkt[i] , level endfor C_stat = c_statistic( y , y_i ) print , i , ' R_flx:' , total( y_i ) , ' Exp_flx:' , total( y ) , ' C_stat:' , c_stat , ' ' str.x = ptr_new(x) str.c_exp = ptr_new(y_i) str.c_stat = C_stat end pro desat::dst_image_synthesis, str , mult_fact , psf ;+ ; NAME: ; dst_image_synthesis ; PURPOSE: ; compute the diffraction component to subtract from the image, the convolution between ; the retrieved intensties, inside the saturation region, and the central core of the ; PSF and the intensities for the bloomed pixels. ; EXPLANATION: ; 1) Computes the convolution product between the intensities obtained from the EM method ; and the central core of the psf. ; 2) The diffraction fringes generated by the retrieved intensities are computed and ; subtracted from the original image ; 3) The bloomed intensities are substituted by the estimated background values. ; CALLING SEQUENCE: ; dst_image_synthesis, str , mult_fact , psf ; INPUTS: ; str = dst_str global structure for the desaturation routine. ; mult_fact = multiplicative factor obtained in desaturation routine necessary to ; renormalize the estimated background, minimizing the interpolation error ; made during the interpolation of the zero-th frequency component in ; dst_bkg_fourier_interp.pro ; psf = structure containing all the psf components, diffraction, dispersion and ; the complete one. This structure is returned as output from dst_psf_gen.pro ; routine. ; OUTPUT: ; It fills the x tag of str ; CALLS: ; CALLED BY: ; DESATURATION ; CALLS: ; CONVOLVE ; PROCEDURE: ;- im = *str.im > 0 s = *str.s b = *str.b opsf = psf.opsf cpsf = psf.cpsf mwav = size(/n_dim, *str.x) - 1 nwav = mwav ? (size(/dim, *str.x))[1] : 1 fin_res = mwav ? fltarr(str.npix,str.npix,nwav) : fltarr(str.npix,str.npix) x_im = mwav ? fltarr(str.npix,str.npix,nwav) : fltarr(str.npix,str.npix) f_im = *str.im check = where(b gt 0 , ct) fq = intarr(long(str.npix)*long(str.npix))+1 fq[where(im gt 15000)] = 0 ;if ct gt 1 then fq[[s,b,*str.z]] = 0 else fq[[s,*str.z]] = 0 f = where(fq eq 1) ijs = array_indices(im , s) for iw = 0, nwav - 1 do begin xtmp = fltarr(str.npix,str.npix) b_im = fltarr(str.npix,str.npix) xtmp[s] = (*str.x)[*,iw] f_im[f] -= (convolve(xtmp, cpsf[*,*,iw])>0)[f] if ct gt 0 then b_im[b] = mult_fact*(convolve(*str.bg_bkup,opsf[*,*,iw])>0)[b] xtmp[s] = (convolve(xtmp, opsf[*,*,iw])>0)[s] x_im[*,*,iw] = xtmp + b_im endfor for iw = 0, nwav - 1 do x_im[*,*,iw] += f_im * reform(fq, str.npix,str.npix) str.x = ptr_new(x_im) end function desat::desaturation, wav, ts, te, path, path_save=path_save, sat_lev=sat_lev, lev=lev, it=it, npix=npix,$ core_dim=core_dim , dwavelength = dwavelength, wavstrategy = wavstrategy, peaklam = peaklam, $ use_prep = use_prep, save_fts=save_fts, aec = aec, loud=loud, psplot = psplot, onewindow = onewindow, $ bkg_filename = bkg_filename ;+ ; NAME: ; desat::desaturation ; PURPOSE: ; Desaturation routine for saturated SDO/AIA images ; EXPLANATION: ; Recover the pixel intensity inside the primary saturation region by means a EM ; based routine for the saturated images inside the user defined time interval. ; Moreover, the fringes generated by the recovered components are subtracted and ; bloomed intensity are substituted by the background intensities. ; CALLING SEQUENCE: ; result = obj -> desaturation( wav, ts, te, path, sat_lev=sat_lev, lev=lev, it=it, core_dim=core_dim ,$ ; save_fts=save_fts, aec = aec, DWAVELENGTH = dwavelength) ; INPUTS: ; file = string array containing the file names for the event taken into account ; ts = user defined analysis starting time ; te = user defined analysis end time ; wav = string array containing the wavelength to be process ; OUTPUT: ; ** Structure RESULTS, 2 tags, length=8, data length=8: ; DATA POINTER ; INFO POINTER ; OPTIONAL: ; sat_level = saturation intensity level (default = 16300) ; lev = EM stopping rule parameter (default = 0.01) (for DESAT method) ; it = maximum iteration number (default = 300) (for DESAT method) ; core_dim = radius for the PSF core (default = 5) ; npix = setting the resulting image to be npix*npix (default = 499) ; peaklam = set the wavelength value to use for the generation of the PSF if wavstrategy = 1 ; wavstrategy = wavelength strategy definition ; 0: generate the PSF using the nominal wavelength value of the passband ; 1: generate the PSF using the wavelength value defined by peaklam ; 2: generate the PSF using the wavelength associated to the BRIGHTEST EMISSIVITY ; value computed by flare_peak_wavelength.pro ; bkg_filename= filename of the image to use as background (for DESAT method) ; KEYWORDS: ; DWAVELENGTH - dwavelength is the fractional change in the wavelength over the nominal wavelength ; so for 131, 131.5 dwavelength would be 0.5/131 = 0.00381679 ; PASSBANDWAV - ; USE_PREP - default 1, if set then use aia_prep to register the data and desat cutout outputs ; path_save = path to folder where .fts file will be saved ; save_fts = if set save fts file at the end of the run (default = 1) ; aec = if set will desaturate also sat. images with short exp. time (default = 1) ; loud = if 1 intermediate and final results are shown in windows (default = 1) ; ONEWINDOW - if set and LOUD is set then plots are written and then overwritten in the same graphic window ; CALLS: ; CALLED BY: - ; CALLS: ; readfts ; indices_analysis ; savefts ; data_restore_2 ; aia_psf_gen ; dst_x_y ; dst_mtrx ; dst_EM ; dst_em_FFT ; dst_image_synthesis ; :History: 28-apr-2015, pass USE_PREP to savefts ; 2-may-2015, G Torre, fix xcen and ycen when cutting a data subset ; :Author: gabriele torre and richard schwartz ; 5-oct-2017, changed sqm to mtrx in all uses because it's a matrix but not a square matrix, sqm is reserved for square matrices ; 5-oct-2017, changed the .fts filename format form aia_xxxxx_xxxxxxxx_xxxxxx_xxx.xxx.fts to aia_xxxxx_xxxxxxxx_xxxxxx_xxxpxxx.fts ; ;- ; ;;;DESAT STRUCTURE DEFINITION default, save_fts, 1 default, loud, 1 default, lev, 1 default, it, 5000 default, npix, 499 default, core_dim, 5 default, sat_lev, 15000 default, aec, 1 default, path_save, './' default, bkg_method, 'quadratic' ;requires two observations prior to ts default, psplot, 0 default, onewindow, 0 default, wavstrategy, 0 default, dwavelength, 0.0 default, peaklam, float(wav) default, bkg_filename,'' loadct, 5 results = replicate({dst_result},size(/n_e, wav)) info_results = PTR_NEW(/ALLOCATE_HEAP) for i_wav = 0, size(/n_e, wav)-1 do begin self_gen = obj_new('desat_gen') ; SG added the object desat_gen range = anytim( [ts,te] ) + [-180.0, 180] file = self_gen->desat_gen::find_file(path, wav[i_wav], range = range) str = {dst_str} self_gen ->desat_gen::readfts, file str.file = ptr_new(file) str.ts = ts str.te = te str.bkg_method = bkg_method ;need this for time selection of background if exist(sat_lev) then str.sat_lev = sat_lev if exist(lev) then str.lev = lev if exist(npix) then str.npix = npix if exist(it) then str.it = it if exist(bkg_filename) then str.bkg_filename = bkg_filename[i_wav] ;GT, added [i_wav], 5-oct-2017 self_gen ->desat_gen::indices_analysis, str, aec = aec ; search the files for background estimation (for DESAT) self_gen ->desat_gen::readfts, file ;;;desaturation wavelength strategy definition if wavstrategy eq 0 then begin peaklam = float(wav[i_wav]) endif if wavstrategy eq 1 and ~exist(peaklam) then begin MESSAGE, 'Error: peaklam not defined for strategy ' +strtrim(wavstrategy,1) endif if wavstrategy eq 2 then begin flare_peak_wavelength, float(wav[i_wav]), peaklam = peaklam;, /plotall endif if wavstrategy ne 0 then begin print, "Desaturation performed at lambda = " + string(peaklam) +" Angstrom" dwavelength =fltarr(size(/n_e,peaklam)) for pw = 0 , size(/n_e,peaklam)-1 do dwavelength[pw] = (peaklam[pw] - float(wav))/float(wav) endif self_gen ->desat_gen::_bld_filenames, *str.time_ind, path_save, peaklam = peaklam self_gen ->desat_gen::ctrl_data, str.npix ;set all the maps to npix x npix pixels if keyword_set(save_fts) then self_gen ->desat_gen::savefts, /original, use_prep = use_prep self_gen->desat_gen::get_data_and_info, str, data, info psf_ptr = self_gen ->desat_gen::psf( info, dwavelength, core_dim ) ;;; BACKGROUND ESTIMATION bg = data * 0.0 self.filenames = ptr_new(*self_gen.filenames) self.psf = ptr_new(*self_gen.psf) if strmatch(str.bkg_filename, '') eq 0 then begin mreadfits,concat_dir(path,str.bkg_filename), info_, bg_ bg_ = [[[bg_]],[[bg_]]] info_ = [[info_],[info_]] self_gen ->desat_gen::ctrl_data, str.npix, data = bg_, info = info_ bg = bg_[*,*,0] endif else begin self -> get_background_estimation, self_gen, str, info, bg ; bg = self -> get_background_estimation(str, info) endelse n_rec = size(/n_e,*str.sat_ind) for i=0, n_rec-1 do begin print, 'Desaturation ' + string(info[i].wavelnth) + ':' + string(i+1) + ' of ' + $ string(n_rec), ' Time:', anytim(/atim,info[i].date_obs) str.im = ptr_new(data[*, *, i]) if strmatch(str.bkg_filename, '') eq 0 then begin str.bg = ptr_new(bg) str.bg_bkup = ptr_new(bg) endif else begin str.bg = ptr_new(convolve(bg[*,*,i],((*self.psf).opsf)[*,*,0])>0) str.bg_bkup = ptr_new(bg[*,*,i]) endelse ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;PLOT;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if loud eq 1 then begin numwin = onewindow ? 31 : i title= 'AIA Desaturation Routine ' text = info[i].INSTRUME+' '+strtrim(info[i].WAVELNTH,1)+' '+anytim(info[i].DATE_OBS, /yoh) if ~psplot then window, numwin, xsize = 900, ysize = 600, title= title + text else $ sps, /colo, /landscape !p.multi=[0,3,2] !p.charsize = 1.5 ; plot original image index2map, info[i], *str.im^0.3, map plot_map, map, /positive, /square ,title = 'Original', thick=1.5 index2map, info[i], *str.bg^0.3, map plot_map, map, /positive, /square ,title = 'Background' endif ;psf_ptr = self -> psf( info, 0, core_dim ) self->identify_saturation_regions, info[i], str, loud ;;;FIND SATURATED PIXEL/DATA POSITION ;psf_ptr = self -> psf( info, dwavelength, core_dim ) mask_im = lonarr(n_elements(*str.im)) + 1 mask_im[ [*str.z, *str.z, *str.g ] ] = 0 outer = where( mask_im ) mult_fact = total((*str.im)[outer])/total((*str.bg)[outer]) *str.bg *= mult_fact ;;; EM FOR SMALL SATURATED REGION Fourier = str.ns lt 1000 ? 0 : 1 print,'***** BEGIN EM *****' if size((*self.psf).cpsf,/n_d) eq 2 then begin case Fourier of 0:begin print, '*** mtrx expectation maximization ****' self->dst_mtrx, str, (*self.psf).cpsf self->dst_em , str, level=level, it=it end 1:begin print, '*** FFT expectation maximization ****' self->dst_em_FFT, str, (*self.psf).cpsf, level=level, it=it end endcase endif print,'***** END EM *****' str.sat_flux = str.x self->dst_image_synthesis, str, mult_fact, *self.psf ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;PLOT;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if loud eq 1 and size((*self.psf).cpsf,/n_d) eq 2 then begin mask_g = *str.im*0 & mask_g[*str.g] = *str.c_exp D = (2./n_elements(*str.g))*(*str.im * alog(f_div(*str.im, mask_g)) + mask_g - *str.im ) ; plot residuals index2map, info[i], d, map plot_map, map, /positive, /square ,title = 'Residuals', thick=1.5 ; plot the de-saturated image index2map, info[i] ,(*str.x>0)^0.3, map plot_map, map, /positive, /square ,title = 'Reconstruction', thick=1.5 !p.multi = 0 !p.thick = 1 endif ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; data_fill = (self_gen->get(/data)) data_fill[*,*,(*str.sat_ind)[i]] = *str.x index2map, (self_gen->get(/index)), data_fill, map self_gen -> set, index=(self_gen->get(/index)), map=map self_gen -> set, grid=30, /limb if keyword_set(save_fts) then begin self_gen ->desat_gen::savefts, reform( info[i].date_obs ), info_str = str, original = 0 endif endfor results[i_wav].data = ptr_new((self_gen->get(/data))[*,*,*str.time_ind]) results[i_wav].info = ptr_new((self_gen->get(/index))[*str.time_ind]) endfor if psplot then begin & device, /close & x & endif return, results end pro desat__define, void void={desat, $ filenames: ptr_new(), $ psf: ptr_new(), $ inherits sdo} return end