#!/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 .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
[docs]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
[docs]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
[docs]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))
[docs]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
[docs]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"]
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,
"""
OpenclProcessing.__init__(self, ctx=ctx, devicetype=devicetype,
platformid=platformid, deviceid=deviceid,
profile=profile)
self._calculate_shapes(sino_shape)
self._init_fft()
self._allocate_memory()
self._compute_filter(filter_name, extra_options)
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 = 2*self.dwidth # TODO nextpow2 ?
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.,
}
if extra_options is not None:
self.extra_options.update(extra_options)
def _init_fft(self):
if __have_clfft__:
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 _compute_filter(self, filter_name, extra_options):
"""
:param str filter_name: filter name
:param dict extra_options: Advanced extra options.
"""
self._init_extra_options(extra_options)
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)
[docs] 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")
[docs] 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])
)
[docs] 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[:]
[docs] 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")
[docs]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)