From 0dbddb03fbd1d75b3c67cb3af072e35db6487f5c Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Sun, 26 Jan 2020 13:23:23 +0100 Subject: [PATCH 01/39] Simplify adjoint method --- odl/tomo/operators/ray_trafo.py | 34 ++++++++++++++++----------------- 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index f5e3346ed09..6e41e293637 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -440,15 +440,14 @@ def adjoint(self): ------- adjoint : `RayBackProjection` """ - if self._adjoint is not None: - return self._adjoint - - kwargs = self._extra_kwargs.copy() - kwargs['domain'] = self.range - self._adjoint = RayBackProjection(self.domain, self.geometry, - impl=self.impl, - use_cache=self.use_cache, - **kwargs) + if self._adjoint is None: + kwargs = self._extra_kwargs.copy() + kwargs['domain'] = self.range + self._adjoint = RayBackProjection(self.domain, self.geometry, + impl=self.impl, + use_cache=self.use_cache, + **kwargs) + return self._adjoint @@ -554,15 +553,14 @@ def adjoint(self): ------- adjoint : `RayTransform` """ - if self._adjoint is not None: - return self._adjoint - - kwargs = self._extra_kwargs.copy() - kwargs['range'] = self.domain - self._adjoint = RayTransform(self.range, self.geometry, - impl=self.impl, - use_cache=self.use_cache, - **kwargs) + if self._adjoint is None: + kwargs = self._extra_kwargs.copy() + kwargs['range'] = self.domain + self._adjoint = RayTransform(self.range, self.geometry, + impl=self.impl, + use_cache=self.use_cache, + **kwargs) + return self._adjoint From 4a4b9342f1821e239091be0acd71eba42cfe3f5b Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Sun, 26 Jan 2020 19:31:19 +0100 Subject: [PATCH 02/39] Make RayTransform operators independent from their implementations --- odl/tomo/backends/__init__.py | 2 + odl/tomo/backends/astra_cpu.py | 36 +++- odl/tomo/backends/astra_cuda.py | 267 +++++++++----------------- odl/tomo/backends/impl.py | 74 +++++++ odl/tomo/backends/skimage_radon.py | 55 +++++- odl/tomo/operators/ray_trafo.py | 298 +++++++++++------------------ 6 files changed, 370 insertions(+), 362 deletions(-) create mode 100644 odl/tomo/backends/impl.py diff --git a/odl/tomo/backends/__init__.py b/odl/tomo/backends/__init__.py index d942667ffdc..6c30076eb75 100644 --- a/odl/tomo/backends/__init__.py +++ b/odl/tomo/backends/__init__.py @@ -13,10 +13,12 @@ from .astra_cpu import * from .astra_cuda import * from .astra_setup import * +from .impl import * from .skimage_radon import * __all__ = () __all__ += astra_cpu.__all__ __all__ += astra_cuda.__all__ __all__ += astra_setup.__all__ +__all__ += impl.__all__ __all__ += skimage_radon.__all__ diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index 7f1c70321eb..cec1a53a267 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -16,6 +16,7 @@ from odl.tomo.backends.astra_setup import ( astra_algorithm, astra_data, astra_projection_geometry, astra_projector, astra_volume_geometry) +from odl.tomo.backends.impl import RayTransformImplBase from odl.tomo.geometry import ( DivergentBeamGeometry, Geometry, ParallelBeamGeometry) from odl.util import writable_array @@ -29,6 +30,7 @@ 'astra_cpu_forward_projector', 'astra_cpu_back_projector', 'default_astra_proj_type', + 'AstraCpuRayTransformImpl', ) @@ -201,9 +203,11 @@ def astra_cpu_back_projector(proj_data, geometry, vol_space, out=None, raise TypeError("`vol_space.impl` must be 'numpy', got {!r}" "".format(vol_space.impl)) if vol_space.ndim != geometry.ndim: - raise ValueError('dimensions {} of reconstruction space and {} of ' - 'geometry do not match'.format( - vol_space.ndim, geometry.ndim)) + raise ValueError( + 'dimensions {} of reconstruction space and {} of ' + 'geometry do not match' + ''.format(vol_space.ndim, geometry.ndim) + ) if out is None: out = vol_space.element() else: @@ -251,6 +255,32 @@ def astra_cpu_back_projector(proj_data, geometry, vol_space, out=None, return out +class AstraCpuRayTransformImpl(RayTransformImplBase): + def __init__(self, geometry, reco_space, proj_space): + super().__init__(geometry, reco_space, proj_space) + + @staticmethod + def can_handle_size(size): + return not size >= 512 ** 2 + + @classmethod + def supports_geometry(cls, geometry): + if geometry.ndim > 2: + return ValueError('`impl` {!r} only works for 2d' + ''.format(cls.__name__)) + + return True + + def call_forward(self, x_real, out_real, **kwargs): + return astra_cpu_forward_projector( + x_real, self.geometry, self.proj_space.real_space, out_real, **kwargs) + + def call_backward(self, x_real, out_real, **kwargs): + return astra_cpu_back_projector( + x_real, self.geometry, self.reco_space.real_space, out_real, **kwargs) + + if __name__ == '__main__': from odl.util.testutils import run_doctests + run_doctests() diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 69627888b23..4323ba7748d 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -10,38 +10,39 @@ from __future__ import absolute_import, division, print_function -from builtins import object from multiprocessing import Lock import numpy as np +import warnings from packaging.version import parse as parse_version from odl.discr import DiscretizedSpace from odl.tomo.backends.astra_setup import ( ASTRA_VERSION, astra_algorithm, astra_data, astra_projection_geometry, - astra_projector, astra_volume_geometry) + astra_projector, astra_volume_geometry, astra_supports, astra_versions_supporting) from odl.tomo.geometry import ( ConeBeamGeometry, FanBeamGeometry, Geometry, Parallel2dGeometry, Parallel3dAxisGeometry) +from odl.tomo.backends.impl import RayTransformImplBase try: import astra + ASTRA_CUDA_AVAILABLE = astra.astra.use_cuda() except ImportError: ASTRA_CUDA_AVAILABLE = False __all__ = ( 'ASTRA_CUDA_AVAILABLE', - 'AstraCudaProjectorImpl', - 'AstraCudaBackProjectorImpl', + 'AstraCudaRayTransformImpl', ) -class AstraCudaProjectorImpl(object): - +class AstraCudaRayTransformImpl(RayTransformImplBase): """Thin wrapper around ASTRA.""" - algo_id = None + algo_forward_id = None + algo_backward_id = None vol_id = None sino_id = None proj_id = None @@ -59,68 +60,37 @@ def __init__(self, geometry, reco_space, proj_space): proj_space : `DiscretizedSpace` Projection space, the space of the result. """ - assert isinstance(geometry, Geometry) - assert isinstance(reco_space, DiscretizedSpace) - assert isinstance(proj_space, DiscretizedSpace) - - self.geometry = geometry - self.reco_space = reco_space - self.proj_space = proj_space + super().__init__(geometry, reco_space, proj_space) self.create_ids() # ASTRA projectors are not thread-safe, thus we need to lock ourselves self._mutex = Lock() - def call_forward(self, vol_data, out=None): - """Run an ASTRA forward projection on the given data using the GPU. - - Parameters - ---------- - vol_data : ``reco_space`` element - Volume data to which the projector is applied. - out : ``proj_space`` element, optional - Element of the projection space to which the result is written. If - ``None``, an element in `proj_space` is created. - - Returns - ------- - out : ``proj_space`` element - Projection data resulting from the application of the projector. - If ``out`` was provided, the returned object is a reference to it. - """ - with self._mutex: - assert vol_data in self.reco_space - if out is not None: - assert out in self.proj_space - else: - out = self.proj_space.element() - - # Copy data to GPU memory - if self.geometry.ndim == 2: - astra.data2d.store(self.vol_id, vol_data.asarray()) - elif self.geometry.ndim == 3: - astra.data3d.store(self.vol_id, vol_data.asarray()) - else: - raise RuntimeError('unknown ndim') - - # Run algorithm - astra.algorithm.run(self.algo_id) - - # Copy result to host - if self.geometry.ndim == 2: - out[:] = self.out_array - elif self.geometry.ndim == 3: - out[:] = np.swapaxes(self.out_array, 0, 1).reshape( - self.proj_space.shape) - - # Fix scaling to weight by pixel size - if (isinstance(self.geometry, Parallel2dGeometry) and - parse_version(ASTRA_VERSION) < parse_version('1.9.9.dev')): - # parallel2d scales with pixel stride - out *= 1 / float(self.geometry.det_partition.cell_volume) - - return out + @classmethod + def supports_geometry(cls, geometry): + # Print a warning if the detector midpoint normal vector at any + # angle is perpendicular to the geometry axis in parallel 3d + # single-axis geometry -- this is broken in some ASTRA versions + if (isinstance(geometry, Parallel3dAxisGeometry) and + not astra_supports('par3d_det_mid_pt_perp_to_axis')): + req_ver = astra_versions_supporting('par3d_det_mid_pt_perp_to_axis') + axis = geometry.axis + mid_pt = geometry.det_params.mid_pt + for i, angle in enumerate(geometry.angles): + if abs(np.dot(axis, + geometry.det_to_src(angle, mid_pt))) < 1e-4: + warnings.warn( + 'angle {}: detector midpoint normal {} is ' + 'perpendicular to the geometry axis {} in ' + '`Parallel3dAxisGeometry`; this is broken in ' + 'ASTRA {}, please upgrade to ASTRA {}' + ''.format(i, geometry.det_to_src(angle, mid_pt), + axis, ASTRA_VERSION, req_ver), + RuntimeWarning) + break + + return True def create_ids(self): """Create ASTRA objects.""" @@ -143,10 +113,8 @@ def create_ids(self): astra_proj_shape = (proj_shape[1], proj_shape[0], proj_shape[2]) astra_vol_shape = self.reco_space.shape - self.in_array = np.empty(astra_vol_shape, - dtype='float32', order='C') - self.out_array = np.empty(astra_proj_shape, - dtype='float32', order='C') + self.vol_array = np.empty(astra_vol_shape, dtype='float32', order='C') + self.proj_array = np.empty(astra_proj_shape, dtype='float32', order='C') # Create ASTRA data structures vol_geom = astra_volume_geometry(self.reco_space) @@ -154,7 +122,7 @@ def create_ids(self): self.vol_id = astra_data(vol_geom, datatype='volume', ndim=self.reco_space.ndim, - data=self.in_array, + data=self.vol_array, allow_copy=False) proj_type = 'cuda' if proj_ndim == 2 else 'cuda3d' @@ -165,70 +133,70 @@ def create_ids(self): self.sino_id = astra_data(proj_geom, datatype='projection', ndim=proj_ndim, - data=self.out_array, + data=self.proj_array, allow_copy=False) # Create algorithm - self.algo_id = astra_algorithm( + self.algo_forward_id = astra_algorithm( 'forward', proj_ndim, self.vol_id, self.sino_id, proj_id=self.proj_id, impl='cuda') - def __del__(self): - """Delete ASTRA objects.""" - if self.geometry.ndim == 2: - adata, aproj = astra.data2d, astra.projector - else: - adata, aproj = astra.data3d, astra.projector3d - - if self.algo_id is not None: - astra.algorithm.delete(self.algo_id) - self.algo_id = None - if self.vol_id is not None: - adata.delete(self.vol_id) - self.vol_id = None - if self.sino_id is not None: - adata.delete(self.sino_id) - self.sino_id = None - if self.proj_id is not None: - aproj.delete(self.proj_id) - self.proj_id = None + # Create algorithm + self.algo_backward_id = astra_algorithm( + 'backward', proj_ndim, self.vol_id, self.sino_id, + proj_id=self.proj_id, impl='cuda') + def call_forward(self, vol_data, out, **kwargs): + """Run an ASTRA forward projection on the given data using the GPU. -class AstraCudaBackProjectorImpl(object): + Parameters + ---------- + vol_data : ``reco_space`` element + Volume data to which the projector is applied. + out : ``proj_space`` element, optional + Element of the projection space to which the result is written. If + ``None``, an element in `proj_space` is created. - """Thin wrapper around ASTRA.""" + Returns + ------- + out : ``proj_space`` element + Projection data resulting from the application of the projector. + If ``out`` was provided, the returned object is a reference to it. + """ + with self._mutex: + assert vol_data in self.reco_space + if out is not None: + assert out in self.proj_space + else: + out = self.proj_space.element() - algo_id = None - vol_id = None - sino_id = None - proj_id = None + # Copy data to GPU memory + if self.geometry.ndim == 2: + astra.data2d.store(self.vol_id, vol_data.asarray()) + elif self.geometry.ndim == 3: + astra.data3d.store(self.vol_id, vol_data.asarray()) + else: + raise RuntimeError('unknown ndim') - def __init__(self, geometry, reco_space, proj_space): - """Initialize a new instance. + # Run algorithm + astra.algorithm.run(self.algo_forward_id) - Parameters - ---------- - geometry : `Geometry` - Geometry defining the tomographic setup. - reco_space : `DiscretizedSpace` - Reconstruction space, the space to which the backprojection maps. - proj_space : `DiscretizedSpace` - Projection space, the space from which the backprojection maps. - """ - assert isinstance(geometry, Geometry) - assert isinstance(reco_space, DiscretizedSpace) - assert isinstance(proj_space, DiscretizedSpace) + # Copy result to host + if self.geometry.ndim == 2: + out[:] = self.proj_array + elif self.geometry.ndim == 3: + out[:] = np.swapaxes(self.proj_array, 0, 1).reshape( + self.proj_space.shape) - self.geometry = geometry - self.reco_space = reco_space - self.proj_space = proj_space - self.create_ids() + # Fix scaling to weight by pixel size + if (isinstance(self.geometry, Parallel2dGeometry) and + parse_version(ASTRA_VERSION) < parse_version('1.9.9.dev')): + # parallel2d scales with pixel stride + out *= 1 / float(self.geometry.det_partition.cell_volume) - # Create a mutually exclusive lock so that two callers cant use the - # same shared resource at the same time. - self._mutex = Lock() + return out - def call_backward(self, proj_data, out=None): + def call_backward(self, proj_data, out, **kwargs): """Run an ASTRA back-projection on the given data using the GPU. Parameters @@ -248,6 +216,7 @@ def call_backward(self, proj_data, out=None): """ with self._mutex: assert proj_data in self.proj_space + if out is not None: assert out in self.reco_space else: @@ -264,10 +233,10 @@ def call_backward(self, proj_data, out=None): astra.data3d.store(self.sino_id, swapped_proj_data) # Run algorithm - astra.algorithm.run(self.algo_id) + astra.algorithm.run(self.algo_backward_id) # Copy result to CPU memory - out[:] = self.out_array + out[:] = self.vol_array # Fix scaling to weight by pixel/voxel size out *= astra_cuda_bp_scaling_factor( @@ -275,56 +244,6 @@ def call_backward(self, proj_data, out=None): return out - def create_ids(self): - """Create ASTRA objects.""" - if self.geometry.motion_partition.ndim == 1: - motion_shape = self.geometry.motion_partition.shape - else: - # Need to flatten 2- or 3-dimensional angles into one axis - motion_shape = (np.prod(self.geometry.motion_partition.shape),) - - proj_shape = motion_shape + self.geometry.det_partition.shape - proj_ndim = len(proj_shape) - - if proj_ndim == 2: - astra_proj_shape = proj_shape - astra_vol_shape = self.reco_space.shape - elif proj_ndim == 3: - # The `u` and `v` axes of the projection data are swapped, - # see explanation in `astra_*_3d_geom_to_vec`. - astra_proj_shape = (proj_shape[1], proj_shape[0], proj_shape[2]) - astra_vol_shape = self.reco_space.shape - - self.in_array = np.empty(astra_proj_shape, - dtype='float32', order='C') - self.out_array = np.empty(astra_vol_shape, - dtype='float32', order='C') - - # Create ASTRA data structures - vol_geom = astra_volume_geometry(self.reco_space) - proj_geom = astra_projection_geometry(self.geometry) - self.sino_id = astra_data(proj_geom, - datatype='projection', - data=self.in_array, - ndim=proj_ndim, - allow_copy=False) - - proj_type = 'cuda' if proj_ndim == 2 else 'cuda3d' - self.proj_id = astra_projector( - proj_type, vol_geom, proj_geom, proj_ndim - ) - - self.vol_id = astra_data(vol_geom, - datatype='volume', - data=self.out_array, - ndim=self.reco_space.ndim, - allow_copy=False) - - # Create algorithm - self.algo_id = astra_algorithm( - 'backward', proj_ndim, self.vol_id, self.sino_id, - proj_id=self.proj_id, impl='cuda') - def __del__(self): """Delete ASTRA objects.""" if self.geometry.ndim == 2: @@ -332,9 +251,12 @@ def __del__(self): else: adata, aproj = astra.data3d, astra.projector3d - if self.algo_id is not None: - astra.algorithm.delete(self.algo_id) - self.algo_id = None + if self.algo_forward_id is not None: + astra.algorithm.delete(self.algo_forward_id) + self.algo_forward_id = None + if self.algo_backward_id is not None: + astra.algorithm.delete(self.algo_backward_id) + self.algo_backward_id = None if self.vol_id is not None: adata.delete(self.vol_id) self.vol_id = None @@ -478,4 +400,5 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry): if __name__ == '__main__': from odl.util.testutils import run_doctests + run_doctests() diff --git a/odl/tomo/backends/impl.py b/odl/tomo/backends/impl.py new file mode 100644 index 00000000000..34ba5996039 --- /dev/null +++ b/odl/tomo/backends/impl.py @@ -0,0 +1,74 @@ +# Copyright 2014-2020 The ODL contributors +# +# This file is part of ODL. +# +# This Source Code Form is subject to the terms of the Mozilla Public License, +# v. 2.0. If a copy of the MPL was not distributed with this file, You can +# obtain one at https://mozilla.org/MPL/2.0/. + +"""Ray transforms.""" + +from abc import abstractmethod, ABC +from __future__ import absolute_import, division, print_function + +from odl.discr import DiscreteLp +from odl.tomo.geometry import Geometry + +__all__ = ('RayTransformImplBase',) + + +class RayTransformImplBase(ABC): + """Base for a RayTransform implementation (a backend)""" + + def __init__(self, geometry, reco_space, proj_space): + """Initialize a new instance. + + Parameters + ---------- + geometry : `Geometry` + Geometry defining the tomographic setup. + reco_space : `DiscreteLp` + Reconstruction space, the space of the images to be forward + projected. + proj_space : `DiscreteLp` + Projection space, the space of the result. + """ + if not isinstance(geometry, Geometry): + raise TypeError('`geometry` must be a `Geometry` instance, got ' + '{!r}'.format(geometry)) + + if not isinstance(reco_space, DiscreteLp): + raise TypeError('`reco_space` must be a `DiscreteLP` instance, got ' + '{!r}'.format(reco_space)) + + if not isinstance(proj_space, DiscreteLp): + raise TypeError('`proj_space` must be a `DiscreteLP` instance, got ' + '{!r}'.format(proj_space)) + + self.geometry = geometry + self.reco_space = reco_space + self.proj_space = proj_space + + @staticmethod + def can_handle_size(size): + """A very general way to check if the implementation is capable + handling reconstruction volumes of a given size.""" + return True # by default no assumptions are made + + @classmethod + def supports_geometry(cls, geom): + """Check if the implementation can handle this geometry.""" + return True + + @classmethod + def supports_reco_space(cls, reco_name, reco_space): + """Check if the implementation can handle the reconstruction space.""" + return True + + @abstractmethod + def call_forward(self, x_real, out_real, **kwargs): + raise NotImplementedError('Needs to be implemented by the subclass.') + + @abstractmethod + def call_backward(self, x_real, out_real, **kwargs): + raise NotImplementedError('Needs to be implemented by the subclass.') diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index 86d5d64ac67..4d3a817ccc5 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -11,19 +11,26 @@ from __future__ import division import numpy as np +from odl.tomo import Parallel2dGeometry from odl.discr import uniform_discr_frompartition, uniform_partition from odl.discr.discr_utils import linear_interpolator, point_collocation +from odl.tomo.backends import RayTransformImplBase from odl.util.utility import writable_array try: import skimage + SKIMAGE_AVAILABLE = True except ImportError: SKIMAGE_AVAILABLE = False -__all__ = ('skimage_radon_forward_projector', 'skimage_radon_back_projector', - 'SKIMAGE_AVAILABLE') +__all__ = ( + 'SKIMAGE_AVAILABLE', + 'skimage_radon_forward_projector', + 'skimage_radon_back_projector', + 'SkimageRayTransformImpl', +) def skimage_proj_space(geometry, volume_space, proj_space): @@ -168,3 +175,47 @@ def skimage_radon_back_projector(sinogram, geometry, vol_space, out=None): out *= scaling_factor return out + + +class SkimageRayTransformImpl(RayTransformImplBase): + def __init__(self, geometry, reco_space, proj_space): + super().__init__(geometry, reco_space, proj_space) + + @staticmethod + def can_handle_size(size): + return not size >= 256 ** 2 + + @classmethod + def supports_geometry(cls, geometry): + if not isinstance(geometry, Parallel2dGeometry): + return TypeError("{!r} backend only supports 2d parallel " + 'geometries'.format(cls.__name__)) + + return True + + @classmethod + def supports_reco_space(cls, reco_name, reco_space): + mid_pt = reco_space.domain.mid_pt + if not np.allclose(mid_pt, [0, 0]): + return ValueError('`{}` must be centered at (0, 0), got ' + 'midpoint {}'.format(reco_name, mid_pt)) + + shape = reco_space.shape + if shape[0] != shape[1]: + return ValueError('`{}.shape` must have equal entries, ' + 'got {}'.format(reco_name, shape)) + + extent = reco_space.domain.extent + if extent[0] != extent[1]: + return ValueError('`{}.extent` must have equal entries, ' + 'got {}'.format(reco_name, extent)) + + return True + + def call_forward(self, x_real, out_real, **kwargs): + return skimage_radon_forward_projector( + x_real, self.geometry, self.proj_space.real_space, out_real) + + def call_backward(self, x_real, out_real, **kwargs): + return skimage_radon_back_projector( + x_real, self.geometry, self.reco_space.real_space, out_real) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 6e41e293637..631d6952cc4 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -18,30 +18,30 @@ from odl.operator import Operator from odl.space.weighting import ConstWeighting from odl.tomo.backends import ( - ASTRA_AVAILABLE, ASTRA_CUDA_AVAILABLE, ASTRA_VERSION, SKIMAGE_AVAILABLE, - AstraCudaBackProjectorImpl, AstraCudaProjectorImpl, - astra_cpu_back_projector, astra_cpu_forward_projector, astra_supports, - astra_versions_supporting, skimage_radon_back_projector, - skimage_radon_forward_projector) -from odl.tomo.geometry import ( - Geometry, Parallel2dGeometry, Parallel3dAxisGeometry) - -ASTRA_CPU_AVAILABLE = ASTRA_AVAILABLE -_SUPPORTED_IMPL = ('astra_cpu', 'astra_cuda', 'skimage') + ASTRA_AVAILABLE, ASTRA_CUDA_AVAILABLE, SKIMAGE_AVAILABLE, + RayTransformImplBase, + SkimageRayTransformImpl, AstraCudaRayTransformImpl, AstraCpuRayTransformImpl) +from odl.tomo.geometry import Geometry + +# Backends that are implemented in ODL and can be specified via `impl` +_SUPPORTED_IMPL = { + 'astra_cpu': AstraCpuRayTransformImpl, + 'astra_cuda': AstraCudaRayTransformImpl, + 'skimage': SkimageRayTransformImpl} + +# Backends that are installed, ordered by preference _AVAILABLE_IMPLS = [] -if ASTRA_CPU_AVAILABLE: - _AVAILABLE_IMPLS.append('astra_cpu') if ASTRA_CUDA_AVAILABLE: _AVAILABLE_IMPLS.append('astra_cuda') +if ASTRA_AVAILABLE: + _AVAILABLE_IMPLS.append('astra_cpu') if SKIMAGE_AVAILABLE: _AVAILABLE_IMPLS.append('skimage') - __all__ = ('RayTransform', 'RayBackProjection') class RayTransformBase(Operator): - """Base class for ray transforms containing common attributes.""" def __init__(self, reco_space, geometry, variant, **kwargs): @@ -111,106 +111,6 @@ def __init__(self, reco_space, geometry, variant, **kwargs): raise TypeError('`geometry` must be a `Geometry` instance, got ' '{!r}'.format(geometry)) - # Handle backend choice - if not _AVAILABLE_IMPLS: - raise RuntimeError('no ray transform back-end available; ' - 'this requires 3rd party packages, please ' - 'check the install docs') - impl = kwargs.pop('impl', None) - if impl is None: - # Select fastest available - if ASTRA_CUDA_AVAILABLE: - impl = 'astra_cuda' - elif ASTRA_AVAILABLE: - impl = 'astra_cpu' - if reco_space.size >= 512 ** 2: - warnings.warn( - "The best available backend ('astra_cpu') may be too " - "slow for volumes of this size. Consider using " - "'astra_cuda' if your machine has an Nvidia GPU. " - "This warning can be disabled by explicitly setting " - "`impl='astra_cpu'`.", - RuntimeWarning) - elif SKIMAGE_AVAILABLE: - impl = 'skimage' - if reco_space.size >= 256 ** 2: - warnings.warn( - "The best available backend ('skimage') may be too " - "slow for volumes of this size. Consider using ASTRA. " - "This warning can be disabled by explicitly setting " - "`impl='skimage'`.", - RuntimeWarning) - else: - raise RuntimeError('no backend') - - impl, impl_in = str(impl).lower(), impl - if impl not in _SUPPORTED_IMPL: - raise ValueError('`impl` {!r} not understood'.format(impl_in)) - if impl not in _AVAILABLE_IMPLS: - raise ValueError('{!r} back-end not available'.format(impl)) - - # Cache for input/output arrays of transforms - self.use_cache = kwargs.pop('use_cache', True) - - # Sanity checks - if impl.startswith('astra'): - if geometry.ndim > 2 and impl.endswith('cpu'): - raise ValueError('`impl` {!r} only works for 2d' - ''.format(impl_in)) - - # Print a warning if the detector midpoint normal vector at any - # angle is perpendicular to the geometry axis in parallel 3d - # single-axis geometry -- this is broken in some ASTRA versions - if ( - isinstance(geometry, Parallel3dAxisGeometry) and - not astra_supports('par3d_det_mid_pt_perp_to_axis') - ): - req_ver = astra_versions_supporting( - 'par3d_det_mid_pt_perp_to_axis' - ) - axis = geometry.axis - mid_pt = geometry.det_params.mid_pt - for i, angle in enumerate(geometry.angles): - if abs(np.dot(axis, - geometry.det_to_src(angle, mid_pt))) < 1e-4: - warnings.warn( - 'angle {}: detector midpoint normal {} is ' - 'perpendicular to the geometry axis {} in ' - '`Parallel3dAxisGeometry`; this is broken in ' - 'ASTRA {}, please upgrade to ASTRA {}' - ''.format(i, geometry.det_to_src(angle, mid_pt), - axis, ASTRA_VERSION, req_ver), - RuntimeWarning) - break - - elif impl == 'skimage': - if not isinstance(geometry, Parallel2dGeometry): - raise TypeError("{!r} backend only supports 2d parallel " - 'geometries'.format(impl)) - - mid_pt = reco_space.domain.mid_pt - if not np.allclose(mid_pt, [0, 0]): - raise ValueError('`{}` must be centered at (0, 0), got ' - 'midpoint {}'.format(reco_name, mid_pt)) - - shape = reco_space.shape - if shape[0] != shape[1]: - raise ValueError('`{}.shape` must have equal entries, ' - 'got {}'.format(reco_name, shape)) - - extent = reco_space.domain.extent - if extent[0] != extent[1]: - raise ValueError('`{}.extent` must have equal entries, ' - 'got {}'.format(reco_name, extent)) - - if reco_space.ndim != geometry.ndim: - raise ValueError('`{}.ndim` not equal to `geometry.ndim`: ' - '{} != {}'.format(reco_name, reco_space.ndim, - geometry.ndim)) - - self.__geometry = geometry - self.__impl = impl - # Generate or check projection space proj_space = kwargs.pop('proj_space', None) if proj_space is None: @@ -285,9 +185,91 @@ def __init__(self, reco_space, geometry, variant, **kwargs): proj_space.dtype, reco_space.dtype)) + if reco_space.ndim != geometry.ndim: + raise ValueError('`{}.ndim` not equal to `geometry.ndim`: ' + '{} != {}'.format(reco_name, reco_space.ndim, + geometry.ndim)) + + # Cache for input/output arrays of transforms + self.use_cache = kwargs.pop('use_cache', True) + + impl = kwargs.pop('impl', None) + if impl is None: # User didn't specify a backend + if not _AVAILABLE_IMPLS: + raise RuntimeError('no ray transform back-end available; ' + 'this requires 3rd party packages, please ' + 'check the install docs') + + # Select fastest available, _AVAILABLE_IMPLS is sorted by speed + impl_type = _SUPPORTED_IMPL[_AVAILABLE_IMPLS[0]] + + # Warn if implementation is not fast enough + if not impl_type.can_handle_size(reco_space.size): + if impl == 'astra_cpu': + warnings.warn( + "The best available backend ('astra_cpu') may be too " + "slow for volumes of this size. Consider using " + "'astra_cuda' if your machine has an Nvidia GPU. " + "This warning can be disabled by explicitly setting " + "`impl='astra_cpu'`.", + RuntimeWarning) + elif impl == 'skimage': + warnings.warn( + "The best available backend ('skimage') may be too " + "slow for volumes of this size. Consider using ASTRA. " + "This warning can be disabled by explicitly setting " + "`impl='skimage'`.", + RuntimeWarning) + else: + # User did specify `impl` + if isinstance(impl, str): + impl, impl_in = str(impl).lower(), impl + + if impl not in _SUPPORTED_IMPL: + raise ValueError('`impl` {!r} not understood'.format(impl_in)) + + if impl not in _AVAILABLE_IMPLS: + raise ValueError('{!r} back-end not available'.format(impl)) + + impl_type = _SUPPORTED_IMPL[impl] + elif isinstance(impl, type): + # User gave the type and leaves instantiation to us + if not issubclass(impl, RayTransformImplBase): + raise ValueError('Type {!r} must be a subclass of ' + '`RayTransformImplBase`.'.format(impl)) + + impl_type = impl + elif isinstance(impl, RayTransformImplBase): + # User gave an object for `impl` + # This is a possible use case, but would clutter this class + raise NotImplementedError( + 'Given for `impl` is an already initialized object of subclass' + ' `RayTransformImplBase`. This is not supported, because' + ' the object cannot be regenerated when `use_cache` is `False`.' + 'Instead give for `impl` the type of your implementation.') + else: + raise ValueError( + 'Given `impl` should be a `str` or the type of a subclass of ' + '`RayTransformImplBase`, now it is a {!r}.'.format(type(impl))) + + # Sanity checks + geometry_support = impl_type.supports_geometry(geometry) + if not geometry_support: + raise geometry_support + + reco_space_support = impl_type.supports_reco_space(reco_space, reco_name) + if not reco_space_support: + raise reco_space_support + + self.__geometry = geometry + + self.impl_type = impl_type + self.__cached_impl = None + self.__reco_space = reco_space + self.__proj_space = proj_space + # Reserve name for cached properties (used for efficiency reasons) self._adjoint = None - self._astra_wrapper = None # Extra kwargs that can be reused for adjoint etc. These must # be retrieved with `get` instead of `pop` above. @@ -304,7 +286,16 @@ def __init__(self, reco_space, geometry, variant, **kwargs): @property def impl(self): """Implementation back-end for the evaluation of this operator.""" - return self.__impl + + # Only skip impl creation (__cached_impl) when `use_cache` is True + if not self.use_cache or self.__cached_impl is None: + # Lazily (re)instantiate the backend + self.__cached_impl = self.impl_type( + self.geometry, + reco_space=self.__reco_space.real_space, + proj_space=self.__proj_space.real_space) + + return self.__cached_impl @property def geometry(self): @@ -338,7 +329,6 @@ def _call(self, x, out=None): class RayTransform(RayTransformBase): - """Discrete Ray transform between L^p spaces.""" def __init__(self, domain, geometry, **kwargs): @@ -396,41 +386,9 @@ def __init__(self, domain, geometry, **kwargs): variant='forward', **kwargs) def _call_real(self, x_real, out_real, **kwargs): - """Real-space forward projection for the current set-up. + """Real-space forward projection for the current set-up.""" - This method also sets ``self._astra_projector`` for - ``impl='astra_cuda'`` and enabled cache. - """ - if self.impl.startswith('astra'): - backend, data_impl = self.impl.split('_') - - if data_impl == 'cpu': - return astra_cpu_forward_projector( - x_real, self.geometry, self.range.real_space, out_real, - **kwargs) - - elif data_impl == 'cuda': - if self._astra_wrapper is None: - astra_wrapper = AstraCudaProjectorImpl( - self.geometry, self.domain.real_space, - self.range.real_space) - if self.use_cache: - self._astra_wrapper = astra_wrapper - else: - astra_wrapper = self._astra_wrapper - - return astra_wrapper.call_forward(x_real, out_real, **kwargs) - else: - # Should never happen - raise RuntimeError('bad `impl` {!r}'.format(self.impl)) - - elif self.impl == 'skimage': - return skimage_radon_forward_projector( - x_real, self.geometry, self.range.real_space, out_real, - **kwargs) - else: - # Should never happen - raise RuntimeError('bad `impl` {!r}'.format(self.impl)) + return self.impl.call_forward(x_real, out_real, **kwargs) @property def adjoint(self): @@ -444,7 +402,7 @@ def adjoint(self): kwargs = self._extra_kwargs.copy() kwargs['domain'] = self.range self._adjoint = RayBackProjection(self.domain, self.geometry, - impl=self.impl, + impl=self.impl_type, use_cache=self.use_cache, **kwargs) @@ -452,7 +410,6 @@ def adjoint(self): class RayBackProjection(RayTransformBase): - """Adjoint of the discrete Ray transform between L^p spaces.""" def __init__(self, range, geometry, **kwargs): @@ -511,39 +468,9 @@ def __init__(self, range, geometry, **kwargs): variant='backward', **kwargs) def _call_real(self, x_real, out_real, **kwargs): - """Real-space back-projection for the current set-up. + """Real-space back-projection for the current set-up.""" - This method also sets ``self._astra_backprojector`` for - ``impl='astra_cuda'`` and enabled cache. - """ - if self.impl.startswith('astra'): - backend, data_impl = self.impl.split('_') - if data_impl == 'cpu': - return astra_cpu_back_projector( - x_real, self.geometry, self.range.real_space, out_real, - **kwargs) - elif data_impl == 'cuda': - if self._astra_wrapper is None: - astra_wrapper = AstraCudaBackProjectorImpl( - self.geometry, self.range.real_space, - self.domain.real_space) - if self.use_cache: - self._astra_wrapper = astra_wrapper - else: - astra_wrapper = self._astra_wrapper - - return astra_wrapper.call_backward(x_real, out_real, **kwargs) - else: - # Should never happen - raise RuntimeError('bad `impl` {!r}'.format(self.impl)) - - elif self.impl == 'skimage': - return skimage_radon_back_projector( - x_real, self.geometry, self.range.real_space, out_real, - **kwargs) - else: - # Should never happen - raise RuntimeError('bad `impl` {!r}'.format(self.impl)) + return self.impl.call_backward(x_real, out_real, **kwargs) @property def adjoint(self): @@ -566,4 +493,5 @@ def adjoint(self): if __name__ == '__main__': from odl.util.testutils import run_doctests + run_doctests() From c4db486b270bbe388ca3d7f261ea4dacfa96f7d7 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Sun, 26 Jan 2020 20:41:09 +0100 Subject: [PATCH 03/39] move __future__ import to beginning --- odl/tomo/backends/impl.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/odl/tomo/backends/impl.py b/odl/tomo/backends/impl.py index 34ba5996039..0f13bf16f71 100644 --- a/odl/tomo/backends/impl.py +++ b/odl/tomo/backends/impl.py @@ -7,10 +7,9 @@ # obtain one at https://mozilla.org/MPL/2.0/. """Ray transforms.""" - -from abc import abstractmethod, ABC from __future__ import absolute_import, division, print_function +from abc import abstractmethod, ABC from odl.discr import DiscreteLp from odl.tomo.geometry import Geometry From b251fbdf7c9c2c6d60359b6251ad9a0ce2d3417e Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Mon, 27 Jan 2020 21:44:18 +0100 Subject: [PATCH 04/39] Keep interface of RayTransform.impl the same, fix lowercasing --- odl/tomo/operators/ray_trafo.py | 45 +++++++++++++++++++++------------ 1 file changed, 29 insertions(+), 16 deletions(-) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 631d6952cc4..296f1222f7c 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -194,6 +194,7 @@ def __init__(self, reco_space, geometry, variant, **kwargs): self.use_cache = kwargs.pop('use_cache', True) impl = kwargs.pop('impl', None) + if impl is None: # User didn't specify a backend if not _AVAILABLE_IMPLS: raise RuntimeError('no ray transform back-end available; ' @@ -205,7 +206,7 @@ def __init__(self, reco_space, geometry, variant, **kwargs): # Warn if implementation is not fast enough if not impl_type.can_handle_size(reco_space.size): - if impl == 'astra_cpu': + if impl_type == AstraCpuRayTransformImpl: warnings.warn( "The best available backend ('astra_cpu') may be too " "slow for volumes of this size. Consider using " @@ -213,25 +214,29 @@ def __init__(self, reco_space, geometry, variant, **kwargs): "This warning can be disabled by explicitly setting " "`impl='astra_cpu'`.", RuntimeWarning) - elif impl == 'skimage': + elif impl_type == SkimageRayTransformImpl: warnings.warn( "The best available backend ('skimage') may be too " "slow for volumes of this size. Consider using ASTRA. " "This warning can be disabled by explicitly setting " "`impl='skimage'`.", RuntimeWarning) + else: + warnings.warn( + "The `impl` backend indicates that it might be too slow " + " for volumes of the input size.", + RuntimeWarning) + else: # User did specify `impl` if isinstance(impl, str): - impl, impl_in = str(impl).lower(), impl + if impl.lower() not in _SUPPORTED_IMPL: + raise ValueError('`impl` {!r} not understood'.format(impl)) - if impl not in _SUPPORTED_IMPL: - raise ValueError('`impl` {!r} not understood'.format(impl_in)) - - if impl not in _AVAILABLE_IMPLS: + if impl.lower() not in _AVAILABLE_IMPLS: raise ValueError('{!r} back-end not available'.format(impl)) - impl_type = _SUPPORTED_IMPL[impl] + impl_type = _SUPPORTED_IMPL[impl.lower()] elif isinstance(impl, type): # User gave the type and leaves instantiation to us if not issubclass(impl, RayTransformImplBase): @@ -263,7 +268,8 @@ def __init__(self, reco_space, geometry, variant, **kwargs): self.__geometry = geometry - self.impl_type = impl_type + self._impl_type = impl_type + self.__impl = impl.lower() if isinstance(impl, str) else impl self.__cached_impl = None self.__reco_space = reco_space self.__proj_space = proj_space @@ -285,12 +291,17 @@ def __init__(self, reco_space, geometry, variant, **kwargs): @property def impl(self): - """Implementation back-end for the evaluation of this operator.""" + """Implementation name string or type.""" + + return self.__impl + + def create_impl(self, from_cache=True): + """Fetches or instantiates implementation backend for evaluation.""" - # Only skip impl creation (__cached_impl) when `use_cache` is True - if not self.use_cache or self.__cached_impl is None: + # Skip impl creation (__cached_impl) when `clear_cache` is True + if not from_cache or self.__cached_impl is None: # Lazily (re)instantiate the backend - self.__cached_impl = self.impl_type( + self.__cached_impl = self._impl_type( self.geometry, reco_space=self.__reco_space.real_space, proj_space=self.__proj_space.real_space) @@ -388,7 +399,8 @@ def __init__(self, domain, geometry, **kwargs): def _call_real(self, x_real, out_real, **kwargs): """Real-space forward projection for the current set-up.""" - return self.impl.call_forward(x_real, out_real, **kwargs) + return self.create_impl(self.use_cache) \ + .call_forward(x_real, out_real, **kwargs) @property def adjoint(self): @@ -402,7 +414,7 @@ def adjoint(self): kwargs = self._extra_kwargs.copy() kwargs['domain'] = self.range self._adjoint = RayBackProjection(self.domain, self.geometry, - impl=self.impl_type, + impl=self.impl, use_cache=self.use_cache, **kwargs) @@ -470,7 +482,8 @@ def __init__(self, range, geometry, **kwargs): def _call_real(self, x_real, out_real, **kwargs): """Real-space back-projection for the current set-up.""" - return self.impl.call_backward(x_real, out_real, **kwargs) + return self.create_impl(self.use_cache) \ + .call_backward(x_real, out_real, **kwargs) @property def adjoint(self): From 468f78ed617f8d5d33d97d74a52da7b0ff797f70 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Mon, 27 Jan 2020 21:57:14 +0100 Subject: [PATCH 05/39] Share cache with adjoint --- odl/tomo/operators/ray_trafo.py | 40 +++++++++++++++++---------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 296f1222f7c..e0c6dee7b52 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -194,6 +194,7 @@ def __init__(self, reco_space, geometry, variant, **kwargs): self.use_cache = kwargs.pop('use_cache', True) impl = kwargs.pop('impl', None) + self.__cached_impl = None if impl is None: # User didn't specify a backend if not _AVAILABLE_IMPLS: @@ -245,16 +246,13 @@ def __init__(self, reco_space, geometry, variant, **kwargs): impl_type = impl elif isinstance(impl, RayTransformImplBase): - # User gave an object for `impl` - # This is a possible use case, but would clutter this class - raise NotImplementedError( - 'Given for `impl` is an already initialized object of subclass' - ' `RayTransformImplBase`. This is not supported, because' - ' the object cannot be regenerated when `use_cache` is `False`.' - 'Instead give for `impl` the type of your implementation.') + # User gave an object for `impl`, meaning to set the + # backend cache to an already initiated object + impl_type = type(impl) + self.__cached_impl = impl else: - raise ValueError( - 'Given `impl` should be a `str` or the type of a subclass of ' + raise TypeError( + 'Given `impl` should be a `str`, or the type or subclass of ' '`RayTransformImplBase`, now it is a {!r}.'.format(type(impl))) # Sanity checks @@ -267,10 +265,8 @@ def __init__(self, reco_space, geometry, variant, **kwargs): raise reco_space_support self.__geometry = geometry - self._impl_type = impl_type self.__impl = impl.lower() if isinstance(impl, str) else impl - self.__cached_impl = None self.__reco_space = reco_space self.__proj_space = proj_space @@ -413,10 +409,13 @@ def adjoint(self): if self._adjoint is None: kwargs = self._extra_kwargs.copy() kwargs['domain'] = self.range - self._adjoint = RayBackProjection(self.domain, self.geometry, - impl=self.impl, - use_cache=self.use_cache, - **kwargs) + + # initiate adjoint with cached implementation if `use_cache` + self._adjoint = RayBackProjection( + self.domain, self.geometry, + impl=self.impl if not self.use_cache else self.create_impl(self.use_cache), + use_cache=self.use_cache, + **kwargs) return self._adjoint @@ -496,10 +495,13 @@ def adjoint(self): if self._adjoint is None: kwargs = self._extra_kwargs.copy() kwargs['range'] = self.domain - self._adjoint = RayTransform(self.range, self.geometry, - impl=self.impl, - use_cache=self.use_cache, - **kwargs) + + # initiate adjoint with cached implementation if `use_cache` + self._adjoint = RayTransform( + self.range, self.geometry, + impl=self.impl if not self.use_cache else self.create_impl(self.use_cache), + use_cache=self.use_cache, + **kwargs) return self._adjoint From 3a3036b4766c5cd659f8bc21bb34bf8b7fafb6b6 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Mon, 27 Jan 2020 22:37:54 +0100 Subject: [PATCH 06/39] Fix self.__impl initialization --- odl/tomo/operators/ray_trafo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index e0c6dee7b52..611862c8b6e 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -266,7 +266,7 @@ def __init__(self, reco_space, geometry, variant, **kwargs): self.__geometry = geometry self._impl_type = impl_type - self.__impl = impl.lower() if isinstance(impl, str) else impl + self.__impl = impl.lower() if isinstance(impl, str) else impl_type self.__reco_space = reco_space self.__proj_space = proj_space From 9b5025edaded939ba3033337f8480712a1f38389 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Tue, 28 Jan 2020 21:27:16 +0100 Subject: [PATCH 07/39] Style changes, change ValueError into TypeError --- odl/tomo/backends/astra_cpu.py | 6 ++++-- odl/tomo/backends/astra_cuda.py | 9 +++++--- odl/tomo/backends/impl.py | 8 +++---- odl/tomo/operators/ray_trafo.py | 38 ++++++++++++++++++++------------- 4 files changed, 37 insertions(+), 24 deletions(-) diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index cec1a53a267..9336a4610dd 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -273,11 +273,13 @@ def supports_geometry(cls, geometry): def call_forward(self, x_real, out_real, **kwargs): return astra_cpu_forward_projector( - x_real, self.geometry, self.proj_space.real_space, out_real, **kwargs) + x_real, self.geometry, self.proj_space.real_space, + out_real, **kwargs) def call_backward(self, x_real, out_real, **kwargs): return astra_cpu_back_projector( - x_real, self.geometry, self.reco_space.real_space, out_real, **kwargs) + x_real, self.geometry, self.reco_space.real_space, + out_real, **kwargs) if __name__ == '__main__': diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 4323ba7748d..1468cab10df 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -19,7 +19,8 @@ from odl.discr import DiscretizedSpace from odl.tomo.backends.astra_setup import ( ASTRA_VERSION, astra_algorithm, astra_data, astra_projection_geometry, - astra_projector, astra_volume_geometry, astra_supports, astra_versions_supporting) + astra_projector, astra_volume_geometry, astra_supports, + astra_versions_supporting) from odl.tomo.geometry import ( ConeBeamGeometry, FanBeamGeometry, Geometry, Parallel2dGeometry, Parallel3dAxisGeometry) @@ -74,7 +75,8 @@ def supports_geometry(cls, geometry): # single-axis geometry -- this is broken in some ASTRA versions if (isinstance(geometry, Parallel3dAxisGeometry) and not astra_supports('par3d_det_mid_pt_perp_to_axis')): - req_ver = astra_versions_supporting('par3d_det_mid_pt_perp_to_axis') + req_ver = astra_versions_supporting( + 'par3d_det_mid_pt_perp_to_axis') axis = geometry.axis mid_pt = geometry.det_params.mid_pt for i, angle in enumerate(geometry.angles): @@ -114,7 +116,8 @@ def create_ids(self): astra_vol_shape = self.reco_space.shape self.vol_array = np.empty(astra_vol_shape, dtype='float32', order='C') - self.proj_array = np.empty(astra_proj_shape, dtype='float32', order='C') + self.proj_array = np.empty(astra_proj_shape, dtype='float32', + order='C') # Create ASTRA data structures vol_geom = astra_volume_geometry(self.reco_space) diff --git a/odl/tomo/backends/impl.py b/odl/tomo/backends/impl.py index 0f13bf16f71..b3baacaf1a8 100644 --- a/odl/tomo/backends/impl.py +++ b/odl/tomo/backends/impl.py @@ -37,12 +37,12 @@ def __init__(self, geometry, reco_space, proj_space): '{!r}'.format(geometry)) if not isinstance(reco_space, DiscreteLp): - raise TypeError('`reco_space` must be a `DiscreteLP` instance, got ' - '{!r}'.format(reco_space)) + raise TypeError('`reco_space` must be a `DiscreteLP` instance, got' + ' {!r}'.format(reco_space)) if not isinstance(proj_space, DiscreteLp): - raise TypeError('`proj_space` must be a `DiscreteLP` instance, got ' - '{!r}'.format(proj_space)) + raise TypeError('`proj_space` must be a `DiscreteLP` instance, got' + ' {!r}'.format(proj_space)) self.geometry = geometry self.reco_space = reco_space diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 611862c8b6e..72d1e311f18 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -20,7 +20,8 @@ from odl.tomo.backends import ( ASTRA_AVAILABLE, ASTRA_CUDA_AVAILABLE, SKIMAGE_AVAILABLE, RayTransformImplBase, - SkimageRayTransformImpl, AstraCudaRayTransformImpl, AstraCpuRayTransformImpl) + SkimageRayTransformImpl, AstraCudaRayTransformImpl, + AstraCpuRayTransformImpl) from odl.tomo.geometry import Geometry # Backends that are implemented in ODL and can be specified via `impl` @@ -118,9 +119,9 @@ def __init__(self, reco_space, geometry, variant, **kwargs): if not reco_space.is_weighted: weighting = None - elif (isinstance(reco_space.weighting, ConstWeighting) and - np.isclose(reco_space.weighting.const, - reco_space.cell_volume)): + elif (isinstance(reco_space.weighting, ConstWeighting) + and np.isclose(reco_space.weighting.const, + reco_space.cell_volume)): # Approximate cell volume # TODO: find a way to treat angles and detector differently # regarding weighting. While the detector should be uniformly @@ -224,8 +225,8 @@ def __init__(self, reco_space, geometry, variant, **kwargs): RuntimeWarning) else: warnings.warn( - "The `impl` backend indicates that it might be too slow " - " for volumes of the input size.", + "The `impl` backend indicates that it might be too " + "slow for volumes of the input size.", RuntimeWarning) else: @@ -235,14 +236,15 @@ def __init__(self, reco_space, geometry, variant, **kwargs): raise ValueError('`impl` {!r} not understood'.format(impl)) if impl.lower() not in _AVAILABLE_IMPLS: - raise ValueError('{!r} back-end not available'.format(impl)) + raise ValueError( + '{!r} back-end not available'.format(impl)) impl_type = _SUPPORTED_IMPL[impl.lower()] elif isinstance(impl, type): # User gave the type and leaves instantiation to us if not issubclass(impl, RayTransformImplBase): - raise ValueError('Type {!r} must be a subclass of ' - '`RayTransformImplBase`.'.format(impl)) + raise TypeError('Type {!r} must be a subclass of ' + '`RayTransformImplBase`.'.format(impl)) impl_type = impl elif isinstance(impl, RayTransformImplBase): @@ -252,15 +254,17 @@ def __init__(self, reco_space, geometry, variant, **kwargs): self.__cached_impl = impl else: raise TypeError( - 'Given `impl` should be a `str`, or the type or subclass of ' - '`RayTransformImplBase`, now it is a {!r}.'.format(type(impl))) + 'Given `impl` should be a `str`, or the type or subclass ' + 'of `RayTransformImplBase`, ' + 'now it is a {!r}.'.format(type(impl))) # Sanity checks geometry_support = impl_type.supports_geometry(geometry) if not geometry_support: raise geometry_support - reco_space_support = impl_type.supports_reco_space(reco_space, reco_name) + reco_space_support = impl_type.supports_reco_space(reco_space, + reco_name) if not reco_space_support: raise reco_space_support @@ -413,7 +417,9 @@ def adjoint(self): # initiate adjoint with cached implementation if `use_cache` self._adjoint = RayBackProjection( self.domain, self.geometry, - impl=self.impl if not self.use_cache else self.create_impl(self.use_cache), + impl=(self.impl + if not self.use_cache + else self.create_impl(self.use_cache)), use_cache=self.use_cache, **kwargs) @@ -499,11 +505,13 @@ def adjoint(self): # initiate adjoint with cached implementation if `use_cache` self._adjoint = RayTransform( self.range, self.geometry, - impl=self.impl if not self.use_cache else self.create_impl(self.use_cache), + impl=(self.impl + if not self.use_cache + else self.create_impl(self.use_cache)), use_cache=self.use_cache, **kwargs) - return self._adjoint + return self._adjoint if __name__ == '__main__': From 894cb073d2bed2f5a79841488fd5fad26cd2aa2b Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Sun, 16 Feb 2020 12:57:00 +0100 Subject: [PATCH 08/39] Remove RayTransformImplBase --- odl/tomo/backends/astra_cpu.py | 52 +++++++++++++++----- odl/tomo/backends/astra_cuda.py | 29 +++++++---- odl/tomo/backends/impl.py | 73 ---------------------------- odl/tomo/backends/skimage_radon.py | 78 ++++++++++++++++++++---------- odl/tomo/operators/ray_trafo.py | 77 +++++++++-------------------- 5 files changed, 134 insertions(+), 175 deletions(-) delete mode 100644 odl/tomo/backends/impl.py diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index 9336a4610dd..2c1fb6105cb 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -11,12 +11,12 @@ from __future__ import absolute_import, division, print_function import numpy as np +import warnings from odl.discr import DiscretizedSpace, DiscretizedSpaceElement from odl.tomo.backends.astra_setup import ( astra_algorithm, astra_data, astra_projection_geometry, astra_projector, astra_volume_geometry) -from odl.tomo.backends.impl import RayTransformImplBase from odl.tomo.geometry import ( DivergentBeamGeometry, Geometry, ParallelBeamGeometry) from odl.util import writable_array @@ -30,7 +30,7 @@ 'astra_cpu_forward_projector', 'astra_cpu_back_projector', 'default_astra_proj_type', - 'AstraCpuRayTransformImpl', + 'AstraCpu', ) @@ -255,21 +255,47 @@ def astra_cpu_back_projector(proj_data, geometry, vol_space, out=None, return out -class AstraCpuRayTransformImpl(RayTransformImplBase): +class AstraCpu: def __init__(self, geometry, reco_space, proj_space): - super().__init__(geometry, reco_space, proj_space) + """Initialize a new instance. + + Parameters + ---------- + geometry : `Geometry` + Geometry defining the tomographic setup. + reco_space : `DiscreteLp` + Reconstruction space, the space of the images to be forward + projected. + proj_space : `DiscreteLp` + Projection space, the space of the result. + """ + if not isinstance(geometry, Geometry): + raise TypeError('`geometry` must be a `Geometry` instance, got ' + '{!r}'.format(geometry)) + + if not isinstance(reco_space, DiscreteLp): + raise TypeError('`reco_space` must be a `DiscreteLP` instance, got' + ' {!r}'.format(reco_space)) + + if not isinstance(proj_space, DiscreteLp): + raise TypeError('`proj_space` must be a `DiscreteLP` instance, got' + ' {!r}'.format(proj_space)) + + if reco_space.size >= 512 ** 2: + warnings.warn( + "The 'astra_cpu' backend may be too slow for volumes of this " + "size. Consider using 'astra_cuda' if your machine has an " + "Nvidia GPU. This warning can be disabled by explicitly " + "setting `impl='astra_cpu'`.", + RuntimeWarning) - @staticmethod - def can_handle_size(size): - return not size >= 512 ** 2 - - @classmethod - def supports_geometry(cls, geometry): if geometry.ndim > 2: - return ValueError('`impl` {!r} only works for 2d' - ''.format(cls.__name__)) + raise ValueError('`impl` {!r} only works for 2d' + ''.format(self.__name__)) - return True + self.geometry = geometry + self.reco_space = reco_space + self.proj_space = proj_space def call_forward(self, x_real, out_real, **kwargs): return astra_cpu_forward_projector( diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 1468cab10df..4cc5c8ef9d6 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -24,7 +24,6 @@ from odl.tomo.geometry import ( ConeBeamGeometry, FanBeamGeometry, Geometry, Parallel2dGeometry, Parallel3dAxisGeometry) -from odl.tomo.backends.impl import RayTransformImplBase try: import astra @@ -35,11 +34,11 @@ __all__ = ( 'ASTRA_CUDA_AVAILABLE', - 'AstraCudaRayTransformImpl', + 'AstraCuda', ) -class AstraCudaRayTransformImpl(RayTransformImplBase): +class AstraCuda: """Thin wrapper around ASTRA.""" algo_forward_id = None @@ -61,15 +60,18 @@ def __init__(self, geometry, reco_space, proj_space): proj_space : `DiscretizedSpace` Projection space, the space of the result. """ - super().__init__(geometry, reco_space, proj_space) + if not isinstance(geometry, Geometry): + raise TypeError('`geometry` must be a `Geometry` instance, got ' + '{!r}'.format(geometry)) - self.create_ids() + if not isinstance(reco_space, DiscreteLp): + raise TypeError('`reco_space` must be a `DiscreteLP` instance, got' + ' {!r}'.format(reco_space)) - # ASTRA projectors are not thread-safe, thus we need to lock ourselves - self._mutex = Lock() + if not isinstance(proj_space, DiscreteLp): + raise TypeError('`proj_space` must be a `DiscreteLP` instance, got' + ' {!r}'.format(proj_space)) - @classmethod - def supports_geometry(cls, geometry): # Print a warning if the detector midpoint normal vector at any # angle is perpendicular to the geometry axis in parallel 3d # single-axis geometry -- this is broken in some ASTRA versions @@ -92,7 +94,14 @@ def supports_geometry(cls, geometry): RuntimeWarning) break - return True + self.geometry = geometry + self.reco_space = reco_space + self.proj_space = proj_space + + self.create_ids() + + # ASTRA projectors are not thread-safe, thus we need to lock ourselves + self._mutex = Lock() def create_ids(self): """Create ASTRA objects.""" diff --git a/odl/tomo/backends/impl.py b/odl/tomo/backends/impl.py deleted file mode 100644 index b3baacaf1a8..00000000000 --- a/odl/tomo/backends/impl.py +++ /dev/null @@ -1,73 +0,0 @@ -# Copyright 2014-2020 The ODL contributors -# -# This file is part of ODL. -# -# This Source Code Form is subject to the terms of the Mozilla Public License, -# v. 2.0. If a copy of the MPL was not distributed with this file, You can -# obtain one at https://mozilla.org/MPL/2.0/. - -"""Ray transforms.""" -from __future__ import absolute_import, division, print_function - -from abc import abstractmethod, ABC -from odl.discr import DiscreteLp -from odl.tomo.geometry import Geometry - -__all__ = ('RayTransformImplBase',) - - -class RayTransformImplBase(ABC): - """Base for a RayTransform implementation (a backend)""" - - def __init__(self, geometry, reco_space, proj_space): - """Initialize a new instance. - - Parameters - ---------- - geometry : `Geometry` - Geometry defining the tomographic setup. - reco_space : `DiscreteLp` - Reconstruction space, the space of the images to be forward - projected. - proj_space : `DiscreteLp` - Projection space, the space of the result. - """ - if not isinstance(geometry, Geometry): - raise TypeError('`geometry` must be a `Geometry` instance, got ' - '{!r}'.format(geometry)) - - if not isinstance(reco_space, DiscreteLp): - raise TypeError('`reco_space` must be a `DiscreteLP` instance, got' - ' {!r}'.format(reco_space)) - - if not isinstance(proj_space, DiscreteLp): - raise TypeError('`proj_space` must be a `DiscreteLP` instance, got' - ' {!r}'.format(proj_space)) - - self.geometry = geometry - self.reco_space = reco_space - self.proj_space = proj_space - - @staticmethod - def can_handle_size(size): - """A very general way to check if the implementation is capable - handling reconstruction volumes of a given size.""" - return True # by default no assumptions are made - - @classmethod - def supports_geometry(cls, geom): - """Check if the implementation can handle this geometry.""" - return True - - @classmethod - def supports_reco_space(cls, reco_name, reco_space): - """Check if the implementation can handle the reconstruction space.""" - return True - - @abstractmethod - def call_forward(self, x_real, out_real, **kwargs): - raise NotImplementedError('Needs to be implemented by the subclass.') - - @abstractmethod - def call_backward(self, x_real, out_real, **kwargs): - raise NotImplementedError('Needs to be implemented by the subclass.') diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index 4d3a817ccc5..94173311d50 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -11,11 +11,13 @@ from __future__ import division import numpy as np -from odl.tomo import Parallel2dGeometry +import warnings -from odl.discr import uniform_discr_frompartition, uniform_partition +from odl.tomo import Parallel2dGeometry, Geometry + +from odl.discr import uniform_discr_frompartition, uniform_partition, \ + DiscretizedSpace from odl.discr.discr_utils import linear_interpolator, point_collocation -from odl.tomo.backends import RayTransformImplBase from odl.util.utility import writable_array try: @@ -29,7 +31,7 @@ 'SKIMAGE_AVAILABLE', 'skimage_radon_forward_projector', 'skimage_radon_back_projector', - 'SkimageRayTransformImpl', + 'Skimage', ) @@ -177,40 +179,66 @@ def skimage_radon_back_projector(sinogram, geometry, vol_space, out=None): return out -class SkimageRayTransformImpl(RayTransformImplBase): +class Skimage: def __init__(self, geometry, reco_space, proj_space): - super().__init__(geometry, reco_space, proj_space) - - @staticmethod - def can_handle_size(size): - return not size >= 256 ** 2 + """Initialize a new instance. + + Parameters + ---------- + geometry : `Geometry` + Geometry defining the tomographic setup. + reco_space : `DiscreteLp` + Reconstruction space, the space of the images to be forward + projected. + proj_space : `DiscreteLp` + Projection space, the space of the result. + """ + if not isinstance(geometry, Geometry): + raise TypeError('`geometry` must be a `Geometry` instance, got ' + '{!r}'.format(geometry)) + + if not isinstance(reco_space, DiscretizedSpace): + raise TypeError( + '`reco_space` must be a `DiscretizedSpace` instance, got {!r}' + ''.format(reco_space) + ) + + if not isinstance(proj_space, DiscretizedSpace): + raise TypeError( + '`proj_space` must be a `DiscretizedSpace` instance, got {!r}' + ''.format(proj_space) + ) + + if reco_space.size >= 256 ** 2: + warnings.warn( + "The 'skimage' backend may be too slow for volumes of this " + "size. Consider using 'astra_cpu', or 'astra_cuda' if your " + "machine has an Nvidia GPU. This warning can be disabled by " + "explicitly setting `impl='astra_cpu'`.", + RuntimeWarning) - @classmethod - def supports_geometry(cls, geometry): if not isinstance(geometry, Parallel2dGeometry): - return TypeError("{!r} backend only supports 2d parallel " - 'geometries'.format(cls.__name__)) - - return True + raise TypeError("{!r} backend only supports 2d parallel " + 'geometries'.format(self.__name__)) - @classmethod - def supports_reco_space(cls, reco_name, reco_space): mid_pt = reco_space.domain.mid_pt if not np.allclose(mid_pt, [0, 0]): - return ValueError('`{}` must be centered at (0, 0), got ' - 'midpoint {}'.format(reco_name, mid_pt)) + raise ValueError('Reconstruction space must be centered at (0, 0),' + 'got midpoint {}'.format(mid_pt)) shape = reco_space.shape if shape[0] != shape[1]: - return ValueError('`{}.shape` must have equal entries, ' - 'got {}'.format(reco_name, shape)) + raise ValueError('`reco_space.shape` must have equal entries, ' + 'got {}'.format(shape)) extent = reco_space.domain.extent if extent[0] != extent[1]: - return ValueError('`{}.extent` must have equal entries, ' - 'got {}'.format(reco_name, extent)) + raise ValueError('`reco_space.extent` must have equal entries, ' + 'got {}'.format(extent)) - return True + self.geometry = geometry + self.reco_space = reco_space + self.proj_space = proj_space def call_forward(self, x_real, out_real, **kwargs): return skimage_radon_forward_projector( diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 72d1e311f18..3aafe72ee43 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -19,16 +19,14 @@ from odl.space.weighting import ConstWeighting from odl.tomo.backends import ( ASTRA_AVAILABLE, ASTRA_CUDA_AVAILABLE, SKIMAGE_AVAILABLE, - RayTransformImplBase, - SkimageRayTransformImpl, AstraCudaRayTransformImpl, - AstraCpuRayTransformImpl) + Skimage, AstraCuda, AstraCpu) from odl.tomo.geometry import Geometry # Backends that are implemented in ODL and can be specified via `impl` _SUPPORTED_IMPL = { - 'astra_cpu': AstraCpuRayTransformImpl, - 'astra_cuda': AstraCudaRayTransformImpl, - 'skimage': SkimageRayTransformImpl} + 'astra_cpu': AstraCpu, + 'astra_cuda': AstraCuda, + 'skimage': Skimage} # Backends that are installed, ordered by preference _AVAILABLE_IMPLS = [] @@ -206,29 +204,6 @@ def __init__(self, reco_space, geometry, variant, **kwargs): # Select fastest available, _AVAILABLE_IMPLS is sorted by speed impl_type = _SUPPORTED_IMPL[_AVAILABLE_IMPLS[0]] - # Warn if implementation is not fast enough - if not impl_type.can_handle_size(reco_space.size): - if impl_type == AstraCpuRayTransformImpl: - warnings.warn( - "The best available backend ('astra_cpu') may be too " - "slow for volumes of this size. Consider using " - "'astra_cuda' if your machine has an Nvidia GPU. " - "This warning can be disabled by explicitly setting " - "`impl='astra_cpu'`.", - RuntimeWarning) - elif impl_type == SkimageRayTransformImpl: - warnings.warn( - "The best available backend ('skimage') may be too " - "slow for volumes of this size. Consider using ASTRA. " - "This warning can be disabled by explicitly setting " - "`impl='skimage'`.", - RuntimeWarning) - else: - warnings.warn( - "The `impl` backend indicates that it might be too " - "slow for volumes of the input size.", - RuntimeWarning) - else: # User did specify `impl` if isinstance(impl, str): @@ -240,33 +215,27 @@ def __init__(self, reco_space, geometry, variant, **kwargs): '{!r} back-end not available'.format(impl)) impl_type = _SUPPORTED_IMPL[impl.lower()] - elif isinstance(impl, type): + elif isinstance(impl, type) or isinstance(impl, object): # User gave the type and leaves instantiation to us - if not issubclass(impl, RayTransformImplBase): - raise TypeError('Type {!r} must be a subclass of ' - '`RayTransformImplBase`.'.format(impl)) - - impl_type = impl - elif isinstance(impl, RayTransformImplBase): - # User gave an object for `impl`, meaning to set the - # backend cache to an already initiated object - impl_type = type(impl) - self.__cached_impl = impl - else: - raise TypeError( - 'Given `impl` should be a `str`, or the type or subclass ' - 'of `RayTransformImplBase`, ' - 'now it is a {!r}.'.format(type(impl))) + forward = getattr(impl, "call_forward", None) + backward = getattr(impl, "call_backward", None) - # Sanity checks - geometry_support = impl_type.supports_geometry(geometry) - if not geometry_support: - raise geometry_support + if not callable(forward) and not callable(backward): + raise TypeError('Type {!r} must be have a `call_forward` ' + 'or `call_backward`.'.format(impl)) - reco_space_support = impl_type.supports_reco_space(reco_space, - reco_name) - if not reco_space_support: - raise reco_space_support + if isinstance(impl, type): + impl_type = impl + else: + # User gave an object for `impl`, meaning to set the + # backend cache to an already initiated object + impl_type = type(impl) + self.__cached_impl = impl + else: + raise TypeError( + '`impl` {!r} should be a `str`, or an object or type ' + 'having a `call_forward()` and/or `call_backward()`. ' + .format(type(impl))) self.__geometry = geometry self._impl_type = impl_type @@ -298,7 +267,7 @@ def impl(self): def create_impl(self, from_cache=True): """Fetches or instantiates implementation backend for evaluation.""" - # Skip impl creation (__cached_impl) when `clear_cache` is True + # Use impl creation (__cached_impl) when `from_cache` is True if not from_cache or self.__cached_impl is None: # Lazily (re)instantiate the backend self.__cached_impl = self._impl_type( From cfdf701b417cd6706f284520afef5fe9291aa513 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Sun, 16 Feb 2020 15:59:50 +0100 Subject: [PATCH 09/39] Remove RayTransformBase; RayBackProjection becomes Operator --- odl/tomo/operators/ray_trafo.py | 251 +++++++++++++------------------- 1 file changed, 105 insertions(+), 146 deletions(-) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 3aafe72ee43..8f08aa529ce 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -10,7 +10,7 @@ from __future__ import absolute_import, division, print_function -import warnings +import copy import numpy as np @@ -40,10 +40,10 @@ __all__ = ('RayTransform', 'RayBackProjection') -class RayTransformBase(Operator): +class RayTransform(Operator): """Base class for ray transforms containing common attributes.""" - def __init__(self, reco_space, geometry, variant, **kwargs): + def __init__(self, reco_space, geometry, **kwargs): """Initialize a new instance. Parameters @@ -54,9 +54,6 @@ def __init__(self, reco_space, geometry, variant, **kwargs): geometry : `Geometry` Geometry of the transform that contains information about the data structure. - variant : {'forward', 'backward'} - Variant of the transform, i.e., whether the ray transform - or its back-projection should be created. Other Parameters ---------------- @@ -90,21 +87,10 @@ def __init__(self, reco_space, geometry, variant, **kwargs): ``dtype='float32'`` and storage order 'C'. Otherwise copies will be needed. """ - variant, variant_in = str(variant).lower(), variant - if variant not in ('forward', 'backward'): - raise ValueError('`variant` {!r} not understood' - ''.format(variant_in)) - - if variant == 'forward': - reco_name = 'domain' - proj_name = 'range' - else: - reco_name = 'range' - proj_name = 'domain' - if not isinstance(reco_space, DiscretizedSpace): - raise TypeError('`{}` must be a `DiscretizedSpace` instance, got ' - '{!r}'.format(reco_name, reco_space)) + raise TypeError( + '`reco_space` must be a `DiscretizedSpace` instance, got ' + '{!r}'.format(reco_space)) if not isinstance(geometry, Geometry): raise TypeError('`geometry` must be a `Geometry` instance, got ' @@ -172,28 +158,55 @@ def __init__(self, reco_space, geometry, variant, **kwargs): else: # proj_space was given, checking some stuff if not isinstance(proj_space, DiscretizedSpace): - raise TypeError('`{}` must be a `DiscretizedSpace` instance, ' - 'got {!r}'.format(proj_name, proj_space)) + raise TypeError( + '`proj_space` must be a `DiscretizedSpace` instance, ' + 'got {!r}'.format(proj_space) + ) if proj_space.shape != geometry.partition.shape: - raise ValueError('`{}.shape` not equal to `geometry.shape`: ' - '{} != {}'.format(proj_name, proj_space.shape, - geometry.partition.shape)) + raise ValueError( + '`proj_space.shape` not equal to `geometry.shape`: ' + '{} != {}' + ''.format(proj_space.shape, geometry.partition.shape) + ) if proj_space.dtype != reco_space.dtype: - raise ValueError('`{}.dtype` not equal to `{}.dtype`: ' - '{} != {}'.format(proj_name, reco_name, - proj_space.dtype, - reco_space.dtype)) + raise ValueError( + '`proj_space.dtype` not equal to `reco_space.dtype`: ' + '{} != {}'.format(proj_space.dtype, reco_space.dtype) + ) if reco_space.ndim != geometry.ndim: - raise ValueError('`{}.ndim` not equal to `geometry.ndim`: ' - '{} != {}'.format(reco_name, reco_space.ndim, - geometry.ndim)) + raise ValueError( + '`reco_space.ndim` not equal to `geometry.ndim`: ' + '{} != {}'.format(reco_space.ndim, geometry.ndim) + ) # Cache for input/output arrays of transforms self.use_cache = kwargs.pop('use_cache', True) + # Check `impl` impl = kwargs.pop('impl', None) - self.__cached_impl = None + impl_type, self.__cached_impl = self._check_impl(impl) + self._impl_type = impl_type + self.__impl = impl.lower() if isinstance(impl, str) else impl_type + + self.geometry = geometry + self.__reco_space = reco_space + self.__proj_space = proj_space + + # Reserve name for cached properties (used for efficiency reasons) + self._adjoint = None + + # Extra kwargs that can be reused for adjoint etc. These must + # be retrieved with `get` instead of `pop` above. + self._extra_kwargs = kwargs + + # Finally, initialize the Operator structure + super(RayTransform, self).__init__( + domain=reco_space, range=proj_space, linear=True + ) + + def _check_impl(self, impl): + cached_impl = None if impl is None: # User didn't specify a backend if not _AVAILABLE_IMPLS: @@ -230,33 +243,15 @@ def __init__(self, reco_space, geometry, variant, **kwargs): # User gave an object for `impl`, meaning to set the # backend cache to an already initiated object impl_type = type(impl) - self.__cached_impl = impl + cached_impl = impl else: raise TypeError( '`impl` {!r} should be a `str`, or an object or type ' 'having a `call_forward()` and/or `call_backward()`. ' - .format(type(impl))) - - self.__geometry = geometry - self._impl_type = impl_type - self.__impl = impl.lower() if isinstance(impl, str) else impl_type - self.__reco_space = reco_space - self.__proj_space = proj_space - - # Reserve name for cached properties (used for efficiency reasons) - self._adjoint = None - - # Extra kwargs that can be reused for adjoint etc. These must - # be retrieved with `get` instead of `pop` above. - self._extra_kwargs = kwargs + ''.format(type(impl)) + ) - # Finally, initialize the Operator structure - if variant == 'forward': - super(RayTransformBase, self).__init__( - domain=reco_space, range=proj_space, linear=True) - elif variant == 'backward': - super(RayTransformBase, self).__init__( - domain=proj_space, range=reco_space, linear=True) + return impl_type, cached_impl @property def impl(self): @@ -264,11 +259,11 @@ def impl(self): return self.__impl - def create_impl(self, from_cache=True): + def create_impl(self, use_cache=True): """Fetches or instantiates implementation backend for evaluation.""" - # Use impl creation (__cached_impl) when `from_cache` is True - if not from_cache or self.__cached_impl is None: + # Use impl creation (__cached_impl) when `use_cache` is True + if not use_cache or self.__cached_impl is None: # Lazily (re)instantiate the backend self.__cached_impl = self._impl_type( self.geometry, @@ -277,11 +272,6 @@ def create_impl(self, from_cache=True): return self.__cached_impl - @property - def geometry(self): - """Geometry of this operator.""" - return self.__geometry - def _call(self, x, out=None): """Return ``self(x[, out])``.""" if self.domain.is_real: @@ -303,68 +293,9 @@ def _call(self, x, out=None): out.imag = result_parts[1] return out - else: raise RuntimeError('bad domain {!r}'.format(self.domain)) - -class RayTransform(RayTransformBase): - """Discrete Ray transform between L^p spaces.""" - - def __init__(self, domain, geometry, **kwargs): - """Initialize a new instance. - - Parameters - ---------- - domain : `DiscretizedSpace` - Discretized reconstruction space, the domain of the forward - projector. - geometry : `Geometry` - Geometry of the transform, containing information about - the operator range (projection/sinogram space). - - Other Parameters - ---------------- - impl : {`None`, 'astra_cuda', 'astra_cpu', 'skimage'}, optional - Implementation back-end for the transform. Supported back-ends: - - - ``'astra_cuda'``: ASTRA toolbox, using CUDA, 2D or 3D - - ``'astra_cpu'``: ASTRA toolbox using CPU, only 2D - - ``'skimage'``: scikit-image, only 2D parallel with square - reconstruction space. - - For the default ``None``, the fastest available back-end is - used, tried in the above order. - range : `DiscretizedSpace`, optional - Discretized projection (sinogram) space, the range of the - forward projector. - Default: Inferred from parameters. - use_cache : bool, optional - If ``True``, data is cached. This gives a significant speed-up - at the expense of a notable memory overhead, both on the GPU - and on the CPU, since a full volume and a projection dataset - are stored. That may be prohibitive in 3D. - Default: True - kwargs - Further keyword arguments passed to the projector backend. - - Notes - ----- - The ASTRA backend is faster if data are given with - ``dtype='float32'`` and storage order 'C'. Otherwise copies will be - needed. - - See Also - -------- - astra_cpu_forward_projector - AstraCudaProjectorImpl - skimage_radon_forward_projector - """ - range = kwargs.pop('range', None) - super(RayTransform, self).__init__( - reco_space=domain, proj_space=range, geometry=geometry, - variant='forward', **kwargs) - def _call_real(self, x_real, out_real, **kwargs): """Real-space forward projection for the current set-up.""" @@ -395,7 +326,7 @@ def adjoint(self): return self._adjoint -class RayBackProjection(RayTransformBase): +class RayBackProjection(Operator): """Adjoint of the discrete Ray transform between L^p spaces.""" def __init__(self, range, geometry, **kwargs): @@ -441,46 +372,74 @@ def __init__(self, range, geometry, **kwargs): The ASTRA backend is faster if data are given with ``dtype='float32'`` and storage order 'C'. Otherwise copies will be needed. - - See Also - -------- - astra_cpu_back_projector - AstraCudaBackProjectorImpl - skimage_radon_back_projector """ domain = kwargs.pop('domain', None) + super(RayBackProjection, self).__init__( - reco_space=range, proj_space=domain, geometry=geometry, - variant='backward', **kwargs) + domain=domain, range=range, linear=True + ) + + # This implementation simply uses a RayTransform to do the job + self._ray_trafo = RayTransform(reco_space=range, + proj_space=domain, + geometry=geometry, + **kwargs) + + # Extra kwargs that can be reused for adjoint etc. These must + # be retrieved with `get` instead of `pop` above. + self._extra_kwargs = kwargs + + @property + def geometry(self): + return self._ray_trafo.geometry + + @property + def use_cache(self): + return self._ray_trafo.use_cache + + def _call(self, x, out=None): + """Return ``self(x[, out])``.""" + if self.domain.is_real: + return self._call_real(x, out, **self._extra_kwargs) + + elif self.domain.is_complex: + result_parts = [ + self._call_real( + x.real, getattr(out, 'real', None), **self._extra_kwargs + ), + self._call_real( + x.imag, getattr(out, 'imag', None), **self._extra_kwargs + ), + ] + + if out is None: + out = self.range.element() + out.real = result_parts[0] + out.imag = result_parts[1] + + return out + + else: + raise RuntimeError('bad domain {!r}'.format(self.domain)) def _call_real(self, x_real, out_real, **kwargs): """Real-space back-projection for the current set-up.""" - return self.create_impl(self.use_cache) \ + return self._ray_trafo.create_impl(self._ray_trafo.use_cache) \ .call_backward(x_real, out_real, **kwargs) @property def adjoint(self): """Adjoint of this operator. + Returns a copy of the internally used `RayTransform`, thus + reusing the cached backend (if `use_cache`). + Returns ------- adjoint : `RayTransform` """ - if self._adjoint is None: - kwargs = self._extra_kwargs.copy() - kwargs['range'] = self.domain - - # initiate adjoint with cached implementation if `use_cache` - self._adjoint = RayTransform( - self.range, self.geometry, - impl=(self.impl - if not self.use_cache - else self.create_impl(self.use_cache)), - use_cache=self.use_cache, - **kwargs) - - return self._adjoint + return copy.copy(self._ray_trafo) if __name__ == '__main__': From d7d42286183f0f3ac6959e2d6aa46755cf7673f9 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Sun, 16 Feb 2020 17:37:58 +0100 Subject: [PATCH 10/39] Move _call logic into backend for potential optimization --- odl/tomo/backends/astra_cpu.py | 23 +++++++----- odl/tomo/backends/astra_cuda.py | 13 ++++++- odl/tomo/backends/skimage_radon.py | 19 +++++++--- odl/tomo/backends/util.py | 28 ++++++++++++++ odl/tomo/operators/ray_trafo.py | 59 +++--------------------------- 5 files changed, 71 insertions(+), 71 deletions(-) create mode 100644 odl/tomo/backends/util.py diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index 2c1fb6105cb..4059441244b 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -16,7 +16,7 @@ from odl.discr import DiscretizedSpace, DiscretizedSpaceElement from odl.tomo.backends.astra_setup import ( astra_algorithm, astra_data, astra_projection_geometry, astra_projector, - astra_volume_geometry) + astra_volume_geometry, complexify) from odl.tomo.geometry import ( DivergentBeamGeometry, Geometry, ParallelBeamGeometry) from odl.util import writable_array @@ -297,15 +297,20 @@ def __init__(self, geometry, reco_space, proj_space): self.reco_space = reco_space self.proj_space = proj_space - def call_forward(self, x_real, out_real, **kwargs): - return astra_cpu_forward_projector( - x_real, self.geometry, self.proj_space.real_space, - out_real, **kwargs) + def call_backward(self, x, out, **kwargs): + def wrapper(x, out, **kwargs): + return astra_cpu_forward_projector( + x, self.geometry, self.proj_space.real_space, out, **kwargs) - def call_backward(self, x_real, out_real, **kwargs): - return astra_cpu_back_projector( - x_real, self.geometry, self.reco_space.real_space, - out_real, **kwargs) + return complexify(wrapper, self.reco_space)(x, out, **kwargs) + + def call_forward(self, x, out, **kwargs): + def wrapper(x_real, out_real, **kwargs): + return astra_cpu_back_projector( + x_real, self.geometry, self.reco_space.real_space, + out_real, **kwargs) + + return complexify(wrapper, self.proj_space)(x, out, **kwargs) if __name__ == '__main__': diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 4cc5c8ef9d6..310b61e4895 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -17,6 +17,7 @@ from packaging.version import parse as parse_version from odl.discr import DiscretizedSpace +from odl.tomo.backends import complexify from odl.tomo.backends.astra_setup import ( ASTRA_VERSION, astra_algorithm, astra_data, astra_projection_geometry, astra_projector, astra_volume_geometry, astra_supports, @@ -158,7 +159,11 @@ def create_ids(self): 'backward', proj_ndim, self.vol_id, self.sino_id, proj_id=self.proj_id, impl='cuda') - def call_forward(self, vol_data, out, **kwargs): + def call_forward(self, x, out, **kwargs): + return complexify(self._call_forward_real, self.reco_space) \ + (x, out, **kwargs) + + def _call_forward_real(self, vol_data, out, **kwargs): """Run an ASTRA forward projection on the given data using the GPU. Parameters @@ -208,7 +213,11 @@ def call_forward(self, vol_data, out, **kwargs): return out - def call_backward(self, proj_data, out, **kwargs): + def call_backward(self, x, out, **kwargs): + return complexify(self._call_backward_real, self.proj_space) \ + (x, out, **kwargs) + + def _call_backward_real(self, proj_data, out, **kwargs): """Run an ASTRA back-projection on the given data using the GPU. Parameters diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index 94173311d50..c52dafefd72 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -18,6 +18,7 @@ from odl.discr import uniform_discr_frompartition, uniform_partition, \ DiscretizedSpace from odl.discr.discr_utils import linear_interpolator, point_collocation +from odl.tomo.backends import complexify from odl.util.utility import writable_array try: @@ -240,10 +241,16 @@ def __init__(self, geometry, reco_space, proj_space): self.reco_space = reco_space self.proj_space = proj_space - def call_forward(self, x_real, out_real, **kwargs): - return skimage_radon_forward_projector( - x_real, self.geometry, self.proj_space.real_space, out_real) + def call_forward(self, x, out, **kwargs): + def wrapper(x_real, out_real, **kwargs): + return skimage_radon_forward_projector( + x_real, self.geometry, self.proj_space.real_space, out_real) - def call_backward(self, x_real, out_real, **kwargs): - return skimage_radon_back_projector( - x_real, self.geometry, self.reco_space.real_space, out_real) + return complexify(wrapper, self.proj_space)(x, out, **kwargs) + + def call_backward(self, x, out, **kwargs): + def wrapper(x, out, **kwargs): + return skimage_radon_back_projector( + x, self.geometry, self.reco_space.real_space, out) + + return complexify(wrapper, self.reco_space)(x, out, **kwargs) diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py new file mode 100644 index 00000000000..e79c28088c9 --- /dev/null +++ b/odl/tomo/backends/util.py @@ -0,0 +1,28 @@ +__all__ = ('complexify',) + + +def complexify(fn, range): + """Wrapper to call a function twice when it is complex""" + + if not range.is_real and not range.is_complex: + raise RuntimeError('Range needs to be real or complex{!r}' + .format(range)) + + # No need to do anything when the range is real + if range.is_real: + return fn + + # Function takes a possibly complex `x` and delivers complex `out` + def complex_fn(x, out, **kwargs): + result_parts = [ + complex_fn(x.real, getattr(out, 'real', None), **kwargs), + complex_fn(x.imag, getattr(out, 'imag', None), **kwargs)] + + if out is None: + out = range.element() + out.real = result_parts[0] + out.imag = result_parts[1] + + return out + + return complex_fn diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 8f08aa529ce..91134b7b099 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -272,35 +272,11 @@ def create_impl(self, use_cache=True): return self.__cached_impl - def _call(self, x, out=None): - """Return ``self(x[, out])``.""" - if self.domain.is_real: - return self._call_real(x, out, **self._extra_kwargs) - - elif self.domain.is_complex: - result_parts = [ - self._call_real( - x.real, getattr(out, 'real', None), **self._extra_kwargs - ), - self._call_real( - x.imag, getattr(out, 'imag', None), **self._extra_kwargs - ), - ] - - if out is None: - out = self.range.element() - out.real = result_parts[0] - out.imag = result_parts[1] - - return out - else: - raise RuntimeError('bad domain {!r}'.format(self.domain)) - - def _call_real(self, x_real, out_real, **kwargs): + def _call(self, x, out, **kwargs): """Real-space forward projection for the current set-up.""" return self.create_impl(self.use_cache) \ - .call_forward(x_real, out_real, **kwargs) + .call_forward(x, out, **kwargs) @property def adjoint(self): @@ -397,36 +373,11 @@ def geometry(self): def use_cache(self): return self._ray_trafo.use_cache - def _call(self, x, out=None): - """Return ``self(x[, out])``.""" - if self.domain.is_real: - return self._call_real(x, out, **self._extra_kwargs) - - elif self.domain.is_complex: - result_parts = [ - self._call_real( - x.real, getattr(out, 'real', None), **self._extra_kwargs - ), - self._call_real( - x.imag, getattr(out, 'imag', None), **self._extra_kwargs - ), - ] - - if out is None: - out = self.range.element() - out.real = result_parts[0] - out.imag = result_parts[1] - - return out - - else: - raise RuntimeError('bad domain {!r}'.format(self.domain)) - - def _call_real(self, x_real, out_real, **kwargs): - """Real-space back-projection for the current set-up.""" + def _call(self, x, out, **kwargs): + """Back-projection for the current set-up.""" return self._ray_trafo.create_impl(self._ray_trafo.use_cache) \ - .call_backward(x_real, out_real, **kwargs) + .call_backward(x, out, **kwargs) @property def adjoint(self): From ad90713c560edd533dd32b1338b7630c9b0a67eb Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Sun, 16 Feb 2020 18:03:07 +0100 Subject: [PATCH 11/39] Remove backends from __all__, update imports --- odl/tomo/backends/__init__.py | 4 ++-- odl/tomo/backends/astra_cpu.py | 4 ++-- odl/tomo/backends/astra_cuda.py | 3 +-- odl/tomo/backends/skimage_radon.py | 3 +-- odl/tomo/backends/util.py | 2 +- odl/tomo/operators/ray_trafo.py | 6 ++++-- 6 files changed, 11 insertions(+), 11 deletions(-) diff --git a/odl/tomo/backends/__init__.py b/odl/tomo/backends/__init__.py index 6c30076eb75..6b3b9391f77 100644 --- a/odl/tomo/backends/__init__.py +++ b/odl/tomo/backends/__init__.py @@ -13,12 +13,12 @@ from .astra_cpu import * from .astra_cuda import * from .astra_setup import * -from .impl import * +from .util import * from .skimage_radon import * __all__ = () __all__ += astra_cpu.__all__ __all__ += astra_cuda.__all__ __all__ += astra_setup.__all__ -__all__ += impl.__all__ +__all__ += util.__all__ __all__ += skimage_radon.__all__ diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index 4059441244b..5c076e5eacf 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -16,7 +16,8 @@ from odl.discr import DiscretizedSpace, DiscretizedSpaceElement from odl.tomo.backends.astra_setup import ( astra_algorithm, astra_data, astra_projection_geometry, astra_projector, - astra_volume_geometry, complexify) + astra_volume_geometry) +from odl.tomo.backends.util import complexify from odl.tomo.geometry import ( DivergentBeamGeometry, Geometry, ParallelBeamGeometry) from odl.util import writable_array @@ -30,7 +31,6 @@ 'astra_cpu_forward_projector', 'astra_cpu_back_projector', 'default_astra_proj_type', - 'AstraCpu', ) diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 310b61e4895..e4fec785019 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -17,7 +17,7 @@ from packaging.version import parse as parse_version from odl.discr import DiscretizedSpace -from odl.tomo.backends import complexify +from odl.tomo.backends.util import complexify from odl.tomo.backends.astra_setup import ( ASTRA_VERSION, astra_algorithm, astra_data, astra_projection_geometry, astra_projector, astra_volume_geometry, astra_supports, @@ -35,7 +35,6 @@ __all__ = ( 'ASTRA_CUDA_AVAILABLE', - 'AstraCuda', ) diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index c52dafefd72..302f1096873 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -18,7 +18,7 @@ from odl.discr import uniform_discr_frompartition, uniform_partition, \ DiscretizedSpace from odl.discr.discr_utils import linear_interpolator, point_collocation -from odl.tomo.backends import complexify +from odl.tomo.backends.util import complexify from odl.util.utility import writable_array try: @@ -32,7 +32,6 @@ 'SKIMAGE_AVAILABLE', 'skimage_radon_forward_projector', 'skimage_radon_back_projector', - 'Skimage', ) diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index e79c28088c9..964aa1396c1 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -1,4 +1,4 @@ -__all__ = ('complexify',) +__all__ = tuple() def complexify(fn, range): diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 91134b7b099..0bf52d1eeb4 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -18,8 +18,10 @@ from odl.operator import Operator from odl.space.weighting import ConstWeighting from odl.tomo.backends import ( - ASTRA_AVAILABLE, ASTRA_CUDA_AVAILABLE, SKIMAGE_AVAILABLE, - Skimage, AstraCuda, AstraCpu) + ASTRA_AVAILABLE, ASTRA_CUDA_AVAILABLE, SKIMAGE_AVAILABLE) +from odl.tomo.backends.astra_cpu import AstraCpu +from odl.tomo.backends.astra_cuda import AstraCuda +from odl.tomo.backends.skimage_radon import Skimage from odl.tomo.geometry import Geometry # Backends that are implemented in ODL and can be specified via `impl` From 2e3c2748abf9cca250ae72302b1f6deb3259b581 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Wed, 19 Feb 2020 19:44:00 +0100 Subject: [PATCH 12/39] Bring `RayBackProjection` class inline in `RayTransform.adjoint` --- odl/tomo/operators/ray_trafo.py | 122 ++++++-------------------------- 1 file changed, 20 insertions(+), 102 deletions(-) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 0bf52d1eeb4..2095a17dbf8 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -10,8 +10,6 @@ from __future__ import absolute_import, division, print_function -import copy - import numpy as np from odl.discr import DiscretizedSpace @@ -39,11 +37,11 @@ if SKIMAGE_AVAILABLE: _AVAILABLE_IMPLS.append('skimage') -__all__ = ('RayTransform', 'RayBackProjection') +__all__ = ('RayTransform',) class RayTransform(Operator): - """Base class for ray transforms containing common attributes.""" + """Linear X-Ray (Radon) transform operator between L^p spaces.""" def __init__(self, reco_space, geometry, **kwargs): """Initialize a new instance. @@ -289,110 +287,30 @@ def adjoint(self): adjoint : `RayBackProjection` """ if self._adjoint is None: - kwargs = self._extra_kwargs.copy() - kwargs['domain'] = self.range - - # initiate adjoint with cached implementation if `use_cache` - self._adjoint = RayBackProjection( - self.domain, self.geometry, - impl=(self.impl - if not self.use_cache - else self.create_impl(self.use_cache)), - use_cache=self.use_cache, - **kwargs) - - return self._adjoint - - -class RayBackProjection(Operator): - """Adjoint of the discrete Ray transform between L^p spaces.""" + # bring `self` into scope to prevent shadowing in inline class + ray_trafo = self - def __init__(self, range, geometry, **kwargs): - """Initialize a new instance. - - Parameters - ---------- - range : `DiscretizedSpace` - Discretized reconstruction space, the range of the - backprojection operator. - geometry : `Geometry` - Geometry of the transform, containing information about - the operator domain (projection/sinogram space). - - Other Parameters - ---------------- - impl : {`None`, 'astra_cuda', 'astra_cpu', 'skimage'}, optional - Implementation back-end for the transform. Supported back-ends: + class RayBackProjection(Operator): + """Adjoint of the discrete Ray transform between L^p spaces.""" - - ``'astra_cuda'``: ASTRA toolbox, using CUDA, 2D or 3D - - ``'astra_cpu'``: ASTRA toolbox using CPU, only 2D - - ``'skimage'``: scikit-image, only 2D parallel with square - reconstruction space. + def _call(self, x, out, **kwargs): + """Back-projection for the current set-up.""" + return ray_trafo.create_impl(ray_trafo.use_cache) \ + .call_backward(x, out, **kwargs) - For the default ``None``, the fastest available back-end is - used, tried in the above order. + @property + def geometry(self): + return ray_trafo.geometry - domain : `DiscretizedSpace`, optional - Discretized projection (sinogram) space, the domain of the - backprojection operator. - Default: Inferred from parameters. - use_cache : bool, optional - If ``True``, data is cached. This gives a significant speed-up - at the expense of a notable memory overhead, both on the GPU - and on the CPU, since a full volume and a projection dataset - are stored. That may be prohibitive in 3D. - Default: True - kwargs - Further keyword arguments passed to the projector backend. + @property + def adjoint(self): + return ray_trafo - Notes - ----- - The ASTRA backend is faster if data are given with - ``dtype='float32'`` and storage order 'C'. Otherwise copies will be - needed. - """ - domain = kwargs.pop('domain', None) - - super(RayBackProjection, self).__init__( - domain=domain, range=range, linear=True - ) - - # This implementation simply uses a RayTransform to do the job - self._ray_trafo = RayTransform(reco_space=range, - proj_space=domain, - geometry=geometry, - **kwargs) - - # Extra kwargs that can be reused for adjoint etc. These must - # be retrieved with `get` instead of `pop` above. - self._extra_kwargs = kwargs - - @property - def geometry(self): - return self._ray_trafo.geometry - - @property - def use_cache(self): - return self._ray_trafo.use_cache - - def _call(self, x, out, **kwargs): - """Back-projection for the current set-up.""" - - return self._ray_trafo.create_impl(self._ray_trafo.use_cache) \ - .call_backward(x, out, **kwargs) - - @property - def adjoint(self): - """Adjoint of this operator. - - Returns a copy of the internally used `RayTransform`, thus - reusing the cached backend (if `use_cache`). + kwargs = self._extra_kwargs.copy() + kwargs['domain'] = self.range + self._adjoint = RayBackProjection(range=self.domain, **kwargs) - Returns - ------- - adjoint : `RayTransform` - """ - return copy.copy(self._ray_trafo) + return self._adjoint if __name__ == '__main__': From bc329b92297f0598c0cd7730223fa3075d1321c2 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Wed, 19 Feb 2020 22:01:22 +0100 Subject: [PATCH 13/39] Fix linear=True kwarg in `RayBackProjection` --- odl/tomo/operators/ray_trafo.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 2095a17dbf8..e7c22689268 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -308,7 +308,8 @@ def adjoint(self): kwargs = self._extra_kwargs.copy() kwargs['domain'] = self.range - self._adjoint = RayBackProjection(range=self.domain, **kwargs) + self._adjoint = RayBackProjection(range=self.domain, linear=True, + **kwargs) return self._adjoint From a1f54eb78fc76ee3b9d3780852cd1635475ffaba Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 20 Feb 2020 22:41:24 +0100 Subject: [PATCH 14/39] Decorate RayTransform backend calls with `_add_default_complex_impl` --- odl/tomo/backends/astra_cpu.py | 19 +++++------- odl/tomo/backends/astra_cuda.py | 14 ++++----- odl/tomo/backends/skimage_radon.py | 18 +++++------ odl/tomo/backends/util.py | 48 ++++++++++++++---------------- 4 files changed, 44 insertions(+), 55 deletions(-) diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index 5c076e5eacf..a4ce0cf8d23 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -17,7 +17,7 @@ from odl.tomo.backends.astra_setup import ( astra_algorithm, astra_data, astra_projection_geometry, astra_projector, astra_volume_geometry) -from odl.tomo.backends.util import complexify +from odl.tomo.backends.util import _add_default_complex_impl from odl.tomo.geometry import ( DivergentBeamGeometry, Geometry, ParallelBeamGeometry) from odl.util import writable_array @@ -297,20 +297,15 @@ def __init__(self, geometry, reco_space, proj_space): self.reco_space = reco_space self.proj_space = proj_space + @_add_default_complex_impl def call_backward(self, x, out, **kwargs): - def wrapper(x, out, **kwargs): - return astra_cpu_forward_projector( - x, self.geometry, self.proj_space.real_space, out, **kwargs) - - return complexify(wrapper, self.reco_space)(x, out, **kwargs) + return astra_cpu_back_projector( + x, self.geometry, self.reco_space.real_space, out, **kwargs) + @_add_default_complex_impl def call_forward(self, x, out, **kwargs): - def wrapper(x_real, out_real, **kwargs): - return astra_cpu_back_projector( - x_real, self.geometry, self.reco_space.real_space, - out_real, **kwargs) - - return complexify(wrapper, self.proj_space)(x, out, **kwargs) + return astra_cpu_forward_projector( + x, self.geometry, self.proj_space.real_space, out, **kwargs) if __name__ == '__main__': diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index e4fec785019..c55676f6e1d 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -17,7 +17,7 @@ from packaging.version import parse as parse_version from odl.discr import DiscretizedSpace -from odl.tomo.backends.util import complexify +from odl.tomo.backends.util import _add_default_complex_impl from odl.tomo.backends.astra_setup import ( ASTRA_VERSION, astra_algorithm, astra_data, astra_projection_geometry, astra_projector, astra_volume_geometry, astra_supports, @@ -64,11 +64,11 @@ def __init__(self, geometry, reco_space, proj_space): raise TypeError('`geometry` must be a `Geometry` instance, got ' '{!r}'.format(geometry)) - if not isinstance(reco_space, DiscreteLp): + if not isinstance(reco_space, DiscretizedSpace): raise TypeError('`reco_space` must be a `DiscreteLP` instance, got' ' {!r}'.format(reco_space)) - if not isinstance(proj_space, DiscreteLp): + if not isinstance(proj_space, DiscretizedSpace): raise TypeError('`proj_space` must be a `DiscreteLP` instance, got' ' {!r}'.format(proj_space)) @@ -158,9 +158,9 @@ def create_ids(self): 'backward', proj_ndim, self.vol_id, self.sino_id, proj_id=self.proj_id, impl='cuda') + @_add_default_complex_impl def call_forward(self, x, out, **kwargs): - return complexify(self._call_forward_real, self.reco_space) \ - (x, out, **kwargs) + return self._call_forward_real(x, out, **kwargs) def _call_forward_real(self, vol_data, out, **kwargs): """Run an ASTRA forward projection on the given data using the GPU. @@ -212,9 +212,9 @@ def _call_forward_real(self, vol_data, out, **kwargs): return out + @_add_default_complex_impl def call_backward(self, x, out, **kwargs): - return complexify(self._call_backward_real, self.proj_space) \ - (x, out, **kwargs) + return self._call_backward_real(x, out, **kwargs) def _call_backward_real(self, proj_data, out, **kwargs): """Run an ASTRA back-projection on the given data using the GPU. diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index 302f1096873..10b183f4f18 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -18,7 +18,7 @@ from odl.discr import uniform_discr_frompartition, uniform_partition, \ DiscretizedSpace from odl.discr.discr_utils import linear_interpolator, point_collocation -from odl.tomo.backends.util import complexify +from odl.tomo.backends.util import _add_default_complex_impl from odl.util.utility import writable_array try: @@ -240,16 +240,12 @@ def __init__(self, geometry, reco_space, proj_space): self.reco_space = reco_space self.proj_space = proj_space + @_add_default_complex_impl def call_forward(self, x, out, **kwargs): - def wrapper(x_real, out_real, **kwargs): - return skimage_radon_forward_projector( - x_real, self.geometry, self.proj_space.real_space, out_real) - - return complexify(wrapper, self.proj_space)(x, out, **kwargs) + return skimage_radon_forward_projector( + x, self.geometry, self.proj_space.real_space, out) + @_add_default_complex_impl def call_backward(self, x, out, **kwargs): - def wrapper(x, out, **kwargs): - return skimage_radon_back_projector( - x, self.geometry, self.reco_space.real_space, out) - - return complexify(wrapper, self.reco_space)(x, out, **kwargs) + return skimage_radon_back_projector( + x, self.geometry, self.reco_space.real_space, out) diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index 964aa1396c1..e9241de0b26 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -1,28 +1,26 @@ __all__ = tuple() -def complexify(fn, range): - """Wrapper to call a function twice when it is complex""" - - if not range.is_real and not range.is_complex: - raise RuntimeError('Range needs to be real or complex{!r}' - .format(range)) - - # No need to do anything when the range is real - if range.is_real: - return fn - - # Function takes a possibly complex `x` and delivers complex `out` - def complex_fn(x, out, **kwargs): - result_parts = [ - complex_fn(x.real, getattr(out, 'real', None), **kwargs), - complex_fn(x.imag, getattr(out, 'imag', None), **kwargs)] - - if out is None: - out = range.element() - out.real = result_parts[0] - out.imag = result_parts[1] - - return out - - return complex_fn +def _add_default_complex_impl(fn): + """Wrapper to call a class method twice when it is complex.""" + + def wrapper(self, x, out, **kwargs): + if self.reco_space.is_real and self.proj_space.is_real: + fn(self, x, out, **kwargs) + return out + elif self.reco_space.is_complex and self.proj_space.is_complex: + result_parts = [ + fn(x.real, getattr(out, 'real', None), **kwargs), + fn(x.imag, getattr(out, 'imag', None), **kwargs) + ] + + if out is None: + out = range.element() + out.real = result_parts[0] + out.imag = result_parts[1] + return out + else: + raise RuntimeError('Domain and range need to be both real ' + 'or both complex.') + + return wrapper From 3e623e18f611126bc1ee55225465220cc0d6140c Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 20 Feb 2020 23:07:18 +0100 Subject: [PATCH 15/39] Fix import --- odl/tomo/backends/astra_cpu.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index a4ce0cf8d23..40a3a5c107b 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -273,13 +273,17 @@ def __init__(self, geometry, reco_space, proj_space): raise TypeError('`geometry` must be a `Geometry` instance, got ' '{!r}'.format(geometry)) - if not isinstance(reco_space, DiscreteLp): - raise TypeError('`reco_space` must be a `DiscreteLP` instance, got' - ' {!r}'.format(reco_space)) - - if not isinstance(proj_space, DiscreteLp): - raise TypeError('`proj_space` must be a `DiscreteLP` instance, got' - ' {!r}'.format(proj_space)) + if not isinstance(reco_space, DiscretizedSpace): + raise TypeError( + '`reco_space` must be a `DiscretizedSpace` instance, got {!r}' + ''.format(reco_space) + ) + + if not isinstance(proj_space, DiscretizedSpace): + raise TypeError( + '`proj_space` must be a `DiscretizedSpace` instance, got {!r}' + ''.format(proj_space) + ) if reco_space.size >= 512 ** 2: warnings.warn( From 25a9b914a6f8014e28a3b0e45a8f0c2d2a31e7da Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 20 Feb 2020 23:09:09 +0100 Subject: [PATCH 16/39] Use `impl_type.__name__` as a return value for custom types --- odl/tomo/operators/ray_trafo.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index e7c22689268..a61ab1dbe54 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -187,10 +187,11 @@ def __init__(self, reco_space, geometry, **kwargs): impl = kwargs.pop('impl', None) impl_type, self.__cached_impl = self._check_impl(impl) self._impl_type = impl_type - self.__impl = impl.lower() if isinstance(impl, str) else impl_type + self.__impl = impl.lower() \ + if isinstance(impl, str) else impl_type.__name__ - self.geometry = geometry - self.__reco_space = reco_space + self._geometry = geometry + self.__vol_space = vol_space self.__proj_space = proj_space # Reserve name for cached properties (used for efficiency reasons) From edeafe12e12f908ee77209d7366b6a60c9588455 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 20 Feb 2020 23:09:47 +0100 Subject: [PATCH 17/39] Add `geometry` with @property --- odl/tomo/operators/ray_trafo.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index a61ab1dbe54..5f57dfcabd5 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -279,6 +279,10 @@ def _call(self, x, out, **kwargs): return self.create_impl(self.use_cache) \ .call_forward(x, out, **kwargs) + @property + def geometry(self): + return self._geometry + @property def adjoint(self): """Adjoint of this operator. From af7cea795b41c9f54643f9a74a5c41c625619624 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 20 Feb 2020 23:11:19 +0100 Subject: [PATCH 18/39] Make `_check_impl` static --- odl/tomo/operators/ray_trafo.py | 1 + 1 file changed, 1 insertion(+) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 5f57dfcabd5..1931bfe7849 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -206,6 +206,7 @@ def __init__(self, reco_space, geometry, **kwargs): domain=reco_space, range=proj_space, linear=True ) + @staticmethod def _check_impl(self, impl): cached_impl = None From b39fa5f7e1c08e43a0e96f9a99a173e3b1844279 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 20 Feb 2020 23:13:00 +0100 Subject: [PATCH 19/39] Change class names of implementations --- odl/tomo/backends/astra_cpu.py | 4 +++- odl/tomo/backends/astra_cuda.py | 2 +- odl/tomo/operators/ray_trafo.py | 20 ++++++++++---------- 3 files changed, 14 insertions(+), 12 deletions(-) diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index 40a3a5c107b..622beaa58db 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -255,7 +255,9 @@ def astra_cpu_back_projector(proj_data, geometry, vol_space, out=None, return out -class AstraCpu: +class AstraCpuImpl: + """Thin wrapper around the ASTRA CPU implementations.""" + def __init__(self, geometry, reco_space, proj_space): """Initialize a new instance. diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index c55676f6e1d..00c08212345 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -38,7 +38,7 @@ ) -class AstraCuda: +class AstraCudaImpl: """Thin wrapper around ASTRA.""" algo_forward_id = None diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 1931bfe7849..ba873ea2a48 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -17,16 +17,16 @@ from odl.space.weighting import ConstWeighting from odl.tomo.backends import ( ASTRA_AVAILABLE, ASTRA_CUDA_AVAILABLE, SKIMAGE_AVAILABLE) -from odl.tomo.backends.astra_cpu import AstraCpu -from odl.tomo.backends.astra_cuda import AstraCuda -from odl.tomo.backends.skimage_radon import Skimage +from odl.tomo.backends.astra_cpu import AstraCpuImpl +from odl.tomo.backends.astra_cuda import AstraCudaImpl +from odl.tomo.backends.skimage_radon import SkimageImpl from odl.tomo.geometry import Geometry # Backends that are implemented in ODL and can be specified via `impl` -_SUPPORTED_IMPL = { - 'astra_cpu': AstraCpu, - 'astra_cuda': AstraCuda, - 'skimage': Skimage} +_IMPL_STR2TYPE = { + 'astra_cpu': AstraCpuImpl, + 'astra_cuda': AstraCudaImpl, + 'skimage': SkimageImpl} # Backends that are installed, ordered by preference _AVAILABLE_IMPLS = [] @@ -217,19 +217,19 @@ def _check_impl(self, impl): 'check the install docs') # Select fastest available, _AVAILABLE_IMPLS is sorted by speed - impl_type = _SUPPORTED_IMPL[_AVAILABLE_IMPLS[0]] + impl_type = _IMPL_STR2TYPE[_AVAILABLE_IMPLS[0]] else: # User did specify `impl` if isinstance(impl, str): - if impl.lower() not in _SUPPORTED_IMPL: + if impl.lower() not in _IMPL_STR2TYPE: raise ValueError('`impl` {!r} not understood'.format(impl)) if impl.lower() not in _AVAILABLE_IMPLS: raise ValueError( '{!r} back-end not available'.format(impl)) - impl_type = _SUPPORTED_IMPL[impl.lower()] + impl_type = _IMPL_STR2TYPE[impl.lower()] elif isinstance(impl, type) or isinstance(impl, object): # User gave the type and leaves instantiation to us forward = getattr(impl, "call_forward", None) From e2e1319da3d908d6be677a41e7a289862f121bfd Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 20 Feb 2020 23:16:40 +0100 Subject: [PATCH 20/39] Change `reco_space` into `vol_space` and formatting --- odl/tomo/backends/astra_cpu.py | 36 +++++++------ odl/tomo/backends/astra_cuda.py | 83 ++++++++++++++++-------------- odl/tomo/backends/skimage_radon.py | 42 +++++++-------- odl/tomo/operators/ray_trafo.py | 51 +++++++++--------- 4 files changed, 109 insertions(+), 103 deletions(-) diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index 622beaa58db..33ebcc8ad1d 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -118,8 +118,10 @@ def astra_cpu_forward_projector(vol_data, geometry, proj_space, out=None, out = proj_space.element() else: if out not in proj_space: - raise TypeError('`out` {} is neither None nor a ' - 'DiscretizedSpaceElement instance'.format(out)) + raise TypeError( + '`out` {} is neither None nor a `DiscretizedSpaceElement` ' + 'instance'.format(out) + ) ndim = vol_data.ndim @@ -204,16 +206,18 @@ def astra_cpu_back_projector(proj_data, geometry, vol_space, out=None, "".format(vol_space.impl)) if vol_space.ndim != geometry.ndim: raise ValueError( - 'dimensions {} of reconstruction space and {} of ' - 'geometry do not match' + 'dimensions {} of reconstruction space and {} of geometry ' + 'do not match' ''.format(vol_space.ndim, geometry.ndim) ) if out is None: out = vol_space.element() else: if out not in vol_space: - raise TypeError('`out` {} is neither None nor a ' - 'DiscretizedSpaceElement instance'.format(out)) + raise TypeError( + '`out` {} is neither None nor a `DiscretizedSpaceElement` ' + 'instance'.format(out) + ) ndim = proj_data.ndim @@ -258,14 +262,14 @@ def astra_cpu_back_projector(proj_data, geometry, vol_space, out=None, class AstraCpuImpl: """Thin wrapper around the ASTRA CPU implementations.""" - def __init__(self, geometry, reco_space, proj_space): + def __init__(self, geometry, vol_space, proj_space): """Initialize a new instance. Parameters ---------- geometry : `Geometry` Geometry defining the tomographic setup. - reco_space : `DiscreteLp` + vol_space : `DiscreteLp` Reconstruction space, the space of the images to be forward projected. proj_space : `DiscreteLp` @@ -275,10 +279,10 @@ def __init__(self, geometry, reco_space, proj_space): raise TypeError('`geometry` must be a `Geometry` instance, got ' '{!r}'.format(geometry)) - if not isinstance(reco_space, DiscretizedSpace): + if not isinstance(vol_space, DiscretizedSpace): raise TypeError( - '`reco_space` must be a `DiscretizedSpace` instance, got {!r}' - ''.format(reco_space) + '`vol_space` must be a `DiscretizedSpace` instance, got {!r}' + ''.format(vol_space) ) if not isinstance(proj_space, DiscretizedSpace): @@ -287,26 +291,24 @@ def __init__(self, geometry, reco_space, proj_space): ''.format(proj_space) ) - if reco_space.size >= 512 ** 2: + if vol_space.size >= 512 ** 2: warnings.warn( "The 'astra_cpu' backend may be too slow for volumes of this " "size. Consider using 'astra_cuda' if your machine has an " - "Nvidia GPU. This warning can be disabled by explicitly " - "setting `impl='astra_cpu'`.", - RuntimeWarning) + "Nvidia GPU.", RuntimeWarning) if geometry.ndim > 2: raise ValueError('`impl` {!r} only works for 2d' ''.format(self.__name__)) self.geometry = geometry - self.reco_space = reco_space + self.vol_space = vol_space self.proj_space = proj_space @_add_default_complex_impl def call_backward(self, x, out, **kwargs): return astra_cpu_back_projector( - x, self.geometry, self.reco_space.real_space, out, **kwargs) + x, self.geometry, self.vol_space.real_space, out, **kwargs) @_add_default_complex_impl def call_forward(self, x, out, **kwargs): diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 00c08212345..b788ce18bb2 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -47,14 +47,14 @@ class AstraCudaImpl: sino_id = None proj_id = None - def __init__(self, geometry, reco_space, proj_space): + def __init__(self, geometry, vol_space, proj_space): """Initialize a new instance. Parameters ---------- geometry : `Geometry` Geometry defining the tomographic setup. - reco_space : `DiscretizedSpace` + vol_space : `DiscretizedSpace` Reconstruction space, the space of the images to be forward projected. proj_space : `DiscretizedSpace` @@ -64,13 +64,17 @@ def __init__(self, geometry, reco_space, proj_space): raise TypeError('`geometry` must be a `Geometry` instance, got ' '{!r}'.format(geometry)) - if not isinstance(reco_space, DiscretizedSpace): - raise TypeError('`reco_space` must be a `DiscreteLP` instance, got' - ' {!r}'.format(reco_space)) + if not isinstance(vol_space, DiscretizedSpace): + raise TypeError( + '`vol_space` must be a `DiscreteLP` instance, got {!r}' + ''.format(vol_space) + ) if not isinstance(proj_space, DiscretizedSpace): - raise TypeError('`proj_space` must be a `DiscreteLP` instance, got' - ' {!r}'.format(proj_space)) + raise TypeError( + '`proj_space` must be a `DiscreteLP` instance, got {!r}' + ''.format(proj_space) + ) # Print a warning if the detector midpoint normal vector at any # angle is perpendicular to the geometry axis in parallel 3d @@ -82,8 +86,9 @@ def __init__(self, geometry, reco_space, proj_space): axis = geometry.axis mid_pt = geometry.det_params.mid_pt for i, angle in enumerate(geometry.angles): - if abs(np.dot(axis, - geometry.det_to_src(angle, mid_pt))) < 1e-4: + if abs( + np.dot(axis, geometry.det_to_src(angle, mid_pt)) + ) < 1e-4: warnings.warn( 'angle {}: detector midpoint normal {} is ' 'perpendicular to the geometry axis {} in ' @@ -95,7 +100,7 @@ def __init__(self, geometry, reco_space, proj_space): break self.geometry = geometry - self.reco_space = reco_space + self.vol_space = vol_space self.proj_space = proj_space self.create_ids() @@ -117,23 +122,23 @@ def create_ids(self): if proj_ndim == 2: astra_proj_shape = proj_shape - astra_vol_shape = self.reco_space.shape + astra_vol_shape = self.vol_space.shape elif proj_ndim == 3: # The `u` and `v` axes of the projection data are swapped, # see explanation in `astra_*_3d_geom_to_vec`. astra_proj_shape = (proj_shape[1], proj_shape[0], proj_shape[2]) - astra_vol_shape = self.reco_space.shape + astra_vol_shape = self.vol_space.shape self.vol_array = np.empty(astra_vol_shape, dtype='float32', order='C') self.proj_array = np.empty(astra_proj_shape, dtype='float32', order='C') # Create ASTRA data structures - vol_geom = astra_volume_geometry(self.reco_space) + vol_geom = astra_volume_geometry(self.vol_space) proj_geom = astra_projection_geometry(self.geometry) self.vol_id = astra_data(vol_geom, datatype='volume', - ndim=self.reco_space.ndim, + ndim=self.vol_space.ndim, data=self.vol_array, allow_copy=False) @@ -167,7 +172,7 @@ def _call_forward_real(self, vol_data, out, **kwargs): Parameters ---------- - vol_data : ``reco_space`` element + vol_data : ``vol_space`` element Volume data to which the projector is applied. out : ``proj_space`` element, optional Element of the projection space to which the result is written. If @@ -180,7 +185,7 @@ def _call_forward_real(self, vol_data, out, **kwargs): If ``out`` was provided, the returned object is a reference to it. """ with self._mutex: - assert vol_data in self.reco_space + assert vol_data in self.vol_space if out is not None: assert out in self.proj_space else: @@ -223,13 +228,13 @@ def _call_backward_real(self, proj_data, out, **kwargs): ---------- proj_data : ``proj_space`` element Projection data to which the back-projector is applied. - out : ``reco_space`` element, optional + out : ``vol_space`` element, optional Element of the reconstruction space to which the result is written. - If ``None``, an element in ``reco_space`` is created. + If ``None``, an element in ``vol_space`` is created. Returns ------- - out : ``reco_space`` element + out : ``vol_space`` element Reconstruction data resulting from the application of the back-projector. If ``out`` was provided, the returned object is a reference to it. @@ -238,9 +243,9 @@ def _call_backward_real(self, proj_data, out, **kwargs): assert proj_data in self.proj_space if out is not None: - assert out in self.reco_space + assert out in self.vol_space else: - out = self.reco_space.element() + out = self.vol_space.element() # Copy data to GPU memory if self.geometry.ndim == 2: @@ -260,7 +265,7 @@ def _call_backward_real(self, proj_data, out, **kwargs): # Fix scaling to weight by pixel/voxel size out *= astra_cuda_bp_scaling_factor( - self.proj_space, self.reco_space, self.geometry) + self.proj_space, self.vol_space, self.geometry) return out @@ -288,7 +293,7 @@ def __del__(self): self.proj_id = None -def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry): +def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): """Volume scaling accounting for differing adjoint definitions. ASTRA defines the adjoint operator in terms of a fully discrete @@ -314,18 +319,18 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry): scaling_factor *= (proj_space.weighting.const / proj_weighting) - scaling_factor /= (reco_space.weighting.const / - reco_space.cell_volume) + scaling_factor /= (vol_space.weighting.const / + vol_space.cell_volume) if parse_version(ASTRA_VERSION) < parse_version('1.8rc1'): # Scaling for the old, pre-1.8 behaviour if isinstance(geometry, Parallel2dGeometry): # Scales with 1 / cell_volume - scaling_factor *= float(reco_space.cell_volume) + scaling_factor *= float(vol_space.cell_volume) elif (isinstance(geometry, FanBeamGeometry) and geometry.det_curvature_radius is None): # Scales with 1 / cell_volume - scaling_factor *= float(reco_space.cell_volume) + scaling_factor *= float(vol_space.cell_volume) # Additional magnification correction src_radius = geometry.src_radius det_radius = geometry.det_radius @@ -333,13 +338,13 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry): elif isinstance(geometry, Parallel3dAxisGeometry): # Scales with voxel stride # In 1.7, only cubic voxels are supported - voxel_stride = reco_space.cell_sides[0] + voxel_stride = vol_space.cell_sides[0] scaling_factor /= float(voxel_stride) elif (isinstance(geometry, ConeBeamGeometry) and geometry.det_curvature_radius is None): # Scales with 1 / cell_volume # In 1.7, only cubic voxels are supported - voxel_stride = reco_space.cell_sides[0] + voxel_stride = vol_space.cell_sides[0] scaling_factor /= float(voxel_stride) # Magnification correction src_radius = geometry.src_radius @@ -349,11 +354,11 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry): # Scaling for the 1.8.x releases if isinstance(geometry, Parallel2dGeometry): # Scales with 1 / cell_volume - scaling_factor *= float(reco_space.cell_volume) + scaling_factor *= float(vol_space.cell_volume) elif (isinstance(geometry, FanBeamGeometry) and geometry.det_curvature_radius is None): # Scales with 1 / cell_volume - scaling_factor *= float(reco_space.cell_volume) + scaling_factor *= float(vol_space.cell_volume) # Magnification correction src_radius = geometry.src_radius det_radius = geometry.det_radius @@ -361,11 +366,11 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry): elif isinstance(geometry, Parallel3dAxisGeometry): # Scales with cell volume # currently only square voxels are supported - scaling_factor /= reco_space.cell_volume + scaling_factor /= vol_space.cell_volume elif (isinstance(geometry, ConeBeamGeometry) and geometry.det_curvature_radius is None): # Scales with cell volume - scaling_factor /= reco_space.cell_volume + scaling_factor /= vol_space.cell_volume # Magnification correction (scaling = 1 / magnification ** 2) src_radius = geometry.src_radius det_radius = geometry.det_radius @@ -377,16 +382,16 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry): # density weighting. det_px_area = geometry.det_partition.cell_volume scaling_factor *= (src_radius ** 2 * det_px_area ** 2 / - reco_space.cell_volume ** 2) + vol_space.cell_volume ** 2) elif parse_version(ASTRA_VERSION) < parse_version('1.9.9.dev'): # Scaling for intermediate dev releases between 1.8.3 and 1.9.9.dev if isinstance(geometry, Parallel2dGeometry): # Scales with 1 / cell_volume - scaling_factor *= float(reco_space.cell_volume) + scaling_factor *= float(vol_space.cell_volume) elif (isinstance(geometry, FanBeamGeometry) and geometry.det_curvature_radius is None): # Scales with 1 / cell_volume - scaling_factor *= float(reco_space.cell_volume) + scaling_factor *= float(vol_space.cell_volume) # Magnification correction src_radius = geometry.src_radius det_radius = geometry.det_radius @@ -394,11 +399,11 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry): elif isinstance(geometry, Parallel3dAxisGeometry): # Scales with cell volume # currently only square voxels are supported - scaling_factor /= reco_space.cell_volume + scaling_factor /= vol_space.cell_volume elif (isinstance(geometry, ConeBeamGeometry) and geometry.det_curvature_radius is None): # Scales with cell volume - scaling_factor /= reco_space.cell_volume + scaling_factor /= vol_space.cell_volume # Magnification correction (scaling = 1 / magnification ** 2) src_radius = geometry.src_radius det_radius = geometry.det_radius @@ -412,7 +417,7 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry): scaling_factor *= (src_radius ** 2 * det_px_area ** 2) else: # Scaling for versions since 1.9.9.dev - scaling_factor /= float(reco_space.cell_volume) + scaling_factor /= float(vol_space.cell_volume) scaling_factor *= float(geometry.det_partition.cell_volume) return scaling_factor diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index 10b183f4f18..a3cd3673ea4 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -15,8 +15,8 @@ from odl.tomo import Parallel2dGeometry, Geometry -from odl.discr import uniform_discr_frompartition, uniform_partition, \ - DiscretizedSpace +from odl.discr import ( + uniform_discr_frompartition, uniform_partition, DiscretizedSpace) from odl.discr.discr_utils import linear_interpolator, point_collocation from odl.tomo.backends.util import _add_default_complex_impl from odl.util.utility import writable_array @@ -179,15 +179,15 @@ def skimage_radon_back_projector(sinogram, geometry, vol_space, out=None): return out -class Skimage: - def __init__(self, geometry, reco_space, proj_space): +class SkimageImpl: + def __init__(self, geometry, vol_space, proj_space): """Initialize a new instance. Parameters ---------- geometry : `Geometry` Geometry defining the tomographic setup. - reco_space : `DiscreteLp` + vol_space : `DiscreteLp` Reconstruction space, the space of the images to be forward projected. proj_space : `DiscreteLp` @@ -197,10 +197,10 @@ def __init__(self, geometry, reco_space, proj_space): raise TypeError('`geometry` must be a `Geometry` instance, got ' '{!r}'.format(geometry)) - if not isinstance(reco_space, DiscretizedSpace): + if not isinstance(vol_space, DiscretizedSpace): raise TypeError( - '`reco_space` must be a `DiscretizedSpace` instance, got {!r}' - ''.format(reco_space) + '`vol_space` must be a `DiscretizedSpace` instance, got {!r}' + ''.format(vol_space) ) if not isinstance(proj_space, DiscretizedSpace): @@ -209,35 +209,33 @@ def __init__(self, geometry, reco_space, proj_space): ''.format(proj_space) ) - if reco_space.size >= 256 ** 2: + if vol_space.size >= 256 ** 2: warnings.warn( "The 'skimage' backend may be too slow for volumes of this " "size. Consider using 'astra_cpu', or 'astra_cuda' if your " - "machine has an Nvidia GPU. This warning can be disabled by " - "explicitly setting `impl='astra_cpu'`.", - RuntimeWarning) + "machine has an Nvidia GPU.", RuntimeWarning) if not isinstance(geometry, Parallel2dGeometry): raise TypeError("{!r} backend only supports 2d parallel " - 'geometries'.format(self.__name__)) + 'geometries'.format(self.__name__)) - mid_pt = reco_space.domain.mid_pt + mid_pt = vol_space.domain.mid_pt if not np.allclose(mid_pt, [0, 0]): raise ValueError('Reconstruction space must be centered at (0, 0),' 'got midpoint {}'.format(mid_pt)) - shape = reco_space.shape + shape = vol_space.shape if shape[0] != shape[1]: - raise ValueError('`reco_space.shape` must have equal entries, ' - 'got {}'.format(shape)) + raise ValueError('`vol_space.shape` must have equal entries, ' + 'got {}'.format(shape)) - extent = reco_space.domain.extent + extent = vol_space.domain.extent if extent[0] != extent[1]: - raise ValueError('`reco_space.extent` must have equal entries, ' - 'got {}'.format(extent)) + raise ValueError('`vol_space.extent` must have equal entries, ' + 'got {}'.format(extent)) self.geometry = geometry - self.reco_space = reco_space + self.vol_space = vol_space self.proj_space = proj_space @_add_default_complex_impl @@ -248,4 +246,4 @@ def call_forward(self, x, out, **kwargs): @_add_default_complex_impl def call_backward(self, x, out, **kwargs): return skimage_radon_back_projector( - x, self.geometry, self.reco_space.real_space, out) + x, self.geometry, self.vol_space.real_space, out) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index ba873ea2a48..b4718aec6c8 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -43,12 +43,12 @@ class RayTransform(Operator): """Linear X-Ray (Radon) transform operator between L^p spaces.""" - def __init__(self, reco_space, geometry, **kwargs): + def __init__(self, vol_space, geometry, **kwargs): """Initialize a new instance. Parameters ---------- - reco_space : `DiscretizedSpace` + vol_space : `DiscretizedSpace` Discretized reconstruction space, the domain of the forward operator or the range of the adjoint (back-projection). geometry : `Geometry` @@ -87,10 +87,10 @@ def __init__(self, reco_space, geometry, **kwargs): ``dtype='float32'`` and storage order 'C'. Otherwise copies will be needed. """ - if not isinstance(reco_space, DiscretizedSpace): + if not isinstance(vol_space, DiscretizedSpace): raise TypeError( - '`reco_space` must be a `DiscretizedSpace` instance, got ' - '{!r}'.format(reco_space)) + '`vol_space` must be a `DiscretizedSpace` instance, got ' + '{!r}'.format(vol_space)) if not isinstance(geometry, Geometry): raise TypeError('`geometry` must be a `Geometry` instance, got ' @@ -99,13 +99,14 @@ def __init__(self, reco_space, geometry, **kwargs): # Generate or check projection space proj_space = kwargs.pop('proj_space', None) if proj_space is None: - dtype = reco_space.dtype + dtype = vol_space.dtype - if not reco_space.is_weighted: + if not vol_space.is_weighted: weighting = None - elif (isinstance(reco_space.weighting, ConstWeighting) - and np.isclose(reco_space.weighting.const, - reco_space.cell_volume)): + elif ( + isinstance(vol_space.weighting, ConstWeighting) + and np.isclose(vol_space.weighting.const, vol_space.cell_volume) + ): # Approximate cell volume # TODO: find a way to treat angles and detector differently # regarding weighting. While the detector should be uniformly @@ -119,9 +120,9 @@ def __init__(self, reco_space, geometry, **kwargs): else: raise NotImplementedError('unknown weighting of domain') - proj_tspace = reco_space.tspace_type(geometry.partition.shape, - weighting=weighting, - dtype=dtype) + proj_tspace = vol_space.tspace_type(geometry.partition.shape, + weighting=weighting, + dtype=dtype) if geometry.motion_partition.ndim == 0: angle_labels = [] @@ -168,16 +169,16 @@ def __init__(self, reco_space, geometry, **kwargs): '{} != {}' ''.format(proj_space.shape, geometry.partition.shape) ) - if proj_space.dtype != reco_space.dtype: + if proj_space.dtype != vol_space.dtype: raise ValueError( - '`proj_space.dtype` not equal to `reco_space.dtype`: ' - '{} != {}'.format(proj_space.dtype, reco_space.dtype) + '`proj_space.dtype` not equal to `vol_space.dtype`: ' + '{} != {}'.format(proj_space.dtype, vol_space.dtype) ) - if reco_space.ndim != geometry.ndim: + if vol_space.ndim != geometry.ndim: raise ValueError( - '`reco_space.ndim` not equal to `geometry.ndim`: ' - '{} != {}'.format(reco_space.ndim, geometry.ndim) + '`vol_space.ndim` not equal to `geometry.ndim`: ' + '{} != {}'.format(vol_space.ndim, geometry.ndim) ) # Cache for input/output arrays of transforms @@ -203,12 +204,12 @@ def __init__(self, reco_space, geometry, **kwargs): # Finally, initialize the Operator structure super(RayTransform, self).__init__( - domain=reco_space, range=proj_space, linear=True + domain=vol_space, range=proj_space, linear=True ) @staticmethod - def _check_impl(self, impl): - cached_impl = None + def _check_impl(impl): + impl_instance = None if impl is None: # User didn't specify a backend if not _AVAILABLE_IMPLS: @@ -245,7 +246,7 @@ def _check_impl(self, impl): # User gave an object for `impl`, meaning to set the # backend cache to an already initiated object impl_type = type(impl) - cached_impl = impl + impl_instance = impl else: raise TypeError( '`impl` {!r} should be a `str`, or an object or type ' @@ -253,7 +254,7 @@ def _check_impl(self, impl): ''.format(type(impl)) ) - return impl_type, cached_impl + return impl_type, impl_instance @property def impl(self): @@ -269,7 +270,7 @@ def create_impl(self, use_cache=True): # Lazily (re)instantiate the backend self.__cached_impl = self._impl_type( self.geometry, - reco_space=self.__reco_space.real_space, + vol_space=self.__vol_space.real_space, proj_space=self.__proj_space.real_space) return self.__cached_impl From 8cf739872d945bd28c36f0e7c261a53025dd9624 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 20 Feb 2020 23:54:42 +0100 Subject: [PATCH 21/39] Change `reco_space` to `vol_space` --- odl/tomo/backends/util.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index e9241de0b26..3f5dbf86048 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -5,10 +5,10 @@ def _add_default_complex_impl(fn): """Wrapper to call a class method twice when it is complex.""" def wrapper(self, x, out, **kwargs): - if self.reco_space.is_real and self.proj_space.is_real: + if self.vol_space.is_real and self.proj_space.is_real: fn(self, x, out, **kwargs) return out - elif self.reco_space.is_complex and self.proj_space.is_complex: + elif self.vol_space.is_complex and self.proj_space.is_complex: result_parts = [ fn(x.real, getattr(out, 'real', None), **kwargs), fn(x.imag, getattr(out, 'imag', None), **kwargs) From c64993fda56145cdc85ed952ff515b4c8a8313c9 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Fri, 21 Feb 2020 00:00:01 +0100 Subject: [PATCH 22/39] Fix `self` in function call --- odl/tomo/backends/util.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index 3f5dbf86048..c2f654b2df1 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -10,8 +10,8 @@ def wrapper(self, x, out, **kwargs): return out elif self.vol_space.is_complex and self.proj_space.is_complex: result_parts = [ - fn(x.real, getattr(out, 'real', None), **kwargs), - fn(x.imag, getattr(out, 'imag', None), **kwargs) + fn(self, x.real, getattr(out, 'real', None), **kwargs), + fn(self, x.imag, getattr(out, 'imag', None), **kwargs) ] if out is None: From 7642de10dad1b4b5de06d0bbe172ed3790a5f9c4 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Sun, 23 Feb 2020 15:07:00 +0100 Subject: [PATCH 23/39] Fix complex spaces, and fix `out` argument --- odl/test/tomo/backends/astra_cuda_test.py | 20 ++++++++++---------- odl/tomo/backends/astra_cuda.py | 23 +++++++++++++---------- odl/tomo/backends/util.py | 7 +++---- odl/tomo/operators/ray_trafo.py | 10 ++++------ 4 files changed, 30 insertions(+), 30 deletions(-) diff --git a/odl/test/tomo/backends/astra_cuda_test.py b/odl/test/tomo/backends/astra_cuda_test.py index ad684d62467..83c1fd69ac8 100644 --- a/odl/test/tomo/backends/astra_cuda_test.py +++ b/odl/test/tomo/backends/astra_cuda_test.py @@ -14,8 +14,7 @@ import pytest import odl -from odl.tomo.backends.astra_cuda import ( - AstraCudaBackProjectorImpl, AstraCudaProjectorImpl) +from odl.tomo.backends.astra_cuda import AstraCudaImpl from odl.tomo.util.testutils import skip_if_no_astra_cuda @@ -85,24 +84,25 @@ def test_astra_cuda_projector(space_and_geometry): """Test ASTRA CUDA projector.""" # Create reco space and a phantom - reco_space, geom = space_and_geometry - phantom = odl.phantom.cuboid(reco_space) + vol_space, geom = space_and_geometry + phantom = odl.phantom.cuboid(vol_space) # Make projection space proj_space = odl.uniform_discr_frompartition(geom.partition, - dtype=reco_space.dtype) + dtype=vol_space.dtype) + + # create RayTransform implementation + astra_cuda = AstraCudaImpl(geom, vol_space, proj_space) # Forward evaluation - projector = AstraCudaProjectorImpl(geom, reco_space, proj_space) - proj_data = projector.call_forward(phantom) + proj_data = astra_cuda.call_forward(phantom) assert proj_data in proj_space assert proj_data.norm() > 0 assert np.all(proj_data.asarray() >= 0) # Backward evaluation - back_projector = AstraCudaBackProjectorImpl(geom, reco_space, proj_space) - backproj = back_projector.call_backward(proj_data) - assert backproj in reco_space + backproj = astra_cuda.call_backward(proj_data) + assert backproj in vol_space assert backproj.norm() > 0 assert np.all(proj_data.asarray() >= 0) diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index b788ce18bb2..c34f6e33a9e 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -164,16 +164,17 @@ def create_ids(self): proj_id=self.proj_id, impl='cuda') @_add_default_complex_impl - def call_forward(self, x, out, **kwargs): + def call_forward(self, x, out=None, **kwargs): return self._call_forward_real(x, out, **kwargs) - def _call_forward_real(self, vol_data, out, **kwargs): + def _call_forward_real(self, vol_data, out=None, **kwargs): """Run an ASTRA forward projection on the given data using the GPU. Parameters ---------- - vol_data : ``vol_space`` element - Volume data to which the projector is applied. + vol_data : ``vol_space.real_space`` element + Volume data to which the projector is applied. Although + ``vol_space`` may be complex, this element needs to be real. out : ``proj_space`` element, optional Element of the projection space to which the result is written. If ``None``, an element in `proj_space` is created. @@ -185,7 +186,8 @@ def _call_forward_real(self, vol_data, out, **kwargs): If ``out`` was provided, the returned object is a reference to it. """ with self._mutex: - assert vol_data in self.vol_space + assert vol_data in self.vol_space.real_space + if out is not None: assert out in self.proj_space else: @@ -218,16 +220,17 @@ def _call_forward_real(self, vol_data, out, **kwargs): return out @_add_default_complex_impl - def call_backward(self, x, out, **kwargs): + def call_backward(self, x, out=None, **kwargs): return self._call_backward_real(x, out, **kwargs) - def _call_backward_real(self, proj_data, out, **kwargs): + def _call_backward_real(self, proj_data, out=None, **kwargs): """Run an ASTRA back-projection on the given data using the GPU. Parameters ---------- - proj_data : ``proj_space`` element - Projection data to which the back-projector is applied. + proj_data : ``proj_space.real_space`` element + Projection data to which the back-projector is applied. Although + ``proj_space`` may be complex, this element needs to be real. out : ``vol_space`` element, optional Element of the reconstruction space to which the result is written. If ``None``, an element in ``vol_space`` is created. @@ -240,7 +243,7 @@ def _call_backward_real(self, proj_data, out, **kwargs): reference to it. """ with self._mutex: - assert proj_data in self.proj_space + assert proj_data in self.proj_space.real_space if out is not None: assert out in self.vol_space diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index c2f654b2df1..467e2508f5d 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -4,10 +4,9 @@ def _add_default_complex_impl(fn): """Wrapper to call a class method twice when it is complex.""" - def wrapper(self, x, out, **kwargs): + def wrapper(self, x, out=None, **kwargs): if self.vol_space.is_real and self.proj_space.is_real: - fn(self, x, out, **kwargs) - return out + return fn(self, x, out, **kwargs) elif self.vol_space.is_complex and self.proj_space.is_complex: result_parts = [ fn(self, x.real, getattr(out, 'real', None), **kwargs), @@ -15,7 +14,7 @@ def wrapper(self, x, out, **kwargs): ] if out is None: - out = range.element() + out = self.proj_space.element() out.real = result_parts[0] out.imag = result_parts[1] return out diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index b4718aec6c8..2040dcdaf0e 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -192,8 +192,6 @@ def __init__(self, vol_space, geometry, **kwargs): if isinstance(impl, str) else impl_type.__name__ self._geometry = geometry - self.__vol_space = vol_space - self.__proj_space = proj_space # Reserve name for cached properties (used for efficiency reasons) self._adjoint = None @@ -270,12 +268,12 @@ def create_impl(self, use_cache=True): # Lazily (re)instantiate the backend self.__cached_impl = self._impl_type( self.geometry, - vol_space=self.__vol_space.real_space, - proj_space=self.__proj_space.real_space) + vol_space=self.domain, + proj_space=self.range) return self.__cached_impl - def _call(self, x, out, **kwargs): + def _call(self, x, out=None, **kwargs): """Real-space forward projection for the current set-up.""" return self.create_impl(self.use_cache) \ @@ -300,7 +298,7 @@ def adjoint(self): class RayBackProjection(Operator): """Adjoint of the discrete Ray transform between L^p spaces.""" - def _call(self, x, out, **kwargs): + def _call(self, x, out=None, **kwargs): """Back-projection for the current set-up.""" return ray_trafo.create_impl(ray_trafo.use_cache) \ .call_backward(x, out, **kwargs) From 05d9fb2dbc407a4a15f158bf6bbf347ee48d41a1 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Sun, 23 Feb 2020 19:43:15 +0100 Subject: [PATCH 24/39] Add properties for `vol_space` and `proj_space` to implementation classes --- odl/tomo/backends/astra_cpu.py | 12 ++++++++++-- odl/tomo/backends/astra_cuda.py | 12 ++++++++++-- odl/tomo/backends/skimage_radon.py | 12 ++++++++++-- 3 files changed, 30 insertions(+), 6 deletions(-) diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index 33ebcc8ad1d..4102f02b599 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -302,8 +302,16 @@ def __init__(self, geometry, vol_space, proj_space): ''.format(self.__name__)) self.geometry = geometry - self.vol_space = vol_space - self.proj_space = proj_space + self._vol_space = vol_space + self._proj_space = proj_space + + @property + def vol_space(self): + return self._vol_space + + @property + def proj_space(self): + return self._proj_space @_add_default_complex_impl def call_backward(self, x, out, **kwargs): diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index c34f6e33a9e..4a9c10a6552 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -100,14 +100,22 @@ def __init__(self, geometry, vol_space, proj_space): break self.geometry = geometry - self.vol_space = vol_space - self.proj_space = proj_space + self._vol_space = vol_space + self._proj_space = proj_space self.create_ids() # ASTRA projectors are not thread-safe, thus we need to lock ourselves self._mutex = Lock() + @property + def vol_space(self): + return self._vol_space + + @property + def proj_space(self): + return self._proj_space + def create_ids(self): """Create ASTRA objects.""" # Create input and output arrays diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index a3cd3673ea4..793361710fd 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -235,8 +235,16 @@ def __init__(self, geometry, vol_space, proj_space): 'got {}'.format(extent)) self.geometry = geometry - self.vol_space = vol_space - self.proj_space = proj_space + self._vol_space = vol_space + self._proj_space = proj_space + + @property + def vol_space(self): + return self._vol_space + + @property + def proj_space(self): + return self._proj_space @_add_default_complex_impl def call_forward(self, x, out, **kwargs): From da3fd97eba0551ee6ea47dca2b58f80c45b840a6 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Sun, 23 Feb 2020 19:43:44 +0100 Subject: [PATCH 25/39] Add docstrings --- odl/tomo/backends/astra_cpu.py | 2 +- odl/tomo/backends/astra_cuda.py | 2 +- odl/tomo/backends/skimage_radon.py | 2 + odl/tomo/backends/util.py | 18 ++++++++- odl/tomo/operators/ray_trafo.py | 64 ++++++++++++++++++++++++++++-- 5 files changed, 81 insertions(+), 7 deletions(-) diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index 4102f02b599..d70ee0de5bd 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -260,7 +260,7 @@ def astra_cpu_back_projector(proj_data, geometry, vol_space, out=None, class AstraCpuImpl: - """Thin wrapper around the ASTRA CPU implementations.""" + """Thin wrapper implementing ASTRA CPU for `RayTransform`.""" def __init__(self, geometry, vol_space, proj_space): """Initialize a new instance. diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 4a9c10a6552..89f668f6c7d 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -39,7 +39,7 @@ class AstraCudaImpl: - """Thin wrapper around ASTRA.""" + """`RayTransform` implementation for CUDA algorithms in ASTRA.""" algo_forward_id = None algo_backward_id = None diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index 793361710fd..fb57284dbcf 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -180,6 +180,8 @@ def skimage_radon_back_projector(sinogram, geometry, vol_space, out=None): class SkimageImpl: + """Scikit-image backend of the `RayTransform` operator.""" + def __init__(self, geometry, vol_space, proj_space): """Initialize a new instance. diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index 467e2508f5d..dfb45bdd472 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -2,7 +2,23 @@ def _add_default_complex_impl(fn): - """Wrapper to call a class method twice when it is complex.""" + """Wrapper to call a class method twice when it is complex. + + This function helps `RayTransform` implementations by splitting + complex-valued forward/backward function calls into two real-valued calls. + + Parameters + ---------- + fn : Callable + Function with signature ``fn(self, x, out=None, **kwargs)``. + ``self`` must be an object instance having ``self.vol_space`` and + ``self.proj_space``. These spaces must be both real or both complex. + + Returns + ------- + Callable + A wrapper function to be executed by a Python decorator. + """ def wrapper(self, x, out=None, **kwargs): if self.vol_space.is_real and self.proj_space.is_real: diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 2040dcdaf0e..509f3be2102 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -207,6 +207,7 @@ def __init__(self, vol_space, geometry, **kwargs): @staticmethod def _check_impl(impl): + """Internal method to verify the validity of the `impl` kwarg.""" impl_instance = None if impl is None: # User didn't specify a backend @@ -256,12 +257,24 @@ def _check_impl(impl): @property def impl(self): - """Implementation name string or type.""" + """Implementation name string. + + If a custom ``impl`` was provided this method returns a ``str`` + of the type.""" return self.__impl def create_impl(self, use_cache=True): - """Fetches or instantiates implementation backend for evaluation.""" + """Fetches or instantiates implementation backend for evaluation. + + Parameters + ---------- + bool : use_cache + If ``True`` returns the cached implementation backend, if it + was generated in a previous call (or given with ``__init__``). + If ``False`` a new instance of the backend will be generated, + freeing up GPU memory and RAM used by the backend. + """ # Use impl creation (__cached_impl) when `use_cache` is True if not use_cache or self.__cached_impl is None: @@ -274,7 +287,24 @@ def create_impl(self, use_cache=True): return self.__cached_impl def _call(self, x, out=None, **kwargs): - """Real-space forward projection for the current set-up.""" + """Forward projection + + Parameters + ---------- + x : DiscreteLpElement + A volume. Must be an element from `RayTransform.domain`. + out : `RayTransform.range` element, optional + A DiscreteLpElement from `RayTransform.range`, to which the result + of the operator evaluation is written. + **kwargs + Arbitrary keyword arguments, passed on to the implementation + backend. + + Returns + ------- + DiscreteLpElement + A sinogram, from the projection space `RayTransform.range`. + """ return self.create_impl(self.use_cache) \ .call_forward(x, out, **kwargs) @@ -287,6 +317,10 @@ def geometry(self): def adjoint(self): """Adjoint of this operator. + The adjoint of the `RayTransform` is the linear `RayBackProjection` + operator. It uses the same geometry and shares the implementation + backend whenever `RayTransform.use_cache` is `True`. + Returns ------- adjoint : `RayBackProjection` @@ -299,7 +333,29 @@ class RayBackProjection(Operator): """Adjoint of the discrete Ray transform between L^p spaces.""" def _call(self, x, out=None, **kwargs): - """Back-projection for the current set-up.""" + """Backprojection + + Parameters + ---------- + x : DiscreteLpElement + A sinogram. Must be an element from + `RayTransform.range` (domain of `RayBackProjection`). + out : `RayBackProjection.domain` element, optional + A volume (an element from the volume space + `RayTransform.domain` or, equivalently, + `RayBackprojection.range`). The result of this + evaluation is written into this variable. + **kwargs + Arbitrary keyword arguments, passed on to the + implementation backend. + + Returns + ------- + DiscreteLpElement + The reconstructed volume, an element in + `RayBackProjection.range`. + """ + return ray_trafo.create_impl(ray_trafo.use_cache) \ .call_backward(x, out, **kwargs) From 3bd56404cb149aeed4f1ecb18582a93817c3e298 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 27 Feb 2020 22:26:58 +0100 Subject: [PATCH 26/39] Update test to include complex adjoint of `RayTransform` --- odl/test/tomo/operators/ray_trafo_test.py | 9 +++++++++ odl/tomo/backends/util.py | 9 ++++++--- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/odl/test/tomo/operators/ray_trafo_test.py b/odl/test/tomo/operators/ray_trafo_test.py index dd7af502cd5..5dad840dc38 100644 --- a/odl/test/tomo/operators/ray_trafo_test.py +++ b/odl/test/tomo/operators/ray_trafo_test.py @@ -339,6 +339,15 @@ def test_complex(impl): assert all_almost_equal(data.real, true_data_re) assert all_almost_equal(data.imag, true_data_im) + # test adjoint for complex data + backproj_r = ray_trafo_r.adjoint + backproj_c = ray_trafo_c.adjoint + true_vol_re = backproj_r(data.real) + true_vol_im = backproj_r(data.imag) + backproj_vol = backproj_c(data) + + assert all_almost_equal(backproj_vol.real, true_vol_re) + assert all_almost_equal(backproj_vol.imag, true_vol_im) def test_anisotropic_voxels(geometry): """Test projection and backprojection with anisotropic voxels.""" diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index dfb45bdd472..46ba9a55144 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -30,9 +30,12 @@ def wrapper(self, x, out=None, **kwargs): ] if out is None: - out = self.proj_space.element() - out.real = result_parts[0] - out.imag = result_parts[1] + range = self.proj_space if x in self.vol_space \ + else self.vol_space + out = range.element() + + out.real = result_parts[0] + out.imag = result_parts[1] return out else: raise RuntimeError('Domain and range need to be both real ' From 2aab34ea1b582382b0000b55ca914d6310eca412 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 27 Feb 2020 22:29:15 +0100 Subject: [PATCH 27/39] Make `_IMPL_STR2TYPE` public as `RAY_TRAFO_IMPLS` --- odl/tomo/backends/skimage_radon.py | 2 +- odl/tomo/operators/ray_trafo.py | 40 +++++++++++++++--------------- 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index fb57284dbcf..7d5bb7af805 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -179,7 +179,7 @@ def skimage_radon_back_projector(sinogram, geometry, vol_space, out=None): return out -class SkimageImpl: +class SkImageImpl: """Scikit-image backend of the `RayTransform` operator.""" def __init__(self, geometry, vol_space, proj_space): diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 509f3be2102..7c5729c8fbf 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -10,6 +10,8 @@ from __future__ import absolute_import, division, print_function +from collections import OrderedDict + import numpy as np from odl.discr import DiscretizedSpace @@ -19,23 +21,21 @@ ASTRA_AVAILABLE, ASTRA_CUDA_AVAILABLE, SKIMAGE_AVAILABLE) from odl.tomo.backends.astra_cpu import AstraCpuImpl from odl.tomo.backends.astra_cuda import AstraCudaImpl -from odl.tomo.backends.skimage_radon import SkimageImpl +from odl.tomo.backends.skimage_radon import SkImageImpl from odl.tomo.geometry import Geometry -# Backends that are implemented in ODL and can be specified via `impl` -_IMPL_STR2TYPE = { - 'astra_cpu': AstraCpuImpl, - 'astra_cuda': AstraCudaImpl, - 'skimage': SkimageImpl} - -# Backends that are installed, ordered by preference -_AVAILABLE_IMPLS = [] -if ASTRA_CUDA_AVAILABLE: - _AVAILABLE_IMPLS.append('astra_cuda') -if ASTRA_AVAILABLE: - _AVAILABLE_IMPLS.append('astra_cpu') +# RAY_TRAFO_IMPLS are used by `RayTransform` when no `impl` is given. +# The last inserted implementation has highest priority. +RAY_TRAFO_IMPLS = OrderedDict() if SKIMAGE_AVAILABLE: - _AVAILABLE_IMPLS.append('skimage') + RAY_TRAFO_IMPLS['skimage'] = SkImageImpl +if ASTRA_AVAILABLE: + RAY_TRAFO_IMPLS['astra_cpu'] = AstraCpuImpl +if ASTRA_CUDA_AVAILABLE: + RAY_TRAFO_IMPLS['astra_cuda'] = AstraCudaImpl + +# Used to warn the user when a backend is not available/installed. +_ALL_IMPLS = ('astra_cuda', 'astra_cpu', 'skimage') __all__ = ('RayTransform',) @@ -211,25 +211,25 @@ def _check_impl(impl): impl_instance = None if impl is None: # User didn't specify a backend - if not _AVAILABLE_IMPLS: + if not RAY_TRAFO_IMPLS: raise RuntimeError('no ray transform back-end available; ' 'this requires 3rd party packages, please ' 'check the install docs') - # Select fastest available, _AVAILABLE_IMPLS is sorted by speed - impl_type = _IMPL_STR2TYPE[_AVAILABLE_IMPLS[0]] + # Select fastest available + impl_type = next(reversed(RAY_TRAFO_IMPLS.values())) else: # User did specify `impl` if isinstance(impl, str): - if impl.lower() not in _IMPL_STR2TYPE: + if impl.lower() not in _ALL_IMPLS: raise ValueError('`impl` {!r} not understood'.format(impl)) - if impl.lower() not in _AVAILABLE_IMPLS: + if impl.lower() not in RAY_TRAFO_IMPLS.keys(): raise ValueError( '{!r} back-end not available'.format(impl)) - impl_type = _IMPL_STR2TYPE[impl.lower()] + impl_type = RAY_TRAFO_IMPLS[impl.lower()] elif isinstance(impl, type) or isinstance(impl, object): # User gave the type and leaves instantiation to us forward = getattr(impl, "call_forward", None) From f952b9e8c5f0066bed911eb56848456865ae67e5 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 27 Feb 2020 23:03:11 +0100 Subject: [PATCH 28/39] Do not reassign to `out` when `out` is not None --- odl/tomo/backends/util.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index 46ba9a55144..6e89b24cb95 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -33,9 +33,9 @@ def wrapper(self, x, out=None, **kwargs): range = self.proj_space if x in self.vol_space \ else self.vol_space out = range.element() + out.real = result_parts[0] + out.imag = result_parts[1] - out.real = result_parts[0] - out.imag = result_parts[1] return out else: raise RuntimeError('Domain and range need to be both real ' From fb2d29d7e6e3ed9e02b89b9d6912d28c19ee14c0 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Mon, 13 Apr 2020 14:34:04 +0200 Subject: [PATCH 29/39] Some formatting updates - Opening and closing parens on separate lines if line needs split - No more backslash line continuation - Sorted imports - Removal of remaining references to `DiscreteLp` --- odl/tomo/backends/__init__.py | 2 +- odl/tomo/backends/astra_cpu.py | 103 ++++++++++++++++++----------- odl/tomo/backends/astra_cuda.py | 94 ++++++++++++++++---------- odl/tomo/backends/skimage_radon.py | 66 ++++++++++-------- odl/tomo/backends/util.py | 12 ++-- odl/tomo/operators/ray_trafo.py | 53 +++++++++------ 6 files changed, 201 insertions(+), 129 deletions(-) diff --git a/odl/tomo/backends/__init__.py b/odl/tomo/backends/__init__.py index 6b3b9391f77..b621dbc4a72 100644 --- a/odl/tomo/backends/__init__.py +++ b/odl/tomo/backends/__init__.py @@ -13,8 +13,8 @@ from .astra_cpu import * from .astra_cuda import * from .astra_setup import * -from .util import * from .skimage_radon import * +from .util import * __all__ = () __all__ += astra_cpu.__all__ diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index d70ee0de5bd..3062433dd4c 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -10,9 +10,10 @@ from __future__ import absolute_import, division, print_function -import numpy as np import warnings +import numpy as np + from odl.discr import DiscretizedSpace, DiscretizedSpaceElement from odl.tomo.backends.astra_setup import ( astra_algorithm, astra_data, astra_projection_geometry, astra_projector, @@ -96,24 +97,34 @@ def astra_cpu_forward_projector(vol_data, geometry, proj_space, out=None, If ``out`` was provided, the returned object is a reference to it. """ if not isinstance(vol_data, DiscretizedSpaceElement): - raise TypeError('volume data {!r} is not a `DiscretizedSpaceElement` ' - 'instance.'.format(vol_data)) + raise TypeError( + 'volume data {!r} is not a `DiscretizedSpaceElement` instance' + ''.format(vol_data) + ) if vol_data.space.impl != 'numpy': - raise TypeError("`vol_data.space.impl` must be 'numpy', got {!r}" - "".format(vol_data.space.impl)) + raise TypeError( + "`vol_data.space.impl` must be 'numpy', got {!r}" + "".format(vol_data.space.impl) + ) if not isinstance(geometry, Geometry): - raise TypeError('geometry {!r} is not a Geometry instance' - ''.format(geometry)) + raise TypeError( + 'geometry {!r} is not a Geometry instance'.format(geometry) + ) if not isinstance(proj_space, DiscretizedSpace): - raise TypeError('`proj_space` {!r} is not a DiscretizedSpace ' - 'instance.'.format(proj_space)) + raise TypeError( + '`proj_space` {!r} is not a DiscretizedSpace instance.' + ''.format(proj_space) + ) if proj_space.impl != 'numpy': - raise TypeError("`proj_space.impl` must be 'numpy', got {!r}" - "".format(proj_space.impl)) + raise TypeError( + "`proj_space.impl` must be 'numpy', got {!r}" + "".format(proj_space.impl) + ) if vol_data.ndim != geometry.ndim: - raise ValueError('dimensions {} of volume data and {} of geometry ' - 'do not match' - ''.format(vol_data.ndim, geometry.ndim)) + raise ValueError( + 'dimensions {} of volume data and {} of geometry do not match' + ''.format(vol_data.ndim, geometry.ndim) + ) if out is None: out = proj_space.element() else: @@ -192,18 +203,23 @@ def astra_cpu_back_projector(proj_data, geometry, vol_space, out=None, 'instance'.format(proj_data) ) if proj_data.space.impl != 'numpy': - raise TypeError('`proj_data` must be a `numpy.ndarray` based, ' - "container got `impl` {!r}" - "".format(proj_data.space.impl)) + raise TypeError( + '`proj_data` must be a `numpy.ndarray` based, container, ' + "got `impl` {!r}".format(proj_data.space.impl) + ) if not isinstance(geometry, Geometry): - raise TypeError('geometry {!r} is not a Geometry instance' - ''.format(geometry)) + raise TypeError( + 'geometry {!r} is not a Geometry instance'.format(geometry) + ) if not isinstance(vol_space, DiscretizedSpace): - raise TypeError('volume space {!r} is not a DiscretizedSpace ' - 'instance'.format(vol_space)) + raise TypeError( + 'volume space {!r} is not a DiscretizedSpace instance' + ''.format(vol_space) + ) if vol_space.impl != 'numpy': - raise TypeError("`vol_space.impl` must be 'numpy', got {!r}" - "".format(vol_space.impl)) + raise TypeError( + "`vol_space.impl` must be 'numpy', got {!r}".format(vol_space.impl) + ) if vol_space.ndim != geometry.ndim: raise ValueError( 'dimensions {} of reconstruction space and {} of geometry ' @@ -226,8 +242,9 @@ def astra_cpu_back_projector(proj_data, geometry, vol_space, out=None, proj_geom = astra_projection_geometry(geometry) # Create ASTRA data structure - sino_id = astra_data(proj_geom, datatype='projection', data=proj_data, - allow_copy=True) + sino_id = astra_data( + proj_geom, datatype='projection', data=proj_data, allow_copy=True + ) # Create projector if astra_proj_type is None: @@ -236,11 +253,13 @@ def astra_cpu_back_projector(proj_data, geometry, vol_space, out=None, # Convert out to correct dtype and order if needed. with writable_array(out, dtype='float32', order='C') as out_arr: - vol_id = astra_data(vol_geom, datatype='volume', data=out_arr, - ndim=vol_space.ndim) + vol_id = astra_data( + vol_geom, datatype='volume', data=out_arr, ndim=vol_space.ndim + ) # Create algorithm - algo_id = astra_algorithm('backward', ndim, vol_id, sino_id, proj_id, - impl='cpu') + algo_id = astra_algorithm( + 'backward', ndim, vol_id, sino_id, proj_id, impl='cpu' + ) # Run algorithm astra.algorithm.run(algo_id) @@ -276,30 +295,32 @@ def __init__(self, geometry, vol_space, proj_space): Projection space, the space of the result. """ if not isinstance(geometry, Geometry): - raise TypeError('`geometry` must be a `Geometry` instance, got ' - '{!r}'.format(geometry)) - + raise TypeError( + '`geometry` must be a `Geometry` instance, got {!r}' + ''.format(geometry) + ) if not isinstance(vol_space, DiscretizedSpace): raise TypeError( '`vol_space` must be a `DiscretizedSpace` instance, got {!r}' ''.format(vol_space) ) - if not isinstance(proj_space, DiscretizedSpace): raise TypeError( '`proj_space` must be a `DiscretizedSpace` instance, got {!r}' ''.format(proj_space) ) + if geometry.ndim > 2: + raise ValueError( + '`impl` {!r} only works for 2d'.format(self.__name__) + ) if vol_space.size >= 512 ** 2: warnings.warn( "The 'astra_cpu' backend may be too slow for volumes of this " "size. Consider using 'astra_cuda' if your machine has an " - "Nvidia GPU.", RuntimeWarning) - - if geometry.ndim > 2: - raise ValueError('`impl` {!r} only works for 2d' - ''.format(self.__name__)) + "Nvidia GPU.", + RuntimeWarning, + ) self.geometry = geometry self._vol_space = vol_space @@ -316,12 +337,14 @@ def proj_space(self): @_add_default_complex_impl def call_backward(self, x, out, **kwargs): return astra_cpu_back_projector( - x, self.geometry, self.vol_space.real_space, out, **kwargs) + x, self.geometry, self.vol_space.real_space, out, **kwargs + ) @_add_default_complex_impl def call_forward(self, x, out, **kwargs): return astra_cpu_forward_projector( - x, self.geometry, self.proj_space.real_space, out, **kwargs) + x, self.geometry, self.proj_space.real_space, out, **kwargs + ) if __name__ == '__main__': diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 89f668f6c7d..499f94c1f4e 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -10,18 +10,18 @@ from __future__ import absolute_import, division, print_function +import warnings from multiprocessing import Lock import numpy as np -import warnings from packaging.version import parse as parse_version from odl.discr import DiscretizedSpace -from odl.tomo.backends.util import _add_default_complex_impl from odl.tomo.backends.astra_setup import ( ASTRA_VERSION, astra_algorithm, astra_data, astra_projection_geometry, - astra_projector, astra_volume_geometry, astra_supports, - astra_versions_supporting) + astra_projector, astra_supports, astra_versions_supporting, + astra_volume_geometry) +from odl.tomo.backends.util import _add_default_complex_impl from odl.tomo.geometry import ( ConeBeamGeometry, FanBeamGeometry, Geometry, Parallel2dGeometry, Parallel3dAxisGeometry) @@ -61,15 +61,15 @@ def __init__(self, geometry, vol_space, proj_space): Projection space, the space of the result. """ if not isinstance(geometry, Geometry): - raise TypeError('`geometry` must be a `Geometry` instance, got ' - '{!r}'.format(geometry)) - + raise TypeError( + '`geometry` must be a `Geometry` instance, got {!r}' + ''.format(geometry) + ) if not isinstance(vol_space, DiscretizedSpace): raise TypeError( '`vol_space` must be a `DiscreteLP` instance, got {!r}' ''.format(vol_space) ) - if not isinstance(proj_space, DiscretizedSpace): raise TypeError( '`proj_space` must be a `DiscreteLP` instance, got {!r}' @@ -79,10 +79,13 @@ def __init__(self, geometry, vol_space, proj_space): # Print a warning if the detector midpoint normal vector at any # angle is perpendicular to the geometry axis in parallel 3d # single-axis geometry -- this is broken in some ASTRA versions - if (isinstance(geometry, Parallel3dAxisGeometry) and - not astra_supports('par3d_det_mid_pt_perp_to_axis')): + if ( + isinstance(geometry, Parallel3dAxisGeometry) + and not astra_supports('par3d_det_mid_pt_perp_to_axis') + ): req_ver = astra_versions_supporting( - 'par3d_det_mid_pt_perp_to_axis') + 'par3d_det_mid_pt_perp_to_axis' + ) axis = geometry.axis mid_pt = geometry.det_params.mid_pt for i, angle in enumerate(geometry.angles): @@ -144,32 +147,46 @@ def create_ids(self): # Create ASTRA data structures vol_geom = astra_volume_geometry(self.vol_space) proj_geom = astra_projection_geometry(self.geometry) - self.vol_id = astra_data(vol_geom, - datatype='volume', - ndim=self.vol_space.ndim, - data=self.vol_array, - allow_copy=False) + self.vol_id = astra_data( + vol_geom, + datatype='volume', + ndim=self.vol_space.ndim, + data=self.vol_array, + allow_copy=False, + ) proj_type = 'cuda' if proj_ndim == 2 else 'cuda3d' self.proj_id = astra_projector( proj_type, vol_geom, proj_geom, proj_ndim ) - self.sino_id = astra_data(proj_geom, - datatype='projection', - ndim=proj_ndim, - data=self.proj_array, - allow_copy=False) + self.sino_id = astra_data( + proj_geom, + datatype='projection', + ndim=proj_ndim, + data=self.proj_array, + allow_copy=False, + ) # Create algorithm self.algo_forward_id = astra_algorithm( - 'forward', proj_ndim, self.vol_id, self.sino_id, - proj_id=self.proj_id, impl='cuda') + 'forward', + proj_ndim, + self.vol_id, + self.sino_id, + self.proj_id, + impl='cuda', + ) # Create algorithm self.algo_backward_id = astra_algorithm( - 'backward', proj_ndim, self.vol_id, self.sino_id, - proj_id=self.proj_id, impl='cuda') + 'backward', + proj_ndim, + self.vol_id, + self.sino_id, + self.proj_id, + impl='cuda', + ) @_add_default_complex_impl def call_forward(self, x, out=None, **kwargs): @@ -220,8 +237,10 @@ def _call_forward_real(self, vol_data, out=None, **kwargs): self.proj_space.shape) # Fix scaling to weight by pixel size - if (isinstance(self.geometry, Parallel2dGeometry) and - parse_version(ASTRA_VERSION) < parse_version('1.9.9.dev')): + if ( + isinstance(self.geometry, Parallel2dGeometry) + and parse_version(ASTRA_VERSION) < parse_version('1.9.9.dev') + ): # parallel2d scales with pixel stride out *= 1 / float(self.geometry.det_partition.cell_volume) @@ -265,7 +284,8 @@ def _call_backward_real(self, proj_data, out=None, **kwargs): shape = (-1,) + self.geometry.det_partition.shape reshaped_proj_data = proj_data.asarray().reshape(shape) swapped_proj_data = np.ascontiguousarray( - np.swapaxes(reshaped_proj_data, 0, 1)) + np.swapaxes(reshaped_proj_data, 0, 1) + ) astra.data3d.store(self.sino_id, swapped_proj_data) # Run algorithm @@ -276,7 +296,8 @@ def _call_backward_real(self, proj_data, out=None, **kwargs): # Fix scaling to weight by pixel/voxel size out *= astra_cuda_bp_scaling_factor( - self.proj_space, self.vol_space, self.geometry) + self.proj_space, self.vol_space, self.geometry + ) return out @@ -328,10 +349,12 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): proj_size = float(proj_space.partition.size) proj_weighting = proj_extent / proj_size - scaling_factor *= (proj_space.weighting.const / - proj_weighting) - scaling_factor /= (vol_space.weighting.const / - vol_space.cell_volume) + scaling_factor *= ( + proj_space.weighting.const / proj_weighting + ) + scaling_factor /= ( + vol_space.weighting.const / vol_space.cell_volume + ) if parse_version(ASTRA_VERSION) < parse_version('1.8rc1'): # Scaling for the old, pre-1.8 behaviour @@ -392,8 +415,9 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): # missing factor src_radius ** 2 in the ASTRA BP with # density weighting. det_px_area = geometry.det_partition.cell_volume - scaling_factor *= (src_radius ** 2 * det_px_area ** 2 / - vol_space.cell_volume ** 2) + scaling_factor *= ( + src_radius ** 2 * det_px_area ** 2 / vol_space.cell_volume ** 2 + ) elif parse_version(ASTRA_VERSION) < parse_version('1.9.9.dev'): # Scaling for intermediate dev releases between 1.8.3 and 1.9.9.dev if isinstance(geometry, Parallel2dGeometry): diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index 7d5bb7af805..cf7d2f1e5e1 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -10,14 +10,14 @@ from __future__ import division -import numpy as np import warnings -from odl.tomo import Parallel2dGeometry, Geometry +import numpy as np from odl.discr import ( - uniform_discr_frompartition, uniform_partition, DiscretizedSpace) + DiscretizedSpace, uniform_discr_frompartition, uniform_partition) from odl.discr.discr_utils import linear_interpolator, point_collocation +from odl.tomo import Geometry, Parallel2dGeometry from odl.tomo.backends.util import _add_default_complex_impl from odl.util.utility import writable_array @@ -189,52 +189,58 @@ def __init__(self, geometry, vol_space, proj_space): ---------- geometry : `Geometry` Geometry defining the tomographic setup. - vol_space : `DiscreteLp` + vol_space : `DiscretizedSpace` Reconstruction space, the space of the images to be forward projected. - proj_space : `DiscreteLp` + proj_space : `DiscretizedSpace` Projection space, the space of the result. """ if not isinstance(geometry, Geometry): - raise TypeError('`geometry` must be a `Geometry` instance, got ' - '{!r}'.format(geometry)) - + raise TypeError( + '`geometry` must be a `Geometry` instance, got {!r}' + ''.format(geometry) + ) if not isinstance(vol_space, DiscretizedSpace): raise TypeError( '`vol_space` must be a `DiscretizedSpace` instance, got {!r}' ''.format(vol_space) ) - if not isinstance(proj_space, DiscretizedSpace): raise TypeError( '`proj_space` must be a `DiscretizedSpace` instance, got {!r}' ''.format(proj_space) ) - - if vol_space.size >= 256 ** 2: - warnings.warn( - "The 'skimage' backend may be too slow for volumes of this " - "size. Consider using 'astra_cpu', or 'astra_cuda' if your " - "machine has an Nvidia GPU.", RuntimeWarning) - if not isinstance(geometry, Parallel2dGeometry): - raise TypeError("{!r} backend only supports 2d parallel " - 'geometries'.format(self.__name__)) - + raise TypeError( + "{!r} backend only supports 2d parallel geometries" + ''.format(self.__name__) + ) mid_pt = vol_space.domain.mid_pt if not np.allclose(mid_pt, [0, 0]): - raise ValueError('Reconstruction space must be centered at (0, 0),' - 'got midpoint {}'.format(mid_pt)) - + raise ValueError( + 'reconstruction space must be centered at (0, 0), ' + 'got midpoint {}'.format(mid_pt) + ) shape = vol_space.shape if shape[0] != shape[1]: - raise ValueError('`vol_space.shape` must have equal entries, ' - 'got {}'.format(shape)) - + raise ValueError( + '`vol_space.shape` must have equal entries, got {}' + ''.format(shape) + ) extent = vol_space.domain.extent if extent[0] != extent[1]: - raise ValueError('`vol_space.extent` must have equal entries, ' - 'got {}'.format(extent)) + raise ValueError( + '`vol_space.extent` must have equal entries, got {}' + ''.format(extent) + ) + + if vol_space.size >= 256 ** 2: + warnings.warn( + "The 'skimage' backend may be too slow for volumes of this " + "size. Consider using 'astra_cpu', or 'astra_cuda' if your " + "machine has an Nvidia GPU.", + RuntimeWarning, + ) self.geometry = geometry self._vol_space = vol_space @@ -251,9 +257,11 @@ def proj_space(self): @_add_default_complex_impl def call_forward(self, x, out, **kwargs): return skimage_radon_forward_projector( - x, self.geometry, self.proj_space.real_space, out) + x, self.geometry, self.proj_space.real_space, out + ) @_add_default_complex_impl def call_backward(self, x, out, **kwargs): return skimage_radon_back_projector( - x, self.geometry, self.vol_space.real_space, out) + x, self.geometry, self.vol_space.real_space, out + ) diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index 6e89b24cb95..2a0cf01487b 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -30,15 +30,19 @@ def wrapper(self, x, out=None, **kwargs): ] if out is None: - range = self.proj_space if x in self.vol_space \ - else self.vol_space + if x in self.vol_space: + range = self.proj_space + else: + range = self.vol_space + out = range.element() out.real = result_parts[0] out.imag = result_parts[1] return out else: - raise RuntimeError('Domain and range need to be both real ' - 'or both complex.') + raise RuntimeError( + 'domain and range need to be both real or both complex' + ) return wrapper diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 7c5729c8fbf..a8482097eb1 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -93,8 +93,10 @@ def __init__(self, vol_space, geometry, **kwargs): '{!r}'.format(vol_space)) if not isinstance(geometry, Geometry): - raise TypeError('`geometry` must be a `Geometry` instance, got ' - '{!r}'.format(geometry)) + raise TypeError( + '`geometry` must be a `Geometry` instance, got {!r}' + ''.format(geometry) + ) # Generate or check projection space proj_space = kwargs.pop('proj_space', None) @@ -105,7 +107,9 @@ def __init__(self, vol_space, geometry, **kwargs): weighting = None elif ( isinstance(vol_space.weighting, ConstWeighting) - and np.isclose(vol_space.weighting.const, vol_space.cell_volume) + and np.isclose( + vol_space.weighting.const, vol_space.cell_volume + ) ): # Approximate cell volume # TODO: find a way to treat angles and detector differently @@ -120,9 +124,11 @@ def __init__(self, vol_space, geometry, **kwargs): else: raise NotImplementedError('unknown weighting of domain') - proj_tspace = vol_space.tspace_type(geometry.partition.shape, - weighting=weighting, - dtype=dtype) + proj_tspace = vol_space.tspace_type( + geometry.partition.shape, + weighting=weighting, + dtype=dtype, + ) if geometry.motion_partition.ndim == 0: angle_labels = [] @@ -188,8 +194,10 @@ def __init__(self, vol_space, geometry, **kwargs): impl = kwargs.pop('impl', None) impl_type, self.__cached_impl = self._check_impl(impl) self._impl_type = impl_type - self.__impl = impl.lower() \ - if isinstance(impl, str) else impl_type.__name__ + if isinstance(impl, str): + self.__impl = impl.lower() + else: + self.__impl = impl_type.__name__ self._geometry = geometry @@ -212,9 +220,10 @@ def _check_impl(impl): if impl is None: # User didn't specify a backend if not RAY_TRAFO_IMPLS: - raise RuntimeError('no ray transform back-end available; ' - 'this requires 3rd party packages, please ' - 'check the install docs') + raise RuntimeError( + 'no ray transform back-end available; this requires ' + '3rd party packages, please check the install docs' + ) # Select fastest available impl_type = next(reversed(RAY_TRAFO_IMPLS.values())) @@ -227,7 +236,8 @@ def _check_impl(impl): if impl.lower() not in RAY_TRAFO_IMPLS.keys(): raise ValueError( - '{!r} back-end not available'.format(impl)) + '{!r} back-end not available'.format(impl) + ) impl_type = RAY_TRAFO_IMPLS[impl.lower()] elif isinstance(impl, type) or isinstance(impl, object): @@ -236,8 +246,10 @@ def _check_impl(impl): backward = getattr(impl, "call_backward", None) if not callable(forward) and not callable(backward): - raise TypeError('Type {!r} must be have a `call_forward` ' - 'or `call_backward`.'.format(impl)) + raise TypeError( + 'type {!r} must be have a `call_forward` ' + 'or `call_backward`.'.format(impl) + ) if isinstance(impl, type): impl_type = impl @@ -306,8 +318,7 @@ def _call(self, x, out=None, **kwargs): A sinogram, from the projection space `RayTransform.range`. """ - return self.create_impl(self.use_cache) \ - .call_forward(x, out, **kwargs) + return self.create_impl(self.use_cache).call_forward(x, out, **kwargs) @property def geometry(self): @@ -356,8 +367,9 @@ def _call(self, x, out=None, **kwargs): `RayBackProjection.range`. """ - return ray_trafo.create_impl(ray_trafo.use_cache) \ - .call_backward(x, out, **kwargs) + return ray_trafo.create_impl( + ray_trafo.use_cache + ).call_backward(x, out, **kwargs) @property def geometry(self): @@ -369,8 +381,9 @@ def adjoint(self): kwargs = self._extra_kwargs.copy() kwargs['domain'] = self.range - self._adjoint = RayBackProjection(range=self.domain, linear=True, - **kwargs) + self._adjoint = RayBackProjection( + range=self.domain, linear=True, **kwargs + ) return self._adjoint From 5c7dad5d13e7e9e84a6d274571867f60def4efc7 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Mon, 13 Apr 2020 14:34:39 +0200 Subject: [PATCH 30/39] Update of README.md in tomo subpackage --- odl/tomo/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/odl/tomo/README.md b/odl/tomo/README.md index 7df9e96565b..263aa203ca4 100644 --- a/odl/tomo/README.md +++ b/odl/tomo/README.md @@ -7,5 +7,5 @@ This directory contains all of the source code related tomographic reconstructio * [analytic](analytic) Analytic reconstruction methods such as filtered back-projection. Also contains various utilities like `parker_weighting`. * [backends](backends) Bindings to external libraries. * [geometry](geometry) Definitions of projection geometries. -* [operators](operators) Defines the `RayTransform` operator and its adjoint ("back-projection"). +* [operators](operators) Defines the `RayTransform` operator. * [util](util) Utilities used internally. \ No newline at end of file From 2fc5d60670a97e31932322769dcef987e403e369 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Mon, 13 Apr 2020 14:35:10 +0200 Subject: [PATCH 31/39] Copyright notice and functool.wraps in backend utils --- odl/tomo/backends/util.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index 2a0cf01487b..fd8ff879e06 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -1,4 +1,16 @@ -__all__ = tuple() +# Copyright 2014-2020 The ODL contributors +# +# This file is part of ODL. +# +# This Source Code Form is subject to the terms of the Mozilla Public License, +# v. 2.0. If a copy of the MPL was not distributed with this file, You can +# obtain one at https://mozilla.org/MPL/2.0/. + +"""Utilities for tomo backends.""" + +from functools import wraps + +__all__ = () def _add_default_complex_impl(fn): @@ -20,6 +32,7 @@ def _add_default_complex_impl(fn): A wrapper function to be executed by a Python decorator. """ + @wraps(fn) def wrapper(self, x, out=None, **kwargs): if self.vol_space.is_real and self.proj_space.is_real: return fn(self, x, out, **kwargs) From d9cfc62712672cc4ce0d4c99c4ae74a1e08371cf Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Mon, 13 Apr 2020 14:38:48 +0200 Subject: [PATCH 32/39] Fix failing import in skimage_radon.py --- odl/tomo/backends/skimage_radon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index cf7d2f1e5e1..fa4ef8d38c7 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -17,7 +17,7 @@ from odl.discr import ( DiscretizedSpace, uniform_discr_frompartition, uniform_partition) from odl.discr.discr_utils import linear_interpolator, point_collocation -from odl.tomo import Geometry, Parallel2dGeometry +from odl.tomo.geometry import Geometry, Parallel2dGeometry from odl.tomo.backends.util import _add_default_complex_impl from odl.util.utility import writable_array From 60a18cbc1e35b2673202690684aec276589a47c1 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 16 Apr 2020 15:03:53 +0200 Subject: [PATCH 33/39] Change of words --- odl/tomo/backends/astra_cuda.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 499f94c1f4e..be8dc1d75ca 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -108,7 +108,7 @@ def __init__(self, geometry, vol_space, proj_space): self.create_ids() - # ASTRA projectors are not thread-safe, thus we need to lock ourselves + # ASTRA projectors are not thread-safe, thus we need to lock manually self._mutex = Lock() @property From 217e33a6eb26809a382de37e18f9f60f2f0b635f Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 16 Apr 2020 15:06:31 +0200 Subject: [PATCH 34/39] Docstring extended with use case --- odl/tomo/backends/util.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index fd8ff879e06..d6bec82f7a4 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -14,17 +14,27 @@ def _add_default_complex_impl(fn): - """Wrapper to call a class method twice when it is complex. + """Wrapper to call a real-valued class method twice when input is complex. - This function helps `RayTransform` implementations by splitting - complex-valued forward/backward function calls into two real-valued calls. + This function helps `RayTransform` implementations with extending methods + that work for real-valued elements to complex-valued elements, by splitting + complex calls into two individual real calls. + + The wrapper will only work for methods of which the class provides a + `vol_space` and `proj_space`. Typically, this will then work as a decorator + on the method, e.g.: + ``` + @_add_default_complex_impl + def call_forward(self, x, out=None, **kwargs): + # Code that should run for real input and output + ``` Parameters ---------- fn : Callable Function with signature ``fn(self, x, out=None, **kwargs)``. ``self`` must be an object instance having ``self.vol_space`` and - ``self.proj_space``. These spaces must be both real or both complex. + ``self.proj_space``. Returns ------- From 204ec7fcfd484a3889751b7745fc37900051b900 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 16 Apr 2020 15:10:24 +0200 Subject: [PATCH 35/39] Removed _ALL_IMPLS, enhanced exception messages, allow duck-typing `impl` string. --- odl/tomo/operators/ray_trafo.py | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index a8482097eb1..7c5a796fb0f 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -23,6 +23,7 @@ from odl.tomo.backends.astra_cuda import AstraCudaImpl from odl.tomo.backends.skimage_radon import SkImageImpl from odl.tomo.geometry import Geometry +from odl.util import is_string # RAY_TRAFO_IMPLS are used by `RayTransform` when no `impl` is given. # The last inserted implementation has highest priority. @@ -34,9 +35,6 @@ if ASTRA_CUDA_AVAILABLE: RAY_TRAFO_IMPLS['astra_cuda'] = AstraCudaImpl -# Used to warn the user when a backend is not available/installed. -_ALL_IMPLS = ('astra_cuda', 'astra_cpu', 'skimage') - __all__ = ('RayTransform',) @@ -194,7 +192,7 @@ def __init__(self, vol_space, geometry, **kwargs): impl = kwargs.pop('impl', None) impl_type, self.__cached_impl = self._check_impl(impl) self._impl_type = impl_type - if isinstance(impl, str): + if is_string(impl): self.__impl = impl.lower() else: self.__impl = impl_type.__name__ @@ -221,8 +219,8 @@ def _check_impl(impl): if impl is None: # User didn't specify a backend if not RAY_TRAFO_IMPLS: raise RuntimeError( - 'no ray transform back-end available; this requires ' - '3rd party packages, please check the install docs' + 'No `RayTransform` back-end available; this requires ' + '3rd party packages, please check the install docs.' ) # Select fastest available @@ -230,13 +228,14 @@ def _check_impl(impl): else: # User did specify `impl` - if isinstance(impl, str): - if impl.lower() not in _ALL_IMPLS: - raise ValueError('`impl` {!r} not understood'.format(impl)) - + if is_string(impl): if impl.lower() not in RAY_TRAFO_IMPLS.keys(): raise ValueError( - '{!r} back-end not available'.format(impl) + 'The {!r} `impl` is not found. This `impl` is either ' + 'not supported, it may be misspelled, or external ' + 'packages required are not available. Consult ' + '`RAY_TRAFO_IMPLS` to find the run-time available ' + 'implementations.'.format(impl) ) impl_type = RAY_TRAFO_IMPLS[impl.lower()] @@ -260,7 +259,7 @@ def _check_impl(impl): impl_instance = impl else: raise TypeError( - '`impl` {!r} should be a `str`, or an object or type ' + '`impl` {!r} should be a string, or an object or type ' 'having a `call_forward()` and/or `call_backward()`. ' ''.format(type(impl)) ) From 6f152463474b82d199eaa31a58cee6f3b0fb5ff2 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 16 Apr 2020 15:11:37 +0200 Subject: [PATCH 36/39] Import sorting order corrected --- odl/tomo/backends/skimage_radon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/odl/tomo/backends/skimage_radon.py b/odl/tomo/backends/skimage_radon.py index fa4ef8d38c7..6ab9ba11d97 100644 --- a/odl/tomo/backends/skimage_radon.py +++ b/odl/tomo/backends/skimage_radon.py @@ -17,8 +17,8 @@ from odl.discr import ( DiscretizedSpace, uniform_discr_frompartition, uniform_partition) from odl.discr.discr_utils import linear_interpolator, point_collocation -from odl.tomo.geometry import Geometry, Parallel2dGeometry from odl.tomo.backends.util import _add_default_complex_impl +from odl.tomo.geometry import Geometry, Parallel2dGeometry from odl.util.utility import writable_array try: From 034494c1ef2a5cf44de57086e757ea8b5218ffa5 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 16 Apr 2020 15:53:25 +0200 Subject: [PATCH 37/39] Renamed `_check_impl` and `create_impl`. Simplified docstrings. --- odl/tomo/operators/ray_trafo.py | 41 +++++++++++++++------------------ 1 file changed, 19 insertions(+), 22 deletions(-) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 7c5a796fb0f..8f66f052dcf 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -190,7 +190,7 @@ def __init__(self, vol_space, geometry, **kwargs): # Check `impl` impl = kwargs.pop('impl', None) - impl_type, self.__cached_impl = self._check_impl(impl) + impl_type, self.__cached_impl = self._initialize_impl(impl) self._impl_type = impl_type if is_string(impl): self.__impl = impl.lower() @@ -212,7 +212,7 @@ def __init__(self, vol_space, geometry, **kwargs): ) @staticmethod - def _check_impl(impl): + def _initialize_impl(impl): """Internal method to verify the validity of the `impl` kwarg.""" impl_instance = None @@ -246,8 +246,8 @@ def _check_impl(impl): if not callable(forward) and not callable(backward): raise TypeError( - 'type {!r} must be have a `call_forward` ' - 'or `call_backward`.'.format(impl) + 'Type {!r} must have a `call_forward()` ' + 'and/or `call_backward()`.'.format(impl) ) if isinstance(impl, type): @@ -275,7 +275,7 @@ def impl(self): return self.__impl - def create_impl(self, use_cache=True): + def get_impl(self, use_cache=True): """Fetches or instantiates implementation backend for evaluation. Parameters @@ -298,26 +298,25 @@ def create_impl(self, use_cache=True): return self.__cached_impl def _call(self, x, out=None, **kwargs): - """Forward projection + """Forward projection. Parameters ---------- x : DiscreteLpElement - A volume. Must be an element from `RayTransform.domain`. + A volume. Must be an element of `RayTransform.domain`. out : `RayTransform.range` element, optional - A DiscreteLpElement from `RayTransform.range`, to which the result - of the operator evaluation is written. + Element to which the result of the operator evaluation is written. **kwargs - Arbitrary keyword arguments, passed on to the implementation + Extra keyword arguments, passed on to the implementation backend. Returns ------- DiscreteLpElement - A sinogram, from the projection space `RayTransform.range`. + Result of the transform, an element of the range. """ - return self.create_impl(self.use_cache).call_forward(x, out, **kwargs) + return self.get_impl(self.use_cache).call_forward(x, out, **kwargs) @property def geometry(self): @@ -343,30 +342,28 @@ class RayBackProjection(Operator): """Adjoint of the discrete Ray transform between L^p spaces.""" def _call(self, x, out=None, **kwargs): - """Backprojection + """Backprojection. Parameters ---------- x : DiscreteLpElement - A sinogram. Must be an element from + A sinogram. Must be an element of `RayTransform.range` (domain of `RayBackProjection`). out : `RayBackProjection.domain` element, optional - A volume (an element from the volume space - `RayTransform.domain` or, equivalently, - `RayBackprojection.range`). The result of this - evaluation is written into this variable. + A volume to which the result of this evaluation is + written. **kwargs - Arbitrary keyword arguments, passed on to the + Extra keyword arguments, passed on to the implementation backend. Returns ------- DiscreteLpElement - The reconstructed volume, an element in - `RayBackProjection.range`. + Result of the transform in the domain + of `RayProjection`. """ - return ray_trafo.create_impl( + return ray_trafo.get_impl( ray_trafo.use_cache ).call_backward(x, out, **kwargs) From 55b0c63841e55f641ac7f221e2f1728f294fb7b2 Mon Sep 17 00:00:00 2001 From: Adriaan Graas Date: Thu, 16 Apr 2020 18:11:19 +0200 Subject: [PATCH 38/39] Whitespace changes --- odl/test/tomo/operators/ray_trafo_test.py | 3 +-- odl/tomo/operators/ray_trafo.py | 3 --- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/odl/test/tomo/operators/ray_trafo_test.py b/odl/test/tomo/operators/ray_trafo_test.py index 5dad840dc38..a5dc3c7eccf 100644 --- a/odl/test/tomo/operators/ray_trafo_test.py +++ b/odl/test/tomo/operators/ray_trafo_test.py @@ -20,7 +20,6 @@ skip_if_no_astra, skip_if_no_astra_cuda, skip_if_no_skimage) from odl.util.testutils import all_almost_equal, simple_fixture - # --- pytest fixtures --- # @@ -106,7 +105,6 @@ def geometry(request): 'par2d skimage half_uniform']) ) - projector_ids = [ " geom='{}' - impl='{}' - angles='{}' ".format(*p.values[0].split()) for p in projectors @@ -349,6 +347,7 @@ def test_complex(impl): assert all_almost_equal(backproj_vol.real, true_vol_re) assert all_almost_equal(backproj_vol.imag, true_vol_im) + def test_anisotropic_voxels(geometry): """Test projection and backprojection with anisotropic voxels.""" ndim = geometry.ndim diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 8f66f052dcf..ca859ba2b60 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -272,7 +272,6 @@ def impl(self): If a custom ``impl`` was provided this method returns a ``str`` of the type.""" - return self.__impl def get_impl(self, use_cache=True): @@ -315,7 +314,6 @@ def _call(self, x, out=None, **kwargs): DiscreteLpElement Result of the transform, an element of the range. """ - return self.get_impl(self.use_cache).call_forward(x, out, **kwargs) @property @@ -362,7 +360,6 @@ def _call(self, x, out=None, **kwargs): Result of the transform in the domain of `RayProjection`. """ - return ray_trafo.get_impl( ray_trafo.use_cache ).call_backward(x, out, **kwargs) From 59829045f1f5c70463ba36902774f892e7e625d6 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Thu, 16 Apr 2020 20:24:36 +0200 Subject: [PATCH 39/39] Turn MD syntax in docstring to rst --- odl/tomo/backends/util.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/odl/tomo/backends/util.py b/odl/tomo/backends/util.py index d6bec82f7a4..17782a56bdc 100644 --- a/odl/tomo/backends/util.py +++ b/odl/tomo/backends/util.py @@ -22,12 +22,11 @@ def _add_default_complex_impl(fn): The wrapper will only work for methods of which the class provides a `vol_space` and `proj_space`. Typically, this will then work as a decorator - on the method, e.g.: - ``` + on the method, e.g. :: + @_add_default_complex_impl def call_forward(self, x, out=None, **kwargs): # Code that should run for real input and output - ``` Parameters ----------