TGV was proposed to preserve high order of smoothness in image processing problems (Bredies et al., 2010). Based on the assumption that the underlying image is piecewise polynomial, TGV exploits sparsity of high order derivatives along the x-axis and the y-axis. For the illustrative purpose, we display in Figure 1 various piecewise polynomials defined on a plane with degree up to two. Given a 2D image u twice continuously differentiable on a bounded set Ω¯⊂ℝ2, the second order TGV of u with the coefficient α = (α₁, α₂) can be defined as the following infimal convolution (Bredies et al., 2010; Guo et al., 2014)

where ∇ is the 2D gradient operator, p is an auxiliary variable, and the operator ℰ~ is defined by

Here the ℓ₁-norm of a matrix treats a matrix as a vector, i.e., ‖X‖1=∑i,j|Xi,j|. Different from the Laplacian operator which only involves all unmixed second partial derivatives, the second-order TGV involves all partial derivatives, similar to Hessian. In (3), when ∇u is equal to p, the first term in the objective function becomes zero and ℰ~ becomes the Hessian of u. Therefore, one can see that TGV(u)≤‖H(u)1 where H(u) is the Hessian of u. This suggests that TGV could yield a faster minimizing sequence than ‖H(u)‖1, therefore it is a better choice as a regularization term for imposing sparsity than the ℓ₁-norm of Hessian in terms of convergence rate.

Since the cortex surface has complicated geometries and topological structures, it is crucial to choose an appropriate regularization tailored to such kind of irregular surfaces. We discretize the cortex surface to be a 3D triangular mesh Ω and define a voxel-based TGV (vTGV) regularization on it. In order to define directional derivatives on triangular mesh, we treat the centroid of each triangular voxel as a dipole. Since each voxel has three voxels connected, three directional derivatives on ℝ³ can be used to define “gradient” of the density function u. Consider a triangular voxel Λ ∈ Ω, which is homeomorphic to ℝ², we assume that q₁, q₂, q₃ are three normal directions along three edges for Λ, where qi∈ℝ3 depends on the shape of the triangle Λ. For instance, Figure 2 illustrates three normal directions associated with a triangular voxel. Although not perpendicular to each other, these three directions can span the tangent plane through each voxel and thereby can be used to fully describe variations of u. The gradient of u restricted on Λ is defined by

Note that this definition is in the local sense and it can be considered as an extension of the gradient operator in ℝ² into the gradient in a 2D manifold. Given a differentiable function p = (p₁, p₂, p₃), the operator ℰ acting on p restricted to Λ is defined by

This operator can be considered as an extension of ℰ~ in Equation (4) tailored to the triangular mesh Ω.

Next, we discuss the discretization of the operators ∇^ and ℰ. On the triangular mesh Ω with M voxels, we first index all voxels and then define a finite difference operator matrix D ∈ ℝ^3M×M as follows. The (i, j)-th entry of D is defined as

where the voxel index is l = ⌈i/3⌉ ∈ {1, …, M}, i.e., the smallest integer no less than i/3, and k_l,1, k_l,2 and k_l,3 are the indices of the voxels adjacent to the l-th voxel. Based on the definition in Equation (6), the discretization of the operator ℰ is defined as E ∈ ℝ^3M×3M of the form

where ⊗ is the Kronecker product of two matrices. Note that each edge is counted twice in Equation (7) so that the operator E can be easily constructed by using D. Moreover, E is symmetrized by taking the average between D^ and its transpose.

One can see that ∇^u is discretized by Du, and ℰ(p) is discretized by Ep. Once D and E are available, TV and the second order vTGV with the coefficients α₁ and α₂ can be defined as

In Equation (10), the parameters α₁ and α₂ balance the first and second order derivative information of the image (Papafitsoros and Valkonen, 2015). It has been proven that for a large ratio α₂/α₁, the second order TGV coincides with TV under certain conditions (Papafitsoros and Valkonen, 2015).

TV is able to well preserve the edges of images, but is known to create piecewise constant result even in regions with smoothly changed intensities (Benning et al., 2013). By considering higher-order derivative information, TGV generalizes TV and is able to reduce staircasing effects by assuming that the image to be reconstructed is piecewise polynomial (including piecewise constant, piecewise linear, piecewise quadratic, etc.)(Bredies and Holler, 2014). In particular, the proposed second order vTGV assumes that the underlying current density distribution is piecewise linear, and thereby this regularization is able to enforce the sparsity of second spatial derivatives. Although a natural image may have higher order smoothness, it is usually sufficient to use the second order vTGV in practice, since performance enhancement is limited but more computations are required for higher order vTGV. Therefore, we only use the second order vTGV regularization in this work.

In order to improve the spatial resolution of the brain image to better separate close sources, we can incorporate sparsity constraint into the model. A natural strategy to impose sparsity is ℓ₀-norm regularization which minimizes the number of non-zero intensity values in the image. However, since the ℓ₀-regularized problem is computational NP-hard, its ℓ₁-norm relaxed version is usually considered in practice. Recently ℓ₁₋₂ regularization has been proposed (Esser et al., 2013; Lou et al., 2014; Yin et al., 2014), and has been shown to provide a sparser result than the widely used ℓ₁-norm regularization.

For a real positive number p, the ℓ_p-norm of u ∈ ℝ^𝕄 is defined as

Different from the ℓ_p-norm, the ℓ₁₋₂ regularization penalty function is defined as

which has shown potential in image processing and compressive sensing reconstruction (Lou et al., 2014; Yin et al., 2015) in terms of sparsity and fast convergence. It promotes sparsity of an image, and achieves the smallest value when only one voxel in the image is non-zero.

We further discuss the sparsity property of ℓ₁₋₂ regularization from the optimization point of view. Consider a minimization problem in 2D minx∈ℝ2R(x) subject to the linear constraint Ax = b where R(x) is a regularization function. To solve the problem graphically, we need to find the level curve of minimum radius to the origin that intersects with the line L : Ax = b. Figure 3 illustrates the solutions when R is ℓ₂, ℓ₁, ℓ_0.001 (used to approximate ℓ₀) and ℓ₁₋₂ when β = 1, respectively. As shown in Figure 3A, the ℓ₂-regularized solution rarely has zero components, indicating that the solution is usually non-sparse. The ℓ₁-regularized solution may not be sparse if the line L is parallel to the level curves. Compared to ℓ_p (0 < p < 1) regularization, the ℓ₁₋₂ regularization is more likely to yield a sparse solution due to the curvature of level curves. Therefore, the ℓ₁₋₂ regularization promotes sparser solutions than the other regularizations being compared. In the EEG inverse problem, the brain images to be reconstructed in general have a sparse structure that the number of sources is limited, which motivates us to apply the ℓ₁₋₂ regularization to solve this problem.

In this paper, we unify the ℓ₁ type and the ℓ₁₋₂ type regularizations by allowing β = 0 in Equation (12), so that the sparsity regularization term could be adjusted by tuning the parameter β.

The following model is proposed to reconstruct the EEG brain image u

where vTGV(α1,α2)2(u) is defined in Equation (10), and ‖u‖_1−2,β is defined in Equation (12). Here α_i > 0 are regularization parameters which control the contribution of each regularization term. Note that if β = 0, the ℓ₁₋₂ regularization reduces to the ℓ₁ regularization. If we require p = 0, then the vTGV regularization reduces to the TV.

Since the dual norm of ‖·‖₂ is itself, i.e., ‖u‖2=‖max‖q2≤1〈u,q〉, the model Equation (13) can be reformulated as

Next we apply the DCA (Tao and An, 1997) to obtain the following two subproblems

In particular, the second subproblem can be solved efficiently using ADMM. By the change of variables, it can be further written as

By introducing the scaled multipliers x~,y~,z~, we have the following augmented Lagrangian function

Note that this version is equivalent to the standard augmented Lagrangian function up to scaling of multipliers. We group the variables u, p, x, y, z into three blocks, i.e., u, p and (x, y, z). Then the ADMM yields the following algorithm

Moreover, u and p can be solved explicitly as follows

In addition, due to the separability of variables, the (x, y, z)-subproblem boils down to three independent subproblems with respect to x, y and z, respectively, each of which has a closed-form solution represented by proximal operators. For example, the z-subproblem can be solved by using the proximal operator of ℓ₁-norm

where prox_γ(x) = sign(x)⊙max{|x|−γ, 0} with componentwise multiplication ⊙, also known as shrinkage operator. Combining DCA for problem Equation (15) and ADMM for the (u, p)-subproblem, we obtain the algorithm summarized in Algorithm 1.

Note that in this study the entire matrix A is scaled by multiplying 10⁵ in order to reduce round-off errors. Algorithm 1 terminates when either the maximal number of iterations or the minimal relative change is reached. Note that there are two loops in the algorithm: outer and inner loop. In our experiments, the maximum number of iterations for each inner loop is set to be 40, and the maximum number of outer loop is set to be 10. The algorithm will also be halted if the relative change of u is smaller than 10⁻³. Here the relative change of u is defined as

In general, ADMM is simple to implement with linear convergence even if part of the objective function is non-differentiable. Our empirical experience shows that the ℓ₁₋₂ regularization further promotes faster convergence of the algorithm than its ℓ₁-regularized counterpart.

In the proposed Algorithm 1, the regularization parameters α₁, α₂, α₃ are selected to make a balance between smoothness and sparsity. Based on our large numbers of experiments, the optimal parameter selection does not change significantly as the source number, size or configuration changes. For different noise levels, the regularization parameters need to be tuned smaller when SNR increases. Table 1 lists all values of α₁ that we use for the synthetic data sets with SNR between 0 and 30 dB. For simplicity, we set α₂ to be equal to α₁. A more detailed discussion about the influence of ratio α₂/α₁ on the reconstruction results can be found in (Papafitsoros and Valkonen, 2015). For α₃, we find that the performance of the proposed method is not sensitive to α₃ as long as it is in the range of α₃ = 0.1 ~ 0.5α₁. Figure 5 illustrates the source reconstruction results with different values of α₃, where we can see that the results look very similar. By taking a careful look at the bottom source, one can see that α₃ = 0.1α₁ yields slightly under-focalized result, while α₃ = 0.4α₁ yields slightly over-focalized result, so α₃ = 0.2 or 0.3 α₁ provides results closest to the ground truth. In our experiment, we fix α₃ to be 0.3α₁ in all test cases. As for the parameter ρ, which controls the convergence speed of Algorithm 1, it is set to 10α₁ by default. For real data sets, we use the same parameters for the same noise level as the synthetic data.

The parameter β in the ℓ₁₋₂ regularization term varies from 0 to 1. When β = 0, the ℓ₁₋₂ regularization becomes the ℓ₁ regularization. In Figure 6, we study the effect of β on the source reconstruction results. Figure 6B shows the change of reconstruction error with different values of β, where we can see that the larger β is, the smaller the reconstruction error will be. When β = 1, the highest reconstruction accuracy is achieved. Figure 6C shows the change of the sparsity term as iteration increases. One can see that comparing to β = 0 (ℓ₁ regularization), β = 1 (ℓ₁₋₂ regularization) helps to promote sparsity. Notice that the sparsity term will decrease rapidly from one inner loop to another since the variable q is redefined in each outer loop. In our experiments, the maximal number of iterations at each inner loop is set to 40. At each inner loop, the solution becomes convergent and stable within the tolerance, so does the sparsity term. Then at the iteration 41, the updated q results in the refinement of the solution and a large drop of the sparsity term (Figure 6C). In sum, β = 1 not only helps reduce the reconstruction error, but also enhances the sparsity term. Therefore, we set β to 1 in the following study.

In our simulation, source is synthesized using the Gaussian-tapered patch. Firstly, a source center is seeded on the cortex surface, then its neighbors are gradually recruited to make a patch. Because the Gaussian function has bell shape, the source intensity distribution reaches a peak at the center and gradually decreases to zero as it moves away from the center. To model different source configurations, we use Gaussian functions with different variations (σ²), illustrated in Figure 4. As σ² goes to infinity, the intensity of the source decays more and more slowly from the center to its neighbors and approximates the constant function.

In addition to various source configurations, we test a variety of sources with different sizes. Specifically, we use the sources containing 100~300 triangular voxels, which corresponds to 1.4~2.2 cm in radius. To study the sensitivity of the result to the measurement noise, we add i.i.d. additive white Gaussian noise to each channel. We also study the influence of the brain noise by adding i.i.d. Gaussian additive noise to the voxel space. As a widely used criterion for noise level measurement, the signal-to-noise ratio (SNR) is defined as

where P_signal and P_noise are the power of the signal and the noise, respectively. In our simulation, SNR is set to 20 dB by default. The effect of different noise levels is also studied by using signals of SNR 0~20 dB. The synthetic signal is normalized to make sure that the amplitude of the signal falls into the range from 10 to 100 μV, which is the typical EEG signal amplitude of an adult human (Aurlien et al., 2004). For synthetic data, we use the head model template provided by Fieldtrip (Oostenveld et al., 2003), where the number of voxels M is equal to 10240.

To evaluate the performance of the proposed method in realistic scenario, we collected two P300 event-related potentials (ERPs) via auditory and visual oddball paradigms, in which a subject detected an occasional target stimulus in a regular train of sensory stimuli. The experiment was conducted with the approval of institutional review board at Hualien Tzu Chi General Hospital, Taiwan (IRB 101-102) with written informed consent from the subject.

P300 is a positive peak occurring about 300 ms or more after a stimulus (Linden, 2005), which reflects information processing associated with attention and memory. In the auditory stimulation setting, two audio signals of 1500 Hz (target, 40 trials) and 1000 Hz frequency (non-target, 160 trials) were randomly presented to the subject. In the visual stimulation setting, two different pictures of a fierce shark (40 trials) and of an old man (160 trials) were randomly presented to the subject. The subject was required to detect the targets by silently counting these events. A 64 channels EEG machine (ANT Neuro, Enschede Netherlands) was used to record the neural signals. The EEG data was sampled at 512 Hz, filtered by a band pass filter of 0.5–30 Hz and was referenced to the average of all channels. In the end, the average was taken across the trials in order to improve the SNR, and the difference between the target and non-target was used for source localization.

In addition to EEG data, high-resolution MRI data (General Electric, Waukesha, WI, USA) were obtained from the subject for realistic head model construction (Oostenveld et al., 2011). We first segmented the head into three layers, i.e., scalp, skull and brain, and then constructed a triangular mesh for each layer (Oostendorp and van Oosterom, 1991; Fuchs et al., 2002). The cortex surface was also triangulated into a fine mesh with 16384 triangles, each corresponding to a potential dipole source. Finally, BEM (Oostendorp and van Oosterom, 1991; Fuchs et al., 2002) was used to calculate the lead field matrix.

Figure 8A illustrates the source localization results of two sources in nearly noiseless (30 dB) and noisy (0 dB) cases. In the nearly noiseless case, MNE successfully locates these two sources but produces a few spurious sources. For TV-ℓ₁ method, although we can see a little magnitude variation in the edge of the sources, the main area of the sources still shows almost uniform current density distribution. Compared to the other two methods, the proposed method shows the closest result to the ground truth, where the magnitude of the current density varies smoothly from the peak to its neighbors. From the noisy case, one can see that the imaging result is sensitive to measurement noise, especially for the bottom source. MNE shows a lot more spurious sources than the nearly noiseless case even after thresholding. The TV-ℓ₁ method shows an enlarged coverage of the bottom source compared to the ground truth. In addition, one can see that the source intensity becomes more flat in the noisy case. The proposed method is more robust to the noise with the coverage of the bottom source shrinks slightly.

To quantify the influence of noise levels on the source reconstruction performance, we test various noise levels and evaluate the results with the criteria defined in Section 2.4. In order to avoid inconsistency due to different noise configurations, we repeat the experiment 50 times by adding random noise and display the averaged result and the standard deviation in Figure 8B. Generally, the performance of all the methods is improved as SNR increases. From the TRE plot, one can see that our method has the smallest total reconstruction error compared to the other two methods. The LE plot shows that the proposed method has the smallest localization error. Compared to the proposed method, the TV-ℓ₁ method has relatively large localization error since it tends to produce an almost uniform current density and thereby has difficulty locating the peak of the source. In the DF plot, both TV-ℓ₁ method and the proposed method show very high focalization degree, this is because they incorporate ℓ₁ or ℓ₁₋₂ regularization to impose sparsity on the source current density. Taken together, the proposed method shows good performance for all three quantitative criteria at every noise level.

In this section, we study the influence of brain noise by adding i.i.d. Gaussian additive noise to each voxel. Figure 9A shows the source imaging results in the nearly noiseless (30 dB) and noisy (0 dB) cases. Note that in this figure the imaging results are not thresholded so as to better visualize the influence of brain noise. In the nearly noiseless case, MNE produces much less spurious sources under the brain noise than under the measurement noise (Figure 8A), indicating that the spurious sources are mainly due to the measurement noise. For the TV-ℓ₁ method, the reconstructed intensity distribution is generally piecewise constant, but we can see that the intensity variation is larger than the result in Figure 8A. The proposed method produces an accurate source intensity distribution that is very close to the ground truth. In the noisy case, generally the performance of all the methods is affected by the noise. The MNE result shows more background activities due to the high level of noise. The TV-ℓ₁ result shows smaller intensity variation than the nearly noiseless case. For example, for the bottom source, we can see four different intensity colors in the nearly noiseless case, but only two different intensity colors in the noisy case. Compared to the TV-ℓ₁ method, the proposed method provides a smoother result. We can see that the source intensity is weakened due to the high noise level.

Figure 9B further quantifies the results using different noise levels. The TRE plot shows that the proposed method has the smallest reconstruction error. In addition, by comparing to the result in Figure 8B with the same noise level, one can see that the reconstruction error under brain noise is smaller, which is consistent with the visualization result. The LE plot shows that the proposed method has the smallest localization error. It is worth noting that the localization errors of all the methods are smaller than those with measurement noise (Figure 8B). Finally, in the DF plot, both the proposed method and the TV-ℓ₁ method achieve high focalization degree. The focalization degree for MNE is much higher than that under measurement noise. In summary, we observe that the brain imaging result is less sensitive to brain noise than to measurement noise. The proposed method demonstrates robust performance under various levels of brain noise.

In addition to noise level, we also investigate the influence of the source size on the reconstruction results. Figure 10A illustrates the reconstructed brain image with two sources of different sizes. In MNE, although it locates these two sources at the approximate locations, however, it is difficult to differentiate the smaller source from the large numbers of spurious sources. TV-ℓ₁ method recovers both sources clearly without spurious sources, but the coverage of the reconstructed sources is enlarged, especially for the small source on the top. Additionally, it fails to recover the intensity variation of the source in space. In contrast, the proposed method accurately reconstructs the size and intensity variation of these two sources.

Figure 10B shows the quantitative results of various source sizes, where the x-axis represents the number of voxels contained in the simulated sources. TRE plot shows that the proposed method has the smallest reconstruction error, which is insensitive to the source size. In the LE plot, the proposed method shows the smallest localization error. As the source size increases, its localization error becomes slightly smaller, which implies that the proposed method has advantages of dealing with larger sources. The TV-ℓ₁, by contrast, shows relatively large localization error due to the uniform intensity of the reconstructed source. In the DF plot, the proposed method demonstrates very high focalization degree. In summary, the proposed method shows consistent outstanding performance over the other two methods regardless of the source size.

We study the performance of the proposed method using sources with different decay speeds (see Figure 4). In Figure 11A we show two sources of different configurations: the top source decays fast as it goes far from the center while the bottom source decays slowly. From the reconstruction results, one can see that the MNE is not able to tell the configuration difference between these two sources. The TV-ℓ₁ method models the source intensity to be piecewise constant, so both of the reconstructed sources decay very slowly. As for the proposed method, we can tell that the bottom source decays more slowly than the top one.

We further evaluate the performance of the methods with different source configurations quantitatively. In Figure 11B, the x-axis represents the variance σ² of the Gaussian function (Figure 4), so the source intensity decays faster and faster from left to right. The TRE plot shows that the proposed method has the smallest reconstruction error among all the methods. By comparing the results of different variance σ², one can see that the proposed method favors smoother sources whose intensity decays faster, i.e., smaller σ². In the LE curve, the proposed method shows much smaller localization error than the other two methods. Again, one can see that the smoother sources have smaller localization errors. Finally, the DF plot shows that the focalization degree does not rely on the source configurations too much. In sum, the proposed method outperforms the other two methods consistently for all three criteria. Compared to constant sources, it favors smoother sources.

To systematically evaluate the performance of the proposed method for different source locations, we randomly select 50 locations in the whole source space, and test its average performance. Figure 12 displays the whisker plot of the quantitative results, where the lower quartile, median and upper quartile are shown. In TRE plot, the proposed method shows the best median reconstruction accuracy. The range of the results is relatively large which indicates the performance varies at different locations. The LE plot shows that the localization error of the proposed method has a median value of around 1 cm, which is the smallest among all the methods. In addition, the range of its localization error is also the smallest. From the DF plot, one can see that the median focalization degree of the proposed method is ~ 97% which is the highest. All in all, the proposed method shows the best average performance for different source locations among all the compared methods.