diff options
Diffstat (limited to 'src/silx/opencl/sinofilter.py')
| -rw-r--r-- | src/silx/opencl/sinofilter.py | 410 |
1 files changed, 410 insertions, 0 deletions
diff --git a/src/silx/opencl/sinofilter.py b/src/silx/opencl/sinofilter.py new file mode 100644 index 0000000..fc447de --- /dev/null +++ b/src/silx/opencl/sinofilter.py @@ -0,0 +1,410 @@ +#!/usr/bin/env python +# /*########################################################################## +# +# Copyright (c) 2016-2023 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ +"""Module for sinogram filtering on CPU/GPU.""" + +__authors__ = ["P. Paleo"] +__license__ = "MIT" +__date__ = "07/06/2019" + +import numpy as np +from math import pi + + +import pyopencl.array as parray +from .common import pyopencl as cl +from .processing import OpenclProcessing +from ..math.fft.clfft import CLFFT, __have_clfft__ +from ..math.fft.npfft import NPFFT +from ..image.tomography import generate_powers, get_next_power, compute_fourier_filter + + +class SinoFilter(OpenclProcessing): + """A class for performing sinogram filtering on GPU using OpenCL. + + This is a convolution in the Fourier space, along one dimension: + + - In 2D: (n_a, d_x): n_a filterings (1D FFT of size d_x) + - In 3D: (n_z, n_a, d_x): n_z*n_a filterings (1D FFT of size d_x) + """ + + kernel_files = ["array_utils.cl"] + powers = generate_powers() + + def __init__( + self, + sino_shape, + filter_name=None, + ctx=None, + devicetype="all", + platformid=None, + deviceid=None, + profile=False, + extra_options=None, + ): + """Constructor of OpenCL FFT-Convolve. + + :param sino_shape: shape of the sinogram. + :param filter_name: Name of the filter. Defaut is "ram-lak". + :param ctx: actual working context, left to None for automatic + initialization from device type or platformid/deviceid + :param devicetype: type of device, can be "CPU", "GPU", "ACC" or "ALL" + :param platformid: integer with the platform_identifier, as given by + clinfo + :param deviceid: Integer with the device identifier, as given by clinfo + :param profile: switch on profiling to be able to profile at the kernel + level, store profiling elements (makes code slightly + slower) + :param dict extra_options: Advanced extra options. + Current options are: cutoff, use_numpy_fft + """ + OpenclProcessing.__init__( + self, + ctx=ctx, + devicetype=devicetype, + platformid=platformid, + deviceid=deviceid, + profile=profile, + ) + + self._init_extra_options(extra_options) + self._calculate_shapes(sino_shape) + self._init_fft() + self._allocate_memory() + self._compute_filter(filter_name) + self._init_kernels() + + def _calculate_shapes(self, sino_shape): + """ + + :param sino_shape: shape of the sinogram. + """ + self.ndim = len(sino_shape) + if self.ndim == 2: + n_angles, dwidth = sino_shape + else: + raise ValueError( + "Invalid sinogram number of dimensions: " "expected 2 dimensions" + ) + self.sino_shape = sino_shape + self.n_angles = n_angles + self.dwidth = dwidth + self.dwidth_padded = get_next_power(2 * self.dwidth, powers=self.powers) + self.sino_padded_shape = (n_angles, self.dwidth_padded) + sino_f_shape = list(self.sino_padded_shape) + sino_f_shape[-1] = sino_f_shape[-1] // 2 + 1 + self.sino_f_shape = tuple(sino_f_shape) + + def _init_extra_options(self, extra_options): + """ + + :param dict extra_options: Advanced extra options. + Current options are: cutoff, + """ + self.extra_options = { + "cutoff": 1.0, + "use_numpy_fft": False, + } + if extra_options is not None: + self.extra_options.update(extra_options) + + def _init_fft(self): + if __have_clfft__ and not (self.extra_options["use_numpy_fft"]): + self.fft_backend = "opencl" + self.fft = CLFFT( + self.sino_padded_shape, + dtype=np.float32, + axes=(-1,), + ctx=self.ctx, + ) + else: + self.fft_backend = "numpy" + print( + "The gpyfft module was not found. The Fourier transforms " + "will be done on CPU. For more performances, it is advised " + "to install gpyfft." + "" + ) + self.fft = NPFFT( + template=np.zeros(self.sino_padded_shape, "f"), + axes=(-1,), + ) + + def _allocate_memory(self): + self.d_filter_f = parray.zeros( + self.queue, (self.sino_f_shape[-1],), np.complex64 + ) + self.is_cpu = self.device.type == "CPU" + # These are already allocated by FFT() if using the opencl backend + if self.fft_backend == "opencl": + self.d_sino_padded = self.fft.data_in + self.d_sino_f = self.fft.data_out + else: + # When using the numpy backend, arrays are not pre-allocated + self.d_sino_padded = np.zeros(self.sino_padded_shape, "f") + self.d_sino_f = np.zeros(self.sino_f_shape, np.complex64) + # These are needed for rectangular memcpy in certain cases (see below). + self.tmp_sino_device = parray.zeros(self.queue, self.sino_shape, "f") + self.tmp_sino_host = np.zeros(self.sino_shape, "f") + + def _compute_filter(self, filter_name): + """ + + :param str filter_name: filter name + """ + self.filter_name = filter_name or "ram-lak" + filter_f = compute_fourier_filter( + self.dwidth_padded, + self.filter_name, + cutoff=self.extra_options["cutoff"], + )[ + : self.dwidth_padded // 2 + 1 + ] # R2C + self.set_filter(filter_f, normalize=True) + + def set_filter(self, h_filt, normalize=True): + """ + Set a filter for sinogram filtering. + + :param h_filt: Filter. Each line of the sinogram will be filtered with + this filter. It has to be the Real-to-Complex Fourier Transform + of some real filter, padded to 2*sinogram_width. + :param normalize: Whether to normalize the filter with pi/num_angles. + """ + if h_filt.size != self.sino_f_shape[-1]: + raise ValueError( + """ + Invalid filter size: expected %d, got %d. + Please check that the filter is the Fourier R2C transform of + some real 1D filter. + """ + % (self.sino_f_shape[-1], h_filt.size) + ) + if not (np.iscomplexobj(h_filt)): + print("Warning: expected a complex Fourier filter") + self.filter_f = h_filt + if normalize: + self.filter_f *= pi / self.n_angles + self.filter_f = self.filter_f.astype(np.complex64) + self.d_filter_f[:] = self.filter_f[:] + + def _init_kernels(self): + OpenclProcessing.compile_kernels(self, self.kernel_files) + h, w = self.d_sino_f.shape + self.mult_kern_args = ( + self.queue, + (int(w), (int(h))), + None, + self.d_sino_f.data, + self.d_filter_f.data, + np.int32(w), + np.int32(h), + ) + + def check_array(self, arr): + if arr.dtype != np.float32: + raise ValueError("Expected data type = numpy.float32") + if arr.shape != self.sino_shape: + raise ValueError( + "Expected sinogram shape %s, got %s" % (self.sino_shape, arr.shape) + ) + if not (isinstance(arr, np.ndarray) or isinstance(arr, parray.Array)): + raise ValueError("Expected either numpy.ndarray or " "pyopencl.array.Array") + + def copy2d(self, dst, src, transfer_shape, dst_offset=(0, 0), src_offset=(0, 0)): + """ + + :param dst: + :param src: + :param transfer_shape: + :param dst_offset: + :param src_offset: + """ + shape = tuple(int(i) for i in transfer_shape[::-1]) + ev = self.kernels.cpy2d( + self.queue, + shape, + None, + dst.data, + src.data, + np.int32(dst.shape[1]), + np.int32(src.shape[1]), + np.int32(dst_offset), + np.int32(src_offset), + np.int32(transfer_shape[::-1]), + ) + ev.wait() + + def copy2d_host( + self, dst, src, transfer_shape, dst_offset=(0, 0), src_offset=(0, 0) + ): + """ + + :param dst: + :param src: + :param transfer_shape: + :param dst_offset: + :param src_offset: + """ + s = transfer_shape + do = dst_offset + so = src_offset + dst[do[0] : do[0] + s[0], do[1] : do[1] + s[1]] = src[ + so[0] : so[0] + s[0], so[1] : so[1] + s[1] + ] + + def _prepare_input_sino(self, sino): + """ + :param sino: sinogram + """ + self.check_array(sino) + self.d_sino_padded.fill(0) + if self.fft_backend == "opencl": + # OpenCL backend: FFT/mult/IFFT are done on device. + if isinstance(sino, np.ndarray): + # OpenCL backend + numpy input: copy H->D. + # As pyopencl does not support rectangular copies, we have to + # do a copy H->D in a temporary device buffer, and then call a + # kernel doing the rectangular D-D copy. + self.tmp_sino_device[:] = sino[:] + if self.is_cpu: + self.tmp_sino_device.finish() + d_sino_ref = self.tmp_sino_device + else: + d_sino_ref = sino + # Rectangular copy D->D + self.copy2d(self.d_sino_padded, d_sino_ref, self.sino_shape) + if self.is_cpu: + self.d_sino_padded.finish() # should not be required here + else: + # Numpy backend: FFT/mult/IFFT are done on host. + if not (isinstance(sino, np.ndarray)): + # Numpy backend + pyopencl input: need to copy D->H + self.tmp_sino_host[:] = sino[:] + h_sino_ref = self.tmp_sino_host + else: + h_sino_ref = sino + # Rectangular copy H->H + self.copy2d_host(self.d_sino_padded, h_sino_ref, self.sino_shape) + + def _get_output_sino(self, output): + """ + :param Union[numpy.dtype,None] output: sinogram output. + :return: sinogram + """ + if output is None: + res = np.zeros(self.sino_shape, dtype=np.float32) + else: + res = output + if self.fft_backend == "opencl": + if isinstance(res, np.ndarray): + # OpenCL backend + numpy output: copy D->H + # As pyopencl does not support rectangular copies, we first have + # to call a kernel doing rectangular copy D->D, then do a copy + # D->H. + self.copy2d( + dst=self.tmp_sino_device, + src=self.d_sino_padded, + transfer_shape=self.sino_shape, + ) + if self.is_cpu: + self.tmp_sino_device.finish() # should not be required here + res[:] = self.tmp_sino_device.get()[:] + else: + if self.is_cpu: + self.d_sino_padded.finish() + self.copy2d(res, self.d_sino_padded, self.sino_shape) + if self.is_cpu: + res.finish() # should not be required here + else: + if not (isinstance(res, np.ndarray)): + # Numpy backend + pyopencl output: rect copy H->H + copy H->D + self.copy2d_host( + dst=self.tmp_sino_host, + src=self.d_sino_padded, + transfer_shape=self.sino_shape, + ) + res[:] = self.tmp_sino_host[:] + else: + # Numpy backend + numpy output: rect copy H->H + self.copy2d_host(res, self.d_sino_padded, self.sino_shape) + return res + + def _do_fft(self): + if self.fft_backend == "opencl": + self.fft.fft(self.d_sino_padded, output=self.d_sino_f) + if self.is_cpu: + self.d_sino_f.finish() + else: + # numpy backend does not support "output=" argument, + # and rfft always return a complex128 result. + res = self.fft.fft(self.d_sino_padded).astype(np.complex64) + self.d_sino_f[:] = res[:] + + def _multiply_fourier(self): + if self.fft_backend == "opencl": + # Everything is on device. Call the multiplication kernel. + ev = self.kernels.inplace_complex_mul_2Dby1D(*self.mult_kern_args) + ev.wait() + if self.is_cpu: + self.d_sino_f.finish() # should not be required here + else: + # Everything is on host. + self.d_sino_f *= self.filter_f + + def _do_ifft(self): + if self.fft_backend == "opencl": + if self.is_cpu: + self.d_sino_padded.fill(0) + self.d_sino_padded.finish() + self.fft.ifft(self.d_sino_f, output=self.d_sino_padded) + if self.is_cpu: + self.d_sino_padded.finish() + else: + # numpy backend does not support "output=" argument, + # and irfft always return a float64 result. + res = self.fft.ifft(self.d_sino_f).astype("f") + self.d_sino_padded[:] = res[:] + + def filter_sino(self, sino, output=None): + """ + + :param sino: sinogram + :param output: + :return: filtered sinogram + """ + # Handle input sinogram + self._prepare_input_sino(sino) + # FFT + self._do_fft() + # multiply with filter in the Fourier domain + self._multiply_fourier() + # iFFT + self._do_ifft() + # return + res = self._get_output_sino(output) + return res + # ~ return output + + __call__ = filter_sino |