boolean_radon_example.py

"""Example of using the boolean ray transform with ODL.

This example is with the 2d ray transform, but you could easily change the
geometry to 3d and get the example you gave. See

github.com/odlgroup/odl/blob/master/examples/tomo/ray_trafo_parallel_3d.py

for an example.

Note that some rounding errors may occur.
"""

import numpy as np
import odl

# Select the type of phantom to use
phantom_type = 'torus'

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

# Create phantom
# Note that to create a phantom from e.g. a numpy array, you would do
# phantom = space.element(my_numpy_array)
if phantom_type == 'shepp_logan':
    # Create a discrete Shepp-Logan phantom (modified version)
    phantom = odl.phantom.shepp_logan(space, modified=True)
elif phantom_type == 'touching_squares':
    phantom = odl.phantom.cuboid(space, min_pt=[-10, -10], max_pt=[-0, -0])
    phantom += odl.phantom.cuboid(space, min_pt=[0, 0], max_pt=[10, 10])
elif phantom_type == 'separate_squares':
    phantom = odl.phantom.cuboid(space, min_pt=[-10, -10], max_pt=[-1, -1])
    phantom += odl.phantom.cuboid(space, min_pt=[1, 1], max_pt=[10, 10])
elif phantom_type == 'torus':
    ellipses = [[1, 0.3, 0.3, -0.5, 0, 0],
                [1, 0.3, 0.3, 0.5, 0, 0]]
    phantom = odl.phantom.ellipse_phantom(space, ellipses)
else:
    raise RuntimeError('unknown phantom_type')

# Make a parallel beam geometry with flat detector
# Angles: uniformly spaced, n = 700, min = 0, max = pi
angle_partition = odl.uniform_partition(0, np.pi, 700)
# Detector: uniformly sampled, n = 700, min = -28, max = 28
detector_partition = odl.uniform_partition(-28, 28, 700)
geometry = odl.tomo.Parallel2dGeometry(angle_partition, detector_partition)

# Ray transform (= forward projection). We use scikit backend
# downloadable by "pip install scikit-image"
# The solver becomes MUCH faster with "impl='astra_cuda'" but this requires
# ASTRA to be installed.
ray_trafo = odl.tomo.RayTransform(space, geometry, impl='scikit')

# Create projection data by calling the ray transform on the phantom
proj_data = ray_trafo(phantom)

# Create boolean projection data by taking positive values
boolean_proj_data = np.greater(proj_data, 0)

# Back-projection can be done by simply calling the adjoint operator on the
# projection data (or any element in the projection space).
backproj = ray_trafo.adjoint(boolean_proj_data)

# Inverse is the data that is in all projections
# Here we subtract a small number to account for numerics.
inverse_result = np.greater_equal(backproj, angle_partition.extent() - 0.0001)

# Shows a slice of the phantom, projections, and reconstruction
phantom.show(title='Phantom')
proj_data.show(title='Projection data (sinogram)')
boolean_proj_data.show(title='Boolean projection data')
backproj.show(title='Back-projected data')
inverse_result.show(title='Inverse')
(phantom - inverse_result).show('Error')