From d4eb486e871a5332087c6ea124d3b69e9390d8cf Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 28 Jun 2016 11:18:57 +0200 Subject: [PATCH] API: move phantoms to odl.phantoms, closes #463 Also did the following fixed to phantoms * Optimized 2d ellipse phantoms in same way as 3d * Fixed bug with submarined in non centered geometries --- doc/source/release_notes.rst | 20 + examples/solvers/chambolle_pock_deconvolve.py | 2 +- examples/solvers/chambolle_pock_tomography.py | 4 +- .../solvers/conjugate_gradient_tomography.py | 4 +- .../filtered_backprojection_parallel_2d.py | 2 +- examples/tomo/ray_trafo_circular_cone.py | 2 +- examples/tomo/ray_trafo_fanflat.py | 2 +- examples/tomo/ray_trafo_helical.py | 2 +- examples/tomo/ray_trafo_parallel_2d.py | 2 +- examples/tomo/ray_trafo_parallel_3d.py | 2 +- examples/tomo/scikit_ray_trafo_parallel_2d.py | 2 +- examples/tomo/stir_project.py | 2 +- examples/tomo/stir_reconstruct.py | 2 +- examples/trafos/fourier_trafo.py | 4 +- examples/visualization/show_2d.py | 2 +- examples/visualization/show_2d_complex.py | 2 +- examples/visualization/show_productspace.py | 6 +- examples/visualization/show_update_2d.py | 2 +- odl/__init__.py | 5 +- odl/phantom/__init__.py | 39 + odl/phantom/emission.py | 151 ++++ odl/phantom/geometric.py | 194 +++++ odl/phantom/misc_phantoms.py | 169 +++++ odl/phantom/noise.py | 49 ++ odl/phantom/phantom_utils.py | 327 ++++++++ odl/phantom/transmission.py | 120 +++ odl/util/__init__.py | 3 - odl/util/phantom.py | 715 ------------------ setup.py | 2 +- test/discr/diff_ops_test.py | 7 +- .../tomo/ray_transform_slow_test.py | 4 +- test/largescale/tomo/stir_slow_test.py | 2 +- test/tomo/backends/astra_cpu_test.py | 4 +- test/tomo/backends/astra_cuda_test.py | 10 +- test/tomo/backends/scikit_test.py | 2 +- test/tomo/operators/ray_trafo_test.py | 6 +- 36 files changed, 1112 insertions(+), 761 deletions(-) create mode 100644 odl/phantom/__init__.py create mode 100644 odl/phantom/emission.py create mode 100644 odl/phantom/geometric.py create mode 100644 odl/phantom/misc_phantoms.py create mode 100644 odl/phantom/noise.py create mode 100644 odl/phantom/phantom_utils.py create mode 100644 odl/phantom/transmission.py delete mode 100644 odl/util/phantom.py diff --git a/doc/source/release_notes.rst b/doc/source/release_notes.rst index 584d072a1c4..6dc6d865e56 100644 --- a/doc/source/release_notes.rst +++ b/doc/source/release_notes.rst @@ -6,6 +6,26 @@ Release Notes ############# +ODL 0.2.4 Release Notes (2016-06-28) +==================================== + +New features +------------ +- Add ``uniform_discr_fromdiscr`` (`PR 467`). +- Add conda build files (`commit 86ff166`). + +Bugfixes +-------- +- Fix bug in submarine phantom with non-centered space (`PR 469`). +- Fix crash when plotting in 1d (`commit 3255fa3`). + +Changes +------- +- Move phantoms to new module odl.phantom (`PR 469`). +- Rename ``RectPartition.is_uniform`` to ``RectPartition.is_uniform`` + (`PR 468`). + + ODL 0.2.3 Release Notes (2016-06-12) ==================================== diff --git a/examples/solvers/chambolle_pock_deconvolve.py b/examples/solvers/chambolle_pock_deconvolve.py index dafdfb4e2ca..9f5be53e78f 100644 --- a/examples/solvers/chambolle_pock_deconvolve.py +++ b/examples/solvers/chambolle_pock_deconvolve.py @@ -57,7 +57,7 @@ # odl.diagnostics.OperatorTest(conv_op).run_tests() # Create phantom -phantom = odl.util.shepp_logan(space, modified=True) +phantom = odl.phantom.shepp_logan(space, modified=True) # Create vector of convolved phantom data = convolution(phantom) diff --git a/examples/solvers/chambolle_pock_tomography.py b/examples/solvers/chambolle_pock_tomography.py index 051e1edcecc..af00cd1261f 100644 --- a/examples/solvers/chambolle_pock_tomography.py +++ b/examples/solvers/chambolle_pock_tomography.py @@ -64,11 +64,11 @@ # Create phantom -discr_phantom = odl.util.shepp_logan(reco_space, modified=True) +discr_phantom = odl.phantom.shepp_logan(reco_space, modified=True) # Create sinogram of forward projected phantom with noise data = ray_trafo(discr_phantom) -data += odl.util.white_noise(ray_trafo.range) * np.mean(data) * 0.1 +data += odl.phantom.white_noise(ray_trafo.range) * np.mean(data) * 0.1 # --- Set up the inverse problem --- # diff --git a/examples/solvers/conjugate_gradient_tomography.py b/examples/solvers/conjugate_gradient_tomography.py index 9566d49fb74..342fe481385 100644 --- a/examples/solvers/conjugate_gradient_tomography.py +++ b/examples/solvers/conjugate_gradient_tomography.py @@ -62,11 +62,11 @@ # Create phantom -discr_phantom = odl.util.shepp_logan(reco_space, modified=True) +discr_phantom = odl.phantom.shepp_logan(reco_space, modified=True) # Create sinogram of forward projected phantom with noise data = ray_trafo(discr_phantom) -data += odl.util.white_noise(ray_trafo.range) * np.mean(data) * 0.1 +data += odl.phantom.white_noise(ray_trafo.range) * np.mean(data) * 0.1 # Optionally pass partial to the solver to display intermediate results partial = (odl.solvers.PrintIterationPartial() & diff --git a/examples/tomo/filtered_backprojection_parallel_2d.py b/examples/tomo/filtered_backprojection_parallel_2d.py index 9ba0d8b44a1..ab410231ac4 100644 --- a/examples/tomo/filtered_backprojection_parallel_2d.py +++ b/examples/tomo/filtered_backprojection_parallel_2d.py @@ -70,7 +70,7 @@ # Create a discrete Shepp-Logan phantom (modified version) -phantom = odl.util.phantom.shepp_logan(reco_space, modified=True) +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) diff --git a/examples/tomo/ray_trafo_circular_cone.py b/examples/tomo/ray_trafo_circular_cone.py index 7bf4a869827..f99600811e0 100644 --- a/examples/tomo/ray_trafo_circular_cone.py +++ b/examples/tomo/ray_trafo_circular_cone.py @@ -45,7 +45,7 @@ ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda') # Create a discrete Shepp-Logan phantom (modified version) -phantom = odl.util.phantom.shepp_logan(reco_space, True) +phantom = odl.phantom.shepp_logan(reco_space, True) # Create projection data by calling the ray transform on the phantom proj_data = ray_trafo(phantom) diff --git a/examples/tomo/ray_trafo_fanflat.py b/examples/tomo/ray_trafo_fanflat.py index 9350084b888..71d570bd5b7 100644 --- a/examples/tomo/ray_trafo_fanflat.py +++ b/examples/tomo/ray_trafo_fanflat.py @@ -44,7 +44,7 @@ ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda') # Create a discrete Shepp-Logan phantom (modified version) -phantom = odl.util.phantom.shepp_logan(reco_space, modified=True) +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) diff --git a/examples/tomo/ray_trafo_helical.py b/examples/tomo/ray_trafo_helical.py index b0da1d4ec71..986056dd773 100644 --- a/examples/tomo/ray_trafo_helical.py +++ b/examples/tomo/ray_trafo_helical.py @@ -46,7 +46,7 @@ ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda') # Create a discrete Shepp-Logan phantom (modified version) -phantom = odl.util.phantom.shepp_logan(reco_space, modified=True) +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) diff --git a/examples/tomo/ray_trafo_parallel_2d.py b/examples/tomo/ray_trafo_parallel_2d.py index 6f33b714b4a..eb93b015880 100644 --- a/examples/tomo/ray_trafo_parallel_2d.py +++ b/examples/tomo/ray_trafo_parallel_2d.py @@ -43,7 +43,7 @@ ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda') # Create a discrete Shepp-Logan phantom (modified version) -phantom = odl.util.phantom.shepp_logan(reco_space, modified=True) +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) diff --git a/examples/tomo/ray_trafo_parallel_3d.py b/examples/tomo/ray_trafo_parallel_3d.py index 8c2a16ba920..67b3a404453 100644 --- a/examples/tomo/ray_trafo_parallel_3d.py +++ b/examples/tomo/ray_trafo_parallel_3d.py @@ -49,7 +49,7 @@ ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda') # Create a discrete Shepp-Logan phantom (modified version) -phantom = odl.util.phantom.shepp_logan(reco_space, modified=True) +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) diff --git a/examples/tomo/scikit_ray_trafo_parallel_2d.py b/examples/tomo/scikit_ray_trafo_parallel_2d.py index b6ee1e2675f..1fb7e2fa7e9 100644 --- a/examples/tomo/scikit_ray_trafo_parallel_2d.py +++ b/examples/tomo/scikit_ray_trafo_parallel_2d.py @@ -37,7 +37,7 @@ ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='scikit') # Create a discrete Shepp-Logan phantom (modified version) -phantom = odl.util.phantom.shepp_logan(reco_space, modified=True) +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) diff --git a/examples/tomo/stir_project.py b/examples/tomo/stir_project.py index 6ef60f80252..6f5860cd671 100644 --- a/examples/tomo/stir_project.py +++ b/examples/tomo/stir_project.py @@ -46,7 +46,7 @@ recon_sp, data_sp, volume, proj_data) # Create shepp-logan phantom -vol = odl.util.shepp_logan(proj.domain, modified=True) +vol = odl.phantom.shepp_logan(proj.domain, modified=True) # Project and show result = proj(vol) diff --git a/examples/tomo/stir_reconstruct.py b/examples/tomo/stir_reconstruct.py index b570f8f8465..a8645a7d0dc 100644 --- a/examples/tomo/stir_reconstruct.py +++ b/examples/tomo/stir_reconstruct.py @@ -35,7 +35,7 @@ volume_file, projection_file) # Create shepp-logan phantom -vol = odl.util.shepp_logan(proj.domain, modified=True) +vol = odl.phantom.shepp_logan(proj.domain, modified=True) # Project data projections = proj(vol) diff --git a/examples/trafos/fourier_trafo.py b/examples/trafos/fourier_trafo.py index 650d13fb037..a5e25f4025c 100644 --- a/examples/trafos/fourier_trafo.py +++ b/examples/trafos/fourier_trafo.py @@ -35,7 +35,7 @@ ft_op = odl.trafos.FourierTransform(discr_space) # Create a phantom and its Fourier transfrom and display them -phantom = odl.util.phantom.shepp_logan(discr_space, modified=True) +phantom = odl.phantom.shepp_logan(discr_space, modified=True) phantom.show(title='Shepp-Logan phantom', show=True) phantom_ft = ft_op(phantom) phantom_ft.show(title='Full Fourier Transform', show=True) @@ -51,7 +51,7 @@ # its complex conjugate. This is faster and more memory efficient. discr_space_real = odl.uniform_discr([-1, -1], [1, 1], (512, 512)) ft_op_halfc = odl.trafos.FourierTransform(discr_space_real, halfcomplex=True) -phantom_real = odl.util.phantom.shepp_logan(discr_space_real, modified=True) +phantom_real = odl.phantom.shepp_logan(discr_space_real, modified=True) phantom_real.show(title='Shepp-Logan phantom, real version', show=True) phantom_real_ft = ft_op_halfc(phantom_real) phantom_real_ft.show(title='Half-complex Fourier Transform', show=True) diff --git a/examples/visualization/show_2d.py b/examples/visualization/show_2d.py index ec9b632dd7d..a582c13bfb3 100644 --- a/examples/visualization/show_2d.py +++ b/examples/visualization/show_2d.py @@ -26,7 +26,7 @@ import odl spc = odl.uniform_discr([0, 0], [1, 1], [100, 100]) -vec = odl.util.shepp_logan(spc, modified=True) +vec = odl.phantom.shepp_logan(spc, modified=True) # Show all data vec.show(show=True) diff --git a/examples/visualization/show_2d_complex.py b/examples/visualization/show_2d_complex.py index cf655ef4742..726c2ff64f7 100644 --- a/examples/visualization/show_2d_complex.py +++ b/examples/visualization/show_2d_complex.py @@ -26,7 +26,7 @@ import odl spc = odl.uniform_discr([0, 0], [1, 1], [100, 100], dtype='complex') -vec = odl.util.shepp_logan(spc, modified=True) * (1 + 0.5j) +vec = odl.phantom.shepp_logan(spc, modified=True) * (1 + 0.5j) # Can also force "instant" plotting vec.show(show=True) diff --git a/examples/visualization/show_productspace.py b/examples/visualization/show_productspace.py index 514325f98dd..3205fb6fbe6 100644 --- a/examples/visualization/show_productspace.py +++ b/examples/visualization/show_productspace.py @@ -30,10 +30,10 @@ pspace = odl.ProductSpace(spc, m) # Making a product space element where each component consists of a -# Shepp-Logan phantom multiplied by the constant (i+1), where i is the +# Shepp-Logan phantom multiplied by the constant i, where i is the # index of the product space component. -vec = pspace.element([odl.util.shepp_logan(spc, modified=True) * i - for i in range(1, m + 1)]) +vec = pspace.element([odl.phantom.shepp_logan(spc, modified=True) * i + for i in range(m)]) # By default 4 uniformly spaced elements are shown. Since there are 7 in # total, the shown components are 0, 2, 4 and 6 diff --git a/examples/visualization/show_update_2d.py b/examples/visualization/show_update_2d.py index da4af16a9b0..3f337cbb21a 100644 --- a/examples/visualization/show_update_2d.py +++ b/examples/visualization/show_update_2d.py @@ -28,7 +28,7 @@ n = 100 m = 20 spc = odl.uniform_discr([0, 0], [1, 1], [n, n]) -vec = odl.util.shepp_logan(spc, modified=True) +vec = odl.phantom.shepp_logan(spc, modified=True) # Create a figure by saving the result of show fig = None diff --git a/odl/__init__.py b/odl/__init__.py index 7310dc77d2f..d4dd2a720d9 100644 --- a/odl/__init__.py +++ b/odl/__init__.py @@ -25,9 +25,9 @@ from __future__ import absolute_import -__version__ = '0.2.3' +__version__ = '0.2.4' __all__ = ('diagnostics', 'discr', 'operator', 'set', 'space', 'solvers', - 'tomo', 'trafos', 'util') + 'tomo', 'trafos', 'util', 'phantom') # Propagate names defined in __all__ of all submodules into the top-level # module @@ -53,6 +53,7 @@ from . import trafos from . import tomo from . import util +from . import phantom from .util import test __all__ += ('test',) diff --git a/odl/phantom/__init__.py b/odl/phantom/__init__.py new file mode 100644 index 00000000000..2a1ec595191 --- /dev/null +++ b/odl/phantom/__init__.py @@ -0,0 +1,39 @@ +# 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 . + +"""Test data for tomography problems.""" + +from __future__ import absolute_import + +__all__ = () + +from . import phantom_utils + +from .emission import * +__all__ += emission.__all__ + +from .geometric import * +__all__ += geometric.__all__ + +from .misc_phantoms import * +__all__ += misc_phantoms.__all__ + +from .noise import * +__all__ += noise.__all__ + +from .transmission import * +__all__ += transmission.__all__ diff --git a/odl/phantom/emission.py b/odl/phantom/emission.py new file mode 100644 index 00000000000..c59490236fe --- /dev/null +++ b/odl/phantom/emission.py @@ -0,0 +1,151 @@ +# 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 . + +"""Phantoms used in emission tomography.""" + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() + +from odl.phantom.phantom_utils import ellipse_phantom, cylinders_from_ellipses + + +__all__ = ('derenzo_sources',) + + +def _derenzo_sources_2d(): + """Return ellipse parameters for a 2d Derenzo sources phantom. + + This is a popular phantom in SPECT and PET. It defines the source + locations and intensities. + """ + return [[1.0, 0.047788, 0.047788, -0.77758, -0.11811, 0.0], + [1.0, 0.063525, 0.063525, -0.71353, 0.12182, 0.0], + [1.0, 0.047788, 0.047788, -0.68141, -0.28419, 0.0], + [1.0, 0.063525, 0.063525, -0.58552, 0.3433, 0.0], + [1.0, 0.047838, 0.047838, -0.58547, -0.45035, 0.0], + [1.0, 0.047591, 0.047591, -0.58578, -0.11798, 0.0], + [1.0, 0.047591, 0.047591, -0.48972, -0.61655, 0.0], + [1.0, 0.047739, 0.047739, -0.48973, -0.28414, 0.0], + [1.0, 0.063747, 0.063747, -0.45769, 0.12204, 0.0], + [1.0, 0.063673, 0.063673, -0.4578, 0.5649, 0.0], + [1.0, 0.04764, 0.04764, -0.39384, -0.45026, 0.0], + [1.0, 0.047591, 0.047591, -0.39381, -0.11783, 0.0], + [1.0, 0.063525, 0.063525, -0.32987, 0.3433, 0.0], + [1.0, 0.03167, 0.03167, -0.31394, -0.7915, 0.0], + [1.0, 0.047591, 0.047591, -0.29786, -0.28413, 0.0], + [1.0, 0.032112, 0.032112, -0.25, -0.68105, 0.0], + [1.0, 0.063488, 0.063488, -0.20192, 0.12185, 0.0], + [1.0, 0.047442, 0.047442, -0.20192, -0.11804, 0.0], + [1.0, 0.079552, 0.079552, -0.15405, 0.59875, 0.0], + [1.0, 0.031744, 0.031744, -0.1862, -0.79155, 0.0], + [1.0, 0.03167, 0.03167, -0.18629, -0.57055, 0.0], + [1.0, 0.031892, 0.031892, -0.12224, -0.68109, 0.0], + [1.0, 0.03167, 0.03167, -0.1217, -0.45961, 0.0], + [1.0, 0.032039, 0.032039, -0.05808, -0.79192, 0.0], + [1.0, 0.031744, 0.031744, -0.058285, -0.57011, 0.0], + [1.0, 0.03167, 0.03167, -0.05827, -0.3487, 0.0], + [1.0, 0.079434, 0.079434, 0.0057692, 0.32179, 0.0], + [1.0, 0.031892, 0.031892, 0.0057692, -0.68077, 0.0], + [1.0, 0.031446, 0.031446, 0.0057692, -0.45934, 0.0], + [1.0, 0.031892, 0.031892, 0.0057692, -0.23746, 0.0], + [1.0, 0.032039, 0.032039, 0.069619, -0.79192, 0.0], + [1.0, 0.031744, 0.031744, 0.069824, -0.57011, 0.0], + [1.0, 0.03167, 0.03167, 0.069809, -0.3487, 0.0], + [1.0, 0.079552, 0.079552, 0.16558, 0.59875, 0.0], + [1.0, 0.031892, 0.031892, 0.13378, -0.68109, 0.0], + [1.0, 0.03167, 0.03167, 0.13324, -0.45961, 0.0], + [1.0, 0.031744, 0.031744, 0.19774, -0.79155, 0.0], + [1.0, 0.03167, 0.03167, 0.19783, -0.57055, 0.0], + [1.0, 0.09533, 0.09533, 0.28269, 0.16171, 0.0], + [1.0, 0.023572, 0.023572, 0.21346, -0.11767, 0.0], + [1.0, 0.032112, 0.032112, 0.26154, -0.68105, 0.0], + [1.0, 0.023968, 0.023968, 0.26122, -0.20117, 0.0], + [1.0, 0.023968, 0.023968, 0.30933, -0.28398, 0.0], + [1.0, 0.023771, 0.023771, 0.30939, -0.11763, 0.0], + [1.0, 0.03167, 0.03167, 0.32548, -0.7915, 0.0], + [1.0, 0.024066, 0.024066, 0.35722, -0.36714, 0.0], + [1.0, 0.023968, 0.023968, 0.35703, -0.20132, 0.0], + [1.0, 0.09538, 0.09538, 0.47446, 0.49414, 0.0], + [1.0, 0.024066, 0.024066, 0.40532, -0.45053, 0.0], + [1.0, 0.024066, 0.024066, 0.40532, -0.28408, 0.0], + [1.0, 0.023671, 0.023671, 0.40537, -0.11771, 0.0], + [1.0, 0.02387, 0.02387, 0.45299, -0.53331, 0.0], + [1.0, 0.02387, 0.02387, 0.45305, -0.36713, 0.0], + [1.0, 0.02387, 0.02387, 0.45299, -0.2013, 0.0], + [1.0, 0.023671, 0.023671, 0.50152, -0.6169, 0.0], + [1.0, 0.023968, 0.023968, 0.50132, -0.45066, 0.0], + [1.0, 0.023968, 0.023968, 0.50132, -0.28395, 0.0], + [1.0, 0.023671, 0.023671, 0.50152, -0.11771, 0.0], + [1.0, 0.024066, 0.024066, 0.54887, -0.69934, 0.0], + [1.0, 0.023771, 0.023771, 0.54894, -0.5333, 0.0], + [1.0, 0.023771, 0.023771, 0.54872, -0.36731, 0.0], + [1.0, 0.023771, 0.023771, 0.54894, -0.20131, 0.0], + [1.0, 0.09533, 0.09533, 0.66643, 0.16163, 0.0], + [1.0, 0.02387, 0.02387, 0.59739, -0.61662, 0.0], + [1.0, 0.023968, 0.023968, 0.59748, -0.45066, 0.0], + [1.0, 0.023968, 0.023968, 0.59748, -0.28395, 0.0], + [1.0, 0.023572, 0.023572, 0.59749, -0.11763, 0.0], + [1.0, 0.023572, 0.023572, 0.64482, -0.53302, 0.0], + [1.0, 0.023671, 0.023671, 0.64473, -0.36716, 0.0], + [1.0, 0.02387, 0.02387, 0.64491, -0.20124, 0.0], + [1.0, 0.02387, 0.02387, 0.69317, -0.45038, 0.0], + [1.0, 0.024066, 0.024066, 0.69343, -0.28396, 0.0], + [1.0, 0.023771, 0.023771, 0.69337, -0.11792, 0.0], + [1.0, 0.023572, 0.023572, 0.74074, -0.36731, 0.0], + [1.0, 0.023671, 0.023671, 0.74079, -0.20152, 0.0], + [1.0, 0.023671, 0.023671, 0.78911, -0.28397, 0.0], + [1.0, 0.02387, 0.02387, 0.78932, -0.11793, 0.0], + [1.0, 0.023572, 0.023572, 0.83686, -0.20134, 0.0], + [1.0, 0.023968, 0.023968, 0.88528, -0.11791, 0.0]] + + +def derenzo_sources(space): + """Create the PET/SPECT Derenzo sources phantom. + + The Derenzo phantom contains a series of circles of decreasing size. + + In 3d the phantom is simply the 2d phantom extended in the z direction as + cylinders. + """ + if space.ndim == 2: + return ellipse_phantom(space, _derenzo_sources_2d()) + if space.ndim == 3: + return ellipse_phantom( + space, cylinders_from_ellipses(_derenzo_sources_2d())) + else: + raise ValueError('dimension not 2, no phantom available') + + +if __name__ == '__main__': + # Show the phantoms + import odl + n = 300 + + # 2D + discr = odl.uniform_discr([-1, -1], [1, 1], [n, n]) + derenzo_sources(discr).show('derenzo_sources 2d') + + # 3D + discr = odl.uniform_discr([-1, -1, -1], [1, 1, 1], [300, 300, 300]) + derenzo_sources(discr).show('derenzo_sources 3d') + + # Run also the doctests + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() diff --git a/odl/phantom/geometric.py b/odl/phantom/geometric.py new file mode 100644 index 00000000000..92223cb6772 --- /dev/null +++ b/odl/phantom/geometric.py @@ -0,0 +1,194 @@ +# 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 . + +"""Phantoms given by simple geometric objects such as cubes or spheres.""" + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() + +import numpy as np + +__all__ = ('cuboid', 'indicate_proj_axis') + + +def cuboid(discr_space, begin=None, end=None): + """Rectangular cuboid. + + Parameters + ---------- + discr_space : `DiscretizedSpace` + Discretized space in which the phantom is supposed to be created. + begin : array-like of length ``discr_space.ndim`` + The lower left corner of the cuboid within the space. + Default: A quarter of the volume from the minimum corner + end : array-like of length ``discr_space.ndim`` + The upper right corner of the cuboid within the space. + Default: A quarter of the volume from the maximum corner + + Returns + ------- + phantom : `LinearSpaceVector` + Returns an element in ``discr_space`` + + Examples + -------- + >>> import odl + >>> space = odl.uniform_discr(0, 1, 6, dtype='float32') + >>> print(cuboid(space, 0.5, 1)) + [0.0, 0.0, 0.0, 1.0, 1.0, 1.0] + >>> space = odl.uniform_discr([0, 0], [1, 1], [4, 6], dtype='float32') + >>> print(cuboid(space, [0.25, 0], [0.75, 0.5])) + [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [1.0, 1.0, 1.0, 0.0, 0.0, 0.0], + [1.0, 1.0, 1.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] + """ + minp = np.asarray(discr_space.domain.min()) + maxp = np.asarray(discr_space.domain.max()) + + if begin is None: + begin = minp * 0.75 + maxp * 0.25 + if end is None: + end = begin * 0.25 + maxp * 0.75 + + begin = np.atleast_1d(begin) + end = np.atleast_1d(end) + + assert begin.size == discr_space.ndim + assert end.size == discr_space.ndim + + def phan(x): + result = True + + for xp, minv, maxv in zip(x, begin, end): + result = (result & + np.less_equal(minv, xp) & np.less_equal(xp, maxv)) + return result + + return discr_space.element(phan) + + +def indicate_proj_axis(discr_space, scale_structures=0.5): + """Phantom indicating along which axis it is projected. + + The number (n) of rectangles in a parallel-beam projection along a main + axis (0, 1, or 2) indicates the projection to be along the (n-1)the + dimension. + + Parameters + ---------- + discr_space : `DiscretizedSpace` + Discretized space in which the phantom is supposed to be created + scale_structures : positive `float` in (0, 1] + Scales objects (cube, cuboids) + + Returns + ------- + phantom : `LinearSpaceVector` + Returns an element in ``discr_space`` + + Examples + -------- + >>> import odl + >>> space = odl.uniform_discr([0] * 3, [1] * 3, [8, 8, 8]) + >>> phan = indicate_proj_axis(space).asarray() + >>> print(np.sum(phan, 0)) + [[ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 3. 3. 0. 0. 0.] + [ 0. 0. 0. 3. 3. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.]] + >>> print(np.sum(phan, 1)) + [[ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 2. 2. 0. 0. 0.] + [ 0. 0. 0. 2. 2. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 1. 1. 0. 0. 0.] + [ 0. 0. 0. 1. 1. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.]] + >>> print(np.sum(phan, 2)) + [[ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 2. 2. 0. 0. 0.] + [ 0. 0. 0. 2. 2. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 0. 2. 0. 0. 0.] + [ 0. 0. 0. 2. 0. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0. 0. 0.]] + """ + if not 0 < scale_structures <= 1: + raise ValueError('`scale_structures` ({}) is not in (0, 1]' + ''.format(scale_structures)) + + assert discr_space.ndim == 3 + + shape = discr_space.shape + phan = np.zeros(shape) + shape = np.array(shape) - 1 + cen = np.round(0.5 * shape) + dx = np.floor(scale_structures * 0.25 * shape) + dx[dx == 0] = 1 + + # cube of size 2 * dx + x0 = (cen - 3 * dx)[0] + x, y, z = cen - 1 * dx + phan[x0:x, y:-y, z:-z] = 1 + + # 1st cuboid of size (dx[0], dx[1], 2 * dx[2]) + x0 = (cen + 1 * dx)[1] + x1 = (cen + 2 * dx)[1] + y0 = cen[1] + z = (cen - dx)[2] + phan[x0:x1, y0:-y, z:-z] = 1 + + # 2nd cuboid of (dx[0], dx[1], 2 * dx[2]) touching the first diagonally + # at a long edge + x0 = (cen + 2 * dx)[1] + x1 = (cen + 3 * dx)[1] + y1 = cen[1] + z = (cen - dx)[2] + phan[x0:x1, y:y1, z:-z] = 1 + + return discr_space.element(phan) + +if __name__ == '__main__': + # Show the phantoms + import odl + + # 1D + discr = odl.uniform_discr(-1, 1, 300) + cuboid(discr).show('cuboid 1d') + + # 2D + discr = odl.uniform_discr([-1, -1], [1, 1], [300, 300]) + cuboid(discr).show('cuboid 2d') + + # 3D + discr = odl.uniform_discr([-1, -1, -1], [1, 1, 1], [300, 300, 300]) + cuboid(discr).show('cuboid 3d') + indicate_proj_axis(discr).show('indicate_proj_axis 3d') + + # Run also the doctests + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() diff --git a/odl/phantom/misc_phantoms.py b/odl/phantom/misc_phantoms.py new file mode 100644 index 00000000000..189fc62ff53 --- /dev/null +++ b/odl/phantom/misc_phantoms.py @@ -0,0 +1,169 @@ +# 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 . + +"""Miscellaneous phantoms that do not fit in other categories.""" + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() + +import numpy as np + + +__all__ = ('submarine',) + + +def submarine(space, smooth=True, taper=20.0): + """Return a 'submarine' phantom consisting in an ellipsoid and a box. + + Parameters + ---------- + space : `DiscreteLp` + Discretized space in which the phantom is supposed to be created. + smooth : `bool`, optional + If `True`, the boundaries are smoothed out. Otherwise, the + function steps from 0 to 1 at the boundaries. + taper : `float`, optional + Tapering parameter for the boundary smoothing. Larger values + mean faster taper, i.e. sharper boundaries. + + Returns + ------- + phantom : ``discr`` element + """ + if space.ndim == 2: + if smooth: + return _submarine_2d_smooth(space, taper) + else: + return _submarine_2d_nonsmooth(space) + else: + raise ValueError('phantom only defined in 2 dimensions, got {}' + ''.format(space.ndim)) + + +def _submarine_2d_smooth(space, taper): + """Return a 2d smooth 'submarine' phantom.""" + + def logistic(x, c): + """Smoothed step function from 0 to 1, centered at 0.""" + return 1. / (1 + np.exp(-c * x)) + + def blurred_ellipse(x): + """Blurred characteristic function of an ellipse. + + If ``space.domain`` is a rectangle ``[0, 1] x [0, 1]``, + the ellipse is centered at ``(0.6, 0.3)`` and has half-axes + ``(0.4, 0.14)``. For other domains, the values are scaled + accordingly. + """ + halfaxes = np.array([0.4, 0.14]) * space.domain.extent() + center = np.array([0.6, 0.3]) * space.domain.extent() + center += space.domain.min() + + # Efficiently calculate |z|^2, z = (x - center) / radii + sq_ndist = np.zeros_like(x[0]) + for xi, rad, cen in zip(x, halfaxes, center): + sq_ndist = sq_ndist + ((xi - cen) / rad) ** 2 + + out = np.sqrt(sq_ndist) + out -= 1 + # Return logistic(taper * (1 - |z|)) + return logistic(out, -taper) + + def blurred_rect(x): + """Blurred characteristic function of a rectangle. + + If ``space.domain`` is a rectangle ``[0, 1] x [0, 1]``, + the rect has lower left ``(0.56, 0.4)`` and upper right + ``(0.76, 0.6)``. For other domains, the values are scaled + accordingly. + """ + xlower = np.array([0.56, 0.4]) * space.domain.extent() + xlower += space.domain.min() + xupper = np.array([0.76, 0.6]) * space.domain.extent() + xupper += space.domain.min() + + out = np.ones_like(x[0]) + for xi, low, upp in zip(x, xlower, xupper): + length = upp - low + out = out * (logistic((xi - low) / length, taper) * + logistic((upp - xi) / length, taper)) + return out + + out = space.element(blurred_ellipse) + out += space.element(blurred_rect) + return out.ufunc.minimum(1, out=out) + + +def _submarine_2d_nonsmooth(space): + """Return a 2d nonsmooth 'submarine' phantom.""" + + def ellipse(x): + """Characteristic function of an ellipse. + + If ``space.domain`` is a rectangle ``[0, 1] x [0, 1]``, + the ellipse is centered at ``(0.6, 0.3)`` and has half-axes + ``(0.4, 0.14)``. For other domains, the values are scaled + accordingly. + """ + halfaxes = np.array([0.4, 0.14]) * space.domain.extent() + center = np.array([0.6, 0.3]) * space.domain.extent() + center += space.domain.min() + + sq_ndist = np.zeros_like(x[0]) + for xi, rad, cen in zip(x, halfaxes, center): + sq_ndist = sq_ndist + ((xi - cen) / rad) ** 2 + + return np.where(sq_ndist <= 1, 1, 0) + + def rect(x): + """Characteristic function of a rectangle. + + If ``space.domain`` is a rectangle ``[0, 1] x [0, 1]``, + the rect has lower left ``(0.56, 0.4)`` and upper right + ``(0.76, 0.6)``. For other domains, the values are scaled + accordingly. + """ + xlower = np.array([0.56, 0.4]) * space.domain.extent() + xlower += space.domain.min() + xupper = np.array([0.76, 0.6]) * space.domain.extent() + xupper += space.domain.min() + + out = np.ones_like(x[0]) + for xi, low, upp in zip(x, xlower, xupper): + out = out * ((xi >= low) & (xi <= upp)) + return out + + out = space.element(ellipse) + out += space.element(rect) + return out.ufunc.minimum(1, out=out) + + +if __name__ == '__main__': + # Show the phantoms + import odl + + space = odl.uniform_discr([-1, -1], [1, 1], [300, 300]) + submarine(space, smooth=False).show('submarine smooth=False') + submarine(space, smooth=True).show('submarine smooth=True') + submarine(space, smooth=True, taper=50).show('submarine taper=50') + + # Run also the doctests + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() diff --git a/odl/phantom/noise.py b/odl/phantom/noise.py new file mode 100644 index 00000000000..47068674276 --- /dev/null +++ b/odl/phantom/noise.py @@ -0,0 +1,49 @@ +# 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 . + +"""Functions to create noise samples of different distributions.""" + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() + +import numpy as np + + +__all__ = ('white_noise',) + + +def white_noise(space): + """Standard gaussian noise in space, pointwise N(0, 1). + + """ + values = np.random.randn(*space.shape) + return space.element(values) + + +if __name__ == '__main__': + # Show the phantoms + import odl + + discr = odl.uniform_discr([-1, -1], [1, 1], [300, 300]) + white_noise(discr).show('white_noise') + + # Run also the doctests + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() diff --git a/odl/phantom/phantom_utils.py b/odl/phantom/phantom_utils.py new file mode 100644 index 00000000000..a97f77e5b75 --- /dev/null +++ b/odl/phantom/phantom_utils.py @@ -0,0 +1,327 @@ +# 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 . + +"""Utilities for creating phantoms.""" + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() + +import numpy as np + + +__all__ = ('ellipse_phantom', 'cylinders_from_ellipses') + + +def cylinders_from_ellipses(ellipses2d): + """Create 3d cylinders from ellipses.""" + ellipses2d = np.asarray(ellipses2d) + ellipses3d = np.zeros((ellipses2d.shape[0], 10)) + ellipses3d[:, [0, 1, 2, 4, 5, 7]] = ellipses2d + ellipses3d[:, 3] = 100000.0 + + return ellipses3d + + +def _getshapes_2d(center, max_radius, shape): + """Calculate indices and slices for the bounding box of a disk.""" + index_mean = shape * center + index_radius = max_radius / 2.0 * np.array(shape) + + min_idx = np.floor(index_mean - index_radius).astype(int) + max_idx = np.ceil(index_mean + index_radius).astype(int) + idx = [slice(minx, maxx) for minx, maxx in zip(min_idx, max_idx)] + shapes = [(idx[0], slice(None)), + (slice(None), idx[1])] + return idx, shapes + + +def ellipse_phantom_2d(space, ellipses): + """Create a phantom of ellipses in 2d space. + + Parameters + ---------- + space : `DiscreteLp` + The space the phantom should be generated in. + ellipses : list of lists + Each row should contain: + 'value', 'axis_1', 'axis_2', 'center_x', 'center_y', 'rotation' + The ellipses should be contained the he rectangle [-1, -1] x [1, 1]. + + Returns + ------- + phantom : `DiscreteLpVector` + The phantom + + Notes + ----- + This function is heavily optimized, achieving runtimes about 20 times + faster than "trivial" implementations. It is therefore recommended to use + in all phantoms. + + The main optimization is that it only considers a subset of all the + points when updating for each ellipse. It does this by first finding + a subset of points that could possibly be inside the ellipse. This + optimization is very good for "spherical" ellipsoids, but not so + much for elongated or rotated ones. + + It also does calculations wherever possible on the meshgrid instead of + individual points. + + See Also + -------- + shepp_logan : The typical use-case for this function. + """ + + # Blank image + p = np.zeros(space.shape, dtype=space.dtype) + + # Create the pixel grid + grid_in = space.grid.meshgrid + minp = space.grid.min() + maxp = space.grid.max() + + # move points to [-1, 1] + grid = [] + for i in range(2): + meani = (minp[i] + maxp[i]) / 2.0 + diffi = (maxp[i] - minp[i]) / 2.0 + grid += [(grid_in[i] - meani) / diffi] + + for ellip in ellipses: + intensity = ellip[0] + a_squared = ellip[1] ** 2 + b_squared = ellip[2] ** 2 + x0 = ellip[3] + y0 = ellip[4] + theta = ellip[5] * np.pi / 180 + + scales = [1 / a_squared, 1 / b_squared] + center = (np.array([x0, y0]) + 1.0) / 2.0 + + # Create the offset x,y and z values for the grid + if theta != 0: + # Rotate the points to the expected coordinate system. + ctheta = np.cos(theta) + stheta = np.sin(theta) + + mat = np.array([[ctheta, stheta], + [-stheta, ctheta]]) + + # Calculate the points that could possibly be inside the volume + # Since the points are rotated, we cannot do anything directional + # without more logic + max_radius = np.sqrt( + np.abs(mat).dot([a_squared, b_squared])) + idx, shapes = _getshapes_2d(center, max_radius, space.shape) + + subgrid = [g[idi] for g, idi in zip(grid, shapes)] + offset_points = [vec * (xi - x0i)[..., np.newaxis] + for xi, vec, x0i in zip(subgrid, + mat.T, + [x0, y0])] + rotated = offset_points[0] + offset_points[1] + np.square(rotated, out=rotated) + radius = np.dot(rotated, scales) + else: + # Calculate the points that could possibly be inside the volume + max_radius = np.sqrt([a_squared, b_squared]) + idx, shapes = _getshapes_2d(center, max_radius, space.shape) + + subgrid = [g[idi] for g, idi in zip(grid, shapes)] + squared_dist = [ai * (xi - x0i) ** 2 + for xi, ai, x0i in zip(subgrid, + scales, + [x0, y0])] + + # Parentheses to get best order for broadcasting + radius = squared_dist[0] + squared_dist[1] + + # Find the pixels within the ellipse + inside = radius <= 1 + + # Add the ellipse intensity to those pixels + p[idx][inside] += intensity + + return space.element(p) + + +def _getshapes_3d(center, max_radius, shape): + """Calculate indices and slices for the bounding box of a ball.""" + index_mean = shape * center + index_radius = max_radius / 2.0 * np.array(shape) + + min_idx = np.floor(index_mean - index_radius).astype(int) + max_idx = np.ceil(index_mean + index_radius).astype(int) + idx = [slice(minx, maxx) for minx, maxx in zip(min_idx, max_idx)] + shapes = [(idx[0], slice(None), slice(None)), + (slice(None), idx[1], slice(None)), + (slice(None), slice(None), idx[2])] + return idx, shapes + + +def ellipse_phantom_3d(space, ellipses): + """Create an ellipse phantom in 3d space. + + Parameters + ---------- + space : `DiscreteLp` + The space the phantom should be generated in. + ellipses : list of lists + Each row should contain: + 'value', 'axis_1', 'axis_2', 'axis_2', + 'center_x', 'center_y', 'center_z', + 'rotation_phi', 'rotation_theta', 'rotation_psi' + The ellipses should be contained the he rectangle + [-1, -1, -1] x [1, 1, 1]. + + Returns + ------- + phantom : `DiscreteLpVector` + The phantom + + Notes + ----- + This function is heavily optimized, achieving runtimes about 20 times + faster than "trivial" implementations. It is therefore recommended to use + in all phantoms. + + The main optimization is that it only considers a subset of all the + points when updating for each ellipse. It does this by first finding + a subset of points that could possibly be inside the ellipse. This + optimization is very good for "spherical" ellipsoids, but not so + much for elongated or rotated ones. + + It also does calculations wherever possible on the meshgrid instead of + individual points. + + See Also + -------- + shepp_logan : The typical use-case for this function. + """ + + # Blank image + p = np.zeros(space.shape, dtype=space.dtype) + + # Create the pixel grid + grid_in = space.grid.meshgrid + minp = space.grid.min() + maxp = space.grid.max() + + # move points to [-1, 1] + grid = [] + for i in range(3): + meani = (minp[i] + maxp[i]) / 2.0 + diffi = (maxp[i] - minp[i]) / 2.0 + grid += [(grid_in[i] - meani) / diffi] + + for ellip in ellipses: + intensity = ellip[0] + a_squared = ellip[1] ** 2 + b_squared = ellip[2] ** 2 + c_squared = ellip[3] ** 2 + x0 = ellip[4] + y0 = ellip[5] + z0 = ellip[6] + phi = ellip[7] * np.pi / 180 + theta = ellip[8] * np.pi / 180 + psi = ellip[9] * np.pi / 180 + + scales = [1 / a_squared, 1 / b_squared, 1 / c_squared] + center = (np.array([x0, y0, z0]) + 1.0) / 2.0 + + # Create the offset x,y and z values for the grid + if any([phi, theta, psi]): + # Rotate the points to the expected coordinate system. + cphi = np.cos(phi) + sphi = np.sin(phi) + ctheta = np.cos(theta) + stheta = np.sin(theta) + cpsi = np.cos(psi) + spsi = np.sin(psi) + + mat = np.array([[cpsi * cphi - ctheta * sphi * spsi, + cpsi * sphi + ctheta * cphi * spsi, + spsi * stheta], + [-spsi * cphi - ctheta * sphi * cpsi, + -spsi * sphi + ctheta * cphi * cpsi, + cpsi * stheta], + [stheta * sphi, + -stheta * cphi, + ctheta]]) + + # Calculate the points that could possibly be inside the volume + # Since the points are rotated, we cannot do anything directional + # without more logic + + max_radius = np.sqrt( + np.abs(mat).dot([a_squared, b_squared, c_squared])) + idx, shapes = _getshapes_3d(center, max_radius, space.shape) + + subgrid = [g[idi] for g, idi in zip(grid, shapes)] + offset_points = [vec * (xi - x0i)[..., np.newaxis] + for xi, vec, x0i in zip(subgrid, + mat.T, + [x0, y0, z0])] + rotated = offset_points[0] + offset_points[1] + offset_points[2] + np.square(rotated, out=rotated) + radius = np.dot(rotated, scales) + else: + # Calculate the points that could possibly be inside the volume + max_radius = np.sqrt([a_squared, b_squared, c_squared]) + idx, shapes = _getshapes_3d(center, max_radius, space.shape) + + subgrid = [g[idi] for g, idi in zip(grid, shapes)] + squared_dist = [ai * (xi - x0i) ** 2 + for xi, ai, x0i in zip(subgrid, + scales, + [x0, y0, z0])] + + # Parentheses to get best order for broadcasting + radius = squared_dist[0] + (squared_dist[1] + squared_dist[2]) + + # Find the pixels within the ellipse + inside = radius <= 1 + + # Add the ellipse intensity to those pixels + p[idx][inside] += intensity + + return space.element(p) + + +def ellipse_phantom(space, ellipses): + """Return a phantom given by ellipses. + + See Also + -------- + ellipse_phantom_2d : The 2d implementation + ellipse_phantom_3d : The 3d implementation + """ + + if space.ndim == 2: + return ellipse_phantom_2d(space, ellipses) + elif space.ndim == 3: + return ellipse_phantom_3d(space, ellipses) + else: + raise ValueError('dimension not 2 or 3, no phantom available') + + +if __name__ == '__main__': + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() diff --git a/odl/phantom/transmission.py b/odl/phantom/transmission.py new file mode 100644 index 00000000000..d639263dc1f --- /dev/null +++ b/odl/phantom/transmission.py @@ -0,0 +1,120 @@ +# 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 . + +"""Phantoms typically used in transmission tomography.""" + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() + +from odl.phantom.phantom_utils import ellipse_phantom + + +__all__ = ('shepp_logan',) + + +def _shepp_logan_ellipse_2d(): + """Return ellipse parameters for a 2d Shepp-Logan phantom.""" + # value axisx axisy x y rotation + return [[2.00, .6900, .9200, 0.0000, 0.0000, 0], + [-.98, .6624, .8740, 0.0000, -.0184, 0], + [-.02, .1100, .3100, 0.2200, 0.0000, -18], + [-.02, .1600, .4100, -.2200, 0.0000, 18], + [0.01, .2100, .2500, 0.0000, 0.3500, 0], + [0.01, .0460, .0460, 0.0000, 0.1000, 0], + [0.01, .0460, .0460, 0.0000, -.1000, 0], + [0.01, .0460, .0230, -.0800, -.6050, 0], + [0.01, .0230, .0230, 0.0000, -.6060, 0], + [0.01, .0230, .0460, 0.0600, -.6050, 0]] + + +def _shepp_logan_ellipse_3d(): + """Return ellipse parameters for a 3d Shepp-Logan phantom.""" + # value axisx axisy axisz, x y z rotation + return [[2.00, .6900, .9200, .810, 0.0000, 0.0000, 0.00, 0.0, 0, 0], + [-.98, .6624, .8740, .780, 0.0000, -.0184, 0.00, 0.0, 0, 0], + [-.02, .1100, .3100, .220, 0.2200, 0.0000, 0.00, -18, 0, 0], + [-.02, .1600, .4100, .280, -.2200, 0.0000, 0.00, 18., 0, 0], + [0.01, .2100, .2500, .410, 0.0000, 0.3500, 0.00, 0.0, 0, 0], + [0.01, .0460, .0460, .050, 0.0000, 0.1000, 0.00, 0.0, 0, 0], + [0.01, .0460, .0460, .050, 0.0000, -.1000, 0.00, 0.0, 0, 0], + [0.01, .0460, .0230, .050, -.0800, -.6050, 0.00, 0.0, 0, 0], + [0.01, .0230, .0230, .020, 0.0000, -.6060, 0.00, 0.0, 0, 0], + [0.01, .0230, .0460, .020, 0.0600, -.6050, 0.00, 0.0, 0, 0]] + + +def _modified_shepp_logan_ellipses(ellipses): + """Modify ellipses to give the modified Shepp-Logan phantom. + + Works for both 2d and 3d. + """ + intensities = [1.0, -0.8, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] + + for ellipse, intensity in zip(ellipses, intensities): + ellipse[0] = intensity + + +def shepp_logan(space, modified=False): + """Create a Shepp-Logan phantom. + + The standard `Shepp-Logan phantom`_ in 2 or 3 dimensions. + + Parameters + ---------- + space : `DiscreteLp` + The 2/3 dimension space that the phantom should be created in. + modified : `bool`, optional + True if the modified Shepp-Logan phantom should be given. + The modified phantom has greatly amplified contrast to aid in + visualization. + + References + ---------- + .. Shepp-Logan phantom: en.wikipedia.org/wiki/Shepp–Logan_phantom + """ + if space.ndim == 2: + ellipses = _shepp_logan_ellipse_2d() + elif space.ndim == 3: + ellipses = _shepp_logan_ellipse_3d() + else: + raise ValueError('dimension not 2 or 3, no phantom available') + + if modified: + _modified_shepp_logan_ellipses(ellipses) + + return ellipse_phantom(space, ellipses) + + +if __name__ == '__main__': + # Show the phantoms + import odl + + # 2D + discr = odl.uniform_discr([-1, -1], [1, 1], [1000, 1000]) + shepp_logan(discr, modified=True).show('shepp_logan 2d modified=True') + shepp_logan(discr, modified=False).show('shepp_logan 2d modified=False') + + # 3D + discr = odl.uniform_discr([-1, -1, -1], [1, 1, 1], [300, 300, 300]) + shepp_logan(discr, modified=True).show('shepp_logan 3d modified=True') + shepp_logan(discr, modified=False).show('shepp_logan 3d modified=False') + + # Run also the doctests + # pylint: disable=wrong-import-position + from odl.util.testutils import run_doctests + run_doctests() diff --git a/odl/util/__init__.py b/odl/util/__init__.py index 9fa61d2988b..8f3b2979c37 100644 --- a/odl/util/__init__.py +++ b/odl/util/__init__.py @@ -30,9 +30,6 @@ from .normalize import * __all__ += normalize.__all__ -from .phantom import * -__all__ += phantom.__all__ - from .graphics import * __all__ += graphics.__all__ diff --git a/odl/util/phantom.py b/odl/util/phantom.py deleted file mode 100644 index d3cb79f1350..00000000000 --- a/odl/util/phantom.py +++ /dev/null @@ -1,715 +0,0 @@ -# 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 . - -"""Some useful phantoms, mostly for tomography tests.""" - -# Imports for common Python 2/3 codebase -from __future__ import print_function, division, absolute_import -from future import standard_library -standard_library.install_aliases() - -import numpy as np - - -__all__ = ('ellipse_phantom_2d', 'ellipse_phantom_3d', - 'cuboid', 'indicate_proj_axis', - 'derenzo_sources', 'shepp_logan', 'submarine_phantom', - 'white_noise') - - -def _shepp_logan_ellipse_2d(): - """Return ellipse parameters for a 2d Shepp-Logan phantom. - - This is the standard phantom in 2d medical imaging. - """ - # value axisx axisy x y rotation - return [[2.00, .6900, .9200, 0.0000, 0.0000, 0], - [-.98, .6624, .8740, 0.0000, -.0184, 0], - [-.02, .1100, .3100, 0.2200, 0.0000, -18], - [-.02, .1600, .4100, -.2200, 0.0000, 18], - [0.01, .2100, .2500, 0.0000, 0.3500, 0], - [0.01, .0460, .0460, 0.0000, 0.1000, 0], - [0.01, .0460, .0460, 0.0000, -.1000, 0], - [0.01, .0460, .0230, -.0800, -.6050, 0], - [0.01, .0230, .0230, 0.0000, -.6060, 0], - [0.01, .0230, .0460, 0.0600, -.6050, 0]] - - -def _shepp_logan_ellipse_3d(): - """Return ellipse parameters for a 3d Shepp-Logan phantom. - - This is the standard phantom in 3d medical imaging. - """ - # value axisx axisy axisz, x y z rotation - return [[2.00, .6900, .9200, .810, 0.0000, 0.0000, 0.00, 0.0, 0, 0], - [-.98, .6624, .8740, .780, 0.0000, -.0184, 0.00, 0.0, 0, 0], - [-.02, .1100, .3100, .220, 0.2200, 0.0000, 0.00, -18, 0, 0], - [-.02, .1600, .4100, .280, -.2200, 0.0000, 0.00, 18., 0, 0], - [0.01, .2100, .2500, .410, 0.0000, 0.3500, 0.00, 0.0, 0, 0], - [0.01, .0460, .0460, .050, 0.0000, 0.1000, 0.00, 0.0, 0, 0], - [0.01, .0460, .0460, .050, 0.0000, -.1000, 0.00, 0.0, 0, 0], - [0.01, .0460, .0230, .050, -.0800, -.6050, 0.00, 0.0, 0, 0], - [0.01, .0230, .0230, .020, 0.0000, -.6060, 0.00, 0.0, 0, 0], - [0.01, .0230, .0460, .020, 0.0600, -.6050, 0.00, 0.0, 0, 0]] - - -def _modified_shepp_logan_ellipses(ellipses): - """Modify ellipses to give the modified shepp-logan phantom. - - Works for both 1d and 2d - """ - intensities = [1.0, -0.8, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] - - for ellipse, intensity in zip(ellipses, intensities): - ellipse[0] = intensity - - -def _derenzo_sources_2d(): - """Return ellipse parameters for a 2d Derenzo sources phantom. - - This is a popular phantom in SPECT and PET. It defines the source - locations and intensities. - """ - return [[1.0, 0.047788, 0.047788, -0.77758, -0.11811, 0.0], - [1.0, 0.063525, 0.063525, -0.71353, 0.12182, 0.0], - [1.0, 0.047788, 0.047788, -0.68141, -0.28419, 0.0], - [1.0, 0.063525, 0.063525, -0.58552, 0.3433, 0.0], - [1.0, 0.047838, 0.047838, -0.58547, -0.45035, 0.0], - [1.0, 0.047591, 0.047591, -0.58578, -0.11798, 0.0], - [1.0, 0.047591, 0.047591, -0.48972, -0.61655, 0.0], - [1.0, 0.047739, 0.047739, -0.48973, -0.28414, 0.0], - [1.0, 0.063747, 0.063747, -0.45769, 0.12204, 0.0], - [1.0, 0.063673, 0.063673, -0.4578, 0.5649, 0.0], - [1.0, 0.04764, 0.04764, -0.39384, -0.45026, 0.0], - [1.0, 0.047591, 0.047591, -0.39381, -0.11783, 0.0], - [1.0, 0.063525, 0.063525, -0.32987, 0.3433, 0.0], - [1.0, 0.03167, 0.03167, -0.31394, -0.7915, 0.0], - [1.0, 0.047591, 0.047591, -0.29786, -0.28413, 0.0], - [1.0, 0.032112, 0.032112, -0.25, -0.68105, 0.0], - [1.0, 0.063488, 0.063488, -0.20192, 0.12185, 0.0], - [1.0, 0.047442, 0.047442, -0.20192, -0.11804, 0.0], - [1.0, 0.079552, 0.079552, -0.15405, 0.59875, 0.0], - [1.0, 0.031744, 0.031744, -0.1862, -0.79155, 0.0], - [1.0, 0.03167, 0.03167, -0.18629, -0.57055, 0.0], - [1.0, 0.031892, 0.031892, -0.12224, -0.68109, 0.0], - [1.0, 0.03167, 0.03167, -0.1217, -0.45961, 0.0], - [1.0, 0.032039, 0.032039, -0.05808, -0.79192, 0.0], - [1.0, 0.031744, 0.031744, -0.058285, -0.57011, 0.0], - [1.0, 0.03167, 0.03167, -0.05827, -0.3487, 0.0], - [1.0, 0.079434, 0.079434, 0.0057692, 0.32179, 0.0], - [1.0, 0.031892, 0.031892, 0.0057692, -0.68077, 0.0], - [1.0, 0.031446, 0.031446, 0.0057692, -0.45934, 0.0], - [1.0, 0.031892, 0.031892, 0.0057692, -0.23746, 0.0], - [1.0, 0.032039, 0.032039, 0.069619, -0.79192, 0.0], - [1.0, 0.031744, 0.031744, 0.069824, -0.57011, 0.0], - [1.0, 0.03167, 0.03167, 0.069809, -0.3487, 0.0], - [1.0, 0.079552, 0.079552, 0.16558, 0.59875, 0.0], - [1.0, 0.031892, 0.031892, 0.13378, -0.68109, 0.0], - [1.0, 0.03167, 0.03167, 0.13324, -0.45961, 0.0], - [1.0, 0.031744, 0.031744, 0.19774, -0.79155, 0.0], - [1.0, 0.03167, 0.03167, 0.19783, -0.57055, 0.0], - [1.0, 0.09533, 0.09533, 0.28269, 0.16171, 0.0], - [1.0, 0.023572, 0.023572, 0.21346, -0.11767, 0.0], - [1.0, 0.032112, 0.032112, 0.26154, -0.68105, 0.0], - [1.0, 0.023968, 0.023968, 0.26122, -0.20117, 0.0], - [1.0, 0.023968, 0.023968, 0.30933, -0.28398, 0.0], - [1.0, 0.023771, 0.023771, 0.30939, -0.11763, 0.0], - [1.0, 0.03167, 0.03167, 0.32548, -0.7915, 0.0], - [1.0, 0.024066, 0.024066, 0.35722, -0.36714, 0.0], - [1.0, 0.023968, 0.023968, 0.35703, -0.20132, 0.0], - [1.0, 0.09538, 0.09538, 0.47446, 0.49414, 0.0], - [1.0, 0.024066, 0.024066, 0.40532, -0.45053, 0.0], - [1.0, 0.024066, 0.024066, 0.40532, -0.28408, 0.0], - [1.0, 0.023671, 0.023671, 0.40537, -0.11771, 0.0], - [1.0, 0.02387, 0.02387, 0.45299, -0.53331, 0.0], - [1.0, 0.02387, 0.02387, 0.45305, -0.36713, 0.0], - [1.0, 0.02387, 0.02387, 0.45299, -0.2013, 0.0], - [1.0, 0.023671, 0.023671, 0.50152, -0.6169, 0.0], - [1.0, 0.023968, 0.023968, 0.50132, -0.45066, 0.0], - [1.0, 0.023968, 0.023968, 0.50132, -0.28395, 0.0], - [1.0, 0.023671, 0.023671, 0.50152, -0.11771, 0.0], - [1.0, 0.024066, 0.024066, 0.54887, -0.69934, 0.0], - [1.0, 0.023771, 0.023771, 0.54894, -0.5333, 0.0], - [1.0, 0.023771, 0.023771, 0.54872, -0.36731, 0.0], - [1.0, 0.023771, 0.023771, 0.54894, -0.20131, 0.0], - [1.0, 0.09533, 0.09533, 0.66643, 0.16163, 0.0], - [1.0, 0.02387, 0.02387, 0.59739, -0.61662, 0.0], - [1.0, 0.023968, 0.023968, 0.59748, -0.45066, 0.0], - [1.0, 0.023968, 0.023968, 0.59748, -0.28395, 0.0], - [1.0, 0.023572, 0.023572, 0.59749, -0.11763, 0.0], - [1.0, 0.023572, 0.023572, 0.64482, -0.53302, 0.0], - [1.0, 0.023671, 0.023671, 0.64473, -0.36716, 0.0], - [1.0, 0.02387, 0.02387, 0.64491, -0.20124, 0.0], - [1.0, 0.02387, 0.02387, 0.69317, -0.45038, 0.0], - [1.0, 0.024066, 0.024066, 0.69343, -0.28396, 0.0], - [1.0, 0.023771, 0.023771, 0.69337, -0.11792, 0.0], - [1.0, 0.023572, 0.023572, 0.74074, -0.36731, 0.0], - [1.0, 0.023671, 0.023671, 0.74079, -0.20152, 0.0], - [1.0, 0.023671, 0.023671, 0.78911, -0.28397, 0.0], - [1.0, 0.02387, 0.02387, 0.78932, -0.11793, 0.0], - [1.0, 0.023572, 0.023572, 0.83686, -0.20134, 0.0], - [1.0, 0.023968, 0.023968, 0.88528, -0.11791, 0.0]] - - -def _make_3d_cylinders(ellipses2d): - """Create 3d cylinders from ellipses.""" - ellipses2d = np.asarray(ellipses2d) - ellipses3d = np.zeros((ellipses2d.shape[0], 10)) - ellipses3d[:, [0, 1, 2, 4, 5, 7]] = ellipses2d - ellipses3d[:, 3] = 100000.0 - - return ellipses3d - - -def ellipse_phantom_2d(space, ellipses): - """Create an ellipse phantom in 2d space. - - Parameters - ---------- - space : `DiscreteLp` - The space the phantom should be generated in. - ellipses : list of lists - Each row should contain: - 'value', 'axis_1', 'axis_2', 'center_x', 'center_y', 'rotation' - The ellipses should be contained the he rectangle [-1, -1] x [1, 1]. - - Returns - ------- - phantom : `DiscreteLpVector` - The phantom - """ - - # Blank image - p = np.zeros(space.size) - - # Create the pixel grid - points = space.points() - minp = space.grid.min() - maxp = space.grid.max() - - # move points to [0, 1] - points = (points - minp) / (maxp - minp) - - # move to [-1, 1] - points = points * 2 - 1 - - for ellip in ellipses: - intensity = ellip[0] - a_squared = ellip[1] ** 2 - b_squared = ellip[2] ** 2 - x0 = ellip[3] - y0 = ellip[4] - phi = ellip[5] * np.pi / 180.0 - - # Create the offset x and y values for the grid - offset_points = points - [x0, y0] - - cos_p = np.cos(phi) - sin_p = np.sin(phi) - - # Find the pixels within the ellipse - scales = [1 / a_squared, 1 / b_squared] - mat = [[cos_p, sin_p], - [-sin_p, cos_p]] - radius = np.dot(scales, np.dot(mat, offset_points.T) ** 2) - inside = radius <= 1 - - # Add the ellipse intensity to those pixels - p[inside] += intensity - - return space.element(p) - - -def _getshapes(center, max_radius, shape): - """Calculate indices and slices for the bounding box of a ball.""" - index_mean = shape * center - index_radius = max_radius / 2.0 * np.array(shape) - - min_idx = np.floor(index_mean - index_radius).astype(int) - max_idx = np.ceil(index_mean + index_radius).astype(int) - idx = [slice(minx, maxx) for minx, maxx in zip(min_idx, max_idx)] - shapes = [(idx[0], slice(None), slice(None)), - (slice(None), idx[1], slice(None)), - (slice(None), slice(None), idx[2])] - return idx, shapes - - -def ellipse_phantom_3d(space, ellipses): - """Create an ellipse phantom in 3d space. - - Parameters - ---------- - space : `DiscreteLp` - The space the phantom should be generated in. - ellipses : list of lists - Each row should contain: - 'value', 'axis_1', 'axis_2', 'axis_2', - 'center_x', 'center_y', 'center_z', - 'rotation_phi', 'rotation_theta', 'rotation_psi' - The ellipses should be contained the he rectangle - [-1, -1, -1] x [1, 1, 1]. - - Returns - ------- - phantom : `DiscreteLpVector` - The phantom - """ - - # Implementation notes: - # This code is optimized compared to the 2d case since it handles larger - # data volumes. - # - # The main optimization is that it only considers a subset of all the - # points when updating for each ellipse. It does this by first finding - # a subset of points that could possibly be inside the ellipse. This - # approximation is accurate for "spherical" ellipsoids, but not so - # accurate for elongated ones. - - # Blank image - p = np.zeros(space.shape) - - # Create the pixel grid - grid_in = space.grid.meshgrid - minp = space.grid.min() - maxp = space.grid.max() - - # move points to [-1, 1] - grid = [] - for i in range(3): - meani = (minp[i] + maxp[i]) / 2.0 - diffi = (maxp[i] - minp[i]) / 2.0 - grid += [(grid_in[i] - meani) / diffi] - - for ellip in ellipses: - intensity = ellip[0] - a_squared = ellip[1] ** 2 - b_squared = ellip[2] ** 2 - c_squared = ellip[3] ** 2 - x0 = ellip[4] - y0 = ellip[5] - z0 = ellip[6] - phi = ellip[7] * np.pi / 180 - theta = ellip[8] * np.pi / 180 - psi = ellip[9] * np.pi / 180 - - scales = [1 / a_squared, 1 / b_squared, 1 / c_squared] - - # Create the offset x,y and z values for the grid - if any([phi, theta, psi]): - # Rotate the points to the expected coordinate system. - cphi = np.cos(phi) - sphi = np.sin(phi) - ctheta = np.cos(theta) - stheta = np.sin(theta) - cpsi = np.cos(psi) - spsi = np.sin(psi) - - mat = np.array([[cpsi * cphi - ctheta * sphi * spsi, - cpsi * sphi + ctheta * cphi * spsi, - spsi * stheta], - [-spsi * cphi - ctheta * sphi * cpsi, - -spsi * sphi + ctheta * cphi * cpsi, - cpsi * stheta], - [stheta * sphi, - -stheta * cphi, - ctheta]]) - - # Calculate the points that could possibly be inside the volume - # Since the points are rotated, we cannot do anything directional - # without more logic - center = (np.array([x0, y0, z0]) + 1.0) / 2.0 - max_radius = np.sqrt( - np.abs(mat).dot([a_squared, b_squared, c_squared])) - idx, shapes = _getshapes(center, max_radius, space.shape) - - subgrid = [g[idi] for g, idi in zip(grid, shapes)] - offset_points = [vec * (xi - x0i)[..., np.newaxis] - for xi, vec, x0i in zip(subgrid, - mat.T, - [x0, y0, z0])] - rotated = offset_points[0] + offset_points[1] + offset_points[2] - np.square(rotated, out=rotated) - radius = np.dot(rotated, scales) - else: - # Calculate the points that could possibly be inside the volume - center = (np.array([x0, y0, z0]) + 1.0) / 2.0 - max_radius = np.sqrt([a_squared, b_squared, c_squared]) - idx, shapes = _getshapes(center, max_radius, space.shape) - - subgrid = [g[idi] for g, idi in zip(grid, shapes)] - squared_dist = [ai * (xi - x0i) ** 2 - for xi, ai, x0i in zip(subgrid, - scales, - [x0, y0, z0])] - - # Parentisis to get best order for broadcasting - radius = squared_dist[0] + (squared_dist[1] + squared_dist[2]) - - # Find the pixels within the ellipse - inside = radius <= 1 - - # Add the ellipse intensity to those pixels - p[idx][inside] += intensity - - return space.element(p) - - -def phantom(space, ellipses): - """Return a phantom given by ellipses.""" - - if space.ndim == 2: - return ellipse_phantom_2d(space, ellipses) - elif space.ndim == 3: - return ellipse_phantom_3d(space, ellipses) - else: - raise ValueError('dimension not 2 or 3, no phantom available') - - -def derenzo_sources(space): - """Create the PET/SPECT Derenzo sources phantom. - - The Derenzo phantom contains a series of circles of decreasing size. - - In 3d the phantom is simply the 2d phantom extended in the z direction as - cylinders. - """ - if space.ndim == 2: - return ellipse_phantom_2d(space, _derenzo_sources_2d()) - if space.ndim == 3: - return ellipse_phantom_3d( - space, _make_3d_cylinders(_derenzo_sources_2d())) - else: - raise ValueError('dimension not 2, no phantom available') - - -def shepp_logan(space, modified=False): - """Create a Shepp-Logan phantom. - - The standard Shepp-Logan phantom in 2 or 3 dimensions. - - References - ---------- - Wikipedia : https://en.wikipedia.org/wiki/Shepp%E2%80%93Logan_phantom - """ - if space.ndim == 2: - ellipses = _shepp_logan_ellipse_2d() - elif space.ndim == 3: - ellipses = _shepp_logan_ellipse_3d() - else: - raise ValueError('dimension not 2 or 3, no phantom available') - - if modified: - _modified_shepp_logan_ellipses(ellipses) - - return phantom(space, ellipses) - - -def submarine_phantom(discr, smooth=True, taper=20.0): - """Return a 'submarine' phantom consisting in an ellipsoid and a box. - - This phantom is used in [Okt2015]_ for shape-based reconstruction. - - Parameters - ---------- - discr : `DiscreteLp` - Discretized space in which the phantom is supposed to be created - smooth : `bool`, optional - If `True`, the boundaries are smoothed out. Otherwise, the - function steps from 0 to 1 at the boundaries. - taper : `float`, optional - Tapering parameter for the boundary smoothing. Larger values - mean faster taper, i.e. sharper boundaries. - - Returns - ------- - phantom : `DiscreteLpVector` - """ - if discr.ndim == 2: - if smooth: - return _submarine_phantom_2d_smooth(discr, taper) - else: - return _submarine_phantom_2d_nonsmooth(discr) - else: - raise ValueError('phantom only defined in 2 dimensions, got {}' - ''.format(discr.dim)) - - -def _submarine_phantom_2d_smooth(discr, taper): - """Return a 2d smooth 'submarine' phantom.""" - - def logistic(x, c): - """Smoothed step function from 0 to 1, centered at 0.""" - return 1. / (1 + np.exp(-c * x)) - - def blurred_ellipse(x): - """Blurred characteristic function of an ellipse. - - If ``discr.domain`` is a rectangle ``[-1, 1] x [-1, 1]``, - the ellipse is centered at ``(0.2, -0.4)`` and has half-axes - ``(0.8, 0.28)``. For other domains, the values are scaled - accordingly. - """ - halfaxes = np.array([0.8, 0.28]) * discr.domain.extent() / 2 - center = np.array([0.2, -0.4]) * discr.domain.extent() / 2 - - # Efficiently calculate |z|^2, z = (x - center) / radii - sq_ndist = np.zeros_like(x[0]) - for xi, rad, cen in zip(x, halfaxes, center): - sq_ndist = sq_ndist + ((xi - cen) / rad) ** 2 - - out = np.sqrt(sq_ndist) - out -= 1 - # Return logistic(taper * (1 - |z|)) - return logistic(out, -taper) - - def blurred_rect(x): - """Blurred characteristic function of a rectangle. - - If ``discr.domain`` is a rectangle ``[-1, 1] x [-1, 1]``, - the rect has lower left ``(0.12, -0.2)`` and upper right - ``(0.52, 0.2)``. For other domains, the values are scaled - accordingly. - """ - xlower = np.array([0.12, -0.2]) * discr.domain.extent() / 2 - xupper = np.array([0.52, 0.2]) * discr.domain.extent() / 2 - - out = np.ones_like(x[0]) - for xi, low, upp in zip(x, xlower, xupper): - length = upp - low - out = out * (logistic((xi - low) / length, taper) * - logistic((upp - xi) / length, taper)) - return out - - out = discr.element(blurred_ellipse) - out += discr.element(blurred_rect) - return out.ufunc.minimum(1, out=out) - - -def _submarine_phantom_2d_nonsmooth(discr): - """Return a 2d nonsmooth 'submarine' phantom.""" - - def ellipse(x): - """Characteristic function of an ellipse. - - If ``discr.domain`` is a rectangle ``[-1, 1] x [-1, 1]``, - the ellipse is centered at ``(0.2, -0.4)`` and has half-axes - ``(0.8, 0.28)``. For other domains, the values are scaled - accordingly. - """ - halfaxes = np.array([0.8, 0.28]) * discr.domain.extent() / 2 - center = np.array([0.2, -0.4]) * discr.domain.extent() / 2 - - sq_ndist = np.zeros_like(x[0]) - for xi, rad, cen in zip(x, halfaxes, center): - sq_ndist = sq_ndist + ((xi - cen) / rad) ** 2 - - return np.where(sq_ndist <= 1, 1, 0) - - def rect(x): - """Characteristic function of a rectangle. - - If ``discr.domain`` is a rectangle ``[-1, 1] x [-1, 1]``, - the rectangle has lower left ``(0.12, -0.2)`` and upper right - ``(0.52, 0.2)``. For other domains, the values are scaled - accordingly. - """ - xlower = np.array([0.12, -0.2]) * discr.domain.extent() / 2 - xupper = np.array([0.52, 0.2]) * discr.domain.extent() / 2 - - out = np.ones_like(x[0]) - for xi, low, upp in zip(x, xlower, xupper): - out = out * ((xi >= low) & (xi <= upp)) - return out - - out = discr.element(ellipse) - out += discr.element(rect) - return out.ufunc.minimum(1, out=out) - - -def cuboid(discr_space, begin, end): - """Rectangular cuboid. - - Parameters - ---------- - discr_space : `DiscretizedSpace` - Discretized space in which the phantom is supposed to be created - begin : array-like or `float` in [0, 1] - The lower left corner of the cuboid within the space grid relative - to the extend of the grid - end : array-like or `float` in [0, 1] - The upper right corner of the cuboid within the space grid relative - to the extend of the grid - - Returns - ------- - phantom : `LinearSpaceVector` - Returns an element in ``discr_space`` - - Examples - -------- - >>> import odl - >>> space = odl.uniform_discr(0, 1, 6, dtype='float32') - >>> print(cuboid(space, 0.5, 1)) - [0.0, 0.0, 0.0, 1.0, 1.0, 1.0] - >>> space = odl.uniform_discr([0, 0], [1, 1], [4, 6], dtype='float32') - >>> print(cuboid(space, [0.25, 0], [0.75, 0.5])) - [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [1.0, 1.0, 1.0, 0.0, 0.0, 0.0], - [1.0, 1.0, 1.0, 0.0, 0.0, 0.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] - """ - ndim = discr_space.ndim - shape = discr_space.shape - - if np.isscalar(begin): - begin = (begin,) * ndim - if np.isscalar(end): - end = (end,) * ndim - - # Create phantom - phan = np.zeros(shape) - - slice1 = [slice(None)] * ndim - - for nn in range(ndim): - start = np.floor(begin[nn] * shape[nn]).astype(int) - stop = np.ceil(end[nn] * shape[nn]).astype(int) - - slice1[nn] = slice(start, stop) - - phan[slice1] = 1 - - return discr_space.element(phan) - - -def indicate_proj_axis(discr_space, scale_structures=0.5): - """Phantom indicating along which axis it is projected. - - The number (n) of rectangles in a parallel-beam projection along a main - axis (0, 1, or 2) indicates the projection to be along the (n-1)the - dimension. - - Parameters - ---------- - discr_space : `DiscretizedSpace` - Discretized space in which the phantom is supposed to be created - scale_structures : positive `float` in (0, 1] - Scales objects (cube, cuboids) - - Returns - ------- - phantom : `LinearSpaceVector` - Returns an element in ``discr_space`` - - Examples - -------- - >>> import odl - >>> space = odl.uniform_discr([0] * 3, [1] * 3, [8, 8, 8]) - >>> phan = indicate_proj_axis(space).asarray() - >>> print(np.sum(phan, 0)) - [[ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 3. 3. 0. 0. 0.] - [ 0. 0. 0. 3. 3. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.]] - >>> print(np.sum(phan, 1)) - [[ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 2. 2. 0. 0. 0.] - [ 0. 0. 0. 2. 2. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 1. 1. 0. 0. 0.] - [ 0. 0. 0. 1. 1. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.]] - >>> print(np.sum(phan, 2)) - [[ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 2. 2. 0. 0. 0.] - [ 0. 0. 0. 2. 2. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 0. 2. 0. 0. 0.] - [ 0. 0. 0. 2. 0. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0. 0. 0.]] - """ - if not 0 < scale_structures <= 1: - raise ValueError('`scale_structures` ({}) is not in (0, 1]' - ''.format(scale_structures)) - - shape = discr_space.shape - phan = np.zeros(shape) - shape = np.array(shape) - 1 - cen = np.round(0.5 * shape) - dx = np.floor(scale_structures * 0.25 * shape) - dx[dx == 0] = 1 - - # cube of size 2 * dx - x0 = (cen - 3 * dx)[0] - x, y, z = cen - 1 * dx - phan[x0:x, y:-y, z:-z] = 1 - - # 1st cuboid of size (dx[0], dx[1], 2 * dx[2]) - x0 = (cen + 1 * dx)[1] - x1 = (cen + 2 * dx)[1] - y0 = cen[1] - z = (cen - dx)[2] - phan[x0:x1, y0:-y, z:-z] = 1 - - # 2nd cuboid of (dx[0], dx[1], 2 * dx[2]) touching the first diagonally - # at a long edge - x0 = (cen + 2 * dx)[1] - x1 = (cen + 3 * dx)[1] - y1 = cen[1] - z = (cen - dx)[2] - phan[x0:x1, y:y1, z:-z] = 1 - - return discr_space.element(phan) - - -def white_noise(space): - """Standard gaussian noise in space, pointwise N(0, 1).""" - values = np.random.randn(*space.shape) - return space.element(values) - - -if __name__ == '__main__': - # Show the phantoms - import odl - n = 300 - - # 2D - discr = odl.uniform_discr([-1, -1], [1, 1], [n, n]) - - shepp_logan(discr, modified=True).show() - shepp_logan(discr, modified=False).show() - derenzo_sources(discr).show() - submarine_phantom(discr, smooth=False).show() - submarine_phantom(discr, smooth=True).show() - submarine_phantom(discr, smooth=True, taper=50).show() - - # Shepp-logan 3d - discr = odl.uniform_discr([-1, -1, -1], [1, 1, 1], [n, n, n]) - with odl.util.Timer(): - shepp_logan_3d = shepp_logan(discr, modified=True) - shepp_logan_3d.show() - - # Run also the doctests - # pylint: disable=wrong-import-position - from odl.util.testutils import run_doctests - run_doctests() diff --git a/setup.py b/setup.py index 20917da7a43..596ee264ed0 100644 --- a/setup.py +++ b/setup.py @@ -89,7 +89,7 @@ def run_tests(self): setup( name='odl', - version='0.2.3', + version='0.2.4', description='Operator Discretization Library', long_description=long_description, diff --git a/test/discr/diff_ops_test.py b/test/discr/diff_ops_test.py index 539f92c2a48..d55eb548c4b 100644 --- a/test/discr/diff_ops_test.py +++ b/test/discr/diff_ops_test.py @@ -15,7 +15,7 @@ # You should have received a copy of the GNU General Public License # along with ODL. If not, see . -"""Unit tests for `tensor_ops`.""" +"""Unit tests for `diff_ops`.""" # Imports for common Python 2/3 codebase from __future__ import print_function, division, absolute_import @@ -30,7 +30,6 @@ finite_diff, PartialDerivative, Gradient, Divergence, Laplacian) from odl.util.testutils import ( all_equal, all_almost_equal, almost_equal, skip_if_no_cuda, never_skip) -from odl.util import cuboid methods = ['central', 'forward', 'backward'] @@ -342,7 +341,7 @@ def test_gradient(method, impl, padding): # DiscreteLp Vector space = odl.uniform_discr([0.] * ndim, [1.] * ndim, [lin_size] * ndim) - dom_vec = cuboid(space, [0.2] * ndim, [0.8] * ndim) + dom_vec = odl.phantom.cuboid(space, [0.2] * ndim, [0.8] * ndim) # gradient grad = Gradient(space, method=method, @@ -413,7 +412,7 @@ def test_divergence(method, impl, padding): div = Divergence(range=space, method=method, padding_method=padding_method, padding_value=padding_value) - dom_vec = cuboid(space, [0.2] * ndim, [0.8] * ndim) + dom_vec = odl.phantom.cuboid(space, [0.2] * ndim, [0.8] * ndim) div([dom_vec] * ndim) diff --git a/test/largescale/tomo/ray_transform_slow_test.py b/test/largescale/tomo/ray_transform_slow_test.py index 3c3f18adb03..63fe06f5fa7 100644 --- a/test/largescale/tomo/ray_transform_slow_test.py +++ b/test/largescale/tomo/ray_transform_slow_test.py @@ -170,8 +170,8 @@ def projector(request): def test_reconstruction(projector): """Test discrete Ray transform using ASTRA for reconstruction.""" - # Create shepp-logan phantom - vol = odl.util.shepp_logan(projector.domain, modified=True) + # Create Shepp-Logan phantom + vol = odl.phantom.shepp_logan(projector.domain, modified=True) # Project data projections = projector(vol) diff --git a/test/largescale/tomo/stir_slow_test.py b/test/largescale/tomo/stir_slow_test.py index 2daa4a2c3ad..e3acced282c 100644 --- a/test/largescale/tomo/stir_slow_test.py +++ b/test/largescale/tomo/stir_slow_test.py @@ -42,7 +42,7 @@ def test_from_file(): proj = stir_projector_from_file(volume_file, projection_file) # Create SPECT phantom - vol = odl.util.derenzo_sources(proj.domain) + vol = odl.phantom.derenzo_sources(proj.domain) # Project data projections = proj(vol) diff --git a/test/tomo/backends/astra_cpu_test.py b/test/tomo/backends/astra_cpu_test.py index 01aad606457..3cd4247b8b0 100644 --- a/test/tomo/backends/astra_cpu_test.py +++ b/test/tomo/backends/astra_cpu_test.py @@ -43,7 +43,7 @@ def test_astra_cpu_projector_parallel2d(): # Create reco space and a phantom reco_space = odl.uniform_discr([-4, -5], [4, 5], (4, 5), dtype='float32') - phantom = odl.util.phantom.cuboid(reco_space, begin=0.5, end=1) + phantom = odl.phantom.cuboid(reco_space, begin=[0, 0], end=[4, 5]) # Create parallel geometry angle_part = odl.uniform_partition(0, 2 * np.pi, 8) @@ -71,7 +71,7 @@ def test_astra_cpu_projector_fanflat(): # Create reco space and a phantom reco_space = odl.uniform_discr([-4, -5], [4, 5], (4, 5), dtype='float32') - phantom = odl.util.phantom.cuboid(reco_space, begin=0.5, end=1) + phantom = odl.phantom.cuboid(reco_space, begin=[0, 0], end=[4, 5]) # Create fan beam geometry with flat detector angle_part = odl.uniform_partition(0, 2 * np.pi, 8) diff --git a/test/tomo/backends/astra_cuda_test.py b/test/tomo/backends/astra_cuda_test.py index e29965922f4..2c19b707f05 100644 --- a/test/tomo/backends/astra_cuda_test.py +++ b/test/tomo/backends/astra_cuda_test.py @@ -43,7 +43,7 @@ def test_astra_cuda_projector_parallel2d(): # Create reco space and a phantom reco_space = odl.uniform_discr([-4, -5], [4, 5], (4, 5), dtype='float32') - phantom = odl.util.phantom.cuboid(reco_space, begin=0.5, end=1) + phantom = odl.phantom.cuboid(reco_space, begin=[0, 0], end=[4, 5]) # Create parallel geometry angle_part = odl.uniform_partition(0, 2 * np.pi, 8) @@ -71,7 +71,7 @@ def test_astra_cuda_projector_fanflat(): # Create reco space and a phantom reco_space = odl.uniform_discr([-4, -5], [4, 5], (4, 5), dtype='float32') - phantom = odl.util.phantom.cuboid(reco_space, begin=0.5, end=1) + phantom = odl.phantom.cuboid(reco_space, begin=[0, 0], end=[4, 5]) # Create fan beam geometry with flat detector angle_part = odl.uniform_partition(0, 2 * np.pi, 8) @@ -104,7 +104,7 @@ def test_astra_cuda_projector_parallel3d(): # Create reco space and a phantom reco_space = odl.uniform_discr([-4, -5, -6], [4, 5, 6], (4, 5, 6), dtype='float32') - phantom = odl.util.phantom.cuboid(reco_space, begin=0.5, end=1) + phantom = odl.phantom.cuboid(reco_space, begin=[0, 0, 0], end=[4, 5, 6]) # Create parallel geometry angle_part = odl.uniform_partition(0, 2 * np.pi, 8) @@ -133,7 +133,7 @@ def test_astra_gpu_projector_circular_conebeam(): # Create reco space and a phantom reco_space = odl.uniform_discr([-4, -5, -6], [4, 5, 6], (4, 5, 6), dtype='float32') - phantom = odl.util.phantom.cuboid(reco_space, begin=0.5, end=1) + phantom = odl.phantom.cuboid(reco_space, begin=[0, 0, 0], end=[4, 5, 6]) # Create circular cone beam geometry with flat detector angle_part = odl.uniform_partition(0, 2 * np.pi, 8) @@ -165,7 +165,7 @@ def test_astra_cuda_projector_helical_conebeam(): # Create reco space and a phantom reco_space = odl.uniform_discr([-4, -5, -6], [4, 5, 6], (4, 5, 6), dtype='float32') - phantom = odl.util.phantom.cuboid(reco_space, begin=0.5, end=1) + phantom = odl.phantom.cuboid(reco_space, begin=[0, 0, 0], end=[4, 5, 6]) # Create circular cone beam geometry with flat detector angle_part = odl.uniform_partition(0, 2 * np.pi, 8) diff --git a/test/tomo/backends/scikit_test.py b/test/tomo/backends/scikit_test.py index a005659b825..631342b6eba 100644 --- a/test/tomo/backends/scikit_test.py +++ b/test/tomo/backends/scikit_test.py @@ -39,7 +39,7 @@ def test_scikit_radon_projector_parallel2d(): # Create reco space and a phantom reco_space = odl.uniform_discr([-5, -5], [5, 5], (5, 5), dtype='float32') - phantom = odl.util.phantom.cuboid(reco_space, begin=0.5, end=1) + phantom = odl.phantom.cuboid(reco_space, begin=[0, 0], end=[5, 5]) # Create parallel geometry angle_part = odl.uniform_partition(0, 2 * np.pi, 8) diff --git a/test/tomo/operators/ray_trafo_test.py b/test/tomo/operators/ray_trafo_test.py index 49adca6edcb..03f7be48af4 100644 --- a/test/tomo/operators/ray_trafo_test.py +++ b/test/tomo/operators/ray_trafo_test.py @@ -162,7 +162,7 @@ def test_projector(projector): # Accept 10% errors places = 1 - # Create shepp-logan phantom + # Create Shepp-Logan phantom vol = projector.domain.one() # Calculate projection @@ -180,8 +180,8 @@ def test_adjoint(projector): # Accept 10% errors places = 1 - # Create shepp-logan phantom - vol = odl.util.shepp_logan(projector.domain, modified=True) + # Create Shepp-Logan phantom + vol = odl.phantom.shepp_logan(projector.domain, modified=True) # Calculate projection proj = projector(vol)