NeuTomPy toolbox, a Python package for tomographic data processing and reconstruction

In this article we present the NeuTomPy Toolbox, a new Python package for tomographic data processing and reconstruction. The toolbox includes pre-processing algorithms, artifacts removal and a wide range of iterative reconstruction methods as well as the Filtered Back Projection algorithm. The NeuTomPy toolbox was conceived primarily for neutron tomography datasets and developed to support the need of users and researchers to compare state-of-the-art reconstruction methods and choose the optimal data processing workflow for their data. In fact, in several cases sparse-view datasets are acquired to reduce scan time during a neutron tomography experiment. Hence, there is great interest in improving quality of the reconstructed images by means of iterative methods and advanced image-processing algorithms. The toolbox has a modular design, multi-threading capabilities and it supports Windows, Linux and Mac OS operating systems. The NeuTomPy toolbox is open source and it is released under the GNU General Public License v3, encouraging researchers and developers to contribute. In this paper we present an overview of the main toolbox functionalities and finally we show a typical usage example. © 2019TheAuthors.PublishedbyElsevierB.


Motivation and significance
Neutron Tomography (NT) has become a routine method at many neutron sources to non-destructively investigate the inner structure of a wide range of objects. The commercial software Octopus [1] by Inside Matters is a well established tool for reconstruction of tomographic data at neutron imaging beamlines. However, this software requires a significant investment and generally users can perform a preliminary data processing with Octopus only at the imaging facility. Data analysis is a crucial step for the output of an experiment, so users usually spend time to optimize the data processing mainly at home. This poses a strong demand of freeware and powerful tools to perform data processing of neutron data.
Image acquisition in NT is very time-consuming with respect to X-ray Computed Tomography (CT) and, in several cases, undersampled datasets are acquired to reduce the scan time and optimize beamtime usage during an experiment. The widely used Filtered Back Projection (FBP) algorithm generates reconstructed images affected by aliasing artifacts when the number of projections does not satisfy the Nyquist-Shannon condition [2]. Iterative reconstruction methods generally outperform analytical methods, such as FBP, to handle under-sampled datasets [3]. Octopus software provides only two reconstruction methods: the FBP and the Simultaneous Algebraic Reconstruction Technique (SART). Modern reconstruction methods are not implemented. On the other hand, several open source tools for tomographic reconstruction are available nowadays but they are mainly developed for X-ray CT and they are not ready to handle neutron data. Some image pre-processing algorithms are mandatory in NT to obtain accurate reconstruction, i.e. the estimation of the rotation axis tilt and the related registration of the projections, the suppression of gammaspots and the data normalization with respect to the radiation dose. Reconstruction tools for X-ray CT generally include some, but not all of such correction algorithms. For example, the ASTRA toolbox [4] is a Matlab and Python package that provides highly efficient implementation of iterative methods for CPUs and GPUs. ASTRA toolbox is only focused on the reconstruction step and it does not include any pre-processing, post-processing algorithms or functions to read and write data. On the other hand, the Python package TomoPy [5] includes several pre-processing and postprocessing algorithms and provides implementation for CPUs of a wide range of iterative reconstruction methods. Moreover TomoPy is not ready to handle neutron data, since it does not include functions to estimate the rotation axis tilt and to compute the related correction on projection data. Furthermore, TomoPy is available only for Linux and Mac OS operating systems. MuhRec [6] is the only free software that was conceived for NT. It includes several filters and pre-processing algorithms and it is currently the main free alternative to Octopus for data processing of neutron data. However, at time of writing, MuhRec does not provide any iterative reconstruction method support.
In this paper we present the NeuTomPy Toolbox, a new Python package for tomographic data processing, that is specifically designed to compensate the shortcomings of the aforementioned software tools. The NeuTomPy toolbox was conceived primarily for NT and developed to support the need of users and researchers to compare state-of-the-art reconstruction methods and choose the optimal data processing workflow for their data. The toolbox has a modular design, multi-threading capabilities and it supports Windows, Linux and Mac OS operating systems. The NeuTomPy toolbox is open source and it is released under the GNU General Public License v3, allowing users to freely use it and encouraging researchers and developers to contribute. Previously, this package has been used for comparative studies [3,7] of reconstruction methods in NT and now is freely distributed to the neutron imaging community.

Software description
Here we describe the architecture of NeuTomPy Toolbox and present its main functionalities.

Software architecture
The NeuTomPy toolbox is written in Python. We chose this programming language because it is open-source, cross-platform, human-readable and allows researchers to use and contribute to it easily. The toolbox is divided into several sub-modules, each of these represents a particular phase of a typical CT reconstruction pipeline. The entire chain is represented in Fig. 1. The NeuTomPy toolbox exploits several Python libraries for scientific computing and image processing, i.e. NumPy [8], NumExpr [9], SciPy [10], scikit-image [11], OpenCV [12] and SimpleITK [13]. In particular, the CT reconstruction step is powered by the ASTRA Toolbox. NeuTomPy combined with ITK-SNAP [14] or 3D Slicer [15] turns out to be a complete open-source software suite for CT.

Software functionalities and sample code snippets
The NeuTomPy toolbox allows to perform the steps of a typical CT reconstruction workflow (Fig. 1). The first task is represented by the reading of a raw dataset. The implemented reader handles TIFF and FITS files and converts a stack of images into a numpy array. A dataset containing raw projections, dark-field, flat-field images and the projection at 180 • can be read by: import neutompy as ntp proj ,dark ,flat , proj_180 = ntp. read_dataset ( proj_180 = True) hence the user can select the data to read from a dialog box. Subsequently, the projection data must be normalized with respect to dark-field and flat-field images to compute the transmission images. If the source intensity is not stable the images can be normalized with respect to the radiation dose [3]. In this case, the user must specify a region of interest (ROI) which corresponds to a background area not covered by the specimen in all the projections (we called it the dose ROI). It can be specified in three different ways: drawing interactively a rectangular selection, specifying the ROI' s coordinates or reading an ImageJ .roi file. For example, to normalize data and select interactively the dose ROI, the Python instruction is: where the function normalize_proj returns a 3D array containing the stack of normalized projections (norm) and a 2D array representing the normalized radiograph at 180 • (norm_180). A common experimental issue in NT is the misalignment of the rotation axis with respect to the vertical axis of the detector.
The function correction_COR evaluates the horizontal offset and the tilt angle by minimizing the squared error between two opposite radiographs computed at different vertical positions, as described in [6], and finally it registers all the projections. The Python instruction for this task is: norm = ntp. correction_COR (norm , proj_0 , proj_180 ) where proj_0 and proj_180 are the projections (raw or normalized) at 0 • and 180 • , respectively. The user selects interactively different ROIs where the sample is visible. Subsequently the results and some information about the evaluation of the rotation axis are shown. We report in Fig. 2 an example for the rotation axis correction: the difference of the projections at 0 • (P 0 ) and the mirrored projection at 180 • (P flipped π ) before and after the correction are shown in the left and right side, respectively.
The NeuTomPy toolbox includes an outlier removal which replaces a pixel value by the median of the neighborhood pixels if it deviates from the median by more than a certain value. This threshold value can be specified by the user as a global value or proportional to the local standard deviation. It is provided also a destriping filter, based on combined wavelet and Fourier analysis, to suppress the ring artifacts [16].
The reconstruction module includes all CPU-and GPU-based algorithms for 2D parallel beam geometry implemented in the ASTRA toolbox and some additional reconstruction methods distributed as ASTRA plugins. The available algorithms are summarized in Table 1. The instruction to perform a CT reconstruction is the following: rec = ntp. reconstruct (norm ,angles ,method , parameters ) where rec is the reconstructed volume, angles is onedimensional array containing the view angles in radians, method is a string which indicates the algorithm to use and parameters is a Python dictionary that contains specific settings of the Fig. 1. Diagram representing the typical CT data processing steps that can be performed by NeuTomPy toolbox. The package has a modular structure that follows the data processing chain. reconstruction algorithm. The allowed values for method and parameters follow the convention of the ASTRA toolbox, reported in the documentation [17]. For example, the following instruction is used to compute with GPU support a FBP reconstruction with the Hamming filter: rec = ntp. reconstruct (norm ,angles , method = " FBP_CUDA " , parameters ={ " FilterType " : " hamming " }) while a SIRT reconstruction with 100 iterations and pixel values limited in the range [0, 2] can be performed by: rec = ntp. reconstruct (norm ,angles , method = " SIRT_CUDA " , parameters ={ " iterations " :100 , " MinConstraint " :0.0 , " MaxConstraint " =2.0}) .
The NeuTomPy toolbox allows to compare and evaluate the performance of different reconstruction algorithms in terms of several image quality indexes. The metrics implemented are the Contrastto-Noise-Ratio (CNR) [3], the Normalized Root Mean Square Error (NRMSE) [3], an edge quality metric [3] and the Structural Similarity Index (SSIM) [18].

Illustrative examples
Here we demonstrate the possibility to perform several reconstruction algorithms and compare them quantitatively using the NeuTomPy toolbox. We used neutron images of a phantom sample Table 1 List of the CT reconstruction methods included in Neu-TomPy Toolbox for two-dimensional parallel-beam geometries.

Method
CPU GPU BP [2] x x FBP [2] x x ART [2] x SART [2] x x CGLS [19] x x SIRT [20] x x NN-FBP [21] x x MR-FBP [22] x x acquired at the IMAT beamline [23,24], ISIS neutron spallation source, UK. The phantom, already analyzed in a previous work [3], is an aluminium cylinder containing four holes of different diameters and filled with iron powder. We used for CT reconstruction an under-sampled dataset with 1/3 of the number of projections required by the Nyquist-Shannon condition. We performed FBP, SIRT and CGLS reconstructions and we compare them in terms of the image quality indexes NRMSE, SSIM and CNR. We consider the SIRT reconstruction (200 iterations) of a full-view dataset, which is sampled to fulfill the Nyquist-Shannon condition, as the reference image for the computation of the NRMSE and SSIM. The CNR was computed considering a ROI that includes one iron rod and with the second ROI outside the sample. The results are shown in Fig. 3. It is clear that the two iterative algorithms outperform the FBP method.
In fact, the CGLS and the SIRT reconstructions have higher CNR and SSIM, and lower NRMSE than the FBP, which indicate better image quality. In general, the under-sampling and the noise in the projection data cause in the reconstructed images a broadening of the attenuation coefficients distribution. However, unlike FBP reconstruction, the CGLS and the SIRT images are characterized by a bimodal distribution of the gray values, which reflects the composition of the sample. The source code of this analysis is omitted here for brevity. However, the source code for this and other examples can be found in the GitHub repository.

Impact
Data processing is the last step of a NT experiment but it is crucial for the interpretation of the results. Advanced image processing algorithms can extract hidden information from data and reduce the tomographic scan time. Hence new software tools, specifically designed for neutron data, are required to compare state-ofthe-art image processing algorithms. Working on robust methods and tools to improve image quality means get better output from NT experiments. However, state-of-the-art iterative reconstruction methods are not implemented in Octopus and MuhRec, which are the leading software for NT reconstruction. The NeuTomPy toolbox solves this shortcoming because it is ready to work with neutron data and allows to perform and compare several iterative reconstruction methods. Researches can define the optimal data processing workflow for their specific problem using the Neu-TomPy toolbox. The code is open-source, hence developers and researchers are invited to contribute.

Conclusions
In this paper we presented the NeuTomPy Toolbox, a new Python package for tomographic data processing. We demonstrated that the toolbox is ready to work with neutron data and allows researchers to state the optimal data processing workflow for their specific investigation. The first release includes preprocessing algorithms, artifacts removal and a wide range of classical and state-of-the-art reconstruction methods. The NeuTomPy toolbox supports Windows, Linux and Mac OS operating systems and it is released as open source. Researchers can freely use it and contribute to the project.
The future development will involve improvement of preprocessing algorithms (e.g. scattering correction), addition of new reconstruction methods and finally the implementation of a Graphical User Interface (GUI).