SIMUS: an open-source simulator for ultrasound imaging. Part I: theory & examples
[[ This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] Abstract
Index Terms —Ultrasonic transducer arrays, Computer simula-tion, Ultrasound imaging, Open-source codes. I. I NTRODUCTION
OMPUTATIONAL ultrasound imaging, which uses nu-merical analysis to solve problems that involve ultrasound wave propagations, has become a standard methodology in the medical and non-destructive ultrasound community. Before considering in vitro or in vivo investigations, computational ul-trasound imaging can be used, for example, to 1) analyze ultra-sound sequences and arrays [1], 2) develop or optimize beam-forming or post-processing algorithms [2], 3) explore multiple configurations through serial tests [3], 4) compare with peers in international challenges [4]. Among the freely available ultra-sound simulators, Field II [5], [6], and k-Wave [7], [8] are likely D. G. is with INSERM at CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image pour la Santé), CNRS UMR 5220 – INSERM U1206 – Université Lyon 1 – INSA Lyon - Université Jean Monnet Saint-Etienne (e-mail: [email protected]; [email protected]). D. G. began this work when he was head of the RUBIC (Research Unit of Biomechanics and Imaging in Cardiology) at the LBUM, CRCHUM (Centre de the most popular. These Matlab toolboxes have widely pro-moted the use of ultrasound simulations for research purposes, and the number of works that use these tools has been increas-ing over the years (Fig. 1). There is a whole range of software packages dedicated to ultrasound imaging, available for free, as open-source or not. A non-exhaustive list of ultrasound-imag-ing programs is available on the k-Wave website recherche du centre hospitalier de l'Université de Montréal), Montreal, QC, Canada. SIMUS: an open-source simulator for ultrasound imaging.
Part I: theory & examples
Damien Garcia C This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] At the time of submission of this paper, only 1-D probes, recti-linear or convex, with elevation focusing, are considered in SIMUS. Although there is a growing interest in ultrasound im-aging with a high number of transducer elements (e.g. 1024) and 2-D matrix arrays [24], [25], it appears that the 1-D config-uration with a limited number of channels (typically 64 to 192) remains by far the most common configuration at present. The main assumptions on which SIMUS relies are 1.
Linearity, 2.
Scatterers acting as monopole sources, 3.
Weak (single) scattering. In essence, these hypotheses are similar to those in Field II. SIMUS, however, works in the Fourier domain. The Fourier domain is indeed more appropriate in many ways. Besides some numerical advantages, the Fourier domain can easily consider the physical aspects that depend on the frequency, such as the directivity of the elements and the attenuation. The challenge was that the function SIMUS could be easily manipulated after a few minutes of training. For the sake of clarity, the syntax of SIMUS has been standardized with that of the functions in-cluded in the MUST toolbox, and the default settings are those commonly used in medical ultrasound. It should be noted that SIMUS calls the PFIELD function, which is the main program. PFIELD calculates acoustic pressure fields and derives the pressure signals reflected by monopole point scatterers and measured by the elements. This article (
Part I: theory ) is completed by a second one (
Part II: comparison with FieldII, k-Wave, and Verasonics , Varray F. and Garcia D.) This first part describes the linear acoustic theory that is in PFIELD. Several approximations and linearizations have been used. It is essential to review them to identify the limits of PFIELD and under which conditions it can be used. I have illustrated Part I with simulations of sound pres-sure fields and ultrasound images. The second article (Part II) is devoted to the comparison of the acoustic pressures generated by PFIELD with those obtained by Field II, k-Wave, and the Verasonics simulator.
Fig. 1. Number of yearly publications that cite [5] (Field II) or [7] (k-Wave). The citation reports are from Web of Science. In this first part, I will outline the theoretical reasoning lead-ing to the equations included in PFIELD. I will first explain how sound pressure fields are simulated and then address the gener-ation of backscattered pressure signals. The last section of Part I will be illustrated with some realistic examples related to med-ical imaging. II.
PFIELD
INSIDE OUT
This section describes the theory inside PFIELD. PFIELD simulates the pressure fields in the frequency domain; i.e. it is assumed that the pressure waves have a harmonic time depend-ence such that the pressure is written as: 𝑝𝑝 ( 𝑿𝑿 , 𝑡𝑡 ) = Re �� 𝑃𝑃 ( 𝑿𝑿 , 𝜔𝜔 , 𝑡𝑡 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 d 𝜔𝜔 +∞−∞ � , (1) with 𝑿𝑿 representing a point in the radiated region of interest, and 𝜔𝜔 being the angular frequency. In the following subsec-tions, I will describe how the pressure component 𝑃𝑃 ( 𝑿𝑿 , 𝜔𝜔 , 𝑡𝑡 ) generated from one array element can be approximated, from which the backscattered echoes will be deduced. Although PFIELD also works for curvilinear arrays, the following sec-tions describe the theory in the context of a rectilinear probe (Fig. 2). The interested reader can refer to the PFIELD code to learn about the slight modifications required for a convex probe. To estimate the waves that are backscattered by a me-dium of point-like scatterers, one must first calculate the acous-tic pressure radiated by a single element transducer and an ul-trasound array. As we will see, the width ( = 2 𝑏𝑏 , see Fig. 2) and height ( ℎ ) of the element transducers are two key parameters. Recalling that only 1-D arrays will be addressed, elevation fo-cusing will be taken into consideration. Fig. 2. Coordinate system for a rectilinear array. In this paper, the height of an element is noted ℎ , and its width is 𝑏𝑏 . 𝑅𝑅 𝑓𝑓 stands for the distance to the elevation focus. This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] A. Overview of the problem
We will use the conventional coordinate system for a recti-linear array (Fig. 2), i.e. 𝑋𝑋 along the azimuthal direction, 𝑌𝑌 for the elevation, and 𝑍𝑍 describing the axial position. A similar sys-tem, noted in lowercase letters ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ), will be applied for in-dividual elements or sub-elements (Fig. 3): i.e. 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 all equal zero at the center of an individual (sub-)element. The model starts with the Rayleigh-Sommerfeld integral, which describes an isolated element behaving like a baffled piston that vibrates on the 𝑥𝑥 - 𝑦𝑦 plane (Fig. 4). From the pressure waves radiated by a single element will follow those of the ultrasound probe. The distance between a point 𝒙𝒙′ = ( 𝑥𝑥 ′ , 𝑦𝑦 ′ , 0) of the piston and a point 𝒙𝒙 = ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) in the field (Fig. 3), is noted 𝑟𝑟 ′ = � ( 𝑥𝑥 − 𝑥𝑥 ′ ) + ( 𝑦𝑦 − 𝑦𝑦 ′ ) + 𝑧𝑧 . (2) The distance between the center of the piston and a point 𝒙𝒙 in the field (Fig. 3), is noted 𝑟𝑟 = �𝑥𝑥 + 𝑦𝑦 + 𝑧𝑧 . (3) Fig. 3. Coordinate system for an individual element. The distance 𝑟𝑟′ is ap-proximated by the expression (6).
Fig. 4. An individual element acts as a baffled piston vibrating in the 𝑧𝑧 -direction. The velocities are delayed to simulate the elevation focusing in-duced by an acoustic lens. Medical ultrasound one-dimensional arrays contain an acoustic lens that focuses the ultrasound waves on the elevation plane. In PFIELD, the incident waves can be focused on the elevation plane at a given distance 𝑅𝑅 𝑓𝑓 (Fig. 2 and Fig. 4), whose value is generally provided by the probe manufacturer. A strat-egy for simulating elevation focusing is to use a large concave element [26] or to position small elements onto a curved surface [27]. Another strategy is to modify the piston velocity delays along the elevational direction (Fig. 4), as described by Eq. (3.27) in [22]: 𝑣𝑣 ( 𝑦𝑦 ′ , 𝜔𝜔 ) = 𝑣𝑣 ( 𝜔𝜔 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖 𝑦𝑦′22𝑅𝑅𝑓𝑓 if | 𝑦𝑦 ′ | ≤ ℎ2 , otherwise, (4) with 𝑣𝑣 ( 𝜔𝜔 ) being the velocity amplitude. B. Acoustic field of a single array element
Let us consider a planar piston embedded in an infinite rigid baffle, and vibrating along the z -direction. The resulting har-monic pressure 𝑃𝑃 at position 𝒙𝒙 = ( 𝑥𝑥 , 𝑦𝑦 , 𝑧𝑧 ) is given by the Ray-leigh-Sommerfeld integral (see e.g. Eq. 1 in [21] or Eq. 6.19 in [22]) 𝑃𝑃 ( 𝒙𝒙 , 𝜔𝜔 , 𝑡𝑡 ) = 𝑖𝑖𝑘𝑘𝑘𝑘2𝑖𝑖𝑖𝑖 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 � � 𝑣𝑣 ( 𝜔𝜔 ) 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 ′ 𝑟𝑟 ′ 𝑑𝑑𝑥𝑥 ′𝑑𝑑𝑦𝑦′ ℎ /2 − ℎ /2 𝑏𝑏−𝑏𝑏 , (5) where 𝑖𝑖 = √− , 𝑡𝑡 is time, 𝜔𝜔 is the angular frequency, 𝜌𝜌 is the medium density, 𝑐𝑐 is the speed of sound, and 𝑘𝑘 = 𝜔𝜔 / 𝑐𝑐 is the wavenumber. Assuming a paraxial propagation with respect to the y -direction, the distance 𝑟𝑟′ [Eq. (2)] can be rewritten, in the far field, as (see Appendix) 𝑟𝑟 ′ ≈ 𝑟𝑟 − 𝑥𝑥 ′ sin 𝜃𝜃 + �𝑦𝑦−𝑦𝑦 ′ � , (6) where the angle 𝜃𝜃 is defined in Fig. 3. The far-field pressure can thus be approximated by inserting (6) into (5) and by sub-stituting 𝑟𝑟 for 𝑟𝑟 ’ in the denominator: 𝑃𝑃 ( 𝒙𝒙 , 𝜔𝜔 , 𝑡𝑡 )= 𝑖𝑖𝑘𝑘𝑘𝑘2𝑖𝑖𝑖𝑖 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 � 𝑣𝑣 ( 𝑦𝑦 ′ , 𝜔𝜔 ) 𝑒𝑒 𝑖𝑖𝑖𝑖�𝑟𝑟−𝑥𝑥 ′ sin 𝜃𝜃+ �𝑦𝑦−𝑦𝑦′�22𝑟𝑟 � 𝑟𝑟 𝑏𝑏 , ℎ / , −ℎ / 𝑑𝑑𝑥𝑥 ′ 𝑑𝑑𝑦𝑦 ′ , (7) which gives 𝑃𝑃 ( 𝒙𝒙 , 𝜔𝜔 , 𝑡𝑡 ) = 𝑖𝑖𝑘𝑘𝑘𝑘2𝑖𝑖𝑖𝑖 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑟𝑟 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 × … �� 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑥𝑥 ′ sin 𝜃𝜃 𝑑𝑑𝑥𝑥 ′𝑏𝑏−𝑏𝑏 � �� 𝑣𝑣 ( 𝑦𝑦 ′ , 𝜔𝜔 ) 𝑒𝑒 𝑖𝑖𝑖𝑖 �𝑦𝑦−𝑦𝑦′�22𝑟𝑟 ℎ / / 𝑑𝑑𝑦𝑦 ′ � . (8) Including the expression of the piston velocity [Eq. (4)] in the integrand of the second integral of Eq. (8) yields 𝑃𝑃 ( 𝒙𝒙 , 𝜔𝜔 , 𝑡𝑡 ) = 𝑖𝑖𝑘𝑘𝑘𝑘2𝑖𝑖𝑖𝑖 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑟𝑟 𝑣𝑣 ( 𝜔𝜔 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 × … (9) This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] �� 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑥𝑥 ′ sin 𝜃𝜃 𝑑𝑑𝑥𝑥 ′𝑏𝑏−𝑏𝑏 � �� 𝑒𝑒 −𝑖𝑖𝑘𝑘 𝑦𝑦′ 𝑅𝑅𝑓𝑓 𝑒𝑒 𝑖𝑖𝑘𝑘 � 𝑦𝑦−𝑦𝑦′ � 𝑟𝑟 𝑑𝑑𝑦𝑦 ′ℎ /2 −ℎ /2 � . The two integrals in the curly brackets will now be deter-mined. The first integral in Eq. (9) yields � 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑥𝑥 ′ sin 𝜃𝜃 𝑑𝑑𝑥𝑥 ′𝑏𝑏−𝑏𝑏 = 2 𝑏𝑏 sinc( 𝑘𝑘𝑏𝑏 sin 𝜃𝜃 ). (10) The second integral of Eq. (9) could be explicitly expressed by using the imaginary error function (erfi). However, the nu-merical estimation of the erfi function, e.g. through estimating the Faddeeva function [28], is computationally expensive. I thus opted for the use of a Gaussian superposition model [20]. The second integral of Eq. (9) is rewritten as � 𝑒𝑒 −𝑖𝑖𝑖𝑖 𝑦𝑦′22𝑅𝑅𝑓𝑓 𝑒𝑒 𝑖𝑖𝑖𝑖 �𝑦𝑦−𝑦𝑦′�22𝑟𝑟 𝑑𝑑𝑦𝑦 ′ℎ / / = � Π � 𝑦𝑦 ′ ℎ � 𝑒𝑒 −𝑖𝑖𝑖𝑖 𝑦𝑦′22𝑅𝑅𝑓𝑓 𝑒𝑒 𝑖𝑖𝑖𝑖 �𝑦𝑦−𝑦𝑦′�22𝑟𝑟 𝑑𝑑𝑦𝑦 ′+∞−∞ . (11) where Π stands for the rectangle function. In the Gaussian superposition model, the rectangle function is approximated by a sum of Gaussians with complex coefficients (Fig. 5): Π � 𝑦𝑦 ′ ℎ � ≈ � 𝐴𝐴 𝑔𝑔 𝑒𝑒 −𝐵𝐵 𝑔𝑔 � 𝑦𝑦′ℎ � 𝐺𝐺𝑔𝑔=1 . (12) Equation (11) thus becomes � Π � 𝑦𝑦 ′ ℎ � 𝑒𝑒 −𝑖𝑖𝑖𝑖 𝑦𝑦′22𝑅𝑅𝑓𝑓 𝑒𝑒 𝑖𝑖𝑖𝑖 �𝑦𝑦−𝑦𝑦′�22𝑟𝑟 𝑑𝑑𝑦𝑦 ′+∞−∞ ≈ � 𝐺𝐺𝑔𝑔=1 𝐴𝐴 𝑔𝑔 � 𝑒𝑒 −𝛼𝛼 𝑔𝑔 𝑦𝑦 ′2 +𝛽𝛽𝑦𝑦 ′ +𝛾𝛾 𝑑𝑑𝑦𝑦 ′+∞−∞ , (13) where 𝛼𝛼 𝑔𝑔 = 𝐵𝐵 𝑔𝑔 ℎ +𝑖𝑖𝑖𝑖2 � 1𝑅𝑅 𝑓𝑓 − ; 𝛽𝛽 = −𝑖𝑖𝑖𝑖𝑦𝑦𝑟𝑟 ; 𝛾𝛾 = 𝑖𝑖𝑖𝑖𝑦𝑦 . (14) Solving the right-hand side Gaussian integral in Eq. (13) yields � Π � 𝑦𝑦 ′ ℎ � 𝑒𝑒 −𝑖𝑖𝑖𝑖 𝑦𝑦′22𝑅𝑅𝑓𝑓 𝑒𝑒 𝑖𝑖𝑖𝑖 �𝑦𝑦−𝑦𝑦′�22𝑟𝑟 𝑑𝑑𝑦𝑦 ′+∞−∞ ≈ � 𝐺𝐺𝑔𝑔=1 𝐴𝐴 𝑔𝑔 � 𝑖𝑖𝛼𝛼 𝑔𝑔 𝑒𝑒 𝛽𝛽24𝛼𝛼𝑔𝑔+𝛾𝛾 . (15) From Eq. (15), the second integral of Eq. (9) thus reduces to � 𝑒𝑒 −𝑖𝑖𝑖𝑖 𝑦𝑦′22𝑅𝑅𝑓𝑓 𝑒𝑒 𝑖𝑖𝑖𝑖 �𝑦𝑦−𝑦𝑦′�22𝑟𝑟 𝑑𝑑𝑦𝑦 ′ℎ / / ≈ � 𝐺𝐺𝑔𝑔 =1 𝐴𝐴 𝑔𝑔 � 𝜋𝜋𝛼𝛼 𝑔𝑔 𝑒𝑒 𝛽𝛽 𝛼𝛼𝑔𝑔 + 𝛾𝛾 . (16) Replacing the two integrals in (9) by their respective expres-sions (10) and (16) provides an estimate of the acoustic pressure generated by a single element: 𝑃𝑃 ( 𝒙𝒙 , 𝜔𝜔 , 𝑡𝑡 ) ≈� 𝑖𝑖𝑏𝑏𝑖𝑖𝑖𝑖 𝜌𝜌𝑐𝑐𝑣𝑣 ( 𝜔𝜔 ) ���������� 𝑃𝑃 Tx ( 𝑖𝑖 ) 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑟𝑟 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 sinc( 𝑘𝑘𝑏𝑏 sin 𝜃𝜃 ) ��������� 𝐷𝐷 ( 𝜃𝜃 , 𝑖𝑖 ) �∑ 𝐴𝐴 𝑔𝑔 � 𝜋𝜋𝛼𝛼𝑔𝑔 𝑒𝑒 𝛽𝛽24𝛼𝛼𝑔𝑔+𝛾𝛾 𝐺𝐺𝑔𝑔=1 �������������� δ ( 𝑦𝑦 , 𝑟𝑟 , 𝑖𝑖 ) . (17) The coefficients 𝐴𝐴 𝑔𝑔 and 𝐵𝐵 𝑔𝑔 can be determined through an op-timization method. We found that four coefficients (Fig. 5) of-fer a good compromise when comparing against Field II and k-Wave (see Part II ): 𝐴𝐴 = 0.187 + 0.275 𝑖𝑖 , 𝐵𝐵 = 4.558 − 𝑖𝑖 , 𝐴𝐴 = 0.288 − 𝑖𝑖 , 𝐵𝐵 = 8.598 − 𝑖𝑖 , 𝐴𝐴 = 𝐴𝐴 ��� , 𝐵𝐵 = 𝐵𝐵 ��� , 𝐴𝐴 = 𝐴𝐴 ��� , 𝐵𝐵 = 𝐵𝐵 ��� . (18) Lists with up to 25 coefficients are provided in [19] and [28]. Fig. 5. The rectangle function can be approximated by the sum of Gaussi-ans, which simplifies the estimation of the integral in Eq. (11).
In practice, it is not the velocity of the element that is known, but the acoustic pressure generated by an element, measured for example with a hydrophone. The first term in brackets (dimen-sionally homogeneous to pressure) in the expression (17) rep-resents the spectrum of the transmit pressure pulse (up to a con-stant multiplier) , noted 𝑃𝑃 𝑇𝑇𝑥𝑥 ( 𝜔𝜔 ) . The sine cardinal sinc term rep-resents the 𝑥𝑥 -directivity of one element, noted 𝐷𝐷 ( 𝜃𝜃 , 𝑘𝑘 ) . The last term in brackets is related to the elevation focusing. It is homo-geneous to a distance and is noted δ ( 𝑦𝑦 , 𝑟𝑟 , 𝑘𝑘 ) . Using these nota-tions, the acoustic pressure of one element finally reduces to 𝑃𝑃 ( 𝒙𝒙 , 𝜔𝜔 , 𝑡𝑡 ) ≈ 𝑃𝑃 Tx ( 𝜔𝜔 ) 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑟𝑟 𝐷𝐷 ( 𝜃𝜃 , 𝑘𝑘 ) δ ( 𝑦𝑦 , 𝑟𝑟 , 𝑘𝑘 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 . (19) https://mathworld.wolfram.com/GaussianIntegral.html This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] More generally, if the transmission is delayed by
Δ𝜏𝜏 , the wave field is given by 𝑃𝑃 ( 𝒙𝒙 , 𝜔𝜔 , 𝑡𝑡 ) ≈ 𝑃𝑃 Tx ( 𝜔𝜔 ) 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑟𝑟 𝐷𝐷 ( 𝜃𝜃 , 𝑘𝑘 ) δ ( 𝑦𝑦 , 𝑟𝑟 , 𝑘𝑘 ) 𝑒𝑒 𝑖𝑖𝑖𝑖Δ𝜏𝜏 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 . (20) It is recalled that the position variables ( 𝒙𝒙 , 𝜃𝜃 , 𝑟𝑟 ) in (20) are relative to the element (Fig. 3). C. Acoustic field of a rectilinear array
The expression (20) models acoustic waves radiated by one element. To derive this expression, the distance 𝑟𝑟 with respect to the element [see Eq. (6)] was simplified by using a Fraunho-fer (far-field) approximation in the azimuthal 𝑥𝑥 -direction, and a Fresnel (paraxial) approximation in the elevation 𝑦𝑦 -direction. In order to fulfill the far-field condition, it may be necessary to split the array elements into 𝜈𝜈 sub-elements, in the azimuthal 𝑥𝑥 -direction, if they are too wide (Fig. 6). In PFIELD, 𝜈𝜈 is chosen so that a sub-element width ( = 2 𝑏𝑏 / 𝜈𝜈 ) is not greater than the minimal wavelength (defined at -6dB): 𝜈𝜈 = � min � , (21) where ⌈ ⌉ being the ceiling function. As an indication, the number 𝜈𝜈 of simulated sub-element(s) per array element is typ-ically 1 for cardiac phased arrays, and 2 for linear arrays. If an array contains 𝑁𝑁 elements, the total number of sub-elements is thus ( 𝜈𝜈𝑁𝑁 ) . Given the properties of linearity, the acoustic wave-field produced by an 𝑁𝑁 -element array can be modeled by super-imposing the ( 𝜈𝜈𝑁𝑁 ) individual sub-element models described by (20): 𝑃𝑃 ( 𝑿𝑿 , 𝜔𝜔 , 𝑡𝑡 ) ≈𝑃𝑃 Tx ( 𝜔𝜔 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 � 𝜈𝜈𝜈𝜈𝑛𝑛=1 𝑊𝑊 𝑛𝑛 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑛𝑛 𝑟𝑟 𝑛𝑛 𝐷𝐷 ( 𝜃𝜃 𝑛𝑛 , 𝑘𝑘 ) δ ( 𝑌𝑌 , 𝑟𝑟 𝑛𝑛 , 𝑘𝑘 ) 𝑒𝑒 𝑖𝑖𝑖𝑖Δ𝜏𝜏 𝑛𝑛 . (22) Fig. 6. To ensure that the far-field condition is met, the 𝑁𝑁 array elements are each divided into 𝜈𝜈 sub-elements. In (22), the position variables ( 𝜃𝜃 𝑛𝑛 , 𝑦𝑦 𝑛𝑛 , 𝑟𝑟 𝑛𝑛 ) are relative to the 𝑛𝑛 th sub-element. One has 𝑦𝑦 𝑛𝑛 = 𝑌𝑌 for a 1-D array. The position 𝑿𝑿 = ( 𝑋𝑋 , 𝑌𝑌 , 𝑍𝑍 ) is relative to the coordinate system of the array depicted in Fig. 2. If 𝑋𝑋 𝑘𝑘 , 𝑛𝑛 stands for the abscissa of the 𝑛𝑛 th sub-element centroid, then 𝑟𝑟 𝑛𝑛 = � ( 𝑋𝑋 − 𝑋𝑋 𝑘𝑘 , 𝑛𝑛 ) + 𝑌𝑌 + 𝑍𝑍 , and sin 𝜃𝜃 𝑛𝑛 = ( 𝑋𝑋 − 𝑋𝑋 𝑘𝑘 , 𝑛𝑛 )/ � ( 𝑋𝑋 − 𝑋𝑋 𝑘𝑘 , 𝑛𝑛 ) + 𝑍𝑍 . (23) We here used 𝑋𝑋 𝑘𝑘 , 𝑛𝑛 = 0 (non-convex array). Let the leftmost to rightmost sub-elements be numbered sequentially, from to 𝑛𝑛 . In the case of a uniform linear array of pitch 𝑝𝑝 (Fig. 6), it can be shown that its centroid abscissa can be written as 𝑋𝑋 𝑘𝑘 , 𝑛𝑛 = 𝑝𝑝2 � � 𝑛𝑛𝜈𝜈 � − 𝑁𝑁 − � + 𝑏𝑏𝜈𝜈 (2( 𝑛𝑛 − 𝜈𝜈 ) − 𝜈𝜈 + 1) (24) The term � 𝑛𝑛𝜈𝜈 � corresponds to the number of the element to which the sub-element 𝑛𝑛 belongs. Note that only 𝑋𝑋 𝑘𝑘 , 𝑛𝑛 (and 𝑍𝑍 𝑘𝑘 , 𝑛𝑛 ≠ ) must be modified for a convex array (see the PFIELD code for details). The transmit delay Δ𝜏𝜏 𝑛𝑛 = Δ𝜏𝜏 ⌈𝑛𝑛 / 𝜈𝜈⌉ in (22) is that of the ⌈𝑛𝑛 / 𝜈𝜈⌉ th element. Each element has been weighted by 𝑊𝑊 𝑛𝑛 = 𝑊𝑊 ⌈𝑛𝑛 / 𝜈𝜈⌉ to consider transmit apodization. The equation (22) is the backbone of PFIELD. Once the transmit pressure 𝑃𝑃 Tx ( 𝜔𝜔 ) is given, it allows simulation of realistic fields of acous-tic pressure produced by an ultrasound array. In PFIELD, the transmit pressure is generated by convolving the one-way re-sponse of the transducer with a windowed sine, as explained in the next section. From the acoustic reciprocity principle, Eq. (22) can also be applied to derive the backscattered echoes, as it will be explained in section IV. Fig. 7.
Left – Pulse-echo response of the transducer. 𝜔𝜔 𝑘𝑘 is the center angular frequency ( 𝜔𝜔 𝑘𝑘 = 2 𝜋𝜋𝑓𝑓 𝑘𝑘 ). 𝜔𝜔 𝑏𝑏 is the angular frequency bandwidth. Right – Example of a one-way transmit pulse (center frequency = 7.6 MHz, frac-tional bandwidth = 77%, excitation pulse = 1 wavelength).
III. S PECTRUM OF THE TRANSMIT PRESSURE
The general expression (22) contains the spectrum of the transmit pressure 𝑃𝑃 Tx ( 𝜔𝜔 ) . In SIMUS and PFIELD, the transmit pressure waveform is obtained by convolving a rectangularly-windowed sinusoid (a “perfect” pulse) with the point spread function (PSF) of the transducer. The angular frequency of the sinusoid is that of the transducer ( 𝜔𝜔 𝑘𝑘 = 2 𝜋𝜋𝑓𝑓 𝑘𝑘 , the center fre-quency 𝑓𝑓 𝑘𝑘 being provided by the manufacturer). The temporal width 𝑇𝑇 of the rectangular window is defined in terms of num-ber of wavelengths 𝑛𝑛 𝜆𝜆 by 𝑇𝑇 = 𝑛𝑛 𝜆𝜆 / 𝑓𝑓 𝑘𝑘 = 2 𝜋𝜋𝑛𝑛 𝜆𝜆 / 𝜔𝜔 𝑘𝑘 . The spectrum This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] of the rectangularly-windowed pulse is given by 𝒮𝒮 P ( 𝜔𝜔 ) = 𝑖𝑖 � sinc �𝑇𝑇 𝑖𝑖−𝑖𝑖 𝑐𝑐 � − sinc �𝑇𝑇 𝑖𝑖+𝑖𝑖 𝑐𝑐 �� (25) The spectrum of the transducer PSF is defined by a general-ized Gaussian window that depends on two positive parameters 𝑝𝑝 and 𝜎𝜎 : 𝒮𝒮 T ( 𝜔𝜔 ) = 𝑒𝑒 −� | 𝑖𝑖−𝑖𝑖 𝑐𝑐 | 𝜎𝜎 𝑖𝑖 𝑐𝑐 �𝑝𝑝 (26) It is designed so that its pulse-echo response has a given bandwidth 𝜔𝜔 𝑏𝑏 at -6 dB (Fig. 7). The pulse-echo fractional band-width is generally given by the manufacturer in percent. For ex-ample, a 65% bandwidth means that the frequency bandwidth of the response is such that 𝜔𝜔 𝑏𝑏 = 0.65 𝜔𝜔 𝑘𝑘 . To determine both 𝑝𝑝 and 𝜎𝜎 , it is postulated that 𝒮𝒮 T (0) = 2 −126 (the smallest positive single-precision floating-point number). It follows that (see Ap-pendix) 𝒮𝒮 T ( 𝜔𝜔 ) = 𝑒𝑒 − ln 2� | 𝜔𝜔−𝜔𝜔𝑐𝑐 | 𝜔𝜔𝑏𝑏 �𝑝𝑝 , with 𝑝𝑝 = ln 126 / ln � 𝑖𝑖 𝑐𝑐 𝑖𝑖 𝑏𝑏 � . (27) Fig. 7 (left panel) depicts an example of the simulated trans-ducer response for a pulse-echo fractional bandwidth of 77%. PFIELD does not include the electrical-acoustic conversions that occur in the piezoelectric elements. The units of the simu-lated acoustic and electrical signals are thus arbitrary (not Pa or V). By using this simplified representation and writing the con-volution in the frequency domain, the spectrum of the transmit pressure 𝑃𝑃 Tx ( 𝜔𝜔 ) is given by this proportionality relationship: 𝑃𝑃 Tx ( 𝜔𝜔 ) ∝ 𝒮𝒮 P ( 𝜔𝜔 ) �𝒮𝒮 T ( 𝜔𝜔 ) . (28) A square root is needed since the transducer response 𝒮𝒮 T ( 𝜔𝜔 ) is two-way (transmit + receive). Fig. 7 (right panel) presents a transmit pulse (one-way) in the temporal domain. IV. SIMUS
INSIDE OUT A. Echoes received by a sub-element
SIMUS uses PFIELD and point-like scatterers to simulate ra-diofrequency (RF) ultrasound signals. These scatterers become individual monopole point sources when an incident wave reaches them. They do not acoustically interact with each other according to the assumption of single weak scattering. Each scatterer is defined by its reflection coefficient ( ℛ 𝓈𝓈 ), which de-scribes how much amplitude of a wave is reflected. Although some tissues, such as blood, are governed by Rayleigh scatter-ing [30], the ℛ 𝓈𝓈 coefficients are assumed constant, i.e. inde-pendent of frequency and incidence angle. From (22), the pres-sure signal received by a scatterer 𝓈𝓈 located at 𝑿𝑿 𝓼𝓼 =( 𝑋𝑋 𝓈𝓈 , 𝑌𝑌 𝓈𝓈 , 𝑍𝑍 𝓈𝓈 ) is 𝑃𝑃 ( 𝑿𝑿 𝓼𝓼 , 𝜔𝜔 , 𝑡𝑡 ) ≈ 𝑃𝑃 Tx ( 𝜔𝜔 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 � 𝜈𝜈𝜈𝜈𝑛𝑛=1 𝑊𝑊 𝑛𝑛 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑛𝑛𝓈𝓈 𝑟𝑟 𝑛𝑛𝓈𝓈 𝐷𝐷 ( 𝜃𝜃 𝑛𝑛𝓈𝓈 , 𝑘𝑘 ) δ ( 𝑌𝑌 𝓈𝓈 , 𝑟𝑟 𝑛𝑛𝓈𝓈 , 𝑘𝑘 ) 𝑒𝑒 𝑖𝑖𝑖𝑖Δ𝜏𝜏 𝑛𝑛 , (29) where 𝑟𝑟 𝑛𝑛𝓈𝓈 = � ( 𝑋𝑋 𝓈𝓈 − 𝑋𝑋 𝑘𝑘 , 𝑛𝑛 ) + 𝑌𝑌 𝓈𝓈 + 𝑍𝑍 𝓈𝓈 and sin 𝜃𝜃 𝑛𝑛𝓈𝓈 =( 𝑋𝑋 𝓈𝓈 − 𝑋𝑋 𝑘𝑘 , 𝑛𝑛 )/ � ( 𝑋𝑋 𝓈𝓈 − 𝑋𝑋 𝑘𝑘 , 𝑛𝑛 ) + 𝑍𝑍 𝓈𝓈 . The principle of acoustic reciprocity dictates that an acoustic response remains identical when the source and receiver are exchanged. Expression (19) can therefore give the pressure received by the 𝑚𝑚 th sub-element from a scatterer 𝓈𝓈 , after accounting for its reflection coefficient ℛ 𝓈𝓈 : 𝑃𝑃 se 𝑚𝑚𝓈𝓈 ( 𝜔𝜔 , 𝑡𝑡 ) ≈ ℛ 𝓈𝓈 𝑃𝑃 ( 𝑿𝑿 𝓼𝓼 , 𝜔𝜔 , 𝑡𝑡 ) ������� 𝑃𝑃 ( 𝑿𝑿 𝓼𝓼 , 𝑖𝑖 ) 𝑒𝑒 −𝑖𝑖𝜔𝜔𝑖𝑖 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑚𝑚𝓈𝓈 𝑟𝑟 𝑚𝑚𝓈𝓈 𝐷𝐷 ( 𝜃𝜃 𝑚𝑚𝓈𝓈 , 𝑘𝑘 ) δ ( 𝑌𝑌 𝓈𝓈 , 𝑟𝑟 𝑚𝑚𝓈𝓈 , 𝑘𝑘 ) , (30) where 𝑟𝑟 𝑚𝑚𝓈𝓈 = � ( 𝑋𝑋 𝓈𝓈 − 𝑋𝑋 𝑘𝑘 , 𝑚𝑚 ) + 𝑌𝑌 𝓈𝓈 + 𝑍𝑍 𝓈𝓈 and sin 𝜃𝜃 𝑚𝑚 𝓈𝓈 =( 𝑋𝑋 𝓈𝓈 − 𝑋𝑋 𝑘𝑘 , 𝑚𝑚 )/ � ( 𝑋𝑋 𝓈𝓈 − 𝑋𝑋 𝑘𝑘 , 𝑚𝑚 ) + 𝑍𝑍 𝓈𝓈 . In (30), the superscript “ se ” means “sub-element”. Inserting (29) in (30), it follows that: 𝑃𝑃 se 𝑚𝑚𝓈𝓈 ( 𝜔𝜔 , 𝑡𝑡 ) ≈ ℛ 𝓈𝓈 𝑃𝑃 Tx ( 𝜔𝜔 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 × �� 𝜈𝜈𝜈𝜈𝑛𝑛=1 𝑊𝑊 𝑛𝑛 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑛𝑛𝓈𝓈 𝑟𝑟 𝑛𝑛𝓈𝓈 𝐷𝐷 ( 𝜃𝜃 𝑛𝑛𝓈𝓈 , 𝑘𝑘 ) δ ( 𝑌𝑌 𝓈𝓈 , 𝑟𝑟 𝑛𝑛𝓈𝓈 , 𝑘𝑘 ) 𝑒𝑒 𝑖𝑖𝑖𝑖Δ𝜏𝜏 𝑛𝑛 � × 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑚𝑚𝓈𝓈 𝑟𝑟 𝑚𝑚𝓈𝓈 𝐷𝐷 ( 𝜃𝜃 𝑚𝑚𝓈𝓈 , 𝑘𝑘 ) δ ( 𝑌𝑌 𝓈𝓈 , 𝑟𝑟 𝑚𝑚𝓈𝓈 , 𝑘𝑘 ). (31) The expression (31) is the acoustic pressure backscattered by a single scatterer and received by the 𝑚𝑚 th sub-element. Assum-ing now that there is a total of 𝑆𝑆 scatterers, the combination of their independent effect (single scattering assumption) gives the total pressure received by the 𝑚𝑚 th sub-element: 𝑃𝑃 se 𝑚𝑚 ( 𝜔𝜔 , 𝑡𝑡 ) ≈ 𝑃𝑃 Tx ( 𝜔𝜔 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 � 𝑆𝑆𝓈𝓈=1 �ℛ 𝓈𝓈 × �� 𝜈𝜈𝜈𝜈𝑛𝑛=1 𝑊𝑊 𝑛𝑛 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑛𝑛𝓈𝓈 𝑟𝑟 𝑛𝑛𝓈𝓈 𝐷𝐷 ( 𝜃𝜃 𝑛𝑛𝓈𝓈 , 𝑘𝑘 ) δ ( 𝑌𝑌 𝓈𝓈 , 𝑟𝑟 𝑛𝑛𝓈𝓈 , 𝑘𝑘 ) 𝑒𝑒 𝑖𝑖𝑖𝑖Δ𝜏𝜏 𝑛𝑛 � × 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑚𝑚𝓈𝓈 𝑟𝑟 𝑚𝑚𝓈𝓈 𝐷𝐷 ( 𝜃𝜃 𝑚𝑚𝓈𝓈 , 𝑘𝑘 ) δ ( 𝑌𝑌 𝓈𝓈 , 𝑟𝑟 𝑚𝑚𝓈𝓈 , 𝑘𝑘 ) � . (32) The expression (32) is for an individual sub-element 𝑚𝑚 . The pressure wave 𝑃𝑃 e ( 𝜔𝜔 , 𝑡𝑡 ) received by one transducer element is the coherent sum of the pressures received by all its sub-ele-ments (Fig. 6). For the element � 𝑛𝑛𝜈𝜈 � : 𝑃𝑃 e � 𝑛𝑛𝜈𝜈 � ( 𝜔𝜔 , 𝑡𝑡 ) = � 𝑃𝑃 se 𝑛𝑛+𝜇𝜇 ( 𝜔𝜔 , 𝑡𝑡 ) 𝜈𝜈−1𝜇𝜇=0 (33) The theoretical pressures were all derived for a single angular This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] frequency 𝜔𝜔 = 2 𝜋𝜋𝑓𝑓 . The full-band waveforms can be obtained by summation in the frequency domain through an inverse fast Fourier transform. B. Radiofrequency signals
The spectrum of the radiofrequency RF signal of the 𝑚𝑚 th sub-element is related to the received acoustic pressure (32) by RF se 𝑚𝑚 ( 𝜔𝜔 , 𝑡𝑡 ) ∝ �𝒮𝒮 T ( 𝜔𝜔 ) 𝑃𝑃 se 𝑚𝑚 ( 𝜔𝜔 , 𝑡𝑡 ) , (34) where it is recalled that �𝒮𝒮 T ( 𝜔𝜔 ) is the one-way transducer response. Alike (33), the RF signal related to one element is the coherent sum of the RF signals of its sub-elements. Fig. 8. Focused pressure fields simulated with PFIELD for a 64-element 2.7-MHz P4-2v cardiac phased array (kerf width = 50 μm, pitch = 0.3 mm, fractional bandwidth = 74%, elevation focus = 6 cm). These RMS (root mean square) acoustic fields illustrate emission sequences such as those used in standard transthoracic echocardiography, with focusing in the axial direction. V. T HE TWO - DIMENSIONAL CASE
In a two-dimensional (2-D) 𝑥𝑥 - 𝑧𝑧 domain, the piston-like ele-ment generates a normal velocity that is constant everywhere (i.e. in [ −∞ , + ∞ ] ) in the 𝑦𝑦 -direction. This situation is obtained when ℎ (element height) and 𝑅𝑅 𝑓𝑓 (distance to elevation focus) both tend to + ∞ . In this limit case, the rectangular function in Eq. (12) becomes for any 𝑦𝑦 ′ , and the coefficients reduce to 𝐺𝐺 = 1 , with 𝐴𝐴 = 1 and 𝐵𝐵 = 0 . Under these conditions, 𝛽𝛽 /(4 𝛼𝛼 ) + 𝛾𝛾 = 0 and √ ( 𝜋𝜋 / 𝛼𝛼 ) = √ (2 𝑖𝑖𝜋𝜋𝑟𝑟 / 𝑘𝑘 ) [see Eq. (14) and (15)]. In 2-D, from (17), the acoustic pressure generated by a single element thus reads 𝑃𝑃 - 𝐷𝐷 ( 𝒙𝒙 , 𝜔𝜔 , 𝑡𝑡 ) ≈� 𝑖𝑖𝑏𝑏𝑖𝑖𝑖𝑖 𝜌𝜌𝑐𝑐𝑣𝑣 ( 𝜔𝜔 ) � 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 𝑟𝑟 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 sinc( 𝑘𝑘𝑏𝑏 sin 𝜃𝜃 ) � , (35) which can be rewritten as 𝑃𝑃 - 𝐷𝐷 ( 𝒙𝒙 , 𝜔𝜔 , 𝑡𝑡 ) ≈√ 𝑏𝑏 �� 𝑖𝑖𝑏𝑏𝑖𝑖𝑖𝑖 𝜌𝜌𝑐𝑐𝑣𝑣 ( 𝜔𝜔 ) ������������ 𝑃𝑃 Tx ( 𝜔𝜔 ) 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 √𝑟𝑟 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 sinc( 𝑘𝑘𝑏𝑏 sin 𝜃𝜃 ) ��������� 𝐷𝐷 ( 𝜃𝜃 , 𝑘𝑘 ) . (36) The acoustic wave (22) radiated by an 𝑁𝑁 -element array then becomes 𝑃𝑃 - 𝐷𝐷 ( 𝑿𝑿 , 𝜔𝜔 , 𝑡𝑡 ) ≈ √ 𝑏𝑏 𝑃𝑃 Tx ( 𝜔𝜔 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 � 𝑊𝑊 𝑛𝑛𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟𝑛𝑛 �𝑟𝑟 𝑛𝑛 𝐷𝐷 ( 𝜃𝜃 𝑛𝑛 , 𝑘𝑘 ) 𝑒𝑒 𝑖𝑖𝑖𝑖Δ𝜏𝜏 𝑛𝑛 𝜈𝜈𝜈𝜈𝑛𝑛=1 . (37) Similarly, the pressure received by the 𝑚𝑚 th sub-element (32) in a 2-D space reduces to 𝑃𝑃 - Dse 𝑚𝑚 ( 𝜔𝜔 , 𝑡𝑡 ) ≈ (2 𝑏𝑏 ) 𝑃𝑃 Tx ( 𝜔𝜔 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 � 𝑆𝑆𝓈𝓈=1 �ℛ 𝓈𝓈 × �� 𝜈𝜈𝜈𝜈𝑛𝑛=1 𝑊𝑊 𝑛𝑛𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟𝑛𝑛𝓈𝓈 �𝑟𝑟 𝑛𝑛𝓈𝓈 𝐷𝐷 ( 𝜃𝜃 𝑛𝑛𝓈𝓈 , 𝑘𝑘 ) 𝑒𝑒 𝑖𝑖𝑖𝑖Δ𝜏𝜏 𝑛𝑛 � × 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟𝑚𝑚𝓈𝓈 �𝑟𝑟 𝑚𝑚𝓈𝓈 𝐷𝐷 ( 𝜃𝜃 𝑚𝑚𝓈𝓈 , 𝑘𝑘 ) � . (38) VI. B AFFLE IMPEDANCE AND D ISPERSIVE MEDIUM A. Finite impedance baffle
The general expressions (22) and (32) were derived by as-suming that the baffle in which the transducer element is em-bedded has an infinite acoustic impedance (rigid baffle). It might be recommended to use a finite baffle impedance to ob-tain more realistic radiation patterns [18], [19]. In such a case, an obliquity factor must be added to the element directivity 𝐷𝐷 ( 𝜃𝜃 , 𝑘𝑘 ) [see Eq. (36)]. The expressions of the obliquity factors This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] for a nonnegative finite impedance are given in [18] (null im-pedance = “soft” baffle) and [19] (positive impedance). The el-ement directivity becomes rigid: 𝐷𝐷 ( 𝜃𝜃 , 𝑘𝑘 ) = sinc( 𝑘𝑘𝑏𝑏 sin 𝜃𝜃 ) . soft: 𝐷𝐷 ( 𝜃𝜃 , 𝑘𝑘 ) = sinc( 𝑘𝑘𝑏𝑏 sin 𝜃𝜃 ) cos 𝜃𝜃 . otherwise: 𝐷𝐷 ( 𝜃𝜃 , 𝑘𝑘 ) = sinc( 𝑘𝑘𝑏𝑏 sin 𝜃𝜃 ) cos 𝜃𝜃cos 𝜃𝜃+𝜁𝜁 , (39) with 𝜁𝜁 standing for the ratio between the medium and baffle impedances. For a baffle of impedance 2.8 MRayl (epoxy) ad-jacent to soft tissues of impedance 1.6 MRayl, 𝜁𝜁 = 1.75 [19]. By default, the baffle is soft in PFIELD and SIMUS. Fig. 9. Plane-wave pressure field simulated with PFIELD for a 128-element 7.6-MHz L11-5v linear array (kerf width = 30 μm, pitch = 0.27 mm, frac- tional bandwidth = 77%, elevation focus = 1.8 cm). This RMS (root mean square) acoustic field illustrates an emission sequence (here, tilt angle = 10 o ) such as that used in “ultrafast” compound plane-wave imaging, without fo-cusing in the axial direction [31]. B. Attenuation
In addition to the distance-dependent decrease in amplitude caused by the inverse law (or inverse square-root law in 2-D), an ultrasound wave propagating in tissues is attenuated through scattering and absorption. In SIMUS, scattering is governed by the reflectivity coefficients ℛ 𝓈𝓈 . In the proposed ultrasound sim-ulator, the medium is assumed non-dispersive to preserve line-arity, which means that waves of different wavelengths travel at the same phase velocity ( = 𝑐𝑐 ). In PFIELD and SIMUS, ab-sorption can be included in the amplitude by adding a fre-quency-dependent multiplicative loss term. For numerical rea-sons (a recursive multiplication is used to calculate the expo-nential terms), this frequency dependence is linear in PFIELD. Amplitude absorption is obtained by substituting 𝑒𝑒 𝑖𝑖 ( 𝑖𝑖+𝑖𝑖𝑖𝑖 𝑎𝑎 ) 𝑟𝑟 for 𝑒𝑒 𝑖𝑖𝑖𝑖𝑟𝑟 , where ( 𝑘𝑘 𝑎𝑎 ) is given by: 𝑘𝑘 𝑎𝑎 = 𝛼𝛼 dB . . (40) The attenuation coefficient 𝛼𝛼 dB is in dB/cm/MHz. A typical value for soft tissues is 0.5 dB/cm/MHz [32]. The accounts for the unit conversion (cm·MHz to m·Hz). Fig. 10. RF signals simulated with SIMUS for a P4-2v phased array trans-mitting a focused wave in a medium that contains five scatterers.
VII.
RMS
PRESSURE FIELDS
Equation (22) gives the acoustic pressure field for a given angular frequency 𝜔𝜔 . By default, the PFIELD Matlab code re-turns the root-mean-square RMS pressure field. If we omit the ( 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 ) term in (22) and define 𝑃𝑃 ( 𝑿𝑿 , 𝜔𝜔 ) by the relationship 𝑃𝑃 ( 𝑿𝑿 , 𝜔𝜔 , 𝑡𝑡 ) ≝ 𝑃𝑃 ( 𝑿𝑿 , 𝜔𝜔 ) 𝑒𝑒 −𝑖𝑖𝑖𝑖𝑡𝑡 , the RMS pressure field is given by 𝑃𝑃 RMS ( 𝑿𝑿 ) = �∫ 𝑃𝑃 ( 𝑿𝑿 , 𝜔𝜔 ) d 𝜔𝜔 𝑐𝑐 . (41) The integral can be estimated by using a midpoint Riemann sum with 𝑁𝑁 𝑓𝑓 equally spaced frequency samples: 𝑃𝑃 RMS ( 𝑿𝑿 ) ≈ �Δ𝜔𝜔 � 𝑃𝑃 �𝑿𝑿 , 2 � 𝑗𝑗𝑁𝑁𝑓𝑓 � 𝜔𝜔 𝑘𝑘 � 𝑓𝑓 𝑗𝑗=0 . (42) This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] The partition width
Δ𝜔𝜔 must be small enough to ensure a proper approximation, but not too small to avoid computational overload with a large 𝑁𝑁 𝑓𝑓 . The angular frequency step is chosen so that the phase increment in (20) be smaller than 𝜋𝜋 every-where in the radiated region of interest ( roi ), i.e. Δ𝜔𝜔 must verify [( Δ𝜔𝜔 𝑐𝑐⁄ ) 𝑟𝑟 + Δ𝜔𝜔Δ𝜏𝜏 ] < 2 𝜋𝜋 for any distance 𝑟𝑟 and transmit delay Δ𝜏𝜏 , i.e.
Δ𝜔𝜔 = min roi {2 𝜋𝜋 /( 𝑟𝑟 / 𝑐𝑐 + Δ𝜏𝜏 )} . Once the transmit delays Δ𝜏𝜏 𝑛𝑛 and the parameters of the array are given, the equations (22) and (42) yield the RMS pressure field at the location specified by 𝑿𝑿 . Fig. 8 and Fig. 9 illustrate two examples of acoustic fields simulated with PFIELD: fo-cused and plane wavefronts with a cardiac and linear array, re-spectively. In these examples, the 3-D equation was applied to take the elevation focusing into account. Fig. 11. Simulation of a three-chamber-view echocardiographic image. RF signals were simulated with SIMUS by using 39,500 scatterers with prede-fined reflection coefficients (top row). 128 focused waves were transmitted to create the 128 scanlines of the B-mode image.
VIII. M EDICAL ULTRASOUND IMAGES WITH
SIMUS Medical ultrasound images, such as B-mode or color Doppler images, can be generated by simulating RF signals through Eq. (32). Once RF signals are obtained, they can be post-processed (e.g. using I/Q demodulation, beamforming, clutter filtering, phase analysis, envelope detection…) by using the Matlab codes available in the MUST toolbox . A series of RF signals given by SIMUS in a simplified configuration (a focused wave that radiates five scatterers) is displayed on Fig. 10. Realistic ultrasound images can be obtained when using a large number of point-like scatterers, as explained in the next paragraphs. A. B-mode imaging
Fig. 11 illustrates the simulation of a long-axis echocardio-graphic image (three-chamber view) by scanning the left heart with a series of focused waves. The fake myocardial tissue con-tained 39,500 randomly distributed point scatterers (density = 10 scatterers per wavelength squared) whose reflection coeffi-cients are shown in Fig. 11. The fan-shaped B-mode image con-sists of 128 scanlines. Each scanline was constructed from one tilted focused wave generated by a 2.7-MHz phased array (64 elements, 0.3-mm pitch, 74% fractional bandwidth, 6-cm ele-vation focus). RF signals were simulated with SIMUS, then I/Q demodulated and beamformed [33]. The B-mode image was ob-tained by log-compressing the real envelopes. As a comparison, Fig. 12 represents the B-mode image sim-ulated with the same scatterers and transducer, but with 32 tilted 50 o -wide diverging waves (tilt angle ranging between -25 o and +25 o ). This cardiac image was obtained by summing the 32 frames coherently after beamforming (compounding). Note that these simulations were not entirely realistic as there was no movement of the scatterers from one frame to the next. Such displacements can highly affect the quality of the compound image in the absence of motion compensation [1]. Fig. 12. Simulation of a three-chamber-view echocardiographic image. RF signals were simulated with SIMUS by using 39,500 scatterers with prede-fined reflection coefficients (see Fig. 11). 32 plane waves were transmitted to create the compound B-mode image. B. Color flow imaging and Vector flow imaging
SIMUS (and the tools available in the MUST toolbox) can also be used to simulate realistic color flow images (color Dop-pler) and vector flow images. Two-dimensional simulations of color and vector Doppler in a 2-D carotid bifurcation were in-troduced in a previous work that combined CFD (computational
This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] fluid dynamics) by SPH (smooth particle hydrodynamics) with SIMUS [23]. In that paper, radiofrequency RF signals were simulated by SIMUS using a plane-wave-imaging sequence (non-tilted waves transmitted by a 128-element linear array). After I/Q demodulation, and beamforming with two distinct sub-apertures (see details in [23]), a one-lag slow-time autocor-relator was applied to recover 2-D velocity fields. Fig. 13 pro-vides two snapshots adapted from [23]. This vector-Doppler methodology was successfully used in vivo with sub-Nyquist RF sampling [34]. Fig. 13. Simulation of color Doppler and vector Doppler in a 2-D carotid bifurcation. The displacements of ~26,000 blood scatterers (top row) were simulated by SPH (smoothed particle hydrodynamics), a Lagrangian CFD method [23]. The inset is at another time. RF signals were simulated with SIMUS by using unsteered plane waves. Color Doppler (bottom row, red-blue patterns) and vector Doppler (bottom row, arrows) were generated by post-processing the I/Q signals (beamforming, autocorrelator, robust smoothing) with tools available in the MUST toolbox. Adapted from [23].
IX. C ONCLUSION
The assumptions and simplifications included in PFIELD and SIMUS make its theory and numerical resolution in the Fourier domain conveninet. The examples show that sufficiently realis-tic ultrasound images can be created for educational and re-search purposes. How PFIELD compares to Field-II, k-Wave and Verasonics is detailed in the second part of this article. The current version of PFIELD (2021), although it includes the 3-D acoustic equation for elevation focusing, is limited to one-di-mensional, linear or convex ultrasonic transducers. A volume version is planned. Nevertheless, based on the equations de-scribed in this paper, and keeping the same hypotheses, an ad-vanced user could adapt SIMUS for the simulation of matrix arrays and volumetric acoustic fields. √𝑎𝑎 + 𝑥𝑥 = √𝑎𝑎 + 𝑥𝑥 � √𝑎𝑎� + O( 𝑥𝑥 ) ⁄ A CKNOWLEDGMENT
I thank my former students from Montreal and Lyon who tested and re-tested the different versions, and contributed to the im-provement of the MUST toolbox. In particular: Jonathan Porée, Daniel Posada, Julia Faurie, Craig Madiena, Vincent Perrot. As well as my colleagues from CREATIS, Lyon, Olivier Bernard, and François Varray for making the most of MUST. A
PPENDIX A. Paraxial and far-field distances
Equation (2) can be rewritten as 𝑟𝑟 ′ = 𝑧𝑧�� ( 𝑥𝑥−𝑥𝑥 ′ ) 𝑧𝑧 � + ( 𝑦𝑦−𝑦𝑦 ′ ) 𝑧𝑧 . (43) In the paraxial (Fresnel) approximation, ( 𝑦𝑦 − 𝑦𝑦 ′ ) / 𝑧𝑧 ≪ [1 + ( 𝑥𝑥 − 𝑥𝑥 ′ ) / 𝑧𝑧 ] . One can thus write 𝑟𝑟 ′ ≈ �𝑧𝑧 + ( 𝑥𝑥 − 𝑥𝑥 ′ ) + �𝑦𝑦−𝑦𝑦 ′ � + ( 𝑥𝑥−𝑥𝑥 ′ ) . (44) In the paraxial approximation, one has also 𝑟𝑟 ≈ 𝑥𝑥 + 𝑧𝑧 . The expression (44) can thus be approximated by 𝑟𝑟 ′ ≈ 𝑟𝑟� 𝑥𝑥 ′2 𝑟𝑟 − 𝑥𝑥𝑥𝑥 ′ 𝑟𝑟 + �𝑦𝑦−𝑦𝑦 ′ � 𝑥𝑥′2𝑟𝑟2 −2 𝑥𝑥𝑥𝑥′𝑟𝑟2 . (45) By keeping the 1 st order of 𝑥𝑥 ′ in the first square root and the 0 th order in the second (far-field approximation), we obtain 𝑟𝑟 ′ ≈ 𝑟𝑟 � − 𝑥𝑥𝑥𝑥 ′ 𝑟𝑟 � + �𝑦𝑦−𝑦𝑦 ′ � . (46) Because sin 𝜃𝜃 = 𝑥𝑥 / √𝑥𝑥 + 𝑧𝑧 ≈ 𝑥𝑥 / 𝑟𝑟 (Fig. 3), (46) reduces to Eq. (6). B. Spectrum of the transducer PSF
We need a 𝜔𝜔 𝑏𝑏 bandwidth at -6 dB. Therefore, from (26): 𝒮𝒮 t �𝜔𝜔 𝑘𝑘 ± 𝑖𝑖 𝑏𝑏 � = 𝑒𝑒 −� 𝑖𝑖 𝑏𝑏 𝑐𝑐 �𝑝𝑝 = (47) One can thus deduce 𝜎𝜎 : 𝜎𝜎 = 𝑝𝑝 𝜔𝜔 𝑏𝑏 𝜔𝜔 𝑘𝑘 . (48) By assuming that 𝒮𝒮 t (0) = 𝒮𝒮 t (2 𝜔𝜔 𝑘𝑘 ) = 𝑒𝑒 − = 2 −126 , (49) the expressions (48) and (49) yield 𝑝𝑝 [Eq. (27)]. This v1 version (2021-02-04) has not yet been reviewed by independent peers. It is for arXiv.org ] R EFERENCES [1] J. Porée, D. Posada, A. Hodzic, F. Tournoux, G. Cloutier, and D. Gar-cia, “High-frame-rate echocardiography using coherent compounding with Doppler-based motion-compensation,”
IEEE Transactions on Medical Imaging , vol. 35, no. 7, pp. 1647–1657, 2016, doi: 10.1109/TMI.2016.2523346. [2] E. Evain, K. Faraz, T. Grenier, D. Garcia, M. De Craene, and O. Ber-nard, “A pilot study on convolutional neural networks for motion esti-mation from ultrasound images,”
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control , pp. 1–9, 2020, doi: 10.1109/TUFFC.2020.2976809. [3] E. Roux, A. Ramalli, P. Tortoli, C. Cachard, M. C. Robini, and H. Liebgott, “2-D Ultrasound sparse arrays multidepth radiation optimi-zation using simulated annealing and spiral-array inspired energy functions,”
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control , vol. 63, no. 12, pp. 2138–2149, 2016, doi: 10.1109/TUFFC.2016.2602242. [4] H. Liebgott, A. Rodriguez-Molares, F. Cervenansky, J. A. Jensen, and O. Bernard, “Plane-wave imaging challenge in medical ultrasound,”
IEEE International Ultrasonics Symposium (IUS) , pp. 1–4, 2016, doi: 10.1109/ULTSYM.2016.7728908. [5] J. A. Jensen, “FIELD: A program for simulating ultrasound systems,” in , 1996, pp. 351–353. [6] J. A. Jensen and P. Munk, “Computer phantoms for simulating ultra-sound B-mode and CFM images,” in
Acoustical Imaging , Springer, Boston, MA, 1997, pp. 75–80. [7] B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the sim-ulation and reconstruction of photoacoustic wave fields,”
JBO, JBOPFO , vol. 15, no. 2, p. 021314, 2010, doi: 10.1117/1.3360308. [8] B. E. Treeby, J. Jaros, D. Rohrbach, and B. T. Cox, “Modelling elastic wave propagation using the k-Wave MATLAB Toolbox,” in , 2014, pp. 146–149, doi: 10.1109/ULTSYM.2014.0037. [9] R. J. McGough, “Rapid calculations of time-harmonic nearfield pres-sures produced by rectangular pistons,”
The Journal of the Acoustical Society of America , vol. 115, no. 5, pp. 1934–1941, Apr. 2004, doi: 10.1121/1.1694991. [10] M. E. Frijlink, H. Kaupang, T. Varslot, and S.-E. Masoy, “Abersim: a simulation program for 3D nonlinear acoustic wave propagation for arbitrary pulses and arbitrary transducer geometries,” in , 2008, pp. 1282–1285, doi: 10.1109/ULTSYM.2008.0310. [11] E. Bossy, M. Talmant, and P. Laugier, “Three-dimensional simula-tions of ultrasonic axial transmission velocity measurement on cortical bone models,”
The Journal of the Acoustical Society of America , vol. 115, no. 5, pp. 2314–2324, Apr. 2004, doi: 10.1121/1.1689960. [12] M. Rao, T. Varghese, and J. A. Zagzebski, “Simulation of ultrasound two-dimensional array transducers using a frequency domain model,”
Medical Physics , vol. 35, no. 7Part1, pp. 3162–3169, 2008, doi: 10.1118/1.2940158. [13] B. Piwakowski and K. Sbai, “A new approach to calculate the field ra-diated from arbitrarily structured transducer arrays,”
IEEE Transac-tions on Ultrasonics, Ferroelectrics, and Frequency Control , vol. 46, no. 2, pp. 422–440, 1999, doi: 10.1109/58.753032. [14] F. Varray, O. Basset, P. Tortoli, and C. Cachard, “CREANUIS: a non-linear radiofrequency ultrasound image simulator,”
Ultrasound in Medicine and Biology , vol. 39, no. 10, pp. 1915–1924, 2013, doi: 10.1016/j.ultrasmedbio.2013.04.005. [15] A. Freedman, “Sound field of a rectangular piston,”
The Journal of the Acoustical Society of America , vol. 32, no. 2, pp. 197–209, 1960, doi: 10.1121/1.1908013. [16] P. R. Stepanishen, “Transient radiation from pistons in an infinite pla-nar baffle,”
The Journal of the Acoustical Society of America , vol. 49, no. 5B, pp. 1629–1638, 1971, doi: 10.1121/1.1912541. [17] J. Marini and J. Rivenez, “Acoustical fields from rectangular ultra-sonic transducers for non-destructive testing and medical diagnosis,”
Ultrasonics , vol. 12, no. 6, pp. 251–256, 1974, doi: 10.1016/0041-624X(74)90132-2. [18] A. R. Selfridge, G. S. Kino, and B. T. Khuri ‐Yakub, “A theory for the radiation pattern of a narrow‐strip acoustic transducer,”
Appl. Phys. Lett. , vol. 37, no. 1, pp. 35–36, 1980, doi: 10.1063/1.91692. [19] P. Pesque and M. Fink, “Effect of the planar baffle impedance in acoustic radiation of a phased array element. Theory and experimenta-tion,” in
IEEE 1984 Ultrasonics Symposium , 1984, pp. 1034–1038, doi: 10.1109/ULTSYM.1984.198462. [20] J. J. Wen and M. A. Breazeale, “A diffraction beam field expressed as the superposition of Gaussian beams,”
The Journal of the Acoustical Society of America , vol. 83, no. 5, pp. 1752–1756, 1988, doi: 10.1121/1.396508. [21] K. B. Ocheltree and L. A. Frizzel, “Sound field calculation for rectan-gular sources,”
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control , vol. 36, no. 2, pp. 242–248, 1989, doi: 10.1109/58.19157. [22] L. W. Schmerr Jr,
Fundamentals of ultrasonic phased arrays . Springer International Publishing, 2015. [23] S. Shahriari and D. Garcia, “Meshfree simulations of ultrasound vec-tor flow imaging using smoothed particle hydrodynamics,”
Phys Med Biol , vol. 63, pp. 1–12, 2018, doi: 10.1088/1361-6560/aae3c3. [24] J. A. Jensen et al. , “SARUS: A synthetic aperture real-time ultrasound system,”
IEEE Transactions on Ultrasonics, Ferroelectrics, and Fre-quency Control , vol. 60, no. 9, pp. 1838–1852, 2013, doi: 10.1109/TUFFC.2013.2770. [25] J. Yu, H. Yoon, Y. M. Khalifa, and S. Y. Emelianov, “Design of a volumetric imaging sequence using a Vantage-256 ultrasound re-search platform multiplexed with a 1024-element fully sampled matrix array,”
IEEE Transactions on Ultrasonics, Ferroelectrics, and Fre-quency Control , vol. 67, no. 2, pp. 248–257, 2020, doi: 10.1109/TUFFC.2019.2942557. [26] B. G. Lucas and T. G. Muir, “The field of a focusing source,”
The Journal of the Acoustical Society of America , vol. 72, no. 4, pp. 1289–1296, 1982, doi: 10.1121/1.388340. [27] J. A. Jensen, “Ultrasound imaging and its modeling,” in
Imaging of complex media with acoustic and seismic waves , M. Fink, W. A. Ku-perman, J.-P. Montagner, and A. Tourin, Eds. Berlin, Heidelberg: Springer, 2002, pp. 135–166. [28] J. A. C. Weideman, “Computation of the complex error function,”
SIAM J. Numer. Anal. , vol. 31, no. 5, pp. 1497–1518, 1994, doi: 10.1137/0731077. [29] W. Liu, P. Ji, and J. Yang, “Development of a simple and accurate ap-proximation method for the Gaussian beam expansion technique,”
The Journal of the Acoustical Society of America , vol. 123, no. 5, pp. 3516–3516, 2008, doi: 10.1121/1.2934433. [30] L. Y. L. Mo and R. S. C. Cobbold, “A stochastic model of the backscattered Doppler ultrasound from blood,”
IEEE Transactions on Biomedical Engineering , vol. BME-33, no. 1, pp. 20–27, 1986, doi: 10.1109/TBME.1986.325834. [31] J. Porée, D. Garcia, B. Chayer, J. Ohayon, and G. Cloutier, “Noninva-sive vascular elastography with plane strain incompressibility assump-tion using ultrafast coherent compound plane wave imaging,”
IEEE Transactions on Medical Imaging , vol. 34, no. 12, pp. 2618–2631, 2015, doi: 10.1109/TMI.2015.2450992. [32] K. J. Parker, R. M. Lerner, and R. C. Waag, “Attenuation of ultra-sound: magnitude and frequency dependence for tissue characteriza-tion.,”
Radiology , vol. 153, no. 3, pp. 785–788, 1984, doi: 10.1148/ra-diology.153.3.6387795. [33] V. Perrot, M. Polichetti, F. Varray, and D. Garcia, “So you think you can DAS? A viewpoint on delay-and-sum beamforming,”
Ultrasonics , vol. 111, p. 106309, 2021, doi: 10.1016/j.ultras.2020.106309. [34] C. Madiena, J. Faurie, J. Porée, and D. Garcia, “Color and vector flow imaging in parallel ultrasound with sub-Nyquist sampling,”