Poisson noise for acquisition data

[epapoutsellis]

There is a problem when using poisson noise method from TestData and acquisitionData. Probably is some scaling issue using np.random.poisson

from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData

import numpy as np 
import numpy                          
import matplotlib.pyplot as plt

from ccpi.optimisation.algorithms import PDHG

from ccpi.optimisation.operators import BlockOperator, Gradient
from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
                      MixedL21Norm, BlockFunction, KullbackLeibler, IndicatorBox

from ccpi.astra.operators import AstraProjectorSimple

import os, sys
import tomophantom
from tomophantom import TomoP2D

model = 1 # select a model number from the library
N = 128 # set dimension of the phantom
path = os.path.dirname(tomophantom.__file__)
path_library2D = os.path.join(path, "Phantom2DLibrary.dat")

phantom_2D = TomoP2D.Model(model, N, path_library2D)    
data = ImageData(phantom_2D)
ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N)

# Create acquisition data and geometry
detectors = N
angles = np.linspace(0, np.pi, 180)
ag = AcquisitionGeometry('parallel','2D',angles, detectors)

Aop = AstraProjectorSimple(ig, ag, device = 'cpu')
sin = Aop.direct(data)

scale = 5
noisy_data1 = AcquisitionData(np.random.poisson( scale * (sin.as_array()))/scale)

from ccpi.framework import TestData

noisy_data2 = TestData.random_noise( sin.as_array()/scale, mode = 'poisson', seed = 5)*scale

plt.imshow(noisy_data1.as_array())
plt.show()

plt.imshow(noisy_data2)
plt.show()

[gfardell]

I think we need to look closely at how we add noise. The TestData implementation is from this package https://kite.com/python/docs/skimage.util.random_noise

I think the input data needs to be between 0-1 (or maybe -1 to 1) as it clips the output for poisson noise. So in your case you'll need to divide the input by np.amax(sin.as_array()) to normalise the data.

But in your example we want to add noise to the intensity and not the logged intensity so it needs a bit more work.

I guess something like:

scale = np.amax(sin.as_array())
exp_data = np.exp(-sin.as_array()/scale)
gain = 400
noise = TestData.random_noise(exp_data/gain, mode = 'poisson')*gain
data_noisy = -scale * np.log(noise)

where the gain determines the noise level.

[gfardell]

I think this is something we want to do in he framework though rather than every time we want to add noise.

[paskino]

I suppose that if we want to do something like @gfardell suggests we'll have to add another staticmethod in TestData that will exp add noise and log.

BTW as of #381 DataContainers do have exp and log unary operations.

[evelinaametova]

I would suggest we distinguish between adding noise to a random image like peppers or cameraman and adding noise to CT sinogram.

for a random image: When adding poisson noise, intensity in each pixel is interpreted as mean of poisson distribution, therefore if image intensity is distributed between 0 and 1, it has to be scaled by some number. In Matlab they overcome this issue by scaling float images by 1e12 for double and 1e6 for single. I would suggest to have the same default numbers, but also allow user to choose the scale.

Poisson distribution cannot have negative mean, image has to be scaled/ clipped. For an image with intensity (min,max), I would suggest to normalise it between (0,1), then scale, add noise, scale back to (0,1), and then back to original range (min, max), all values below min and above max, has to be clipped.

for sinogram For sinogram I would suggest very similar to Gemma's implementation but add a bit of a physical meaning to gain and sin. In general, CT detectors outputs unit values. Poisson noise level depends on it's mean: smaller mean -> higher noise, larger mean -> less pronounced noise. Meaning that noise level for uint8 detector and uint16 detector will be considerably different.

N = 512
n_angles = 360

# initialise loader
loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N, N))
# get ImageGeometry
ig = data.geometry
ig.voxel_size_x = 0.005
ig.voxel_size_y = ig.voxel_size_x

# Create Acquisition data
angles = np.linspace(0, np.pi, n_angles, dtype=np.float32)

ag = AcquisitionGeometry(geom_type = 'parallel',
                         dimension = '2D', 
                         angles = angles,
                         pixel_num_h = N,
                         pixel_size_h = ig.voxel_size_x)

dev = 'gpu'

Aop = AstraProjectorSimple(ig, ag, dev)    
sino = Aop.direct(data).as_array()

n_bits = 8
I_0 = 2 ** n_bits - 1

sino_tmp = np.uint32(np.round(I_0 * np.exp(-sino)))
# here we treat each pixel as an independent implementation of poisson process with a mean given with pixel intensity
sino_noisy = -np.log(np.float32(np.random.poisson(sino_tmp )) / I_0)

plt.imshow(sino_noisy )
plt.show()

note, that in this implementation intensity values in data also have physical meaning, i.e. attenuation coefficient. If we keep pixel/ voxel size equal to 1, then there will be simply not enough penetration, i.e. X-ray penetration length will be too long for given attenuation values. Alternatively, data intensity has to be scaled to have enough X-ray penetration.

[epapoutsellis]

@gfardell , @evelinaametova Is it possible to create a utility for Poisson noise for tomodata, so we all using the same code? We can add to CIL-demos or directly to CCPi-Framework?

[paskino]

@evelinaametova @gfardell can you see if this can be closed given the current implementation?

[jakobsj]

@evelinaametova We discussed at tcon, is this something that can be quickly added, in time for the Pubs deadline?

[gfardell]

@evelinaametova is this used in your notebooks, otherwise I think we should bump it to the next milestone (again)

[MargaretDuff]

Had a chat with @paskino. We thought about creating a demo notebook discussing poisson noise and how to add it correctly. We would utilise some of the code by @evelinaametova in this issue.