Neutron source reconstruction from pinhole imaging at National Ignition Facility P. Volegov, C. R. Danly, D. N. Fittinghoff, G. P. Grim, N. Guler, N. Izumi, T. Ma, F. E. Merrill, A. L. Warrick, C. H. Wilde, and D. C. Wilson Citation: Review of Scientific Instruments 85, 023508 (2014); doi: 10.1063/1.4865456 View online: http://dx.doi.org/10.1063/1.4865456 View Table of Contents: http://scitation.aip.org/content/aip/journal/rsi/85/2?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Iterative reconstruction of Fourier-rebinned PET data using sinogram blurring function estimated from point source scans Med. Phys. 37, 5530 (2010); 10.1118/1.3490711 Noise and resolution in images reconstructed with FBP and OSC algorithms for CT Med. Phys. 34, 585 (2007); 10.1118/1.2409481 Neutron source images and spectra to guide neutron diagnostics specifications for the National Ignition Facility Rev. Sci. Instrum. 77, 10E722 (2006); 10.1063/1.2351885 Simultaneous iterative reconstruction of emission and attenuation images in positron emission tomography from emission data only Med. Phys. 29, 1962 (2002); 10.1118/1.1500400 Plasma ignition schemes for the Spallation Neutron Source radio-frequency driven H − source Rev. Sci. Instrum. 73, 1017 (2002); 10.1063/1.1430521
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
REVIEW OF SCIENTIFIC INSTRUMENTS 85, 023508 (2014)
Neutron source reconstruction from pinhole imaging at National Ignition Facility P. Volegov,1 C. R. Danly,1 D. N. Fittinghoff,2 G. P. Grim,1 N. Guler,1 N. Izumi,2 T. Ma,2 F. E. Merrill,1 A. L. Warrick,2 C. H. Wilde,1 and D. C. Wilson1 1 2
Los Alamos National Laboratory, Los Alamos, New Mexico 87544, USA Livermore National Laboratory, Livermore, California 94550, USA
(Received 9 October 2013; accepted 30 January 2014; published online 18 February 2014) The neutron imaging system at the National Ignition Facility (NIF) is an important diagnostic tool for measuring the two-dimensional size and shape of the neutrons produced in the burning deuterium-tritium plasma during the ignition stage of inertial confinement fusion (ICF) implosions at NIF. Since the neutron source is small (∼100 μm) and neutrons are deeply penetrating (>3 cm) in all materials, the apertures used to achieve the desired 10-μm resolution are 20-cm long, single-sided tapers in gold. These apertures, which have triangular cross sections, produce distortions in the image, and the extended nature of the pinhole results in a non-stationary or spatially varying point spread function across the pinhole field of view. In this work, we have used iterative Maximum Likelihood techniques to remove the non-stationary distortions introduced by the aperture to reconstruct the underlying neutron source distributions. We present the detailed algorithms used for these reconstructions, the stopping criteria used and reconstructed sources from data collected at NIF with a discussion of the neutron imaging performance in light of other diagnostics. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4865456] I. INTRODUCTION
The National Ignition Facility (NIF) team is currently attempting to achieve thermonuclear ignition of a deuteriumtritium (DT) plasma.1 The strategy being used is to generate the required densities and pressures through indirect ablative compression of plastic (CH) capsules containing DT ice layers, which surround DT gas. In this design, implosion process forms a region of less dense fuel surrounded by higher density DT fuel. The inner region reaches temperatures of ∼10 keV. All of the DT reaches an areal density, ρR, of greater than ∼1 g/cm2 . Fusion reactions and resulting energy deposition from alpha particles further heat this hot spot which then initiates an ignition front that propagates through the surrounding “cold” high density DT fuel. This reaction could release over 106 joules of fusion energy and more than 1017 neutrons. Tuning the implosion to achieve ignition requires the measurement of plasma conditions at very small length and time scales. The hot spot is predicted to have a 30–50 μm diameter, while the surrounding cold fuel is expected to be ∼100 μm in diameter, and the fusion conditions must be maintained for ∼100 ps.1 Establishing the appropriate conditions for ignition, such as symmetry, hot-spot volume, and cold-fuel volume, depends on iterative tuning, and until ignition is achieved, measurements must be made at non-ignition conditions. A neutron imaging diagnostic, known as the Neutron Imaging System (NIS), has been built at NIF and is being used to provide data on the size and shape of the fusion hotspot and the surrounding cold fuel during the initial tuning campaigns.2 This aperture imaging system is composed of two basic pieces: an aperture array that is used to form the neutron images and a scintillator-based detector system that 0034-6748/2014/85(2)/023508/12/$30.00
is used to measure the neutron flux passing through the aperture array. The BCF-99-55 scintillating fiber bundle is viewed with two MCP gated image collection systems, which allows the collection of two independently timed images, providing images from two neutron energy intervals. The detector system is capable of collecting two fastgated and independently timed images. Because the detector array is positioned 28 m from the neutron source, the neutron arrival time at the detector is correlated to the neutron energies, allowing measurement of the neutron source distributions from two energy ranges by gating the detectors at two different times. Typically, one detector is gated from 13 to 17 MeV to view the 14-MeV neutrons that are generated from DT fusion processes and provides information on the size and shape of the hot spot. The second detector is gated to measure the source distribution of lower energy neutrons, typically in the range from 6 to 12 MeV. These lower energy neutrons are predominantly DT-fusion neutrons that have scattered in the surrounding cold fuel and, therefore, provide information on the distribution of the cold fuel surrounding the hotspot. Because the scintillation light from BCF-99-55 decays quickly (∼3 ns) the light from the primary neutrons is essentially gone by the time the scattered neutrons arrive at the scintillator, allowing the two images to be measured with the same scintillator. Since the images are produced in the same scintillator the relative positions of the hot spot and cold fuel images have been determined by collecting the two images simultaneously with a calibration experiment. This relative image registration is accurate to within 5 μm in the source plane. We note that the absolute registration relative to the original source position is not known accurately due to alignment precision of ∼100 μm in the source plane. However, we are
85, 023508-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-2
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
able to extract this information from the data as described in Sec. III. Because of the long mean free path of neutrons within aperture materials, the neutron imaging system forms images using extended apertures or “pinholes” in materials with short interaction lengths. The aperture array currently in use at NIF is composed of triangular pinholes and mini-penumbral apertures machined along the 20 cm length of wedged gold layers. Gold has a mean free path for elastic scattering of ∼3 cm for 14 MeV neutrons. These equilateral triangular pinholes are 226-μm high at the downstream end of the aperture array and taper to 5-μm high at the surface facing the neutron source. The triangular shape of the pinholes results in image distortion while the extended nature of the pinhole results in a nonstationary point spread function generating distortions which vary across the pinhole field of view.3–8 To provide the size and shape information required by the tuning campaigns at the National Ignition Facility, two Maximum-Likelihood reconstruction algorithms have been developed to remove these distortions. In this work, we present details of these algorithms, the stopping criteria used and present reconstructed neutron sources from NIF measurements along with a discussion of the performance in light of other diagnostics. In Sec. II, we present an overview of the NIS system and give a mathematical description of the formation of images with non-stationary point-spread functions. In Sec. III, we describe the employed image processing techniques: the data reduction prior to reconstruction, determination of the aperture array pointing, and the iterative algorithm that we have used to reconstruct the source images including discussion of the stopping criterion. The results of the reconstructions and conclusions will then be discussed in Sec. IV. II. NEUTRON IMAGING SYSTEM OVERVIEW
As shown in Figure 1, the NIS aperture array has 20 triangular pinhole apertures and three mini-penumbral apertures.2 The array was designed so that the fields of view of all apertures overlap at 26.5 cm from the front face. In its usage for the inertial confinement fusion (ICF) implosions at the NIF, the front surface of the aperture array is positioned at 32.5 cm from the source. This results in each aperture pointing to a different location at the source plane, resulting in a combined 0.5 mm field of view. A. Image formation
For a neutron aperture imager, neglecting the neutrons that scatter in the aperture material or secondary scattering of neutrons in the detector, the mean neutron flux at the detection plane can be written as ˆ (nˆ · ) r ˆ 3r , S(r , E, )d F (r, E) = e− ∫r μt (ξ ,E)|dξ | 4π |r − r |2 S
(1) where F(r, E) is the neutron flux at the detection plane, i.e., F(r, E)AE is a number of neutrons crossed differential area A of detection plane about r, with differential energies E energy about E; μt (ξ , E) is the total linear attenuation co-
FIG. 1. Photograph of the rear-surface of the NIS aperture array. The array has 20 pinhole apertures with triangular cross sections arranged in four rows, with three penumbral apertures in the center row. The triangular apertures shown here are 226-μm high and taper to 5-μm high at the front surface of the array, which faces the source.
efficient for neutrons with energy E at the point ξ ; nˆ is a unit ˆ ≡ (r, ˆ r ) ≡ r−r vector normal to the detection surface; |r−r | is a unit vector in the direction from the point r to the ˆ is the differential neutron source, i.e., point r; S r , E, ˆ V E is a number of neutrons emitted in S r, E, the volume V about the point r , with differential energies E energy about E, in the solid angle about direction deˆ and integral over S is an integral fined by the unit vector ; over the source volume. In Eq. (1) and after we use the following conventions: bold font denotes vectors and matrices, italic denotes scalars, a hat sign denotes a unit vector, a capital italic superscript T denotes a transposed matrix/vector, and the vectors are column vectors. In general, a linear, position-sensitive detector can be described by a spatial sensitivity function, also called a point spread function of a detector, P(r, r , E), such that P(r, r , E)AE is the signal generated in the sensor at the position r by a unit neutron flux through the differential area A of the detection plane about r , with differential energies E about a central energy, E. In this case, we can write the recorded image as ∞ I (r) = 0
ˆ 3 r dE, K(r, r , E)S(r , E, )d
(2)
S
where the total sensitivity function of the neutron imaging system is defined by the following formula: ˆ (nˆ · ) r K(r, r , E) = P (r, r , E)e− ∫r μt (ξ ,E)|dξ | d 3 r . 4π |r −r |2 S
(3)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
Volegov et al.
023508-3
Rev. Sci. Instrum. 85, 023508 (2014)
Here the integration is over the detector sensitive volume, for example, a scintillator surface. It should be noted that the function K(r, r , E) models the blurring of the imaging system and is called the point spread function (PSF) of the image collection system. The last two equations can be simplified assuming that the extent of the neutron source and the detector are small compared to the distance between them. In this case we can reduce the dimensionality of the problem: ∞ I (u) =
K(u, v, E)S(v, E) d 2 vdE. 0
(4)
S
Here the vectors u = r − rsens and v = r − rsrc define a position on the detector plane and the source plane, respectively, and S (v, E) is a projection of the neutron source on the plane orthogonal to the axis connecting the center of the source, rsrc , and the center of the detector, rsens , which is the central axis of the imaging system: ∞ S(v, E) =
ˆ 0 l, E)dl, S(r +
(5)
−∞ rsens −rsrc . |rsens −rsrc |
ˆ0 ≡ Assuming that both the where the unit vector detector plane and the source plane are orthogonal to the central axis of the imaging system, the sensitivity function can be written as 1 u P (u, u , E)e− ∫v μt (ξ ,E)|dξ | d 2 u . K(u, v, E) = 2 4π L Ds
(6) Here Ds denotes the integration over the detector plane, L = |rsens − rsrc | is the distance between the source center and the detector plane. This approximation is called the paraxial approximation or a small-angle approximation. Assuming that the energy distribution of the source is independent of the position so that the function S (v, E) is separable and S (v, E) = S (v) U (E), we can further simplify the imaging equation: (7) I (u) = K(u, v)S(v)d 2 v, S
where energy distribution function U(E) is considered to be known. In practice, using time-of-flight techniques, the images are recorded in several energy intervals, assuming that spatial and energy distributions of the neutrons are separable inside each interval. The kernel function K(u, v) in this case is defined by K(u, v) ∞ 1 u U (E) P (u, u , E)e− ∫v μt (ξ ,E)|dξ | d 2 u dE. = 2 4π L 0
Ds
(8)
FIG. 2. A combination of images, i.e., point-spread functions, of point sources displaced in x and y in steps of 100-μm from the central axis of a NIF triangular pinhole. As it can be seen, the shape of a point-spread function strongly depends on position of the source relative to the central axis of an aperture.
Typically a point spread function of a detector P(u, u , E) is modeled by a Gaussian: P (u, u , E) = √
1 2π σ
e−
|u−u |2 2σ 2
,
(9)
where σ is the RMS of the Gaussian blur. The spatial dependence of the point-spread function for a triangular pinhole aperture used at NIF is illustrated in Figure 2, which shows an image corresponding to a combination of point sources, each located on a 100-μm grid. The equations that have been considered so far connecting the neutron source with the recorded image are idealized imaging equations because they do not take into account various imperfections of the imaging system, such as the contribution of the scattered neutrons or imperfections in the detection system. III. IMAGE PROCESSING A. Data reduction
Data collected from the two imaging systems are normalized before source reconstructions are performed. In the process of collecting data for a NIF experiment, three dark field images are collected to measure the baseline levels of the CCD camera system. These dark fields are averaged together and subtracted from the data images. During image collection some neutrons are scattered from the flight path and interact directly with elements of the CCD array. Each direct interaction generates a few bright pixels, often called “stars.” After the dark fields are subtracted these stars are differentiated from statistical deviations of the local signal levels and set to the average of the neighboring pixels. At the present yields (∼1014 neutrons) less than 0.1% of the pixels are affected by this procedure.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-4
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
FIG. 3. The image produced by the aperture array is shown with the same contrast range before (left) and after (right) the scatter background subtraction. The images are 521 by 521 pixels with each pixel being 341 μm.
Once the average dark field is subtracted and stars are removed, the images are divided by flat-field images.9 The recording system has been calibrated and is operated in a linear regime both during pinhole imaging and flat-field collection. The flat-field images are measurements of the light generated by a uniform distribution of neutrons interacting with the scintillating fiber array. By dividing the data images by the flat fields, the variation in response across the scintillator and image collection systems is removed. The image is then typically re-binned from the substantial pixel oversampling of 37 μm pixels at the scintillator plane to 341 μm pixels, which maps to ∼4 μm pixels at the source location. The remainder of the analysis has been typically performed at this pixel size. Of the 40% of neutrons that interact in the scintillator, some small fraction, ∼10%, interact a second time before exiting the scintillator material. This results in a long-range scatter background, which extends to more than 3 cm. Because of the extent of this long-range scatter background, it is of little consequence to the resolution for a single pinhole, where peak scatter background is 0, pi, j > 0, xj(0) > 0 and ensures that the solution is also positive, which is consistent with Poisson statistics.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-8
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
For the Gaussian noise model, the likelihood function is defined by L(S) = (2π )− 2 |Cεε |− 2 e− 2 (I−KS) m
1
1
T
C−1 εε (I−KS)
,
(21)
where Cεε is a m × m covariance matrix of a random vector ε, |Cεε | denotes the determinant of the covariance matrix, C−1 εε is inverse of the covariance matrix, and m is the total number of pixels in the recorded image I. The log-likelihood function in this case is 1 l(S) = − (I − KS)T C−1 εε (I − KS) + const. 2
(22)
Using this likelihood function the MLE solution of the problem (13) is simply −1 T −1 K Cεε I, (23) SMLE = KT C−1 εε K which is exactly a weighted least squares solution. If the noise is uncorrelated, i.e., the image pixels are independent the matrix Cεε is a diagonal matrix with the diagonal elements σi2 , i = 1, .., m being a variance of the noise in the ith pixel. Direct application of Eq. (23) in the context of reconstruction of neutron images is somewhat impractical due to difficulties in inverting a large and ill-conditioned matrix KT C−1 εε K. An iteration scheme, which we will refer to as MLGAUSS, similar to algorithm for the Poisson noise described above, see Eq. (19), was developed for Gaussian noise distributions:15 m yi − y˜i(t) (t+1) (t) xj = xj 1 + h pi,j , σi2 i=1 (24) j = 1, .., n. The step size h is computed from 1D maximization of the likelihood function, Eq. (22). The details of the step size computation are described in Ref. 15. It should be noted that the formula (24) does not require that quantities yi , pi, j , and xj(t) be positive, and does not ensure that the solution is positive. In the case of Poisson noise, where a variance of the noise in a pixel is equal to the mean value of the signal in this pixel, i.e., σi2 = y˜i , i = 1, ..., m, and assuming the step size h ≡ 1, formula (24) corresponds to the iteration scheme shown in Eq. (19). An important question for practical application of the iteration scheme (24), or (19), is when to stop the iterations. It is known that when the described algorithm passes a certain point the image starts to deteriorate. This optimal point of truncation depends on statistics (noise): the higher the statistics, the later in the iteration the degradation begins. The number of iterations, therefore, becomes the regularization parameter. The goal then is to stop the iterations as soon as the solution is consistent with the data within the experimental accuracy of the measurement. The agreement of the solution with the data is usually characterized by normalized chi-square (or reduced chi-squared) criteria, defined by the following formula: 1 yi − y˜i . m i=1 σi2 m
χn2 =
(25)
For a solution that is consistent with the data, the expectation value of χn2 is 1, and the variance is 2/m. Thus the observed value of χn2 is easy to interpret: if χn2 is about 1 within, say, a variance, then our solution is consistent with the data, if χn2 1, then the solution captures not only the neutron distribution but also the noise (or the noise is overestimated), and finally if χn2 1, then the solution is not consistent with the data (or the noise is underestimated). These simple rules may be formalized by introducing P- and Q-values. Given a particular value, we can calculate the probability to obtain χn2 that is less than the P-value and the probability to obtain χn2 that is greater than the Q-Value. The most widely accepted critical value is Q = 0.05, meaning that only in 5% of the trials would result in a higher value of χn2 . This leads to the simple stopping rule: continue iterations until the Q-Value is higher than the threshold value. To apply the stopping rule described in the previous paragraph, the variance of the noise should be reliably estimated. In practice, this is not always possible. An alternative is to compute the relative change of the residual from one iteration to the next iteration: ρ ≡ (ρ (t) − ρ (t + 1) /ρ (t) ), where ρ (t) is some measure of the difference, and stop the iterations as soon as this relative difference is below some value. Typically, reduced chi-square is used as a measure of the difference: 1 yi − y˜i(t) . m i=1 σi2 m
ρ (t) =
(26)
However, one may employ a different measure such as the Kullback-Leibler divergence: m yi (t) ρ = yi ln . (27) y˜i(t) i=1 Empirically it was found that critical value 10−3 provides a robust stopping criterion in most of the cases. More complicated and sophisticated stopping criteria exist; see, for example, Refs. 15–18. These criteria will be investigated in future work. IV. APPLICATION OF MLE ALGORITHMS TO DATA
Many neutron images have been collected during the ignition tuning experiments at NIF over the last three years. The neutron source distributions have each been reconstructed with the algorithms described above. Typically these reconstructed sources are 50–100 μm elliptical sources with very little features. This type of reconstruction has proved little opportunity to assess the performance of these reconstruction algorithms. Recently, however, a new set of experiments have been performed to study hydro-dynamic implosion physics at NIF. These experiments are designed to be driven with a stronger shock and thus implode along a higher adiabat of the thermonuclear fuel, resulting in a larger stagnation radius and larger neutron images. In addition, the initial experiments, which were executed to tune the shape of the implosion, resulted in a source with a shape more complicated than the typical simple ellipsoid, and with record neutron yields. This more complicated shape with high
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-9
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
FIG. 9. Primary neutron image (a) and down-scattered neutrons image (b) for a high-yield NIF shot (NIF shot number N130710-002-999).
statistics has provided a feature rich source for tests of these reconstruction algorithms. The raw images from this experiment (NIF shot number N130710-002-999) for the primary neutrons and for the down-scattered neutrons are shown in Figure 9. The reconstructed primary source distributions using two different algorithms (EMML and MLGAUSS) are shown in Figure 10. Convergence for both of these algorithms is illustrated by
the χ 2 evolutions with respect to the number of iterations shown in Figure 11. In spite of apparent difference in convergence behavior, after 50 iterations both algorithms demonstrated almost exactly the same convergence: EMML – χ 2 = 1.31, χ 2 /χ 2 = 0.8K × 10−3 , and MLGAUSS – χ 2 = 1.32, χ 2 /χ 2 = 1.2 × 10−3 . The reconstructed source distributions, as it can be seen from top raw images in Figure 10, are similar but differ in size: MLGAUSS produces bigger and
FIG. 10. Primary neutron image reconstruction of the source by using EMML (left: (a) 50 iterations and (c) 262 iterations) and MLGAUSS (right: (b) 50 iterations and (d) 262 iteration) algorithms (NIF shot number N130710-002-999).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-10
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
FIG. 13. Overlaid images filtered in red and cyan colors for primary and downscattered sources, respectively (NIF shot number N130710-002-999).
FIG. 11. Reduced χ 2 as a function of iteration (NIF shot number N130710002-999, primary source reconstruction).
smother source. In this example EMML algorithm changes χ 2 very little after 50 iterations, specifically after 262 iterations: χ 2 = 1.28, χ 2 /χ 2 = 3.5 × 10−5 , while perceptibly changing the reconstructed source (see Figures 10(a) and 10(c)) making it implausible from the physical considerations. MLGAUSS, however, exhibits an opposite behavior: the residual changes perceptibly – after 262 iterations: χ 2 = 1.19, χ 2 /χ 2 = 7.0 × 10−7 – while change of the reconstructed source is not that drastic. Explanation of such a behavior lays in the background noise. As it was mentioned in Sec. III A, dark field images are subtracted from the data images with the result that the noise in the background pixels is no longer consistent with a Poisson noise model: the pixels can take negative values, and generally the average of the background pixels is not necessarily a zero, it could be negative. Such an inconsistency, even with non-positive pixels being discarded, drives fast corruption of the reconstructed image using EMML algorithm which assumes a Poisson noise model. The MLGAUSS algorithm on the other hand does not require non-negativity of the data and thus exhibits more robust behavior with respect to the
background subtraction procedure. Summarizing, in selecting between those two algorithms – EMML and MLGAUSS – one must take into account that EMML requires that the data be all positive and ensures that the reconstructed image is also positive, which is a very strong regularization provision. On the other hand, MLGAUSS does not impose nonnegativity on the reconstructed image, but does not require non-negativity of the data and handles negative values from background subtraction, making this the preferred algorithm for this application. Application of the described algorithm to the low statistic data is illustrated using down-scattered neutrons image (see Figure 9(b)). Here intensity of the down-scattered image is about 20 times lower than the intensity of the primary image. However, even in this case, both algorithms reconstruct a plausible neutron source as illustrated in Figure 12, which shows the image and the reconstructed source for the down-scattered neutrons, again with the two algorithms. Note that reconstruction of low statistic data requires significantly less iterations: EMML – 12 iterations and MLGAUSS – 28 iterations. It is important to compare and validate the reconstructed source distributions against physical models of the process and relevant experimental results obtained using a different technique. Figure 13 shows in red the reconstructed
FIG. 12. Down-scattered neutron image reconstruction of the source by using EMML ((a) 12 iterations) and MLGAUSS ((b) 28 iterations) algorithms (NIF shot number N130710-002-999).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-11
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
demonstrating that these algorithms are effective in removing the unwanted image distortions. The algorithm based on Gaussian noise model – MLGAUSS – has been found to be tolerant to background subtraction and, therefore, used in these applications. Algorithms based on a Poisson noise model, such as EMML, maybe better for other applications. These algorithms are applicable to a wide range of imaging modalities; from pinhole imaging to penumbral imaging and we are investigating how this new analysis capability can be used in the future to enable new imaging systems with higher resolution at low statistics, such as neutron imaging with coded apertures.
ACKNOWLEDGMENTS
FIG. 14. The X-ray image polar view, the NIS LOS is shown with the arrow (NIF shot number N130710-002-999).
distribution of the primary fusion neutrons generated in the hot DT plasma overlaid on the distribution of the cold fuel in cyan, which has been reconstructed from the down scattered neutron image. The reconstructed primary image clearly shows a source distribution with two connected lobes, consistent with a toroid of burning nuclear fuel. The cold fuel distribution is consistent with this toroid of burning fuel with cold fuel penetrating through the vertical axis of the burning plasma. An x-ray image collected along the vertical axis shows the orthogonal view from the top of the toroid shown in Figure 14. One would expect this high energy x-ray image to be consistent with a toroidal primary neutron image, as only the hottest regions of the nuclear fuel will generate fusion neutrons and this region will also be the brightest emitter of x-rays.19, 20 A comparison of the raw image to the reconstructed source shows that a significant level of distortions has been removed to form the source neutron distribution. The fact that the resulting source reconstruction is consistent with the orthogonal x-ray image demonstrates that this method of reconstruction is effective for this application. V. CONCLUSIONS
Iterative Maximum Likelihood Estimation algorithms were used to develop image reconstruction tools for the neutron images that have been recently collected at NIF. Accurate pointing of the aperture axis to the neutron source has proven to be an important parameter in accurately reconstructing neutron source distributions. The array of multiple pinhole apertures has allowed us to extract accurate pointing information, which is necessary in the calculation of the proper non-stationary point spread function for reconstructions. These tools have been used to analyze dozens of neutron images. Most recently these algorithms have been used to reconstruct a toroidal neutron source, which required the removal of significant image distortions due to the pinhole aperture and detector system. This reconstructed source shape has been independently validated with an orthogonal x-ray view,
Additional credit goes to the dedicated staff and technicians of NIF, whose hard work and operational expertise resulted in the data that are shown here. The authors wish to acknowledge D. Jedlovec, M. A. Talison, O. Drury, D. Kalantar, and R. Wood for their hard work and expertise in fielding and aligning the neutron imaging system. This work has been performed under the auspices of the U.S. Department of Energy for NNSA Campaign 10 (Inertial Confinement Fusion) with Steve Batha as program manager. 1 S. H. Glenzer, D. A. Callahan, A. J. MacKinnon, J. L. Kline, G. Grim, E. T.
Alger, R. L. Berger, L. A. Bernstein, R. Betti, D. L. Bleuel, T. R. Boehly, D. K. Bradley, S. C. Burkhart, R. Burr, J. A. Caggiano, C. Castro, D. T. Casey, C. Choate, D. S. Clark, P. Celliers, C. J. Cerjan, G. W. Collins, E. L. Dewald, P. DiNicola, J. M. DiNicola, L. Divol, S. Dixit, T. Döppner, R. Dylla-Spears, E. Dzenitis, M. Eckart, G. Erbert, D. Farley, J. Fair, D. Fittinghoff, M. Frank, L. J. A. Frenje, S. Friedrich, M. Gatu Johnson, C. Gibson, E. Giraldez, V. Glebov, S. Glenn, N. Guler, S. W. Haan, B. J. Haid, B. A. Hammel, A. V. Hamza, C. A. Haynam, G. M. Heestand, M. Hermann, H. W. Hermann, D. G. Hicks, D. E. Hinkel, J. P. Holder, D. M. Holunda, J. B. Horner, W. W. Hsing, H. Huang, N. Izumi, M. Jackson, O. S. Jones, D. H. Kalantar, R. Kauffman, J. D. Kilkenny, R. K. Kirkwood, J. Klingmann, T. Kohut, J. P. Knauer, J. A. Koch, B. Kozioziemki, G. A. Kyrala, A. L. Kritcher, J. Kroll, K. La Fortune, L. Lagin, O. L. Landen, D. W. Larson, D. LaTray, R. J. Leeper, S. Le Pape, J. D. Lindl, R. Lowe-Webb, T. Ma, J. McNaney, A. G. MacPhee, T. N. Malsbury, E. Mapoles, C. D. Marshall, N. B. Meezan, F. Merrill, P. Michel, J. D. Moody, A. S. Moore, M. Moran, K. A. Moreno, D. H. Munro, B. R. Nathan, A. Nikroo, R. E. Olson, C. D. Orth, A. E. Pak, P. K. Patel, T. Parham, R. Petrasso, J. E. Ralph, H. Rinderknecht, S. P. Regan, H. F. Robey, J. S. Ross, M. D. Rosen, R. Sacks, J. D. Salmonson, R. Saunders, J. Sater, C. Sangster, M. B. Schneider, F. H. Séguin, M. J. Shaw, B. K. Spears, P. T. Springer, W. Stoeffl, L. J. Suter, C. A. Thomas, R. Tommasini, R. P. J. Town, C. Walters, S. Weaver, S. V. Weber, P. J. Wegner, P. K. Whitman, K. Widmann, C. C. Widmayer, C. H. Wilde, D. C. Wilson, B. Van Wonterghem, B. J. MacGowan, L. J. Atherton, M. J. Edwards, and E. I. Moses, Phys. Plasmas 19(5), 056318 (2012). 2 F. E. Merrill, D. Bower, R. Buckles, D. D. Clark, C. R. Danly, O. B. Drury, J. M. Dzenitis, V. E. Fatherley, D. N. Fittinghoff, R. Gallegos, G. P. Grim, N. Guler, E. N. Loomis, S. Lutz, R. M. Malone, D. D. Martinson, D. Mares, D. J. Morley, G. L. Morgan, J. A. Oertel, I. L. Tregillis, P. L. Volegov, P. B. Weiss, C. H. Wilde, and D. C. Wilson, Rev. Sci. Instrum. 83(10), 10D317 (2012). 3 D. C. Wilson, G. P. Grim, I. L. Tregillis, M. D. Wilke, M. V. Patel, S. M. Sepke, G. L. Morgan, R. Hatarik, E. N. Loomis, C. H. Wilde, J. A. Oertel, V. E. Fatherley, D. D. Clark, D. N. Fittinghoff, D. E. Bower, M. J. Schmitt, M. M. Marinak, D. H. Munro, F. E. Merrill, M. J. Moran, T. S. F. Wang, C. R. Danly, R. A. Hilko, S. H. Batha, M. Frank, and R. Buckles, Rev. Sci. Instrum. 81(10), 10D335 (2010). 4 E. Loomis, G. Grim, C. Wilde, D. C. Wilson, M. Wilke, J. Finch, G. Morgan, I. Tregillis, and D. Clark, Sixth International Conference on Inertial Fusion Sciences and Applications, Parts 1–4, Vol. 244, 2010.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-12 5 G.
Volegov et al.
P. Grim, R. D. Day, D. D. Clark, V. E. Fatherley, F. P. Garcia, S. A. Jaramillo, A. J. Montoya, G. L. Morgan, J. A. Oertel, T. A. Ortiz, J. R. Payton, P. D. Pazuchanics, D. A. W. Schmidt, A. C. Valdez, C. H. Wilde, and M. D. Wilke, Proc. SPIE 6707, H7070 (2007). 6 C. A. Barrera, E. C. Morse, and M. J. Moran, Rev. Sci. Instrum. 77(10), 10E716 (2006). 7 M. D. Wilke, S. H. Batha, P. A. Bradley, R. D. Day, D. D. Clark, V. E. Fatherley, J. P. Finch, R. A. Gallegos, F. P. Garcia, G. P. Grim, S. A. Jaramillo, A. J. Montoya, M. J. Moran, G. L. Morgan, J. A. Oertel, T. A. Ortiz, J. R. Payton, P. Pazuchanics, D. W. Schmidt, A. C. Valdez, C. H. Wilde, and D. C. Wilson, Rev. Sci. Instrum. 79(10), 10E529 (2008). 8 V. Y. Glebov, D. D. Meyerhofer, T. C. Sangster, C. Stoeckl, S. Roberts, C. A. Barrera, J. R. Celeste, C. J. Cerjan, L. S. Dauffy, D. C. Eder, R. L. Griffith, S. W. Haan, B. A. Hammel, S. P. Hatchett, N. Izumi, J. R. Kimbrough, J. A. Koch, O. L. Landen, R. A. Lerche, B. J. MacGowan, M. J. Moran, E. W. Ng, T. W. Phillips, P. M. Song, R. Tommasini, B. K. Young, S. E. Caldwell, G. P. Grim, S. C. Evans, J. M. Mack, T. J. Sedillo, M. D. Wilke, D. C. Wilson, C. S. Young, D. Casey, J. A. Frenje, C. K. Li, R. D. Petrasso, F. H. Seguin, J. L. Bourgade, L. Disdier, M. Houry, I. Lantuejoul, O. Landoas, G. A. Chandler, G. W. Cooper, R. J. Leeper, R. E. Olson, C. L. Ruiz, M. A. Sweeney, S. P. Padalino, C. Horsfield, and B. A. Davis, Rev. Sci. Instrum. 77(10), 10E715 (2006).
Rev. Sci. Instrum. 85, 023508 (2014) 9 N.
Guler, P. Volegov, C. R. Danly, G. P. Grim, F. E. Merrill, and C. H. Wilde, Rev. Sci. Instrum. 83(10), 10D316 (2012). 10 M. Z. Tarasko, Preprint of the Physics and Energy Institute, PEI-156 (1969). 11 A. J. Rockmore and A. Macovski, IEEE Trans. Nucl. Sci. 23(4), 1428–1432 (1976). 12 A. P. Dempster, N. M. Laird, and D. B. Rubin, J. R. Stat. Soc. Ser. B 39(1), 1–38 (1977). 13 W. H. Richardson, J. Opt. Soc. Am. 62(1), 55–59 (1972). 14 L. B. Lucy, Astron. J. 79(6), 745–754 (1974). 15 V. I. Gelfgat, E. L. Kosarev, and E. R. Podolyak, Comput. Phys. Commun. 74(3), 335–348 (1993). 16 E. Veklerov and J. Llacer, IEEE Trans. Med. Imaging 6(4), 313–319 (1987). 17 C. Falcon, I. Juvells, J. Pavia, and D. Ros, Phys. Med. Biol. 43(5), 1271– 1283 (1998). 18 A. Gaitanis, G. Kontaxakis, G. Spyrou, G. Panayiotakis, and G. Tzanakos, Comput. Med. Imaging Graph. 34(2), 131–141 (2010). 19 P. A. Bradley, D. C. Wilson, F. J. Swenson, and G. L. Morgan, Rev. Sci. Instrum. 74(3), 1824–1827 (2003). 20 S. M. Glenn, L. R. Benedetti, D. K. Bradley, B. A. Hammel, N. Izumi, S. F. Khan, G. A. Kyrala, T. Ma, J. L. Milovich, A. E. Pak, V. A. Smalyuk, R. Tommasini, and R. P. Town, Rev. Sci. Instrum. 83(10), 10E519 (2012).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
REVIEW OF SCIENTIFIC INSTRUMENTS 85, 023508 (2014)
Neutron source reconstruction from pinhole imaging at National Ignition Facility P. Volegov,1 C. R. Danly,1 D. N. Fittinghoff,2 G. P. Grim,1 N. Guler,1 N. Izumi,2 T. Ma,2 F. E. Merrill,1 A. L. Warrick,2 C. H. Wilde,1 and D. C. Wilson1 1 2
Los Alamos National Laboratory, Los Alamos, New Mexico 87544, USA Livermore National Laboratory, Livermore, California 94550, USA
(Received 9 October 2013; accepted 30 January 2014; published online 18 February 2014) The neutron imaging system at the National Ignition Facility (NIF) is an important diagnostic tool for measuring the two-dimensional size and shape of the neutrons produced in the burning deuterium-tritium plasma during the ignition stage of inertial confinement fusion (ICF) implosions at NIF. Since the neutron source is small (∼100 μm) and neutrons are deeply penetrating (>3 cm) in all materials, the apertures used to achieve the desired 10-μm resolution are 20-cm long, single-sided tapers in gold. These apertures, which have triangular cross sections, produce distortions in the image, and the extended nature of the pinhole results in a non-stationary or spatially varying point spread function across the pinhole field of view. In this work, we have used iterative Maximum Likelihood techniques to remove the non-stationary distortions introduced by the aperture to reconstruct the underlying neutron source distributions. We present the detailed algorithms used for these reconstructions, the stopping criteria used and reconstructed sources from data collected at NIF with a discussion of the neutron imaging performance in light of other diagnostics. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4865456] I. INTRODUCTION
The National Ignition Facility (NIF) team is currently attempting to achieve thermonuclear ignition of a deuteriumtritium (DT) plasma.1 The strategy being used is to generate the required densities and pressures through indirect ablative compression of plastic (CH) capsules containing DT ice layers, which surround DT gas. In this design, implosion process forms a region of less dense fuel surrounded by higher density DT fuel. The inner region reaches temperatures of ∼10 keV. All of the DT reaches an areal density, ρR, of greater than ∼1 g/cm2 . Fusion reactions and resulting energy deposition from alpha particles further heat this hot spot which then initiates an ignition front that propagates through the surrounding “cold” high density DT fuel. This reaction could release over 106 joules of fusion energy and more than 1017 neutrons. Tuning the implosion to achieve ignition requires the measurement of plasma conditions at very small length and time scales. The hot spot is predicted to have a 30–50 μm diameter, while the surrounding cold fuel is expected to be ∼100 μm in diameter, and the fusion conditions must be maintained for ∼100 ps.1 Establishing the appropriate conditions for ignition, such as symmetry, hot-spot volume, and cold-fuel volume, depends on iterative tuning, and until ignition is achieved, measurements must be made at non-ignition conditions. A neutron imaging diagnostic, known as the Neutron Imaging System (NIS), has been built at NIF and is being used to provide data on the size and shape of the fusion hotspot and the surrounding cold fuel during the initial tuning campaigns.2 This aperture imaging system is composed of two basic pieces: an aperture array that is used to form the neutron images and a scintillator-based detector system that 0034-6748/2014/85(2)/023508/12/$30.00
is used to measure the neutron flux passing through the aperture array. The BCF-99-55 scintillating fiber bundle is viewed with two MCP gated image collection systems, which allows the collection of two independently timed images, providing images from two neutron energy intervals. The detector system is capable of collecting two fastgated and independently timed images. Because the detector array is positioned 28 m from the neutron source, the neutron arrival time at the detector is correlated to the neutron energies, allowing measurement of the neutron source distributions from two energy ranges by gating the detectors at two different times. Typically, one detector is gated from 13 to 17 MeV to view the 14-MeV neutrons that are generated from DT fusion processes and provides information on the size and shape of the hot spot. The second detector is gated to measure the source distribution of lower energy neutrons, typically in the range from 6 to 12 MeV. These lower energy neutrons are predominantly DT-fusion neutrons that have scattered in the surrounding cold fuel and, therefore, provide information on the distribution of the cold fuel surrounding the hotspot. Because the scintillation light from BCF-99-55 decays quickly (∼3 ns) the light from the primary neutrons is essentially gone by the time the scattered neutrons arrive at the scintillator, allowing the two images to be measured with the same scintillator. Since the images are produced in the same scintillator the relative positions of the hot spot and cold fuel images have been determined by collecting the two images simultaneously with a calibration experiment. This relative image registration is accurate to within 5 μm in the source plane. We note that the absolute registration relative to the original source position is not known accurately due to alignment precision of ∼100 μm in the source plane. However, we are
85, 023508-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-2
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
able to extract this information from the data as described in Sec. III. Because of the long mean free path of neutrons within aperture materials, the neutron imaging system forms images using extended apertures or “pinholes” in materials with short interaction lengths. The aperture array currently in use at NIF is composed of triangular pinholes and mini-penumbral apertures machined along the 20 cm length of wedged gold layers. Gold has a mean free path for elastic scattering of ∼3 cm for 14 MeV neutrons. These equilateral triangular pinholes are 226-μm high at the downstream end of the aperture array and taper to 5-μm high at the surface facing the neutron source. The triangular shape of the pinholes results in image distortion while the extended nature of the pinhole results in a nonstationary point spread function generating distortions which vary across the pinhole field of view.3–8 To provide the size and shape information required by the tuning campaigns at the National Ignition Facility, two Maximum-Likelihood reconstruction algorithms have been developed to remove these distortions. In this work, we present details of these algorithms, the stopping criteria used and present reconstructed neutron sources from NIF measurements along with a discussion of the performance in light of other diagnostics. In Sec. II, we present an overview of the NIS system and give a mathematical description of the formation of images with non-stationary point-spread functions. In Sec. III, we describe the employed image processing techniques: the data reduction prior to reconstruction, determination of the aperture array pointing, and the iterative algorithm that we have used to reconstruct the source images including discussion of the stopping criterion. The results of the reconstructions and conclusions will then be discussed in Sec. IV. II. NEUTRON IMAGING SYSTEM OVERVIEW
As shown in Figure 1, the NIS aperture array has 20 triangular pinhole apertures and three mini-penumbral apertures.2 The array was designed so that the fields of view of all apertures overlap at 26.5 cm from the front face. In its usage for the inertial confinement fusion (ICF) implosions at the NIF, the front surface of the aperture array is positioned at 32.5 cm from the source. This results in each aperture pointing to a different location at the source plane, resulting in a combined 0.5 mm field of view. A. Image formation
For a neutron aperture imager, neglecting the neutrons that scatter in the aperture material or secondary scattering of neutrons in the detector, the mean neutron flux at the detection plane can be written as ˆ (nˆ · ) r ˆ 3r , S(r , E, )d F (r, E) = e− ∫r μt (ξ ,E)|dξ | 4π |r − r |2 S
(1) where F(r, E) is the neutron flux at the detection plane, i.e., F(r, E)AE is a number of neutrons crossed differential area A of detection plane about r, with differential energies E energy about E; μt (ξ , E) is the total linear attenuation co-
FIG. 1. Photograph of the rear-surface of the NIS aperture array. The array has 20 pinhole apertures with triangular cross sections arranged in four rows, with three penumbral apertures in the center row. The triangular apertures shown here are 226-μm high and taper to 5-μm high at the front surface of the array, which faces the source.
efficient for neutrons with energy E at the point ξ ; nˆ is a unit ˆ ≡ (r, ˆ r ) ≡ r−r vector normal to the detection surface; |r−r | is a unit vector in the direction from the point r to the ˆ is the differential neutron source, i.e., point r; S r , E, ˆ V E is a number of neutrons emitted in S r, E, the volume V about the point r , with differential energies E energy about E, in the solid angle about direction deˆ and integral over S is an integral fined by the unit vector ; over the source volume. In Eq. (1) and after we use the following conventions: bold font denotes vectors and matrices, italic denotes scalars, a hat sign denotes a unit vector, a capital italic superscript T denotes a transposed matrix/vector, and the vectors are column vectors. In general, a linear, position-sensitive detector can be described by a spatial sensitivity function, also called a point spread function of a detector, P(r, r , E), such that P(r, r , E)AE is the signal generated in the sensor at the position r by a unit neutron flux through the differential area A of the detection plane about r , with differential energies E about a central energy, E. In this case, we can write the recorded image as ∞ I (r) = 0
ˆ 3 r dE, K(r, r , E)S(r , E, )d
(2)
S
where the total sensitivity function of the neutron imaging system is defined by the following formula: ˆ (nˆ · ) r K(r, r , E) = P (r, r , E)e− ∫r μt (ξ ,E)|dξ | d 3 r . 4π |r −r |2 S
(3)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
Volegov et al.
023508-3
Rev. Sci. Instrum. 85, 023508 (2014)
Here the integration is over the detector sensitive volume, for example, a scintillator surface. It should be noted that the function K(r, r , E) models the blurring of the imaging system and is called the point spread function (PSF) of the image collection system. The last two equations can be simplified assuming that the extent of the neutron source and the detector are small compared to the distance between them. In this case we can reduce the dimensionality of the problem: ∞ I (u) =
K(u, v, E)S(v, E) d 2 vdE. 0
(4)
S
Here the vectors u = r − rsens and v = r − rsrc define a position on the detector plane and the source plane, respectively, and S (v, E) is a projection of the neutron source on the plane orthogonal to the axis connecting the center of the source, rsrc , and the center of the detector, rsens , which is the central axis of the imaging system: ∞ S(v, E) =
ˆ 0 l, E)dl, S(r +
(5)
−∞ rsens −rsrc . |rsens −rsrc |
ˆ0 ≡ Assuming that both the where the unit vector detector plane and the source plane are orthogonal to the central axis of the imaging system, the sensitivity function can be written as 1 u P (u, u , E)e− ∫v μt (ξ ,E)|dξ | d 2 u . K(u, v, E) = 2 4π L Ds
(6) Here Ds denotes the integration over the detector plane, L = |rsens − rsrc | is the distance between the source center and the detector plane. This approximation is called the paraxial approximation or a small-angle approximation. Assuming that the energy distribution of the source is independent of the position so that the function S (v, E) is separable and S (v, E) = S (v) U (E), we can further simplify the imaging equation: (7) I (u) = K(u, v)S(v)d 2 v, S
where energy distribution function U(E) is considered to be known. In practice, using time-of-flight techniques, the images are recorded in several energy intervals, assuming that spatial and energy distributions of the neutrons are separable inside each interval. The kernel function K(u, v) in this case is defined by K(u, v) ∞ 1 u U (E) P (u, u , E)e− ∫v μt (ξ ,E)|dξ | d 2 u dE. = 2 4π L 0
Ds
(8)
FIG. 2. A combination of images, i.e., point-spread functions, of point sources displaced in x and y in steps of 100-μm from the central axis of a NIF triangular pinhole. As it can be seen, the shape of a point-spread function strongly depends on position of the source relative to the central axis of an aperture.
Typically a point spread function of a detector P(u, u , E) is modeled by a Gaussian: P (u, u , E) = √
1 2π σ
e−
|u−u |2 2σ 2
,
(9)
where σ is the RMS of the Gaussian blur. The spatial dependence of the point-spread function for a triangular pinhole aperture used at NIF is illustrated in Figure 2, which shows an image corresponding to a combination of point sources, each located on a 100-μm grid. The equations that have been considered so far connecting the neutron source with the recorded image are idealized imaging equations because they do not take into account various imperfections of the imaging system, such as the contribution of the scattered neutrons or imperfections in the detection system. III. IMAGE PROCESSING A. Data reduction
Data collected from the two imaging systems are normalized before source reconstructions are performed. In the process of collecting data for a NIF experiment, three dark field images are collected to measure the baseline levels of the CCD camera system. These dark fields are averaged together and subtracted from the data images. During image collection some neutrons are scattered from the flight path and interact directly with elements of the CCD array. Each direct interaction generates a few bright pixels, often called “stars.” After the dark fields are subtracted these stars are differentiated from statistical deviations of the local signal levels and set to the average of the neighboring pixels. At the present yields (∼1014 neutrons) less than 0.1% of the pixels are affected by this procedure.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-4
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
FIG. 3. The image produced by the aperture array is shown with the same contrast range before (left) and after (right) the scatter background subtraction. The images are 521 by 521 pixels with each pixel being 341 μm.
Once the average dark field is subtracted and stars are removed, the images are divided by flat-field images.9 The recording system has been calibrated and is operated in a linear regime both during pinhole imaging and flat-field collection. The flat-field images are measurements of the light generated by a uniform distribution of neutrons interacting with the scintillating fiber array. By dividing the data images by the flat fields, the variation in response across the scintillator and image collection systems is removed. The image is then typically re-binned from the substantial pixel oversampling of 37 μm pixels at the scintillator plane to 341 μm pixels, which maps to ∼4 μm pixels at the source location. The remainder of the analysis has been typically performed at this pixel size. Of the 40% of neutrons that interact in the scintillator, some small fraction, ∼10%, interact a second time before exiting the scintillator material. This results in a long-range scatter background, which extends to more than 3 cm. Because of the extent of this long-range scatter background, it is of little consequence to the resolution for a single pinhole, where peak scatter background is 0, pi, j > 0, xj(0) > 0 and ensures that the solution is also positive, which is consistent with Poisson statistics.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-8
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
For the Gaussian noise model, the likelihood function is defined by L(S) = (2π )− 2 |Cεε |− 2 e− 2 (I−KS) m
1
1
T
C−1 εε (I−KS)
,
(21)
where Cεε is a m × m covariance matrix of a random vector ε, |Cεε | denotes the determinant of the covariance matrix, C−1 εε is inverse of the covariance matrix, and m is the total number of pixels in the recorded image I. The log-likelihood function in this case is 1 l(S) = − (I − KS)T C−1 εε (I − KS) + const. 2
(22)
Using this likelihood function the MLE solution of the problem (13) is simply −1 T −1 K Cεε I, (23) SMLE = KT C−1 εε K which is exactly a weighted least squares solution. If the noise is uncorrelated, i.e., the image pixels are independent the matrix Cεε is a diagonal matrix with the diagonal elements σi2 , i = 1, .., m being a variance of the noise in the ith pixel. Direct application of Eq. (23) in the context of reconstruction of neutron images is somewhat impractical due to difficulties in inverting a large and ill-conditioned matrix KT C−1 εε K. An iteration scheme, which we will refer to as MLGAUSS, similar to algorithm for the Poisson noise described above, see Eq. (19), was developed for Gaussian noise distributions:15 m yi − y˜i(t) (t+1) (t) xj = xj 1 + h pi,j , σi2 i=1 (24) j = 1, .., n. The step size h is computed from 1D maximization of the likelihood function, Eq. (22). The details of the step size computation are described in Ref. 15. It should be noted that the formula (24) does not require that quantities yi , pi, j , and xj(t) be positive, and does not ensure that the solution is positive. In the case of Poisson noise, where a variance of the noise in a pixel is equal to the mean value of the signal in this pixel, i.e., σi2 = y˜i , i = 1, ..., m, and assuming the step size h ≡ 1, formula (24) corresponds to the iteration scheme shown in Eq. (19). An important question for practical application of the iteration scheme (24), or (19), is when to stop the iterations. It is known that when the described algorithm passes a certain point the image starts to deteriorate. This optimal point of truncation depends on statistics (noise): the higher the statistics, the later in the iteration the degradation begins. The number of iterations, therefore, becomes the regularization parameter. The goal then is to stop the iterations as soon as the solution is consistent with the data within the experimental accuracy of the measurement. The agreement of the solution with the data is usually characterized by normalized chi-square (or reduced chi-squared) criteria, defined by the following formula: 1 yi − y˜i . m i=1 σi2 m
χn2 =
(25)
For a solution that is consistent with the data, the expectation value of χn2 is 1, and the variance is 2/m. Thus the observed value of χn2 is easy to interpret: if χn2 is about 1 within, say, a variance, then our solution is consistent with the data, if χn2 1, then the solution captures not only the neutron distribution but also the noise (or the noise is overestimated), and finally if χn2 1, then the solution is not consistent with the data (or the noise is underestimated). These simple rules may be formalized by introducing P- and Q-values. Given a particular value, we can calculate the probability to obtain χn2 that is less than the P-value and the probability to obtain χn2 that is greater than the Q-Value. The most widely accepted critical value is Q = 0.05, meaning that only in 5% of the trials would result in a higher value of χn2 . This leads to the simple stopping rule: continue iterations until the Q-Value is higher than the threshold value. To apply the stopping rule described in the previous paragraph, the variance of the noise should be reliably estimated. In practice, this is not always possible. An alternative is to compute the relative change of the residual from one iteration to the next iteration: ρ ≡ (ρ (t) − ρ (t + 1) /ρ (t) ), where ρ (t) is some measure of the difference, and stop the iterations as soon as this relative difference is below some value. Typically, reduced chi-square is used as a measure of the difference: 1 yi − y˜i(t) . m i=1 σi2 m
ρ (t) =
(26)
However, one may employ a different measure such as the Kullback-Leibler divergence: m yi (t) ρ = yi ln . (27) y˜i(t) i=1 Empirically it was found that critical value 10−3 provides a robust stopping criterion in most of the cases. More complicated and sophisticated stopping criteria exist; see, for example, Refs. 15–18. These criteria will be investigated in future work. IV. APPLICATION OF MLE ALGORITHMS TO DATA
Many neutron images have been collected during the ignition tuning experiments at NIF over the last three years. The neutron source distributions have each been reconstructed with the algorithms described above. Typically these reconstructed sources are 50–100 μm elliptical sources with very little features. This type of reconstruction has proved little opportunity to assess the performance of these reconstruction algorithms. Recently, however, a new set of experiments have been performed to study hydro-dynamic implosion physics at NIF. These experiments are designed to be driven with a stronger shock and thus implode along a higher adiabat of the thermonuclear fuel, resulting in a larger stagnation radius and larger neutron images. In addition, the initial experiments, which were executed to tune the shape of the implosion, resulted in a source with a shape more complicated than the typical simple ellipsoid, and with record neutron yields. This more complicated shape with high
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-9
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
FIG. 9. Primary neutron image (a) and down-scattered neutrons image (b) for a high-yield NIF shot (NIF shot number N130710-002-999).
statistics has provided a feature rich source for tests of these reconstruction algorithms. The raw images from this experiment (NIF shot number N130710-002-999) for the primary neutrons and for the down-scattered neutrons are shown in Figure 9. The reconstructed primary source distributions using two different algorithms (EMML and MLGAUSS) are shown in Figure 10. Convergence for both of these algorithms is illustrated by
the χ 2 evolutions with respect to the number of iterations shown in Figure 11. In spite of apparent difference in convergence behavior, after 50 iterations both algorithms demonstrated almost exactly the same convergence: EMML – χ 2 = 1.31, χ 2 /χ 2 = 0.8K × 10−3 , and MLGAUSS – χ 2 = 1.32, χ 2 /χ 2 = 1.2 × 10−3 . The reconstructed source distributions, as it can be seen from top raw images in Figure 10, are similar but differ in size: MLGAUSS produces bigger and
FIG. 10. Primary neutron image reconstruction of the source by using EMML (left: (a) 50 iterations and (c) 262 iterations) and MLGAUSS (right: (b) 50 iterations and (d) 262 iteration) algorithms (NIF shot number N130710-002-999).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-10
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
FIG. 13. Overlaid images filtered in red and cyan colors for primary and downscattered sources, respectively (NIF shot number N130710-002-999).
FIG. 11. Reduced χ 2 as a function of iteration (NIF shot number N130710002-999, primary source reconstruction).
smother source. In this example EMML algorithm changes χ 2 very little after 50 iterations, specifically after 262 iterations: χ 2 = 1.28, χ 2 /χ 2 = 3.5 × 10−5 , while perceptibly changing the reconstructed source (see Figures 10(a) and 10(c)) making it implausible from the physical considerations. MLGAUSS, however, exhibits an opposite behavior: the residual changes perceptibly – after 262 iterations: χ 2 = 1.19, χ 2 /χ 2 = 7.0 × 10−7 – while change of the reconstructed source is not that drastic. Explanation of such a behavior lays in the background noise. As it was mentioned in Sec. III A, dark field images are subtracted from the data images with the result that the noise in the background pixels is no longer consistent with a Poisson noise model: the pixels can take negative values, and generally the average of the background pixels is not necessarily a zero, it could be negative. Such an inconsistency, even with non-positive pixels being discarded, drives fast corruption of the reconstructed image using EMML algorithm which assumes a Poisson noise model. The MLGAUSS algorithm on the other hand does not require non-negativity of the data and thus exhibits more robust behavior with respect to the
background subtraction procedure. Summarizing, in selecting between those two algorithms – EMML and MLGAUSS – one must take into account that EMML requires that the data be all positive and ensures that the reconstructed image is also positive, which is a very strong regularization provision. On the other hand, MLGAUSS does not impose nonnegativity on the reconstructed image, but does not require non-negativity of the data and handles negative values from background subtraction, making this the preferred algorithm for this application. Application of the described algorithm to the low statistic data is illustrated using down-scattered neutrons image (see Figure 9(b)). Here intensity of the down-scattered image is about 20 times lower than the intensity of the primary image. However, even in this case, both algorithms reconstruct a plausible neutron source as illustrated in Figure 12, which shows the image and the reconstructed source for the down-scattered neutrons, again with the two algorithms. Note that reconstruction of low statistic data requires significantly less iterations: EMML – 12 iterations and MLGAUSS – 28 iterations. It is important to compare and validate the reconstructed source distributions against physical models of the process and relevant experimental results obtained using a different technique. Figure 13 shows in red the reconstructed
FIG. 12. Down-scattered neutron image reconstruction of the source by using EMML ((a) 12 iterations) and MLGAUSS ((b) 28 iterations) algorithms (NIF shot number N130710-002-999).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-11
Volegov et al.
Rev. Sci. Instrum. 85, 023508 (2014)
demonstrating that these algorithms are effective in removing the unwanted image distortions. The algorithm based on Gaussian noise model – MLGAUSS – has been found to be tolerant to background subtraction and, therefore, used in these applications. Algorithms based on a Poisson noise model, such as EMML, maybe better for other applications. These algorithms are applicable to a wide range of imaging modalities; from pinhole imaging to penumbral imaging and we are investigating how this new analysis capability can be used in the future to enable new imaging systems with higher resolution at low statistics, such as neutron imaging with coded apertures.
ACKNOWLEDGMENTS
FIG. 14. The X-ray image polar view, the NIS LOS is shown with the arrow (NIF shot number N130710-002-999).
distribution of the primary fusion neutrons generated in the hot DT plasma overlaid on the distribution of the cold fuel in cyan, which has been reconstructed from the down scattered neutron image. The reconstructed primary image clearly shows a source distribution with two connected lobes, consistent with a toroid of burning nuclear fuel. The cold fuel distribution is consistent with this toroid of burning fuel with cold fuel penetrating through the vertical axis of the burning plasma. An x-ray image collected along the vertical axis shows the orthogonal view from the top of the toroid shown in Figure 14. One would expect this high energy x-ray image to be consistent with a toroidal primary neutron image, as only the hottest regions of the nuclear fuel will generate fusion neutrons and this region will also be the brightest emitter of x-rays.19, 20 A comparison of the raw image to the reconstructed source shows that a significant level of distortions has been removed to form the source neutron distribution. The fact that the resulting source reconstruction is consistent with the orthogonal x-ray image demonstrates that this method of reconstruction is effective for this application. V. CONCLUSIONS
Iterative Maximum Likelihood Estimation algorithms were used to develop image reconstruction tools for the neutron images that have been recently collected at NIF. Accurate pointing of the aperture axis to the neutron source has proven to be an important parameter in accurately reconstructing neutron source distributions. The array of multiple pinhole apertures has allowed us to extract accurate pointing information, which is necessary in the calculation of the proper non-stationary point spread function for reconstructions. These tools have been used to analyze dozens of neutron images. Most recently these algorithms have been used to reconstruct a toroidal neutron source, which required the removal of significant image distortions due to the pinhole aperture and detector system. This reconstructed source shape has been independently validated with an orthogonal x-ray view,
Additional credit goes to the dedicated staff and technicians of NIF, whose hard work and operational expertise resulted in the data that are shown here. The authors wish to acknowledge D. Jedlovec, M. A. Talison, O. Drury, D. Kalantar, and R. Wood for their hard work and expertise in fielding and aligning the neutron imaging system. This work has been performed under the auspices of the U.S. Department of Energy for NNSA Campaign 10 (Inertial Confinement Fusion) with Steve Batha as program manager. 1 S. H. Glenzer, D. A. Callahan, A. J. MacKinnon, J. L. Kline, G. Grim, E. T.
Alger, R. L. Berger, L. A. Bernstein, R. Betti, D. L. Bleuel, T. R. Boehly, D. K. Bradley, S. C. Burkhart, R. Burr, J. A. Caggiano, C. Castro, D. T. Casey, C. Choate, D. S. Clark, P. Celliers, C. J. Cerjan, G. W. Collins, E. L. Dewald, P. DiNicola, J. M. DiNicola, L. Divol, S. Dixit, T. Döppner, R. Dylla-Spears, E. Dzenitis, M. Eckart, G. Erbert, D. Farley, J. Fair, D. Fittinghoff, M. Frank, L. J. A. Frenje, S. Friedrich, M. Gatu Johnson, C. Gibson, E. Giraldez, V. Glebov, S. Glenn, N. Guler, S. W. Haan, B. J. Haid, B. A. Hammel, A. V. Hamza, C. A. Haynam, G. M. Heestand, M. Hermann, H. W. Hermann, D. G. Hicks, D. E. Hinkel, J. P. Holder, D. M. Holunda, J. B. Horner, W. W. Hsing, H. Huang, N. Izumi, M. Jackson, O. S. Jones, D. H. Kalantar, R. Kauffman, J. D. Kilkenny, R. K. Kirkwood, J. Klingmann, T. Kohut, J. P. Knauer, J. A. Koch, B. Kozioziemki, G. A. Kyrala, A. L. Kritcher, J. Kroll, K. La Fortune, L. Lagin, O. L. Landen, D. W. Larson, D. LaTray, R. J. Leeper, S. Le Pape, J. D. Lindl, R. Lowe-Webb, T. Ma, J. McNaney, A. G. MacPhee, T. N. Malsbury, E. Mapoles, C. D. Marshall, N. B. Meezan, F. Merrill, P. Michel, J. D. Moody, A. S. Moore, M. Moran, K. A. Moreno, D. H. Munro, B. R. Nathan, A. Nikroo, R. E. Olson, C. D. Orth, A. E. Pak, P. K. Patel, T. Parham, R. Petrasso, J. E. Ralph, H. Rinderknecht, S. P. Regan, H. F. Robey, J. S. Ross, M. D. Rosen, R. Sacks, J. D. Salmonson, R. Saunders, J. Sater, C. Sangster, M. B. Schneider, F. H. Séguin, M. J. Shaw, B. K. Spears, P. T. Springer, W. Stoeffl, L. J. Suter, C. A. Thomas, R. Tommasini, R. P. J. Town, C. Walters, S. Weaver, S. V. Weber, P. J. Wegner, P. K. Whitman, K. Widmann, C. C. Widmayer, C. H. Wilde, D. C. Wilson, B. Van Wonterghem, B. J. MacGowan, L. J. Atherton, M. J. Edwards, and E. I. Moses, Phys. Plasmas 19(5), 056318 (2012). 2 F. E. Merrill, D. Bower, R. Buckles, D. D. Clark, C. R. Danly, O. B. Drury, J. M. Dzenitis, V. E. Fatherley, D. N. Fittinghoff, R. Gallegos, G. P. Grim, N. Guler, E. N. Loomis, S. Lutz, R. M. Malone, D. D. Martinson, D. Mares, D. J. Morley, G. L. Morgan, J. A. Oertel, I. L. Tregillis, P. L. Volegov, P. B. Weiss, C. H. Wilde, and D. C. Wilson, Rev. Sci. Instrum. 83(10), 10D317 (2012). 3 D. C. Wilson, G. P. Grim, I. L. Tregillis, M. D. Wilke, M. V. Patel, S. M. Sepke, G. L. Morgan, R. Hatarik, E. N. Loomis, C. H. Wilde, J. A. Oertel, V. E. Fatherley, D. D. Clark, D. N. Fittinghoff, D. E. Bower, M. J. Schmitt, M. M. Marinak, D. H. Munro, F. E. Merrill, M. J. Moran, T. S. F. Wang, C. R. Danly, R. A. Hilko, S. H. Batha, M. Frank, and R. Buckles, Rev. Sci. Instrum. 81(10), 10D335 (2010). 4 E. Loomis, G. Grim, C. Wilde, D. C. Wilson, M. Wilke, J. Finch, G. Morgan, I. Tregillis, and D. Clark, Sixth International Conference on Inertial Fusion Sciences and Applications, Parts 1–4, Vol. 244, 2010.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
023508-12 5 G.
Volegov et al.
P. Grim, R. D. Day, D. D. Clark, V. E. Fatherley, F. P. Garcia, S. A. Jaramillo, A. J. Montoya, G. L. Morgan, J. A. Oertel, T. A. Ortiz, J. R. Payton, P. D. Pazuchanics, D. A. W. Schmidt, A. C. Valdez, C. H. Wilde, and M. D. Wilke, Proc. SPIE 6707, H7070 (2007). 6 C. A. Barrera, E. C. Morse, and M. J. Moran, Rev. Sci. Instrum. 77(10), 10E716 (2006). 7 M. D. Wilke, S. H. Batha, P. A. Bradley, R. D. Day, D. D. Clark, V. E. Fatherley, J. P. Finch, R. A. Gallegos, F. P. Garcia, G. P. Grim, S. A. Jaramillo, A. J. Montoya, M. J. Moran, G. L. Morgan, J. A. Oertel, T. A. Ortiz, J. R. Payton, P. Pazuchanics, D. W. Schmidt, A. C. Valdez, C. H. Wilde, and D. C. Wilson, Rev. Sci. Instrum. 79(10), 10E529 (2008). 8 V. Y. Glebov, D. D. Meyerhofer, T. C. Sangster, C. Stoeckl, S. Roberts, C. A. Barrera, J. R. Celeste, C. J. Cerjan, L. S. Dauffy, D. C. Eder, R. L. Griffith, S. W. Haan, B. A. Hammel, S. P. Hatchett, N. Izumi, J. R. Kimbrough, J. A. Koch, O. L. Landen, R. A. Lerche, B. J. MacGowan, M. J. Moran, E. W. Ng, T. W. Phillips, P. M. Song, R. Tommasini, B. K. Young, S. E. Caldwell, G. P. Grim, S. C. Evans, J. M. Mack, T. J. Sedillo, M. D. Wilke, D. C. Wilson, C. S. Young, D. Casey, J. A. Frenje, C. K. Li, R. D. Petrasso, F. H. Seguin, J. L. Bourgade, L. Disdier, M. Houry, I. Lantuejoul, O. Landoas, G. A. Chandler, G. W. Cooper, R. J. Leeper, R. E. Olson, C. L. Ruiz, M. A. Sweeney, S. P. Padalino, C. Horsfield, and B. A. Davis, Rev. Sci. Instrum. 77(10), 10E715 (2006).
Rev. Sci. Instrum. 85, 023508 (2014) 9 N.
Guler, P. Volegov, C. R. Danly, G. P. Grim, F. E. Merrill, and C. H. Wilde, Rev. Sci. Instrum. 83(10), 10D316 (2012). 10 M. Z. Tarasko, Preprint of the Physics and Energy Institute, PEI-156 (1969). 11 A. J. Rockmore and A. Macovski, IEEE Trans. Nucl. Sci. 23(4), 1428–1432 (1976). 12 A. P. Dempster, N. M. Laird, and D. B. Rubin, J. R. Stat. Soc. Ser. B 39(1), 1–38 (1977). 13 W. H. Richardson, J. Opt. Soc. Am. 62(1), 55–59 (1972). 14 L. B. Lucy, Astron. J. 79(6), 745–754 (1974). 15 V. I. Gelfgat, E. L. Kosarev, and E. R. Podolyak, Comput. Phys. Commun. 74(3), 335–348 (1993). 16 E. Veklerov and J. Llacer, IEEE Trans. Med. Imaging 6(4), 313–319 (1987). 17 C. Falcon, I. Juvells, J. Pavia, and D. Ros, Phys. Med. Biol. 43(5), 1271– 1283 (1998). 18 A. Gaitanis, G. Kontaxakis, G. Spyrou, G. Panayiotakis, and G. Tzanakos, Comput. Med. Imaging Graph. 34(2), 131–141 (2010). 19 P. A. Bradley, D. C. Wilson, F. J. Swenson, and G. L. Morgan, Rev. Sci. Instrum. 74(3), 1824–1827 (2003). 20 S. M. Glenn, L. R. Benedetti, D. K. Bradley, B. A. Hammel, N. Izumi, S. F. Khan, G. A. Kyrala, T. Ma, J. L. Milovich, A. E. Pak, V. A. Smalyuk, R. Tommasini, and R. P. Town, Rev. Sci. Instrum. 83(10), 10E519 (2012).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitationnew.aip.org/termsconditions. Downloaded to IP: 137.189.170.231 On: Tue, 23 Dec 2014 15:07:36
