I Introduction
The need for fast and accurate us simulation tools is as great as ever. Such tools have been and continue to be extensively used for the design and characterization of us transducers [23, 15]. The use of end-to-end us scanner simulators is also increasing for the training and evaluation of physicians and sonographers [53, 25], as they eliminate the need for volunteers or patients and can provide on-demand exposure to specific care and diagnosis scenarios. Simulation tools are also crucial for the development, assessment, and validation of image analysis and image reconstruction methods. This includes, for example, the analysis and characterization of speckle patterns [50, 10], the optimization of application-specific acquisition sequences [16], the fine-tuning and benchmarking of image reconstruction parameters and methods [32], or the development of image quality metrics [33]. The recent breakthrough of deep learning-based methods in medical image analysis [12], medical image reconstruction [51, 52], and ultrasound imaging [48] comes with a critical need for large-scale datasets to feed the data-intensive algorithms involved. In this context, the generation of synthetic data is of great interest [11], with the potential to generate (infinitely) large, highly diverse, and unbiased datasets.
With the myriad of us applications in which simulations can be leveraged, it is no wonder that many us simulators are available.^{1}^{1}1A fairly exhaustive list of available software can be found on the kwave website: http://www.k-wave.org/acousticsoftware.php. Among them, and in the context of linear acoustics, fieldii [18, 21] is considered a reference, thanks to its ability to calculate acoustic fields of arbitrarily shaped, excited, and apodized transducers, as well as simulating the acquisition of ultrasound pulse-echo signals in the presence of (weak) tissue inhomogeneities [17]. It is based on the sir method developed by Tupholme [43] and Stepanishen [38, 37], which is an analytic (mesh-free) method dedicated to the evaluation of transient acoustic fields that are of particular interest for pulse-echo us imaging.
The sir approach amounts to evaluating a time-dependent surface integral derived from the Rayleigh-Sommerfeld equations. A detailed review of the method can be found in [13]. Once the sir is known, the resulting velocity potential can be computed by the time convolution of the surface excitation and the sir. Physical quantities such as pressure or particle velocity can eventually be obtained from the velocity potential. Much attention and efforts were devoted to the derivation of analytic sir expressions for uniformly excited radiators of specific shapes and boundary (baffle) conditions, resulting in complete expressions for the circular piston [13], the slit [24], the spherically focused radiator [2], the triangular piston [19], and the rectangular piston [34]. Semi-analytic expressions were also derived for more complex geometries such as the cylindrical shell [40] and the toroidal shell [3]. Non-uniform excitations were also investigated in [14, 39, 42, 49], but analytic expressions were restricted to specific transducer shapes with specific excitation amplitude distributions.
Analytic expressions only exist for a restricted set of transducer shapes, excitation distributions, and baffle conditions. To cope with this limitation, numerical methods were proposed. The main principle consists of representing the transducer surface as a set of characteristic sub-elements for which (simple) analytic expressions of the sir exist. The sir of the transducer is obtained by summing all the sub-element sir (superposition principle). Such a discretization also enables approximating surface apodization and delay by appropriately weighting and delaying the sir of each sub-element before summation. A popular characteristic sub-element type is the rectangle, as it enables reasonably good approximations of arbitrary surfaces, and more importantly because there exist computationally efficient far-field approximations for both rigid and soft baffle conditions [20]. Another relevant approach [30] consists of evaluating the surface integral numerically by discretizing the surface into ideal points. The sir of the transducer can then be computed by a simple weighted sum of Dirac delta functions. This strategy has the advantage of enabling exact representation of surfaces, at the cost of requiring a much larger number of sub-elements than the strategy deployed in fieldii [31].
A major numerical difficulty of the sir method comes from the signal properties of such responses. They are typically of very short duration and are characterized by abrupt slope changes induced by the edges and vertices of the transducer aperture, resulting in very high frequencies [26, 2].^{2}^{2}2A typical example of a sir and corresponding frequency spectrum is shown in supplementary material (Fig. S1). Yet, because the electromechanical impulse responses of conventional transducer elements are bandlimited, the excitation waveform, and most importantly the field signal of interest (eg, pressure field), are also bandlimited. Ideally, one would want to sample the sir at the same rate as the excitation waveform before time convolution of the two signals. This is obviously not possible in most cases and sampling at the rate of the sir (ie, several additional orders of magnitude) is computationally prohibitive. Different strategies based on the time integration of the sir were proposed to tackle this computational difficulty, either numerically using an adaptive sampling of the sir [2] or from analytic expressions [9, 20]. Such strategies enabled for more efficient sampling rates, but still require oversampling ratios of several factors to prevent from detrimental aliasing.
In summary, the most important needs for a computational method relying on the sir are: an accurate representation of arbitrary transducer shapes; an efficient way of evaluating the surface integral at each time instant; and an efficient sampling strategy such that the sir can be sampled at rates similar to those required for the excitation waveform. To address these needs, we propose to represent the shape of transducers as nurbs surfaces, to evaluate the surface integral numerically using high-order Gaussian quadrature rules, and to express the sir in bspline bases for efficient sampling of the time axis. The use of nurbs enables accurate representation of complex surfaces and exact representations of common transducer elements. High-order Gaussian quadrature rules (eg, Gauss-Legendre), require much fewer quadrature points to achieve a desired accuracy [1, Sec. 25.4]. Finally, relying on numerical integration allowed us to express the calculation of the sir in high-order bspline bases [46, 44, 41]. This enabled sampling the sir at low sampling rates, as required by the excitation waveform, without introducing additional errors on simulated field signals.
Ii Proposed Approach
Let us recall [38, 8] that the sir of a radiating surface (assumed to be embedded in an infinite planar baffle) can be expressed at a field point and at a time instant as
(1) |
where is the mean sound speed in the medium (assumed homogeneous), is the (inward) surface normal, represents the spatial distribution of surface velocity or pressure over the radiating surface, and is a term depending on the boundary conditions, expressed as
(2) |
Let us consider a generic excitation waveform imposed as boundary condition on the surface (eg, normal velocity). The velocity potential can then be expressed as from which physical quantities can be obtained. For instance, the pressure can be evaluated as . Thus, the acoustic propagation from a radiating surface can be interpreted as the convolution of the excitation waveform derivative and the sir. Note that is bandlimited as it includes the electromechanical impulse response of the radiator.
ii.1 Numerical quadrature of the spatial impulse response
At each time instant and field point , the integrand of Eq. 1 is a real-valued function that can define as for the purpose of derivation. Let us assume that the radiating surface is a smooth surface parametrized by a mapping whose jacobian determinant , , where the jacobian matrix is defined as . Thus, the surface integral of onto can be rewritten as
(3) |
Assuming that is a well-behaved function,^{3}^{3}3This is generally true as there is no realistic interest for the evaluation of the sir on the radiating surface, where the integrand of Eq. 1 is obviously singular. one can approximate Eq. 3 by means of Gaussian quadrature as the finite sum
(4) |
where are the quadrature coordinates (parametric space), for which , , and , are the corresponding quadrature weights, jacobian determinants, and quadrature points (physical space), respectively.
Using Eq. 4, we can rewrite Eq. 1 as
(5) |
where and . The normal vector at some surface point can be computed from the surface parametrization (mapping) as
(6) |
where are the coordinates defining the parametric space of the surface mapping. The jacobian determinants, represented as in Eq. 5, can then be computed from the normal vectors directly as
(7) |
Finally, by grouping all spatially dependent weighting terms in Eq. 5 under a global weighting term , we can obtain a compact expression for the evaluation of the sir by means of Gaussian quadrature as
(8) |
Thus, provided that such a surface mapping exists and that a suitable Gaussian quadrature rule is deployed with a sufficient amount of quadrature points to achieve a desired accuracy, the calculation of the sir can be well approximated by a sum of shifted-and-weighted Dirac delta functions.
ii.2 Non-uniform rational B-spline surface representations
We first review some basic principles of nurbs (surface) representations. For a detailed reference, the reader is referred to the book by Piegl and Tiller [29]. Let be a nondecreasing sequence of real numbers representing a nonperiodic, clamped, or open knot vector, with knots, defined as
(9) |
. The corresponding (nonnegative) bspline basis functions of degree are defined recursively as [5, 6, 7]
(10) | ||||
(11) |
Note that we restrict ourselves to nonperiodic knot vectors onto which bspline basis functions are interpolating at the endpoints of such knot vectors but are (in general) non-interpolating at interior knots. Also, definitions are not strictly limited to the
interval. Yet, it is so common to the nurbs community that it is also adopted in the present work. Typical examples of bspline basis functions on knot vectors of the form are shown in Fig. 1 for different degrees. This type of knot vectors results in bspline basis functions that are bernstein polynomials, because bspline representations are a generalization of bezier representations [29, Sec. 2.2].A nurbs surface of degree in directions can be represented by a bivariate vector-valued piecewise rational function (mapping) defined as [29, Sec. 4.4]
(12) |
where are the control points (forming a bidirectional net), are the corresponding weights, and are the (nonrational) bspline basis functions of degrees and that are defined on the nondecreasing and nonperiodic (ie, clamped) knot vectors
(13) | ||||
(14) |
respectively, where and . Expressions for the derivatives of nurbs surfaces also exist and can be found in [29, Sec. 4.5]. They are essential tools for the calculations of the surface normal vectors and the jacobian determinants defined in Eq. 6 and Eq. 7, respectively. Note that if both and are defined with no interior points, such as in the examples shown in Fig. 1, the nurbs surface is a (smooth) rational bezier surface.
An important property of nurbs representations is that any nurbs surface can be decomposed into a union of (smooth) rational bezier surfaces (or patches), that share at most a common edge or a common vertex. Such a decomposition can be performed by a simple procedure known as knot refinement [29, Sec. 5.3]. Each (decomposed) rational bezier patch is actually also a nurbs surface such that it can be readily represented by Eq. 12. Rational bezier patches are of particular interest as they are smooth surfaces as opposed to “general” nurbs surfaces that may contain breakpoints. This is an essential property to allow the use of Gaussian quadrature rules onto such surfaces (Sec. II.1).
Also, most transducer elements composing conventional us transducer arrays are (at most) quadric surfaces and nurbs surfaces can represent quadric surfaces exactly (suppl. mat., Sec. S-II.1).
ii.3 Spatial impulse response in B-spline bases
Thanks to the nurbs representation of the radiating surface, we can evaluate the global weigthing term from the sum of (shifted-and-weighted) Dirac delta functions to approximate the sir [Eq. 8]. Before diving into the well-known difficulty of sampling Dirac delta functions, it is important to keep in mind that the sir is mainly a physical concept as such a quantity cannot be measured in physical conditions due to the electromechanical impulse response of transducer elements that are bandlimited. Thus, signals of interest that will be eventually measured are also bandlimited. Such a signal can be expressed generically as
(15) |
where represents some bandpass excitation waveform. Using Eq. 8 we obtain the corresponding approximation of such a signal as
(16) |
Equation 16 tells us that the signal is composed of shifted-and-weighted replicas of the excitation waveform . The excitation waveform is typically measured in physical conditions or approximated using a model pulse. Thus, it can be expressed as a sampled signal of uniformly spaced samples , , with a sampling interval . Equation 16 can then be rewritten as the discrete convolution
(17) |
where is some (interpolating) basis function (sometimes referred to as sampling kernel), and where it is assumed that boundary conditions are handled properly (via signal extension of ). Since is a bandpass signal, a natural and error-free choice for would be the (normalized) sinc function, . But because the sinc function is of infinite support, it needs to be truncated by multiplying it by some finite-support window. Generally, truncated sinc functions require rather large supports to achieve high accuracy [41]. Other options are standard (finite-support) interpolating basis functions such as those deployed for nearest-neighbor and linear interpolation, which are computationally efficient but of low accuracy.
Instead, as the discrete convolution of Eq. 17 consists of shifting and interpolating , we propose to rely on the concept of generalized interpolation [41]. By doing so, we can rewrite Eq. 17 as
(18) |
where is some basis function and are the corresponding basis coefficients defined such that
(19) |
As a result, one can rely on potentially non-interpolating basis functions, possessing better properties than interpolating ones, with an additional operation that consists of finding the basis coefficients. A particularly elegant and efficient way of finding these coefficients is yet another (pre-filtering) convolution operation with the convolution-inverse^{4}^{4}4Not to be confused with the inverse function of . as
(20) |
provided that exists. In the case of interpolating basis functions, the sequence of coefficients is equal to the sequence of samples .
We will restrict ourselves to the family of bspline basis functions as they benefit from maximal approximation orders for given supports [46, 44, 41]. For the purpose of interpolation, we can express a bspline basis function of degree as [35, 36]
(21) |
, , where the one-sided power function is defined as [41]
(22) |
We are typically interested in bspline basis functions of degree greater than one, since is almost identical to the nearest-neighbor basis function and
corresponds to the basis function of linear interpolation. Also, bspline basis functions of even degrees are usually not computationally interesting because they have the same support as bspline basis functions of the following (odd) degrees. Since bspline basis functions are symmetric, the pre-filtering operation defined in
Eq. 20 can be efficiently performed by a series of consecutive causal and anti-causal iir filters with poles and , respectively [45, 46, 47]. The pairs of poles can be derived from the z-transform of the convolution-inverse , many of which are tabulated in [45, 46, 47].So far we have assumed that the samples and corresponding coefficients were extended properly for the purpose of convolution operations. Because the corresponding excitation waveform is bandlimited and can be modeled as a windowed sinusoidal rf signal of finite support,^{5}^{5}5Acoustic pulses are theoretically infinite responses with (rapidly) decaying trailing oscillations. Thus, they may be truncated when their trailing pulse envelope falls below some level to achieve a desired accuracy. we consider zero boundary conditions. Also, since we want to work with finite support sequences to perform the discrete convolution in Eq. 18, we impose zero boundary conditions to the coefficients directly, namely , .
Once these coefficients are computed, they can be used readily for all quadrature points to evaluate the signal of interest at any field point using Eq. 18, namely as a discrete (full) convolution of the (pre-filtered) coefficients and the summation of shifted-and-weighted basis functions wrt the propagation times and weights , respectively. An illustration of the convolution involved in the evaluation of Eq. 18 is described in supplementary material (Sec. S-II.2). For the purpose of interpretation (and compactness), let us define the basis sir as
(23) |
such that Eq. 18 can be interpreted as the convolution of the (pre-filtered) coefficients and the basis sir, namely
(24) |
Even if is not a quantity that can be measured, and as such is less important to approximate accurately than , it remains interesting to be able to evaluate it.^{6}^{6}6Note that it is also possible to obtain the field response to a continuous-wave excitation directly from the sir by evaluating its fourier transform. Since is a linear combination of basis functions [Eq. 23], it is possible to obtain an approximation of from with the convolution-inverse as
(25) |
This is similar to expressing a cardinal spline basis function from its corresponding (non-interpolating) bspline basis function. Thus, the approximation of the sir obtained from Eq. 25 when considering non-interpolating basis functions will be of infinite support with rapidly decaying oscillations [41].
ii.4 Complete process and implementation details
Figure 2 summarizes the complete process of the proposed approach developed in Secs. II.3, II.2 and II.1 for the approximation of (bandpass) field signals using the sir of radiating surfaces. The cylindrical shape of the transducer element [Fig. ((a))(a)] is represented (exactly) as a nurbs surface using Eq. 12. As this nurbs surface is also a smooth bezier surface, a Gauss-Legendre quadrature rule can be defined in the parametric space directly [Fig. ((b))(b)], The corresponding eight quadrature points can be mapped into the physical space using the nurbs representation. The normal vectors and jacobian determinants can be computed using Eq. 6 and Eq. 7, respectively. From the distances between the quadrature points and the field point, some basis function (eg, cubic bspline) can be evaluated at the corresponding time instants on a sampled time axis of minimum support [Fig. ((c))(c)], and weighted accordingly. Their summation [Eq. 23] results in the basis sir [Fig. ((d))(d)]. The (pre-filtered) coefficients [Fig. ((f))(f)] are obtained by a series of causal and anti-causal iir filters applied to the excitation waveform. Finally, the field signal [Fig. ((g))(g)] is obtained by the convolution of the (pre-filtered) coefficients [Fig. ((f))(f)] and the basis sir [Fig. ((d))(d)]. This process must be repeated for all field points of interest, except for the coefficients that only need to be computed once.
The extension of the proposed approach to transducer arrays and pulse-echo acquisitions is straightforward and can be efficiently implemented (suppl. mat., Sec. S-II.3).
Iii Experiments and Results
To validate the proposed approach, we performed two numerical experiments (Secs. III.2 and III.1). The goal of the first one is to validate the core of the proposed approach, as described in Eq. 24, namely a convolution of (pre-filtered) coefficients and a signal (representing the basis sir) composed of shifted-and-weighted basis functions. The second experiment consists of evaluating the accuracy of the complete proposed approach on field signals radiated by transducer elements with specific shapes allowing analytic expressions for the sir.
For the two experiments, the pulse model considered was an analytic expression of the time derivative of a log-normal-modulated sinusoidal rf pulse (suppl. mat., Sec. S-III
). For the log-normal distribution
, we used the parameters and . The frequency of the sinusoidal was set to , resulting in a waveform centered at with a bandwidth of at [Fig. ((b))(b)]. The waveform was truncated at a trailing pulse envelope level of (double-precision floating-point format), resulting in a duration of . The fwhm of the resulting waveform corresponds to .In addition to the bspline basis functions, we also considered o-moms functions [4], that are derived from bspline basis functions. In general, we compared the classical interpolating basis functions for nearest-neighbor, linear, and quadratic^{7}^{7}7By quadratic, it is implied that the optimal parameter was used for the keys basis function, resulting in a quadratic interpolation method of order three. keys [22], as well as bspline and o-moms basis functions of different degrees. For the purpose of quantifying the accuracy in the following experiments, we relied on the relative two-norm error defined between an estimated (rf signal) vector and its true counterpart as
(26) |
iii.1 Convergence order of various basis functions
We considered a reference signal consisting of the time convolution between the analytic waveform considered [Fig. ((a))(a)] and a stream of Dirac delta functions at random-uniform times and of random-normal amplitudes . The corresponding reference signal can be expressed as
(27) |
which can be seen as a generic expression for Eq. 16, namely signals we want to approximate using the proposed approach Eq. 24. Note that Eq. 27 can be evaluated exactly using the analytic expression for the excitation waveform.
We considered a signal duration corresponding to times the time resolution cell of the waveform considered and populated it with an average of random Diracs per time resolution cell. We compared three interpolating basis functions, namely nearest-neighbor (degree zero), linear (degree one), and quadratic keys (degree two), four (non-interpolating) bspline basis functions of degrees , , , and , and the (non-interpolating) o-moms basis function of degree three. To validate that the theoretical approximation order of these basis functions (ie, degree plus one) is achieved by the proposed method, we compared each approximated signal with the analytic one using Eq. 26 at log-linearly spaced sampling rates ranging from to .
The resulting convergence curves are depicted in Fig. 3. One can see that the theoretical approximation orders are accurately validated for all basis functions and that the use of high-order non-interpolating bspline basis functions provides a major advantage. It can also be mentioned that the o-moms of degree three seems to be “over-performing” at low frequencies, namely from to . This low-frequency range represents a “rough” Nyquist-rate range for the excitation waveform considered [Fig. ((b))(b)], especially below .
iii.2 Validation against analytic solutions
We considered two transducer-element shapes and corresponding baffle conditions for which analytic expressions are available for evaluating the sir at any field point: the spherically focused element (ie, spherical cap) with a rigid baffle boundary [2] and the rectangular plane with a soft baffle condition [34]. In all cases, we considered the same excitation waveform [in Fig. S5], characterized by a (center) wavelength for . We evaluated the field signal radiated by the transducer element at three characteristic field points (A, B, and C) positioned relatively to the transducer element center. The three field points were positioned such that their projections onto the surface lie on an axis of symmetry, an edge, and outside the surface projection.
We restricted the comparison to the following basis functions: nearest-neighbor, linear, quadratic keys, bspline of degrees three and five, and o-moms of degree three. Two sampling rates of and were considered, namely a rather low sampling rate for the excitation waveform considered, and a rather high one [Fig. ((b))(b)]. Both shapes considered were represented exactly by nurbs surfaces. For the Gaussian quadrature rule, we relied on the well-known Gauss-Legendre one. As the regularity of the integrand involved in the evaluation of the sir [Eq. 1
] has not been studied in depth, we relied on a heuristic strategy consisting of selecting the number of quadrature points in each
direction of the nurbs surface to obtain a spatial sampling rate equivalent to the sampling rate considered for the time dimension. By doing so, it was observed that the resulting accuracy was not bound by the quadrature rule in the cases studied.Because no analytic expressions could be derived for the field signals, we relied on a very high sampling rate to evaluate the analytic sir, the excitation waveform, and their discrete convolution. We used a reference sampling rate of for the - case and of for the - case to guarantee integer downsampling factors. We also compared the sir obtained using Eq. 25.
iii.2.1 Spherically focused element with a rigid baffle condition
The geometry of the spherically focused transducer element considered is defined by an active diameter and a spherical radius defined such that , namely a similar ratio to the one studied in [2, Fig. 4]. The nurbs representation and corresponding four smooth bezier patches is identical to the one described in Fig. S2. Our heuristic strategy to define the number of quadrature points led to Gauss-Legendre quadrature points per bezier patch for the sampling rate of , and for the sampling rate of . The three field points (A, B, and C) at which field signals were evaluated all lie in the same plane of revolution at a depth of . The lateral coordinate of the first one (A) is . For the second one (B), the lateral coordinate was computed such that its projection onto the surface lies on an edge, resulting in . The lateral coordinate of the last one (C) was simply set to such that its projection onto the surface lies outside.
The relative two-norm errors of the field signals at each field point for both sampling rates and all basis functions considered are reported in Table 1. These results indicate that the basis function of the highest order performs best (ie, bspline5), with a relative error of approximately and at a sampling rate of and , respectively. The field signals and sir for the bspline basis function of degree five is depicted in Fig. 4. Despite the very different sir at the three field points considered, the relative errors are similar at a given sampling rate. One can note that the approximated sir typically contains ripple artifacts because the very high frequencies cannot be accounted for at such low sampling rates. This is especially visible for the field point A at a sampling rate of . Yet, the field signals do not suffer from such artifacts. We can also observe that the approach would tend to an accurate approximation of the sir at much higher rates, should such a quantity be of interest. Field signals obtained with the other basis functions considered can be found in Fig. S11, S10, S9, S8 and S7.
Freq. | Point | Nearest | Linear | Keys | bspline3 | o-moms3 | bspline5 |
origin=c]90 | A | ||||||
B | |||||||
C | |||||||
origin=c]90 | A | ||||||
B | |||||||
C |
iii.2.2 Rectangular element with a soft baffle condition
The geometry of the rectangular transducer element considered is defined by a width of and a height of , chosen to reflect typical width-to-height ratios of transducer elements composing linear arrays. The nurbs representation is simply a bilinear surface, resulting in a single smooth bezier patch onto which the Gauss-Legendre quadrature rule was deployed. Our heuristic strategy to define the number of quadrature points led to Gauss-Legendre quadrature points for the sampling rate of , and for the sampling rate of . The three field points (A, B, and C) at which field signals were evaluated all lie in a plane parallel to the element. They were positioned at a depth of and an elevation of , with lateral coordinates of , , and , respectively. Note that this is quite an extreme case that was selected on purpose as this results in sir with very high frequencies.
The relative two-norm errors of the field signals at each field point for both sampling rates and all basis functions considered are reported in Table 2. As for the spherically focused cased, the higher-order basis function performs best. In general, the order of relative errors are also comparable between the two element shapes. The field signals and sir for the bspline basis function of degree five are depicted in Fig. 5. Similarly to the spherically focused case, ripples due to the high-frequency components of the sir can be observed in the approximations. Field signals obtained with the other basis functions considered can be found in Fig. S16, S15, S14, S13 and S12.
The same experiment was also carried out for a rigid baffle condition. All results are consistent for the two different baffle conditions. The detailed results are reported in supplementary material (Sec. S-III.1.1).
Freq. | Point | Nearest | Linear | Keys | bspline3 | o-moms3 | bspline5 |
origin=c]90 | A | ||||||
B | |||||||
C | |||||||
origin=c]90 | A | ||||||
B | |||||||
C |
Iv Discussion
The results obtained from the two experiments carried out demonstrate that the proposed approach can attain a high accuracy wrt analytic reference signals. The comparison of the relative errors obtained in all cases of the second experiment on realistic transducer element shapes (Tables 2, S-I and 1) with those of the convergence-order study on a random stream of Dirac delta functions (Fig. 3) shows that the latter provides an upper bound on the relative error. As such, the first experiment provides a robust way for selecting appropriate basis functions to guarantee a desired accuracy, assuming that suitable numerical quadrature is deployed.
iv.1 Benefit of high-order B-spline basis functions
In pulse-echo us imaging, it is typically acceptable to consider a sampling rate that enable preserving frequencies with a relative spectrum magnitude of approximately . In the experiments carried out, this corresponds to a sampling rate of [Fig. ((b))(b)]. One can note (Fig. 3) that using a bspline basis function of degree five results in a relative error of approximately , and should therefore not induce additional errors in the simulated signals at that sampling rate. On the other hand, using a nearest-neighbor basis function would require a sampling rate of to achieve the same relative error. This is an important observation and a major advantage of the proposed method as large-scale simulations highly benefit from low sampling rates (because of the many discrete convolutions). Thus, the use of high-order basis functions is of primary interest, even if it requires basis functions of slightly greater supports.
iv.2 Non-uniform rational B-spline representation of surfaces
Our choice of representing radiating surfaces as nurbs surfaces was primarily motivated by the fact that they enable representing quadric surfaces exactly, and that most us transducer elements are at most quadric surfaces. Other representations could be used in the proposed approach (eg, analytic ones), provided that the parametrization is smooth such that Gaussian quadrature can be deployed. The fact that any nurbs surface can be decomposed into smooth bezier patches guarantees that this requirement is always met.
Also, nurbs representations are heavily used in cad software. Hence, the design of new transducer element shapes could be evaluated directly from these nurbs definitions without additional processing (such as meshing), similarly to the concept of iga. The optimization of transducer element shapes could also be performed directly from their nurbs definitions with direct evaluation of the field quantities of interest.
iv.3 Gaussian quadrature rules
We only considered the well-known Gauss-Legendre quadrature rule in the present study. Since there is no realistic scenario in which field signals would need to be computed onto the radiating surface (ie, non-singular integrand), the Gauss-Legendre quadrature rule guarantees high-order accuracies. Even though not reported here, we also evaluated the performance of typical (low-order) quadrature rules such as the midpoint, trapezoidal, and simpson ones. As expected, they all performed much less efficiently than the Gauss-Legendre one. We also conducted some preliminary evaluations using the Gauss-Legendre-Lobatto rule. This quadrature rule may be promising in the context of radiating surfaces as it includes the endpoints of the integration interval (ie, the surface edges), at the cost of a slightly reduced accuracy than the Gauss-Legendre rule. However, we could not yet conclude on which one performs generally best. An in-depth study on the regularity of the integrand should be performed to obtain an optimal selection of the quadrature point number.
iv.4 Comparison with other strategies
It is interesting to note that the proposed approach can be used to interpret previously proposed discretized approaches, as different surface representations, quadrature rules, and basis functions can be deployed. For instance, Piwakowski and Delannoy [30] proposed an approach to compute the sir by discretizing the radiating surface into many ideal points and to average the Dirac delta functions between each time samples. This strategy can be considered in the proposed formulation, namely by using a midpoint quadrature rule together with a nearest-neighbor basis function. This would of course result in much worse accuracy than using a high-order Gaussian quadrature rule together with a high-order basis function. In general, the strategy of averaging sir values was proposed and used in many approaches [2, 30, 9, 20]. Thus, these approaches could also benefit from the formulation in basis functions proposed here, in particular to benefit from higher-order basis functions.
From a pure computational perspective, simulation methods based on the sir are particularly suited for parallel implementations such as gpu-based ones. A critical point when it comes to gpu-based implementations is the complexity of arithmetic operations. Analytic expressions for the sir of specific radiating surfaces involve complex operations such as hyperbolic and inverse hyperbolic operations, with many cases depending on the relative positioning of field points wrt radiating surfaces. Such operations would typically not be ideally suited for gpu-based implementations. The proposed formulation relies on many more operations of much lower arithmetic complexity,^{8}^{8}8The most complicated arithmetic operations are the square root and the cosine in the case of soft baffle conditions. and is as such better suited for gpu-based implementations. So far, we did not perform an exhaustive benchmark of our current implementation against well-known software such as fieldii. Our initial experiments in the context of pulse-echo sa imaging indicated that we could reduce the computing time by approximately two orders magnitude, considering consumer-level cpu (for fieldii) and gpu (for the proposed approach). Note that the initial goal of developing this approach was to enable us generating a sufficient amount of data for the purpose of training cnn-based image reconstruction methods in the context of us imaging [27, 28].
V Conclusion
We proposed a spline-based sir approach for the simulation of field signals radiated by arbitrary shapes embedded in both rigid and soft baffles and excited by bandpass waveforms. This approach consists of representing a transducer surface as a nurbs surface and decomposing it into smooth bezier patches onto which high-order Gaussian quadrature rules can be deployed to approximate the surface integral involved in the computation of the sir. Using high-order bspline bases to express the sir, the basis sir amounts to a sum of shifted-and-weighted basis functions that depend on the positions and weights of the quadrature points. The resulting field signal is then obtained by the convolution of the basis coefficients, derived from the excitation waveform, and the basis sir. The use of nurbs enables accurate representations of complex surfaces, and common transducer shapes can be represented exactly. High-order Gaussian quadrature rules enable using fewer quadrature points to attain a desired accuracy than low-order ones commonly used (eg, midpoint or trapezoidal). High-order bspline basis functions enables using a simulation sampling rate identical to the one required to represent the excitation waveform accurately. The numerical experiments demonstrated that the proposed approach can attain errors as low as the sampling errors of the excitation waveform. The extension to transducer arrays and pulse-echo settings is straightforward. This approach is also well-suited to parallel implementations. An initial gpu implementation enabled us to reduce the computing time by up to two orders of magnitude compared with the well-known fieldii simulator.
Acknowledgments
This work was supported in part by the Swiss National Science Foundation under Grant 205320_175974.
References
- Abramowitz and Stegun [1964] Abramowitz, M., and Stegun, I. A., eds. (1964). number 55 in Applied Mathematics Series Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (US Government Printing Office, Washington, D.C.).
- Arditi [1981] Arditi, M. (1981). “Transient fields of concave annular arrays,” Ultrason. Imaging 3(1), 37–61.
- Bæk et al. [2012] Bæk, D. B., Jensen, J. A., and Willatzen, M. (2012). “Spatial impulse response of a rectangular double curved transducer,” J. Acoust. Soc. Am. 131(4), 2730–2741.
- Blu et al. [2001] Blu, T., Thevenaz, P., and Unser, M. (2001). “Moms: maximal-order interpolation of minimal support,” IEEE Trans. Image Process. 10(7), 1069–1080.
- Cox [1972] Cox, M. G. (1972). “The numerical evaluation of B-splines,” IMA J. Appl. Math. 10(2), 134–149.
- de Boor [1972] de Boor, C. (1972). “On calculating with B-splines,” J. Approx. Theory 6(1), 50–62.
- de Boor [1987] de Boor, C. (1987). Applied Mathematical Sciences A Practical Guide to Splines, 4 ed. (Springer-Verlag, New York).
- Delannoy et al. [1979] Delannoy, B., Lasota, H., Bruneel, C., Torguet, R., and Bridoux, E. (1979). “The infinite planar baffles problem in acoustic radiation and its experimental verification,” J. Appl. Phys. 50(8), 5189–5195.
- D’hooge et al. [1997] D’hooge, J., Nuyts, J., Bijnens, B., De Man, B., Suetens, P., Thoen, J., Herregods, M.-C., and Van de Werf, F. (1997). “The calculation of the transient near and far field of a baffled piston using low sampling frequencies,” J. Acoust. Soc. Am. 102(1), 78–86.
- Foster et al. [1983] Foster, D. R., Arditi, M., Foster, F. S., Patterson, M. S., and Hunt, J. W. (1983). “Computer simulations of speckle in B-scan images,” Ultrason. Imaging 5(4), 308–330.
- Frangi et al. [2018] Frangi, A. F., Tsaftaris, S. A., and Prince, J. L. (2018). “Simulation and synthesis in medical imaging,” IEEE Trans. Med. Imaging 37(3), 673–679.
- Greenspan et al. [2016] Greenspan, H., van Ginneken, B., and Summers, R. M. (2016). “Guest editorial deep learning in medical imaging: Overview and future promise of an exciting new technique,” IEEE Trans. Med. Imaging 35(5), 1153–1159.
- Harris [1981a] Harris, G. R. (1981a). “Review of transient field theory for a baffled planar piston,” J. Acoust. Soc. Am. 70(1), 10–20.
- Harris [1981b] Harris, G. R. (1981b). “Transient field of a baffled planar piston having an arbitrary vibration amplitude distribution,” J. Acoust. Soc. Am. 70(1), 186–204.
- Hunt et al. [1983] Hunt, J. W., Arditi, M., and Foster, F. S. (1983). “Ultrasound transducers for pulse-echo medical imaging,” IEEE Trans. Biomed. Eng. BME-30(8), 453–481.
- Jensen et al. [2016] Jensen, J., Stuart, M. B., and Jensen, J. A. (2016). “Optimized plane wave imaging for fast and high-quality ultrasound imaging,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 63(11), 1922–1934.
- Jensen [1991] Jensen, J. A. (1991). “A model for the propagation and scattering of ultrasound in tissue,” J. Acoust. Soc. Am. 89(1), 182–190.
- Jensen [1996a] Jensen, J. A. (1996a). “Field: A program for simulating ultrasound systems,” in 10th Nord. Conf. Biomed. Imaging, Vol. 4, pp. 351–353.
- Jensen [1996b] Jensen, J. A. (1996b). “Ultrasound fields from triangular apertures,” J. Acoust. Soc. Am. 100(4), 2049–2056.
- Jensen [2001] Jensen, J. A. (2001). “Speed-accuracy trade-offs in computing spatial impulse responses for simulating medical ultrasound imaging,” J. Comput. Acoust. 09(03), 731–744.
- Jensen and Svendsen [1992] Jensen, J. A., and Svendsen, N. B. (1992). “Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 39(2), 262–267.
- Keys [1981] Keys, R. (1981). “Cubic convolution interpolation for digital image processing,” IEEE Trans. Acoust. 29(6), 1153–1160.
- Kino [1987] Kino, G. S. (1987). Acoustic Waves: Devices, Imaging, and Analog Signal Processing (Prentice-Hall, Englewood Cliffs, NJ).
- Lasota et al. [1984] Lasota, H., Salamon, R., and Delannoy, B. (1984). “Acoustic diffraction analysis by the impulse response method: A line impulse response approach,” J. Acoust. Soc. Am. 76(1), 280–290.
- Lewiss et al. [2014] Lewiss, R. E., Hoffmann, B., Beaulieu, Y., and Phelan, M. B. (2014). “Point-of-care ultrasound education,” J. Ultrasound Med. 33(1), 27–32.
- Lockwood and Willette [1973] Lockwood, J. C., and Willette, J. G. (1973). “High‐speed method for computing the exact solution for the pressure variations in the nearfield of a baffled piston,” J. Acoust. Soc. Am. 53(3), 735–741.
- Perdios et al. [2020] Perdios, D., Vonlanthen, M., Martinez, F., Arditi, M., and Thiran, J.-P. (2020). “CNN-based image reconstruction method for ultrafast ultrasound imaging,” https://arxiv.org/abs/2008.12750.
- Perdios et al. [2021] Perdios, D., Vonlanthen, M., Martinez, F., Arditi, M., and Thiran, J.-P. (2021). “CNN-based ultrasound image reconstruction for ultrafast displacement tracking,” IEEE Trans. Med. Imag. 40(3), 1078–1089.
- Piegl and Tiller [1997] Piegl, L., and Tiller, W. (1997). Monographs in Visual Communication The NURBS Book, 2 ed. (Springer Berlin Heidelberg, Berlin, Heidelberg).
- Piwakowski and Delannoy [1989] Piwakowski, B., and Delannoy, B. (1989). “Method for computing spatial pulse response: Time‐domain approach,” J. Acoust. Soc. Am. 86(6), 2422–2432.
- Piwakowski and Sbai [1999] Piwakowski, B., and Sbai, K. (1999). “A new approach to calculate the field radiated from arbitrarily structured transducer arrays,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 46(2), 422–440.
- Rindal et al. [2019] Rindal, O. M. H., Austeng, A., Fatemi, A., and Rodriguez-Molares, A. (2019). “The effect of dynamic range alterations in the estimation of contrast,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 66(7), 1198–1208.
- Rodriguez-Molares et al. [2020] Rodriguez-Molares, A., Rindal, O. M. H., D’hooge, J., Masoy, S.-E., Austeng, A., Lediju Bell, M. A., and Torp, H. (2020). “The generalized contrast-to-noise ratio: A formal definition for lesion detectability,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 67(4), 745–759.
- San Emeterio and Ullate [1992] San Emeterio, J. L., and Ullate, L. G. (1992). “Diffraction impulse response of rectangular transducers,” J. Acoust. Soc. Am. 92(2), 651–662.
- Schoenberg [1946a] Schoenberg, I. J. (1946a). “Contributions to the problem of approximation of equidistant data by analytic functions. part a. on the problem of smoothing or graduation. a first class of analytic approximation formulae,” Q. Appl. Math. 4(1), 45–99.
- Schoenberg [1946b] Schoenberg, I. J. (1946b). “Contributions to the problem of approximation of equidistant data by analytic functions. part b. on the problem of osculatory interpolation. a second class of analytic approximation formulae,” Q. Appl. Math. 4(2), 112–141.
- Stepanishen [1971a] Stepanishen, P. R. (1971a). “The time‐dependent force and radiation impedance on a piston in a rigid infinite planar baffle,” J. Acoust. Soc. Am. 49(3B), 841–849.
- Stepanishen [1971b] Stepanishen, P. R. (1971b). “Transient radiation from pistons in an infinite planar baffle,” J. Acoust. Soc. Am. 49(5B), 1629–1638.
- Stepanishen [1981] Stepanishen, P. R. (1981). “Acoustic transients from planar axisymmetric vibrators using the impulse response approach,” J. Acoust. Soc. Am. 70(4), 1176–1181.
- Theumann et al. [1990] Theumann, J., Arditi, M., Meister, J., and Jaques, E. (1990). “Acoustic fields of concave cylindrical transducers,” J. Acoust. Soc. Am. 88(2), 1160–1169.
- Thevenaz et al. [2000] Thevenaz, P., Blu, T., and Unser, M. (2000). “Interpolation revisited,” IEEE Trans. Med. Imaging 19(7), 739–758.
- Tjøtta and Tjøtta [1982] Tjøtta, J. N., and Tjøtta, S. (1982). “Nearfield and farfield of pulsed acoustic radiators,” J. Acoust. Soc. Am. 71(4), 824–834.
- Tupholme [1969] Tupholme, G. E. (1969). “Generation of acoustic pulses by baffled plane pistons,” Mathematika 16(2), 209–224.
- Unser [1999] Unser, M. (1999). “Splines: a perfect fit for signal and image processing,” IEEE Signal Process. Mag. 16(6), 22–38.
- Unser et al. [1991] Unser, M., Aldroubi, A., and Eden, M. (1991). “Fast b-spline transforms for continuous image representation and interpolation,” IEEE Trans. Pattern Anal. Mach. Intell. 13(3), 277–285.
- Unser et al. [1993a] Unser, M., Aldroubi, A., and Eden, M. (1993a). “B-spline signal processing: Part I—Theory,” IEEE Trans. Signal Process. 41(2), 821–833.
- Unser et al. [1993b] Unser, M., Aldroubi, A., and Eden, M. (1993b). “B-spline signal processing: Part II—Efficiency design and applications,” IEEE Trans. Signal Process. 41(2), 834–848.
- van Sloun et al. [2020] van Sloun, R. J. G., Cohen, R., and Eldar, Y. C. (2020). “Deep learning in ultrasound imaging,” Proc. IEEE 108(1), 11–29.
- Verhoef et al. [1984] Verhoef, W. A., Cloostermans, M. J. T. M., and Thijssen, J. M. (1984). “The impulse response of a focused source with an arbitrary axisymmetric surface velocity distribution,” J. Acoust. Soc. Am. 75(6), 1716–1721.
- Wagner et al. [1983] Wagner, R. F., Smith, S. W., Sandrik, J. M., and Lopez, H. (1983). “Statistics of speckle in ultrasound B-scans,” IEEE Trans. Sonics Ultrason. 30(3), 156–163.
- Wang [2016] Wang, G. (2016). “A perspective on deep imaging,” IEEE Access 4, 8914–8924.
- Wang et al. [2018] Wang, G., Ye, J. C., Mueller, K., and Fessler, J. A. (2018). “Image reconstruction is a new frontier of machine learning,” IEEE Trans. Med. Imaging 37(6), 1289–1296.
- Ziv et al. [2006] Ziv, A., Wolpe, P. R., Small, S. D., and Glick, S. (2006). “Simulation-based medical education: An ethical imperative,” Simul. Healthc. J. Soc. Simul. Healthc. 1(4), 252–256.
S-I Introduction
A typical example of a sir and corresponding frequency spectrum is shown in Fig. S1 for a rectangular element designed to work at a driving frequency of , from which it is clear that a sampling rate of a few is required.
S-Ii Proposed Approach
S-ii.1 Non-uniform rational B-spline surface representations
As most transducer elements composing conventional us transducers are (at most) quadric surfaces, the assumption that the radiating surface can be represented exactly by a nurbs surface is generally valid in us imaging since nurbs surfaces can represent quadric surfaces exactly [29, Chap. 8]. For instance, transducer elements composing linear and phased arrays are cylindrical shells (ie, elliptic cylinders), those composing convex arrays are toroidal shells (ie, hyperbolic paraboloids), and those composing 2-D matrix arrays are simple rectangular planes. Spherically focused transducer elements are spherical caps. An example of the latter is shown in Fig. S2, for which the nurbs surface definition can be obtained by revolving a circular arc (ie, represented by nurbs curve). The nurbs surface can then be decomposed into four rational biquadratic bezier smooth patches. Each of these smooth patches is defined by Eq. 12 onto the characteristic unit square space
. As they are smooth surfaces, Gaussian quadrature rules can be deployed onto the 2-D parametric space (also referred to as parent element) by means of the tensor product. Resulting quadrature coordinates can then be mapped onto the radiating surface exactly using
Eq. 12. This process is also illustrated in Fig. S2, in which the Gauss-Legendre quadrature rule was considered to obtain the quadrature coordinates of a grid. Note that Gaussian quadrature rules are typically defined on a characteristic interval [1, Sec. 25.4], but can be mapped to arbitrary intervals.Thanks to their small (or nonexistent) curvatures,all other aforementioned transducer element shapes can be represented by nurbs surfaces that are themselves single rational bezier patches already, and thus do not require any further decomposition. An example of such as case for the cylindrical shell is shown in Fig. S3.
S-ii.2 Complete process and implementation details
An illustration of the convolution involved in the evaluation of Eq. 18 for each quadrature point is depicted in Fig. S4 for four different basis functions, namely nearest-neighbor, linear, quadratic keys, and cubic bspline. As the first three basis functions are interpolating, their corresponding (pre-filtered) coefficients are equal to the excitation waveform samples. The coefficients of the cubic bspline have larger amplitudes than the excitation waveform to compensate for its non-interpolating property. One can already note the major issue associated with near-neighbor and linear interpolation, namely a shift error and an underestimated response, respectively. The advantage of using a cubic bspline over a quadratic Keys, both having a support of four samples, can already be noticed. This will be confirmed by the numerical validation of the expected convergence orders (Sec. III.1).
S-ii.3 Extension to arrays
The extension to us transducers composed of multiple transducer elements, typically arranged as arrays, is straightforward. Let us consider a generic definition of a transducer array composed of a set of transducer elements. Each element can be represented by a nurbs surface and decomposed into a set of smooth bezier patches onto which a Gaussian quadrature rule can be defined. These elements are generally of the same shape, which would result in similar nurbs representations, although this is not required. Each -th element can be driven by a different excitation waveform , from which the corresponding basis coefficients needs only to be evaluated once [Eq. 20]. Using the principle of superposition (linear acoustics) and Eq. 24, the field signal at any field point radiated by such an array can be expressed as
(S1) |
where each element basis sir can be computed using Eq. 23. It is common to define explicitly the delays and apodization weights applied across the array elements, rather than implicitly including them in the different excitation waveforms. Such delays and apodization weights are typically used to shape the transmit beam (beamforming), for instance to focus at a desired position or to steer an unfocused wavefront. Thus, and without loss of generality, Eq. S1 can be rewritten as
(S2) |
In pulse-echo imaging, it is conventionally assumed that all (identical) transducer elements composing the array have the same electromechanical impulse response. An identical electric excitation is also typically used, such that the characteristic excitation waveform is identical. Only the delays and apodization weights applied to the elements may therefore differ such that Eq. S2 can be simplified as
(S3) |
Finally, the echo signal scattered by an ideal reflector at can be easily obtained by the convolution of the field signal of the array [Eq. S3] with the field signal of the receive element [Eq. 24], multiplied by the scattering amplitude [17]. Note that the excitation waveform on transmit is generally different from the one on receive, because the first contains both the electric excitation and the electromechanical impulse response, whereas the second contains the electromechanical impulse response only.
S-Iii Experiments and Results
The analytic expression of the time derivative of a log-normal-modulated sinusoidal rf pulse considered for the experiments (Secs. III.2 and III.1) is shown in Fig. S5. It is a fairly good model for the electromechanical impulse response of transducer elements. Note that the time derivative guarantees a zero dc component [Fig. ((b))(b)], a physical property of such electromechanical impulse response.
S-iii.1 Validation against analytic solutions
S-iii.1.1 Rectangular element with a rigid baffle condition
The same experiment as described in Sec. III.2.2 was carried out for the rectangular transducer element with a rigid baffle condition. The relative two-norm errors of the field signals at each field point for both sampling rates and all basis functions considered are reported in Table S-I. All results are consistent for the two different baffle conditions, even at field point C where the sir of the rigid and soft baffle conditions differ greatly. The field signals and sir for the bspline basis function of degree five are depicted in Fig. S6. By comparing Fig. 5 and Fig. S6, one can see that the effect of the soft baffle condition is visible for all field points, and in particular at field point C. Field signals obtained with the other basis functions considered can be found in Fig. S21, S20, S19, S18 and S17.
Freq. | Point | Nearest | Linear | Keys | bspline3 | o-moms3 | bspline5 |
origin=c]90 | A | ||||||
B | |||||||
C | |||||||
origin=c]90 | A | ||||||
B | |||||||
C |