Next Article in Journal / Special Issue
The Neutron Imaging Instrument CONRAD—Post-Operational Review
Previous Article in Journal
Imaging Meets Cytometry: Analyzing Heterogeneous Functional Microscopic Data from Living Cell Populations
Previous Article in Special Issue
Neutron Imaging Using a Fine-Grained Nuclear Emulsion
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Improved Acquisition and Reconstruction for Wavelength-Resolved Neutron Tomography †

1
Multimodal Sensor Analytics Group, Oak Ridge National Lab, Oak Ridge, TN 37831, USA
2
Neutron Scattering Division, Oak Ridge National Lab, Oak Ridge, TN 37831, USA
*
Author to whom correspondence should be addressed.
This report has been authored by UT-Battelle, LLC, under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (https://www.energy.gov/downloads/doe-public-access-plan).
J. Imaging 2021, 7(1), 10; https://doi.org/10.3390/jimaging7010010
Submission received: 31 October 2020 / Revised: 17 December 2020 / Accepted: 22 December 2020 / Published: 15 January 2021
(This article belongs to the Special Issue Neutron Imaging)

Abstract

:
Wavelength-resolved neutron tomography (WRNT) is an emerging technique for characterizing samples relevant to the materials sciences in 3D. WRNT studies can be carried out at beam lines in spallation neutron or reactor-based user facilities. Because of the limited availability of experimental time, potential imperfections in the neutron source, or constraints placed on the acquisition time by the type of sample, the data can be extremely noisy resulting in tomographic reconstructions with significant artifacts when standard reconstruction algorithms are used. Furthermore, making a full tomographic measurement even with a low signal-to-noise ratio can take several days, resulting in a long wait time before the user can receive feedback from the experiment when traditional acquisition protocols are used. In this paper, we propose an interlaced scanning technique and combine it with a model-based image reconstruction algorithm to produce high-quality WRNT reconstructions concurrent with the measurements being made. The interlaced scan is designed to acquire data so that successive measurements are more diverse in contrast to typical sequential scanning protocols. The model-based reconstruction algorithm combines a data-fidelity term with a regularization term to formulate the wavelength-resolved reconstruction as minimizing a high-dimensional cost-function. Using an experimental dataset of a magnetite sample acquired over a span of about two days, we demonstrate that our technique can produce high-quality reconstructions even during the experiment compared to traditional acquisition and reconstruction techniques. In summary, the combination of the proposed acquisition strategy with an advanced reconstruction algorithm provides a novel guideline for designing WRNT systems at user facilities.

1. Introduction

WRNT is a technique used to characterize samples in 3D based on the time-of-flight (TOF) technique that is capable of identifying the time of arrival of each neutron, thus its wavelength. Each WRNT scan involves collecting WR projection data of an object from multiple orientations and then using an algorithm to obtain a hyper-spectral 3D reconstruction. Figure 1 shows a schematic of a WRNT instrument where a spectral source of neutrons, typically produced by the spallation of neutrons from the interaction of a proton beam with a heavy element target such as tungsten or mercury, impinges on the sample of interest which is then rotated to obtain a full set of tomographic measurements. WRNT is performed at spallation neutron sources based facilities [1,2,3,4] by using a pulsed (spectral) source combined with a time-of-flight detector [5] that allows events of different wavelengths, which are inversely proportional to their kinetic energy, to be recorded according to arrival time at the detector. Wavelength-resolved neutron imaging has been used to study the phase composition of samples [6,7] and to investigate the crystalline nature of engineering materials [8] using the principles of Bragg scatter i.e., making use of the fact that crystals diffract strongly when they satisfy the Bragg condition at a certain orientation and neutron wavelength [9].
While WRNT beam lines are in operation, there are challenges that still have to be addressed in order to produce high-fidelity tomographic reconstructions at these facilities. Specifically, the first phase of innovations were mainly focused on the instrument design for producing the highest flux possible, and efficient detectors to get high-quality measurements. As a result, there has been relatively less focus on the joint design of acquisition and reconstruction algorithms to further improve the performance of the tomographic system. A typical approach for the tomography involves rotating the sample sequentially in a 180 or 360 range and obtaining a dataset corresponding to hundreds of wavelength-resolved projection measurements at each position (see Figure 1)—a sequential step-and-scan approach. The data from all the measurements are then processed at the end of the experiment using an analytic reconstruction algorithm such as filtered back projection (FBP) to obtain a tomographic reconstruction [7,10]. However, this approach has a few drawbacks. Firstly, the measured data can be of very low signal-to-noise (SNR) in each wavelength bin due to the limited time available for an experiment. Next, the limited availability of beam time at facilities often forces users to acquire a sparse set of projection data. This results in a reconstruction with significant artifacts when FBP-like algorithms are used for this sparse-view and noisy dataset. Finally, the time taken to acquire a single tomographic scan can be long with experiments spanning several days during which the user may not have accurate feedback on the current state of the reconstruction from the thus far acquired data. In summary, while there has been significant improvements in WRNT systems, most of these have focused on the hardware, leaving open the possibility of further improvements via the design of new acquisition and reconstruction algorithms.
In this paper, we propose a new framework for WRNT acquisition and reconstruction based on the use of an interlaced scanning strategy combined with a model-based image reconstruction algorithm. A key benefit of our approach is that we can produce intermediate reconstructions which significantly suppress noise and artifacts compared to traditional approaches. This also implies that the final reconstructions are available at the end of the experiment time and are of higher quality (lower noise for a given resolution) compared to the FBP methods in [7,10]. Specifically, we adapt the scanning technique of the time-interlaced model-based iterative reconstruction (TIMBIR) algorithm [11] designed in the context of time-resolved synchrotron based X-ray tomography in order to improve the spatio-temporal resolution. We note that an approach which is intuitively similar has also been explored in the context of white-beam time-resolved neutron tomography where a acquisition based on the golden ratio [12] has been used to improve the spatio-temporal resolution of the reconstructions [13,14,15]. While the two acquisition schemes are similar in spirit, the interlaced scanning algorithm of [11] is more general than the golden section method and can be made to be very similar to the golden section method by setting the appropriate parameters during acquisition. In contrast to the sequential step-and-scan approach, the interlaced scan ensures that subsequent measurements are further apart so that there is sufficient diversity in the data at any given time during data acquisition. A benefit of such a strategy is that if there are unexpected downtimes in the course of an experiment, it allows for a more uniform set of measurements to be available for reconstruction compared to traditional sequential scanning protocols. Secondly, we propose to use of a model-based image reconstruction (MBIR) algorithm on the acquired wavelength-resolved/hyper-spectral data instead of the standard FBP techniques. The MBIR algorithm works by minimizing a cost function that balances a data-fitting term that incorporates a physics based model for the imaging system and noise characteristics of the detector, and a regularization term that incorporates a model for the underlying 3D object itself (such as local smoothness). MBIR methods have enabled significant improvements in performance for several tomography applications including for neutron tomography [16,17,18,19] especially when dealing with sparse and low signal-to-noise ratio (SNR) data; a scenario that is also common in WRNT as we seek to provide high-quality real-time feedback to users in the course of an experiment. We demonstrate the utility of our proposed method by implementing this system at the SNAP beam line of the Spallation Neutron Source (SNS) at the Oak Ridge National Laboratory (ORNL) and demonstrate real-time feedback capability during the course of a WRNT of a magnetite sample.
The rest of this paper is organized as follows. In Section 2 we present specific details of the proposed acquisition algorithm and the MBIR method. In Section 3 we provide experimental results from the measurement of a magnetite sample highlighting the utility of the proposed method for WRNT. Finally, in Section 4 we present our conclusions and discuss future research directions.

2. Interlaced Scanning and Model-Based Image Reconstruction

Our objective in the design of WRNT systems is to be able to make measurements such that we can provide real-time feedback to users via intermediate high-quality reconstructions and also ensure the same high-quality reconstruction upon completion of the experiment. We achieve these objectives via the design of a new acquisition protocol and combine it with a MBIR algorithm that can produce high-quality reconstructions from sparse, noisy and in-complete data sets.

2.1. Interlaced Scanning

The conventional approach to a tomographic scan at a WRNT beam line is to rotate the sample about a single axis, and at each position measure a wavelength-resolved projection dataset. Typically, the sample is rotated using a fixed step size which is chosen based on the total number of measurements to be made in the course of the experiment. The sequence of angles usually follows a linear increasing pattern—starting from 0 and progressing to 180 using a fixed step size about the axis of rotation (see Figure 2a). However, the disadvantage of such a scanning protocol is that at any given intermediate time in the course of the scan, the object is only viewed from a very limited angular range; making the reconstruction obtained from such partial data to have significant streak artifacts. This can also lead to challenges if, due to some unforeseen circumstance, beam time is lost during an experiment thus producing a dataset which is heavily sampled only from a small angular range. In order to address this challenge, we require that the measurements are spread apart and process the sparse set of measurements using appropriate algorithms so that we can provide a continuous high-quality reconstruction to the end users (see Figure 2b).
2
Here, we propose to use an interlaced scanning strategy [11] to acquire the data. Based on the time available for a given scan (with a reasonable SNR), we can decide on the total number of measurements we plan to measure ( N θ ). Once N θ is determined, the conventional data acquisition generates a uniform set of N θ measurements between 0 and 180 (see Figure 2a,c). In contrast to the conventional data acquisition approach described above, we propose to use the method described in [11] which utilizes a sequence of angles that are interlaced so that subsequent angles are far apart. In this algorithm, in addition to N θ , we choose a number, K, which represents the number of “half-rotations” over which to acquire the data. This K ( N θ ) can be chosen by the user who desires to obtain roughly K reconstructions in the course of an experiment (though other heuristics can be also used). Once these numbers are chosen, the list of angles, θ n , (in the range 0 to 180 ) is given by
θ n = n mod N θ K K + B r n K N θ mod K 180 N θ
n 0 , . . , N θ 1 , mod is the standard modulo operator, and  B r is the bit-reversal operator that takes the integer representation of a number and reverses the order of the bits [20]. An example of this sequence for N = 36 and K = 3 is shown in Figure 2d and is contrasted with the standard sequential scanning protocol in Figure 2c. Notice that both the sampling methods acquire the same set of views, but the main difference is in the order in which these views are measured. When K = N θ the resulting sequence will be maximally separated i.e., each subsequent measurement is farthest apart from the previous, which is in complete contrast to the conventional scanning where each subsequent measurement is closest to the previous one for a given N θ . Hence, at the end of the experiment, both protocols will have measured the same set of views but the interlaced scan acquires data in a manner that covers the full angular range better even at intermediate times in the experiment. The advantage of interlaced scanning is that if there is a disruption, say mid-way through an experiment, we can obtain a reconstruction with fewer artifacts than from a conventional scan. Furthermore, if we desire to provide intermediate reconstructions to the end-users during an experiment, this scanning protocol can enable higher-quality reconstructions for a generic sample compared to the traditional scan as we will demonstrate in the results section. We note that such sampling patterns are best used when the time to rotate the object from one position to the next is smaller than the time to make the measurement itself, which is typical for WRNT systems where it can take of the order of an hour or more to acquire one projection measurement compared to a few seconds to align the sample in a particular direction with respect to the beam. Finally, we note that for this case the interlaced scanning strategy is similar but not exactly the same as the golden ratio based scan of [13] that has been investigated for time-resolved neutron CT. One of the advantages of the proposed strategy compared to the golden section method is that the parameter K can be adjusted to control how separated successive angles are.

2.2. Model-Based Image Reconstruction for Streaming Hyper-Spectral Data

Model-based image reconstruction [21] approaches have been developed for the last few decades and are well established to be the method of choice when dealing with sparse, limited-view and noisy tomographic datasets. In the scientific user facility community, such methods have been developed for electron tomography [22,23,24], synchrotron based X-ray CT [11,25,26] and neutron laminography and tomography [16,17,18,19]. In this section, we briefly explain the core ideas behind MBIR and how we have adapted it for the WRNT systems in order to perform hyper-spectral tomography. We skip some of the mathematical details since these are well established in the literature.
At any given stage in the course of the experiment, when n views have been measured, the MBIR is obtained as
f ^ ( n ) argmin f c ( n ) ( f )
where
c ( n ) ( f ) = 1 2 t = 1 n g ( t ) A ¯ ( t ) f W ( t ) 2 + R ¯ ( f ; ψ )
where g ( t ) is the vector of normalized hyper-spectral projection measurements at view t organized as g = g 1 ( t ) g S ( t ) , g s t is a vector of all normalized measurements at view t and wavelength bin s, f is the vector containing all the unknown voxels of hyper-spectral attenuation coefficients organized as f = f 1 f S , A ¯ ( t ) is a matrix that represents that forward-projection operator at view t at all the wavelengths (a block diagonal matrix with entries A ( t ) representing the forward projection at view angle θ t 1 ), W ( t ) is a diagonal weight matrix containing the inverse noise-variance of the measurements g ( t ) , and R ¯ ( f ; ψ ) is a regularization term that enforces certain desirable properties in the reconstructed volume (like local smoothness, sharp edges etc.). Notice that this is a very high-dimensional optimization problem over all unknown voxels across all wavelength bins. The entries of the weight matrix ( W ( t ) ) are typically set to be the measured (un-normalized) count data in the case of transmission tomography [27] and have the intuition that if this number is very small, the overall term weights less relative to other terms in the cost-function. The inclusion of this term allows for a dose-weighting of the acquired data and can be a useful way to reject some of the individual measurements. In summary, the MBIR approach involves formulating the tomographic reconstruction as minimizing a function that balances a fitting term that enforces a low discrepancy between the projection of the final reconstruction and the true measurements; and a regularization term that enforces certain desirable properties in the reconstruction.
There are many potential choices for the regularization term R ¯ in Equation (3); some of which can enforce desirable properties in each wavelength bin and others that can model the correlations across bins (i.e., adjacent bins will be correlated). In this paper, we choose an R ¯ which regularizes each wavelength bin independently.
R ¯ ( f ; ψ ) = s = 1 S R ( f s ; ψ s )
where f s is a vectorized version of all the voxels in the sth wavelength bin, and ψ s is the set of parameters associated with that bin. For R ( ; ) we choose a q-generalized Markov-random field (qGGMRF)-based function [28]. It is given by
R ( x ; ψ s ) = { j , k } N w j k ρ ( x j x k ; ψ s ) ρ ( x j x k ; ψ s ) = x j x k σ x , s 2 c + x j x k σ x , s 2 p s
where ψ s = [ p s , σ x , s , c ] are the parameters of the function, N is the set of pairs of neighboring voxels (e.g., a 26 point neighborhood), 1 p 2 , c and σ x are qGGMRF parameters. c is typically set to a very small number and is used to make the overall function differentiable. The weights, w j k , are inversely proportional to the distance between voxels j and k, normalized to 1. This model provides a greater degree of flexibility in the quality of reconstructions compared to an algorithm specifically designed for a total-variation regularizer that may force the reconstructions to appear “waxy” [21]. In particular, when  p = 1 we get a behavior similar to a total-variation model and when p = 2 the regularizer is a quadratic function allowing for smoother reconstructions. Finally, we note that the choice of regularizer can be significantly improved further compared to the above choice. Exploiting the correlations across wavelength channels in the MBIR framework can help further improve the reconstruction quality as has been observed in several hyper-spectral imaging applications [25,29,30,31]. Furthermore, if we are measuring samples with known materials and spectral responses, this can be used as a part of the model in the MBIR framework to reduce the dimensionality of the problem dramatically (instead of reconstructing every voxel at each channel we only need to find the proportion of each material at each voxel) [31]. However, in this proof-of-concept study we implement the regularizer (5) for its simplicity and the fact that it enables parallel reconstruction of each of the individual hyper-spectral bins.
With these choices of functions, the cost-function (3) can be re-written as
c ( n ) ( f ) = 1 2 t = 1 n s = 1 S g s ( t ) A ( t ) f s W s ( t ) + s = 1 S R ( f s ; ψ s ) .
Because we have chosen a regularizer that is separable, this can be re-written as
c ( n ) ( f ) = s = 1 S 1 2 t = 1 n g s ( t ) A ( t ) f s + R ( f s ; ψ s )
Therefore the hyper-spectral reconstruction after measuring n views can be computed as
f ^ s ( n ) argmin f s 1 2 t = 1 n g s ( t ) A ( t ) f s + R ( f s ; ψ s )
s 1 , . . . , S . Finally, we note that even though
c ( n ) ( f ) = c ( n 1 ) ( f ) + 1 2 g ( n ) A ¯ ( n ) f W ( n ) 2
n 1 , the MBIR method involves processing all the data in order to obtain a reconstruction. However, we can attain an accelerated convergence to the solution by initializing the optimization routine at the nth stage with the solution from the ( n 1 ) th stage. In contrast to this, linear methods used for tomography like the filtered-back projection can be implemented by keeping track of the current reconstruction and simply adding the contribution from the new measurement.
We use the optimized gradient method (OGM) [32] to find a minimum of the above cost function for each wavelength channel in equation (8). The algorithm involves a standard gradient computation combined with a step size determined using Nesterov’s method. In order to prevent the notation from getting cumbersome, we drop the indices corresponding to the wavelength channel (s) and the view index (t) and merely highlight what forms the core of each iteration. Specifically, for each iteration k,
h ( k + 1 ) x ( k ) 1 L c ( x ( k ) )
t ( k + 1 ) 1 + 1 + 4 ( t ( k ) ) 2 2
x ( k + 1 ) h ( k + 1 ) + t ( k ) 1 t ( k + 1 ) ( h ( k + 1 ) h ( k ) ) + t ( k ) t ( k + 1 ) ( h ( k + 1 ) x ( k ) )
where t ( 0 ) = 1 , L is the Lipschitz constant of the gradient of c ( . ) , h ( 0 ) = x ( 0 ) is an initial estimate for the reconstruction. The gradient of the cost-function c ( . ) is given by
c ( x ) = A T W ( g A x ) + R ( x ; ψ ) .
We use the ASTRA toolbox [33,34] to implement GPU accelerated forward (A) and back-projection ( A T ) operators. We also accelerate the computation of the expression for R ( x ; ψ ) using a parallel implementation. Overall, this leads to a algorithm which can be implemented rapidly for each wavelength channel and across all channels in parallel.
Our overall framework for acquisition and reconstruction is summarized in Algorithm 1. We note that it is not necessary to reconstruct the data in order to proceed to the next measurement in the proposed framework. Furthermore, MBIR methods are computationally expensive and obtaining real-time reconstructions depends on the acquisition time, size of the detectors and the number of wavelength bins. In practice, for the type of samples and detectors used in this paper we have observed that it is possible to obtain a useful (partial) reconstruction from the data in lesser time than it takes to measure the next projection using our implementation of the algorithm. We believe that using even more optimized implementations [35,36] we can obtain large hyper-spectral reconstructions in near-real time for WRNT systems using fairly small compute clusters. We plan to open-source our reconstruction code as a part of the pyMBIR package [37] upon publication of the current work.
Algorithm 1: The proposed framework for wavelength-resolved neutron tomography systems. Prior to the scan we decide a fixed number of views to measure N θ and the number of half-rotations over which this will occur K. After each measurement the object can be reconstructed using the MBIR algorithm. The reconstruction from the previous stage is used as a initial condition for the next stage in order to accelerate the convergence i.e., time to reconstruct for the MBIR method.
1: function [ f ^ ( N θ ) ] Reconstruct( N θ ,K, ψ )
2:    %Inputs: Number of views N θ , Number of half-rotations K, Regularization parameters ψ
3:    %Outputs: Wavelength-resolved reconstruction f ^ ( N θ )
4:    Compute θ n according to (1) n 0 , . . , N θ 1 ▹ Interlaced scanning angles
5:    Initialize f ^ ( 0 ) 0
6:    Measure open-beam data λ o , s s 1 , . . , S ▹ Normalization data for transmission tomography
7:     n 1
8:    while n N θ do
9:         λ ( n ) MEASURE P R O J E C T I O N ( θ n 1 )
10:         g ( n ) NORMALIZE( λ ( n ) , λ o )▹ Can include further corrections
11:         f ^ ( n ) MBIR ( g ( 1 ) , . . , g ( n ) ; ψ , f ^ ( n 1 ) )▹ Equation (8)
12:         n n + 1
13:    end while
14:    Return f ^ ( N θ )
15: end function

3. Results

In order to illustrate our approach, we acquired a dataset consisting of two magnetite samples (see Figure 3) using the SNAP instrument at the Spallation Neutron Source at Oak Ridge National Laboratory. The samples were taped to a Aluminum post and separated by a piece of Cadmium. The acquisition time for each projection image was set so that the accumulated proton charge was set to 4 Coulombs and a 512 × 512 pixel micro-channel plate detector (MCP) [5] was used to acquire the data. The acquisition hardware was programmed in order to acquire an interlaced scanning pattern using K = 3 and N = 36 (see Figure 2). We note that due to a change in beam power, the acquisitions took longer than anticipated and we only acquired 30 projection measurements in the course of the allotted beam time. Furthermore, for the exposure time we had chosen (approximately 1.5 h corresponding to the accumulated charge of 4 Coulombs) the signal-to-noise ratio of the data was low (see Figure 3) with counts in the range of approximately 10 to 50 in a region of the sample in each wavelength bin. For each projection measurement, we binned the data into a total of 1400 time-of-flight/wavelength bins and used the first 1200 bins for the reconstructions. As the data was acquired, we pre-processed it using a median filter with a window size of 3 to reduce the impact of outliers, combined four adjacent wavelength bins by summing them in order to boost the SNR, normalized the data by applying the log operation to the ratio of the open-beam measurements to the acquired measurements with the sample, and then applied a stripe removal filter (Fourier wavelet algorithm) using TomoPy [38] to suppress ring-artifacts in the reconstruction. While we accumulated the counts in the four adjacent bins for these results, in the future we anticipate being able to do so during acquisition i.e., choosing the size of the wavelength/time-of-flight bins so that the acquired data has a sufficient signal-to-noise ratio and satisfies the wavelength resolution needed for the specific scientific use case. We then reconstructed the resulting normalized and binned data using the FBP and MBIR methods and compared the reconstructions obtained using two different scanning protocols— (1) The conventional method where we sequentially measure the sample from 0 to 180 degrees in the course of an experiment and (2) The interlaced scanning protocol proposed in this work. We note that since we only performed the experiment once using the interlaced scanning protocol, the results from the conventional scanning protocol are obtained by collecting all the data and retrospectively processing it as if it were acquired using the conventional method. The filter cut-off parameter for FBP was set to 0.35 for all wavelength bins and the Ram-Lak filter was used. The p s was set to 1.2 for MBIR and the σ x , s was set to 0.4 for all the wavelength-bins. These parameter values produced reconstructions of reasonable visual quality i.e., a reasonable balance between noise levels and visual resolution. The question of how to choose the parameters is an open question for WRNT with our observation that this is often done in a heuristic manner to attain a reasonable visual quality. For FBP there is a filter parameter for each reconstruction channel, and the same holds for the type of regularizer used in the presented MBIR approach resulting in a large number of parameters if they have to be objectively chosen for each channel. The time required to reconstruct a 3D sub-volume comprising 4 slices of size 460 × 460 from one wavelength bin using the MBIR approach and all the projection views was approximately 1 min (maximum of 250 MBIR iterations) using a machine with 40 cores and 2 Tesla P40 GPUs. Therefore depending on the available computing resources and desired rate at which we seek feedback from the system we can select the appropriate number of slices and wavelength bins to reconstruct over. We note that this compute can be dramatically accelerated with highly optimized MBIR implementations [35,36]; in this paper we focus on demonstrating a proof-of-concept that highlights the quality improvements that can be achieved using the proposed system.
Figure 4 show the reconstructions of a single cross section at various stages of the two acquisition protocols when processed using FBP and MBIR algorithms. We emphasize that these results compare reconstructions using different scanning strategies over the 180 angular range and different reconstruction algorithms corresponding to the acquired data; and not a comparison of how different algorithms perform when we have a sparse projection dataset as might be typical. Notice that even after about 22.5 h of experiment time (15 views), we can start seeing the features of the underlying sample clearly using the MBIR method combined with an interlaced scanning strategy. This illustrates that merely acquiring the data using a standard method combined with a MBIR algorithm is not the best strategy. This result highlights that the interlaced scanning strategy can pave the way for instrument users to be able to view reconstructions in real time and start getting meaningful feedback much prior to the completion of the experiment. This will be crucial for WRNT systems since beam time is at a premium and the making the instrument robust to disruptions can be highly impactful. In Figure 5 we quantify the accuracy of the proposed method by comparing the normalized root mean-squared error (NRMSE) and structural-similarity (SSIM) index [39] of a single cross-section between the current reconstruction and the final reconstruction upon completion of the experiment using a given algorithm. This gives a measure of how similar the “current” reconstruction is to the final solution and how fast it is converging to the desired solution when using the MBIR method but different acquisition protocols. Once again notice that these metrics are consistent with the visual inspection, further highlighting the value of combining an interlaced scan with a MBIR method for WRNT systems.
Figure 6 show the overall qualitative improvements in the final reconstruction between the MBIR and the FBP method for WRNT. Notice that the MBIR method can significantly suppress noise while preserving the resolution of the reconstruction compared to the FBP algorithm. The edges are clearly reconstructed and even the region of the cadmium plate which is highly attenuating is clearly reconstructed compared to the FBP algorithm. Furthermore, a plot of the average reconstructed value in a central patch across the wavelength dimension (Figure 7) shows the significant reduction in noise that can be obtained using the proposed MBIR technique. We believe that these results can be further improved by the use of more sophisticated regularization techniques that exploit correlations across adjacent wavelength bins. Finally, we note that the minor difference in qualitative and quantitative results of MBIR using all the data but in the conventional scanning mode and interlaced mode is due to the effect of the strip-removal filter in TomoPy, whose output varies according to the ordering of the input views.

4. Conclusions

In this paper we presented a new method for acquisition and reconstruction for wavelength-resolved neutron tomography beam lines. The core idea behind our approach is to combine an interlaced scanning protocol with the use of a model-based reconstruction algorithm that can produce high-quality reconstructions from low SNR, sparse and incomplete datasets. Our method is ideal for step-and-shoot type systems where the time to acquire the data is significantly higher than the time to position the sample in a certain orientation. An advantage of our approach is that the combination of a new acquisition protocol combined with a model-based reconstruction algorithm can produce high-quality reconstructions even during the course of the experiment; making it more robust to potential disruptions in the experiment and also enabling users to visualize their sample in real-time paving the way for them to modify/terminate the experiment based on this feedback. While this method is ideal for samples with unknown shapes and distribution of features, in the next phase of research we foresee designing sample-adaptive scanning protocols that can enable even more efficient acquisition of the data along the lines of [40].

Author Contributions

All the authors worked to conceive the overall idea. S.V. and P.B. developed the acqusition protocol. S.V. developed the MBIR algorithm and wrote the code. Y.Z. and H.B. designed the experiment at the beam line including the implementation of the new scanning protocol. C.H. contributed to the choice of sample for the demonstration experiment. L.D. advised with the choice of experimental parameters such as beam-energies by conducting simulations based on the magnetite material. All authors have read and agreed to the published version of the manuscript.

Funding

S.V.Venkatakrishnan was supported by ORNL’s Eugene P. Wigner Distinguished Fellowship program.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study is available on request from the corresponding author.

Acknowledgments

We thank Erik Stringfellow for helping with preparing the sample for measurement at the SNAP beam line. This research used resources at the Spallation Neutron Source, a DOE Office of Science User Facility operated by the Oak Ridge National Laboratory.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Strobl, M. Future prospects of imaging at spallation neutron sources. Nucl. Instrum. Methods Phys. Res. Sect. A Accel. Spectrometers Detect. Assoc. Equip. 2009, 604, 646–652. [Google Scholar] [CrossRef]
  2. Shinohara, T.; Kai, T.; Oikawa, K.; Segawa, M.; Harada, M.; Nakatani, T.; Ooi, M.; Aizawa, K.; Sato, H.; Kamiyama, T.; et al. Final design of the energy-resolved neutron imaging system RADEN at J-PARC. J. Phys. Conf. Ser. 2016, 746, 012007. [Google Scholar] [CrossRef]
  3. Kockelmann, W.; Minniti, T.; Pooley, D.E.; Burca, G.; Ramadhan, R.; Akeroyd, F.A.; Howells, G.D.; Moreton-Smith, C.; Keymer, D.P.; Kelleher, J.; et al. Time-of-Flight Neutron Imaging on IMAT at ISIS: A New User Facility for Materials Science. J. Imaging 2018, 4, 47. [Google Scholar] [CrossRef] [Green Version]
  4. Nelson, R.O.; Vogel, S.C.; Hunter, J.F.; Watkins, E.B.; Losko, A.S.; Tremsin, A.S.; Borges, N.P.; Cutler, T.E.; Dickman, L.T.; Espy, M.A.; et al. Neutron Imaging at LANSCE—From Cold to Ultrafast. J. Imaging 2018, 4, 45. [Google Scholar] [CrossRef] [Green Version]
  5. Tremsin, A.S.; Vallerga, J.V.; McPhate, J.B.; Siegmund, O.H.; Raffanti, R. High Resolution Photon Counting With MCP-Timepix Quad Parallel Readout Operating at >1 KHz Frame Rates. IEEE Trans. Nucl. Sci. 2013, 60, 578–585. [Google Scholar] [CrossRef]
  6. Myhre, K.G.; Zhang, Y.; Bilheux, H.Z.; Johnson, J.A.; Bilheux, J.C.; Miskowiec, A.J.; Hunt, R.D. Nondestructive Tomographic Mapping of Uranium and Gadolinium Using Energy-Resolved Neutron Imaging; Technical Report; Oak Ridge National Lab. (ORNL): Oak Ridge, TN, USA, 2018. [Google Scholar]
  7. Woracek, R.; Penumadu, D.; Kardjilov, N.; Hilger, A.; Boin, M.; Banhart, J.; Manke, I. Neutron Bragg-edge tomography for phase mapping. Phys. Procedia 2015, 69, 227–236. [Google Scholar] [CrossRef] [Green Version]
  8. Song, G.; Lin, J.Y.; Bilheux, J.C.; Xie, Q.; Santodonato, L.J.; Molaison, J.J.; Skorpenske, H.D.; M Dos Santos, A.; Tulk, C.A.; An, K.; et al. Characterization of Crystallographic Structures Using Bragg-Edge Neutron Imaging at the Spallation Neutron Source. J. Imaging 2017, 3, 65. [Google Scholar] [CrossRef] [Green Version]
  9. Cereser, A.; Strobl, M.; Hall, S.A.; Steuwer, A.; Kiyanagi, R.; Tremsin, A.S.; Knudsen, E.B.; Shinohara, T.; Willendrup, P.K.; Fanta, A.B.S.; et al. Time-of-flight three dimensional neutron diffraction in transmission mode for mapping crystal grain structures. Sci. Rep. 2017, 7, 9561. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Vogel, S.C.; Bourke, M.A.; Craft, A.E.; Harp, J.M.; Kelsey, C.T.; Lin, J.; Long, A.M.; Losko, A.S.; Hosemann, P.; McClellan, K.J.; et al. Advanced Postirradiation Characterization of Nuclear Fuels Using Pulsed Neutrons. JOM 2020, 72, 187–196. [Google Scholar] [CrossRef]
  11. Mohan, K.; Venkatakrishnan, S.; Gibbs, J.; Gulsoy, E.; Xiao, X.; De Graef, M.; Voorhees, P.; Bouman, C. TIMBIR: A Method for Time-Space Reconstruction from Interlaced Views. Comput. Imaging IEEE Trans. 2015, 1, 96–111. [Google Scholar] [CrossRef]
  12. Kohler, T. A projection access scheme for iterative reconstruction based on the golden section. In Proceedings of the IEEE Symposium Conference Record Nuclear Science 2004, Rome, Italy, 16–22 October 2004; Volume 6, pp. 3961–3965. [Google Scholar]
  13. Kaestner, A.P.; Munch, B.; Trtik, P. Spatiotemporal computed tomography of dynamic processes. Opt. Eng. 2011, 50, 123201. [Google Scholar]
  14. Kaestner, A.; Lehmann, E.; Hovind, J.; Radebe, M.; De Beer, F.; Sim, C. Verifying neutron tomography performance using test objects. Phys. Procedia 2013, 43, 128–137. [Google Scholar] [CrossRef] [Green Version]
  15. Kaestner, A.P.; Trtik, P.; Zarebanadkouki, M.; Kazantsev, D.; Snehota, M.; Dobson, K.J.; Lehmann, E.H. Recent developments in neutron imaging with applications for porous media research. Solid Earth 2016, 7, 1281–1292. [Google Scholar] [CrossRef] [Green Version]
  16. Abir, M.; Islam, F.; Wachs, D.; Lee, H.K. Sparse-view neutron CT reconstruction of irradiated fuel assembly using total variation minimization with Poisson statistics. J. Radioanal. Nucl. Chem. 2016, 307, 1967–1979. [Google Scholar] [CrossRef]
  17. Kazantsev, D.; Bleichrodt, F.; van Leeuwen, T.; Kaestner, A.; Withers, P.J.; Batenburg, K.J.; Lee, P.D. A novel tomographic reconstruction method based on the robust Student’s t function for suppressing data outliers. IEEE Trans. Comput. Imaging 2017, 3, 682–693. [Google Scholar] [CrossRef] [Green Version]
  18. Venkatakrishnan, S.; Cakmak, E.; Billheux, H.; Bingham, P.; Archibald, R.K. Model-based iterative reconstruction for neutron laminography. In Proceedings of the 2017 51st Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, CA, USA, 29 October–1 November 2017; pp. 1864–1869. [Google Scholar]
  19. Barnard, R.C.; Bilheux, H.; Toops, T.; Nafziger, E.; Finney, C.; Splitter, D.; Archibald, R. Total variation-based neutron computed tomography. Rev. Sci. Instruments 2018, 89, 053704. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  20. Cramblitt, R.M.; Allebach, J.P. Analysis of time-sequential sampling with a spatially hexagonal lattice. JOSA 1983, 73, 1510–1517. [Google Scholar] [CrossRef]
  21. Bouman, C.A. Model Based Image Processing. 2013. Available online: https://engineering.purdue.edu/~bouman/publications/pdf/MBIP-book.pdf (accessed on 30 October 2020).
  22. Saghi, Z.; Holland, D.J.; Leary, R.; Falqui, A.; Bertoni, G.; Sederman, A.J.; Gladden, L.F.; Midgley, P.A. Three-Dimensional Morphology of Iron Oxide Nanoparticles with Reactive Concave Surfaces. A Compressed Sensing-Electron Tomography (CS-ET) Approach. Nano Lett. 2011, 11, 4666–4673. [Google Scholar] [CrossRef]
  23. Goris, B.; den Broek, W.V.; Batenburg, K.; Mezerji, H.H.; Bals, S. Electron tomography based on a total variation minimization reconstruction technique. Ultramicroscopy 2012, 113, 120–130. [Google Scholar] [CrossRef]
  24. Venkatakrishnan, S.; Drummy, L.; Jackson, M.; De Graef, M.; Simmons, J.; Bouman, C. Model based iterative reconstruction for Bright-Field Electron Tomography. IEEE Trans. Comput. Imaging 2015, 1, 1–15. [Google Scholar] [CrossRef]
  25. Gürsoy, D.; Biçer, T.; Lanzirotti, A.; Newville, M.G.; De Carlo, F. Hyperspectral image reconstruction for x-ray fluorescence tomography. Opt. Express 2015, 23, 9014–9023. [Google Scholar] [CrossRef] [PubMed]
  26. Mohan, K.A.; Venkatakrishnan, S.V.; Drummy, L.F.; Simmons, J.; Parkinson, D.Y.; Bouman, C.A. Model-Based Iterative Reconstruction for Synchrotron X-ray Tomography. In Proceedings of the 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Florence, Italy, 4–9 May 2014. [Google Scholar]
  27. Sauer, K.; Bouman, C. Bayesian Estimation of Transmission Tomograms Using Segmentation Based Optimization. IEEE Trans. Nucl. Sci. 1992, 39, 1144–1152. [Google Scholar] [CrossRef]
  28. Thibault, J.B.; Sauer, K.; Bouman, C.; Hsieh, J. A three-dimensional statistical approach to improved image quality for multislice helical CT. Med. Phys. 2007, 34, 4526–4544. [Google Scholar] [CrossRef] [PubMed]
  29. Golbabaee, M.; Vandergheynst, P. Hyperspectral image compressed sensing via low-rank and joint-sparse matrix recovery. In Proceedings of the 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Kyoto, Japan, 25–30 March 2012; pp. 2741–2744. [Google Scholar]
  30. Zhuang, L.; Bioucas-Dias, J.M. Fast hyperspectral image denoising and inpainting based on low-rank and sparse representations. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2018, 11, 730–742. [Google Scholar] [CrossRef]
  31. Monier, E.; Oberlin, T.; Brun, N.; Tencé, M.; de Frutos, M.; Dobigeon, N. Reconstruction of Partially Sampled Multiband Images—Application to STEM-EELS Imaging. IEEE Trans. Comput. Imaging 2018, 4, 585–598. [Google Scholar] [CrossRef] [Green Version]
  32. Kim, D.; Fessler, J.A. An optimized first-order method for image restoration. In Proceedings of the 2015 IEEE International Conference on Image Processing (ICIP), Quebec City, QC, Canada, 27–30 September 2015; pp. 3675–3679. [Google Scholar] [CrossRef]
  33. Palenstijn, W.; Batenburg, K.; Sijbers, J. Performance improvements for iterative electron tomography reconstruction using graphics processing units (GPUs). J. Struct. Biol. 2011, 176, 250–253. [Google Scholar] [CrossRef] [Green Version]
  34. van Aarle, W.; Palenstijn, W.J.; Beenhouwer, J.D.; Altantzis, T.; Bals, S.; Batenburg, K.J.; Sijbers, J. The ASTRA Toolbox: A platform for advanced algorithm development in electron tomography. Ultramicroscopy 2015, 157, 35–47. [Google Scholar] [CrossRef] [Green Version]
  35. Wang, X.; Sridhar, V.; Ronaghi, Z.; Thomas, R.; Deslippe, J.; Parkinson, D.; Buzzard, G.T.; Midkiff, S.P.; Bouman, C.A.; Warfield, S.K. Consensus equilibrium framework for super-resolution and extreme-scale CT reconstruction. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, Denver, CO, USA, 9–19 November 2019; pp. 1–23. [Google Scholar]
  36. Marchesini, S.; Trivedi, A.; Enfedaque, P.; Perciano, T.; Parkinson, D. Sparse Matrix-Based HPC Tomography. arXiv 2020, arXiv:2003.12677. [Google Scholar]
  37. Venkatakrishnan, S.V. pyMBIR; Technical Report; Oak Ridge National Lab. (ORNL): Oak Ridge, TN, USA, 2019. [Google Scholar]
  38. Gürsoy, D.; De Carlo, F.; Xiao, X.; Jacobsen, C. TomoPy: A framework for the analysis of synchrotron tomographic data. J. Synchrotron Radiat. 2014, 21, 1188–1193. [Google Scholar] [CrossRef] [Green Version]
  39. Wang, Z.; Bovik, A.C.; Sheikh, H.R.; Simoncelli, E.P. Image quality assessment: From error visibility to structural similarity. IEEE Trans. Image Process. 2004, 13, 600–612. [Google Scholar] [CrossRef] [Green Version]
  40. Batenburg, K.J.; Palenstijn, W.J.; Balázs, P.; Sijbers, J. Dynamic angle selection in binary tomography. Comput. Vis. Image Underst. 2013, 117, 306–318. [Google Scholar] [CrossRef] [Green Version]
Figure 1. A block diagram of the wavelength-resolved neutron tomography setup at a spallation neutron source. An incident hyper-spectral neutron beam, depicted here with multiple color arrows each corresponding to a neutron wavelength, provides the transmission signal through a sample placed in front of the detector. A time-of-flight (TOF) imaging detector measures the number of neutrons detected as a function of time, which is equivalent to wavelength. These measurements can be binned into wavelengths to obtain a projection measurement corresponding to each wavelength and rotation angle.
Figure 1. A block diagram of the wavelength-resolved neutron tomography setup at a spallation neutron source. An incident hyper-spectral neutron beam, depicted here with multiple color arrows each corresponding to a neutron wavelength, provides the transmission signal through a sample placed in front of the detector. A time-of-flight (TOF) imaging detector measures the number of neutrons detected as a function of time, which is equivalent to wavelength. These measurements can be binned into wavelengths to obtain a projection measurement corresponding to each wavelength and rotation angle.
Jimaging 07 00010 g001
Figure 2. Illustration of two ways to acquire tomographic datasets. (a) Representation of the conventional scanning protocol using 6 equally spaced angles between 0 and 180 . (b) Example of one realization of the proposed interlaced scanning protocol that ensures subsequent measurements are far apart ( K = 2 , N θ = 6 in Equation (1)). List of view angles using (c) the conventional approach and (d) one realization of the interlaced scan that was used in the experiments in this paper ( K = 3 , N θ = 36 ). The advantage of the interlaced scan is that at any given point during the course of the experiment, the object is scanned fairly uniformly, thereby enabling reconstructions with fewer artifacts that the end users can utilize to make decisions. The parameters of the interlaced scan can be further adjusted so that subsequent measurements are maximally separated (not shown here).
Figure 2. Illustration of two ways to acquire tomographic datasets. (a) Representation of the conventional scanning protocol using 6 equally spaced angles between 0 and 180 . (b) Example of one realization of the proposed interlaced scanning protocol that ensures subsequent measurements are far apart ( K = 2 , N θ = 6 in Equation (1)). List of view angles using (c) the conventional approach and (d) one realization of the interlaced scan that was used in the experiments in this paper ( K = 3 , N θ = 36 ). The advantage of the interlaced scan is that at any given point during the course of the experiment, the object is scanned fairly uniformly, thereby enabling reconstructions with fewer artifacts that the end users can utilize to make decisions. The parameters of the interlaced scan can be further adjusted so that subsequent measurements are maximally separated (not shown here).
Jimaging 07 00010 g002
Figure 3. Illustration of the experimental set-up. Two magnetite crystals are mounted on a holder, separated by a piece of cadmium and held by a tape. Each projection image consisted of measurements using a 512 × 512 micro-channel plate detector which natively binned the images into energy-bins. Each such projection measurement took approximately 1.5 h of experiment time and the counts per bin in a region of the sample were approximately in the range of 0 to 35 counts as shown in the plot above (average signal strength in a region indicated by the red box).
Figure 3. Illustration of the experimental set-up. Two magnetite crystals are mounted on a holder, separated by a piece of cadmium and held by a tape. Each projection image consisted of measurements using a 512 × 512 micro-channel plate detector which natively binned the images into energy-bins. Each such projection measurement took approximately 1.5 h of experiment time and the counts per bin in a region of the sample were approximately in the range of 0 to 35 counts as shown in the plot above (average signal strength in a region indicated by the red box).
Jimaging 07 00010 g003
Figure 4. Cross section from the reconstruction of a single wavelength bin using filtered back projection (FBP) and model-based image reconstruction (MBIR) with different scanning protocols. In each inset, the top row is the FBP reconstruction, the bottom is the MBIR. The left column is the result from the conventional scan, while the right column is from the interlaced scanning strategy. Notice that the MBIR reconstructions are of significantly higher visual quality than FBP. Furthermore, the interlaced scanning provides a clearer image of the sample in a more uniform manner thereby providing useful feedback very early in the scan compared to existing protocols. All images are displayed in the viewing window 0.001 , 0.0035 .
Figure 4. Cross section from the reconstruction of a single wavelength bin using filtered back projection (FBP) and model-based image reconstruction (MBIR) with different scanning protocols. In each inset, the top row is the FBP reconstruction, the bottom is the MBIR. The left column is the result from the conventional scan, while the right column is from the interlaced scanning strategy. Notice that the MBIR reconstructions are of significantly higher visual quality than FBP. Furthermore, the interlaced scanning provides a clearer image of the sample in a more uniform manner thereby providing useful feedback very early in the scan compared to existing protocols. All images are displayed in the viewing window 0.001 , 0.0035 .
Jimaging 07 00010 g004
Figure 5. Comparison of normalized root mean squared error (NRMSE) and structural-similarity (SSIM) metrics between the current reconstruction and the final reconstruction in Figure 4 using the MBIR algorithm. Notice that the SSIM and NRMSE converge faster for the interlaced scanning approach compared to the conventional scanning protocol.
Figure 5. Comparison of normalized root mean squared error (NRMSE) and structural-similarity (SSIM) metrics between the current reconstruction and the final reconstruction in Figure 4 using the MBIR algorithm. Notice that the SSIM and NRMSE converge faster for the interlaced scanning approach compared to the conventional scanning protocol.
Jimaging 07 00010 g005
Figure 6. XY, XZ and YZ cross sections of the final reconstruction corresponding to the first wavelength bin. Notice that the features of the sample are clearly visible in MBIR even though the input data has only 30 projections. The filter and regularization parameters are chosen for approximately similar visual resolution. Notice that the MBIR approach results in a significantly lower noise reconstruction compared to the FBP approach.
Figure 6. XY, XZ and YZ cross sections of the final reconstruction corresponding to the first wavelength bin. Notice that the features of the sample are clearly visible in MBIR even though the input data has only 30 projections. The filter and regularization parameters are chosen for approximately similar visual resolution. Notice that the MBIR approach results in a significantly lower noise reconstruction compared to the FBP approach.
Jimaging 07 00010 g006
Figure 7. Line profile through the center of the rock. Notice that the FBP produces a very noisy reconstruction as a function of wavelength index.
Figure 7. Line profile through the center of the rock. Notice that the FBP produces a very noisy reconstruction as a function of wavelength index.
Jimaging 07 00010 g007
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Venkatakrishnan, S.; Zhang, Y.; Dessieux, L.; Hoffmann, C.; Bingham, P.; Bilheux, H. Improved Acquisition and Reconstruction for Wavelength-Resolved Neutron Tomography. J. Imaging 2021, 7, 10. https://0-doi-org.brum.beds.ac.uk/10.3390/jimaging7010010

AMA Style

Venkatakrishnan S, Zhang Y, Dessieux L, Hoffmann C, Bingham P, Bilheux H. Improved Acquisition and Reconstruction for Wavelength-Resolved Neutron Tomography. Journal of Imaging. 2021; 7(1):10. https://0-doi-org.brum.beds.ac.uk/10.3390/jimaging7010010

Chicago/Turabian Style

Venkatakrishnan, Singanallur, Yuxuan Zhang, Luc Dessieux, Christina Hoffmann, Philip Bingham, and Hassina Bilheux. 2021. "Improved Acquisition and Reconstruction for Wavelength-Resolved Neutron Tomography" Journal of Imaging 7, no. 1: 10. https://0-doi-org.brum.beds.ac.uk/10.3390/jimaging7010010

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop