ODL projection example

import odl
import numpy as np

# Discrete reconstruction space: discretized functions on the rectangle
# [-20, 20]^2 with 300 samples per dimension.
reco_space = odl.uniform_discr(
    min_corner=[-1, -1], max_corner=[1, 1], nsamples=[256, 256],
    dtype='float32')

# Make a parallel beam geometry with flat detector

# Angles: uniformly spaced, n = 180, min = 0, max = pi
angle_partition = odl.uniform_partition(0, np.pi, 180)

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

# Create a discrete Shepp-Logan phantom (modified version)
phantom = odl.util.phantom.shepp_logan(reco_space, modified=True)

# ray transform aka forward projection. We use ASTRA CUDA/CPU backend.
ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda')
ray_trafo_cpu = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cpu')

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

# 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(proj_data)
backproj_cpu = ray_trafo_cpu.adjoint(proj_data_cpu)

# Shows a slice of the phantom, projections, and reconstruction
phantom.show(title='Phantom')
proj_data.show(title='cuda sinogram')
backproj.show(title='cuda backprojection')
proj_data_cpu.show(title='cpu sinogram')
backproj_cpu.show(title='cpu backprojection')
(proj_data_cpu - proj_data).show(title='diff sinogram')
(backproj_cpu - backproj).show(title='diff backprojection')