Add Poisson noise

AddNoise illustrates how one can add pre-log Poisson noise to RTK-simulated projections. This simulation has been used in previous articles, e.g.[Rit et al, Med Phys, 2016] and described as The same simulations were repeated with Poisson noise. The Shepp Logan densities were weighted by 0.01879 mm\(^{−1}\), i.e., the linear attenuation coefficient of water at 75 keV. The number of photons received per detector pixel without object in the beam was constant (…) and equal to \(10^4\).

#include <itkImageFileWriter.h>
#include <itkMultiplyImageFilter.h>
#include <itkExpImageFilter.h>
#include <itkShotNoiseImageFilter.h>
#include <itkThresholdImageFilter.h>
#include <itkLogImageFilter.h>

#include <rtkConstantImageSource.h>
#include <rtkSheppLoganPhantomFilter.h>

using ImageType = itk::Image<float, 3>;

// Constant parameters
const float I0 = 1e4;     // Number of photons before attenuation
const float mu = 0.01879; // mm^-1, attenuation coefficient of water at 75 keV

int
main()
{
  // Simulate a Shepp Logan projection
  auto geometry = rtk::ThreeDCircularProjectionGeometry::New();
  geometry->AddProjection(1000., 0., 0.);

  auto source = rtk::ConstantImageSource<ImageType>::New();
  source->SetSize(itk::MakeSize(64, 64, 1));
  source->SetSpacing(itk::MakeVector(2., 2., 2.));
  source->SetOrigin(itk::MakePoint(-63., -63., 0.));

  auto shepp = rtk::SheppLoganPhantomFilter<ImageType, ImageType>::New();
  shepp->SetInput(source->GetOutput());
  shepp->SetGeometry(geometry);
  shepp->SetPhantomScale(70.);

  // Use ITK to add pre-log Poisson noise
  auto multiply = itk::MultiplyImageFilter<ImageType>::New();
  multiply->SetInput(shepp->GetOutput());
  multiply->SetConstant(-mu);

  auto exp = itk::ExpImageFilter<ImageType, ImageType>::New();
  exp->SetInput(multiply->GetOutput());

  auto multiply2 = itk::MultiplyImageFilter<ImageType>::New();
  multiply2->SetInput(exp->GetOutput());
  multiply2->SetConstant(I0);

  auto poisson = itk::ShotNoiseImageFilter<ImageType>::New();
  poisson->SetInput(multiply2->GetOutput());

  auto threshold = itk::ThresholdImageFilter<ImageType>::New();
  threshold->SetInput(poisson->GetOutput());
  threshold->SetLower(1.);
  threshold->SetOutsideValue(1.);

  auto multiply3 = itk::MultiplyImageFilter<ImageType>::New();
  multiply3->SetInput(threshold->GetOutput());
  multiply3->SetConstant(1. / I0);

  auto log = itk::LogImageFilter<ImageType, ImageType>::New();
  log->SetInput(multiply3->GetOutput());

  auto multiply4 = itk::MultiplyImageFilter<ImageType>::New();
  multiply4->SetInput(log->GetOutput());
  multiply4->SetConstant(-1. / mu);

  itk::WriteImage(multiply4->GetOutput(), "projection.mha");

  return 0;
}

import itk
from itk import RTK as rtk

# Constant values of the script
I0 = 1e4  # Number of photons before attenuation
mu = 0.01879  # mm^-1 (attenuation coefficient of water at 75 keV)

# Simulate a Shepp Logan projection
ImageType = itk.Image[itk.F, 3]
geometry = rtk.ThreeDCircularProjectionGeometry.New()
geometry.AddProjection(1000.0, 0.0, 0)
projection = rtk.constant_image_source(
    size=[64, 64, 1], spacing=[2.0] * 3, ttype=[ImageType], origin=[-63.0, -63.0, 0]
)
projection = rtk.shepp_logan_phantom_filter(
    projection, geometry=geometry, phantom_scale=70
)
itk.imwrite(projection, "projection.mha")

# Use ITK to add pre-log Poisson noise
noisy_projection = itk.multiply_image_filter(projection, constant=-mu)
noisy_projection = itk.exp_image_filter(noisy_projection)
noisy_projection = itk.multiply_image_filter(noisy_projection, constant=I0)
noisy_projection = itk.shot_noise_image_filter(noisy_projection)
noisy_projection = itk.threshold_image_filter(
    noisy_projection, lower=1.0, outside_value=1.0
)
noisy_projection = itk.multiply_image_filter(noisy_projection, constant=1.0 / I0)
noisy_projection = itk.log_image_filter(noisy_projection)
noisy_projection = itk.multiply_image_filter(noisy_projection, constant=-1.0 / mu)
itk.imwrite(noisy_projection, "noisy_projection.mha")

# Alternative NumPy implementation
# projection = rtk.constant_image_source(size=[64,64,1],spacing=[2.]*3, ttype=[ImageType], origin=[-63.,-63.,0])
# projection = rtk.shepp_logan_phantom_filter(projection, geometry=geometry, phantom_scale=70)
# import numpy as np
# proj_array = itk.array_view_from_image(projection)
# proj_array = I0*np.exp(-1.*mu*proj_array)
# proj_array = np.maximum(np.random.poisson(proj_array), 1)
# proj_array = np.log(I0/proj_array)/mu
# projection_noisy = itk.image_view_from_array(proj_array)
# projection_noisy.CopyInformation(projection)
# itk.imwrite(projection_noisy, 'projection.mha')

import matplotlib.pyplot as plt

plt.figure(figsize=(14, 4))
plt.subplot(131)
plt.title("Noiseless projection")
plt.imshow(
    itk.array_view_from_image(projection)[0, :, :],
    cmap=plt.cm.gray,
    origin="lower",
    clim=[0, 140],
)
plt.xlabel("u pixel index")
plt.ylabel("v pixel index")
plt.subplot(132)
plt.title("Noisy projection")
plt.imshow(
    itk.array_view_from_image(noisy_projection)[0, :, :],
    cmap=plt.cm.gray,
    origin="lower",
    clim=[0, 140],
)
plt.xlabel("u pixel index")
plt.subplot(133)
plt.title("Central vertical profile")
plt.plot(itk.array_view_from_image(projection)[0, :, 32], label="Noiseless")
plt.plot(itk.array_view_from_image(noisy_projection)[0, :, 32], label="Noisy")
plt.xlabel("v pixel index")
plt.ylabel("Pixel intensity")
plt.legend()
plt.savefig("AddNoise.png", bbox_inches="tight")
plt

The plot resulting from the Python version is

AddNoise