Significant advances have been made in inertial confinement fusion (ICF) ignition techniques such as laser ignition1,2 and Z pinches3–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 Laser7 (x-ray), the National Ignition Facility (NIF)8 (deuterium and tritium neutrons), and the Qin-1 facility6 [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. Pecora13 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.

Expanding the 3D intensity distribution of a pulsed radiation source in cylindrical harmonic function, we getSr=∑m=−MMsm(ρ,z)eimθ,(1)where r = (ρ, θ, z) is a point in the cylindrical coordinate system, M denotes the maximum expansion order of the cylindrical harmonic functions, and s_m are the expansion coefficients. Considering that the source intensities are real numbers, we have the following equation:sm(ρ,z)=s−m*(ρ,z).(2)

Expanding the source with non-negative m values, we haveS(r)=s0(ρ,z)+2∑m=1MResm(ρ,z)cos⁡mθ−2∑m=1MImsm(ρ,z)sin⁡mθ.(3)

Letting sma(ρ,z)=Resm(ρ,z) and smb(ρ,z)=Imsm(ρ,z), we getS(r)=s0(ρ,z)+2∑m=1Msma(ρ,z)cos⁡mθ−2∑m=1Msmb(ρ,z)sin⁡mθ⋅(4)In the iterative method, the expansion coefficients of the cylindrical harmonic function decomposition are obtained assm=arg minsmL(Irec,I),(5)where I denotes the projections of the experimental measurements and I_rec denotes the projections calculated by the reconstructed source. The projections can be expressed asIreci=Pp̂i{S},(6)where Pp̂i{} denotes the projection operation, which represents the projection of the source along the direction p̂i.

The objective function takes the formL(Irec,I)=λ2VNxNy∑i=1VIreci(x,y)−Ii(x,y)22+KsTV,(7)where V is the number of imaging views, N_xN_y is the pixel number of a projection image, ⋅2 denotes the l₂ norm and ⋅TV the total variation,20s is the vector of expansion coefficients, and K is the transform matrix from coefficients to image.

s can be solved by the following optimization problem:mins,dλ2I−As22+d1,s.t.d=∇⃗Ks≡φ(s),(8)where A is the product of the projection matrix P and the transformation matrix K. We convert this into an unconstrained problem:mins,dλ2I−As22+d1+μ2d−φ(s)22.(9)

The iteration steps are as follows, using the alternating direction method of multiplier (ADMM) algorithm:21sk+1=minsλ2I−As22+μk2dk−φ(s)−bk22,(10)dk+1=shrinkφ(sk+1)+bk,1μk,(11)bk+1=bk+φ(sk+1)−dk+1.(12)

The problem (10) is solved by a quasi-Newton method called L-BFGS.22,23 The parameter μ is updated adaptively during the iteration:24μk+1=γμk,Δk+1≥ηΔk,μk,Δk+1<ηΔk,(13)Δk+1≡1n1sk+1−sk22+1n2dk+1−dk22+bk+1−bk22,(14)where γ > 0 and η ∈ [0, 1).

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.ALGORITHM 1.

Iterative method using cylindrical harmonic decomposition.

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-Net25 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.

The loss function is given by the following formula:Lt=α1MSE(St,S)+α2L1(P{St},I)+α3L1(St,St−1)+α4∑z=1Nz−1MSE(St(z),St(z+1)),where α_i ∈ [0, 1] are the parameters, MSE(·) denotes the mean-square error, and L₁(·) denotes the l₁ norm. S_t = f_Θ(η) is the output 3D image of the neural network at the tth iteration, η is the random noise image, f_Θ is the CNN with parameters Θ, S is the image of the iterative reconstructed result, S_t(z) is the 2D image of the output image at layer z, and N_z is the number of layers. The algorithm for reducing artifacts is shown in Algorithm 2.ALGORITHM 2.

Artifact reduction via deep image prior.

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.01, d₀ = 0, b₀ = 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 α₁ = 0.9, α₂ = 0.1, α₃ = 0.1, and α₄ = 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.

The results are shown in Fig. 3. The perspective images show a significant difference between the 3D distribution of the analytical results and the simulation data. The results of the iterative method have vertical artifacts, while results of the iteration with DIP have fewer artifacts and are closer to the simulated data. To evaluate the impact of noisy projection data on reconstruction algorithms, 5% Gaussian noise is added to the projection data, and the reconstruction results of the three algorithms in this case are shown in Fig. 4. It can be seen that both the analytical and iterative methods are subject to strong noise interference, resulting in significant reconstruction errors, while the iteration with the DIP algorithm can maintain a high-quality result. The projections of the reconstructed results are shown in Fig. 5. The analytical results exhibit large errors between the projections of the reconstruction and simulation data, while the results of the iteration with DIP exhibits the smallest projection errors in both data cases. The convergence curve of the iterativeprocess and the loss curve of the DIP are shown in Fig. 6. The iterative and DIP processes converge in both cases. The convergence value is larger with noisy data. Numerical evaluations are used to calculate the indices of the reconstruction results, including structural similarity (SSIM),26 peak signal-to-noise ratio (PSNR), root mean square error between 3D reconstruction results and 3D simulation data [RMSE(S)], and root mean square error between reconstruction result projections and simulation projections [RMSE(I)]. The SSIM index of two 3D images X, Y is defined as:SSIM = [l(X,Y)]α[c(X,Y)]β[s(X,Y)]γwhere l(X,Y), c(X,Y), and s(X,Y) are the luminance comparison function, the contrast comparison function, and the structure comparison function, respectively.  α,β,γ>0 are parameters to adjust the relative importance of these three terms. The results are shown in Table I. The evaluations show that in both noiseless and noisy cases, the results of the iteration with the DIP algorithm are better than those of the iterative method, and remarkably outperform the analytical algorithm.

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₁ = 0 ns, T₂ = 40 ns, T₃ = 80 ns, and T₄ = 120 ns. The starting times of the main pulse current in shots 2 023 052 405 and 2 023 052 406 were T₀₁ = −100 ns and T₀₂ = −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³. 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 Adam30 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₁. 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₂, 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₃ and T₄ 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.

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.

Acknowledgment. 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.