Skip to content

Instantly share code, notes, and snippets.

@adler-j
Last active December 21, 2016 11:10
Show Gist options
  • Save adler-j/b3c6ee6a9fb917b9c31c5c6ec3532977 to your computer and use it in GitHub Desktop.
Save adler-j/b3c6ee6a9fb917b9c31c5c6ec3532977 to your computer and use it in GitHub Desktop.
# Copyright 2014-2016 The ODL development group
#
# This file is part of ODL.
#
# ODL is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# ODL is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ODL. If not, see <http://www.gnu.org/licenses/>.
"""
Example using a filtered back-projection (FBP) in cone-beam 3d using `fbp_op`.
Note that the FBP is only approximate in this geometry, but still gives a
decent reconstruction that can be used as an initial guess in more complicated
methods.
"""
import numpy as np
import odl
from odl.tomo.util.utility import axis_rotation_matrix
# --- Set-up geometry of the problem --- #
full_reco_space = odl.uniform_discr(
min_pt=[-20, -20, -20], max_pt=[20, 20, 20], shape=[500, 500, 500],
dtype='float32')
angle_partition = odl.uniform_partition(0, 2 * np.pi, 360)
detector_partition = odl.uniform_partition([-40, -40], [40, 40], [500, 500])
full_geometry = odl.tomo.CircularConeFlatGeometry(
angle_partition, detector_partition, src_radius=40, det_radius=40,
axis=[0, 0, 1])
# --- Create ray transform and filtering operator --- #
full_ray_trafo = odl.tomo.RayTransform(full_reco_space, full_geometry,
impl='astra_cuda')
# %% --- Create raw data --- #
# Create a discrete Shepp-Logan phantom (modified version)
phantom = odl.phantom.shepp_logan(full_reco_space, modified=True)
phantom.show('phantom')
# Create projection data by calling the ray transform on the phantom
proj_data = full_ray_trafo(phantom)
# %% --- Reconstruct single slice --- #
slice_normal = np.array([1, 1, 0], dtype=float)
slice_shift = np.array([0, 0, 5 * np.sqrt(2)], dtype=float)
slice_reco_space = odl.uniform_discr(
min_pt=[-20 - slice_shift[0], -20 - slice_shift[1], -full_reco_space.cell_sides[2] / 2 - slice_shift[2]],
max_pt=[20 - slice_shift[0], 20 - slice_shift[1], full_reco_space.cell_sides[2] / 2 - slice_shift[2]],
shape=[500, 500, 1],
dtype='float32')
rot_geometry = odl.tomo.CircularConeFlatGeometry(
angle_partition, detector_partition, src_radius=40, det_radius=40,
axis=slice_normal)
slice_ray_trafo = odl.tomo.RayTransform(slice_reco_space, rot_geometry,
impl='astra_cuda')
slice_fbp_op = odl.tomo.fbp_op(slice_ray_trafo)
fbp_reconstruction = slice_fbp_op(proj_data)
fbp_reconstruction.show(title='FBP, slice normal = [1, 1, 0], slice shift = [0, 0, 5 * sqrt(2)]')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment