FBP_slice_example.py

# 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)]')