diff --git a/examples/tomo/filtered_backprojection_parallel_2d_complex.py b/examples/tomo/filtered_backprojection_parallel_2d_complex.py index 5d443fdc525..b37f17c7798 100644 --- a/examples/tomo/filtered_backprojection_parallel_2d_complex.py +++ b/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) diff --git a/odl/tomo/analytic/filtered_back_projection.py b/odl/tomo/analytic/filtered_back_projection.py index 1b55192e077..dd356d0766a 100644 --- a/odl/tomo/analytic/filtered_back_projection.py +++ b/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] diff --git a/odl/trafos/util/ft_utils.py b/odl/trafos/util/ft_utils.py index 2e7bb31d8f5..2e621155fa5 100644 --- a/odl/trafos/util/ft_utils.py +++ b/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) diff --git a/odl/util/utility.py b/odl/util/utility.py index 45297db4bc8..27dc7a11fb7 100644 --- a/odl/util/utility.py +++ b/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()