diff options
Diffstat (limited to 'silx/opencl/sinofilter.py')
-rw-r--r-- | silx/opencl/sinofilter.py | 559 |
1 files changed, 559 insertions, 0 deletions
diff --git a/silx/opencl/sinofilter.py b/silx/opencl/sinofilter.py new file mode 100644 index 0000000..3e4d69c --- /dev/null +++ b/silx/opencl/sinofilter.py @@ -0,0 +1,559 @@ +#!/usr/bin/env python +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2016 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.""" + +from __future__ import absolute_import, print_function, with_statement, division + +__authors__ = ["P. Paleo"] +__license__ = "MIT" +__date__ = "01/02/2019" + +import numpy as np +from math import pi +from itertools import product +from bisect import bisect + +from .common import pyopencl as cl +import pyopencl.array as parray +from .processing import OpenclProcessing +from ..math.fft import FFT +from ..math.fft.clfft import __have_clfft__ +from ..utils.deprecation import deprecated + + +def nextpow2(N): + p = 1 + while p < N: + p *= 2 + return p + +def compute_ramlak_filter(dwidth_padded, dtype=np.float32): + """ + Compute the Ramachandran-Lakshminarayanan (Ram-Lak) filter, used in + filtered backprojection. + + :param dwidth_padded: width of the 2D sinogram after padding + :param dtype: data type + """ + #~ dwidth_padded = dwidth * 2 + L = dwidth_padded + h = np.zeros(L, dtype=dtype) + L2 = L//2+1 + h[0] = 1/4. + j = np.linspace(1, L2, L2//2, False).astype(dtype) # np < 1.9.0 + # h[2::2] = 0 + h[1:L2:2] = -1./(pi**2 * j**2) + # h[-1:L2-1:-2] = -1./(pi**2 * j**2) + h[L2:] = np.copy(h[1:L2-1][::-1]) + return h + + +def tukey(N, alpha=0.5): + """ + Compute the Tukey apodization window. + + :param int N: Number of points. + :param float alpha: + """ + apod = np.zeros(N) + x = np.arange(N)/(N-1) + r = alpha + M1 = (0 <= x) * (x < r/2) + M2 = (r/2 <= x) * (x <= 1 - r/2) + M3 = (1 - r/2 < x) * (x <= 1) + apod[M1] = (1 + np.cos(2*pi/r * (x[M1] - r/2)))/2. + apod[M2] = 1. + apod[M3] = (1 + np.cos(2*pi/r * (x[M3] - 1 + r/2)))/2. + return apod + + +def lanczos(N): + """ + Compute the Lanczos window (truncated sinc) of width N. + + :param int N: window width + """ + x = np.arange(N)/(N-1) + return np.sin(pi*(2*x-1))/(pi*(2*x-1)) + + +def compute_fourier_filter(dwidth_padded, filter_name, cutoff=1.): + """ + Compute the filter used for FBP. + + :param dwidth_padded: padded detector width. As the filtering is done by the + Fourier convolution theorem, dwidth_padded should be at least 2*dwidth. + :param filter_name: Name of the filter. Available filters are: + Ram-Lak, Shepp-Logan, Cosine, Hamming, Hann, Tukey, Lanczos. + :param cutoff: Cut-off frequency, if relevant. + """ + Nf = dwidth_padded + #~ filt_f = np.abs(np.fft.fftfreq(Nf)) + rl = compute_ramlak_filter(Nf, dtype=np.float64) + filt_f = np.fft.fft(rl) + + filter_name = filter_name.lower() + if filter_name in ["ram-lak", "ramlak"]: + return filt_f + + w = 2 * pi * np.fft.fftfreq(dwidth_padded) + d = cutoff + apodization = { + # ~OK + "shepp-logan": np.sin(w[1:Nf]/(2*d))/(w[1:Nf]/(2*d)), + # ~OK + "cosine": np.cos(w[1:Nf]/(2*d)), + # OK + "hamming": 0.54*np.ones_like(filt_f)[1:Nf] + .46 * np.cos(w[1:Nf]/d), + # OK + "hann": (np.ones_like(filt_f)[1:Nf] + np.cos(w[1:Nf]/d))/2., + # These one is not compatible with Astra - TODO investigate why + "tukey": np.fft.fftshift(tukey(dwidth_padded, alpha=d/2.))[1:Nf], + "lanczos": np.fft.fftshift(lanczos(dwidth_padded))[1:Nf], + } + if filter_name not in apodization: + raise ValueError("Unknown filter %s. Available filters are %s" % + (filter_name, str(apodization.keys()))) + filt_f[1:Nf] *= apodization[filter_name] + return filt_f + + +def generate_powers(): + """ + Generate a list of powers of [2, 3, 5, 7], + up to (2**15)*(3**9)*(5**6)*(7**5). + """ + primes = [2, 3, 5, 7] + maxpow = {2: 15, 3: 9, 5: 6, 7: 5} + valuations = [] + for prime in primes: + valuations.append(range(maxpow[prime]+1)) + powers = product(*valuations) + res = [] + for pw in powers: + res.append(np.prod(list(map(lambda x : x[0]**x[1], zip(primes, pw))))) + return np.unique(res) + + + +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 = self.get_next_power(2*self.dwidth) + 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., + "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 = FFT( + self.sino_padded_shape, + dtype=np.float32, + axes=(-1,), + backend="opencl", + 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 = FFT( + template=np.zeros(self.sino_padded_shape, "f"), + axes=(-1,), + backend="numpy", + ) + + 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 get_next_power(self, n): + """ + Given a number, get the closest (upper) number p such that + p is a power of 2, 3, 5 and 7. + """ + idx = bisect(self.powers, n) + if self.powers[idx-1] == n: + return n + return self.powers[idx] + + 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, + 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, + np.int32(self.d_sino_f.shape[::-1]), + 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: + """ + self.kernels.cpy2d( + self.queue, + np.int32(transfer_shape[::-1]), + 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]) + ) + + 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() + 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() + res[:] = self.tmp_sino_device[:] + 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() + 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. + self.kernels.inplace_complex_mul_2Dby1D( + *self.mult_kern_args + ) + if self.is_cpu: + self.d_sino_f.finish() + 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 + + + + +# ------------------- +# - Compatibility - +# ------------------- + +@deprecated(replacement="Backprojection.sino_filter", since_version="0.10") +def fourier_filter(sino, filter_=None, fft_size=None): + """Simple np based implementation of fourier space filter. + This function is deprecated, please use silx.opencl.sinofilter.SinoFilter. + + :param sino: of shape shape = (num_projs, num_bins) + :param filter: filter function to apply in fourier space + :fft_size: size on which perform the fft. May be larger than the sino array + :return: filtered sinogram + """ + assert sino.ndim == 2 + num_projs, num_bins = sino.shape + if fft_size is None: + fft_size = nextpow2(num_bins * 2 - 1) + else: + assert fft_size >= num_bins + if fft_size == num_bins: + sino_zeropadded = sino.astype(np.float32) + else: + sino_zeropadded = np.zeros((num_projs, fft_size), + dtype=np.complex64) + sino_zeropadded[:, :num_bins] = sino.astype(np.float32) + + if filter_ is None: + h = np.zeros(fft_size, dtype=np.float32) + L2 = fft_size // 2 + 1 + h[0] = 1 / 4. + j = np.linspace(1, L2, L2 // 2, False) + h[1:L2:2] = -1. / (np.pi ** 2 * j ** 2) + h[L2:] = np.copy(h[1:L2 - 1][::-1]) + filter_ = np.fft.fft(h).astype(np.complex64) + + # Linear convolution + sino_f = np.fft.fft(sino, fft_size) + sino_f = sino_f * filter_ + sino_filtered = np.fft.ifft(sino_f)[:, :num_bins].real + + return np.ascontiguousarray(sino_filtered.real, dtype=np.float32) |