Noninvasive Epicardial and Endocardial Extracellular Potentials Imaging with Low-Rank and Non-local Total Variation Regularization

Abstract

Epicardial and endocardial extracellular potentials imaging has significant implications for diagnosing cardiac diseases. In this paper, we propose a novel noninvasive epicardial and endocardial potentials (EEP) reconstruction method based on low rank and non-local total variation regularization (LRNLTV). The low rank constraint retains the spatiotemporal correlation of EEP and the nonlocal TV constraint make use of the non-local similarity. The constrained optimization problem is solved by augmented Lagrangian multiplier (ALM) method. The simulated PVC data, intervention premature data and clinical PVC data while both EEP map and activation map together with the correlation coefficient (CC), structural similarity (SSIM) and locating error are used to evaluate the proposed method. We also compare LRNLTV method with TV method (L1-norm based) and low rank method, and the proposed framework achieves better results both in reconstruction accuracy and the boundary.

1 Introduction

Reconstructing epicardial and endocardial extracellular potentials from the records of the cardiac electric events by body surface ECG has significant values for better understanding and diagnosis of cardiac disease. However, even with many efforts from noninvasive electrophysiological imaging community, it remains one of open key issues comparing with locating the arrhythmia site accurately by place the catheter on the surface of the heart [1], while is invasive and poses a certain threat to human health.

Improving the reconstruction quality of epicardial and endocardial potentials (EEP) imaging is essential, especially in applications like location of ventricular tachycardia (VT) [2] and premature ventricular contraction (PVC). With the wide adoption of regularizations to achieve unique solution, various approaches have been proposed, including the use of L2-norm based method, L1-norm based method and the use of method based on the rank of the matrix and L0-norm.

Generalized minimum residuals (GMRes) which is a L2-based method has been employed in [3] to reconstruct the epicardial potential during activation and repolarization of ventricular pacing. However, the shortcoming of L2-norm that the reconstruction results are too smooth and the accuracy is low have been pointed out in [4], and the total variation regularization (L1-based method) was proposed to improve the reconstruction accuracy. However, TV method only take use of the local similarity of the heart surface potential and compute frame by frame, the computing time is too long, the spatiotemporal correlation get lost as well, and a low rank and sparse decomposition method based on sparsity theory was presented in [5], the spatiotemporal correlation is aggregated in the method to improve the reconstruction accuracy, but the noise in ECG still have a great impact on the result.

In this study, we introduce a novel low rank and non-local total variation regularization (LRNLTV) method to solve the ill-posed inverse problem of EEP imaging. Base on the sparsity of EEP, it can be decomposed into two parts: the low rank background part and the sparse time-variant part, and the non-local total variation regularization is applied to the low rank background part to retain the boundary information. The innovations of the method is as follows: (1) spatiotemporal correlation of the EEP can be extract by low rank and sparse decomposition to get the dynamic potential of each node on heart surface; (2) non-local TV constraint take advantage of the nonlocal similarity to obtain better performance on the boundary.

2 Methodology

In this section, a low rank with non-local total variation (LRNLTV) framework is present for solving the inverse problem of EEP imaging. At first, we use the boundary element model (BEM) to obtain the relationship between the EEP and body surface potential (BSP) in the forward model of EEP imaging problem [6]. Then, we introduce the LRNLTV framework to reconstruct the EEP from BSP. At last, the optimization problem is solved by the augmented lagrangian multiplier method, singular value thresholding [7] (SVT) and the gradient descent algorithm.

2.1 BEM-Based Dynamic Electrocardiography Model

In this part, we introduce the dynamic forward model for the relationship between BSP and EEP. Different from the traditional standard “capped” full heart epicardial model, our heart surface model is a union of the endocardium and epicardium, where the heart surface denotes both the endocardium and epicardium surface. According to the previous work, we assume that the potential of cardiac electrical source distributed on a closed surface, and the potential distribution can be considered as an equivalent source. Thus, the relationship between cardiac electrical source and body surface can be expressed by,

$$ \upsigma_{e} \nabla^{2} \phi_{e} \left( r \right) = \nabla \cdot \left( { - D_{i} \left( r \right)\nabla \varphi \left( r \right)} \right)\quad r \,in\,\Omega _{h} $$

(1)

In which, $\upsigma_{e}$ denotes the body conductivity in cardiothoracic cavity, $\phi_{e} \left( r \right)$ represent the value of extracellular potential at spatial coordinate r, $D_{i}$ is the intracellular conductivity tensor of myocardial cells, $\varphi \left( r \right)$ is transmembrane potential of cardiac myocytes and $\Omega _{h}$ denotes the interior space of myocardial tissue. In order to simplify the forward model, we suppose that there is no other power electrical source in the volume between the heart surface and the body surface and it is an isotropic homogeneous medium, then the cardiac electric field can be represented as a Laplace equation, unify the boundary condition, it can be solved by BEM as,

$$ \upphi = {\text{Hu}} $$

(2)

where $\upphi \in {\text{R}}^{{{\text{J}} \times 1}}$ denotes the potentials measured on the body surface with J-lead electrodes at a particular moment, u ∈ R^I×1 is the EEP to be reconstructed and H ∈ R^I×J is a time invariant transfer matrix for specific patients which can be get from the heart-torso model built with the thoracic CT scan or MRI scan data.

However, electrocardiogram signal is a time sequence signal include temporal and spatial information, and the above ϕ is only one frame of the ECG data with spatial information, the single frame forward model can be extended to a multi-frame one,

$$ \Phi = {\text{HU}} $$

(3)

$$ {\text{U}} \in {\text{R}}^{{{\text{I}} \times {\text{C}}}} , \,\upphi \in {\text{R}}^{{{\text{J}} \times {\text{C}}}} $$

where $\Phi = \left[ {\begin{array}{*{20}c} {\upphi_{1} } & {\upphi_{2} } & {\begin{array}{*{20}c} \ldots & {\upphi_{c} } \\ \end{array} } \\ \end{array} } \right]$, $\upphi_{c}$ is the c-th frame of body surface and ${\text{U}} = \left[ {\begin{array}{*{20}c} {{\text{u}}_{1} } & {{\text{u}}_{2} } & {\begin{array}{*{20}c} \ldots & {{\text{u}}_{c} } \\ \end{array} } \\ \end{array} } \right]$, ${\text{U}}_{c}$ is the c-th frame of EEP.

2.2 Low Rank and Sparse Decomposition

The epicardial and endocardial potential is inherent sparse. Therefore, EEP can be decomposed into two parts [8], while the first one is the background part that almost never change with time. And the second part is the time-variant part that changes rapidly over time. The background part represents the highly relevant information between frames and the time-variant part refers to the potential or activation propagation. Thus, we apply the low rank and sparse decomposition approach to reconstruct the EEP, and the EEP matrix can be written as,

$$ {\text{U}} = {\text{L}} + {\text{S}} $$

(4)

where L is the background part and S is the time-variant part. When measuring the body surface potential information, noise will inevitably be introduced while small noise might cause large error in U. However, if the rank of background part L is as small as possible and the time-variant part is as sparse as possible, the effect of the noise will be effectively eliminated, and the inverse problem becomes a constrained minimization problem,

$$ \hbox{min} rank\left( {\text{L}} \right) + \lambda \left\| {\text{S}} \right\|_{0} + \mu \left\| {\Phi - {\text{HU}}} \right\|_{F} \quad s.t. \, {\text{U}} = {\text{L}} + {\text{S}} $$

(5)

In which $rank\left( {\text{L}} \right)$ is the rank of L, $\left\| {{\text{S}}_{0} } \right\|$ indicates the number of nonzero elements in S, $\mu$ is the weight parameter of the fidelity item, and $\left\| {\Phi - {\text{HU}}} \right\|_{F}$ is Frobenius norm of the matrix, which is defined quadratic sum of the elements in the matrix. However, the minimization problem in Eq. (5) is a nonconvex problem, we need to relax the above equation to a convex form as nuclear norm is a convex approximation of rank and L1 norm is a convex relaxation of L0 norm [9], and the objective function can be expressed by,

$$ { \hbox{min} }_{{{\text{L}},{\text{S}}}} \left\| {\text{L}} \right\|_{*} + \lambda \left\| {\text{S}} \right\|_{1} + \mu \left\| {\Phi - {\text{HU}}} \right\|_{F} \quad s.t. \,{\text{U}} = {\text{L}} + {\text{S}} $$

(6)

where $\left\| {\text{L}} \right\|_{*}$ denotes the nuclear norm of L, which is the sum of singular values of L, and $\left\| {\text{S}} \right\|_{1}$ represents the L1 norm of S, which is defined as the sum of all the element in matrix S. The background part L and the time-variant part S can be obtained by solving the Eq. (6), and then we can get the solution U which we are interested. The detailed solving process will be introduced in Sect. 2.4.

2.3 Nonlocal TV Constraint

Although the presented low rank and sparse decomposition can eliminate the influence of noise to a certain extent, the remaining noise in the background part L will cause the unsmooth solution, therefore, non-local total variation constraint is adopted to background part L to ensure the non-local smoothness of the solution and better details. According to the work in [10], we define the nonlocal gradient as,

$$ \nabla_{NLTV} {\text{L}}\left( {a,b} \right) = \left( {{\text{L}}\left( b \right) - {\text{L}}\left( a \right)} \right)\sqrt {w\left( {a,b} \right)} $$

(7)

where ${\text{L}}\left( a \right)$ and ${\text{L}}\left( b \right)$ denotes the value of position a and b in the background matrix L, and $w\left( {a,b} \right)$ indicates the weights between voxel a and voxel b that is non-negative and symmetric. $w\left( {a,b} \right)$ can be formulated as,

$$ w\left( {a,b} \right) = { \exp }\left( { - \frac{{\mathop \sum \nolimits_{k = - i}^{i} G\left( k \right) \cdot \left| {{\text{L}}\left( {a + i} \right) - {\text{L}}\left( {b + i} \right)} \right|^{2} }}{{2h^{2} }}} \right) $$

(8)

where $G\left( k \right)$ is Gaussian kernel with patch size ((2i + 1) × (2i + 1) and $h$ denotes the filtering parameter. Then, the NLTV constraint term for background matrix L can be given as,

$$ J\left( {\text{L}} \right) = \left\| {\nabla_{NLTV} {\text{L}}} \right\|_{1} = \sum\nolimits_{a} {\sqrt {\sum\nolimits_{{b \in\Omega }} {\left( {{\text{L}}\left( b \right) - {\text{L}}\left( a \right)} \right)^{2} w\left( {a,b} \right)} } } $$

(9)

Thus, the low rank and non-local TV framework for the EEP imaging problem is present below,

$$ { \hbox{min} }_{{{\text{L}},{\text{S,}}\,{\text{U}}}} \left\| {\text{L}} \right\|_{ * } + \lambda \left\| {\text{S}} \right\|_{1} + \mu \left\| {\Phi - {\text{HU}}} \right\|_{F} + \alpha J\left( {\text{E}} \right) \quad s.t. \, {\text{U}} = {\text{L}} + {\text{S}}, {\text{L}} = {\text{E}} $$

(10)

2.4 Optimization Method

We apply the Augmented Lagrange Multiplier method (ALM) [11] to the above optimization problems so that the problem becomes an unconstrained problem, and the augmented Lagrange function is present below,

$$ \begin{aligned} {\mathcal{L}}\left( {{\text{L}},{\text{S}},{\text{U}},{\text{E}}} \right) & = \left\| L \right\|_{*} + \lambda \left\| {\text{S}} \right\|_{1} - \left\langle {Z,{\text{U}} - \left( {{\text{L}} + {\text{S}}} \right)} \right\rangle + \frac{\beta }{2}\left\| {{\text{U}} - \left( {{\text{L}} + {\text{S}}} \right)} \right\|_{F}^{2} + \\ & \frac{\mu }{2}\left\| {{\text{HU}} -\Phi } \right\|_{F}^{2} + \alpha J\left( E \right) + \left\langle {Z_{L} ,{\text{L}} - {\text{E}}} \right\rangle + \frac{{\beta_{L} }}{2}\left\| {{\text{L}} - {\text{E}}} \right\|_{F} \\ \end{aligned} $$

(11)

where $\beta$ and $\beta_{L}$ are penalty parameters, $Z$ and $Z_{L}$ denotes the Lagrange multipliers. However, it is difficult to solve the above function directly, thus we divide it into three sub-problems and solve the sub problems separately.

L, S Subproblem

Separate all the terms with L, we can formulate the L subproblem as,

$$ { \hbox{min} }_{\text{L}} \left\| {\text{L}} \right\|_{ *} + \frac{{\beta + \beta_{L} }}{2}\left\| {{\text{L}} - \frac{{\beta \left( {{\text{U}} - {\text{S}} - \frac{Z}{\beta }} \right) + \frac{{\beta_{L} }}{2}\left( {E + \frac{{Z_{L} }}{{\beta_{L} }}} \right)}}{{\beta + \beta_{L} }}} \right\|_{F}^{2} $$

(12)

The problem in (12) can be solved by singular value thresholding (SVT) as,

$$ {\text{L}} = {\text{F}}_{{\frac{1}{{\beta + \beta_{L} }}}} \left( X \right)\quad s.t. \, X = \frac{{\beta \left( {{\text{U}} - {\text{S}} - \frac{Z}{\beta }} \right) + \frac{{\beta_{L} }}{2}\left( {E + \frac{{Z_{L} }}{{\beta_{L} }}} \right)}}{{\beta + \beta_{L} }} $$

(13)

where ${\text{F}}_{\varepsilon } \left( X \right) = U{\text{S}}_{\varepsilon } \left( s \right)V^{T}$, and $UsV^{T}$ is the singular value decomposition of X. The S subproblem can be solved in the same way,

$$ {\text{S}} = {\text{S}}_{{\frac{\lambda }{\beta }}} \left( {Y_{S} } \right) \quad s.t. \quad Y_{S} = {\text{U}} - {\text{L}} + \frac{\beta }{Z} $$

(14)

U Subproblem

Combine the terms related to U from the augmented Lagrange function, we have the U subproblem,

$$ { \hbox{min} }_{\text{U}} \frac{\beta }{2}\left\| {{\text{U}} - \left( {{\text{L}} + {\text{S}} + \frac{Z}{\beta }} \right)} \right\|_{F}^{2} + \frac{\mu }{2}\left\| {{\text{HU}} -\Phi } \right\|_{F}^{2} $$

(15)

Equation (15) can be solved directly as it is a convex problem,

$$ {\text{U}} = \left( {\mu {\text{H}}^{T} {\text{H}} + \beta } \right)\left[ {\beta \left( {{\text{L}} + {\text{S}} + \frac{Z}{\beta }} \right) + \mu {\text{H}}^{T}\Phi } \right] $$

(16)

E Subproblem

Similarly, combine all the E-related term, we can formulate the E subproblem,

$$ { \hbox{min} }_{\text{E}} \,\alpha J\left( {\text{E}} \right) + \frac{{\beta_{L} }}{2}\left\| {{\text{E}} - \left( {{\text{L}} - \frac{{Z_{L} }}{{\beta_{L} }}} \right)} \right\|_{F}^{2} $$

(17)

The above problem can be solved by gradient descent algorithm as shown below,

$$ {\text{E}} = {\text{E}} - dx\left( {\alpha \xi_{NLTV} \left( {\text{E}} \right) + \beta_{L} \left( {{\text{E}} - {\text{L}} + \frac{{Z_{L} }}{{\beta_{L} }}} \right)} \right) . $$

(18)

where $\xi_{NLTV}$ denotes the Euler-Lagrange of $J\left( {\text{E}} \right)$ as shown in Eq. (19), $dx$ is the step length in the gradient descent algorithm,

$$ \xi_{NLTV} \left( {\text{E}} \right) = \sum\nolimits_{a,b \in \varOmega } {\left( {{\text{E}}\left( a \right) - {\text{E}}\left( b \right)} \right)w\left( {a,b} \right)} \left[ {\frac{1}{{\left| {\nabla_{NLTV} {\text{E}}\left( a \right)} \right|}} + \frac{1}{{\left| {\nabla_{NLTV} {\text{E}}\left( b \right)} \right|}}} \right] $$

(19)

The Lagrange multipliers are updated by,

$$ Z^{k + 1} = Z^{k} - \beta \left( {{\text{U}} - {\text{S}} - {\text{L}}} \right) $$

(20)

$$ Z_{L}^{k + 1} = Z_{L}^{k} - \beta \left( {{\text{L}} - {\text{E}}} \right) $$

(21)

3 Experiments

Our experiment consists of three parts based on three types of the datasets. Firstly, we apply the proposed algorithm to the simulated ventricular pacing dataset provided by Karlsruhe Institute of Technology (KIT), the simulation dataset include pacing at eight different pacing positions on the left ventricle (LV) and the right ventricle (RV). After that, the LRNLTV algorithm is used to reconstruct the real intervention pace data including two healthy hearts paced with a catheter device which have 21 cases of 551 different heartbeats shared in the online database EDGAR. Finally, we tried the LRNLTV algorithm on 5 real patients with PVC, and the data was collected by the ethical committee of Zhejiang Provincial Hospital. In the three experiment, we compare the reconstruction result of EEP map and activation sequence map by low rank algorithm, TV algorithm and LRNLTV algorithm.

3.1 Simulation Experiments

As we can locate the disease from the potential distribution on EEP map and the activation time distribution on activation map, and accurately determining the location of heart disease is very important for subsequent treatment, we reconstruct the EEP map and activation map in this experiment. The EEP ground truth and transfer matrix is given in the simulated dataset and we compute the EEP by $\Phi = {\text{HU}}$. However, the measured body surface potentials are accompanied by a certain degree of noise in real situation, we add 30 dB Gaussian white noise M to the BSP as $\Phi = {\text{HU}} + {\text{M}}$.

3.1.1 Spatial EEP Reconstruction

In order to compare the reconstruction performance of the three methods (LR, TV and LRNLTV), we use these three methods to reconstruct EEP maps respectively. Take a pacing site at left ventricle apex as an example, we reconstructed the spatial EEP distribution at four moments (20 ms, 60 ms, 100 ms, 160 ms) as shown in Fig. 1, all three methods reconstructed the pacing point at the apex of left ventricle and spreading to the remote area of the lateral left ventricle. However, for the boundary region of low potential and high potential, the reconstruction result of TV for edge information is poor, while LR algorithm have a better result but the boundary is not particularly clear yet. It is worth mentioning that the proposed LRNLTV algorithm have the best performance that we can clearly see the boundary in the reconstructed EEP map at four moment.

Fig. 1.

The reconstruction of EEP at four moments (20 ms, 60 ms, 100 ms, 160 ms) using LRNLTV, LR and TV, the pacing site is at apex LV.

Correlation coefficient (CC) [12] and structural similarity (SSIM) is used to analyze the results quantitatively, while CC and SSIM is defined as,

$$ {\text{CC}}_{i} = \frac{{Cov\left( {u_{ri} ,u_{ti} } \right)}}{{\sqrt {D\left( {u_{ri} } \right)} \sqrt {D\left( {u_{ti} } \right)} }} $$

(22)

$$ {\text{SSIM}}_{i} \left( {{\text{x}},{\text{y}}} \right) = \frac{{\left( {2v_{x} v_{y} + c_{1} } \right)\left( {2Cov\left( {x,y} \right) + c_{2} } \right)}}{{\left( {v_{x}^{2} + v_{y}^{2} + c_{1} } \right)\left( {D\left( x \right)^{2} + D\left( y \right)^{2} + c_{2} } \right)}} $$

(23)

where $u_{ri}$ and $u_{ti}$ are the i-th frame of the truth and the reconstructed result, $Cov\left( {u_{ri} ,u_{ti} } \right)$ is the covariance, $D\left( { \ldots \ldots } \right)$ denotes the variance, $v_{x}$ and $v_{y}$ represent the mean value which are constant avoiding denominator be 0. Error histogram of 8 pacing sites is shown in Fig. 2, the proposed LRNLTV algorithm has the highest CC and SSIM which represent the best reconstruction performance in three methods.

Fig. 2.

The CC and SSIM of simulated EEP reconstruction

3.1.2 Activation Reconstruction

Based on the result of EEP map, we get the activation map which the pacing site is at apex LV illustrated in Fig. 3, while the activation time of a node on the heart is the time when the potential drop at maximum. And the earliest activated point is considered the pacing site, we also calculate the coordinate of the pacing site of ground truth and the result reconstructed by three methods as Truth (81.82, −93.53, 1.19), LRNLTV(81.5, −94.6, 0.19), LR(81, −92.1, 0.09), TV(80.13, −92.13, 0.28), and the locating error of the methods are 1.5 mm, 1.98 mm and 2.37 mm.

Fig. 3.

The ground truth and the reconstruction of activation using LRNLTV, LR and TV, the earliest activated point is at apex LV.

As illustrated in Fig. 3, all three methods suggest that the earliest activated point at apex LV, but LR and TV are less effective in reconstructing details and boundary information, and Table 1 also indicates that the proposed LRNLTV have the highest CC, SSIM and the lowest locating error at different pacing sites.

Table 1. CC, SSIM and locating error between reconstructed activation results and ground truth

3.2 Real Experiments

3.2.1 Interventional Premature Data

In this part of the experiment, we reconstruct the activation map of interventional premature data which include BSP measured by 120 electrodes sampled at 2 kHz, transfer matrix and the gold standard Carto System measured pacing site coordination (x, y, z) of 21 cases of 551 heartbeats. We apply the aforementioned 3 methods to the reconstruction process, two cases of reconstructed EEP map and 14 locating errors of LV and 7 locating errors of RV are displayed in Fig. 4.

Fig. 4.

Location error of LRNLTV, LR, TV compared with Carto System. (Color figure online)

As illustrated in Fig. 4, the position of red point is the ground truth record by Carto System, the white point is the pacing site reconstructed by LRNLTV method, the black point is the location result of LR method and the yellow point is the result of TV method. In the two examples shown in Fig. 4, all the three method indicates that the pacing site at RVOT anterior in case 1 and the LV apex in case 2 while the LRNLTV method has the smallest locating error of 3.78 mm and 6.62 mm. Besides, the error histogram also points out that the proposed LRNLTV method have the highest locating accuracy while the mean value of the locating error is 8.6 mm Owning to the nonlocal similarity in the whole map, details are retained, as TV and LR have larger locating errors.

3.2.2 Clinical PVC Data

To verify the feasibility of the proposed LRNLTV algorithm for clinical data, we apply the algorithm to data of 5 real patients with PVC which consist of the thoracic CT scan whose axial spatial resolution is 0.6–1 mm and the 64-lead ECG signal with a sample frequency of 2 kHz, while CT scan also recorded the location of 64-lead electrodes. Before reconstruction, the relationship between torso and heart is needed, the heart model can be obtained by the segmentation of CT slices, and we can also get the coordinates of the 64-lead ECG electrodes from CT scan to build the torso model. By matching the torso model and the heart model, we can get the transfer matrix.

Figure 5 shows the activation map reconstruction of a male PVC patient aged 33. ECG of 3000 ms to 4000 ms is used to the reconstruction. Because the ECG signal contains a lot of noise, the signal is preprocessed by wavelet transform before reconstruction.

Fig. 5.

Activation reconstruction of a male PVC patient aged 33 compared with Ensite 3000

As illustrated in Fig. 6, the earliest activated site is located at the RVOT septum near the ventricular septal side, it is consistent with the position measured by the gold standard Ensite 3000 System which is invasive. For the other 4 patients, the pacing site is located at RVOT anterior, posterior, free wall and LV apex, which is in concert with Ensite 3000 results and the proposed method is effective in locating the origin of PVC.

Fig. 6.

Locating error and SSIM of different $\mu$ and $\beta = 0.02, \beta_{L} = 0.01$

4 Discussion

We have presented a novel low rank and non-local TV method for EEP reconstruction. By taking advantage of the spatiotemporal information and the nonlocal similarity of EEP, LRNLTV perform better both in locating accuracy and the boundary compared with the LR and TV method, while the background part of LR method have noise inevitably, and the TV method only uses the local similarity of EEP, the temporal correlation between frames get lost. ALM is applied to solve the optimization problem, many parameters are introduced due to the number of constraints, the parameter adjustment is quite important for results, and the computing time is also significant for application of the proposed method.

4.1 Parameter Adjustment

In this paper, the main parameters are $\lambda , \mu , \alpha , \beta , \beta_{L}$. $\lambda$ is the weight coefficient fixed at $1/\sqrt {{ \hbox{max} }\left( {I,J} \right)}$ where (I, J) is the size of solution [13]. $\mu , \beta , \beta_{L}$ are the weight coefficient of fidelity term, Fig. 7 shows the locating error and mean value of SSIM change for a case of LV lateral-epi pace, and $\mu$ is fixed at [0.03, 0.1], $\beta , \beta_{L}$ is in the range of [0.0002, 0.03], and $\alpha$ is the weight coefficient of NLTV constraint fixed at 5.

4.2 Program Running Time

All the process of optimization and solution is computed by MATLAB 2018a of windows version with 32 GB RAM and a 10 core 3.3 GHz processor. As shown in Table 2, we record the computing time of 3 different datasets, LR method have the lowest computing time. In particular, LRNLTV method have the highest reconstruction accuracy, while the time cost is acceptable compared with the other two methods.

Table 2. The computing time of three different datasets with LRNLTV, LR and TV method

5 Conclusion

In this paper, we have presented an innovative LRNLTV method for reconstructing EEP from BSP while both spatiotemporal correlation and nonlocal similarity of EEP are considered. The proposed method can improve the reconstruction accuracy with moderate computing time compared with LR and TV method when applied to simulated and real data, and the reconstruction of clinic data proves the effectiveness of the method.