Multiaxial neutron/x-ray imaging and three-dimensional (3D) reconstruction techniques play a crucial role in gaining valuable insights into the generation and evolution mechanisms of pulsed radiation sources. Owing to the short emission time (∼200 ns) and drastic changes of the pulsed radiation source, it is necessary to acquire projection data within a few nanoseconds in order to achieve clear computed tomography 3D imaging. As a consequence, projection data that can be used for computed tomography image reconstruction at a certain moment are often available for only a few angles. Traditional algorithms employed in the process of reconstructing 3D images with extremely incomplete data may introduce significant distortions and artifacts into the final image. In this paper, we propose an iterative image reconstruction method using cylindrical harmonic decomposition and a self-supervised denoising network algorithm based on the deep image prior method. We augment the prior information with a 2D total variation prior and a 3D deep image prior. Single-wire Z-pinch imaging experiments have been carried out at Qin-1 facility in five views and four frames, with a time resolution of 3 ns for each frame and a time interval of 40 ns between adjacent frames. Both numerical simulations and experiments verify that our proposed algorithm can achieve high-quality reconstruction results and obtain the 3D intensity distribution and evolution of extreme ultraviolet and soft x-ray emission from plasma.
I. INTRODUCTION
Significant advances have been made in inertial confinement fusion (ICF) ignition techniques such as laser ignition^{1,2} and Z pinches^{3–5} in recent years. Spatial structure diagnosis of the burning plasma during the ignition stage of ICF implosions is important for evaluating the ignition state, verifying physical models, and improving experimental parameters. Multiaxial passive imaging is generally adopted to acquire three-dimensional (3D) information of pulsed radiation sources. Only several two-dimensional (2D) line-integral projections obtained at a few views are available, owing to the short durations of pulsed sources and the limitations of current technology. Projection images from up to four views are captured at the Omega Laser^{7} (x-ray), the National Ignition Facility (NIF)^{8} (deuterium and tritium neutrons), and the Qin-1 facility^{6} [extreme ultraviolet (XUV)/soft x-ray (SR)]. The reconstruction of a 3D source with a minimal number of 2D projections is a severely ill-posed problem. Some methods developed in computed tomography, such as the algebraic reconstruction technique and the expectation maximization algorithm, have been applied to the reconstruction of pulsed radiation sources,^{7,9,10} but these introduce severe artifacts into the reconstructed 3D image.
An approach based on basis function expansion has been proposed that increases the possibility of accurate reconstruction of pulsed radiation source with few-view data. Chen et al.^{11} represented the radiation source as a sum of multiple 3D Gaussian functions to perform 3D reconstruction with x-ray pinhole projection data from three angles at the Gekko XII facility. Li et al.^{12} leveraged the zeroth-order L shell to reconstruct the plasma density distribution of the plasmasphere in the geomagnetic equatorial plane from extreme ultraviolet sensor data from the IMAGE satellite. Pecora^{13} derived an analytical reconstruction algorithm by decomposing the radiation source into a linear combination of spherical harmonic functions. Volegov et al.^{7} applied this analytical algorithm to reconstruct the x-ray source at the Omega Laser and the deuterium and tritium neutron source at the NIF. They later developed an analytical reconstruction method based on cylindrical harmonic decomposition and used two-axis imaging data to reconstruct the neutron source at the NIF.^{8} Analytical methods based on basis function expansions are susceptible to the impact of noise in the measured projections, and, owing to the limited number of views, the maximum expansion order is limited, resulting in insufficient function expression capabilities, which are prone to introducing artifacts in the reconstructed image.^{14}
Following the realization that convolutional neural networks (CNNs) have the capacity for application to image processing,^{15–18} deep learning methods have been considered for use in the post-processing procedure. Many effective algorithms based on CNNs rely on training sets containing large amounts of data to adjust the parameters. Nevertheless, it is difficult to have a large amount of labeled data in pulsed radiation imaging, owing to the lack of real 3D distribution data. Ulyanov et al.^{19} proposed the concept of a deep image prior (DIP). According to this, the CNN structure itself has universal image prior information. The image can be recovered by updating parameters of the generator network iteratively without a training set, which makes it possible to suppress artifacts in the image.
In this paper, we propose an iterative 3D image reconstruction algorithm based on cylindrical harmonic decomposition. It can increase the expansion order of the basis functions and boost the representation ability of the cylindrical harmonic functions. Furthermore, we propose an image post-processing method based on the DIP, which reduces artifacts in the reconstructed results through unsupervised learning.
II. METHOD
The initial value will affect the result of the iteration under the condition of incomplete data. In the algorithm, we take the result of the analytical method of cylindrical harmonic function decomposition as the initial value of the iterative method. The maximum expansion order of the analytical method is M_{a} = V. The iterative method is shown in Algorithm 1.
Parameters: M, maximum expansion order; K, iteration number; ɛ, difference between the projection images of the results in two adjacent iterations; λ, μ_{0}, d_{0}, b_{0}, parameters of ADMM algorithm; γ, η, adaptive updating parameters. |
1: Calculate s_{m} by analytical algorithm, t ← 0, s^{(0)} ←s_{m}. |
2: Compute $S0r$, $Irec0\u2190Pp\u0302S0$. |
3: for k = 1, …, K do |
4: t ← t + 1 |
5: for n_{z} = 1: N_{z} do |
6: Obtain s^{(t)} by L-BFGS algorithm. |
7: Compute $S(t)z$, set $S(t)zj=0$ if $S(t)zj<0$. |
8: $dt\u2190shrink\phi st+bt\u22121,1\mu t\u22121$ |
9: $bt\u2190bt\u22121+\phi st\u2212dt$ |
10: Update μ_{t} |
11: end for |
12: if $1NxNyIrect\u2212Irec2t\u221212\u2264\epsilon $, break; |
13: end for |
Parameters: M, maximum expansion order; K, iteration number; ɛ, difference between the projection images of the results in two adjacent iterations; λ, μ_{0}, d_{0}, b_{0}, parameters of ADMM algorithm; γ, η, adaptive updating parameters. |
1: Calculate s_{m} by analytical algorithm, t ← 0, s^{(0)} ←s_{m}. |
2: Compute $S0r$, $Irec0\u2190Pp\u0302S0$. |
3: for k = 1, …, K do |
4: t ← t + 1 |
5: for n_{z} = 1: N_{z} do |
6: Obtain s^{(t)} by L-BFGS algorithm. |
7: Compute $S(t)z$, set $S(t)zj=0$ if $S(t)zj<0$. |
8: $dt\u2190shrink\phi st+bt\u22121,1\mu t\u22121$ |
9: $bt\u2190bt\u22121+\phi st\u2212dt$ |
10: Update μ_{t} |
11: end for |
12: if $1NxNyIrect\u2212Irec2t\u221212\u2264\epsilon $, break; |
13: end for |
To enhance the quality of reconstructed images, we reprocess the iterative results of the cylindrical harmonic decomposition. We utilize a CNN with coding structure similar to U-Net^{25} to reduce artifacts from the reconstructed 3D image by unsupervised learning. The structure of the neural network is illustrated in Fig. 1. The channels are gradually doubled from 8 to 128 in the contracting path with downsampling convolutions, while they are halved in the expansive path. There are two concatenations in the middle of the network structure. A sigmoid activation function is appended at the final layer to map the output values to the interval of (0,1). The length and width of the neural network input are the same as the image size of the reconstruction result. The input value is random noise with a Gaussian distribution. The output is an image of the same size as the input.
Parameters: σ, standard deviation of Gaussian noise; K_{e}, epochs; |
l_{r}, learning rate |
1 η ∼ U(0, 0.1): |
2: for t = 1, …, K_{e} do |
3: ξ ∼N(0, σ) η ←η + ξ, |
4: calculate loss L_{t} |
5: update Θ using Adam |
6: end for |
7: $S\u2190f\Theta \eta $ |
Parameters: σ, standard deviation of Gaussian noise; K_{e}, epochs; |
l_{r}, learning rate |
1 η ∼ U(0, 0.1): |
2: for t = 1, …, K_{e} do |
3: ξ ∼N(0, σ) η ←η + ξ, |
4: calculate loss L_{t} |
5: update Θ using Adam |
6: end for |
7: $S\u2190f\Theta \eta $ |
III. SIMULATION
We reconstructed the simulated data of an aluminum single-wire Z pinch. Owing to technological limitations, we can currently only simulate 2D cross-sectional data of the plasma. Axisymmetric rotation on the sectional data is performed to obtain 3D data. Then, the 3D data are projected in five horizontal angles. Figure 2 shows the distribution of these five angles.
We apply the analytical method, iterative method, and the iterative and DIP post-processing algorithm, respectively, to reconstruct the source. In the iterative method, the maximum expansion order M = 5 and the number of iterations K = 20. The other parameters are set as follows: λ = 1, μ_{0} = 0.01, d_{0} = 0, b_{0} = 0, γ = 2, and η = 0.9. In the DIP procedure, the standard deviation σ = 1/30, the epochs K_{e} = 2400, and the learning rate l_{r} = 0.01. The parameters of the loss function are set as α_{1} = 0.9, α_{2} = 0.1, α_{3} = 0.1, and α_{4} = 0.1. The code in this paper consists of MATLAB and Python code, running on MATLAB R2022a and Python3.9 software, respectively. The hardware platform is an Intel Xeon Platinum 9242 CPU@2.30 GHz and a 24 GB NVIDIA RTX3090 GPU. The analytical method and iterative method take ∼7 and 30 min, respectively, to get the solutions. The running time of the iterative and DIP post-processing algorithm is relatively long, with the iterations of the DIP taking most of the time (2.5 h for 2400 iterations) and a complete reconstruction taking about 3 h.
. | Algorithm . | SSIM . | PSNR . | RMSE(S) . | RMSE(I) . |
---|---|---|---|---|---|
Without noise | AM | 0.811 15 | 27.596 | 0.041 71 | 7.7679 |
IM | 0.955 87 | 39.872 | 0.010 15 | 1.0324 | |
IM+DIP | 0.979 49 | 42.454 | 0.007 54 | 0.2957 | |
With noise | AM | 0.750 02 | 26.971 | 0.044 82 | 8.1900 |
IM | 0.828 43 | 32.365 | 0.024 09 | 2.3296 | |
IM+DIP | 0.955 21 | 39.607 | 0.010 46 | 0.5068 |
. | Algorithm . | SSIM . | PSNR . | RMSE(S) . | RMSE(I) . |
---|---|---|---|---|---|
Without noise | AM | 0.811 15 | 27.596 | 0.041 71 | 7.7679 |
IM | 0.955 87 | 39.872 | 0.010 15 | 1.0324 | |
IM+DIP | 0.979 49 | 42.454 | 0.007 54 | 0.2957 | |
With noise | AM | 0.750 02 | 26.971 | 0.044 82 | 8.1900 |
IM | 0.828 43 | 32.365 | 0.024 09 | 2.3296 | |
IM+DIP | 0.955 21 | 39.607 | 0.010 46 | 0.5068 |
IV. EXPERIMENTAL RESULTS
We carried out a five-axis imaging experiment of a single-wire Z pinch on the double pulse current generator Qin-1 facility, which couples a ∼10 kA 20 ns prepulse generator and a ∼450 kA 400 ns main current generator.^{27,28} Details of the Qin-1 facility can be found in Ref. 29. We measured the extreme ultraviolet/soft x-ray (XUV/SR) photons emitted from the Z-pinch plasma at five angles. Figure 7 illustrates the components of the XUV/SR pinhole imaging system. The values of the optical image captured by the CMOS camera are approximately proportional to the number of photons emitted by the plasma during the gating time (∼3 ns) of the microchannel plate. The principle of the pinhole imaging system is described in Ref. 6. The system can obtain 2D integral projections of the x-ray source in five views at four moments in a one-shot experiment.
A 30 µm silver single wire with a height of 20 mm was used for Z-pinch imaging in shots 2 023 052 405 and 2 023 052 406. The object distance was 300 mm, and the image distance was 224 mm. The four different imaging moments were T_{1} = 0 ns, T_{2} = 40 ns, T_{3} = 80 ns, and T_{4} = 120 ns. The starting times of the main pulse current in shots 2 023 052 405 and 2 023 052 406 were T_{01} = −100 ns and T_{02} = −80 ns, respectively. The raw projection data and the preprocessed image (pseudocolor image) of angle 1 from shot 2 023 052 405 are shown in Fig. 8. The projection image has four divisions, with an interval of 40 ns between adjacent divisions. Starting from Fig. 8(b), time increases in the clockwise direction.
The number of source voxels is 200 × 200 × 460 and the size of each voxel is about 45 × 45 × 45 µm^{3}. We employ the analytical method and the iterative method with DIP to perform 3D reconstruction from the experimental data. The results of the analytical method, with maximum expansion order M_{a} = 5, are shown in Fig. 9. The analytical algorithm can roughly reconstruct the contour of the radiation area, and some hot spots are visible along the z axis. However, there remain massive burr-like artifacts in the reconstruction results, with the shape and distribution of the hot spots greatly disturbed. The internal 3D radiation distribution is severely degraded. The analytical method is unable to obtain high-quality reconstructed images.
In the iterative algorithm, the maximum expansion order M = 30 and the number of iterations K = 20. The standard deviation of random noise with Gaussian distribution was set to 1/30. The parameters of the network were updated by 2400 iterations using the Adam^{30} optimization method with learning rate set to 0.01. Figure 10 shows 3D spatial distributions of XUV/SR emission from the two shots at four moments reconstructed by the iterative method with DIP. Slice images at z = −47 from shot 2 023 052 405 and at z = 113 from shot 2 023 052 405 are shown in Fig. 11. The plasma is in the Z-pinch run-in phase of accelerating inward implosion at T_{1}. At this time, XUV/SR radiation has already been generated, and there are some hot spots with strong radiation along the z axis. At time T_{2}, which is near the stagnation phase, the plasma is compressed to a small size and forms a thin cylindrical structure, with its radiation reaching maximum. T_{3} and T_{4} are in the plasma collapse phase after stagnation, the radiation weakens, and the 3D distribution gradually becomes irregular. Compared with shot 2 023 052 405, the hotspot in shot 2 023 052 406 is more evenly distributed on the z axis and has a higher overall compression rate. It can be seen from the z-slice images at the strongest radiation spot that near the stagnation phase, the hotspot in shot 2 023 052 406 is closer to a circular shape and is more uniform than that in shot 2 023 052 405.
V. DISCUSSION
The projection angle is not limited to a distribution in the same plane by the proposed algorithm. The iterative method still works when the imaging angle is not horizontally distributed. It is also applicable for other projection modes, such as cone beam. Our algorithm has a good suppressive effect on noisy data. Prior knowledge about plasma and imaging could be added in iteration, such as noise characteristics and image smoothness.
The inclination angle of pinhole imaging in our experiment is ∼1/15, and in the single-wire Z-pinch experiment, the radial size of the plasma pinch process before the stagnation phase with which we are concerned is also very small. Therefore, we adopt the parallel beam imaging mode to approximate the pinhole imaging.^{31} The alignment of multi-axis imaging systems directly affects the registration of imaging data from various angles. Currently, we manually register data from five angles based on the features in projections, which introduces a certain degree of error to 3D reconstruction. Self-absorption on the photons emitted by the plasma in our Z-pinch experiment occurs when the wavelength is in the XUV band. The algorithm currently does not take account of the influence of this attenuation, which is also a possible factor leading to reconstruction errors.
In the future, we will conduct experimental calibration on the center of multi-axis imaging, improve the synchronization accuracy of the system, and calibrate the response of different axis imaging devices to obtain high-precision imaging data. We will consider the impact of attenuation to reduce reconstruction errors and obtain more accurate and reliable information of the reconstructed source.
ACKNOWLEDGMENTS
This paper was supported partially by a grant from NNSFC No. 12027811. We would like to thank the research group of Professor Wu Jian from Xi’an Jiaotong University for providing support and assistance.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Jianpeng Gao: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal). Liang Sheng: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Xinyi Wang: Data curation (equal); Methodology (equal); Visualization (equal). Yanhong Zhang: Data curation (equal); Methodology (equal). Liang Li: Data curation (equal); Funding acquisition (equal); Methodology (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Baojun Duan: Funding acquisition (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Mei Zhang: Data curation (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Yang Li: Formal analysis (equal); Validation (equal). Dongwei Hei: Funding acquisition (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding authors upon reasonable request.