Reconstruction problem with statistical solver: mlem

[ginsanity15]

I was trying to perform CT reconstruction with statistical solver: mlem. It works fine when the projection is noise-free. However, the reconstruction result will not converge after I add bit possion noise to the projection data. I have pasted my code below, please let me know if I have missed something here.

import numpy as np
import odl

#Reconstruction space: discretized functions on the rectangle
#[-20, 20]^2 with 200 samples per dimension.
reco_space = odl.uniform_discr(
    min_pt=[-20, -20], max_pt=[20, 20], shape=[200, 200])

#Make a parallel beam geometry with flat detector
#Angles: uniformly spaced, n = 400, min = 0, max = pi
angle_partition = odl.uniform_partition(0, 2*np.pi, 1000)

#Detector: uniformly sampled, n = 400, min = -30, max = 30
detector_partition = odl.uniform_partition(-30, 30, 400)
geometry = odl.tomo.Parallel2dGeometry(angle_partition, detector_partition)

#Create the forward operator
ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda')

#Create phantom
discr_phantom = odl.phantom.shepp_logan(reco_space, modified=True)

#Create sinogram of forward projected phantom with noise
data = ray_trafo(discr_phantom)
data += odl.phantom.white_noise(ray_trafo.range) * np.mean(data) * 0.1

#Optionally pass callback to the solver to display intermediate results
callback = (odl.solvers.CallbackPrintIteration(step=1) &
            odl.solvers.CallbackShow(step=1))

#Select iteration number
maxiter = 5

#Choose a starting point
#Start with a uniform disk in FOV
x = ray_trafo.domain.one()

#Start with reconstruction of FBP
#fbp = odl.tomo.fbp_op(ray_trafo)
#x = fbp(data)
#recon = x.asarray()
#nag_flag = recon<0
#recon[nag_flag] = 0
#x = ray_trafo.domain.element(recon)
#x.show(title='Reconstruction outcome of Filtered Backprojection')

#Run the algorithm
odl.solvers.iterative.statistical.mlem(ray_trafo, x, data, maxiter, sensitivities = 1,
                                           callback = callback,noise = 'poisson')

#Display images
discr_phantom.show(title='original image')
data.show(title='sinogram')
x.show(title='reconstructed image', force_show=True)

Edit: @adler-j typeset code

[adler-j]

The noise you add is white according to this line:

data += odl.phantom.white_noise(ray_trafo.range) * np.mean(data) * 0.1

and the MLEM solver is only defined for positive values, but this will give you negative values, causing it to diverge. Perhaps we should have a check for this in the MLEM solver.

If you actually add poisson noise according to e.g.

data = odl.phantom.poisson_noise(data)  # data is the "lambda" parameter of the noise

the solver should work.

Please also note that you typically don't want to manually set the sensitivities parameter. It is computed automatically.

[ginsanity15]

Okay I see, thank you for pointing out my problem. I just reconstructed projection with the real poisson noise and the solver works fine. Thanks for the recommendation on the sensitivities as well, I was just wondering if miracle would happen if I tune "sensitivities" a little bit.

[adler-j]

Great to hear! In general the sensitivites parameter is given since in applications they are typically measured separately. For simulation studies you should just compute them automatically.

[spcebn]

hi,do you know how to simulate the low-dose data

[kohr-h]

hi,do you know how to simulate the low-dose data

You need to check the physics to get realistic values. Here are the steps:

• For a given phantom composition (let's say water, bone and fat) get the realistic values for the attenuation coefficients, for instance from NIST. There are libraries for this, e.g., physdata.
• Define a realistic phantom size in a unit of your choice. You need to transform the attenuation coefficients according to that unit.
• Compute the ray transform of the phantom, multiplied by the source intensity (for monochromatic X-rays).
• Sample a Poisson distribution in each point, with the expected value being the result of the previous step above.
• "Low dose" means that the source intensity (expected number of photons per unit area, without object in between) is low, whatever "low" means for your application.

Here's some Python pseudocode:

intensity = 10000  # or whatever value makes sense
attenuation_coeff =  # <- code for phantom creation goes here
proj_data = ray_trafo(attenuation_coeff)
noise_free_data = intensity * np.exp(-proj_data)
noisy_data = odl.phantom.poisson_noise(noise_free_data)

A complete example that does all this (but more complex, for spectral CT) is here -- check out the "A complete example" section.

[kohr-h]

For a given phantom composition (let's say water, bone and fat) get the realistic values for the attenuation coefficients, for instance from NIST. There are libraries for this, e.g., physdata.

Of course, you can also just play around with values, e.g., take an odl.phantom.shepp_logan phantom and multiply it with some constant that gives you a workable result, like np.exp(-proj_data) is neither extremely close to 1 everywhere (too little contrast) nor very close to 0 anywhere (too dense object, "photon starvation").

[adler-j]

Just as a general pointer, the mass attenuation coefficient of water at typical x-ray energies is about 0.2 cm^2 g^-1.

As an example, this generates data with that coefficient and reconstructs it using FBP, with units carefully written out

import odl
import numpy as np

attenuation_coefficient = 0.2  # unit is cm^2/g
intensity = 50000  # number of photons per pixel

space = odl.uniform_discr([-12.8, -12.8], [12.8, 12.8], [256, 256])  # units in cm, 256 pixels

phantom = odl.phantom.shepp_logan(space)  # units is g/cm^3
attenuation = phantom * attenuation_coefficient  # Unit is cm^-1

# The ray transform respects the units given to it. 
# Since our input has units cm^-1 and we take the line integral of it, the output is unitless
ray_trafo = odl.tomo.RayTransform(space, odl.tomo.parallel_beam_geometry(space))
proj_data = ray_trafo(attenuation)

noise_free_data = intensity * np.exp(-proj_data)  # unit is photons per pixel
noisy_data = odl.phantom.poisson_noise(noise_free_data)  # unit is photons per pixel

# Convert g / cm^2
log_noisy_data = -np.log(noisy_data / intensity) / attenuation_coefficient

# FBP is "inverse" of line integral, hence result has units g / cm^3 as expected
fbp_op = odl.tomo.fbp_op(ray_trafo, filter_type='Hann')
fbp = fbp_op(log_noisy_data)

log_noisy_data.show('log_noisy_data')
fbp.show('fbp window', clim=[0.95, 1.15])