adler-j · December 11, 2016 22:56 · Dec 11, 2016 · Dec 11, 2016 · Dec 11, 2016
diff --git a/boolean_radon_example.py b/boolean_radon_example.py
@@ -14,14 +14,16 @@
 import odl
 
 # Select the type of phantom to use
-phantom_type = 'touching_squares'
+phantom_type = 'torus'
 
 # Discrete reconstruction space: discretized functions on the rectangle
-# [-20, 20]^2 with 100 samples per dimension.
+# [-20, 20]^2 with 300 samples per dimension.
 space = odl.uniform_discr(
     min_pt=[-20, -20], max_pt=[20, 20], shape=[300, 300])
 
 # Create phantom
+# Note that to create a phantom from e.g. a numpy array, you would do
+# phantom = space.element(my_numpy_array)
 if phantom_type == 'shepp_logan':
     # Create a discrete Shepp-Logan phantom (modified version)
     phantom = odl.phantom.shepp_logan(space, modified=True)

diff --git a/boolean_radon_example.py b/boolean_radon_example.py
@@ -19,7 +19,7 @@
 # Discrete reconstruction space: discretized functions on the rectangle
 # [-20, 20]^2 with 100 samples per dimension.
 space = odl.uniform_discr(
-    min_pt=[-20, -20], max_pt=[20, 20], shape=[100, 100])
+    min_pt=[-20, -20], max_pt=[20, 20], shape=[300, 300])
 
 # Create phantom
 if phantom_type == 'shepp_logan':

diff --git a/boolean_radon_example.py b/boolean_radon_example.py
@@ -0,0 +1,74 @@
+"""Example of using the boolean ray transform with ODL.
+
+This example is with the 2d ray transform, but you could easily change the
+geometry to 3d and get the example you gave. See
+
+github.com/odlgroup/odl/blob/master/examples/tomo/ray_trafo_parallel_3d.py
+
+for an example.
+
+Note that some rounding errors may occur.
+"""
+
+import numpy as np
+import odl
+
+# Select the type of phantom to use
+phantom_type = 'touching_squares'
+
+# Discrete reconstruction space: discretized functions on the rectangle
+# [-20, 20]^2 with 100 samples per dimension.
+space = odl.uniform_discr(
+    min_pt=[-20, -20], max_pt=[20, 20], shape=[100, 100])
+
+# Create phantom
+if phantom_type == 'shepp_logan':
+    # Create a discrete Shepp-Logan phantom (modified version)
+    phantom = odl.phantom.shepp_logan(space, modified=True)
+elif phantom_type == 'touching_squares':
+    phantom = odl.phantom.cuboid(space, min_pt=[-10, -10], max_pt=[-0, -0])
+    phantom += odl.phantom.cuboid(space, min_pt=[0, 0], max_pt=[10, 10])
+elif phantom_type == 'separate_squares':
+    phantom = odl.phantom.cuboid(space, min_pt=[-10, -10], max_pt=[-1, -1])
+    phantom += odl.phantom.cuboid(space, min_pt=[1, 1], max_pt=[10, 10])
+elif phantom_type == 'torus':
+    ellipses = [[1, 0.3, 0.3, -0.5, 0, 0],
+                [1, 0.3, 0.3, 0.5, 0, 0]]
+    phantom = odl.phantom.ellipse_phantom(space, ellipses)
+else:
+    raise RuntimeError('unknown phantom_type')
+
+# Make a parallel beam geometry with flat detector
+# Angles: uniformly spaced, n = 700, min = 0, max = pi
+angle_partition = odl.uniform_partition(0, np.pi, 700)
+# Detector: uniformly sampled, n = 700, min = -28, max = 28
+detector_partition = odl.uniform_partition(-28, 28, 700)
+geometry = odl.tomo.Parallel2dGeometry(angle_partition, detector_partition)
+
+# Ray transform (= forward projection). We use scikit backend
+# downloadable by "pip install scikit-image"
+# The solver becomes MUCH faster with "impl='astra_cuda'" but this requires
+# ASTRA to be installed.
+ray_trafo = odl.tomo.RayTransform(space, geometry, impl='scikit')
+
+# Create projection data by calling the ray transform on the phantom
+proj_data = ray_trafo(phantom)
+
+# Create boolean projection data by taking positive values
+boolean_proj_data = np.greater(proj_data, 0)
+
+# Back-projection can be done by simply calling the adjoint operator on the
+# projection data (or any element in the projection space).
+backproj = ray_trafo.adjoint(boolean_proj_data)
+
+# Inverse is the data that is in all projections
+# Here we subtract a small number to account for numerics.
+inverse_result = np.greater_equal(backproj, angle_partition.extent() - 0.0001)
+
+# Shows a slice of the phantom, projections, and reconstruction
+phantom.show(title='Phantom')
+proj_data.show(title='Projection data (sinogram)')
+boolean_proj_data.show(title='Boolean projection data')
+backproj.show(title='Back-projected data')
+inverse_result.show(title='Inverse')
+(phantom - inverse_result).show('Error')