Features overview

Nabu is a tomography processing software. The term “processing” stands for various things: pre-processing the radios/sinograms, reconstruction, and image/volume post-processing.

Pre-processing

Pre-processing means steps done before reconstruction.

API : nabu.preproc

Flat-field normalization

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. A linear interpolation is done between flats to estimate the “current” flat to use.

It’s also possible to normalize all frames with the synchrotron current during this step.

Configuration file: section [preproc]: flatfield_enabled = 1.

API: FlatField and CudaFlatField

CCD hotspots removal

Some pixels might have unusually high values. These are called “hotspots”. A simple an efficient way to remove them is to apply a thresholded median filter: a pixel is replaced by the median of its neighborhood if its values exceeds a certain relative threshold.

Configuration file: section [preproc]: ccd_filter_enabled = 1 and ccd_filter_threshold = 0.04.

API: CCDCorrection and CudaCCDCorrection

Double flat-field

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.

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.

Configuration file key: section [preproc]: double_flatfield_enabled = 1.

API: DoubleFlatField and CudaDoubleFlatField

You can pre-compute the double flatfield and provide it to the nabu configuration file, in order to avoid recomputing it each time. To this end, a CLI tool nabu-double-flatfield is available. Once you generated a double flat-field file, eg. dff.h5, you can provide processes_file = dff.h5 in [preproc].

Logarithm

When doing an end-to-end volume reconstruction, the -log() step has to be explicitly specified.

This step enables the conversion of the projections to the usual linear model setting, assuming a Beer–Lambert–Bouguer transmission law: I = I_0 exp(-mu(x, y)) to mu(x, y) = -log(I/I_0)

Numerical issues might occur when performing this step. For this reason, the projections values can be “clipped” to a certain range before applying the logarithm.

Configuration file: section [preproc]: take_logarithm = 1, log_min_clip = 1e-6, log_max_clip = 10.0.

API: CCDProcessing.Log and CudaLog

Projections rotation and detector tilt auto-detection

Each projection can be rotated by a user-defined angle. This can be useful to correct the effects of a detector tilt, or rotation axis tilt in some cases.
A command-line interface tool nabu-rotate is also available.

Configuration file: [preproc]: rotate_projections (value in degrees).

API: Rotation and CudaRotation

Nabu can also attempt at estimating the detector tilt, if any. If a tilt is detected, then the projection images will be rotated by the found angle value.

Configuration file: [preproc]: tilt_correction and autotilt_options.
Parameter tilt_correction can be the following:

  • A scalar value: tilt correction angle in degrees. It is equivalent to rotate_projections parameter.

  • 1d-correlation: auto-detect tilt with the 1D correlation method (fastest, but works best for small tilts)

  • fft-polar: auto-detect tilt with polar FFT method (slower, but works well on all ranges of tilts)

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

Sinogram-based rings artefacts removal

There are three deringer acting on sinogram in nabu.

Fourier-Wavelets

Nabu implements the following paper:

Beat Münch, Pavel Trtik, Federica Marone, and Marco Stampanoni, “Stripe and ring artifact removal with combined wavelet - Fourier filtering,” Opt. Express 17, 8567-8591 (2009)

The library has two implementations:

In the configuration file, the relevant keys are sino_rings_correction = munch and sino_rings_options.

The options are the following:

  • levels: the number of wavelets decomposition levels. A high number indicates that the stripes will be filtered far down the wavelets pyramid, which can yield artefacts.

  • sigma: standard deviation of the Gaussian filter applied on the (Fourier transform of) wavelets coefficients

  • padding, if not empty, should be a tuple of integers indicating how much to pad the images horizontally before performing the wavelets transform. Padding helps to avoid boundary artefacts in the edges of the image.

Example: sino_rings_options = sigma=1.0 ; levels=10 ; padding=(100, 100)

Algotom/tomocupy

Nabu offers the rings artefacts removal techniques from Nghia Vo described in

Nghia T. Vo, Robert C. Atwood, and Michael Drakopoulos, “Superior techniques for eliminating ring artifacts in X-ray micro-tomography,” Opt. Express 26, 28396-28412 (2018)

There are again two implementations:

In the configuration file, the relevant keys are sino_rings_correction = vo and sino_rings_options.

The options are the following: sino_rings_options = snr=3.0; la_size=51; sm_size=21; dim=1 - please refer to algotom remove_all_stripe() for their meaning.

Mean subtraction

A very simple deringer, yet often efficient, is to subtract a curve from each line of the sinogram. In the simplest case, this curve can be the mean of all projections. In short: sinogram = sinogram - sinogram.mean(axis=0). Usually it’s good to high-pass filter this curve to avoid low-frequency effects in the new singoram.

Nabu implements this principle in SinoMeanDeringer and CudaSinoMeanDeringer

In the configuration file, the relevant keys are sino_rings_correction = mean-subtraction and sino_rings_options.

The options are the following: sino_rings_options = filter_cutoff=(0, 30)

Flats distortion correction

For highly-structured wavefronts (phase contrast on nano-tomography beamlines), flats distortion estimation/correction is available. If activated, each radio is correlated with its corresponding flat, in order to determine and correct the flat distortion.

Warning

This is currently quite slow (many correlations), and only supported by the “full radios pipeline”. Also the implementation is CPU-only.

Configuration file: [preproc]: flat_distortion_params and flat_distortion_correction_enabled

API: estimate_flat_distortion

Phase retrieval

Phase retrieval is still part of the “pre-processing”, although it has a dedicated section in the configuration file.

This step enables to retrieve the image phase information from intensity. After that, reconstruction consists in obtaining the (deviation from unity of the real part of the) refractive index from the phase (instead of absorption coefficient from the attenuation).

In general, phase retrieval is a non-trivial problem. Nabu currently implements a simple phase retrieval method (single-distance Paganin method). More methods are planned to be integrated in the future.

See also: Phase retrieval

Paganin phase retrieval

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

API : PaganinPhaseRetrieval and CudaPaganinPhaseRetrieval

Contrast Transfer Function (CTF) phase retrieval

A single distance CTF retrieval method is also available.

Configuration file: [phase]: method = CTF, delta_beta, ctf_geometry, ctf_translations_file and ctf_advanced_params.

API: CTFPhaseRetrieval

Unsharp Mask

Although it is a general-purpose image processing utility, unsharp mask is often used after Paganin phase retrieval as a mean to sharpen the image (high pass filter).

Each radio is processed as UnsharpedImage = (1 + coeff)*Image - coeff * ConvolvedImage where ConvolvedImage is the result of a Gaussian blur applied on the image.

Configuration file: section [phase]: unsharp_coeff = 1., unsharp_sigma = 1., unsharp_method = gaussian.

Setting coeff to zero (default) disables unsharp masking.

Reconstruction

Tomographic reconstruction is the process of reconstructing the volume from projections/sinograms.

Configuration file: section [reconstruction]. Related keys: angles_file, angle_offset, rotation_axis_position, enable_halftomo, fbp_filter_type, padding_type, start_x, end_x, start_y, end_y, start_z, end_z.

FBP

FBP (Filtered Back-Projection) is the workhorse method for reconstructing (parallel) tomography data.

Configuration file: [reconstruction]: method = FBP.

Other relevant parameters: padding_type (= edges), fbp_filter_type (= ramlak).

From the API, nabu provides two implementations of FBP:

The backprojector has many features that are not commonly found in other libraries: ROI reconstruction, easy tuning of the center of rotation, custom padding type and filtering, native half-tomography support (without sinogram stitching).

HBP

HBP (Hierarchical Back-Projection) is an accelerated (but approximate) version of FBP. It works by doing FBP recursively on smaller images, and merging the results. On huge slices, it’s substantially faster (5-10X) than FBP.

Configuration file: [reconstruction]: method = HBP.

Other relevant parameters: hbp_reduction_steps (= 2), hbp_legs (= 4).

API: HierarchicalBackprojector

Cone

Nabu provides the FDK (Feldkamp-Davis-Kress) method for reconstructing cone-beam data. It uses the BP3D_CUDA backprojector of the astra toolbox. The filtering and pre-weighting are done by nabu.

Configuration file: [reconstruction]: method = cone.

Other relevant parameters: source_sample_dist, sample_detector_dist (if not found in the data file ; will overwrite them if found).

From the API: ConebeamReconstructor has an API that offers most of the features of FBP implementations above (ROI reconstruction, custom filtering/padding, horizontal corrections, …)

MLEM

Maximum Likelihood Expectation Maximization is commonly used in emission tomography, or generally when the noise model follows a Poisson distribution (rather than Gaussian for limit case of many-photons in X-ray tomography).

For now nabu offers an implementation from corrct, aiming at providing advanced corrections (eg. absorption).

Configuration file: [reconstruction]: method = mlem ; implementation = corrct.

API: MLEMReconstructor

Reconstructor

This utility can only be used from the API (Reconstructor and CudaReconstructor).

The purpose of this class is to quickly reconstruct slices over the three main axes, possibly in a Region Of Interest. For example: “I want to reconstruct 50 vertical slices along the Y axis”, or “I want to reconstruct 10 vertical slices along the X axis, with each slice in the subregion (1000-1200, 2000-2200)”.

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 and VolumeHistogram

File formats

HDF5

By default, the reconstructions are saved in HDF5 format, following the Nexus NXTomo convention. The output data is saved along with the software configuration needed to obtain it.

When a multi-stage reconstruction is performed, the volume is reconstructed by chunks. Each chunk of reconstructed slices is saved to a HDF5 file. Then, a “HDF5 master file” is created to stitch together the reconstructed sub-volumes.

Tiff

Reconstruction can be saved as tiff files by specifying file_format = tiff in the configuration [output] section. Mind that tiff support currently has the following limitations:

  • Data is saved as float32 data type, no normalization

  • No metadata (configuration used to obtain the reconstruction, date, version,…)

It’s possible to use one single big file for the whole volume (instead of the default one-file-per-slice) using tiff_single_file = 1 in the configuration file. Mind that the generated tiff file cannot be natively read by all software - for example imagej needs an extension (bio-formats)

EDF

Reconstruction can be saved as tiff files by specifying file_format = edf in the configuration [output] section.

Jpeg2000

Reconstruction can be saved as jpeg2000 files by specifying file_format = jpeg2000 in the configuration [output] section. 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.

  • No metadata (configuration used to obtain the reconstruction, date, version,…)

The following can be configured through the configuration file:

  • jpeg2000_compression_ratio: Compression ratio for Jpeg2000 output

  • float_clip_values: Lower and upper bounds to use when converting from float32 to int. Floating point values are clipped to these (min, max) values before being cast to integer.

Note

Jpeg2000 needs the openjpeg library (libopenjp2 in Debian-like distributions, openjpeg=2.4.0 in conda) and the glymur python package. Multi-threaded encoding needs glymur >= 0.9.3 and openjpeg >= 2.4.0.

Save/resume from processing steps

It is possible to save arbitrary processing steps, and to resume the processing pipeline from a saved step. It speeds up the process of “reconstructing - editing configuration - reconstructing - etc”.
In the section [pipeline], there are two corresponding options save_steps and resume_from_step.
For example, to set a check point to the “sinogram” generation, set save_steps = sinogram.

Both save_steps and resume_from_step can be specified at the same time. If the file is not found, then all the processing is done in normal mode, and the steps are saved to the file. Otherwise, the pipeline resumes from the step. The file where data is dumped can be specified with the parameter steps_file.

Parameters estimation

Nabu provides some parameter estimation tools. Please see the dedicated page.