An image reconstruction method for accurately identifying intracranial lesion areas in stroke

Abstract

The present invention discloses an image reconstruction method for accurately identifying intracranial lesion areas of stroke, which is used for electrical impedance tomography to reconstruct the conductivity distribution of hemorrhagic stroke. The method regards the electrical impedance tomography problem as solving a linear ill-posed problem and determines the objective function; the process of the reconstruction method is: obtaining object field information, solving the relative boundary measurement value vector and sensitivity matrix; setting initialization parameters; using a weighted matrix to transform the data fidelity term; solving the objective function in an iterative manner; judging whether the iteration is finished; and imaging according to the obtained grayscale value. The present invention effectively overcomes the edge over-smoothing effect and suppresses staircase artifacts, improves the stability, applicability and spatial resolution of the algorithm to reconstruct the image, and expands the application scope of the total variation regularization algorithm in the field of medical detection.

Description

Technical Field

[0001] The present invention belongs to the technical field of electrical impedance tomography, and in particular relates to an image reconstruction method for accurately identifying intracranial lesion areas of cerebral stroke. Background Art

[0002] Currently, hemorrhagic stroke is considered a serious brain disease that exposes brain tissue to the irritation of blood, leading to disability or death in patients. In the diagnosis of hemorrhagic stroke, computed tomography (CT) and magnetic resonance imaging (MRI) are two commonly used imaging techniques with high spatial resolution. However, both technologies are unable to detect the early process of hemorrhagic stroke, hindering timely treatment. In addition, CT and MRI have low temporal resolution, long imaging time, bulky equipment, and are relatively expensive.

[0003] Electrical Impedance Tomography (EIT) is a technology that can reconstruct the conductivity distribution in the detection area and has attracted widespread attention. Physiological or pathological changes in biological tissues can cause changes in conductivity, which can be reflected by the reconstructed images obtained with EIT. At the same time, compared with CT and MRI, EIT technology is portable, safe, economical and non-invasive. In addition, a higher temporal resolution can be obtained. Due to its superiority, EIT technology has been used in medical imaging of lung function, breast cancer, liver fat and brain function. In particular, when hemorrhagic stroke occurs, the distribution of conductivity changes due to the higher conductivity of blood than that of normal brain tissue, making EIT possible even in the early monitoring of hemorrhagic stroke. However, recovering the conductivity distribution from the measurement data of EIT is usually an ill-posed inverse problem, which is not enough to guarantee the existence, uniqueness and stability of the solution, resulting in a relatively low spatial resolution of the reconstructed image.

[0004] For ill-posed problems, in order to obtain an effective and stable solution to the inverse problem, a regularization method is usually used to approximate the true solution. At present, typical regularization methods include the Tikhonov regularization (TR) method based on the L2 norm and the total variation (TV) regularization method based on the L1 norm. The Tikhonov regularization algorithm has been widely used due to its simplicity, but the image reconstruction quality depends on the selection of regularization parameters, such as M. Rezghi et al. published in 2009 in the Journal of Computational and Applied Mathematics, Vol. 231, pp. 914-924, entitled "A new variant of L-curve for Tikhonov regularization". The Tikhonov regularization algorithm has good performance for continuous distribution of the measured medium. However, when the measured medium is discontinuously distributed, excessive smoothness is imposed on the boundary, thereby reducing the resolution of the reconstructed image. In order to maintain sharp edges, the total variation regularization algorithm provides a solution. This method preserves the discontinuity of the boundary and allows reconstruction of sharp edges to produce clearer images. It has received widespread attention in the field of EIT, such as G González et al., published in Computers & Mathematics with Applications, Volume 74, Pages 564-576 in 2017, entitled "Isotropic and anisotropic total variation regularization in electrical impedance tomography". However, since the total variation method often produces piecewise constant solutions, this method will produce a "staircase effect" in smooth areas while reconstructing sharp edges, reducing the overall resolution of the reconstructed image, which seriously limits the application of the total variation method in EIT.

[0005] In order to reduce the step effect of the total variation regularization method and improve the spatial resolution of its reconstructed image, the present invention proposes an image reconstruction method for accurately identifying the intracranial lesion area of ​​stroke to solve the inverse problem of electrical impedance tomography. By comparing the image effects of absolute reconstruction and relative reconstruction in the measured area by this method, the stability and applicability of the algorithm proposed in the present invention are verified. Summary of the invention

[0006] The technical problem solved by the present invention is to provide an image reconstruction method with good stability and applicability for accurately identifying intracranial lesion areas of stroke, which is used for EIT reconstruction of the conductivity distribution of hemorrhagic stroke. In order to reduce the difficulty of solving the objective function, a weighted matrix is ​​proposed to transform the data fidelity term; in order to further remove artifacts and improve robustness, a soft threshold operator is introduced; in order to solve the minimized objective function, an effective iterative algorithm is adopted; in order to verify the resolution of the reconstructed image by the proposed method, numerical simulation is used to study the image effects of absolute reconstruction and relative reconstruction. Compared with the Tikhonov regularization method and the total variation regularization method, the image reconstruction method for accurately identifying intracranial lesion areas of stroke proposed in the present invention shows its advantages in effectively suppressing step artifacts and accurately identifying object field media in electrical impedance tomography reconstruction images.

[0007] The present invention adopts the following technical solution to solve the above technical problems: an image reconstruction method for accurately identifying intracranial lesion areas of stroke, which treats electrical impedance tomography as a linear ill-posed problem Ag=b meas , where A is the sensitivity matrix, b meas is the relative boundary measurement value vector, g is the required imaging grayscale value, and the established minimization objective function is:

[0008]

[0009] in

[0010]

[0011]

[0012] Where S(g) is the objective function; λ is the regularization parameter, which is a positive scalar used to maintain the balance between the fidelity term and the mixed penalty term; ω∈[0,1] is a weighting factor that controls the trade-off between the first-order TV and higher-order TV penalty terms; ||·|| p and ||·||1 are the Lp norm and L1 norm of the vector respectively; represents the gradient; l x and l y Respectively represent the finite element decomposition in the x-direction and the y-direction; and are first-order and second-order discrete differentials, respectively.

[0013] The process of the reconstruction algorithm is as follows:

[0014]

[0015] The beneficial effects of the present invention are as follows: when calculating the objective function, the image reconstruction method for accurately identifying the intracranial lesion area of ​​stroke described in the present invention introduces a weighting matrix to transform the objective function in order to ensure the sparsity of the algorithm and reduce the difficulty of solving the objective function; in order to further remove artifacts and improve robustness, a soft threshold operator is introduced; finally, in order to solve the minimized objective function, an effective iterative algorithm is adopted. The present invention improves the stability and applicability of the algorithm, expands the application scope of the total variation regularization algorithm, effectively overcomes the edge over-smoothing effect of the traditional Tikhonov regularization algorithm and the serious image step artifacts of the total variation regularization algorithm, and has a good effect in improving the spatial resolution of the image reconstructed by the electrical impedance tomography technology. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 A flowchart of an image reconstruction method for accurately identifying intracranial lesion areas of stroke according to the present invention;

[0017] Figure 2 The circular single-section measured field, the mode of the excitation current and the measurement voltage, and the electrode distribution of the electrical resistance tomography system of the present invention;

[0018] Figure 3 The finite element mesh of the present invention is (a) finite element mesh of the forward problem and (b) finite element mesh of the inverse problem;

[0019] Figure 4 It is a schematic diagram of the image results of absolute reconstruction when the present invention selects the real distribution of six models (models 1-6), a typical TR algorithm, a TV algorithm and an image reconstruction method for accurately identifying the intracranial lesion area of ​​stroke;

[0020] Figure 5 It is a schematic diagram of the relatively reconstructed image results of a typical TR algorithm, a TV algorithm and an image reconstruction method for accurately identifying intracranial lesion areas of stroke when selecting the real distribution of two groups of models (models 7-12) in the present invention;

[0021] Figure 6 It is a quantitative evaluation of the spatial resolution of the relative reconstruction model.

[0022] In the figure: 1- measured field, 2- electrode, 3- excitation current, 4- measurement voltage. DETAILED DESCRIPTION

[0023] An image reconstruction method for accurately identifying intracranial lesion areas of stroke according to the present invention is described in conjunction with the accompanying drawings and embodiments.

[0024] The present invention provides an image reconstruction method for accurately identifying intracranial lesion areas of stroke, which is used for EIT reconstruction of conductivity distribution of hemorrhagic stroke. In order to reduce the difficulty of solving the objective function, a weighted matrix is ​​proposed to transform the data fidelity term; at the same time, in order to further remove artifacts and improve robustness, a soft threshold operator is introduced; finally, an effective iterative algorithm is used to complete the final inverse problem solution.

[0025] like Figure 1 The figure shows a flow chart of an image reconstruction method for accurately identifying intracranial lesion areas of stroke according to the present invention.

[0026] like Figure 2 The figure shows a circular single-section measured field of a resistance tomography system of electrical tomography, the modes of excitation current 3 and measurement voltage 4, and the distribution of electrodes 2, wherein 16 electrodes 2 are evenly distributed on the outer wall of the measured field 1.

[0027] The finite element method is used to solve the EIT, such as Figure 3 As shown in Figure 2. For the EIT forward problem, triangular meshing is used, as shown in Figure 2. Figure 3 (a); the EIT inverse problem is divided into square grids, as shown in Figure 3 As shown in (b); in order to avoid the impact of inverse violations on imaging quality during image reconstruction, the number of subdivided grids for the forward problem should be greater than the number of subdivided grids for the inverse problem.

[0028] In order to better reflect the difference between the image reconstruction method for accurately identifying intracranial lesion areas of stroke in the present invention and the typical TR algorithm and TV regularization algorithm, the solution results of six absolute reconstruction models and two groups of relative reconstruction models under these three regularization algorithms are given respectively.

[0029] The embodiment includes the following specific steps: A method for accurately identifying the intracranial lesion area of ​​cerebral stroke by image reconstruction, which treats electrical tomography as a linear ill-posed problem Ag = b meas In the formula, A is the sensitivity matrix, b meas is the relative boundary measurement value vector, and g is the imaging grayscale value.

[0030] The inverse problem can be expressed using the objective function in the form of least squares optimization: Where f(g) is the objective function. In order to solve the uncertainty problem of the inverse problem, the regularization method provides a new method to improve the stability of the solution. Add a regularization term, also called a penalty term, to the objective function of the least squares. The general form of the regularization method can be described as: Here, λ is a regularization factor that controls the least squares term (often called the data fidelity term) and the regularization term R(g).

[0031] The TR algorithm can replace the regularization term R(g) with the regularization term It can be described as: However, the TR algorithm imposes excessive smoothness on the boundary when the measured medium is discontinuously distributed, thereby reducing the resolution of the reconstructed image. In order to maintain sharp edges, the TV regularization algorithm provides a solution. The form of the TV regularization algorithm can be described as: Although the reconstructed image maintains the discontinuity on the edge, a piecewise constant solution will be produced, so that this method will produce a "staircase effect" in smooth areas while reconstructing sharp edges, thereby reducing the overall resolution of the reconstructed image.

[0032] In order to suppress the step effect of the TV regularization method and improve the resolution of the reconstructed image, the present invention proposes an image reconstruction method for accurately identifying the intracranial lesion area of ​​stroke to solve the inverse problem of electrical impedance tomography based on the data fidelity term and the regularization term. The objective function to be minimized is:

[0033]

[0034] in

[0035]

[0036]

[0037] Where S(g) is the objective function; λ is the regularization parameter, which is a positive scalar used to maintain the balance between the fidelity term and the mixed penalty term; ω∈[0,1] is a weighting factor that controls the trade-off between the first-order TV and higher-order TV penalty terms; ||·|| p and ||·||1 are the Lp norm and L1 norm of the vector respectively; represents the gradient; l x and l y Respectively represent the finite element decomposition in the x-direction and the y-direction; and are the first and second order discrete differentials respectively. The reconstruction steps are as follows:

[0038] (1) Obtaining physical field information:

[0039] For the reconstruction model, the boundary measurement values ​​and sensitivity matrices required for reconstruction are obtained respectively. The boundary measurement value is to place the measured object in the electrical tomography measurement system, and 16 electrodes (such as Figure 2As shown in the figure, the boundary voltage under the cycle excitation cycle measurement is collected by adopting the mode of current excitation voltage measurement and no measurement of the excitation electrode, and a total of 208 measurement values ​​are obtained; the sensitivity matrix is ​​calculated based on the boundary measurement voltage of the empty field without inclusions, combined with the sensitivity theory, and the calculation formula is:

[0040]

[0041] Where A ij is the sensitivity coefficient of the jth electrode pair to the ith electrode pair, φ i ,φ j The excitation current is I for the i-th electrode pair and the j-th electrode pair respectively. i , I j The electric potential distribution in the time domain.

[0042] (2) Set initialization parameters:

[0043] Regularization parameter λ = 75, weighting factor ω = 0.3, norm p = 0.9, positive constant ε = 10 -6 , the iteration ends θ min =10 -4 , the maximum number of iterations k max =50, penalty parameters β1 = 100, β2 = 6, β3 = 6, soft threshold coefficient κ = 0.08, Lagrange multiplier and the initial solution g 0 =0 n×1 .

[0044] (3) Use the weighted matrix to transform the data fidelity term:

[0045] In order to reduce the difficulty of solving the objective function, a weighted matrix is ​​used to transform the data fidelity term, so the minimized objective function can be transformed into:

[0046]

[0047] in

[0048]

[0049] Where W k is the weight matrix calculated in the previous iteration and updated in each iteration, g k+1 and g k is the column vector of the (k+1)th solution and the column vector of the kth solution, ε is a small positive constant used to prevent the denominator from being zero, and is chosen to be 10 -6 .

[0050] (4) Iterate and solve the objective function:

[0051] By introducing multiple auxiliary variables, the complex constraint problem is decomposed into a series of simple constraint sub-problems. At the same time, each sub-problem has a closed-form solution so that it can be solved using existing numerical methods. To this end, three auxiliary variables f, q and r are introduced, which can be expressed as:

[0052]

[0053] Its augmented Lagrangian function can be expressed as:

[0054]

[0055] Where γ1, γ2 and γ3 represent augmented Lagrange multipliers, and β1, β2, β3>0 are penalty parameters. Therefore, the solution of auxiliary variables and gray values ​​can be described as:

[0056]

[0057] The f subproblem can be expressed as:

[0058]

[0059] In the formula, is the dot-wise product, "sgn" represents the function of signum;

[0060] The q subproblem can be expressed as:

[0061]

[0062] The r subproblem can be expressed as:

[0063]

[0064] The subproblem for the gray value g can be expressed as:

[0065]

[0066] The conversion yields:

[0067]

[0068] In the formula, the superscript T represents the transposition operator. and A T A is a circulant matrix.

[0069] To further reduce the staircase effect, a region of interest is introduced, defined as:

[0070]

[0071] In the formula, R refers to the region of interest, and κ is the threshold coefficient for distinguishing the background from the medium. Therefore, the gray value solution g can be converted to:

[0072] g k+1 =R.*g k+1

[0073] In the formula, . represents the corresponding multiplication between elements;

[0074] Finally, the calculation and update of Lagrange multipliers can be expressed as:

[0075]

[0076] (5) Iteration termination criteria:

[0077] Determine whether the iteration meets the iteration termination condition Or whether the maximum number of iterations k≤k is reached max If yes, the iteration is terminated and the next step is performed; if not, set k = k + 1 and jump back to step (3) to continue iterating and solving.

[0078] (6) Output results:

[0079] Imaging is performed based on the grayscale value g obtained by the final solution.

[0080] Figure 4 Six typical medium models are selected for the present invention, and the results of the absolute reconstruction image are shown as follows: the traditional TR algorithm, the TV regularization algorithm and the present method. It can be seen from the reconstructed image that when the image is reconstructed using the TR algorithm, excessive smoothness is observed on the boundary of the reconstructed image, making it difficult to accurately identify the size and shape of the image. Compared with the TR algorithm, the TV method avoids excessive smoothness in the TR algorithm during the reconstruction process and has a better reconstruction effect. However, a serious staircase effect is found in the background of the reconstructed image. Compared with the TR and TV methods, the present method eliminates excessive smoothing and staircase effects well during the reconstruction process. The reconstructed image has the highest quality, with a clear background and clear boundaries, and can accurately detect the size and position of the hemorrhagic tissue in the brain.

[0081] In clinical hemorrhagic stroke research, it is of great significance to reconstruct the relative changes of conductivity distribution as the bleeding process develops. The purpose of relative reconstruction is to obtain the changes in conductivity distribution between two time intervals (model 9 and model 12). The absolute reconstruction of model 7, model 8, model 10 and model 11 is also depicted to clearly understand the performance of the algorithm, such as Figure 5As shown. It is obvious that the TR method has difficulty in distinguishing the boundaries of the medium in absolute reconstruction, while blurring is observed on the boundaries of the reconstructed image in relative reconstruction. For the TV method, the boundaries can be well identified in both absolute and relative reconstruction, and a serious stair-step effect is observed. In addition, for both TR and TV methods, the reconstructed medium is larger than the actual change. In the absolute and relative reconstructed images, the method proposed in the present invention can effectively suppress blurring and eliminate the stair-step effect. In addition, the reconstructed medium size is almost the same as the actual size. At the same time, spatial resolution (SpatialResolution, SRE) is used to quantitatively evaluate the image reconstruction quality, and its expression can be expressed as: In the formula, Among them, P d and P0 represent the area of ​​the medium and the entire detection area, respectively. i is the area of ​​the ith pixel, [g d ] i =1(if g i greater than one quarter of the maximum amplitude of g), otherwise select [g d ] i =0.

[0082] from Figure 6 It can be seen that compared with the typical TR algorithm and TV regularization algorithm, an image reconstruction method for accurately identifying intracranial lesion areas of stroke has the smallest spatial resolution, can effectively suppress edge blur, eliminate the step effect, and accurately depict the internal medium distribution of the measured field, which significantly improves the accuracy of solving the inverse problem of electrical impedance tomography and the quality of image reconstruction.

[0083] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. An image reconstruction method for accurately identifying intracranial lesion areas of stroke, characterized by: Electrical impedance tomography is formulated as a linear ill-posed inverse problem Ag = b meas , where A is the sensitivity matrix, b meas is the relative boundary measurement value vector, g is the imaging gray value, and the minimization objective function is established as: in Where S(g) is the objective function; λ is the regularization parameter, which is a positive scalar used to maintain the balance between the fidelity term and the mixed penalty term; ω∈[0,1] is a weighting factor that controls the trade-off between first-order TV and higher-order TV penalties; ||·|| p and ||·||1 are the Lp norm and L1 norm of the vector respectively; represents the gradient; l x and l y Respectively represent the finite element decomposition in the x-direction and the y-direction; and They are first-order and second-order discrete differentials respectively. The specific reconstruction steps are as follows: (1) Obtaining physical field information: The boundary measurement values ​​and sensitivity matrix required for reconstruction are obtained for the reconstruction model. The boundary measurement value is to place the measured object in the electrical tomography measurement system, evenly distribute 16 electrodes outside the measured field, adopt the mode of current excitation voltage measurement and no measurement of the excitation electrode, collect the boundary voltage under cyclic excitation cyclic measurement, and obtain a total of 208 measurement values; the sensitivity matrix is ​​calculated based on the boundary measurement voltage of the empty field without medium, combined with the sensitivity theory, and the calculation formula is: In the formula, A ij is the sensitivity coefficient of the jth electrode pair to the ith electrode pair, φ i ,φ j The excitation current is I for the i-th electrode pair and the j-th electrode pair respectively. i , I j Time field potential distribution; (2) Set initialization parameters: Regularization parameter λ = 75, weighting factor ω = 0.3, norm p = 0.9, positive constant ε = 10 -6 , the iteration ends θ min =10 -4 , the maximum number of iterations k max =50, penalty parameters β1 = 100, β2 = 6, β3 = 6, soft threshold coefficient κ = 0.08, Lagrange multiplier and the initial solution g 0 =0 n×1 ; (3) Use the weighted matrix to transform the data fidelity term: In order to reduce the difficulty of solving the objective function, a weighted matrix is ​​used to transform the data fidelity term, and the minimized objective function is transformed into: in, Where W k is the weight matrix calculated in the previous iteration and updated in each iteration, g k+1 and g k is the column vector of the (k+1)th solution and the column vector of the kth solution, ε is a small positive constant used to prevent the denominator from being zero, and is chosen to be 10 -6 ; (4) Iterate and solve the objective function: By introducing multiple auxiliary variables, the complex constraint problem is decomposed into a series of simple constraint sub-problems. At the same time, each sub-problem has a closed-form solution so that it can be solved by existing numerical methods. For this purpose, three auxiliary variables f, q and r are introduced, which are specifically expressed as: Its augmented Lagrangian function is expressed as: Where γ1,γ2 and γ3 are augmented Lagrange multipliers, β1,β2,β3>0 are penalty parameters; The solution for the auxiliary variables and grayscale values ​​is described as: The f subproblem is expressed as: Where "°" is the point-by-point product, "sgn" represents the function of signum; for the q subproblem, it is expressed as: The r subproblem is expressed as: The subproblem for the gray value g is expressed as: Convert to: The superscript T represents the transposition operator. and A T A is a circulant matrix; In order to further reduce the staircase effect, a region of interest is introduced and defined as: Where R is the region of interest, κ is the threshold coefficient for distinguishing the background from the medium, and the gray value solution g is converted to: g k+1 =R.*g k +1 , where . represents the corresponding multiplication between elements; the calculation and update method of Lagrange multipliers is expressed as: (5) Iteration termination criteria: Determine whether the iteration meets the iteration termination condition θ=||g k+1 -g k ||2 / ||g k ||2≤θ min Or whether the maximum number of iterations k≤k is reached max If yes, the iteration is terminated and the next step is performed; if no, set k = k + 1 and jump back to step (3) to continue iterating and solving; (6) Output results: Imaging is performed according to the gray value g obtained by the final solution.

Images