How to add Poisson noise to projection data？

[Hznnn]

dataset = dicom.read_file('./L067_FD_1_1.CT.0001.0001.2015.12.22.18.09.40.840353.358074219.IMA')
img = dataset.pixel_array.astype(np.float32)
RescaleSlope = dataset.RescaleSlope
RescaleIntercept = dataset.RescaleIntercept
CT_img = img * RescaleSlope + RescaleIntercept
reco_space = odl.uniform_discr(min_pt=[-128, -128], max_pt=[128, 128], shape=[512, 512], dtype='float32')
angle_partition = odl.uniform_partition(0, 2 * np.pi, 1000)
detector_partition = odl.uniform_partition(-360, 360, 1000)
geometry = odl.tomo.FanBeamGeometry(angle_partition, detector_partition,
                                            src_radius=500, det_radius=500)
ray_trafo = odl.tomo.RayTransform(reco_space, geometry)
fbp = odl.tomo.fbp_op(ray_trafo, filter_type='Hann', frequency_scaling=0.8)
photons_per_pixel = 1e5
mu_water = 0.02
epsilon = 0.0001
nonlinear_operator = odl.ufunc_ops.exp(ray_trafo.range) * (- mu_water * ray_trafo)
phantom = img/1000.0
proj_data = nonlinear_operator(phantom)
proj_data = odl.phantom.poisson_noise(proj_data * photons_per_pixel) / photons_per_pixel
proj_data = -np.log(epsilon + proj_data) / mu_water
degrade = fbp(proj_data)

The result degrade looks lower than 1e5. I am looking forward to hear some better suggestions from you.

[adler-j]

What is the scale of proj_data? It needs to be at scaled to [0, 1] for the above code to work.

[Hznnn]

proj_data in proj_data = nonlinear_operator(phantom)?

[adler-j]

Yes indeed.

[adler-j]

Another option is that you reproduce this using e.g. the shepp-logan phantom, that way others can test your code and much more easily help.

[Hznnn]

reco_space = odl.uniform_discr(min_pt=[-128, -128], max_pt=[128, 128], shape=[512, 512], dtype='float32')
angle_partition = odl.uniform_partition(0, 2 * np.pi, 1000)
detector_partition = odl.uniform_partition(-360, 360, 1000)
geometry = odl.tomo.FanBeamGeometry(angle_partition, detector_partition,
                                            src_radius=500, det_radius=500)
ray_trafo = odl.tomo.RayTransform(reco_space, geometry)

How to set min_pt and max_pt inreco_space, detector_partition, even if I read the documentation, it is difficult to understand these parameters.

[adler-j]

Try using the more simple odl.tomo.cone_beam_geometry constructor which has fewer options.

[Hznnn]

I want to add noise to the mayo data to get simulated noise data. The configuration of this process in IRT Library in http://web.eecs.umich.edu/~fessler/code/ is

param.NProj = 720;
param.DSD   = 1270;
param.DSO   = 870;
param.nu    = 848;
param.du    = 0.6;

%%  image params
param.nx    = 512;
param.ny    = 512; 
param.dx    = 0.5;

ig = image_geom('nx', param.nx, 'ny', param.ny, 'dx', param.dx);
sg = sino_geom('fan', 'nb', param.nu, 'na',param.NProj, 'ds',param.du,'dsd',param.DSD,'dso',param.DSO ,'units','cm','strip_width','d');
systemMatrix = Gtomo_nufft_new(sg, ig);
tmp = fbp2(sg, ig);

[Hznnn]

How to make the same settings with odl？

[adler-j]

Sorry for the late answer. I'm not very well versed with fesslers implementation. Which of the parameters indicate the photon fluence?

[sunchang2017]

Where is the document that explain the mechanism of adding Poisson Noise? I find it difficult to read documents, too... I also curious about how to make the same settings in IRT library with odl without considering the photon fluence.

[ozanoktem]

Where is the document that explain the mechanism of adding Poisson Noise? I find it difficult to read documents, too... I also curious about how to make the same settings in IRT library with odl without considering the photon fluence.

Check out https://odlgroup.github.io/odl/odl.phantom.noise.html for documentation on various ODL functions to create noise samples of different distributions, including Poisson.

[Binganyuan]

o_spa

cool!