Can we use MLEM solver for Cone Beam CT?

[ginsanity15]

I was trying to perform reconstruction of the Shepp-Logan Phantom with MLEM solver and I noticed the reconstructed Shepp-Logan Phantom doesn't have homogeneous intensity on its outer boundary. For me it seems that part of the boundary will be much brighter than the other, while precisely which part of the boundary will be brighter depend on the setting of "axis" in geometry of the tomography (In my case, it is the setting of "axis" in "CircularConeFlatGeometry"). To make my idea clear, I have attached my code as well as two reconstructed image with different different axes below.

I have found paper: "EM reconstruction algorithms for emission and transmission tomography" published by Lange and Carson, which has given different update equation for emission and transmission tomography and it seems the equation implemented in ODL is only suitable for emission tomography. Could that be the cause of this inhomogeneous boundary?

import numpy as np
import odl

#--- Set up the forward operator (ray transform) --- #

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

#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, 500)

#Detector: uniformly sampled, n = 400, min = -30, max = 30
detector_partition = odl.uniform_partition([-30, -30], [30, 30], [400, 400])
geometry = odl.tomo.CircularConeFlatGeometry(
    angle_partition, detector_partition, src_radius=140, det_radius=40,
    axis=(1, 1, 1))

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

#--- Generate artificial data --- #

#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.poisson_noise(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 = 15

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

#Create norm volume
norm_proj = ray_trafo.range.one()
#norm_volume = fbp(norm_proj)
norm_volume = ray_trafo.adjoint(norm_proj)

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

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

Edit: @adler-j improved typesetting

[ginsanity15]

Since no one had answered this question, I will give some of my thoughts here.

After reading couple of papers related to apply statistical iterative algorithms in transmission computed tomography, I found the update equation for transmission CT always requires knowledge of the number of photons that would be detected in the absence of phantom (for reference, see section II, part A, B, C of the paper in this link: https://web.eecs.umich.edu/~fessler/papers/files/jour/95/pre/tip,lange.pdf), which hasn't been included in ODL, because when reconstructing with functions in ODL, we are actually working on the log-converted sinogram (sum of attenuation coefficients along each directions), instead of the raw sinogram, which is the photon counting along different direction.

Another way to address this missed information is: we have assumed the above mentioned variable is uniform. However, in reality the detected photon in the absence of phantom is not homogeneous. More precisely, it is higher in the center of detector (given the assumption that the perpendicular line across the x-ray source hit the detector at the central point) and lower in the edges of detector, and that is why when working with CT scanner in reality, we will not only record projection image of the sample, but also flat field image (projection image without the sample). As a result of this ill modelling, we are actually reconstructing image from projections with lower than reality count in the center and higher than reality counts in the edges. That's why the reconstruction from those projections will have a higher than reality intensity along the direction of the incident x-ray, which means to a higher attenuation coefficients.

What do you think of this explanation? Have I missed some functions in ODL? Please let me know if I haven't made my point clear.

[ozanoktem]

The way ML-EM iterates are formulated have noting to do with whether the forward model is given by the ray transform (as in CT) or some other variant as in SPECT. It is related to the model for the statistical properties of data. Each statistical model yields a different data likelihood, which in turn yields a different ML-EM iterative scheme. The optional parameter noise for the mlem solver is used to specify the statistical model. If the available options does not cover your case, you need to work out the the loglikelihood-function for your case and then add it to the mlem solver with an appropriate value for noise .

[ginsanity15]

Well that is precisely what I have done, I implemented update equations other than the one currently given in mlem solver. And my real concern is that ODL alone is not enough to perform statistical reconstruction of transmission computed tomography, we need to import data outside ODL as well. The comment above was trying to proof why ODL isn't enough. Anyway, I will explore more in ODL and it is very likely that this library is more powerful than I thought.

By the way, does the current implementation of loglikelihood in ODL works? I noticed that you are trying to perform multiplication between op.domain element and op.range element. Is that possible for us to do so?

import odl
import numpy as np

reco_space = odl.uniform_discr(
    min_pt=[-20, -20, -20], max_pt=[20, 20, 20], shape=[256, 256, 256])
angle_partition = odl.uniform_partition(0, 2*np.pi, 500)
detector_partition = odl.uniform_partition([-30, -30], [30, 30], [400, 400])

geometry = odl.tomo.CircularConeFlatGeometry(
        angle_partition, detector_partition, src_radius=140, det_radius=40,
    axis=(1, 1, 1))

ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda')

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

data = ray_trafo(discr_phantom)

x = ray_trafo.domain.one()

output = odl.solvers.iterative.statistical.loglikelihood(x, data)

print(output)

I was testing with the code above and it returns the error as follow: Traceback (most recent call last):

  File "<ipython-input-52-45e744791cc0>", line 1, in <module>
    odl.solvers.iterative.statistical.loglikelihood(x, data)

  File "C:\Users\bobgao\AppData\Local\Continuum\anaconda2\lib\site-packages\odl\solvers\iterative\statistical.py", line 234, in loglikelihood
    return np.sum(data * np.log(x + 1e-8) - x)

TypeError: unsupported operand type(s) for *: 'DiscreteLpElement' and 'DiscreteLpElement'

edit: @adler-j typeset

[adler-j]

[original question]

This is basically because the MLEM algorithm has a slight positive bias (simply because it can only give positive answers, even if the ground truth is zero.) However, as the number of iterations go to infinity and noise goes to zero, this should be reduced. If you try e.g. 100 or 1000 iterations, you'll find that this problem mostly goes away.

After reading couple of papers related to apply statistical iterative algorithms in transmission computed tomography, I found the update equation for transmission CT always requires knowledge of the number of photons that would be detected in the absence of phantom

This is correct!

which hasn't been included in ODL, because when reconstructing with functions in ODL, we are actually working on the log-converted sinogram (sum of attenuation coefficients along each directions), instead of the raw sinogram, which is the photon counting along different direction.

Working with the log-converted sinogram is a choice you make as a user. If you want a physically correct forward model for CT imaging, you can achieve this by doing e.g.

attenuation_coefficient = ?
exp = odl.ufunc_ops.exp(ray_trafo.range)
nonlinear_ray_trafo = exp * (-attenuation_coefficient * ray_trafo)

with this model, you would get the correct number of counts.

What do you think of this explanation? Have I missed some functions in ODL? Please let me know if I haven't made my point clear.

This explanation, while very valid for actual data, fails in your case since you generate the data yourself and as such anything related to physical reality (e.g. nonlinearity of the forward model, etc) cannot impact your result. The true explanation for this behaviour in this case should be the one I gave above.

By the way, does the current implementation of loglikelihood in ODL works? I noticed that you are trying to perform multiplication between op.domain element and op.range element. Is that possible for us to do so?

As you correctly observe, this does (and should) not work. Instead you should do e.g.

odl.solvers.iterative.statistical.loglikelihood(ray_trafo(x), data)

[ginsanity15]

I just tried to run the first code with 150 iterations and the problem of inhomogeneous boundary intensity has indeed gone away. Thanks for your help as well as the discussion of ML-EM.

[adler-j]

I'll close this issue then in order to keep discussions coherent. If you have further questions please open a new issue! We'll try to be faster in answering next time :)

[ginsanity15]

Sure, I will contact you if I get lost in ODL again. I will also close the issue now.