A method of X-ray tomography image reconstruction based on fuzzy PID control

Abstract

The present invention designs an X-ray tomography image reconstruction method based on fuzzy PID control, which mainly includes the following steps: first, using an X-ray tomography system whose ray source is a fan beam, the projection value at different projection angles is obtained by measurement, the projection value is mapped to a parallel beam, and then the initial reconstruction value is calculated by Radon inverse transform; then, the above reconstruction algorithm is used as the controlled object, and the density distribution inverted by the algorithm is used as the feedback quantity, and the parameters of the PID controller are adjusted by fuzzy rules to correct the deviation of the reconstructed density value, so that the entire closed loop converges; finally, if the number of iterations reaches the preset number of iterations, the iteration is terminated and the reconstruction result is output, otherwise the iteration continues. This method proposes a method for mapping fan beam projection data into a parallel beam and combining it with the closed-loop control principle, which can effectively reduce image reconstruction errors and significantly improve image reconstruction quality. It has practical value and good application prospects in the field of X-ray imaging.

Description

Technical Field

[0001] The present invention relates to an X-ray tomography image reconstruction method based on fuzzy PID control, which belongs to the field of X-ray tomography. The fan beam projection data is directly reconstructed by using the inverse Radon transform and then mapped to the parallel beam point by point, and combined with the fuzzy PID closed-loop control. The method can effectively reduce the reconstruction error through the closed-loop control and significantly improve the image reconstruction quality. Background Art

[0002] The advent of X-ray Computer Tomography (XCT) technology is another epoch-making breakthrough in the engineering field after the discovery of radiation, and is recognized as one of the ten major technologies that promoted human development in the 20th century. XCT performs multi-angle radiation scanning on the object to be measured from the outside, thereby accurately and intuitively reconstructing the distribution of the internal tissue of the object's cross section, greatly enhancing the ability of humans to observe the internal structure of the object. This technology is non-invasive, fast, and highly accurate. Since the first XCT machine was introduced in 1971, it has been widely used in medical diagnosis, industrial testing, safety inspection, aerospace, and national defense.

[0003] The two most important components of the CT detection system are the hardware system and the software algorithm. The most important components of the hardware system include the radiation source, the rotating stage and the detector for receiving the radiation. The software part mainly includes three aspects: CT scanning data acquisition, CT image reconstruction and CT image application, among which CT image reconstruction is the core. The reconstruction algorithm has the most direct and huge impact on the reconstruction result, so the main way to improve the CT reconstruction effect is to optimize and improve the reconstruction algorithm. Currently, the commonly used reconstruction algorithms include analytical algorithms (AR) and algebraic algorithms (IR).

[0004] Scholars at home and abroad have conducted a lot of research on image reconstruction algorithms. The analytical method is mainly based on the filtered back projection algorithm. For example, in 2014, Pelt, Batenburg et al. published an article titled "Improving filtered backprojection reconstruction by data-dependent filtering" in the 11th issue of Volume 23 of IEEE Transactions on Image Processing. They used custom data-dependent filters to improve the filtered back projection method to minimize the reconstruction error. The final reconstruction accuracy is similar to that of the algebraic method, but the computational complexity is smaller and the generation speed is faster. It can also be extended to utilize certain forms of prior knowledge to improve the reconstruction accuracy in specific cases. X-ray scanning methods can be divided into fan beam, cone beam and parallel beam. In practice, the commonly used fan beam projection data reconstruction methods can be roughly divided into two categories. One is the direct reconstruction of fan beam projection. This method performs filtered back projection reconstruction after appropriate weighting. For example, Chen GH published "A new framework of image reconstruction from fan beam projections" in "Medical Physics" Vol. 30, No. 6, pp. 1151-1161 in 2003. A new framework was proposed to reconstruct fan beam projection data, and more stringent requirements were set for the reconstruction scanning trajectory. The other is the rearrangement algorithm, which recombines fan beam projection data into parallel beam projection data and reconstructs according to the parallel beam filtered back projection algorithm. For example, Deepak Bharkhada et al. published "A comparative study on interpolation methods for controlled cardiac CT" in "International Journal of Imaging Systems and Technology" Vol. 17, No. 2, pp. 91-98 in 2007. In the article "Controlled Cardiac CT", the fan beam data was rearranged based on the nearest neighbor interpolation method of triangles, and then reconstructed using the filtered back projection algorithm. The proposed method has good noise resistance. However, rearranging the fan beam into parallel beams requires interpolation of the original data, which will introduce errors and result in loss of spatial resolution.

[0005] Radon transform and inverse Radon transform are the basis of tomography, and propose the relationship between the distribution of two-dimensional objects and one-dimensional projection. In 1987, G. Beylkin et al. published a paper titled "Discrete Radon transform" in IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. 35, No. 2, pp. 162-172, introducing the discrete Radon transform and its precise inversion algorithm. The precise reversibility, flexibility of definition, and fast calculation of the discrete Radon transform make it possible to become a new application. In 2012, Amir Averbuch et al. published an article titled "CT Reconstruction From Parallel and Fan-Beam Projections by a 2-D Discrete Radon Transform" in IEEE Transactions on Image Processing, Vol. 21, No. 2, pp. 733-741. They verified that for both parallel beams and equidistant fan beams, the original distribution can be accurately reconstructed by using a 2-D discrete Radon transform on the projection data, and that the reconstruction accuracy improves accordingly as the amount of projection data increases.

[0006] The fuzzy PID control system has a closed-loop control system with a negative feedback channel, and realizes real-time control of the control target by adaptively selecting the control coefficient through fuzzy rules. In 2016, Zafer Civelek et al. published an article entitled "Control of Pitch Angle of Wind Turbine by Fuzzy Pid Controller" in Intelligent Automation & Soft Computing, Vol. 22, No. 3, pp. 1-9. They used conventional PID, fuzzy control and fuzzy PID to control the pitch angle of wind turbines. By evaluating the steady-state time and steady-state error of the output power, the results clearly show that the performance of the fuzzy PID controller is better than that of the ordinary PID and fuzzy controllers. In 2017, Yuzhong Wang et al. published an article entitled "Improved fuzzy PID controller design using predictivefunctional control structure" in ISA Transactions, Vol. 71, No. 2, pp. 354-363. They proposed a new PID control method that inherits the advantages of fuzzy PID control and predictive functional control. First, the current process input signal is used to obtain a prediction of future process behavior, and then multi-step predictive fuzzy PID control is introduced to obtain optimal control. The effectiveness of the control scheme is verified.

[0007] When using fan beam inverse Radon transform for image reconstruction, the loss of fan beam projection data converted into parallel beam data by rearrangement algorithm is relatively large, so the number of projection angles and the number of projection data are required to be high. When the projection data is insufficient, the reconstruction result has large artifacts. Based on the above background, the present invention proposes a new method for mapping fan beam data to parallel beams, and combines it with a fuzzy PID closed-loop control X-ray density reconstruction method. The fuzzy PID closed-loop control is combined with the improved fan beam inverse Radon transform, and the PID controller parameters are adaptively adjusted by the fuzzy controller. The closed-loop is iterated multiple times to continuously reduce the reconstruction error, and finally improve the accuracy of the reconstruction result. Summary of the invention

[0008] The purpose of the present invention is to propose an X-ray tomography based on the fuzzy PID closed-loop control principle. The method has high reconstruction accuracy and can effectively improve the reconstruction quality of X-ray tomography images. The reconstruction system includes an X-ray source, a stage for the object to be measured, and a sensor for detection and reception. Among them, the sensor is a regular s-gon, and s≥3, the X-ray source is located at the s vertices of the s-gon, the detection plate is located at the s sides of the s-gon, and d detectors are evenly distributed on each side. The rays of each X-ray source are fan-shaped beams, covering only the detection plates on the opposite sides, and the sensor arranges a total of Q=sd rays. The proposed reconstruction method first calculates the intersection set between the fan-shaped beam rays and the intersection set between the imaginary parallel beams, calculates the reconstruction value of the fan-shaped beam intersection set by Radon inverse transform, and then maps it to the parallel beam intersection set, and then calculates the density value of the reconstructed area, and finally introduces the fuzzy PID closed-loop model for error control, thereby obtaining the final reconstructed two-dimensional density distribution value. Specifically including the following steps:

[0009] Step 1: Select an initial X-ray source. The projection angle is set to 0. The corresponding projection data on the detection board is:

[0010]

[0011] The projection angle of the next ray source along the clockwise direction is Similarly, the corresponding projection data on the detection board is p2. Similarly, the projection angle of the last ray source, that is, the sth ray source, is The projection data on the detection board is p s , so the projection data matrix of all angles is:

[0012]

[0013] At a certain projection angle, the fan-shaped ray beams with equal angle distribution are parallel in the imaginary space. Corresponding to the imaginary parallel beams with equal interval distribution, when the projection angle is β1=β, the nth parallel ray l p (β,n) corresponds to the nth fan beam ray l when the angle between the ray source and the X-axis is β2=β f (β,n), and 1≤n≤d, calculate the intersection set of parallel beams and the intersection set of fan beams, and use the intersection points between parallel beams and fan beams to describe this mapping relationship. The projection angle is β1=β i The nth i Parallel rays and projection angle β2 = β j The nth j The intersection of the parallel rays represents a p (β i ,n i ,β j ,n j), and 1≤i≤s, 1≤j≤s, 1≤n i ≤d,1≤n j ≤d; projection angle is β1=β i The nth i The fan-shaped beam rays and the projection angle are β2 = β j The nth j The intersection point of the fan-shaped beam rays is denoted by a f (β i ,n i ,β j ,n j ), and 1≤i≤s, 1≤j≤s, 1≤n i ≤d,1≤n j ≤d; F is used to represent the reconstruction value, then

[0014] F(a p (β i ,n i ,β j ,n j ))=F(a f (β i ,n i ,β j ,n j )) (3)

[0015] The reconstruction area is a square with a side length of D, which is divided into N×N square grids, denoted as F' N×N , the intersection point set of the fan beam located in the reconstruction area is screened as S f , the corresponding intersection point set of parallel beams is S p At this time, the fan beam is assumed to be parallel, so the fan beam intersection set S can be obtained by using the reconstruction formula of the parallel beam Radon transform f The reconstruction value F f :

[0016]

[0017] Where p represents the projection data matrix, (r,θ) represents the polar coordinates of the reconstruction point, (φ,l) represents the integration direction, and F f Represents the reconstruction value at the point (r, β) parallel beam.

[0018] Discretizing formula (4) yields:

[0019]

[0020] Where L is the length of the detection plate, t is the projection angle index, and 0≤t≤s-1, Δφ is the interval of the projection angle, and p is fitted three times. There are a total of dp.pieces fitting intervals. The range of the wth fitting interval is (l min,w ,l max,w ), a 2,w 、a 1,w 、a 0,w are the coefficients of the quadratic term, the linear term and the constant term of the fitting, respectively, 0,w represents an integral quantity, we have:

[0021]

[0022] in

[0023] R t =rcos(θ-φ t ) (7)

[0024] From formula (5), we can get F′ f After that, according to the mapping relationship S f →S p Get F′ f →F p , divide the reconstruction area into uniform grid points S N×N , put F p At the corresponding scatter point coordinates S p Interpolate to S N×N Got F.

[0025] Step 2: According to the required density distribution, a new projection data matrix p′ is obtained from formula (8) as negative feedback information.

[0026] p'=W·F (8)

[0027] Where W represents the sensitivity matrix and F represents the reconstructed density distribution.

[0028] Step 3: The negative feedback projection data matrix at the kth iteration is recorded as p′(k). The kth deviation p(k)=pp′(k) can be calculated based on the feedback data matrix. The kth output information q(k) of the PID controller is obtained using formula (9). Then, the reconstructed distribution is updated using the formula and the described hypothetical parallel method. The selected digital PID controller is:

[0029] q(k)=q(k-1)+(K p +K i +K d )·p(k)-(K p +2·K d )·p(k-1)+K d ·p(k-2) (9)

[0030] Among them, K p ,K i , K d are proportional, integral, and differential coefficients respectively, and a fuzzy controller is introduced to establish the second norm of the deviation p(k) and the derivative of the deviation The fuzzy inference rules between the second norm of and the PID controller parameters are used to adjust K p ,K i , K d Changes in three parameters; q(k), q(k-1) are the outputs of the k-th and (k-1)-th iteration calculations of the PID controller, respectively; p(k), p(k-1), p(k-2) are the inputs of the k-th, (k-1)-th and (k-2)-th iteration calculations of the PID controller, respectively; k = 1, 2, 3, ..., m, where m is the preset number of iterations; the initial value is set to q N×N (0)=0,p N×N (0)=0,p N×N (-1)=0,p′ N×N (1)=0,0 N×N represents an N×N all-zero matrix;

[0031] Step 4: Determine whether the number of iterations k is greater than the preset number of iterations m. If yes, execute step 5. If not, set k=k+1 and return to step 2.

[0032] Step 5: End the iteration and output the imaging results.

[0033] The beneficial effects of the present invention are as follows: based on Radon transform, a new fan beam to parallel beam mapping method is proposed, the closed-loop control principle is introduced into the X-ray tomography image reconstruction algorithm, the inversion error is used as the control object, and the reconstructed image algorithm is used as the controlled object, and the quality of image reconstruction is significantly improved in the case of low-angle projection. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 It is an implementation flow chart.

[0035] Figure 2 It is a specific implementation of the original model.

[0036] Figure 3 is the hypothetical parallel beam coverage case.

[0037] Figure 4 This is the actual fan beam coverage

[0038] Figure 5 It uses the hypothetical parallel algorithm to reconstruct the effect diagram.

[0039] Figure 6 It combines the closed-loop control principle with the hypothetical parallel algorithm to reconstruct the effect diagram. DETAILED DESCRIPTION

[0040] The present invention will be further described below in conjunction with examples.

[0041] This example uses 5 X-ray sources and 71 detectors are evenly set on each side of the sensor. The reconstruction flow chart is as follows Figure 1 As shown in the figure, it is a hypothetical parallel X-ray tomography image reconstruction method based on the closed-loop control principle. The density distribution of the original distribution model is shown in the attached figure. Figure 4 As shown, the density value of the white part is 1, the density value of the black part is 0, and it is composed of 49×49 pixels. The specific steps include:

[0042] Step 1: Use the imaginary parallel method to calculate the density value of the measured area under the current sensor setting. The 7 X-ray sources used in this example are located at the seven vertices of the regular heptagon, emitting fan-shaped beams that irradiate the 71 sensors on the detection board on the opposite sides. The entire sensor is arranged with a total of 497 rays. Select an initial X-ray source, and the projection angle is set to 0 at this time. The corresponding projection data on the detection board is:

[0043]

[0044] The projection angle of the next ray source along the clockwise direction is Similarly, the corresponding projection data on the detection board is p2. Similarly, the projection angle of the last ray source is The projection data on the detection board is p7, so the projection data at all angles is:

[0045]

[0046] The fan-shaped beams of rays with medium angle distribution in real space correspond to parallel beams with medium distance in imaginary space, as shown in the following figure. Figure 2 As shown, at the projection angle In the case of Corresponding to, at the projection angle When the fifth ray in the equidistant parallel beam The projection angle is The fifth ray and projection angle of the fan beam are The intersection point of the 45th ray of the fan beam is expressed as Similarly, the projection angle is When the fifth ray of the parallel beam and the projection angle are The intersection point of the 45th ray of the parallel beam is expressed as Let F represent the reconstruction value, then:

[0047]

[0048] The intersection mapping relationship of other rays is similar. This mapping relationship is described by the intersection set between parallel beams and the intersection set of fan beams. The intersection set S of fan beams at all projection angles within the reconstruction area is selected. f , the corresponding intersection point set of parallel beams is S p At this time, since the fan beam is assumed to be parallel, the reconstruction formula of the parallel beam is directly used to obtain the fan beam intersection set S by Radon transformation. f The reconstruction value F f :

[0049]

[0050] Where p represents the projection data matrix, (r,θ) represents the polar coordinates of the reconstruction point, (φ,l) represents the integration direction, and F f represents the reconstruction value at the parallel beam of point (r,θ). Discretizing formula (13) yields:

[0051]

[0052] Where L is the length of the detection plate, t is the projection angle index, and 0≤t≤s-1, Δφ is the interval of the projection angle, and p is fitted three times. There are a total of dp.pieces fitting intervals. The range of the wth fitting interval is (l min,w ,l max,w ), a 2,w 、a 1,w 、a 0,w are the coefficients of the quadratic term, the linear term and the constant term of the fitting, respectively, 0,w represents an integral quantity, we have:

[0053]

[0054] in

[0055] R t =rcos(θ-φ t ) (16)

[0056] From formula (14), we can get F′ f After that, according to the mapping relationship S f →S p Get F′ f →F p , divide the reconstruction area into uniform grid points S 49×49 , put F p At the corresponding scatter point coordinates S pInterpolate to S 49×49 Got F.

[0057] Step 2: According to the required density distribution, a new projection data matrix p′ is obtained from formula (18) as negative feedback information.

[0058] p'=W·F (17)

[0059] Where W represents the sensitivity matrix and F represents the reconstructed density distribution.

[0060] Step 3: The negative feedback projection data matrix at the kth iteration is recorded as p′(k). The kth deviation p can be calculated based on the feedback data matrix. 49×49 (k) = p 49×49 -p′ 49×49 (k), use formula (19) to get the k-th output information q(k) of the PID controller, and then use the formula and the described hypothetical parallel method to update the reconstructed distribution. The selected digital PID controller is:

[0061] q(k)=q(k-1)+(K p +K i +K d )p(k)-(K p +2K d )p(k-1)+K d p(k-2) (18)

[0062] Among them, K p ,K i , K d are proportional, integral, and differential coefficients respectively, and a fuzzy controller is introduced to establish the second norm of the deviation p(k) and the derivative of the deviation The fuzzy inference rules between the second norm of and the PID controller parameters are used to adjust K p ,K i , K d The change of three parameters, taking the second norm of the deviation p(k) as the input error e(k), and fixing K i =2, select a two-input two-output fuzzy controller to achieve K p , K d Adaptive adjustment of two parameters, the left boundary of the fuzzy control function uses the Z-type membership function, the right boundary uses the S-type membership function, and the middle uses the triangle membership function. The function expression is:

[0063]

[0064] In the formula, parameters a and c determine the zero position of the triangular function, and parameter b determines the peak position of the function.

[0065] q(k), q(k-1) are the outputs of the k-th and (k-1)-th iterations of the PID controller, respectively; p(k), p(k-1), p(k-2) are the inputs of the k-th, (k-1)-th and (k-2)-th iterations of the PID controller, respectively; k = 1, 2, 3, ..., m, where m is the preset number of iterations; the initial value is set to q 49×49 (0) = O 49×49 ,p 49×49 (0) = O 49×49 ,p 49×49 (-1) = O 49×49 , p′ 49×49 (1) = O 49×49 ,O 49×49 Represents a 49×49 all-zero matrix.

[0066] Step 4: Set the preset number of iterations to m=100, and determine whether the number of iterative calculations k is greater than the preset number of iterations m. If yes, execute step 5; if not, set k=k+1 and return to step 2.

[0067] Step 5: End the iteration and output the imaging results. Use the hypothetical parallel algorithm to invert the density distribution inside the field as shown in the attached figure. Figure 4 As shown in the figure, the inversion results obtained by combining the closed-loop control principle with the hypothetical parallel algorithm are shown in the attached figure. Figure 5 shown.

[0068] By comparing Figure 3 , Attachment Figure 4 and attached Figure 5 ,Combining the closed-loop control principle and the imaginary parallel algorithm can ,effectively improve the image reconstruction quality, and the reconstructed image is closer to the ,original distribution and has clearer edges.

[0069] The above description of the present invention and its embodiments is not limited thereto, and the accompanying drawings are only one embodiment of the present invention. Without departing from the purpose of the present invention, any structure or embodiment similar to the technical solution designed without creativity shall fall within the protection scope of the present invention.

Claims

1. A method for X-ray tomography image reconstruction based on fuzzy PID control, It is characterized in that For X-ray scanning imaging with few angles, a new mapping of fan-beam projection data to parallel beams is established. The density distribution of the measured area is reconstructed using the classic reconstruction method, Radon inverse transform. The projection data is calculated from the reconstructed density distribution and the deviation is compared with the original projection data. Fuzzy PID control is used to adaptively adjust the PID controller parameters according to the deviation and the deviation change. The deviation is continuously reduced through multiple closed loops, so that the density distribution of the measured area can be quickly and accurately reconstructed, thereby obtaining the density distribution of the measured area. The proposed X-ray tomography image reconstruction method based on fuzzy PID control is characterized by reducing data loss by constructing a direct mapping from fan beam to parallel beam, and directly reconstructing the density distribution of the measured area using the Radon inverse transform method, which specifically includes the following steps: In X-ray tomography, the sensor is a regular s-gon with s≥3, the X-ray source is located at the s vertices of the s-gon, the detection plate is located at the s sides of the s-gon, and d detectors are evenly distributed on each side; the rays of each X-ray source are fan-shaped beams, covering only the detection plates on the opposite sides, and the sensor arranges a total of Q=sd rays; the proposed reconstruction method first calculates the intersection set between the fan-shaped beam rays and the intersection set between the imaginary parallel beams, calculates the reconstruction value of the fan-shaped beam intersection set by Radon inverse transform, and then maps it to the parallel beam intersection set, and then calculates the density value of the reconstructed area, and finally introduces a fuzzy PID closed-loop model for error control, thereby obtaining the final reconstructed two-dimensional density distribution value; specifically, the following steps are included: Step 1: Select an initial X-ray source. The projection angle is set to 0. The corresponding projection data on the detection board is: The projection angle of the next ray source along the clockwise direction is Similarly, the corresponding projection data on the detection board is p 2 ; Similarly, the projection angle of the last ray source, that is, the sth ray source, is The projection data on the detection board is p s , so the projection data matrix of all angles is: At a certain projection angle, the fan-shaped ray beams with equal angle distribution are parallel in the imaginary space, corresponding to the imaginary parallel beams with equal interval distribution. When the projection angle is β 1 =β, the nth parallel ray l p (β,n) corresponds to the angle between the ray source and the X-axis is β 2 = β when the nth fan beam ray l f (β,n), and 1≤n≤d, calculate the intersection set of parallel beams and the intersection set of fan beams, and use the intersection points between parallel beams and fan beams to describe this mapping relationship; the projection angle is β 1 =β i The nth i Parallel rays and projection angle β 2 =β j The nth j The intersection of the parallel rays represents a p (β i ,n i ,β j ,n j ), and 1≤i≤s, 1≤j≤s, 1≤n i ≤d,1≤n j ≤d; projection angle is β 1 =β i The nth i The fan-shaped beam rays and the projection angle are β 2 =β j The nth j The intersection point of the fan-shaped beam rays is denoted by a f (β i ,n i ,β j ,n j ), and 1≤i≤s, 1≤j≤s, 1≤n i ≤d,1≤n j ≤d; F is used to represent the reconstruction value, then F(a p (b i ,n i ,b j ,n j ))=F(a f (b i ,n i ,b j ,n j )) (3) The reconstruction area is a square with a side length of D, which is divided into N×N square grids, denoted as F' N×N , the intersection point set of the fan beam located in the reconstruction area is screened as S f , the corresponding intersection point set of parallel beams is S p ; At this time, the fan beam is assumed to be parallel, so the fan beam intersection set S can be obtained by using the reconstruction formula of the parallel beam Radon transform f The reconstruction value F f : Where p represents the projection data matrix, (r,θ) represents the polar coordinates of the reconstruction point, (φ,l) represents the integration direction, and F f represents the reconstruction value at the parallel beam of point (r, θ); discretizing formula (4) yields: Where L is the length of the detection plate, t is the projection angle index, and 0≤t≤s-1, Δφ is the interval of the projection angle, and p is fitted three times. There are a total of dp.pieces fitting intervals. The range of the wth fitting interval is (l min,w ,l max,w ), a 2,w 、a 1,w 、a 0,w are the coefficients of the quadratic term, the linear term and the constant term of the fitting, respectively, 0,w represents an integral quantity, we have: in R t =rcos(θ-φ t ) (7) From formula (5), we can get F′ f After that, according to the mapping relationship S f →S p Get F′ f →F p , divide the reconstruction area into uniform grid points S N×N , put F p At the corresponding scatter point coordinates S p Interpolate to S N×N Get F; Step 2: According to the required density distribution, the new projection data matrix p′ is obtained from formula (8) as negative feedback information: p'=W·F (8) Where W represents the sensitivity matrix and F represents the reconstructed density distribution; Step 3: The negative feedback projection data matrix at the kth iteration is recorded as p′(k). The kth deviation p(k)=pp′(k) can be calculated based on the feedback data matrix. The kth output information q(k) of the PID controller is obtained using formula (9). Then, the reconstructed distribution is updated using the formula and the described hypothetical parallel method. The selected digital PID controller is: q(k)=q(k-1)+(K p +K i +K d )·p(k)-(K p +2·K d )·p(k-1)+K d ·p(k-2) (9) Among them, K p ,K i , K d are proportional, integral, and differential coefficients respectively, and a fuzzy controller is introduced to establish the second norm of the deviation p(k) and the derivative of the deviation The fuzzy inference rules between the second norm of and the PID controller parameters are used to adjust K p ,K i , K d Changes in three parameters; q(k), q(k-1) are the outputs of the k-th and (k-1)-th iteration calculations of the PID controller, respectively; p(k), p(k-1), p(k-2) are the inputs of the k-th, (k-1)-th and (k-2)-th iteration calculations of the PID controller, respectively; k = 1, 2, 3, ..., m, where m is the preset number of iterations; the initial value is set to q N×N (0)=0,p N×N (0)=0,p N×N (-1)=0,p′ N×N (1)=0,0 N×N represents an N×N all-zero matrix; Step 4: Determine whether the number of iterations k is greater than the preset number of iterations m. If yes, execute step 5. If not, set k=k+1 and return to step 2. Step 5: End the iteration and output the imaging results.

Images