Commit 04807474 authored by Pierre Paleo's avatar Pierre Paleo

Merge branch 'release_2020.4' into 'master'

Release 2020.4.0

Closes #183, #186, and #184

See merge request !79
parents 4a159813 fd59f606
Pipeline #35269 passed with stages
in 12 minutes and 52 seconds
# Change Log
## 2020.4.0
This is a version adds a number of features.
### Application
- Automatic center of rotation estimation, working for both half-acquisition and regular scans
- Command line tool for splitting a NXTomo file by "z" series
- Volume histogram. Command line tool for merging histogram of multiple volumes.
- Enable to perform flat-field normalization with darks/flats from another dataset
- Sinogram normalization (baseline subtraction)
- Nabu does not need `tomwer_processes.h5` to get the "final" darks/refs anymore.
- Fix `fbp_filter_type = none`: now it will actually disable any filtering
- Fix auto-CoR not working when flatfield is disabled
### Library
- Add `misc.histogram` for computing partial histograms and merging them
- Add `preproc.sinogram.SinoNormalization`
- Fix double flat-field when `dff_sigma > 0`  which was giving nonsense results.
- Half-tomography: add support for center of rotation on the left side
## 2020.3.0
This is a release for ESRF User Service Mode restart.
nabu.misc.histogram module
.. automodule:: nabu.misc.histogram
......@@ -9,6 +9,7 @@ Submodules
nabu.resources.cli.histogram module
.. automodule:: nabu.resources.cli.histogram
nabu.resources.cli.nx\_z\_splitter module
.. automodule:: nabu.resources.cli.nx_z_splitter
......@@ -9,6 +9,8 @@ Submodules
nabu.resources.nxflatfield module
.. automodule:: nabu.resources.nxflatfield
......@@ -23,6 +23,7 @@ Submodules
......@@ -22,7 +22,7 @@ copyright = '2020, ESRF'
author = 'Pierre Paleo'
# The full version, including alpha/beta/rc tags
release = '2020.3.0'
#release = '2020.4.0'
# -- General configuration ---------------------------------------------------
......@@ -56,6 +56,18 @@ Configuration file: section `[preproc]`: `take_logarithm = 1`, `log_min_clip = 1
API: [CCDProcessing.Log](apidoc/nabu.preproc.ccd.rst#nabu.preproc.ccd.Log) and [CudaLog](apidoc/nabu.preproc.ccd_cuda.rst#nabu.preproc.ccd_cuda.CudaLog)
### Sinogram normalization
Sinograms can sometimes be "messy" for various reasons. For example, a synchrotron beam refill can take place during the acquisition, and not be compensated properly by flats.
In this case, you can "normalize" the sinogram to get rid of defects. Currently, only a baseline removal is implemented.
Mind that you probably lose quantativeness by using this additional normalization !
Configuration file: `[preproc]` : `sino_normalization = chebyshev`.
API: [SinoNormalization](apidoc/nabu.preproc.sinogram.rst#nabu.preproc.sinogram.SinoNormalization)
## Phase retrieval
Phase retrieval is still part of the "pre-processing", although it has a dedicated section in the [configuration file](nabu_config_file).
......@@ -85,6 +97,7 @@ Configuration file: section `[phase]`: `unsharp_coeff = 1.` and `unsharp_sigma
Setting `coeff` to zero (default) disables unsharp masking.
## Reconstruction
Tomographic reconstruction is the process of reconstructing the volume from projections/sinograms.
......@@ -107,6 +120,16 @@ The purpose of this class is to quickly reconstruct slices over the three main a
The Reconstructor enables to reconstruct slices/region of interest without reconstructing the whole volume.
## Post-processing
### Histogram
Nabu can compute the histogram of the reconstucted (sub-) volume. As the volume usually does not fit in memory, the histogram is computed by parts, and the final histogram is obtained by merging partial histograms.
Configuration file: section `[postproc]`: `output_histogram = 1`, `histogram_bins = 1000000`.
API : [PartialHistogram](apidoc/nabu.misc.histogram.rst#nabu.misc.histogram.PartialHistogram) and [VolumeHistogram](apidoc/nabu.misc.histogram.rst#nabu.misc.histogram.VolumeHistogram)
## File formats
### HDF5
......@@ -118,7 +141,7 @@ When a [multi-stage reconstruction]( is performed, the volume is rec
### Tiff
Reconstruction can be saved as tiff files by specifying `file_format = tiff` in the configuration `[output]` section.
In the current version (2020.3), tiff support has the following limitations:
Mind that tiff support currently has the following limitations:
- One file per slice
- Data is saved as `float32` data type, no normalization
- No metadata (configuration used to obtain the reconstruction, date, version,...)
......@@ -126,7 +149,7 @@ In the current version (2020.3), tiff support has the following limitations:
### Jpeg2000
Reconstruction can be saved as jpeg2000 files by specifying `file_format = jpeg2000` in the configuration `[output]` section.
In the current version (2020.3), jpeg2000 support has the following limitations:
Mind that jpeg2000 support currently has the following limitations:
- When exporting to `uint16` data type (mandatory for jpeg2000), the normalization from `float32` to `uint16` is done slice-wise instead of volume-wise. This is slightly less accurate.
- Only lossless compression is supported. In the future, compression will be tunable through Nabu configuration.
- No metadata (configuration used to obtain the reconstruction, date, version,...)
......@@ -138,6 +161,14 @@ In the current version (2020.3), jpeg2000 support has the following limitations:
Nabu provides a method to find the half-shift between two images. The center of axial vertical rotation is obtained when the fist image is a radiography at the rotation angle 0 and the second image is given by the radiography at the rotation angle 180 after flipping the image horizontally. The rotation axis position is the center of the image plus the found shift.
Configuration file: section `[reconstruction]`, key `rotation_axis_position`. Values can be:
- A number (known CoR)
- Empty: the CoR is set to the middle of the detector
- `centered` (or `auto`): this searches for a CoR around the center, but does not work for half-acquisition
- `global`: new method, should return a CoR in any setting.
API: [CenterOfRotation](apidoc/nabu.preproc.alignment.rst#nabu.preproc.alignment.CenterOfRotation)
### Detector Translation Along the Beam
......@@ -16,6 +16,7 @@ Quick links:
Latest news:
* October 2020: `Version 2020.4.0 released <>`_
* August 2020: `Version 2020.3.0 released <>`_
* June 2020: `Version 2020.2.0 released <>`_
......@@ -82,6 +83,7 @@ Release history
.. toctree::
:maxdepth: 3
......@@ -14,18 +14,8 @@ pip install git+
The above `pip install` command should automatically install the nabu dependencies.
For the record, the mandatory dependencies are the following:
- numpy
- pytest
- psutil
- [silx](
- [tomoscan](
All of them can be installed with `pip`.
There are optional dependencides enabling various features:
- Computations acceleration (GPU/manycore CPU): `pycuda`, `pyopencl`
- Computations distribution: `distributed`, `dask_jobqueue`
Please note that Nabu supports Python >= 3.5.
......@@ -52,6 +52,7 @@ Each section describe a usual processing step. In the current version, the avail
- `preproc`: phase retrieval, CCD corrections, ...
- `phase`: phase retrieval
- `reconstruction`: tomography reconstruction
- `postproc`: post-processing: histogram
- `output`: output data description
- `resources`: computing resources description
- `about`: extra information about the current configuration file
......@@ -71,6 +72,6 @@ Yes. Nabu is foremost a library, meaning that all its component can be accessed
### Compatibility policy
At this early development stage, it is not entirely clear which keys should be in the configuration file. In the future, some keys might be removed, added, or have their name changed. Also, with new features coming, new keys/sections might appear in the configuration file.
During the development of Nabu, some features will be added, leading to new keys in the configuration file. Besides, some keys might be renamed or even deleted if deemed necessary.
Thus, the configuration file might change according to the version of nabu. For this reason, there is a section `[about]` in the configuration file, which defines a `nabu_version` and a `nabu_config_version`. A given version of the nabu software (`nabu_version`) will be compatible with a certain range of nabu configuration files (`nabu_config_version`).
In any case, **a configuration file from an "old" version of nabu will be supported by newer versions** for some time, with a deprecation warning when using obsolete keys.
\ No newline at end of file
# Version 2020.4.0
Version 2020.4.0 adds a number of features, notably command line tools for manipulating datasets.
``` note:: The changelog is available at
## Highlights
This section highlights some of the available [features](
### Fully automatic center of rotation estimation
Various methods exist for estimating automatically the center of rotation (CoR). Nabu had one, but it did not work for half tomography acquisitions. Another method is now available and should work with any kind of acquisition.
The paramater is `rotation_axis_position` (in section`[reconstruction]`), the value can be:
- A number (known CoR)
- Empty: the CoR is set to the middle of the detector
- `centered` (was `auto`): this searches for a CoR around the center, but does not work for half-acquisition
- `global`: new method, should return a CoR in any setting.
### Histograms
Nabu can now compute a histogram of the reconstructed volume. This histogram can serve for further volume processing (conversion to uint16 for example). The computation is done while the data is still in memory to avoid extraneous disk reads/writes.
The option is available in a new section `[postproc]` with the parameter `output_histogram = 1`.
A command line tool is also available to merge the histogram of multiple volumes: `nabu-histogram <file1> <file2> ... <output>`.
### Use darks/flats from another dataset
Sometimes, flats cannot be acquired properly during an acquisition. In order to still have the ability to reconstruct the data, a common workaround is to use flats from another dataset.
Suppose you want to reconstruct a dataset which does not have flats/darks. In this case, nabu will detect that these images are missing and won't do flat-field correction. However you can force the flat-field normalization with `flatfield_enabled = forced` in the `[preproc]` section, and by providing `processes_file = /path/to/nabu_processes.h5`. This `nabu_processes.h5` file (or `tomwer_processes.h5`) can be generated either by tomwer or nabu.
### Sinogram normalization
A sinogram normalization is now available in the section `[preproc]` : `sino_normalization = chebyshev`. If enabled, each line of the sinogram undergoes a baseline removal. This can be useful to correct a "refill" taking place during the scan. See [this discussion]( for more details.
### Split a NXTomo file by "z series"
Z-series is the name given to a multi-stage scan, i.e scanning a volume in several stages by moving the sample vertically after each stage (this is different from helical tomography).
For the time being, the CLI tool `nxtomomill` merges everything in a single file.
Projections have a varying "z" identified by the `entry/sample/z_translation` key (in conformity to the Nexus-Tomo standard).
The command line tool `nabu-zsplit` enables to split a file into distinct "z" to reconstruct each volume individually.
``` warning:: This is a temporary solution. This will be natively be handled by nxtomomill in the future.
\ No newline at end of file
__version__ = "2020.3.1"
__version__ = "2020.4.0"
__nabu_modules__ = [
......@@ -358,12 +358,13 @@ class FullFieldPipeline:
@use_options("double_flatfield", "double_flatfield")
def _init_double_flatfield(self):
options = self.processing_options["double_flatfield"]
avg_is_on_log = (options["sigma"] is not None)
self.double_flatfield = self.DoubleFlatFieldClass(
......@@ -395,7 +396,7 @@ class FullFieldPipeline:
if "unsharp_mask" not in self.processing_steps:
self.register_callback("phase", self._reshape_radios_after_phase)
self.register_callback("phase", FullFieldPipeline._reshape_radios_after_phase)
@use_options("unsharp_mask", "unsharp_mask")
def _init_unsharp(self):
......@@ -405,7 +406,7 @@ class FullFieldPipeline:
options["unsharp_sigma"], options["unsharp_coeff"],
mode="reflect", method="gaussian"
self.register_callback("unsharp_mask", self._reshape_radios_after_phase)
self.register_callback("unsharp_mask", FullFieldPipeline._reshape_radios_after_phase)
@use_options("take_log", "mlog")
def _init_mlog(self):
......@@ -458,6 +459,7 @@ class FullFieldPipeline:
self._rec_roi = (x_s, x_e, y_s, y_e)
self._allocate_recs(y_e - y_s, x_e - x_s)
@use_options("reconstruction", "reconstruction")
def _init_reconstruction(self):
options = self.processing_options["reconstruction"]
......@@ -471,6 +473,8 @@ class FullFieldPipeline:
rot_center = options["rotation_axis_position"]
if options.get("cor_estimated_auto", False):"Estimated center of rotation: %.2f" % rot_center)
if self.sino_builder._halftomo_flip:
rot_center = self.sino_builder.rot_center
self.reconstruction = self.FBPClass(
......@@ -484,7 +488,7 @@ class FullFieldPipeline:
"axis_correction": options["axis_correction"],
if options["fbp_filter_type"] == "none":
if options["fbp_filter_type"] is None:
self.reconstruction.fbp = self.reconstruction.backproj
@use_options("histogram", "histogram")
......@@ -91,11 +91,11 @@ class CudaFullFieldPipeline(FullFieldPipeline):
def _register_callbacks(self):
self.register_callback("read_chunk", self._read_data_callback)
self.register_callback("read_chunk", CudaFullFieldPipeline._read_data_callback)
if self.reconstruction is not None:
self.register_callback("reconstruction", self._rec_callback)
self.register_callback("reconstruction", CudaFullFieldPipeline._rec_callback)
if self.writer is not None:
self.register_callback("save", self._saving_callback)
self.register_callback("save", CudaFullFieldPipeline._saving_callback)
......@@ -33,7 +33,7 @@ def pipeline_step(step_attr, step_desc):
self.logger.debug("End " + step_desc)
callback = self._callbacks.get(self._steps_component2name[step_attr], None)
if callback is not None:
return res
return wrapper
return decorator
......@@ -32,18 +32,18 @@ class MedianFilter(CudaProcessing):
Size of the median filter, in the format (y, x).
mode: str
Boundary handling mode. Available modes are:
"reflect": cba|abcd|dcb
"nearest": aaa|abcd|ddd
"wrap": bcd|abcd|abc
"constant": 000|abcd|000
- "reflect": cba|abcd|dcb
- "nearest": aaa|abcd|ddd
- "wrap": bcd|abcd|abc
- "constant": 000|abcd|000
Default is "reflect".
threshold: float, optional
Threshold for the "thresholded median filter".
A thresholded median filter only replaces a pixel value by the median
if this pixel value is greater or equal than median + threshold.
Cuda Parameters
Please refer to the documentation of the CudaProcessing class for
the other parameters.
......@@ -86,3 +86,20 @@ __global__ void nlog(float* array, int Nx, int Ny, int Nz, float clip_min, float
array[pos] = -logf(val);
// Reverse elements of a 2D array along "x", i.e:
// arr = arr[:, ::-1]
// launched with grid (Nx/2, Ny)
__global__ void reverse2D_x(float* array, int Nx, int Ny) {
uint x = blockDim.x * blockIdx.x + threadIdx.x;
uint y = blockDim.y * blockIdx.y + threadIdx.y;
if ((x >= Nx/2) || (y >= Ny)) return;
uint pos = y*Nx + x;
uint pos2 = y*Nx + (Nx - 1 - x);
float tmp = array[pos];
array[pos] = array[pos2];
array[pos2] = tmp;
......@@ -105,6 +105,11 @@ def get_first_hdf5_entry(fname):
return entry
def hdf5_entry_exists(fname, entry):
with HDF5File(fname, "r") as fid:
res = fid.get(entry, None) is not None
return res
def get_h5_value(fname, h5_path, default_ret=None):
with HDF5File(fname, "r") as fid:
......@@ -63,8 +63,12 @@ class CudaDoubleFlatField(DoubleFlatField, CudaProcessing):
return o
def _proc_mlog(x, o):
cumath.log(x, out=o)
def _proc_mlog(x, o, min_clip=None):
if min_clip is not None:
garray.maximum(x, min_clip, output=o)
cumath.log(o, out=o)
cumath.log(x, out=o)
o *= -1
return o
......@@ -71,7 +71,15 @@ class SinoProcessing(object):
self.rot_center = rot_center
self._rot_center_int = int(round(self.rot_center))
# If CoR is on the left: "flip" the logic
self._halftomo_flip = False
rc = self._rot_center_int
if rc < (self.n_x - 1)/2:
rc = self.n_x - 1 - rc
self._rot_center_int = rc
self.rot_center = self.n_x - self.rot_center
self._halftomo_flip = True
def _set_halftomo(self, halftomo):
self.halftomo = halftomo
......@@ -123,7 +131,12 @@ class SinoProcessing(object):
sinos = np.zeros(self.sinos_halftomo_shape, dtype=np.float32)
for i in range(n_z):
sinos[i][:] = convert_halftomo(radios[:, i, :], self._rot_center_int)
radio = radios[:, i, :]
if self._halftomo_flip:
radio = radio[:, ::-1]
sinos[i][:] = convert_halftomo(radio, self._rot_center_int)
if self._halftomo_flip:
sinos[i][:] = sinos[i][:, ::-1]
return sinos
......@@ -42,6 +42,20 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing):
d = self.n_x - rc # will have to be adapted for varying axis pos
self.halftomo_weights = np.linspace(0, 1, d, endpoint=True, dtype="f")
self.d_halftomo_weights = garray.to_gpu(self.halftomo_weights)
if self._halftomo_flip:
self.xflip_kernel = CudaKernel(
blk = (32, 32, 1)
self._xflip_blksize = blk
self._xflip_gridsize_1 = (
updiv(self.n_x, blk[0]),
updiv(self.n_angles, blk[1]),
self._xflip_gridsize_2 = self._halftomo_gridsize
# Overwrite parent method
......@@ -80,15 +94,28 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing):
sinos = garray.zeros(self.sinos_halftomo_shape, dtype=np.float32)
# need to use a contiguous 2D, array otherwise kernel does not work
d_radio = radios[:, 0, :].copy()
for i in range(n_z):
d_radio[:] = radios[:, i, :].copy()
if self._halftomo_flip:
d_radio, n_x, n_a,
grid=self._xflip_gridsize_1, block=self._xflip_blksize
radios[:, i, :].copy(), # need to copy this 2D, array otherwise kernel does not work
n_a, n_x, rc,
if self._halftomo_flip:
sinos[i], 2*rc, n_a2,
grid=self._xflip_gridsize_2, block=self._xflip_blksize
return sinos
......@@ -125,7 +152,3 @@ class CudaSinoNormalization(SinoNormalization, CudaProcessing):
sino, *self._cuda_kernel_args, **self._cuda_kernel_kwargs
return sino
......@@ -5,12 +5,13 @@ import h5py
from nabu.testutils import compare_arrays, utilstest
from nabu.preproc.sinogram import SinoProcessing, convert_halftomo
from nabu.cuda.utils import __has_pycuda__
if __has_pycuda__:
import pycuda.gpuarray as garray
from nabu.preproc.sinogram_cuda import CudaSinoProcessing
def bootstrap(request):
cls = request.cls
radio, sino_ref, cor = get_data_h5("halftomo.h5")
......@@ -26,37 +27,87 @@ def bootstrap(request):
def get_data_h5(*dataset_path):
dataset_relpath = os.path.join(*dataset_path)
dataset_path = utilstest.getfile(dataset_relpath)
with h5py.File(dataset_path, 'r') as hf:
with h5py.File(dataset_path, "r") as hf:
radio = hf["entry/radio/results/data"][()]
sino = hf["entry/sino/results/data"][()]
cor = hf["entry/sino/configuration/configuration/rotation_axis_position"][()]
return radio, sino, cor
class TestHalftomo:
def test_halftomo(self):
sino_processing = SinoProcessing(
radios_shape=self.radios.shape, rot_center=self.rot_center, halftomo=True
sinos_halftomo = sino_processing.radios_to_sinos(self.radios)
_, err = compare_arrays(sinos_halftomo[0], self.sino_ref, self.tol, return_residual=True)
assert err < self.tol, "Something wrong with SinoProcessing.radios_to_sino, halftomo=True"
_, err = compare_arrays(
sinos_halftomo[0], self.sino_ref, self.tol, return_residual=True
assert (
err < self.tol
), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True"
@pytest.mark.skipif(not(__has_pycuda__), reason="Need pycuda for this test")
@pytest.mark.skipif(not (__has_pycuda__), reason="Need pycuda for this test")
def test_cuda_halftomo(self):
sino_processing = CudaSinoProcessing(
radios_shape=self.radios.shape, rot_center=self.rot_center, halftomo=True
d_radios = garray.to_gpu(self.radios)
d_sinos = garray.zeros(sino_processing.sinos_halftomo_shape, "f")
sino_processing.radios_to_sinos(d_radios, output=d_sinos)
sino_halftomo = d_sinos.get()[0]
_, err = compare_arrays(sino_halftomo, self.sino_ref, self.tol, return_residual=True)
assert err < self.tol, "Something wrong with SinoProcessing.radios_to_sino, halftomo=True"
_, err = compare_arrays(
sino_halftomo, self.sino_ref, self.tol, return_residual=True
assert (
err < self.tol
), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True"
def _flip_array(arr):
if arr.ndim == 2:
return np.fliplr(arr)
res = np.zeros_like(arr)
for i in range(arr.shape[0]):
res[i] = np.fliplr(arr[i])
return res
def test_halftomo_left(self):
na, nz, nx = self.radios.shape
left_cor = nx - 1 - self.rot_center
radios = self._flip_array(self.radios)
sino_processing = SinoProcessing(
radios_shape=radios.shape, rot_center=left_cor, halftomo=True
sinos_halftomo = sino_processing.radios_to_sinos(radios)
_, err = compare_arrays(
assert (
err < self.tol
), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True"
@pytest.mark.skipif(not (__has_pycuda__), reason="Need pycuda for this test")
def test_cuda_halftomo_left(self):
na, nz, nx = self.radios.shape
left_cor = nx - 1 - self.rot_center
radios = self._flip_array(self.radios)
sino_processing = CudaSinoProcessing(
radios_shape=radios.shape, rot_center=left_cor, halftomo=True
d_radios = garray.to_gpu(radios)
d_sinos = garray.zeros(sino_processing.sinos_halftomo_shape, "f")
sino_processing.radios_to_sinos(d_radios, output=d_sinos)