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

Merge branch 'sino_normalization_url' into 'master'

Sino normalization url

Closes #255

See merge request !161
parents 0667b851 50edc0af
Pipeline #73822 passed with stages
in 7 minutes and 40 seconds
__version__ = "2022.1.3"
__version__ = "2022.1.4"
__nabu_modules__ = [
"app",
"cuda",
......
......@@ -2,6 +2,18 @@
typedef pycuda::complex<float> complex;
// Generic operations
#define OP_ADD 0
#define OP_SUB 1
#define OP_MUL 2
#define OP_DIV 3
//
#ifndef GENERIC_OP
#define GENERIC_OP OP_ADD
#endif
// arr2D *= arr1D (line by line, i.e along fast dim)
__global__ void inplace_complex_mul_2Dby1D(complex* arr2D, complex* arr1D, int width, int height) {
int x = blockDim.x * blockIdx.x + threadIdx.x;
......@@ -17,6 +29,23 @@ __global__ void inplace_complex_mul_2Dby1D(complex* arr2D, complex* arr1D, int w
arr2D[i]._M_im = a._M_im * b._M_re + a._M_re * b._M_im;
}
__global__ void inplace_generic_op_2Dby2D(float* arr2D, float* arr2D_other, int width, int height) {
int x = blockDim.x * blockIdx.x + threadIdx.x;
int y = blockDim.y * blockIdx.y + threadIdx.y;
if ((x >= width) || (y >= height)) return;
int i = y*width + x;
#if GENERIC_OP == OP_ADD
arr2D[i] += arr2D_other[i];
#elif GENERIC_OP == OP_SUB
arr2D[i] -= arr2D_other[i];
#elif GENERIC_OP == OP_MUL
arr2D[i] *= arr2D_other[i];
#elif GENERIC_OP == OP_DIV
arr2D[i] /= arr2D_other[i];
#endif
}
// arr3D *= arr1D (along fast dim)
__global__ void inplace_complex_mul_3Dby1D(complex* arr3D, complex* arr1D, int width, int height, int depth) {
......
......@@ -598,6 +598,7 @@ class ChunkedPipeline:
self.sino_normalization = self.SinoNormalizationClass(
kind=options["method"],
radios_shape=self._get_shape("sino_normalization"),
normalization_array=options["normalization_array"]
)
@use_options("build_sino", "sino_builder")
......
......@@ -108,10 +108,16 @@ nabu_config = {
},
"sino_normalization": {
"default": "",
"help": "Sinogram normalization method. Available methods are: chebyshev, none. Default is none (no normalization)",
"help": "Sinogram normalization method. Available methods are: chebyshev, subtraction, division, none. Default is none (no normalization)",
"validator": sino_normalization_validator,
"type": "advanced",
},
"sino_normalization_file": {
"default": "",
"help": "Path to the file when sino_normalization is either 'subtraction' or 'division'. To specify the path within a HDF5 file, the syntax is /path/to/file?path=/entry/data",
"validator": no_validator,
"type": "advanced",
},
"processes_file": {
"default": "",
"help": "Path to the file where some operations should be stored for later use. By default it is 'xxx_nabu_processes.h5'",
......
......@@ -2,6 +2,7 @@ import os
import posixpath
import numpy as np
from silx.io import get_data
from silx.io.url import DataUrl
from ...utils import copy_dict_items, compare_dicts
from ...io.utils import hdf5_entry_exists, get_h5_value
from ...io.reader import import_h5_to_dict
......@@ -33,6 +34,7 @@ class ProcessConfig(ProcessConfigBase):
self._update_rotation_angles()
self._get_translation_file("reconstruction", "translation_movements_file", "translations")
self._get_translation_file("phase", "ctf_translations_file", "ctf_translations")
self._get_user_sino_normalization()
def _get_translation_file(self, config_section, config_key, dataset_info_attr):
......@@ -143,6 +145,22 @@ class ProcessConfig(ProcessConfigBase):
self.dataset_info.axis_position = cor
def _get_user_sino_normalization(self):
self._sino_normalization_arr = None
norm = self.nabu_config["preproc"]["sino_normalization"]
if norm not in ["subtraction", "division"]:
return
norm_path = "silx://" + self.nabu_config["preproc"]["sino_normalization_file"].strip()
url = DataUrl(norm_path)
try:
norm_array = get_data(url)
self._sino_normalization_arr = norm_array.astype("f")
except (ValueError, OSError) as exc:
error_msg = "Could not load sino_normalization_file %s. The error was: %s" % (norm_path, str(exc))
self.logger.error(error_msg)
raise ValueError(error_msg)
@property
def do_halftomo(self):
"""
......@@ -312,8 +330,10 @@ class ProcessConfig(ProcessConfigBase):
if nabu_config["preproc"]["sino_normalization"] is not None:
tasks.append("sino_normalization")
options["sino_normalization"] = {
"method": nabu_config["preproc"]["sino_normalization"]
"method": nabu_config["preproc"]["sino_normalization"],
"normalization_array": self._sino_normalization_arr
}
#
# Sinogram-based rings artefacts removal
#
......
......@@ -100,6 +100,8 @@ sino_normalizations = {
"none": None,
"": None,
"chebyshev": "chebyshev",
"subtraction": "subtraction",
"division": "division",
}
cor_methods = {
......
......@@ -188,6 +188,7 @@ class ProcessConfigBase:
"""
raise ValueError("Base class")
def _get_rotation_axis_position(self):
self.dataset_info.axis_position = self.nabu_config["reconstruction"]["rotation_axis_position"]
......
......@@ -118,7 +118,7 @@ class PaganinPhaseRetrieval:
- If an positive integer `n` is provided, then `n` threads are used
- If a negative integer `n` is provided, then `n_max + n + 1` threads are used (so -1 to use all available threads)
- If `True` is provided, then `n_max - 2` threads are used
- If `False`or `None` is provided (default), then FFTW is not used
- If `False` or `None` is provided (default), then FFTW is not used
Important
----------
......
......@@ -273,10 +273,16 @@ class SinoNormalization:
"""
kinds = [
"chebyshev"
"chebyshev",
"subtraction",
"division",
]
operations = {
"subtraction": np.subtract,
"division": np.divide
}
def __init__(self, kind="chebyshev", sinos_shape=None, radios_shape=None):
def __init__(self, kind="chebyshev", sinos_shape=None, radios_shape=None, normalization_array=None):
"""
Initialize a SinoNormalization class.
......@@ -286,19 +292,48 @@ class SinoNormalization:
Normalization type. They can be the following:
- chebyshev: Each sinogram line is estimated by a Chebyshev polynomial
of degree 2. This estimation is then subtracted from the sinogram.
- subtraction: Each sinogram is subtracted with a user-provided array.
The array can be 1D (angle-independent) and 2D (angle-dependent)
- division: same as previously, but with a division operation.
Default is "chebyshev"
sinos_shape: tuple, optional
Shape of the sinogram or sinogram stack.
Either this parameter or 'radios_shape' has to be provided.
radios_shape: tuple, optional
Shape of the projections or projections stack.
Either this parameter or 'sinos_shape' has to be provided.
normalization_array: numpy.ndarray, optional
Normalization array when kind='subtraction' or kind='division'.
"""
self._get_shapes(sinos_shape, radios_shape)
self._set_kind(kind)
self._set_kind(kind, normalization_array)
_get_shapes = SinoBuilder._get_shapes
def _set_kind(self, kind):
def _set_kind(self, kind, normalization_array):
check_supported(kind, self.kinds, "sinogram normalization kind")
self.normalization_kind = kind
self._normalization_instance_method = self._normalize_chebyshev # default
if kind in ["subtraction", "division"]:
if not isinstance(normalization_array, np.ndarray):
raise ValueError(
"Expected 'normalization_array' to be provided as a numpy array for normalization kind='%s'" % kind
)
if normalization_array.shape[-1] != self.sinos_shape[-1]:
n_a, n_x = self.sinos_shape[-2:]
raise ValueError(
"Expected normalization_array to have shape (%d, %d) or (%d, )"
% (n_a, n_x, n_x)
)
self.norm_operation = self.operations[kind]
self._normalization_instance_method = self._normalize_op
self.normalization_array = normalization_array
#
# Chebyshev normalization
#
def _normalize_chebyshev_2D(self, sino):
output = sino # inplace
......@@ -325,21 +360,32 @@ class SinoNormalization:
def _normalize_chebyshev(self, sino):
if sino.ndim == 2:
res = self._normalize_chebyshev_2D(sino)
self._normalize_chebyshev_2D(sino)
else:
res = self._normalize_chebyshev_3D(sino)
return res
self._normalize_chebyshev_3D(sino)
return sino
def _get_norm_func_name(self, prefix=""):
return prefix + "normalize_" + self.normalization_kind
#
# Array subtraction/division
#
def _normalize_op(self, sino):
if sino.ndim == 2:
self.norm_operation(sino, self.normalization_array, out=sino)
else:
for i in range(sino.shape[0]):
self.norm_operation(sino[i], self.normalization_array, out=sino[i])
return sino
#
# Dispatch
#
def normalize(self, sino):
"""
Normalize a sinogram or stack of sinogram.
The process is done in-place, meaning that the sinogram content is overwritten.
"""
norm_func = getattr(self, self._get_norm_func_name(prefix="_"))
return norm_func(sino)
return self._normalization_instance_method(sino)
......@@ -155,9 +155,9 @@ CudaSinoProcessing = deprecated_class(
class CudaSinoNormalization(SinoNormalization):
def __init__(self, kind="chebyshev", sinos_shape=None, radios_shape=None, cuda_options=None):
def __init__(self, kind="chebyshev", sinos_shape=None, radios_shape=None, normalization_array=None, cuda_options=None):
super().__init__(
kind=kind, sinos_shape=sinos_shape, radios_shape=radios_shape,
kind=kind, sinos_shape=sinos_shape, radios_shape=radios_shape, normalization_array=normalization_array
)
self._get_shapes(sinos_shape, radios_shape)
self.cuda_processing = CudaProcessing(**(cuda_options or {}))
......@@ -166,41 +166,93 @@ class CudaSinoNormalization(SinoNormalization):
_get_shapes = SinoBuilder._get_shapes
def _init_cuda_normalization(self):
self._cuda_kernel = CudaKernel(
self._get_norm_func_name(),
filename=get_cuda_srcfile("normalization.cu"),
signature="Piii",
)
self._cuda_kernel_args = [
np.int32(self.n_x), np.int32(self.n_angles), np.int32(self.n_z)
]
blk = (1, 64, 16) # TODO tune ?
self._cuda_kernel_kwargs = {
"block": blk,
"grid": (1, int(updiv(self.n_angles, blk[1])), int(updiv(self.n_z, blk[2]))),
}
#
# Chebyshev normalization
#
def _init_cuda_normalization(self):
self._d_tmp = garray.zeros(self.sinos_shape[-2:], "f")
if self.normalization_kind == "chebyshev":
self._chebyshev_kernel = CudaKernel(
"normalize_chebyshev",
filename=get_cuda_srcfile("normalization.cu"),
signature="Piii",
)
self._chebyshev_kernel_args = [
np.int32(self.n_x), np.int32(self.n_angles), np.int32(self.n_z)
]
blk = (1, 64, 16) # TODO tune ?
self._chebyshev_kernel_kwargs = {
"block": blk,
"grid": (1, int(updiv(self.n_angles, blk[1])), int(updiv(self.n_z, blk[2]))),
}
elif self.normalization_array is not None:
normalization_array = self.normalization_array
# If normalization_array is 1D, make a 2D array by repeating the line
if normalization_array.ndim == 1:
normalization_array = np.tile(normalization_array, (self.n_angles, 1))
self._d_normalization_array = garray.to_gpu(normalization_array.astype("f"))
if self.normalization_kind == "subtraction":
generic_op_val = 1
elif self.normalization_kind == "division":
generic_op_val = 3
self._norm_kernel = CudaKernel(
"inplace_generic_op_2Dby2D",
filename=get_cuda_srcfile("ElementOp.cu"),
signature="PPii",
options=["-DGENERIC_OP=%d" % generic_op_val]
)
self._norm_kernel_args = [
self._d_normalization_array, np.int32(self.n_angles), np.int32(self.n_x)
]
blk = (32, 32, 1)
self._norm_kernel_kwargs = {
"block": blk,
"grid": (int(updiv(self.n_angles, blk[0])), int(updiv(self.n_x, blk[1])), 1)
}
def normalize(self, sinos):
def _normalize_chebyshev(self, sinos):
if sinos.flags.c_contiguous:
self._cuda_kernel(
sinos, *self._cuda_kernel_args, **self._cuda_kernel_kwargs
self._chebyshev_kernel(
sinos, *self._chebyshev_kernel_args, **self._chebyshev_kernel_kwargs
)
else:
d_tmp = garray.zeros(sinos.shape[1:], "f")
# This kernel seems to have an issue on arrays that are not C-contiguous.
# We have to process image per image.
nz = np.int32(1)
nthreadsperblock = (1, 32, 1) # TODO tune
nblocks = (1, int(updiv(self.n_angles, nthreadsperblock[1])), 1)
for i in range(sinos.shape[0]):
d_tmp[:] = sinos[i][:]
self._cuda_kernel(
d_tmp,
self._d_tmp[:] = sinos[i][:]
self._chebyshev_kernel(
self._d_tmp,
np.int32(self.n_x), np.int32(self.n_angles), np.int32(1),
grid=nblocks,
block=nthreadsperblock
)
sinos[i][:] = d_tmp[:]
sinos[i][:] = self._d_tmp[:]
return sinos
#
# Array subtraction/division
#
def _normalize_op(self, sino):
if sino.ndim == 2:
# Things can go wrong if "sino" is a non-contiguous 2D array
# But this should be handled outside this function, as the processing is in-place
self._norm_kernel(sino, *self._norm_kernel_args, **self._norm_kernel_kwargs)
else:
if sino.flags.forc:
# Contiguous 3D array. But pycuda wants the same shape for both operands.
for i in range(sino.shape[0]):
self._norm_kernel(sino[i], *self._norm_kernel_args, **self._norm_kernel_kwargs)
else:
# Non-contiguous 2D array. Make a temp. copy
for i in range(sino.shape[0]):
self._d_tmp[:] = sino[i][:]
self._norm_kernel(self._d_tmp, *self._norm_kernel_args, **self._norm_kernel_kwargs)
sino[i][:] = self._d_tmp[:]
return sino
......@@ -14,6 +14,8 @@ def bootstrap(request):
cls = request.cls
cls.sino = get_data("sino_refill.npy")
cls.tol = 1e-7
cls.norm_array_1D = np.arange(cls.sino.shape[-1]) + 1
cls.norm_array_2D = np.arange(cls.sino.size).reshape(cls.sino.shape) + 1
@pytest.mark.usefixtures("bootstrap")
......@@ -28,7 +30,7 @@ class TestSinoNormalization:
@pytest.mark.skipif(not(__has_pycuda__), reason="Need pycuda for sinogram normalization with cuda backend")
def test_dff_cuda(self):
def test_sino_normalization_cuda(self):
sino_proc = SinoNormalization(
kind="chebyshev", sinos_shape=self.sino.shape
)
......@@ -44,3 +46,54 @@ class TestSinoNormalization:
assert np.max(np.abs(res - ref)) < self.tol
def get_normalization_reference_result(self, op, normalization_arr):
# Perform explicit operations to compare with numpy.divide, numpy.subtract, etc
if op == "subtraction":
ref = self.sino - normalization_arr
elif op == "division":
ref = self.sino / normalization_arr
return ref
def test_sino_array_subtraction_and_division(self):
with pytest.raises(ValueError):
SinoNormalization(kind="subtraction", sinos_shape=self.sino.shape)
def compare_normalizations(normalization_arr, op):
sino_normalization = SinoNormalization(
kind=op, sinos_shape=self.sino.shape,
normalization_array=normalization_arr
)
sino = self.sino.copy()
sino_normalization.normalize(sino)
ref = self.get_normalization_reference_result(op, normalization_arr)
assert np.allclose(sino, ref), "operation=%s, normalization_array dims=%d" % (op, normalization_arr.ndim)
compare_normalizations(self.norm_array_1D, "subtraction")
compare_normalizations(self.norm_array_1D, "division")
compare_normalizations(self.norm_array_2D, "subtraction")
compare_normalizations(self.norm_array_2D, "division")
@pytest.mark.skipif(not(__has_pycuda__), reason="Need pycuda for sinogram normalization with cuda backend")
def test_sino_array_subtraction_cuda(self):
with pytest.raises(ValueError):
CudaSinoNormalization(kind="subtraction", sinos_shape=self.sino.shape)
def compare_normalizations(normalization_arr, op):
sino_normalization = CudaSinoNormalization(
kind=op, sinos_shape=self.sino.shape,
normalization_array=normalization_arr
)
sino = garray.to_gpu(self.sino)
sino_normalization.normalize(sino)
ref = self.get_normalization_reference_result(op, normalization_arr)
assert np.allclose(sino.get(), ref)
compare_normalizations(self.norm_array_1D, "subtraction")
compare_normalizations(self.norm_array_2D, "subtraction")
compare_normalizations(self.norm_array_1D, "division")
compare_normalizations(self.norm_array_2D, "division")
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