**A Regularized Weighted Smoothed** *L***<sup>0</sup> Norm Minimization Method for Underdetermined Blind Source Separation**

### **Linyu Wang, Xiangjun Yin, Huihui Yue \* and Jianhong Xiang**

College of Information and Communication Engineering, Harbin Engineering University, Harbin 150001, China; wanglinyu@hrbeu.edu.cn (L.W.); yinxiangjun@hrbeu.edu.cn (X.Y.); xiangjianhong@hrbeu.edu.cn (J.X.) **\*** Correspondence: yuehuihui@hrbeu.edu.cn

Received: 22 October 2018; Accepted: 30 November 2018; Published: 4 December 2018

**Abstract:** Compressed sensing (CS) theory has attracted widespread attention in recent years and has been widely used in signal and image processing, such as underdetermined blind source separation (UBSS), magnetic resonance imaging (MRI), etc. As the main link of CS, the goal of sparse signal reconstruction is how to recover accurately and effectively the original signal from an underdetermined linear system of equations (ULSE). For this problem, we propose a new algorithm called the weighted regularized smoothed *L*0-norm minimization algorithm (WReSL0). Under the framework of this algorithm, we have done three things: (1) proposed a new smoothed function called the compound inverse proportional function (CIPF); (2) proposed a new weighted function; and (3) a new regularization form is derived and constructed. In this algorithm, the weighted function and the new smoothed function are combined as the sparsity-promoting object, and a new regularization form is derived and constructed to enhance de-noising performance. Performance simulation experiments on both the real signal and real images show that the proposed WReSL0 algorithm outperforms other popular approaches, such as SL0, BPDN, NSL0, and L*p*-RLSand achieves better performances when it is used for UBSS.

**Keywords:** image reconstruction; nullspace measurement matrix; regularized least squares problem; smoothed *L*0-norm; sparse signal recovery; UBSS; weighted function

### **1. Introduction**

The problem that UBSS [1,2] needs to address is how to separate multiple signals from a small number of sensors. The essence of this problem is to solve the optimal solution of the undetermined linear system of equations (ULSE). Fortunately, as a new undersampling technique, compressed sensing (CS) [3–5] is an effective way to solve ULSE, which makes it possible to apply CS to UBSS.

The model of CS is shown in Figure 1. According to this figure, it can be see that CS boils down to the form,

$$\mathbf{y} = \Phi \mathbf{x} + \mathbf{b},\tag{1}$$

where **<sup>Φ</sup>** = [*φ*1,*φ*2, ...,*φn*] <sup>∈</sup> <sup>R</sup>*m*×*<sup>n</sup>* is a sensing matrix with the condition of *<sup>m</sup> <sup>n</sup>* and *<sup>φ</sup><sup>i</sup>* <sup>∈</sup> <sup>R</sup>*m*, *i* = 1, 2, ..., *n*, which can be further represented as **Φ** = *ψϕ*, while *ψ* is a random matrix, and *ϕ* is the sparse basis matrix. **<sup>y</sup>** <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* is the vector of measurements. Moreover, **<sup>b</sup>** <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* denotes the additive noise.

**Figure 1.** Frame of compressed sensing (CS).

To solve the ULSE in Equation (1), we try to recover the sparse signal **x** from the given {**y**, **Φ**} by CS. According to CS, this problem is transformed into solving the *L*0-norm minimization problem.

$$\|(P\_0)\quad\underset{\mathbf{x}\in\mathbb{R}^n}{\text{arg min}}\|\|\mathbf{x}\|\|\_{\mathcal{O}^\times} \quad \text{s.t.} \ ||\boldsymbol{\Phi}\mathbf{x} - \mathbf{y}||\_2^2 \le \epsilon. \tag{2}$$

where  denotes error. This rather wonderful attempt is actually supported by a brilliant theory [6]. Based on this theory, in the noiseless case, it is proven that the sparsest solution is indeed a real signal when x is sufficiently sparse and **Φ** satisfies the restricted isometry property (RIP) [7]:

$$1 - \delta\_K \le \frac{||\Phi \mathbf{x}||\_2^2}{||\mathbf{x}||\_2^2} \le 1 + \delta\_{K'} \tag{3}$$

where *K* is the sparsity of signal **x** and *δ<sup>K</sup>* ∈ (0, 1) is a constant. In Equation (2), the *L*0-norm is nonsmooth, which leads an NP-hard problem. In practice, two alternative approaches are usually employed to solve the problem [8]:


For greedy search, the main methods are based on greedy matching pursuit (GMP) algorithms, such as orthogonal matching pursuit (OMP) [9,10], stage-wise orthogonal matching pursuit (StOMP) [11], regularized orthogonal matching pursuit (ROMP) [12], compressive sampling matching pursuit (CoSaMP) [13], generalized orthogonal matching pursuit (GOMP) [14,15], and subspace pursuit (SP) [16,17] algorithms. The objective function of these algorithms is given by:

$$\mathop{\arg\min}\_{\mathbf{x}\in\mathbb{R}^{n}}\frac{1}{2}\parallel\Phi\mathbf{x}-\mathbf{y}\parallel\_{2}^{2},\ \text{s.t.}\parallel\mathbf{x}\parallel\_{0}\leq K.\tag{4}$$

As shown in the above equation, the features of GMP algorithms can be concluded as:


The advantage of GMP algorithms is that the computational complexity is low, but the reconstruction accuracy is not high in the noise case.

At present, the relaxation method for *P*<sup>0</sup> is widely used. The relaxation method is mainly divided into two categories: the constraint-type algorithm and the regularization method. The constraint-type algorithm can also be divided into *L*1-norm minimization methods and smoothed *L*0-norm minimization methods. The representative algorithm of the former is the BPalgorithm [18], and the latter is the smoothed *L*0-norm minimization (SL0) algorithm. For the SL0 algorithm, the objective function can be expressed as:

$$\begin{aligned} \left(P\_F\right) \underset{\mathbf{x}\in\mathbb{R}^n}{\arg\min} \, F\_\sigma(\mathbf{x}), \quad \text{s.t.} \; ||\Phi\mathbf{x} - \mathbf{y}||\_2^2 \le \epsilon. \\\ \lim\_{\sigma\to 0} F\_\sigma(\mathbf{x}) = \lim\_{\sigma\to 0} \sum\_{i=1}^n f\_\sigma(\mathbf{x}\_i) \approx ||\mathbf{x}||\_0. \end{aligned} \tag{5}$$

where *Fσ*(**x**) is a smoothed function, which approximates the *L*0-norm when *σ* → 0. Compared with *L*<sup>1</sup> or *Lp*, a small *σ* is selected to make the function close to *L*0-norm [8]; therefore, *Fσ*(**x**) are closer to the optimal solution.

Based on the idea of approximation, Mohimani used a Gauss function to approximate the *L*0-norm [19], which is described as:

$$f\_{\sigma}(\mathbf{x}\_{i}) = 1 - \exp(-\frac{\mathbf{x}\_{i}^{2}}{2r^{2}}).\tag{6}$$

According to the equation, we can know:

$$f\_{\mathcal{F}}(\mathbf{x}\_{i}) \approx \begin{cases} 1 & \text{if } \mathbf{x}\_{i} \gg \sigma \\ 0 & \text{if } \mathbf{x}\_{i} \ll \sigma \end{cases} \tag{7}$$

when *σ* is a small enough positive value, the Gauss function is almost equal to the *L*0-norm. Furthermore, the Gauss function is differentiable and smoothed; hence, it can be optimized by optimization methods such as the gradient descent (GD) method. Zhao proposed another smoothed function: the hyperbolic tangent (tanh) [20],

$$f\_{\mathcal{T}}(\mathbf{x}\_{i}) = \frac{\exp(\frac{\mathbf{x}\_{i}^{2}}{2\sigma^{2}}) - \exp(-\frac{\mathbf{x}\_{i}^{2}}{2\sigma^{2}})}{\exp(\frac{\mathbf{x}\_{i}^{2}}{2\sigma^{2}}) + \exp(-\frac{\mathbf{x}\_{i}^{2}}{2\sigma^{2}})}. \tag{8}$$

This smoothed function makes a closer approximation to the *L*0-norm than the Gauss function, as shown in [19], with the same *σ*; hence, it performs better in sparse signal recovery. Indeed, a large number of simulation experiments confirmed this view.

Another relaxation method is the regularization method. For CS, sparse signal recovery in the noise case is a very practical and unavoidable problem. Fortunately, the regularization method makes the solution of this problem possible [21,22]. The regularization method can be described as a "relaxation" approach that tries to solve the following unconstrained recovery problem:

$$\|(P\_{\upsilon})\mathop{\arg\min}\_{\mathbf{x}\in\mathbb{R}^n} \frac{1}{2} \|\,\boldsymbol{\Phi}\mathbf{x} - \mathbf{y}\|\|\_{2}^{2} + \lambda\upsilon(\mathbf{x}),\tag{9}$$

where *<sup>λ</sup>* <sup>&</sup>gt; 0 is the parameter that balances the trade-off between the deviation term **<sup>Φ</sup><sup>x</sup>** <sup>−</sup> **<sup>y</sup>** <sup>2</sup> <sup>2</sup> and the sparsity regularizer *υ*(**x**). The sparse prior information is enforced via the regularizer *υ*(**x**), and a proper *υ*(**x**) is crucial to the success of the sparse signal recovery task: it should favor sparse solutions and make sure the problem *Pυ* can be solved efficiently in the meantime.

For regularization, various sparsity regularizers have been proposed as the relaxation of the *L*0-norm. The most popular algorithms are the convex *L*1-norm [22,23] and the nonconvex *Lp*-norm to the *p*th power [24,25]. In the noiseless case, the *L*1-norm is equivalent to the *L*0-norm, and the *L*1-norm is the only norm with sparsity and convexity. Hence, it can be optimized by convex optimization methods. However, according to [8], in the noisy case, the *L*1-norm is not exactly equivalent to the *L*0-norm, so the effect of promoting sparsity is not obvious. Compared to the *L*1-norm, the nonconvex *Lp*-norm

to the *p*th power makes a closer approximation to the *L*0-norm; therefore, *Lp*-norm minimization has a better sparse recovery performance [8].

In view of the above explanation, in this paper, a compound inverse proportional function (CIPF) function is proposed as a new smoothed function, and a new weighted function is proposed to promote sparsity. For the noise case, a new regularization form is derived and constructed to enhance de-noising performance. The experimental simulation verifies the superior performance of this algorithm in signal and image recovery, and it has achieved good results when applied to UBSS.

This paper is organized as follows: Section 2 introduces the main work of this paper. The steps of the ReRSL0algorithm and the selection of related parameters are described in Section 3. Experimental results are presented in Section 4 to evaluate the performance of our approach. Section 5 verifies the effect of the proposed weighted regularized smoothed *L*0-norm minimization (WReSL0) algorithm in UBSS. Section 6 concludes this paper.

### **2. Main Work of This Paper**

In this paper, based on the *PF* in Equation (9), we propose a new objective function, which is given by:

$$\underset{\mathbf{x}\in\mathbb{R}^n}{\arg\min} \,\mathbf{W}H\_{\sigma}(\mathbf{x}),\text{ s.t. } ||\Phi\mathbf{x}-\mathbf{y}||\_2^2 \le \epsilon.\tag{10}$$

According to this equation, We not only propose a smoothed function approximating the *L*0-norm, but also propose a weighted function to promote sparsity. This section focuses on the relevant contents of **W** = [*w*1, *w*2, ...*wn*] *<sup>T</sup>* and *Hσ*(**x**).

### *2.1. New Smoothed Function: CIPF*

According to [26], some properties of the smoothed functions are summarized in the following: **Property**: Let *<sup>f</sup>* : <sup>R</sup> <sup>→</sup> [−∞, <sup>+</sup>∞] and, define *<sup>f</sup>σ*(*r*) <sup>≈</sup> *<sup>f</sup>σ*(*r*/*σ*) for any *<sup>σ</sup>* <sup>&</sup>gt; 0. The function *<sup>f</sup>* has the property, if:


It follows immediately from **Property** that { *<sup>f</sup>σ*(*r*)} converges to the *<sup>L</sup>*0-norm as *<sup>σ</sup>* <sup>→</sup> <sup>0</sup>+, i.e.,

$$\lim\_{\sigma \to 0^{+}} f\_{\sigma}(r) = \begin{cases} 0 & \text{if } r = 0 \\ 1 & \text{otherwise.} \end{cases} \tag{11}$$

Based on **Property**, this paper proposes a new smoothed function model called CIPF, which satisfies **Property** and better approximates the *L*0-norm. The smoothed function model is given as:

$$f\_{\sigma}(r) = 1 - \frac{\sigma^2}{ar^2 + \sigma^2}. \tag{12}$$

In Equation (12), *α* denotes a regularization factor, which is a large constant. By experiments, the factor *α* is determined to be 10, which is a good result of the simulation. *σ* represents a smoothed factor, and when it is smaller, it will make the proposed model closer to the *L*0-norm. Obviously,

$$\lim\_{r \to 0} f\_{\sigma}(r) = \left\{ \begin{array}{ll} 0, & r = 0 \\ 1, & r \neq 0 \end{array} \text{ or approximately } f\_{\sigma}(r) \approx \begin{cases} 0, & |\, r| \ll \sigma \\ 1, & |\, r| \gg \sigma \end{cases} \text{ is satisfied. Let:} $$

$$H\_{\sigma}(\mathbf{x}) = \sum\_{i=1}^{n} f\_{\sigma}(\mathbf{x}\_{i}) = n - \sum\_{i=1}^{n} \frac{\sigma^{2}}{a \mathbf{x}\_{i}^{2} + \sigma^{2}} \tag{13}$$

where *Hσ*(**x**) ≈ ||**x**||<sup>0</sup> for small values of *σ*, and the approximation tends to equality when *σ* → 0.

Figure 2 shows the effect of the CIPF model approximating the L0-norm. Obviously, the CIPF model makes a better approximation.

In conclusion, the merits of the CIPF model can be summarized as follows:


These merits make it possible to reduce the computational complexity on the premise of ensuring the accuracy of sparse signal reconstruction, which is of practical significance for sparse signal reconstruction.

**Figure 2.** Different functions used in the literature to approximate the *L*0-norm; some of them are plotted in this figure, and the *L*0.5-norm is displayed for comparison. CIPF, compound inverse proportional function.

### *2.2. New Weighted Function*

Candès et al. [27] proposed the weighted *L*1-norm minimization method, which employs the weighted norm to enhance the sparsity of the solution. They provided an analytical result of the improvement in the sparsity recovery by incorporating the weighted function with the objective function. Pant et al. [28] applied another weighted smoothed *L*0-norm minimization method, which uses a similar weighted function to promote sparsity. The weighted function can be summarized as follows:


From the two weighted functions, we can find a phenomenon: a large signal entry *xi* is weighted with a small *wi*; on the contrary, a small signal entry *xi* is weighted with a large value *wi*. By analysis, the large *wi* forces the solution **x** to concentrate on the indices where *wi* is small, and by construction, these correspond precisely to the indices where **x** is nonzero.

Combined with the above idea, we propose a new weighted function, which is given by:

$$w\_i = \mathbf{e}^{-\frac{|x\_i|}{\sigma}}, \text{ s.t. } i = 1, 2, \dots, n. \tag{14}$$

As for Candès et al., when the signal entry is zero or close to zero, the value of *wi* will be very large, which is not suitable for computation by a computer. Although Pant et al. noticed the problem and improved the weighted function to avoid it, the constant *ζ* depends on experience. Actually, the proposed weighted function can avoid the two problems. Moreover our weighted function can

be satisfied with the phenomenon. When the small signal entry *xi* can be weighted with a large *wi* and a large signal entry *xi* can be weighted with a small *wi*, this can make the large signal entry and small signal entry closer. In this way, the direction of optimization can be kept as consistent as possible, and the optimization process tends to be more optimal. Therefore, the proposed weighted function can have a better effect.

### **3. New Algorithm for CS: WReSL0**

### *3.1. WReSL0 Algorithm and Its Steps*

Here, in order to analyze the problem more clearly, we rewrite Equation (10) as follows:

$$\underset{\mathbf{x}\in\mathbb{R}^n}{\arg\min} \mathbf{W}H\_{\mathcal{F}}(\mathbf{x})\_{\prime} \text{ s.t. } ||\Phi\mathbf{x} - \mathbf{y}||\_2^2 \le \epsilon.$$

where *<sup>H</sup>σ*(**x**) = **<sup>I</sup>** <sup>−</sup> *<sup>σ</sup>*<sup>2</sup> *<sup>α</sup>***x**2+*σ*<sup>2</sup> (**<sup>I</sup>** <sup>∈</sup> <sup>R</sup>*<sup>N</sup>* is a unit vector) is a differentiable smoothed accumulated function. The weighted function **W** = e<sup>−</sup> <sup>|</sup>**x**<sup>|</sup> *<sup>σ</sup>* . Therefore, we can obtain the gradient of CIPF, which is written as:

$$\mathbf{G} = \frac{\partial H\_{\mathcal{C}}(\mathbf{x})}{\partial \mathbf{x}} = \frac{2a\sigma^2 \mathbf{x}}{\left(a\mathbf{x}^2 + \sigma^2\right)^2} \tag{15}$$

According to Equation (15), as in [28], we can obtain:

$$\mathbf{WG} = \left(\mathbf{e}^- \frac{|\mathbf{x}|}{\sigma}\right)^T \frac{2a\sigma^2 \mathbf{x}}{\left(a\mathbf{x}^2 + \sigma^2\right)^2} \tag{16}$$

Solving the problem of ULSE is to solve the optimization problem in Equation (10). As for this problem, there are many methods, such as split Bregman methods [29–31], FISTA [32], alternating direction methods [33], gradient descent (GD) [34], etc. In order to reduce the computational complexity, this paper adopts the GD method to optimize the proposed objective function.

Given *σ*, a small target value *σ*min, and a sufficiently large initial value *σ*max, after referring to the annealing mechanism in simulated annealing [35], this paper proposes a monotonically-decreasing sequence {*σt*|*t* = 2, 3, ..., *T*}, which is generated as:

$$
\sigma\_t = \sigma\_{\text{max}} \theta^{-\gamma(t-1)}, \text{ s.t. } t = 1, 2, 3, \dots, T. \tag{17}
$$

where *<sup>γ</sup>* = log*<sup>θ</sup>* (*σ*max/*σ*min) *<sup>T</sup>*−<sup>1</sup> , *<sup>θ</sup>* is a constant that is larger than one, and *<sup>T</sup>* is the maximum number of iterations. Using such a monotonically-decreasing sequence can avoid the case of too small of a *σ* leading to the local optimum.

Similar to SL0, WReSL0 also consists of two nested iterations: the external loop, which begins with a sufficiently large value of *σ*, i.e, *σ*max, responsible for the gradually decreasing strategy in Equation (17), and the internal loop, which for each value of *σ*, finds the maximizer of *Hσ*(**x**) on {**x**|||*A***x** − **y**||<sup>2</sup> ≤ }.

According to the GD algorithm, the internal loop consists of the gradient descent step, which is given by:

$$
\hat{\mathbf{x}} = \mathbf{x} + \mu \mathbf{d},\tag{18}
$$

where *d* = *g* and *μ* denotes a step size factor. This part is similar to SL0, followed by solving the problem:

$$\underset{\mathbf{x}^\* \in \mathbb{R}^n}{\arg\min} ||\mathbf{x}^\* - \hat{\mathbf{x}}||\_{2^\prime}^2 \quad \text{s.t.} ||\Phi \mathbf{x}^\* - \mathbf{y}||\_2^2 \le \varepsilon \tag{19}$$

where **x**∗ denotes the optimal solution. By regularization, this form can be converted to another form as follows,

$$\underset{\mathbf{x}^\* \in \mathbb{R}^n}{\arg\min} ||\mathbf{x}^\* - \hat{\mathbf{x}}||\_2^2 + \lambda ||\Phi \mathbf{x}^\* - \mathbf{y}||\_2^2. \tag{20}$$

where *λ* is the regularization parameter, which is adapted to balance the fit of the solution to the data *y* and the approximation of the solution to the maximizer of *Hσ*(**x**). Weighted least squares (WLS) can be used to solve this problem, and the solution is:

$$\mathbf{x}^\* = \left[ \begin{bmatrix} \mathbf{I}\_{\mathrm{n}} \\ \mathbf{\Phi} \end{bmatrix}^H \begin{bmatrix} \mathbf{I}\_{\mathrm{n}} & \mathbf{0} \\ \mathbf{0} & \lambda \mathbf{I}\_{\mathrm{m}} \end{bmatrix} \begin{bmatrix} \mathbf{I}\_{\mathrm{n}} \\ \mathbf{\Phi} \end{bmatrix} \right]^{-1} \begin{bmatrix} \mathbf{I}\_{\mathrm{n}} \\ \mathbf{\Phi} \end{bmatrix}^H \begin{bmatrix} \mathbf{I}\_{\mathrm{n}} & \mathbf{0} \\ \mathbf{0} & \lambda \mathbf{I}\_{\mathrm{m}} \end{bmatrix} \begin{bmatrix} \mathbf{\hat{x}} \\ \mathbf{y} \end{bmatrix} . \tag{21}$$

By calculation, Equation (21) is equivalent to:

$$\mathbf{x}^\* = \left(\mathbf{I}\_{\mathrm{il}} + \lambda \boldsymbol{\Phi}^H \boldsymbol{\Phi}\right)^{-1} \left(\hat{\mathbf{x}} + \lambda \boldsymbol{\Phi}^H \mathbf{y}\right) \tag{22}$$

where **I***<sup>n</sup>* and **I***<sup>m</sup>* are both identity matrices of size *n* × *n* and *m* × *m*, respectively. Therefore, we can obtain:

$$\begin{array}{l} \mathbf{x}^{\*} - \dot{\mathbf{x}} = \left(\mathbf{I}\_{n} + \lambda \boldsymbol{\Phi}^{H} \boldsymbol{\Phi}\right)^{-1} \left(\dot{\mathbf{x}} + \lambda \boldsymbol{\Phi}^{H} \mathbf{y}\right) - \dot{\mathbf{x}} \\ = \left(\mathbf{I}\_{\mathrm{H}} + \lambda \boldsymbol{\Phi}^{H} \boldsymbol{\Phi}\right)^{-1} \left(\dot{\mathbf{x}} + \lambda \boldsymbol{\Phi}^{H} \mathbf{y} - \left(\mathbf{I}\_{n} + \lambda \boldsymbol{\Phi}^{H} \boldsymbol{\Phi}\right) \dot{\mathbf{x}}\right) \\ = \left(\mathbf{I}\_{n} + \lambda \boldsymbol{\Phi}^{H} \boldsymbol{\Phi}\right)^{-1} \left(\dot{\mathbf{x}} + \lambda \boldsymbol{\Phi}^{H} \mathbf{y} - \dot{\mathbf{x}} - \lambda \boldsymbol{\Phi}^{H} \boldsymbol{\Phi} \dot{\mathbf{x}}\right) \\ = -\left(\lambda^{-1} \mathbf{I}\_{n} + \boldsymbol{\Phi}^{H} \boldsymbol{\Phi}\right)^{-1} \boldsymbol{\Phi}^{H} \left(\boldsymbol{\Phi} \dot{\mathbf{x}} - \mathbf{y}\right) \end{array}$$

According to the above analysis and derivation, we can get:

$$\mathbf{x}^\* = \hat{\mathbf{x}} - \left(\lambda^{-1}\mathbf{I}\_n + \Phi^H \Phi\right)^{-1} \Phi^H \left(\Phi \hat{\mathbf{x}} - \mathbf{y}\right) \tag{23}$$

The initial value of the internal loop is the maximizer of *Hσ*(**x**) obtained for *σ*max. To increase the speed, the internal loop is repeated a fixed and small number of times (L). In other words, we do not wait for the GD method to converge in the internal loop.

According to the explanation above, we can conclude the steps of the proposed WReSL0 algorithm, which are given in Table 1. As for *σ*, it can be shown that function *Hσ*(**x**) remains convex in the region where the largest magnitude of the component of **x** is less than *σ*. As the algorithm starts at the original value **x**(0) = Φ*H*(ΦΦ*H*)−1**y**, the above choice of *σ*<sup>1</sup> ensures that the optimization starts in a convex region. This greatly facilitates the convergence of the WReSL0 algorithm.

**Table 1.** Weighted regularized smoothed *L*0-norm minimization (WReSL0) algorithm using the GD method.

• Initialization:

(1) Set *L*, *μ* = *σ*/(2*α*), **x**ˆ(0) = **Φ***H*(**ΦΦ***H*)−1**y**.

(2) Set *<sup>σ</sup>*max <sup>=</sup> <sup>√</sup>*<sup>α</sup>* max <sup>|</sup>**x**|, *<sup>σ</sup>*min <sup>=</sup> 0.01, and *<sup>σ</sup><sup>t</sup>* <sup>=</sup> *<sup>σ</sup>*max*θ*−*γ*(*t*−1), where *<sup>γ</sup>* <sup>=</sup> log*<sup>θ</sup>* (*σ*max/*σ*min) *<sup>T</sup>*−<sup>1</sup> , and *<sup>T</sup>* is the maximum number of iterations.

• **while** *t* < *T*, **do** (1) Let *σ* = *σt*. (2) Let **x** = **x**ˆ(*t*−1). **for** *l* = 1, 2, ..., *L* (a) **x** ← **x** − *μ* e<sup>−</sup> <sup>|</sup>**x**<sup>|</sup> *σ <sup>T</sup>* <sup>2</sup>*ασ*2**<sup>x</sup>** (*α***x**<sup>2</sup>+*σ*<sup>2</sup>) 2 . (b) **<sup>x</sup>** <sup>←</sup> **<sup>x</sup>** <sup>−</sup> *<sup>λ</sup>*−1**I***<sup>n</sup>* <sup>+</sup> **<sup>Φ</sup>***H***Φ**−<sup>1</sup> **<sup>Φ</sup>***<sup>H</sup>* (**Φx**<sup>ˆ</sup> <sup>−</sup> **<sup>y</sup>**) (3) Set **x**ˆ(*t*−1) = **x**. • The estimated value is **<sup>x</sup>**<sup>ˆ</sup> <sup>=</sup> **<sup>x</sup>**ˆ(*t*).

### *3.2. Selection of Parameters*

The selection of parameters *μ* and *σ* will affect the performance of the WReSL0 algorithm; thus, this paper discusses the selection of these two above parameters in this section.

### 3.2.1. Selection of Parameter *μ*

According to the algorithm, each iteration consists of a descent step *xi* ← *xi* − *μ* e<sup>−</sup> <sup>|</sup>*xi*<sup>|</sup> *σ* <sup>2</sup>*ασ*<sup>2</sup>*xi* (*αxi* <sup>2</sup>+*σ*2) <sup>2</sup> , 1 ≤ *i* ≤ *n*, followed by a projection step. If for some values of *i*, we have |*xi*| *σ*, then the algorithm does not change the value of *xi* in that descent step; however, it might be changed in the projection step. If we are looking for a suitably large *μ*, a suitable choice is to make the algorithm force all those values of **x** satisfying |*xi*| *σ* toward zero. Therefore, we can get:

$$
\mu \mathbf{x}\_i - \mu \left(\mathbf{e}^- \frac{|\mathbf{x}\_i|}{\sigma}\right) \frac{2a\sigma^2 \mathbf{x}\_i}{\left(a\mathbf{x}\_i^2 + \sigma^2\right)^2} \approx 0\tag{24}
$$

and:

$$\left(\mathbf{e}^{-}\left|\frac{\mathbf{x}\_{i}}{\sigma}\right|\right)\xrightarrow{\mathbf{x}\_{i}\to\mathbf{0}}\mathbf{1}\tag{25}$$

Combining Equations (24) and (25), we can further obtain:

$$\mathbf{x}\_{i} - \mu \frac{2a\sigma^{2}\mathbf{x}\_{i}}{\left(a\mathbf{x}\_{i}^{2} + \sigma^{2}\right)^{2}} \approx 0\tag{26}$$

By calculation, we can obtain:

$$
\mu \approx \frac{\left(\alpha x\_i^2 + \sigma^2\right)^2}{2\alpha \sigma^2} \xrightarrow{x\_i \to 0} \frac{\sigma^2}{2\alpha} \tag{27}
$$

According to the above derivation, we have come to the conclusion that *<sup>μ</sup>* <sup>≈</sup> *<sup>σ</sup>*<sup>2</sup> <sup>2</sup>*<sup>α</sup>* . Therefore, we can set *μ* = *<sup>σ</sup>*<sup>2</sup> 2*α* .

### 3.2.2. Selection of Parameter *σ*

According to Equation (17), the descending sequence of *<sup>σ</sup>* is generated by *<sup>σ</sup><sup>t</sup>* <sup>=</sup> *<sup>σ</sup>*max *<sup>σ</sup>*min *<sup>σ</sup>*max *<sup>t</sup>*−<sup>1</sup> *<sup>T</sup>*−<sup>1</sup> (it is obtained through simplification of Equation (17)). Parameter *σ*min and parameter *σ*max should be appropriately selected. The selection of *σ*min and *σ*max is discussed below.

For the initial value of *σ*, i.e., *σ*max, here, let *x*˜ = max{|**x**|}; suppose there is a constant *b*, in order to make the algorithm converge quickly; let parameter *σ*max satisfy:

$$H\_{\mathcal{V}}(\tilde{\mathbf{x}}) = 1 - \frac{\sigma\_{\text{max}}^2}{a\tilde{\mathbf{x}}^2 + \sigma\_{\text{max}}^2} \le b \Rightarrow \sigma\_{\text{max}} \ge \left(\sqrt{\frac{1 - b}{b}a}\right)\tilde{\mathbf{x}}.\tag{28}$$

From the equation, we can see that constant *b* satisfies <sup>1</sup>−*<sup>b</sup> <sup>b</sup>* ≥ 0; thus 0 < *b* ≤ 1, and here, we define constant *<sup>b</sup>* as 0.5. Hence, *<sup>σ</sup>*max <sup>=</sup> <sup>√</sup>*<sup>α</sup>* max{|**x**|}.

For the final value *σ*min, when *σ*min → 0, *Hσ*min (**x**) → ||**x**||0. That is, the smaller *σ*min, the more *Hσ*min (**x**) can reflect the sparsity of signal **x**, but at the same time, it is also more sensitive to noise; therefore, the value *σ*min should not be too small. Combining [19], we choose *σ*min = 0.01.

### **4. Performance Simulation and Analysis**

The numerical simulation platform is MATLAB 2017b, which is installed on a computer with a Windows 10, 64-bit operating system. The CPU of the simulation computer is the Intel (R) Core (TM) i5-3230M, and the frequency is 2.6 GHz. In this section, the performance of the WReSL0 algorithm is verified by signal and image recovery in the noise case.

Here, some state-of-the-art algorithms are selected for comparison. The parameters are selected to obtain the best performance for each algorithm: for the BPDNalgorithm [36], the regularization parameter *λ* = *σ<sup>N</sup>* 2log(*n*); for the SL0 algorithm [19], the initial value of smoothed factor *δ*max = 2 max{|**x**|}, the final value of smoothed factor *δ*min = 0.01, scale factor is set as step size *L* = 5, and the attenuation factor *ρ* = 0.8; for the NSL0algorithm [20], the initial value of smoothed factor *δ*max = 4 max{**x**}, the final value of smoothed factor *δ*min = 0.01, the step size *L* = 10, and the attenuation factor *ρ* = 0.8; for L*p*-RLSalgorithm [24], the number of iterations *T* = 80, the norm initial value *p*<sup>1</sup> = 1, the norm final value *pT* = 0.1, the initial value of regularization factor <sup>1</sup> = 1, the final value of regularization factor  *<sup>T</sup>* = 0.01, and the algorithm termination threshold *Et* = 10<sup>−</sup>25; for the WReSL0 algorithm, the initial value of smoothed factor *<sup>σ</sup>*max <sup>=</sup> <sup>√</sup>*<sup>c</sup>* max{|**x**|}, the final value of smoothed factor *σ*min = 0.01, the iterations *T* = 30, the step size *L* = 5, and the regularization parameter *λ* = 0.1. All experiments are based on 100 trials.

### *4.1. Signal Recovery Performance in the Noise Case*

In this part, we discuss signal recovery performance in the noise case. We add noise **b** to the measurement vector **y**; moreover, **b** = *δN***Ω**, **Ω** is randomly formed and follows the Gaussian distribution of N (0, 1). For signal recovery under noise conditions, we evaluate the performance of algorithms by the normalized mean squared error (NMSE) and the CPU running time (CRT). NMSE is defined as ||*<sup>x</sup>* <sup>−</sup> *<sup>x</sup>*||2/||*x*||2. CRT is measured with *tic* and *toc*. In order to analyze the de-noising performance of the WReSL0 algorithm in context closer to the real situation, we constructed a certain signal as an experimental object in the experiments in this section. The signal is given by:

$$\begin{cases} \begin{aligned} \mathbf{x}\_1 &= a\_1 \sin(2\pi f\_1 T\_s \mathbf{t}) \\\\ \mathbf{x}\_2 &= \beta\_1 \cos(2\pi f\_2 T\_s \mathbf{t}) \\\\ \mathbf{x}\_3 &= a\_2 \sin(2\pi f\_3 T\_s \mathbf{t}) \\\\ \mathbf{x}\_4 &= \beta\_2 \cos(2\pi f\_4 T\_s \mathbf{t}) \\\\ \mathbf{x}\_5 &= \mathbf{x}\_1 + \mathbf{x}\_2 + \mathbf{x}\_3 + \mathbf{x}\_4 \end{aligned} \end{cases} \tag{29}$$

where *α*<sup>1</sup> = 0.2, *α*<sup>2</sup> = 0.1, *β*<sup>1</sup> = 0.3, and *β*<sup>2</sup> = 0.4. *f*<sup>1</sup> = 50 Hz; *f*<sup>2</sup> = 100 Hz; *f*<sup>3</sup> = 200 Hz; and *f*<sup>4</sup> = 300 Hz. Here, **t** is a sequence with **t** = [1, 2, 3, ..., *n*], and *Ts* is sampling interval with the value of <sup>1</sup> *fs* . *fs* is the sampling frequency with the value of 800 Hz. The object that needs to be reconstructed can be expressed as:

$$\mathbf{y} = \boldsymbol{\Phi}\mathbf{x} + \delta\_N \mathbf{\Omega}.\tag{30}$$

where *<sup>x</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* is a sparse signal in the frequency domain, and it is the Fourier transform expression of <sup>X</sup> , *<sup>y</sup>* <sup>∈</sup> <sup>R</sup>*m*. Here, let *<sup>n</sup>* <sup>=</sup> 128, *<sup>m</sup>* <sup>=</sup> 64. Moreover, **<sup>Φ</sup>** can be represented as **<sup>Φ</sup>** <sup>=</sup> *ψϕ*; here, *<sup>ψ</sup>* is a randn matrix generated by a Gaussian distribution, and *ϕ* is a sparse basis matrix generated by Fourier transform. Here, *ϕ* can be given by Fourier **I***n*×*<sup>n</sup>* , and **I***n*×*<sup>n</sup>* is a unit matrix. This target signal X is sparse in Fourier space; hence, the signal X can be recovered from given {*y*, **Φ**} by CS recovery methods.

Figure 3 shows the signal recovery effect. Obviously, BPDN and SL0 do not perform well, while NSL0, L*p*-RLS and the proposed WReSL0 perform quite well. This verifies that the regularization mechanism has a good de-noising effect. Figure 4 shows the frequency spectrum of the recovered signal by the selected algorithms. The spectrum of the signal recovered by our proposed WReSL0 algorithm is almost the same as the original signal, while other algorithms fail to achieve this effect.

**Figure 3.** Signal recovery effect by BPDN, SL0, NSL0, L*p*-RLS, and weighted regularized smoothed *L*0-norm minimization (WReSL0) when noise intensity *δ<sup>N</sup>* = 0.2. (**a**) signal recovery by the BPDN algorithm; (**b**) signal recovery by the SL0 algorithm; (**c**) signal recovery by NSL0 algorithm; (**d**) signal recovery by the L*p*-RLS algorithm; (**e**) signal recovery by the WReSL0 algorithm.

**Figure 4.** Frequency spectrum analysis of the original signal and the signal recovered by BPDN, SL0, NSL0, L*p*-RLS, and WReSL0 when noise intensity *δ<sup>N</sup>* = 0.2. (**a**) original signal; (**b**) signal recovery by the BPDN algorithm; (**c**) signal recovery by the SL0 algorithm; (**d**) signal recovery by the NSL0 algorithm; (**e**) signal recovery by the L*p*-RLS algorithm; (**f**) signal recovery by the WReSL0 algorithm.

Table 2 shows the CRT of all algorithms. The *n* changes according to a given sequence [170, 220, 270, 320, 370, 420, 470, 520]. From the table, for any *n*, SL0 has the shortest computation time, followed by WReSL0, NSL0, and L*p*-RLS, and BPDN has the longest computation time. The BPDN algorithm is generally implemented by the quadratic programming method, and the computational complexity of this method is very high, thus resulting in a large increase in the overall computation time of the algorithm. Furthermore, in L*p*-RLS, the iterative process adopts the conjugate

gradient method with high complexity, while NSL0 and WReSL0 do not. Compared with NSL0, WReSL0 is more prominent in the decrease of computation time.

**Table 2.** Signal CPU running time (CRT) analysis for BPDN, SL0, NSL0, L*p*-RLS, and the proposed WReSL0 with signal length changes according to the sequence [170,220,270,320,370,420,470,520] when *δ<sup>N</sup>* = 0.2.


The performance of each algorithm under different noise intensities is shown in Figure 5. When *δ<sup>N</sup>* = 0, SL0 outperforms other algorithms, but with the increase of *δN*, the effect of SL0 becomes worse and worse. This result further illustrates that the traditional constrained sparse recovery algorithm does not have the performance of anti-noising. For BPDN, NSL0, L*p*-RLS, and WReSL0, they all applied the regularization mechanism, and they are indeed superior to SL0 in the noise case. Therefore, the proposed WReSL0 in this paper has the best de-noising performance.

**Figure 5.** NMSE analysis by BPDN, SL0, NSL0, L*p*-RLS, and WReSL0 when noise intensity *δ<sup>N</sup>* changes according to the sequence [0, 0.1, 0.2, 0.3, 0.4, 0.5].

### *4.2. Image Recovery Performance in the Noise Case*

Real images are considered to be approximately sparse under some proper basis, such as the DCT basis, DWT basis, etc. Here, we choose the DWT basis to recover these images. We compare the recovery performances based on the four real images in Figure 6: boat, Barbara, peppers, and Lena. The size of these images is 256 × 256; the compression ratio (CR; defined as *m*/*n*) is 0.5; and the noise *δ<sup>N</sup>* equals 0.01. We still choose SL0, BPDN, NSL0, and L*p*-RLS to make comparisons. For image recovery, the object of image processing is given by:

$$
\mathbf{Y} = \boldsymbol{\Phi}\boldsymbol{\lambda} + \mathbf{B} \tag{31}
$$

Here, Y, *<sup>X</sup>*, <sup>B</sup> are matrices, and among these, Y,B <sup>∈</sup> <sup>R</sup>*m*×*n*, <sup>X</sup> <sup>∈</sup> <sup>R</sup>*n*×*n*. In order to meet the basic requirements of CS, we perform the following processing:

$$\mathbf{Y}\_{i} = \Phi \mathbf{X}\_{i} + \mathbf{B}\_{i} \text{ s.t. } i = 1, 2, \dots, n. \tag{32}$$

where Y*i*, *Xi*, B*<sup>i</sup>* are the column vectors of Y, *X*, B, respectively. B*<sup>i</sup>* = *δN*Ω, Ω obeys the Gaussian distribution N (0, 1).

To perform image recovery, we valuate it by the peak signal to noise ratio (PSNR) and the structural similarity index (SSIM). PSNR is defined as:

$$\text{PSNR} = 10 \log(255^2/\text{MSE}) \tag{33}$$

where MSE <sup>=</sup> ||*<sup>x</sup>* <sup>−</sup> *<sup>x</sup>*||<sup>2</sup> <sup>2</sup>, and SSIM is defined as:

$$\text{SSIM}(p,q) = \frac{(2\mu\_p + \mu\_q + c\_1)(2\sigma\_{pq} + c\_2)}{(\mu\_p^2 + \mu\_q^2 + c\_1)(\sigma\_p^2 + \sigma\_q^2 + c\_2)}.\tag{34}$$

Among these, *μ<sup>p</sup>* is the mean of image *p*, *μ<sup>q</sup>* is the mean of image *q*, *σ<sup>p</sup>* is the variance of image *p*, *σ<sup>q</sup>* is the variance of image *q*, and *σpq* is the covariance between image *p* and image *q*. Parameters *c*<sup>1</sup> = *z*1*L* and *c*<sup>2</sup> = *z*2*L*, for which *z*<sup>1</sup> = 0.01, *z*<sup>2</sup> = 0.03, and *L* is the dynamic range of pixel values. The range of SSIM is [−1, 1], and when these two images are the same, SSIM equals one.

(**a**) Original Boat (**b**) Original Barbara (**c**) Original Peppers (**d**) Original Lena

**Figure 6.** Original images: (**a**) boat; (**b**) Barbara; (**c**) peppers; (**d**) Lena.

Figure 7 shows the recovery effect of boat and Barbara with noise intensity *δ<sup>N</sup>* = 0.01. For boat and Barbara, the recovered images by SL0 and BPDN have obvious water ripples, while recovered images by other algorithms have no such water ripples. Similarly, for peppers and Lena, the recovered images by SL0 and BPDN are blurred compared with the recovered images by other algorithms. The NSL0, L*p*-RLS, and WReSL0 algorithms are also effective at noisy image recovery. For the NSL0, L*p*-RLS, and WReSL0 algorithms, their recovery effects are very similar. In order to further analyze the advantages and disadvantages of the algorithms, we analyze the PSNR and SSIM of the images recovered by these algorithms, and the results are shown in Tables 3 and 4. By observation and analysis, L*p*-RLS performs

better than NSL0, and at the same time, WReSL0 outperforms L*p*-RLS. Hence, the WReSL0 proposed by this paper is superior to the other selected algorithms in image processing.

(**d**) Recovered Lena

**Figure 7.** Image recovery effect by the BPDN, SL0, NSL0, L*p*-RLS, and WReSL0 algorithms with noise intensity *δ<sup>N</sup>* = 0.01. In (**a**–**d**), from left to right, are: image recovered by the BPDN, SL0, NSL0, L*p*-RLS, and WReSL0 algorithms.


**Table 3.** PSNR and SSIM analysis of recovered images (boat and Barbara) by SL0, BPDN, NSL0, L*p*-RLS, and WReSL0.


**Table 4.** PSNR and SSIM analysis of recovered images (peppers and Lena) by SL0, BPDN, NSL0, L*p*-RLS, and WReSL0.

### **5. Application in Underdetermined Blind Source Separation**

The problem of UBSS stems from cocktail reception, which is shown in Figure 8. Suppose the source signal matrix **<sup>S</sup>**(*t*)=[*s*1(*t*),*s*2(*t*), ...,*sm*(*t*)]*T*, the mixed matrix (Sensors) **<sup>A</sup>** is *<sup>m</sup>* × *<sup>n</sup>* (*<sup>m</sup> <sup>n</sup>*) matrix, the Gaussian noise **G**(*t*)=[*g*1(*t*), *g*2(*t*), ..., *gm*(*t*)]*<sup>T</sup>* is generated by Gaussian distribution, and the observed mixed signal matrix **X**(*t*)=[*x*1(*t*), *x*2(*t*), ..., *xn*(*t*)]*T*; therefore, the general mathematical models of UBSS can be summarized as:

$$\mathbf{X}(t) = \mathbf{A}\mathbf{S}(t) + \mathbf{G}(t) \tag{35}$$

**Figure 8.** Schematic diagram of cocktail reception signal mixing.

In fact, each signal has *<sup>L</sup>* data collected; therefore, **<sup>X</sup>** <sup>∈</sup> <sup>R</sup>*m*×*L*, **<sup>A</sup>** <sup>∈</sup> <sup>R</sup>*m*×*<sup>n</sup>* (*<sup>m</sup> <sup>n</sup>*), **<sup>S</sup>**(*t*) <sup>∈</sup> <sup>R</sup>*n*×*L*, and **<sup>G</sup>** <sup>∈</sup> <sup>R</sup>*m*×*L*, and **<sup>G</sup>** can be represented as *<sup>δ</sup>N***<sup>W</sup>** (**<sup>W</sup>** obeys <sup>N</sup> (0, 1)). The purpose of UBSS is to use the mixed signal matrix **x**(*t*) to estimate the sof the source signal matrix **s**(*t*). In fact, this is the process of solving the underdetermined linear system of equations (ULSE). For this problem, we can use the two-step method to solve it, which is shown in Figure 9.

**Figure 9.** Schematic diagram of two-step method for UBSS.

From Figure 9, firstly, we get the mixed matrix by the clustering method and then use CS technology to separate the signal, so as to restore the original signal.

### *5.1. Process Analysis of CS Applied to UBSS*

### 5.1.1. Solving the Mixed Matrix by the Potential Function Method

In this section, we choose the potential function method to solve the mixed matrix **A**. To verify the performance of the proposed WReSL0 algorithm better, we choose four simulated signals and four real images to organize experiments in this section.

Suppose there are four source signals, which are:

$$\begin{cases} s\_1(t) = 5\sin(2\pi f\_1 t) \\\\ s\_2(t) = 5\sin(2\pi f\_2 t) \\\\ s\_3(t) = 5\sin(2\pi f\_3 t) \\\\ s\_4(t) = 5\sin(2\pi f\_4 t) \\\\ \mathbf{S} = [s\_1(t), s\_2(t), s\_3(t), s\_4(t)]^T \end{cases} \tag{36}$$

where *f*<sup>1</sup> = 310 Hz, *f*<sup>2</sup> = 210 Hz, *f*<sup>3</sup> = 110 Hz, and *f*<sup>4</sup> = 10 Hz. The length of each source signal *si* (*i* = 1, 2, 3, 4) is 1024, and the sample frequency is 1024 Hz. These four signals are shown in Figure 10.

The four source images are the classic standard test images: boat, Barbara, peppers, and Lena, which are in Figure 6.

Suppose there are two sensors that receive signals and another two sensors that receive images. Mixed matrices **A** and **B** are set as:

$$\begin{aligned} \mathbf{A} = \begin{bmatrix} A\_1 \\ A\_2 \end{bmatrix} = \begin{bmatrix} 0.9930 & 0.9941 & 0.1092 & 0.9304 \\ 0.2116 & 0.0757 & 0.9647 & 0.3837 \end{bmatrix} \\\\ \mathbf{B} = \begin{bmatrix} B\_1 \\ B\_2 \end{bmatrix} = \begin{bmatrix} 0.9354 & 0.9877 & -0.6730 & 0.1097 \\ 0.3535 & 0.07846 & 0.7396 & 0.9940 \end{bmatrix} \end{aligned} \tag{37}$$

By this mixed matrix and added Gaussian noise (*δ<sup>N</sup>* = 0.1), we can get the two mixed signals, which are shown in Figure 11, and the two mixed images, which are shown in Figure 12. Then, we can get the estimated mixed matrix **Aˆ** and **Bˆ** by clustering by the potential function method [37]. As shown in Figure 13, the potential function method can cluster well. By clustering, we get the estimated values of **A** and **B**, as follows:

$$\begin{aligned} \hat{\mathbf{A}} &= \begin{bmatrix} \hat{A}\_1 \\ \hat{A}\_2 \end{bmatrix} = \begin{bmatrix} 0.9792 & 0.9969 & 0.1097 & 0.9239 \\ 0.2028 & 0.0785 & 0.9940 & 0.3827 \end{bmatrix} \\ \mathbf{B} &= \begin{bmatrix} \hat{B}\_1 \\ \hat{B}\_2 \end{bmatrix} = \begin{bmatrix} 0.9478 & 0.9431 & -0.6483 & 0.1130 \\ 0.3476 & 0.0765 & 0.7075 & 0.9979 \end{bmatrix} \end{aligned} \tag{38}$$

(**a**) Mixed image *I*<sup>1</sup> (**b**) Mixed image *I*<sup>2</sup>

**Figure 12.** Mixed image by sensors.

(**b**) Mixed signals' polar coordinate scatter plot

(**d**) Mixed images' polar coordinate scatter plot

**Figure 13.** Clustering analysis.

By calculation, the error of solving the mixed matrix is ||**A**−**Aˆ** ||*<sup>F</sup>* ||**A**||*<sup>F</sup>* <sup>×</sup> 100% <sup>=</sup> 1.763% and ||**B**−**Bˆ** ||*<sup>F</sup>* ||**B**||*<sup>F</sup>* × 100% = 3.64%. This error range is much smaller than the classical k-means and fuzzy c-means, thus laying a foundation for the reconstruction of compressed sensing.

### 5.1.2. Using CS to Separate Source Signals

The next problem is to get **S**(*t*) from known **A**(*t*) and**X**(*t*). Here, we solve this problem by CS. The solution process is similar to the image reconstruction process. The difference is that the sparse basis used here is the Fourier basis. Then, we apply the proposed RWeSL0 algorithm to this process. First, we transform the obtained **x**(*t*) into column vectors:

$$\mathbf{x}(t) = [\mathbf{x}\_1(t), \mathbf{x}\_2(t)]^T \Rightarrow \tilde{\mathbf{x}}(t) = \begin{bmatrix} \mathbf{x}\_1(t) \\ \mathbf{x}\_2(t) \end{bmatrix} \tag{39}$$

Then, we use the Fourier (for the sparse signal) or DWT (for the image) basis for sparse representation and extend the matrix and the valuated mixed matrix to obtain the sensing matrix.

$$\begin{aligned} \mathbf{\tilde{A}} &= \mathbf{\tilde{A}} \otimes \mathbf{I}\_{L \times L} \text{ or } \mathbf{\tilde{B}} = \mathbf{\tilde{B}} \otimes \mathbf{I}\_{L \times L} \\\\ \mathbf{\tilde{Y}} &= Fourier(\mathbf{I}\_{L \times L}) / \sqrt{L}, \text{ or } \mathbf{Y} = DWT(\mathbf{I}\_{L \times L}) / \sqrt{L} \end{aligned}$$

$$\boldsymbol{\Psi} = \begin{bmatrix} \mathbf{\tilde{Y}} & \mathbf{0} & \dots & \mathbf{0} \\ \mathbf{0} & \mathbf{Y} & \dots & \dots \\ \dots & \dots & \dots & 0 \\ \mathbf{0} & \dots & \mathbf{0} & \mathbf{Y} \end{bmatrix} \tag{40}$$

$$\boldsymbol{\Phi} = \mathbf{\tilde{A}} \boldsymbol{\Psi}, \text{ or } \boldsymbol{\Phi} = \mathbf{\tilde{B}} \boldsymbol{\Psi}$$

For this equation, ⊗ denotes the Kronecker product sign, *Fourier*(·) represents the Fourier transform, and DWT represents the discrete wavelet transform. Therefore, the CS-UBSS model can be described as:

$$\begin{aligned} \dot{\mathbf{X}}(t) &= \tilde{\mathbf{A}}(t)\mathbf{S}(t) + \mathbf{G}(t) \\\\ &= \tilde{\mathbf{A}}(t)\tilde{\mathbf{P}}\Theta(t) + \mathbf{G}(t) \\\\ &= \Phi\Theta(t) + \mathbf{G}(t) \\\\ \text{or} \\ \tilde{\mathbf{X}}(t) &= \tilde{\mathbf{B}}(t)\mathbf{S}(t) + \mathbf{G}(t) \\\\ &= \tilde{\mathbf{B}}(t)\tilde{\mathbf{Y}}\Theta(t) + \mathbf{G}(t) \\\\ &= \Phi\Theta(t) + \mathbf{G}(t) \end{aligned} \tag{41}$$

where **Θ** is the Fourier transform or DWT of **S**(*t*), so **Θ** is a sparse signal. As for UBSS in the images, firstly, each image matrix needs to be transformed into a row vector, then the four row vectors form a matrix **S**(*t*). At the same time, the sparse basis in Equation (40) needs to be replaced by DWT.

Then, we can recover the source signal by CS. In summary, the above can be described as the flowchart in Figure 14.

**Figure 14.** Flowchart of UBSS by CS.

### *5.2. Performance Analysis of the WReSL0 Algorithm Applied to UBSS*

### 5.2.1. The Effect of the WReSL0 Algorithm Applied to UBSS

In this section, we evaluate the effect of the WReSL0 algorithm applied to UBSS by the separation of signals and spectrum analysis.

The effect of the separation of signals is shown in Figure 15: the source signals are well separated, and the separation signals and the original signals are very similar. Figure 16 displays the error between the original source signal and the recovered source signal. It indicates that the error between the original source signal and the recovered source signal is fairly small, and the WReSL0 algorithm can better deal with the problem of UBSS. In addition, We get the time-frequency diagram of the restored signal by short-time Fourier transform. Figure 17 is the time-frequency diagram. From this figure, we find that each signal has the same frequency as the original signal, and it also validates the rationality of the proposed algorithm for UBSS.

**Figure 16.** Separation signal error analysis.

**Figure 17.** Separation signals' frequency spectrum. Subfigures (**a**–**d**) show the frequency spectrums of separation signals *s*ˆ1, *s*ˆ2, *s*ˆ3, and *s*ˆ4.

### 5.2.2. Performance Comparisons of the Selected Algorithms

Here, we use the SL0, NSL0, and L*p*-RLS algorithms and the classical shortest path method (SPM) [38] to make a comparison in different noise cases. In order to analyze the situation of signal recovery clearly, we apply average SNR (ASNR) (for the signal) and average peak SNR (APSNR) (for the image) to evaluate. Let the original source signal be *si* and the recovered source signal be *s*ˆ*i*, so ANSR is defined as:

$$\begin{aligned} \text{ASNR} &= \frac{1}{n} \sum\_{i=1}^{n} \text{SNR}\_i \\\\ \text{SNR}\_i &= 20 \log \frac{||\beta\_i - s\_i||\_2}{||s\_i||\_{2\text{-}}} \end{aligned} \tag{42}$$

and PSNR is defined as:

APSNR = <sup>1</sup> *n n* ∑ *i*=1 PSNR*<sup>i</sup>* (43)

$$\text{PSNR}\_{i} = 10 \log \frac{255^2 \times M \times N}{||\beta\_i - s\_i||\_2}$$

where *M* and *N* are the width and height of the image.

The ASNR comparisons are shown in Table 5. From the table, we can see that ASNR attenuates sharply when *δ<sup>N</sup>* increases from 0.15–0.2. The reason is that the error of the valuated mixed matrix **A**ˆ increases obviously, which leads those CS recovery algorithms to perform poorly. In fact, from this table, our proposed RWeSL0 algorithm performs well when *δ<sup>N</sup>* is less than 0.15, and when *δ<sup>N</sup>* is greater than 0.15, the L*p*-RLS algorithm performs best, followed by our proposed RWeSL0 algorithm.

The APSNR comparisons are shown in Table 6. In this table, It is clear that APSNR is not high, and it drops greatly when *δ<sup>N</sup>* increases from 0.15–0.2. From Figure 18, we can see that these separated images seem to be enveloped in mist, which leads to a low APSNR. Therefore, we will try our best to improve this problem in the future.

In summary, the CS technique can be used in UBSS and performs well especially for the signal recovery. Our proposed WReSL0 algorithm can perform well in UBSS for the signal recoverywhen the noise is small; and regarding image recovery, we will develop this in the future.


**Table 5.** Average SNR (ASNR) analysis for separated signals by SPM, SL0, NSL0, L*p*-RLS, and the proposed WReSL0 with *δ<sup>N</sup>* changing according to sequence [0,0.1,0.15,0.18,0.2] with 100 runs.

**Table 6.** APSNR analysis for separated images by SPM, SL0, NSL0, L*p*-RLS, and the proposed WReSL0 with *δ<sup>N</sup>* changing according to the sequence [0,0.1,0.15,0.18,0.2] with 100 runs.


(**a**) Separated Boat (**b**) Separated Barbara (**c**) Separated Peppers (**d**) Separated Lena

**Figure 18.** Separated images: (**a**) boat; (**b**) Barbara; (**c**) peppers; (**d**) Lena.

### **6. Conclusions**

In this paper, we propose the WReSL0 algorithm to recover the sparse signal from given {**y**, **Φ**} in the noise case. The WReSL0 algorithm is constructed under the GD method, in which the update process of **x** in the inner loop adopts the regularization mechanism to enhance the de-noising performance. As a key part of the WReSL0 algorithm, a weighted smoothed function **W***THσ*(**x**) is proposed to promote sparsity and provide the guarantee of robust and accurate signal recovery. Furthermore, We deduced the value of *μ* and the initial value *σ*max to ensure the optimization performance of the algorithm. Performance simulation experiments on both real signals and real images show that the proposed WReSL0 algorithm performs better than the *L*<sup>1</sup> or *Lp* regularization methods and the classical *L*<sup>0</sup> regularization methods. Finally, we apply the proposed WReSL0 algorithm to solve the problem of UBSS and also make comparisons with the classical SPM, SL0, NSL0, and Lp-RLS algorithms. Experiments show that this algorithm has some advanced performance. In addition, we would also like to apply the the proposed algorithm to other CS applications such as the RPCA [39], SAR imaging [40], and other de-noising methods [41].

**Author Contributions:** All authors have made great contributions to the work. L.W., X.Y., H.Y., and J.X. conceived of and designed the experiments; X.Y. and H.Y. performed the experiments and analyzed the data; X.Y. gave insightful suggestions for the work; X.Y. and H.Y. wrote the paper.

**Funding:** This research received no external funding.

**Acknowledgments:** This paper is supported by the National Key Laboratory of Communication Anti-jamming Technology.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


c 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Detail Preserved Surface Reconstruction from Point Cloud**

### **Yang Zhou 1,2, , Shuhan Shen 1,2\* and Zhanyi Hu 1,2**


Received: 18 February 2019; Accepted: 11 March 2019; Published: 13 March 2019

**Abstract:** In this paper, we put forward a new method for surface reconstruction from image-based point clouds. In particular, we introduce a new visibility model for each line of sight to preserve scene details without decreasing the noise filtering ability. To make the proposed method suitable for point clouds with heavy noise, we introduce a new likelihood energy term to the total energy of the binary labeling problem of Delaunay tetrahedra, and we give its *s*-*t* graph implementation. Besides, we further improve the performance of the proposed method with the dense visibility technique, which helps to keep the object edge sharp. The experimental result shows that the proposed method rivalled the state-of-the-art methods in terms of accuracy and completeness, and performed better with reference to detail preservation.

**Keywords:** computer vision; 3D reconstruction; point cloud

### **1. Introduction**

Image-based scene reconstruction is a fundamental problem in Computer Vision. It has many practical applications in fields such as entertainment industry, robotics, cultural heritage digitalization and geographic systems. Image-based scene reconstruction has been studied for decades due to its low cost data acquisition and various usages. In recent years, researchers have made tremendous progress in this field. As far as small objects under controlled conditions are concerned, the performance of current scene reconstruction methods could achieve results comparable to those generated by laser scans or structured-light based methods [1,2]. However, when it comes to large scale scenes with multi-scale objects, current reconstruction methods have some problems with the completeness and accuracy, especially when concerning scene details [3].

Scene details such as small scale objects and object edges are an essential part of scene surfaces. Figure 1 shows an example of preserving scene details in reconstructing an ancient Chinese architecture. In general, cultural heritage digitalization projects, representing scene details such as the brackets in Figure 1, are among the most important tasks. The point cloud representation is often redundant and noisy, and the mesh representation is concise but it sometimes lose some information. Therefore, preserving scene details in reconstructing multi-scale scenes has been a difficult problem in surface reconstruction. The existing surface reconstruction methods [4–7] either ignore the scene details or rely on further refinement to restore them. Firstly, this is because, compared with noise, the supportive points in such part of the scene are sparse, making it difficult to distinguish true surface points from false ones. Secondly, the visibility models and associated parameters employed in existing methods are not particularly suitable for large scale ranges, where scene details are usually compromised for overall accuracy and completeness. While the first case seems to be unsolvable due to the lack of sufficient information, in this work, we focus on the second case. In particular, we extend the work of our conference paper [8], and suggest a new method with a new visibility model for surface reconstruction.

**Figure 1.** An example of the application of the proposed method in the field of cultural heritage digitalization, (**a**) point cloud; (**b**) mesh.

In many previous surface reconstruction methods [4–6,9–13], visibility information that records a 3D point is seen by the views used to help to generate accurate surface meshes. To use the visibility information, assumptions of the visibility model are made so that the space between camera center and the 3D point is free-space and the space behind the point along the line of sight is full-space. However, the points are often contaminated with noise and the full-space scales are often hard to determinate. To preserve scene details without decreasing the noise filtering ability, we propose a new visibility model with error tolerance and adaptive end weights. We also introduce a new likelihood energy representing the punishment of wrongly classifying a part of space as free-space or full-space, which helps to improve the ability of the proposed method to efficiently filter noise. Moreover, we further improve the performance of the proposed method with the dense visibility technique, which helps to keep the object edge sharp. Experimental results show that the proposed method rivals the state-of-the-art methods in terms of accuracy and completeness, and performs better with reference to detail preservation.

### **2. Related Work**

In recent years, various works have been done to advance the image-based scene reconstruction. Referring to small objects, silhouette based methods [14–18] are proposed. The silhouettes provide proper bounds for the objects, which help reduce the computing cost and yield a good model for the scene. However, good silhouettes rely on effective image segmentation, which remains a difficult task. Furthermore, silhouettes can hardly be used for large scale scenes. Volumetric methods such as space carving [19–21], level sets [22,23] and volumetric graph cut [9–11] often yield good results for small objects. However, in the case of large scale scenes, the computational and memory costs increase rapidly as scene scale grows. Consequently, this makes them unsuitable for large scale scene reconstruction. Vogiatzis et al. [9] proposed a volumetric graph cut based method to reconstruct an object by labeling voxels as inside or outside, in which a photo-consistency term is introduced to enhance the final result. Tran and Davis [10] and Lempitsky et al. [11] also exploited the same idea, the former by adding predetermined locations of possible surface as surface constraints and the latter by estimating visibility based on position and orientation of local surface patches, then optimizing on a CW-complex.

As far as outdoor scenes, uncontrollable imaging conditions and multiple scale structures make it hard to reconstruct scene surfaces. A common process of reconstructing large scale scenes is to generate the dense point cloud from calibrated images first, and then extract the scene surface. The dense point cloud can be generated through the depth fusion method [24–28] or through the feature expansion methods [29]. In depth-map fusion methods, depth-maps are usually computed

independently, and then merged into one point cloud; in feature expansion methods instead, the sparse point cloud is generated first, and then expanded to points near the seed points. Usually, the depth-map fusion based methods yield a relatively denser but noisier point cloud. Once the dense point cloud is generated, scene surface can be reconstructed through poisson surface reconstruction [30] or through graph cut based methods [4,5,12,13]. In [12], firstly the dense point cloud is used to generate Delaunay tetrahedra, then a visibility model is introduced to weight the facets in Delaunay tetrahedra, and finally the inside–outside binary labeling of tetrahedra is solved and the surface is extracted; it consists of triangles between tetrahedra with different labels. The basic assumption of the visibility model in [12] is that the space between camera center and 3D point is free-space, while the space behind 3D point is full-space. Then, this typical visibility model is promoted by a refined version in [13], namely soft-visibility, to cope with noisy point clouds. Apart from using dense point clouds, Bódis-Szomorú et al. [31] also proposed a method to produce a surface mesh by fitting the meshes reconstructed by single views and a sparse point cloud. The method is parallelizable and the quality of the final meshes is comparable to that of a state-of-the-art pipeline [32].

Besides the point-cloud based methods, large scale scene reconstruction can also be achieved through volume-based methods. Häne et al. [6,33] proposed a method for scene reconstruction and object classification. The scene space is represented by voxels and each one is given a class label. The accuracy of the classification relies on the decision tree, which is made up of labeled images. The extensive works in [34,35] reduce the heavy memory cost of the method in [33] by introducing octree structure and block scheme, respectively. Savinov et al. [36] exploited an idea similar to that in [6,33], and used full multilabel ray potential and continuously inspired anisotropic surface regularization together to yield 3D semantic models. Ummenhofer and Brox [37] proposed a method that is capable of handling a billion points. This method uses an octree structure to manage the points and then reconstruct the scene using a level-set alike method.

In this paper, we extend the work of our conference paper [8], follow the idea in [5,13] and exploit the visibility model for the line of sight in the binary labeling problem of Delaunay tetrahedra. The main contributions of the proposed method are: (1) two new visibility models for the visibility information of the points on the vertices of Delaunay tetrahedra and inside them, respectively; and (2) the dense visibility technique, which exploits the visibility information of unmerged points for better performance of preserving scene details. Experimental comparison results show that the proposed method rivals the state-of-the-art methods [5,24,27,29,38] in terms of accuracy and completeness, and performs better in detail preservation.

### **3. Visibility Models and Energies**

The pipeline of the proposed method is shown in Figure 2. The input of the proposed method is a dense point cloud generated by multi-view stereo (MVS) methods. Each point in the point cloud is attached with the visibility information recording that from which views the point is seen. Next, Delaunay tetrahedra are constructed from given input point cloud, and the scene reconstruction problem is formulated as a binary labeling problem to label tetrahedra as inside or outside. Then, an *s*-*t* graph is constructed with tetrahedra as the vertices and facets of tetrahedra as the edges. Finally, by minimizing the energy defined on the *s*-*t* graph, the vertices are separated into two parts, i.e., inside and outside, and the scene surface is extracted which consists of triangle facets lying between tetrahedra with different labels. The key issue of this process is to find a proper energy. In doing so, we introduce a new visibility model to formulate the energy of the binary labeling problem. In the following subsections, we detail the visibility models in [12,13] and our new one. Before that, we give the meaning of the symbols used in this work in Table 1 for better understanding.

**Figure 2.** Pipeline of the proposed method in 2D. From left to right are the input point cloud, Delaunay tetrahedra, the *s*-*t* graph, the energy minimization result and the final surface mesh.


**Table 1.** Symbols used in this work.

### *3.1. Existing Visibility Models*

The typical visibility model [12] assumes that, for each line of sight *v*, the space between camera center *c* and point *p* is free-space; the space behind point *p* along the line of sight is full-space. For the Delaunay tetrahedra of the point cloud, the tetrahedra intersected by segment (*c*, *p*) should be labeled as outside; the tetrahedron right behind the point *p* should be labeled as inside. For example, in Figure 3a, according to the above model, tetrahedra *T*1–*T*<sup>5</sup> should be labeled as outside, while tetrahedron *T*<sup>6</sup> should be labeled as inside.

For a single line of sight, the above label assignment is desirable. While taking all lines of sight into consideration, some surface part might not be handled properly, for example, that between *T*<sup>2</sup> and *T*<sup>3</sup> in Figure 3a. This scenario violates the label assignment principle described above. To punish the conflicts, the facets intersected by a line of sight are given weights (energy) of *W*(*lTi* , *lTj* ). Similarly, weights for punishing bad label assignments of the first tetrahedron and last tetrahedron are *D*(*lT*<sup>1</sup> ) and *D*(*lTNv*<sup>+</sup><sup>1</sup> ), respectively. Therefore, the visibility energy is the sum of the penalties of all the bad label assignments in all the lines of sight, as

$$E\_{vis}^{typical} = \sum\_{v \in \mathcal{V}} \left[ D(l\_{T\_1}) + \sum\_{i=1}^{N\_v - 1} \mathcal{W}(l\_{T\_i \prime} l\_{T\_{i+1}}) + D(l\_{T\_{N\_v + 1}}) \right] \tag{1}$$

where

$$D(l\_{T\_1}) = \begin{cases} a\_v & \text{if } l\_{T\_1} = 1 \\ 0 & \text{otherwise} \end{cases}, \\ D(l\_{T\_{\overline{\mathbb{N}}-1}}) = \begin{cases} a\_v & \text{if } l\_{T\_{\overline{\mathbb{N}}-1}} = 0 \\ 0 & \text{otherwise} \end{cases}, \\ \mathcal{W}(l\_{T\_i}, l\_{T\_{\overline{\mathbb{N}}}}) = \begin{cases} a\_v & \text{if } l\_{\overline{\mathbb{N}}} = 0 \land l\_{\overline{\mathbb{N}}} = 1 \\ 0 & \text{otherwise} \end{cases}$$

where V denotes the visibility information set, containing all lines of sight; *Nv* is the number of tetrahedra intersected by a single line of sight *v*, indexed from the camera center *c* to the point *p*; *Nv* + 1 denotes the tetrahedron behind the point *p*; *lT* is the label of tetrahedron *T*, with 1 stands for inside and 0 outside; *α<sup>v</sup>* is the weight of a line of sight *v*, which can be the photo consistency score of point *p* in current view.

**Figure 3.** Typical visibility model and soft visibility model in 2D. (**a**) Typical visibility model; and (**b**) soft visibility model for a single line of sight *v*; how to assign weight (energy) to the tetrahedron containing camera center, the end tetrahedron and the facets intersected by (*c*, *p*) or its extension are shown.

The typical visibility model as well as the energy formulation in Equation (1) described above is effective in most cases, but it has several flaws in relation to dense and noisy point clouds [13]. The surfaces reconstructed by [12] tend to be overly complex, with bumps on the surface and handles inside the model. One possible solution to these problems is the soft visibility proposed in [13], as shown in Figure 3b. In the soft visibility model, the basic assumption is similar to the previous typical visibility model. However, the edge weights are multiplied by a weight factor (<sup>1</sup> <sup>−</sup> *<sup>e</sup>*−*d*2/2*σ*<sup>2</sup> ), in which *d* represents the distance between the point *p* and the intersecting point. In addition, the end tetrahedron in the soft visibility model is shifted to a distance of 3*σ* along the line of sight.

### *3.2. Our Proposed Visibility Model*

Although the soft visibility model is effective to filter noise points and helps to yield visually smoothed models, it sometimes performs poorly in preserving details, especially in a large scene containing some relatively small scale objects (see Experimental Results). According to our observations, this happens mainly because of the improperly chosen relaxation parameter *σ* and the strong constraint imposed on the end of line of sight in the tetrahedron *kσ* from the point *p* along the line of sight. In some cases, such end tetrahedra would be free-space even though the point *p* is a true surface point.

To balance noise filtering and detail preserving, we propose a new visibility model, which is shown in Figure 4a. In our visibility model, we also use the relaxed visibility constraints in the soft visibility model in the space between the camera center *c* and the point *p*, i.e., the weight factor (<sup>1</sup> <sup>−</sup> *<sup>e</sup>*−*d*2/2*σ*<sup>2</sup> ) is kept to ensure that the final model is not overly complex. Then, we set the end of line of sight in the tetrahedron just right behind the point *p*, to avoid the wrong end in the soft visibility model in the case of small scale objects. To determine the weight of the *t*-edge of the end tetrahedron, we compare the end tetrahedra of noisy points and true surface points on datasets with quasi-truth. Figure 4b shows a typical end tetrahedron of noisy points and that of true surface points on densely sampled surfaces in 2D space. Noise points tend to appear in somewhere a bit away from true surface, which makes the end tetrahedra (triangles in 2*D*) thin and long, and true surface points are often surrounded by other true surface points, which makes their end tetrahedra flat and wide. Based on the above observations, we set a weight of *<sup>α</sup>v*(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*−*r*2/2*σ*<sup>2</sup> ) to the *t*-edge of the end tetrahedron, where *r* is the radius of the circumsphere of the end tetrahedron.

**Figure 4.** Our visibility model and end tetrahedra comparison in 2D. (**a**) In our visibility model, for a single line of sight *v*, how to assign weight (energy) to the tetrahedron containing camera center, the end tetrahedron and the facets intersected by (*c*, *p*) is shown. (**b**) From left to right: Typical end tetrahedron of noise points and that of true surface points on densely sampled surfaces.

With our new visibility model, our visibility energy is formulated as

$$E\_{\rm vis} = \sum\_{v \in \mathcal{V}} \left[ D(l\_{T\_1}) + \sum\_{i=1}^{N\_v - 1} \mathcal{W}(l\_{T\_i \prime} l\_{T\_{i+1}}) + D(l\_{T\_{\rm Nv+1}}) \right] \tag{2}$$

where

$$D(l\_{T\_1}) = \begin{cases} a\_{\upsilon} & \text{if } l\_{T\_1} = 1 \\ 0 & \text{otherwise} \end{cases}, D(l\_{T\_{N\_{\upsilon}+1}}) = \begin{cases} a\_{\upsilon}(1 - e^{-r^2/2r^2}) & \text{if } l\_{T\_{N\_{\upsilon}+1}} = 0 \\ 0 & \text{otherwise} \end{cases}$$

$$W(l\_{T\_i}, l\_{T\_j}) = \begin{cases} a\_{\upsilon}(1 - e^{-d^2/2r^2}) & \text{if } l\_{T\_i} = 0 \land l\_{T\_j} = 1 \\ 0 & \text{otherwise} \end{cases}$$

Instead of setting the tolerance parameter *σ* to a constant manually, we set *σ* adaptively within 0.5–1% of the length of line of sight in our visibility model. The underlying reason is that, typically, in depth-map fusion methods, the error bound used to filter bad correspondences while generating dense point clouds is adaptively set to 0.5–1% of the depth of the point in the current view, as in [28]. This gives each point a confidence interval along the line of sight.

### **4. Likelihood Energy for Efficient Noise Filtering**

In both the typical visibility model and our proposed one, the end of each line of sight is set in the tetrahedra right behind the point *p*. This practice sometimes could weaken the ability of noise filtering. When the surface is sampled very densely, unexpected handles could appear inside the model [13], or part of the surface could fail to be reconstructed, as shown in Figure 5. This is mainly due to the unbalanced links of *s*-edges and *t*-edges in the *s*-*t* graph, i.e., the *s*-edges are too strong to be cut as the camera centers are consistent for all lines of sight, while the *t*-edges are weak because their weights

are scattered by the varying locations of the points. The noisier the point cloud is, the greater the gap between *s*-edges and *t*-edges becomes.

**Figure 5.** Surface reconstruction without and with the likelihood energy. From left to right: (**a**) point cloud with heavy noise; (**b**) reconstructed meshes without; and (**c**) with the likelihood energy.

### *4.1. Likelihood Energy*

To solve these problems, we introduce a likelihood energy *Elike* to the total energy of the binary labeling problem. *Elike* is defined as

$$E\_{likr} = \sum\_{i=1}^{N} D(l\_{T\_i})\_\prime \text{ where } D(l\_{T\_i}) = \begin{cases} l\mathcal{U}\_{\text{out}}(T\_i) & \text{if } l\_{T\_i} = 1\\ l\mathcal{U}\_{\text{in}}(T\_i) & \text{otherwise} \end{cases} \tag{3}$$

where *N* is the total number of Delaunay tetrahedra. *Elike* measures the penalties of a wrong label assignment. For each tetrahedron, it is attached with two attributes that describe how likely it is to be outside or inside. If it is mistakenly labeled, a penalty is introduced, i.e., *Uout*(*Ti*) or *Uin*(*Ti*).

To evaluate the likelihood of the label assignment of a tetrahedron, we employ the measure *free-space support* [5], which is used to measure the emptiness of a compact space. For each line of sight *v* that intersects tetrahedra *T*, it contributes to the emptiness of *Ti*, thus increasing the probability of *T* to be outside. The free-space support *f*(*T*) of tetrahedron *T* is computed as

$$f(T) = \sum\_{\upsilon \in \mathcal{V}\_T} a\_{\upsilon\_\tau} \text{ with } \mathcal{V}\_T = \{\upsilon | \upsilon \cap T \neq \mathcal{O}\} \tag{4}$$

where V*<sup>T</sup>* is the set of lines of sight *v* that intersect with tetrahedron *T*. To evaluate *f*(*T*) to adequately describe the likelihood energy, we set *Uout*(*T*) = *λ f*(*T*) and *Uin*(*T*) = *λ*(*β* − *f*(*T*)), where *λ* is a constant used to scale the range of *f*(*T*), and *β* is a constant greater than *f*(*T*) for all tetrahedra *T*.

### *4.2. Implementation of the Likelihood Energy*

For the likelihood term *Elike*, we link two edges for each vertex *i* in *s*-*t* graph. One is from source *s* to vertex *i*, with weight *Uout*(*Ti*); the other is from vertex *i* to sink *t*, with weight *Uin*(*Ti*). However, to reduce the complexity of the graph, we cross out the *s*-edges of all vertices, and only link the vertices to sink *t* whose correspondent tetrahedra have lower *f*(*T*) than the 75th percentile of all *f*(*T*)s. Since some of the vertices in *s*-*t* graph have a heavy-weight edge linked to source *s*, we rely on the visibility constraints to label truly outside tetrahedra, instead of linking them with extra *s*-edges.

Note that the free space support threshold of 75th percentile was empirically set, as shown in Figures 6 and 7. We generated dense point clouds of four datasets, and applyiedDelaunay tetrahedralization to them. The four datasets were Fountain-P11 [3], Herz-Jesu-P8 [3], Temple-P312 [1] and scan23-P49 [2]. Then, we labeled the tetrahedra with method described in [5], and took the result as the quasi-truth. Finally, we evaluated the ratios of number of outside and inside tetrahedra in different proportions of all free-space support scores, as shown in Figure 6. In Figure 6, we can see that, generally, tetrahedra with lower *f*(*T*) (*f*(*T*) lower than the 75th percentile of all free-space support scores) had a higher probability to be truly inside. As shown in Figure 7, the true positive rates and false positive rates with different free-space support thresholds were evaluated. It is noteworthy that, when the free-space support threshold was set as the 75th percentile of each dataset, both true positive rate and false positive rate were reasonable. Therefore, we only linked those tetrahedra to sink *t* whose free-space support was lower than 75th percentile of all free-space support scores. In Figure 5, we can see that the likelihood energy is helpful for filtering noise.

**Figure 6.** Free-space support analysis. The four graphs show the ratios of number of outside and inside tetrahedra in different percentiles of all free-space support scores. The four datasets are: (**a**) Fountain-P11 [3]; (**b**) Herz-Jesu-P8 [3]; (**c**) Temple-P312 [1]; and (**d**) scan23-P49 [2].

**Figure 7.** Free-space support threshold evaluation. The true positive rates and false positive rates in the same four datasets as in Figure 6 are evaluated by setting different free-space support thresholds.

### **5. Surface Reconstruction with Energy Minimization**

With the likelihood energy and the proposed visibility model, the total energy of the binary labeling problem of the Delaunay tetrahedra is formulated as

$$E\_{total} = E\_{vis} + \lambda\_{likc} E\_{likc} + \lambda\_{qual} E\_{qual} \tag{5}$$

where *λlike* and *λqual* are two constant balancing factors; *Evis* and *Elike* are defined in Equations (2) and (3); and *Equal* is the surface quality energy introduced in [13] as

$$E\_{qual} = \sum\_{f} w\_f \tag{6}$$

where

$$w\_f = 1 - \min\{\cos(\phi), \cos(\psi)\}, \text{ if } l\_{T\_1^f} \neq l\_{T\_2^f}$$

The total energy *Etotal* defined in Equation (5) could be represented as an *s*-*t* graph and minimized by the Maxflow/Mincut algorithm [39]. In *Etotal*, the graph construction for the visibility energy *Evis* and the surface quality energy *Equal* is straightforward, and the likelihood *Elike* is implemented as described in Section 4.2. Then, the energy minimization problem is solved using the Maxflow/Mincut algorithm on the *s*-*t* graph, and the optimal label assignment of Delaunay tetrahedra is yielded. Ultimately, a triangle mesh, which consists of triangles lying between tetrahedra with different labels, is extracted. A further optional refinement can be applied as described in [7].

### **6. Dense Visibility for Edge Preservation**

Although the proposed visibility model in Figure 4 and the energy formulation in Equation (5) are carefully designed for preserving the scene details, the object edges sometimes still appear to be inaccurate concerning bumps and dents. When referring to the original depth maps, the depths of object edges are quite smooth. This scenario is shown in Figure 8. Figure 8a,b shows the object edges in 3D meshes reconstructed by the method in [5] and the proposed method described in the previous sections, as well as the original image and the corresponding depth map with a similar view point. We can easily see that the object edges in the reconstructed meshes failed to keep the smoothness as in the depth map. This could be due to the error of either point locations or visibility information which is introduced in the depth map fusion process. To circumvent such problem, instead of fusing the depths in matched cameras, we generate the dense point cloud simply by joining all of the 3D points recovered with depth maps and camera parameters. This could result in a much denser and more redundant point cloud, but it also contains more useful information for better surface reconstruction. To alleviate the memory and the computational cost, points are sampled and then Delaunay tetrahedra are constructed from the sampled point cloud. Instead of discarding those points unused for tetrahedralization and their visibility information, we apply a modified version of our visibility model to use their visibility information, as shown in Figure 9. To keep the ability to filter noise and select true surface points, we keep most of the visibility model in Figure 4 along the line of sight and only modify the part near the 3D point. The difference between the visibility models in Figure 9 and the one in Figure 4 is that, for a point *p* that lies in a tetrahedron, the end of the line of sight is set in the tetrahedron right behind the one that contains *p*, and the facet between them is punished with a weight multiplied by the same weight factor as in the soft visibility model. In this way, we keep the end of the line of sight close to the 3D point and the tolerance to noise. In Figure 8c, we show the surface meshes containing object edges reconstructed by the method in [5] and the proposed method with the dense visibility technique. We can infer from the results in Figure 8 that the dense visibility technique is helpful for preserving object edges.

**Figure 8.** Object edges in the reconstructed meshes, the original image and the depth map. From left to right: (**a**) the object edges in surface meshes reconstructed by the method in [5] (top) and our method (bottom) without the dense visibility technique; (**b**) the original image and the corresponding depth map with a similar view point; and (**c**) the object edges in surface meshes reconstructed by the method in [5] (top) and our method (bottom) with the dense visibility technique.

**Figure 9.** The modified version of our visibility model in 2D. In this visibility model, for a single line of sight *v*, how to assign weight (energy) to the tetrahedron containing camera center, the end tetrahedron and the facets intersected by (*c*, *p*) or its extension is shown.

### **7. Experimental Results**

In our experiments, the input dense point cloud was generated from images with the open source library OpenMVG (http://imagine.enpc.fr/~moulonp/openMVG/) and OpenMVS (http: //cdcseacave.github.io/openMVS/). Sparse point cloud was generated by OpenMVG, then densified with OpenMVS. Delaunay tetrahedralization was computed using CGAL (http://www.cgal.org/) library. Maxflow/Mincut algorithm [39] aws used. The proposed method was tested on public benchmark MVS dataset [2] and Tanks and Temples dataset [40].

We first tested the proposed method on MVS dataset [2]. MVS dataset [2] contains over one hundred scenes consisting of images depicting compact objects under controlled lighting conditions. Figure 10 shows the result on the MVS dataset [2]. The reference model is in the first column. From the second column to the last column, there are the models of Tola et al. [27], Furukawa and Ponce [29], Campbell et al. [24], Jancosek and Pajdla [5] and the proposed method, respectively. The final meshes given by Tola et al. [27], Furukawa and Ponce [29] and Campbell et al. [24] were generated by Poisson surface reconstruction method [30] and trimmed, which were provided in the benchmark MVS dataset [2]; the meshes of Jancosek and Pajdla [5] were generated by OpenMVS, which contains an reimplementation of the method in [5]. Figure 10 shows that the proposed method could reconstruct complex scenes as well as regular scenes, and it had great potential for preserving scene details. Figure 11 shows the accuracy and completeness of reconstructed meshes over twenty scenes, which are scans 1–6, 9–10, 15, 21, 23–24, 29, 36, 44, 61, 110, 114, 118 and 122. The detailed information of the 3D models evaluated on MVS dataset [2] is presented in Table 2. The evaluation method is described in [2], in which accuracy is measured as the distance from the MVS reconstruction to the reference model, and the completeness is measured from the reference model to the MVS reconstruction. In addition to the evaluation of the surface, we also evaluated the point clouds generated by the methods in [24,27,29] as well as OpenMVS, with the point clouds in [24,27,29] provided by MVS dataset [2]. We can see in Figure 11 that generally the point clouds achieved better scores than the surface meshes, which complies with the results in [2]. The underlying reason could be that the mesh representation discarded points that are redundant but close to the true surface, and simultaneously fixed unexpected gaps. Comparing the result within the surface meshes, the proposed method without the dense visibility technique was not outstanding, in terms of both accuracy and completeness. By applying the dense visibility technique, the proposed method achieved the best median accuracy and median completeness, and the second best mean accuracy and mean completeness. In odd rows of Figure 12 are the local point cloud and the surface meshes of Jancosek and Pajdla [5], the proposed method without the dense visibility technique and the proposed method; in the even rows of Figure 12 are the evaluation result (lower is better) of the corresponding local models through the method in [41]. We show the ability of the three methods to preserve thin objects and object edges. In some cases, the method in [5] failed completely in reconstructing them; even the complete ones were less accurate, both visually and quantitatively, than those provided by the proposed method with or without the dense visibility technique. In addition, compared with the method in [5], the proposed method showed a tremendous capability of preserving sharp object edges.

To better visualize the differences between the models generated by Jancosek and Pajdla [5] and the proposed method, we enlarged some parts of the meshes generated by the two methods. The enlarged views are shown in Figure 12. We also tested the proposed method on benchmark Tanks and Temples dataset [40]. Scenes in Tanks and Temples dataset [40] are realistic and contain plenty of objects with different scales in both outdoor conditions and indoor environments. We evaluated the proposed method as well as two other methods on four outdoor scenes (Barn, Ignatius, Truck and Courthouse) of the training set and the eight scenes of the intermediate set of Tanks and Temples dataset [40]. Figure 13 shows the result of the proposed method on four scenes of the training set of Tanks and Temples dataset [40]. Looking at Figure 13 from left to right, we can see the input images, the precision and recall of the model generated by the proposed method and the F-scores of the models generated by Colmap [38], Jancosek and Pajdla [5], and the proposed method with and

without the dense visibility technique. The evaluation method is described in [40], in which precision is measured as the distance from the MVS reconstruction to the reference model, the completeness is measured from the reference model to the MVS reconstruction, and the F-score is the harmonic mean of precision and recall with a given threshold. Table 3 presents the evaluation result of the 3D models through the method in [42]. Since the evaluation method in [40,42] takes point clouds as the input, the meshes generated by Jancosek and Pajdla [5] and the proposed method were sampled to acquire point clouds. The point clouds yielded by Colmap [38] are provided in Tanks and Temples dataset [40]. From the evaluation result of the four scenes in Figure 13 and Table 3, we concluded similarly as for the MVS dataset [2] that the proposed method outperformed the method in [5] with the dense visibility technique, while it performd slightly worse than the method in [5] without it. The proposed method performed better than the other three methods in Barn, Truck and Courthouse and rivalled Colmap [38] in Ignatius when evaluated through the method in [40], while it achieved the best result in the four scenes when evaluated through the method in [42]. Figure 14 shows the detailed views of the proposed method and the method of Jancosek and Pajdla [5] on four scenes of the training set. As in Figure 12, we also present the local models and the evaluation result using the method in [41] in odd rows and even rows, respectively. In Figure 14, Jancosek and Pajdla's method [5] inaccurately reconstructed the edge of the roof, wrongly fixed the gap between the arm and the chest of the statue, and failed to reconstruct the rearview mirror of the truck and the light stand of the lamp, while the proposed method performed well in these parts. It is noteworthy that the proposed method also had a great ability of deducing the close form of the real surface when there were no supportive points due to the occlusion, a common phenomenon in 3D reconstruction. A good example is the arm and the chest of the statue in Figure 14. We can see that the proposed method fixed the vacant part of the chest reasonably, even though it was occluded by the arm of the statue, while the method in [5] gave a less satisfactory solution. The underlying reason is that, in the proposed visibility models, the end of a line of sight lies in the tetrahedra right behind the 3D point, which is good for separating small scale objects from other objects. Table 4 and Figure 15 present the result of the proposed method on the intermediate set of Tanks and Temples dataset [40]. The detailed information of the 3D models evaluated on the intermediate set of Tanks and Temples dataset [40] is presented in Table 5. In Table 4, we can see that the result of the proposed method and that in [5] are not outstanding. However, compared to the evaluation result of the point clouds of OpenMVG + OpenMVS, the result of both the proposed method and the method in [5] achieves a substantial boost in all scenes except M60. The underlying reason could be that the two methods simultaneously filtered most noise points and fixed unexpected gaps during the surface reconstruction process. These two behaviors had opposite effects on the evaluation result, since the former one increased the F-score while the latter one decreased it. In Table 4, we can also find that generally with the dense visibility technique the proposed method outperformed the method in [5], while without the dense visibility technique it performed worse than the method in [5]. Figure 15 shows the detailed views of the proposed method and the method of Jancosek and Pajdla [5] on eight scenes of the intermediate set. Table 4 shows the F-scores of several public methods and the proposed method on eight scenes of the intermediate set. In dealing with the thin objects such as the human legs in Francis, the horse ear in Horse, the wire in M60, the barrel in Panther, the bar in Playground and the handle in Train, Jancosek and Pajdla's method [5] either failed to reconstruct it or wrongly fixed the gaps, while the proposed method performed well in these parts. Compared to the meshes of the method in [5], the proposed method also kept distinguishing object edges such as the cloth folds in Family, the horse mouth in Horse and the doorframe in Lighthouse. Therefore, we can infer that the proposed method has a strong ability to preserve the scene details.

**Figure 10.** Result of the five methods on MVS dataset [2]. GT is the reference model; Tol is Tola et al. [27]; Fur is Furukawa and Ponce [29]; Cam is Campbell et al. [24]; Jan is Jancosek and Pajdla [5]; Our is the proposed method; and ∗\_Sur is the surface generated by poisson surface reconstruction method [30].

**Figure 11.** Quantitative evaluation (lower is better) of the methods on MVS dataset [2]. The abbreviations of the methods are given in Figure 10. In addition, OMVS is OpenMVS; Our∗ is the proposed method without the dense visibility technique; and ∗\_Pts is the point cloud generated by method ∗.

**Figure 12.** Detailed views of three methods on MVS dataset [2]. From left to right are: (**a**) the point cloud generated by OpenMVS; (**b**) the mesh of Jancosek and Pajdla [5]; (**c**) the mesh of the proposed method without the dense visibility technique; and (**d**) the mesh of the proposed method. In the even rows are the evaluation result (lower is better) of the corresponding local models in the odd rows through the method in [41]. The unit is mm for all numbers.

**Figure 13.** Result of the proposed method on four scenes of the training set of Tanks and Temples dataset [40]. From left to right are: (**a**) the input images; (**b**) the precision of the model generated by the proposed method; (**c**) the recall of the model generated by the proposed method; and (**d**) the evaluation result (higher is better) of the models generated by Colmap [38] and three other methods depicted in Figure 10.

**Table 2.** Information of the 3D models evaluated on MVS dataset [2]. ∗\_Pts is the number of points in a 3D point cloud. ∗\_Vtx and ∗\_Fcs are the number of vertices and facets of a 3D mesh, respectively. The unit of all numbers is million.


**Table 3.** Evaluation result of the 3D models evaluated on the training set of Tanks and Temples dataset [40] through the method in [42]. M stands for million. Precision is expressed as a proportion 1:k, where k is the size of the scene divided by the standard error. The unit of all other numbers is mm.



**Table 3.** *Cont.*

**Figure 14.** Detailed views of three methods on four scenes of the training set of Tanks and Temples dataset [40]. From left to right are: (**a**) the point cloud generated by OpenMVS; (**b**) the mesh of Jancosek and Pajdla [5]; (**c**) the mesh of the proposed method without the dense visibility technique; and (**d**) the mesh of the proposed method. In the even rows are the evaluation result (lower is better) of the corresponding local models in the odd rows through the method in [41]. The unit is m for all numbers.

**Figure 15.** Detailed views of three methods on the intermediate set of Tanks and Temples dataset [40]. From top to bottom are Family, Francis, Horse, Lighthouse, M60, Panther, Playground and Train. From left to right are: (**a**) the point cloud generated by OpenMVS; (**b**) the mesh of Jancosek and Pajdla [5]; (**c**) the mesh of the proposed method without the dense visibility technique; and (**d**) the mesh of the proposed method.


**Table 4.** Leaderboard <sup>1</sup> of the methods and the result of the proposed method with respect to F-score on the intermediate set of Tanks and Temples dataset [40].

<sup>1</sup> https://www.tanksandtemples.org/leaderboard/.

**Table 5.** Information of the 3D models evaluated on the intermediate set of Tanks and Temples dataset [40]. ∗\_Pts is the number of points in a 3D point cloud. ∗\_Vtx and ∗\_Fcs are the number of vertices and facets of a 3D mesh, respectively. The unit of all numbers is million.


### **8. Conclusions**

In this paper, we present a new surface reconstruction method. The proposed method is designed to preserve scene details while keeping the ability to filter noise. To make the proposed method efficient to filter out noise and to select true surface points, we introduce a new visibility model with error tolerance and adaptive end weights. Along with the proposed visibility model, a new likelihood energy term is added to the total energy of the binary labeling problem to promote the robustness of the proposed method to noise. Moreover, we further improve the performance of the proposed method with the dense visibility technique, which avoids the error introduced in the point cloud generation process and provide denser visibility information. We tested the proposed method on two publicly available benchmark datasets. Experimental results on different datasets show that the proposed method rivalled the state-of-the-art methods in terms of accuracy and completeness, while it could preserve scene details such as thin objects and sharp edges. Our future work will consist in segmenting the model and adding semantic knowledge to each part of the model.

**Author Contributions:** Conceptualization, Y.Z.; Data curation, Y.Z. and S.S.; Formal analysis, Y.Z.; Funding acquisition, S.S. and Z.H.; Investigation, Y.Z.; Methodology, Y.Z. and S.S.; Project administration, S.S. and Z.H.; Resources, Y.Z. and S.S.; Software, Y.Z. and S.S.; Supervision, S.S. and Z.H.; Validation, Y.Z.; Visualization, Y.Z.; Writing—original draft, Y.Z.; and Writing—review and editing, Y.Z., S.S. and Z.H.

**Funding:** This research was funded by National Natural Science Foundation of China grant numbers 61632003, 61873265 and 61573351.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Efficient Fiducial Point Detection of ECG QRS Complex Based on Polygonal Approximation**

### **Seungmin Lee, Yoosoo Jeong, Daejin Park \* , Byoung-Ju Yun \* and Kil Houm Park \***

School of Electronics Engineering, Kyungpook National University, Daegu 41566, Korea; lsm1106@knu.ac.kr (S.L.); ysjung@ee.knu.ac.kr (Y.J.)

**\*** Correspondence: boltanut@knu.ac.kr (D.P.); bjisyun@ee.knu.ac.kr (B.-J.Y.); khpark@ee.knu.ac.kr (K.H.P.); Tel.: +82-53-950-5548 (D.P.)

Received: 20 November 2018; Accepted: 13 December 2018; Published: 19 December 2018

**Abstract:** Electrocardiogram signal analysis is based on detecting a fiducial point consisting of the onset, offset, and peak of each waveform. The accurate diagnosis of arrhythmias depends on the accuracy of fiducial point detection. Detecting the onset and offset fiducial points is ambiguous because the feature values are similar to those of the surrounding sample. To improve the accuracy of this paper's fiducial point detection, the signal is represented by a small number of vertices through a curvature-based vertex selection technique using polygonal approximation. The proposed method minimizes the number of candidate samples for fiducial point detection and emphasizes these sample's feature values to enable reliable detection. It is also sensitive to the morphological changes of various QRS complexes by generating an accumulated signal of the amplitude change rate between vertices as an auxiliary signal. To verify the superiority of the proposed algorithm, error distribution is measured through comparison with the QT-DB annotation provided by Physionet. The mean and standard deviation of the onset and the offset were stable as −4.02 ± 7.99 ms and −5.45 ± 8.04 ms, respectively. The results show that proposed method using small number of vertices is acceptable in practical applications. We also confirmed that the proposed method is effective through the clustering of the QRS complex. Experiments on the arrhythmia data of MIT-BIH ADB confirmed reliable fiducial point detection results for various types of QRS complexes.

**Keywords:** electrocardiogram; QRS complex; fiducial point; polygonal approximation; dynamic programming; QT-database; MIT-BIH arrhythmia database

### **1. Introduction**

Electrocardiogram (ECG) signals are electronically converted signals from the depolarization and repolarization of the atria and ventricle [1]. Generally, signals are composed of a P-wave, QRS complex, and T-wave, which occurred in the depolarization of the atrium and ventricle, and the repolarization of the ventricle, respectively [2,3]. Signals are determined by cardiac activity, so the signal instantly shows an arrhythmia. Signal analysis to diagnose arrhythmia has been widely used to recognize the deformation of the signal and analyze the feature value when arrhythmia occurs [4], and it is used for monitoring such as mental stress [5] and fear [6]. The study of signal analysis is subdivided according to various techniques and their purposes. In general, the signal analysis system can be divided into the noise removal step [7], fiducial point detection step [8–10], and feature value acquisition and arrhythmia classification step [11–13]. Various applications related to signal analysis include a signal monitoring system [14], heart rate measurement, signal compression [15,16], personal authentication [17,18], and so on.

The fiducial point detection step identifies the most important fiducial points of the ECG signal, which are represented by the onset, the peak, and the offset of the P-wave, the QRS complex, and the T-wave [8,19–23]. With accurate detection of these fiducial points, the width and the interval information of ECG waveforms used as feature values can be accurately measured. Therefore, the detection of accurate fiducial points is an important research field that greatly affects all subsequent ECG signal analysis.

The R-peak, which is the peak of the QRS complex, is the easiest to detect because it has the largest amplitude value. It is used not only for measuring heart rate, but also for detecting other fiducial points. Typical R-peak detection methods are based on the Pan's method, which is approximately 99% accurate. However, accurate detection methods of the fiducial points of the QRS complex other than the R-peak have been not clearly determined. Their low detection rate and inaccuracy are problematic due to signal deformation caused by various arrhythmias. These difficulties stem from the ambiguity of the reference regions that serve as a starting point.

Among various signal compression techniques, a method using polygonal approximation compresses a signal through a small number of vertices. In this case, since the onset and offset of the waveform represent the boundary between the waveform region and the baseline region, they are well-preserved as vertices.

Figure 1 illustrates the difference between the existing fiducial point detection method and the polygonal approximation method.

**Figure 1.** Motivation of the proposed work.

As shown in Figure 1a, since the existing technique uses all samples, the position of the fiducial point is ambiguous due to the samples having similar feature values nearby. Therefore, there is a high possibility of error in the threshold value, and the detection result is unreliable.

On the other hand, since the polygonal approximation uses the atomic vertices, there is a large difference in the feature values of the vertices. As shown in Figure 1b, even if the tolerance used in the approximation and the number of vertices changes, the fiducial point can be represented as a vertex of the same or similar position, so that stable detection thereof is possible.

In this paper, we propose a curvature-based vertex selection technique to solve the ambiguity of the fiducial points in the QRS complex. Our approach is roughly divided into three stages. The first stage consists of initial vertex selection using curvature-based polygonal approximation. Since the curvature value of the fiducial point is large, most of the fiducial points are represented as vertices. The second stage is an incremental vertex selection using repetitive sequential polygonal approximation, and the third stage performs an additional vertex optimization step using dynamic programming. These steps are applied for a missing case due to ambiguous curvature value.

After applying the polygonal approximation, unwanted variations of the QRS complex—such as the presence of Q-waves, S-waves, and the polarity of the waveform—remain a major problem. To mitigate the side effects of this problem, we have prepared an auxiliary signal by accumulating data from the polygonal approximation signal. This cumulative signal monotonically increases by

accumulating the absolute value of the rate of change of the amplitude between the vertices of the polygonal signal. Thus, by effectively expressing and emphasizing the feature of the fiducial point, we could effectively represent changes in the QRS complex's shape and polarity. From the cumulative signal, we analyze feature values for each vertex, such as the amplitude difference between R-peak and vertex, the time difference between reference point and vertex and the angles with neighbor vertices. Then, we determine the vertex with the largest sum of these feature values as a fiducial point.

This paper is organized as follows. Section 2 briefly reviews the ECG signal composition and explains why detecting the fiducial point is relatively difficult. Section 3 introduces the curvature-based vertex selection technique proposed in this paper and shows the expected benefits when applied to ECG signals. Section 4 details our algorithm for generating the cumulative signal from the polygonal approximated signal and detecting the fiducial point therefrom. In Section 5, the performance of the proposed algorithm is verified through experiments on QT-DB [24] and MIT-BIH ADB [25], and Section 6 concludes the paper.

### **2. Composition of ECG Signal**

The ECG signal, which consists of the P-wave, QRS complex, and T-wave, includes the corresponding onset, offset, and peak points, which are referred as fiducial points. Figure 2 shows the fiducial points and feature values of the P-wave, QRS complex, and T-wave of the ECG signal.

**Figure 2.** Composition of ECG signal.

As shown in Figure 2, the ECG signal is divided into a waveform region in which the amplitude is changed by depolarization and repolarization, and a baseline region in which no amplitude change occurs. The boundary point is a fiducial point of each waveform. However, the actual input signal contains various noise, not ideal forms, as shown in Figure 2. Typical types of noise are as follows [26–29].


1 and 2 are typical high- and low-frequency ECG signal noises, respectively. Since noise complicates the baseline, it is difficult to estimate baseline and thus detect the boundary with an ambiguous waveform. Most of study use the bandpass filter for suppress noises and in some case, it uses a notch filter [30] for aiming to suppress the power line interference, such as high frequency of 50 Hz or 60 Hz.

Figure 3 shows the result of applying a high-pass and low-pass filter to suppress baseline deviation and power line interference.

**Figure 3.** Noises of ECG signal and filtering results.

Peaks, such as the R-peak, can be easily detected because the amplitude has a local maximum or minimum and the rate of amplitude change is large enough. In contrast, the boundary between the baseline and the waveform is still ambiguous, even when filtering is applied. This is because the amplitude change occurs slowly, and the features are similar to the surrounding samples. In this paper, we propose an effective fiducial point detection technique by emphasizing ambiguous fiducial points based on polygonal approximation.

### **3. Polygonal Approximation of ECG Signal**

Signal approximation techniques have been widely studied by using polygonal approximation, such as sequential polygonal and cyclic polygonal approximation, and polynomial approximation, such as B-spline. However, these techniques are problematic, since they select too many vertices and errors are not minimized.

In the curvature-based vertex selection technique [31], curvature-based polygonal approximation [32], sequential polygonal approximation [33], and dynamic programming [34] are proposed as methods for minimizing vertices and resultant errors. This algorithm selects an initial vertex using a curvature-based polygonal approximation. A fiducial point with a large curvature is selected as a vertex. However, there is a problem when a fiducial point having an ambiguous curvature value is not selected as an initial vertex. To solve this problem, the sequential polygonal approximation method is adopted to select additional vertices, and the dynamic programming technique can optimize the position of the additional vertices by minimizing errors.

The algorithm flow of the curvature-based vertex selection method for the input ECG signal (*S*) is summarized as follows.


$$V^I = \{v\_1^I, v\_{2'}^I \cdots, v\_S^I\} \tag{1}$$

3. We apply the sequential polygonal approximation method to the interval between each initial vertex to select additional vertices. Equation (2) represents a set of *NVi* − 1 additional vertices between the *i*-th initial vertex and the *i* + 1-th initial vertex, and both end vertices coincide with the two initial vertices.

$$\begin{aligned} V\_i &= \{ v\_{i,0\prime} \cdot \cdot \cdot \cdot v\_{i,N\_{V\_i}} \} \\ v\_{i,0} &= v\_{i\prime}^I \qquad v\_{i,N\_{V\_i}} = v\_{i+1}^I \end{aligned} \tag{2}$$

4. Dynamic programming is applied to the additional vertices to optimize their position. Equation (3) is a set of corrected vertices for the additional vertex set *Vi*.

$$\begin{aligned} V\_i^{Opt} &= \{ v\_{i,0}^{Opt}, \cdot, \cdot, v\_{i,N\_{V\_i}}^{Opt} \} \\ v\_{i,0}^{Opt} &= v\_{i,0} = v\_{i'}^I \qquad v\_{i,N\_{V\_i}}^{Opt} = v\_{i,N\_{V\_i}} = v\_{i+1}^I \end{aligned} \tag{3}$$

5. Repeat steps 2–4 to proceed with polygonal approximation for the entire input signal. Equation (4) represents the set of *NV* vertices as the result of vertex selection.

$$V = \{v\_1, \dots, v\_{N\_V}\}, v\_i = (v\_{x\_i}, v\_{y\_i}) \tag{4}$$

Figure 4 shows the result of each step of the polygonal approximation.

**Figure 4.** Additional vertex calibration results using dynamic programming.

In general, when the curvature-based polygonal approximation is applied to a pole with a large curvature value, it appears as an initial vertex, as shown in Figure 4a. However, with a smooth transition of the amplitude value near the fiducial point, such as the onset of the QRS complex and the offset of the P-wave in Figure 4b, there is a side effect wherein the fiducial point is not selected as the initial vertex. By applying additional vertex selection and correction to solve this problem, the fiducial points are efficiently represented by the vertex as shown in region A, and the similarity between the original signal and the approximated signal is preserved, as shown in B and C.

### **4. Fiducial Point Detection Based on Polygonal Approximation**

The curvature-based vertex selection technique represents the ECG signal as a small number of vertices, and then detects the onset and offset of the QRS complex by analyzing the characteristic values of each vertex. However, it is not easy to express the characteristic value from the fiducial point because the QRS complex has various shapes based on its polarity, as well as the presence of the Qand S-peaks.

To resolve the difficulty of extracting features from the QRS complex's ambiguous shape, robust fiducial point detection using various auxiliary signals has been proposed. In Pan's method, the R-peak detection is assisted by an auxiliary signal generated by the derivative of the signal and an average filter. In Manriquez's method, the fiducial point is detected by a threshold value of the auxiliary signal generated from the Hilbert transform. This paper is also based on the auxiliary signal, for which we propose the cumulative signal of the polygonal approximation to preserve morphological features of the vertex that represent the fiducial point.

### *4.1. Generate the Cumulative Signal*

To acquire the cumulative signal, we first obtain the amplitude difference (*VD*) for the vertex, as shown in Equation (4).

$$\begin{aligned} V^D &= \{v\_1^D, \dots, v\_{N\_V}^D\}, & v\_i^D &= (v\_{x\_i}^D v\_{x\_i}^D) \\ v\_{x\_i}^D &= v\_{x\_{i'}} & v\_{y\_i}^D &= v\_{y\_i} - v\_{y\_{i-1'}} & v\_{y\_1}^D &= 0 \end{aligned} \tag{5}$$

With the absolute value of the amplitude difference obtained using Equation (5), that value is accumulated as shown in Equation (6) to generate the cumulative signal.

$$\begin{aligned} \boldsymbol{V}^{D'} &= \{\boldsymbol{v}^{D'}\_1, \dots, \boldsymbol{v}^{D'}\_{N\_V}\}, & \boldsymbol{v}^{D'}\_i &= \left(\boldsymbol{v}^{D'}\_{x\_i}, \boldsymbol{v}^{D'}\_{x\_i}\right) \\ \boldsymbol{v}^{D'}\_{x\_i} &= \boldsymbol{v}^{D}\_{x\_i \prime}, & \boldsymbol{v}^{D'}\_{y\_i} &= \sum\_{k=1}^i |\boldsymbol{v}^{D}\_{y\_k}| \end{aligned} \tag{6}$$

This simplifies the signal as monotonically increased, even if the QRS complex appears as a downward wave or includes Q-peaks and S-peaks that appear as downward or upward waves. The vertex corresponding to the fiducial point also maintains the feature of dividing the baseline and waveform regions.

Figure 5 shows the cumulative signal results for various shapes of the polygonal approximation signals.

**Figure 5.** Comparison of cumulative signal results according to waveform type.

As shown in Figure 5 even if the signal includes Q- or S-peaks, or the QRS complex shows a downward wave, it can be expressed as a cumulative signal of similar shape. In this case, the features are also similar, so that effective fiducial point detection is possible.

### *4.2. Algorithm of Fiducial Point Detection*

The feature value of each vertex of the accumulated signal is analyzed to determine the fiducial point. In this paper, we propose three types of features for each vertex to determine the fiducial point: the amplitude difference between R-peak and vertex, the time difference between the reference point and vertex, and the angles with neighbor vertices.

### 4.2.1. Amplitude Difference between R-Peak and Vertex

Figure 6, which is magnified from red-dotted box in Figure 5b, show amplitude difference between R-peak and vertex from this cumulative signal.

**Figure 6.** The amplitude difference between R-peak and the vertex in the cumulative signal.

The onset and offset of the QRS complex are the boundary points between it and the baseline region. The amplitude difference between the R-peak and the vertex in the accumulated signal is close to the maximum value near the fiducial point. Therefore, the Q-peak with a largest amplitude difference in Figure 5a become a smallest amplitude difference in the cumulative signal, which makes it easier to determine the fiducial point.

Equation (7) represents the feature value obtained by using the amplitude difference between the R-peak and the vertex.

$$\begin{aligned} A\_{i\_L} &= v\_{y\_i}^{D'} - v\_{y\_1}^{D'} \\ A\_{i\_R} &= v\_{y\_{N\_V}}^{D'} - v\_{y\_i}^{D'} \end{aligned} \tag{7}$$

*AiL* denotes an amplitude difference between the previous R-peak and vertex, which is used to detect the offset. Similarly, *AiR* is used to detect the onset.

### 4.2.2. Time Difference between Reference Point and Vertex

Figure 7 shows the time difference between the vertex which is 0.3 s away from the R-peak, and the reference point.

**Figure 7.** The time difference between the reference point and vertex in the cumulative signal.

In this paper, we suggest the time difference as a second feature value for excluding the case of detecting the onset of the P-wave. Generally, the normal width of the QRS complex is about 0.08 to 0.12 s. We use a reference point based on the point 0.3 s away from the R-peak for estimating the time difference and consider that the larger the time difference is, the more likely it is to be the fiducial point.

If the ventricular arrhythmia occurred and the width of the QRS complex is increased, fiducial point may have lower time difference feature value when there is a vertex at notch in QRS complex. However, this problem can be easily solved from the amplitude difference between previous R-peak and vertex, since the amplitude of notch and fiducial point is similar to previous R-peak and baseline, respectively.

Equation (8) represents the feature value obtained by using the time difference between the references and vertex.

$$\begin{aligned} T\_{\dot{i}\_L} &= \left(\upsilon\_{\mathbf{x}\_1}^{D'} + 0.3 \times F\right) - \upsilon\_{\mathbf{x}\_i}^{D'} \\ T\_{\dot{i}\_R} &= \upsilon\_{\mathbf{x}\_i}^{D'} - \left(\upsilon\_{\mathbf{x}\_{N\_V}}^{D'} - 0.3 \times F\right) \end{aligned} \tag{8}$$

*F* denotes a sampling frequency and *TiL* denotes a time difference between the previous R-peak and vertex, which is used to detect the offset. Similarly, *TiR* is used to detect the onset.

### 4.2.3. Angles with Neighbor Vertices

Since the fiducial point is boundary between the waveform region and the baseline region, most significant feature of vertices of fiducial point is angle between the horizontal line and straight line connecting the neighbor vertex.

Figure 8 shows the angles of the vertices corresponding to the fiducial points with the left and right vertices in the cumulative signal.

In the case of the onset, the angle with the left vertex is close to 0 degrees, and the angle with the right vertex is close to 90 degrees. In the case of the offset, it is reversed. *θiL* and *θiR* mean the left and right angles of the *i*-th vertex, respectively.

**Figure 8.** The angle of fiducial point in cumulative signal.

### 4.2.4. Detecting the Fiducial Point

Based on the feature values of the fiducial point, we can calculate the feature value of each vertex in the searching interval. Our approach provides a method to select the point with the highest probabilities of being the fiducial point by summarizing all feature values. Equations (9) and (10) represent equations used to detect the onset and offset of the QRS complex, respectively.

$$w\_{Q\_{\rm av}} = \arg\max \{ \omega\_A(A\_{i\_R}) + \omega\_T(T\_{i\_R}) + \omega\_{\theta\_C}(\theta\_{i\_L}) + \omega\_{\theta\_S}(\theta\_{i\_R}) \},\tag{9}$$

$$v\_{S\_{eff}} = \operatorname\*{argmax} \{ \omega\_A(A\_{i\_L}) + \omega\_T(T\_{i\_L}) + \omega\_{\theta\_S}(\theta\_{i\_L}) + \omega\_{\theta\_C}(\theta\_{i\_R}) \},\tag{10}$$

where,

$$\begin{aligned} \omega\_A(A\_i) &= \frac{A\_i}{\max(A\_i)'} & \omega\_T(T\_i) &= \frac{T\_i}{0.3 \times F}, \\ \omega\_{\theta\_C}(\theta\_i) &= \sqrt{\cos \theta\_i} & \omega\_{\theta\_S}(\theta\_i) &= \sin^2 \theta\_i \end{aligned}$$

*ωA*, *ωT*, and *ωθ* are weight functions for normalizing each feature value. *ω<sup>A</sup>* uses the maximum value in the search interval because sensitive amplitude changes in the QRS complex may be affected according to each heartbeat. *ω<sup>T</sup>* uses 0.3 s, which means the reference, and the *ωθ<sup>C</sup>* and *ωθ<sup>S</sup>* are used to have higher feature value when the baseline and waveform angles are close to 0 and 90 degree, respectively.

In the case of the baseline direction angle, the square root is added to consider that the vertex corresponding to the fiducial point may have a high value of about 30 to 40 degrees due to noise or signal distortion. On the other hand, since the waveform angle has a low possibility of distortion because of the large amplitude change, a square is added to have a low feature value for the vertices except the fiducial point.

### **5. Experiment and Analysis of Results**

Figure 9 is the flowchart of our proposed algorithm.

In the preprocessing step, noise suppression and R-peak detection are performed. Breathing and muscle movements cause the low-frequency noise, and power noise of 30 Hz or 60 Hz causes the high-frequency noise. In this paper, these noises are suppressed using a Butterworth bandpass filter of 1–25 Hz and the most widely known Pan's method is applied to R-peak detection. The QT-DB and MIT-BIH ADB provided by Physionet are used to evaluate the performance of the proposed algorithm.

**Figure 9.** Algorithm flowchart.

### *5.1. Experiment in QT-DB*

The QT-DB contains a total of 105 fifteen-minute excerpts of two-channel ECGs, which were carefully selected to avoid significant baseline wander or other artifacts. Within each record, around 30 numbers of beats were manually annotated by cardiologists, who can identify the onset, peak, and offset of the QRS complex. The proposed algorithm works on a single-channel signal, while cardiologist manually recorded one annotation considering both channel of signal simultaneously. Therefore, to compare our approach with the manual annotations on the QT-DB for each of the two single-channels, it is reasonable to choose the annotation result of the channel with less error [10]. After all the errors are obtained, the mean of total error *μ* and the standard deviation of total error *σ* are computed by averaging the intrarecording mean and standard deviation of each set of data. The standard deviation of the total error is used to measure the criterion for the algorithm's stability. In the CSE working party [35], tolerances for standard deviation of the error for the onset and offset of QRS complex are suggested as 6.5 ms and 11.6 ms, respectively.

We summarized the experimental results by our proposed method in Figure 10 according to the types of DB constituting the QT-DB.


**Figure 10.** Distribution of standard deviation for data.

Most of the data is regarded as having satisfied the tolerance or shown on the detection result similar to the tolerance. Especially in the case of MIT-BIH normal sinus rhythm DB (represented as a red marker), it can be confirmed that the standard deviation of the total error satisfies tolerance. On the other hand, MIT-BIH arrhythmia DB and sudden-death patient DB, represented by a black marker, contains various arrhythmia heartbeats and caused large errors compared to another DB.

Table 1 shows the performance of the proposed algorithm compared to the existing detection algorithm.


**Table 1.** QRS segmentation performance comparison in the QT-DB.

As with other algorithms, the standard deviation of error of the onset is out of tolerance. On the other hand, the offset satisfies the tolerance and confirmed a lower error than other algorithms.

To analyze experimental results in detail, we tried to cluster the experimental results for each data set. The experimental results of some data according to the types of QRS complex are shown in Figure 11.

**Figure 11.** Clustering results.

Clustering was performed around the R-peak, and the position of the R-peak was shown as the origin. As shown in Figure 11a,b, not only the detection results of the normal type of QRS complex are stable, but also the fiducial point detection is well-performed for various types of QRS complex, such as absence of Q-peak, widen QRS complex and downward QRS complex (see Figure 11c–f).

### *5.2. Experiment in MIT-BIH ADB*

Since MIT-BIH ADB only annotate the arrhythmia type of each heartbeat, statistical analysis for fiducial point detection is hard. We use the mean length of QR section, which start from the onset to R-peak, and RS section, which start from R-peak to the offset, as reference for onset and offset detection, respectively. According to the type of heartbeat, we separate into normal heartbeats and abnormal heartbeats, and calculate the mean and standard deviation of each type of heartbeats, respectively.

However, in the case of arrhythmia, there are various types of abnormal heartbeats. Most occur with premature atrial contraction (PAC; annotated as 'A') or premature ventricular contraction (PVC; annotated as 'V'), with the remainder very rare. For example, supraventricular premature or ectopic beat (SVP; annotated as 'S') only appears twice in record 208 over the total 107,265 heartbeats in MIT-BIH ADB. Thus, the test for abnormal heartbeat only proceeded for one of the most typical abnormal heartbeat types for records with more than 30 PACs or PVCs.

Figure 12 shows the distribution of standard deviation of fiducial points in MIT-BIH ADB.

**Figure 12.** Distribution of standard deviation for MIT-BIH ADB.

Figure 12a is distribution of fiducial points for normal heartbeats (annotated as 'N') or other representable type of heartbeat, such as left bundle branch block (LBBB; annotated as 'L'), right bundle branch block (RBBB; annotated as 'R') or pacemaker (PM; annotated as '/'). Considering that the length of each record is 30 min, which is 60 times longer than the length of QT-DB used in Section 5.1, and most signals are unstable with various arrhythmia, it can be confirmed that the detection result of fiducial point is excellent.

Generally, since the PAC does not affect the QRS complex, the stability of the fiducial point detection is high as shown in Figure 12b. On the other hand, since PVC deforms QRS complex into various shapes, the standard deviation of width is increased.

Table 2 represents the detailed distribution of stable data in Figure 12b. In the case of PVC, the results for uniform PVC are shown.


**Table 2.** Detailed results for stable data in Figure 12.

Figure 13 shows the results of the proposed algorithm for the part of the MIT-BIH ADB data in which normal heartbeat and PVC occur consecutively.

**Figure 13.** The result of fiducial point detection for arrhythmia data.

The proposed algorithm detected the fiducial points regardless of whether Q- or S-waves are present or not in the normal heartbeat. In addition, we could reliably detect the fiducial points for the various types of QRS complex with arrhythmia.

From this experiment, we can confirm that the proposed polygonal approximation method can effectively detect fiducial points compared with other algorithms. In addition, it obtains stable results for various types of QRS including arrhythmia. In some cases, detection error is high due to the acquisition of signals such as the V2 and V5 channel, rather than signal acquisition by the modified lead II (MLII) channel, which emphasizes the QRS complex. The distortion in the preprocessing process causes an error. Therefore, we consider that the signal acquisition clearly reflecting the shape of the QRS complex—such as the MLII channel—or the stable noise suppression technique can improve the performance of the proposed algorithm.

### **6. Conclusions**

In this paper, the proposed algorithm focuses on mitigating the ambiguity of the fiducial point by polygonal approximation to extract the feature points of ECG signals, motivated from the characteristic that the polygonal approximation preserves the fiducial point as a vertex. Therefore, our approach resulted in better performance compared to other signal compression techniques. In addition, we propose an effective auxiliary signal for stable detection results in various types of QRS complex using these features. Experimental results show that the proposed auxiliary signal-based technique enables stable detection for various applications in real ECG databases equipped with QT-DB and MIT-BIH ADB.

In future research, the proposed method is required not only to expand the detection of P- and T-waves, but also to apply post-processing to improve the existing technique. This method can also be extended to study adaptive signal compression technique according to importance of intervals by reapplying fiducial point detection results. The trade-off between the accuracy and energy consumption could be achieved by adjusting the fiducial points compression ratio, so the extremely long-time ECG monitoring services are capable in the battery-operated wearable systems.

**Author Contributions:** S.L. designed the entire core architecture and performed the hardware/software implementation and experiments; Y.J. performed the simulation and measurements; D.P. has his responsibility in writing entire paper as the main corresponding author; B.-J.Y. presented the initial concept and has his role as co-corresponding author; K.H.P. is the principle investigator serving co-corresponding author.

**Funding:** This study was supported by the BK21 Plus project funded by the Ministry of Education, Korea (21A20131600011). This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2014R1A6A3A04059410, NRF-2018R1A6A1A03025109, 2018R1A6A3A01011035). The APC was funded by NRF-2018R1A6A1A03025109.

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

### **Abbreviations**

The following abbreviations are used in this manuscript:


### **References**


c 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Signal Amplification Gains of Compressive Sampling for Photocurrent Response Mapping of Optoelectronic Devices**

### **George Koutsourakis 1,\* , James C. Blakesley <sup>1</sup> and Fernando A. Castro 1,2**


### Received: 17 May 2019; Accepted: 25 June 2019; Published: 28 June 2019

**Abstract:** Spatial characterisation methods for photodetectors and other optoelectronic devices are necessary for determining local performance, as well as detecting local defects and the non-uniformities of devices. Light beam induced current measurements provide local performance information about devices at their actual operating conditions. Compressed sensing current mapping offers additional specific advantages, such as high speed without the use of complicated experimental layouts or lock-in amplifiers. In this work, the signal amplification advantages of compressed sensing current mapping are presented. It is demonstrated that the sparsity of the patterns used for compressive sampling can be controlled to achieve significant signal amplification of at least two orders of magnitude, while maintaining or increasing the accuracy of measurements. Accurate measurements can be acquired even when a point-by-point scan yields high noise levels, which distort the accuracy of measurements. Pixel-by-pixel comparisons of photocurrent maps are realised using different sensing matrices and reconstruction algorithms for different samples. The results additionally demonstrate that such an optical system would be ideal for investigating compressed sensing procedures for other optical measurement applications, where experimental noise is included.

**Keywords:** non-destructive testing; current mapping; digital micromirror device; compressed sensing

### **1. Introduction**

The non-uniformities of material structure and local defects can have an influence on the overall performance of optoelectronic devices, such as solar cells and photodiodes. Therefore, it is important to develop methods that provide spatially resolved information on the defects and inhomogeneities of such semiconductor devices. Light/laser beam induced current (LBIC) methods have been established for the spatial characterisation of solar cells [1], photodiodes [2,3], and other sensors and *p*–*n* junction devices [4,5]. For the realisation of current mapping, a light beam scans the device being tested, and the induced current is measured for every point. A variety of different system approaches have been proposed, making the LBIC measurement systems able to deliver spatial maps of electrical properties [6], local reflectivity [7], performance parameters [8], and material properties of optoelectronic devices. Recent implementations utilise multiple laser wavelengths that enable measurements on a larger range of samples and for different energy ranges [1,9].

Although useful and sometimes necessary, photoresponse mapping measurements are usually time-consuming, since a small spot size has to scan the entire active area of the device for a total area current map, which means the smaller the spot size, the lengthier the measurements. Focusing the laser beam on a small spot often requires elaborate optical elements and accurate alignment. A very frequent solution is to use a microscope objective lens to achieve a spot size of several micrometres [10]. The point-by-point scan is realised by an *x*–*y* stage, which means that there is always a time delay

from moving from one point to the next. Continuous acquisition methods have been reported in order to accelerate this process [11]; nevertheless, this can result in small distortions of the current. The alternative option to maximise scan speed is to use piezo-electric mirror systems to guide the beam on the sample [12,13]. Spot sizes of several micrometres also mean very weak signals. For this reason, lock-in techniques for accurate current readings were introduced, even in very early systems [14], and have been used in almost every LBIC system implementation ever since. Combining all of the above features into one system is not trivial, and LBIC systems can become very complicated to realise.

The first attempt to utilise digital light processing (DLP) for current mapping was by a fast tomographic current mapping method for photovoltaic (PV), based on a digital micromirror device (DMD) for implementing the scan [15]. DLP devices utilise a DMD to create light projections [16]. Photocathode quantum efficiency mapping using a digital micromirror device has been reported, where the DMD implements the laser raster scan [17]. A DLP projector has been utilised for low-resolution spatial uniformity characterisation of solar cells [18]. A DMD-based system has also been reported for fast spectral response measurements of PV devices [19], where additional frequency modulation for each wavelength band has been reported, in order to accelerate measurements [20]. High-frequency light modulation of more than 40 GHz has also been introduced recently, with an Si light emitter embedded in a *p*-channel, metal oxide, semiconductor field effect transistor (PMOSFET) structure [21].

Using a DMD to apply compressed sensing (CS) current mapping of PV devices has been demonstrated in recent work [22–24]. CS current mapping has also been demonstrated by utilising an LCD (liquid crystal display) monitor to project the necessary patterns for compressive sampling [25]. The CS current mapping method is based on the CS sampling theory [26,27]. According to this theory, it is possible to reconstruct a signal from highly incomplete or inaccurate information. Compression of signals is something very common in everyday life. For instance, in JPEG image compression, most of the signal information is thrown away at the transform compression stage. Only the necessary elements for describing the image in the transform domain are kept (*K* elements). The image is reconstructed using these very few *K* elements, which provide a sparse representation of the image. The aim of CS imaging is to directly acquire the *K* coefficients necessary for an almost exact reconstruction of a signal. This is achieved by only acquiring *M* < *N* measurements for capturing an *N* pixel image, where *K* < *M*. There are a large number of compressive sampling applications, such as CS Magnetic resonance imaging (MRI) [28], the single-pixel camera [29], CS radar imaging [30], CS confocal microscopy [31], and many more.

In previous work, we have presented the CS current mapping methodology and a design for a CS current mapping measurement system for solar cells [22,23]. In this work, the signal amplification aspects of the sampling process and technical approaches for optimised sampling are presented. Although the performance of different CS aspects (algorithms, transforms, and matrices) can be investigated by simulations, aberrations of compressive sampling due to instrumentation and optics only show in experimental investigations, such as the one presented in this work. The significant signal amplification gains of CS current mapping for optoelectronic devices are demonstrated and discussed. The optimum ways to achieve such amplification for the measured signal and the impact of sensing matrix sparsity (defined later on) on the accuracy of measurements are investigated for the first time for such an application. Three types of devices are used to illustrate that the choice of sensing matrix sparsity is dependent of sample and measurement instrumentation. The robustness of CS current mapping against long-term measurement noise is studied, and a pixel-by-pixel comparison of compressive and point-by-point sampling for current mapping is realised. Different measurement settings and samples are tested using the DMD optical system. This comparison aids in determining the most suitable occasions in which each sampling method should be realised, and presents a realistic performance evaluation of compressive sampling for this specific application. In addition, it is demonstrated that this optical setup is ideal for realistic experimental comparisons of reconstruction algorithms for optical measurement applications of compressed sensing.

### **2. Methodology**

### *2.1. Experimental Layout*

The optical current mapping system used in this work is based on a DMD kit and is presented in Figure 1. A single mode fibre-coupled laser source of 40 mW at a 637 nm wavelength is used. The light output of the fibre is collimated such that the beam overfills the DMD micro-mirror area. The DMD is a V-7000 module, consisting of a 1024 × 768 pixel micromirror array, each micromirror having a pixel size of 13.7 × 13.7 μm. A spatial filter is used to reject the diffracted and non-collimated components of the beam. Finally, a mirror is used for guiding the beam onto the sample, which is placed horizontally on a *z*-stage platform. A National Instruments PXIe-4139 source measure unit (SMU) is used for measuring the current for both cases of sampling (raster scanning and patterns). The sample is placed at the focal plane of the last lens, so that the scanning spot or the patterns are projected onto the sample. In order to apply a compressive or a point-by-point scan, a number of micromirrors are grouped together to form one pixel, and the number of grouped micromirrors depends on the selected optical resolution. The spot shape is square. The sampling methods are presented in Figure 2. As can be seen in the picture of the DMD on the right of Figure 2, not all of the active area of the DMD (1024 × 768 pixels) is used. A square 700 × 700 pixel area is used to project the patterns, in order to create a square projection. Groups of 7 × 7 micromirrors are binned together, creating projections of 100 × 100 pixels. This results in a 100 × 100 resolution of the final current maps. The sampling rate that can be achieved is 30 points or patterns per second and this sampling rate is used for all silicon samples of this work. For the organic device measured, a slower sampling rate was selected (5 samples/s), due to the slower response of the specific organic photovoltaic device [32].

**Figure 1.** On the left, schematic of the optical system used for compressed sensing (CS) current mapping. On the right, a photo of the system.

**Figure 2.** Schematics of the two different sampling modes used, and photos of how each is implemented on the digital micromirror device (DMD): on the left, point-by-point sampling; on the right, compressive sampling.

### *2.2. Compressed Sensing Current Mapping*

For the application of compressive sampling, a series of binary patterns are projected onto the sample's area to be measured, and the photocurrent response of the sample is measured for each pattern. The patterns are generated by the DMD, assigning pixels (binned groups of micromirrors) as either "on" or "off", illuminating or shading different points of the sample, as can be seen in Figure 2. Similar to JPEG image compression, the sequence of patterns measures and compresses the necessary

information, in order to successfully reconstruct the photocurrent response map. This is a standard procedure for optical CS imaging systems [33] that is analytically described for CS current mapping of photovoltaic (PV) devices in [34]. Compared to a point-by-point scan, fewer measurements are required in order to produce a current map of a sample.

In summary, for the application of CS current mapping, a series of binary patterns Φ = {ϕ*m*} *M m*=1 are projected onto the sample, in order to acquire a compressed representation of the signal **x**, which has *N* elements, using *M* < *N* linear measurements. Each row of the sensing matrix Φ is a binary pattern expressed as a vector, which makes Φ an *N* × *M* matrix. The current response of the PV device is measured for each pattern, populating the measurement vector **y**. An underdetermined problem **x**: **y** = Φ**x** is created, since **y** has fewer elements than **x**. Random binary matrices are used in this work to produce the sensing matrix, as it has been shown that they possess the necessary properties needed for compressive sampling [35]. The discrete cosine transform (DCT) is applied as a basis to provide the sparse representation of the signal. Two different algorithms are used in this work for solving this underdetermined problem and reconstructing the current map. The first is the 1 norm minimisation basis pursuit algorithm, included in the 1 magic toolbox in MatLab developed by Candès and Romberg [36]. A second algorithm used is the orthogonal matching pursuit (OMP) algorithm [37]. Using one of the above algorithms, the underdetermined problem is solved, and the current map is reconstructed. 

Although it is not within the scope of this work to investigate different reconstruction algorithms, the right choice of algorithm can be crucial for the successful reconstruction of the final current map. Nevertheless, given a specific algorithm, the choice of sensing matrix sparsity does not significantly affect the reconstructed image, as will be demonstrated in this work. This is shown by acquiring similar current mapping results when using different sensing matrices, for each of the two algorithms. In addition, it is demonstrated that this simple optical experimental setup is ideal for comparing CS reconstruction algorithms under real measurement conditions, and not just in simulations. Although there are a large number of reconstruction algorithms reported in the literature for use with CS, the algorithms of this work are selected due to their simplicity and known theoretical performance.

Three samples are used in this work, and an area of 1 cm by 1 cm of each sample is always measured, as well as a monocrystalline silicon (*c*-Si) reference cell, with an active area of 2 cm by 2 cm; an organic photovoltaic (OPV) cell, with an area of 1 cm by 1 cm, non-uniform performance, and a weak current; and a large multicrystalline silicon (*mc*-Si) solar cell with an area of 8 cm by 8 cm, which yields noisy measurements due to its large area. Photocurrent response measurements are acquired at short circuit conditions for all the samples. The samples are presented in Figure 3, with a random pattern projected on them using the DMD optical system. A series of such random patterns are used for compressive sampling.

**Figure 3.** The three samples used in this work, with a random pattern for compressive sampling projected on them. On the left is the c-Si reference cell, in the middle is the organic photovoltaic (OPV) device, and on the right is the large mc-Si cell.

### *2.3. Sensing Matrix Sparsity*

The impact of sensing matrix sparsity on the measurement process and on measurement accuracy can be significant. In this work, 100 × 100 pixel random binary sensing matrices are used, with different levels of sparsity, which means that they have a proportion of pixels in the "on" state between 1% and 99%. In this scenario, 50% means that half of the elements of the sensing matrix are in the "on" state, and the rest are in the "off" state. As a result, the projected patterns on the sample have half of their pixels bright ("on") and the other half dark ("off"). A proportion of 1% simply means that only 1% of the pixels are in the "on" state, resulting in 100 illuminated pixels for a 10,000 pixel projection. As a result, the amplitude of the current response measured, when a series of patterns (rows in a sensing matrix) is projected onto the sample ,will depend on the sparsity of the sensing matrix. This influences signal levels, and so has an impact on the measurement signal-to-noise ratio (SNR). It should be noted that "measurement SNR" is the SNR at the sampling level—the final image SNR of the reconstructed current maps will be lower, and will also depend on the artefacts inserted by the reconstruction procedure. In reality, initial sampling SNR is only one of the factors that influences the final image SNR of the reconstructed image, but it is still a very significant factor for compressive sampling, as will be shown below. Increased sparsity will mean fewer pixels in the "on" state, while reduced sparsity will mean more pixels in the "on" state. Although one could argue that regarding sparsity, 1% and 99% can be the same thing, for the sake of clarity the above convention is adopted throughout this work. This is explained in Figure 4.

It has been demonstrated in CS microscopy that sensing matrix sparsity can have an influence on CS imaging applications [38]. When using very sparse sensing matrices, the probability of having two adjacent pixels in the "on" state at the same time is small. If, in a projected pattern, there are two adjacent pixels in the "on" state simultaneously, the result may be an overlapping excited area in the sample. In CS application cases, as in the optical system of this work, due to light scattering and the diffusion of charge carriers, it may be uncertain to which of the two adjacent pixels the additional measured signal, which contributes to the global current reading of the specific pattern, is generated. Consequently, there may eventually be increased measurement noise in the final reconstructed current map, because of this uncertainty. On the other hand, with very sparse matrices the measured signal is significantly reduced. When using less sparse sensing matrices, many more pixels are in the "on" state, which results in a significant signal amplification, especially when compared with the point-by-point sampling case. The cases when sparser or less-sparse matrices are most appropriate for CS current mapping can eventually depend on the sample to be measured or the background noise of measurements.

**Figure 4.** Visualisation of the sparsity of individual patterns of sensing matrices. Sparser patterns have more dark pixels than bright pixels, which is equivalent to a larger proportion of micromirrors being in the "off" state.

### **3. Results**

### *3.1. Signal Amplification*

The photocurrent signal level that a conventional LBIC system has to accurately measure in order to produce the current map can be in the range of nA. In our case, when the optical system is used to implement a point-by-point photocurrent scan, the current values are indeed in the nA range, as can be seen on the right in Figure 5, for the c-Si reference sample. In the same figure, the values of compressive sampling are also presented. All the values are in the range of 0.45–0.50 mA, which means that the current signal is enhanced by at least three orders of magnitude. This is an important feature that can be highly advantageous in cases where the signal level of individual pixel points is very weak to measure with a point-by-point process without a lock-in system.

**Figure 5.** On the left, the current measurements for 10,000 patterns for compressive sampling—all the values are within a very small value range. On the right, the current measurements for a 10,000-pixel point-by-point current map, where the values have a range of 1 order of magnitude.

All values are within a very narrow range (0.45 mA to 0.50 mA), and this never changes during measurements for a specific sample or a sensing matrix sparsity choice. All the necessary information for reconstructing the current map is within the scatter of the measurements. This means that when acquiring measurements, the minimum and maximum instrument reading range can be set easily in a way that provides a very high dynamic range for the sampling procedure, which can increase accuracy of measurements. This specific feature of compressive sampling is utilised in the next section to correct long-term noise during measurements. Additionally, problematic measurements, such as spikes or zero values, will appear as outliers, and can be excluded easily from the reconstruction process, along with their corresponding pattern. Although the signal levels are greatly enhanced with compressive sampling, actual measurements will still contain noise as any measurement, which always influences the reconstruction process.

In practice, while signal levels are significantly enhanced by using compressive sampling, the background measurement noise levels are kept relatively stable, depending on the measurement settings of the instrument. In order to show the influence of measurement SNR on the method's performance, the SNR was calculated for all samples and cases of sensing matrix sparsity. The results are presented in Table 1. The SNR for every projected pattern during compressive sampling is calculated using 30 samples for each measurement (pattern). The measurement SNR is calculated with the Formula (1)

$$\text{SNR (signal-to-noise ratio)} = \frac{\text{Mean Value}}{\text{Standard Deviation}} \tag{1}$$



SNR is calculated for each projected pattern, and the measurement SNR is the average for all the patterns. The values of average current and measurement SNR for all three samples, and for different cases of sampling procedures, are presented in Table 1. The difference between compressive sampling and the raster scans (point-by-point scans) regarding signal amplitude and SNR is significant. In particular, for the large area mc-Si cell, the dark current present results in high levels of noise for the point-by-point scan. In all cases of different samples, the signal is amplified at least two orders of magnitude compared to the raster scan, as can be observed in Table 1. The sparsity of sensing matrices also has an effect on SNR of measurements. As it can be observed in Table 1 and in Figure 6a, the SNR increases significantly for all of the samples, with decreasing sparsity levels of sensing matrices. Specific "falls" of the SNR trend in Figure 6a (for example, at 30% and 95% for the c-Si reference cell) are due to changes in the measurement range for the photocurrent reading of the instrument as the signal increases. This results in slightly higher background noise levels when a higher value is selected for the range of the instrument, when the measured signal reaches the limit for the previous range. For a specific choice of range, the SNR increases steadily, until it saturates before the range changes. This behaviour can be observed from 30% to 90% of sparsity levels for the c-Si reference cell, for 50% to 90% for the OPV, and from 0% to 90% for the large mc-Si cell. This shows that the choice of sensing matrix sparsity for optimising SNR would also depend on the specific instrument used for measurements.

**Figure 6.** (**a**) Measurement SNR for the three different samples in the case of compressive sampling, with reducing levels of sensing matrix sparsity. (**b**) Signal amplitude (current reading) for the three different samples in the case of compressive sampling, with reducing levels of sensing matrix sparsity. (**c**) SNR as a function of average current of samples, while reducing sparsity levels.

In Figure 6b, the average current measured for each sparsity level is presented. As expected, all the cells demonstrate a linear response, with current increasing while sparsity is decreasing. The silicon devices both demonstrate similar trends, since they have similar efficiencies, while the OPV low-efficiency device produces a much lower current. The falls observed in Figure 6a due to changes in range are not observed in Figure 6b, since what is affected when the range changes is the background noise levels and not the measured signal. The correlation between SNR and measured current is presented in Figure 6c. The same behaviour seen in Figure 6a can be observed for the silicon devices. It is also clear that the OPV device demonstrates the same SNR levels as the c-Si reference device for given measured current values; although the efficiency of the OPV sample is low, and the current gains are not as high as for the silicon devices, its SNR still increases significantly for less sparse sensing matrices.

The influence of measurement SNR can be observed in Figure 7. While a raster scan is possible with this optical system for the two smaller samples, the high noise levels of the large mc-Si sample result in a very noisy current map. On the other hand, with the signal amplification when using compressive sampling, the acquisition of a current map is possible even with such high noise levels. It is clear that the number of pixels in the "on" state of the sensing matrix can be increased in order to amplify the measured signal. This does not affect the reconstruction performance, as will be demonstrated in a following section. Thus, CS current mapping can provide reliable results, even in cases of very weak signals or high noise levels, when a raster scan is not possible. It has to be noted

that the SNR analyzed here is the measurement SNR at the sampling stage, and not the final current map SNR. The final current map SNR will also depend on the choice of sensing matrix, transform, reconstruction algorithm, undersampling level, and of course, the initial measurement SNR that is discussed in this work.

**Figure 7.** Current maps of the three samples used in this work, with the c-Si reference cell on the left, the OPV in the middle, and the large mc-Si cell on the right. In the top row, point-by-point current maps of the three samples. In the bottom row, CS current maps of the three samples, using the orthogonal matching pursuit (OMP) algorithm for reconstruction and for 50% undersampling.

### *3.2. Low-Frequency Noise Correction*

Although a reference measurement for the laser light source has been implemented into the optical system using a photodiode, there is a more convenient and practical way of removing long-term noise during measurements in the case of compressive sampling. Low frequency noise/drift of signal that is independent of the sample's instantaneous performance can be due to laser instability or temperature changes of the sample. Such changes can be easily filtered out when compressive sampling is applied. As described in the previous section, when compressive sampling is applied, the complete measurement set spans within a very small range of values. This range is constant for a specific current map and sensing matrix, and any changes due to long-term noise will appear as drifts from this range. In addition, any spikes or other significant instantaneous changes of laser power are visible as outliers, and can be removed from the measurement set without losing any information. This is because fewer measurements than the pixels of the current map are applied, removing one more measurement, and the corresponding pattern will have no effect on the reconstruction.

This feature is demonstrated in Figure 8, for a case of compressed sensing current mapping of a small area of the large mc-Si sample with drifting measurement data. The laser source power changed slightly over time, simulating the potential effects of temperature or light source instability. This created a drift of the measured signal, which affected the reconstruction process and resulted in a very noisy current map. The signal was unstable and increased slightly over time. The OMP algorithm is used for reconstruction in this case. As can be observed in Figure 8, this small drift results in a noisy reconstruction of the current map. Nevertheless, since the compressively sampled measurements are expected to always be within a short range of values, this noise can be corrected. Even in this case of more intense deformation of the signal, the actual average signal difference due to this drift is around 2.5%. Still, this affects the reconstruction process if there is no correction of the sampled data. Using a generated curve to normalise the data, the drift is completely removed from the sampled data, after using a polynomial fitting to generate a curve on the corrupted sampled data. Although in most cases it is not necessary, this correction procedure is used for all cases in this work, in order to ensure

that any drift of the signal is not affecting reconstruction performance. Such a correction would not be possible with a raster scan, as there would be a chance that real information would be removed.

**Figure 8.** A case of drift correction of sampled data, using 80% pixels "on" patterns. On the top row, the uncorrected map and sampled data, at the bottom, the corrected map and data.

### *3.3. Reconstruction Performance*

For a quantitative evaluation of the performance of the method, depending on sensing matrix sparsity, the point-by-point and reconstructed current maps were compared at a pixel-by-pixel level. This is straightforward to achieve using the DMD optical system, as it includes no moving parts and in both sampling cases, the coordinates of the current maps coincide accurately. Pearson's correlation coefficient ρ(**x**ˆ, **x**) was calculated for different levels of undersampling used for reconstruction, for different levels of sensing matrix sparsity and for the two different algorithms. The correlation coefficient is calculated by dividing the covariance of the point-by-point and reconstructed current map by the product of their standard deviations:

$$\rho(\mathbf{\hat{x}}, \mathbf{x}) = \frac{\text{cov}(\mathbf{\hat{x}}, \mathbf{x})}{\sigma\_{\mathbf{\hat{x}}} \sigma\_{\mathbf{x}}} \tag{2}$$

where **x** is the point-by-point current map, and **x**ˆ is the CS-reconstructed current map, both in vector form. For the case of the large mc-Si sample, a pixel-by-pixel comparison is not possible with this optical system. Due to the large area of this sample, and since no lock-in is used, the raster-scanned current map is very noisy due to high dark current, and cannot be used as a reference for the reconstructed current maps for a pixel-by-pixel comparison.

In Figure 9, the reconstructed current maps of the c-Si reference cell and the OPV cell are presented along with the raster scan, using the same DMD optical system. By using compressive sampling, the current maps were acquired with half the number of measurements that the raster scan required. A number of sensing matrices with different sparsity levels were used, from 1% of pixels in the "on" state up to 99%. As it can be observed in Figure 9, for this sample, and for a given algorithm, all sensing matrices with different sparsity levels exhibit similar reconstruction performance. In all cases, defects

like broken fingers in the silicon device and non-uniformities in the OPV device are clearly imaged. In the case of sensing matrices with 99% of the pixels "on", it is almost as if the whole sample is illuminated, significantly increasing the measured signal and measurement SNR without affecting the accuracy of the reconstructed current map.

**Figure 9.** CS current mapping with sensing matrices with different sparsity levels (number of pixels in the "on" state). On the top left, the point-by-point scan is also included for comparison. It can be observed that the differences in reconstruction performance for different sensing matrices are negligible.

On the other hand, when using different reconstruction algorithms, the reconstruction performance can vary. In Figure 10, the correlation coefficient between the point-by-point and the CS current maps for the two samples is presented as a function of measurements used for reconstruction, for sensing matrices with different sparsity levels, and for two different reconstruction algorithms ( 1, OMP). Although for the c-Si reference cell the differences between different algorithms are not significant, for the OPV sample the reconstruction performance varies significantly between the two algorithms. This shows that some algorithms can have a different performance for different samples, depending on the features of the current map. This has to be taken into consideration when choosing a reconstruction algorithm. Nevertheless, as can be observed in the graphs of Figure 10, sensing matrix sparsity does not affect reconstruction performance for a given reconstruction algorithm. This shows that the right sparsity level of the sensing matrices can be set each time, considering background noise, signal levels, and equipment sensitivity for acquiring the current map of a specific sample.

**Figure 10.** Correlation coefficient as a function of the number of measurements used for reconstruction, for sensing matrices with different sparsity levels, and for the two different algorithms. On the left, the graph for the c-Si reference cell; on the right, the results of the OPV cell.

When approaching 100% of measurements used for reconstruction, the performance of the algorithms, especially that of the 1 algorithm, declines. This is because when measurement noise is included, the optimisation algorithm fails to find a solution for 100% of measurements used, as has been previously demonstrated in [22]. This is because the optimisation algorithm is increasingly constrained as we reach 100%, and has fewer degrees of freedom to filter out noise. There are algorithms available in the literature that can expect some noise in the measurements, and would not have such issues when approaching 100%. For the graphs in Figure 10, reconstruction was implemented for 99.0% as well as 99.9% of sampling, in order to accurately draw these curves. In reality, this area of undersampling is meaningless for compressive sampling, and such problems will not arise in real applications.

### **4. Conclusions**

Compressed sensing photocurrent mapping provides fast and reliable measurements with simple experimental layouts. In this work, the signal amplification gains and the ability to optimise CS current mapping by controlling sensing matrix sparsity levels is demonstrated. By setting the right sparsity levels of sensing matrices, a significant increase in the SNR of measurements can be achieved. This provides the means to acquire current maps of samples with very weak currents or high dark currents, where a point-by-point scan would fail. In addition, current mapping systems can be put together without the need of a lock-in amplifier, allowing measurements when the application of lock-in techniques is not possible. The nature of compressive sampling allows long-term noise correction to be applied without the need of a reference measurement of the light source. For this experimental application of CS, the selected sensing matrix sparsity for optimum signal amplification does not affect current map reconstruction performance for a given reconstruction algorithm. As a result, sensing matrix sparsity can be a crucial setting that can be controlled in order to optimise measurement accuracy of CS current mapping.

It is apparent from the results of this work that in CS current mapping, different reconstruction algorithms behave differently for different samples. A future investigation of different algorithms and transforms for this CS application is necessary in order to fully optimise CS current mapping. Since a direct pixel-by-pixel comparison with a raster scan is possible, the DMD-based optical current mapping system of this work offers the opportunity to investigate the performance of different algorithms and transforms for compressive sampling. In this way, tools for this CS application can be evaluated experimentally in a realistic way, including instrument noise and system specific features. Such an evaluation of CS tools can be useful for other optical CS applications, where a comparison with a point-by-point scan is not always possible.

**Author Contributions:** Conceptualization, G.K.; methodology, G.K.; software, G.K.; formal analysis, G.K.; investigation, G.K.; data curation, G.K.; writing—original draft preparation, G.K.; writing—review and editing, G.K., J.C.B., and F.A.C.; visualization, G.K.; funding acquisition, J.C.B. and F.A.C.

**Funding:** This research was funded by the European Metrology Programme for Innovation and Research (EMPIR) 16ENG02 PV-Enerate, co-financed by the participating states and from the European Union's Horizon 2020 research and innovation programme.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
