# Optimized Virtual Sources Distributions for 3-D Ultrafast Diverging Wave Compounding Imaging: A Simulation Study

Goulven Le Moign, Patrice Masson, Olivier Basset, Hervé Liebgott and Nicolas Quaegebeur

**Abstract**— Ultrafast ultrasound imaging allows observing rapid phenomena; combined with 3-D imaging it has the potential to provide more accurate analysis on organs which leads, at the end, to better diagnosis. Coherent compounding using diverging waves is commonly used to reconstruct high quality images on large volumes while keeping frame rate high enough to allow dynamic analysis. In practice, the virtual sources that drive the diverging waves are often distributed by a deterministic way: following a regular grid, concentric rings and spirals. Even though those deterministic distributions can offer various trade-off in terms of imaging performance, other distributions can be considered to improve imaging performance. It is herein suggested to look at alternative virtual sources distributions for optimizing the lateral resolution and the secondary lobes level on several PSFs by means of a multiobjective genetic algorithm. The optimization framework leads to seven pseudo-irregular distributions of virtual source distributions that have not yet been found in the literature. An analysis on the imaging performance with a simulated phantom shows that these new distributions offer different trade-offs between lateral resolution and contrast, respectively measured on point-like reflectors and anechoic cysts. As an example, one of these optimized distributions improves the lateral resolution by 16% and gives equivalent contrast values on cysts, when compared to a concentric-rings-based distribution.

**Index Terms**— 3-D ultrafast ultrasound imaging, 2-D ultrasound probe, diverging wave compounding, virtual sources, multiobjective genetic algorithm

## I. INTRODUCTION

In medical imaging, echography has become a standard tool in many specialties like cardiology [1], vascular medicine [2] or kidney investigation [3]. Echography is evolving and new applications, still at the research stage, have recently emerged like functional brain imaging [4] or vector flow mapping [5]. These advances are mainly enabled by the use of innovative broad insonifications, that can achieve frame rates of several kHz depending on the depth of the investigation [6]. Although these recent advances have allowed 2-D fast imaging, exploring the third spatial dimension is still a challenging task for fast ultrasound imaging. Probe technologies have evolved to perform 3-D imaging by using 2-D arrays that can explore a whole volume without having to move the probe. Such arrays typically contain thousands of elements,

This work was supported by the Region Auvergne Rhone Alpes within the program SCUSI in France, by the LABEX CELYA (ANR-10-LABX-0060) of Université de Lyon in France, by the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon in France and by the Natural Sciences and Engineering Research Council of Canada (NSERC).

G. Le Moign was affiliated with GAUS, Dept. Mechanical Engineering, Université de Sherbrooke, Sherbrooke, QC, Canada J1K 2R1 and with CREATIS, Univ. Lyon, INSA Lyon, UCBL1, CNRS UMR 5220, Inserm U1206, UJM Saint-Etienne. 69100, Villeurbanne, France. (e-mail: goulven.le.moign@gmail.com)

P. Masson and N. Quaegebeur are affiliated with GAUS, Dept. Mechanical Engineering, Université de Sherbrooke, Sherbrooke, QC, Canada J1K 2R1.

O. Basset and H. Liebgott are affiliated with CREATIS, Univ. Lyon, INSA Lyon, UCBL1, CNRS UMR 5220, Inserm U1206, UJM Saint- Etienne. 69100, Villeurbanne, France.

which requires numerous channels to be addressed. Hence, implementing 2-D arrays in 3-D fast ultrasound imaging has been suggested by limiting the number of channels, by using either sparse array [7], [8], Row Column Addressing (RCA) [9], [10] or μ-beamforming strategies [11], [12]. Another aspect to consider in such implementation is the way the probe is driven to explore a given volume. Multiline transmit [13], fan-beam [14] and broad insonifications by plane waves or diverging waves using Coherent Compounding [15], [16] also, appear as good candidates to do so. Moreover, recent works allow performing Coherent Compounding with RCA and sparse probe technologies [17], [18]. Several studies have tested 3-D plane or diverging wave sequences [16], [18] and they have shown, like in 2-D imaging, that the choice of the number of virtual sources relies on a compromise between imaging volume rate and quality. Improving the image quality by adding more virtual sources is limited at a certain point, as a plateau becomes evident beyond a certain number of virtual sources [16], [19].

Even though the number of virtual sources is fixed accord- ing to a desired imaging volume rate, they can be distributed in many ways. Regular-grid-distributions are often used as a generic approach. It has been reported in [20] and [21] for 2-D ultrasound imaging that the spatial periodicity of the virtual sources can generate grating lobes. Thus, a trade-off has to be made between the lateral resolution and the grating lobes. Additionally, employing aperiodic ring-based and spiral-based plane wave distributions enhances contrast in simulations [22]
and experiments [23], while maintaining frame rate and spatial resolution.

In this paper, a generic framework is proposed for searching and assessing other geometries of virtual sources distributions that can offer competitive imaging performance. This framework is inspired by previous studies on sparse array distributions [24], radar [25] and acoustic source localization [26] that have been suggested in the literature; these studies generally aim at improving the spatial resolution and contrast of an imaging system by reducing both the main beam width and the secondary lobes of an array. Several algorithms have been suggested in the literature to handle such optimization, like Simulated Annealing [27], Genetic Algorithm [28] or Particle Swarm Optimization [29]. The present study also focuses on optimizing the lateral resolution and the side lobes level. A first optimization approach would be to define a single objective function by aggregating both metrics into a single one. This requires to choose a *priori* weighting coefficients to set up the priority between the two metrics. However, these coefficients need to be set carefully to avoid a quality metric to be overwhelmed by the other one. An optimization based on Non-Dominated Solutions (NDS) [30] is suggested to bypass this problem, and also because such approach allows to find a set of solutions with a single optimization. Some evolutionary algorithms [30], [31] have been adapted to handle such multiobjective approach: it is the case for multiobjective particle swarm optimization (MOPSO) [32] and non-dominated sorting genetic algorithm (NSGA-II, or MOGA for MultiObjective Genetic Algorithm) [33]. A MOPSO [34] have been tried, as well as the MOGA from the Global Optimization Toolbox in Matlab [35]. The MOGA has been chosen, following better preliminary results. In general, the side area of a reconstructed image in ultrasound imaging has lower levels for the quality metrics. In order to better distribute the levels of the quality metrics outside the center of the probe, a special attention is given to the homogeneity of the quality metrics. The global framework, objective functions and multiobjective optimization algorithm are detailed in Section II. The optimized distribution that have been found are presented in III-A. The imaging performance when using optimized distributions are compared to nine deterministic distributions in Section III-B.

## II. METHOD

### A. Configuration for the simulation

The overall configuration that includes the probe, the inspected volume and virtual sources is illustrated in Fig. 1 and described thereafter. A fully populated 2-D probe of 1,024 (32 × 32) elements is considered and simulated using the Matlab based toolbox FIELD II [36], [37]. The pulse-echo bandwidth of the elements is 40%. The inter-element pitch is set to $λ/2 \approx 257$ µm, with $λ = c_p/f_0$ where $c_p = 1540$ m/s denotes the propagation velocity of the compressional waves, and $f_0 = 3$ MHz the center frequency of the probe. The elements are separated by a $λ/20$ kerf. The total aperture of the probe (footprint) is 8.4 mm × 8.4 mm.

A total of 25 virtual sources (VS) are used in the proposed study. This number is considered high enough for capturing

<figure><img src="image_1.png" /><figcaption>Fig. 1: Configuration of the study, from top to bottom: positioning of the 25 virtual sources (VS), 2-D probe of 32 × 32 elements, nine scatterers used for the optimization and a pyramidal region used for a volumetric analysis in the results section</figcaption></figure>

geometrical tendencies in the distributions while keeping a high relative imaging volume rate. This should safely result in a sufficiently high imaging volume rate for 3-D cardiac applications such as 3-D ultrafast echocardiography. The optimization process, presented later in this section, is responsible for distributing the VS inside a “box”. This box bounds the VS from −10 mm to 10 mm in the *x* and *y* directions and from −40 mm to −10 mm in depth. These positions are continuous variables. In the case of an ideal probe, we assume that a distribution with a given symmetry will guarantee the same symmetry in computed metrics over a given volume. On the contrary, a non-symmetrical VS distribution will lead to non-symmetrical quality metrics over the volume, which indicates a preferential direction with the probe. As we impose the homogeneity of quality metrics to be of primary concern, only symmetrical distributions are considered. Considering this, a central symmetry of order *N* over the *z* axis with *N* → ∞ would be required. However, in practice, homogeneity of quality metrics is limited by the probe square shape and element discretization. Thus, a bi-planar symmetry system is selected in the present study. The two planes of symmetry are the *xz* plan on *y* = 0 and the *yz* plan on *x* = 0. This symmetry
can be seen in Fig. 1 by looking at the “replicated VS” (in black) that are distributed according to the “VS” (in red).

The considered media is homogeneous, isotropic and without attenuation. To take into account a whole imaging volume, a total of 27 scatterers are distributed (Fig. 1) and used to compute 3-D Point Spread Functions (PSF). They are located at [20; 40; 60] (mm) in depth (z) and at [−5; 0; 5], [−10; 0; 10] and [−15; 0; 15] (mm) in lateral directions x and y for the indicated depth location respectively. Regarding the symmetry system imposed on the VS distribution and the previous assumption on the quality metrics, only 9 of the 27 scatterers are used for the optimization. The 25 VS are used for each PSF to construct. Time-domain signals received by the elements, usually called Radio Frequency (RF) signals, are simulated using FIELD II with a sampling frequency of 48 MHz. Each PSF is calculated individually in order to decrease the computational load, instead of putting all the scatterers and then construct a whole imaging volume. For this reason and with a given VS distribution, all the RF signals are computed for each VS, probe element and scatterer. So, for a given VS distribution, it results in the calculation of 230,400 RF signals (number of VS × number of probe elements × number of individual scatterers). A standard time-domain delay-and-sum beamforming is used to construct envelope volumes using a Graphics Processing Unit [38] with Nvidia CUDA [39]; these signals are log-compressed afterward. No apodization has been used for reconstruction, in order to avoid favorizing or deteriorating the imaging performance of a given fraction of the imaging volume.

## B. Imaging criteria and extracted metrics

In order to assess other geometries of VS with respect to quality metrics and their homogeneity over an imaging volume, two quantitative metrics are defined on PSF. The first one measures the width of the main lobe at −6 dB (Full Width at Half Maximum, FWHM) on log-compressed volume; the second measures the amplitude of the highest secondary lobe (Secondary Lobes Level, SLL) on log-compressed volume.

A total of 9 individual scatterers are used in order to construct 9 PSF in order to represent a given volume; each PSF is used to measure a FWHM and SLL values. The metric that measures the FWHM is based on a measurement of the position of the voxels that exceed −6 dB (volumes are log-compressed). All voxels above −6 dB are projected along an axis being the direction between the center of the probe and the position of the scatterer. The mean diameter of that 2-D projection is then computed and used as the lateral resolution metric, similarly to FWHM in 2-D echography that is often used to assess the lateral resolution of an imaging system [40]. Regarding the geometry of the probe and the beamforming described previously, a natural expansion of the lateral resolution nearly proportional to the depth is usually observed. As the study focuses on the homogeneity of the quality metrics over the volume, the measurement of the FWHM in mm is converted into an angle ($\theta_{FWHM}$) viewed

from the center of the probe:

$$ \theta_{\text{FWHM}} = 2 \tan^{-1} \left( \frac{\text{FHWM}}{2z} \right). \qquad (1) $$

An SLL detection is performed on each PSF (log-compressed) volumes in order to extract the second metric. High SLL reduces the contrast of an ultrasound imaging system [41], that is why the latter is considered as a second quality metric. This estimation is done using a 3-D lobes' detection toolbox (MinimaMaxima3D [42]) in Matlab, and is displayed in dB relatively to the amplitude of the main lobe.

## C. Multiobjective optimization

For each of the 9 scatterers, the values of $\theta_{FWHM}$ and SLL are concatenated into two datasets referred as $\mathcal{D}_{FWHM}$ and $\mathcal{D}_{SLL}$:

$$ \begin{aligned} \mathcal{D}_{\text{FWHM}} &= [\theta_{\text{FWHM1}}, \dots, \theta_{\text{FWHM9}}] (^{\circ}), \\ \mathcal{D}_{\text{SLL}} &= [\text{SLL}_1, \dots, \text{SLL}_9] (\text{dB}). \end{aligned} \qquad (2) $$

For a first set of objective functions in the optimization, the assessment of the homogeneity of these datasets is done with a dispersion measurement by computing the standard deviation $\sigma$:

$$ \begin{aligned} \sigma_1 &= \sigma(\mathcal{D}_{\text{FWHM}}), \\ \sigma_2 &= \sigma(\mathcal{D}_{\text{SLL}}). \end{aligned} \qquad (3) $$

Concurrently, a second set of objective functions is suggested, based on the maximum value on the datasets:

$$ \begin{aligned} \max_1 &= \max(\mathcal{D}_{\text{FWHM}}), \\ \max_2 &= \max(\mathcal{D}_{\text{SLL}}). \end{aligned} \qquad (4) $$

For the first set of objective functions described in Eq. (3), the spread across the different scatterers is minimized by reducing the standard deviation, allowing a direct reduction of the inhomogeneity across the quality metrics of the 9 PSF. However, when considering the second set of objective functions, the idea is to increase the lower values of FWHM and SLL in the dataset to guarantee a maximal threshold for FWHM and SLL. Finally, the two objectives for FWHM and SLL homogeneity are concatenated into multiobjective vectors of two values, $[\sigma_1, \sigma_2]$ and $[\max_1, \max_2]$.

The decision space (also referred to as the search space [31]) is restricted in order to reduce the complexity of the objective functions and to accelerate the optimization process [43]. Here, 75 input variables are required so that 25 VS can freely move into a 3-D Cartesian space. A first restriction is applied by limiting the search space. Physical boundaries have been set up on the possible VS localization ($-10 \le x, y \le 10$, $-40 \le z \le -10$ (mm)) as indicated in Fig. 1. In the present study, one VS moves in the z direction and is set at the center of the probe ($x = y = 0$) as suggested in the literature [15], [16], [18], [21]. This leads to 24 3-D Cartesian coordinates for the 24 remaining VS. Applying the bi-planar symmetry system, as illustrated in Fig.1, reduces this number from 24 to 6. Taking all these considerations into account reduces the amount of input variables from 75 to 19. As mentioned in the previous section, a vector of two objectives describing FWHM
and SLL homogeneity ([σ₁, σ₂] or [max₁, max₂]) is considered for optimization.

A schematic view of the optimization procedure is presented in Fig. 2, illustrating an application of MOGA to this problem. First, an initial population is created by the algorithm: several

<figure><img src="image_2.png" /><figcaption>Fig. 2: Overview of a multiobjective genetic algorithm applied to the VS distribution</figcaption></figure>

VS distributions are generated randomly. The population is set at 1,000 individuals. From the population, each individual (VS distribution) is evaluated to extract either the dispersion or the maximum value of the quality metrics datasets. Depending on how good the individuals perform, they might be selected to create children (new solutions) for the next population.

In the study, this selection is made by a tournament of 4 individuals. The selected individuals (parents) are used to generate new individuals (children) by means of crossover. In this study, the crossover fraction is set at 40% (40% of the new population is produced by crossover). The remaining individuals of the new population in this algorithm are created by means of a mutation function. Mutation function is here operated with the mutationadaptfeasible function from the MOGA, that randomly selects a variable and modifies its value by adding a step, which depends on the boundaries set on variables and the current search space. The optimization is stopped according to the evolution of non-dominated solutions, generation after generation. Non-dominated solutions are here considered optimal enough when only a few of them have been improved after 25 generations, and when this improvement

is under 5% when compared to the previous one. A total of 250 generations have thus been reached when respecting these conditions. At the end of the optimization, non-dominated solutions are picked up from the last generated population to be the optimized solutions of the problem. They are thus shown and analyzed in the next section.

## D. Deterministic distributions

Three shapes of deterministic VS distributions with a constant height (height = -z) are implemented to further assess the imaging performance obtained using the optimized VS distributions. All the present deterministic distributions are, as well as the optimized ones, constituted by a total of 25 VS. The first deterministic distribution to be investigated in the section III is a 5 × 5 regular grid of VS. The second distribution contains two concentric rings and a central VS (x = y = 0). Each configuration is then declined into three different configurations in terms of aperture size (5 mm, 10 mm and 20 mm) and height (10 mm, 20 mm and 40 mm). These declination have been chosen to have a set of deterministic distributions that can be compared to the optimized distribution, regarding to the limit set on the position of the VS and also to the 9 scatterers taken into account (Fig. 1).

## E. Analysis of the imaging performance

1) PSF simulation: PSFs are first used to validate the performance of the optimized distributions of VS. To this end, a single scatterer is located at [0; 0; 40] mm (x; y; z) in an homogeneous and anechoic media. PSFs are displayed in planar surface C-scans with a constant depth z = 40 mm with -7 ≤ x, y ≤ 7 (mm).

2) Phantom with anechoic cysts and reflectors: All VS distributions are used to reconstruct images from a simulated phantom that contains point-like reflectors and spherical anechoic cysts. Reflectors are used to assess lateral resolution with the full width at -6 dB on log-compressed images. Reflectors are scatterers which have a reflection coefficient of 5. The cyst phantom is a composition of scatterers which values are uniformly distributed among the phantom, with a scatterer density of 4 scatterers/mm³. The reflection coefficients of these scatterers are normally distributed (σ² = 1, μ = 0), and scaled to be between -1 and 1. The anechoic cysts have no scatterers and they are used to assess Contrast Ratio (CR) on envelop images, here expressed by:

$$CR = 20 \log_{10} \left( \frac{\mu_{in}}{\mu_{out}} \right), \qquad (5)$$

where μ and σ are respectively the mean and the standard deviation of the pixels' values inside (in) or outside (out) cysts. Different cysts radius are used to avoid the edges of the cysts when taking the inside region while assessing CR values. To that end, we choose a minimum cyst radius of 2 mm plus an added value, corresponding to an approximated lateral resolution value (1.2λF [44]). All images are simulated with the same configuration introduced in section II.
3) **Volumetric analysis based on FWHM and SLL:** All VS distributions are used for imaging over a given volume. The motivation for this analysis is to give hints on volumetric imaging performance that includes out-of-plane locations in the volume. This consists in evaluating the two quality metrics (FWHM and SLL) used for the optimization, explained in section II-B. These metrics come from PSFs that are calculated from a single scatterer which is moved at several locations in order to have numerous PSFs that cover the whole pyramidal region shown in Fig. 1. The sector angle of the pyramidal region is fixed at 40° (20° from the normal to the probe center); the height is 45 mm with a base starting at $z = 65$ mm and a top ending at $z = 20$ mm. This leads to a $47.3 \text{ mm} \times 47.3 \text{ mm}$ base and a $14.6 \text{ mm} \times 14.6 \text{ mm}$ truncated apex. The single scatterer is moved inside this volume with a step $\delta_{x,y,z} = 2.5$ mm; a total of 3,131 PSFs are thus generated, allowing to extract 3,131 FWHM and SLL values. Because large bi-variant datasets need to be analyzed, bagplots [45], [46] have been chosen to summarize information. They are based on half-space depth [47] which acts as a median value for bi-variant data.

### III. RESULTS AND DISCUSSION

#### A. Optimization results

A total number of 127 and 143 optimal solutions have been generated respectively for the optimization based on standard deviation [$\sigma_1, \sigma_2$] and maximum value [$\max_1, \max_2$]. The two sets of non-dominated solutions are plotted in Fig. 3 in their respective objective space, which represents the outputs of the two objective functions, [$\sigma_1, \sigma_2$] and [$\max_1, \max_2$] (3). Fig. 3 shows the compromise between the two objectives: when $\sigma_1$ or $\max_1$ decreases, $\sigma_2$ or $\max_2$ increases and vice versa. All the solutions for the standard-deviation-based ($\sigma$-based) optimization ensure a standard deviation of, at worst, 0.34° in $D_{\text{FWHM}}$ and 1.19 dB in $D_{\text{SLL}}$. The solutions based on maximum value ensure at least 3.89° in FWHM and -17.4 dB in SLL, for the 9 selected scatterers. The two fronts of non-dominated solutions, often called Pareto-front [32], [33], are discontinuous. Such discontinuities are usual in many two-objective problems: they exist because of the definition of non-dominated solutions [30] and the landscape of feasible solutions of the multiobjective function. Specific solution groups are then plotted in different colors based on a visual check of similar VS distribution. Solutions have then been classified into 3 distinct groups for the $\sigma$-based optimization and 4 distinct groups for the max-based optimization (Fig. 3). Due to the number of optimized solutions and the similarity observed on the early stage that is based on a visual check, a single solution is selected for each group. These selected solutions are encircled in Fig. 3 and represent their group. These solutions are systematically the best for FWHM ($\sigma(D_{\text{FWHM}})$ and $\max(D_{\text{FWHM}})$) and the worst for SLL ($\sigma(D_{\text{SLL}})$ and $\max(D_{\text{SLL}})$). These specific VS distributions are presented in Fig. 4, when looking toward the z direction.

Fig. 4 shows the deterministic distributions from Fig. 4.a to Fig. 4.i, previously described in II-D. The optimized distributions that have been identified in Fig. 3 are displayed from

<figure><img src="image_3.png" /><figcaption>Fig. 3: (a) Standard-deviation ($\sigma$) and (b) maximum value (max) optimized solutions on their respective objective space. Solutions are categorized into groups and displayed in different colors, from top left to bottom right on the front, respectively to the shape of the corresponding VS distributions. A single solution from each group color is encircled for further analysis.</figcaption></figure>

Fig. 4.j to Fig. 4.l for the optimization based on the standard deviation (3) and from Fig. 4.m to Fig. 4.p for the optimization based on the maximum (4).

Fig 4 shows how optimized VS distributions are different from each other and from deterministic distributions. For example in Fig. 4.j, the first distribution of the standard-deviation-based objective function (named $\sigma$-n1) has half the VS located close to the x and y boundaries and close to the probe constituting a pseudo-circle; the second half of VS are located close to the center of the probe in the (x, y) plane. However, distribution max_n4 in Fig. 4.p has two distinguishable layers: the first layer contains 9 VS that forms a cross and the second one has the 16 remaining VS that nearly form a regular grid close to the probe.

#### B. Imaging results

1) **PSF simulation:** The planar surfaces PSF are displayed in Fig. 5. The highest secondary lobe level observed among deterministic distributions is -22 dB for distribution
<figure><img src="image_4.png" /><figcaption>Fig. 4: Shapes of the VS distribution: deterministic distributions are displayed from Fig. 4.a to Fig. 4.i, and optimized distribution from Fig. 4.j to Fig. 4.p. Color and size of each source indicate the height (−z) of the VS: the smaller and darker, the closer to the probe (−10 mm).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5: PSFs over a planar surface (z = 40 mm). The −6 dB isoline is highlighted by a white line.</figcaption></figure>

spiral_10_mm; by comparison, the highest secondary lobe
level observed among optimized distributions is −19 dB with
distribution max_n1. Each PSF in Fig. 5 has its −6 dB isoline
highlighted by a white line. Deterministic distributions offer
in most cases a more isotropic PSF than optimized distribu-
tions: the shape of the −6 dB isoline is circular when using
deterministic distributions (from Fig 5.a to from Fig 5.i), while
this shape is more elliptical for most optimized distributions,

expect for max_n3 and max_n4.

2) Phantom with anechoic cysts and reflectors: Images from the simulated phantom are presented in Fig. 6 for the xz plan. Fig. 6 let us see that the visibility of point-like reflectors and cysts depend on their location as well as the VS distribution that is being used. FWHM (Eq. (1)) measured on reflectors and CR (Eq. (5)) measured on cysts are displayed as boxplots in Fig. 7. The lowest FWHM value observed is 2.3° for
<figure><img src="image_6.png" /><figcaption>Fig. 6: Images in the *xz* plan from the simulated phantom with cysts and reflectors</figcaption></figure>

distribution `max_n1` and the highest is 6.9° for distribution `spiral_10_mm`. Distributions `max_n2` and `max_n3` are the only two that have all their FWHM values below 4°. Distributions `max_n3` and `max_n4` have a mean FWHM of 16% and 9% lower (3.24° and 3.50° respectively), in comparison to distribution `rings_20_mm` (3.86°) while keeping equivalent mean contrasts values measured in cysts (–19.6 dB, –19.5 dB and –19.3 dB respectively for distributions `rings_20_mm`, `max_n3` and `max_n4`). The lowest CR value measured when using a deterministic VS distribution is with distribution `spiral_20_mm` (–23.7 dB) while the lowest CR value measured when using an optimized VS distribution is with distribution `max_n1` (–29.8 dB). The latter also have the highest CR value as well as distribution `reg_40_mm` (–12.7 dB and –12.5 dB respectively).

Fig. 7 shows an improvement of the FWHM with deterministic distributions as the aperture of the VS distribution increases. This tendency has been reported in [21]: the wider the aperture of the VS distribution, the lower the FWHM. Utilizing spirals as a VS distribution instead of a regular grid improves CR; this is consistent with findings in [23], where CR was enhanced through the use of spiral-based and aperiodic ring-based distributions for plane waves.

3) **Volumetric analysis based on FWHM and SLL:** Bagplots, based on the measurement of FWHM and SLL for 3,131

different PSFs, are plotted in Fig. 8 for all distributions. Optimized VS distributions (from Fig. 8.j to Fig. 8.p,) always offer a lower FWHM value of the half-space depth than deterministic distributions (from Fig. 8.a to Fig. 8.i). The half-space depth with the lowest FWHM value between deterministic distributions stands for `reg_10_mm`, with a FWHM of 3.68°; in comparison, the half-space depth with the highest FWHM among optimized distribution is 3.59° for the distribution `σ_n1` and this value decreases down to 3.26° for both distributions `max_n1` and `max_n2`. Fig. 8.a and Fig. 8.p show half-space depths with the lowest SLL value (–23.3 dB). Sizes of bags and fences of optimized distributions (from Fig. 8.j to Fig. 8.p) are often wider than bags and fences constructed with deterministic distributions with the height of the VS being fixed at 10 mm (from Fig. 8.a to Fig. 8.c), but lower in comparison to deterministic distributions when the height is fixed at 40 mm (from Fig. 8.g to Fig. 8.i). Distributions `rings_10_mm` and `spiral_10_mm` have the lowest sizes of both bags and fences while looking at SLL. To us, the flattening of bags and fences reflects how quality metrics are homogeneous. By following this logic, distributions `rings_10_mm` and `spiral_10_mm` offer the most homogeneous distribution of SLL when imaging the pyramidal region indicated in Fig. 1. Then, optimized distribution offer then a more homogeneous distribution of SLL than deterministic distribution when the height of the VS
<figure><img src="image_7.png" /><figcaption>**Fig. 7:** Boxplots: CR measured on the 7 cysts and FWHM measured on the 12 reflectors from the phantom images (Fig. 6) in the xz and yz plan for all VS distributions. Mean values are indicated by a cross mark.</figcaption></figure>

is 20 mm or 40 mm. Optimized distribution also offer a more homogeneous distribution of FWHM than these deterministic distributions in most cases. Bagplots can also help assess the choice of the 9 scatterers that have been used for the optimization. To give an example, distributions σ_n1 and σ_n3 are used in the following. Bags and fences that come from Fig. 8.j and Fig. 8.l are plotted in Fig. 9. The FWHM and SLL value that comes from the 9 scatterers of the optimization are indicated by crosses, to make a comparison with bags and fences. When looking at bagplots, the bag and fence of distribution σ_n1 are flattening in the FWHM direction than SLL direction, in opposition with distribution σ_n3. By looking at these two examples, the 9 scatterers used for the optimization follow these flattening tendencies: they have a small variation in FWHM when FWHM's standard deviation is prioritized, as well as SLL when the SLL's standard deviation is prioritized. To us, this is an indicator that selecting those 9 scatterers was a suitable option to drive the optimization and has an impact on the whole pyramidal region presented in Fig. 1.

### C. Limitations and perspectives

The study has been conducted through simulations; further *in-vitro* and *in-vivo* experiments might validate the optimized distributions that have been found. Nonetheless, several studies [21], [48] have demonstrated agreement between *in-silico* and *in-vitro/in-vivo* conditions for both 2-D and 3-D ultrafast diverging wave compounding. Such agreement is also observed for 3-D plane wave compounding [23]. The enhancements that

the optimized distributions have shown in our study should be comparable between *in-silico* simulations and potential upcoming *in-vitro/in-vivo* experiments. In the presence of moving objects, imaging performances could be enhanced by applying a motion compensation (MoCo) technique as to the one presented in [19], which demonstrated optimality with 25 virtual sources. Moreover, a more recent 3-D MoCo technique shown in [48] could possibly be applied. These optimized virtual sources distributions are anticipated to be suitable for cardiac applications including 3-D echocardiography.

Optimizations and *in-silico* results have been conducted without noise. In [21] and [18], 5 and 25 virtual sources have been used, respectively, in 2-D and 3-D diverging wave compounding imaging *in-vitro*, showing a decent signal-to-noise ratio (SNR). Additionally, we employed all the elements of the probe for both transmission and reception in our simulations without any apodization. Therefore, the optimized distributions should provide a high SNR within an angular aperture of 30° (Fig. 1). Beyond 30°, the SNR could decrease rapidly; the lateral positions of the virtual sources may be adjusted to achieve a desired angular aperture. In cases of suboptimal SNR, coded excitation [49] could be employed to enhance the SNR and mitigate strong attenuation.

One of the goals of the bi-planar symmetry system described section II-A was to guarantee the optimized distributions of virtual sources to offer the same level of symmetry, for the distribution of the quality metrics over a given volume. The scope of the study has thus excluded the search for distributions like concentric rings if it does not respect this symmetry, as well as spirals. The constraints imposed to the current framework could be adapted to search for spiral-based or concentric-based distributions, that should offer a all-roundly symmetrical distribution of quality metrics. In addition, this bi-planar symmetry system has resulted in optimized distributions that yield an anisotropic PSF, except for distributions max_n3 and max_n4 (as shown in Fig 5.o and Fig 5.p). This has led to anisotropic resolutions and could potentially lead to inconsistencies in imaging performance across various lateral directions. In order to maintain a PSF that is as isotropic as possible, we recommend utilizing distributions max_n3, max_n4, or a deterministic distribution.

Section II-A indicated that no apodization has been used in transmission and reception. In [22], [23] the imaging point is reconstructed only if the main beam of the wave propagates through it; employing a beamforming strategy akin to this could further reduce SLL. Integrating this into the optimization framework may yield additional optimized distributions.

The multiobjective functions used in the optimization framework has focused on isolated PSFs. This could also be adapted to take into account different metrics, like the homogeneity of PSF intensity.

## IV. CONCLUSION

The study has introduced novel distributions of virtual sources for 3-D ultrafast diverging wave compounding imaging, discovered through the application of a multiobjective optimization framework. Some of these new distributions have
<figure><img src="image_8.png" /><figcaption>Fig. 8: Bagplot from all VS distributions for 3,131 FWHM and SLL values measured in the pyramidal region indicated in Fig. 1.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 9: Comparison between the total volume performance from two optimized distributions and the 9 scatterers used during the optimization, using bagplots.</figcaption></figure>

exhibited an enhancement in lateral resolution, reducting it down to 16% in comparison to deterministic distributions that possess similar contrast and PSF isotropy properties. These optimized distributions could be employed in cardiac applications to provide an improvement in imaging performances.

## REFERENCES

[1] M Cheitlin, W Armstrong, G Aurigemma, G Beller, F Bierman, J Davis, P Douglas, D Faxon, L Gillam, and T Kimball. ACC/AHA/ASE

2003 Guideline Update for the Clinical Application of Echocardiography: Summary Article*1A report of the American College of Cardiology/American Heart Association Task Force on Practice Guidelines (ACC/AHA/ASE Committee to update the 1997 guidelines for the clinical application of echocardiography). *Journal of the American Society of Echocardiography*, 16(10):1091–1110, October 2003.

[2] Anne Long. Guide pratique d'écho-Doppler vasculaire. ELSEVIER, 2017.

[3] Gerhard H. Mostbeck, Thomas Zontsich, and Karl Turetschek. Ultrasound of the kidney: obstruction and medical diseases. *European Radiology*, 11(10):1878–1889, October 2001.

[4] Claudia Errico, Bruno-Félix Osmanski, Sophie Pezet, Olivier Couture, Zsolt Lenkei, and Mickael Tanter. Transcranial functional ultrasound imaging of the brain using microbubble-enhanced ultrasensitive Doppler. *NeuroImage*, 124:752–761, January 2016.

[5] Jørgen A. Jensen, Svetoslav Nikolov, Alfred C.H. Yu, and Damien Garcia. Ultrasound Vector Flow Imaging: I: Sequential Systems. *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, pages 1–1, 2016.

[6] Mickael Tanter and Mathias Fink. Ultrafast imaging in biomedical ultrasound. *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, 61(1):102–119, January 2014.

[7] A. Austeng and S. Holm. Sparse 2-D arrays for 3-D phased array imaging - design methods. *IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control*, 49(8):1073–1086, August 2002.

[8] Alessandro Ramalli, Enrico Boni, Emmanuel Roux, Hervé Liebgott, and Piero Tortoli. Design, implementation, and medical applications of 2-d ultrasound sparse arrays. *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, 69(10):2739–2755, 2022.

[9] H. Ermert, O. Keitmann, R. Oppelt, B. Granz, A. Pesavento, M. Vester, B. Tillig, and V. Sander. A new concept for a real-time ultrasound transmission camera. In *2000 IEEE Ultrasonics Symposium. Proceedings. An International Symposium (Cat. No.00CH37121)*, volume 2, pages 1611–1614, San Juan, Puerto Rico, 2000. IEEE.

[10] Morten F. Rasmussen and Jørgen A. Jensen. 3D ultrasound imaging performance of a row-column addressed 2D array transducer: a simulation study. In Johan G. Bosch and Marvin M. Doyley, editors, *2013 IEEE International Ultrasonics Symposium (IUS)*, page 86750C, Lake Buena Vista (Orlando Area), Florida, USA, mar 2013.

[11] Pedro Santos, Geir Ultveit Haugen, Lasse Lovstakken, Eigil Samset, and Jan D'hooge. Diverging Wave Volumetric Imaging Using Subaperture Beamforming. *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, 63(12):2114–2124, December 2016.
[12] C. Chen, Z. Chen, D. Bera, S. B. Raghunathan, M. Shabanimotlagh, E. Noothout, Z.Y. Chang, J. Ponte, C. Prins, H.J. Vos, J.G. Bosch, M.D. Verweij, N. de Jong, and M.A.P. Pertijs. A front-end ASIC with receive sub-array beamforming integrated with a 32 × 32 PZT matrix transducer for 3-D transesophageal echocardiography. In *2016 IEEE Symposium on VLSI Circuits (VLSI-Circuits)*, pages 1–2, Honolulu, HI, USA, June 2016. IEEE.

[13] Ling Tong, Alessandro Ramalli, Ruta Jasaityte, Piero Tortoli, and Jan D'hooge. Multi-transmit beam forming for fast cardiac imaging—experimental validation and in vivo application. *IEEE Transactions on Medical Imaging*, 33(6):1205–1219, 2014.

[14] Mustafa Karaman, Ira O. Wygant, Ömer Oralkan, and Butrus T. Khuri-Yakub. Minimally redundant 2-d array designs for 3-d medical ultrasound imaging. *IEEE Transactions on Medical Imaging*, 28(7):1051–1061, 2009.

[15] G. Montaldo, M. Tanter, J. Bercoff, N. Benech, and M. Fink. Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography. *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, 56(3):489–506, March 2009.

[16] Jean Provost, Clement Papadacci, Juan Esteban Arango, Marion Imbault, Mathias Fink, Jean-Luc Gennisson, Mickael Tanter, and Mathieu Pernot. 3D ultrafast ultrasound imaging in vivo. *Physics in Medicine and Biology*, 59(19):L1–L13, October 2014.

[17] J Sauvage, M Flesch, G Férin, A Nguyen-Dinh, J Porée, M Tanter, M Pernot, and T Deffieux. A large aperture row column addressed probe for *in vivo* 4D ultrafast doppler ultrasound imaging. *Physics in Medicine & Biology*, 63(21):215012, October 2018.

[18] Emmanuel Roux, François Varray, Lorena Petrusca, Christian Cachard, Piero Tortoli, and Hervé Liebgott. Experimental 3-D Ultrasound Imaging with 2-D Sparse Arrays using Focused and Diverging Waves. *Scientific Reports*, 8(1), December 2018.

[19] Yinran Chen, Jan D'hooge, and Jianwen Luo. Doppler-based motion compensation strategies for 3-d diverging wave compounding and multiplane-transmit beamforming: A simulation study. *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, 65(9):1631–1642, 2018.

[20] Bastien Denarie, Thor Andreas Tangen, Ingvild Kinn Ekroll, Natale Rolim, Hans Torp, Tore Bjästad, and Lasse Lovstakken. Coherent plane wave compounding for very high frame rate ultrasonography of rapidly moving targets. *IEEE Transactions on Medical Imaging*, 32(7):1265–1276, 2013.

[21] Clement Papadacci, Mathieu Pernot, Mathieu Couade, Mathias Fink, and Mickael Tanter. High-Contrast Ultrafast Imaging of the Heart. *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, 61(2):15, 2014.

[22] Sua Bae, Jiwon Park, and Tai-Kyong Song. Contrast and Volume Rate Enhancement of 3-D Ultrasound Imaging Using Aperiodic Plane Wave Angles: A Simulation Study. *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, 66(11):1731–1748, 2019.

[23] Sua Bae, Bae-Hyung Kim, Azra Alizad, Mostafa Fatemi, and Tai-Kyong Song. Experimental study of aperiodic plane wave imaging for ultrafast 3-d ultrasound imaging. *IEEE Transactions on Biomedical Engineering*, 69(8):2679–2690, 2022.

[24] Emmanuel Roux, Alessandro Ramalli, Piero Tortoli, Christian Cachard, Marc C. Robini, and Herve Liebgott. 2-D Ultrasound Sparse Arrays Multidepth Radiation Optimization Using Simulated Annealing and Spiral-Array Inspired Energy Functions. *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, 63(12):2138–2149, December 2016.

[25] L Haupt. Thinned Arrays Using Genetic Algorithms. *IEEE Transactions on Antennas and Propagation*, page 7, 1994.

[26] Florent Le Courtois, Jean-Hugh Thomas, Franck Poisson, and Jean-Claude Pascal. Genetic optimisation of a plane array geometry for beamforming. Application to source localisation in a high speed train. *Journal of Sound and Vibration*, 371:78–93, June 2016.

[27] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by Simulated Annealing. *Science*, 220(4598):671–680, May 1983.

[28] John H. Holland. *Adaptation in natural and artificial systems: An introductory analysis with applications to biology, control, and artificial intelligence*. Adaptation in natural and artificial systems: An introductory analysis with applications to biology, control, and artificial intelligence. U Michigan Press, Oxford, England, 1975.

[29] James Kennedy and Russell Eberhart. Particle Swarm Optimization. In *Proceedings of ICNN'95 - International Conference on Neural Networks*, page 7, 1995.

[30] Kalyanmoy Deb. Multi-Objective Optimization Using Evolutionary

Algorithms: An Introduction. In *Multi-objective optimization using evolutionary algorithms*, page 25, 2001.

[31] Carlos A. Coello, Gary B Lamont, and David A Van Veldhuizen. *Evolutionary Algorithms for Solving Multi-Objective Problems Second Edition*. Springer US, 2007.

[32] Carlos A. Coello, Gregorio Toscano Pulido, and M Salazar Lechuga. Handling multiple objectives with particle swarm optimization. *IEEE Transactions on evolutionary computation*, 8(3):256–279, 2004.

[33] Kalyanmoy Deb, Samir Agrawal, Amrit Pratap, and T Meyarivan. A Fast Elitist Non-dominated Sorting Genetic Algorithm for Multi-objective Optimization: NSGA-II. In G. Goos, J. Hartmanis, J. van Leeuwen, Marc Schoenauer, Kalyanmoy Deb, Günther Rudolph, Xin Yao, Evelyne Lutton, Juan Julian Merelo, and Hans-Paul Schwefel, editors, *Parallel Problem Solving from Nature PPSN VI*, volume 1917, pages 849–858. Springer Berlin Heidelberg, Berlin, Heidelberg, 2000.

[34] Víctor Martínez-Cagigal. Multi-objective particle swarm optimization (mopso), 2019.

[35] MathWorks ©. Global optimization toolbox, 2021.

[36] J.A. Jensen and N.B. Svendsen. Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers. *IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control*, 39(2):262–267, March 1992.

[37] Jørgen Jensen. Field: A program for simulating ultrasound systems. *Medical and Biological Engineering and Computing*, 34:351–352, 01 1996.

[38] NVIDIA. nVidia V100SXM2. https://www.nvidia.com/en-us/data-center/v100/, 2018. nVidia. Retrieved October 21, 2020.

[39] NVIDIA, Péter Vingelmann, and Frank H.P. Fitzek. CUDA, release: 10.2.89. https://developer.nvidia.com/cuda-toolkit, 2020. nVidia CUDA. Retrieved October 21, 2020.

[40] H. Liebgott, A. Rodriguez-Molares, F. Cervenansky, J.A. Jensen, and O. Bernard. Plane-Wave Imaging Challenge in Medical Ultrasound. In *2016 IEEE International Ultrasonics Symposium (IUS)*, pages 1–4, Tours, France, September 2016. IEEE.

[41] Thomas L. Szabo. *Diagnostic Ultrasound Imaging: Inside Out*. Academic Press, 2004.

[42] Samuel Pichardo. MinimaMaxima3D. https://www.mathworks.com/matlabcentral/fileexchange/17997-minimamaxima3d, 2020. MATLAB Central File Exchange. Retrieved October 21, 2020.

[43] Marc C. Robini. Theoretically Grounded Acceleration Techniques for Simulated Annealing. In Ivan Zelinka, Václav Snášel, and Ajith Abraham, editors, *Handbook of Optimization*, volume 38, pages 311–335. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.

[44] F.J. Harris. On the use of windows for harmonic analysis with the discrete fourier transform. *Proceedings of the IEEE*, 66(1):51–83, 1978.

[45] Peter J. Rousseeuw, Ida Ruts, and John W. Tukey. The Bagplot: A Bivariate Boxplot. *The American Statistician*, 53(4):382–387, November 1999.

[46] Mark Geurts. bagplot.m [source code]. https://github.com/mwgeurts/libra/blob/master/bagplot.m, 2018. [Online; accessed March-2021].

[47] J. W. Tukey. Mathematics and the picturing of data. *Proceedings of the International Congress of Mathematicians, Vancouver, 1975*, 2:523–531, 1975.

[48] Yinran Chen, Zichen Zhuang, Jianwen Luo, and Xiongbiao Luo. Doppler and pair-wise optical flow constrained 3d motion compensation for 3d ultrasound imaging. *IEEE Transactions on Image Processing*, pages 1–1, 2023.

[49] Feifei Zhao, Ling Tong, Qiong He, and Jianwen Luo. Coded excitation for diverging wave cardiac imaging: A feasibility study. *Physics in Medicine & Biology*, 62(4):1565, 2017.