A Modified Regularization Approach to MIT Image Reconstruction Based on PSO-SA

Background: Magnetic induction tomography (MIT) is a promising imaging modality of electrical tomography. The precision of image reconstruction algorithms in MIT places the key role in its application in industrial and biomedical fields. Image reconstruction in MIT is an ill-posed inverse problem. Regularization algorithms are effective to handle the ill-posed problem. However it is difficult to select its proper regularization parameter, larger or smaller regularization parameter will both effect on the quality of the reconstructed images. Methods: In this paper, firstly, the regularization principle of image reconstruction of MIT is analyzed. Then, we take the dimension of the Hessian matrix as the prior knowledge to obtain a novel penalty term, and add it to hybrid regularization algorithm, including Tikhonov and NOSER regularization, to propose a modified hybrid regularization algorithm. Secondly, selection for regularization parameters of image reconstruction algorithm is considered as a typical optimization problem, Particle swarm optimization with simulated annealing (PSO-SA) algorithm is used to solve the optimal solution. Results: Numberial experiments with six typical models and cerebral hemorrhage model were carried out. The correlation coefficient (CC), relative error (RE) and condition number of Hessian matrix are used as evaluation metrics to evaluate the effectiveness of the proposed method. The results obtained from the proposed method have lower REs and higher CCs compared to the other two methods. And the condition number of Hessian matrix in the proposed method is greatly reduced. Conclusions: The experimental results can prove the validity of the proposed method. The proposed method is capable of enhancing accuracy and noise immunity, which most closely coincides with the true conductivity distributions. Our method has better reconstructed quality compared with the other two methods, and it provides a theoretical reference for the development of application of the MIT technique in clinical decision application .

Many image reconstruction schemes have been proposed in recent decades, and they can be classified into non-iterative and iterative algorithms. Noniterative algorithms mainly include linear back project (LBP), Tikhonov, NOSER, Total variation (TV), and hybrid regularization algorithm. Li et al. [8] presented an improved back projection image reconstruction algorithm based on the magnetic field lines distribution, which performs well in reflecting location and shape information of the perturbation. Sun et al. [9] also made some improvement in linear back projection algorithm by decomposing the existing measuring data and reformulating the backprojecting way. Hao et al. [10] presented an improved Tikhonov algorithm to alleviate the ill-posed characteristics of the MIT inverse problem. The two weighted matrices were estimated using the solution of the Tikhonov regularization algorithm. Lei et al. [11] proposed an extended Tikhonov regularization by a combination estimation method and an extended stabilizing item. Hyaric et al. [12] extended the NOSER regularization algorithm in more realistic electrode models to reduce much of the work needed in the pre-calculations. Andrea et al. [13] proposed TV function to realize in vivo imaging of physiological data for EIT, which was solved the lagged diffusivity method and the primal dual-interior point method (PD-IPM) and preserves well discontinuities in reconstruction profiles. Gong et al. [14] adapted Total generalized variation (TGV) regularization to the finite element model (FEM) framework for EIT reconstruction, which promotes more realistic images compared with TV regularization. Song et al. [15] proposed an adaptive total variation (SATV) regularization algorithm to optimize the regularization term and the regularization factor according to the spatial feature. Chen et al. [16] proposed a hybrid algorithm, which firstly produced the preliminary image region using Tikhonov regularization algorithm, and then obtained the final reconstruction image using variation regularization algorithm. He et al. [17] proposed a hybrid regularization algorithm including Tikhonov and NOSER regularization, which improved the spatial resolution and prevented distinguishable sketch of disturbances from the background in 3D imaging. Liu et al. [18] proposed a hybrid regularization algorithm by combining the advantages of the Tikhonov and TV regularization, which enhances spatial resolution effectively. Among Iterative algorithms, Newton-Raphson (NR), the Landweber, conjugate gradient (CG) and simultaneous iterative reconstruction technique (SITR) algorithm received considerable attention. Han et al. [19] introduced weighting matrix and L1norm regularization in iterative NR algorithm. This can reduce the drawbacks of large estimation errors and improve reconstruction image algorithm stability. He et al. [20] proposed a modified Newton-Raphson method to realize the EMT image reconstruction for solid-gas two-phase flow measurement. Yan et al. [21] proposed the Landweber iterative algorithm with fuzzy thresholding, in which the threshold value of each iteration was generated by minimizing the measure of fuzziness of current reconstruction image. Liu et al. [22] proposed an improved iterative Landweber algorithm, which uses a Tikhonov regularization reconstruction image as the initial iterative value for the iterative. Zhang et al. [23] used a wavelet function to fuse reconstruction images by iterative Landweber and conjugate gradient least square (CGLS) algorithms. Wang et al. [24] used the conjugate gradient algorithm for EIT image reconstruction based on pulmonary prior information. Yan et al. [25] presented a method of image reconstruction using the CG algorithm with regularized pre-processing factors. Following analyzing a pre-iterative method for SIRT algorithm, a preconditioning for the projected SIRT is proposed by Hao et al. [26] to accelerate the image reconstruction speed and alleviate the ill-posed nature of the EMT inverse problem. Considering extreme changes in conductivity causing the instability in the reconstruction, a penalty term is added to iterative schemes based on optimization methods.
Therefore, in the actual image reconstruction processes, regularization parameter selection has gradually become a research focus. Ando et al. [27] developed the empirical likelihood information criterion to select regularization parameter of the penalized empirical likelihood estimator. Dario et al. [28] proposed an automatic parameter selection for Tikhonov regularization in ECT by analyzing the energy evolution of the current density data as a function of the regularization parameter.An extended L-curve method proposed by Xu et al. [29], which determined the regularization parameter associated with either global corner or the new corner. Liao et al. [30] exploited generalized cross-validation (GCV) method for an automatic regularization parameter selection scheme. Wen et al. [31] used GCV to select regularization parameter for TV regularization problems, and determined the optimal regularization parameter in each iteration of the first-order primal-dual method. Scherzer et al. [32] chosen the regularization parameter for Tikhonov regularization algorithm by using Morozov's discrepancy principle. Recently, more efforts have been done in multiple regularization parameter selection for improving the quality of reconstruction images. Zhang et al. [33] presented particle swarm optimization (PSO) to draw inspiration from analyzing the dynamic characteristics of particle positions. Gong et al. [34] proposed an adaptive parameter technique for strategy adaptation in differential evolution (DE). However, there is a problem with image reconstruction algorithms all mentioned: most of the current iterative algorithms only rely on experience to determine penalty term parameters and matrix construction, but these factors have a decisive effect on the image reconstruction quality. As a result, the final scheme effect cannot reach the best. Inspired by the dimension of the Hessian matrix, a novel penalty term is determined and added to the combined regularization of Tikhonov and NOSER regularization, constructing a modified hybrid regularization algorithm. We compare it with Tikhonov regularization and other regularization process as the benchmark.
The combined particle swarm optimization (PSO) and simulated annealing (SA) optimizers have been adopted for optimal parameters selection. The correlation coefficient (CC), relative error (RE) and condition number of Hessian matrix are used as evaluation metrics, verifying the effectiveness of the proposed method.

MIT basic principle
There are two major computational problems in MIT, the forward problem and the inverse problem. The forward problem is to obtain the measurements voltages from the conductivity distribution. The inverse problem is to convert the measurements into images through suitable image reconstruction algorithms. The forward problem of MIT is defined as computing the electromagnetic fields for a given geometry, distribution of dielectric properties, operating frequency, and source current, subjected to known boundary conditions. Assuming time-harmonic fields with angular frequency  , the governing Maxwell's equations [35] are expressed as Where B is the magnetic flux density, E is the electric field intensity  is conductivity, 0  is permeability, s J is the current source in excitation coil, and the displacement currents are negligible,  is the Hamilton operator.
Magnetic vector potential A is introduced to solve the forward problem of MIT, which satisfies the following equation: By introducing the Coulomb standard condition, namely Equation (4) can be solved by solving finite-element method as a solution for A . Besides, the induced voltage of the sensing coil can be presented as follows: Where V is the induced voltage, coil n represents the turns of the coil and l d represents the line integral, l denotes the line along the sensing coils.
For image reconstruction, it is to find a solution (  ) to the MIT conductivity problem by constructing a discrete linear problem of the form: Where  denotes the conductivity distribution, V denotes the vector of the voltage which is induced in the sensing coils by the conductivity, n m R S   represents the sensitivity matrix, and m being the number of excitation coils and n that of receiver coils.
The MIT inverse problem It is ill-conditioned and under-determined for limited measure data. For achieving good MIT image, researchers have paid more attention to develop image reconstruction algorithms for MIT.

Regularization reconstruction
The image reconstruction process of MIT system involves the inversion of conductivity distribution as measured by a certain algorithm performance, using the measured voltage value V . According to the least-squares method [36], Equation (7) is expressed as the least sum of squares error of the mathematical model of the MIT system. This can make the calculated voltage value  S and the measured voltage value V match as closely as possible, and the mathematical model: In order to minimize the objective function, we calculated its partial derivative to obtain the vector of the reconstruction conductivity, evaluating k iterate estimates k  , using the formulas given by: The m V denotes the vector of the measured voltage, the process is iterated until the maximum number of iterations is reached or certain criteria are met. To provide stability to the ill-posed least squares problem by damping or filtering the unwanted low singular values in an implicit fashion, a penalty term, also referred to as regularization term is usually added into Equation (8). The general form of regularization algorithm can be described as: the data fidelity, weighing how close between V and  S , ) ( R  is penalty term, which impose constraints to the solution. Two typical regularization algorithms are described in detail below.

Tikhonov regularization
Tikhonov regularization algorithm is one of the representative methods for solving the ill-posed MIT inverse problem. The solution can be obtained through minimizing the following optimization objective functional For the standard Tikhonov regularization algorithm, the penalty term is 2 I , I denotes regularization matrix, which is able to be calculated by the Equation (12).
Then the solution of standard Tikhonov regularization algorithm is given by: Tikhonov regularization algorithms makes the numerical solution stable, but it will cause excessive smoothness to the edge of the reconstruction images.

NOSER regularization
Newton-one-step error reconstruction (NOSER) is the development based on satisfactory results for minimize Newton algorithm execution time. It takes only one step of a Newton's method, starting from a constant conductivity distribution and approximates linearization. For NOSER regularization algorithm, the objective functional The NOSER regularization algorithm uses for the constraining matrix and  is a constant known as a regularization parameter. Given an initial conductivity, the NOSER regularization algorithm solves the inverse problem by ： This algorithm can locate the perturbed regions accurately and quickly according to the Equation (15), but it is sensitive to the signal noise of inverse problem.

Hybrid regularization algorithm
Based on the previous analysis, adding both Tikhonov and NOSER penalty term to Equation (8) can decrease the condition number of the Hessian matrix. However, the high singular values of the Hessian matrix reduce rapidly, which will remove the detailed information due to the fact that the high singular values play key roles in image reconstruction. To get better reconstruction results, we take the dimension of the Hessian matrix as the prior information to obtain a novel penalty term, and add it into hybrid regularization algorithm including Tikhonov and NOSER regularization, to propose a modified hybrid algorithm. The objective function can be defined as follows: In Equation (16), the first term is called the data fidelity term, the second term is the NOSER penalty term, applied to enhance the stability, the third term represents Tikhonov penalty, filtering the contribution of lower singular values to the solution of Equation (16).The fourth term represents novel penalty term, which can dampen the quick attenuation of high singular values effectively.
L can be calculated as Where n is the dimension of the Hessian matrix and all the elements in the n-dimensional matrix B are 1. The Newton-Raphson iteration algorithm is used for solving the minimization of the objective function Equation (16). The first-order gradient function of the objective function of the modified hybrid regularization algorithm can be obtained The second-order gradient function (Hessian matrix) of the objective function can be calculated as follows: The final iterative formula is shown as: The singular value decomposition was incorporated to decompose three regularization Hessian matrix Equation (19).

Hybrid regularization parameters selection
The multi-parameter regularization algorithm offers a certain degree of flexibility and enhances the image quality. In this work, the particle swarm optimization (PSO) and simulated annealing (SA) optimizers are considered to select better regularization parameters for improving the performance of image reconstruction.
Particle swarm optimization (PSO) PSO is a global optimization algorithm inspired by the behavior of bird flocks and fish schools searching for food [37]. Every potential solution is called particle in PSO algorithm. Each particle i has its own position i    . In order to prevent the i th particle missing the optimal solution, the velocity of each particle is limited between . The search and optimization process is carried out by Algorithm 1.
Simulated annealing (SA) The concept of simulated annealing (SA) is proposed for finding optimized problems solutions in 1983 [38]. The SA mainly simulates the modeling of molecular movement in materials during annealing. In the process of annealing, it will accept a solution which is worse than the current solution with a certain probability, and jump out of the local optima. The method simulates the physical process of annealing a solid material and then cooling it naturally after heating it to a high enough temperature. Algorithm 2 gives specific solving process of SA. calculate particle's fitness 6: if fitness is better than that of the best particle then 7: update best particle and save fitness 8: end if 9: end for 10: for all particles in swarm do 11: retrieve best particle from neighborhood 12： update position and velocity 13: update position and velocity 14: end for 15:until reaching termination condition 1 16:return solution of best particle in swarm  Regularization parameters selection based on PSO-SA Regularization parameters selection in the inverse problem can be conducted as a typical optimization problem. Particle swarm optimization with simulated annealing behavior (PSO-SA) algorithm combines the advantages of PSO (fast calculation and easy mechanism) and the advantages of SA (ability to jump away from local optimum solutions and converge to the global optimum solution). In PSO-SA, particle swarms are randomly created according to the PSO structure. Then, the next position of each particle is determined by using the Equation (24) of position movement of the PSO. At this stage the SA metropolis acceptance rule is introduced for the new position of the particle. The rule determines whether to accept the new position or recalculate another candidate position according to the fitness function difference between the new and old positions. This enables the solution to jump out of local optima, and decreases the vibration near the end of locating a solution. For MIT inverse problem, the fitness function is defined as Where S is the sensitivity matrix,  ,  and  are the

Stimulation setting and steps
To evaluate the performance of the proposed method (named as the modified hybrid regularization algorithm), six typical models are carried out in COMSOL and MATLAB environment. The finite element method and single coil excitation strategy are used to solve the forward problem of 8-coil MIT system. The material of the coils is set to copper. The diameter and cross-sectional area of the coils are set to 50 and 4e-6 2 m , respectively. The radius of the object space is 100 mm . The conductivity of the background and the inclusion in the object are set as 0.001 and 2 , respectively. And we compared the proposed method with Tikhonov regularization and hybrid regularization (Tikhonov and NOSER regularization) process as the benchmark. The combined particle swarm optimization (PSO) and simulated annealing (SA) optimizers is applied for optimal parameters selection. In figures and tables, TK, HR and MHR represent the Tikhonov, hybrid regularization and the proposed method, respectively.

Evalution metrics
Three metrics were used to quantitatively evaluate the performance, including correlation coefficient (CC), relative error (RE) and condition number. Condition number is used to evaluate the degree of Hessian matrix. The larger the condition number, the more illconditioned the Hessian matrix. It can be defined as:

Reconstruction results and analysis of typical models Numerical stability
The Hessian matrix affects the stability of the reconstruction algorithms. And singular values distribution reflects the ill-conditioned degree of Hessian matrix. The larger the ratio between the highest and lowest singular values is, the more illconditioned the Hessian matrix is. In this section, we conduct a singular values decomposition simulation about Hessian matrix of three algorithms, as shown in Fig. 1. In addition, Hessian matrix's condition number was utilized as evaluation metric for the stability of three image reconstruction algorithms, which is shown in Table 1. From Fig. 1, the singular values of the original Hessian matrix are basically close to 0. The singular values of Hessian matrix obtained by hybrid regularization algorithm is higher than those of the Tikhonov regularization algorithm. But its high singular values decay rapidly. The attenuation speed of the high singular values with the proposed method relatively slow. And the proposed method produces the higher singular values than the other two algorithms.
The results above demonstrate weak ill-posedness of Hessian matrix with the proposed method. The condition number of the Hessian matrix decreases from 1.4675e+24 to 4.7963e+13 in Table 1, verifying that the proposed method is a potential algorithm with good stability.

Precision
We conduct six typical models simulations to analyze the precision of the proposed method. The results are shown in Table 2. And the interface of the real targets in models (a) and (b) are marked by a black solid circle to study on edges characteristics. The first column of Table 2 shows the original simulation models, the red denotes conductive materials namely copper, and the blue refers to air. The second column of Table 2 shows the reconstruction images obtained by Tikhonov regularization algorithm. The third column of Table 2 shows the reconstruction images obtained by hybrid regularization algorithm. The fourth column of Table 2 shows the reconstruction images obtained by the proposed method. The initial value of the proposed method was obtained by using the linear back project (LBP) algorithm. From Table 2, for the reconstruction results of the Tikhonov regularization algorithm, the edges are oversmoothed, resulting in the serious loss of edges feature information. And the degree of edges information losing in the complex models (c)-(f) are more serious than those in simple models (a) and (b).
In particular, the information losing in the model (f) is the most serious, so that the target can't be found. For the reconstruction results of hybrid regularization algorithm, the surrounding region of the target has stronger artifacts. The true shape of the target is hard to determine in the model (f). Compared with the Tikhonov and hybrid regularization algorithm, the proposed method has a cleaner background and better imaging performance. The location and size of the target is clearly revealed in the background area, and the artifacts in the images can be effectively eliminated.
Although the model (f) is the most complex model among the selected models, its shape can be still be well preserved by the proposed method. It is obvious that the proposed method has the best performance in terms of the six typical models in comparison with the other two algorithms The CC and RE values are used to evaluate the precision of the proposed method, as shown in Table 3. From Table 3 The average RE value of the Tikhonov, hybrid regularization and the proposed method are 1.090, 1.0303 and 0.832, respectively. The results above indicate that the quality of the reconstruction images of simple models is better than that of complex models.

Tolerance with different noise levels
To validate the noise tolerance of the proposed method, three models (b), (c) and (f) are chosen as examples. Gaussian noises with a signal-to-noise ratio (SNR) of 55, 45 and 35dB are added to original voltage data, respectively. The reconstruction images with three algorithms are shown in Tables 4, 5 and 6. The reconstruction images obtained by the Tikhonov and hybrid regularization algorithms are quickly degraded with the increase of noise levels, especially in the interface region of the targets. However, with a different degree of interference, the proposed method can obviously distinguish the target from the background, which means that the proposed method will significantly enhance the quality the reconstruction results when different the SNR of noise is added. For model (b), the reconstruction results of three regularization algorithms maintain well with low noise levels added. However, with the increase of noise levels, the reconstruction images obtained by using the Tikhonov regularization algorithm has more artifacts, which increases the blurring the true features of target. The reconstruction image of hybrid regularization algorithm also has stronger artifacts, especially in the surrounding region of the target. In contrast to hybrid regularization algorithm, the proposed method still can clearly distinguish the interface between the target and the background. For model (c), as the noiselevels increase, the interface of the target obtained by the Tikhonov regularization algorithm has serious degradation. The target shape with hybrid regularization algorithm is deformed to some extent. Compared with the other two algorithms, the proposed method obtains better effects with relatively high noise levels added. For model (f), when relatively high noise levels added are high, the position and size of the target obtained by the Tikhonov regularization algorithm can't be revealed in the background area, and hybrid regularization algorithms can hardly determine the true shape of the target. However, the proposed method still reconstructs the true shape of target. These results show that simple model (b) is less sensitive to noise than complex models (c) and (f), and the noise tolerance performance of the proposed method is better than the other two algorithms.

Reconstruction results and analysis of typical models
To further verify the performance of the proposed method, a brain CT 3D model is used in this paper. A head model was designed as the imaging area based on the size of 99% of the national standard for Chinese adults. The specific size parameters of the head model are shown in Table 10. The model comprised six materials (tissue types) listed in Table 11. The triangular mesh is used for the forward problem calculations. Figure 2a shows distribution of mesh quality with head model. The conductivity distribution of section (z=0) is selected as research area, as shown in Figure 2b. Three different left cerebrum hemorrhagic length with 5cm, 3cm and 2cm are taken as the lesions to investigate, respectively.    Three cerebral hemorrhage length of 5cm, 3cm and 2cm are set in the simulation, as shown in Table 12.
The first column shows the original distributions. The reconstruction images of three algorithms are shown in the other three columns of Table 12, respectively. And the lesion area is marked with black solid box. The lesion edges obtained by the Tikhonov regularization algorithm are over-smooth, increasing the blurring the true features. When the length of cerebral hemorrhage is 2cm, a part of the lesion area is reconstruction outside the black solid box. The construction lesions obtained by hybrid regularization algorithm have more artifacts, which effects the determination of the edge of lesion area. The proposed method obtains the better results. The lesion and the normal tissue are clearly distinguished. The disturbance of normal tissue significantly reduces and the quality of reconstruction images is obviously improved. Compared with the left and right brain, it can absolutely show the lesion area of the left brain. Moreover, the quality of the model with cerebral hemorrhage length 5 cm is superior than that of the model with cerebral hemorrhage length 2 cm. Table 13 gives the CC values of three cerebral hemorrhage length models. From Table 13, the average CC value of the Tikhonov, hybrid regularization and the proposed method are 0.7694, 0.8332 and 0.8858, respectively. Compared with Tikhonov and hybrid regularization algorithms, the proposed method produces larger CC values. For the results of the proposed method, the CC values of the model with cerebral hemorrhage length 5 cm, 3cm and 2cm are 0.8985, 0.8826 and 0.8764, respectively, verifying that the quality of the model with cerebral hemorrhage length 5 cm is better.  To test the noise immunity of the proposed method, noise with a signal-to-noise ratio (SNR) of 55, 45 and 35dB is added to simulated voltages of cerebral hemorrhage length (5cm), respectively, as shown in Table 14. From Table 14, the quality of reconstruction images obtained by the Tikhonov and hybrid regularization algorithms degrades quickly with the increase of noise. Compared with the Tikhonov and hybrid regularization algorithms, the proposed method is less sensitive to the change in the noise level. As the low noise levels are added, the reconstruction results obtained by the Tikhonov regularization algorithm can't clearly show the difference between the lesion and the normal tissue. With relatively high noise levels added, the lesion shape of hybrid regularization algorithm is difficult to discern, while the proposed method still can distinguish the interface between the lesion and the background. Table 15 gives the values of the three algorithms with different noise levels for cerebral hemorrhage length 5cm. From Table 15, compared with the other algorithms, the CC values of the proposed method are the largest. As the SNR changes from 55dB to 35dB, the CC values of Tikhonov regularization algorithm and hybrid regularization algorithm decrease by 0.2181 and 0.1519, respectively, while the CC values of the proposed method only decrease by 0.1034. The results show that the proposed method can endure a relatively high level of noise.

Conclusion
The image quality in MIT greatly depends on factors such as the penalty term and the chosen regularization parameter. This study attempts to promote the performance of hybrid regularization algorithm including Tikhonov and NOSER regularization via adding a novel penalty term. The penalty term is determined based on the dimension of the Hessian matrix. Then, the combined particle swarm optimization (PSO) and simulated annealing (SA) optimizers have been adopted for optimal hybrid parameters selection. Numerical simulations with six typical models and the cerebral hemorrhage model are performed to prove the effectiveness of the proposed method and to compare it with Tikhonov and hybrid regularization algorithm. The CC, RE and Hessian matrix's condition number are utilized as evaluation metrics. The results demonstrate that the proposed method has obvious advantages in terms of stability, accuracy and noise immunity.