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)