Commit 55bbe459 authored by Pierre Paleo's avatar Pierre Paleo
Browse files

Merge branch 'dff_cuda' into 'master'

Double flatfield: cuda backend

Closes #117

See merge request !42
parents 70d4b118 6a1b3af6
Pipeline #26985 passed with stages
in 2 minutes and 52 seconds
nabu.preproc.double\_flat\_field module
=======================================
nabu.preproc.double\_flatfield module
=====================================
.. automodule:: nabu.preproc.double_flat_field
.. automodule:: nabu.preproc.double_flatfield
:members:
:undoc-members:
:show-inheritance:
nabu.preproc.double\_flatfield\_cuda module
===========================================
.. automodule:: nabu.preproc.double_flatfield_cuda
:members:
:undoc-members:
:show-inheritance:
......@@ -10,7 +10,8 @@ Submodules
nabu.preproc.alignment
nabu.preproc.ccd
nabu.preproc.ccd_cuda
nabu.preproc.double_flat_field
nabu.preproc.double_flatfield
nabu.preproc.double_flatfield_cuda
nabu.preproc.phase
nabu.preproc.phase_cuda
nabu.preproc.rings
......
......@@ -29,19 +29,19 @@ API: [CCDCorrection](apidoc/nabu.preproc.ccd.rst#nabu.preproc.ccd.CCDCorrection)
### Double flat-field
This is a projections-based rings artefacts removal.
This is a projections-based rings artefacts removal.
Some defects might remain in the projections even after flat-fielding. These defects are visible as structured noise in the projections or sinograms.
Some defects might remain in the projections even after flat-fielding. These defects are visible as structured noise in the projections or sinograms.
By doing the average of all projections, the genuine features are canceled out and only the defects remain. This "average image" is used to remove the defects from the radios. Schematically, this method does `radio = radio / mean(radios)` (or other operations if the logarithm of radios was already taken).
Be aware that by doing so, you might lose the quantitativeness of the reconstructed data.
This method assumes that when averaging all the radios, the genuine feature will cancel and only spurious artefacts will remain. This assumption can fail if genuine features are (more or less) independent from the projection angle, ex. ring-shaped.
This method assumes that when averaging all the radios, the genuine feature will cancel and only spurious artefacts will remain. This assumption can fail if genuine features are (more or less) independent from the projection angle, ex. ring-shaped.
Configuration file key: section `[preproc]`: `double_flatfield_enabled = 1`.
API: [DoubleFlatField](apidoc/nabu.preproc.double_flat_field.rst#nabu.preproc.double_flat_field.DoubleFlatField)
API: [DoubleFlatField](apidoc/nabu.preproc.double_flatfield.rst#nabu.preproc.double_flatfield.DoubleFlatField) and [CudaDoubleFlatField](apidoc/nabu.preproc.double_flatfield_cuda.rst#nabu.preproc.double_flatfield_cuda.CudaDoubleFlatField)
### Logarithm
......@@ -67,7 +67,7 @@ See also: [Phase retrieval](phase.md)
### Paganin phase retrieval
The Paganin method consists in applying a band-pass filter on the radios.
The Paganin method consists in applying a band-pass filter on the radios.
It depends on the δ/β ratio (assumed to be constant in all the image) and the incoming beam energy.
Configuration file: section `[phase]`: `method = Paganin`, `delta_beta = 1000.0`
......@@ -88,7 +88,7 @@ Setting `coeff` to zero (default) disables unsharp masking.
Tomographic reconstruction is the process of reconstructing the volume from projections/sinograms.
Configuration file: section `[reconstruction]`.
Configuration file: section `[reconstruction]`.
Related keys: `angles_file`, `angle_offset`, `rotation_axis_position`, `enable_halftomo`
`start_x`, `end_x`, `start_y`, `end_y`, `start_z`, `end_z`.
......
......@@ -4,8 +4,8 @@ from .logger import LoggerOrPrint
from ..utils import check_supported
from ..io.reader import ChunkReader
from ..preproc.ccd import FlatField, Log
from ..preproc.double_flatfield import DoubleFlatField
from ..preproc.phase import PaganinPhaseRetrieval
from ..preproc.double_flat_field import DoubleFlatField
from ..preproc.sinogram import SinoProcessing
from ..misc.unsharp import UnsharpMask
from ..resources.utils import is_hdf5_extension
......@@ -20,11 +20,11 @@ class FullFieldPipeline:
"""
FlatFieldClass = FlatField
DoubleFlatFieldClass = DoubleFlatField
PaganinPhaseRetrievalClass = PaganinPhaseRetrieval
UnsharpMaskClass = UnsharpMask
SinoProcessingClass = SinoProcessing
MLogClass = Log
DoubleFlatFieldClass = DoubleFlatField
FBPClass = None # For now we don't have a plain python/numpy backend for reconstruction
def __init__(self, process_config, sub_region, logger=None, extra_options=None):
......
from math import ceil
import numpy as np
from ..preproc.ccd_cuda import CudaFlatField, CudaLog
from ..preproc.double_flatfield_cuda import CudaDoubleFlatField
from ..preproc.phase_cuda import CudaPaganinPhaseRetrieval
from ..preproc.sinogram_cuda import CudaSinoProcessing
from ..preproc.sinogram import SinoProcessing
......@@ -19,6 +20,7 @@ class CudaFullFieldPipeline(FullFieldPipeline):
"""
FlatFieldClass = CudaFlatField
DoubleFlatFieldClass = CudaDoubleFlatField
PaganinPhaseRetrievalClass = CudaPaganinPhaseRetrieval
UnsharpMask = CudaUnsharpMask
SinoProcessingClass = CudaSinoProcessing
......
import numpy as np
from os.path import dirname
from silx.image.utils import gaussian_kernel
import pycuda.gpuarray as garray
from pycuda.compiler import SourceModule
from ..utils import updiv, get_cuda_srcfile
from ..cuda.utils import __has_pycuda__
from ..misc.utils import ConvolutionInfos
from .processing import CudaProcessing
if __has_pycuda__:
import pycuda.gpuarray as garray
from pycuda.compiler import SourceModule
class Convolution(CudaProcessing):
"""
......
import numpy as np
from .utils import get_cuda_context
import pycuda.driver as cuda
import pycuda.gpuarray as garray
dev_attrs = cuda.device_attribute
from .utils import get_cuda_context, __has_pycuda__
if __has_pycuda__:
import pycuda.driver as cuda
import pycuda.gpuarray as garray
dev_attrs = cuda.device_attribute
# NB: we must detach from a context before creating another context
class CudaProcessing(object):
......
import pycuda.gpuarray as garray
from ..cuda.utils import __has_pycuda__
from ..cuda.convolution import Convolution
from pycuda.elementwise import ElementwiseKernel
from ..cuda.processing import CudaProcessing
from .unsharp import UnsharpMask
if __has_pycuda__:
import pycuda.gpuarray as garray
from pycuda.elementwise import ElementwiseKernel
class CudaUnsharpMask(UnsharpMask, CudaProcessing):
def __init__(self, shape, sigma, coeff, mode="reflect", method="gaussian",
......
......@@ -115,6 +115,8 @@ class DoubleFlatField(CCDProcessing):
self.output_is_mlog = output_is_mlog
self.average_is_on_log = average_is_on_log
self.sigma_filter = sigma_filter
if self.sigma_filter is not None and abs(float(self.sigma_filter)) < 1e-4:
self.sigma_filter = None
self.filter_mode = filter_mode
proc = lambda x,o: np.copyto(o, x)
if self.input_is_mlog:
......@@ -158,7 +160,7 @@ class DoubleFlatField(CCDProcessing):
acc /= radios.shape[0]
if self.sigma_filter is not None and abs(float(self.sigma_filter)) > 1e-4:
if self.sigma_filter is not None:
acc = acc - gaussian_filter(acc, self.sigma_filter, mode=self.filter_mode)
self.doubleflatfield = self.postproc(acc)
# Handle small values to avoid issues when dividing
......
import numpy as np
from .double_flatfield import DoubleFlatField
from ..utils import check_shape
from ..cuda.utils import get_cuda_context, __has_pycuda__
from ..cuda.processing import CudaProcessing
from ..misc.unsharp_cuda import CudaUnsharpMask
if __has_pycuda__:
import pycuda.gpuarray as garray
import pycuda.cumath as cumath
class CudaDoubleFlatField(DoubleFlatField, CudaProcessing):
def __init__(
self,
shape,
result_url=None,
sub_region=None,
input_is_mlog=True,
output_is_mlog=False,
average_is_on_log=False,
sigma_filter=None,
filter_mode="reflect",
cuda_options={},
):
"""
Init double flat field with Cuda backend.
"""
DoubleFlatField.__init__(
self,
shape,
result_url=result_url,
sub_region=sub_region,
input_is_mlog=input_is_mlog,
output_is_mlog=output_is_mlog,
average_is_on_log=average_is_on_log,
sigma_filter=sigma_filter,
filter_mode=filter_mode
)
CudaProcessing.__init__(self, **cuda_options)
self._init_gaussian_filter()
def _init_gaussian_filter(self):
if self.sigma_filter is None:
return
self._unsharp_mask = CudaUnsharpMask(
self.shape, self.sigma_filter, -1.,
mode=self.filter_mode, method="log"
)
@staticmethod
def _proc_copy(x, o):
o[:] = x[:]
return o
@staticmethod
def _proc_expm(x, o):
o[:] = x[:]
o[:] *= -1
cumath.exp(o, out=o)
return o
@staticmethod
def _proc_mlog(x, o):
cumath.log(x, out=o)
o *= -1
return o
def _init_processing(self, input_is_mlog, output_is_mlog, average_is_on_log, sigma_filter, filter_mode):
self.input_is_mlog = input_is_mlog
self.output_is_mlog = output_is_mlog
self.average_is_on_log = average_is_on_log
self.sigma_filter = sigma_filter
if self.sigma_filter is not None and abs(float(self.sigma_filter)) < 1e-4:
self.sigma_filter = None
self.filter_mode = filter_mode
# proc = lambda x,o: np.copyto(o, x)
proc = self._proc_copy
if self.input_is_mlog:
if not self.average_is_on_log:
# proc = lambda x,o: np.exp(-x, out=o)
proc = self._proc_expm
else:
if self.average_is_on_log:
# proc = lambda x,o: -np.log(x, out=o)
proc = self._proc_mlog
# postproc = lambda x: x
postproc = self._proc_copy
if self.output_is_mlog:
if not self.average_is_on_log:
# postproc = lambda x: -np.log(x)
postproc = self._proc_mlog
else:
if self.average_is_on_log:
# postproc = lambda x: np.exp(-x)
postproc = self._proc_expm
self.proc = proc
self.postproc = postproc
def compute_double_flatfield(self, radios, recompute=False):
"""
Read the radios and generate the "double flat field" by averaging
and possibly other processing.
Parameters
----------
radios: array
Input radios chunk.
recompute: bool, optional
Whether to recompute the double flatfield if already computed.
"""
if not(isinstance(radios, garray.GPUArray)):
raise ValueError("Expected pycuda.gpuarray.GPUArray for radios")
if self._computed:
return self.doubleflatfield
acc = garray.zeros(radios[0].shape, "f")
tmpdat = garray.zeros(radios[0].shape, "f")
for i in range(radios.shape[0]):
self.proc(radios[i], tmpdat)
acc += tmpdat
acc /= radios.shape[0]
if self.sigma_filter is not None:
# acc = acc - gaussian_filter(acc, self.sigma_filter, mode=self.filter_mode)
self._unsharp_mask.unsharp(acc, tmpdat)
acc[:] = tmpdat[:]
self.postproc(acc, tmpdat)
self.doubleflatfield = tmpdat
# Handle small values to avoid issues when dividing
# self.doubleflatfield[np.abs(self.doubleflatfield) < self._small] = 1.
cumath.fabs(self.doubleflatfield, out=acc)
acc -= self._small # acc = abs(doubleflatfield) - _small
garray.if_positive(
acc,
self.doubleflatfield,
garray.zeros_like(acc) + self._small,
out=self.doubleflatfield
)
if self.writer is not None:
self.writer.write(self.doubleflatfield.get(), "double_flatfield")
self._computed = True
return self.doubleflatfield
def get_double_flatfield(self, radios=None, compute=False):
"""
Get the double flat field or a subregion of it.
Parameters
----------
radios: array, optional
Input radios chunk
compute: bool, optional
Whether to compute the double flatfield anyway even if a dump file
exists.
"""
if self.reader is None:
if radios is None:
raise ValueError(
"result_url was not provided. Please provide 'radios' to this function"
)
return self.compute_double_flatfield(radios)
if radios is not None and compute:
res = self.compute_double_flatfield(radios)
else:
res = self._load_dff_dump()
res = garray.to_gpu(res)
self._computed = True
return res
def apply_double_flatfield(self, radios):
"""
Apply the "double flatfield" filter on a chunk of radios.
The processing is done in-place !
"""
check_shape(radios.shape, self.radios_shape, "radios")
dff = self.get_double_flatfield(radios=radios)
for i in range(self.n_angles):
radios[i] /= dff
return radios
......@@ -6,8 +6,11 @@ import pytest
from silx.io.url import DataUrl
from tomoscan.esrf.mock import MockHDF5
from nabu.io.reader import HDF5Reader
from nabu.preproc.double_flat_field import DoubleFlatField, __have_scipy__
from nabu.preproc.double_flatfield import DoubleFlatField, __have_scipy__
from nabu.preproc.double_flatfield_cuda import CudaDoubleFlatField, __has_pycuda__
if __has_pycuda__:
import pycuda.gpuarray as garray
@pytest.fixture(scope="class")
def bootstrap(request):
......@@ -35,7 +38,10 @@ def bootstrap(request):
file_path=path.join(cls.dname, "dff.h5"),
data_path="/entry/double_flatfield/results/data"
)
cls.ff_cuda_dump_url = DataUrl(
file_path=path.join(cls.dname, "dff_cuda.h5"),
data_path="/entry/double_flatfield/results/data"
)
golden = 0
for i in range(10):
golden += exp(-i)
......@@ -64,3 +70,26 @@ class TestDoubleFlatField:
assert np.max(np.abs(mydf2 - mydf)) < self.tol
assert np.max(np.abs(mydf - self.golden)) < self.tol
@pytest.mark.skipif(not(__has_pycuda__), reason="Need pycuda for double flatfield with cuda backend")
def test_dff_cuda(self):
import pycuda.autoinit
dff = CudaDoubleFlatField(
self.radios.shape,
result_url=self.ff_cuda_dump_url
)
d_radios = garray.to_gpu(self.radios)
mydf = dff.get_double_flatfield(radios=d_radios).get()
assert path.isfile(dff.result_url.file_path())
dff2 = CudaDoubleFlatField(
self.radios.shape,
result_url=self.ff_cuda_dump_url
)
mydf2 = dff2.get_double_flatfield(radios=d_radios).get()
assert np.max(np.abs(mydf2 - mydf)) < self.tol
assert np.max(np.abs(mydf - self.golden)) < self.tol
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