Recommended way of adding noise to sinogram data?

[liutianlin0121]

Hi!

I am wondering if there is a recommended way of adding various levels of noise to sinogram data? A naive way seems to be the following:

from dival import get_standard_dataset
import numpy as np
import matplotlib.pyplot as plt

dataset = get_standard_dataset('lodopab')
val_data = dataset.get_data_pairs('validation')

sample_at = 0
sinogram, single_gt = val_data[sample_at]
poisson_lam = 0.00001 # noise level
noisy_sinogram = (sinogram.data + np.random.poisson(poisson_lam, sinogram.shape))

fig, ax = plt.subplots(1, 2, figsize=(10, 10))
ax[0].imshow(sinogram.data)
ax[1].imshow(noisy_sinogram)

However, this naive way is perhaps not ideal: If I understood correctly, the sinogram in val_data[sample_at] has already been corrupted once during the construction of lodopab dataset, as it was mentioned that "Poisson noise corresponding to 4096 incident photons per pixel before attenuation is applied to the projection data". Thus, the superposition (sinogram.data + np.random.poisson(poisson_lam, sinogram.shape)) seems to be another round of noise adding upon the noisy data. My questions are the following:

• Is there a convenient method of adding various levels of noise on clean sinogram data?
• If the above method exists, it will also be great if I could use this method to recover the default noisy sinogram currently used in lodopab. To recover them, I am wondering what is the current default level of noise added to the lodopab CT sinogram? I've tried to track down the noise level in here but somehow couldn't find it.

Thanks a lot!! Also, big congrats for the paper published in the Journal of Imaging!

[MBaltz]

Hello @liutianlin0121 I'll try help you with the little I know.

About the method to adding noise in the sinogram, I think what you proposed is correct. Like used here: example

Parameters and others informations about LoDoPaB, I believe that what will find in this paper about the LoDoPaB: The LoDoPaB-CT Dataset

Anyway, I am also curious about the reconstruction of the noisy sinograms to their original form.

[liutianlin0121]

@MBaltz Thanks for pointing out the example!

I think my implementation above could be different from the example you pointed out. In the example you pointed out, the proj_data = ray_trafo(phantom) should be a clean sinogram, so adding noise on that clean sinogram makes sense. In my naive implementation above, the sinogram in sinogram, single_gt = val_data[sample_at] should already be a noisy signogram, and that's why I think that it is probably a bit problematic to add noise on the noisy signogram one more time.

In light of that example you pointed out, it will be great to have a similar script that exactly generates the current standard lodopab sinograms. For example, in that example, the poisson parameter is set to be 0.3 (together with other parameters). I am not sure if this is the case for the standard lodopab dataset.

[jleuschn]

Hi @liutianlin0121 , hi @MBaltz ,

you can have a look at the script used for generating the LoDoPaB dataset here. For resimulating with a different setup or noise level, this script (and maybe also this dival example) may be useful resources.

The Poisson noise in this context is not used in an additive manner. Instead we define the expected value lam (in unit "photons") for each sinogram pixel and then draw a sample from the Poisson distribution with this parameter lam, like done by np.random.poisson. One can of course reformulate it as an additive noise which is dependent on the ground truth, for this see e.g. eq. (3) of this paper you mentioned.

The Poisson noise level essentially is controlled by choosing an expected photon count. For LoDoPaB it is PHOTONS_PER_PIXEL = 4096 per detector pixel without attenuation (i.e. when the ground truth image is all zero).

Since we usually deal with log-transformed sinograms, and for LoDoPaB the values have been divided by a factor dival.util.constants.MU_MAX, we need to undo those transforms in order to obtain the expected photon counts we want to pass to np.random.poisson:

1. obs = ray_trafo(gt)  # clean sinogram (post-log values)---if the "inverse crime" should be avoided, upscale gt and use a different ray_trafo operating on the higher image resolution
2. obs_pre_log = np.exp(-MU_MAX * obs)  # expected intensity ratio I1/I0
3. obs_pre_log_noisy = np.random.poisson(obs_pre_log * PHOTONS_PER_PIXEL) / PHOTONS_PER_PIXEL

Now we can apply the log-transform again (whereby we need to deal with possible 0-photon counts, for which the logarithm would be undefined):

4. obs_pre_log_noisy_pos = np.maximum(obs_pre_log_noisy, 0.1 / PHOTONS_PER_PIXEL)
5. obs_noisy = np.log(obs_pre_log_noisy_pos) / (-MU_MAX)

I think doing this directly in numpy is the easiest option. (In theory the class ObservationGroundTruthPairDataset via GroundTruthDataset.create_pair_dataset could be used for this (like in ct_simulate_fan_beam_from_lodopab_ground_truth.py), with forward_op implementing steps 1 and 2, with noise_type='poisson' and noise_kwargs={'scaling_factor': PHOTONS_PER_PIXEL}, and with post_processor implementing steps 4 and 5.)

Hope this helps, and feel free to ask any further questions.

Thanks for your contributions!

Best regards, Johannes

[liutianlin0121]

Thanks so much!! @jleuschn

[MBaltz]

Thank you for your attention @jleuschn :)