1. Introduction

A major drawback of a parallel-beam tomography system is that the reconstructable area [1] is limited by the beam size of the X-ray source and the size of the detector, which together determine the field of view (FOV). In practice, it is not always possible to prepare samples to fit to the FOV such as valuable fossils [2,3] or cultural heritage specimens [4]. In such cases, different scanning schemes to extend the FOV are needed. There are three basic data acquisition schemes which can be used: a grid scan [5,6], a 360-degree scan with offset center-of-rotation (COR), known as a half-acquisition scan [7,8], and a helical scan [9,10]. Each of these scanning schemes requires different pre-processing methods before the data can be reconstructed using the traditional approach, which reconstructs a slice from a 180-degree sinogram. In a grid scan and half-acquisition scan, image stitching methods are required. In a helical scan, data need to be converted to an equivalent circular scan dataset. With high-flux X-ray sources at synchrotron facilities, acquiring a tomographic dataset can be done in minutes, seconds or even faster [11] making in-situ and time-resolved experiments possible. Scientists therefore need to reconstruct a few slices quickly for assessment and decision making about the next step of their experiment. In some cases, users would like to reconstruct slices from a whole sample to locate the region-of-interest (ROI) for subsequent, higher resolution scans. This requires a processing pipeline that is easy-to-use, reliable, and automated.

Although the aforementioned data acquisition schemes are easy to implement, not many applications employing these techniques at synchrotron facilities can be found in the literature [12,13]. One of the main reasons is that there is a lack of data processing methods and open-source software to process such data. Some open-source software [14,15] can process a half-acquisition scan but users must input the COR manually based on visual inspection. Tomosaic software [16] can be used for a grid scan but is currently limited to processing of 180-degree scans. Additionally, metadata about spatial positions of the scans are required for auto-stitching methods to work. For helical parallel-beam scans, we found no Python open-source software to reconstruct this type of data.

The scope of this work is to provide the research community with data processing methods and open-source software for advanced scanning schemes, which are: easy-to-use, automatable, and capable of the highest possible reconstruction quality. Specifically, we introduce:

Considering that high-speed experiments and large datasets are possible using synchrotron-based tomography systems, the proposed methods and methods implemented in Algotom have been developed for working in the sinogram space to reduce the computational cost and to parallelize the computational work. To process sinograms independently, the information recorded in each sinogram must be independent. This requires that a parallel-beam tomography system must be aligned, i.e the rotation axis of a sample stage is parallel to the plane of the detector and perpendicular to the rows of the pixel-grid detector [18]. Alignment is a standard routine at synchrotron facilities [19,20]. If the independence condition cannot be met due to image rotation or radial lens distortion of a detecting system, Algotom provides methods for extracting tilted sinograms and unwarped sinograms with low computational cost.

Details of the methods are presented in Section 2. It is divided into two sub-sections: a section about methods used for processing a half-acquisition scan, a grid scan, and a helical scan; and a section about methods in the data processing pipeline implemented in Algotom. Section 3 show results of applying these methods to process two datasets: one acquired by using a grid scan with offset COR; and one acquired by using a helical scan with offset COR. Data used for this paper were collected at Beamline I12 and Beamline I13-2 of DLS [21].

2. Methods

2.1 Pre-processing techniques for data acquisition schemes extending the FOV

2.1.1 Methods for a 360-degree scan with offset center-of-rotation

2.1.1.1 Determination of the overlap-side and the overlap-area between sinograms

In a half-acquisition scan, a common routine is that a 360-degree sinogram is converted to an equivalent 180-degree sinogram (Fig. 1(a) and 1(c)), by stitching two halves of the 360-degree sinogram, before reconstruction. Image stitching is a well-established research area in computational photography. There are many choices of algorithms where feature-based approaches have been proven to be fast, robust, and fully automated [22–24]. However, applying these methods to sinogram images of X-ray tomography is not straightforward due to three major challenges. Firstly, the feature space of sinogram images is often not rich enough to find similar features as can be seen in Fig. 1(b) (the blue frames). Secondly, overlap areas are often chosen to be small to maximize the FOV, e.g., 10% to 20% of the area of a sinogram. This strongly affects the reliability of both feature-based and Fourier transform-based methods [23]. Lastly, there are regions of free space beyond the object (Fig. 1(d)), the features in the beam profile and optical artifacts are often very similar and thus highly correlated (Fig. 2(b)). They can cause incorrectly determined direction of stitching of the sinograms by automatic correlation-based methods (Fig. 1(e)).

 

Fig. 1. Demonstration of the routine of stitching sinograms. (a) 360-degree sinogram of a half-acquisition scan. (b) Two halves of the 360-degree sinogram used for stitching where the second half is flipped left-to-right. Blue frames indicate the selection areas. In the first search, the region at the outermost left side of the 2^nd sinogram is fixed while the same size region is shifted across the 1^st sinogram. (c) Correctly stitched 180-sinogram. (d) In the second search, the region at the outermost right side of the 2^nd sinogram is fixed. Correlation between regions in the solid red frames (Fig. 2(b)) is higher than the correlation between the solid blue frames (Fig. 2(a)) in (b). (e) Wrongly stitched sinogram if the determination is based on correlation correlation only.

Download Full Size | PDF

 

Fig. 2. Demonstration the routine of determining the overlap-side between sinograms. (a) Graph of correlation metrics corresponding to Fig. 1(b). (b) Graph of correlation metrics corresponding to Fig. 1(d) where the minimum metric is smaller than the one in (a). (c) Graph of points around the minimum metric point, the red dot in (a). (d) Graph of points around the minimum-metric point, the blue dot in (b).

Download Full Size | PDF

Here we use a correlation-based method with a technique proposed to solve the problem of wrong-side image-stitching caused by free space. A routine for determining the overlap-side and overlap-area between sinograms is as follows, where the first sinogram is used as the reference; the overlap-side, left or right, is referred to this sinogram:

The routine described above performs 1D searching, which has low computational cost, across the full width of a sinogram. This enables automated processing. The reason for using a quadratic coefficient to determine the overlap-side is that it represents the sharpness of the true overlap correlation (Fig. 2(c)); a small mismatch in the overlap destroys the correlation whereas accidental strong correlation due to low contrast features within space region persists over a broader range (Fig. 2(d)). The robustness of the algorithm can be improved by applying some processing steps such as denoising using a smoothing filter and normalizing intensities between two regions before calculating the correlation metric. Details of these implementations are in Code 1 (Ref. [17]).

For automation to work there must be a certain level of overlap between recording areas which can be coarsely adjusted, given the FOV and the sample stage position, before data are acquired. In our experience, an overlap of 10% to 20% of an image width of 2k pixels works for most cases. In the example [25] shown in Fig. 1 and Fig. 2, the width of the sinogram is 2.8k pixels. The calculated overlap-area (or overlap-width) is 418 pixels. The width of the search window is 100 pixels but insensitive. A larger window size increases the reliability of the result. However, it must be smaller than the overlap-area which can be estimated from data acquisition settings.

Here the correlation metric used is the inversion of the Pearson coefficient [26], a common metric to measure linear correlation between two images. Other metrics available in Python open-source software such as Spearman, Kendall’s tau, or point-biserial correlation coefficient [27], can be used. However, we found that the Pearson coefficient is reliable enough for our data. The inversion means that the minimum metric corresponds to the highest correlation. This is done to be consistent with other metrics used by our software (Ref. [17]) where the minimum values are the best values. The metric is calculated as

If the overlap-side is the left side of the 1^st sinogram,

If the overlap-side is the right side of the 1^st sinogram,

From the determined overlap-side and the overlap-area, a linear ramp to blend and stitch sinograms is used. Details of these implementations are shown in Code 1 (Ref. [17]). Although the above routine is demonstrated on sinograms, it certainly can be used on projection images. However, this may reduce the reliability of the method because projection images may exhibit repeated patterns, e.g., projections of a rectangular grid.

2.1.1.2 Reconstruction without converting to a 180-degree sinogram

A common approach to processing data from a half-acquisition scan is that two halves of a 360-degree sinogram are stitched to form a 180-degree sinogram. Then, it can be reconstructed using a standard reconstruction method. This approach works for most of the cases. However, stitching a sinogram suffering from intensity-fluctuation problems can result in sharp changes of intensities in the overlap-area which may be enhanced by the ramp-filter in the FBP (filtered back-projection) reconstruction method [28]. This can give rise to ring artifacts or streak artifacts. Here we introduce a very simple approach where there is no need to use image stitching or conversion of a 360-degree sinogram to a 180-degree sinogram. It needs to be used with back-projection-based reconstruction methods such as the FBP method or the BPF (back-projection and filtering) method [29].

The pre-processing steps apply linear ramp weighting row-by-row across the overlap area of a 360-degree sinogram (Fig. 3(a) and 3(b)), where the negative of the logarithm of the intensities is used; then zero-padding is applied to the result to expand the width of the sinogram to size of (2×W −OA) (Fig. 3(c)). Details of these implementations are in Code 1 (Ref. [17]). Back-projection in the full range of 360 degrees is applied to the expanded sinogram to reconstruct a slice as can be seen in Fig. 4.

 

Fig. 3. Pre-processing steps applied to the 360-degree sinogram for reconstruction without prior conversion to a 180-degree sinogram. (a) Flat-field-corrected and logarithm-transformed sinogram. (b) Linear-ramp-weighted overlap-area (arrows) of sinogram (a). (c) Zero-padded result of sinogram (b).

Download Full Size | PDF

 

Fig. 4. Demonstration of the back-projection process from the sinogram in Fig. 3(c). (a) Back-projection using the angle of [0°; 90°]. (b) Back-projection using the angle of [0°; 180°]. (c) Back-projection using the angle of [0°; 270°]. (d) Back-projection using the full range of 360 degree resulting in a slice.

Download Full Size | PDF

2.1.2 Methods for grid scans

2.1.2.1 Combination of grid scans and half-acquisition scans with image stitching applied in the sinogram domain

Using a grid scanning technique to scan a sample larger than the FOV of a detecting system is an easy implementation from the data acquisition perspective. However, to reconstruct the whole sample efficiently in term of storage usage and computational cost is challenging. Data processing steps commonly used in X-ray micro-tomography are flat-field correction, zinger removal, ring artefact removal, distortion correction, denoising, and tomographic reconstruction [30]. They need to be combined properly with an image stitching step to produce a full reconstructed sample.

In most synchrotron-based tomography systems, the rotation speed of a sample stage is much higher than the translation speed, so a common data collection routine is to perform a full tomographic fly-scan at each translation position instead of doing translation scanning at each angle position. Different grid-scanning schemes [6,13,16] can be used depending on how the rotation axis is positioned with respect to the FOV and a sample. This dictates how a data processing pipeline should be constructed.

In the first approach, sub-volumes of a sample are scanned and reconstructed independently then the stitching is performed on the reconstruction space. Scanning each sub-volume can be done using a traditional 180-degree scan with the rotation axis at the middle of the FOV or with a half-acquisition scan [13]. A disadvantage of this approach is that the intermediate reconstruction results of each scan need to be stored for later use in the stitching step. More importantly, stitching reconstructed images is not straightforward. Reconstructing from truncated projections introduces artifacts [31–33]. Additionally, the low-frequency components and the DC-component of slices between overlapped regions may vary depending on the contribution of sample parts outside the FOV of these slices. This affects the quality of the blended images.

In a synchrotron-based tomography system, the detector position is often fixed relative to the source. In the previous approach, each sub-volume of the sample is translated in 3D to the reconstructable volume, while the position of the rotation axis is unchanged with respect to the FOV. In the second approach, the rotation axis is translated horizontally in 1D relative to the FOV and the sample is translated vertically in 1D to scan at different heights (Fig. 5(a)). Thus, the translation dimension is 1D less than the first approach. To process data acquired by this type of scanning, stitching is performed in projection space or sinogram space before reconstruction. Data can be acquired in the range of [0; 180°] with the rotation axis positioned at the horizontal center of the grid as shown in Fig. 5(a) and used in [16]. Here we propose a modification of this approach by acquiring data in the range of [0;360°] (Fig. 5(b)) to reconstruct the same volume but with a half number of translations. This scanning technique is simply derived from a half-acquisition scan where the rotation axis is positioned with respect to the FOV of the grid as a whole.

 

Fig. 5. Demonstration of the grid-scanning schemes where stitching is performed in the projection space or the sinogram space. The rectangular areas represent the FOV of a detecting system and the red line indicates the position of the rotation axis. (a) Rotation axis is at the center of the grid where data is collected in the angle range of [0°; 180°]. (b) Position of the rotation axis is the same as in (a) but data is collected in the angle range of [0°; 360°]. This approach reduces a half of the number of stage translations compared to the approach in (a).

Download Full Size | PDF

To avoid storing intermediate results, i.e., stitched data, the data processing pipeline is conducted in the sinogram domain. If there is radial lens distortion [34] it needs to be corrected prior to the auto-stitching method or the results will not be reliable. We developed a correction routine, accompanied by Python software [35,36], where a reference pattern is imaged to characterize the distortion produced by the detector system. Our approach is to extract information from this pattern to calculate parameters of a distortion correction model. From the parameters determined, an undistorted sinogram is generated using a few adjacent rows of projection images (Code 1, Ref. [17]).

A demonstration of how to reconstruct a slice of a full sample acquired using our approach is shown in Fig. 6. 360-degree sinograms of each tomograph (Fig. 6(a) and 6(b)) are stitched to form a single 360-degree sinogram. This sinogram is converted to a 180-degree sinogram (Section 2.1.1.1), then the FBP reconstruction method is applied resulting in the slice in Fig. 6(e). It can also be reconstructed directly from the sinogram in Fig. 6(c) using the method described in Section 2.1.1.2. Python implementations of these steps are in Code 1 (Ref. [17]). A practical advantage of the proposed scanning technique is that the pre-processing steps of COR determination and image stitching can all use the same method described in Section 2.1.1.1. This method allows search of the whole width of a sinogram with a low computational cost. It thus enables a fully automated data processing pipeline. Details of the data processing pipeline to reconstruct a whole sample using our proposed data acquisition method are presented in Section 3.

 

Fig. 6. Processing steps to reconstruct a slice from data acquired using a combined grid scan and half-acquisition scan. (a) 360-degree sinogram from a tomograph of a grid scan. (b) 360-degree sinogram from a neighboring tomograph of the tomograph in (a). (c) 360-degree sinogram stitched from sinograms (a) and (b). (d) 180-degree sinogram converted from sinogram (c). (e) Reconstructed image of sinogram (d).

Download Full Size | PDF

2.1.2.2 Sample detection

Depending on the shape of a sample there are parts of a grid scan where no portions of the sample are present (e.g., the 2^nd row of the 1^st column of the grid in Fig. 5(a)). For automating a processing pipeline, one needs to detect these areas and avoid using them for determining the overlap-side and overlap-area between cells of the grid. A method to indicate if there is a sample in an image is needed. It can happen that a part of a sample may stay inside the projection images at certain angles but outside at other angles, so using sinogram images for detection is more reliable. Here we propose a simple detection technique by exploiting an interesting feature of a sinogram is that its 2D Fourier transform exhibits a double-wedge region of significant coefficients (Fig. 7(a) and 7(c)). This feature shows not only in a 360-degree sinogram but also in a 180-degree sinogram or a sinogram with the COR out of the FOV. A mathematical derivation and thorough discussion about this feature are fully presented in Ref. [37].

 

Fig. 7. Demonstration of the sample detection algorithm. (a) Image with a part of a sample inside the FOV. (b) Image with no part of the sample inside the FOV. (c) Logarithmic scale of the 2D Fourier transform of image (a). (d) Logarithmic scale of the 2D Fourier transform of image (b). (e) Smoothed image of the top-left quarter of the image (c). (f) Smoothed image of the top-left quarter of image (d) shown no distinct boundary compared to (e).

Download Full Size | PDF

Knowing that significant Fourier coefficients of a free space image do not exhibit the same feature as a sinogram (Fig. 7(b) and 7(d)), our technique aims to detect the boundary of the double-wedge region. If such a boundary exists, the sinogram contains a sample. Due to the symmetry of the Fourier space, using a quarter of an image is enough for detecting the boundary. The detection algorithm works as follows:

Details of these implementations are in the Algotom package, Code 1 (Ref. [17]).

2.1.3 Methods for a helical scan

2.1.3.1 Calibration and data acquisition

In a helical scan, a sample is rotated and translated along the rotation axis during the scanning process. To process data acquired using the helical scan, a proposed approach is to convert the data to an equivalent circular scan data before reconstruction. It is important here to determine the pixel size precisely with respect to the translation increment; taking into account that the X-ray beam is slightly divergent in a synchrotron-based tomography system, and the mechanical error of the translation stage needs to be smaller than the pixel size. To do that, we scan a dense sphere as shown in Fig. 8. The sphere is segmented and its center of mass (COM) is calculated. Using the positions of the COM of the sphere (Fig. 8(c)), referred to the coordinate of an image, in correlation to the translation positions of the stage we can determine the pixel size and evaluate the accuracy of the translation stage.

 

Fig. 8. Point-like object for characterizing a translation stage used for a helical scan. (a). Chromium sphere. (b) Overlay of a small number of projections of the sphere during a helical scan. (c) Trajectory of the COM of the sphere.

Download Full Size | PDF

The routine for collecting helical data is demonstrated in Fig. 9 and described as follows:

 

Fig. 9. Demonstration of the data-acquisition routine for the helical scan. (a) Translating the stage to see the top of the sample, recording y_s. (b) Translating the stage to see the bottom of the sample, recording y_e. (c) Moving the stage to the starting point of (y_s - ½P_t). Part of the image above the top of the sample has no use. (d) Moving the stage to the stopping point of (y_e + ½P_t). Part of the image below the bottom of the sample has no use.

Download Full Size | PDF

Suppose that one would like to scan a sample with a length of (y_e−y_s) where y_e and y_s are determined by moving the translation stage to see the bottom and the top of the sample, respectively. To be able to reconstruct the top and the bottom of the sample, it needs to be inside the FOV at least in the angle range of [0, 180°]. This is done by shifting the translation stage to the position of

2.1.3.2 Converting from helical acquisition data to equivalent circular acquisition data

Figure 10 demonstrates the process of generating a circular sinogram from helical scan data. As the sample is translated, to extract 1D-projections of the same slice of the sample one must extract data at different rows of the collected images. Details of the conversion routine are as follows.

 

Fig. 10. (a) 1D-projections of the same slice of the sample located at different rows across images. (b) Incorrect sinogram generated by extracting the same row (along the blue line in (a)) across images. (c) Correct sinogram generated by extracting rows along a tilted plane (indicated by the red lines in (a)).

Download Full Size | PDF

Suppose that one needs to generate a circular sinogram at the height of

Knowing that at each movement the stage is translated by

In this projection, the row position of the slice y_k is given by

If j₀ is not an integer, the 1D-projection is extracted by interpolation between the nearest rows. The routine is repeated for all projections in the range of [i₀; i₀ + N₁₈₀] resulting in the circular sinogram as shown in Fig. 10(c). Note that the corresponding angle range must be used for reconstruction instead of resetting to the range of [0°; 180°]. Otherwise, the reconstructed images are rotated between slices. In general, one can extract all projections of the slice y_k which is in the range of [i₀; i₀ +FOV/y_step] and use them all for reconstruction. This helps to improve the signal-to-noise ratio of the reconstructed image. Python implementations of the mentioned steps are in Code 1 (Ref. [17]).

Equation (12) shows that the precise determination of the pixel size, Δx, in relation to the stage translation is crucial to the quality of the extracted sinogram. If the pixel size cannot be determined using the above calibration routine, an alternative way is to calculate the total number of vertically shifted pixels between a projection at its counterpart at 180° or 360°. This requires that the pitch needs to be smaller than the FOV, to be able to record the same part of the sample after 180° or 360° rotation. If this approach cannot yield a good result, one can generate multiple sinograms using estimated pixel sizes and measure the quality metrics of reconstructed images [39] to find the best value.

2.1.3.3 Ring artifact removal and blob removal

A significant advantage of using a helical scan is that it helps to clean ring artifacts caused by partially defective pixels and scintillator defects of a detecting system (Fig. 11) [12,38]. It can also remove low-frequency ring artifacts caused by the flat-field correction from an X-ray source having an unstable beam profile, or from a strongly absorbing sample recorded by a low dynamic range camera. However, in some cases where flat-field images are noisy there are still some ring artifacts around the middle of the reconstructed image (Fig. 11(c)). These ring artifacts are caused by the flat-field correction process and can be cleaned (Fig. 11(d)) using the sorting-based method, algorithm 3, in [40].

 

Fig. 11. Demonstration of the advantage of the helical scan in removing ring artifacts. (a) Reconstructed image using a circular scan. (b) Reconstructed image using the helical scan. (c) Magnified view from the red box in (b). (d) Same area as in (c) after the ring-artifact removal.

Download Full Size | PDF

Another problem is that unresponsive pixels can give rise to streak artifacts in a reconstructed image as can be seen in Fig. 12. To remove these artifacts we adapt algorithm 4 and algorithm 6 as used in [40] for removing ring artifacts caused by unresponsive pixels in a circular scan. First, a binary blob map (Fig. 13(a)) is generated by applying algorithm 4 to each row of a flat-field image (Fig. 12(a)) where the intensity profile of each row is normalized by dividing with the median-filtered version of itself. Then, the intensities inside the blobs are interpolated from neighboring pixels using the blob mask in combination with Eq. (12). Details of these implementations are in Code 1 (Ref. [17]). As can be seen in Fig. 13(b) and 13(c), all the streak artifacts are removed.

 

Fig. 12. (a) Defective pixels (arrowed) visible in the flat field image. (b) Impact of the defective pixels on the sinogram extracted. Note that the locations of the affected areas (arrowed) depend on the pitch used and Eq. (12). (c) Streak artifacts (arrowed) in the reconstructed image caused by the defective areas in (b).

Download Full Size | PDF

 

Fig. 13. Demonstration of the method used for removing streak artifacts in the helical scan caused by the unresponsive pixels. (a) Binary blob map generated from the flat-field image in Fig. 12(a). (b) Corrected sinogram. (c) Reconstructed image from sinogram (b).

Download Full Size | PDF

2.2 Data processing pipeline and the Algotom package

In the tomographic reconstruction process, a raw dataset goes through a series of data processing steps [30], commonly called a data processing pipeline, before the final result can be used for post analysis. The pipeline is designed to satisfy two basic requirements: the reconstructed images show the least artifacts and the measured refractive indices are closest to those calculated. As shown in Fig. 14, the reconstruction process is divided into three stages: pre-processing, reconstruction, and post-processing. In each stage, many data processing methods can be used where each of them may work in different spaces: the projection space, the sinogram space, or the reconstruction space. In the pre-processing stage, there are different ways of building a pipeline to optimize the computational cost and the quality of the reconstructed images. It is important that pre-processing methods are chosen to work in the same space, because switching between the projection space and the sinogram space requires intermediate storage of data.

 

Fig. 14. Three stages of the tomographic reconstruction process where many data processing methods can be added to the pipeline depending on the set-up and the hardware issues of a tomography system.

Download Full Size | PDF

Algotom development was started at the I12-JEEP beamline in 2014 as Python codes to process data acquired by the beamline’s large FOV detector, which uses two imaging sensors to cover a rectangular FOV [11]. Images from these cameras must be stitched, using the same approach as shown in Section 2.1.1.1, before tomographic reconstruction can take place. Data processing methods for improving the quality of tomographic data; removing artifacts caused by imperfections of hardware components; making use the beamline capabilities; processing data acquired by non-traditional scanning techniques; and automating data processing pipeline have been actively developed at I12 over the years. These methods have been used internally by I12’s users and refined further for publication and sharing with the research community through open-source software [36,41]. Specifically, methods for automated determination of the center-of-rotation [18], distortion correction [34], and removing ring artifacts [40] have been contributed to Tomopy [14] and Savu [15], which are software frameworks for high-throughput processing of tomographic data.

In contrast to Savu and Tomopy, which are optimized for speed, Algotom is a package of data processing algorithms and tools which are designed to be easy-to-use and easy-to-deploy prototype methods. The development of Algotom has focused on pre-processing methods which work in the sinogram space to reduce computational cost. Methods working in the projection space such as phase filter, distortion correction, or rotation correction have been adapted to work in the sinogram space. Figure 14 shows the variety of data processing methods which can be added to the pipeline depending on the set-up and the hardware issues of a tomography system. The following sections present data processing methods available in the Algotom package and commonly used at the I12 beamline. The software is the implements of the methods presented in the previous section, methods developed by us earlier [18,34,40–42], methods published in the literature [28,43], and some of our methods not published elsewhere. It provides users with a complete set of tools to process data acquired by the data acquisition schemes for extending the FOV, aiming for the highest quality possible.

2.2.1 Pre-processing methods

2.2.1.1 MTF (modulation transfer function) deconvolution

The loss of resolution caused by the scattering of scintillation photons in a scintillator-based detector system is reduced by determining the MTF of the scintillator using a calibration image and following the routine described in [42]. Then, raw images are deconvolved by the division in the Fourier domain.

2.2.1.2 Flat-field correction

The intensity ratio of the image to the incoming beam at each pixel is found simply by the division of an image-with-sample by an image-without-sample, i.e., flat-field image, where a dark-field image (camera’s dark noise) is subtracted from each image before the division.

2.2.1.3 Zinger removal

Zingers are prominent bright dots in projection images caused by scattered X-rays hitting the detector system CCD or CMOS chip (Fig. 15(a) and 15(b)). They produce line artifacts across a reconstructed image (Fig. 15(c)). There are some ways to remove them such as normalization-based techniques [44,45] or statistical methods [15,30]. They may work better in the projection space than the sinogram space or vice versa. Here we propose a unified method which can work in both spaces. It can be used on raw images or flat-field-corrected images. Zingers in flat-field image, if not removed, cause stripe artifacts in a sinogram and ring artifacts in a reconstructed image due to the flat-field correction process. These artifacts can be removed using the corresponding removal methods [40].

 

Fig. 15. Artifacts caused by zingers. (a) Zingers in the sinogram space. (b) Zingers in the projection space. (c) Line artifacts caused by the zingers.

Download Full Size | PDF

Our approach is a normalization-based method. However, in contrast to other techniques [14,45], which use a smoothing filter for producing a smoothed version of the original image, we use the average of the images shifted around the original position but not including the original image (Fig. 16(a) and 16(b)). The shifting step is the same size as a zinger. In the normalized image, obtained by dividing the original image by the average of the shifted versions, the zingers are located using a threshold parameter (Fig. 16(c)). Then, their values in the original image are replaced by the values in the averaged image, only at these locations with size equal to or smaller than the given size (Fig. 16(d)). A major advantage of our approach is that it does not significantly alter the image, even when used to remove low-contrast zingers.

 

Fig. 16. Demonstration of the zinger removal method. (a) Original image. (b) Average of the shifted versions of image (a) but excluding image (a). (c) Normalized image. (d) Image with the zingers removed.

Download Full Size | PDF

2.2.1.4 Ring artifact removal

Ring artifacts may come from the irregular response of a detecting system caused by the non-linear response of the sensor chip, dust on the optics components, and defects in the scintillator. Our work on classifying types of ring artifacts and removing them efficiently were published in [40] and available in open source software such as Sarepy [41], Tomopy [14], and Savu [15]. In Algotom, including with these methods we also provide supplementary tools for customizing ring removal methods (Code 1, Ref. [17]).

Another source of ring artifacts is introduced by the flat-field correction process. It is caused by features in a flat-field image which are uncorrelated to projection images. In Fig. 17, there are features (arrowed) in the flat-field image caused by the X-ray optics [20]. During tomographic scanning, these features change over time, so that they cannot be removed using the flat-field correction technique (Fig. 17(c)). As a result, there are very large ring artifacts (Fig. 17(d)) which are difficult to remove without altering the image. In such cases, using a helical scan is a good solution.

 

Fig. 17. Ring artifacts caused by the change of the flat-field image. (a) Flat-field image. (b) Image with the sample. (c) Flat-field corrected image. (d) Ring artifacts (arrowed) in the reconstructed image.

Download Full Size | PDF

2.2.1.5 Distortion correction

The distortion correction routine is used to correct radial lens distortion of a lens-coupled detector system. Distortion parameters are calculated from a reference pattern using our Vounwarp software [35,36]. From the parameters determined, projection images are undistorted by a backward mapping process. Directly generating an undistorted sinogram using a few adjacent rows of projection images is also available in Algotom.

2.2.1.6 Tilted sinogram generation

The tilt correction routine is used to correct the problem of projection-image rotation but working in sinogram space using a few adjacent sinograms.

2.2.1.7 Center-of-rotation determination

A method developed by us and published in [18] is used to calculate the center-of-rotation from a 180-sinogram. A method for a half-acquisition sinogram presented in Section 2.1.1 is also implemented in Algotom.

2.2.1.8 Double-wedge filter

The previously described double-wedge property of the 2D Fourier transform of a sinogram (Section 2.1.2.2) [37] can used to remove sample parts outside the FOV of a sinogram. In fast tomography experiments [11,46], the number of projections collected is often insufficient for direct reconstruction methods to work well [1]. As a result, iterative reconstruction methods, mainly linear algebra-based, are preferred. Because these methods construct and solve a set of linear equations from a sinogram, they require that the sums of intensities of each sinogram-row are the same. This condition is not satisfied in the case of a sample [46] larger than the FOV as can be seen in Fig. 18(a). It shows a sinogram of the sample where there are parts coming in and out the field of few of a single row of a detector. Applying iterative reconstruction methods [47–50] on this sinogram results in various types of artifacts as shown in Fig. 18(b)–18(d). The slice from the FBP method (Fig. 18(e)) is not affected by the problem. However, as the number of projections is insufficient, the slice is noisy.

 

Fig. 18. Artifacts caused by an extended object to iterative reconstruction methods. (a) Sinogram where there are parts come in and out the FOV (arrows). (b) Reconstructed image from the SART method with 500 iterations (a high iteration number was used due to the slow convergence of the SART method caused by the extended object). (c) Result of the CGLS method with 100 iterations. (d) Result of the SIRT method with 100 iterations. (e) Result of the FBP method.

Download Full Size | PDF

We can reduce the effect of extended objects on the performance of linear algebra-based methods by removing parts of a sample larger than the width of its sinogram. This can be done by masking the 2D Fourier transform of this sinogram. The double-wedge region in the Fourier space is symmetric and well-defined if the sinogram used is 360-degree (Fig. 19(a)). Because a sinogram is often collected in a 180-degree range in parallel-beam tomography, forming a 360-degree sinogram from this sinogram can be done [18] by: flipping the 180-degree sinogram from left to right; shifting this image to the COR; then appending the shifted image vertically to the original sinogram (Fig. 19(a)).

 

Fig. 19. Demonstration of the steps of the double-wedge filter. (a) 360-degree sinogram constructed from the 180-degree sinogram in Fig. 18(a). (b) Logarithmic scale of the 2D Fourier transform of image (a). (c) Double-wedge binary mask. (d) Result of applying the double-wedge filter for 10 times. (e) Difference between image (d) and the image in Fig. 18(a).

Download Full Size | PDF

The double-wedge property has been exploited to calculate the COR [18], sinogram restoration [51], or sample detection as shown in Section 2.1.2.2. A mathematical explanation of the formation of a double-wedge region is fully presented in Ref. [37]. Here we only use a formula which defines the wedge-angle of a double-wedge region to design the mask. The angle of a wedge (Fig. 19(c)) is calculated by

The routine of removing parts of a sample larger the width of its sinogram, dubbed the double-wedge filter (Code 1, Ref. [17]), is as follows:

The result of applying the double-wedge filter to the sinogram in Fig. 18(a) is shown in Fig. 19(d). As can be seen the external parts of the sample are removed. However, there is a side effect that the low-frequency components are altered as seen by comparing the difference between the sinograms before and after the filtering (Fig. 19(e)). We can reduce this effect by applying normalization to the processed sinogram (Fig. 19(d)) before reconstruction. Figure 20 shows improvements of the iterative reconstruction methods after the double-wedge filter is applied. Strong artifacts disappear as compared to Fig. 18.

 

Fig. 20. Improvement of the reconstructed results after using the double-wedge filter. (a) From the SART method with 500 iterations. (c) From the CGLS method with 100 iterations. (d) From the SIRT method with 100 iterations.

Download Full Size | PDF

The double-wedge filter can also be used to remove horizontal stripes in a sinogram which may be caused by the fluctuation of an X-ray source (Fig. 21(a) and 21(d)) or by the extinction contrast of crystalline samples (Fig. 21(b) and 21(e)) [52]. Note that these horizontal stripes do not have significant impact on the reconstruction quality of the FBP method but they may affect iterative reconstruction methods.

 

Fig. 21. Applications of the double-wedge filter. (a) Horizontal stripes caused by the fluctuation of an X-ray source. (b) Partially horizontal stripes caused by a crystalline sample. (c) Magnified view from the blue frame in (b). (d) Result of applying the double-wedge filter to image (a). (e) Result of applying the double-wedge filter to image (b). (f) Magnified view from the blue frame in (e).

Download Full Size | PDF

2.2.1.9 Sinogram stitching

This module of methods is used for processing data acquired using a grid scan, a half-acquisition scan, or their combination. It includes the methods of determining the overlap-side and the overlap-area between sinograms, as shown in Section 2.1.1.1, and the method for stitching them. In practice, it may happen that there are no overlaps between images. In such a case, Algotom also provides methods for joining sinograms based on interpolation.

2.2.1.10 Helical sinogram transformation

This module is used for a helical scan as presented in Section 2.1.3.1.

2.2.1.11 Denoising or contrast enhancement

In a synchrotron-based tomography system, the partial coherence of the X-ray source introduces a propagation-based phase contrast effect resulting in edge-enhanced images as shown in Fig. (1), Fig. (3), or Fig. (6). For applications which only use absorption-contrast information, this effect is unwanted as it causes streak artifacts in reconstructed images (e.g., Fig. (4)). To reduce the effect, a low-pass filter can be used to remove the high-frequency edges in the images. One of the well-known filters was published by Paganin et al. [53], implemented in many tomographic software packages [14,15]. As a low-pass filter it reduces noise and the dynamic range of an image. This helps to enhance the contrast between low-contrast features which can be confused if this enhancement comes from the phase effect. To properly retrieve phase shift from edge-enhanced images, a class of methods known as the propagation-based methods need at least two measurements at different sample-detector distances [54,55]. Paganin’s method is considered to be a phase retrieval method if it is applied to the image of a single homogeneous material [53]. However, the advantage of phase-contrast imaging over absorption-contrast imaging is only shown in a multi-material sample with small-density differences [56,57]. Considering Paganin’s formula which is rewritten as [53,58]

As R, i.e., δ/β, is undetermined for a multi-material sample, users must manually select this value as a trade-off between the contrast enhancement and the sharpness of a reconstructed image. In practice, R can be arbitrarily adjusted in the range of 10 to 10⁴. Because of that, the constant of (λz)/(4πΔx²) has no meaning and can be omitted which gives a much simpler formula

The adjustable range of the R parameter is not changed much in Eq. (18) compared to Eq. (17) because the value of (λz)/(4πΔx²) is about between 0.1 and 1 in the near-field hard X-ray regime with a micro-resolution detector. The 1D version of this filter for working in the sinogram space is given as

This 1D filter is useful to denoise or improve the contrast of reconstructed images from under-sampled sinograms in grid scans because applying Eq. (18) to all projection images before sinogram generation is very computationally expensive. Implementations of the filters in Eqs. (18) and (19) are shown in Code 1 (Ref. [17]).

2.2.2 Reconstruction methods

2.2.2.1 Direct reconstruction methods and edge-padding for removing cupping artifacts

Large-size images and high number of datasets limit the choices of reconstruction methods when users have to process their data in a reasonable timescale. Here, we mainly use the FBP method [28] or the DFI (direct Fourier inversion) method [43], depending on the availability of GPU resources. The under-sampling problem, caused by applying these methods to data having a low ratio between the number of projections and the image width, can be improved using denoising methods before or after reconstruction [53,59].

In FFT (fast Fourier transform)-based reconstruction methods, such as the FBP or DFI method, the 1D FFT is applied to a sinogram row-by-row before another processing step, e.g., interpolation, is used. FFT requires that an input is periodic, which is not the case for a sinogram of a sample larger than the FOV. This non-periodicity leads to artifacts in the processed sinogram. Zero-padding is often used to reduce these artifacts [1]. However, the sharp changes in intensities at padded edges introduce errors in the Fourier domain which give rise to cupping artifacts in reconstructed images (Fig. 22). Published techniques can be used to tackle the problem, known as the truncated projection problem [31–33]. They mainly utilize extrapolation techniques to compensate for missing data which are not straightforward to implement. Here we use a simple edge-padding technique to avoid the problem caused by sharp changes in intensities at the edges. It removes the cupping artifacts in the slices of both the FBP method and the DFI method as can be seen in Fig. 23. This technique may alter the DC-component of a sinogram. However, this variance is not significant and can be corrected by normalization. The technique is already used elsewhere for the FBP method [60] but not for the DFI method. Here, we demonstrate this good practice because some tomographic open-source software currently uses built-in zero-padding for the FFT-based reconstruction methods [14,47,61]. To show the advantage of using edge-padding, we implement the FBP method and the DFI method ourselves in the Algotom package (Ref. [17]).

 

Fig. 22. Cupping artifacts from a sample larger than the FOV. (a) Sinogram. (b) Reconstructed image using the FBP method in Astra toolbox [47]. (c) Reconstructed image using the DFI method in Tomopy [14]. (d) Line profile along the red line in (b) and the blue line in (c).

Download Full Size | PDF

 

Fig. 23. Demonstration of the cupping artifact removal. (a) Edge-padding sinogram. (b) Slice from our implementation of the FBP method (Ref. [17]). (c). Slice from our implementation of the DFI method (Ref. [17]). (d). (d) Line profile along the red line in (b) and the blue line in (c).

Download Full Size | PDF

2.2.3 Post-processing methods

2.2.3.1 Ring artifact removal

Ring artefacts still present in reconstructed images may be reduced by transforming reconstructed images between the polar coordinates and the Cartesian coordinates and applying FFT-based filters to remove ring artifacts. Routines for these procedures are included here.

2.2.3.2 Image downsampling and image rescaling

Methods for downsampling a 2D or 3D image dataset include the mean, the median, the maximum, or the minimum value of a pixel grid or a voxel grid. Methods for rescaling 32-bit data to 16-bit data or 8-bit data are also provided. It is important to be aware that rescaling causes information loss. The global extrema or user-chosen percentile of a 3D data or a 4D data (time-series tomography) need to be used for rescaling to limit the loss.

3. Results

3.1 Using a grid scan in combination with half-acquisition scans

Datasets used for demonstrating this combined scanning technique were collected at DLS, Beamline I12 with the following conditions: Monochromatic X-ray beam energy 76 keV. Imaging module 3 was used with a PCO.edge 5.5 sCMOS camera having a 2,560 × 2,160 pixels sensor. The optics module gives a pixel size of 3.2 µm resulting in a FOV of 8.2 mm wide and 6.9-mm high. To cover a whole sample, which is a pebble having a size of ∼3.6 cm × 5.1 cm × 2.2 cm (width × height × depth), a grid scan of 3 columns by 8 rows was performed giving 24 tomographic datasets in total. Each tomographic dataset contains 3601 projections evenly distributed in the range of [0°; 360°]. The rotation axis was offset to increase the width of the stitched FOV to approximately 6 times (2 × 3 image widths) of the FOV of the detector. Linear translations of the sample were selected to get the overlap between views of ∼0.8 mm, which is equivalent to ∼250 pixels. The scanning was done row-by-row to avoid the mismatch between sinograms of the same sample slice, which may be caused by the mechanical error of the vertical translation stage. Raw datasets were stored in HDF (Hierarchical Data Format) files accompanying by NeXus metadata files [62]. A critical advantage of using the HDF file format for tomography is that subsets of a 3D dataset can be extracted quickly and efficiently. This is important for a data processing pipeline because different methods in a pipeline may work in different slice dimensions.

Figure 24 shows how the data are processed. The following methods are used for removing artifacts: flat-field correction, zinger removal, and ring artifact removal methods which are the combination of algorithm 3, 4, 5, and 6 in [40]. Sinograms are stitched and converted to 180-degree sinograms using the method described in Section 2.1.1.1. The CORs are determined using Eq. (5). The denoising method with R = 300 is applied (Eq. (18)). Finally, the images are reconstructed using the FBP method.

 

Fig. 24. Demonstration of the data processing routine used for tomographic data acquired using the grid scan in combination with the half-acquisition scan.

Download Full Size | PDF

In practice, synchrotron users often reconstruct a few slices of a sample for quick decision-making during an experiment, or when optimizing parameters of data processing methods prior to performing full reconstruction of the sample. As shown in Fig. 24, reconstructing a few slices from data acquired using the combination of a grid scan and half-acquisition scans is straightforward and fully automated using the methods described in Section 2.1.1. Specifically, reconstructing a slice of the sample shown in Fig. 24 from the raw data takes around 2.5 minutes on a normal workstation at Beamline I12 having 16GB RAM memory and Intel Xeon E-2124 CPU. The computational cost of the auto-stitching and auto-centering steps is most expensive taking nearly 2 minutes. The pre-processing step of flat-field correction and the gridrec reconstruction [14] cost only 0.5 minutes, giving the reconstructed slice with a size of ∼14k × 14k pixels. One needs to calculate the overlap-areas and overlap-sides for stitching sinograms only once. After that time for reconstructing another slice is below 1 minute.

For reconstructing the whole sample, the approach is slightly different. Firstly, overlap-areas and overlap-sides for stitching and converting sinograms are calculated column-by-column and row-by-row between cells in a grid. This is necessary to determine the global width of reconstructed images. Because the overlap-areas may change between the rows of the grid, possibly caused by the mechanical error of the translation stage, the widths of the reconstructed images may vary correspondingly. The global width chosen is the maximum value among the rows of the grid. Sinograms are then padded before reconstruction to compensate for the mismatch. For stitching in the vertical direction, the overlap-sides and overlap- areas between the grid rows are calculated using the method in Section 2.1.1.1, which is applied to the projection images at the same angle between the grid rows. Sinograms at the top and bottom of two grid cells at the same column can also be used for determining the overlaps, but this approach is more costly due to the overhead of extracting sinograms. Knowing the overlaps, sinograms of the same slice are combined using the distance-weighted average before reconstruction as shown in Fig. 24.

Although reconstructing the whole sample into a massive single HDF file is straightforward, it is not practical for post-analysis or transferring the file. Here, we divided the full sample into 12 sub-sections and reconstructed them independently using 12 GPU nodes of the DLS Hamilton cluster. It took nearly a day for processing all the data using all the processing methods described above. The shape of the 3D reconstructed image (32-bit) is of ∼15k × 14k × 14k pixels occupying 11.3 TB of the storage system. The 3D rendering of the 3D reconstructed image, from the Avizo software, is shown in Fig. 25(a) and 25(b) where the image volume was downsampled by a factor of 8 to make the visualization possible.

 

Fig. 25. (a) Vertical slice of the 3D reconstructed image. (b) 3D rendering of the 3D reconstructed image.

Download Full Size | PDF

3.2 Using a helical scan in combination with a half-acquisition scan

A dataset used for demonstration this type of scanning was also collected at the beamline I12 using the same setup as in Section 3.1. The monochromatic X-ray beam energy was 54 keV. Parameters of the helical scan are: y_start = 3.75 mm; y_stop = 19.5 mm; P_t= 6.5 mm; and N₁₈₀ = 1,801. The pixel size determined using the routine in Section 2.1.3.1 is of 3.24 µm. The rotation axis was offset to approximately double the FOV. The slice thickness, Δy, is chosen to be equal the pixel size. Figure 26 shows some projections of the tomographic data where Fig. 26(c) is created by combining the last columns of all projection images. It helps to show the direction of extracting sinograms for reconstruction and the areas (the blue overlays) where the data have no use or are not reconstructable. Figure 26(c) also shows another possible way of generating sinograms which includes two steps [15]: applying a shearing transform in the same plane as in Fig. 26(c); extracting columns of the resulting images. This approach, however, requires storing of the intermediate data and its performance relies on the IO performance of a storage system.

 

Fig. 26. Some projections acquired using the helical scan in combination with the half-acquisition scan. (a) Projection at 0°. (b) Projection at 360°. (c) 2D image generated by combining the last columns of all projection images. The red line indicates the direction of sinogram extraction.

Download Full Size | PDF

Data processing methods used in the pipeline (Fig. 27) for this type of scan are: flat-field correction, circular sinogram generation including the blob removal (Section 2.1.3), partial ring artifact removal (algorithm 3 in [40]), zinger removal, 180-degree sinogram conversion using the auto-centering and auto-stitching methods in Section 2.1.1.1, and FBP method for reconstruction. A bottleneck in the processing pipeline is the step for generating the circular sinograms where the program walks through the 3D dataset and extract subsets of the data. Here, the time cost was around 10s per sinogram for this dataset on the DLS Hamilton cluster.

 

Fig. 27. Demonstration of the data processing routine used for tomographic data acquired using the helical scan in combination with the half-acquisition scan.

Download Full Size | PDF

4. Discussion and conclusion

Although there are many data processing methods presented in the paper and they are all important in each scanning technique and data processing pipeline, we would like to highlight a few important contributions of our work.

Firstly, there is a method for determining the overlap-side between sinograms. This method leads to the determination of the overlap-area between sinograms of a grid scan and/or the center of rotation in a half-acquisition scan. Practically, the method searches the full width of a sinogram in both left and right directions with low computational cost. This enables automated data processing. Explanations of the feasibility of the proposal method, in comparison to other approaches, for tomographic imaging are shown in Section 2.1.1.1.

Secondly, there is a combination of a grid scan with half-acquisition scans. This scanning technique allows a multiple doubling of the FOV of a parallel-beam tomography system with the least number of stage translations. Practically, it makes use of the proposed method for both sinogram stitching and 360-degree sinogram conversion. This enables to reconstruct a slice of the grid scan fully automatically with low computational cost.

Thirdly, it is a detailed procedure of calibration, data acquisition, and data conversion for helical scans. Although it is easy to perform a helical scan, processing the acquired data is very cumbersome without a proper calibration and acquisition. The pixel size of a detector in relation to the translation of a stage needs to be calibrated as X-ray beams are not perfectly parallel. The start and stop position of the translation stage need to be calculated properly to be able to reconstruct the region of interest. Details of the conversion algorithm and the methods of removing characteristic artifacts are also provided. All these contributions allow users or beamline scientists to perform helical scans with ease.

Lastly, it is software which not only implements all the proposed methods but also common methods in a data processing pipeline. It provides users with a complete set of tools to reconstruct data acquired by the scanning techniques described in this paper with ease and the highest possible quality. They all help to extend applications of parallel-beam tomography systems.

As stated about the scope of this work, we have demonstrated data processing methods and software (Ref. [17]) developed for reconstructing data acquired by the scanning techniques, which extend the FOV of a parallel-beam tomography system, including half-acquisition scans, grid scans, helical scans, and their combinations. Methods in a full data processing pipeline are provided in the Algotom package to give users a complete set of tools for tomographic reconstruction with the highest possible quality. Considering the importance of decision making during an in-situ experiment, a very large data can be produced in a single beamtime, and a high demand of delivering a high-throughput system of data acquisition and data processing at a synchrotron facility; our approaches were designed to be easy-to-use, automatable, and efficient; e.g., a large slice of a grid scan can be reconstructed fully automatically in a few minutes using a normal workstation computer as shown in Section 3.1.

Although the methods and the software have been developed and tested on parallel-beam X-ray tomography systems, they should be applicable to other types of parallel-beam tomography such as neutron tomography, X-ray fluorescence tomography, or electron tomography. Certainly, characteristic pre-processing methods of image alignment for electron tomography, scattering correction for neutron tomography, or self-absorption correction for fluorescence tomography need to be added to the software. To distribute the work and further development, the software is also made available at: https://github.com/algotom/algotom. Details of how to use the Algotom package to reconstruct images from the raw data presented in Section 3.1 can be found in the “examples” folder of Code 1 (Ref. [17]). Links to the raw data uploaded to zenodo.org are provided inside the codes.