{'fraction_non_alphanumeric': 0.03165511545229289, 'fraction_numerical': 0.008025999209968758, 'mean_word_length': 4.846016584444211, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Compressed Ultrasound Imaging:from Sub-Nyquist Rates to Super-Resolution

The multi-billion dollar, worldwide medical ultrasound (US) market continues to grow annually. Its non-ionizing nature, real-time capabilities and relatively low cost, compared to other imaging modalities, have led to significant applications in many different fields, including cardiology, angiology, obstetrics and emergency medicine. Facilitated by ongoing innovations, US continues to change rules and norms regarding patient screening, diagnosis and surgery. This huge and promising market is constantly driven by new imaging and processing techniques. From 3D images to sophisticated software, hardware and portability improvements, it is clear that the status of US as one of the leading medical imaging technologies is ensured for many years ahead. However, as imaging systems evolve, new engineering challenges emerge. Acquisition, transmission and processing of huge amounts of data are common for all ultrasound-based imaging modalities. Moreover, achieving higher resolution is constantly on demand, as improved diagnosis could be achieved by better visualization of organs and blood vessels deep within tissues. In this article, our goal is to motivate further interest and research in emerging processing techniques, as well as their applications in medical ultrasound, enabled by recent advancements in signal processing algorithms and deep learning. We address some of the primary challenges and potential remedies from a signal processing perspective, by exploiting the inherent structure of the received US signal.


I. INTRODUCTION
The multi-billion dollar, worldwide medical ultrasound (US) market continues to grow annually. Its non-ionizing nature, real-time capabilities and relatively low cost, compared to other imaging modalities, have led to significant applications in many different fields, including cardiology, angiology, obstetrics and emergency medicine. Facilitated by ongoing innovations, US continues to change rules and norms regarding patient screening, diagnosis and surgery. This huge and promising market is constantly driven by new imaging and processing techniques. From 3D images to sophisticated software, hardware and portability improvements, it is clear that the status of US as one of the leading medical imaging technologies is ensured for many years ahead. However, as imaging systems evolve, new engineering challenges emerge. Acquisition, transmission and processing of huge amounts of data are common for all ultrasound-based imaging modalities. Moreover, achieving higher resolution is constantly on demand, as improved diagnosis could be achieved by better visualization of organs and blood vessels deep within tissues.
In this article, our goal is to motivate further interest and research in emerging processing techniques, as well as their applications in medical ultrasound, enabled by recent advancements in signal processing algorithms and deep learning [1]. We address some of the primary challenges and potential remedies from a signal processing perspective, by exploiting the inherent structure of the received US signal [2], [3]. In particular, we discuss the mathematical models which underlie the received signals and show how prior information can be exploited to achieve unprecedented improvements in both data acquisition trade-offs and resolution. We begin by detailing the processing steps necessary to obtain an ultrasound image in the most commonly used schemes today. Once these are established, we turn to the challenges inherent to this process, focusing on data rates, access to the actual signal returns, and resolution.
The first challenge we address is that of the data-rate bottleneck in US imaging. The formation of a US image requires large amounts of data due mostly to two factors: High sampling rates, multiples or even magnitudes of order beyond the Nyquist rate to ensure high resolution beamforming, and a large number of data channels (and corresponding receivers), typically tens to hundreds, to enable high spatial resolution. These high data rates constrain advanced processing algorithms and portable hardware [3]. Typically, to enable real time imaging, a significant amount of the processing is done in hardware, with only reduced-dimensional digital images being exported from the device. This inability to access the original, raw signals ('channel data') limits the ability to take advantage of advanced, digital processing and machine learning methods [4]. Additionally, processing in hardware means that bulky and costly systems must be used, housing all required electronics, limiting portability and increasing cost. Currently, whenever reducing system costs becomes a priority, an unnecessary reduction in performance follows.
To address the problem of data-rate reduction, we review methods brought upon by the ideas of compressed sensing [2], [3], [5], where prior knowledge about the frequency content and structure of ultrasound signals may be used to represent it using a reduced set of samples. In particular, we exploit the fact that the recived signal is sparse, centered around the transmitted frequency, and that the beamformed signal used to form an image corresponds to time-shifted and scaled reflections of the transmitted pulse, thus obeying the Finite Rate of Innovation model [6]. The result is an ability to reconstruct an equivalent image using fewer samples, reducing required sampling rates by orders of magnitude, even below the Nyquist limit of the transmitted pulse [6].
To reduce the number of channels used in US probes, we turn to the field of sparse array design [7], [8]. With an appropriate choice of sparse array and post processing, the performance of a large physical array, as characterized by the beam pattern, can be achieved using fewer elements enabling a reduction in the number of elements without loss of performance. Algorithms such as Convolutional Beamforming [9], can be used to reconstruct an equivalent image with four to ten times fewer channels.
Addressing the data bottleneck leads to the availability of US channel data, enabling the application of emerging deep learning techniques directly on the raw channel returns. This approach can improve a seemingly limitless array of applications, including noise reduction, enhanced resolution, advanced diagnostic tools, addressing technician dependence, and rapid image formation, to name a few [4], [10]. We demonstrate two such deep learning algorithms: the LISTA algorithm [11], used to perform efficient recovery of an ultrasound signal from a small set of samples given its sparse representation [12], and the ABLE algorithm [13], used to form high resolution US images from their channel data. Moreover, reducing the data size and computational load needed in today's traditional US machines, paves the way for over-WIFI, efficient, wireless US imaging, bringing US machines closer to the patients, independent of their geographical location.
Finally, we turn to address the challenge of increasing US image resolution [14]. As a wave-based imaging modality, conventional spatial resolution for US signals is on the order of its wavelength, ranging from a few tenths of a millimeter to a few millimeters, depending on application [15]. This prohibits the use of traditional US imaging in applications where the features being imaged are of microscopic scale, such as the microscopic blood vessels developing around malignant cancerous tumors, or around inflamed regions such as the abdomen in Crohn's disease. Adressing this challenge may lead to unprecedented, non-invasive, non-radiative diagnostic capabilities for early detection and/or treatment of many newly detectable pathologies. This challenge is addressed by reviewing techniques in the field of Contrast-Enhanced UltraSound (CEUS) imaging, where the injection of designated agents into the blood stream allows for the exploitation of prior knowledge about their frequency domain behaviour, as well as their sparsity in several domains [16]. This results in an ability to map subwavelength features such as the microscopic vessels (microvasculature) created near malignant tumors. We review the use of sparsity and deep learning in this context, and demonstrate the improved resolution that can be obtained by algorithms utilizing them, such as the SUSHI method [17].
We conclude our review by discussing several outstanding challenges related to signal processing and image processing, and the exciting benefits their solution may enable. We hope this review will inspire further work by the signal processing community to bring ultrasound closer to the patient side and pave the way to high quality and portable ultrasound imaging at the home and primary care.

II. BASICS OF ULTRASOUND IMAGING
We begin by detailing the stages taken to perform ultrasound imaging. The most common mode of US imaging is referred to as Brightness mode, or B-mode. In B-mode, a gray-scale image of a desired plane (or volume in the case of 3D) is generated, with the brightness at each point representing its acoustic reflectivity. Scanning is performed using a designated system with a probe (transducer) moved around to determine the imaging region.
The typical stages of B-mode formation include: Transmission, where the imaging plane is insonified by ultrasonic pulses emitted by the probe; Reception, where their echoes are received in the transducer elements; Delay, where geometric calculations are used to match the timing of the recorded signals across the elements; and Sum, where the signals of all elements are summed (or averaged) to yield a brightness value at each point. The Delay and Sum stages are often lumped together in an algorithm termed Delay-and-Sum (DAS) Beamforming.
Once summed, the signal is referred to as RF data and undergoes further processing to enhance the resulting image, such as envelope detection, up-sampling, dynamic range compression, contrast adjustment and edge enhancement [18], before being stored as digital image data. In most commercially available US systems, the sheer size of data before summing usually means that only RF data, or even image data, is exported digitally, with prior stages performed in hardware.
In transmit stage, the imaging region is insonified with acoustic energy, most commonly by the formations of thin, focused, axial beams, transmitted one at a time and swept laterally across the imaging region. In receive stage, the signals acquired from each such beam are referred to as an image 'line', since they are each used to form the brightness along a line in the beam axis. To mitigate the effect of echoes from one beam to the reception of another, a beam is transmitted only after the previous line's echoes are acquired. Once a line is received, its data may be described as a set of signals ϕ m (t; θ), where m denotes the channel number, t denotes time and θ the lateral coordinate of the beam. In this form, the data is referred to as channel data.
In the delay stage, an appropriate time delay is applied to each channel so that the signals are aligned in time to yield the delayed signalφ whereφ m (t, θ) is the delayed signal for channel m at time t and τ m (t; θ) is the time delay for the same channel, given by Here c = 1540m/sec is the effective speed of sound in tissue and δ m is the distance between the beam origin and receiving element m. In the Sum stage, the signals are summed to yield the brightness-mode value at each point in the line: where Φ(t; θ) denotes the beamformed data. Repeating the process for each line yields a brightness value for each point in the imaging region. In practice, the delay is performed in the digital domain after sampling on a sufficiently dense grid, to enable implementation of (2) in digital. Many improvements upon this basic DAS scheme have been developed over the years, such as apodization and adaptive beamforming. See box "Advanced Beamforming" for more information. In practice, DAS remains the most commonly used beamforming technique, thanks to its simplicity and applicability for real-time imaging.
Although the method of focused-beam, line-based beamforming described above is the most common in US imaging, other techniques have been developed for various settings. Plane wave and diverging wave schemes insonify the entire scan region with a single pulse allowing for high frame rate [19], [20]. M-Mode imaging narrows the scan region to a single line, monitoring axial movements at a high refresh rate [21]. Finally, various versions of Doppler ultrasound use designated transmit schemes to improve spectral resolution and estimate flow velocities [22].
III. DATA-RATE REDUCTION As previously described, the standard approach to digital beamforming is through DAS, thanks to its relatively low computational load. The delays applied in DAS typically are on the order of nanoseconds, which results in high sampling rate requirements [23]. For practical reasons, US signals are sampled at lower rates, in the order of tens of megahertz,  with rounding or interpolation used to perform fine delaying. This reduced rate is still higher than the actual Nyquist rate of the signal, which is twice its bandwidth [3]. For fine delay resolution, a common rule of thumb is to sample at a rate 4-10 times higher than the central frequency of the transducer (which typically is 1-20 MHz, depending on application). Sampling at such high rates leads to huge amounts of data that need to be stored and processed.
Here we review techniques for reduction in sampling rate, based on the combination of compressed sensing (CS) [5], [6] and sub-Nyquist sampling [3] and by exploiting the finite rate of innovation (FRI) [24] structure of the beamformed signal. We also consider methods based on algorithm unfolding in the recovery stage [12], [25]. We begin by describing an approach of equivalently implementing DAS beamforming in the Fourier domain, referred to as frequency domain beamforming (FDBF) [6]. We then show how FDBF can pave the way to sub-Nyquist sampling, leading to rates that are much lower than those imposed by time-domain considerations. The suggested sub-Nyquist system is implemented based on the ideas of Xampling presented in [5] and further elaborated on in the box "Xampling for US imaging".

A. Frequency-Domain Beamforming
In [6] the authors suggested transforming the process of beamforming into the frequency domain in order to reduce the sampling rate, showing that beamforming can be performed in frequency, using a small set of Fourier coefficients. Let c[k] denote the kth Fourier series coefficient of the beamformed signal in (3) andĉ m [k] be the kth Fourier series coefficient of the delayed signal at channel m. Due to the linearity of the Fourier transform and following (3), where M is the number of channels andĉ m [k] are the delayed Fourier coefficients, given bŷ Here, I [a,b) is the indicator function, equal to 1 when a ≤ t < b and zero otherwise. The beam is supported on [0, T B (θ)), where T B (θ) = min m τ −1 m (t, θ) and T is the pulse penetration depth [6].
In [6], it is shown that the Fourier coefficients of the delayed signal can be expressed aŝ where c m [k] are the Fourier coefficients of the received signals in each channel before the delay is applied. The coefficients Q k,m;θ [n], illustrated in Fig. 3, are the Fourier coefficients of a distortion function, determined solely by the the imaging setup and geometry, essentially transforming the delay stage of beamforming from the temporal to the frequency domain. These coefficients may be computed offline and stored in memory for real time applications. Due to the decay properties of {Q k,m;θ [n]}, the terms in the last sum may be neglected outside a finite set around n = 0. The summation in (9) can thus be approximated by a relatively small finite sum

XAMPLING FOR US IMAGING
Xampling is a sampling method designed to sample signals below their Nyquist rate, in a way that enables perfect reconstruction [26]. In the context of FRI signals, which is the case in US imaging, Xampling can be used to obtain a sub-Nyquist representation by sampling only a small subset of frequency values. In US imaging, each signal has one main band of energy in the frequency domain, around the central frequency of the probe. Using the Xampling mechanism a subset of the main band of energy is obtained. Specifically, let κ be an arbitrary set comprised of K consecutive frequency components, and let ϕ m (t) be the received signal at channel m. The set of subsamples of the received signal is obtained using an appropriate sampling filter [5] where ϕ m (n s ) is the n s th sample of the analog US signal and where r is a constant determined by the support of the transmitted pulse, T is the duration of the received signals, and T sN is the sub-Nyquist sampling period. The filter f (t), is designed according to the required set κ, and effectively zeros out the frequency components that are not included in κ. Explicitly, the frequency response of f (t) satisfies Therefore, the number of samples is equal to the number of Fourier coefficients of interest. After filtering the signal, it is sampled at its effective Nyquist rate. Finally, applying a Fourier transform on the sampled signal results in the Fourier coefficients of the US signal, enabling the application of Fourier domain beamforming [6]. 2 illustrates Xampling for US imaging. Xampling framework for Fourier domain beamforming. In the first stage, signals received at all channels are filtered using an analog filter that obtains the desired set of frequency components κ, as described in (6). The signals are then sampled at their effective Nyquist rate, followed by an FFT and application of FDBF, resulting in low rate Fourier coefficients of the beamformed signal. where N 1 , N 2 are chosen empirically. The combination of (10) and (7) yields the FDBF formula: By applying an inverse Fourier transform on c[k], the beamformed signal is obtained.
The main advantage of this method is that it requires only a limited set of significant Fourier coefficients of the received signals. Each such coefficient is calculated using FFT, as a linear combination of samples of the signal. Thus, for fewer Fourier coefficients, less samples are required, and a lower sampling rate may be used.

B. CS Based Beam Reconstruction
A further reduction in sampling rate may be obtained by acquiring only a subset of the non-negligible Fourier coefficients. However, since in this case the signal is sampled below its effective Nyquist rate, applying an inverse FFT results in aliasing in the time domain. To overcome aliasing effects, we exploit the structure of the beamformed signal, relying on techniques from compressed sensing (CS) and finite rate of innovation (FRI) sampling [3], [5], [6]. More specifically, due to the way a beamformed signal is formed, it can be expressed as a shifted and scaled sum of replicas of the (known) transmitted pulse where h(t) is the transmitted pulse, L is the number of scattering elements in direction θ, {b l } L l=1 are the unknown amplitudes of the reflections and {t l } L l=1 denote the time of arrival of the reflection from the lth element. Such a signal is completely defined by the 2L unknown amplitudes and delays.
By quantizing the delays {t l } L l=1 with quantization step where T s is the sampling period and f s is the sampling frequency, we may write the Fourier coefficients of the beamformed signal as: where N = T /T s , h[k] are the Fourier coefficients of the transmitted pulse, and Defining an M -length measurement vector c with kth entry c[k], (13) may be rewritten as where H is an M ×M diagonal matrix with h[k] as its entries, D is an M × N matrix formed by taking a set of rows from an N ×N Fourier matrix, and b is a vector of length N st with lth entry b l , where N st is the number of samples required for standard DAS beamforming. Once formulated this way, the problem is that of recovering a sparse vector b, given measurements c. A typical beamformed ultrasound signal is comprised of a relatively small number of strong reflections and many scattered echoes. Thus, the vector b defined in (15), is typically not strictly sparse, but rather compressible. This property can be captured by using the l 1 norm as an objective function in an optimization problem: with a parameter representing a noise level. Problem (16) can be solved using second-order methods such as interior point methods [27], or first-order methods, based on iterative shrinkage ideas such as the NESTA algorithm [6], which was shown to be highly suitable for US signals. The techniques discussed above, are evaluated in Fig. 7 using in-vivo liver data of a healthy volunteer. Acquisition was performed using the Verasonics Vantage 256 System, using the 64-elements phased array transducer P4-2v. The frequency response of this probe is centered at 2.72 MHz, and a sampling rate of 10.8 MHz was used, leading to 1920 samples per image line. In (a), DAS beamformed US image is presented, using 1920 samples. Sub Figs. (b) and (d) show the potential usage of FDBF and CS reconstruction using 230 and 130 samples, leading to 8, and 15 fold reduction in data size and sampling rate, respectively. A description of a specific imaging setup that was used to implement the suggested methods above is given in the "Experimental BF demonstration" box. The described methods can also be extended to 3D US imaging [28]; an example of implementing these methods for 3D imaging is given in Fig. 8, presenting a 12-fold rate reduction.
Another powerful means of solving the recovery problem in real time is through the use of deep learning networks, specifically using algorithm unrolling [11], [25], [29]. In [12], it was suggested to use the unrolled LISTA architecture, which is a learned version of the well known ISTA method for sparse recovery. For further details, see "Sparse Recovery Using LISTA" box. As opposed to complicated deep learning approaches for beamforming, LISTA consists of fast and efficient training and testing steps since only a small fraction of the samples of the beamformed signal needs to be stored for training the network and for full image recovery. In addition, the computational complexity of LISTA, reduces to only O(n), when compared to standard iterative approaches like NESTA, which suffers from a O(n 2 ) computational load. In Fig. 7, the usage of LISTA for recovery of subsampled data is depicted in subfigures (c), and (e). Using LISTA, only 230, 130 samples are needed and the resulting images yield better visual and quantitative results compared to NESTA [12].
The results, shown in Fig. 7, demonstrate that a significant reduction in sampling rate, while fully recovering the US image is possible. Moreover, the usage of deep unfolding methods, such as LISTA, provides an efficient recovery method for subsampled US signals. The following section, suggests another approach for reduction in data rates by reducing the number of receiving channels.

IV. SPARSE ARRAY DESIGN
We now turn to designing arrays which can enable a reduction in the number of required elements, without compromising on image quality. Image contrast and resolution are limited by the number of elements in the probe, leading to designs of hundreds of elements in typical transducers. Thus, techniques that reduce the number of receive channels are of great importance in paving the way to portable low cost devices.

A. Convolutional Beamforming Algorithm (COBA)
Recall that (3) with ∠ϕ m (t) and |ϕ m (t)| denoting the phase and magnitude of ϕ m (t), respectively. Defining s(t) and u(t) as the vectors of size 2N − 1 whose entries are s n (t) and u n (t), respectively, s(t) can be written as where * s denotes a linear convolution in the lateral direction, which can be implemented efficiently using the convolution theorem and FFT. This beamforming procedure is referred to as convolutional beamforming algorithm (COBA).

EXPERIMENTAL BF DEMONSTRATION
A low-rate frequency domain beamforming mechanism is implemented on an ultrasound imaging system. The setup is shown in Fig. 4. A state of the art GE ultrasound machine, phantom and ultrasound probe with 64 acquisition channels were used for the scan. The radiated depth r = 15.7 cm and speed of sound c = 1540 m/sec yield a signal of duration T = 2r/c 204 µsec. The acquired signal is characterized by a narrow band-pass bandwidth of 1.77 MHz, centered at a carrier frequency f 0 ≈ 3.4 MHz. The signals are sampled at the rate of 50 MHz and are then digitally demodulated and down-sampled to the demodulated processing rate of f p ≈ 2.94 MHz, resulting in 1224 samples per transducer element. Linear interpolation is applied in order to improve beamforming resolution, leading to 2448 samples used to perform beamforming in time. Fig. 6(Left) presents a schematic block diagram of the system's front end. As illustrated in Fig. 6(Right), in-phase and quadrature components of the received signals were used to obtain the desired set of their Fourier coefficients. Using this set, beamforming in frequency was performed according to (11), yielding the Fourier coefficients of the beamformed signal. Frequency domain beamforming was performed at a low rate: M = 100 Fourier coefficients of the beamformed signal were calculated, using K = 120 Fourier coefficients of each one of the received signals. This corresponds to 240 real-valued samples used for beamforming in frequency. Hence, beamforming in frequency is performed at a rate corresponding to 240/2448 ≈ 1/10 of the demodulated processing rate. This implementation was done on a commercial system which samples each channel at a high rate. Data and processing rate reduction took place following the Fourier transform, in the frequency domain. The set of 120 Fourier coefficients of the received signals, required for frequency domain beamforming, are obtained from only 120 low-rate samples. The performance of the proposed method was verified by scanning a healthy consenting volunteer. Fig. 5 compares low-rate beamforming in frequency and standard time-domain beamforming of a cardiac scan of the volunteer. As can be seen, sufficient image quality is retained, despite the significant reduction in processing rate.
where ω 0 is the central frequency of the transducer and δ is the distance between two adjacent array elements. Based on (31), the beam pattern of COBA can be written as the product of two DAS beam patterns, leading to where a n is the intrinsic apodization, calculated by a = I M * s I M , with I M denoting a binary vector whose mth entry is 1 if m ∈ M . The significance of the last expression is in demonstrating that COBA yields a beam pattern equivalent to a DAS beamformer using a ULA twice as large as the original array used for imaging. As (35) implies, effectively the array used for imaging is larger, and corresponds to the convolution of the original array with itself. Therefore, [9] suggested using a sparse array, which is formed by removing some of the elements from the original array such that a desired beam pattern is formed, in a method called sparse COBA (SCOBA). The essential idea is to remove elements such that after convolution, the original beampattern is preserved.
More specifically, the behaviour of the sparse array is determined by the sum co-array [36]. Given an arbitrary array U , the sum co-array S U is the array that includes all distinct elements of the form n+m where n, m ∈ U . The beam pattern of the sparse convolutionally beamformed signal is where θ is the angular direction, c is the speed of sound and ω 0 is the transducer central frequency. By appropriately choosing the sparse array U ⊆ M in a way that its sum co-array yields a full ULA, the reconstructed US image after COBA will have the same image quality as DAS beamforming. One can also use sparse COBA for super resolution (SCOBAR), when choosing a sparse array such that its sum co-array is larger than the original ULA, which leads to enhanced image contrast and resolution. Further reduction in the number of elements was considered by using a fractal array geometry [7], [8], [37]. The box, "The Sum Co-Array" describes several of these array geometries. 10 compares in-vivo cardiac image reconstructions of a standard DAS beamformer and the convolutional beamformers presented here, using different array geometries for SCOBA, and SCOBAR. The data was captured using a GE breadboard, phased array probe with 64 channels. The radiated depth was 16 cm, probe carrier frequency 3.4 MHz and sampling frequency 16 MHz. The number of elements used was 21 for SCOBA based on the geometry given by (40), 27 elements for SCOBAR using (43), and 15 elements for SCOBA with fractal array geometry of order 3 and generator array G = {0, 1}. The DAS beamformed image is shown in 10 (a), while (b) demonstrates the improvement of using COBA with the full ULA. SCOBA, and fractal geometry based SCOBA are presented in (c), (e) and achieve similar image quality to DAS with a fourth of the number of elements, while SCOBAR, presented in (d), yields even higher resolution and improved image quality.

B. Combining Low-Rate Sampling and Sparse Arrays
The ideas of temporal (sampling) and spatial (array design) dilution, discussed in the last two sections, may be combined, achieving further reduction in data size. This technique is referred to as the Compressed Frequency domain Convolutional Beamforming Algorithm (CFCOBA) [38].The idea in CFCOBA, is to calculate the sub-sampled Fourier coefficients of the convolutionally beamformed signal, and to reconstruct it using NESTA or LISTA, based on the FRI and CS frameworks. In [38], it is shown that the convolutionally beamformed signal obeys an FRI model, using the square of the known transmitted SPARSE RECOVERY USING LISTA Learned ISTA (LISTA) [11], [25] is a data-driven, learned counterpart of the widely used iterative shrinkage/thresholding algorithm (ISTA) [30] for sparse recovery. Consider a sparse recovery problem where y is the measured signal, A is the measurement matrix, both known, λ a regularization parameter and x is the underlying sparse vector to be recovered. ISTA solves (17) iteratively through the iterations where x k+1 , x k are subsequent approximations of the solution, W e = µA T , W t = I − µA T A for some chosen step size parameter µ, and S λ is a soft thresholding operator with threshold λ.
Adapting ISTA to a learning-based framework is performed by unrolling: cascading layers of a deep network, each mimicking an iteration in the iterated algorithm [25]. In the unfolded version of LISTA, the layers consist of trainable convolutions W k e , W k t , and a trainable shrinkage parameter λ k , corresponding to the threshold in ISTA, for 1 ≤ k ≤ K, where K is the total number of layers. The thresholding operation in LISTA is replaced by a smooth approximation (though other smooth approximations are possible) with the operations performed element-wise. Training LISTA is done in a supervised manner, by first recovering a set of sparse codes x i , 1 ≤ i ≤ N , given a set of measurements y i and the known measurement matrix A. This recovery is typically performed using an iterative sparse solver, e.g. ISTA, FISTA, and NESTA [30], [31]. Training is typically performed by minimizing the meansquared-error (MSE) between the unrolled network's output, given the measurements y i and their respective latent codes using standard loss minimization methods, such as stochastic gradient descent. The basic architecture of LISTA is given in Fig. 9 within the yellow dashed line. Among the advantages of unrolled networks over their iterative counterparts are fixed computational complexity (a feedforward network requires only a single pass), which leads to typically much faster execution times, improved performance, and no need to know the matrix A. Fig. 9 describes the architecture suggested in [12], where the last layer G, transfers the sparse representation of the output, into a beamformed signal. Unrolled networks may also be used to recover the positions of strong reflectors, such as contrast agents in an US frame, given that their positions in space and time are sparse [4], [32], [33]. Fig. 9. LISTA block diagram, with each layer applying an iteration step (rectangles represent trainable convolutional layers and S λ i represents soft thresholding). Φ stands for the beamformed signal. Image adapted from [12].
pulse, hence it can be recovered based on the same approach discussed in Section III-B. In Fig. 7 (f) the combination of FDBF and COBA is demonstrated, using 230 samples out of 1920 used for DAS and only 15 receiving channels, based on the fractal array geometry, leading to an overall 36-fold reduction in data size for a liver image. Another example of the huge reduction is depicted in Fig. 10 (f), where invivo cardiac images were acquired with the setup described in Section IV-A, using a GE breadboard US machine. It can be seen, that using CFCOBA, only 100 samples and 15 elements out of 3328 samples and 64 elements, respectively, are needed for high quality reconstruction, yielding a rate reduction by a factor of 142. Reconstruction of combined time and space sub sampled data using the LISTA approach [12], [25], as described in Section III-B and further elaborated on in the "Sparse Recovery Using LISTA" box, is shown in Fig. 7 (g), which depicts 36-fold reduction in data rates for reconstruction of the beamformed signal using the efficient LISTA algorithm. This huge reduction in data size, makes wireless US systems possible, as illustrated in Fig. 12.

V. SUPER-RESOLUTION ULTRASOUND IMAGING
We now review US imaging performed at a resolution finer than its wavelength, and discuss the methods to obtain it.
As a wave based imaging modality, the resolution of ultrasound imaging is bound to the diffraction limit, meaning ADVANCED BEAMFORMING In US beamforming, channel signals are delayed and summed to yield brightness values in an imaging region, with the objective of spatially focusing and maximizing the energy in the direction of interest. Apodization is the process of assigning weights to each channel before summing, to reduce the contribution of certain channels, for example due to low SNR. Classically, the weights are represented using a window function (e.g. Hamming or Tukey), centered around the element at the center of the beam. To formulate apodization mathematically, we consider a ULA of M transducer elements. Assume that all recieved signals have been delayed correctly for a given point in the imaging region. Thus, the delayed signal at channel m can be written as y m = x + n m , where x is a signal of interest that represents the reflection of the imaged region and n m ∼ N (0, σ 2 n ) is additive interference, including off target reflections, clutter and system noise [34]. In vector form where y is the delayed data, 1 is a vector of ones and n is the vector of interference. The classic DAS beamformer, averaging the delayed data over all elements, can be written as With a windowing function w m , we obtain Adaptive beamformers assign weights using the data itself. One such algorithm is the minimum variance (MV) beamformer [35], which maintains a unity gain in a chosen focusing direction (the foresight), while minimizing the energy received from other directions. It can be written as the solution to the following minimization problem: The solution is where R y = E[yy H ] is the received correlation matrix. This requires the inversion of R y , creating computational issues and prohibiting it from being used in practice in US applications. Approaches such as Eigenvalue-Based Minimum Variance, spatial averaging and diagonal loading [13] have been developed to remedy this, but they still suffer from a high computational load.
To overcome these issues, Iterative Maximum a-Posteriori (iMAP) beamforming, was suggested in [34]. Assuming the beamformed signal is a gaussian random variable x ∼ N (0, σ 2 x ), and that the noise at each element is an uncorrelated Gaussian variable, the iMAP estimator of x is given by The solution is where the weights are based on parameters of the prior distribution of the signal of interest and the interference. In most realistic scenarios, the parameters are unknown and need to be estimated from the data, hence the use for iMAP. The statistics are calculated based on the previous estimate in an iterative fashion, with x DAS used as the initial estimate: Iterating between (27) and (29) yields increasingly better estimates of the distribution of the parameters. Deep neural networks may also be used to perform adaptive beamforming. One such network, ABLE [13], uses MV beamformed images, computed offline, as targets, to train adaptive weights that minimize the loss function: where P ABLE and P M V are the pixel values of the network output and target, respectively. A comparison between different adaptive beamforming methods is given in Fig. 1, where iMAP2 indicates two iterations of (27) and (29).
that it cannot accurately image objects on a smaller scale than its wavelength. The wavelength used during scanning is often determined by considerations of imaging depth and transducer constraints, and thus cannot be chosen to be arbitrarily small. This practically limits the capabilities of US imaging at roughly a tenth of a millimetre.
The vast majority of modern super-resolution ultrasound methods rely on Contrast Enhanced UltraSound (CEUS) [39], [40], where an ultrasound contrast agent (UCA), is injected into the bloodstream and then used for localization. The UCA is composed of microscopic bubbles (Microbubbles or MBs). Prior knowledge about the differing behaviour of the MBs and tissue in the scanning region is then exploited in order to separate their respective signals. Processing of the MB signal is performed by relying on further prior knowledge such as their sparsity or shape (as observed in the US image) used to pinpoint locations of MBs over multiple frames. Stacking the frames yields a vascularity map of a finer resolution than achieved in B-Mode imaging.
We first discuss methods relying on localization -the tracking of, and subsequent replacement of MBs with point objects. Then, we show how sparsity in different domains can be exploited to yield non-localization based techniques.

A. Ultrasound Localization Microscopy
Ultrasound Localization Microscopy (ULM) is a family of algorithms which achieve super resolution by a simple principle: Objects of sub-wavelength size undergo diffraction, yielding a blurred shape in the US image. Thus, given an observed blurred shape at a certain point in space, it may be replaced by a point source representing a MB. The set of all points where an MB is detected is then used to construct a super resolved vascularity map [16]. The stages of ULM consist of: Acquisition, where multiple frames are acquired; Detection, where MB signal is separated from tissue; Isolation, where nearby or overlapping MB signals are discarded; Localization, where the blurred shape created due to the diffraction of each MB is replaced by a point reflector; Tracking, where estimation of flow is done based on MB positional differences between consecutive frames; and Mapping, where the data is accumulated and post-processed for a visualization of the arterial map. An illustration of these stages is given in 13.
During acquisition, a large number of frames is acquired. Hundreds, or even thousands of frames may be needed in order to produce a satisfactory vascular map. To this end, ultrafast imaging methods have been developed, often trading off nominal image resolution in favor of frame rate [33]. Using frame rates in excess of 500 frames per second, an analogue to optical localization microscopy is obtained by capturing the transient signal decorrelation of MBs using plane wave imaging [20]. Furthermore, a designated transmit-receive scheme may be used in order to enable improved separation of MB from tissue, such as harmonic imaging techniques exploiting their pulsating behaviour; for further details, see "Harmonic Imaging in Contrast Enhanced Ultrasound" box.
In the detection stage, the positions of microbubbles within each frame are determined, and they are separated from the tissue. Classically, a centroid localization approach was used, typically by fitting a two-dimensional gaussian brightness curve per MB. More recently, assuming the MBs are sparse i.e. only a few are present per frame, their detection has been treated as a sparse recovery problem [33]. The formulation and example methods for sparse recovery are given in the "Sparse Recovery Using LISTA" box. For example, one can define x as the (vectorized) map of MB locations and seek a super resolved image by solving for a sparse x using Here y is the vectorized image frame, and A is the measurement matrix, with each column being a shifted version of the point spread function (PSF), an estimated blurring filter representing the system's response to a point object. Using LISTA to solve (46) enables super resolution without requiring exact knowledge of the system PSF and yields high resolution recovery techniques [4]. Despite its high spatial resolution, ULM's reliance on MB sparsity for the sake of detection and isolation implies that MB concentration must be limited during scanning. This, in turn, means that only a small number of MBs may be detected per frame, leading to a requirement for a large number of frames, with long associated acquisition periods, dictating stringent restrictions on patients and forcing frames with significant overlap to be discarded, lowering efficacy.
The advent of deep learning in recent years has led to the development of methods that can circumvent the requirement for sparsity by analyzing interference patterns between adjacent MBs. One such example is the frame-based Deep ULM method presented in [33]. Deep ULM generates simulated ground-truth data representing microbubbles at known locations on a grid denser by orders of magnitude than those of the ultrasound data. The MB positions are convolved with point spread functions estimated from empirical measurements, representing the blurred shape of the MBs due to diffraction. The network is then trained to minimize the difference between THE SUM CO-ARRAY Assume a uniform linear array (ULA) with 2N − 1 elements, at positions I = {− (N − 1), . . . , (N − 1)}, on a grid with spacing δ. This array can be thinned while exhibiting a fully sampled equivalent array after COBA. This can be utilized for sparse beamforming, without compromising image quality. The sum co-array for a ULA with position set I is defined as the setS This includes the positions of all array elements, as well positions obtained by their summation/difference. A sparse array with position set J can be obtained from the full ULA with position set I by removing some of its elements. Ideal removal of the elements would be such that they lie on the sum co-array of J: In (38), I, J and S j are integer sets such that I represents the element locations of the fully sampled ULA, while J and S j represent the element locations for the thinned array and its resulting sum co-array. Now, assume that N is not prime and can be factored as N = AB, where A, B = N + (natural numbers including zero). Define Let U A + U B = {n + m : n ∈ U A , m ∈ U B }. Then Denoting U ⊂ I the array geometry defined as it holds that I ⊂ S U , where S U is the sumset of U . Therefore, the family of sets (41) satisfies (38), and the number of elements in each set is 2A + 2B − 3. This configuration does not achieve the set S I which corresponds to the sum co-array of the full ULA I. To achieve this, another set is introduced and now the sparse array geometry is defined as The array V is in fact obtained by adding two small ULAs to U , each of size A − 1 at its edges. It can be verified that This implies that elements positioned on the sum co-array of V can produce the same images as those with the sum co-array of the full ULA I, despite having fewer elements. Further reduction in the size of the array is obtained using fractal arrays [8], which are defined recursively by where r is the array order, G is the generator array in fractal terminology, and min(G) = 0. The translation factor L is given by L = 2 max(G) + 1. By choosing the generator array and the fractal order appropriately, desired beampatterns can be obtained. Other approaches for separating MBs from the tissue include spatio-temporal filtering based on the singular value decomposition (SVD) [4]. SVD filtering includes collecting a series of consecutive frames, stacking them as vectors in a matrix, performing SVD of the matrix and removing the largest singular values, assumed to be related to the tissue.
To overcome difficulties of SVD filtering, the task of clutter removal was formulated as a convex optimization problem by leveraging a low-rank-and-sparse decomposition in [32]. Then, the authors proposed an efficient deep learning solution to this optimization problem through unfolding a robust PCA algorithm, referred to as Convolutional rObust pRincipal cOmpoNent Analysis (CORONA). This approach harnesses the power of both deep learning and model-based frameworks, and leads to improved separation performance. Fig. 12. Cloud-based US processing (illustration). Probe is connected to the cloud for sampling, on cloud computation, and display tablet. Detaching the processing task from data acquisition will reduce machine size and increase user flexibility, as well as encourage new algorithm development. This also makes data available outside the confines of the specific machine it was acquired with, regardless of its vendor.  .

B. Non-localization methods for super-resolution
Several approaches have been developed to mitigate the reliance of super-resolution imaging on sparsity alone. Here we review two of them: SUSHI, based on exploitation of 2ndorder statistics of a pixel between frames [17], and Triple-SAT, which exploits flow estimation [40].
Sparsity-based super-resolution ultrasound hemodynamic imaging (SUSHI) makes use of the statistical independence between signal fluctuations, over time, of different vessels [17]. This is done through estimation of the per-pixel auto-correlation, which can be expressed as where h is the system PSF, ∆ xL is the pixel size, g n [τ ] is the autocorrelation function of the temporal fluctuations of a pixel n and g il [τ ] is the cross correlation of pixels i, l. Here, τ stands for discrete pre-determined delay of the autocorrelation function, and i, l are indices of dependent volume cells located in the same streamline. The significance of this expression is that MBs flowing independently in different vessels are uncorrelated, and therefore yield only expressions of the form of the first term, where the squaring of the signal creates a squaring of their PSF, making it narrower and thus improving resolution.
A dense imaging grid with microscopic resolution is then created, where the vascular map can be assumed spatially sparse. SUSHI exploits this to perform sparse recovery in the spatial domain (e.g. using LISTA). Note that the sparsity assumption holds, even if the microbubbles densely populate a region of the frame, making use of MBs even if they themselves are not sparsely distributed locally. 16 demonstrates SUSHI on a Rabbit Kidney, illustrating the blurring caused by diffraction and the resulting, super-resolved arterial map. Due to it's reliance on temporal statistics of the signal, SUSHI's performance depends on the achievable frame rate, and becomes more effective with the advent of ultra-fast imaging modes [17].
When frame rates are low compared with vascular flow velocities the simultaneous, sparsity-based super-resolution and tracking (triple-SAT) may be used. Triple-SAT uses optical flow estimation to first estimate the velocities of MBs detected in a frame. The calculated velocities of subsequent frames are then used to estimate 'tracks' within the frame; an example of the method is depicted in Fig. 17. As MBs flow inside blood vessels, they are more likely to be found in certain positions in the imaging plane. These tracks are then weighted in, per frame, and increase the likelihood for signal in their region to be detected as MBs.

C. Current limitations of super-resolution CEUS
We now discuss some challenges of super-resolution US. Complete imaging of vascular maps, and capillaries more so, requires that they be populated in some frames during scanning, yielding trade-offs between several desired properties: For high probability of population, it is beneficial to take a large number of frames. However, for conventional imaging schemes, a number of frames in the order of thousands, which is typically needed to populate organs such as the liver or the heart, requires several minutes of acquisition, creating discomfort, logistic complications and huge amounts of data, as well as increasing the probability of patient and organ motion mid-scan. Thus, a higher frame-rate is desired. This may be achieved by using plane-wave or diverging-wave

HARMONIC IMAGING IN CONTRAST ENHANCED ULTRASOUND
Contrast-enhanced ultrasound imaging is enabled by injection of a contrast agent, consisting of microscopic, gas-filled bubbles, encapsulated in a phospholipid or protein shell [41], into the blood stream. These MBs have a typical mean diameter between 1 − 3µm, and circulate in the patient's bloodstream for a few minutes before dissolving. MBs exhibit nonlinear vibration when insonified in their resonance frequency, even in low pressure, meaning that their spectrum exhibits higher order harmonics, where linear tissue typically will only reflect the same frequency transmitted by the probe. Tissue signal interferes with the signal received from blood vessels, reducing image contrast. Several schemes for tissue suppression are currently in practice, and depend on two key conditions: Transmitted energy in the MB resonance frequency, and MB motion relative to tissue. Where a transmit frequency close to that of the MB resonance frequency is used, significant sub-harmonics and higher order harmonics are present in the MB signal, and harmonic imaging methods may be used. Harmonic imaging involves two key methods: Pulse inversion (illustrated in 15), where consecutive pulses have inverse polarity, and amplitude modulation, where the pulses have differing amplitudes. In pulse inversion, the sum of received signals cancels the linear (fundamental frequency) part of the signal, whereas in amplitude modulation, the difference between received signals, with scaling to equalize their amplitudes, will do the same. High-pass filtering may also be applied to further reduce residual tissue signal. After processing, the data contains mostly harmonic signals and can be processed in a normal Bmode beamforming (where imaging quality might be improved depending on application) or in a CEUS imaging scheme.   imaging, which are faster and allow insonification of the entire scan region using a single pulse, allowing for ultrafast imaging with frame rates of up to hundreds of frames per second. However, these schemes suffer from poor resolution, causing an increased blurring of MBs and difficulty in localization. The inherent 3D motion of organs and 3D structure of blood vessels also places a challenge and opportunity for 3D superresolution imaging schemes to be developed. Currently, 3D imaging modes typically rely on sequential acquisition of 2D slices of a given volume, effectively forcing a projection of the 3D image on the imaging plane, one at a time, as is shown in 8. Methods relying on a complete data-set acquired using designated, 2D or 3D arrays, may be able to exploit the complex structure in the data in similar manners to those proposed in this section, while taking into account further complexities introduced by the higher dimensionality. One such complexity is the increase in data sizes, which will challenge the device's ability to acquire at reasonable framerate.

VI. CONCLUSIONS AND OUTSTANDING CHALLENGES
In this review, we detailed common approaches for ultrasound imaging, and how their limitations can be overcome by modern signal processing and learning methods. In particular, we focused on super-resolution, reduction in channel-data rates and sparse-array design. For each, a mathematical model was described, as well as how it may exploit prior knowledge of the data structure to improve results. For example, using Fourierdomain knowledge of the signal and the pulse strucutre yields a reduction in sampling rate and thus data size. Modelling contrast agent sparsity in several domains, leads to superresolved images in conditions previously not enabled such as high concentrations and low frame rates. Further advances in signal processing techniques, as well as the increasing popularity of unfolding algorithms using deep neural networks, may lead to leaps in performance in upcoming years, and also help address issues such as technician dependence and improved diagnostics.
There are, however, still several major challenges and opportunities, all of which can benefit from a signal processing point of view.
Super-resolution imaging faces challenges in the fields of tissue suppression, organ motion detection and microbubble localization, with image and RF signal processing techniques constantly being developed and tested for different organs and scenarios. A development towards the direction of 3D superresolution imaging may be necessary in order to overcome several limitations, such as the 3D nature of microvascular flow and organ motion. Moreover, complex imaging schemes designed to increase frame rate e.g. by compounding several, rapidly transmitted plane waves, are continuously developed and tested, paving the way to early, noninvasive diagnosis of inflammatory and cancerous pathologies.
3D scanning capability is a key topic of research in the US community even regardless of super-resolution. True, realtime, 3D imaging is considered one of the ultimate goals in ultrasound imaging and is challenging due to the shear amount of data needed to image typical volumes of interest such as a heart or bladder (see Fig. 8), and the large number of RF chains needed to convert the received pulses to digital. Combining sub-Nyquist techniques and 2D or 3D sparse arrays, may reduce the amount of data to feasible dimensions, enabling sufficiently high frame rates at sufficient imaging quality. This would in turn reduce the dependence on the operator, as the current, limited performance of 3D scanning requires great expertise in order to yield valuable images. Data size reduction will pave the way towards wireless, cloud-based medicine, where experts from around the world may provide real-time analysis to data gathered remotely by technicians, in real time, using low-cost, portable systems, as shown in 12. Deep learning methods, developed directly on channel data, to target specific clinical applications and extraction of further clinical insights from channel data directly will allow unprecedented inference, regardless of data visualization constraints.
With advances in hardware, structure-based processing techniques, and the advent of deep learning methods, ultrasound has the potential to offer imaging on the edge to patients worldwide at an affordable price.