Clutter suppression in ultrasound: performance evaluation and review of low-rank and sparse matrix decomposition methods

Abstract

Vessel diseases are often accompanied by abnormalities related to vascular shape and size. Therefore, a clear visualization of vasculature is of high clinical significance. Ultrasound color flow imaging (CFI) is one of the prominent techniques for flow visualization. However, clutter signals originating from slow-moving tissue are one of the main obstacles to obtain a clear view of the vascular network. Enhancement of the vasculature by suppressing the clutters is a significant and irreplaceable step for many applications of ultrasound CFI. Currently, this task is often performed by singular value decomposition (SVD) of the data matrix. This approach exhibits two well-known limitations. First, the performance of SVD is sensitive to the proper manual selection of the ranks corresponding to clutter and blood subspaces. Second, SVD is prone to failure in the presence of large random noise in the dataset. A potential solution to these issues is using decomposition into low-rank and sparse matrices (DLSM) framework. SVD is one of the algorithms for solving the minimization problem under the DLSM framework. Many other algorithms under DLSM avoid full SVD and use approximated SVD or SVD-free ideas which may have better performance with higher robustness and less computing time. In practice, these models separate blood from clutter based on the assumption that steady clutter represents a low-rank structure and that the moving blood component is sparse. In this paper, we present a comprehensive review of ultrasound clutter suppression techniques and exploit the feasibility of low-rank and sparse decomposition schemes in ultrasound clutter suppression. We conduct this review study by adapting 106 DLSM algorithms and validating them against simulation, phantom, and in vivo rat datasets. Two conventional quality metrics, signal-to-noise ratio (SNR) and contrast-to-noise ratio (CNR), are used for performance evaluation. In addition, computation times required by different algorithms for generating clutter suppressed images are reported. Our extensive analysis shows that the DLSM framework can be successfully applied to ultrasound clutter suppression.

Background

Angiology, which concerns vessel-related diseases, is one of the most important branches of medical science since vascular diseases are very common and cause death to a large number of people every year [1]. Vascular diseases can primarily be divided into several categories based on the type of vessel. Arterial diseases include aneurysms, thrombosis, vasculitides, and vasospastic disorders. Venous diseases include venous thrombosis, chronic venous insufficiency, and varicose veins. There are also diseases associated with capillaries. One such example is the capillary hemangioma. Currently, the most accepted classification of vascular abnormalities is tumors and deformities which were adopted in 1996 by the International Society for the Study of Vascular Anomalies [2]. Therefore, many major clinical diseases have been shown to cause vascular growth abnormalities. For example, many cardiovascular diseases are related to aneurysms or other vascular variations [3, 4]; the growth of many tumors in cancer is also highly dependent on angiogenesis [5, 6]. Similarly, angiogenesis is also an important feature of diabetes-related diseases [7,8,9] and endometriosis [10]. Therefore, blood vessel imaging is indispensable in clinical fields and medical research [11], including but not limited to diagnosis, treatment planning, surgery, and follow-up treatment results.

Several medical imaging modalities such as duplex ultrasound (DUS), computed tomography (CT), magnetic resonance imaging (MRI), and digital subtraction angiography (DSA) have been employed thus far to ensure a proper visualization of blood vessels. Among different vascular imaging modalities, ultrasound has become the primary choice, for it is safe, economical, easy-to-use, and most importantly, real-time [11]. Duplex ultrasound is the combination of color flow imaging (CFI) and grayscale/brightness mode (B-mode) imaging, whereas the CFI is used to observe the blood flow direction and velocities, and the B-mode ultrasound is used to visualize two-dimensional anatomy images simultaneously. By simultaneous processing frequency, phase, and amplitude of the backscattered ultrasound signal, CFI can rapidly identify the flow direction and velocities in the region of interest. Moreover, CFI can be used to mark flow abnormalities, including stenoses and occlusions [12]. The comparison between ultrasound and other vascular imaging methods is shown in Table 1.

Table 1 A comparison of vessel imaging methods [11, 13]

Due to the excellent performance, ultrasound CFI has been increasingly used for the diagnosis of vessel-related diseases [14]. However, as one of the most promising and widely applicable methods with low cost and no risk, CFI still has some obvious disadvantages. First, due to the tissue scattering of the ultrasound beam, the intensity of the blood backscatter is several orders of magnitude less than that of the tissue backscatter, which makes it hard to image blood flow clearly [12, 15]. Second, more than three pulses are needed to estimate the velocity because of the stochastic behavior of blood signals and the impact of tissue clutters [12]. The requirement for multiple pulses limits frame rates and the number of scan lines. Third, CFI is limited by the insonation angle which is the angle between the ultrasound beam and the flow direction [16, 17]. Generally, an accurate measurement requires Doppler angles ranging from 30° to 60°, where smaller angles will result in lower speeds and greater angles will produce a significant overestimation of the velocity [16, 17]. Last but not least, blood signals and clutter signals will possess a significant overlap, that is, when the blood flow rate is very slow (such as in small blood vessels) or when the tissue movement is obvious. The overlap will be harmful to blood vessel visualization [18, 19]. Most of these disadvantages are caused by clutter, as a consequence, clutter suppression is particularly important in ultrasound blood flow imaging. Figure 1 shows the clutter in two CFI images and illustrates the importance of clutter filtering.

Fig. 1

A set of comparison images showing CFI with and without clutter filters. a CFI raw data in Brightness mode. b The same data after clutter suppression by SVD. In the upper right window, the raw CFI data contain a lot of tissue clutter in the background, which is suppressed by SVD in the second image

The main purpose of clutter filtering is to suppress gross-moving tissue clutter and beam sidelobe leakages [19]. An efficient clutter suppression is a prerequisite for CFI to present accurate and clear blood flow maps. The most significant effect of clutter reduction is an increase in the signal-to-noise ratio (SNR) of the blood signal, which enables clearer blood flow maps and reduces erroneous moving tissue signals. Meanwhile, pure blood flow signals also help reduce the number of pulses needed to estimate the speed, thereby increasing the frame rate. In addition, the overlapping frequency spectra of slow blood flow and fast-moving tissue will no longer hinder the microvascular flow detection or add bias to high-velocity flow [19, 20].

However, the perfect removal of clutter signals is still impossible for now since clutter signals are 40 to 100 dB stronger than blood signals and they exhibit similar properties [15].

In early development of CFI clutter filtering, tissue signals and blood signals were assumed to have completely different frequency characteristics. This assumption holds that tissue and blood signals exhibit non-overlapping frequency spectra since the tissue is considered to be nearly stationary whereas red blood cells are rapidly moving [18]. Based on this assumption, finite impulse response (FIR) and infinite impulse response (IIR) high-pass filters were used to filter clutter signals and enhance the sensitivity of blood flow [15, 18]. Nowadays, it is recognized that FIR and IIR filters have distinct drawbacks. FIR filters require a high order to separate blood from clutter, whereas IIR filters take a long time to settle [18, 21]. Furthermore, both types of high-pass filters suffer from the insufficient number of slow time samples, which leads to inefficient suppression of clutter [19, 22]. Another clutter removal approach introduces linear regression filters [23,24,25]. The regression filter eliminates clutter signals by taking the least square fitting of signals from the signal model [18]. Studies suggest that polynomial regression filters and IIR filters have better performance than FIR filters. In the case of contrast-enhanced ultrasound vascular imaging, pulse inversion technique has been introduced toward the end of clutter rejection [26,27,28]. In this approach, the linearity property of tissue echo is exploited for distinguishing tissue from blood [26, 29, 30]. Although these methods significantly improve the SNR of blood signals, they are not considered in this paper because of their invasiveness.

The aforementioned traditional clutter suppression algorithms, such as FIR and IIR, have at least one of the following issues: (1) long settling time, (2) inability to adaptively suppress the clutter based on data property, (3) inadequate temporal sample or resolution. Besides, two main reasons are resulting in the imperative innovation of ultrasound clutter filtering. First, new ultrasound technologies like plane wave ultrasound have brought a higher frame rate and imaging speed. Traditional filters cannot meet the higher clutter filtering performance requirements, though they do not suffer from the settling time due to the high frame rate. Second, the underlying assumption of traditional filters does not hold in the presence of significant tissue motion stemming from the sonographer’s sinusoidal hand movement or the patient’s breathing and heartbeat [31, 32]. In such a scenario, the frequency bands corresponding to tissue and blood overlap with each other without a definite boundary between them. Hence, high-pass filters fail to separate blood from tissue when the clutter signal dominates with non-zero Doppler frequency caused by substantial tissue movements.

To resolve these issues, eigen-based filters [33,34,35] have been proposed which take both spatial and temporal samples into consideration to develop an adaptive clutter suppression scheme. The techniques related to these eigen-based filters have been widely applied in the field of computer science which is mainly used for processing high-dimensional data. Meanwhile, these techniques are not based on incomplete traditional assumptions. Matrix decomposition is the principal idea behind these algorithms and it is assumed that clutter and blood signals lie in different subspaces. Therefore, eigen-based filters are considered adaptive to gross motions induced by the sonographer or the subject being examined. Based on different assumptions, research proves that eigen-based filters perform better than traditional methods [20, 22].

Most of the eigen-based filters for ultrasound clutter suppression are based on singular value decomposition (SVD) or eigenvalue decomposition and improve upon it [36,37,38,39]. To perform the subspace separation task, slow-time temporal ultrasound frames are stacked as columns of data matrix, known as the Casorati matrix [40]. The SVD of this Casorati matrix provides the opportunity to distinguish blood from clutter. It has been reported in the literature that the most dominant singular values and vectors correspond to clutter, the next few represent blood and the least significant ones correspond to noise [41]. In these eigen-based approaches, the eigen or singular values representing clutter and noise are set to zero to find the blood component of the echo signal [41, 42].

Many SVD-based techniques have been proposed which work with conventional line-by-line scanning [20, 43,44,45]. These methods suffer from lacking an adequate number of temporal samples due to low frame rate associated with focused ultrasound imaging [21]. Recent clutter suppression algorithms [18, 42, 46,47,48,49] have resolved this issue by incorporating ultrafast plane-wave imaging. However, the blood signal in plane-wave ultrasound is even weaker than normal ultrasound due to the unfocused wave [50, 51]. The sidelobe in plane-wave imaging is also much higher than that in conventional imaging due to the same reason. Therefore, plane-wave ultrasound has a higher and more urgent filtering requirement than traditional CFI. Recent methods have extended SVD-based clutter suppression to a higher order by analyzing a data tensor instead of a two-dimensional matrix [42, 48, 52]. Since the first few singular values do not necessarily correspond to the clutter signal in the presence of a large temporal misalignment among the frames, the motion correction step has been introduced in SVD-based clutter rejection [53]. Since SVD was initially combined with plane-wave imaging in 2015, almost all the clutter suppression research have been based on plane-wave ultrasound since SVD can reach its full potential on large datasets [18].

Although SVD-based techniques are promising for suppressing clutter optimally, they have two major drawbacks. First, there is still no uniform and efficient standard for rank selection which presents boundaries of tissue and blood flow [42]. Proper subspace rank selection, which is done by extensive manual intervention, is crucial for the optimality of clutter rejection. Recent methods suggest different criteria for selecting the optimal ranks [54]. In addition, [21] proposes K-means clustering of the decomposed components as a criterion for selecting singular values and vectors corresponding to clutter and blood. Though different ideas are proposed for automatic rank selection [55], there is still no efficient and standard method. An example that briefly explains the problem of SVD threshold selection is shown in Fig. 2. The selected rank will affect blood signals. A large threshold range cannot effectively filter clutter and noise, and a small range will lose part of the blood signals. The second drawback is that SVD is sensitive to noise. It fails to obtain the optimal result while processing data with large random noise [56].

Fig. 2

A set of pictures showing the threshold selection of SVD. a The original simulation data in brightness mode. b–e The processed images by SVD with different thresholds. Parameters b and e represent the selected rank of blood and noise signal, respectively. The full rank of the data is 20

The aforementioned issues can potentially be resolved by taking the framework called decomposition into low-rank and sparse matrices (DLSM) [21] into account. SVD is one of the calculation methods in DLSM framework and there are also approximate SVD or SVD-free algorithms. This is a well-established framework in the field of computer vision due to its robustness to large noise and information corruption [56]. The underlying assumption of this approach is that steady tissue is a low-rank component and moving blood exhibits sparsity [50]. It has been noticed that both temporal and spatial information can be used to separate tissue and blood signals since tissue signals have a higher temporal–spatial coherence than blood signals (e.g., the blood scatterers are unique and constantly changing). A convex optimization problem is usually solved to decompose the data matrix into low-rank clutter and sparse blood components. A recent technique has used this model for the concurrent removal of clutter and noise [57]. Furthermore, recent work has incorporated deep learning with low-rank and sparse decomposition for improved clutter suppression performance [21].

The main purpose of this work is to demonstrate the feasibility of 106 established low-rank and sparse decomposition algorithms in ultrasound clutter suppression and to provide suggestions for most suitable DLSM models, optimization methods, and algorithms for ultrasound clutter suppression. The paper is organized as follows. “Decomposition into low-rank and sparse matrices (DLSM) framework” section illustrates DLSM framework including decomposition types, loss functions, and relationships with tensor decomposition. “Experiments” section includes detailed experimental settings and results on simulation, phantom, and in vivo rat datasets. Finally, the paper reveals our discoveries in “Discussion” and “Conclusion” sections including the appropriate DLSM algorithms for clutter suppression and the shortcomings of the remaining algorithms. Comments and possible solutions are also proposed in response to different shortcomings.

Decomposition into low-rank and sparse matrices (DLSM) framework

Low-rank and sparse structures are attractive since they usually represent part of the large and high-dimensional data which we are most interested in. Noise and data corruption can be fixed when decomposing matrices into low-rank and sparse components. Methods like sparse representation and low-rank modeling have achieved great success in computer vision, natural language processing, system identification, bioinformatics, etc. [58,59,60]. So far, many different models, optimization methods, and algorithms are proposed aiming at solving the low-rank and sparse matrix recovery problems. Meanwhile, many classifications have been proposed [58, 61,62,63] according to linearity, convexity, number of subspaces, or number of addition matrices.

Decomposition into low-rank and sparse matrices (DLSM) is one of the relatively detailed and comprehensive frameworks [61] which classifies various models of matrix decomposition according to the number of constrained component matrices. DLSM framework provides a suitable framework for signal processing, system identification, computer vision, machine learning, etc. This decomposition idea is becoming more popular and widely used in recent years, especially after the robust principal component pursuit (RPCP) was purposed in papers of Candes et al. [56], and Chandrasekharan et al. in 2009 [64]. In the beginning, these algorithms are designed to deal with high-dimensional data which are often regarded as an extremely high-dimensional data matrix. Since many dimensions are usually independent, it is possible to recover the matrix from corruption or noise. These ideas are based on the assumption that the uncorrupted information matrix is highly correlated within the observing time-window and therefore lies in the low-rank subspace. At the same time, the moving foreground objects, noise, or other special signals constitute the correlated sparse outliers.

Based on similar assumptions, several algorithms under DLSM framework have been validated that they can be successfully applied to ultrasound clutter suppression [18, 22, 36, 38, 55, 65]. In medical ultrasound, tissue and blood flow also lie in different subspace. In terms of temporal information, tissue signals and blood signals have different spectral features due to the different movement patterns of blood and tissue. As for spatial features, the blood signal has an extremely lower spatial coherence than tissue signal because the irregular movement and arrangement of red blood cells produce constantly changing scatterers, whereas the tissue movement is overall patterned. Therefore, they gain a low rank and sparsity characteristics, respectively, and lie in different subspaces. Due to the robust and efficient performance of DLSM frameworks in separating low-rank and sparse components, it can show great potential in the field of ultrasound clutter suppression.

Overall, DLSM framework is divided into decomposition problems, minimization problems, loss function and solvers (algorithms used to solve the optimization problems) [63] as Fig. 3 shows. The permutations and combinations of models and optimization methods and solvers lead to various algorithms, which is the origin of the DLSM framework. DLSM framework and its application in the ultrasound clutter suppression will be briefly illustrated in the following subsections.

Fig. 3

The schematic diagram of DLSM framework. DLSM framework contains 5 branches, which are models (or called math formulations), decomposition problems, minimization problems, loss functions, solvers (or called algorithms). Examples are shown beside the branches

Preprocessing and notations

Preprocessing of ultrasound data is necessary for integration into an input matrix or tensor in a special shape when applying DLSM algorithms. In general, the input of the DLSM algorithm consists of a sequence of n consecutive ultrasound data ($F_1\dots F_n$) with the original size of $F\in \mathfrak {R}^{i_1\times i_2}$. For a two-dimensional DLSM algorithm, the input M ($M\in \mathfrak {R}^{m\times n}, m=i_1\times i_2$) is in matrix form in most cases which consists of n resized ultrasound data frames ($F\in \mathfrak {R}^{m\times 1}$) arranged in order. In terms of high-order DLSM algorithms, the input is typically an N-order tensor T ($T\in \mathfrak {R}^{t_1\times t_2\dots t_n}$). T is generally third order and concatenated by original size ultrasound frames, where $T=[F_1^{i_1\times i_2},\dots F_n^{i_1\times i_2}], T\in \mathfrak {R}^{i_1\times i_2\times n}$. Next, the input M (or T) is decomposed into several components through the DLSM algorithm as follows:

$$\begin{aligned} M = \sum \limits _{x=1}^X K_x \end{aligned}$$

(1)

where $0\le X \le 3$ and $K_1,K_2,K_3$ typically represent low-rank L, sparse S, and noise components E, respectively. The specific components $K_x$ and the number of X depend on the purpose (interested in sparse or low rank components) and the decomposition formulation.

Decomposition formulations

Implicit decomposition

Implicit decomposition $(X=1)$: Under the condition that x is equal to 1, matrix M is approximately equal to a target low-rank matrix L under the constraint condition, because the information that people interested in mainly lies in the low rank component in most cases. Sparse matrix S can be obtained from the difference between M and S (e.g., $S=M-L$). However, this processing is the opposite in the application of ultrasound clutter suppression because the blood signal is relatively sparse. The formulation of this problem is as follows:

$$\begin{aligned} \min \, f(M,L) \quad \text {s.t.}\; L \end{aligned}$$

(2)

where $M\approx L$, f(.) is a loss function used for the minimization term which depends on specific solvers or algorithms. Models like principal component analysis (PCA), non-negative matrix factorization (NMF), and matrix completion (MC) are in this category.

For the applications targeted to sparse components, implicit decomposition sets the target matrix $K_1$ as a sparse matrix S. Then low-rank matrix L is the difference between M and S which can be calculated as $L=M-S$. Sparse dictionary learning [66], sparse linear approximation, and compressive sensing [66,67,68], etc. are built under the same idea.

$$\begin{aligned} \min \, f(M,S) \quad \text {s.t.}\; S \end{aligned}$$

(3)

where $M\approx S$, and the difference contains noise and other information. In this case, implicit decomposition can be used in the compressed sensing and signal recovery similar to unsupervised clustering [69] and image recognition [70].

Before more robust explicit decomposition method was proposed, the main development of ultrasound clutter suppression was based on PCA or SVD or eigenvalues, which belong to implicit decomposition [19, 22, 35, 38, 39, 45]. Although many experiments have proved that these eigen-based filters greatly improve the performance than traditional IIR and regression filters, many authors realize that the filtering method based on implicit decomposition is not robust to accelerated tissue movements and different kinds of noise [22, 35, 39]. Moreover, their expensive computational complexity is not suitable for real-time processing.

Explicit decomposition

Explicit decomposition $(X=2)$: Under this condition, M is usually decomposed into a low-rank matrix $K_1=L$ and a sparse matrix $K_2=S$ ($M\approx L+S$). This is called explicit decomposition because there are two constraints. One is sparse constraint over S and the other is low-rank constraint over L. Therefore, explicit decomposition is more robust than implicit decomposition. The formulation of explicit decomposition is as follows:

$$\begin{aligned} \min \, f(L)+f(S) \quad \text {s.t.}\; L,S \end{aligned}$$

(4)

where $M\approx L+S$ and f(.) represents loss function. The explicit decomposition includes robust principal component analysis (RPCA), robust non-negative matrix factorization (RNMF), robust matrix completion (RMC), and robust subspace tracking (RST) [59, 71].

These methods generally work better and are more robust than implicit decomposition because of the additional constraints [71]. In this way, RPCA has been used as a powerful tool in MRI, CT, and ultrasound imaging [72,73,74]. Many optimization algorithms have been proposed for cluster suppression in ultrasound imaging using RPCA, RMC [21, 55, 75].

Stable decomposition

Stable decomposition $(X=3)$: Due to the fact that there are always noise and corruption caused by special cases in the real world, an additional matrix $K_3$ is added to represent unexpected components. $K_3$ could represent distortion, shadows, and noise according to special situations ($M\approx S+L+N$). It is more stable than the explicit decomposition since more detailed information is separated and taken into account. The stable decomposition can handle more complex situations in the real life such as dynamic videos and maritime monitoring videos which are corrupted by complicated noise.

$$\begin{aligned} \min \, f(S)+f(L)+f(N) \quad \text {s.t.}\; L, S \end{aligned}$$

(5)

Stable decomposition methods include Stable Principal Component Analysis (Stable PCA) or Stable Non-negative Matrix Factorization (Stable NMF) and Three Term Decomposition (TTD). These methods can deal with more complex situations. In terms of US imaging, it is usually assumed that signal M contains clutter signals L (low-rank), blood signals S (sparse), and noise N. Since ultrasound signals have complex noises and dynamic clutter signals, this assumption $M=S+L+N$ is more acceptable when there are meticulous requirements such as microvascular imaging. Although some literature mentions the stable decomposition of blood (L, S, E respectively represent blood flow signals, clutter signals, and noise), they do not illustrate whether constraints are added to noise component. Therefore, the stable decomposition formulation is still a promising research area for ultrasound clutter suppression.

Models under DLSM framework

As of today, many models, also called problem formulations, have been proposed. According to different math formulations and features, methods are usually classified under families such as robust principal component analysis (RPCA), non-negative matrix factorization (NMF), matrix completion (MC), and Subspace Tracking (ST). Different models have different functions and aims. However, it has been proved that the solutions of many robust models can be mutually expressed in closed forms [76]. For instance, RPCA via principal component pursuit [56] can be considered as MC models using $l_1$-norm loss function [63]. In addition, these models can be flexibly generated in any decomposition formulations. For example, adding constraints on noise components based on the RPCA will change it from explicit decomposition to stable decomposition.

Robust principal component analysis

Principal component analysis (PCA) generates a set of linearly uncorrelated variables which is called principal components, from a set of observations by orthogonal transformation. Similar mathematical tools include SVD and eigenvalue decomposition. RPCA is based on the extension of PCA (expansion from implicit decomposition to explicit decomposition), which aims to recover low-rank components and reduce the impact of grossly corrupted data. RPCA can be approached by principal component pursuit (PCP) [56, 64], Bayesian RPCA [77,78,79], and so on. RPCA problem is generally expressed as follows:

$$\begin{aligned} M = L + S \end{aligned}$$

(6)

where L is low-rank matrix and S is sparse matrix. According to the nature of L and S, the most intuitive way to solve the RPCA problem is to minimize the rank of L and the $l_0$-norm of S:

$$\begin{aligned} \min _{L,S}\, rank(L) + \lambda \Vert S\Vert _{l_0}\quad \text {s.t.}\;M-L-S=0 \end{aligned}$$

(7)

where $\lambda$ is a balanced parameter. However, this formulation is NP-hard. Therefore, optimization problems like PCP are needed.

The convex optimization principal component pursuit (PCP) was first proposed by Candes et al. [56, 63, 80] to address the RPCA problem. It becomes one of the most famous methods of face recognition and background modeling in recent years. PCP uses the following formula to convexly optimize RPCA problem:

$$\begin{aligned} \min _{L,S}\, \Vert L\Vert _{*} + \lambda \Vert S\Vert _{l_1}\quad \text {s.t.}\;M-L-S=0 \end{aligned}$$

(8)

where $\Vert .\Vert _{*}$ and $\Vert .\Vert _{l_1}$ are the nuclear norm and $l_1$-norm, respectively. Although this method excels in computer vision, there are still some limitations to sparse components recovery. First, it requires expensive computational algorithms. Second, it is a batch method which is not suitable for real-time applications, especially for plane-wave ultrasound with high frame rates. Third, it has very high requirements for low rank and sparse properties; however, the complex blood flow or noise may make it difficult for ultrasound data to meet such requirements. To accelerate the algorithms and achieve higher precision, different solvers have been developed [81,82,83]. Solvers for real-time implementations have also been proposed [84, 85].

The stable principal component pursuit (SPCP) is a stable expanded form based on PCP, which mainly aims at reducing the impact of noise. SPCP adds noise term E based on PCP and constrains it by Frobenius norm.

Matrix completion

The main purpose of matrix completion (MC) is to recover low-rank observation matrix of its missing entries. The Netflix movie rating matrix recover problem is one of the most classic examples. The classic low-rank matrix completion problem can be seen as finding the lowest rank of the matrix L which matches the matrix M, for all the measured entries in set $\Omega$. The mathematical formulation of MC problem is as follows:

$$\begin{aligned} \min _{L}\, rank(L)\quad \text {s.t.}\;L_{m,n}=M_{m,n}\;\forall i,j\in \Omega \end{aligned}$$

(9)

Due to the implicit decomposition of MC is not robust to noise which only affects a small-scale data [86, 87], MC is generally extended to explicit decomposition by adding restrictions, which is called robust matrix completion (RMC). The common RMC obtains stronger robustness than MC by adding sparse constraints, and its formulation after convex optimization is as follows:

$$\begin{aligned} \min _{L,S}\, \Vert L\Vert _{*} + \lambda \Vert S\Vert _{l_1}\quad \text {s.t.}\;P_\Omega (L+S)=P_\Omega (M) \end{aligned}$$

(10)

where $P_\Omega (M)$ is the projection of the complete dataset on the measured entries $\Omega$. Although the form of decomposition is the same as PCP, the unique constraints of RMC make it supervised while the PCP is unsupervised learning [63], which is consistent with the purpose of RMC.

Nonnegative matrix factorization

The nonnegative matrix factorization (NMF) is also a widely used matrix factorization and dimension reduction model under DLSM framework. The main unique feature of NMF is that low-rank factor matrix is subject to nonnegative constraints consistent with the physically natural features in many fields [88, 89]. The NMF problem is generally expressed as follows:

$$\begin{aligned} M \approx WH^\top \end{aligned}$$

(11)

where $W\in \mathfrak {R}^{m\times k}$ and $H\in \mathfrak {R}^{n\times k}$ are two nonnegative matrices, and $k<\min \{m,n\}$ due to the goal of dimension reduction. The most common formulation for the optimization problem of NMF is as following:

$$\begin{aligned} \min _{W,H}\, f(W,H) = \Vert M-WH^\top \Vert _F^2\quad \text {s.t.}\;W\ge 0,H\ge 0 \end{aligned}$$

(12)

where $\Vert .\Vert _F^2$ is the Frobenius norm. The problem (14) is a non-convex problem and it is NP-hard to find its global minimum [88, 90]. Consequently, optimization algorithms and solvers are developed for the local minimum.

Subspace tracking

The subspace tracking (ST) can be regarded as the dynamic RPCA designed to handle increasing new data or dynamic subspaces. The data at each moment t are processed as the increments and then discarded. This idea addresses the problem when new observations come in asynchronously in online streaming environments. It makes subspace tracking more efficient and less computationally expensive on extremely long data sequences [91]. Since ST can recover subspaces from incomplete frame vectors, it has the potential to further improve efficiency by downsampling the input frames [63]. The general formulation for the ST problem is as follows:

$$\begin{aligned} m_t = \sum _{x=1}^X k_x = l_t+s_t+e_t \quad for\; t=1,2,\dots ,n;\; X\in {1,2,3} \end{aligned}$$

(13)

where $m_t$ is input frame data at time t, and $l_t$, $s_t$, $e_t$ are low-rank, sparse, and noise components of $m_t$. The number of k is determined according to different decomposition forms, and the constraint conditions and approximate approximations on each component are determined according to different optimization methods.

Low-rank representation

Low-rank representation (LRR) can also be called low-rank optimization or low-rank minimization. Other unclassified models can be regarded as LRR. LRR is a minimization problem in mathematics. In LRR, the cost function measures the fit between the input matrix M and the approximation matrix L [63]. The mathematical formulation of LRR problem is as follows:

$$\begin{aligned} \min \, \Vert M-\hat{M}\Vert _F\quad \text {s.t.}\;rank(L)\le r \end{aligned}$$

(14)

where M is the input matrix, $\hat{M}$ the approximate matrix, $\Vert .\Vert _F$ the Frobenius norm, and r the rank. The basic form of LRR is similar to other models; therefore, most of the other unclassified models can be regarded as a category in LRR. For instance, RPCA and NMF can be obtained by similar architectures. Constraints other than rank constraints can be added for specific applications. LRR can be extended into an explicit or stable form by adding constraints on the sparse and noise components.

The extension to tensor

In DLSM framework, only some of the single-dimensional information is used when images are pretreated into data matrix M as vectors. This means that some multidimensional information is not taken into account in the process of decomposition. To improve the results, the tensor decomposition is proposed.

Tensor DLSM

When it comes to tensor, the most intuitive idea is to change all matrices to tensors directly since a tensor can be seen as a combination of several matrices. It is very similar to DLSM framework which subjects to $T=L+S+E$. The tensor DLSM extends all components to a tensor form as Fig. 4.

$$\begin{aligned} T=L+S+E \end{aligned}$$

(15)

where T, L, S and N represent the data tensor, low-rank tensor, sparse tensor and noise tensor, respectively. Similar to the matrix DLSM framework, it can be optimized and solved by a minimization problem. Some other classic matrix decomposition optimization methods have also been extended to tensor. The tensor robust principal component method [92] has been proposed based on tensor singular value decomposition (t-SVD) [93]. It has been demonstrated the effectiveness of image denoising. Another robust low-rank tensor recovery model based on RPCA has also been published for complex multilinear data analysis [94]. Rank sparsity tensor decomposition (RSTD) [95] and some other ideas based on stable principal component pursuit (PCP) also have been utilized in image processing.

Fig. 4

The illustration of tensor decomposition

Tensor decomposition

There are two classical tensor decomposition forms which are CANDECOMP/PARAFAC (CP) decomposition and Tucker decomposition [96]. Given a tensor $T\in \mathfrak {R}^{t_1\times t_2\times \cdots t_n}$, the CP decomposition and Tucker decomposition can be modeled as follows:

• Tucker decomposition

  $$\begin{aligned} T=g\times \prod _{i=1}^N U_i+\varepsilon \end{aligned}$$

  (16)

where $g\in \mathfrak {R}^{r_1\times r_2\times \cdots r_n}$ is the core tensor, r is the rank of factor matrix $U_i\in \mathfrak {R}^{t_i\times r_i}$, and $\varepsilon$ represents the residuals. Figure 5 is a schematic representation of the Tucker decomposition. The Tucker decomposition is usually regarded as a non-convex optimization problem [63]. Two most famous and widely used solvers for Tucker decomposition are Tucker-ALS based on alternating least squares [96] and Tucker-ADAL based on alternating direction augmented Lagrangian [94]. SVD based on Tucker decomposition is generally called high-order singular value decomposition (HOSVD) [97, 98], which calculates the singular values of the three expansions $U_1,U_2,U_3$ of a three-dimensional tensor under Tucker decomposition. HOSVD-based ultrasound clutter optimization has been proposed [52, 99] and proved to be more robust to low sampling rates than SVD.

• CP decomposition

  $$\begin{aligned} T=U_1\circ U_2 \cdots \circ U_R+\varepsilon \end{aligned}$$

  (17)

where R is the number of the components, $U_i\in \mathfrak {R}^{t_i\times r_i}$, $\varepsilon$ represents the residuals, and $U_1\circ U_2 \cdots \circ U_R$ is the CP model [71]. Figure 6 is a schematic representation of the CP decomposition. CP decomposition is similar in form to Tucker decomposition since the number of components in the factor matrices is the same [96]. The original CP problem is NP-hard. Therefore, the Frobenius norm is generally used to relax the low-rank constraint. Similar to Tucker decomposition, CP decomposition problem can also be solved by alternating least squares, called CP-ALS. To the best of our knowledge, there is currently no well-known article applying CP decomposition to ultrasound clutter filtering.

Fig. 5

The illustration of Tucker decomposition

Fig. 6

The illustration of CANDECOMP/PARAFAC (CP) decomposition

Minimization problems

The decomposition problems generally turn into minimization problems or optimization problems in its original form or its Lagrangian form [63].

$$\begin{aligned} \min _{K_i} \, \sum _{i=1}^{x}{\lambda _if_i(K_i)}\quad \text {s.t.} \;C_i \end{aligned}$$

(18)

where $\lambda _i$ are the regularization parameters, $f_i(.)$ are the loss functions for low-rank, sparse, and noise components, $C_i$ are the constraints on $K_i$. Consistent with the decomposition problems, the minimization problems can be divided into three categories according to the number of constraints and loss functions imposed.

• $x=1$ is the case of implicit decomposition: $\min _L \, \lambda _1f_1(L) \quad \text {s.t.} \; C_1$

  where $C_1$ is $\Vert M-L\Vert _2=0$ or other forms. For sparse decomposition, the low-rank components are replaced by sparse components. This problem can be NP-hard, non-convex, or under specific constraints. Therefore, other formats of the loss functions are applied to relax the constraints when the problem is NP-hard or non-convex. For example, the loss function f is rank loss function in original MC model as $\min _{L}\, rank(L) \quad \text {s.t.} \; \Vert M-L\Vert _2=0$.

• $x=2$ is the case of explicit decomposition: $\min _{L,S} \; \lambda _1f_1(L)+\lambda _2f_2(S) \quad \text {s.t.} \; C_2$

  where $C_2$ is $\Vert M-L-S\Vert _2=0$ or other forms. For example, the $f_1$ and $f_2$ loss functions are rank and $l_0-norm$ loss functions in original RPCA model as $\min _{L,S}\, rank(L)+\lambda \Vert S\Vert _{l_0} \quad \text {s.t.} \;\Vert M-L-S\Vert _2=0$.

• $x=3$ is the case of stable decomposition: $\min _{L,S} \, \lambda _1f_1(L)+\lambda _2f_2(S)+\lambda _3f_3(N) \quad \text {s.t.} \; C_3$

  where $C_3$ is $\Vert M-L-S-E\Vert _2=0$ or other forms. For example, the $f_1$ and $f_2$ loss functions are rank and $l_0-norm$ in original RPCA model as $\min _{L,S}\, rank(L)+\lambda \Vert S\Vert _{l_0} \quad \text {s.t.} \; \Vert M-L-S\Vert _2=0$. The stable decomposition is generally adding constraints on the noise component based on the robust decomposition. The Frobenius norm loss function ($\lambda \Vert M-L-S\Vert _F^2=0$) is used in most cases.

Although there are some algorithms that can solve non-convex problems through mathematical approximation [100], in general, non-convex problems are difficult to solve with weak convergence. This is also an important role that minimization problems play.

Loss functions

The loss function can be seen as a constraint of the decomposed matrices. In DLSM framework, loss functions are used on the minimization matrices as norm formats. For example, in implicit decomposition, explicit decomposition, and stable decomposition, the functions f(S), f(L), f(E), represent the loss functions or norms on sparse component, low-rank component, and noise component, respectively. However, in most cases, the original loss function will be replaced by other forms of the loss function to optimize and solve the problem. The common loss function forms (or norm forms) can be listed as follows:

• $l_0 \,$ norm loss function $(\Vert M\Vert _{l_0})$ is the number of non-zero entries.

• $l_1 \,$norm loss function $(\Vert M\Vert _{l_1} = \sum _{i,j}|M_{i,j}|)$ is the Manhattan distance.

• $l_2 \,$norm loss function $(\Vert M\Vert _{l_2} = \sqrt{\sum _{i,j}M^2_{i,j}})$ is also called the Frobenius norm ($l_{F} \,$norm loss function $(\Vert M\Vert _{l_{F}} = \sqrt{\sum _{i,j}M^2_{i,j}})$ ).

• $l_{\infty } \,$norm loss function $(\Vert M\Vert _{l_{\infty }} = max_{i,j}|M_{i,j}|)$ is also called the max norm ($(\Vert M\Vert _{max} = max_{i,j}|M_{i,j}|)$).

• $l_{*} \,$norm loss function $(\Vert M\Vert _{l_{*}})$ is the sum of singular values.

Solvers

The models are solved by specific algorithms, which are called solvers in DLSM [63, 71] framework. Solvers are generally applied to the models after the minimization problem has been optimized and the loss function has been relaxed. Solvers can be broadly divided into regularization-based approaches and statistical-based approaches [101]. As for regularization approaches, the data matrices are regularized by convex surrogates with different features [63]. Typical regularization approaches include singular value thresholding (SVT) [102], accelerated proximal gradient (APG) [103], and augmented Lagrange multiplier (ALM) [83]. In terms of statistical-based approaches, prior distributions are used to capture low-rank or sparse properties and predict the joint distribution of the measured entries and unknown variables. Meanwhile, posterior distributions of the unknown variables can be approximated by Bayesian-based methods [63].

Although many solvers are proposed to solve the optimization problems under DLSM framework, most of the mainstream algorithms for ultrasound clutter suppression are based on SVD. SVD-based clutter suppression algorithms that are proposed and reviewed [19, 20, 22, 38, 43] based on traditional CFI before 2011. In these algorithms, SVD is used as one of the steps or iterations within many of the algorithms we evaluated. After 2015, with the rapid development of ultrasound technologies like plane-wave ultrasound, SVD was combined with ultrafast plane-wave imaging, which can provide a huge amount of data at a high frame rate, to improve the effectiveness of SVD and overcome the limitation of low frame rate [18, 51, 104, 105]. Due to the excellent performance of SVD on large datasets, SVD-based clutter suppression algorithms based on the plane-wave ultrasound have become a popular and mainstream research area. The SVD-based algorithms have been used in functional ultrasound [106, 107], super resolution ultrasound localization microscopy [104, 108] and high-sensitivity microvessel perfusion imaging [18, 51] due to its excellent performance in the ultrasound clutter suppression and the microvascular imaging [21].

Due to the obvious disadvantages of SVD, DLSM framework contains many approximate SVD and non-SVD algorithms for higher efficiency and lower computational cost, which have the potential for real-time ultrasound clutter suppression.

Experiments

DLSM framework has been successfully utilized to video surveillance, face recognition, texture modeling, video inpainting, audio separation, and latent semantic indexing [109]. However, only a few algorithms under DLSM framework have been applied to ultrasound clutter suppression. Herein, we apply DLSM algorithms as the clutter filter for CFI. To that end, we test if DLSM algorithms can be used for clutter suppression and conduct simulation experiments, phantom experiments, and in vivo experiments. Finally, we will conclude a list of algorithms that are suitable for ultrasound clutter suppression.

Experiment data

Three datasets are used in this experiment which are simulation data, phantom data, and in vivo rat data. For each dataset, raw RF-data, complex envelope data, and B-mode data formats are used for analysis. The specific parameters and obtaining process of three datasets and a brief introduction of three data formats are given in the following subsections.

Simulation data

The simulation data include a set of ultrasound simulation frames as Fig. 7 shows. The ultrasound simulation data are generated by the Field II simulation program implemented in MATLAB [110, 111]. A cube $A\in \mathfrak {R}^{60 \times 60 \times 60 }$ is built to represent the tissue filled with scatterers given the fact that each voxel is $1\, \text{mm}^3$. A vessel through the middle of the cube with a radius of $20\, \text{mm}$ is generated by scatterers flowing to the right. The max velocity in the center of the vessel is $15\, \text{mm}/\text{s}$. Assuming sound waves travel from the top to the bottom and focus on the center. Probe frequency and sampling frequency are set to $7.27\, \text{mHz}$ and $40\, \text{mHz}$, respectively. The frame rate is set to $1000\, \text{fps}$ and 64 active elements are used for beamforming.

Fig. 7

The illustration of the simulation data. a The simulation cube with tissue scatterers and blood scatterers. The red blood scatterers are in the middle and moving to the right. The simulated sound waves focus in the center. b A series of simulation data frames obtained from simulation experiments

Phantom data

The phantom was created to simulate a cube of tissue including one blood vessel which travels across the cube in the middle. Knox unflavored gelatin, water, and sugar-free Metamucil psyllium fiber supplement were gently heated and mixed to prepare the phantom gel which represents soft tissue. An intra-venous tube simulating a venous structure model runs through the gel cube. Probe frequency and sampling frequency are set to 10 MHz and 40 MHz, respectively. The Alpinion E-Cube R12 ultrasound system is used in ultrasound data collection with an L3-12H linear array probe. Figure 8 briefly illustrates the phantom experiment.

Fig. 8

The illustration of the phantom experiments. a The illustration of phantom data collection experiment. b The B-mode image of the first frame in phantom data

Rat data

The acquisition of the rat data was under the supervision of the Animal Care Facility of Concordia University. A 27-week-old Sprague–Dawley male rat was anesthetized for ultrasound scanning. The experiment followed the guidelines of the Canadian Council on Animal Care and was approved by the Animal Ethics Committee of Concordia University $(\#30000259)$. The probe frequency and the sampling frequency were set to 10 MHz and 40 MHz, respectively. Similarly, as with phantom data, the Alpinion E-Cube R12 research ultrasound system with an L3-12H linear array probe was used. The schematic diagram of the in vivo rat experiment is shown in Fig. 9.

Fig. 9

The illustration of the in vivo rat experiments. a The illustration of the in vivo rat data collection experiment. b A schematic representation of sparse component of the in vivo rat data

Data formats

Both real and simulated ultrasound data are available in three formats, which are raw RF data, complex envelope data, and B-mode data. Common ultrasound probes generally consist of a piezoelectric transducer array that emits and receives signals. The backscatter signal which is processed by the preamplifier and the time gain compensation is referred as radio-frequency (RF) signal. The RF signal then processed by an envelope detector becomes complex envelope data. Last, the complex envelope data are log-compressed into a grayscale format. And the data are further passed through intensity mapping and post-processing filtering. The final readable image is commonly called brightness mode (B-mode) image. RF frames generally have a very large size since the sampling rate of the RF data is usually extremely high. This high sampling rate is not necessary for envelope data as it does not have high-frequency contents. Therefore, envelope and B-mode images can be downsampled by a large factor. RF data may also be downsampled by a small factor, but the Nyquist sampling rate should be considered to avoid aliasing.

Experiment methods

In “Decomposition into low-rank and sparse matrices (DLSM) framework” section, DLSM framework is introduced and built as Fig. 3 shows. The DLSM algorithms are classified in four groups which are implicit decomposition, explicit decomposition, stable decomposition, and tensor decomposition. In this experiment, all algorithms are selected from LRSLibrary [61, 63, 71] which provides a group of low-rank and sparse matrix decomposition algorithms in moving object detection. In LRSLibrary, these algorithms are further subdivided into robust PCA (RPCA), subspace tracking (ST), matrix completion (MC), three-term decomposition (TTD), low-rank representation (LRR), nonnegative matrix factorization (NMF), nonnegative tensor factorization (NTF), and standard tensor decomposition (TD) according to the models. Due to the flexible conversion between models and their similar mathematical formulations, in this paper, TTD can be a subcategory in stable PCA under stable decomposition. Similarly, NTF belongs to the subcategory under the tensor decomposition (TD) model.

In the first step, the DLSM algorithms are applied to three formats of simulation data to verify the performance of all algorithms compared to sparse component with ground truth and give a computing time contrast. Then, all algorithms are used on phantom data to find out if DLSM suits ultrasound data with real ultrasound features. In the third step, rat data are used for verification and comparison. The acquired data have three formats which are RF data format, complex envelope data format, and B-mode data format. The results of different data formats and different datasets are grouped for comparison to find the optimal conditions of ultrasound clutter suppression.

All experiments are processed by a normal desktop computer with an i7-4770 CPU @ 3.40 GHz and 16.0 GB RAM.

Evaluation metrics

Two main indicators are used to evaluate the performance of various algorithms, which are signal-to-noise ratio (SNR) [112] and contrast-to-noise ratio (CNR). The SNR and CNR are calculated as follows:

$$\begin{aligned} \text{SNR}=\frac{\mu _1}{\sigma _1},\quad \text{CNR}=\frac{|\mu _1-\mu _2|}{\sqrt{\sigma _1^2+\sigma _2^2}/2} \end{aligned}$$

(19)

where $\mu _1$ and $\sigma _1$ are the mean intensity value and the standard deviation of the background window, $\mu _2$ and $\sigma _2$ are the mean intensity value and the standard deviation of the target window.

Experiment results

The results of 106 DLSM algorithms on three datasets and their three formats are reported in this section. The results of all algorithms include the SNR, CNR, calculation times, and images for visual observation. Since all the output images are sparse components of the same data and are very similar, we classify the results according to their performance and report the number of algorithms in each category instead of SNR and CNR of all algorithms.

The results of all algorithms are divided into several categories. The results which fall in the first category are considered to be good results as they give the correct sparse matrix with a pure blank background which means high robustness to noise and dynamic background and strong decomposition ability. The cases when the output sparse component is more than 100 times higher than background pixel values are also regarded as good results. The results which fall in the second category are considered to be defective. These results either contain background noise which is supposed to be part of the low-rank components, or are noisy and algorithms failed to decompose. The results in the third category are not considered because some algorithms failed to run due to some limitations like non-negative limitations or real input limitations. Algorithms with this type of results are called restricted algorithms in this section.

The information that all algorithms, including their model classifications, are from LRSLibrary [61, 63, 71], and they have all been proved to be successfully applied to moving object detection on traffic video.

Simulation experiments

The experiments first applied simulation data to verify the availability and approximate performance of all algorithms. Meanwhile, the computation cost and time of these algorithms on ultrasound clutter suppression are also tested. The first experiment applied all 106 DLSM algorithms to the RF simulation data. Among 106 algorithms, 11 of them were out of memory and failed to run. These algorithms cannot deal with the large size of simulation data because they use the full singular value decomposition or QR decomposition and require a huge memory to initialize (7.9 GB). Meanwhile, there are 6 algorithms that require non-negative input and cannot take RF data as input. Consequently, a total of 17 of these two kinds of algorithms are classified as restricted algorithms.

In terms of the remaining 89 algorithms, only 19 of them are able to output relatively pure sparse components that match the ground truth without any processes of RF simulation data. To be precise, only three algorithms (abbreviation: LRR-ROSL, RPCA-IALM, RPCA-IALM-BLWS) give a truly pure background as zero matrices (all entries in sparse matrices except the ones presenting simulated vessel are 0). The other 16 results highlight simulated vessel with a non-zero background. Since the value of the background pixels is 1000–10,000 times less than the value of sparse component, we consider it to be a pure result without low-rank components. The possible reason is the particular small values of RF simulation data and low dynamic range. In general, these results with the CNR values above 1.6 are classified as good results in Table 2. The results of the other 44 algorithms are very noisy with the CNR values less than 1.1. As for these algorithms, the sparse parts in simulation data are not clearly determined and the clutter is not well suppressed. The remaining 24 algorithms give blank output due to low dynamic range and other reasons. Almost all the DLSM algorithms give an SNR of about 0.759, so SNR is not reported in detail here.

Table 2 The 19 algorithms with the CNR values above 1.6

Due to the extremely small data values and dynamic ranges, a large number of algorithms are invalidated. Therefore, in the second step, the order of magnitude and dynamic range of RF simulation data are expanded to re-examine the performance of all algorithms.

After processing RF data, 56 algorithms show good results. Among these algorithms, 16 of them give correct sparse components with a zero-valued background, others give sparse components 1000–10,000 times greater than background pixel values. The results of the remaining 33 algorithms are noisy. These algorithms either do not correctly isolate sparse components or contain inseparable background noise with similar values. Most good results have a CNR greater than 1.3, while noisy results generally have a CNR less than 1. Similarly, almost all the DLSM algorithms with good results give an SNR of about 0.759. There are a few good results with a CNR less than 1. The algorithms with such results only highlight the sparsest parts which reduce the mean intensity values of the target window. However, these results are considered to be good because the unhighlighted sparse components still have higher intensities than backgrounds. The results after increasing dynamic range are listed in Table 3. Examples of different kinds of results in simulation exp

Table 3 The 16 algorithms with pure background after increasing dynamic range

eriments are shown in Fig. 10.

Fig. 10

The output result images of simulation data. a The output of sparse component obtained by the IALM algorithm on original simulated RF data. It is a typical good result representing correct decomposition and pure sparse components. b The output of sparse component obtained by the ADM algorithm on original simulated RF data. It is a typical noisy result with background noise as sparse components. c The output of sparse component obtained by the OSTD algorithm on processed simulated RF data with larger dynamic range. The algorithms with a CNR less than 1 in Table 3 give such results with pure background because they only show the sparsest parts

The complex envelope simulation data are obtained by Hilbert transform based on the RF data. For this reason, the complex envelope data do not have the problem of miniature pixel values and low dynamic range. However, the SNR and CNR of the complex envelope simulation data are lower than the SNR and CNR of RF simulation data. Except for the 11 algorithms that are limited by frame size, 24 of the remaining algorithms show good results. In addition, 13 algorithms are affected by complex numbers generated by the Hilbert transform and thus failed to run. The results of other algorithms are noisy. After enlarging the dynamic range of the complex envelope simulation data, 8 algorithms which failed on original simulation data give good results on the preprocessed data. These algorithms are sensitive to the changes of dynamic ranges. The results on complex envelope simulation data are shown in Table 4. Obviously, the CNR in the results of complex envelope simulation data is far less than the CNR on the RF data.

Table 4 The algorithms with good results on complex envelope simulation data

The third step of the simulation experiment is using B-mode data. As for the results of B-mode simulation data, 40 DLSM algorithms have successfully detected the simulate vessel on original B-mode simulation data. Meanwhile, 12 algorithms are affected by high peak values in the background and keep static peaks into sparse components. These algorithms give pure sparse matrices after suppressing peak values. After enlarging the dynamic range of the original B-mode data, another 10 algorithms successfully detect correct sparse components. Therefore, 62 algorithms can successfully separate the correct sparse components. The other algorithms which give very noisy results may need parameter adjustment and threshold process. The results of simulation experiment on B-mode data are reported in Table 5.

Table 5 The algorithms with good results on B-mode simulation data

Phantom experiments

The next set of experimental data used for testing is phantom data. The phantom data are used to test whether DLSM framework is suitable for ultrasound clutter suppression with real ultrasound noise and other ultrasound features. The phantom data also consist of three formats, which are RF phantom data, complex envelope phantom data, and B-mode phantom data. Except for 11 inapplicable algorithms due to size limitation, the other 95 algorithms have a huge difference in computing time ranging from less than 0.1 s to more than 500 s.

As for RF phantom data, the order of magnitude of all pixels is first adjusted into the range of $10^{\pm 3}$. However, the structured peak pixels that are caused by bright structure generated at the rebound reflection interface still affect many algorithms. 36 algorithms clearly show the simulated vessel with a pure background with an average CNR of 3.5. Meanwhile, 43 algorithms only highlight bright edges as the sparse components with an average CNR of 0.4. The structured peak pixels of RF phantom data can compromise the calculation of some algorithms when these bright edges have a pixel value $10^3$ times larger than the remaining pixel values. Therefore, the peak values are processed logarithmically to achieve the gray balance and reduce the dynamic range. After logarithmic processing, 22 additional algorithms are able to display sparse components correctly excluding bright and static edges. Among them, 19 algorithms are previously affected by peaks, and 3 algorithms are defective on the original data. The results of remaining algorithms are still noisy, and parameter adjustment should be applied to these algorithms for better performance. The results of RF phantom experiments are shown in Table 6.

Table 6 The algorithms with good results in RF phantom experiments

The complex envelope phantom data are then used for experiments. The number of good results of the complex envelope phantom data is less than the number of good results of RF phantom data. 26 algorithms successfully detected the simulated vessel and 33 algorithms only showed bright edges, which is an intra-venous (IV) tube representing the vessel wall. The CNR of all results is less than 0.4. After suppressing the edge brightness logarithmically, 11 of these algorithms that have been affected by edges can separate pure sparse components. It shows that the extremely high bright structures can affect the sensitivity to sparse components in many algorithms. However, there are still some algorithms that give noisy results. At the same time, 19 algorithms cannot take complex numbers as input. The results on complex envelope phantom data are shown in Table 7.

Table 7 The algorithms with good results on complex envelope phantom data

The last format of data to be applied is B-mode data. As for B-mode phantom data, 49 algorithms successfully give good results with a CNR higher than 2. However, the results of four of these algorithms contain bright edges which are considered to be low-rank components. The results of 35 algorithms only show the edges. Among them, a few sensitive algorithms can also partly detect sparse partition with obvious motion. However, only several bright pixels with motion can be detected. After reducing the dynamic range, no algorithm is affected by edges and 75 algorithms give sparse components with pure background. Examples of good results and noisy results in phantom experiments are shown in Fig. 11. The results of B-mode phantom data are shown in Table 8.

Fig. 11

Three typical output results of phantom experiments. a A typical good result showing pure sparse components without noise. This image is obtained by ALM algorithm on original phantom data. b A typical output affected by bright edge structures. This image is obtained by APG algorithm on original phantom data. Because the pixel values of bright edges are 1000 times larger than the pixel values in the rest of the image, the flow sparse component in the middle of the tube cannot be observed. c A typical noisy result showing sparse components with indivisible noise. This image is obtained by RSTD algorithm on original phantom data

Table 8 The algorithms with good results on B-mode phantom data

In vivo experiments

In the third step, rat data are used to test the performance of these algorithms on real ultrasound data with small vessels-like tissues. The RF rat data, complex envelope rat data, and B-mode rat data are used to be compared.

Fig. 12

The examples of the results of rat experiments. a The B-mode image of rat data for comparison. b is obtained by ALM algorithm on original rat data. The dynamic background and noise are filtered out relatively well. c is obtained by APG algorithm on original rat data. Large areas of dynamic tissue are classified as sparse components. Since there is no ground truth for in vivo rat data, the results are described using relatively good and relatively noisy

In terms of RF rat data, 82 algorithms give very good and similar results, and 7 algorithms show noisy and meaningless results. The other 17 algorithms are restricted due to the size limitation or non-negative limitation. The complex envelope rat data remains to share similar results with RF rat data. As for B-mode rat data, 92 algorithms successfully detected vessel-like tissues and only 3 algorithms failed to show any part of the sparse components. Examples of good results and noisy results in in vivo rat experiments are shown in Fig. 12. Since most algorithms give results with similar SNR and CNR, the evaluation of results combine subjective observations and numerical analysis. Due to the unknown in vivo structure, we lack ground truth for the accuracy of the assessment results. Only algorithms with pure backgrounds are shown in Table 9 due to similar results and limited space.

Table 9 The algorithms with pure backgrounds on in vivo data

Discussion

A total of 106 algorithms were tested in this paper. Analyzing the results obtained from simulation, phantom, and in vivo experiments, we found that 11 algorithms require huge memory (7.9 GB for frame size $250\times 125$, 20 frames) due to the singular value decomposition or QR decomposition process. Since typical ultrasound frames are large in size, the left unitary matrix in full singular value decomposition demands an excessive amount of memory, e.g., ADM. There are two possible solutions to this problem. First, the approximate SVD can be calculated and stored in every iteration instead of full SVD [55, 113]. Second, small overlapping patches from the ultrasound frames can be considered to formulate the data matrix which will substantially reduce the size of the Casorati matrix and eventually the memory footprint. Another advantage of using this windowing technique is that it can automatically equalize uneven noise distribution by normalizing the power locally [51]. The 11 algorithms with size limitation are listed in Table 10.

Table 10 The 11 algorithms with size limitation

Moreover, there are 6 algorithms that require non-negative input. Since ultrasound RF data usually contain both positive and negative values, these algorithms are not suitable for working with RF data for clutter suppression. 20 algorithms giving good results on RF simulation data are tested with the absolute value of RF data to confirm the impact of non-negative requirements on ultrasound clutter suppression. Although all of these 53 algorithms are still capable of showing show high-contrast vessel structures, the SNR obtained with absolute value RF data (0.81) is slightly greater than the original SNR (0.76) showing a significant increase of background noise in sparse components. At the same time, the CNR obtained with absolute value RF data (1.69) is slightly lower than the original CNR (1.73), which proves that nonnegative requirement has only limited effect on the accuracy of the DLSM decomposition. The 6 algorithms with size limitations are listed in Table 11.

Table 11 The algorithms with non-negative requirement

Another type of restricted algorithms is affected by complex inputs. From the experiment results, it is obvious that complex envelope data are not suitable for ultrasound clutter suppression since it takes longer calculation time and gives poor performance. Also, 13 algorithms are affected by complex value and cannot separate low-rank and sparse components well. Among them, 13 algorithms cannot take complex numbers as input, and some algorithms are stuck in a longer loop that requires more than 300 s. Algorithms which failed due to complex numbers are listed below in Table 12. In addition, the extremely small CNR obtained from the envelope data is only one-hundredth of the ones obtained from other datasets which indicates that envelope data are not suitable as an input form of ultrasound clutter suppression.

Table 12 The 13 algorithms that cannot take complex numbers as input

Among the remaining algorithms, 17 algorithms are easily affected by outliers. These algorithms cannot denoise the peak values when the dynamic range is roughly greater than $11.5\, \text{bits}$, which is the natural logarithm of difference between maximum and minimum. These 17 algorithms have performed well on pre-processed data and showed good results in simulation experiments and phantom experiments. However, they are not robust to outliers. In the simulation experiments, these algorithms divided the background peak pixels into sparse components, resulting in a noisy background. Similarly, they are not robust to large-shaped structured outliers and divide bright static edges into sparse components in phantom experiments. These 17 algorithms that are susceptible to outliers are listed in Table 13. Since complex envelope data are not suitable for ultrasound clutter suppression, the performance of the algorithms on complex envelope data has not been considered.

Table 13 The algorithms not robust to the outliers

In addition, a pure background (0 dB) is of great significance for vascular image segmentation and process and analysis of other medical images [46, 114]. However, this is a difficult goal due to the probe jitter, dynamic backgrounds, noise, shadows, and many other reasons. Therefore, only a few results have pure background on simulation data and phantom data. Furthermore, no result has pure background on in vivo rat data because of the complex tissue motions and the harsh conditions. Some algorithms have a strong ability dealing with these challenges and give pure backgrounds on simulation data and phantom data. These algorithms are listed in Table 14.

Table 14 The algorithms with the potential to give a pure background

Overall, in terms of calculation time, DLSM algorithms take the longest time to run complex envelope data in comparison with RF data and B-mode data. Due to its large amount of calculations, complex envelope data take twice as long as RF data do to run. This confirms again that complex envelope data are not suitable for ultrasound clutter suppression. Meanwhile, RF data require slightly less computation time than B-mode data. This may be caused by the extra information RF data contain. Meanwhile, we can find that DLSM algorithms use a slightly longer time on preprocessed data than on original data. However, the algorithms separate sparse components more accurately. Table 15 lists the average time taken by the fastest 20 algorithms on different datasets and different data formats.

Table 15 The average time taken by the fastest 20 algorithms

The experimental results prove that ultrasound data are very different from ordinary video surveillance frames. All DLSM algorithms can be successfully applied to surveillance images. However, some of them are not suitable for ultrasound data. There are quite a few algorithms that are not suitable for ultrasound RF data and complex envelope data, which may due to the complexity of the RF data and the complex space of complex envelope data. The simulation experiments prove that some algorithms are still not robust to ultrasonic clutter and are not sensitive to the data with overall small pixel values ($<10^{-3}$). As for these algorithms, the low-rank components of the results often contain inseparable background flicker, noise, and tiny motion. These algorithms have been listed in Table 13. Meanwhile, the phantom experiment results prove that some DLSM algorithms are not robust and stable with a high dynamic range greater than $10\, \text{bits}$. For ultrasound data, an area with small values often exists in a uniform tissue. Edges that are much brighter than other tissues are also common due to the strong reflections at the interface. The CNR after preprocessing the ultrasound data is generally higher than the CNR of raw data. The result of the data that removed the peak is also significantly better than the results of raw data. Therefore, it is necessary to preprocess the ultrasound image when applying the DLSM algorithm. Moreover, parameter adjustment or other math improvements are necessary when applying some DLSM algorithms on ultrasound data to get the best filtering performance.

On the other hand, in terms of ultrasound data formats, experiments show that B-mode ultrasound data can make more algorithms successful for vascular detection. The B-mode ultrasound data may lose information. However, the outliers that may affect DLSM algorithms may also be weakened by Hilbert transform an absolute process. This might be the reason why more DLSM algorithms work for B-mode data. Although B-mode data has more good results than RF data, RF data require slightly less average calculation time and are more suitable for real-time requirements. The algorithms in Table 16 are relatively stable in all three datasets. These algorithms all require less than 1 s for computation while giving the correct sparse components. Experiments show that they may be more suitable for ultrasound clutter suppression.

Table 16 The algorithms require less than 1 s calculation time

Finally, 22 algorithms which are most robust to noise with the best performance in all the experiments are listed in Table 17. These algorithms may have a strong ability for ultrasound clutter suppression.

Table 17 The most robust algorithms with the best performance

In this paper, we adapted different techniques originally proposed for natural images in the field of computer vision for ultrasound color flow imaging. As ultrasound images have unique characteristics due to the physics of sound propagation, these images have the so-called “speckle noise”. We believe that the results of this paper can be generalized to other imaging modalities that are affected by diffraction, such as optical coherence tomography (OCT).

Conclusion

The performance of 106 established low-rank and sparse decomposition algorithms for clutter filtering has been tested in this work. Our results show that few robust matrix decomposition techniques are suitable for solving the limitations of SVD-based ultrasound clutter suppression methods such as sensitivity to large noise. In addition, several matrix decomposition techniques show the potential for real-time implementation on commercial ultrasound machines due to their low computational complexity. Furthermore, some preprocessing is necessary when applying this framework to ultrasound data. Finally, some of the algorithms studied in this work can automatically estimate the optimal power Doppler images without requiring extensive manual tuning, which may pave the way for easier commercial and clinical translation of ultrasound clutter suppression (Additional file 1).