# 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](apidoc/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](apidoc/nabu.preproc.ccd) and [CudaFlatField](apidoc/nabu.preproc.ccd_cuda)


### 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](apidoc/nabu.preproc.ccd) and [CudaCCDCorrection](apidoc/nabu.preproc.ccd_cuda)


### 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](apidoc/nabu.preproc.double_flatfield) and [CudaDoubleFlatField](apidoc/nabu.preproc.double_flatfield_cuda)

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](apidoc/nabu.preproc.ccd) and [CudaLog](apidoc/nabu.preproc.ccd_cuda)

### 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](apidoc/nabu.misc.rotation) and [CudaRotation](apidoc/nabu.misc.rotation_cuda)

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](apidoc/nabu.reconstruction.sinogram)

### 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:
  - numpy-based using `pywavelets`: [MunchDeringer](apidoc/nabu.reconstruction.rings) and an implementation from Francesco Brun (see `nabu/thirdparty`).
  - cuda-based using [pycudwt](https://pypi.org/project/pycudwt): [CudaMunchDeringer](apidoc/nabu.reconstruction.rings_cuda)

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:
  - numpy-based, using the package [algotom](https://github.com/algotom/algotom): [VoDeringer](apidoc/nabu.reconstruction.rings)
  - cuda/cupy-based, using the package [tomocupy](https://github.com/tomography/tomocupy): [CudaVoDeringer](apidoc/nabu.reconstruction.rings_cuda)

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](apidoc/nabu.reconstruction.rings) and [CudaSinoMeanDeringer](apidoc/nabu.reconstruction.rings_cuda)

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

A 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"
```

Configuration file: `[preproc]`: `flat_distortion_params` and `flat_distortion_correction_enabled`

API: [estimate_flat_distortion](apidoc/nabu.estimation.distortion)


## Phase retrieval

Phase retrieval is still part of the "pre-processing", although it has a dedicated section in the [configuration file](nabu_config_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](phase.md)

### 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](apidoc/nabu.preproc.phase) and [CudaPaganinPhaseRetrieval](apidoc/nabu.preproc.phase_cuda)

### 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](apidoc/nabu.preproc.ctf)

### 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`
`start_x`, `end_x`, `start_y`, `end_y`, `start_z`, `end_z`.

### FBP

Configuration file: section `[reconstruction]`: `method = FBP`, `padding_type = edges`, `fbp_filter_type = ramlak`

This is the Filtered Back Projection reconstruction method.

### Reconstructor

This utility can only be used from the API ([Reconstructor](apidoc/nabu.reconstruction.reconstructor) and [CudaReconstructor](apidoc/nabu.reconstruction.reconstructor_cuda)).

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](apidoc/nabu.misc.histogram) and [VolumeHistogram](apidoc/nabu.misc.histogram)

## File formats

### HDF5

By default, the reconstructions are saved in HDF5 format, following the [Nexus NXTomo convention](https://manual.nexusformat.org/classes/applications/NXtomo.html). The output data is saved along with the software configuration needed to obtain it.

When a [multi-stage reconstruction](nabu_cli.md) 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)](https://docs.openmicroscopy.org/bio-formats/5.8.2/users/imagej/installing.html)

### 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](estimation.md).