Skip to content

Instantly share code, notes, and snippets.

@adler-j
Created February 22, 2017 10:17
Show Gist options
  • Save adler-j/948196c0f95222df86cabe56934a6232 to your computer and use it in GitHub Desktop.
Save adler-j/948196c0f95222df86cabe56934a6232 to your computer and use it in GitHub Desktop.
import numpy as np
import odl
# --- Set-up geometry of the problem --- #
# Discrete reconstruction space: discretized functions on the rectangle
# [-20, 20]^2 with 300 samples per dimension.
reco_space = odl.uniform_discr(
min_pt=[-20, -20], max_pt=[20, 20], shape=[300, 300], dtype='float32')
# Angles: uniformly spaced, n = 1000, min = 0, max = pi
angle_partition = odl.uniform_partition(0, np.pi, 1000)
# Detector: uniformly sampled, n = 500, min = -30, max = 30
detector_partition = odl.uniform_partition(-30, 30, 500)
# Make a parallel beam geometry with flat detector
geometry = odl.tomo.Parallel2dGeometry(angle_partition, detector_partition)
# --- Create Filtered Back-Projection (FBP) operator --- #
# Ray transform (= forward projection).
ray_trafo = odl.tomo.RayTransform(reco_space, geometry)
# Fourier transform in detector direction
fourier = odl.trafos.FourierTransform(ray_trafo.range, axes=[1])
# Create ramp in the detector direction
ramp_function = fourier.range.element(lambda x: np.abs(x[1]) / (2 * np.pi))
# Create ramp filter via the convolution formula with fourier transforms
ramp_filter = fourier.inverse * ramp_function * fourier
# Create filtered back-projection by composing the back-projection (adjoint)
# with the ramp filter.
fbp = ray_trafo.adjoint * ramp_filter
# --- Show some examples --- #
# Create a discrete Shepp-Logan phantom (modified version)
phantom = odl.phantom.shepp_logan(reco_space, modified=True)
# Create projection data by calling the ray transform on the phantom
proj_data = ray_trafo(phantom)
# Calculate filtered back-projection of data
fbp_reconstruction = fbp(proj_data)
# Shows a slice of the phantom, projections, and reconstruction
fig = phantom.show(coords=[None, 0])
fig = fbp_reconstruction.show(coords=[None, 0], fig=fig)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment