Commit c139317d authored by Pierre Paleo's avatar Pierre Paleo
Browse files

Implement edge padding in cuda sinofilter

parent 1a843fb2
......@@ -9,8 +9,8 @@ from ..cuda.processing import CudaProcessing
from .unsharp import UnsharpMask
class CudaUnsharpMask(UnsharpMask, CudaProcessing):
def __init__(self, shape, sigma, coeff, mode="reflect",
device_id=None, cleanup_at_exit=True):
def __init__(self, shape, sigma, coeff, mode="reflect", method="gaussian",
ctx=None, device_id=None, cleanup_at_exit=True):
"""
NB: For now, this class is designed to use the lowest amount of GPU memory
as possible. Therefore, the input and output image/volumes are assumed
......@@ -18,8 +18,10 @@ class CudaUnsharpMask(UnsharpMask, CudaProcessing):
"""
if not(__have_cuda__):
raise ImportError("Need cuda and pycuda")
CudaProcessing.__init__(self, device_id=device_id, cleanup_at_exit=cleanup_at_exit)
UnsharpMask.__init__(self, shape, sigma, coeff, mode=mode)
CudaProcessing.__init__(
self, device_id=device_id, ctx=ctx, cleanup_at_exit=cleanup_at_exit
)
UnsharpMask.__init__(self, shape, sigma, coeff, mode=mode, method=method)
self._init_convolution()
self._init_mad_kernel()
......@@ -48,5 +50,8 @@ class CudaUnsharpMask(UnsharpMask, CudaProcessing):
assert isinstance(image, garray.GPUArray)
assert isinstance(output, garray.GPUArray)
self.convolution(image, output=output)
self.mad_kernel(output, -self.coeff, image, 1. + self.coeff)
if self.method == "gaussian":
self.mad_kernel(output, -self.coeff, image, 1. + self.coeff)
else: # log
self.mad_kernel(output, self.coeff, image, 1.)
return output
......@@ -10,7 +10,7 @@ except ImportError:
from .unsharp import UnsharpMask
class OpenclUnsharpMask(UnsharpMask, OpenclProcessing):
def __init__(self, shape, sigma, coeff, mode="reflect",
def __init__(self, shape, sigma, coeff, mode="reflect", method="gaussian",
ctx=None, devicetype="all", platformid=None, deviceid=None,
block_size=None, memory=None, profile=False):
"""
......@@ -57,5 +57,8 @@ class OpenclUnsharpMask(UnsharpMask, OpenclProcessing):
assert isinstance(image, parray.Array)
assert isinstance(output, parray.Array)
self.convolution(image, output=output)
self.mad_kernel(output, -self.coeff, image, 1. + self.coeff)
if self.method == "gaussian":
self.mad_kernel(output, -self.coeff, image, 1. + self.coeff)
else:
self.mad_kernel(output, self.coeff, image, 1.)
return output
......@@ -38,6 +38,7 @@ class Backprojector(CudaProcessing):
self._axis_array = None
self.extra_options = {
"filter_name": "ramlak",
"padding_mode": "zeros",
"axis_correction": None,
}
self.extra_options.update(extra_options)
......@@ -79,7 +80,12 @@ class Backprojector(CudaProcessing):
def compile_kernels(self):
# Configure sinogram filter
self.sino_filter = SinoFilter(self.sino_shape)
self.sino_filter = SinoFilter(
self.sino_shape,
filter_name=self.extra_options["filter_name"],
padding_mode=self.extra_options["padding_mode"],
ctx=self.ctx,
)
# Configure backprojector
fname = get_cuda_srcfile("backproj.cu")
self.gpu_projector = CudaKernel("backproj", filename=fname)
......
......@@ -7,32 +7,32 @@ from silx.math.fft import FFT
from silx.image.tomography import compute_fourier_filter, get_next_power
from ..cuda.kernel import CudaKernel
from ..cuda.processing import CudaProcessing
from ..utils import get_cuda_srcfile
from ..utils import get_cuda_srcfile, check_supported, updiv
class SinoFilter(CudaProcessing):
available_padding_modes = ["zeros", "edges"]
def __init__(
self,
sino_shape,
filter_name=None,
padding_mode="zeros",
extra_options=None,
**cuda_options
):
"""
Build a sinogram filter process.
"""
super().__init__(**cuda_options)
self._init_extra_options(extra_options)
self.calculate_shapes(sino_shape)
self.init_fft()
self.allocate_memory(sino_shape)
self.compute_filter(filter_name)
self.init_kernels()
self._calculate_shapes(sino_shape)
self._init_fft()
self._allocate_memory(sino_shape)
self._compute_filter(filter_name)
self._set_padding_mode(padding_mode)
self._init_kernels()
def _init_extra_options(self, extra_options):
......@@ -42,8 +42,12 @@ class SinoFilter(CudaProcessing):
if extra_options is not None:
self.extra_options.update(extra_options)
def _set_padding_mode(self, padding_mode):
check_supported(padding_mode, self.available_padding_modes, "padding mode")
self.padding_mode = padding_mode
def calculate_shapes(self, sino_shape):
def _calculate_shapes(self, sino_shape):
self.ndim = len(sino_shape)
if self.ndim == 2:
n_angles, dwidth = sino_shape
......@@ -64,9 +68,12 @@ class SinoFilter(CudaProcessing):
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)
#
self.pad_left = (self.dwidth_padded - self.dwidth)//2
self.pad_right = self.dwidth_padded - self.dwidth - self.pad_left
def init_fft(self):
def _init_fft(self):
self.fft = FFT(
self.sino_padded_shape,
dtype=np.float32,
......@@ -75,7 +82,7 @@ class SinoFilter(CudaProcessing):
)
def allocate_memory(self, sino_shape):
def _allocate_memory(self, sino_shape):
self.d_filter_f = garray.zeros((self.sino_f_shape[-1],), np.complex64)
self.d_sino_padded = self.fft.data_in
self.d_sino_f = self.fft.data_out
......@@ -108,8 +115,7 @@ class SinoFilter(CudaProcessing):
self.d_filter_f[:] = self.filter_f[:]
def compute_filter(self, filter_name):
def _compute_filter(self, filter_name):
self.filter_name = filter_name or "ram-lak"
filter_f = compute_fourier_filter(
self.dwidth_padded,
......@@ -119,7 +125,7 @@ class SinoFilter(CudaProcessing):
self.set_filter(filter_f, normalize=True)
def init_kernels(self):
def _init_kernels(self):
fname = get_cuda_srcfile("ElementOp.cu")
if self.ndim == 2:
kernel_name = "inplace_complex_mul_2Dby1D"
......@@ -134,15 +140,47 @@ class SinoFilter(CudaProcessing):
)
self.kern_args = (self.d_sino_f, self.d_filter_f)
self.kern_args += (self.d_sino_f.shape[::-1])
self._pad_edges_kernel = CudaKernel(
"padding_edge",
filename=get_cuda_srcfile("padding.cu"),
signature="Piiiiiiii"
)
self._pad_block = (32, 32, 1)
self._pad_grid = tuple([updiv(n, p) for n, p in zip(self.sino_padded_shape[::-1], self._pad_block)])
def check_array(self, arr):
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))
def _pad_sino(self, sino):
if self.padding_mode == "edges":
self.d_sino_padded[:, :self.dwidth] = sino[:]
self._pad_edges_kernel(
self.d_sino_padded,
self.dwidth,
self.n_angles,
self.dwidth_padded,
self.n_angles,
self.pad_left,
self.pad_right,
0,
0,
grid=self._pad_grid,
block=self._pad_block
)
else: # zeros
self.d_sino_padded.fill(0)
if self.ndim == 2:
self.d_sino_padded[:, :self.dwidth] = sino[:]
else:
self.d_sino_padded[:, :, :self.dwidth] = sino[:]
def filter_sino(self, sino, output=None, no_output=False):
"""
Perform the sinogram siltering.
......@@ -157,15 +195,9 @@ class SinoFilter(CudaProcessing):
If set to True, no copy is be done. The resulting data lies
in self.d_sino_padded.
"""
self.check_array(sino)
self._check_array(sino)
# copy2d/copy3d
self.d_sino_padded.fill(0)
if self.ndim == 2:
#~ self._copy2d_sino(sino)
self.d_sino_padded[:, :self.dwidth] = sino[:]
else:
# TODO same workaround as above ?
self.d_sino_padded[:, :, :self.dwidth] = sino[:]
self._pad_sino(sino)
# FFT
self.fft.fft(self.d_sino_padded, output=self.d_sino_f)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment