Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

So you think you can DAS? A viewpoint on delay-and-sum beamforming

So you think you can DAS? A viewpoint on delay-and-sum beamforming Perrot et al. So you think you can DAS? A viewpoint on delay-and-sum beamforming Vincent Perrot, Maxime Polichetti, François Varray, Damien Garcia compared to the size of the receiver. A telecommunication an- Abstract – Delay-and-sum (DAS) is the most widespread digital tenna then receives quasi-planar waves only. An array of beamformer in high-frame-rate ultrasound imaging. Its implementa- equally spaced unit antennas can be steered by phase-shifting tion is simple and compatible with real-time applications. In this view- to collect a radio wavefront from one direction. The radio waves point article, we describe the fundamentals of DAS beamforming. The of the different unit antennas are time-shifted with predefined underlying theory and numerical approach are detailed so that users delays (the greater the delays, the greater the receive angle) so can be aware of its functioning and limitations. In particular, we dis- that they can be summed to increase the signals in the desired cuss the importance of the 𝒇𝒇 -number and speed of sound on image quality, and propose one solution to set their values from a physical receive direction while attenuating those from undesired direc- viewpoint. We suggest determining the 𝒇𝒇 -number from the directivity tions [3]. This describes the basic concept of delay-and-sum. of the transducer elements and the speed of sound from the phase dis- Computer-controlled antenna technology based on DAS has persion of the delayed signals. Simplified Matlab codes are provided been the subject of intensive research for military radars and for the sake of clarity and openness. The effect of the 𝒇𝒇 -number and sonars [4], as it allowed antennas to be quickly steered to detect speed of sound on the lateral resolution and contrast-to-noise ratio was airplanes, submarines, and missiles. The DAS principle had investigated in vitro and in vivo. If not properly preset, these two fac- tors had a substantial negative impact on standard metrics of image also been one of the most useful techniques in seismological quality (namely and 𝐅𝐅𝐅𝐅 ). When beamforming with DAS in research. It allowed observation of small seismic events by con- vitro or in vivo, it is recommended to optimize these parameters in or- necting a large number of horizontally distributed seismic sen- der to use it wisely and prevent image degradation. sors [5]. In geophysics, DAS is often referred to as “diffraction summation”. It was the first computer implementation of seis- Keywords – beamforming, delay-and-sum, 𝒇𝒇 -number, speed of mic data analysis [6]. A similar technique that incorporates am- sound plitude correction and frequency-dependent phase correction, called “Kirchhoff migration”, was then proposed and is still in common use [7]. Amplitude correction due to wave spreading 1. Introduction is not a major concern in ultrasound imaging since the recorded ELAY-AND-SUM is the most basic digital beamformer signals can be time-gain compensated [8]. In addition, unlike D for medical ultrasound imaging. Because of its simplicity geophysical signals, medical ultrasound signals are narrow- and efficiency, it is very likely the most widely used in high- band, so there is little need for frequency-dependent phase cor- frame-rate (ultrafast) ultrasound. Before being used for ultra- rection. An important aspect of geophysics is the wide range of sound imaging [1], the technique of delay-and-sum (DAS) has propagation speeds, which has promoted the development of historically been linked to ground-based and airborne radar as advanced reconstruction techniques [6]. Although the study of well as telecommunication [2]. It originates from steerable ar- speed heterogeneity may be valuable in some ultrasound appli- ray antennas in shortwave communication [3]. Similar to a cations [9], [10], the speed of sound is in most situations con- medical ultrasound transducer, an antenna array is defined as a sidered homogeneous in soft tissue and assumed to be equal to group of connected individual antennas (or elements) that oper- 1540 m/s. As we will see, this fixed value can occasionally be ate together. In radar and ultrasonic imaging, a phased array suboptimal. refers to a computer-controlled array that is electronically steered to transmit or receive in different directions without DAS position in ultrafast ultrasound imaging – With the ac- moving the unit elements. Although the term is generally re- cess to open-architecture ultrasound scanners and their raw served for cardiac application, most transducers used in ultra- data [11], ultrafast ultrasound imaging has quickly gained pop- sound medical imaging are phased arrays since the delays can ularity in our community. Although the content of the manu- be modified electronically in transmission or reception. script could also be derived for focused beams (see Appendix 6.A), we opted for diverging [12] and plane waves to be in line Multidisciplinary DAS – In telecommunication, in contrast to with the most recent literature. “Ultrafast imaging” might be a ultrasound imaging, the involved travel distances are very large misnomer; however, we keep this abusive terminology because This work was supported in part by the LABEX CELYA (ANR-10-LABX- All the authors are with CREATIS, CNRS UMR 5220, INSERM U1206, 0060) and LABEX PRIMES (ANR-10-LABX-0063) of Université de Lyon, Université Lyon 1, INSA Lyon, Université Jean Monnet Saint-Etienne (e-mail: program “Investissements d'Avenir” (ANR-16-IDEX-0005). damien.garcia@inserm.fr; garcia.damien@gmail.com). 𝐅𝐅𝐅𝐅 𝐂𝐂𝐂𝐂𝐂𝐂 2 Perrot et al. unfocused beams are generally used for high-frame-rate pur- waves). Forming an image from sound requires estimating the poses. Even if DAS remains a widespread beamformer in high- round-trip traveltimes of the wavefronts towards and from all frame-rate ultrasound imaging, there is a continuing interest in scatterers. They can be explicitly expressed when transmitting the development of alternative beamforming methods. In par- circular and plane waves, as now explained. ticular, a number of recent studies have focused on adaptive beamformers that aim to improve contrast and spatial resolution provided by DAS. A series of popular adaptive beamforming techniques are discussed in [13]. In such investigations, the DAS is most often chosen as the substandard reference method. However, as it will be explained and illustrated, care must be taken to implement it correctly. In particular, special attention must be paid to the speed of sound and the receive aperture used during beamforming. A well-implemented DAS can have a sig- nificant impact on image quality. It is highly recommended to use it properly when comparing it with other beamformers. Our objective is to dissect the DAS by providing a detailed theoret- ical and pedagogical overview and suggesting some tricks for appropriate and efficient utilization. The advantages of the DAS are many: 1) DAS is based on Fig. 1. Parameters that define a circular-wave transmit. It becomes a 𝜃𝜃 - basic concepts of wave propagation (linearity, straight-ray tilted plane-wave transmit when 𝛽𝛽 → 0 . Figure adapted from [20]. propagation, weak backscattering); 2) its implementation is simple and can be parallelized; 3) it is numerically robust, fast To make things simpler, we focus on rectilinear arrays (for and compatible with real-time applications; and 4) because it is which the positions of the elements satisfy 𝑧𝑧 = 0) because the data-independent, it preserves the temporal coherence and sta- use of unfocused waves is still marginal with convex transduc- tistical properties of the real envelopes [14], [15]. The imple- ers [21]. The Cartesian coordinate system associated with a uni- mentation of DAS in the space domain, as opposed to fre- form linear array (ULA) is conventionally defined as follows: quency-based approaches [16], [17], enables the notion of 𝑓𝑓 - the 𝑥𝑥 -axis is parallel to the transducer and points from the left- number and directivity, both of which will be discussed in this most to the rightmost element (Fig. 1), and 𝑥𝑥 = 0 at the center article. In this paper, we put DAS-based beamforming in the of the ULA. The 𝑧𝑧 -axis is perpendicular to the ULA, points context of high-frame-rate imaging before detailing and exem- downward, and 𝑧𝑧 = 0 at the level of the ULA. Note here that plifying its specificities. We then show the effects of the receive what we refer to as a ULA is not necessarily the complete trans- aperture and speed of sound on image quality and propose po- ducer array (i.e. when using a full aperture), but can be a sub- tential solutions to adjust these parameters. In particular, we ex- array (i.e. when using a sub-aperture). Let 𝑝𝑝 denote the pitch of plain how to determine a proper receive aperture (𝑓𝑓 -number) the ULA, i.e. the center-to-center distance between two adja- from the directivity of the array elements. A method is also pro- cent elements, and 𝑁𝑁 the number of elements. The center-to- posed to optimize the speed of sound. It is shown that these two e center distance from the first to the last element of the ULA fundamental aspects can significantly improve the quality of therefore is 𝐿𝐿 = (𝑁𝑁 − 1)𝑝𝑝 . The coordinates of the centers DAS-derived images. Simplified short Matlab codes are pro- (𝑥𝑥 ,𝑧𝑧 ) of the elements are given by vided in the appendix for the sake of clarity and pedagogy. 𝑖𝑖 𝑖𝑖 Complete open-source Matlab codes can be found in the MUST 𝑝𝑝 ( ) 𝑥𝑥 = 2𝑖𝑖 −𝑁𝑁 − 1 𝑖𝑖 e Matlab UltraSound Toolbox released by D. Garcia and down- � with 𝑖𝑖 = 1 …𝑁𝑁 . (1) loadable from www.biomecardio.com/MUST. 𝑧𝑧 = 0 𝑖𝑖 We assume that the speed of sound (= 𝑐𝑐 ) in the medium is 2. DAS in depth uniform. The two-way traveltimes of the wavefronts, from the A. Pulse-echo with a virtual point source transducer to a scatterer of coordinate 𝑿𝑿 = (𝑥𝑥 ,𝑧𝑧 ), and back 𝑠𝑠 𝑠𝑠 𝑠𝑠 High-frame-rate ultrasound is most often used in conjunction from the scatterer to an element #𝑖𝑖 of the transducer, are then with the emission of circular or plane waves. The purpose is to defined by insonify a large region of interest with a wide wavefront so that 𝑑𝑑 (𝑿𝑿 ) +𝑑𝑑 (𝑿𝑿 ,𝑥𝑥 ) TX 𝑠𝑠 RX 𝑠𝑠 𝑖𝑖 a complete image can be obtained with a single transmission. In (2) ( ) 𝜏𝜏 𝑿𝑿 = −𝑡𝑡 . 𝑖𝑖 𝑠𝑠 0 𝑐𝑐 practice, however, a number of transmissions are necessary to allow the corresponding backscattered signals to be combined and 𝑑𝑑 are the transmit and receive distances that are 𝑑𝑑 TX RX to get a high-quality ultrasound image [19]. Assuming that the described below. The parameter 𝑡𝑡 is the start time of the acqui- insonified medium is essentially composed of pointlike Ray- sition. It can be used to reduce the volume of the acquired data leigh backscatters, the latter behave as monopole sources when by considering only the region of interest (often useful in color they are reached by the wavefront. These secondary sources re- Doppler, when the region of interest is a few centimeters from emit the signals quasi-uniformly in all directions (spherical 3 Perrot et al. the probe). A factor that 𝑡𝑡 can take into account is also the dif- shortest distance between the virtual source and the transducer ferent speed of sound through the acoustic lens of the trans- (Fig. 2). The H term reduces this distance to 𝑧𝑧 if |𝑥𝑥 | ≤ 𝐿𝐿 /2 0 0 ducer. The lens adds a supplementary propagation time that can (Fig. 2, left), otherwise it is the distance to the closest end of the be corrected by an additional 𝑡𝑡 ; for instance, Verasonics scan- array segment (Fig. 2, right). The limiting value of circular- ners perform this correction by default for non-custom probe if wave imaging, when 𝛽𝛽 tends towards zero, is plane-wave im- the lens correction is defined by the manufacturer. Finally, 𝑡𝑡 aging. Substituting (4) into (5) and taking the limit as 𝛽𝛽 → 0 can correct delays due to the pulse length. yields the transmit distance for plane-wave imaging (see appen- dix 6.B): The receive distance 𝑑𝑑 represents the distance traveled by RX the spherical wavefront generated by the point scatterer, from 𝐿𝐿 ( ) ( ) ( ) (6) lim 𝑑𝑑 𝑿𝑿 = � sgn 𝜃𝜃 −𝑥𝑥 � sin 𝜃𝜃 +𝑧𝑧 cos(𝜃𝜃 ). TX 𝑠𝑠 𝑠𝑠 𝑠𝑠 the scatterer to element #𝑖𝑖 . It is given by 𝛽𝛽 →0 2 2 (3) ( ) ( ) 𝑑𝑑 𝑿𝑿 ,𝑥𝑥 = �𝑥𝑥 −𝑥𝑥 +𝑧𝑧 . RX 𝑠𝑠 𝑖𝑖 𝑖𝑖 𝑠𝑠 𝑠𝑠 The transmit distance 𝑑𝑑 represents the distance traveled by TX the wavefront from the ULA to a given point scatterer. In high- frame-rate (“ultrafast”) ultrasound imaging, it is common to use circular or plane waves as they can be easily designed. As we will see below, a plane wave is a limiting case of a circular −𝜋𝜋 𝜋𝜋 wave. We define a circular-wave transmit by its tilt 𝜃𝜃 ∈ �, � 2 2 and its angular width 𝛽𝛽 ∈ ]0,𝜋𝜋 [, both represented in Fig. 1. The angle 𝛽𝛽 is the angle between the two lines passing through the virtual source and the two centers of the edge elements. The tilt angle 𝜃𝜃 is measured counterclockwise with respect to the 𝑧𝑧 - axis. Note that this angle 𝜃𝜃 is > 𝜋𝜋 for a focused beam. Trigono- Fig. 3. Hyperbolic signatures in RF signals in the presence of three scatter- metric manipulations can provide the coordinates of the corre- ers. They have an eccentricity of ≈ 𝑐𝑐 . The hyperbolas flatten as the speed of sponding virtual point source: sound increases. The RF signals were simulated by using SIMUS [22]. 𝐿𝐿 sin(2𝜃𝜃 ) 𝑥𝑥 = , 2 sin(𝛽𝛽 ) B. Hyperbolic signatures and DAS (4) 𝐿𝐿 cos(𝛽𝛽 ) + cos(2𝜃𝜃 ) ⎨ The analog signals returned by the piezoelectric elements are 𝑧𝑧 = − . ⎪ 0 2 sin(𝛽𝛽 ) commonly referred to as RF (radio frequency) signals. They are amplified and sampled before being processed by a beam- former. We call 𝓈𝓈 (𝑥𝑥 ,𝑡𝑡 ) the digital signal recorded by the ele- 𝑖𝑖 ment #𝑖𝑖 . The time input 𝑡𝑡 is sometimes called fast-time. The argument 𝑥𝑥 represents the (lateral) 𝑥𝑥 -position of element #𝑖𝑖 𝑖𝑖 given by (1). In the following, we will define 𝓈𝓈 (𝑡𝑡 ) ≡ 𝓈𝓈 (𝑥𝑥 ,𝑡𝑡 ). 𝑖𝑖 𝑖𝑖 These signals 𝓈𝓈 are bandpass modulated signals. The complex 𝑖𝑖 envelope of 𝓈𝓈 , which can be obtained by I/Q demodulation 𝑖𝑖 (downmixing + low-pass filtering; see 6.B in the appendix for . The signal 𝐼𝐼 𝐼𝐼 (𝑡𝑡 ) is com- a short Matlab script), is noted 𝐼𝐼𝐼𝐼 𝑖𝑖 𝑖𝑖 plex, with the real and imaginary parts (𝐼𝐼 (𝑡𝑡 ) and 𝐼𝐼 (𝑡𝑡 )) being Fig. 2. Illustration of the transmit distance 𝑑𝑑 for circular wave imaging, 𝑇𝑇𝑇𝑇 𝑖𝑖 𝑖𝑖 as given by Eq. (5). 𝐿𝐿 stands for the array aperture. (𝑥𝑥 ,𝑧𝑧 ) are the coordi- 𝑠𝑠 𝑠𝑠 the in-phase and quadrature components, respectively. Its mod- nates of one scatterer. 2 2 ulus |𝐼𝐼 𝐼𝐼 (𝑡𝑡 )| = � 𝐼𝐼 (𝑡𝑡 ) +𝐼𝐼 (𝑡𝑡 ) is the real envelope. A B- 𝑖𝑖 𝑖𝑖 𝑖𝑖 mode ultrasound image is obtained by log-compressing real en- In circular-wave imaging, the transmit distance 𝑑𝑑 can be velopes. Additional post-processing is generally used to en- TX written as a function of the virtual point source as follows: hance the images (e.g. speckle filtering [23]). When the signals 𝓈𝓈 (𝑡𝑡 ) recorded by the array elements (𝑖𝑖 = 𝑑𝑑 (𝑿𝑿 ) 𝑖𝑖 TX 𝑠𝑠 1 …𝑁𝑁 ) are stacked side by side in a 𝑡𝑡𝑥𝑥 -plane (𝑥𝑥 − fast time 2 2 = � (𝑥𝑥 −𝑥𝑥 ) + (𝑧𝑧 −𝑧𝑧 ) (5) 𝑠𝑠 0 𝑠𝑠 0 plane), the presence of a scatterer is indicated by a hyperbolic 2 2 −� H(|𝑥𝑥 |−𝐿𝐿 /2 )(|𝑥𝑥 |−𝐿𝐿 /2) +𝑧𝑧 , 0 0 0 signature (Fig. 3). The equation of the hyperbola ℋ related to 𝑠𝑠 ( ) a given scatterer of coordinate 𝑿𝑿 = 𝑥𝑥 ,𝑧𝑧 can be derived by where H represents the Heaviside step function. A generalized 𝑠𝑠 𝑠𝑠 𝑠𝑠 ( ) combining (2) and (3) with 𝑡𝑡 = 𝜏𝜏 𝑿𝑿 +𝑡𝑡 and 𝑥𝑥 = 𝑥𝑥 : equation that includes focused imaging is given in the Appendix 𝑖𝑖 𝑠𝑠 0 𝑖𝑖 6.A. This equation is only valid for points located in the shadow ℋ (𝑿𝑿 ) 𝑠𝑠 𝑠𝑠 of the line of sensor elements (shaded area in Fig. 1). This equa- 𝑑𝑑 (𝑿𝑿 ) TX 𝑠𝑠 �𝑡𝑡 − � (7) ( ) 𝑥𝑥 − 𝑥𝑥 tion also has the condition that the minimum transmission delay 𝑐𝑐 𝑠𝑠 = � (𝑥𝑥 ,𝑡𝑡 ) ∈ ℝ ×ℝ � − − 1 = 0� . 2 2 2 𝑧𝑧 ⁄𝑐𝑐 𝑧𝑧 is zero. The square root term that includes H represents the 𝑠𝑠 𝑠𝑠 4 Perrot et al. All the hyperbolas are geometrically similar, i.e. one hyper- The center frequency 𝑓𝑓 is also the frequency used for 𝑐𝑐 bola can be rescaled and relocated to coincide with another. downmixing during I/Q demodulation. Once beamformed I/Q (𝐼𝐼𝐼𝐼 ) are obtained: 1) a B-mode image can be generated by Their eccentricity is bf taking their log-compressed moduli, or 2) a Doppler image can (8) 𝑒𝑒 = 1 +𝑐𝑐 ≈ 𝑐𝑐 . be produced by analyzing the temporal phase shifts from one frame to another. Only a few correctly selected need to be 𝑖𝑖 Eq. (8) shows that the higher the speed of sound, the flatter summed along the hyperbolas in Eq. (10) to obtain beamformed the hyperbolas are (Fig. 3). DAS is a simplified inverse problem data that can yield high-quality images. We will now see how that consists in determining an amplitude image in the 𝑧𝑧𝑥𝑥 -plane to select an appropriate subset of ⟦1,𝑁𝑁 ⟧, i.e. how to choose a from a set of ultrasonic signals of the 𝑡𝑡𝑥𝑥 -plane. Assuming that proper receive aperture. the medium consists of a very large number of randomly dis- tributed point scatterers, and neglecting multiple scattering, DAS produces an amplitude image as if snapshotting the con- tribution of each reflector emitting simultaneously. This is achieved by adding the amplitudes along the hyperbolas ℋ of 𝑠𝑠 the 𝑡𝑡𝑥𝑥 -plane described by (7). The DAS-beamformed signal value at 𝑿𝑿 = (𝑥𝑥 ,𝑧𝑧 ) is written as 𝑠𝑠 𝑠𝑠 𝑠𝑠 𝓈𝓈 (𝑿𝑿 ) = � 𝓈𝓈 �𝜏𝜏 (𝑿𝑿 )� . bf 𝑠𝑠 𝑖𝑖 𝑖𝑖 𝑠𝑠 (9) 𝑖𝑖 ⊆⟦1,𝑁𝑁 ⟧ The subscript “bf” stands for “beamformed”. As we will see ⟦ ⟧ later, a subset of 1,𝑁𝑁 may be preferred to consider the ele- ment directivity (by using an 𝑓𝑓 -number). In other words, it is often advantageous to discard some signals (i.e. use only some 𝑖𝑖 in 1 …𝑁𝑁 ) when using Eq. (9). 𝑒𝑒 Because the signals are obtained by sampling, the term Fig. 4. The directivity (Eq. 11) of one array element (thick closed curve) depends on the width of the element in relation to the wavelength. The con- ( ) 𝓈𝓈 �𝜏𝜏 𝑿𝑿 � are generally unknown and must be estimated from 𝑖𝑖 𝑖𝑖 𝑠𝑠 ical sectors show areas where the directivity is greater than -3dB. the closest discrete values by using, for example, a 𝑞𝑞 -point in- terpolation; e.g. 𝑞𝑞 = 1 for nearest-neighbor interpolation, 𝑞𝑞 = 2 for linear interpolation, … 𝑞𝑞 = 6 for 5-lobe Lanczos (win- dowed sinc) interpolation [24]. As a recall, a 𝑞𝑞 -point interpola- tor can be written as a weighted arithmetic mean of 𝑞𝑞 data points. These interpolating weights can be conveniently in- cluded in a beamforming sparse matrix, as explained in the sec- tion 2.F entitled “DAS as a matrix product”. C. Beamforming I/Q signals Note that Eq. (9) is valid only for RF signals. In many situa- tions, however, the signals 𝓈𝓈 that the elements record are digi- 𝑖𝑖 tally I/Q demodulated before beamforming. I/Q signals are in- Fig. 5. An appropriate 𝑓𝑓 -number 𝑓𝑓 can be determined from the directivity deed low-frequency signals and are thus easier to handle than of the array elements. The 𝑓𝑓 -number is related to the so-called angle-of- RF. More importantly, the magnitudes of the I and Q signals view (= 2𝛼𝛼 ). The directivity curve and conical sector are those of Fig. 4. contain amplitude and phase information that are classically used to generate B-mode or Doppler images [25]. I/Q signals D. Receive f-number can therefore be beamformed directly onto the image grid (e.g. The signal amplitude along a hyperbolic signature ℋ is not 𝑠𝑠 256×256 size), which cannot be done with RF signals in gen- uniform. It is maximal at its vertex whose abscissa is 𝑥𝑥 = 𝑥𝑥 . 𝑠𝑠 eral. The beamforming of RF signals indeed requires fine axial What essentially governs the signal amplitude along the hyper- sampling to ensure correct envelope detection. bolas is the directivity of the elements. The elements are indeed Summation and demodulation are non-commutating opera- not omnidirectional. As a result, they do not receive uniformly tors. To preserve the relative phases when delaying-and-sum- in all directions. If the ratio (𝑊𝑊 /𝜆𝜆 ) between the element width ming I/Q signals after a demodulation, phase rotators are nec- (𝑊𝑊 ) and the wavelength (𝜆𝜆 ) is relatively large, the element has essary (e.g. see Eq. 4 in [26]). Beamformed I/Q signals can thus a high directivity: it receives mainly in front of it. The smaller be written as: this ratio is, the more omnidirectional it becomes. This ratio is ~1 in most vascular linear arrays and ~0.5 in most cardiac 2 𝑓𝑓 𝜏𝜏 (𝑿𝑿 ) 𝑐𝑐 𝑖𝑖 𝑠𝑠 𝐼𝐼𝐼𝐼 (𝑿𝑿 ) = � 𝐼𝐼𝐼𝐼 �𝜏𝜏 (𝑿𝑿 )� 𝑒𝑒 . bf 𝑠𝑠 𝑖𝑖 𝑖𝑖 𝑠𝑠 (10) 𝑖𝑖 ⊆⟦1,𝑁𝑁 ⟧ 𝑖𝑖𝑖𝑖 𝐼𝐼𝐼𝐼 5 Perrot et al. phased arrays (note that in clinical practice, cardiac phased ar- angle-of-view 𝛼𝛼 in Eq. (13) can thus be determined by using: rays for transthoracic echocardiography are uniform linear ar- ( ) (14) 𝜆𝜆 = 𝜆𝜆 = 𝑐𝑐 /𝑓𝑓 = 𝑐𝑐 / 𝑓𝑓 +𝐵𝐵 /2 . min max 𝑐𝑐 rays). The directivity of an element depends on the wave-prop- −𝜋𝜋 𝜋𝜋 agation angle 𝜗𝜗 ∈ �, � with respect to the 𝑧𝑧 -axis (Fig. 4). For Taking into account the 𝑓𝑓 -number, the DAS equation (10) at 2 2 𝑿𝑿 = (𝑥𝑥 ,𝑧𝑧 ) then becomes: a piston-like element in a soft baffle [27], it is given by 𝑠𝑠 𝑠𝑠 𝑠𝑠 𝑊𝑊 2 𝑓𝑓 𝜏𝜏 (𝑿𝑿 ) 𝑐𝑐 𝑖𝑖 𝑠𝑠 𝐷𝐷 (𝜗𝜗 ) = cos𝜗𝜗 sinc�𝜋𝜋 sin𝜗𝜗 � . (11) 𝐼𝐼𝐼𝐼 (𝑿𝑿 ) = � 𝐼𝐼𝐼𝐼 �𝜏𝜏 (𝑿𝑿 )� 𝑒𝑒 , bf 𝑠𝑠 𝑖𝑖 𝑖𝑖 𝑠𝑠 𝜆𝜆 𝑖𝑖 ⊆⟦1,𝑁𝑁 𝑒𝑒 𝑧𝑧 𝑠𝑠 (15) with 𝑖𝑖 subject to ≥ 𝑓𝑓 . 2 |𝑥𝑥 −𝑥𝑥 | 𝑠𝑠 𝑖𝑖 E. Speed of sound As shown by Eq. (8), the eccentricity of the hyperbolas is al- most equal to the speed of sound 𝑐𝑐 . In other words, the speed of sound governs the shape of the hyperbolas. If the 𝑐𝑐 parameter in DAS does not match the actual speed of sound, the DAS al- gorithm adds the signal amplitudes along incorrect hyperbolas, which distorts the output point spread function. It follows that an error on 𝑐𝑐 may affect image quality significantly (Fig. 6). It is generally assumed that the speed of sound in soft tissues is 1540 m/s on average. In vivo, this may be untrue if the medium is relatively heterogeneous [28]. In vitro, temperature or storage conditions can modify the mechanical properties of a phantom and therefore the speed of sound [29]. To complicate matters, Fig. 6. Effect of the speed of sound on DAS beamforming. Inadequate speeds can lead to distorted (gamma-compressed) B-mode images (left col- several factors might tend to reduce the average speed of sound: umn) because the received signals are summed on erroneous hyperbolic fat tissues, presence of bubbles in the coupling gel, air interface, paths (top right). The proposed 𝒬𝒬 metric analyzes the phase distribution on 𝑝𝑝 curved or bent rays due to aberrations. If the average speed of these hyperbolas (center column). An appropriate speed of sound tends to uniformize the distribution of the phase and maximizes 𝒬𝒬 (bottom right). sound in the insonified medium is not known, the hyperbolas 𝑝𝑝 used during DAS may deviate significantly from the actual hy- perbolas (Fig. 6). As the signal phase is not uniform along in- The sound pressure received by an element is directly propor- correct hyperbolas, this inaccuracy can have negative effects on tional to 𝐷𝐷 (𝜗𝜗 ) and thus has a maximum amplitude when 𝜗𝜗 = 0. the quality of ultrasound images. We propose one solution to It decreases as 𝜗𝜗 tends towards ±𝜋𝜋 /2 (Fig. 3). It follows that estimate an optimal average speed of sound. The idea is to de- the signal-to-noise ratio along a hyperbola ℋ also decreases 𝑠𝑠 termine the speed of sound that returns the hyperbolas with when |𝑥𝑥 −𝑥𝑥 | increases, i.e. when one moves away from the 𝑠𝑠 minimal phase dispersion. The technique that we introduce is vertex (Fig. 3). For this reason, it is recommended to use the top inspired by Yoon et al. [30]. We define the following dimen- part only of the hyperbolas during receive beamforming by us- sionless phase-based quality metric for a given speed of sound: ing an 𝑓𝑓 -number. The receive 𝑓𝑓 -number is defined by the ratio between the depth 𝑧𝑧 and the width of the receive aperture. It is 𝑠𝑠 ( ) � 𝐼𝐼𝐼𝐼 �𝑿𝑿 𝑘𝑘 �� bf 𝑠𝑠 𝒬𝒬 (𝑐𝑐 ) = � . (16) related to the angle-of-view 2𝛼𝛼 (defined in Fig. 5) as follows: 𝑝𝑝 Var(𝜑𝜑 ) 𝑘𝑘 𝑘𝑘 =1…𝑀𝑀 𝑓𝑓 = . (12) # 𝑀𝑀 is the number of spatial points where the signals are beam- ( ) 2 tan 𝛼𝛼 formed. The numerator represents the squared real envelope By using the directivity expression (11), the angle-of-view (intensity) of the signal, after DAS beamforming, at point should satisfy: 𝑿𝑿 (𝑘𝑘 ). The denominator is the sample variance of the un- 𝑠𝑠 th wrapped phase along the 𝑘𝑘 hyperbola ℋ associated with the 𝑠𝑠 ( ) 𝐷𝐷 𝛼𝛼 = 𝐷𝐷 ⟹ 𝑡𝑡 ℎ𝑟𝑟𝑒𝑒𝑠𝑠 ℎ scatterer located at 𝑿𝑿 (𝑘𝑘 ). Note that 𝒬𝒬 can also be calculated (13) 𝑖𝑖 𝑊𝑊 𝑠𝑠 𝑝𝑝 𝛼𝛼 =�𝜗𝜗 ∈ � 0, � � cos𝜗𝜗 sinc�𝜋𝜋 sin𝜗𝜗 � = 𝐷𝐷 � 𝑡𝑡 ℎ𝑟𝑟 ℎ 2 𝜆𝜆 from compound data to return better estimates. In this case, the variances are estimated on compound hyperbolas obtained by On the basis of image quality and SNR (signal-to-noise ratio), coherent summation. We define the expected speed of sound 𝑐𝑐 ̂ a rule-of-thumb compromise is to discard amplitudes below as the one that uniformizes the phases along the hyperbolas, i.e. values of -3dB (i.e. 𝐷𝐷 = 0.71). If 𝑊𝑊 /𝜆𝜆 = 1, the expres- 𝑡𝑡 ℎ𝑟𝑟 ℎ 𝑐𝑐 ̂ maximizes the 𝒬𝒬 metric (Fig. 6): 𝑝𝑝 sions (12) and (13) yield 𝑓𝑓 = 1.2 (see 6.B in the appendix for a Matlab script). As shown by Eq. (11), the element directivity 𝑐𝑐 ̂ = arg max�𝒬𝒬 (𝑐𝑐 )� . 𝑝𝑝 (17) is wavelength-dependent. Consequently, it is preferable to con- 𝑐𝑐 sider the smallest significant wavelength of the signals. For a In [30], only the sum of the phase variances – the denominator signal of center frequency 𝑓𝑓 and bandwidth 𝐵𝐵 (both in Hz), the 𝑐𝑐 𝑒𝑒𝑠𝑠 𝑒𝑒𝑠𝑠 𝑖𝑖𝑖𝑖 6 Perrot et al. in (16) – was minimized, without taking into account the re- envelop detection. For this reason, unless it is strictly necessary spective intensities. The 𝒬𝒬 metric considers these intensities – to post-process RF signals (e.g. in RF-based motion tracking 𝑝𝑝 [31]), this approach is not recommended as it unnecessarily bur- the numerator in (16) – to prioritize the contribution of the dens data storage and calculations. As mentioned earlier, fine bright speckles. The appendix 6.D provides a simplified Matlab axial sampling is indeed required to allow subsequent envelope function (called ezsos) to estimate the speed of sound by max- detection in beamformed RF signals (e.g. by demodulation or imizing (16). Note that we propose to estimate an average speed through a Hilbert transform). This is not the case when beam- of sound (i.e. a constant), not a mapping. The region of interest, forming I/Q. In an extreme situation, it is possible to beamform and its average speed of sound, remain more or less time-invar- I/Q onto a single pixel whereas it is not with RF (how to detect iant when scanning an organ. A single transmit is thus needed. the envelope with a single RF spot?). In short, the original I/Q The I/Q signals can be beamformed at a limited number of can be easily decimated to reduce their number of samples (𝑛𝑛 ), points (say a few hundred). This process has to be done once. 𝑠𝑠 and they require far fewer beamforming points (= 𝑀𝑀 ). This sig- nificantly reduces the number of computational steps (propor- tional to 𝑀𝑀 ×𝑛𝑛 , Fig. 7) when DASing. 𝑠𝑠 It is to be noted that several frames (all with the same virtual source) can be beamformed simultaneously [32]. In such case, becomes a matrix whose each column contains the 𝑛𝑛 𝑁𝑁 I/Q s 𝑒𝑒 data values of each frame. Now, to obtain the interpolated I/Q along the hyperbolas required for estimating the metric 𝒬𝒬 (16), 𝑝𝑝 the signals 𝐼𝐼𝐼𝐼 must be arranged in a matrix of size 𝑖𝑖 , 𝑖𝑖 =1…𝑁𝑁 𝑒𝑒 ( ) 𝑛𝑛 𝑁𝑁 ×𝑁𝑁 , as represented in Fig. 14 of the appendix (see also s e 𝑒𝑒 the Matlab function ezsos in the appendix 6.C). 𝐅𝐅 is large but very sparse since it mostly contains zero- DAS valued elements (Fig. 7). Let 𝑛𝑛𝑛𝑛 𝑧𝑧 denote its number of nonzero elements. If a 𝑞𝑞 -point interpolation is used to estimate the sig- nals at 𝜏𝜏 (𝑿𝑿 ), we have 𝑛𝑛𝑛𝑛 𝑧𝑧 ≤ (𝑀𝑀 𝑁𝑁 𝑞𝑞 ). The right-hand term is 𝑖𝑖 𝑠𝑠 𝑒𝑒 an upper bound since limited apertures (i.e. 𝑓𝑓 > 0) can be used Fig. 7. Example of a sparse DAS matrix. DAS beamforming can be written during receive beamforming, which reduces 𝑛𝑛𝑛𝑛 𝑧𝑧 . The sparsity as a sparse matrix-vector multiplication. For example, the generation of a 256×256 image by beamforming 128 signals each containing 1000 samples [33] of the DAS matrix verifies: leads to a DAS matrix of size (256×256) × (128×1000) = 65536 × 128000. 𝑛𝑛𝑛𝑛 𝑧𝑧 𝑞𝑞 Note that the matrix shown is for illustration purpose only (its aspect can be 𝑠𝑠𝑠𝑠𝑠𝑠𝑝𝑝𝑠𝑠 𝑠𝑠 (𝐅𝐅 ) = 1− ≥ 1− . (19) DAS quite different). 𝑀𝑀 = number of beamformed data (is generally chosen equal 𝑀𝑀 𝑛𝑛 𝑁𝑁 𝑛𝑛 s e 𝑠𝑠 to the number of pixels of the image to be formed); 𝑛𝑛 = number of time 𝑠𝑠 samples per raw signals; 𝑁𝑁 = number of array elements. To give an example, if each array element acquires 1000 I/Q 𝑒𝑒 samples and if the signals 𝐼𝐼𝐼𝐼 are interpolated linearly (i.e. 𝑞𝑞 = 𝑖𝑖 2), then the DAS matrix has a sparsity greater than 1-2/1000 = F. DAS as a matrix product 99.8%, i.e. a density smaller than 0.2%. Sparse matrix-vector As can be seen by Eq. (15), the DAS is a linear operator that multiplication (SpMV) can be computed on GPUs [34]. Its transforms the I/Q signals recorded by the array elements (𝐼𝐼𝐼𝐼 𝑖𝑖 computational complexity is 𝑂𝑂 (𝑛𝑛𝑛𝑛 𝑧𝑧 ) ≤ 𝑂𝑂 (𝑀𝑀 𝑁𝑁 𝑞𝑞 ) [33]. A con- 𝑒𝑒 with 𝑖𝑖 = 1 …𝑁𝑁 ) to beamformed I/Q data (𝐼𝐼𝐼𝐼 ). If the time se- 𝑒𝑒 bf struction of the complex DAS matrix (for I/Q beamforming) is ries 𝐼𝐼𝐼𝐼 contain 𝑛𝑛 samples each, and are stacked in a column 𝑖𝑖 s given in the Matlab function (called ezdas) included in the ap- ( ) vector = �𝐼𝐼 𝐼𝐼 , … ,𝐼𝐼 𝐼𝐼 � of size 𝑛𝑛 𝑁𝑁 × 1 , DAS beam- pendix. The proposed function ezdas uses a linear interpolation, 1 𝑁𝑁 s 𝑒𝑒 which is recommended in most situations. For the sake of com- forming can be written synthetically as: pleteness, the dasmtx function of the MUST toolbox = 𝐅𝐅 × , (18) 𝑏𝑏 𝑓𝑓 DAS (www.biomecardio.com/MUST/) includes several interpola- tion methods: nearest neighbor, linear, quadratic, 5-point least- where 𝐅𝐅 is the DAS matrix (Fig. 7) that contains the in- DAS squares parabolic, and 3-or-5-lobe Lanczos. terpolation weights (see the end of paragraph 2.B), the phase rotators (see equation 10), and the sum truncation induced by a Fig. 7 illustrates the sparsity of such a matrix. If the ultrasound positive f-number (see equation 15). If the data (e.g. B-mode or sequence (array, transmit, receive) and the 𝑀𝑀 beamforming Doppler) to be constructed contain 𝑀𝑀 pixels, then the column point locations remain unchanged, 𝐅𝐅 only needs to be cal- DAS vector that contains the beamformed I/Q signals is of 𝑏𝑏 𝑓𝑓 culated once. Once the DAS matrix is created (or loaded), length 𝑀𝑀 . The matrix 𝐅𝐅 is of size (𝑀𝑀 ×𝑛𝑛 𝑁𝑁 ). It is complex DAS s e SpMV-based beamforming is very fast and compatible with due to the phase rotations present in Eq. (15). The matrix-vector real-time visualization [32]. product for DAS beamforming is schematized in Fig. 7. It is understood that RF signals can also be processed by this matrix product, in which case the DAS matrix is real. Nevertheless, sufficiently fine axial sampling is required to allow subsequent 𝑰𝑰𝑰𝑰 𝑰𝑰𝑰𝑰 𝑰𝑰𝑰𝑰 𝑰𝑰𝑰𝑰 𝑖𝑖𝑡𝑡 𝑰𝑰𝑰𝑰 7 Perrot et al. 3. In vitro and in vivo examples were transmitted at a pulse repetition frequency >6 kHz to cre- ate one coherently compound image. The parameter 𝒬𝒬 was cal- 𝑝𝑝 The respective effects of the speed of sound and 𝑓𝑓 -number on culated by using the compound hyperbolas whose vertices be- image quality were evaluated using the PICMUS dataset avail- longed to a 128×256 beamforming grid. able online [35]. We analyzed RF data that were acquired in an in vitro CIRS phantom (040GSE) and sampled at four times the center frequency. These data were obtained by transmitting steered plane waves with a 128-element 5-MHz linear array. The RF data were IQ-demodulated then beamformed using the DAS matrix of Eq. (18) on the Cartesian image grid specified in PICMUS. A series of real envelopes were constructed from one unsteered plane wave (no transmission delay) for a large range of speed of sound (1400 to 1700 m/s) or 𝑓𝑓 -number (0 to 4). For reference, high-quality real envelopes were also gener- ated by coherent compounding with the whole dataset (75 planes waves with steering ranging from -16 to +16 degrees). The optimal 𝑓𝑓 -number was estimated by using Eq. (12) and (13) with 𝐷𝐷 = 0.71 (-3dB threshold). The optimal speed of 𝑡𝑡 ℎ𝑟𝑟 ℎ Fig. 9. In vitro results obtained in a CIRS phantom using single plane wave sound was determined by maximizing the phase-based metric imaging: effect of the speed of sound and 𝑓𝑓 -number on FWHM and CNR. using Eq. (16). Contrast and lateral resolution were quantified The errors in % are expressed relative to the optimum values (shown by the using the contrast-to-noise ratio (CNR) and full width at half dot). The horizontal and vertical lines represent the optimal speed of sound and 𝑓𝑓 -number estimated by the equations introduced in the manuscript. maximum (FWHM), respectively, as detailed in [35]. The CNR and FWHM reported in our work (Fig. 9) were an average from two cysts and seven nylon wires, respectively (see Figures 2b A. In vitro results and 2d in [35]). From Eq. (14), the minimum significant wavelength trans- mitted by the elements was 𝜆𝜆 ≈ 0.23 mm. It followed that min 𝑊𝑊 /𝜆𝜆 ≈ 1.17. From Eq. (12) and (13), with 𝐷𝐷 = 0.71 min 𝑡𝑡 ℎ𝑟𝑟 ℎ (i.e. -3dB directivity threshold), the directivity-based 𝑓𝑓 -number was thus 𝑓𝑓 ≈ 1.4. Maximizing 𝒬𝒬 (16) by emitting one plane # 𝑝𝑝 wave returned an optimal speed of sound 𝑐𝑐 ̂ = 1570 m/s. Fig. 8. In vitro results obtained in a CIRS phantom (experimental data from PICMUS). The top and bottom rows show the effect of the speed of sound and 𝑓𝑓 -number on CNR (contrast-to-noise ratio) and lateral resolution, re- spectively. The thick vertical bars represent the optimal speed of sound and 𝑓𝑓 -number estimated by the equations introduced in the manuscript. The gray dots represent the metrics after compounding 75 plane waves (reference val- ues), while the black dots are for a single plane wave transmit. Fig. 10. In vivo results obtained in 8 carotids. The curves represent the 𝒬𝒬 𝑝𝑝 metric as a function of the speed of sound. The boxplot shows the distribu- The speed of sound was also measured in the carotids of eight th th tion of the estimated speeds of sound (median, range, 25 , 75 percentiles). healthy volunteers by maximizing the 𝒬𝒬 metric. A cardiologist 𝑝𝑝 used a linear-array transducer (ATL L7-4, center frequency = 5.2 MHz, element width = 0.245 mm, fractional bandwidth As anticipated, both the speed of sound and the 𝑓𝑓 -number at -6dB = 65%) connected to a Verasonics scanner (V-1-128, had a significant impact on the contrast-to-noise ratio (CNR) Verasonics Inc.) to scan the common carotid artery (protocol and lateral resolution (Fig. 8 and Fig. 9). In the CIRS phantom, similar to that described in [36]). Five steered plane waves CNR and lateral resolution were both optimal for 𝑐𝑐 ≈ 1550- (transmit beam angles evenly spaced between -10° and 10°) 𝑒𝑒𝑠𝑠 𝑒𝑒𝑠𝑠 8 Perrot et al. 1600 m/s and 𝑓𝑓 ≈ 1.1-1.6. These values determined using a respective effects are well known, how to set these parameters brute-force search were consistent with those estimated by had not been clearly explained in the literature. We have ad- physical reasoning (phase homogeneity and element directiv- dressed this issue by proposing two simple approaches that can be summarized by equations (12)-(13) and (16)-(17). As we ity) using Eq. (12), (13) and (16), (17); see Fig. 8 and Fig. 9. have argued, for a relatively homogeneous medium, the 𝑓𝑓 -num- ber is essentially related to the directivity of the array elements. B. In vivo results It depends on the transducer (element width, center frequency) The calculated 𝑓𝑓 -number was also 𝑓𝑓 ≈ 1.4 (same transducer and pulse-echo bandwidth. A modification is needed when an- as in vitro). The speeds of sound estimated through maximizing gled receiving is performed, as in vector Doppler [36], [37]. In 𝒬𝒬 in the eight carotids were 1510 ± 41 m/s (Fig. 10). Fig. 11 𝑝𝑝 this case, the angle variable 𝜗𝜗 in Eq. (13) must be replaced by shows a carotid whose speed of sound was estimated at 1400 ( | |), with 𝜃𝜃 being the receive steering angle. 𝜗𝜗 + 𝜃𝜃 𝑇𝑇𝑅𝑅 𝑇𝑇𝑅𝑅 m/s. This value may seem abnormally low. However, this speed An appropriately chosen 𝑓𝑓 -number helps to optimize the bias- of sound probably does not reflect the average speed of sound variance tradeoff when summing the signals along the hyperbo- in the carotid artery. This value has somewhat improved the im- las. An 𝑓𝑓 -number that is too large (too small an aperture) in- age by maximizing the metric we introduced. This apparent duces a summation based on a small sample of elements, thus speed of sound probably reduced the deterioration of the image generating a significant bias. In the case of an 𝑓𝑓 -number that is caused by a combination of independent effects (fat, air, aber- too small (too large an aperture), the summation involves low rations...). The 𝑓𝑓 -number had an obvious positive impact on im- st nd directivities (and therefore low SNRs), thus causing a high var- age quality (Fig. 11, 1 row vs. 2 row). In comparison, the iance. Too high bias or variance affects the contrast and lateral effect of the speed of sound was less noticeable, except for the st nd resolution (Fig. 8). We heuristically found that a -3dB-directiv- change of scale in the 𝑧𝑧 -direction (Fig. 11, 1 column vs. 2 ity threshold gives a good compromise (i.e. Eq. (13) with column). 𝐷𝐷 = 0.71). Similar conclusions could be reached with 𝑡𝑡 ℎ𝑟𝑟 ℎ Doppler or vector Doppler. Fig. 12 illustrates one vector Dop- pler example recomputed from [36] (plane wave imaging, re- ceive angle 𝜃𝜃 = ±15°, 3 cm-diameter disk rotating at 300 𝑇𝑇𝑅𝑅 rpm). The smallest velocity-vector errors coincided with the di- rectivity-derived 𝑓𝑓 -number (≈ 2.6) obtained from Eq. (13) with (𝜗𝜗 + |𝜃𝜃 |) instead of 𝜗𝜗 . 𝑇𝑇𝑅𝑅 Fig. 11. In vivo example in one carotid. The directivity-derived 𝑓𝑓 -number was 1.4 and the estimated speed of sound was 1400 m/s (instead of the as- Fig. 12. Vector Doppler obtained in a rotating disk (data from [36]): 17 un- sumed 1540 m/s). 𝑓𝑓 -number = 0 means “full aperture”. The artefacts related steered plane waves were transmitted; I/Q signals were DAS-beamformed to a full aperture are well visible in the lumen (top row). with ±15° receive angles. The curves show the median relative errors on the 𝑥𝑥 - and 𝑧𝑧 -velocity components as a function of the 𝑓𝑓 -number. The thick ver- tical bars represent the optimal f-number estimated by Eq. (12)-(13). 4. Discussion In this essentially educational article, we have considered the fundamentals of the delay-and-sum (DAS) for ultrasound image beamforming. The underlying theory has been described in depth so that users can be aware of the functioning and limita- tions of this technique. In particular, we have explained that the DAS can be written numerically as a sparse matrix-vector prod- Fig. 13. Knowing the speed of sound (𝑐𝑐 ) allows one to determine which uct and we have clarified the choice of the 𝑓𝑓 -number and speed hyperbola must be chosen. The 𝑓𝑓 -number allows one to decide which part of sound from a physical and practical viewpoint. of the hyperbola should be considered. A. The 𝑓𝑓 -number and speed of sound in DAS To put it simply, the 𝑓𝑓 -number allows one to decide which We have seen that the 𝑓𝑓 -number and speed of sound both in- part of the hyperbola should be considered. Knowing the speed fluence image quality returned by DAS (Fig. 8). Although their 𝑒𝑒𝑠𝑠 9 Perrot et al. of sound, on the other hand, allows one to determine which hy- different stages, for example at the beamforming [48] or com- perbola should be chosen (Fig. 13). According to the in vitro pounding [49] level. If properly trained, DL systems are ex- results obtained in a calibration phantom, the speed of sound pected to provide high-quality images with fewer input data. has a significant impact on the standard CNR and FWHM met- Data-driven DL algorithms, however, are highly dependent on rics (Fig. 8). However, when considering the in vivo results, the accurate and clean training datasets to learn from: for DL sys- overall physiological appearance appears almost unchanged tems to learn to produce high-quality images, they must be (Fig. 11), except for the axial dimensions. This would probably heavily trained with high-quality images. The question then have little or no diagnostic impact in a clinical setting. It then arises as to how to provide such images of optimal quality for a seems that adjusting the speed of sound would be impactful wide variety of organs. If DAS is the chosen beamformer to mainly for calibration or in vitro comparisons. Although this is train DL systems, it is essential to ensure that its parameters are only a speculation, there might nevertheless be some interest in correctly set. super-resolution localization of microbubbles, as each mi- crobubble is expected to be located at the peak or centroid of 5. Conclusion the bubble-PSF (point spread function) [38]. Optimizing the DAS is highly popular in ultrasound imaging. If properly de- speed of sound and thus refining the PSF could make localiza- signed and compound, it can be very effective in terms of in tion more robust. On a more clinical level, it has been shown vivo image quality, even in challenging situations such as the that the longitudinal speed of sound has a potential diagnostic beating heart [50]. Yet DAS is sometimes considered mediocre value in hepatic steatosis, a common liver disease [39], [40]. In by investigations based mainly on in vitro experiments because particular, Imbault et al. [39] estimated the speed of sound by it is often suboptimally addressed. This opinion is often biased locally analyzing the spatial coherence [41]. The measurement by its misuse. For example, the 𝑓𝑓 -number is sometimes ne- of the proposed dimensionless metric 𝒬𝒬 (16) based on I/Q sig- glected or wrongly chosen. It is imperative that the receive 𝑓𝑓 - 𝑝𝑝 nals could provide an alternative method. number be properly selected. In addition, the choice of the speed of sound can have a significant impact on the metrics B. DAS vs. advanced beamformers measured in vitro. To get the most out of DAS and to prevent As stated earlier, DAS is the best-known and simplest beam- any misguided discredit, it is suggested to optimize its parame- forming technique in medical ultrasound imaging. It is fast, ters in order to use it wisely. When it comes to DAS, one must easy to program and well efficient in most situations. Moreover, not DAS out of step. DAS allows for low-complexity real-time beamforming [42]. Yet, several studies have focused on adaptive beamformers 6. Appendix whose intended purpose is to improve image quality for medical A. Generalized transmit distance (focused and circular beams) ultrasound. One typical approach for adaptive beamforming is If (𝑥𝑥 ,𝑧𝑧 ) is a focus point, a reasoning similar to that used to derive to perform a weighted sum and adjust the weights to the data. 0 0 Eq. (5) leads to: This is similar to an ordinary DAS, except that instead of each of the signal samples (inside the receive aperture) contributing 2 2 𝑑𝑑 (𝑿𝑿 ) = sign(𝑧𝑧 −𝑧𝑧 )� (𝑥𝑥 −𝑥𝑥 ) + (𝑧𝑧 −𝑧𝑧 ) TX 𝑠𝑠 𝑠𝑠 0 𝑠𝑠 0 𝑠𝑠 0 (20) equally to the final sum, some data samples contribute more 2 2 + � H(|𝑥𝑥 | +𝐿𝐿 /2 )(|𝑥𝑥 | +𝐿𝐿 /2) +𝑧𝑧 . 0 0 0 than others do, and the relative contributions depend on these Equations (5) and (20) can be combined to yield a generalized transmit samples. Minimum variance (MV) and its variant (eigen-based distance (focusing or diverging wavefronts): MV) as well as 𝑝𝑝 -DAS are ones of these techniques [43]–[45]. (21) Adaptive methods based on the coherence factor [46] also be- 2 2 𝑑𝑑 (𝑿𝑿 ) = sign(𝑧𝑧 −𝑧𝑧 )� (𝑥𝑥 −𝑥𝑥 ) + (𝑧𝑧 −𝑧𝑧 ) TX 𝑠𝑠 𝑠𝑠 0 𝑠𝑠 0 𝑠𝑠 0 long in this category, though the weights are pixelwise. Because 2 2 ( )� (| | ( ) )(| | ( ) ) + sign 𝑧𝑧 H 𝑥𝑥 + sign 𝑧𝑧 𝐿𝐿 /2 𝑥𝑥 + sign 𝑧𝑧 𝐿𝐿 /2 +𝑧𝑧 . they are data-dependent, some adaptive beamforming tech- 0 0 0 0 0 0 niques are computationally expensive and cannot be used in B. Plane wave as a limit of diverging wave real-time. More importantly, although they can improve com- mon metrics (e.g. CNR, FWHM) in specific anechoic-cyst- or Taking the limit of the virtual source coordinates, given by Eq. (4), as 𝛽𝛽 → 0 yields wire-based calibration phantoms, it is not clear if they have sig- nificant value in a clinical context, as to whether they can help ( ) 𝐿𝐿 sin(2𝜃𝜃 ) cos 𝜃𝜃 sin(𝜃𝜃 ) lim 𝑥𝑥 = = 𝐿𝐿 , in a better diagnosis at the patient’s bedside. Another downside ⎪ + 𝛽𝛽 →0 2 𝛽𝛽 𝛽𝛽 (22) of this adaptability is that local phase information, from one 𝐿𝐿 1 + cos(2𝜃𝜃 ) cos (𝜃𝜃 ) lim 𝑧𝑧 = − = −𝐿𝐿 , emission to the next, can be distorted. This may have a detri- ⎪ 0 𝛽𝛽 →0 2 𝛽𝛽 𝛽𝛽 mental impact on velocity estimation by color Doppler: the ( ) ( ) ( ) ( ) Let 𝑓𝑓 𝜃𝜃 = 𝐿𝐿 cos 𝜃𝜃 sin(𝜃𝜃 ) and 𝑔𝑔 𝜃𝜃 = 𝐿𝐿 cos 𝜃𝜃 . The transmit dis- more data-dependent the beamformer, the higher the Doppler tance (5) as 𝛽𝛽 → 0 is thus reduced to error [47]. To avoid biased conclusions when comparing alter- native beamformers against DAS through these metrics, it is lim 𝑑𝑑 (𝑥𝑥 ,𝑧𝑧 ) TX 𝑠𝑠 𝑠𝑠 𝛽𝛽 →0 therefore essential to ensure that the speed of sound and/or the 2 2 (23) = � (𝑥𝑥 −𝑓𝑓 (𝜃𝜃 )⁄𝛽𝛽 ) + (𝑧𝑧 +𝑔𝑔 (𝜃𝜃 )⁄𝛽𝛽 ) 𝑠𝑠 𝑠𝑠 𝑓𝑓 -number have been correctly set for DAS beamforming. 2 2 ( ( ) ( )⁄ ) ( ( )⁄ ) −�sgn 𝜃𝜃 𝑓𝑓 𝜃𝜃 𝛽𝛽 −𝐿𝐿 /2 + 𝑔𝑔 𝜃𝜃 𝛽𝛽 . Another increasingly growing direction for ultrasound image formation is deep learning. Deep learning (DL) can intervene at Expanding the four square terms then factoring yields: 10 Perrot et al. % See also DAS, DASMTX. lim 𝑑𝑑 (𝑥𝑥 ,𝑧𝑧 ) TX 𝑠𝑠 𝑠𝑠 𝛽𝛽 →0 2 2 % -- Damien Garcia -- 2019/11 �𝑓𝑓 (𝜃𝜃 ) +𝑔𝑔 (𝜃𝜃 ) % www.biomecardio.com 𝛽𝛽 2 siz0 = size(x); 2(𝑧𝑧 𝑔𝑔 (𝜃𝜃 )−𝑥𝑥 𝑓𝑓 (𝜃𝜃 )) 𝑥𝑥 +𝑧𝑧 ² 𝑠𝑠 𝑠𝑠 𝑠𝑠 𝑠𝑠 2 (24) � [nl,nc] = size(SIG); �1 + 𝛽𝛽 + 𝛽𝛽 2 2 2 2 𝑓𝑓 (𝜃𝜃 ) +𝑔𝑔 (𝜃𝜃 ) 𝑓𝑓 (𝜃𝜃 ) +𝑔𝑔 (𝜃𝜃 ) x = x(:); z = z(:); % ULA (uniform linear array): sgn(𝜃𝜃 )𝑓𝑓 (𝜃𝜃 )𝐿𝐿 𝐿𝐿 − 1− 𝛽𝛽 + 𝛽𝛽 � . % x-coordinates of the elements 2 2 2 2 ( ) ( ) ( ) ( ) 𝑓𝑓 𝜃𝜃 +𝑔𝑔 𝜃𝜃 4(𝑓𝑓 𝜃𝜃 +𝑔𝑔 𝜃𝜃 ) xe = ((0:nc-1)-(nc-1)/2)*param.pitch; L = xe(end)-xe(1); % length of the array 1/2 ( ) at the first order yields: Using the Taylor series of 1 +𝑥𝑥 % coordinates of the virtual source 𝑓𝑓 (𝜃𝜃 )(sgn(𝜃𝜃 )𝐿𝐿 − 2𝑥𝑥 ) + 2𝑔𝑔 (𝜃𝜃 )𝑧𝑧 x0 = vsource(1); z0 = vsource(2); 𝑠𝑠 𝑠𝑠 lim 𝑑𝑑 (𝑥𝑥 ,𝑧𝑧 ) = (25) TX 𝑠𝑠 𝑠𝑠 𝛽𝛽 →0 2 2 2�𝑓𝑓 (𝜃𝜃 ) +𝑔𝑔 (𝜃𝜃 ) % transmit & receive distances dTX = hypot(x-x0,z-z0)-... −𝜋𝜋 𝜋𝜋 Noting that cos(𝜃𝜃 ) is always positive (since 𝜃𝜃 ∈ �, � ), Eq. (25) can hypot((abs(x0)-L/2)*(abs(x0)>L/2),z0); 2 2 be simplified to obtain Eq. (6) after replacing 𝑓𝑓 (𝜃𝜃 ) and 𝑔𝑔 (𝜃𝜃 ). dRX = hypot(x-xe,z); % two-way travel times C. Determine the f-number in Matlab tau = (dTX+dRX)/param.c; % Note: in Matlab, sinc(x) = sin(pi*x)/(pi*x) f = @(th,width,lambda)... % fast-time indices abs(cos(th)*sinc(width/lambda*sin(th))-0.71); idxt = tau*param.fs + 1; alpha = fminbnd(@(th) f(th,width,lambda),0,pi/2); fnumber = 1/2/tan(alpha); % boolean vectors I = idxt>=1 & idxt<=nl-1; Iaperture = abs(x-xe)<=(z/2/param.fnumber); I = I&Iaperture; D. A simple Matlab code for DAS beamforming % linear indices function [bfSIG,M] = ezdas(SIG,x,z,vsource,param) idx = idxt + (0:nc-1)*nl; idx = idx(I); %EZDAS Delay-And-Sum beamforming idxf = floor(idx); % (Easy version of DAS) idx = idxf-idx; % BFSIG = EZDAS(SIG,X,Z,VSOURCE,PARAM) beamforms % DAS matrix % the RF or I/Q signals stored in the array SIG, [i,~] = find(I); % and returns the beamformed signals BFSIG. The s = [idx+1;-idx]; % (for linear interpolation) % signals are beamformed at the points specified if ~isreal(SIG) % if IQ: phase rotations % by X and Z. s = s.*exp(2i*pi*param.fc*[tau(I);tau(I)]); end % 1) SIG must be a 2-D array. The first dimension M = sparse([i;i],[idxf;idxf+1],s,numel(x),nl*nc); % (i.e. each column) corresponds to single RF % or I/Q signals over (fast-) time, with the % DAS beamforming % 1st column corresponding to the 1st element. bfSIG = reshape(M*SIG(:),siz0); % 2) VSOURCE contains the coordinates [x0,z0] of % the virtual point source. Use large x0,z0 for % plane waves. % 3) PARAM is a structure that contains the % parameter values required for the delay-and- % sum (see below for details). % Note: SIG must be complex for I/Q data % (i.e. SIG = complex(I,Q) = I + 1i*Q). % [~,M] = EZDAS(...) also returns the DAS matrix. % PARAM must contain the following fields: % --------------------------------------- % 1) PARAM.fs: sampling frequency (Hz) % 2) PARAM.pitch: element pitch (m) % 3) PARAM.fc: center frequency (Hz) % 4) PARAM.c: longitudinal velocity (m/s) % 5) PARAM.fnumber: receive f-number Fig. 14. DAS matrix product to obtain interpolated I/Q along the M hyper- bolas ℋ associated to the 𝑀𝑀 pixels of the beamformed image. 𝑛𝑛 = number 𝑠𝑠 𝑠𝑠 % --- of time samples per raw signals; 𝑁𝑁 = number of array elements. 𝑒𝑒 % NOTE #1: Interpolation method: EZDAS uses a % linear interpolation to generate the DAS matrix. % --- E. A simple Matlab code for estimating the speed of sound % NOTE #2: EZDAS is for pedagogical purpose. Use % DAS of the MUST toolbox for more options. function c = ezsos(IQ,x,z,vsource,param) % --- % 11 Perrot et al. %EZSOS Speed-of-sound estimation if wasrow, RF = RF(:); end % (Easy version of SOS) % %-- Time vector % c = EZSOS(IQ,X,Z,VSOURCE,PARAM) returns the nl = size(RF,1); % speed of sound that yields an "optimal" real- t = (0:nl-1)'/Fs; % envelope image by DAS beamforming. The optimal % speed of sound is estimated by analyzing the %-- Downmixing of the RF signals % phases along the diffraction hyperbolas whose IQ = double(RF).*exp(-1i*2*pi*Fc*t); % vertices are located at (X,Z). % %-- Low-pass filter % The input arguments are the same as those of [b,a] = butter(5,.5); % EZDAS. Except that the signals MUST be complex IQ = filtfilt(b,a,IQ)*2; % (I/Q data). Enter "help EZDAS" for details. % %-- Recover the initial size (if was a row vector) % See also SOS, EZDAS. if wasrow, IQ = IQ.'; end % -- Damien Garcia & Vincent Perrot -- 2019/11 % www.biomecardio.com References assert(~isreal(IQ)) [nl,nc] = size(IQ); c0 = param.c; [1] R. E. McKeighen and M. P. Buchin, “New techniques for dynamically variable electronic delays for real time ultrasonic imaging,” in 1977 IQ0 = IQ; Ultrasonics Symposium, 1977, pp. 250–254, doi: IQ = sparse(1:nl*nc,kron(1:nc,ones(1,nl)),... 10.1109/ULTSYM.1977.196834. IQ(:),nl*nc,nc); [2] R. J. Mailloux, “Phased array theory and technology,” Proceedings of the IEEE, vol. 70, no. 3, pp. 246–291, 1982, doi: c = round(fminbnd(@PBC,1200,1700)); 10.1109/PROC.1982.12285. [3] H. T. Friis and C. B. Feldman, “A multiple unit steerable antenna for function Qp = PBC(c) short-wave reception,” Proceedings of the Institute of Radio Engi- param.c = c; neers, vol. 25, no. 7, pp. 841–917, 1937, doi: [~,M] = ezdas(IQ0,x,z/c0*c,vsource,param); 10.1109/JRPROC.1937.228354. [4] P. J. Kahrilas, “HAPDAR—An operational phased array radar,” Pro- hyperb = full(M*IQ); ceedings of the IEEE, vol. 56, no. 11, pp. 1967–1975, 1968, doi: % (contains the diffraction hyperbolas) 10.1109/PROC.1968.6773. [5] P. E. Green, R. A. Frosch, and C. F. Romney, “Principles of an experi- A = angle(hyperb); mental large aperture seismic array (LASA),” Proceedings of the A(A==0) = NaN; IEEE, vol. 53, no. 12, pp. 1821–1833, 1965, doi: A = unwrap(A,[],2); 10.1109/PROC.1965.4453. Qp = abs(mean(hyperb,2,'omitnan')).^2./... [6] Ö. Yilmaz, O. Yilmaz, and S. M. Doherty, “Migration,” in Seismic std(A,0,2,'omitnan').^2; data analysis: processing, inversion, and interpretation of seismic Qp = -mean(Qp,'omitnan'); data, SEG Books, 2001, pp. 463–653. [7] S. H. Gray, J. Etgen, J. Dellinger, and D. Whitmore, “Seismic migra- end tion problems and solutions,” Geophysics, vol. 66, no. 5, pp. 1622– 1640, 2001, doi: 10.1190/1.1487107. end [8] S. D. Pye, S. R. Wild, and W. N. McDicken, “Adaptive time gain compensation for ultrasonic imaging,” Ultrasound in Medicine & Bi- F. A simple Matlab code for I/Q demodulation ology, vol. 18, no. 2, pp. 205–212, 1992, doi: 10.1016/0301- function IQ = ezrf2iq(RF,Fs,Fc) 5629(92)90131-S. [9] A. Kyriakou, E. Neufeld, B. Werner, M. M. Paulides, G. Szekely, and %EZRF2IQ I/Q demodulation of RF data N. Kuster, “A review of numerical and experimental compensation % (easy version of RF2IQ) techniques for skull-induced phase aberrations in transcranial focused ultrasound,” International Journal of Hyperthermia, vol. 30, no. 1, pp. % IQ = RF2IQ(RF,Fs,Fc) demodulates the radio- 36–46, 2014, doi: 10.3109/02656736.2013.861519. % frequency (RF) bandpass signals and returns the [10] S. J. Sanabria, E. Ozkan, M. Rominger, and O. Goksel, “Spatial do- % Inphase/Quadrature (I/Q) components. IQ is a main reconstruction for imaging speed-of-sound with pulse-echo ul- % complex whose real (imaginary) part contains the trasound: simulation and in vivo study,” Phys. Med. Biol., vol. 63, no. % inphase (quadrature) component. 21, p. 215015, 2018, doi: 10.1088/1361-6560/aae2fb. [11] E. Boni, A. C. H. Yu, S. Freear, J. A. Jensen, and P. Tortoli, “Ultra- % 1) Fs is the sampling frequency (in Hz), sound open platforms for next-generation imaging technique develop- % 2) Fc represents the center frequency (in Hz). ment,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Fre- quency Control, vol. 65, no. 7, pp. 1078–1092, 2018, doi: % EZRF2IQ demodulates along columns for 2-D and 10.1109/TUFFC.2018.2844560. % 3-D RF data. Each column corresponds to a single [12] J. Faurie, M. Baudet, J. Poree, G. Cloutier, F. Tournoux, and D. Gar- % RF signal over (fast-) time. cia, “Coupling myocardium and vortex dynamics in diverging-wave echocardiography,” IEEE Trans Ultrason Ferroelectr Freq Control, % EZRF2IQ is for pedagogical purpose. You may use vol. 66, no. 3, pp. 425–432, 2019, doi: % RF2IQ for more options. 10.1109/TUFFC.2018.2842427. [13] O. M. H. Rindal, A. Austeng, A. Fatemi, and A. Rodriguez-Molares, % See also RF2IQ. “The effect of dynamic range alterations in the estimation of contrast,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency % -- Damien Garcia -- 2019/11 Control, vol. 66, no. 7, pp. 1198–1208, 2019, doi: % www.biomecardio.com 10.1109/TUFFC.2019.2911267. [14] F. Destrempes and G. Cloutier, “A critical review and uniformized %-- Convert to column vector (if RF is a row vector) representation of statistical distributions modeling the ultrasound echo wasrow = isrow(RF); 12 Perrot et al. envelope,” Ultrasound in Medicine & Biology, vol. 36, no. 7, pp. [34] M. Maggioni and T. Berger-Wolf, “Optimization techniques for 1037–1051, 2010, doi: 10.1016/j.ultrasmedbio.2010.04.001. sparse matrix–vector multiplication on GPUs,” Journal of Parallel and Distributed Computing, vol. 93–94, pp. 66–86, 2016, doi: [15] S. Salles, H. Liebgott, O. Basset, C. Cachard, D. Vray, and R. Lavarello, “Experimental evaluation of spectral-based quantitative ul- 10.1016/j.jpdc.2016.03.011. trasound imaging using plane wave compounding,” IEEE Transac- [35] H. Liebgott, A. Rodriguez-Molares, F. Cervenansky, J. A. Jensen, and tions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 61, O. Bernard, “Plane-wave imaging challenge in medical ultrasound,” no. 11, pp. 1824–1834, Nov. 2014, doi: IEEE International Ultrasonics Symposium (IUS), pp. 1–4, 2016, doi: 10.1109/ULTSYM.2016.7728908. 10.1109/TUFFC.2014.006543. [16] Jian-Yu Lu, “Experimental study of high frame rate imaging with lim- [36] C. Madiena, J. Faurie, J. Porée, and D. Garcia, “Color and vector flow ited diffraction beams,” IEEE Transactions on Ultrasonics, Ferroelec- imaging in parallel ultrasound with sub-Nyquist sampling,” IEEE trics, and Frequency Control, vol. 45, no. 1, pp. 84–97, 1998, doi: Trans Ultrason Ferroelectr Freq Control, vol. 65, no. 5, pp. 795–802, 10.1109/58.646914. 2018. [37] B. Y. S. Yiu and A. C. H. Yu, “Least-squares multi-angle Doppler es- [17] D. Garcia, L. Le Tarnec, S. Muth, E. Montagnon, J. Porée, and G. Cloutier, “Stolt’s f-k migration for plane wave ultrasound imaging,” timators for plane-wave vector flow imaging,” IEEE Transactions on IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Ultrasonics, Ferroelectrics, and Frequency Control, vol. 63, no. 11, Control, vol. 60, no. 9, pp. 1853–1867, 2013, doi: pp. 1733–1744, 2016, doi: 10.1109/TUFFC.2016.2582514. 10.1109/TUFFC.2013.2771. [38] J. Yu, L. Lavery, and K. Kim, “Super-resolution ultrasound imaging method for microvasculature in vivo with a high temporal accuracy,” [18] D. Garcia, L. Kadem, D. Savéry, P. Pibarot, and L.-G. Durand, “Ana- lytical modeling of the instantaneous maximal transvalvular pressure Scientific Reports, vol. 8, no. 1, pp. 1–11, 2018, doi: 10.1038/s41598- gradient in aortic stenosis,” Journal of Biomechanics, vol. 39, no. 16, 018-32235-2. pp. 3036–3044, 2006, doi: 10.1016/j.jbiomech.2005.10.013. [39] M. Imbault et al., “Robust sound speed estimation for ultrasound- [19] G. Montaldo, M. Tanter, J. Bercoff, N. Benech, and M. Fink, “Coher- based hepatic steatosis assessment,” Phys Med Biol, vol. 62, no. 9, pp. ent plane-wave compounding for very high frame rate ultrasonogra- 3582–3598, 2017, doi: 10.1088/1361-6560/aa6226. phy and transient elastography,” IEEE Transactions on Ultrasonics, [40] R. E. Zubajlo et al., “Experimental validation of longitudinal speed of Ferroelectrics, and Frequency Control, vol. 56, no. 3, pp. 489–506, sound estimates in the diagnosis of hepatic steatosis (part II),” Ultra- 2009, doi: 10.1109/TUFFC.2009.1067. sound in Medicine & Biology, vol. 44, no. 12, pp. 2749–2758, 2018, [20] J. Porée, D. Posada, A. Hodzic, F. Tournoux, G. Cloutier, and D. Gar- doi: 10.1016/j.ultrasmedbio.2018.07.020. cia, “High-frame-rate echocardiography using coherent compounding [41] R. Mallart and M. Fink, “Adaptive focusing in scattering media with Doppler-based motion-compensation,” IEEE Transactions on through sound‐speed inhomogeneities: The van Cittert Zernike ap- Medical Imaging, vol. 35, no. 7, pp. 1647–1657, 2016, doi: proach and focusing criterion,” The Journal of the Acoustical Society 10.1109/TMI.2016.2523346. of America, vol. 96, no. 6, pp. 3721–3732, 1994, doi: [21] S. Bae, P. Kim, and T. Song, “Ultrasonic sector imaging using plane 10.1121/1.410562. wave synthetic focusing with a convex array transducer,” The Journal [42] B. Y. S. Yiu, I. K. H. Tsang, and A. C. H. Yu, “GPU-based beam- of the Acoustical Society of America, vol. 144, no. 5, pp. 2627–2644, former: fast realization of plane wave compounding and synthetic ap- 2018, doi: 10.1121/1.5065391. erture imaging,” IEEE Transactions on Ultrasonics, Ferroelectrics, [22] S. Shahriari and D. Garcia, “Meshfree simulations of ultrasound vec- and Frequency Control, vol. 58, no. 8, pp. 1698–1705, 2011, doi: tor flow imaging using smoothed particle hydrodynamics,” Phys Med 10.1109/TUFFC.2011.1999. Biol, vol. 63, pp. 1–12, 2018, doi: 10.1088/1361-6560/aae3c3. [43] J. F. Synnevag, A. Austeng, and S. Holm, “Adaptive beamforming ap- [23] C. P. Loizou and C. S. Pattichis, “An overview of despeckle-filtering plied to medical ultrasound imaging,” IEEE Transactions on Ultra- techniques,” in Handbook of speckle filtering and tracking in cardio- sonics, Ferroelectrics, and Frequency Control, vol. 54, no. 8, pp. vascular ultrasound imaging and video, IET Digital Library, 2018, pp. 1606–1613, 2007, doi: 10.1109/TUFFC.2007.431. 95–109. [44] B. M. Asl and A. Mahloojifar, “Eigenspace-based minimum variance [24] W. Burger and M. J. Burge, “Interpolation,” in Principles of digital beamforming applied to medical ultrasound imaging,” IEEE Transac- image processing: core algorithms, London: Springer-Verlag, 2009, tions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 57, pp. 210–237. no. 11, pp. 2381–2390, 2010, doi: 10.1109/TUFFC.2010.1706. [25] T. L. Szabo, Diagnostic ultrasound imaging: inside out. Academic [45] M. Polichetti, F. Varray, J.-C. Béra, C. Cachard, and B. Nicolas, “A Press, 2004. nonlinear beamformer based on p-th root compression—application to [26] D. C. M. Horvat, J. S. Bird, and M. M. Goulding, “True time-delay plane wave ultrasound imaging,” Applied Sciences, vol. 8, no. 4, p. bandpass beamforming,” IEEE Journal of Oceanic Engineering, vol. 599, 2018, doi: 10.3390/app8040599. 17, no. 2, pp. 185–192, 1992, doi: 10.1109/48.126975. [46] J. Camacho, M. Parrilla, and C. Fritsch, “Phase coherence imaging,” [27] A. R. Selfridge, G. S. Kino, and B. T. Khuri‐Yakub, “A theory for the IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency radiation pattern of a narrow‐strip acoustic transducer,” Appl. Phys. Control, vol. 56, no. 5, pp. 958–974, 2009, doi: Lett., vol. 37, no. 1, pp. 35–36, 1980, doi: 10.1063/1.91692. 10.1109/TUFFC.2009.1128. [28] D. Napolitano et al., “Sound speed correction in ultrasound imaging,” [47] M. Polichetti, V. Perrot, H. Liebgott, B. Nicolas, and F. Varray, “In- Ultrasonics, vol. 44, pp. e43–e46, 2006, doi: 10.1016/j.ul- fluence of beamforming methods on velocity estimation: in vitro ex- tras.2006.06.061. periments,” in 2018 IEEE International Ultrasonics Symposium (IUS), [29] E. L. Madsen, J. A. Zagzebski, R. A. Banjavie, and R. E. Jutila, “Tis- 2018, pp. 1–4, doi: 10.1109/ULTSYM.2018.8580186. sue mimicking materials for ultrasound phantoms,” Medical Physics, [48] A. C. Luchies and B. C. Byram, “Deep neural networks for ultrasound vol. 5, no. 5, pp. 391–394, 1978, doi: 10.1118/1.594483. beamforming,” IEEE Transactions on Medical Imaging, vol. 37, no. [30] C. Yoon, Y. Lee, J. H. Chang, T. Song, and Y. Yoo, “In vitro estima- 9, pp. 2010–2021, 2018, doi: 10.1109/TMI.2018.2809641. tion of mean sound speed based on minimum average phase variance [49] M. Gasse, F. Millioz, E. Roux, D. Garcia, H. Liebgott, and D. Fribou- in medical ultrasound imaging,” Ultrasonics, vol. 51, no. 7, pp. 795– let, “High-quality plane wave compounding using convolutional neu- 802, 2011, doi: 10.1016/j.ultras.2011.03.007. ral networks,” IEEE Transactions on Ultrasonics, Ferroelectrics, and [31] F. Vignon et al., “Fast frame rate 2D cardiac deformation imaging Frequency Control, vol. 64, no. 10, pp. 1637–1639, 2017, doi: based on RF data: what do we gain?,” in 2017 IEEE International Ul- 10.1109/TUFFC.2017.2736890. trasonics Symposium (IUS), 2017, pp. 1–4, doi: [50] J. Porée, M. Baudet, F. Tournoux, G. Cloutier, and D. Garcia, “A dual 10.1109/ULTSYM.2017.8092648. tissue-Doppler optical-flow method for speckle tracking echocardiog- [32] G. Y. Hou et al., “Sparse matrix beamforming and image reconstruc- raphy at high frame rate,” IEEE Trans Med Imaging, vol. 37, no. 9, tion for 2-D HIFU monitoring using harmonic motion imaging for fo- pp. 2022–2032, 2018, doi: 10.1109/TMI.2018.2811483. cused ultrasound (HMIFU) with in vitro validation,” IEEE Transac- tions on Medical Imaging, vol. 33, no. 11, pp. 2107–2117, 2014, doi: 10.1109/TMI.2014.2332184. [33] J. R. Gilbert, C. Moler, and R. Schreiber, “Sparse matrices in Matlab: design and implementation,” SIAM J. Matrix Anal. Appl., vol. 13, no. 1, pp. 333–356, 1992, doi: 10.1137/0613024. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Electrical Engineering and Systems Science arXiv (Cornell University)

So you think you can DAS? A viewpoint on delay-and-sum beamforming

Loading next page...
 
/lp/arxiv-cornell-university/so-you-think-you-can-das-a-viewpoint-on-delay-and-sum-beamforming-0LmS3lsK5T

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

ISSN
0041-624X
eISSN
ARCH-3348
DOI
10.1016/j.ultras.2020.106309
Publisher site
See Article on Publisher Site

Abstract

Perrot et al. So you think you can DAS? A viewpoint on delay-and-sum beamforming Vincent Perrot, Maxime Polichetti, François Varray, Damien Garcia compared to the size of the receiver. A telecommunication an- Abstract – Delay-and-sum (DAS) is the most widespread digital tenna then receives quasi-planar waves only. An array of beamformer in high-frame-rate ultrasound imaging. Its implementa- equally spaced unit antennas can be steered by phase-shifting tion is simple and compatible with real-time applications. In this view- to collect a radio wavefront from one direction. The radio waves point article, we describe the fundamentals of DAS beamforming. The of the different unit antennas are time-shifted with predefined underlying theory and numerical approach are detailed so that users delays (the greater the delays, the greater the receive angle) so can be aware of its functioning and limitations. In particular, we dis- that they can be summed to increase the signals in the desired cuss the importance of the 𝒇𝒇 -number and speed of sound on image quality, and propose one solution to set their values from a physical receive direction while attenuating those from undesired direc- viewpoint. We suggest determining the 𝒇𝒇 -number from the directivity tions [3]. This describes the basic concept of delay-and-sum. of the transducer elements and the speed of sound from the phase dis- Computer-controlled antenna technology based on DAS has persion of the delayed signals. Simplified Matlab codes are provided been the subject of intensive research for military radars and for the sake of clarity and openness. The effect of the 𝒇𝒇 -number and sonars [4], as it allowed antennas to be quickly steered to detect speed of sound on the lateral resolution and contrast-to-noise ratio was airplanes, submarines, and missiles. The DAS principle had investigated in vitro and in vivo. If not properly preset, these two fac- tors had a substantial negative impact on standard metrics of image also been one of the most useful techniques in seismological quality (namely and 𝐅𝐅𝐅𝐅 ). When beamforming with DAS in research. It allowed observation of small seismic events by con- vitro or in vivo, it is recommended to optimize these parameters in or- necting a large number of horizontally distributed seismic sen- der to use it wisely and prevent image degradation. sors [5]. In geophysics, DAS is often referred to as “diffraction summation”. It was the first computer implementation of seis- Keywords – beamforming, delay-and-sum, 𝒇𝒇 -number, speed of mic data analysis [6]. A similar technique that incorporates am- sound plitude correction and frequency-dependent phase correction, called “Kirchhoff migration”, was then proposed and is still in common use [7]. Amplitude correction due to wave spreading 1. Introduction is not a major concern in ultrasound imaging since the recorded ELAY-AND-SUM is the most basic digital beamformer signals can be time-gain compensated [8]. In addition, unlike D for medical ultrasound imaging. Because of its simplicity geophysical signals, medical ultrasound signals are narrow- and efficiency, it is very likely the most widely used in high- band, so there is little need for frequency-dependent phase cor- frame-rate (ultrafast) ultrasound. Before being used for ultra- rection. An important aspect of geophysics is the wide range of sound imaging [1], the technique of delay-and-sum (DAS) has propagation speeds, which has promoted the development of historically been linked to ground-based and airborne radar as advanced reconstruction techniques [6]. Although the study of well as telecommunication [2]. It originates from steerable ar- speed heterogeneity may be valuable in some ultrasound appli- ray antennas in shortwave communication [3]. Similar to a cations [9], [10], the speed of sound is in most situations con- medical ultrasound transducer, an antenna array is defined as a sidered homogeneous in soft tissue and assumed to be equal to group of connected individual antennas (or elements) that oper- 1540 m/s. As we will see, this fixed value can occasionally be ate together. In radar and ultrasonic imaging, a phased array suboptimal. refers to a computer-controlled array that is electronically steered to transmit or receive in different directions without DAS position in ultrafast ultrasound imaging – With the ac- moving the unit elements. Although the term is generally re- cess to open-architecture ultrasound scanners and their raw served for cardiac application, most transducers used in ultra- data [11], ultrafast ultrasound imaging has quickly gained pop- sound medical imaging are phased arrays since the delays can ularity in our community. Although the content of the manu- be modified electronically in transmission or reception. script could also be derived for focused beams (see Appendix 6.A), we opted for diverging [12] and plane waves to be in line Multidisciplinary DAS – In telecommunication, in contrast to with the most recent literature. “Ultrafast imaging” might be a ultrasound imaging, the involved travel distances are very large misnomer; however, we keep this abusive terminology because This work was supported in part by the LABEX CELYA (ANR-10-LABX- All the authors are with CREATIS, CNRS UMR 5220, INSERM U1206, 0060) and LABEX PRIMES (ANR-10-LABX-0063) of Université de Lyon, Université Lyon 1, INSA Lyon, Université Jean Monnet Saint-Etienne (e-mail: program “Investissements d'Avenir” (ANR-16-IDEX-0005). damien.garcia@inserm.fr; garcia.damien@gmail.com). 𝐅𝐅𝐅𝐅 𝐂𝐂𝐂𝐂𝐂𝐂 2 Perrot et al. unfocused beams are generally used for high-frame-rate pur- waves). Forming an image from sound requires estimating the poses. Even if DAS remains a widespread beamformer in high- round-trip traveltimes of the wavefronts towards and from all frame-rate ultrasound imaging, there is a continuing interest in scatterers. They can be explicitly expressed when transmitting the development of alternative beamforming methods. In par- circular and plane waves, as now explained. ticular, a number of recent studies have focused on adaptive beamformers that aim to improve contrast and spatial resolution provided by DAS. A series of popular adaptive beamforming techniques are discussed in [13]. In such investigations, the DAS is most often chosen as the substandard reference method. However, as it will be explained and illustrated, care must be taken to implement it correctly. In particular, special attention must be paid to the speed of sound and the receive aperture used during beamforming. A well-implemented DAS can have a sig- nificant impact on image quality. It is highly recommended to use it properly when comparing it with other beamformers. Our objective is to dissect the DAS by providing a detailed theoret- ical and pedagogical overview and suggesting some tricks for appropriate and efficient utilization. The advantages of the DAS are many: 1) DAS is based on Fig. 1. Parameters that define a circular-wave transmit. It becomes a 𝜃𝜃 - basic concepts of wave propagation (linearity, straight-ray tilted plane-wave transmit when 𝛽𝛽 → 0 . Figure adapted from [20]. propagation, weak backscattering); 2) its implementation is simple and can be parallelized; 3) it is numerically robust, fast To make things simpler, we focus on rectilinear arrays (for and compatible with real-time applications; and 4) because it is which the positions of the elements satisfy 𝑧𝑧 = 0) because the data-independent, it preserves the temporal coherence and sta- use of unfocused waves is still marginal with convex transduc- tistical properties of the real envelopes [14], [15]. The imple- ers [21]. The Cartesian coordinate system associated with a uni- mentation of DAS in the space domain, as opposed to fre- form linear array (ULA) is conventionally defined as follows: quency-based approaches [16], [17], enables the notion of 𝑓𝑓 - the 𝑥𝑥 -axis is parallel to the transducer and points from the left- number and directivity, both of which will be discussed in this most to the rightmost element (Fig. 1), and 𝑥𝑥 = 0 at the center article. In this paper, we put DAS-based beamforming in the of the ULA. The 𝑧𝑧 -axis is perpendicular to the ULA, points context of high-frame-rate imaging before detailing and exem- downward, and 𝑧𝑧 = 0 at the level of the ULA. Note here that plifying its specificities. We then show the effects of the receive what we refer to as a ULA is not necessarily the complete trans- aperture and speed of sound on image quality and propose po- ducer array (i.e. when using a full aperture), but can be a sub- tential solutions to adjust these parameters. In particular, we ex- array (i.e. when using a sub-aperture). Let 𝑝𝑝 denote the pitch of plain how to determine a proper receive aperture (𝑓𝑓 -number) the ULA, i.e. the center-to-center distance between two adja- from the directivity of the array elements. A method is also pro- cent elements, and 𝑁𝑁 the number of elements. The center-to- posed to optimize the speed of sound. It is shown that these two e center distance from the first to the last element of the ULA fundamental aspects can significantly improve the quality of therefore is 𝐿𝐿 = (𝑁𝑁 − 1)𝑝𝑝 . The coordinates of the centers DAS-derived images. Simplified short Matlab codes are pro- (𝑥𝑥 ,𝑧𝑧 ) of the elements are given by vided in the appendix for the sake of clarity and pedagogy. 𝑖𝑖 𝑖𝑖 Complete open-source Matlab codes can be found in the MUST 𝑝𝑝 ( ) 𝑥𝑥 = 2𝑖𝑖 −𝑁𝑁 − 1 𝑖𝑖 e Matlab UltraSound Toolbox released by D. Garcia and down- � with 𝑖𝑖 = 1 …𝑁𝑁 . (1) loadable from www.biomecardio.com/MUST. 𝑧𝑧 = 0 𝑖𝑖 We assume that the speed of sound (= 𝑐𝑐 ) in the medium is 2. DAS in depth uniform. The two-way traveltimes of the wavefronts, from the A. Pulse-echo with a virtual point source transducer to a scatterer of coordinate 𝑿𝑿 = (𝑥𝑥 ,𝑧𝑧 ), and back 𝑠𝑠 𝑠𝑠 𝑠𝑠 High-frame-rate ultrasound is most often used in conjunction from the scatterer to an element #𝑖𝑖 of the transducer, are then with the emission of circular or plane waves. The purpose is to defined by insonify a large region of interest with a wide wavefront so that 𝑑𝑑 (𝑿𝑿 ) +𝑑𝑑 (𝑿𝑿 ,𝑥𝑥 ) TX 𝑠𝑠 RX 𝑠𝑠 𝑖𝑖 a complete image can be obtained with a single transmission. In (2) ( ) 𝜏𝜏 𝑿𝑿 = −𝑡𝑡 . 𝑖𝑖 𝑠𝑠 0 𝑐𝑐 practice, however, a number of transmissions are necessary to allow the corresponding backscattered signals to be combined and 𝑑𝑑 are the transmit and receive distances that are 𝑑𝑑 TX RX to get a high-quality ultrasound image [19]. Assuming that the described below. The parameter 𝑡𝑡 is the start time of the acqui- insonified medium is essentially composed of pointlike Ray- sition. It can be used to reduce the volume of the acquired data leigh backscatters, the latter behave as monopole sources when by considering only the region of interest (often useful in color they are reached by the wavefront. These secondary sources re- Doppler, when the region of interest is a few centimeters from emit the signals quasi-uniformly in all directions (spherical 3 Perrot et al. the probe). A factor that 𝑡𝑡 can take into account is also the dif- shortest distance between the virtual source and the transducer ferent speed of sound through the acoustic lens of the trans- (Fig. 2). The H term reduces this distance to 𝑧𝑧 if |𝑥𝑥 | ≤ 𝐿𝐿 /2 0 0 ducer. The lens adds a supplementary propagation time that can (Fig. 2, left), otherwise it is the distance to the closest end of the be corrected by an additional 𝑡𝑡 ; for instance, Verasonics scan- array segment (Fig. 2, right). The limiting value of circular- ners perform this correction by default for non-custom probe if wave imaging, when 𝛽𝛽 tends towards zero, is plane-wave im- the lens correction is defined by the manufacturer. Finally, 𝑡𝑡 aging. Substituting (4) into (5) and taking the limit as 𝛽𝛽 → 0 can correct delays due to the pulse length. yields the transmit distance for plane-wave imaging (see appen- dix 6.B): The receive distance 𝑑𝑑 represents the distance traveled by RX the spherical wavefront generated by the point scatterer, from 𝐿𝐿 ( ) ( ) ( ) (6) lim 𝑑𝑑 𝑿𝑿 = � sgn 𝜃𝜃 −𝑥𝑥 � sin 𝜃𝜃 +𝑧𝑧 cos(𝜃𝜃 ). TX 𝑠𝑠 𝑠𝑠 𝑠𝑠 the scatterer to element #𝑖𝑖 . It is given by 𝛽𝛽 →0 2 2 (3) ( ) ( ) 𝑑𝑑 𝑿𝑿 ,𝑥𝑥 = �𝑥𝑥 −𝑥𝑥 +𝑧𝑧 . RX 𝑠𝑠 𝑖𝑖 𝑖𝑖 𝑠𝑠 𝑠𝑠 The transmit distance 𝑑𝑑 represents the distance traveled by TX the wavefront from the ULA to a given point scatterer. In high- frame-rate (“ultrafast”) ultrasound imaging, it is common to use circular or plane waves as they can be easily designed. As we will see below, a plane wave is a limiting case of a circular −𝜋𝜋 𝜋𝜋 wave. We define a circular-wave transmit by its tilt 𝜃𝜃 ∈ �, � 2 2 and its angular width 𝛽𝛽 ∈ ]0,𝜋𝜋 [, both represented in Fig. 1. The angle 𝛽𝛽 is the angle between the two lines passing through the virtual source and the two centers of the edge elements. The tilt angle 𝜃𝜃 is measured counterclockwise with respect to the 𝑧𝑧 - axis. Note that this angle 𝜃𝜃 is > 𝜋𝜋 for a focused beam. Trigono- Fig. 3. Hyperbolic signatures in RF signals in the presence of three scatter- metric manipulations can provide the coordinates of the corre- ers. They have an eccentricity of ≈ 𝑐𝑐 . The hyperbolas flatten as the speed of sponding virtual point source: sound increases. The RF signals were simulated by using SIMUS [22]. 𝐿𝐿 sin(2𝜃𝜃 ) 𝑥𝑥 = , 2 sin(𝛽𝛽 ) B. Hyperbolic signatures and DAS (4) 𝐿𝐿 cos(𝛽𝛽 ) + cos(2𝜃𝜃 ) ⎨ The analog signals returned by the piezoelectric elements are 𝑧𝑧 = − . ⎪ 0 2 sin(𝛽𝛽 ) commonly referred to as RF (radio frequency) signals. They are amplified and sampled before being processed by a beam- former. We call 𝓈𝓈 (𝑥𝑥 ,𝑡𝑡 ) the digital signal recorded by the ele- 𝑖𝑖 ment #𝑖𝑖 . The time input 𝑡𝑡 is sometimes called fast-time. The argument 𝑥𝑥 represents the (lateral) 𝑥𝑥 -position of element #𝑖𝑖 𝑖𝑖 given by (1). In the following, we will define 𝓈𝓈 (𝑡𝑡 ) ≡ 𝓈𝓈 (𝑥𝑥 ,𝑡𝑡 ). 𝑖𝑖 𝑖𝑖 These signals 𝓈𝓈 are bandpass modulated signals. The complex 𝑖𝑖 envelope of 𝓈𝓈 , which can be obtained by I/Q demodulation 𝑖𝑖 (downmixing + low-pass filtering; see 6.B in the appendix for . The signal 𝐼𝐼 𝐼𝐼 (𝑡𝑡 ) is com- a short Matlab script), is noted 𝐼𝐼𝐼𝐼 𝑖𝑖 𝑖𝑖 plex, with the real and imaginary parts (𝐼𝐼 (𝑡𝑡 ) and 𝐼𝐼 (𝑡𝑡 )) being Fig. 2. Illustration of the transmit distance 𝑑𝑑 for circular wave imaging, 𝑇𝑇𝑇𝑇 𝑖𝑖 𝑖𝑖 as given by Eq. (5). 𝐿𝐿 stands for the array aperture. (𝑥𝑥 ,𝑧𝑧 ) are the coordi- 𝑠𝑠 𝑠𝑠 the in-phase and quadrature components, respectively. Its mod- nates of one scatterer. 2 2 ulus |𝐼𝐼 𝐼𝐼 (𝑡𝑡 )| = � 𝐼𝐼 (𝑡𝑡 ) +𝐼𝐼 (𝑡𝑡 ) is the real envelope. A B- 𝑖𝑖 𝑖𝑖 𝑖𝑖 mode ultrasound image is obtained by log-compressing real en- In circular-wave imaging, the transmit distance 𝑑𝑑 can be velopes. Additional post-processing is generally used to en- TX written as a function of the virtual point source as follows: hance the images (e.g. speckle filtering [23]). When the signals 𝓈𝓈 (𝑡𝑡 ) recorded by the array elements (𝑖𝑖 = 𝑑𝑑 (𝑿𝑿 ) 𝑖𝑖 TX 𝑠𝑠 1 …𝑁𝑁 ) are stacked side by side in a 𝑡𝑡𝑥𝑥 -plane (𝑥𝑥 − fast time 2 2 = � (𝑥𝑥 −𝑥𝑥 ) + (𝑧𝑧 −𝑧𝑧 ) (5) 𝑠𝑠 0 𝑠𝑠 0 plane), the presence of a scatterer is indicated by a hyperbolic 2 2 −� H(|𝑥𝑥 |−𝐿𝐿 /2 )(|𝑥𝑥 |−𝐿𝐿 /2) +𝑧𝑧 , 0 0 0 signature (Fig. 3). The equation of the hyperbola ℋ related to 𝑠𝑠 ( ) a given scatterer of coordinate 𝑿𝑿 = 𝑥𝑥 ,𝑧𝑧 can be derived by where H represents the Heaviside step function. A generalized 𝑠𝑠 𝑠𝑠 𝑠𝑠 ( ) combining (2) and (3) with 𝑡𝑡 = 𝜏𝜏 𝑿𝑿 +𝑡𝑡 and 𝑥𝑥 = 𝑥𝑥 : equation that includes focused imaging is given in the Appendix 𝑖𝑖 𝑠𝑠 0 𝑖𝑖 6.A. This equation is only valid for points located in the shadow ℋ (𝑿𝑿 ) 𝑠𝑠 𝑠𝑠 of the line of sensor elements (shaded area in Fig. 1). This equa- 𝑑𝑑 (𝑿𝑿 ) TX 𝑠𝑠 �𝑡𝑡 − � (7) ( ) 𝑥𝑥 − 𝑥𝑥 tion also has the condition that the minimum transmission delay 𝑐𝑐 𝑠𝑠 = � (𝑥𝑥 ,𝑡𝑡 ) ∈ ℝ ×ℝ � − − 1 = 0� . 2 2 2 𝑧𝑧 ⁄𝑐𝑐 𝑧𝑧 is zero. The square root term that includes H represents the 𝑠𝑠 𝑠𝑠 4 Perrot et al. All the hyperbolas are geometrically similar, i.e. one hyper- The center frequency 𝑓𝑓 is also the frequency used for 𝑐𝑐 bola can be rescaled and relocated to coincide with another. downmixing during I/Q demodulation. Once beamformed I/Q (𝐼𝐼𝐼𝐼 ) are obtained: 1) a B-mode image can be generated by Their eccentricity is bf taking their log-compressed moduli, or 2) a Doppler image can (8) 𝑒𝑒 = 1 +𝑐𝑐 ≈ 𝑐𝑐 . be produced by analyzing the temporal phase shifts from one frame to another. Only a few correctly selected need to be 𝑖𝑖 Eq. (8) shows that the higher the speed of sound, the flatter summed along the hyperbolas in Eq. (10) to obtain beamformed the hyperbolas are (Fig. 3). DAS is a simplified inverse problem data that can yield high-quality images. We will now see how that consists in determining an amplitude image in the 𝑧𝑧𝑥𝑥 -plane to select an appropriate subset of ⟦1,𝑁𝑁 ⟧, i.e. how to choose a from a set of ultrasonic signals of the 𝑡𝑡𝑥𝑥 -plane. Assuming that proper receive aperture. the medium consists of a very large number of randomly dis- tributed point scatterers, and neglecting multiple scattering, DAS produces an amplitude image as if snapshotting the con- tribution of each reflector emitting simultaneously. This is achieved by adding the amplitudes along the hyperbolas ℋ of 𝑠𝑠 the 𝑡𝑡𝑥𝑥 -plane described by (7). The DAS-beamformed signal value at 𝑿𝑿 = (𝑥𝑥 ,𝑧𝑧 ) is written as 𝑠𝑠 𝑠𝑠 𝑠𝑠 𝓈𝓈 (𝑿𝑿 ) = � 𝓈𝓈 �𝜏𝜏 (𝑿𝑿 )� . bf 𝑠𝑠 𝑖𝑖 𝑖𝑖 𝑠𝑠 (9) 𝑖𝑖 ⊆⟦1,𝑁𝑁 ⟧ The subscript “bf” stands for “beamformed”. As we will see ⟦ ⟧ later, a subset of 1,𝑁𝑁 may be preferred to consider the ele- ment directivity (by using an 𝑓𝑓 -number). In other words, it is often advantageous to discard some signals (i.e. use only some 𝑖𝑖 in 1 …𝑁𝑁 ) when using Eq. (9). 𝑒𝑒 Because the signals are obtained by sampling, the term Fig. 4. The directivity (Eq. 11) of one array element (thick closed curve) depends on the width of the element in relation to the wavelength. The con- ( ) 𝓈𝓈 �𝜏𝜏 𝑿𝑿 � are generally unknown and must be estimated from 𝑖𝑖 𝑖𝑖 𝑠𝑠 ical sectors show areas where the directivity is greater than -3dB. the closest discrete values by using, for example, a 𝑞𝑞 -point in- terpolation; e.g. 𝑞𝑞 = 1 for nearest-neighbor interpolation, 𝑞𝑞 = 2 for linear interpolation, … 𝑞𝑞 = 6 for 5-lobe Lanczos (win- dowed sinc) interpolation [24]. As a recall, a 𝑞𝑞 -point interpola- tor can be written as a weighted arithmetic mean of 𝑞𝑞 data points. These interpolating weights can be conveniently in- cluded in a beamforming sparse matrix, as explained in the sec- tion 2.F entitled “DAS as a matrix product”. C. Beamforming I/Q signals Note that Eq. (9) is valid only for RF signals. In many situa- tions, however, the signals 𝓈𝓈 that the elements record are digi- 𝑖𝑖 tally I/Q demodulated before beamforming. I/Q signals are in- Fig. 5. An appropriate 𝑓𝑓 -number 𝑓𝑓 can be determined from the directivity deed low-frequency signals and are thus easier to handle than of the array elements. The 𝑓𝑓 -number is related to the so-called angle-of- RF. More importantly, the magnitudes of the I and Q signals view (= 2𝛼𝛼 ). The directivity curve and conical sector are those of Fig. 4. contain amplitude and phase information that are classically used to generate B-mode or Doppler images [25]. I/Q signals D. Receive f-number can therefore be beamformed directly onto the image grid (e.g. The signal amplitude along a hyperbolic signature ℋ is not 𝑠𝑠 256×256 size), which cannot be done with RF signals in gen- uniform. It is maximal at its vertex whose abscissa is 𝑥𝑥 = 𝑥𝑥 . 𝑠𝑠 eral. The beamforming of RF signals indeed requires fine axial What essentially governs the signal amplitude along the hyper- sampling to ensure correct envelope detection. bolas is the directivity of the elements. The elements are indeed Summation and demodulation are non-commutating opera- not omnidirectional. As a result, they do not receive uniformly tors. To preserve the relative phases when delaying-and-sum- in all directions. If the ratio (𝑊𝑊 /𝜆𝜆 ) between the element width ming I/Q signals after a demodulation, phase rotators are nec- (𝑊𝑊 ) and the wavelength (𝜆𝜆 ) is relatively large, the element has essary (e.g. see Eq. 4 in [26]). Beamformed I/Q signals can thus a high directivity: it receives mainly in front of it. The smaller be written as: this ratio is, the more omnidirectional it becomes. This ratio is ~1 in most vascular linear arrays and ~0.5 in most cardiac 2 𝑓𝑓 𝜏𝜏 (𝑿𝑿 ) 𝑐𝑐 𝑖𝑖 𝑠𝑠 𝐼𝐼𝐼𝐼 (𝑿𝑿 ) = � 𝐼𝐼𝐼𝐼 �𝜏𝜏 (𝑿𝑿 )� 𝑒𝑒 . bf 𝑠𝑠 𝑖𝑖 𝑖𝑖 𝑠𝑠 (10) 𝑖𝑖 ⊆⟦1,𝑁𝑁 ⟧ 𝑖𝑖𝑖𝑖 𝐼𝐼𝐼𝐼 5 Perrot et al. phased arrays (note that in clinical practice, cardiac phased ar- angle-of-view 𝛼𝛼 in Eq. (13) can thus be determined by using: rays for transthoracic echocardiography are uniform linear ar- ( ) (14) 𝜆𝜆 = 𝜆𝜆 = 𝑐𝑐 /𝑓𝑓 = 𝑐𝑐 / 𝑓𝑓 +𝐵𝐵 /2 . min max 𝑐𝑐 rays). The directivity of an element depends on the wave-prop- −𝜋𝜋 𝜋𝜋 agation angle 𝜗𝜗 ∈ �, � with respect to the 𝑧𝑧 -axis (Fig. 4). For Taking into account the 𝑓𝑓 -number, the DAS equation (10) at 2 2 𝑿𝑿 = (𝑥𝑥 ,𝑧𝑧 ) then becomes: a piston-like element in a soft baffle [27], it is given by 𝑠𝑠 𝑠𝑠 𝑠𝑠 𝑊𝑊 2 𝑓𝑓 𝜏𝜏 (𝑿𝑿 ) 𝑐𝑐 𝑖𝑖 𝑠𝑠 𝐷𝐷 (𝜗𝜗 ) = cos𝜗𝜗 sinc�𝜋𝜋 sin𝜗𝜗 � . (11) 𝐼𝐼𝐼𝐼 (𝑿𝑿 ) = � 𝐼𝐼𝐼𝐼 �𝜏𝜏 (𝑿𝑿 )� 𝑒𝑒 , bf 𝑠𝑠 𝑖𝑖 𝑖𝑖 𝑠𝑠 𝜆𝜆 𝑖𝑖 ⊆⟦1,𝑁𝑁 𝑒𝑒 𝑧𝑧 𝑠𝑠 (15) with 𝑖𝑖 subject to ≥ 𝑓𝑓 . 2 |𝑥𝑥 −𝑥𝑥 | 𝑠𝑠 𝑖𝑖 E. Speed of sound As shown by Eq. (8), the eccentricity of the hyperbolas is al- most equal to the speed of sound 𝑐𝑐 . In other words, the speed of sound governs the shape of the hyperbolas. If the 𝑐𝑐 parameter in DAS does not match the actual speed of sound, the DAS al- gorithm adds the signal amplitudes along incorrect hyperbolas, which distorts the output point spread function. It follows that an error on 𝑐𝑐 may affect image quality significantly (Fig. 6). It is generally assumed that the speed of sound in soft tissues is 1540 m/s on average. In vivo, this may be untrue if the medium is relatively heterogeneous [28]. In vitro, temperature or storage conditions can modify the mechanical properties of a phantom and therefore the speed of sound [29]. To complicate matters, Fig. 6. Effect of the speed of sound on DAS beamforming. Inadequate speeds can lead to distorted (gamma-compressed) B-mode images (left col- several factors might tend to reduce the average speed of sound: umn) because the received signals are summed on erroneous hyperbolic fat tissues, presence of bubbles in the coupling gel, air interface, paths (top right). The proposed 𝒬𝒬 metric analyzes the phase distribution on 𝑝𝑝 curved or bent rays due to aberrations. If the average speed of these hyperbolas (center column). An appropriate speed of sound tends to uniformize the distribution of the phase and maximizes 𝒬𝒬 (bottom right). sound in the insonified medium is not known, the hyperbolas 𝑝𝑝 used during DAS may deviate significantly from the actual hy- perbolas (Fig. 6). As the signal phase is not uniform along in- The sound pressure received by an element is directly propor- correct hyperbolas, this inaccuracy can have negative effects on tional to 𝐷𝐷 (𝜗𝜗 ) and thus has a maximum amplitude when 𝜗𝜗 = 0. the quality of ultrasound images. We propose one solution to It decreases as 𝜗𝜗 tends towards ±𝜋𝜋 /2 (Fig. 3). It follows that estimate an optimal average speed of sound. The idea is to de- the signal-to-noise ratio along a hyperbola ℋ also decreases 𝑠𝑠 termine the speed of sound that returns the hyperbolas with when |𝑥𝑥 −𝑥𝑥 | increases, i.e. when one moves away from the 𝑠𝑠 minimal phase dispersion. The technique that we introduce is vertex (Fig. 3). For this reason, it is recommended to use the top inspired by Yoon et al. [30]. We define the following dimen- part only of the hyperbolas during receive beamforming by us- sionless phase-based quality metric for a given speed of sound: ing an 𝑓𝑓 -number. The receive 𝑓𝑓 -number is defined by the ratio between the depth 𝑧𝑧 and the width of the receive aperture. It is 𝑠𝑠 ( ) � 𝐼𝐼𝐼𝐼 �𝑿𝑿 𝑘𝑘 �� bf 𝑠𝑠 𝒬𝒬 (𝑐𝑐 ) = � . (16) related to the angle-of-view 2𝛼𝛼 (defined in Fig. 5) as follows: 𝑝𝑝 Var(𝜑𝜑 ) 𝑘𝑘 𝑘𝑘 =1…𝑀𝑀 𝑓𝑓 = . (12) # 𝑀𝑀 is the number of spatial points where the signals are beam- ( ) 2 tan 𝛼𝛼 formed. The numerator represents the squared real envelope By using the directivity expression (11), the angle-of-view (intensity) of the signal, after DAS beamforming, at point should satisfy: 𝑿𝑿 (𝑘𝑘 ). The denominator is the sample variance of the un- 𝑠𝑠 th wrapped phase along the 𝑘𝑘 hyperbola ℋ associated with the 𝑠𝑠 ( ) 𝐷𝐷 𝛼𝛼 = 𝐷𝐷 ⟹ 𝑡𝑡 ℎ𝑟𝑟𝑒𝑒𝑠𝑠 ℎ scatterer located at 𝑿𝑿 (𝑘𝑘 ). Note that 𝒬𝒬 can also be calculated (13) 𝑖𝑖 𝑊𝑊 𝑠𝑠 𝑝𝑝 𝛼𝛼 =�𝜗𝜗 ∈ � 0, � � cos𝜗𝜗 sinc�𝜋𝜋 sin𝜗𝜗 � = 𝐷𝐷 � 𝑡𝑡 ℎ𝑟𝑟 ℎ 2 𝜆𝜆 from compound data to return better estimates. In this case, the variances are estimated on compound hyperbolas obtained by On the basis of image quality and SNR (signal-to-noise ratio), coherent summation. We define the expected speed of sound 𝑐𝑐 ̂ a rule-of-thumb compromise is to discard amplitudes below as the one that uniformizes the phases along the hyperbolas, i.e. values of -3dB (i.e. 𝐷𝐷 = 0.71). If 𝑊𝑊 /𝜆𝜆 = 1, the expres- 𝑡𝑡 ℎ𝑟𝑟 ℎ 𝑐𝑐 ̂ maximizes the 𝒬𝒬 metric (Fig. 6): 𝑝𝑝 sions (12) and (13) yield 𝑓𝑓 = 1.2 (see 6.B in the appendix for a Matlab script). As shown by Eq. (11), the element directivity 𝑐𝑐 ̂ = arg max�𝒬𝒬 (𝑐𝑐 )� . 𝑝𝑝 (17) is wavelength-dependent. Consequently, it is preferable to con- 𝑐𝑐 sider the smallest significant wavelength of the signals. For a In [30], only the sum of the phase variances – the denominator signal of center frequency 𝑓𝑓 and bandwidth 𝐵𝐵 (both in Hz), the 𝑐𝑐 𝑒𝑒𝑠𝑠 𝑒𝑒𝑠𝑠 𝑖𝑖𝑖𝑖 6 Perrot et al. in (16) – was minimized, without taking into account the re- envelop detection. For this reason, unless it is strictly necessary spective intensities. The 𝒬𝒬 metric considers these intensities – to post-process RF signals (e.g. in RF-based motion tracking 𝑝𝑝 [31]), this approach is not recommended as it unnecessarily bur- the numerator in (16) – to prioritize the contribution of the dens data storage and calculations. As mentioned earlier, fine bright speckles. The appendix 6.D provides a simplified Matlab axial sampling is indeed required to allow subsequent envelope function (called ezsos) to estimate the speed of sound by max- detection in beamformed RF signals (e.g. by demodulation or imizing (16). Note that we propose to estimate an average speed through a Hilbert transform). This is not the case when beam- of sound (i.e. a constant), not a mapping. The region of interest, forming I/Q. In an extreme situation, it is possible to beamform and its average speed of sound, remain more or less time-invar- I/Q onto a single pixel whereas it is not with RF (how to detect iant when scanning an organ. A single transmit is thus needed. the envelope with a single RF spot?). In short, the original I/Q The I/Q signals can be beamformed at a limited number of can be easily decimated to reduce their number of samples (𝑛𝑛 ), points (say a few hundred). This process has to be done once. 𝑠𝑠 and they require far fewer beamforming points (= 𝑀𝑀 ). This sig- nificantly reduces the number of computational steps (propor- tional to 𝑀𝑀 ×𝑛𝑛 , Fig. 7) when DASing. 𝑠𝑠 It is to be noted that several frames (all with the same virtual source) can be beamformed simultaneously [32]. In such case, becomes a matrix whose each column contains the 𝑛𝑛 𝑁𝑁 I/Q s 𝑒𝑒 data values of each frame. Now, to obtain the interpolated I/Q along the hyperbolas required for estimating the metric 𝒬𝒬 (16), 𝑝𝑝 the signals 𝐼𝐼𝐼𝐼 must be arranged in a matrix of size 𝑖𝑖 , 𝑖𝑖 =1…𝑁𝑁 𝑒𝑒 ( ) 𝑛𝑛 𝑁𝑁 ×𝑁𝑁 , as represented in Fig. 14 of the appendix (see also s e 𝑒𝑒 the Matlab function ezsos in the appendix 6.C). 𝐅𝐅 is large but very sparse since it mostly contains zero- DAS valued elements (Fig. 7). Let 𝑛𝑛𝑛𝑛 𝑧𝑧 denote its number of nonzero elements. If a 𝑞𝑞 -point interpolation is used to estimate the sig- nals at 𝜏𝜏 (𝑿𝑿 ), we have 𝑛𝑛𝑛𝑛 𝑧𝑧 ≤ (𝑀𝑀 𝑁𝑁 𝑞𝑞 ). The right-hand term is 𝑖𝑖 𝑠𝑠 𝑒𝑒 an upper bound since limited apertures (i.e. 𝑓𝑓 > 0) can be used Fig. 7. Example of a sparse DAS matrix. DAS beamforming can be written during receive beamforming, which reduces 𝑛𝑛𝑛𝑛 𝑧𝑧 . The sparsity as a sparse matrix-vector multiplication. For example, the generation of a 256×256 image by beamforming 128 signals each containing 1000 samples [33] of the DAS matrix verifies: leads to a DAS matrix of size (256×256) × (128×1000) = 65536 × 128000. 𝑛𝑛𝑛𝑛 𝑧𝑧 𝑞𝑞 Note that the matrix shown is for illustration purpose only (its aspect can be 𝑠𝑠𝑠𝑠𝑠𝑠𝑝𝑝𝑠𝑠 𝑠𝑠 (𝐅𝐅 ) = 1− ≥ 1− . (19) DAS quite different). 𝑀𝑀 = number of beamformed data (is generally chosen equal 𝑀𝑀 𝑛𝑛 𝑁𝑁 𝑛𝑛 s e 𝑠𝑠 to the number of pixels of the image to be formed); 𝑛𝑛 = number of time 𝑠𝑠 samples per raw signals; 𝑁𝑁 = number of array elements. To give an example, if each array element acquires 1000 I/Q 𝑒𝑒 samples and if the signals 𝐼𝐼𝐼𝐼 are interpolated linearly (i.e. 𝑞𝑞 = 𝑖𝑖 2), then the DAS matrix has a sparsity greater than 1-2/1000 = F. DAS as a matrix product 99.8%, i.e. a density smaller than 0.2%. Sparse matrix-vector As can be seen by Eq. (15), the DAS is a linear operator that multiplication (SpMV) can be computed on GPUs [34]. Its transforms the I/Q signals recorded by the array elements (𝐼𝐼𝐼𝐼 𝑖𝑖 computational complexity is 𝑂𝑂 (𝑛𝑛𝑛𝑛 𝑧𝑧 ) ≤ 𝑂𝑂 (𝑀𝑀 𝑁𝑁 𝑞𝑞 ) [33]. A con- 𝑒𝑒 with 𝑖𝑖 = 1 …𝑁𝑁 ) to beamformed I/Q data (𝐼𝐼𝐼𝐼 ). If the time se- 𝑒𝑒 bf struction of the complex DAS matrix (for I/Q beamforming) is ries 𝐼𝐼𝐼𝐼 contain 𝑛𝑛 samples each, and are stacked in a column 𝑖𝑖 s given in the Matlab function (called ezdas) included in the ap- ( ) vector = �𝐼𝐼 𝐼𝐼 , … ,𝐼𝐼 𝐼𝐼 � of size 𝑛𝑛 𝑁𝑁 × 1 , DAS beam- pendix. The proposed function ezdas uses a linear interpolation, 1 𝑁𝑁 s 𝑒𝑒 which is recommended in most situations. For the sake of com- forming can be written synthetically as: pleteness, the dasmtx function of the MUST toolbox = 𝐅𝐅 × , (18) 𝑏𝑏 𝑓𝑓 DAS (www.biomecardio.com/MUST/) includes several interpola- tion methods: nearest neighbor, linear, quadratic, 5-point least- where 𝐅𝐅 is the DAS matrix (Fig. 7) that contains the in- DAS squares parabolic, and 3-or-5-lobe Lanczos. terpolation weights (see the end of paragraph 2.B), the phase rotators (see equation 10), and the sum truncation induced by a Fig. 7 illustrates the sparsity of such a matrix. If the ultrasound positive f-number (see equation 15). If the data (e.g. B-mode or sequence (array, transmit, receive) and the 𝑀𝑀 beamforming Doppler) to be constructed contain 𝑀𝑀 pixels, then the column point locations remain unchanged, 𝐅𝐅 only needs to be cal- DAS vector that contains the beamformed I/Q signals is of 𝑏𝑏 𝑓𝑓 culated once. Once the DAS matrix is created (or loaded), length 𝑀𝑀 . The matrix 𝐅𝐅 is of size (𝑀𝑀 ×𝑛𝑛 𝑁𝑁 ). It is complex DAS s e SpMV-based beamforming is very fast and compatible with due to the phase rotations present in Eq. (15). The matrix-vector real-time visualization [32]. product for DAS beamforming is schematized in Fig. 7. It is understood that RF signals can also be processed by this matrix product, in which case the DAS matrix is real. Nevertheless, sufficiently fine axial sampling is required to allow subsequent 𝑰𝑰𝑰𝑰 𝑰𝑰𝑰𝑰 𝑰𝑰𝑰𝑰 𝑰𝑰𝑰𝑰 𝑖𝑖𝑡𝑡 𝑰𝑰𝑰𝑰 7 Perrot et al. 3. In vitro and in vivo examples were transmitted at a pulse repetition frequency >6 kHz to cre- ate one coherently compound image. The parameter 𝒬𝒬 was cal- 𝑝𝑝 The respective effects of the speed of sound and 𝑓𝑓 -number on culated by using the compound hyperbolas whose vertices be- image quality were evaluated using the PICMUS dataset avail- longed to a 128×256 beamforming grid. able online [35]. We analyzed RF data that were acquired in an in vitro CIRS phantom (040GSE) and sampled at four times the center frequency. These data were obtained by transmitting steered plane waves with a 128-element 5-MHz linear array. The RF data were IQ-demodulated then beamformed using the DAS matrix of Eq. (18) on the Cartesian image grid specified in PICMUS. A series of real envelopes were constructed from one unsteered plane wave (no transmission delay) for a large range of speed of sound (1400 to 1700 m/s) or 𝑓𝑓 -number (0 to 4). For reference, high-quality real envelopes were also gener- ated by coherent compounding with the whole dataset (75 planes waves with steering ranging from -16 to +16 degrees). The optimal 𝑓𝑓 -number was estimated by using Eq. (12) and (13) with 𝐷𝐷 = 0.71 (-3dB threshold). The optimal speed of 𝑡𝑡 ℎ𝑟𝑟 ℎ Fig. 9. In vitro results obtained in a CIRS phantom using single plane wave sound was determined by maximizing the phase-based metric imaging: effect of the speed of sound and 𝑓𝑓 -number on FWHM and CNR. using Eq. (16). Contrast and lateral resolution were quantified The errors in % are expressed relative to the optimum values (shown by the using the contrast-to-noise ratio (CNR) and full width at half dot). The horizontal and vertical lines represent the optimal speed of sound and 𝑓𝑓 -number estimated by the equations introduced in the manuscript. maximum (FWHM), respectively, as detailed in [35]. The CNR and FWHM reported in our work (Fig. 9) were an average from two cysts and seven nylon wires, respectively (see Figures 2b A. In vitro results and 2d in [35]). From Eq. (14), the minimum significant wavelength trans- mitted by the elements was 𝜆𝜆 ≈ 0.23 mm. It followed that min 𝑊𝑊 /𝜆𝜆 ≈ 1.17. From Eq. (12) and (13), with 𝐷𝐷 = 0.71 min 𝑡𝑡 ℎ𝑟𝑟 ℎ (i.e. -3dB directivity threshold), the directivity-based 𝑓𝑓 -number was thus 𝑓𝑓 ≈ 1.4. Maximizing 𝒬𝒬 (16) by emitting one plane # 𝑝𝑝 wave returned an optimal speed of sound 𝑐𝑐 ̂ = 1570 m/s. Fig. 8. In vitro results obtained in a CIRS phantom (experimental data from PICMUS). The top and bottom rows show the effect of the speed of sound and 𝑓𝑓 -number on CNR (contrast-to-noise ratio) and lateral resolution, re- spectively. The thick vertical bars represent the optimal speed of sound and 𝑓𝑓 -number estimated by the equations introduced in the manuscript. The gray dots represent the metrics after compounding 75 plane waves (reference val- ues), while the black dots are for a single plane wave transmit. Fig. 10. In vivo results obtained in 8 carotids. The curves represent the 𝒬𝒬 𝑝𝑝 metric as a function of the speed of sound. The boxplot shows the distribu- The speed of sound was also measured in the carotids of eight th th tion of the estimated speeds of sound (median, range, 25 , 75 percentiles). healthy volunteers by maximizing the 𝒬𝒬 metric. A cardiologist 𝑝𝑝 used a linear-array transducer (ATL L7-4, center frequency = 5.2 MHz, element width = 0.245 mm, fractional bandwidth As anticipated, both the speed of sound and the 𝑓𝑓 -number at -6dB = 65%) connected to a Verasonics scanner (V-1-128, had a significant impact on the contrast-to-noise ratio (CNR) Verasonics Inc.) to scan the common carotid artery (protocol and lateral resolution (Fig. 8 and Fig. 9). In the CIRS phantom, similar to that described in [36]). Five steered plane waves CNR and lateral resolution were both optimal for 𝑐𝑐 ≈ 1550- (transmit beam angles evenly spaced between -10° and 10°) 𝑒𝑒𝑠𝑠 𝑒𝑒𝑠𝑠 8 Perrot et al. 1600 m/s and 𝑓𝑓 ≈ 1.1-1.6. These values determined using a respective effects are well known, how to set these parameters brute-force search were consistent with those estimated by had not been clearly explained in the literature. We have ad- physical reasoning (phase homogeneity and element directiv- dressed this issue by proposing two simple approaches that can be summarized by equations (12)-(13) and (16)-(17). As we ity) using Eq. (12), (13) and (16), (17); see Fig. 8 and Fig. 9. have argued, for a relatively homogeneous medium, the 𝑓𝑓 -num- ber is essentially related to the directivity of the array elements. B. In vivo results It depends on the transducer (element width, center frequency) The calculated 𝑓𝑓 -number was also 𝑓𝑓 ≈ 1.4 (same transducer and pulse-echo bandwidth. A modification is needed when an- as in vitro). The speeds of sound estimated through maximizing gled receiving is performed, as in vector Doppler [36], [37]. In 𝒬𝒬 in the eight carotids were 1510 ± 41 m/s (Fig. 10). Fig. 11 𝑝𝑝 this case, the angle variable 𝜗𝜗 in Eq. (13) must be replaced by shows a carotid whose speed of sound was estimated at 1400 ( | |), with 𝜃𝜃 being the receive steering angle. 𝜗𝜗 + 𝜃𝜃 𝑇𝑇𝑅𝑅 𝑇𝑇𝑅𝑅 m/s. This value may seem abnormally low. However, this speed An appropriately chosen 𝑓𝑓 -number helps to optimize the bias- of sound probably does not reflect the average speed of sound variance tradeoff when summing the signals along the hyperbo- in the carotid artery. This value has somewhat improved the im- las. An 𝑓𝑓 -number that is too large (too small an aperture) in- age by maximizing the metric we introduced. This apparent duces a summation based on a small sample of elements, thus speed of sound probably reduced the deterioration of the image generating a significant bias. In the case of an 𝑓𝑓 -number that is caused by a combination of independent effects (fat, air, aber- too small (too large an aperture), the summation involves low rations...). The 𝑓𝑓 -number had an obvious positive impact on im- st nd directivities (and therefore low SNRs), thus causing a high var- age quality (Fig. 11, 1 row vs. 2 row). In comparison, the iance. Too high bias or variance affects the contrast and lateral effect of the speed of sound was less noticeable, except for the st nd resolution (Fig. 8). We heuristically found that a -3dB-directiv- change of scale in the 𝑧𝑧 -direction (Fig. 11, 1 column vs. 2 ity threshold gives a good compromise (i.e. Eq. (13) with column). 𝐷𝐷 = 0.71). Similar conclusions could be reached with 𝑡𝑡 ℎ𝑟𝑟 ℎ Doppler or vector Doppler. Fig. 12 illustrates one vector Dop- pler example recomputed from [36] (plane wave imaging, re- ceive angle 𝜃𝜃 = ±15°, 3 cm-diameter disk rotating at 300 𝑇𝑇𝑅𝑅 rpm). The smallest velocity-vector errors coincided with the di- rectivity-derived 𝑓𝑓 -number (≈ 2.6) obtained from Eq. (13) with (𝜗𝜗 + |𝜃𝜃 |) instead of 𝜗𝜗 . 𝑇𝑇𝑅𝑅 Fig. 11. In vivo example in one carotid. The directivity-derived 𝑓𝑓 -number was 1.4 and the estimated speed of sound was 1400 m/s (instead of the as- Fig. 12. Vector Doppler obtained in a rotating disk (data from [36]): 17 un- sumed 1540 m/s). 𝑓𝑓 -number = 0 means “full aperture”. The artefacts related steered plane waves were transmitted; I/Q signals were DAS-beamformed to a full aperture are well visible in the lumen (top row). with ±15° receive angles. The curves show the median relative errors on the 𝑥𝑥 - and 𝑧𝑧 -velocity components as a function of the 𝑓𝑓 -number. The thick ver- tical bars represent the optimal f-number estimated by Eq. (12)-(13). 4. Discussion In this essentially educational article, we have considered the fundamentals of the delay-and-sum (DAS) for ultrasound image beamforming. The underlying theory has been described in depth so that users can be aware of the functioning and limita- tions of this technique. In particular, we have explained that the DAS can be written numerically as a sparse matrix-vector prod- Fig. 13. Knowing the speed of sound (𝑐𝑐 ) allows one to determine which uct and we have clarified the choice of the 𝑓𝑓 -number and speed hyperbola must be chosen. The 𝑓𝑓 -number allows one to decide which part of sound from a physical and practical viewpoint. of the hyperbola should be considered. A. The 𝑓𝑓 -number and speed of sound in DAS To put it simply, the 𝑓𝑓 -number allows one to decide which We have seen that the 𝑓𝑓 -number and speed of sound both in- part of the hyperbola should be considered. Knowing the speed fluence image quality returned by DAS (Fig. 8). Although their 𝑒𝑒𝑠𝑠 9 Perrot et al. of sound, on the other hand, allows one to determine which hy- different stages, for example at the beamforming [48] or com- perbola should be chosen (Fig. 13). According to the in vitro pounding [49] level. If properly trained, DL systems are ex- results obtained in a calibration phantom, the speed of sound pected to provide high-quality images with fewer input data. has a significant impact on the standard CNR and FWHM met- Data-driven DL algorithms, however, are highly dependent on rics (Fig. 8). However, when considering the in vivo results, the accurate and clean training datasets to learn from: for DL sys- overall physiological appearance appears almost unchanged tems to learn to produce high-quality images, they must be (Fig. 11), except for the axial dimensions. This would probably heavily trained with high-quality images. The question then have little or no diagnostic impact in a clinical setting. It then arises as to how to provide such images of optimal quality for a seems that adjusting the speed of sound would be impactful wide variety of organs. If DAS is the chosen beamformer to mainly for calibration or in vitro comparisons. Although this is train DL systems, it is essential to ensure that its parameters are only a speculation, there might nevertheless be some interest in correctly set. super-resolution localization of microbubbles, as each mi- crobubble is expected to be located at the peak or centroid of 5. Conclusion the bubble-PSF (point spread function) [38]. Optimizing the DAS is highly popular in ultrasound imaging. If properly de- speed of sound and thus refining the PSF could make localiza- signed and compound, it can be very effective in terms of in tion more robust. On a more clinical level, it has been shown vivo image quality, even in challenging situations such as the that the longitudinal speed of sound has a potential diagnostic beating heart [50]. Yet DAS is sometimes considered mediocre value in hepatic steatosis, a common liver disease [39], [40]. In by investigations based mainly on in vitro experiments because particular, Imbault et al. [39] estimated the speed of sound by it is often suboptimally addressed. This opinion is often biased locally analyzing the spatial coherence [41]. The measurement by its misuse. For example, the 𝑓𝑓 -number is sometimes ne- of the proposed dimensionless metric 𝒬𝒬 (16) based on I/Q sig- glected or wrongly chosen. It is imperative that the receive 𝑓𝑓 - 𝑝𝑝 nals could provide an alternative method. number be properly selected. In addition, the choice of the speed of sound can have a significant impact on the metrics B. DAS vs. advanced beamformers measured in vitro. To get the most out of DAS and to prevent As stated earlier, DAS is the best-known and simplest beam- any misguided discredit, it is suggested to optimize its parame- forming technique in medical ultrasound imaging. It is fast, ters in order to use it wisely. When it comes to DAS, one must easy to program and well efficient in most situations. Moreover, not DAS out of step. DAS allows for low-complexity real-time beamforming [42]. Yet, several studies have focused on adaptive beamformers 6. Appendix whose intended purpose is to improve image quality for medical A. Generalized transmit distance (focused and circular beams) ultrasound. One typical approach for adaptive beamforming is If (𝑥𝑥 ,𝑧𝑧 ) is a focus point, a reasoning similar to that used to derive to perform a weighted sum and adjust the weights to the data. 0 0 Eq. (5) leads to: This is similar to an ordinary DAS, except that instead of each of the signal samples (inside the receive aperture) contributing 2 2 𝑑𝑑 (𝑿𝑿 ) = sign(𝑧𝑧 −𝑧𝑧 )� (𝑥𝑥 −𝑥𝑥 ) + (𝑧𝑧 −𝑧𝑧 ) TX 𝑠𝑠 𝑠𝑠 0 𝑠𝑠 0 𝑠𝑠 0 (20) equally to the final sum, some data samples contribute more 2 2 + � H(|𝑥𝑥 | +𝐿𝐿 /2 )(|𝑥𝑥 | +𝐿𝐿 /2) +𝑧𝑧 . 0 0 0 than others do, and the relative contributions depend on these Equations (5) and (20) can be combined to yield a generalized transmit samples. Minimum variance (MV) and its variant (eigen-based distance (focusing or diverging wavefronts): MV) as well as 𝑝𝑝 -DAS are ones of these techniques [43]–[45]. (21) Adaptive methods based on the coherence factor [46] also be- 2 2 𝑑𝑑 (𝑿𝑿 ) = sign(𝑧𝑧 −𝑧𝑧 )� (𝑥𝑥 −𝑥𝑥 ) + (𝑧𝑧 −𝑧𝑧 ) TX 𝑠𝑠 𝑠𝑠 0 𝑠𝑠 0 𝑠𝑠 0 long in this category, though the weights are pixelwise. Because 2 2 ( )� (| | ( ) )(| | ( ) ) + sign 𝑧𝑧 H 𝑥𝑥 + sign 𝑧𝑧 𝐿𝐿 /2 𝑥𝑥 + sign 𝑧𝑧 𝐿𝐿 /2 +𝑧𝑧 . they are data-dependent, some adaptive beamforming tech- 0 0 0 0 0 0 niques are computationally expensive and cannot be used in B. Plane wave as a limit of diverging wave real-time. More importantly, although they can improve com- mon metrics (e.g. CNR, FWHM) in specific anechoic-cyst- or Taking the limit of the virtual source coordinates, given by Eq. (4), as 𝛽𝛽 → 0 yields wire-based calibration phantoms, it is not clear if they have sig- nificant value in a clinical context, as to whether they can help ( ) 𝐿𝐿 sin(2𝜃𝜃 ) cos 𝜃𝜃 sin(𝜃𝜃 ) lim 𝑥𝑥 = = 𝐿𝐿 , in a better diagnosis at the patient’s bedside. Another downside ⎪ + 𝛽𝛽 →0 2 𝛽𝛽 𝛽𝛽 (22) of this adaptability is that local phase information, from one 𝐿𝐿 1 + cos(2𝜃𝜃 ) cos (𝜃𝜃 ) lim 𝑧𝑧 = − = −𝐿𝐿 , emission to the next, can be distorted. This may have a detri- ⎪ 0 𝛽𝛽 →0 2 𝛽𝛽 𝛽𝛽 mental impact on velocity estimation by color Doppler: the ( ) ( ) ( ) ( ) Let 𝑓𝑓 𝜃𝜃 = 𝐿𝐿 cos 𝜃𝜃 sin(𝜃𝜃 ) and 𝑔𝑔 𝜃𝜃 = 𝐿𝐿 cos 𝜃𝜃 . The transmit dis- more data-dependent the beamformer, the higher the Doppler tance (5) as 𝛽𝛽 → 0 is thus reduced to error [47]. To avoid biased conclusions when comparing alter- native beamformers against DAS through these metrics, it is lim 𝑑𝑑 (𝑥𝑥 ,𝑧𝑧 ) TX 𝑠𝑠 𝑠𝑠 𝛽𝛽 →0 therefore essential to ensure that the speed of sound and/or the 2 2 (23) = � (𝑥𝑥 −𝑓𝑓 (𝜃𝜃 )⁄𝛽𝛽 ) + (𝑧𝑧 +𝑔𝑔 (𝜃𝜃 )⁄𝛽𝛽 ) 𝑠𝑠 𝑠𝑠 𝑓𝑓 -number have been correctly set for DAS beamforming. 2 2 ( ( ) ( )⁄ ) ( ( )⁄ ) −�sgn 𝜃𝜃 𝑓𝑓 𝜃𝜃 𝛽𝛽 −𝐿𝐿 /2 + 𝑔𝑔 𝜃𝜃 𝛽𝛽 . Another increasingly growing direction for ultrasound image formation is deep learning. Deep learning (DL) can intervene at Expanding the four square terms then factoring yields: 10 Perrot et al. % See also DAS, DASMTX. lim 𝑑𝑑 (𝑥𝑥 ,𝑧𝑧 ) TX 𝑠𝑠 𝑠𝑠 𝛽𝛽 →0 2 2 % -- Damien Garcia -- 2019/11 �𝑓𝑓 (𝜃𝜃 ) +𝑔𝑔 (𝜃𝜃 ) % www.biomecardio.com 𝛽𝛽 2 siz0 = size(x); 2(𝑧𝑧 𝑔𝑔 (𝜃𝜃 )−𝑥𝑥 𝑓𝑓 (𝜃𝜃 )) 𝑥𝑥 +𝑧𝑧 ² 𝑠𝑠 𝑠𝑠 𝑠𝑠 𝑠𝑠 2 (24) � [nl,nc] = size(SIG); �1 + 𝛽𝛽 + 𝛽𝛽 2 2 2 2 𝑓𝑓 (𝜃𝜃 ) +𝑔𝑔 (𝜃𝜃 ) 𝑓𝑓 (𝜃𝜃 ) +𝑔𝑔 (𝜃𝜃 ) x = x(:); z = z(:); % ULA (uniform linear array): sgn(𝜃𝜃 )𝑓𝑓 (𝜃𝜃 )𝐿𝐿 𝐿𝐿 − 1− 𝛽𝛽 + 𝛽𝛽 � . % x-coordinates of the elements 2 2 2 2 ( ) ( ) ( ) ( ) 𝑓𝑓 𝜃𝜃 +𝑔𝑔 𝜃𝜃 4(𝑓𝑓 𝜃𝜃 +𝑔𝑔 𝜃𝜃 ) xe = ((0:nc-1)-(nc-1)/2)*param.pitch; L = xe(end)-xe(1); % length of the array 1/2 ( ) at the first order yields: Using the Taylor series of 1 +𝑥𝑥 % coordinates of the virtual source 𝑓𝑓 (𝜃𝜃 )(sgn(𝜃𝜃 )𝐿𝐿 − 2𝑥𝑥 ) + 2𝑔𝑔 (𝜃𝜃 )𝑧𝑧 x0 = vsource(1); z0 = vsource(2); 𝑠𝑠 𝑠𝑠 lim 𝑑𝑑 (𝑥𝑥 ,𝑧𝑧 ) = (25) TX 𝑠𝑠 𝑠𝑠 𝛽𝛽 →0 2 2 2�𝑓𝑓 (𝜃𝜃 ) +𝑔𝑔 (𝜃𝜃 ) % transmit & receive distances dTX = hypot(x-x0,z-z0)-... −𝜋𝜋 𝜋𝜋 Noting that cos(𝜃𝜃 ) is always positive (since 𝜃𝜃 ∈ �, � ), Eq. (25) can hypot((abs(x0)-L/2)*(abs(x0)>L/2),z0); 2 2 be simplified to obtain Eq. (6) after replacing 𝑓𝑓 (𝜃𝜃 ) and 𝑔𝑔 (𝜃𝜃 ). dRX = hypot(x-xe,z); % two-way travel times C. Determine the f-number in Matlab tau = (dTX+dRX)/param.c; % Note: in Matlab, sinc(x) = sin(pi*x)/(pi*x) f = @(th,width,lambda)... % fast-time indices abs(cos(th)*sinc(width/lambda*sin(th))-0.71); idxt = tau*param.fs + 1; alpha = fminbnd(@(th) f(th,width,lambda),0,pi/2); fnumber = 1/2/tan(alpha); % boolean vectors I = idxt>=1 & idxt<=nl-1; Iaperture = abs(x-xe)<=(z/2/param.fnumber); I = I&Iaperture; D. A simple Matlab code for DAS beamforming % linear indices function [bfSIG,M] = ezdas(SIG,x,z,vsource,param) idx = idxt + (0:nc-1)*nl; idx = idx(I); %EZDAS Delay-And-Sum beamforming idxf = floor(idx); % (Easy version of DAS) idx = idxf-idx; % BFSIG = EZDAS(SIG,X,Z,VSOURCE,PARAM) beamforms % DAS matrix % the RF or I/Q signals stored in the array SIG, [i,~] = find(I); % and returns the beamformed signals BFSIG. The s = [idx+1;-idx]; % (for linear interpolation) % signals are beamformed at the points specified if ~isreal(SIG) % if IQ: phase rotations % by X and Z. s = s.*exp(2i*pi*param.fc*[tau(I);tau(I)]); end % 1) SIG must be a 2-D array. The first dimension M = sparse([i;i],[idxf;idxf+1],s,numel(x),nl*nc); % (i.e. each column) corresponds to single RF % or I/Q signals over (fast-) time, with the % DAS beamforming % 1st column corresponding to the 1st element. bfSIG = reshape(M*SIG(:),siz0); % 2) VSOURCE contains the coordinates [x0,z0] of % the virtual point source. Use large x0,z0 for % plane waves. % 3) PARAM is a structure that contains the % parameter values required for the delay-and- % sum (see below for details). % Note: SIG must be complex for I/Q data % (i.e. SIG = complex(I,Q) = I + 1i*Q). % [~,M] = EZDAS(...) also returns the DAS matrix. % PARAM must contain the following fields: % --------------------------------------- % 1) PARAM.fs: sampling frequency (Hz) % 2) PARAM.pitch: element pitch (m) % 3) PARAM.fc: center frequency (Hz) % 4) PARAM.c: longitudinal velocity (m/s) % 5) PARAM.fnumber: receive f-number Fig. 14. DAS matrix product to obtain interpolated I/Q along the M hyper- bolas ℋ associated to the 𝑀𝑀 pixels of the beamformed image. 𝑛𝑛 = number 𝑠𝑠 𝑠𝑠 % --- of time samples per raw signals; 𝑁𝑁 = number of array elements. 𝑒𝑒 % NOTE #1: Interpolation method: EZDAS uses a % linear interpolation to generate the DAS matrix. % --- E. A simple Matlab code for estimating the speed of sound % NOTE #2: EZDAS is for pedagogical purpose. Use % DAS of the MUST toolbox for more options. function c = ezsos(IQ,x,z,vsource,param) % --- % 11 Perrot et al. %EZSOS Speed-of-sound estimation if wasrow, RF = RF(:); end % (Easy version of SOS) % %-- Time vector % c = EZSOS(IQ,X,Z,VSOURCE,PARAM) returns the nl = size(RF,1); % speed of sound that yields an "optimal" real- t = (0:nl-1)'/Fs; % envelope image by DAS beamforming. The optimal % speed of sound is estimated by analyzing the %-- Downmixing of the RF signals % phases along the diffraction hyperbolas whose IQ = double(RF).*exp(-1i*2*pi*Fc*t); % vertices are located at (X,Z). % %-- Low-pass filter % The input arguments are the same as those of [b,a] = butter(5,.5); % EZDAS. Except that the signals MUST be complex IQ = filtfilt(b,a,IQ)*2; % (I/Q data). Enter "help EZDAS" for details. % %-- Recover the initial size (if was a row vector) % See also SOS, EZDAS. if wasrow, IQ = IQ.'; end % -- Damien Garcia & Vincent Perrot -- 2019/11 % www.biomecardio.com References assert(~isreal(IQ)) [nl,nc] = size(IQ); c0 = param.c; [1] R. E. McKeighen and M. P. Buchin, “New techniques for dynamically variable electronic delays for real time ultrasonic imaging,” in 1977 IQ0 = IQ; Ultrasonics Symposium, 1977, pp. 250–254, doi: IQ = sparse(1:nl*nc,kron(1:nc,ones(1,nl)),... 10.1109/ULTSYM.1977.196834. IQ(:),nl*nc,nc); [2] R. J. Mailloux, “Phased array theory and technology,” Proceedings of the IEEE, vol. 70, no. 3, pp. 246–291, 1982, doi: c = round(fminbnd(@PBC,1200,1700)); 10.1109/PROC.1982.12285. [3] H. T. Friis and C. B. Feldman, “A multiple unit steerable antenna for function Qp = PBC(c) short-wave reception,” Proceedings of the Institute of Radio Engi- param.c = c; neers, vol. 25, no. 7, pp. 841–917, 1937, doi: [~,M] = ezdas(IQ0,x,z/c0*c,vsource,param); 10.1109/JRPROC.1937.228354. [4] P. J. Kahrilas, “HAPDAR—An operational phased array radar,” Pro- hyperb = full(M*IQ); ceedings of the IEEE, vol. 56, no. 11, pp. 1967–1975, 1968, doi: % (contains the diffraction hyperbolas) 10.1109/PROC.1968.6773. [5] P. E. Green, R. A. Frosch, and C. F. Romney, “Principles of an experi- A = angle(hyperb); mental large aperture seismic array (LASA),” Proceedings of the A(A==0) = NaN; IEEE, vol. 53, no. 12, pp. 1821–1833, 1965, doi: A = unwrap(A,[],2); 10.1109/PROC.1965.4453. Qp = abs(mean(hyperb,2,'omitnan')).^2./... [6] Ö. Yilmaz, O. Yilmaz, and S. M. Doherty, “Migration,” in Seismic std(A,0,2,'omitnan').^2; data analysis: processing, inversion, and interpretation of seismic Qp = -mean(Qp,'omitnan'); data, SEG Books, 2001, pp. 463–653. [7] S. H. Gray, J. Etgen, J. Dellinger, and D. Whitmore, “Seismic migra- end tion problems and solutions,” Geophysics, vol. 66, no. 5, pp. 1622– 1640, 2001, doi: 10.1190/1.1487107. end [8] S. D. Pye, S. R. Wild, and W. N. McDicken, “Adaptive time gain compensation for ultrasonic imaging,” Ultrasound in Medicine & Bi- F. A simple Matlab code for I/Q demodulation ology, vol. 18, no. 2, pp. 205–212, 1992, doi: 10.1016/0301- function IQ = ezrf2iq(RF,Fs,Fc) 5629(92)90131-S. [9] A. Kyriakou, E. Neufeld, B. Werner, M. M. Paulides, G. Szekely, and %EZRF2IQ I/Q demodulation of RF data N. Kuster, “A review of numerical and experimental compensation % (easy version of RF2IQ) techniques for skull-induced phase aberrations in transcranial focused ultrasound,” International Journal of Hyperthermia, vol. 30, no. 1, pp. % IQ = RF2IQ(RF,Fs,Fc) demodulates the radio- 36–46, 2014, doi: 10.3109/02656736.2013.861519. % frequency (RF) bandpass signals and returns the [10] S. J. Sanabria, E. Ozkan, M. Rominger, and O. Goksel, “Spatial do- % Inphase/Quadrature (I/Q) components. IQ is a main reconstruction for imaging speed-of-sound with pulse-echo ul- % complex whose real (imaginary) part contains the trasound: simulation and in vivo study,” Phys. Med. Biol., vol. 63, no. % inphase (quadrature) component. 21, p. 215015, 2018, doi: 10.1088/1361-6560/aae2fb. [11] E. Boni, A. C. H. Yu, S. Freear, J. A. Jensen, and P. Tortoli, “Ultra- % 1) Fs is the sampling frequency (in Hz), sound open platforms for next-generation imaging technique develop- % 2) Fc represents the center frequency (in Hz). ment,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Fre- quency Control, vol. 65, no. 7, pp. 1078–1092, 2018, doi: % EZRF2IQ demodulates along columns for 2-D and 10.1109/TUFFC.2018.2844560. % 3-D RF data. Each column corresponds to a single [12] J. Faurie, M. Baudet, J. Poree, G. Cloutier, F. Tournoux, and D. Gar- % RF signal over (fast-) time. cia, “Coupling myocardium and vortex dynamics in diverging-wave echocardiography,” IEEE Trans Ultrason Ferroelectr Freq Control, % EZRF2IQ is for pedagogical purpose. You may use vol. 66, no. 3, pp. 425–432, 2019, doi: % RF2IQ for more options. 10.1109/TUFFC.2018.2842427. [13] O. M. H. Rindal, A. Austeng, A. Fatemi, and A. Rodriguez-Molares, % See also RF2IQ. “The effect of dynamic range alterations in the estimation of contrast,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency % -- Damien Garcia -- 2019/11 Control, vol. 66, no. 7, pp. 1198–1208, 2019, doi: % www.biomecardio.com 10.1109/TUFFC.2019.2911267. [14] F. Destrempes and G. Cloutier, “A critical review and uniformized %-- Convert to column vector (if RF is a row vector) representation of statistical distributions modeling the ultrasound echo wasrow = isrow(RF); 12 Perrot et al. envelope,” Ultrasound in Medicine & Biology, vol. 36, no. 7, pp. [34] M. Maggioni and T. Berger-Wolf, “Optimization techniques for 1037–1051, 2010, doi: 10.1016/j.ultrasmedbio.2010.04.001. sparse matrix–vector multiplication on GPUs,” Journal of Parallel and Distributed Computing, vol. 93–94, pp. 66–86, 2016, doi: [15] S. Salles, H. Liebgott, O. Basset, C. Cachard, D. Vray, and R. Lavarello, “Experimental evaluation of spectral-based quantitative ul- 10.1016/j.jpdc.2016.03.011. trasound imaging using plane wave compounding,” IEEE Transac- [35] H. Liebgott, A. Rodriguez-Molares, F. Cervenansky, J. A. Jensen, and tions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 61, O. Bernard, “Plane-wave imaging challenge in medical ultrasound,” no. 11, pp. 1824–1834, Nov. 2014, doi: IEEE International Ultrasonics Symposium (IUS), pp. 1–4, 2016, doi: 10.1109/ULTSYM.2016.7728908. 10.1109/TUFFC.2014.006543. [16] Jian-Yu Lu, “Experimental study of high frame rate imaging with lim- [36] C. Madiena, J. Faurie, J. Porée, and D. Garcia, “Color and vector flow ited diffraction beams,” IEEE Transactions on Ultrasonics, Ferroelec- imaging in parallel ultrasound with sub-Nyquist sampling,” IEEE trics, and Frequency Control, vol. 45, no. 1, pp. 84–97, 1998, doi: Trans Ultrason Ferroelectr Freq Control, vol. 65, no. 5, pp. 795–802, 10.1109/58.646914. 2018. [37] B. Y. S. Yiu and A. C. H. Yu, “Least-squares multi-angle Doppler es- [17] D. Garcia, L. Le Tarnec, S. Muth, E. Montagnon, J. Porée, and G. Cloutier, “Stolt’s f-k migration for plane wave ultrasound imaging,” timators for plane-wave vector flow imaging,” IEEE Transactions on IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Ultrasonics, Ferroelectrics, and Frequency Control, vol. 63, no. 11, Control, vol. 60, no. 9, pp. 1853–1867, 2013, doi: pp. 1733–1744, 2016, doi: 10.1109/TUFFC.2016.2582514. 10.1109/TUFFC.2013.2771. [38] J. Yu, L. Lavery, and K. Kim, “Super-resolution ultrasound imaging method for microvasculature in vivo with a high temporal accuracy,” [18] D. Garcia, L. Kadem, D. Savéry, P. Pibarot, and L.-G. Durand, “Ana- lytical modeling of the instantaneous maximal transvalvular pressure Scientific Reports, vol. 8, no. 1, pp. 1–11, 2018, doi: 10.1038/s41598- gradient in aortic stenosis,” Journal of Biomechanics, vol. 39, no. 16, 018-32235-2. pp. 3036–3044, 2006, doi: 10.1016/j.jbiomech.2005.10.013. [39] M. Imbault et al., “Robust sound speed estimation for ultrasound- [19] G. Montaldo, M. Tanter, J. Bercoff, N. Benech, and M. Fink, “Coher- based hepatic steatosis assessment,” Phys Med Biol, vol. 62, no. 9, pp. ent plane-wave compounding for very high frame rate ultrasonogra- 3582–3598, 2017, doi: 10.1088/1361-6560/aa6226. phy and transient elastography,” IEEE Transactions on Ultrasonics, [40] R. E. Zubajlo et al., “Experimental validation of longitudinal speed of Ferroelectrics, and Frequency Control, vol. 56, no. 3, pp. 489–506, sound estimates in the diagnosis of hepatic steatosis (part II),” Ultra- 2009, doi: 10.1109/TUFFC.2009.1067. sound in Medicine & Biology, vol. 44, no. 12, pp. 2749–2758, 2018, [20] J. Porée, D. Posada, A. Hodzic, F. Tournoux, G. Cloutier, and D. Gar- doi: 10.1016/j.ultrasmedbio.2018.07.020. cia, “High-frame-rate echocardiography using coherent compounding [41] R. Mallart and M. Fink, “Adaptive focusing in scattering media with Doppler-based motion-compensation,” IEEE Transactions on through sound‐speed inhomogeneities: The van Cittert Zernike ap- Medical Imaging, vol. 35, no. 7, pp. 1647–1657, 2016, doi: proach and focusing criterion,” The Journal of the Acoustical Society 10.1109/TMI.2016.2523346. of America, vol. 96, no. 6, pp. 3721–3732, 1994, doi: [21] S. Bae, P. Kim, and T. Song, “Ultrasonic sector imaging using plane 10.1121/1.410562. wave synthetic focusing with a convex array transducer,” The Journal [42] B. Y. S. Yiu, I. K. H. Tsang, and A. C. H. Yu, “GPU-based beam- of the Acoustical Society of America, vol. 144, no. 5, pp. 2627–2644, former: fast realization of plane wave compounding and synthetic ap- 2018, doi: 10.1121/1.5065391. erture imaging,” IEEE Transactions on Ultrasonics, Ferroelectrics, [22] S. Shahriari and D. Garcia, “Meshfree simulations of ultrasound vec- and Frequency Control, vol. 58, no. 8, pp. 1698–1705, 2011, doi: tor flow imaging using smoothed particle hydrodynamics,” Phys Med 10.1109/TUFFC.2011.1999. Biol, vol. 63, pp. 1–12, 2018, doi: 10.1088/1361-6560/aae3c3. [43] J. F. Synnevag, A. Austeng, and S. Holm, “Adaptive beamforming ap- [23] C. P. Loizou and C. S. Pattichis, “An overview of despeckle-filtering plied to medical ultrasound imaging,” IEEE Transactions on Ultra- techniques,” in Handbook of speckle filtering and tracking in cardio- sonics, Ferroelectrics, and Frequency Control, vol. 54, no. 8, pp. vascular ultrasound imaging and video, IET Digital Library, 2018, pp. 1606–1613, 2007, doi: 10.1109/TUFFC.2007.431. 95–109. [44] B. M. Asl and A. Mahloojifar, “Eigenspace-based minimum variance [24] W. Burger and M. J. Burge, “Interpolation,” in Principles of digital beamforming applied to medical ultrasound imaging,” IEEE Transac- image processing: core algorithms, London: Springer-Verlag, 2009, tions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 57, pp. 210–237. no. 11, pp. 2381–2390, 2010, doi: 10.1109/TUFFC.2010.1706. [25] T. L. Szabo, Diagnostic ultrasound imaging: inside out. Academic [45] M. Polichetti, F. Varray, J.-C. Béra, C. Cachard, and B. Nicolas, “A Press, 2004. nonlinear beamformer based on p-th root compression—application to [26] D. C. M. Horvat, J. S. Bird, and M. M. Goulding, “True time-delay plane wave ultrasound imaging,” Applied Sciences, vol. 8, no. 4, p. bandpass beamforming,” IEEE Journal of Oceanic Engineering, vol. 599, 2018, doi: 10.3390/app8040599. 17, no. 2, pp. 185–192, 1992, doi: 10.1109/48.126975. [46] J. Camacho, M. Parrilla, and C. Fritsch, “Phase coherence imaging,” [27] A. R. Selfridge, G. S. Kino, and B. T. Khuri‐Yakub, “A theory for the IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency radiation pattern of a narrow‐strip acoustic transducer,” Appl. Phys. Control, vol. 56, no. 5, pp. 958–974, 2009, doi: Lett., vol. 37, no. 1, pp. 35–36, 1980, doi: 10.1063/1.91692. 10.1109/TUFFC.2009.1128. [28] D. Napolitano et al., “Sound speed correction in ultrasound imaging,” [47] M. Polichetti, V. Perrot, H. Liebgott, B. Nicolas, and F. Varray, “In- Ultrasonics, vol. 44, pp. e43–e46, 2006, doi: 10.1016/j.ul- fluence of beamforming methods on velocity estimation: in vitro ex- tras.2006.06.061. periments,” in 2018 IEEE International Ultrasonics Symposium (IUS), [29] E. L. Madsen, J. A. Zagzebski, R. A. Banjavie, and R. E. Jutila, “Tis- 2018, pp. 1–4, doi: 10.1109/ULTSYM.2018.8580186. sue mimicking materials for ultrasound phantoms,” Medical Physics, [48] A. C. Luchies and B. C. Byram, “Deep neural networks for ultrasound vol. 5, no. 5, pp. 391–394, 1978, doi: 10.1118/1.594483. beamforming,” IEEE Transactions on Medical Imaging, vol. 37, no. [30] C. Yoon, Y. Lee, J. H. Chang, T. Song, and Y. Yoo, “In vitro estima- 9, pp. 2010–2021, 2018, doi: 10.1109/TMI.2018.2809641. tion of mean sound speed based on minimum average phase variance [49] M. Gasse, F. Millioz, E. Roux, D. Garcia, H. Liebgott, and D. Fribou- in medical ultrasound imaging,” Ultrasonics, vol. 51, no. 7, pp. 795– let, “High-quality plane wave compounding using convolutional neu- 802, 2011, doi: 10.1016/j.ultras.2011.03.007. ral networks,” IEEE Transactions on Ultrasonics, Ferroelectrics, and [31] F. Vignon et al., “Fast frame rate 2D cardiac deformation imaging Frequency Control, vol. 64, no. 10, pp. 1637–1639, 2017, doi: based on RF data: what do we gain?,” in 2017 IEEE International Ul- 10.1109/TUFFC.2017.2736890. trasonics Symposium (IUS), 2017, pp. 1–4, doi: [50] J. Porée, M. Baudet, F. Tournoux, G. Cloutier, and D. Garcia, “A dual 10.1109/ULTSYM.2017.8092648. tissue-Doppler optical-flow method for speckle tracking echocardiog- [32] G. Y. Hou et al., “Sparse matrix beamforming and image reconstruc- raphy at high frame rate,” IEEE Trans Med Imaging, vol. 37, no. 9, tion for 2-D HIFU monitoring using harmonic motion imaging for fo- pp. 2022–2032, 2018, doi: 10.1109/TMI.2018.2811483. cused ultrasound (HMIFU) with in vitro validation,” IEEE Transac- tions on Medical Imaging, vol. 33, no. 11, pp. 2107–2117, 2014, doi: 10.1109/TMI.2014.2332184. [33] J. R. Gilbert, C. Moler, and R. Schreiber, “Sparse matrices in Matlab: design and implementation,” SIAM J. Matrix Anal. Appl., vol. 13, no. 1, pp. 333–356, 1992, doi: 10.1137/0613024.

Journal

Electrical Engineering and Systems SciencearXiv (Cornell University)

Published: Jul 23, 2020

There are no references for this article.