LINEAR OVER RANGES ITERATIVE ALGORITHMS FOR IMAGE RECONSTRUCTION IN ELECTRICAL CAPACITANCE TOMOGRAPHY

. The paper concerns the non-linear algorithms for image reconstruction in electrical capacitance tomography for which Jacobi matrix computation time is very long. The paper presents the idea of an iterative linearization in nonlinear problems, which leads to a reduction in the number of steps calculating Jacobi matrix. The linear Landweber algorithm with sensitivity matrix updating and non-linear Levenberg-Marquardt algorithm with Jacobi matrix updating in selected steps only were presented.


Introduction
Electrical capacitance tomography is an imaging technique which allows to visualize spatial (three-dimensional or crosssectional) and temporal distribution of electric permittivity inside tomographic sensor. The image is reconstructed from measurements of mutual capacitance of sensor electrodes surrounding examined space. This method was first time proposed by Maurice Beck and Andrzej Plaskowski from University of Manchester for monitoring of industrial processes [6,14].
Image reconstruction in electrical capacitance tomography is a non-linear problem [10]. The measured capacitances are nonlinear function of the spatial distribution of permittivity. The inverse problem consists of determining the distribution of electric permittivity inside the probe using the measured capacitance. Lack of knowledge of analytical form of non-linear inverse transform forces calculation of permittivity distribution by solving an optimization problem, for example by minimizing the mean square norm.
The Newton-Raphson method could be applied for leastsquare minimization. Because this method requires the calculation of Hessian matrix which is time consuming, the Gauss-Newton method is preferable in which the approximation of Hessian matrix is used. The application of Gauss-Newton algorithm in the electrical capacitance tomography is presented in [5]. The weakness of the Gauss-Newton method is bad behavior when Jacobi matrix is ill-conditioned, what is the case of electrical tomography. Even small measurement errors cause a significant change of step direction of iterative algorithm searching a functional minimum. For regularized non-linear least-square optimization the Levenberg-Marquardt method can be used. However, the selection of the value of the regularization parameter used in Lavenberg-Marquard method is not a simple task. The value could be selected experimentally or set manually. Some automatic methods of regularization parameter value were proposed like generalized cross validation (GCV) or L-curve.
The application of the Lavenberg-Marquard algorithm in electrical tomography is described in [1]. The application of GCV method for automatic selection of the value of regularization parameter in the Levenberg-Marquardt algorithm did not bring satisfactory results [1]. In contrast to the GCV, L curve method gave good results of image reconstruction [3].
Practical application of non-linear optimization algorithm is limited by the calculation time, because these algorithms need to determine the Jacobian matrix (solve a forward problem) in each step. The forward problem consists of determining the electric field distribution in tomographic probe by numerical solution of the generalized Poisson equation, which is a complex computational problem. For this reason, linearization and iterative algorithms for solving the linear problem are preferred in electrical capacitance tomography.
An interesting idea presented in this paper is to use an iterative linearization in non-linear problem. The iterative linearization reduces the number of steps of the algorithm in which the Jacobian matrix is determined. Transformation of linear and nonlinear algorithm in a linear over ranges algorithm was presented in the papers [15,16]. The linear Landweber algorithm with sensitivity matrix updating and Levenberg-Marquardt algorithm with the calculation of Jacobi matrix only in selected steps were discussed.

Landweber algorithm with sensitivity matrix updating
In the linear approach, the image reconstruction in capacitance tomography can be defined as the minimization of a quadratic residue of the standard linear equation with respect to the vector ε , which can be written as: ε εε 2 2 1 arg min 2 =-Sc (1) where the vector ε describes the spatial distribution of electric permittivity inside the sensor, S is the sensitivity matrix describing the impact of changes in the distribution of permittivity on the value of the mutual capacitance of electrodes in the tomographic sensor, and c is the vector of measured capacitances. The Landweber iterative algorithm, which searches a minimum of quadratic norm functional is given by the formula: The Landweber algorithm is steepest descent gradient type method with very small value of step length a which guarantees the convergence of the algorithm [8].
Because the linear model of the tomographic system (sensitivity matrix) depends on the distribution of permittivity in the tomographic sensor, the approximation of the sensitivity matrix is used in the algorithm. Most often, because there is no a priori information about the examined object, the sensitivity matrix calculated for a uniform distribution of the permittivity is used in the linear model. For the same reason, uniform distribution of permittivity is taken as the first stage of an iterative algorithm. So, taken the form of sensitivity matrix, which is a good option for the first steps of the algorithm, it is becoming worse approximation model of the system with the next steps. Suppose that in subsequent iterations of the algorithm, getting better estimate of the electrical permittivity distribution is found. This information can be used to determine a better approximation model of linear system at given step. If after a certain number of steps, the permittivity distribution changes significantly, a new approximation of the sensitivity matrix can be determined.
In the simplest version, the sensitivity matrix can be updated after a given number of iterations. In this way, ranges linear algorithm can be obtained, which solves the non-linear problem over linear ranges [13,17] . The step of the algorithm can be written using the following equation: and the function ( ) gi describes the moments of sensitivity matrix updating.
In the first approach it can be assumed that through a number of steps the algorithm is moving in a small range, in which the function does not change rapidly and can be approximated by linear relationship. If the sensitivity matrix has to be updated every certain number of iterations, the function has the form of a stepped function given by the formula: where ceil is a rounding up to the nearest integer, mod is the remainder of the division. In the case of observing changes in the estimated vector, the function defining moment of renovation sensitivity matrix may have the form: d is a certain small value. The value of the parameter that determine when to upgrade the sensitivity matrix, such as the number of iterations or threshold by which must change estimate the distribution permittivity is not obvious. This issue was the subject of research. Some papers [9,17] show that the modified Landweber algorithm with sensitivity matrix updating achieves far better results than the linear version. Even the use of two, three modifications of the sensitivity matrix, with the total number of steps equal to several hundreds, brings a great improvement in the quality of the reconstructed image. The additional modifications of the above-described ranges linear iterative algorithm are possible [12]. The first modification is adaptive control of step length parameter ( ) i a . The standard approach is a relaxation strategy, involving reduction of step length during the calculation, which allows the algorithm to reach a stabilized solution.
The proposed modification consists of checking whether the value of the residue decreases in the next step. If in the next step the condition of decreasing the value of the residuum is not reached, the algorithm returns to the previous estimate of the solution, and the value of the parameter ( ) i a is decreased. This strategy allows to take greater step length at the beginning of calculations to accelerate the search and prevents divergence of the algorithm in the case of too long step. Another modification of the Landweber algorithm is to use the normalized sensitivity matrix [12]. The method of normalization of the sensitivity matrix was taken from a linear back projection (LBP) algorithm [7]. The matrix is normalized in columns, so that the sum of the coefficients in the columns is equal to unity. The matrix normalized in columns is given by the formula: is column unit vector, and diag is an operator that creates diagonal matrix from a vector. If the measurements are normalized according to the following formula: where c , max c , min c are vectors of capacitance measurements for the object, for the uniform distribution of electric permittivity with the minimum min e and maximum max e value of the measurement range respectively, the reconstructed permittivity values are normalized and bounded in the range 0;1 . The iterative algorithm with the application of normalized sensitivity matrix which searches for a minimum of a quadratic norm is given by the formula: If the initial solution is the uniform permittivity distribution with , the first step of the algorithm is given by the equation (9) The first step of this algorithm is equivalent to a linear back projection algorithm LBP when step length ( ) The application of normalized sensitivity matrix compensates for the lack of sensitivity in the center of the field of view of the tomographic sensor. In our forward problem solver the sensitivity matrix is calculated inaccurately at the edges of the field of view because of the used Cartesian grid. It turned out that sensitivity matrix normalization minimized the impact of the errors of the sensitivity matrix coefficients for the edges of the field of view.
This approach is preferred in the case of image reconstruction with the sensitivity matrix calculated using a regular Cartesian discretization grid because it is less sensitive to the errors of calculation of the sensitivity coefficients on the edges of the sensor.

Levenberg-Marquardt algorithm with the calculation of Jacobi matrix only in selected steps
The problem of image reconstruction in electric capacitance tomography described by a nonlinear equation: can be defined as a search for roots of the function: In the electric capacitance tomography, it can be assumed that the problem is underdetermined and there are many solutions. The task is to find such N dimensional vector ε 12 , ,..., T N e e e éů = ęú ëű , for which the condition: is satisfied. If the function f is continuous and differentiable in the vicinity of the point e , it can be expanded in a Taylor series: where ( ) ε J is a partial derivative matrix and ( ) O ε represents the components of the higher order. Neglecting higher order components, the following relationship is obtained: where e D is the increase defined by the tangent to the function at a point ε or in other words the distance from this point to the intersection of the tangent with the abscissa. From the above equation and using linear approximation of function h , the following iterative formula can be obtained: The matrix 1 -

J
does not exist in most cases in electric capacitance tomography because the problem is undetermined. Instead of the direct inverse the left pseudoinverse can be used what leads to the formula known as the Gauss-Newton method: If we introduce regularization parameter 2 m I , where m is a positive number, and I is the identity matrix, we get regularized Gauss-Newton algorithm called Levenberg-Marquardt (LM) algorithm [11]. Additionally, if we take into account the step length ( ) i a , we obtain the following iterative formula: If the function is highly non-linear, Jacobi matrix (sensitivity matrix) should be calculated in each step of the algorithm. However, if a function within a certain range is approximately linear, Jacobi matrix does not change much. If we assume that the function domain could be divided into several ranges where the function is approximately linear, it is possible to modify the algorithm, so that the Jacobi matrix will be constant in each range and it will change only when moving to an adjacent range [18]. Then the algorithm takes the form: gi is the function describing moments of updating of sensitivity matrix. The problem is to determine linear ranges in this method. There are different approaches, including simple, assuming that the linear approximation is valid for a certain number of steps, during which the Jacobi matrix is not updated. In this approach, the function has the form of a step function given by the formula (4), as in the modified algorithm Landweber with sensitivity matrix updating. A more sophisticated method for determining the moment of Jacobi matrix calculation can consist in observation of changes of the estimated vector ε . If the change of the estimated vector from the last modification of Jacobi matrix exceeds a certain preset threshold, the Jacobi matrix should be recalculated. In this method, the function describing the update moments of the sensitivity matrix may be in the form given by (5).

Results of image reconstruction
The described above linear over ranges algorithms were verified experimentally using numerical simulated data and real measurements. Here, exemplary results of the experiments will be shown only for modified Levenberg-Marquardt algorithm. The sensitivity matrix updating scheme used in both algorithms is similar. The Levenberg-Marquardt algorithm switches between the Gauss-Newton algorithm (when the regularization parameter is very small) and steepest descent algorithm (where the regularization parameter is large).
Measurements were performed using ET3 tomograph build at our Division [2]. The multichannel device worked with different gains for pairs of adjacent (low amplification) and opposite electrodes (high amplification). The cylindrical sensor with one ring of 16 electrodes was used. The diameter of the cylinder was equal to 160 mm. The height of the electrodes was equal to the cylinder diameter (the dimension in the direction of the axis of the cylinder). The physical object used in the experiments consisted of 6 Perspex rods with relative permittivity value equal about 2.9. The diameter of the rods was equal 20 mm. The length of rods was much greater than electrode height. The shape of the test object was selected arbitrary to allow two-dimensional modeling of the electric field distribution in a tomographic sensor. The object placed in the sensor is shown in Fig. 1. The collected data for the object were normalized using the measurements for empty and fully filled with Perspex sensor.
Two-dimensional simulations of measurements were performed using Matlab toolbox -ECTsim‖ developed by our group. The numerical modeling of the sensor and the object, electric field calculation and sensitivity matrix calculation is performed using regular Cartesian grid. The size of discretization matrix of sensor model used in that simulation was 5252 elements. The shape of the object does not change in the Z-axis so that the two-dimensional approximation of the electric field distribution is close to the three-dimensional distribution. The numerical representation of the object used in the real measurements is shown in Fig. 2. The capacitance measurements were simulated for the object. To simulate a measurement process in electrical capacitance tomography so called forward problem has to be solved. The potential distribution in the sensor with boundary conditions imposed by voltage applications on the electrodes has to be calculated. Knowing the voltage distribution the mutual capacitance formed by two electrodes selected from the set of electrodes surrounding the object can be calculated using the Gauss's law. An alternative method of capacitance calculation is application of discrete linear approximation = cS ε .
(19) This method was used in ECTsim toolbox. The potential distribution was calculated for all application electrodes for the specified sensor geometry and the permittivity distribution inside the sensor. The mutual impedance principle is used for sensitivity matrix calculation. The sensitivity distributions for all electrode pairs are calculated. The calculated two-dimensional sensitivity maps constitute the rows of the sensitivity matrix. Because the measured value of capacitance formed by two electrodes is the same independently of the fact which electrode is an application or measuring electrode, there are The Gaussian noise was added to the generated data. The measurement error (standard deviation) corresponds to around 0.1 % of the measurement range limited by the maximum value measured for the sensor full filled with the material of high permittivity.
The modified LM algorithm was implemented using the equation (18). The regularization parameter was calculated automatically by the L-curve method [4]. The value of the regularization parameter is a compromise between the value of residual error and the norm of the regularized solution. The norm of a regularized solution versus the norm of the residual norm for different value of regularization parameter is plotted in log-log scale. The value of the regularization parameter is selected as a point that corresponds to the corner of L-curve. The curves for reconstruction from simulated data and real measurement are shown in Fig. 3 and Fig. 4 respectively. The value of the regularization parameter is calculated after every update of the sensitivity matrix. The step length ( ) i a was constant for all iterations. The value was set arbitrarily and equal 0.15 for reconstruction from simulated data and 0.1 for reconstruction from real measurements.
The sensitivity matrix (Jacobian) was calculated only a few times in the iterative process. In the shown example the sensitivity matrix was updated every 100 steps. The moments of sensitivity matrix recalculation are visible on the plot of reconstruction error. The algorithm was stopped arbitrarily after 650 steps when the algorithm reached a steady state. The reconstructed distribution of permittivity for simulated and real data are shown in Fig. 5 and Fig. 6.
To evaluate the reconstruction quality quantitatively a square norm was used. The norm was calculated in image space, capacitance measurements space and sensitivity matrix space. In image space, the norm is a discrepancy between the reconstructed IAPGOŚ 1/2017 p-ISSN 2083-0157, e-ISSN 2391-6761 and true permittivity distribution. In case of real measurement data, the numerical representation of the object is used as true permittivity. In capacitance space, the capacitance residual error is the normalized discrepancy between the measured capacitances and the values calculated for the given estimate of permittivity. In sensitivity matrix space the error is the norm between the sensitivity matrix computed for a true permittivity distribution and the sensitivity matrix computed for given estimate of a permittivity distribution. In case of real data the sensitivity matrix for a true permittivity distribution is calculated using the numerical model of a physical object. The presented in this paper nonlinear iterative algorithm updates the sensitivity matrix. The norm in sensitivity space shows how the approximated model fits to the true model. The square norm is normalized using the formula: For example in image space t x denotes the true permittivity distribution. 0 x is the start solution and i x is the permittivity estimate at i-th step of iterative algorithm.
The plots of reconstruction error for simulated and real data are shown in Fig. 7 and Fig. 8. The points of discontinuity in sensitivity matrix error corresponds to the moments of sensitivity matrix recalculation. The sensitivity 2D map calculated for uniform distribution and obtained after 6 updates are shown in Fig. 9 and Fig. 10 respectively.

Conclusions
The linear over ranges algorithms were introduced in this paper. The results of image reconstruction using modified Levenberg-Marquardt algorithm were presented. The linear over ranges LM algorithm was convergent. The reconstruction quality norms decrease in function of iteration number with one exception. The discrepancy norm for a sensitivity matrix increases after first update of sensitivity matrix but decreases after next updates. This behavior was observed for both cases: reconstruction from simulated data and reconstruction from real measurements.
Only a few updates of the sensitivity matrix were enough to obtain acceptable image quality. The limited number of updates of Jacobi matrix up to six for dozens of iterations of the modified LM algorithm did not degrade image quality, and significantly reduces the number of calculations. Moreover, in the experiments the simplest method of selection of the moments when the sensitivity matrix was recalculated was used.
The carried out experiments had shown that the linear over ranges algorithm are very effective in terms of picture quality relative to the number of calculations. The application of the custom finite element method previously elaborated for sensitivity matrix calculation [19] accelerates the computations but the forward problem computations are still very time consuming comparing to the computation of the inverse transform.