Neutron source three-dimensional image reconstruction method for constraining EM iteration initial value based on column harmonic function decomposition method

Abstract

The invention relates to a pulse neutron source image diagnosis method, in particular to a neutron source three-dimensional image reconstruction method for constraining an EM iteration initial value based on a column harmonic function decomposition method, and solves the technical problem that a high-precision and high-accuracy neutron source three-dimensional image is difficult to obtain by the existing column harmonic function decomposition method and iteration method. The neutron source three-dimensional image reconstruction method for constraining the EM iteration initial value based on the column harmonic function decomposition method comprises the following steps that 1, N neutron source two-dimensional projection images I0 are obtained, the projection angles corresponding to the N neutron source two-dimensional projection images I0 are determined, and N is larger than 1; 2) obtaining N neutron source two-dimensional projection images I; the method comprises the steps of (1) obtaining a column harmonic function reconstruction result, (2) obtaining a column harmonic function reconstruction result, (4) carrying out iteration updating by utilizing an EM algorithm until the difference between a kth iteration result and a (k-1) th iteration result is smaller than 10 <-4 >, stopping iteration and taking a last iteration result as an output result (k is greater than or equal to 1), and (5) recombining the output result into a three-dimensional matrix and taking the recombined three-dimensional matrix as a neutron source reconstruction result.

Description

Technical Field

[0001] The invention relates to a pulse neutron source image diagnosis method, in particular to a neutron source three-dimensional image reconstruction method based on a column harmonic function decomposition method constraining an EM iteration initial value. Background Art

[0002] In the diagnosis of Z-pinch or laser-driven inertial confinement fusion, the three-dimensional information acquisition technology of neutrons produced by the fusion target (neutron source) is a cutting-edge technology developed in recent years. With the help of a multi-axis pinhole camera system, the technology gives a three-dimensional spatial image of the neutron source through image reconstruction algorithm inversion, which can comprehensively and intuitively give the neutron implosion state. Usually, the number of optical axes of a multi-axis pinhole camera system is very limited, resulting in a strong ill-posed integral projection equation.

[0003] The methods for neutron source reconstruction usually use column harmonic decomposition method and iterative methods, such as EM algorithm, SART algorithm, etc.; the column harmonic decomposition method fully utilizes the structural characteristics of inertial confinement fusion target capsule that can be expanded using column harmonic functions, but the expansion order of column harmonic functions is limited by the projection angle, so that its expression ability cannot fully characterize the actual three-dimensional spatial structure of the neutron source, and its reconstruction results cannot ensure that the forward projection is consistent with the actual projection; the iterative method is more universal in image reconstruction, and its effective solution range is wider, so it is difficult for the above two methods to obtain high-precision and high-accuracy three-dimensional images of neutron sources. Summary of the invention

[0004] The purpose of the present invention is to solve the technical problem that both the existing column harmonic function decomposition method and the iteration method are difficult to obtain high-precision and high-accuracy three-dimensional images of neutron sources, and to provide a neutron source three-dimensional image reconstruction method based on the column harmonic function decomposition method to constrain the initial value of EM iteration.

[0005] In order to achieve the above object, the present invention adopts the following technical solution:

[0006] A neutron source three-dimensional image reconstruction method based on column harmonic function decomposition method to constrain EM iteration initial value, which is special in that it includes the following steps:

[0007] 1] Use a multi-axis pinhole camera system located on the same horizontal plane to obtain N neutron source two-dimensional projection images I0, and determine their corresponding projection angles, N>1;

[0008] 2] Preprocess all neutron source two-dimensional projection images I0 to obtain N neutron source two-dimensional projection images I;

[0009] 3] Generate and solve a linear equation system based on the projection angle and N neutron source two-dimensional projection images I to obtain the column harmonic function reconstruction result;

[0010] 4】Use the reconstruction result of cylindrical harmonic functions as the initial value of the EM algorithm, and iteratively update it using the EM algorithm until the difference between the k-th iteration result and the (k - 1)-th iteration result is less than 10 -4 and then stop the iteration, and use the last iteration result as the output result, where k ≥ 1;

[0011] 5】Recombine the output result into a three-dimensional matrix, and use the recombined three-dimensional matrix as the reconstruction result of the neutron source to complete the three-dimensional image reconstruction of the neutron source.

[0012] Further, step 2】is specifically as follows:

[0013] 2.1. Respectively perform pixel size correction on each two-dimensional projection image I0 of the neutron source to obtain N pixel-corrected projection images I1:

[0014] 2.2. Respectively perform shearing on the N pixel-corrected projection images I1 so that the image centers coincide with the neutron source center to obtain N sheared images I2;

[0015] 2.3. Respectively perform median filtering on the N sheared images I2 to obtain N two-dimensional projection images I of the neutron source.

[0016] Further, step 3】is specifically as follows:

[0017] 3.1. Generate a system of linear equations based on the projection angles obtained in step 1】and the N two-dimensional projection images I of the neutron source obtained in step 2.3, and solve this system of linear equations using the analytical method to obtain the Hankel transform value of the cylindrical harmonic function expansion coefficient s m

[0018]

[0019] where M is the maximum value of the cylindrical harmonic function expansion order, m is the cylindrical harmonic function expansion order, and m = 1, 2,..., M; is the Hankel transform value of the expansion coefficient s0 corresponding to m = 0 of the cylindrical harmonic function expansion; ξ is the spatial frequency after Fourier transform; z is a variable; is the projection angle of the eigenfunction; is the Fourier transform value of the two-dimensional projection image I of the neutron source; Re{...} represents taking the real part, and Im{...} represents taking the imaginary part;

[0020] 3.2. Perform inverse Hankel transform on the Hankel transform value m of the cylindrical harmonic function expansion coefficient s to obtain the cylindrical harmonic function expansion coefficient s m ;

[0021] ​3.3. Substitute the expansion coefficient s of the cylindrical harmonic function m into the expansion formula of the cylindrical harmonic function to obtain the reconstruction result of the cylindrical harmonic function S(ρ, θ, z); then convert the reconstruction result of the cylindrical harmonic function S(ρ, θ, z) into the representation in the rectangular coordinate system to obtain the reconstruction result of the cylindrical harmonic function S(x, y, z); where ρ is the projection radius of the pixel point on the neutron source on the z-plane, and θ is the angle between the projection radius and the x-axis.

[0022] Furthermore, step 4 is specifically as follows:

[0023] 4.1. In the Matlab software, express the N two-dimensional projection images I of the neutron source by the integral equation WX = Y, and convert the reconstruction result of the cylindrical harmonic function S(x, y, z) obtained in step 3.3 into a one-dimensional vector and assign it to the initial value X of the iteration of the EM algorithm 0 :

[0024] WX = Y

[0025] X 0 = reshape(S(x, y, z), 1, []);

[0026] In the formula, X is the one-dimensional representation form of the neutron source; W is the projection matrix generated by the projection angle and the pixel size of the N two-dimensional projection images I of the neutron source; Y is the one-dimensional vector converted from the N two-dimensional projection images I of the neutron source; reshape is the matrix conversion code;

[0027] 4.2. Iteratively update X through the following formula. When the difference between the k-th iteration result X k and the (k - 1)-th iteration result X k-1 is less than 10 -4 , stop the iteration, and take X k as the final output result;

[0028]

[0029] In the formula, represents the j-th value of the X k vector; w t,j represents the value of the t-th row and j-th column of the projection matrix W; w l,j represents the value of the l-th row and j-th column of the projection matrix W, where l = 1, 2,..., T, and T is the number of rows of the projection matrix W; y l is the l-th value of the Y vector, and W l is the l-th row array of the projection matrix W.

[0030] Furthermore, it also includes step 6:

[0031] Use a 3D program to draw the 3D structure of the neutron source reconstruction result.

[0032] Further, step 2.1 is specifically as follows:

[0033] Use the program code in Matlab software to adjust the pixel size of the two-dimensional projection images of N neutron sources according to the standard that the area corresponding to each pixel point to the neutron source is 1 micron × 1 micron, and obtain the pixel-corrected projection image I1 through the following formula:

[0034] I1 = imresize(I n0 , scale, 'bilinear')

[0035]

[0036] In the formula, imresize is the pixel adjustment statement; scale is the pixel correction ratio between the pixel-corrected projection image and the two-dimensional projection image of the neutron source; 'bilinear' is the pixel intensity difference method, and linear interpolation is used in this embodiment; I n0 represents the two-dimensional projection image I0 of the neutron source on the nth optical axis, where n is 1, 2,... N; a n represents the area corresponding to the pixel of the two-dimensional projection image I0 of the neutron source obtained on the nth optical axis.

[0037] Further, step 2.3 is specifically as follows:

[0038] Use the program code in Matlab software to perform median filtering on the N shear images I2 respectively through the following formula to obtain the two-dimensional projection images I of N neutron sources:

[0039] I n = medfilt2(I n2 , [3 3])

[0040] In the formula, medfilt2 is the two-dimensional image median filtering statement; [3 3] represents that the median filtering window size is 3×3; I n2 represents the shear image I2 on the nth optical axis; I n represents the two-dimensional projection image of the neutron source on the nth optical axis.

[0041] Further, in step 1], N = 3.

[0042] Further, in step 2.2, the pixel sizes of the N shear images I2 are the same, and both the row and column pixels are even numbers.

[0043] Advantages of the present invention:

[0044] 1) The three-dimensional image reconstruction method of neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method combines the advantages of the cylindrical harmonic function decomposition method and the EM iteration method, can well eliminate the rotational artifacts of the cylindrical harmonic function and the striped artifacts of the EM algorithm, and the obtained three-dimensional image of the neutron source has higher accuracy.

[0045] 2) The three-dimensional image reconstruction method of neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method can overcome the ill-posedness and reconstruct a three-dimensional map with high precision in the case of only a few projection images.

[0046] 3) The three-dimensional image reconstruction method of neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method combines the advantages of two reconstruction methods, namely the cylindrical harmonic function decomposition method and the iterative method. It proposes a reconstruction method that uses the cylindrical harmonic function decomposition method as the initial value of the EM iteration. First, the three-dimensional space function of the neutron source is preliminarily given by the cylindrical harmonic function decomposition method, and then the EM algorithm is used to optimize the reconstruction result of the cylindrical harmonic function to find the optimal solution in the neighborhood of the reconstruction result of the cylindrical harmonic function, further improving the accuracy of the three-dimensional structure of the neutron radiation area of the fusion pellet.

[0047] 4) The three-dimensional image reconstruction method of neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method can still obtain good reconstruction results with very few projection angles, and the number of projection angles can be as few as 2 to 5.

[0048] 5) The three-dimensional image reconstruction method of neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method has a wide range of applications. It can be applied to self-luminous radiators with a certain symmetric structure, such as Z-pinch implosion experiments, etc.

[0049] 6) The three-dimensional image reconstruction method of neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method uses a multi-optical axis pinhole imaging system located on the same horizontal plane to obtain the three-dimensional intensity distribution information of the neutron source, and can be applied to verifying the compression symmetry of the pellet, verifying physical models, evaluating the success or failure of ignition, etc. in the diagnostic research of inertial confinement fusion devices.

[0050] 7) The three-dimensional image reconstruction method of neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method takes the reconstruction result of the cylindrical harmonic function of the neutron source as the initial value in the iterative process of the EM algorithm, and uses the integral projection equation to obtain the three-dimensional distribution of the neutron source. This method constrains the reconstruction result of the EM iteration algorithm near the analytical result of the cylindrical harmonic function, and fully utilizes the prior knowledge that the inertial confinement fusion target is a ray source with a certain spherical symmetry. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 is a flowchart of an embodiment of the three-dimensional image reconstruction method of neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method;

[0052] Figure 2 It is a schematic diagram of a two-dimensional projection image of a neutron source obtained by using a multi-optical axis pinhole imaging system;

[0053] Figure 3 It is a schematic diagram of the two-dimensional projection image I of the neutron source corresponding to the nth optical axis in the embodiment of the three-dimensional image reconstruction method of the neutron source based on the constraint of the EM iteration initial value by the cylindrical harmonic function decomposition method;

[0054] Figure 4(a) is a perspective view of the simulated three-dimensional structure of the neutron source in the embodiment of the three-dimensional image reconstruction method of the neutron source based on the constraint of the EM iteration initial value by the cylindrical harmonic function decomposition method;

[0055] Figure 4(b) is a perspective view of the reconstructed three-dimensional structure of the neutron source in the embodiment of the three-dimensional image reconstruction method of the neutron source based on the constraint of the EM iteration initial value by the cylindrical harmonic function decomposition method;

[0056] Figure 5(a) is a three-dimensional structure diagram of the simulated three-dimensional structure of the neutron source in the embodiment of the three-dimensional image reconstruction method of the neutron source based on the constraint of the EM iteration initial value by the cylindrical harmonic function decomposition method;

[0057] Figure 5(b) is a three-dimensional structure diagram of the reconstructed three-dimensional structure of the neutron source in the embodiment of the three-dimensional image reconstruction method of the neutron source based on the constraint of the EM iteration initial value by the cylindrical harmonic function decomposition method;

[0058] Figure 6(a) is a schematic diagram of the Z-axis section of the simulated three-dimensional structure of the neutron source in the embodiment of the three-dimensional image reconstruction method of the neutron source based on the constraint of the EM iteration initial value by the cylindrical harmonic function decomposition method;

[0059] Figure 6(b) is a schematic diagram of the Z-axis section of the reconstructed three-dimensional structure of the neutron source in the embodiment of the three-dimensional image reconstruction method of the neutron source based on the constraint of the EM iteration initial value by the cylindrical harmonic function decomposition method. Detailed implementation manners

[0060] As Figure 1 shown, a three-dimensional image reconstruction method of a neutron source based on the constraint of the EM iteration initial value by the cylindrical harmonic function decomposition method includes the following steps:

[0061] 1. Use a multi-optical axis pinhole imaging system located on the same horizontal plane to obtain N two-dimensional projection images I0 of the neutron source, and determine their corresponding projection angles N>1; the multi-optical axis pinhole imaging system determines the number of optical axes, which is the same as the number of two-dimensional projection images I0 of the neutron source. The pixel size of the two-dimensional projection image of the neutron source obtained by each optical axis and the area corresponding to the pixel to the neutron source;

[0062] As Figure 2As shown, a neutron source in an asymmetric Gaussian-like three-dimensional structure is established, N = 3. Using three pinhole imaging systems on the same horizontal plane, three two-dimensional projection images of the neutron source are obtained; each pixel is 100*100, and the pixel size of the two-dimensional projection image of the neutron source corresponds to an area of 1 square micron of the neutron source. The projection angles of each optical axis are respectively (90, 0), (90, 120), and (90, 240).

[0063] 2. Preprocess all the two-dimensional projection images I0 of the neutron source to obtain N two-dimensional projection images I of the neutron source. Specifically:

[0064] 2.1. Perform pixel size correction on each two-dimensional projection image I0 of the neutron source respectively to obtain N pixel-corrected projection images I1:

[0065] Use the program code in Matlab software to adjust the pixel size of the N two-dimensional projection images of the neutron source according to the standard that the area corresponding to each pixel point to the neutron source is 1 micron × 1 micron, and obtain the pixel-corrected projection image I1 through the following formula:

[0066] I1 = imresize(I n0 , scale, 'bilinear')

[0067]

[0068] In the formula, imresize is the pixel adjustment statement; scale is the pixel correction ratio between the pixel-corrected projection image and the two-dimensional projection image of the neutron source, and the value in this embodiment is 1; 'bilinear' is the pixel intensity difference method, and linear interpolation is adopted in this embodiment; I n0 represents the two-dimensional projection image I0 of the neutron source on the nth optical axis, where n is 1, 2,... N; a n represents the area corresponding to the pixel of the two-dimensional projection image I0 of the neutron source obtained on the nth optical axis;

[0069] 2.2. Perform shearing on the N pixel-corrected projection images I1 respectively to make the image center coincide with the neutron source center, and obtain N sheared images I2;

[0070] Preferably, it is required that the pixel sizes of the N sheared images I2 are the same, and the number of rows and columns of pixels are both even;

[0071] 2.3. Use the program code in Matlab software to perform median filtering on the N sheared images I2 respectively through the following formula to obtain N two-dimensional projection images I of the neutron source:

[0072] I n = medfilt2(I n2 , [3 3])

[0073] In the formula, medfilt2 is the median filtering statement for a two-dimensional image; [3 3] indicates that the median filtering window size is 3×3; I n2 represents the sheared image I2 of the nth optical axis;

[0074] 3】 Obtain the linear system of cylindrical harmonic functions of the two-dimensional projection images I of N neutron sources and solve it to obtain the cylindrical harmonic function reconstruction result, specifically:

[0075] 3.1. According to the projection angles obtained in step 1】 and the N two-dimensional projection images I of neutron sources obtained in step 2.3, generate a linear system of equations, and solve this linear system of equations using the analytical method to obtain the Hankel transform value m of the cylindrical harmonic function expansion coefficient s

[0076]

[0077] In the formula, M is the maximum value of the cylindrical harmonic function expansion order, m is the cylindrical harmonic function expansion order, m = 1, 2,... M; is the Hankel transform value of the expansion coefficient s0 corresponding to m = 0 of the cylindrical harmonic function expansion order; ξ is the spatial frequency after Fourier transform; z is a variable; is the projection angle of the eigenfunction; is the Fourier transform value of the two-dimensional projection image I of the neutron source; Re{...} represents taking the real part, and Im{...} represents taking the imaginary part;

[0078] The representation method of the neutron source in the cylindrical coordinate system is S(ρ, θ, z), where ρ is the projection radius of the pixel point on the neutron source on the Z plane, θ is the angle between the projection radius and the x-axis, that is, the projection angle of the nth optical axis z is a variable, and S(ρ, θ, z) is expressed by the cylindrical harmonic function expansion formula as:

[0079]

[0080] Among them, e imθ is the eigenfunction of the angle θ between the projection radius and the x-axis. Since the neutron source is an unknown quantity, the Hankel transform value m of the cylindrical harmonic function expansion coefficient s is the parameter to be solved in this system of equations.

[0081] Such as Figure 3As shown, the neutron source can be formed by stacking two-dimensional images of the Z-section layer by layer. A rectangular coordinate system with O as the origin of coordinates is established. The representation method of the neutron source in the rectangular coordinate system is S(x, y, z); a plane coordinate system (u, z) with O' as the origin of coordinates is established on the projection plane. Then, the two-dimensional projection image I of the neutron source corresponding to this optical axis is represented on the plane coordinate axes as I n (u, z), is the one-dimensional Fourier transform of I n (u, z) along the u direction. Using the N two-dimensional projection images I of all optical axes to form a system of linear equations, which is formula (1); when the maximum value M of the order of the cylindrical harmonic function expansion of the neutron source is N - 1, the system of equations in formula (1) is full rank, and the Hankel transform value of the cylindrical harmonic function expansion coefficient s m has only one valid solution.

[0082] 3.2. Perform the inverse Hankel transform on the Hankel transform value of the cylindrical harmonic function expansion coefficient s m to obtain the cylindrical harmonic function expansion coefficient s ; m

[0083] 3.3. Substitute the cylindrical harmonic function expansion coefficient s m into the cylindrical harmonic function expansion formula to obtain the reconstruction result of the cylindrical harmonic function S(ρ, θ, z); then convert the reconstruction result of the cylindrical harmonic function S(ρ, θ, z) into the representation form in the rectangular coordinate system to obtain the reconstruction result of the cylindrical harmonic function S(x, y, z);

[0084] 4. Take the reconstruction result S(x, y, z) of the cylindrical harmonic function as the initial value of the EM algorithm iteration, and use the EM algorithm to iteratively update until the difference between the k-th iteration result and the (k - 1)-th iteration result is less than 10 -4 and stop the iteration, and take the last iteration result as the output result;

[0085] 4.1. In the Matlab software, express the N two-dimensional projection images I of the neutron source with integral equations, and convert the reconstruction result S(x, y, z) of the cylindrical harmonic function obtained in step 3.3 into a one-dimensional vector and assign it to the initial value X of the EM algorithm iteration 0 , specifically:

[0086] WX = Y

[0087] X 0 = reshape(S(x, y, z), 1, []);

[0088] In the formula, X is the one-dimensional expression form of the neutron source, which is the unknown; W is the projection angle composed of known parameters The projection matrix generated from the pixel sizes of the two-dimensional projection images I of N neutron sources; Y is the one-dimensional vector converted from the two-dimensional projection images I of N neutron sources; X 0 is the initial value of the iteration of the EM algorithm; reshape is the matrix conversion code;

[0089] 4.2. Iteratively update X through formula (2). When the result X of the k-th iteration k differs from the result X of the (k - 1)-th iteration k-1 by less than 10 -4 , stop the iteration and take X k as the final output result;

[0090]

[0091] In the formula, represents the j-th value of the X k vector; w t,j represents the value at the t-th row and j-th column of the projection matrix W; w l,j represents the value at the l-th row and j-th column of the projection matrix W, where l = 1, 2, …, T, and T is the number of rows of the projection matrix W; y l is the l-th numerical value of the Y vector, and W l is the l-th row array of the projection matrix W;

[0092] 5】 Recombine the final output result X k into a three-dimensional matrix Q(x, y, z), and take the recombined three-dimensional matrix Q(x, y, z) as the neutron source reconstruction result to complete the three-dimensional image reconstruction of the neutron source.

[0093] 6】 Use a three-dimensional program to draw the three-dimensional structure of the recombined three-dimensional matrix Q(x, y, z).

[0094] As Figures 4(a) to 6(b) shown, Figure 4(a) , 4(b) are the perspective view of the simulated three-dimensional structure of the neutron source and the perspective view of the reconstructed three-dimensional structure of the neutron source in the embodiment of the three-dimensional image reconstruction method of the neutron source based on the constraint of the cylindrical harmonic function decomposition method for the EM iteration initial value; Figure 5(a) , 5(b) are the three-dimensional structure diagram of the simulated three-dimensional structure of the neutron source and the three-dimensional structure diagram of the reconstructed three-dimensional structure of the neutron source in the embodiment of the three-dimensional image reconstruction method of the neutron source based on the constraint of the cylindrical harmonic function decomposition method for the EM iteration initial value; Figure 6(a) , 6(b) are the schematic diagrams of the Z-axis cross-sections of the simulated three-dimensional structure of the neutron source and the reconstructed three-dimensional structure of the neutron source; It can be seen from Figures 4(a) to 6(b) that the similarity between the simulated neutron source and the reconstructed neutron source is 0.9816.

[0095] In the research on neutron pinhole image diagnosis in the deuterium-tritium combustion region of inertial confinement fusion, the acquisition mechanism of the two-dimensional projection image of the neutron source is the integration result of the neutron rays in the combustion region along the pinhole direction. It is similar to the principle of computer tomography. By establishing an integral equation using N two-dimensional projection images I of the neutron source and solving this equation, the three-dimensional spatial result of the neutron source can be obtained.

[0096] When a multi-axis pinhole camera system located on the same horizontal plane (assuming the Z-axis) obtains two-dimensional projection images of the neutron source in different directions at the same time, according to the structural characteristics of the Z-pinch or laser-driven inertial confinement fusion target pellet (i.e., the neutron source) with a Gaussian-like distribution, the Z-sectional image of the three-dimensional structure neutron source can be expanded into a set of orthogonal cylindrical harmonic function expressions. By introducing this cylindrical harmonic function expression into the integral equation of the neutron two-dimensional projection image and using Fourier transform and the central slice theorem, the integral equation of the neutron two-dimensional projection image is transformed from a linear integral equation system in the spatial domain to a linear equation system in the frequency domain. Among them, the two-dimensional projection image of the neutron source is Fourier-transformed into an expression in the spatial domain, and the unknown to be solved is the Hankel transform of the expansion coefficients of the cylindrical harmonic function of the neutron source. This method aims at the situation where the number of optical axes of the multi-axis pinhole camera system is small and the uncertainty of the reconstruction result of the linear integral equation system in the spatial domain is strong. By converting it into a linear equation system in the frequency domain, the effective solution obtained by the analytical method of this equation is unique. Performing an inverse Hankel transform on the obtained effective solution results in the expansion coefficients of the cylindrical harmonic function at different expansion orders of the cylindrical harmonic function, and linearly combining the expansion coefficients of the cylindrical harmonic function gives the reconstruction result of the cylindrical harmonic function. Using this reconstruction result of the cylindrical harmonic function as the initial value of the EM algorithm, the integral equation of the neutron two-dimensional projection image is iteratively solved, and finally the optimal result is obtained as the three-dimensional intensity distribution data of the neutron source.

[0097] In this embodiment, the diagnostic object is the three-dimensional spatial distribution information of the neutron radiation intensity of the neutron source. All the optical axes of the multi-axis pinhole camera system need to be located on the same horizontal plane, and the angles are expressed in polar coordinates (pitch angle, projection angle), where the opening times of the two-dimensional projection images of the neutrons obtained in different directions are the same. In the method of the present invention, the number of pinhole imaging angles is extremely small, being 2 - 5; the order of the cylindrical harmonic function expansion is 1 less than the number of projection angles; when the difference between the k-th and (k - 1)-th iterative reconstruction results is less than 10 -4 the EM algorithm stops iterating.

Claims

1. A three-dimensional image reconstruction method for neutron sources that constrains the initial value of EM iteration based on the cylindrical harmonic function decomposition method, characterized in that, It includes the following steps:

1. Use a multi-optical-axis pinhole imaging system on the same horizontal plane to obtain N two-dimensional projection images I0 of the neutron source, and determine their corresponding projection angles, where N > 1; 2. Preprocess all the two-dimensional projection images I0 of the neutron source to obtain N two-dimensional projection images I of the neutron source; 3. Generate and solve a system of linear equations based on the projection angles and the N two-dimensional projection images I of the neutron source to obtain the reconstruction result of the cylindrical harmonic function; 4】Use the reconstruction result of the cylindrical harmonic function as the initial value of the EM algorithm, and iteratively update it using the EM algorithm until the difference between the k-th iteration result and the (k - 1)-th iteration result is less than 10 -4 to stop the iteration, and use the last iteration result as the output result, where k ≥ 1; 5. Recombine the output result into a three-dimensional matrix, and use the recombined three-dimensional matrix as the reconstruction result of the neutron source to complete the three-dimensional image reconstruction of the neutron source.

2. The three-dimensional image reconstruction method of neutron source based on the constraint of EM iteration initial value by column harmonic function decomposition method according to claim 1, characterized in that, Step 2 is specifically as follows: 2.

1. Correct the pixel size of each two-dimensional projection image I0 of the neutron source respectively to obtain N pixel-corrected projection images I1: 2.

2. Perform shearing on the N pixel-corrected projection images I1 respectively to make the image center coincide with the neutron source center to obtain N sheared images I2; 2.

3. Perform median filtering on the N sheared images I2 respectively to obtain N two-dimensional projection images I of the neutron source.

3. The three-dimensional image reconstruction method of neutron source based on the constraint of the initial value of EM iteration by column harmonic function decomposition method according to claim 2, characterized in that, Step 3 is specifically as follows: 3.

1. Generate a system of linear equations based on the projection angles obtained in step 1 and the N two-dimensional projection images I of the neutron sources obtained in step 2.3, and solve this system of linear equations using the analytical method to obtain the expansion coefficients s of the cylindrical harmonic function m Hankel transform value of In the formula, M is the maximum value of the order of the cylindrical harmonic function expansion, m is the order of the cylindrical harmonic function expansion, and m = 1, 2, … M; is the Hankel transform value of the expansion coefficient s0 corresponding to m = 0 of the cylindrical harmonic function expansion; ξ is the spatial frequency after Fourier transform; z is a variable; e imφ is the projection angle is the eigenfunction of is the Fourier transform value of the two-dimensional projection image I of the neutron source; Re{...} represents taking the real part, and Im{...} represents taking the imaginary part; 3.

2. Hankel transform value of the expansion coefficient s of cylindrical harmonic functions m Perform the inverse Hankel transform on the Hankel transform value to obtain the expansion coefficient s of cylindrical harmonic functions m ; 3.

3. Substitute the expansion coefficient s of the cylindrical harmonic function m into the expansion formula of the cylindrical harmonic function to obtain the reconstruction result of the cylindrical harmonic function S(ρ, θ, z); then convert the reconstruction result of the cylindrical harmonic function S(ρ, θ, z) into the representation in the rectangular coordinate system to obtain the reconstruction result of the cylindrical harmonic function S(x, y, z); where ρ is the projection radius of the pixel point on the neutron source on the z-plane, and θ is the angle between the projection radius and the x-axis.

4. The three-dimensional image reconstruction method of a neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method according to claim 3, characterized in that Step 4 is specifically as follows: 4.

1. In Matlab software, represent the two-dimensional projection images I of N neutron sources by the integral equation WX = Y, and after converting the reconstructed result of the cylindrical harmonic function S(x, y, z) obtained in step 3.3 into a one-dimensional vector, assign it to the initial value X for iteration of the EM algorithm 0 : WX = Y X 0 = reshape(S(x,y,z),1,[]); where X is the one-dimensional representation of the neutron source; W is the projection angle and the projection matrix generated by the pixel sizes of the two-dimensional projection images I of N neutron sources; Y is the one-dimensional vector converted from the two-dimensional projection images I of N neutron sources; reshape is the matrix conversion code; 4.

2. Update X iteratively using the following formula. When the result X of the k-th iteration k differs from the result X of the (k - 1)-th iteration k-1 by less than 10 -4 , stop the iteration and use X k as the final output result; In the formula, represents the j-th value of the X k vector; w t,j represents the value of the t-th row and j-th column of the projection matrix W; w l,j represents the value of the l-th row and j-th column of the projection matrix W, where l = 1, 2, …, T and T is the number of rows of the projection matrix W; y l is the l-th numerical value of the Y vector, and W l is the l-th row array of the projection matrix W.

5. The three-dimensional image reconstruction method of a neutron source based on constraining the initial value of EM iteration by the cylindrical harmonic function decomposition method according to claim 4, characterized in that It further includes step 6: Use a three-dimensional program to draw the three-dimensional structure of the neutron source reconstruction result.

6. The three-dimensional image reconstruction method of a neutron source based on the constraint of the initial value of EM iteration by the column harmonic function decomposition method according to claim 5, characterized in that Step 2.1 is specifically as follows: Use the program code in Matlab software to adjust the pixel size of the N two-dimensional projection images of the neutron source according to the standard that the area corresponding to each pixel point to the neutron source is 1 micron × 1 micron through the following formula to obtain the pixel-corrected projection image I1: I1 = imresize(I n0 , scale, 'bilinear') In the formula, imresize is the pixel adjustment statement; scale is the pixel correction ratio between the pixel-corrected projection image and the two-dimensional projection image of the neutron source; 'bilinear' is the pixel intensity difference method, and linear interpolation is adopted in this embodiment; I n0 represents the two-dimensional projection image I0 of the neutron source on the nth optical axis, where n = 1, 2,..., N; a n represents the area of the neutron source corresponding to the pixels of the two-dimensional projection image I0 obtained on the nth optical axis.

7. The three-dimensional image reconstruction method of neutron source based on the constraint of the initial value of EM iteration by column harmonic function decomposition method according to claim 6, characterized in that Step 2.3 is specifically as follows: Use the program code in Matlab software to perform median filtering on the N sheared images I2 respectively through the following formula to obtain N two-dimensional projection images I of the neutron source: I n = medfilt2(I n2 , [3 3]) In the formula, medfilt2 is the median filtering statement for two-dimensional images; [3 3] indicates that the median filtering window size is 3×3; I n2 represents the sheared image I2 of the nth optical axis; I n represents the two-dimensional projection image of the neutron source of the nth optical axis.

8. The three-dimensional image reconstruction method of the neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method according to claim 7, characterized in that: In step 1, N = 3.

9. The three-dimensional image reconstruction method of the neutron source based on the constraint of the initial value of EM iteration by the cylindrical harmonic function decomposition method according to claim 8, characterized in that: In step 2.2, the pixel sizes of the N sheared images I2 are the same, and the number of rows and columns of pixels are both even.

Images