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

Merge branch 'release_2020.2' into 'master'

Release 2020.2

See merge request !47
parents 5b605233 1d53ff4b
Pipeline #28499 passed with stages
in 3 minutes and 40 seconds
# Change Log
# 2020.1.0: 2020/04/29
## 2020.2.0
This is a "preview release", preparing an important release foreseen for September 2020.
### Library features
- Vertical translations
- Double flat-field (radios-based rings artefacts removal)
- Support for half tomography
- `preproc.alignment` module
- Support for flat-field normalization with more than 2 flats
### Application
- Take configuration file `start_z` and `end_z` into account for `nabu` CLI
- Phase retrieval: `method = none` is now default (no phase retrieval by default)
- Add `overwrite_results`
- Add support for `rotation_axis_position = auto` (except for half-tomo)
### Internal
- Re-wrote internal processing pipeline with a much simpler design
- NXresults: config is now a H5 dataset. "image stack" default interpretation.
## 2020.1.0: 2020/04/29
This is the first official release of Nabu. Please see the [features overview](www.silx.org/pub/nabu/doc/features.html)
\ No newline at end of file
# About Nabu
Nabu is a tomography processing software being developed at [ESRF](https://www.esrf.eu) by the Data Analysis Unit. It is part of the [new ESRF tomography software suite](https://tomotools.gitlab-pages.esrf.fr/tomo-bridge-docs/data-proc/projects).
## Why Nabu
The European Synchrotron has several tomography beamlines. Each of them use dedicated software, which over the years led a variety of different tools spread over the beamlines with poor maintainability. This is summarized in [ESRF current situation for tomography software](https://tomotools.gitlab-pages.esrf.fr/tomo-bridge-docs/data-proc/situation/#current-situation-limitations-and-opportunities).
Nabu is an effort to unify tomography software in a new toolkit with the following requirements:
- Library of tomography processing, with "applications" built on top of it, usable by both non-experts and power-users
- High performance processing (parallelization with Cuda/OpenCL, computations distribution, memory re-use)
- Support of [multiple techniques](https://tomotools.gitlab-pages.esrf.fr/tomo-bridge-docs/#tomography-techniques-at-esrf), not only absorption and phase contrast
- Extensively documented
- Focus on maintainability with a [bus factor](https://en.wikipedia.org/wiki/Bus_factor) greater than one
- Compatible with ESRF legacy software, progressively replacing it
Nabu does not aim at being the new universal tomography reconstruction software. Well-established software like [Astra](http://www.astra-toolbox.com), [tomopy](http://tomopy.readthedocs.io), [Savu](https://github.com/DiamondLightSource/Savu) and [UFO](https://github.com/ufo-kit) have an extensive set of features. Nabu foremost focuses on ESRF needs, while being designed so that it can be re-used in other projects.
## Project management
### Development spot
The various projects are hosted on [gitlab.esrf.fr](https://gitlab.esrf.fr/tomotools).
### Members and meeting minutes
The [new ESRF tomography software suite](https://tomotools.gitlab-pages.esrf.fr/tomo-bridge-docs/data-proc/projects) is developed by ESRF Data Analysis Unit.
The weekly meeting minutes [are available here](https://gitlab.esrf.fr/tomotools/minutes/-/tree/master).
### Roadmap
[Roadmap for Nabu developments](https://tomotools.gitlab-pages.esrf.fr/tomo-bridge-docs/data-proc/roadmap)
......@@ -9,8 +9,8 @@ Submodules
nabu.app.fullfield
nabu.app.fullfield_cuda
nabu.app.logger
nabu.app.simple_pipeline
nabu.app.utils
Module contents
---------------
......
nabu.app.logger module
======================
nabu.app.utils module
=====================
.. automodule:: nabu.app.logger
.. automodule:: nabu.app.utils
:members:
:undoc-members:
:show-inheritance:
nabu.misc.fourier_filters module
========================
nabu.misc.fourier\_filters module
=================================
.. automodule:: nabu.misc.fourier_filters
:members:
......
......@@ -8,6 +8,7 @@ Submodules
:maxdepth: 4
nabu.misc.binning
nabu.misc.fourier_filters
nabu.misc.unsharp
nabu.misc.unsharp_cuda
nabu.misc.unsharp_opencl
......
......@@ -17,6 +17,7 @@ Submodules
nabu.preproc.rings
nabu.preproc.rings_cuda
nabu.preproc.shift
nabu.preproc.shift_cuda
nabu.preproc.sinogram
nabu.preproc.sinogram_cuda
......
nabu.preproc.shift\_cuda module
===============================
.. automodule:: nabu.preproc.shift_cuda
:members:
:undoc-members:
:show-inheritance:
nabu.resources.cor module
=========================
.. automodule:: nabu.resources.cor
:members:
:undoc-members:
:show-inheritance:
nabu.resources.logger module
============================
.. automodule:: nabu.resources.logger
:members:
:undoc-members:
:show-inheritance:
......@@ -16,9 +16,11 @@ Submodules
:maxdepth: 4
nabu.resources.computations
nabu.resources.cor
nabu.resources.dataset_analyzer
nabu.resources.dataset_validator
nabu.resources.gpu
nabu.resources.logger
nabu.resources.machinesdb
nabu.resources.nabu_config
nabu.resources.params
......
......@@ -22,7 +22,7 @@ copyright = '2020, Pierre Paleo'
author = 'Pierre Paleo'
# The full version, including alpha/beta/rc tags
release = '2020.1.0'
release = '2020.2.0'
# -- General configuration ---------------------------------------------------
......
......@@ -12,7 +12,7 @@ API : [nabu.preproc](apidoc/nabu.preproc)
This is usually the first step done on the radios. Each radio is normalized by the incoming beam intensity and CCD dark current as `radio = (radio - dark)/(flat - dark)`.
Usually, several "flats" images are acquired. Nabu automatically determines which flat image to use depending on the radio angle.
Usually, several "flats" images are acquired. Nabu automatically determines which flat image to use depending on the radio angle. A linear interpolation is done between flats to estimate the "current" flat to use.
Configuration file: section `[preproc]`: `flatfield_enabled = 1`.
......@@ -27,6 +27,9 @@ Configuration file: section `[preproc]`: `ccd_filter_enabled = 1` and `ccd_filte
API: [CCDCorrection](apidoc/nabu.preproc.ccd.rst#nabu.preproc.ccd.CCDCorrection) and [CudaCCDCorrection](apidoc/nabu.preproc.ccd_cuda.rst#nabu.preproc.ccd_cuda.CudaCCDCorrection)
``` warning:: In version 2020.2.0, this feature is not available through the command line interface, meaning that you cannot use it via the configuration file yet. It is only available in the library.
```
### Double flat-field
This is a projections-based rings artefacts removal.
......
......@@ -14,6 +14,10 @@ Quick links:
* `Nabu development place <https://gitlab.esrf.fr/tomotools/nabu/>`_
* Project maintainer: `Pierre Paleo <mailto:pierre.paleo@esrf.fr>`_
Latest news:
* June 2020: `Version 2020.2.0 released <v2020_2_0.md>`_
Getting started
---------------
......@@ -25,6 +29,7 @@ Getting started
quickstart.md
definitions.md
nabu_config_file.md
about.md
Features
......@@ -69,6 +74,14 @@ API reference
apidoc/nabu.rst
Release history
----------------
.. toctree::
:maxdepth: 3
v2020_2_0.md
Indices and tables
......
......@@ -22,6 +22,11 @@ nabu nabu.conf
More options are available, see `nabu --help`.
For example, to reconstruct slices 1000 to 1100: `nabu nabu.conf --slice 1000-1100`
``` warning:: In version 2020.2.0, running "nabu nabu.conf" is likely to cause memory errors as nabu will attempt at reconstructing all the dataset at once. You should use always use the slice parameter or specify start_z and end_z in the configuration file.
```
## See also
......
# Version 2020.2.0
Version 2020.2.0 is a "preview" of a major release foreseen for September 2020. It brings new features since 2020.1.0.
``` note:: The changelog is available at https://gitlab.esrf.fr/tomotools/nabu/-/blob/master/CHANGELOG.md
```
## Highlights
This section highlights some of the available [features](features.md).
### Binning
For a quick reconstruction, you can perform binning on the data. There are three parameters in the configuration file:
- `binning`: binning in the horizontal direction. Each reconstructed slice will have its dimensions divided by this factor. For now only factors 2 and 3 are supported.
- `binning_z`: binning in the vertical direction. The number of reconstructed slices will be divided by this number. For now only factors 2 and 3 are supported.
- `projections_subsampling`: read one projection out of 2 or 3. This accelerates the processing but results in lower quality data.
### Automatic center of rotation estimation
The default `rotation_axis_position` value (if left blank) is half of the detector width. The value `rotation_axis_position = auto` can now be provided, so that the value will be automatically estimated.
In the current version, automatic estimation lacks robustness and does not work for half tomography.
### Half tomography
Half tomography reconstruction is supported by setting the parameter `enable_halftomo` to 1.
### Double flat-field
[Double flat-field](features.md#double-flat-field) is a projections-based rings artefacts removal. The relevant configuration key is `double_flatfield_enabled`.
## Foreseen developments
Most of the following foreseen developments should materialize for "September Milestone".
### Input/output
- "merge" HDF5 reconstructions by creating a master file
- Allow other formats than HDF5 (ex. tiff, jp2k, npy)
- Enable to dump intermediate results (pre-processed radios/sinograms) and run the processing from a previous intermediate result
- Generate flats/darks instead of relying on a `tomwer_processes.h5` file
### Reconstruction
- Determine the minimum number of detector rows to read ("Paganin margin")
- Reconstruction with overlapping vertical regions for phase retrieval
- More robust Center of Rotation estimation methods (notably for half-tomography)
- Sample/detector shifts. This is theoretically already implemented, but real datasets for tests are missing.
- Sinogram normalization (background subtraction)
- Iterative reconstruction (need a forward projector implementation)
### Computations distribution
Distributing the computation will be a crucial feature of Nabu.
- Distribution on the local machine
- Distribution on the computing cluster (via the SLURM task dispatcher)
......@@ -9,7 +9,7 @@
/**
* In-place flat-field normalization.
* In-place flat-field normalization with linear interpolation.
* This kernel assumes that all the radios are loaded into memory
* (although not necessarily the full radios images)
* and in radios[x, y z], z in the radio index
......@@ -20,7 +20,8 @@
* Nx: number of pixel horizontally in the radios
* Nx: number of pixel vertically in the radios
* Nx: number of radios
* flats_indices: indices of flats, in sorted order
* flats_indices: indices of flats to fetch for each radio
* flats_weights: weights of flats for each radio
* darks_indices: indices of darks, in sorted order
**/
__global__ void flatfield_normalization(
......@@ -31,43 +32,30 @@ __global__ void flatfield_normalization(
int Ny,
int Nz,
int* flats_indices,
int* darks_indices,
int* radios_indices
float* flats_weights
) {
uint x = blockDim.x * blockIdx.x + threadIdx.x;
uint y = blockDim.y * blockIdx.y + threadIdx.y;
uint z = blockDim.z * blockIdx.z + threadIdx.z;
if ((x >= Nx) || (y >= Ny) || (z >= Nz)) return;
uint pos = (z*Ny+y)*Nx + x;
int radio_idx = radios_indices[z];
float dark_val = 0.0f, flat_val = 1.0f;
#if N_FLATS == 1
flat_val = flats[y*Nx + x];
#else
// interpolation between 2 flats
for (int i = 0; i < N_FLATS-1; i++) {
int ind_prev = flats_indices[i];
int ind_next = flats_indices[i+1];
if (ind_prev >= radio_idx) {
flat_val = flats[(i*Ny+y)*Nx + x];
break;
}
else if (ind_prev < radio_idx && radio_idx < ind_next) {
// Linear interpolation
// TODO nearest interpolation
int delta = ind_next - ind_prev;
float w1 = 1.0f - (radio_idx*1.0f - ind_prev) / delta;
float w2 = 1.0f - (ind_next*1.0f - radio_idx) / delta;
flat_val = w1 * flats[(i*Ny+y)*Nx + x] + w2 * flats[((i+1)*Ny+y)*Nx + x];
break;
}
else if (ind_next <= radio_idx) {
flat_val = flats[((i+1)*Ny+y)*Nx + x];
break;
}
int prev_idx = flats_indices[z*2 + 0];
int next_idx = flats_indices[z*2 + 1];
float w1 = flats_weights[z*2 + 0];
float w2 = flats_weights[z*2 + 1];
if (next_idx == -1) {
flat_val = flats[(prev_idx*Ny+y)*Nx + x];
}
else {
flat_val = w1 * flats[(prev_idx*Ny+y)*Nx + x] + w2 * flats[(next_idx*Ny+y)*Nx + x];
}
#endif
#if (N_DARKS == 1)
dark_val = darks[y*Nx + x];
......
......@@ -387,7 +387,7 @@ class ChunkReader(object):
self.binning = binning
binning_function = get_binning_function(binning[::-1]) # user provides (binning_x, binning_y)
if binning_function is None:
raise NotImplementedError("Binning factor %s is not implemented yet" % binning)
raise NotImplementedError("Binning factor %s is not implemented yet" % str(binning))
self.binning_function = binning_function
......
......@@ -142,7 +142,7 @@ class FlatField(CCDProcessing):
def _set_parameters(self, radios_shape, flats, darks, radios_indices, interpolation):
self._set_radios_shape(radios_shape)
self._check_frames(flats, "flats", 1, 2)
self._check_frames(flats, "flats", 1, 9999)
self._check_frames(darks, "darks", 1, 1)
self.flats = flats
self.darks = darks
......@@ -171,6 +171,7 @@ class FlatField(CCDProcessing):
for flat_idx, flat_url in self.flats.items():
self.flats_arr[flat_idx] = self.flats_reader.get_data(flat_url)
self._sorted_flat_indices = sorted(self.flats.keys())
self._flat2arrayidx = {k: v for k, v in zip(self._sorted_flat_indices, np.arange(self.n_flats))}
def _load_darks(self):
self.n_darks = len(self.darks)
......
......@@ -33,6 +33,7 @@ class CudaFlatField(FlatField):
)
self._set_cuda_options(cuda_options)
self._init_cuda_kernels()
self._precompute_flats_indices_weights()
self._load_flats_and_darks_on_gpu()
def _set_cuda_options(self, user_cuda_options):
......@@ -51,7 +52,7 @@ class CudaFlatField(FlatField):
self.cuda_kernel = CudaKernel(
"flatfield_normalization",
self._cuda_fname,
signature="PPPiiiPPP",
signature="PPPiiiPP",
options=[
"-DN_FLATS=%d" % self.n_flats,
"-DN_DARKS=%d" % self.n_darks,
......@@ -60,26 +61,41 @@ class CudaFlatField(FlatField):
self._nx = np.int32(self.shape[1])
self._ny = np.int32(self.shape[0])
def _precompute_flats_indices_weights(self):
flats_idx = np.zeros((self.n_radios, 2), dtype=np.int32)
flats_weights = np.zeros((self.n_radios, 2), dtype=np.float32)
for i, idx in enumerate(self.radios_indices):
prev_next = self.get_previous_next_indices(self._sorted_flat_indices, idx)
if len(prev_next) == 1: # current index corresponds to an acquired flat
weights = (1, 0)
f_idx = (self._flat2arrayidx[prev_next[0]], -1)
else: # interpolate
prev_idx, next_idx = prev_next
delta = next_idx - prev_idx
w1 = 1 - (idx - prev_idx) / delta
w2 = 1 - (next_idx - idx) / delta
weights = (w1, w2)
f_idx = (self._flat2arrayidx[prev_idx], self._flat2arrayidx[next_idx])
flats_idx[i] = f_idx
flats_weights[i] = weights
self.flats_idx = flats_idx
self.flats_weights = flats_weights
def _load_flats_and_darks_on_gpu(self):
# Flats
self.d_flats = garray.zeros((self.n_flats,) + self.shape, np.float32)
# ~ for i, flat_arr in enumerate(self.flats_arr.values()):
for i, flat_idx in enumerate(self._sorted_flat_indices):
self.d_flats[i].set(np.ascontiguousarray(self.flats_arr[flat_idx], dtype=np.float32))
self.d_flats_indices = garray.to_gpu(
np.array(self._sorted_flat_indices, dtype=np.int32)
)
# Darks
self.d_darks = garray.zeros((self.n_darks,) + self.shape, np.float32)
# ~ for i, dark_arr in enumerate(self.darks_arr.values()):
for i, dark_idx in enumerate(self._sorted_dark_indices):
self.d_darks[i].set(np.ascontiguousarray(self.darks_arr[dark_idx], dtype=np.float32))
self.d_darks_indices = garray.to_gpu(
np.array(self._sorted_dark_indices, dtype=np.int32)
)
# Radios indices
self.d_radios_indices = garray.to_gpu(self.radios_indices)
# Indices
self.d_flats_indices = garray.to_gpu(self.flats_idx)
self.d_flats_weights = garray.to_gpu(self.flats_weights)
def normalize_radios(self, radios):
"""
......@@ -105,8 +121,7 @@ class CudaFlatField(FlatField):
self._ny,
np.int32(self.n_radios),
self.d_flats_indices,
self.d_darks_indices,
self.d_radios_indices
self.d_flats_weights,
)
return radios
......
Markdown is supported
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