diff --git a/help/examples/deconvolution_example.pro b/help/examples/deconvolution_example.pro new file mode 100644 index 0000000..2c479e2 --- /dev/null +++ b/help/examples/deconvolution_example.pro @@ -0,0 +1,112 @@ +;;; +;;; Examples of using both FOXSI deconvolution methods. +;;; + +; +; The "simple" method (just a wrapper for max_likelihood.pro) +; + +foxsi, 2014 + +;;;;; Create FOXSI map + +dat = data_lvl2_d6 +cen_map = cen1_pos0 +trange = [t1_pos0_start, t1_pos0_end] +t1 = trange[0] + t_launch +t2 = trange[1] + t_launch +i = where(dat.wsmr_time gt t1 and dat.wsmr_time lt t2 and dat.error_flag eq 0) +map = foxsi_image_map(dat[i],erange=[4.,15.],cen_map,/xy) + +; Find centroid of FOXSI data +center = [33,-233] ; this is the centroid for RHESSI data (microflare 1) +sub_map, map, sub, xrange=[center[0]-100.,center[0]+100.], yrange=[center[1]-100.,center[1]+100.] +cenmap = map_centroid(sub,threshold=0.05*max(sub.data)) + +; Refine centroid of FOXSI data +sub_map, map, sub, xrange=[cenmap[0]-100.,cenmap[0]+100.], yrange=[cenmap[1]-100.,cenmap[1]+100.] +cenmap = map_centroid(sub,threshold=0.05*max(sub.data)) + +; Use centroid to shift map to align with RHESSI data +shiftmap = center - cenmap +map = shift_map(map,shiftmap[0],shiftmap[1]) + +; Restore whichever PSF you want. +restore, GETENV('FOXSIDB') + '/../calibration_data/foxsi_psf_targ1_pos0_367.sav' + +deconv = deconv_foxsi_simple( map, psf_in=psf ) +movie_map, deconv, /noscale + + +; +; The FOXSI custom-built method that uses a maximum likelihood method to find a +; probable source given data in various detector strips +; + +foxsi, 2014 + +pitch = 7.735 ; strip pitch in um +npix = 12. ; number of strips on one side of an image +fov = pitch*(npix-1)/60. ; field of view of the image + +; Choose some FOXSI data. +tr = [t1_pos2_start,t1_pos2_end] ; Time range (example here is first target, pos 2) +d0 = data_lvl2_d0[ where(data_lvl2_d0.error_flag eq 0)] +d1 = data_lvl2_d1[ where(data_lvl2_d1.error_flag eq 0)] +d4 = data_lvl2_d4[ where(data_lvl2_d4.error_flag eq 0)] +d5 = data_lvl2_d5[ where(data_lvl2_d5.error_flag eq 0)] +d6 = data_lvl2_d6[ where(data_lvl2_d6.error_flag eq 0)] + +; This method requires that the data be in their native form - i.e. the counts collected +; in each detector strip (NOT maps made by rebinning the data to rotate all the detectors +; to the same orientation). So make the maps in detector coordinates. Note that the maps +; will all be rotated differently. +map0 = make_map( foxsi_image_det( d0, trange = tr ), dx=pitch, dy=pitch ) +map1 = make_map( foxsi_image_det( d1, trange = tr ), dx=pitch, dy=pitch ) +map4 = make_map( foxsi_image_det( d4, trange = tr ), dx=pitch, dy=pitch ) +map5 = make_map( foxsi_image_det( d5, trange = tr ), dx=pitch, dy=pitch ) +map6 = make_map( foxsi_image_det( d6, trange = tr ), dx=pitch, dy=pitch ) + +; Select a submap with a small FOV (because large ones take forever to run!) +map0 = make_submap( map0, cen=map_centroid(map0, thresh=0.3*max(map0.data)), fov=fov ) +map1 = make_submap( map1, cen=map_centroid(map1, thresh=0.3*max(map1.data)), fov=fov ) +map4 = make_submap( map4, cen=map_centroid(map4, thresh=0.3*max(map4.data)), fov=fov ) +map5 = make_submap( map5, cen=map_centroid(map5, thresh=0.3*max(map5.data)), fov=fov ) +map6 = make_submap( map6, cen=map_centroid(map6, thresh=0.3*max(map6.data)), fov=fov ) + +; This avoids issues like some maps being slightly different sizes. +map1 = coreg_map( map1, map0, /same, /resc ) +map4 = coreg_map( map4, map0, /same, /resc ) +map5 = coreg_map( map5, map0, /same, /resc ) +map6 = coreg_map( map6, map0, /same, /resc ) + +; Append maps into one array. +; Detectors 2 and 3 are just zero arrays; we will only use the Si detectors for now. +maps = replicate( map0, 7 ) +maps[0].data = map0.data +maps[1].data = map1.data +maps[2].data = 0. +maps[3].data = 0. +maps[4].data = map4.data +maps[5].data = map5.data +maps[6].data = map6.data + +; Now the data are prepared! +; Next is to compute the transformation matrix that captures FOXSI's imaging response. +; This step takes a lot of computation! (Several minutes for this example) +; Once you have the matrix file saved, you can use it for any measured image set that has +; the same number of strips, same detectors, and same PSF. +matrix_file = 'workdir/matrix.sav' +measured_dim = npix ; dimensions of (square) measured image +det = [1,1,0,0,1,1,1] +restore, 'calibration_data/foxsi_psf_on_axis.sav', /v ; on-axis PSF + +; Here's the matrix computation. +matrix = deconv_define_matrix( matrix=matrix_file, psf=psf, $ + measured_dim=measured_dim, detector=det ) + +; This step actually does the deconvolution. The output is a set of maps after various +; numbers of deconvolution iterations. It is up to the user to decide what iteration +; number to go to. +deconv = deconv_foxsi( maps, matrix_file, max=20 ) +movie_map, deconv, /nosc diff --git a/help/examples/deconvolution_example_foxsi1.pro b/help/examples/deconvolution_example_foxsi1.pro new file mode 100644 index 0000000..a76e006 --- /dev/null +++ b/help/examples/deconvolution_example_foxsi1.pro @@ -0,0 +1,269 @@ +;;; +;;; Examples of using both FOXSI deconvolution methods. +;;; + +; +; EXAMPLE 1: The "simple" method (just a wrapper for max_likelihood.pro) +; For Detector 6 only. +; See next example for a multi-detector case. +; + +foxsi, 2012 + +;;;;; Create FOXSI map + dat = data_lvl2_d6 + cen_map = cen4 + i = where((dat.wsmr_time ge t4_start and dat.wsmr_time le t4_end) or (dat.wsmr_time ge t6_start and dat.wsmr_time le t6_end) and dat.error_flag eq 0) + map = foxsi_image_map(dat[i],erange=[4.,15.],cen_map,/xy) + + ; Find centroid of FOXSI data + cenmap = map_centroid(map,threshold=0.05*max(map.data)) + + ; Refine centroid of FOXSI data + sub_map, map, sub, xrange=[cenmap[0]-100.,cenmap[0]+100.], yrange=[cenmap[1]-100.,cenmap[1]+100.] + cenmap = map_centroid(sub,threshold=0.05*max(sub.data)) + +; Get the PSF. Because the foxsi_psf routine was written specifically for the FOXSI-1 flare, it is at the correct off-axis angle and rotation. +psf=foxsi_psf( pix = sub.dx, fov=sub.dx*sqrt(n_elements(sub.data))/60. ) + +deconv = deconv_foxsi_simple( map, psf_in=psf, iter=[1,2,5,10,25,50,100,200] ) +movie_map, deconv, /noscale + + +; +; EXAMPLE 2: The "simple" method (just a wrapper for max_likelihood.pro) +; This one uses all the detector data except for Det3. +; It also plots the results if you have an AIA context image. +; + +foxsi, 2012 + +;;;;; Create FOXSI map +cen_map = cen4 +d0 = data_lvl2_d0[ where(data_lvl2_d0.error_flag eq 0)] +d1 = data_lvl2_d1[ where(data_lvl2_d1.error_flag eq 0)] +d2 = data_lvl2_d1[ where(data_lvl2_d2.error_flag eq 0)] +d4 = data_lvl2_d4[ where(data_lvl2_d4.error_flag eq 0)] +d5 = data_lvl2_d5[ where(data_lvl2_d5.error_flag eq 0)] +d6 = data_lvl2_d6[ where(data_lvl2_d6.error_flag eq 0)] + +; Select the data from EITHER target 4 OR target 6. (Both targets are the same position.) +tr4 = [t4_start,t4_end] +tr6 = [t6_start,t6_end] +d0 = d0[ where( (d0.wsmr_time ge tr4[0] and d0.wsmr_time le tr4[1]) or (d0.wsmr_time ge tr6[0] and d0.wsmr_time le tr6[1]) ) ] +d1 = d1[ where( (d1.wsmr_time ge tr4[0] and d1.wsmr_time le tr4[1]) or (d1.wsmr_time ge tr6[0] and d1.wsmr_time le tr6[1]) ) ] +d2 = d2[ where( (d2.wsmr_time ge tr4[0] and d2.wsmr_time le tr4[1]) or (d2.wsmr_time ge tr6[0] and d2.wsmr_time le tr6[1]) ) ] +d4 = d4[ where( (d4.wsmr_time ge tr4[0] and d4.wsmr_time le tr4[1]) or (d4.wsmr_time ge tr6[0] and d4.wsmr_time le tr6[1]) ) ] +d5 = d5[ where( (d5.wsmr_time ge tr4[0] and d5.wsmr_time le tr4[1]) or (d5.wsmr_time ge tr6[0] and d5.wsmr_time le tr6[1]) ) ] +d6 = d6[ where( (d6.wsmr_time ge tr4[0] and d6.wsmr_time le tr4[1]) or (d6.wsmr_time ge tr6[0] and d6.wsmr_time le tr6[1]) ) ] + +map0 = foxsi_image_map( d0, erange=[4.,15.],cen_map,/xy ) +map1 = foxsi_image_map( d1, erange=[4.,15.],cen_map,/xy ) +map2 = foxsi_image_map( d2, erange=[4.,15.],cen_map,/xy ) +map4 = foxsi_image_map( d4, erange=[4.,15.],cen_map,/xy ) +map5 = foxsi_image_map( d5, erange=[4.,15.],cen_map,/xy ) +map6 = foxsi_image_map( d6, erange=[4.,15.],cen_map,/xy ) + +; Select a submap with a small FOV (because large ones take forever to run!) +fov=3. +map0 = make_submap( map0, cen=map_centroid(map0, thresh=0.3*max(map0.data)), fov=fov ) +map1 = make_submap( map1, cen=map_centroid(map1, thresh=0.3*max(map1.data)), fov=fov ) +map2 = make_submap( map2, cen=map_centroid(map1, thresh=0.3*max(map2.data)), fov=fov ) +map4 = make_submap( map4, cen=map_centroid(map4, thresh=0.3*max(map4.data)), fov=fov ) +map5 = make_submap( map5, cen=map_centroid(map5, thresh=0.3*max(map5.data)), fov=fov ) +map6 = make_submap( map6, cen=map_centroid(map6, thresh=0.3*max(map6.data)), fov=fov ) + +; This avoids issues like some maps being slightly different sizes. +map1 = coreg_map( map1, map0, /same, /resc ) +map2 = coreg_map( map2, map0, /same, /resc ) +map4 = coreg_map( map4, map0, /same, /resc ) +map5 = coreg_map( map5, map0, /same, /resc ) +map6 = coreg_map( map6, map0, /same, /resc ) + +temp = [map0,map1,map2,map4,map5,map6] +map = map0 +map.data = total(temp.data,3) +print, 'Total photons: ', total(map.data) +print, 'Photons per detector: ', total(map.data)/6. +; IDL> print, 'Total photons: ', total(map.data) +; Total photons: 7166.12 +; IDL> print, 'Photons per detector: ', total(map.data)/6. +; Photons per detector: 1194.35 + +; Get the PSF. Because the foxsi_psf routine was written specifically for the FOXSI-1 +; flare, it is at the correct off-axis angle and rotation. +psf=foxsi_psf( pix = map0.dx, fov=map0.dx*sqrt(n_elements(map0.data))/60. ) + +deconv = deconv_foxsi_simple( map, psf_in=psf, iter=[1,2,5,10,25,50,50,75,100] ) +movie_map, deconv, /noscale + +; The next part overplots the result on AIA. This requires you to have an AIA map +; already. Replace the directory and file with yours. +i = 5 ; which AIA frame to use. +j = 9 ; which FOXSI frame to use. +restore, getenv('FOXSIPKG')+'/data_2012/aia-maps-flare.sav', /v +cen = map_centroid(aia[i],thresh=0.3*max(aia[i].data)) +fx_map = shift_map( deconv[j], cen[0]-deconv[j].xc+12, cen[1]-deconv[j].yc+2 ) +ref_map = shift_map( deconv[0], cen[0]-deconv[0].xc+7, cen[1]-deconv[0].yc+2 ) +fx_map.data = gauss_smooth(fx_map.data,0.5) +;;;ref_map.data = gauss_smooth(ref_map.data,0.5) ; image doesn't really need smoothing +aia[i].id = 'AIA 94A' +fraction = [30,50,70,90]/100. ; fractions of encircled energy for contour levels + +popen, 'foxsi1_deconv_simple_on_aia', xsi=7, ysi=7, /port + !p.multi=0 + aia_lct, r,g,b, wave=94, /load + reverse_ct + plot_map, aia[i], fov=1.8, cen=cen, col=255, charsi=1.5, /limb, grid=2, gcol=255, $ + lcol=255 + loadct, 7 + plot_map, fx_map, /over, col=128, th=6, $ + levels=levels_encircled_energy( fx_map.data, fraction) + xyouts, 910, -170, 'FOXSI 4-15 keV ', col=128, charsi=1.8 + xyouts, 910, -255, 'Deconvolved '+fx_map.id, col=128, charsi=1.8 + xyouts, 0.2, 0.94, 'Targets 4+6, Det 0,1,2,4,5,6', /norm, charsi=2 + xyouts, 0.3, 0.0, 'Gauss smooth = 0.6', /norm, charsi=2 +pclose +spawn, 'open foxsi1_deconv_simple_on_aia.ps' +spawn, 'ps2eps -s legal foxsi1_deconv_simple_on_aia.ps' +popen, 'foxsi1_raw_on_aia', xsi=7, ysi=7, /port + !p.multi=0 + aia_lct, r,g,b, wave=94, /load + reverse_ct + plot_map, aia[i], fov=1.8, cen=cen, col=255, charsi=1.5, /limb, grid=2, gcol=255, $ + lcol=255 + loadct, 7 + plot_map, ref_map, /over, col=128, th=6, $ + levels=levels_encircled_energy( ref_map.data, fraction) + xyouts, 910, -170, 'FOXSI 4-15 keV ', col=128, charsi=1.8 + xyouts, 910, -255, 'Raw map', col=128, charsi=1.8 + xyouts, 0.2, 0.94, 'Targets 4+6, Det 0,1,2,4,5,6', /norm, charsi=2 + xyouts, 0.38, 0.0, 'No smooth', /norm, charsi=2 +pclose +spawn, 'open foxsi1_raw_on_aia.ps' +spawn, 'ps2eps -s legal foxsi1_raw_on_aia.ps' + + + +; +; The FOXSI custom-built method that uses a maximum likelihood method to find a +; probable source given data in various detector strips +; + +foxsi, 2012 + +pitch = 7.735 ; strip pitch in um +npix = 12. ; number of strips on one side of an image +fov = pitch*(npix-1)/60. ; field of view of the image + +; Choose some FOXSI data. +d0 = data_lvl2_d0[ where(data_lvl2_d0.error_flag eq 0)] +d1 = data_lvl2_d1[ where(data_lvl2_d1.error_flag eq 0)] +d2 = data_lvl2_d1[ where(data_lvl2_d2.error_flag eq 0)] +d4 = data_lvl2_d4[ where(data_lvl2_d4.error_flag eq 0)] +d5 = data_lvl2_d5[ where(data_lvl2_d5.error_flag eq 0)] +d6 = data_lvl2_d6[ where(data_lvl2_d6.error_flag eq 0)] + +; Select the data from EITHER target 4 OR target 6. (Both targets are the same position.) +tr4 = [t4_start,t4_end] +tr6 = [t6_start,t6_end] +d0 = d0[ where( (d0.wsmr_time ge tr4[0] and d0.wsmr_time le tr4[1]) or (d0.wsmr_time ge tr6[0] and d0.wsmr_time le tr6[1]) ) ] +d1 = d1[ where( (d1.wsmr_time ge tr4[0] and d1.wsmr_time le tr4[1]) or (d1.wsmr_time ge tr6[0] and d1.wsmr_time le tr6[1]) ) ] +d2 = d2[ where( (d2.wsmr_time ge tr4[0] and d2.wsmr_time le tr4[1]) or (d2.wsmr_time ge tr6[0] and d2.wsmr_time le tr6[1]) ) ] +d4 = d4[ where( (d4.wsmr_time ge tr4[0] and d4.wsmr_time le tr4[1]) or (d4.wsmr_time ge tr6[0] and d4.wsmr_time le tr6[1]) ) ] +d5 = d5[ where( (d5.wsmr_time ge tr4[0] and d5.wsmr_time le tr4[1]) or (d5.wsmr_time ge tr6[0] and d5.wsmr_time le tr6[1]) ) ] +d6 = d6[ where( (d6.wsmr_time ge tr4[0] and d6.wsmr_time le tr4[1]) or (d6.wsmr_time ge tr6[0] and d6.wsmr_time le tr6[1]) ) ] + + +; This method requires that the data be in their native form - i.e. the counts collected +; in each detector strip (NOT maps made by rebinning the data to rotate all the detectors +; to the same orientation). So make the maps in detector coordinates. Note that the maps +; will all be rotated differently. +map0 = make_map( foxsi_image_det( d0, erange=[4.,15.] ), dx=pitch, dy=pitch ) +map1 = make_map( foxsi_image_det( d1, erange=[4.,15.] ), dx=pitch, dy=pitch ) +map2 = make_map( foxsi_image_det( d2, erange=[4.,15.] ), dx=pitch, dy=pitch ) +map4 = make_map( foxsi_image_det( d4, erange=[4.,15.] ), dx=pitch, dy=pitch ) +map5 = make_map( foxsi_image_det( d5, erange=[4.,15.] ), dx=pitch, dy=pitch ) +map6 = make_map( foxsi_image_det( d6, erange=[4.,15.] ), dx=pitch, dy=pitch ) + +; Select a submap with a small FOV (because large ones take forever to run!) +map0 = make_submap( map0, cen=map_centroid(map0, thresh=0.3*max(map0.data)), fov=fov ) +map1 = make_submap( map1, cen=map_centroid(map1, thresh=0.3*max(map1.data)), fov=fov ) +map2 = make_submap( map2, cen=map_centroid(map1, thresh=0.3*max(map2.data)), fov=fov ) +map4 = make_submap( map4, cen=map_centroid(map4, thresh=0.3*max(map4.data)), fov=fov ) +map5 = make_submap( map5, cen=map_centroid(map5, thresh=0.3*max(map5.data)), fov=fov ) +map6 = make_submap( map6, cen=map_centroid(map6, thresh=0.3*max(map6.data)), fov=fov ) + +; This avoids issues like some maps being slightly different sizes. +map1 = coreg_map( map1, map0, /same, /resc ) +map2 = coreg_map( map2, map0, /same, /resc ) +map4 = coreg_map( map4, map0, /same, /resc ) +map5 = coreg_map( map5, map0, /same, /resc ) +map6 = coreg_map( map6, map0, /same, /resc ) + +; Append maps into one array. +; Detectors 2 and 3 are just zero arrays; we will only use the Si detectors for now. +maps = replicate( map0, 7 ) +maps[0].data = map0.data +maps[1].data = map1.data +maps[2].data = map2.data +maps[3].data = 0. +maps[4].data = map4.data +maps[5].data = map5.data +maps[6].data = map6.data + +; Now the data are prepared! +; Next is to compute the transformation matrix that captures FOXSI's imaging response. +; This step takes a lot of computation! (Several minutes for this example) +; Once you have the matrix file saved, you can use it for any measured image set that has +; the same number of strips, same detectors, and same PSF. +matrix_file = 'matrix_foxsi1_flare.sav' +measured_dim = npix ; dimensions of (square) measured image +det = [1,1,1,0,1,1,1] +psf=foxsi_psf( pix = 1, fov=5. ) + +; Here's the matrix computation. +; THIS IS SLOW! IF YOU HAVE ALREADY DONE IT, DON'T REPEAT. +matrix = deconv_define_matrix( matrix=matrix_file, psf=psf, $ + measured_dim=measured_dim, detector=det ) + +; This step actually does the deconvolution. The output is a set of maps after various +; numbers of deconvolution iterations. It is up to the user to decide what iteration +; number to go to. +deconv = deconv_foxsi( maps, matrix_file, max=100 ) +movie_map, deconv, /nosc + +print, 'Total photons: ', total(maps.data) +print, 'Photons per detector: ', total(maps.data)/6. + +;IDL> print, 'Total photons: ', total(maps.data) +;Total photons: 6181.00 +;IDL> print, 'Photons per detector: ', total(maps.data)/6. +;Photons per detector: 1030.17 + +; The next part overplots the result on AIA. This requires you to have an AIA map +; already. Replace the directory and file with yours. +i = 5 ; which AIA frame to use. +j = 20 ; which FOXSI frame to use. +restore, getenv('FOXSIPKG')+'/data_2012/aia-maps-flare.sav', /v +cen = map_centroid(aia[i],thresh=0.3*max(aia[i].data)) +fx_map = shift_map( deconv[j], cen[0]-deconv[j].xc+5, cen[1]-deconv[j].yc+3 ) +aia[i].id = 'AIA 94A' +fraction = [30,50,70,90]/100. ; fractions of encircled energy for contour levels + +popen, 'foxsi1_deconv_on_aia', xsi=7, ysi=7, /port + !p.multi=0 + aia_lct, r,g,b, wave=94, /load + reverse_ct + plot_map, aia[i], fov=1.8, cen=cen, col=255, charsi=1.5, /limb, grid=2, gcol=255, lcol=255 + loadct, 7 + plot_map, fx_map, /over, col=128, th=6, $ + levels=levels_encircled_energy( fx_map.data, fraction) + xyouts, 910, -170, 'FOXSI 4-15 keV ', col=128, charsi=1.8 + xyouts, 910, -255, 'Deconvolved '+fx_map.id, col=128, charsi=1.8 + xyouts, 0.2, 0.94, 'Targets 4+6, Det 0,1,2,4,5,6', /norm, charsi=2 + xyouts, 0.38, 0.0, 'No smooth', /norm, charsi=2 +pclose +spawn, 'open foxsi1_deconv_on_aia.ps' +spawn, 'ps2eps -s legal foxsi1_deconv_on_aia.ps' + diff --git a/img/deconv_define_matrix.pro b/img/deconv_define_matrix.pro new file mode 100644 index 0000000..7662bd9 --- /dev/null +++ b/img/deconv_define_matrix.pro @@ -0,0 +1,172 @@ +;+ +; NAME: DECONV_DEFINE_MATRIX +; +; PURPOSE: +; This routine computes a transformation matrix describing FOXSI's imaging instrument +; response. It is used for deconvolving FOXSI images. It defines a matrix that +; computes the counts in each measured detector strip i due to each source pixel j. +; Source and expected image pixels can be sized differently. No effective area or +; spectral response is included. +; +; Notes: +; -- The PSF must be >=3x the source image dimensions, and square. +; -- No actual source is needed at this point. This computes the transformation +; matrix that can be used for any source. This matrix can be applied to any +; detector image with the specified pixel pitch and dimension. +; -- It's not fast! Computing the matrix with defaults takes several minutes +; on Lindsay's laptop. There are several ways the routine could be developed to +; compute the matrix more quickly. +; -- The matrix roughly conserves flux for a centered source, so that most source and +; measured images on each detector will roughly have the same number of photons. +; (100k photons/source => 100k photons on each detector.) This might not be the +; physical case -- for off-axis angles about half the flux might fall of the small +; region of the detector considered here, but getting it right requires better +; knowledge of the PSF wings than we currently have. Please note that the effective +; area of the instrument is NOT included in this routine. +; +; KEYWORDS: +; PSF: plot_map structure containing the point spread function. The deconvolved +; image will have a pixel size equal to the PSF pixel size. If no PSF is +; supplied then the default is an on-axis PSF with 1" pixels. +; The PSF FOV should be >3x the image FOV (MEASURED_DIM x STRIP PITCH) in order +; to allow for rotations and "full reach." +; across the detector. +; MATRIX_FILE: Save the matrix variable in an IDL save file of this name. +; Default is 'matrix.sav' in current directory. +; To use this routine for anything useful, the matrix must be saved. +; MEASURED_DIM: Number of strips across the image. If this number is n, the 2D image +; has n^2 strip crossings. Default: 10 strips +; This MUST match the size of the image you want to deconvolve. +; Source map dimensions are calculated from this variable (and pitch +; and PSF pixel size) and so runtime goes as MEASURED_DIM^4. +; PITCH: Detector strip size in arcsec. Default: 7.735 +; This MUST match the strip size of the image you want to +; deconvolve. +; DETECTOR_MASK: Flags to show which detectors you want to include in the transform. +; Detectors with zero flags just get zeroes for their transforms. +; +; HISTORY: +; 2016-oct-30 LG Created, based on my code for FOXSI-SMEX. +; 2020-mar-02 LG Cleaned up code for use by others. +;- + + +FUNCTION DECONV_DEFINE_MATRIX, psf=psf, matrix_file=matrix_file, $ + measured_dim=measured_dim, pitch=pitch, $ + detector_mask = detector_mask, stop=stop + + COMMON FOXSI_PARAM + + default, matrix_file, 'matrix.sav' + default, measured_dim, long(10) ; expected image dimensions in detector pixels + default, pitch, 7.73493 ; detector strip pitch (for expected image) + default, detector_mask, [0,0,0,0,0,0,1] ; which of the 7 detectors to use. + + ; Pull the detector rotations (defined in common file) + rotation = [rot0,rot1,rot2,rot3,rot4,rot5,rot6] + + ; If a PSF was supplied, do some basic checks on it and return if errors are found. + ; If no PSF input then get a default one with 0.8 arcsec pixels, 7' off-axis. + if exist( PSF ) then begin + ; Check basic properties of PSF + psf_size = size( PSF.data ) + if psf_size[1] ne psf_size[2] or psf.dx ne psf.dy then begin + print, 'User-supplied PSF must be square and have square pixels." + return, -1 + endif + source_pix = psf.dx + source_FOV = measured_dim*pitch + source_dim = fix(source_FOV/source_pix) + if psf_size[1]*source_pix lt 3*source_fov then begin + print, 'User-supplied PSF must be greater than three times the source FOV' + return, -1 + endif + endif else begin + source_pix = 0.8 ; default source pixel size in arcsec + source_FOV = measured_dim*pitch + source_dim = fix(source_FOV/source_pix) + psf = foxsi_psf( pix=source_pix, fov=3*source_fov/60. ) + endelse + + ; to match variable names in legacy code sections. + w_dim = long( source_dim ) + h_dim = long( measured_dim ) + + ; Get normalization in order to conserve flux for a source map that is concentrated + ; at its center. Reference is a FOV-size submap of the PSF, which should + ; integrate to 1. + testmap = make_submap( psf, cen=[0.,0.], fov=h_dim*pitch/60. ) + psf.data = psf.data / total(testmap.data) + + ; Define the elements of the transformation matrix. + ; Dimensions are source basis, measurement basis, detectors + matrix = fltarr( w_dim^2, h_dim^2, 7 ) + + print + print, 'Computing transformation matrix.' + print, w_dim, '^2 source pixels, ', source_pix, ' arcsec' + print, h_dim, '^2 strip crossings, ', pitch, ' arcsec' + print + + for det=0, 6 do begin ; loop through detectors + if detector_mask[ det ] eq 0 then continue + + for col=0, w_dim-1 do begin + for row=0, w_dim-1 do begin + + undefine, list + if row eq 0 and col mod 5 eq 0 then print, ' Progress: Det ', det, $ + fix(100*float(col)/w_dim), ' percent.' + shift_psf = shift_map( psf, col*source_pix, row*source_pix ) + sub_map, shift_psf, small_psf, xr=[-0.5*w_dim,1.5*w_dim], $ + yr=[-0.5*w_dim,1.5*w_dim] + + ; This next section does the following: + ; (1) Account for detector rotation by rotating the shifted PSF around the + ; center of the positive quadrant (corresponding to detector area). + ; (2) Take the positive quadrant of the rotated image as the true + ; detector area. + ; (3) Rebin to detector strip pitch. + ; (4) Append the values for each strip for a given detector to a list of + ; all strip values. + ; Note that the PSF 2D array size is larger than the source 2D array size. + ; (Positive quadrant is 1.5^2 times the source array size so that sources + ; at detector corners are not eliminated.) + + center = w_dim*source_pix/2.*[1,1] + rot_psf0 = rot_map( small_psf, -rotation[det], rcen=center ) + rot_psf0.roll_angle=0 + sub_map, rot_psf0, sub0, xr=center[0]+h_dim*pitch*[-1.,1.]/2, $ + yr=center[1]+h_dim*pitch*[-1.,1.]/2 + ; Find the number of source pixels that will divide evenly into the number + ; of detector strips. IDL integer math makes this line work! + dim = fix(h_dim*pitch/source_pix)/h_dim*h_dim + half = (size( sub0.data ))[1]/2 + img = sub0.data[ half-dim/2:half+dim/2, half-dim/2:half+dim/2 ] + frebin = frebin( img, h_dim, h_dim, /total ) + push, list, reform( frebin, h_dim^2) + + ;;; Geometry checks show the uncommented one below to be correct! + ;;; This was determined empirically. + ; matrix[w_dim-col*w_dim-row,*,det] = list + ; matrix[col*w_dim+row,*,det] = list + matrix[row*w_dim+col,*,det] = list + ; matrix[w_dim-row*w_dim-col,*,det] = list + + if keyword_set( stop ) then stop + + endfor + endfor + + endfor ; end detector loop + + ; Save matrix and the source pixel size to file. + ; If the filename wasn't given with the extension, add it. + if strpos( matrix_file, '.sav') lt 0 then matrix_file += '.sav' + save, matrix, source_pix, file= matrix_file + print, 'Transformation matrix saved in ', matrix_file, '.' + + return, matrix ; This is no longer technically needed because the matrix is + ; saved to file, but it is useful for debugging. + +END diff --git a/img/deconv_foxsi.pro b/img/deconv_foxsi.pro index 7e775b3..ae359c6 100644 --- a/img/deconv_foxsi.pro +++ b/img/deconv_foxsi.pro @@ -1,432 +1,119 @@ -; Function to perform deconvolution of FOXSI images using a measured or calculated PSF. - -FUNCTION DECONV_FOXSI, DETECTOR, TIME_RANGE, PIX=PIX, FOV=FOV, ERANGE=ERANGE, ALL=ALL, $ - PLATE_SCALE=PLATE_SCALE, PSF_SMOOTH=PSF_SMOOTH, $ - IMG_SMOOTH=IMG_SMOOTH, ROTATION=ROTATION, $ - MEASURED_PSF = MEASURED_PSF, STOP = STOP, FIRST = FIRST, $ - PSF_map=PSF_map, FIX4=FIX4, iter=iter, $ - RECONV_map = reconv_map, CSTAT=cstat, year=year - - ; Detector should be an index array saying which detectors we're using. - ; Example [0,0,0,0,0,1] for D6. - - default, pix, 1. - default, fov, 2. - default, erange, [5.,12.] - default, plate_scale, 1.3 - ; 2. is a good value for D6, 2. is a good value for D4. - ; 1.3 is nominal value. - default, psf_smooth, 15 ; only applies to measured PSF - default, img_smooth, 15 - default, rotation, 79 - default, measured_psf, 0 - default, iter, [1,2,3,4,5,10,20,40,80,100] - - - if year eq 2012 then flare=[967,-207] $ ; from RHESSI flarelist - else flare=[0,-200] - - ; usable times - if year eq 2012 then t_launch = 64500 else t_launch = 69060 - t4_start = t_launch + 340 ; Target 4 (flare) - t4_end = t_launch + 420. ; slightly altered from nominal 421.2 - t5_start = t_launch + 423.5 ; Target 5 (off-pointing) - t5_end = t_launch + 435.9 - t6_start = t_launch + 438.5 ; Target 6 (flare) - t6_end = t_launch + 498.3 - - default, time_range, [t4_start, t6_end] +;+ +; NAME: DECONV_FOXSI +; +; PURPOSE: +; Perform maximum likelihood deconvolution on a real or simulated FOXSI image (where +; "image" means a subset of strip values for one or more detectors), using a +; Lucy-Richardson style technique. A transformation matrix is needed that computes the +; flux in detector strip j due to source pixel i. Because the transformation matrix +; is time-consuming to compute and can be used for deconvolution of any images that +; share dimensions, pixel size, and desired source resolution, that matrix is +; computed in a separate function and saved to file. That filename must be input to +; this routine using the MATRIX_FILE keyword. +; (See DECONV_DEFINE_MATRIX.) +; +; The matrix should include the transformations for all 7 detectors (but some can be +; zeros if you don't want to use those detectors). +; +; Note that the input maps need to be in detector coordinates (NOT solar coordinates) +; since this is the actual measurement basis. This means that the elements of IMAGE +; contain detector images that should look rotated with respect to each other. +; +; INPUT: +; IMAGE: The FOXSI image map to deconvolve, in a plot_map structure with 7 +; elements for the 7 detectors. These must be in detector coordinates. +; Must be square. +; MATRIX_FILE: IDL save file containing the transformation matrix, named MATRIX, +; and the SOURCE_PIXEL var (source pixel size in arcsec) +; +; OPTIONAL KEYWORDS: +; MAX_ITER: Stop after this many iterations. Default 40. +; +; NOTES: +; - Would be nice, but not essential, to adapt the code to handle non-square images. +; +; +; HISTORY: +; 2016-oct-30 LG New version that pulls in many features of the FOXSI-SMEX +; version. +; 2020-mar-02 LG Cleaning up the code to be usable by others. +;- + + +FUNCTION DECONV_FOXSI, image, matrix_file, max_iter=max_iter, stop=stop + + default, max_iter, 40 - dim = fov*60./pix - - loadct, 5 - - ; Step 1: Define a PSF map normalized to unity - ; Use either the X5 measured values for 7' off-axis or else the 3-component - ; 2D Gaussian PSF that came from a fit to the measurement. - - print - print, 'Defining PSF...' - print - - if keyword_set(measured_psf) then begin - - if year eq 2012 then $ - f=file_search('data_2012/45az*.fits') $ ; file containing measured PSF 7' offaxis - else f=file_search('data_2014/atfocus_1s_10times.fits') - fits_read, f, data, ind - - m0 = make_map( float(data[*,*,0]), dx=plate_scale, dy=plate_scale ) - m1 = make_map( float(data[*,*,1]), dx=plate_scale, dy=plate_scale ) - m2 = make_map( float(data[*,*,2]), dx=plate_scale, dy=plate_scale ) - m3 = make_map( float(data[*,*,3]), dx=plate_scale, dy=plate_scale ) - m4 = make_map( float(data[*,*,4]), dx=plate_scale, dy=plate_scale ) - m5 = make_map( float(data[*,*,5]), dx=plate_scale, dy=plate_scale ) - m6 = make_map( float(data[*,*,6]), dx=plate_scale, dy=plate_scale ) - m7 = make_map( float(data[*,*,7]), dx=plate_scale, dy=plate_scale ) - m8 = make_map( float(data[*,*,8]), dx=plate_scale, dy=plate_scale ) - m9 = make_map( float(data[*,*,9]), dx=plate_scale, dy=plate_scale ) - m=[m0,m1,m2,m3,m4,m5,m6,m7,m8,m9] - m_sum = m0 - m_sum.data = total(m.data, 3) - med = median( m_sum.data ) ; subtract a constant value. - m_sum.data = m_sum.data - med - for i=0, n_elements(m)-1 do begin - med = median(m[i].data) - m[i].data = m[i].data - med - endfor - ; Switch PSF to center of a map. - cen = [250., -326.]/1.3*plate_scale - psf = shift_map( make_submap( m_sum, cen=cen, fov=fov ), -cen[0], -cen[1] ) - psf = rot_map( psf, rotation ) - psf = make_map( frebin( psf.data, dim, dim, /total ), dx=pix, dy=pix ) - psf.data = smooth(psf.data, psf_smooth) - - endif else begin - - ; or, use a gaussian PSF - params = [ 0.0, 0.0, 0.0, 0.9875, 0.218387, 0.0762158, 1.27836, 1.77492, 4.36214,$ - 7.21397, 47.5, 240.314, 0.0 ] ; parameters from Steven - - ;psf1 = psf_gaussian( npix=[120,120], /double, st_dev=[1.27836, 1.77492]/pix ) - ;psf2 = psf_gaussian( npix=[120,120], /double, $ - ; st_dev=2.*[4.36214, 7.21397]/pix )*params[4] - ;psf3 = psf_gaussian( npix=[120,120], /double, $ - ; st_dev=2.*[47.5, 240.314]/pix )*params[4]*params[5] - psf1 = psf_gaussian( npix=[dim,dim], /double, st_dev=[1.27836, 1.77492]/pix ) - psf2 = psf_gaussian( npix=[dim,dim], /double, $ - st_dev=2.*[4.36214, 7.21397]/pix )*params[4] - psf3 = psf_gaussian( npix=[dim,dim], /double, $ - st_dev=2.*[47.5, 240.314]/pix )*params[4]*params[5] - psftest = make_map( psf1+psf2+psf3, dx=pix, dy=pix ) - psf = rot_map( psftest, -45 ) - psf.roll_angle=0. - - endelse - - psf.data = psf.data / total(psf.data) ; Renormalize - - ; Step 2: Retrieve flare data and prep it. - - print - print, 'Retrieving flare data...' - print - - t0 = time_range[0] - t_launch - t1 = time_range[1] - t_launch - - tr = [t0,t1] - - n_det=n_elements( detector[ where( detector gt 0 ) ] ) - - if detector[0] gt 0 then $ - img0=foxsi_image_det( data_lvl2_d0, erange=erange, trange=tr,thr_n=4., year=year) - if detector[1] gt 0 then $ - img1=foxsi_image_det( data_lvl2_d1, erange=erange, trange=tr,thr_n=4., year=year) - if detector[2] gt 0 then $ - img2=foxsi_image_det( data_lvl2_d2, erange=erange, trange=tr,thr_n=4., year=year) - if detector[3] gt 0 then $ - img3=foxsi_image_det( data_lvl2_d3, erange=erange, trange=tr,thr_n=4., year=year) - if detector[4] gt 0 then $ - img4=foxsi_image_det( data_lvl2_d4, erange=erange, trange=tr,thr_n=4., year=year) - if detector[5] gt 0 then $ - img5=foxsi_image_det( data_lvl2_d5, erange=erange, trange=tr,thr_n=4., year=year) - if detector[6] gt 0 then $ - img6=foxsi_image_det( data_lvl2_d6, erange=erange, trange=tr,thr_n=4., year=year) - - ; Rotation angles for all detectors (as designed, no tweaks yet). - rot0 = 82.5 - rot1 = 75. - rot2 = -67.5 - rot3 = -75. - rot4 = 97.5 - rot5 = 90. - rot6 = -60. - -; imdim = 1000/fix(pix) - imdim=128*7.78/pix - - if n_det eq 1 then begin - - i = where( detector ne 0 ) - case i of - 0: begin - img0 = foxsi_image_det( data_lvl2_d0, erange=erange, trange=tr, thr_n=4., year=2012) - map = rot_map( make_map( img0, dx=7.78, dy=7.78 ), rot0 ) - end - 1: begin - img1 = foxsi_image_det( data_lvl2_d1, erange=erange, trange=tr, thr_n=4., year=2012) - map = rot_map( make_map( img1, dx=7.78, dy=7.78 ), rot1 ) - end - 2: begin - img2 = foxsi_image_det( data_lvl2_d2, erange=erange, trange=tr, thr_n=4., year=2012) - map = rot_map( make_map( img2, dx=7.78, dy=7.78 ), rot2 ) - end - 3: begin - img3 = foxsi_image_det( data_lvl2_d3, erange=erange, trange=tr, thr_n=4., year=2012) - map = rot_map( make_map( img3, dx=7.78, dy=7.78 ), rot3 ) - end - 4: begin - img4 = foxsi_image_det( data_lvl2_d4, erange=erange, trange=tr, thr_n=4., year=2012) - if keyword_set( fix4 ) then img4[18,*] = (img4[17,*] + img4[19,*])/2. - map = rot_map( make_map( img4, dx=7.78, dy=7.78 ), rot4 ) - end - 5: begin - img5 = foxsi_image_det( data_lvl2_d5, erange=erange, trange=tr, thr_n=4., year=2012) - map = rot_map( make_map( img5, dx=7.78, dy=7.78 ), rot5 ) - end - 6: begin - img6 = foxsi_image_det( data_lvl2_d6, erange=erange, trange=tr, thr_n=4., year=2012) - map = rot_map( make_map( img6, dx=7.78, dy=7.78 ), rot6 ) - end - endcase - - if keyword_set(stop) then stop - - first = map - first_cen = map_centroid( first, thr=0.1*max(first.data) ) - first = make_submap( first, cen=first_cen, fov=3 ) - ; this gets returned in keyword FIRST to use as sample unprocessed image. - - flare_new=map_centroid( map, thresh=0.3*max(map.data) ) - 45 - raw = shift_map( make_submap( map, cen=flare_new, fov=fov+1), -flare_new[0], -flare_new[1] ) - raw.xc=0 - raw.yc=0 - raw = rebin_map( raw, dim, dim ) - raw.dx = pix - raw.dy = pix - centr = map_centroid( raw ) - raw = make_submap( raw, cen=centr, fov=fov ) - raw.dx = pix - raw.dy = pix - raw.data = smooth(raw.data,img_smooth/pix) - - endif else begin - img0 = foxsi_image_det( data_lvl2_d0, erange=erange, trange=tr, thr_n=4.) - img1 = foxsi_image_det( data_lvl2_d1, erange=erange, trange=tr, thr_n=4.) - img2 = foxsi_image_det( data_lvl2_d2, erange=erange, trange=tr, thr_n=4.) - img3 = foxsi_image_det( data_lvl2_d3, erange=erange, trange=tr, thr_n=4.) - img4 = foxsi_image_det( data_lvl2_d4, erange=erange, trange=tr, thr_n=4.) - if keyword_set( fix4 ) then img4[18,*] = (img4[17,*] + img4[19,*])/2. - img5 = foxsi_image_det( data_lvl2_d5, erange=erange, trange=tr, thr_n=4.) - img6 = foxsi_image_det( data_lvl2_d6, erange=erange, trange=tr, thr_n=4.) - map0_raw = make_map( img0, dx=7.78, dy=7.78 ) - map1_raw = make_map( img1, dx=7.78, dy=7.78 ) - map2_raw = make_map( img2, dx=7.78, dy=7.78 ) - map3_raw = make_map( img3, dx=7.78, dy=7.78 ) - map4_raw = make_map( img4, dx=7.78, dy=7.78 ) - map5_raw = make_map( img5, dx=7.78, dy=7.78 ) - map6_raw = make_map( img6, dx=7.78, dy=7.78 ) - map0 = rebin_map( map0_raw, imdim, imdim ) - map1 = rebin_map( map1_raw, imdim, imdim ) - map2 = rebin_map( map2_raw, imdim, imdim ) - map3 = rebin_map( map3_raw, imdim, imdim ) - map4 = rebin_map( map4_raw, imdim, imdim ) - map5 = rebin_map( map5_raw, imdim, imdim ) - map6 = rebin_map( map6_raw, imdim, imdim ) - map0.dx = pix - map0.dy = pix - map1.dx = pix - map1.dy = pix - map2.dx = pix - map2.dy = pix - map3.dx = pix - map3.dy = pix - map4.dx = pix - map4.dy = pix - map5.dx = pix - map5.dy = pix - map6.dx = pix - map6.dy = pix - map0.data = map0.data/total(map0.data)*total(img0) - map1.data = map1.data/total(map1.data)*total(img1) - map2.data = map2.data/total(map2.data)*total(img2) - map3.data = map3.data/total(map3.data)*total(img3) - map4.data = map4.data/total(map4.data)*total(img4) - map5.data = map5.data/total(map5.data)*total(img5) - map6.data = map6.data/total(map6.data)*total(img6) - ;map0.data = smooth( map0.data, img_smooth/pix ) - ;map1.data = smooth( map1.data, img_smooth/pix ) - ;map2.data = smooth( map2.data, img_smooth/pix ) - ;map3.data = smooth( map3.data, img_smooth/pix ) - ;map4.data = smooth( map4.data, img_smooth/pix ) - ;map5.data = smooth( map5.data, img_smooth/pix ) - ;map6.data = smooth( map6.data, img_smooth/pix ) - map0 = rot_map( map0, rot0 ) - map1 = rot_map( map1, rot1 ) - map2 = rot_map( map2, rot2 ) - map3 = rot_map( map3, rot3 ) - map4 = rot_map( map4, rot4 ) - map5 = rot_map( map5, rot5 ) - map6 = rot_map( map6, rot6 ) - map0_raw = rot_map( map0_raw, rot0 ) - map1_raw = rot_map( map1_raw, rot1 ) - map2_raw = rot_map( map2_raw, rot2 ) - map3_raw = rot_map( map3_raw, rot3 ) - map4_raw = rot_map( map4_raw, rot4 ) - map5_raw = rot_map( map5_raw, rot5 ) - map6_raw = rot_map( map6_raw, rot6 ) - centr0 = map_centroid(map0,th=0.2*max(map0.data)) - 45 - centr1 = map_centroid(map1,th=0.2*max(map1.data)) - 45 - centr2 = map_centroid(map2,th=0.2*max(map2.data)) - 45 - centr3 = map_centroid(map3,th=0.2*max(map3.data)) - 45 - centr4 = map_centroid(map4,th=0.2*max(map4.data)) - 45 - centr5 = map_centroid(map5,th=0.2*max(map5.data)) - 45 - centr6 = map_centroid(map6,th=0.2*max(map6.data)) - 45 - raw0 = shift_map( make_submap( map0, cen=centr0, fov=fov+1), -centr0[0], -centr0[1] ) - raw1 = shift_map( make_submap( map1, cen=centr1, fov=fov+1), -centr1[0], -centr1[1] ) - raw2 = shift_map( make_submap( map2, cen=centr2, fov=fov+1), -centr2[0], -centr2[1] ) - raw3 = shift_map( make_submap( map3, cen=centr3, fov=fov+1), -centr3[0], -centr3[1] ) - raw4 = shift_map( make_submap( map4, cen=centr4, fov=fov+1), -centr4[0], -centr4[1] ) - raw5 = shift_map( make_submap( map5, cen=centr5, fov=fov+1), -centr5[0], -centr5[1] ) - raw6 = shift_map( make_submap( map6, cen=centr6, fov=fov+1), -centr6[0], -centr6[1] ) - - ; choose which dets to include in the image! - raw=raw6 - raw.data = detector[0]*raw0.data + $ - detector[1]*raw1.data + $ - detector[2]*raw2.data + $ - detector[3]*raw3.data + $ - detector[4]*raw4.data + $ - detector[5]*raw5.data + $ - detector[6]*raw6.data - - ; save a raw, totally unprocessed image. - first=map6_raw - first.data = detector[0]*map0_raw.data + $ - detector[1]*map1_raw.data + $ - detector[2]*map2_raw.data + $ - detector[3]*map3_raw.data + $ - detector[4]*map4_raw.data + $ - detector[5]*map5_raw.data + $ - detector[6]*map6_raw.data - first_cen = map_centroid( first, thr=0.1*max(first.data) ) - first = make_submap( first, cen=first_cen, fov=3 ) - ; this gets returned in keyword FIRST to use as sample unprocessed image. - -; raw.xc=0 -; raw.yc=0 -; raw.dx = pix -; raw.dy = pix - - ;raw.data = smooth( raw.data, img_smooth/pix ) - raw = make_submap( raw, cen=[0.,0.], fov=fov ) - - endelse + size = size( image.data ) + if size[1] ne size[2] then begin + print, 'Measured image must be square.' + return, -1 + endif + measured_dim = size[1] + h = reform( double(image.data), measured_dim^2, 7 ) - ; a kludge to fix the problem that im dimensions may one-off. - data_psf = psf.data - data_raw = raw.data - nDimPsf = size( psf.data ) - nDimRaw = size( raw.data ) - if n_elements(raw.data) gt n_elements(psf.data) then begin - sub_map, raw, new, dim=nDimPsf[1:2], ref=psf - raw = new + ; Must supply a transformation matrix file. + if not exist(matrix_file) then begin + print, 'User must supply a matrix file.' + return, -1 endif - if n_elements(psf.data) gt n_elements(raw.data) then begin - sub_map, psf, new, dim=nDimRaw[1:2], ref=raw - psf = new + + ; Restore matrix file and check that contents are appropriate. + restore, matrix_file + if not exist(source_pix) then begin + print, 'Matrix file must also contain source pixel size.' + return, -1 endif - ; Step 3: Do the deconvolution! :D - - print - print, 'Performing deconvolution...' - print - - ;iter = [1,2,3,4,5,10,20,40,80,100] - n = n_elements(iter) - deconv_map = replicate( raw, n+1 ) ; leave an extra at beginning to contain raw image. - reconv_map = replicate( raw, n+1 ) - -; for j=0, n-1 do begin -; undefine, deconv -; print, j, iter[j] -; for i=0, iter[j] do max_likelihood, raw.data, psf.data, deconv, reconv -; deconv_map[j+1].data = deconv -; reconv_map[j+1].data = reconv -; deconv_map[j+1].id=strtrim(iter[j],2)+' iter' -; endfor - - undefine, deconv - undefine, reconv - for j=0, iter[n-1] do begin - max_likelihood, raw.data, psf.data, deconv, reconv - if total( j eq iter ) gt 0 then begin - i = where( j eq iter ) - deconv_map[i+1].data = deconv - reconv_map[i+1].data = reconv - deconv_map[i+1].id=strtrim(iter[i],2)+' iter' - endif - endfor - - print - print, 'Done.' - print - - deconv_map.roll_angle = 0. - - if keyword_set(stop) then stop - - cstat = fltarr(n) - print, 'Cash statistics:' - for j=0, n-1 do print, 'Iteration ', iter[j], ': ', $ - c_statistic( reconv_map[j].data, raw.data ) - for j=0, n-1 do cstat[j] = c_statistic( reconv_map[j].data, raw.data ) - - ch=1.2 - popen, 'deconv-result', xsi=8, ysi=11 - !p.multi=[0,3,4] - if measured_psf eq 1 then tit='D5 measured PSF' else tit='D5 fit PSF' - plot_map, psf, tit=tit+' smooth='+strtrim(psf_smooth,2), charsi=ch - plot_map, raw, tit='Raw FOXSI image', charsi=ch - for j=1, n do plot_map, deconv_map[j], tit=deconv_map[j].id, charsi=ch - - ; Plot on top of AIA 131 image. -; restore, 'data_2012/aia-maps-flare.sav' -; !p.multi=[0,2,2] -; for j=7, 10 do begin -; plot_map, aia[0], fov=1.5 -; plot_map, shift_map(deconv_map[j], flare[0]-6, flare[1]+3), /over, /per,$ -; color=255, thick=4, lev=[10,30,50,70,90] -; endfor - - ; Plot on top of AIA diff image. Find AIA image closest in time. - - restore, 'data_2012/aia_avg.sav' - iavg = closest( anytim(avg.time), anytim('2012-nov-2')+time_range[0] ) - -; ch=1.2 -; !p.multi=[0,2,2] -; for j=7, 10 do begin -; rotmap = rot_map( deconv_map[j], (rotation-90) ) -; rotmap.roll_angle=0. -; plot_map, avg[iavg], fov=1.5, dmin=0., $ -; title = strtrim(iter[j-1],2)+ ' iter' -; plot_map, shift_map(rotmap, flare[0]-6, flare[1]+3), /over, 4, $ -; col=255, thick=4, lev=[20,40,60,80], /per -; endfor - - ch=1.1 - fov=1.5 - !p.multi=[0,2,2] - plot_map, psf, tit='PSF', charsi=ch, fov=fov - plot_map, psf, /over, lev=[50], /per, thick=5, col=255 - xyouts, -40, 35, '50% contour', size=1.2, col=255 - plot_map, raw, tit='Raw image', charsi=ch, fov=fov - plot_map, raw, /over, lev=[50], /per, thick=5, col=255 - xyouts, -40, 35, '50% contour', size=1.2, col=255 - plot_map, deconv_map[8], tit=deconv_map[8].id+', deconvolved', charsi=ch, fov=fov - plot_map, deconv_map[8], /over, lev=[50], /per, thick=5, col=255 - xyouts, -40, 35, '50% contour', size=1.2, col=255 - plot_map, reconv_map[8], tit=deconv_map[8].id+', reconvolved', charsi=ch, fov=fov - plot_map, reconv_map[8], /over, lev=[50], /per, thick=5, col=255 - xyouts, -40, 35, '50% contour', size=1.2, col=255 - - pclose + ; Check that matrix is appropriate for the given data. + size = size( matrix ) + if size[2] ne measured_dim^2 then begin + print, 'Matrix does not match image dimensions.' + return, -1 + endif + + source_dim = sqrt(size[1]) - psf_map = psf + ; Renormalize the detector layers of the matrix so that its total is proportional + ; to the number of counts in each detector. This allows for differing effective + ; areas, but it assumes that differing count numbers are due to that reason only. + for det=0, 6 do matrix[*,*,det] *= total(h[*,det]) + matrix /= total(matrix) +;; stop + + ; Some useful matrices + ; 2D matrix that reshapes the detector dimension + matrix_2D = fltarr( (size(matrix))[1], measured_dim^2*7 ) + for i=0, (size(matrix))[1]-1 do matrix_2D[i,*] = reform( matrix[i,*,*], $ + measured_dim^2*7 ) + inv_matrix = transpose( matrix_2D ) ; this is S_ji (inverse matrix) + + ; + ; Do the deconvolution + ; + + ; Grey-scale initial guess for source W + W0 = dblarr( source_dim^2 )+1./source_dim^2 + ; Let the first map in the output be the grey-scale map. + map = make_map( reform( w0, source_dim, source_dim ), dx=source_pix, dy=source_pix ) + map.id = '0 Iterations' + + w = w0 ; current best guess + C = dblarr( measured_dim^2, 7 ) + for i=1, max_iter do begin + if i mod (max_iter/20) eq 0 then $ + print, 'Iter ', i, ' of', max_iter, ' , C-state = ', c_statistic(h, w) + for j=0,6 do C[*,j] = reform( w#matrix[*,*,j] ) ; 7 reconvolved maps + temp1 = h/c + temp1[ where( finite(temp1) eq 0 ) ] = 0. + temp2 = reform(temp1,measured_dim^2*7)#inv_matrix + w = w*reform(temp2) + map = [map, make_map( reform( w, source_dim, source_dim ), dx=source_pix, $ + dy=source_pix )] + map[i].id = strtrim(i,2)+' Iterations' + if keyword_set( stop ) then stop + endfor - return, deconv_map + return, map -END \ No newline at end of file +END \ No newline at end of file