A Low-Rank and Joint-Sparse Model for Ultrasound Signal Reconstruction

TL;DR

This paper introduces a low-rank and joint-sparse model for ultrasound signal reconstruction that effectively reduces data requirements while maintaining high image quality, addressing data transfer and storage bottlenecks in dense sensor arrays.

Contribution

The paper proposes a novel low-rank and joint-sparse modeling approach that leverages inter-channel correlations for efficient ultrasound signal reconstruction from limited samples.

Findings

Can recover high-quality images from 10% of the samples

Effective exploitation of inter-channel correlations

Reduces data transfer and storage needs

Abstract

With the introduction of very dense sensor arrays in ultrasound (US) imaging, data transfer rate and data storage became a bottleneck in ultrasound system design. To reduce the amount of sampled channel data, we propose to use a low-rank and joint-sparse model to represent US signals and exploit the correlations between adjacent receiving channels. Results show that the proposed method is adapted to the ultrasound signals and can recover high quality image approximations from as low as 10% of the samples.

Equations

$$\text{X}=\text{YD}$$

$$\text{rank}(\text{X})\le k\text{when}k<N$$

$$\text{B}={\text{P}}_{\Omega }(\text{X})+{\text{N}}_e$$

$$\hat{\text{X}}=arg\underset{\text{X}}{\min }{\left|\left|\text{X}\right|\right|}_{\ast }+\alpha {\left|\left|{\text{Y}}_t\text{X}\right|\right|}_{2,1}+\frac{1}{2\mu }{\left|\left|\text{B}-{\text{P}}_{\Omega }(\text{X})\right|\right|}_F^2$$

$$\hat{\text{D}}=arg\underset{\text{D}}{\min }{\left|\left|\text{D}\right|\right|}_{\ast }+\alpha {\left|\left|\text{D}\right|\right|}_{2,1}+\frac{1}{2\mu }{\left|\left|\text{B}-{\text{P}}_{\Omega }(\text{YD})\right|\right|}_F^2$$

$${\text{D}}^{s+1}=\arg \underset{\text{D}}{\min }\frac{1}{2\gamma }{\left|\left|\left(\begin{array}{c}{\text{b}}_1^s\\ {\text{b}}_2^s\\ {\text{b}}_3^s\end{array}\right)+\left(\begin{array}{c}\text{I}\\ \text{I}\\ \text{Y}\end{array}\right)\text{D}-\left(\begin{array}{c}{\text{w}}_1^s\\ {\text{w}}_2^s\\ {\text{w}}_3^s\end{array}\right)\right|\right|}_F^2$$

$$\begin{gathered} \left(\begin{array}{c}{\text{w}}_1^{s+1}\\ {\text{w}}_2^{s+1}\\ {\text{w}}_3^{s+1}\end{array}\right)=\arg \underset{{\text{w}}_1,{\text{w}}_2,{\text{w}}_3}{\min }\{\frac{1}{2\gamma }||\left(\begin{array}{c}{\text{b}}_1^s\\ {\text{b}}_2^s\\ {\text{b}}_3^s\end{array}\right)+\left(\begin{array}{c}\text{I}\\ \text{I}\\ \text{Y}\end{array}\right){\text{D}}^{s+1} \\ {-\left(\begin{array}{c}{\text{w}}_1\\ {\text{w}}_2\\ {\text{w}}_3\end{array}\right)||}_F^2+\sum _{i=1}^3{g}_i({\text{w}}_i)\} \end{gathered}$$

$$\left(\begin{array}{c}{\text{b}}_1^{s+1}\\ {\text{b}}_2^{s+1}\\ {\text{b}}_3^{s+1}\end{array}\right)=\left(\begin{array}{c}{\text{b}}_1^s\\ {\text{b}}_2^s\\ {\text{b}}_3^s\end{array}\right)+\left(\begin{array}{c}\text{I}\\ \text{I}\\ \text{Y}\end{array}\right){\text{D}}^{s+1}-\left(\begin{array}{c}{\text{w}}_1^{s+1}\\ {\text{w}}_2^{s+1}\\ {\text{w}}_3^{s+1}\end{array}\right)$$

Full text

1 Introduction

Ultrasound echography is one of the most used diagnostic imaging techniques as it is real-time, safe, low-cost and portable. Conventional ultrasound imaging is usually performed by scanning a medium using sequential focused beams, each firing allowing the reconstruction of one line of the final image, i.e. Single-Line-Transmission (SLT) imaging. A high-resolution image requires sufficient numbers of transmissions as well as a sampling rate that is significantly higher than the Nyquist rate of the signal [1]. Consequently, with such high sampling rate, and taking into account the number of transmissions and the number of transducer elements, the amount of sampled data can become enormous, which restricts the data storage and transportation in most of the commercial systems today. In this context, several strategies based on compressed sensing (CS) have been proposed [2]. However, the reconstruction accuracy of CS-based methods is highly dependent on the basis and the number of measurements for each channel cannot be lower than the sparsity thereby limiting the data reduction rate. The aim of the current study was therefore to further reduce data rates by exploiting the low-rank property of US signals.

2 Methodology

Let us rearrange all the received pre-beamformed data in a 2D matrix $\textbf{X}\in\,\mathbb{R}^{M\times\,N}$. According to the fact that the pre-beamformed US signals from different transducer elements are joint-sparse in the Fourier domain  [3], X can be expressed in a matrix form as:

[TABLE]

where $\textbf{Y}\in\,\mathbb{R}^{M\times\,k}$ is a partial 1D Fourier matrix with $k(k\ll M)$ frequencies inside the bandwidth. This hypothesis is reasonable, since the US radio frequency (RF) signals are bandlimited by the impulse response of the transducer. Assuming that the signal bandwidth is 1, the $k=M(2f_{c}/f_{s})$ frequencies are the frequencies inside $[f_{c}/2,3f_{c}/2]$ and $[-f_{c}/2,-3f_{c}/2]$, where $f_{c}$ and $f_{s}$ are the center and sampling frequency of the RF signal. Thus $\textbf{D}\in\,\mathbb{R}^{k\times\,N}$ is the corresponding 1D Fourier coefficient matrix of the signal at those $k$ frequencies. In addition, since the maximum bandwidth of the RF signal is 1, the real joint-sparsity $K_{k}$ of the signal in Fourier domain is no more than $k$, i.e. $K_{k}\leq k$.

The factorization form in (1) implies a low-rank structure of X with:

[TABLE]

where $\text{rank}(\textbf{X})=k$ if and only if D is full-rank. In practice, thanks to the strong correlation between US signals, we have $rank(\textbf{X})\ll k$. Thus, X has a low-rank and joint-sparse structure when we have enough number of US channel signals, i.e. $N>\text{rank}(\textbf{X})$. Fortunately, in the field of US, the above condition is generally satisfied with $N\gg\,\text{rank}(\textbf{X})$, implying that the low-rank and joint-sparse property can be applied in US signal reconstruction.

Assuming that the data is contaminated by random noise, the measurements B can be modeled as :

[TABLE]

where $\Omega$ is a set of locations where the signal $\textbf{X}(m,n)$ is observed, i.e. $\textbf{X}(m,n)$ known if $(m,n)\in\Omega$. $\textbf{P}_{\Omega}(\textbf{X})$ represents the corresponding values of X in the locations of $\Omega$ and $\textbf{N}_{e}$ is the additive noise term. In these settings, the reconstruction problem thus amounts to solve (3) for X, under the constraint that X is low-rank and joint-sparse. This problem could be reformulated as an unconstrained optimization problem as in  [4]:

[TABLE]

where $\left|\left|\textbf{X}\right|\right|_{\ast}=\sum_{i}\sigma_{i}$ is the sum of the singular values (i.e. the nuclear norm) that aims at imposing the low-rank property of the RF signal X. This assumption has been extensively used in matrix completion, see [5, 6]; $\textbf{Y}_{t}$ is the adjoint operator of Y with the relation $\textbf{Y}_{t}\textbf{X}=\textbf{D}$. In our case, $\textbf{Y}_{t}$ and Y are the Fourier and Inverse Fourier matrix with effective frequencies $f_{i}$; $\left|\left|\textbf{D}\right|\right|_{2,1}=\sum_{q=1}^{k}\left|\left|\textbf{D}^{q\rightarrow}\right|\right|_{2}$ (i.e. the $l_{2,1}$ norm) that is used to explore the joint sparsity property of D and $q\rightarrow$ denotes the $q$-th row; $\left|\left|\cdot\right|\right|_{F}$ is the Frobenius norm. The parameters $\alpha$ and $\mu$ give the trade-off among the nuclear norm term $\left|\left|\textbf{X}\right|\right|_{\ast}$, the $l_{2,1}$ norm term $\left|\left|\textbf{Y}_{t}\textbf{X}\right|\right|_{2,1}$ and the data consistency term $\left|\left|\textbf{B}-\textbf{P}_{\Omega}(\textbf{X})\right|\right|_{F}^{2}$. Fundamentally, we are looking for a matrix with minimum rank and joint-sparsity subject to the acquired data. Since X has the same rank as D, this model could be reformulated as:

[TABLE]

The objective function in (5) aims at estimating D instead of X directly from the acquired samples. It is worth to notice that D is a $k\times\,N$ matrix with $k\ll M$, which means the number of variables to be estimated in D is $\frac{M}{k}$ times less than X and the problem (4) is much simplified. To solve the optimization problem in (5), we proposed in the following section an algorithm based on simultaneous direction method of multipliers (SDMM) [7].

3 Algorithm

In this section, we adapted the SDMM optimization framework to solve the problem in (5), which iteratively solves the above optimization problem as follows:

Step 1: Updating D, given $\textbf{b}_{i}$ and $\textbf{w}_{i}$:

[TABLE]

where $\textbf{w}_{1}=\textbf{w}_{2}=\textbf{D}$ and $\textbf{w}_{3}=\textbf{YD}$ are setting to indicate the three convex objective functions, $\textbf{b}_{1}$, $\textbf{b}_{2}$ and $\textbf{b}_{3}$ are Lagrangian parameters that have the same dimensions as $\textbf{w}_{1}$, $\textbf{w}_{2}$ and $\textbf{w}_{3}$, respectively. $(\cdot)^{s}$ means the updated value of $(\cdot)$ from the $s$-th iteration. $\gamma>0$ is a penalty parameter. As stated in [7], (6) is a classical $l_{2}$ norm minimization problem and can be efficiently solved.

Step 2: Solving $\textbf{w}_{i}$ using $\textbf{b}_{i}$ and D:

[TABLE]

where $g_{1}(\textbf{w}_{1})=\left|\left|\textbf{w}_{1}\right|\right|_{\ast}$, $g_{2}(\textbf{w}_{2})=\alpha\left|\left|\textbf{w}_{2}\right|\right|_{2,1}$ and $g_{3}(\textbf{w}_{3})=\frac{1}{2\mu}\left|\left|\textbf{B}-\textbf{P}_{\Omega}(\textbf{w}_{3})\right|\right|_{F}^{2}$ are the three convex terms of the optimization problem. Due to the separate structure of (7), it can be solved by minimizing the three subproblems that corresponding to the update of $\textbf{w}_{1}$, $\textbf{w}_{2}$ and $\textbf{w}_{3}$, respectively.

Step 3: Updating $\textbf{b}_{i}$, given D and $\textbf{w}_{i}$:

[TABLE]

4 Results

To validate the proposed methodology, in vivo cardiac pre-beamformed RF data ($f_{c}$-3.5 MHz; $f_{s}$-25 MHz, matrix size $4838\times 4544$) was captured using an experimental ultrasound system and reconstructed from a subset of samples obtained at different sampling rates (SR). The parameters in the optimization algorithm were chosen as $\{\gamma,\alpha,\mu\}=\{1,0.01,1e-6\}$ by cross validation. The contrast-to-noise-ratio (CNR) was calculated on the 8-bit grayscale B-mode image.

Figure 1 and Figure 2 show the B-mode images reconstructed for different sampling rates together with the corresponding CNR. Strong artifacts appear for the reconstructed image with a sampling rate of 5%, corresponding to a reduction in CNR of 2.7 dB. However, when the sampling rate increases, the artifacts disappear and the CNR improves, resulting in images visually very similar to the reference image. Image quality remained good even with only 10% of the samples. In comparison, given the sparsity of the problem (i.e. $2f_{c}/f_{s}$), using solely the sparsity promoting prior maximally achieve a SR of 28%.

5 Conclusion

This paper introduces a new data model based on the low-rank and joint-sparse priors. This regularization strategy allows for estimating missing data and thus reducing drastically the sampling rate. In vivo experiment was performed to validate and evaluate the proposed method. The results demonstrate that the proposed approach is capable of reconstructing the whole image from a sparse set of samples (e.g. 10% of samples) while keeping adequate image quality. However, the reconstruction algorithm is time consuming, thus real-time imaging is not realistic. Fast algorithm will be devised in future work.