Structurally random Fourier domain compressive sampling and frequency domain beamforming for ultrasound imaging

Advances in ultrasound technology have fueled the emergence of Point-Of-Care Ultrasound (PoCU) imaging, including improved ease-of-use, superior image quality, and lower cost ultrasound. One of the approaches that can make the adoption of PoCU universal is to make the data acquisition module as simple as a "stethoscope" while further processing and image construction can be done using cloud-based processors. Toward this goal, we use Structurally Random Matrices (SRM) for compressive sensing of ultrasound data, Fourier sparsifying matrix for recovery in 1D, and frequency domain approach for 2D ultrasound image reconstruction. This approach is demonstrated through wire phantom and in vivo carotid arteries data from ultrasound system using 25%, 12.5%, and 6.25% of the full data rate and ultrasound images of similar perceived quality quantified by Structural Similarity Index Metric (SSIM).


INTRODUCTION
With a portable ultrasound device, trained clinicians can identify more than forty medical conditions with a high degree of accuracy in two minuets [1]. They can detect collapsed longs, gallstones, blood clots, heart problems and blockages in the stomach, intestines, and kidneys. Using a portable high resolution ultrasound system, medics and firefighters can make critical patient care decisions in the "golden hour"-the important 60-minutes window after Traumatic Brain Injury when effective medical treatment is most crucial. Point-Of-Care Ultrasound (PoCU) Imaging is an emerging technology that provides accessibility to fast and accurate imaging system when it is needed. The goal in these systems is to move the hardware and software complexity from the data acquisition module to cloud-based processing units where the final image can be displayed on a hand-held device of a physician, for example.
Compressed sensing (CS) methods give us a fresh look at data acquisition, in general. The theory says one can design sensing or sampling protocols that captures the useful information content embedded in a sparse signal and condense it into a small amount of data [2]. It is a very simple and efficient data acquisition method which acquires data at a rate lower than the data needed in classical Nyquist sampling framework and uses computational power in the decoder from what it appears to be an incomplete set of measurements. The key ingredient to this recepie is sparsityhow much small coefficients can be ignored in a signal when compressed in a known basis. Sparsity is a fundamental modeling tool which allows accurate statistical estimation and classification, efficient data compression, and so on [2,3]. In the context of CS, sparsity has significant influence on the data acquisition process; it determines if one can acquire signals nonadaptively and efficiently. In recent years, the area of CS has branched out to a number of new fronts and has worked its way into several application areas, such as radar, communications, and ultrasound imaging. Eldar and her colleagues published extensively in applying CS theory to ultrasound imaging [4][5][6][7]. The key idea in all of these works is that the ultrasound signal can be modeled as Finite Rate of Innovation (FRI) signals which is a sum of a weighted and delayed known signal [8]. Filtered through a sampling kernel, the unknown parameters (delays and weights) can be estimated from the observations at much lower sampling rate as dictated by the Nyquist-Shannon theorem. Reconstruction of the original signals from the sub-Nyquist signals are carried out by using high resolution spectral estimation techniques such as the annihilating filter [9], matrix-pencil method [10], minimumnorm method [11], and other variants. Other papers follow general CS approach without any model [12,13]. The wave atom decomposition [14] is shown to have significantly sparse solution in ultrasound. Wave atom for ultrasound has been recently studied in [15] but the method was tested using simulations and not experimental data. In [16], we used the wave atom basis, a Bernouli random sensing matrix augmented with Robust Capon Beamformer to reconstruct ultrasound image of in-vivo cardiac data. When it comes to practical implementation, random sensing matrices require significant memory buffering for storage of matrix elements and high computational complexity due to their unstructured nature [2]. Also, with wave atom basis, there is no fast and efficient computation method like Fast Fourier Transform (FFT) to get the multiplication of sensing and sparsifying matrices as needed in recovery. Thus, the computational efficiency of the system in practice heavily depends on the structure of these matrices.
In this work, we use Structurally Random Matrices (SRM) for sensing, the Fourier basis for sparsifying matrix, a regularized-1 optimization technique for signal reconstruction in 1D and a frequency domain beamforming method for 2D image reconstruction. We use the SRM sensing matrices due to their low complexity for real time ultrasound imaging, hardware implementation simplicity, and optimal or near optimal performance for recovery [17]. The Fourier basis is used for signal reconstruction using the FFT or other fast form Fourier implementations (e.g., Sparse Fast Fourier Transform [18]). The optimization problem based on regularized-`1-norm is solved using the YALL1 package [19] from Rice university to reconstruct the signals in 1D. Finally, the CS-based beamforming is applied in the frequency domain to get the image of a wire phantom as well as in vivo carotid arteries image using the Verasonics Vantage system.

SYSTEM MODEL FOR COMPRESSIVE SENSING
We use a small bold letter for a vector, a capital bold letter for a matrix, a capital bold letter with subindex for a row or column of a matrix, and a capital letter with two subindices for an entry of a matrix. It is assumed that the radio frequency (RF) signals received by the i'th transducer (for 1  i  L) in our setup are sampled using an (M ⇥ N ) (with M ⌧ N ) sensing matrix as follows: The signal xi = (Xi1, ..., XiM ) T 2 C M is the noisy observation vector, yi = (Yi1, ..., YiN ) T 2 C N is the signal to recover, and vi = (Vi1, ..., ViM ) T 2 C M in (1) represents any noise including instrumentation errors and errors due to the signal being only approximately sparse. In CS, the received signals can be reconstructed from a relatively small number of linear incoherent measurements, given the fact that the signals have sparse representation in a known basis (i.e., yi = zi ) which is incoherent with the measurement matrix . The transform coefficient vector zi is compressible (or K-sparse) in domain. It has been shown that when the sensing matrix is orthonormal, the minimum number of measurements needed to recover the overwhelming majority of yi is given by [2]: where C is a constant, K is the number of nonzero components in zi, and the mutual coherence µ( , ) is the largest entry in : (2), the number of measurements required for recovery is proportional to the coherency between the sensing and the sparsifying matrices [2].
The incoherency between the sensing and sparsifying matrices does not guarantee a good recovery. The other important condition is the stochastic independence of the compressed measurements [20], which in case of sub-Gaussian matrices, rows of the random matrix are needed to generate this stochastic independecy. There are various reconstruction techniques based on`0and`1-minimization methods.`0-minimization techniques are NP-hard while`1-based techniques can be solved in polynomial time using convex optimization methods [21]. For a certain sensing system, the stability of sparse approximation techniques must be evaluated in the presence of measurement noise and imprecise sparsity. For this purpose, one of the most accepted conditions is the Restricted Isometry Property (RIP). RIP suggests that the K-sparse vectors are not in the null space of the sensing matrix, hence, the vectors can be reconstructed successfully. A matrix has the RIP of order k and constant k 2 (0, 1) if it satisfies the following condition.
In other words, matrices satisfy the RIP condition can be constructed using random distributions such as Gaussian or Bernoulli as their entries [3]. Then, yi can be recovered from xi with high probability if the elements of are independent realizations of a Gaussian random distribution or are following a Bernoulli distribution of ±1.

RANDOM BASED SENSING MATRICES
To date, the proposed random based sensing matrices fall into three types of scenarios [2,17,22].
• Random matrices from a sub-Gaussian distribution: The sensing matrix is generated by drawing each element independently from a sub-Gaussian distribution such as Gaussian or Bernoulli [22][23][24]. This family of random matrices is incoherent with all other recovery dictionaries, therefore, there is no need to have a prior knowledge about the sparsifying matrix; making this solution universal. Also, the number of measurements needed to reconstruct the signal is optimum as However the drawback of using random matrices in practical applications is that they need large memory buffering for storage of the matrix elements and have high computational complexity due to their unstructured nature.
• Random matrices from an orthonormal basis: In this scenario, the measurements are taken from rows of an orthogonal matrix. This selection is done uniformly at random from among all subsets of rows of size M . One example is Fourier matrix (or partial FFT) as given in [25].
Then, the minimum number of samples required for exact reconstruction of signal depends on the coherency between the sensing and sparsifying matrices as shown in (2). It has a fast and efficient implementation which reduces the complexity both at encoder and decoder, but the sparsifying basis is limited to the identity matrix.
• Structurally random matrices: The SRM uses an orthogonal matrix whose columns are permuted randomly or the sign of its entries in each column are reversed simultaneously with the same probability [20]. Three steps are involved in making the SRM matrices: (i) Prerandomization, which randomizes the signal either by uniformly permuting the locations or flipping its sample sign (matrix R 2 C N ); (ii) Transformation, which applies an FFT (matrix F 2 C N ) to the signals from step (i); and (iii) Subsampling, which randomly picks up M rows out of N from matrix of step (ii) (matrix D 2 C M ). In summary, the SRM is defined as a product of these matrices with a normalization of energy as [17] = r N M DFR, where R is a uniform permutation matrix which scrambles the signal's sample locations globally, F is the FFT matrix as given in (5), and D is simply a random subset of M rows of the FR matrix. In this paper, we use the SRM sensing matrices due to their low complexity for real time ultrasound imaging, hardware implementation simplicity, and optimal or near optimal performance for recovery.

RECONSTRUCTION TECHNIQUES IN 1D AND 2D
In this section, we introduce our signal recovery method as well as frequency domain beamforming approach to produce the ultrasound image. First, we review the recovery method we used based on regularized-`1 [26], and then the Delay-and-Sum frequency domain beamforming to reconstruct the image.

Fourier Domain Signal Reconstruction
To recover the signal at the decoder, a non-linear reconstruction algorithm such as basis pursuit, orthogonal matching pursuit, iterative thresholding with projection onto convex sets and lots of variants are proposed in the literature [27,28]. The quality of recovery of the RF signals depends mostly on (i) choice of the sensing matrix (ii) choice of the best basis in which the signals have the most sparse representation, (iii) the incoherency between the sensing and the basis in which the signals are sparse ; and (iv) the ratio of the number of compressed measurements acquired to the number of information bearing (non-zero) components of the signal in that basis. Recently, Demanent and Ying [14] showed that wrapped oscillatory patterns like ultrasound waves have sparse representation in the wave atom basis. In our previous work [16], we followed [15,29] and used wave atoms as the sparsifying basis as shown in the third block of Fig.1 of [16]. The problem with wave atom transform is that there is no fast and efficient computation method like the FFT to get the multiplication as needed in recovery. Thus, the computational efficiency of the system in practice heavily depends on the structure of these matrices.
In this work, we use a Fourier dictionary ( ) in the recovery phase due to two reasons: (i) Existence of fast and efficient implementation of this dictionary; and (ii) No need to transform the signals to time domain as long as the beamforming can be implemented in the frequency domain as well. The issue with using this sparsifying matrix is the coherency with the sensing matrix we used. But given the random structure of the sensing matrix , we expect that this will not have a significant effect with experimentally acquired ultrasound data and the performance in practical situations. In recovery, instead of using the traditional 1-norm or basis pursuit, the optimization problem used is based on regularized-`1 [16,26] which can be written as with ⌧ > 0 being the regularization parameter. The`1regularization term in (7) promotes -domain sparsity in the solution and the`2-norm term keeps the solution close to the measurements. The YALL1 package [19] from Rice University is used to solve the optimization in (7) suitable for ultrasound signals with high dynamic range. The regularization parameter ⌧ is selected empirically based on a trade off between the quality of reconstruction and the speed of convergence. We have done a sensitivity analysis to see the effect of this parameter with all the other ones defined in the YALL1 package and select the best values to recover the signals.

Frequency Domain Beamforming
In the frequency domain, the recordings at transducers are denoted by Yi(!q), (1  i  L), for L transducers and Q frequency bins (1  q  Q) which can be presented as [16] Yi(!q) = Hi(!q)P (ri, !q) + Ni(!q), where Hi(!q) is the frequency response of transducer i, P (ri, !q) is the pressure field at the receiver i in response to the transmitting wave, and Ni(!q) is the observation noise. In this model, the receivers are assumed to be point transducers. In practice, however, due to the limited size of the transducer, Hi(!q) is also space dependent. The pressure P (ri, !q) is the multiplication of frequency domain Green's function of the medium and the source field generated from the scatterers at location rs denoted by F (rs, !q). Defining the frequency dependent near field array steering vector a(rs, !q) as a collection of Green's functions of the medium for the receive array as a(rs, !q) , [ e j(!q /c)|rs r 1 | 4⇡|rs r1| · · · , e j(!q /c)|rs r L | 4⇡|rs rL| with c being the speed of sound, the (L ⇥ 1) vector of received signals can be represented as with y(!q) = [Y1(!q), · · · , YL(!q)] T . Similarly, the vector n(!q) is the stack vector of all the L observation noise Ni(!q), and the notation is the Hadamard product of the two vectors. The (L⇥1) vector h(!q) = [↵(r1, !q)H1(!q), · · · , ↵(rL, !q)HL!q)] T takes care of both the frequency dependent attenuation factors ↵(ri, !q) of the medium as well as the frequency response of the transducers. For any observation point x inside the ROI, the search steering vector is denoted by a(x, !q). Before applying the beamforming method, we use a wavelet filter to suppress the noise in the received signals. We use the Matlab wavelet denoising tool based on non-parametric function estimation to efficiently filter the signals. To be consistent with the original signals (without compressive sensing), we use the same tool to filter SSIM out the noise. Withỹ(!q) being the filtered vector of y(!q), the non-coherent DAS beamforming results in the following image reconstructed from the recorded signals The following section summarizes our experimentation with an ultrasound imaging device to compare the image reconstructed from frequency domain compressive sampling and beamforming with the original data reconstructed by our research based ultrasound system.

EXPERIMENTATION
This section explains our experimental setup to study the proposed approach with SRM sensing matrix at the encoder, the Fourier transform as the sparsifying basis, and frequency domain beamforming as our image reconstruction method. Our database includes ultrasound data from ultrasound phantoms (wire and breast phantoms), carotid arteries, bladder, liver, breast, spleen, kidney, and liver. The ultrasound system used to acquire data is Verasonics (http://verasonics.com). The Vantage system in our lab, uses state-of-the-art hardware and software technologies to provide access to raw ultrasound data from each channel, while preserving the ability to perform high quality real-time imaging at clinically useful frame rates. The L7-4v, a 128-element linear array transducer with pitch of 0.298(mm), sensitivity of -64.5(dB), and center frequency of 5MHz is used to acquire the data. In this paper, we present two sets of data on our wire phantom, as well as on in vivo carotid arteries data. These data are acquired in ideal mode which means that each transducer sends a pulse sequentially and the received data is acquired by all the transducers, producing a (128 ⇥ 128 ⇥ number of samples) data matrix. Fig. 1 (a) shows the original image at full sampling frequency of 10MHz acquired by the system. The subsequent figures show our results after using 25%, 12.5%, and 6.25% of the original number of samples acquired by the system. The SRM sensing matrices is used based on uniformly permuting the locations of a Fourier matrix and subsampling to M measurements with ratio M/N given before. The regularized-`1 solver using the YALL1 package with ⇢ = 0.001, ⌧ = 0.5 ⇤ ⇢ is used to recover the signals in frequency domain. Before applying image reconstruction method, we use the Wavelet Toolbox TM of Mathworks to suppress the noise. For quantitative comparisons, we use the structural similarity (SSIM) index [30] in the resulted image to measure similarity between two images [5]. We present the SSIM for subplots (b-d) in Fig. 1, in table 1. Since a higher value of SSIM implies a greater level of similarity between the images, using CS results in some degradation in the quality of the image reconstruction but with such limited computation in the encoder, this degradation is acceptable. Finally, the same method is applied to the carotid arteries data obtained from the ultrasound imaging device and presented in Fig. 2. Similarly, Fig. 2 (b-d) show the image reconstructed from the data subsampled again at M/N = 25%, 12.5%, 6.25% respectively. Visual observations show that most of tissue boundaries data including the arteries and walls of the vessels are preserved even after removal of almost 83% of the data. However, 1/16 sample rate reduction

SUMMARY
In this work, we applied the SRM compressive sensing matrices to experimental ultrasound data from a Verasonics system to reduce the number of samples acquired in current ultrasound machines. The Fourier transform was used as the sparsifying matrix followed by the frequency domain beamforming. The results showed that after intensively reducing the sampling rate by the SRM matrices and recovering the signals using FFT, we can get similar image quality as the one reconstructed from the original samples at full data rate.