Image Deconvolution for 511 keV Extended Sources

updated on 2025-04-11

This notebook introduces COSI’s image deconvolution with the Compton data space (CDS) in the Galactic coordinate system. An example of the image analysis will be presented using the 511keV thin disk 3-month simulation data created for DC2.

We have two options on the coordinate system to describe the Compton scattering direction (\(\psi\chi\)) in the image deconvolution, namely Galactic coordinates or detector coordinates. Using Galactic coordinates is intuitive, and the spectral fitting adopts this convention. However, image deconvolution in Galactic coordinates take significant computation time since the detector response needs to be converted into Galactic coordinates for each sky pixel. To reduce the burden on the end-user, a pre-computed converted response is provided in DC2 and DC3 for several main sources (511 keV, Al-26, Ti-44, continuum). These pre-computed responses assume that we analyze 3-month data without extracting time intervals, and the pixel resolution of the model map is already fixed. While there is less flexibility in binning/modeling, it is relatively fast to perform the image deconvolution since the most computationally heavy part, the coordinate conversion of the response, can be skipped.

You can also check out the related notebook 511keV-ScAtt-DataReduction.ipynb for imaging in detector coordinates.

[ ]:

import logging
import sys
logger = logging.getLogger('cosipy')
logger.setLevel(logging.INFO)
logger.addHandler(logging.StreamHandler(sys.stdout))

[ ]:

import matplotlib.pyplot as plt

# cosipy uses astropy units
import astropy.units as u

import healpy as hp
from mhealpy import HealpixMap

from histpy import Histogram

from cosipy.data_io import BinnedData
from cosipy.image_deconvolution import  DataIF_COSI_DC2, ImageDeconvolution
from cosipy.util import fetch_wasabi_file

%matplotlib inline

0. Files needed for this notebook

From wasabi

cosi-pipeline-public/COSI-SMEX/develop/data/Responses/PointSourceResponse/psr_gal_511_DC2.h5
- a pre-computed 511 keV line response file converted into the Galactic coordinate system
cosi-pipeline-public/COSI-SMEX/DC2/Data/Sources/511_thin_disk_3months_unbinned_data.fits.gz
cosi-pipeline-public/COSI-SMEX/DC2/Data/Backgrounds/albedo_photons_3months_unbinned_data.fits.gz
- In this notebook, only the albedo gamma-ray background is considered for a tutorial.
- If you want to consider all of the background components, please replace it with cosi-pipeline-public/COSI-SMEX/DC2/Data/Backgrounds/total_bg_3months_unbinned_data.fits.gz
- Note that total_bg_3months_unbinned_data.fits.gz is 14.15 GB.

From docs/tutorials/image_deconvolution/511keV/GalacticCDS

inputs_511keV_DC2.yaml
imagedeconvolution_parfile_gal_511keV.yml

You can download the data and detector response from wasabi. You can skip this cell if you already have downloaded the files.

[ ]:

# Response file:
# wasabi path: COSI-SMEX/develop/Data/Responses/PointSourceResponse/psr_gal_511_DC2.h5
# File size: 3.82 GB
fetch_wasabi_file('COSI-SMEX/develop/Data/Responses/PointSourceResponse/psr_gal_511_DC2.h5', checksum='50eb36c2bc0089ba230af457ec768fa0')

[ ]:

# Source file (511 keV thin disk model):
# wasabi path: COSI-SMEX/DC2/Data/Sources/511_thin_disk_3months_unbinned_data.fits.gz
# File size: 202.45 MB
fetch_wasabi_file('COSI-SMEX/DC2/Data/Sources/511_thin_disk_3months_unbinned_data.fits.gz', checksum='e0db59ba35aabd0f3b9371b85986841d')

[ ]:

# Background file (albedo gamma):
# wasabi path: COSI-SMEX/DC2/Data/Backgrounds/albedo_photons_3months_unbinned_data.fits.gz
# File size: 2.69 GB
fetch_wasabi_file('COSI-SMEX/DC2/Data/Backgrounds/albedo_photons_3months_unbinned_data.fits.gz', checksum='4c125410d8f127d7e12682f008d5651d')

1. Create binned event/background files in the Galactic coordinate system

please modify “path_data” corresponding to your environment.

[ ]:

path_data = "./"
# example
# path_data = "/Users/ckierans/Software/COSItools/COSItools/cosipy/docs/tutorials/image_deconvolution/511keV/GalacticCDS/"

Source

[ ]:

%%time

signal_filepath = path_data + "511_thin_disk_3months_unbinned_data.fits.gz"

binned_signal = BinnedData(input_yaml = "inputs_511keV_DC2.yaml")

binned_signal.get_binned_data(unbinned_data = signal_filepath, psichi_binning="galactic")

Background

[ ]:

%%time

bkg_filepath = path_data + "albedo_photons_3months_unbinned_data.fits.gz"

binned_bkg = BinnedData(input_yaml = "inputs_511keV_DC2.yaml")

binned_bkg.get_binned_data(unbinned_data = bkg_filepath, psichi_binning="galactic")

Convert the data into sparse matrices & add the signal to the background

[ ]:

signal = binned_signal.binned_data.to_dense()
bkg = binned_bkg.binned_data.to_dense()
event = signal + bkg

Save the binned histograms

[ ]:

signal.write("511keV_dc2_galactic_signal.hdf5", overwrite = True)
bkg.write("511keV_dc2_galactic_bkg.hdf5", overwrite = True)
event.write("511keV_dc2_galactic_event.hdf5", overwrite = True)

Load the saved files

[ ]:

signal = Histogram.open("511keV_dc2_galactic_signal.hdf5")
bkg = Histogram.open("511keV_dc2_galactic_bkg.hdf5")
event = Histogram.open("511keV_dc2_galactic_event.hdf5")

In DC2, the number of time bins should be 1 when you perform the image deconvolution using the Galactic CDS. It is because the pre-computed response files in the galactic coordinate have no time axis, and all of the events are assumed to be projected into a single Galactic CDS. In the future, we plan to introduce more flexible binning.

[ ]:

bkg.axes['Time'].edges

2. Load the response matrix

[ ]:

%%time

response_path = path_data + "psr_gal_511_DC2.h5"

image_response = Histogram.open(response_path)

[ ]:

image_response.axes.labels

[ ]:

image_response.contents.shape

3. Imaging deconvolution

Brief overview of the image deconvolution

Basically, we have to maximize the following likelihood function

\[\log L = \sum_i X_i \log \epsilon_i - \sum_i \epsilon_i\]

\(X_i\): detected counts at \(i\)-th bin ( \(i\) : index of the Compton Data Space)

\(\epsilon_i = \sum_j R_{ij} \lambda_j + b_i\) : expected counts ( \(j\) : index of the model space)

\(\lambda_j\) : the model map (basically gamma-ray flux at \(j\)-th pixel)

\(b_i\) : the background at \(i\)-th bin

\(R_{ij}\) : the response matrix

Since we have to optimize the flux in each pixel, and the number of parameters is large, we adopt an iterative approach to find a solution of the above equation. The simplest one is the ML-EM (Maximum Likelihood Expectation Maximization) algorithm. It is also known as the Richardson-Lucy algorithm.

\[\lambda_{j}^{k+1} = \lambda_{j}^{k} + \delta \lambda_{j}^{k}\]

\[\delta \lambda_{j}^{k} = \frac{\lambda_{j}^{k}}{\sum_{i} R_{ij}} \sum_{i} \left(\frac{ X_{i} }{\epsilon_{i}} - 1 \right) R_{ij}\]

We refer to \(\delta \lambda_{j}^{k}\) as the delta map.

As for now, the two improved algorithms are implemented in COSIpy.

Accelerated ML-EM algorithm (Knoedlseder+99)

\[\lambda_{j}^{k+1} = \lambda_{j}^{k} + \alpha^{k} \delta \lambda_{j}^{k}\]

\[\alpha^{k} < \mathrm{max}(- \lambda_{j}^{k} / \delta \lambda_{j}^{k})\]

Practically, in order not to accelerate the algorithm excessively, we set the maximum value of \(\alpha\) (\(\alpha_{\mathrm{max}}\)). Then, \(\alpha\) is calculated as:

\[\alpha^{k} = \mathrm{min}(\mathrm{max}(- \lambda_{j}^{k} / \delta \lambda_{j}^{k}), \alpha_{\mathrm{max}})\]

Noise damping using gaussian smoothing (Knoedlseder+05, Siegert+20)

\[\lambda_{j}^{k+1} = \lambda_{j}^{k} + \alpha^{k} \left[ w_j \delta \lambda_{j}^{k} \right]_{\mathrm{gauss}}\]

\[w_j = \left(\sum_{i} R_{ij}\right)^\beta\]

\(\left[ ... \right]_{\mathrm{gauss}}\) means that the differential image is smoothed by a gaussian filter.

3-1. Prepare DataInterface containing all necessary datasets

[ ]:

%%time

data_interface = DataIF_COSI_DC2.load(name = "511keV",
                                      event_binned_data = event.project(['Em', 'Phi', 'PsiChi']),
                                      dict_bkg_binned_data = {"albedo": bkg.project(['Em', 'Phi', 'PsiChi'])},
                                      rsp = image_response,
                                      coordsys_conv_matrix=None)

3-2. Initialize the instance of the image deconvolution class

First, we prepare an instance of the ImageDeconvolution class and then register the dataset and parameters for the deconvolution. After that, you can start the calculation.

please modify this parameter_filepath corresponding to your environment.

[ ]:

parameter_filepath = "imagedeconvolution_parfile_gal_511keV.yml"

[ ]:

image_deconvolution = ImageDeconvolution()

# set data_interface to image_deconvolution
image_deconvolution.set_dataset([data_interface])

# set a parameter file for the image deconvolution
image_deconvolution.read_parameterfile(parameter_filepath)

3-3. Initialize image_deconvolution

In this process, a model map is defined following the input parameters, and it is initialized. Also, it prepares ancillary data for the image deconvolution, e.g., the expected counts with the initial model map, gaussian smoothing filter etc.

I describe parameters in the parameter file.

model_definition:class

Specify the name of the image class that will be used in the image deconvolution. As for now, it should be “AllSkyImage”.

model_definition:property

Name	Unit	Description	Notes
coordinate	str	the coordinate system of the model map	As for now, it must be ‘galactic’
nside	int	NSIDE of the model map	it must be the same as NSIDE of ‘lb’ axis of the coordinate conversion matrix
scheme	str	SCHEME of the model map	As for now, it must be ‘ring’
energy_edges:value	list of float	The definition of the energy bins of the model map	As for now, it must be the same as that of the response matrix
energy_edges:unit	str	The physical unit of the above list	In most of cases, it is keV.
unit	str	The physical unit of the image	As for now, it should be “”cm-2 s-1 sr-1”

model_definition:initialization

Name	Unit	Description	Notes
algorithm	str	the method name to initialize the model map	As for now, only ‘flat’ can be used
parameter_flat:value	list of float	the list of photon fluxes for each energy band	the length of the list should be the same as the length of “energy_edges” - 1
parameter_flat:unit	str	The physical unit of the above list	As for now, it should be “”cm-2 s-1 sr-1”

deconvolution

Name	Unit	Description	Notes
algorithm	str	the name of the image deconvolution algorithm	As for now, only ‘RL’ is supported

parameter_RL:iteration_max	int	The maximum number of the iteration
parameter_RL:acceleration:activate	bool	whether the accelerated ML-EM algorithm (Knoedlseder+99) is used
parameter_RL:acceleration:accel_factor_max	float	the maximum value for the acceleration parameter
parameter_RL:response_weighting:activate	bool	whether a delta map is renormalized based on the exposure time on each pixel, namely \(w_j = (\sum_{i} R_{ij})^{\beta}\) (see Knoedlseder+05, Siegert+20)
parameter_RL:response_weighting:index	float	\(\beta\) in the above equation
parameter_RL:smoothing:activate	bool	whether a Gaussian filter is used (see Knoedlseder+05, Siegert+20)
parameter_RL:smoothing:FWHM:value	float	the FWHM of the Gaussian in the filter
parameter_RL:smoothing:FWHM:unit	str	the unit of the above value. deg. or rad.
parameter_RL:background_normalization:activate	bool	whether the background normalization factor is optimized at each iteration	As for now, the single background normalization factor is used in all of the bins
parameter_RL:background_normalization:range	list of float	the range of the normalization factor	should be positive
parameter_RL:stopping_criteria:statistics	str	the statistics to be checked	In most of cases, it is “log-likelihood”
parameter_RL:stopping_criteria:threshold	float	the iteration is terminated when the increase of the selected statistics is below this value
parameter_RL:save_results:activate	bool	whether an updated model map, detal map, likelihood etc. are saved after the image deconvolution
parameter_RL:save_results:directory	bool	the directory where the files are saved
parameter_RL:save_results:only_final_result	bool	If True, only the result after the last iteration is saved

[ ]:

image_deconvolution.initialize()

(You can change the parameters as follows)

Note that when you modify the parameters, do not forget to run “initialize” again!

[ ]:

image_deconvolution.override_parameter("deconvolution:parameter:iteration_max = 50")
image_deconvolution.override_parameter("deconvolution:parameter:background_normalization_optimization:activate = True")
image_deconvolution.override_parameter("deconvolution:parameter:acceleration:accel_factor_max = 2.0")
image_deconvolution.override_parameter("deconvolution:parameter:smoothing:activate = False")
#image_deconvolution.override_parameter("deconvolution:parameter:smoothing:activate = True")
image_deconvolution.override_parameter("deconvolution:parameter:response_weighting:activate = False")

image_deconvolution.initialize()

3-4. Start the image deconvolution

With MacBook Pro with M1 Max and 64 GB memory, it takes about 1.5 minutes for 50 iterations.

[ ]:

%%time

image_deconvolution.run_deconvolution()

[ ]:

image_deconvolution.results

4. Analyze the results

Examples to see/analyze the results are shown below.

Log-likelihood

Plotting the log-likelihood vs the number of iterations

[ ]:

x, y = [], []

for result in image_deconvolution.results:
    x.append(result['iteration'])
    y.append(result['log-likelihood'])

plt.plot(x, y)
plt.grid()
plt.xlabel("iteration")
plt.ylabel("log-likelihood")
plt.show()

Alpha (the factor used for the acceleration)

Plotting \(\alpha\) vs the number of iterations. \(\alpha\) is a parameter to accelerate the EM algorithm (see the beginning of Section 4). If it is too large, reconstructed images may have artifacts.

[ ]:

x, y = [], []

for result in image_deconvolution.results:
    x.append(result['iteration'])
    y.append(result['accel_factor'])

plt.plot(x, y)
plt.grid()
plt.xlabel("iteration")
plt.ylabel("accel_factor")
plt.show()

Background normalization

Plotting the background nomalization factor vs the number of iterations. If the backgroud model is accurate and the image is reconstructed perfectly, this factor should be close to 1.

[ ]:

x, y = [], []

for result in image_deconvolution.results:
    x.append(result['iteration'])
    y.append(result['background_normalization']['albedo'])

plt.plot(x, y)
plt.grid()
plt.xlabel("iteration")
plt.ylabel("background_normalization")
plt.show()

The reconstructed images

[ ]:

def plot_reconstructed_image(result, source_position = None): # source_position should be (l,b) in degrees
    iteration = result['iteration']
    image = result['model']

    for energy_index in range(image.axes['Ei'].nbins):
        map_healpxmap = HealpixMap(data = image[:,energy_index], unit = image.unit)

        _, ax = map_healpxmap.plot('mollview')

        _.colorbar.set_label(str(image.unit))

        if source_position is not None:
            ax.scatter(source_position[0]*u.deg, source_position[1]*u.deg, transform=ax.get_transform('world'), color = 'red')

        plt.title(label = f"iteration = {iteration}, energy_index = {energy_index} ({image.axes['Ei'].bounds[energy_index][0]}-{image.axes['Ei'].bounds[energy_index][1]})")

Plotting the reconstructed images in all of the energy bands at the 50th iteration

[ ]:

iteration = 49

plot_reconstructed_image(image_deconvolution.results[iteration])

An example to plot the image in the log scale

[ ]:

iteration_idx = 49

result = image_deconvolution.results[iteration_idx]

iteration = result['iteration']
image = result['model']

data = image[:,0]
data[data <= 0 * data.unit] = 1e-12 * data.unit

hp.mollview(data, min = 1e-5, norm ='log', unit = str(data.unit), title = f'511 keV image at {iteration}th iteration', cmap = 'magma')

plt.show()

Note: when you use a smoothing filter, you may see arfacts like many rings around the galactic center as seen in the next image. This is likely caused by the calculation errors of hp.smoothing and the high contrast of the 511 keV image. It is expected that these artifacts will be mitigated by using finer pixel resolution.

For reference, the parameters used for the image below is:

algorithm: "RL"
parameter:
  iteration_max: 50
  acceleration:
    activate: True
    accel_factor_max: 2.0
  response_weighting:
    activate: False
    index: 0.5
  smoothing:
    activate: False
    FWHM:
      value: 2.0
      unit: "deg"
  stopping_criteria:
    statistics: "log-likelihood"
    threshold: 0.01
  background_normalization_optimization:
    activate: True
    range: {"albedo": [0.01, 10.0]}

[ ]:

iteration_idx = 49

result = image_deconvolution.results[iteration_idx]

iteration = result['iteration']
image = result['model']

data = image[:,0]
data[data <= 0 * data.unit] = 1e-12 * data.unit

hp.mollview(data, min = 1e-5, norm ='log', unit = str(data.unit), title = f'511 keV image at {iteration}th iteration', cmap = 'magma')

plt.show()