Medical Ultrasound Imaging is the most widespread real-time non-invasive imaging system: it exploits the ability of human tissue to reflect ultrasound signals In particular, different ultrasound modalities process the reflected echoes, generating either morphological images or measures of tissue properties. The best known biomedical ultrasound modality, called B-mode, produces the standard black and white images by means of several transmission and reception phases, resulting in a map of acoustic tissue impedance. During the last decades the machine performances, image quality and computational power have increased generating new strategies, techniques and methods. In particular, several quantitative modalities have been developed, allowing to measure physical parameters of tissues with diagnostic significance (Doppler modalities [CW] and elastography [ARFI3] can be cited among them). The performance of each modality is deeply dependent on how the transmission and reception phases are tuned reflecting both on the image quality and the accuracy of measurements. The two phases are characterized by some parameters such as pulse shape, central frequency, transmission focus, transmit and receive apodization. Concerning the transmission phase, the effect of the parameters is summarized by the transmit beam pattern which encodes the information about energy distribution across the image field. While the large majority of ultrasound modalities deals with wideband beam patterns, two specific quantitative applications, Continuous Wave (CW) Doppler and Acoustic Radiation Force Impulse (ARFI) elastography employ narrowband beam patterns. The CW modality consists in a continuous emission of a narrow band pulse from some probe channels while echoes are acquired from a region of interest by using a different set of channels [CW1]. This technique can produce a time frequency representation of blood flow and the measures accuracy is of paramount important [CW2]
. The ARFI is a specific elastography technique in which the probe emits a narrow band shock pulse to generate a shear wave in the tissue; then several tracking pulses are emitted to sample the propagation of the shear wave estimating the local shear wave velocity; finally, such an estimation allows to compute the elasticity map of the tissue which is displayed as a color-coded image[ARFI3, ARFI, ARFI1, ARFI2, ARFI4, ARFI5]. Generally, the ideal transmit beam pattern, both narrowband or wideband, should have a main lobe width as uniform as possible along the depth of interest and a side lobe level as low as possible. Moreover the main lobe intensity along the depth should be as uniform as possible, since its maximum peak is bounded by both mechanical and thermal limits assuring the safety of the insonification. These requirements are in contrast with the typical hourglass-shaped beam obtained with a fixed focalization depth [Jen1, Harris1], in which the main lobe width has his minimum nearby the focal depth and enlarges before and above. Furthermore, the near field region below the focal depth is often characterized by side lobes whose intensity is even higher than the value along the central line. Finally, the combined effect of geometric and thermoviscous attenuation causes a quick energy drop-out immediately after the focal depth. These drawbacks have a negative impact on both image quality and measurement accuracy. For example, in ARFI elastography a non uniform main lobe width generates shear waves of non uniform temporal length along the depth and this causes bias in the estimation of shear wave times of arrival. Moreover, the presence of side lobes with high intensity generates spurious shear waves that hinder the effect of the desired one. Finally the quick energy drop-out after the focal depth causes a rapid decrease of shear wave energy along the depth, thus dramatically decreasing the SNR and the accuracy of the estimated shear wave time of arrival [ARFIbias]. Similar considerations hold also for CW [CW2]. The tuning of the high level parameters of a standard beam pattern, for instance focal depth, aperture, transmit frequency, may improve some beam pattern features but often at the expense of others. For example, increasing the aperture and moving the focal depth forward reduces the energy drop out but decreases the main lobe width uniformity. A different technique, named synthetic transmit beam [MO] improves the transmit beam uniformity retrospectively, i.e. recombining different laterally spaced received signals. Though quite effective, this technique is not feasible for CW and ARFI methods in which only one narrowband transmission is performed. A different approach addressed in the literature [ku, cu, cu2, tr, Opt, Opt1, Opt2]
, consists in optimizing a set of low level parameters, such as the transmit intensity at each channel, i.e. the apodization window, by minimizing a suitable cost function in the wideband case. In this framework the solution of this problem can admit a closed form or may be easily computed as the cost function is convex. To the best of our knowledge the optimization of transmit delays has received much less attention. In this work we propose to optimize the narrowband beam pattern, suited for CW Doppler and ARFI elastography, considering each single transmit delay as an independent variable and keeping fixed the transmit intensity. The choice of optimizing apodization-free beam patterns, in which all the active elements transmit at the same voltage, enlarges the scope of application of this method also to low-price commercial apparati, equipped with simple two level transmit pulsers instead of linear amplifiers. The much higher number of degrees of freedom, in comparison with the classical single-focus delay parametrization, allows the energy to be distributed over a greater depth range, thus improving the beam uniformity along the depth. However, the optimization problem turns out to be non-convex, thus requiring specific strategies to avoid local minima. To compare quantitatively the optimized beam patterns with a large set of standard ones, we make use of quality metrics. Existing beam pattern metrics, like peak side lobe level, integrated side lobe level or full width half maximum of main lobe are quite effective for far field beam patterns or to evaluate near field beam pattern at a single depth. However, they fail to capture the global beam pattern quality along a depth range. For this reason, we introduced a novel set of metrics specifically designed to measure how much a beam pattern approaches the ideal behaviour required by the applications described above. The paper is organized as follows: section 2 briefly describes the optimization problem, in section 3 we define the simulation setup, we present the results obtained with a non-gradient based optimization algorithm and we propose some metrics to evaluate the results. Finally, in section 4 some conclusions are drawn.
2. The Optimization Problem
We proceed in describing the simulation of the transmit (TX) beam pattern (BP) and the corresponding optimization problem. We focus on the linear probe case, as for convex probes the treatment is analogous. In particular, we consider a probe composed by piezoelectric elements of width
that convert the electric pulse into a pressure wave and vice-versa. We consider a Cartesian coordinate system whose origin corresponds to the geometrical center of the probe (see Figure1), so that the center of each probe element is in the position with . Each point of the image plane is identified by a positive value of the coordinate.
2.1. The BP equation
The transmit BP is the pattern of radiation generated from a beam and it can be described as the energy or the power, for a periodic signal, of the propagating wavefront, evaluated at every spatial point of the field. The propagating wavefront can be computed from the signal emitted by the probe and the spatial-temporal impulse response characterising the probe itself. In particular, the impulse response function of each -th element is the integral over the element surface of a spherical wave starting from each point of the element surface [Harris1]
where is the ultrasound speed propagation. We assume that every active element emits the same temporal TX waveform . Thus, the resulting signal is obtained in every point by convolving with the sum of all the impulse responses of the active elements and with a term accounting for the tissue thermoviscous attenuation . The waveform can be trasmitted with a different delay for each active element. The BP shape is dependent from the set of transmission times thanks to the various destructive and constructive interferences which are generated during the propagation. We introduce the transmission delay , associated with the -th element, and we describe the delay in the signal emission as a translation in time . By defining
we can describe the propagating wavefront in time as:
In the narrow band case we can take as an infinite pulse with given frequency , i.e. . Thus, the power of the BP is
where . Moving to the time frequency domain and relying on Plancherel’s theorem, stating the equivalence of signal power in time and frequency domain, the BP power can be written as the sum of the squared absolute values of the Fourier series coefficients:
where is the Fourier Series of . Since, we have assumed to be a pure tone, only one coefficient of the series is different from zero. Therefore, the equation 2 reduces to:
where we denote . In particular, we can exploit the modulus properties obtaining:
Finally, making use of the trigonometric form for the exponential and using the cosine addition formula, the BP takes the form
is obtained as the product between the Fourier Series of
and the conjugate of the Fourier transform ofat frequency .
2.2. Least Squares problem formulation
The aim of the optimization problem is to make as similar as possible to a prescribed BP shape by varying the frequency , the number of active elements and their corresponding delays. We are interested in BP shapes symmetric with respect to the probe central transmit line, so that we impose . Moreover, we want to consider the delays as free variables so that the BP changes according to parameters:
In order to compute the BP in equation 3 we sampled by taking a grid of points with where and (see Fig.1). From a computational point of view we approximated , the integral over a single element surface, according with [Jen1, Jen2]. Moreover, by fixing the central frequency , we have pre-computed the attenuation coefficients and the impulse response function for each element. An example of impulse response function for a linear probe with pulse central frequency MHz is given in the left panel of Figure 1. Then we compute by using the discrete Fourier transform over time at the frequency
. The sampled BP is a vectordepending non-linearly on , whose components are given by where is an unrolled index, namely . Therefore, we can formulate a Least Squares (LS) problem as follows
where is the discretization of the prescribed BP shape we want to approximate. Since the vector is separately periodic w.r.t. the delays , the optimization problem is non linear and non convex as the parameter domain contains an
dimensional torus which makes it not suitable for being solved by means of standard optimization techniques. Among the many methods developed to deal with non convex optimization, e.g. genetic algorithms, simulated annealing, manifold optimization to mention a few, we propose to use the particle swarm optimization. This method has been proved to be particularly efficient in solving this problem, after a number of attempts with different methods.
3. Simulation Setting and Results
Let us consider a linear probe with a depth range of cm. Our purpose is to obtain a beam shape as uniform as possible along depth, with the lowest possible level of sidelobes. The goal is to depart from the usual hour-glass shaped beam pattern related to a fixed focal depth, by focalizing each element at a different depth to achieve a more stretched shape. We define the desired BP as a rectangular shape centered at depth cm with length of cm, with value one inside and zero outside the rectangle (right panel of Fig. 1). This shape aims to describe the uniformity along depth disregarding the behavior before the first depth of interest. Moreover, the zero value at the sides of the rectangle enforces a low level of side lobes. We decided to optimize the frequency of narrowband pulse and the aperture with a greedy approach and the vector of delays with Particle Swarm Optimization Algorithm [PSO], with the following rationale. By choosing a finite set of frequencies suitable for the probe, we can generate the corresponding impulse response map set and rapidly simulate a beam. The aperture is an integer and is reasonable to assume the active surface length is greater than the smallest dimension of the rectangle. The choice of PSO algorithm is justified by the nature of the problem: its domain contains a N-torus so every gradient-based algorithm would work in a fixed local chart, thus being very sensitive to the initialization. On the contrary, PSO adopts multiple random initialization hence being less affected by the manifold domain structure. The simulator used to generate the beam images is developed in Python™. We have validated it by a comparison with the well-known ultrasound simulator FIELD II [Jen2]. Finally, we have used PySwarms 444https://pyswarms.readthedocs.io/en/latest/ free software research toolkit to apply PSO algorithm in Python™. At each iteration the algorithm generates a BP for each particle corresponding to a vector of TX delays and measures the adherence between the actual and prescribed beam, than it updates each particle taking in account its adherence value and the overall direction of the swarm.
3.1. The standard simulated beam pattern
First, we analyze a set of narrow-band BP obtained with a standard focalization law as a function of focal depth and aperture (as displayed in Fig. 2), at MHz, equal to the frequency found by the optimization algorithm for a fair comparison. Looking at the panel, the well-known drawbacks of standard delay profile can be observed. For small aperture the main lobe is quite large and its intensity profile over depth rapidly decreases. Increasing the aperture, the main lobe width shrinks around the focal depth but it enlarges considerably in the near field region, where the intensity peak is even not aligned with the main lobe axis. These drawbacks can severely affect the accuracy of Doppler frequency shift measurements in CW and shear wave propagation estimation in ARFI elastography [ARFIbias].
3.2. Optimized beam pattern
In Fig. 3 we display some examples of the BP obtained optimizing the aperture, frequency and delay profile. Focusing the analysis on the depth range of interest given by the rectangle of the prescribed beam pattern (right panel in Figure 1), it can be seen that the intensity peak is located on the main lobe axis for every depth. Moreover, the main lobe width is considerably uniform over all the depth range especially looking at the near field in comparison with the BP obtained with the same aperture ( elements) and a standard delay profile.
Furthermore an improvement is achieved also in terms of intensity uniformity along the BP main lobe axis (Fig. 4).
These improvements are achieved thanks to the optimized delay profile displayed in Fig. 5, together with the standard profiles for three different focal depths. The departure of the optimized profile from the family of standard profiles, confirms the effectiveness of the proposed approach in optimizing the narrow band TX beam pattern.
Considering each delay as a free variable has increased significantly the number of delays combination therefore augmenting the possible BP shapes. In particular, the optimized profile can be seen as a combination of standard delay profiles at different focal depths.
3.3. Evaluation Metrics
In order to allow a quantitative evaluation of optimized BPs with respect to standard ones, we devised a set of metrics aimed at providing a synthetic evaluation of the overall BP quality. The first metric is aimed at measuring the main lobe width average value and uniformity along the depth. To do so, we have made the assumption that the Main Lobe (ML) at each depth consists in the BP portion centered on the transmission line that differs from the central value less than dB (Figure 6). The choice of this threshold is coherent with the usual definition of Full Width at Half Maximum. Moreover, this assumption penalizes those patterns whose maximum value is not on the central line.
Therefore, we define Main Lobe Width at a certain fixed depth as:
where is the azimuth coordinate of the beam central line and .
The second metric evaluates the overall side lobe level. Let us define the Side Lobe Level at a fixed depth as the average over of all the BP values outside the region of the main lobe as defined in above.
where denotes the measure of the complementary set of . Notice that in the definition the is normalized with respect to the central line value at the each depth
, in order to cope with the overall decrease trend of the BP along the depth. By computing mean and standard deviation ofand over a set of depths of interest, we can evaluate the main lobe width uniformity and narrowness and the side lobes uniformity and average level.
Finally, in order to evaluate the drop out of main lobe power along depth we compute the mean integral of the central line values.
where denotes a range of depths of interest. This quantity measures how much the line is uniform along the depth and gives information about the power drop: it results in a summary index to have a quantitative evaluation corresponding to Fig. 4.
3.3.1. Quantitative Analysis
|Depths of Interest||Depths|
|BP||MLW (mm)||CLP||SLL (dB)||MLW (mm)||CLP||SLL (dB)|
|f (cm)||A. E.||Mean||STD||Mean||STD||Mean||STD||Mean||STD|
In Table 1 we reported the measurements for all the standard BPs displayed in Fig. 2 and for the optimized ones displayed in Fig. 3. The first block, highlighted in blue, takes into account the depths in corresponding to the rectangle in the prescribed BP as displayed in right panel of Fig. 1. The second block, highlighted in red, takes into account all the depths . This second block allows to quantify the BP behavior also beyond the optimization region and to evaluate the sensitivity of the optimization with respect to the choice of the prescribed BP. All the three optimized BPs have a thinner and more uniform MLW with respect to all the 16 evaluated standard BPs. Moreover, two out of three optimized BP have a higher than all the standard BPs, while the one obtained by activating 22 elements has a comparable with the best cases in standard BPs set. More significantly, the best values of for standard BP are in general achieved with large apertures, but this condition is associated with worse values of . Looking at the we can observe that the optimized BPs achieve values in the middle between the best and the worst values of standard BP, both for mean and standard deviation. However, the lowest average SLL for standard BP is just 1.5 dB lower than the lowest average for optimized BP and is achieved at the price of a much higher . Finally, the above considerations are still valid for the results in the extended range of depths (red block), with a further improvement of the optimized BP (22 elements) with respect to standard ones in terms of standard deviation of . These results support the thesis that the proposed optimization method enforces the prescribed beam pattern behaviour also outside the optimization zone, thus making the choice of prescribed BP not critical.
Our work provides an effective approach in optimizing the narrow band transmit beam pattern, therefore providing a method to improve accuracy and reduce artifacts in important medical ultrasound applications such as CW Doppler and ARFI elastography imaging. The adopted strategy makes it easy to extend the calculations to convex probes, by a simple change from Cartesian to polar coordinates, as well as to steered BP. Future developments will include the apodization weights in the set of variables to be optimized, compensating in this way the slight worsening of side lobe level observed with the present method with respect to standard delay laws [Sz]. Moreover, the extension to wide band beam patterns will be investigated as well.