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.trafos.util import next_fast_len
 
 
 __all__ = ('fbp_op', 'fbp_filter_op', 'tam_danielson_window',
@@ -379,7 +380,8 @@ def fourier_filter(x):
         if padding:
             # Define padding operator
             ran_shp = (ray_trafo.range.shape[0],
-                       ray_trafo.range.shape[1] * 2 - 1)
+                       next_fast_len(
+                           ray_trafo.range.shape[1] * 2 - 1, impl, rfft=True))
             resizing = ResizingOperator(ray_trafo.range, ran_shp=ran_shp)
 
             fourier = FourierTransform(resizing.range, axes=1, impl=impl)
@@ -435,12 +437,14 @@ 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 = next_fast_len(
+                    ray_trafo.range.shape[1] * 2 - 1, impl, rfft=True)
             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 = next_fast_len(
+                    ray_trafo.range.shape[2] * 2 - 1, impl, rfft=True)
             else:
                 padded_shape_v = ray_trafo.range.shape[2]
 

odl/trafos/util/ft_utils.py

@@ -20,12 +20,13 @@
     is_real_dtype, is_numeric_dtype, is_real_floating_dtype,
     is_complex_floating_dtype, complex_dtype, dtype_repr,
     is_string,
-    normalized_scalar_param_list, normalized_axes_tuple)
+    normalized_scalar_param_list, normalized_axes_tuple, nextpow2)
 
 
 __all__ = ('reciprocal_grid', 'realspace_grid',
            'reciprocal_space',
-           'dft_preprocess_data', 'dft_postprocess_data')
+           'dft_preprocess_data', 'dft_postprocess_data',
+           'next_fast_len')
 
 
 def reciprocal_grid(grid, shift=True, axes=None, halfcomplex=False):
@@ -649,6 +650,19 @@ def reciprocal_space(space, axes=None, halfcomplex=False, shift=True,
     return recip_spc
 
 
+def next_fast_len(n, impl, rfft=False):
+    """The smallest size for which a fast fft implementation is available."""
+    n = int(n)
+    if impl == 'pyfftw':
+        if rfft:
+            return 2 ** nextpow2(n)
+        else:
+            import pyfftw
+            return pyfftw.next_fast_len(n)
+    elif impl == 'numpy':
+        return 2 ** nextpow2(n)
+
+
 if __name__ == '__main__':
     from doctest import testmod, NORMALIZE_WHITESPACE
     testmod(optionflags=NORMALIZE_WHITESPACE)

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 0
+    return int(np.ceil(np.log2(n)))
+
+
 if __name__ == '__main__':
     from odl.util.testutils import run_doctests
     run_doctests()