Merge branch 'flatfield' into 'master'

Flatfield

See merge request paleo/nabu!4

.gitignore

 .vscode/
 *.geany
-doc/_build/
 
 # Created by https://www.gitignore.io/api/c,c++,python
 

nabu/misc/unsharp.py

@@ -16,18 +16,46 @@ class UnsharpMask(object):
     A helper class for unsharp masking.
     """
 
-    def __init__(self, shape, sigma, coeff, mode="reflect"):
+    avail_methods = ["gaussian", "log"]
+
+    def __init__(self, shape, sigma, coeff, mode="reflect", method="gaussian"):
         """
         Initialize a Unsharp mask.
-        UnsharpedImage =  (1 + coeff)*Image - coeff * ConvolutedImage
+        `UnsharpedImage =  (1 + coeff)*Image - coeff * ConvolutedImage`
+
+        If method == "log":
+        `UnsharpedImage = Image + coeff*ConvolutedImage`
+
+        Parameters
+        -----------
+        shape: tuple
+            Shape of the image.
+        sigma: float
+            Standard deviation of the Gaussian kernel
+        coeff: float
+            Coefficient in the linear combination of unsharp mask
+        mode: str, optional
+            Convolution mode. Default is "reflect"
+        method: str, optional
+            Method of unsharp mask. Can be "gaussian" (default) or "log" for
+            Laplacian of Gaussian.
         """
         self.shape = shape
         self.ndim = len(self.shape)
         self.sigma = sigma
         self.coeff = coeff
+        self._set_method(method)
         self._set_padding_mode(mode)
         self._compute_gaussian_kernel()
 
+    def _set_method(self, method):
+        if method not in self.avail_methods:
+            raise ValueError(
+                "Unknown unsharp method '%s'. Available are %s"
+                % (method, str(self.avail_methods))
+            )
+        self.method = method
+
     def _set_padding_mode(self, mode):
         self.mode = mode
         if mode not in PaddingMode.values():
@@ -48,8 +76,9 @@ class UnsharpMask(object):
         if not(__have_scipy__):
             raise ValueError("Need scipy to use self.unsharp()")
         image_b = self._blur2d(image)
-        return (1 + self.coeff) * image - self.coeff * image_b
-
-
-
+        if self.method == "gaussian":
+            res = (1 + self.coeff) * image - self.coeff * image_b
+        else: # LoG
+            res = image + self.coeff * image_b
+        return res
 

nabu/misc/unsharp_cuda.py

@@ -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

nabu/misc/unsharp_opencl.py

@@ -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

nabu/preproc/ccd_cuda.py

 import numpy as np
 from typing import Union, Any
-from ..preproc.ccd import CCDCorrection, FlatField
+from ..preproc.ccd import CCDProcessing, CCDCorrection, FlatField
 from ..cuda.kernel import CudaKernel
 from ..cuda.medfilt import MedianFilter
-from ..utils import get_cuda_srcfile
+from ..utils import get_cuda_srcfile, updiv
 
 try:
     import pycuda.gpuarray as garray
@@ -200,11 +200,15 @@ class CudaLog(CCDProcessing):
     def _init_kernels(self):
         self._do_clip_min = int(self.clip_min is not None)
         self._do_clip_max = int(self.clip_max is not None)
+        self.clip_min = np.float32(self.clip_min or 0)
+        self.clip_max = np.float32(self.clip_max or 1)
         self._nlog_srcfile = get_cuda_srcfile("ElementOp.cu")
-        nz, ny, nx = radios.shape
+        nz, ny, nx = self.radios.shape
         self._nx = np.int32(nx)
         self._ny = np.int32(ny)
         self._nz = np.int32(nz)
+        self._nthreadsperblock = (16, 16, 4) # TODO tune ?
+        self._nblocks = tuple([updiv(n, p) for n, p in zip([nx, ny, nz], self._nthreadsperblock)])
 
         self.nlog_kernel = CudaKernel(
             "nlog",
@@ -244,7 +248,9 @@ class CudaLog(CCDProcessing):
             self._ny,
             self._nz,
             self.clip_min,
-            self.clip_max
+            self.clip_max,
+            grid=self._nblocks,
+            block=self._nthreadsperblock
         )
         return output
 

nabu/reconstruction/fbp.py

@@ -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)

nabu/reconstruction/filtering.py

@@ -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)