diff options
Diffstat (limited to 'silx/opencl/sinofilter.py')
| -rw-r--r-- | silx/opencl/sinofilter.py | 205 |
1 files changed, 41 insertions, 164 deletions
diff --git a/silx/opencl/sinofilter.py b/silx/opencl/sinofilter.py index 3e4d69c..ad17acb 100644 --- a/silx/opencl/sinofilter.py +++ b/silx/opencl/sinofilter.py @@ -29,135 +29,21 @@ from __future__ import absolute_import, print_function, with_statement, division __authors__ = ["P. Paleo"] __license__ = "MIT" -__date__ = "01/02/2019" +__date__ = "07/06/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 .common import pyopencl as cl from .processing import OpenclProcessing from ..math.fft import FFT from ..math.fft.clfft import __have_clfft__ +from ..image.tomography import generate_powers, get_next_power, compute_fourier_filter 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): """ @@ -213,10 +99,10 @@ class SinoFilter(OpenclProcessing): self.sino_shape = sino_shape self.n_angles = n_angles self.dwidth = dwidth - self.dwidth_padded = self.get_next_power(2*self.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 + 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): @@ -268,16 +154,6 @@ class SinoFilter(OpenclProcessing): 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): """ @@ -286,9 +162,9 @@ class SinoFilter(OpenclProcessing): self.filter_name = filter_name or "ram-lak" filter_f = compute_fourier_filter( self.dwidth_padded, - filter_name, + self.filter_name, cutoff=self.extra_options["cutoff"], - )[:self.dwidth_padded//2+1] # R2C + )[:self.dwidth_padded // 2 + 1] # R2C self.set_filter(filter_f, normalize=True) def set_filter(self, h_filt, normalize=True): @@ -313,22 +189,18 @@ class SinoFilter(OpenclProcessing): print("Warning: expected a complex Fourier filter") self.filter_f = h_filt if normalize: - self.filter_f *= pi/self.n_angles + 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) - ) + 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: @@ -350,18 +222,16 @@ class SinoFilter(OpenclProcessing): :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]) - ) + 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)): @@ -376,11 +246,10 @@ class SinoFilter(OpenclProcessing): 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]] + 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) @@ -401,7 +270,7 @@ class SinoFilter(OpenclProcessing): # 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() + 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)): @@ -415,7 +284,6 @@ class SinoFilter(OpenclProcessing): def _get_output_sino(self, output): """ - :param Union[numpy.dtype,None] output: sinogram output. :return: sinogram """ @@ -433,14 +301,14 @@ class SinoFilter(OpenclProcessing): src=self.d_sino_padded, transfer_shape=self.sino_shape) if self.is_cpu: - self.tmp_sino_device.finish() - res[:] = self.tmp_sino_device[:] + 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() + 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 @@ -467,11 +335,12 @@ class SinoFilter(OpenclProcessing): def _multiply_fourier(self): if self.fft_backend == "opencl": # Everything is on device. Call the multiplication kernel. - self.kernels.inplace_complex_mul_2Dby1D( + ev = self.kernels.inplace_complex_mul_2Dby1D( *self.mult_kern_args ) + ev.wait() if self.is_cpu: - self.d_sino_f.finish() + self.d_sino_f.finish() # should not be required here else: # Everything is on host. self.d_sino_f *= self.filter_f @@ -508,7 +377,7 @@ class SinoFilter(OpenclProcessing): # return res = self._get_output_sino(output) return res - #~ return output + # ~ return output __call__ = filter_sino @@ -519,6 +388,14 @@ class SinoFilter(OpenclProcessing): # - Compatibility - # ------------------- + +def nextpow2(N): + p = 1 + while p < N: + p *= 2 + return p + + @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. |