Improve fbp speed

examples/tomo/filtered_backprojection_parallel_2d_complex.py

@@ -34,7 +34,7 @@
 
 
 # Ray transform (= forward projection).
-ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda')
+ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cpu')
 
 # Create filtered back-projection operator
 fbp = odl.tomo.fbp_op(ray_trafo)

odl/tomo/analytic/filtered_back_projection.py

@@ -12,6 +12,7 @@
 
 from odl.discr import ResizingOperator
 from odl.trafos import FourierTransform, PYFFTW_AVAILABLE
+from odl.util import nextpow2
 
 
 __all__ = ('fbp_op', 'fbp_filter_op', 'tam_danielson_window',
@@ -379,7 +380,7 @@ def fourier_filter(x):
         if padding:
             # Define padding operator
             ran_shp = (ray_trafo.range.shape[0],
-                       ray_trafo.range.shape[1] * 2 - 1)
+                       nextpow2(ray_trafo.range.shape[1] * 2 - 1))
             resizing = ResizingOperator(ray_trafo.range, ran_shp=ran_shp)
 
             fourier = FourierTransform(resizing.range, axes=1, impl=impl)
@@ -435,12 +436,12 @@ def fourier_filter(x):
         if padding:
             # Define padding operator
             if used_axes[0]:
-                padded_shape_u = ray_trafo.range.shape[1] * 2 - 1
+                padded_shape_u = nextpow2(ray_trafo.range.shape[1] * 2 - 1)
             else:
                 padded_shape_u = ray_trafo.range.shape[1]
 
             if used_axes[1]:
-                padded_shape_v = ray_trafo.range.shape[2] * 2 - 1
+                padded_shape_v = nextpow2(ray_trafo.range.shape[2] * 2 - 1)
             else:
                 padded_shape_v = ray_trafo.range.shape[2]
 

odl/util/utility.py

@@ -47,6 +47,7 @@
     'repr_string',
     'attribute_repr_string',
     'method_repr_string',
+    'nextpow2',
     'run_from_ipython',
     'npy_random_seed',
     'unique',
@@ -1607,6 +1608,28 @@ def unique(seq):
         return unique_values
 
 
+def nextpow2(n):
+    """
+    Compute the integer which is a power of two.
+
+    Parameters
+    ----------
+    n : int
+
+    Examples
+    ========
+    >>> odl.util.nextpow2(0)
+    1
+    >>> odl.util.nextpow2(7)
+    8
+    >>> odl.util.nextpow2(513)
+    1024
+    """
+    if n == 0:
+        return 1
+    return 2 ** int(np.ceil(np.log2(n)))
+
+
 if __name__ == '__main__':
     from odl.util.testutils import run_doctests
     run_doctests()