Optimization of array encoding for ultrasound imaging

Techniques in machine learning have gained a significant amount of research attention in the field of ultrasound imaging, particularly in the application of deep neural networks (DNNs) (Liu et al., 2019; Hyun et al., 2021). Following the unprecedented successes of the deep convolutional neural network (CNN) AlexNet in non-medical imaging tasks (Krizhevsky et al., 2017), similar networks are frequently trained and utilized for ultrasound-adjacent tasks as an additional post-processing step to a conventional imaging pipeline by performing classification, segmentation, etc. on image data (Mischi et al., 2020).

Other deep learning approaches perform image formation directly. In a prototypical setup, a DNN is fed the received echoes of a transducer array to produce a B-mode image, aiding or even supplanting the use of a conventional beamformer (van Sloun et al., 2020; Goudarzi and Rivaz, 2022). Often the goal is to create a network that can be evaluated faster than classical imaging techniques without degrading image quality (Chen et al., 2021). Other times, the DNN is designed to directly improve image quality. For example, Hyun et al. (2019) utilize a CNN as a beamformer that reduces speckle, citing considerable advantages over conventional delay-and-sum techniques.

While a well-trained DNN can offer these and other benefits, their use necessarily introduces complications. Most notably, neural networks are susceptible to hallucinations, in which it unexpectedly generates features that are not present in the underlying data, an issue that is uniquely consequential in medical imaging (Bhadra et al., 2021). Our approach to machine learning entirely circumvents this and other issues by using a custom architecture whose forward operation is entirely acoustically motivated. Specifically, the proposed ML method improves imaging quality exclusively through an optimized parameterization of experimental hardware via the transmit sequence. After data is acquired with an optimized sequence, the B-mode image to be evaluated is generated through the REFoCUS imaging framework, which itself is based solely on the linear nature pulse-echo responses (We describe the decoding and imaging processes in detail in Section 2.1 and Section 2.2 respectively).

A similar strategy to our own is employed by Chen et al. (2021), in which a machine learning model is trained to optimize the apodization weights of planewave transmissions, although a DNN is used to decode the received echoes. Specifically, a stacked denoising autoencoder architecture (a type of DNN whose first trainable layer is taken to be the weights themselves) is trained to optimize the reconstruction of the multistatic data set. In contrast, our ML model allows full control of the encoding sequence through a continuous treatment of apodization weights and time delays, with no other trainable parameters. At the same time, our training workflow optimizes the encoding sequence according to direct measures of image quality. In doing so, we provide notable counterexamples to the principle that a better reconstruction of the multistatic data set leads to higher quality images.

By detaching our machine learning model from conventional neural network architectures, our approach obtains theoretical and practical advantages beyond avoiding hallucinations. First, an optimized encoding sequence generated by our ML model can be analyzed separately from the model itself, such that its beneficial features can be identified and further studied. This is in contrast to a trained DNN, whose per-layer weights and biases have no meaning outside the context of the DNN. Second, the smaller parameter space of our ML model dramatically reduces the training time needed to generate an optimized encoding sequence. This allows for greater flexibility in the overall ML workflow, as a different configuration of the ML model can be retrained at a minimal cost when the experimental setting changes (e.g., targeting a different viewing depth or prioritizing different image features via an alternate loss function). Finally, optimization of data acquisition through the transmit sequence makes our imaging process very robust to different targets. While conventional DNNs are typically trained for very specific imaging tasks on a constrained set of training data, we will show that the proposed approach not only generalizes well to out-of-sample training data, but also to data with dramatically different features. This is because our ML model is fundamentally optimized based on the underlying acoustic principles of ultrasound imaging, where the processing of backscattered echoes is influenced primarily through manipulation of the transmission sequence.

We consider a linear transducer array of $N_{E}$ physical elements, of which $N_{T}$ are capable of transmitting a diverging wave and $N_{R}$ can receive backscattered echoes. Moreover, it is typical to have elements that perform both functions, such that $N_{T}=N_{R}=N_{E}$. Taken across all pairs of transmit and receive elements, collected RF (radio frequency) signals form a time dependent $N_{T}\times N_{R}$ matrix $\mathbf{U}(t)$, within which the matrix element $u_{TR}(t)$ is the signal observed by transmitting on element $T$ and receiving on element $R$. In the synthetic transmit aperture model, the complete collection of pulse-echo response signals between pairs of transmit and receive elements make up the multistatic data set. Because these signals can be delayed and summed across the receive channel to produce transmit focus throughout the domain, as in delay-and-sum (DAS) beamforming, these data can be considered the mathematical basis of our imaging (Jensen et al., 2006).

As previously discussed, it is inadvisable to construct $\mathbf{U}(t)$ directly by firing each transmit element in sequence across a total of $N_{T}$ transmissions (Chiao et al., 1997). Instead, the REFoCUS technique seeks to construct an approximation of $\mathbf{U}(t)$ using the responses from an arbitrary scanning sequence. We consider a scanning sequence of $N_{M}$ separate transmissions, each of which is parameterized by $N_{T}$ time delays $t_{MT}$ and apodization weights $w_{MT}$. The response observed by receive element $R$ from transmit $M$ can be described as a weighted sum of time delayed matrix elements from the multistatic data set,

which are similarly collected into a time dependent $N_{M}\times N_{R}$ matrix $\mathbf{S}(t)$ of focused responses. Collected across each transmit, these time delays and apodization weights make up the encoding sequence, denoted $\boldsymbol{\sigma}=(t,w)_{MT}$.

It is now convenient to work in the frequency domain, where each time delay of $u_{TR}(t)$ becomes a complex phase shift of $u_{TR}(\omega)$ (Ali et al., 2020). For simplicity, we represent the frequency dependent Fourier transform of $\mathbf{U}(t)$ and $\mathbf{S}(t)$ as $\mathbf{U}(\omega)$ and $\mathbf{S}(\omega)$ respectively. In doing so, we can express the transmission response of receive element $R$ to transmission $M$ as

This describes a linear relationship between the multistatic data set $\mathbf{U}(\omega)$ and the transmission responses $\mathbf{S}(\omega)$ that is fully parameterized by $\boldsymbol{\sigma}$, given by

where $\mathbf{H}(\omega)$ is a frequency dependent $N_{M}\times N_{T}$ encoding matrix given by

In practice, we collect only discrete, equispaced samples from $\mathbf{S}(t)$, from which we have a total of $N_{\omega}$ angular frequencies $\omega$. Because calculations involving the collection of matrices $\mathbf{S}(\omega)$ are typically performed in parallel across frequencies, we use calligraphic letters to denote the collection of all frequencies. As an example, we denote $\boldsymbol{\mathcal{S}}_{\omega}=\mathbf{S}(\omega)$, and represent the above equation more compactly as $\boldsymbol{\mathcal{S}}=\boldsymbol{\mathcal{H}}\boldsymbol{\mathcal{U}}$. When we wish to emphasize the dependence of our encoding matrices $\boldsymbol{\mathcal{H}}$ on the sequence $\boldsymbol{\sigma}$, we add $\boldsymbol{\sigma}$ as a subscript.

The objective of the REFoCUS framework is to create an approximation of the multistatic data set $\widehat{\boldsymbol{\mathcal{U}}}$ from the recorded echoes in $\boldsymbol{\mathcal{S}}$. This is done by applying a frequency dependent decoder $\mathbf{H}^{\dagger}(\omega)$ to $\mathbf{S}(\omega)$, resulting in

This decoder is theoretically arbitrary, with options explored in the literature ranging from a per-frequency conjugate transpose (Bottenus, 2018) to a fully trained neural network (Chen et al., 2021). In this work, we follow Ali et al. (2020); Bottenus et al. (2023) and apply a Tikhonov regularized pseudoinverse that has been shown to work well experimentally. This choice of decoder depends only on the original encoding matrix, which is itself dependent only on our sequence $\boldsymbol{\sigma}$. The Tikhonov decoder is given explicitly by

where $\boldsymbol{\mathcal{H}}^{*}$ is the conjugate transpose of $\boldsymbol{\mathcal{H}}$ and $\gamma$ is a regularization constant. This regularization both suppresses the contribution of noise to the recovered multistatic channel data and produces a more stable inversion of $\boldsymbol{\mathcal{H}}$, which together results in a more accurate and robust recovery of $\widehat{\boldsymbol{\mathcal{U}}}$. This is particularly important in the underdetermined system of interest, in which the number of transmits $N_{M}$ is significantly smaller than the number of array elements $N_{E}$. In all usages of regularization herein we take a value of $\gamma=0.1\sigma_{\text{max}}(\omega)$, where $\sigma_{\text{max}}(\omega)$ is the spectral norm of each per-frequency encoding matrix $\mathbf{H}(\omega)$. This particular value has been observed empirically to avoid issues of over- or under-regularization in the presence of noise (Bottenus et al., 2023).

Given a multistatic data set $\boldsymbol{\mathcal{U}}$, we produce the 2D B-mode image of interest with conventional DAS beamforming. This technique maps ground-truth or approximated multistatic IQ data (in-phase and quadrature data derived from the experimentally acquired RF data) at each frequency to a single IQ data matrix, taking a sum of the time delayed data in $\widehat{\boldsymbol{\mathcal{U}}}$. By manipulating these post-acquisition time delays, an image is created that is focused at individual pixels in the imaging domain. Mathematically, we denote this beamforming operation as $\mathcal{B}:\mathbb{C}^{N_{M}\times N_{R}\times N_{\omega}}\to\mathbb{C}^{N_{x}\times N_{z}}$, where $(N_{x},N_{z})$ are the number of pixels in the lateral and axial directions of our image.

To make the resulting image suitable for interpretation by a human observer, we apply a number of non-linear post-processing steps to this re-focused IQ data. These include the computation of the signal envelope, log-scaling the dynamic range, and clipping the image to a minimum decibel value. In our notation, we represent these operations collectively as $|\cdot|_{\text{Im}}:~{}\mathbb{C}^{N_{x}\times N_{z}}\to\mathbb{R}^{N_{x}\times N_{z}}$. When applied in the experimental context, the total forward operation of our ML model takes in transmission responses collected in $\boldsymbol{\mathcal{S}}$ and an encoding sequence $\boldsymbol{\sigma}$, and produces an image represented by the real matrix $|\mathcal{B}(\boldsymbol{\mathcal{H}}_{\boldsymbol{\sigma}}^{\dagger}\boldsymbol{\mathcal{S}})|_{\text{Im}}$.

To ensure that the generated encoding sequences is useful in a clinical setting, we train our ML model with a loss function $\mathcal{L}:~{}\mathbb{R}^{N_{x}\times N_{z}}\to~{}\mathbb{R}^{+}$ that directly measures image quality. While there is flexibility in the specific choice of loss function (Hyun, 2022), achieving a lower value for the loss should correspond to higher resolution, higher SNR, wider FOV, fewer artifacts, etc.

A clear choice is to use a loss function that quantifies these qualities directly. For example, the generalized contrast to noise ratio (gCNR) is a measure of histogram overlap between the brightness of predefined speckle and anechoic regions, such that a higher value indicates greater contrast in the image (Rodriguez-Molares et al., 2019). We can in principle improve the gCNR by minimizing the loss function

However, complications with implementation make this loss function impractical. In particular, the histogram operator is not differentiable and standard continuous approximations are unstable, which makes the metric incompatible with backpropagation during training.

On the other hand, a natural and easily implementable choice is a normalized $\ell^{2}$ comparison to a reference image that has desirable qualities. This introduces further flexibility in the choice of imaging target, which we discuss in Section 2.4 and Section 4.3. With a given image $\widehat{\mathbf{X}}$ and reference image $\mathbf{X}$, this loss is then defined by

Importantly, we will show that training the ML model to minimize $\mathcal{L}_{\ell^{2}}$ still makes a meaningful improvement to the gCNR of produced images.

We define an optimized encoding sequence $\boldsymbol{\sigma}^{*}$ as one that, for a given loss function $\mathcal{L}$, beamformer $\mathcal{B}$, and decoder $\boldsymbol{\mathcal{H}}^{\dagger}_{\boldsymbol{\sigma}}$, minimizes the non-convex optimization problem

where this expectation is taken over all possible transmission responses $\boldsymbol{\mathcal{S}}$. As it is impossible to evaluate this expectation directly, we instead approach the minimization problem from a machine learning perspective, where we instead minimize over a representative sample of training data $\{\boldsymbol{\mathcal{U}}_{i}\}$. Each piece of ground-truth multistatic data is related to a transmission response by our encoding matrix $\boldsymbol{\mathcal{H}}_{\boldsymbol{\sigma}}$ by $\boldsymbol{\mathcal{S}}_{i}=\boldsymbol{\mathcal{H}}_{\boldsymbol{\sigma}}\boldsymbol{\mathcal{U}}_{i}$.

Similarly, we must restrict the space of encoding sequences to those that are practically realizable on a physical system, which we denote as $\Sigma$. For example, we limit the total transmission power by restricting the set of apodization weights $\{w_{MT}\}$ to an $\ell^{\infty}$ ball, effectively clipping each value to the range $[-1,1]$. Other restrictions stem from quantization of numerical values to match machine clock cycles. These small discretizations are much more dependent on the specific acquisition system, but also less impactful to the optimized result.

Altogether, we are left to solve

Numerically, we implement and train the ML model that minimizes Equation 10 in PyTorch, a Python package for building deep learning and other machine learning models (Paszke et al., 2019). This platform is flexible to our unique architecture, which lacks the abstract layers of a neural network and instead propagates information only through acoustically motivated operations, namely the encoding/decoding of the multistatic data and the generation of a B-mode image with DAS beamforming. In all other respects, we mirror the conventional ML pipeline by minimizing our loss function over batches of training data. We perform the numerical optimization using the Adam algorithm, an adaptive version of standard stochastic gradient descent (Kingma and Ba, 2017), and accommodate the constraint set $\Sigma$ for physically realizable encoding sequences with a simple projection of our sequence during each update. The update loop of our ML model is depicted in Figure 1, for which we emphasize the portion that is deployed in an experimental context.

To effectively update the parameters of any ML model requires knowledge of first order derivatives, which we obtain using the reverse-mode automatic differentiation, or backpropagation, capabilities of PyTorch (Griewank and Walther, 2008). In the context of ultrasound imaging, standard techniques are often sufficient when the loss function depends primarily on multistatic data (Chen et al., 2021) or on individual pieces of RF data (Noda et al., 2020). This is because platforms like PyTorch have built-in analytic derivatives for operations that are common in deep learning applications, such as matrix multiplication and Fourier transforms. Gradients for all other operations can be constructed at runtime at the cost of additional computation time and memory overhead. However, this process is made difficult when an ML model directly incorporates beamforming, which we consider a critical component to produce high-quality encoding sequences (See Section 4.4). Backpropagation through the DAS beamformer rapidly becomes a computational bottleneck when it is implemented directly as a sum of linear interpolations of IQ data into per-pixel focused channel data, as the number of intermediate values that must be stored in memory grows rapidly. This is especially relevant while training the proposed ML model, as we must simultaneously perform this expensive calculation over entire batches of training data.

In our application, we circumvent this issue through a novel PyTorch implementation of the derivative of the DAS beamforming operator $\mathcal{B}$. Because the beamformed image $\mathcal{I}$ is linear as a function of the multistatic data $\mathcal{U}$ (Perrot et al., 2021), the analytic derivative needed for backpropagation is simply its adjoint operator. We provide PyTorch the DAS adjoint for arbitrary data, which eliminates the additional cost and memory overhead.

We do not directly construct and transpose the DAS beamforming matrix, as it is impractically large for reasonably sized multistatic data, even under a sparse representation (Perrot et al., 2021). Instead, we apply this adjoint matrix-free, recognizing that $\mathcal{B}$ can be described as a sum of linear interpolations, or

From this, we derive from simple analytic formulas that

These two component adjoints have known, albeit obscure formulas (Claerbout, 2014), and so we provide pseudocode for both operations for clarity. The code used to perform the optimization, along with the RF data that support the findings of this study, are available at github.com/jcs15c/optimal_ultrasound_encoding (Spainhour et al., 2024).

While training on true in vivo data would best reflect the experimental setting, it is understood that large collections of general purpose, labeled clinical data are rare (Liu et al., 2019). Instead, we utilize simulated multistatic data created with Field II (Jensen, 2004) to train and evaluate our proposed ML model. The simulated transducer configuration is described in Table 1. We have found that multistatic data depicting randomly placed anechoic lesions in an underdeveloped speckle pattern is a particularly effective class of training data, as encoding sequences trained on it generalize well to other classes of data (See Section 4.2).

For training, we generate with Field II a collection of 500 instances of underdeveloped speckle data, of which 20% form a validation set for by-hand tuning of our (reasonably few) hyperparameters. Notably, the proposed model is exposed only to this type of data during training. For thorough testing of the optimized encoding sequence, we use additional samples of underdeveloped speckle data, as well as collections of other types of simulated data. These include conventional imaging targets, namely isolated point scatterers and anechoic lesions in fully developed background speckle. These data are also used to evaluate the optimized encoding sequence according to standard ultrasound imaging quality metrics, such as the cystic resolution and gCNR respectively.

To explore performance on more complex RF data, we follow a procedure similar to Hyun et al. (2019) and generate arrangements of scatterers with amplitudes weighted according a grayscale image. These images are derived from a set of 160 samples from the publicly available validation set of an ImageNet competition (Howard et al., 2018), where we have extracted and smoothed the middle 256$\times$256 pixels from each and converted them to grayscale. The first type of image-derived data interpolates the pixel values of each grayscale image to nearby scatterers that are distributed throughout the spatial domain, allowing us to consider data with scatterers of arbitrary amplitude. The second type of image-derived completely removes scatterers near pixels of sufficiently low brightness, resulting in anechoic regions in background speckle that are more arbitrary than standard circular cysts.

We also test our optimized encoding sequences on experimental hardware with both a tissue-mimicking phantom and a wire target, as described in Section 3.4.

We must similarly consider the imaging target used by the loss function in Equation 8 during the supervised learning. A natural choice is an “unencoded” image, where each simulated ground-truth multistatic data set is directly focused with DAS beamforming, thereby removing artifacts induced by encoding. Alternatively, one can compare to a more “idealized” target that represents the ground-truth contrast, as this also reduces the influence of any particular speckle pattern during training. Each image of ground-truth contrast is generated alongside the associated simulated multistatic data, where spatial masks are used to create homogeneous anechoic and scattering regions. An alternative approach to describing ground-truth echogenicity would be to directly convolve a high-quality point spread function (PSF) with each scatterer (Goudarzi and Rivaz, 2022), but we find our approach using a spatial mask to be simpler computationally. For image-derived data, the ground-truth contrast is simply the original grayscale image. During training, we use ground-truth contrast targets, and discuss the consequences of this decision in Section 4.3.

In Figure 2 we show an example of each target type for each category of simulated data, highlighting the particular configuration of data that is exclusively used during training.

We now demonstrate our capability to effectively train the proposed ML model, thereby producing encoding sequences that meaningfully improve image quality. In the most general case, the available hardware allows us to manipulate both time delays and apodization weights for each array element. This allows the greatest freedom when designing encoding sequences, and predictably the greatest improvement over conventional sequences. As a prototypical configuration of the training workflow, we train the ML model using 400 simulated instances of the previously described underdeveloped speckle background with anechoic lesions, processed in batches of size 8. We evaluate the model with the $\ell^{2}$ loss function in Equation 8 using ground-truth contrast as the imaging target. The training is performed using the Adam optimizer using a learning rate of 0.1 over 25 epochs. The available Verasonics hardware is limited to a minimum duty cycle of two clock cycles, which effectively nullifies low values for our weights. As a result, we appropriately restrict the parameter space of encoding sequences $\Sigma$ during training (See Section 2.3).

While this is one of many possible configurations of the proposed training workflow (e.g., the choice loss function can vary according to the specific imaging task) we will see that the optimized encoding sequence created in the above configuration improves image quality fairly broadly across a number of different metrics.

We first demonstrate these improvements by evaluating the optimized encoding sequence on several types of simulated ultrasound data in Figure 3, for which we uniformly observe significantly improved contrast and reduction of scattering artifacts. For comparison, we use two conventional planewave encodings of varying maximum angle extent, as such sequences characterize a common tradeoff between FOV and resolution at depth (Bottenus et al., 2023). For example, the first generates 15 beams with 1 degree of separation, which results in a higher resolution along a narrower FOV. This is particularly notable when imaging point targets, where some scatterers are essentially undetectable. Steering these 15 planewaves sufficiently far apart to cover the entire imaging domain results in shallower contrast for anechoic regions, as there is more significant scattering throughout the domain.

By contrast, the proposed ML model has generated a sequence that recovers this wider FOV simultaneously with improved contrast. Importantly, we observe through the improvements on the image-derived data that the optimized encoding sequence is performant regardless of the specific arrangement of scatterers in the field. Instead, the point spread function for each scatterer is improved more generally and consistently throughout the FOV.

We show these improvements quantitatively in Table 2, where this prototypical optimized encoding sequence leads to improved $\ell^{2}$ error against the ground-truth contrast. This is the same metric that is minimized over the training set of anechoic lesions in underdeveloped speckle (See Equation 8), and it is therefore natural that our ML model improves this metric for this type of data. However, we see that these improvements in the $\ell^{2}$ loss are consistent across several other types of imaging targets that are visually quite distinct.

We can also demonstrate that improvements in this $\ell^{2}$ loss are strongly associated with improvements in more conventional image quality metrics, namely the gCNR and cystic resolution, thereby justifying the decision to select an $\ell^{2}$ loss for training of the ML model.

The cystic resolution quantifies detectability of an anechoic cyst in background speckle with a given level of contrast (Ranganathan and Walker, 2007). We compute this metric through cystic contrast, which is the concentration of energy within a certain distance from a point target in an anechoic field. When cystic contrast is considered as a function of target radius, the cystic resolution is the minimum radius that achieves a desired level of contrast. For our examples we use a contrast of -20 dB, a common threshold for lesion detectability that represents an order of magnitude difference from the surrounding speckle. By explicitly plotting the cystic contrast for each encoding sequence in Figure 4(a), we can see that the optimized encoding sequence produces lesion detectability on par with narrowly concentrated planewaves, and that this conclusion is not sensitive to the exact choice of threshold.

In contrast to these narrow planewaves, however, our optimal encoding sequence maintains this level of resolution throughout the imaging domain. To see this, we use Field II to simulate the responses from single, isolated scatterers distributed throughout the viewing window, and compute the cystic resolution for each when imaged with the optimized sequence. Shown in Figure 5, this procedure generates a map of lesion detectability as a function of target position. This strongly suggests that the benefits of our optimized encoding sequence stem from more general acoustic principles, rather than overfitting to the training data, as the ML model is never exposed to these isolated point targets during training.

We additionally measure contrast using the gCNR, which also varies according to the spatial position of the lesion target. To account for this variance, we consider simulated anechoic lesion data for which the position and radius of each lesion is random. We categorize these lesions based on position, and tabulate the average gCNR for each category in Figure 5. To consider more complicated lesion shapes, we also measure the gCNR of the binarized, image-derived data of Figure 2(b), which we record in Table 3. Altogether, we see that our optimized encoding sequence outperforms conventional planewave and Hadamard encodings in terms of both FOV and contrast.

We also compare our optimized sequence to standard techniques in apodization, which similarly seek to improve the PSF by reducing side lobes (Jeong and Kwon, 2020). To do this, we consider planewaves with $10^{\circ}$ spacing, which offers good resolution at depth while imaging the full FOV, and apply a Tukey window to the receive channel. As expected, this improves cyst detectability in the image considerably. However, it does so at the cost of resolution, as shown in Figure 6. Specifically, there are still clear artifacts around point scatterers, and the striation in the cystic resolution map for the unapodized case is still present after apodization. These effects are not present in our optimized encoding sequence.

We can also apply the proposed ML model and training framework to more restrictive hardware specifications. In doing so, we demonstrate that there is still considerable improvement present, although such sequences do not perform as well as a fully arbitrary system, as expected. To this end, we train two separate configurations of the ML model. In the first, the apodization weights are fixed to some uniform value, and the time delays are the only trainable parameters of the ML model. In the second, only the apodization weights can be changed, and each time delay is fixed to zero in the ML model. As in Section 3.1, we train the model on anechoic lesions in underdeveloped speckle and test the resulting sequence on data far outside the training set. We compare the two resulting optimized sequences to the appropriate conventional alternative: wide spanning planewaves for delay-only optimization, and a truncated Hadamard encoding for weight-only optimization.

The results of this comparison are compiled in Table 4. We see that in each case the optimized sequence performs better according to the $\ell^{2}$ loss metric than its conventional alternative, which we have now shown correlates with high quality as measured by gCNR and cystic resolution. Furthermore, we observe that while there is still considerable improvement when only delays are optimized through the training of our ML model, the magnitude of improvements when only optimized weights are generated is comparable to that of a fully optimized encoding sequence. This demonstrates that in this imaging scenario, optimization of the weights has a more significant influence on the quality of the resulting sequence than delay optimization alone.

These simulated results demonstrate the theoretical utility of applying an ML optimization strategy for data acquisition and imaging within the REFoCUS framework. However, we can also verify that our optimized sequences produce meaningful improvements in several practical scenarios. In the first, we consider a separate configuration of the ML model that is trained to emulate an imaging system under the influence of electronic noise. We recreate this during training by the addition of random noise to the encoded channel data prior to imaging. In the current example, this additive noise is sampled from a random normal distribution, filtered to have a bandwidth equal to that of the transducer, and then scaled to produce an average of 20 dB channel signal-to-noise ratio. It is in this context that restricting the set of apodization weights to an $\ell^{\infty}$ ball is most critical, as otherwise all noise would be suppressed simply by arbitrarily increasing the apodization weights.

By accounting for this noise during training, the ML model generates a new encoding sequence that specifically suppresses this noise. We see this in comparison to the optimized encoding sequence evaluated throughout Section 3.1. In Figure 7(a), we consider these two sequences, trained with and without the presence of noise, along with a conventional planewave configuration for comparison. We then apply each sequence to data to which noise is, or is not added. In each of the four cases, we see considerable improvement over the conventional planewave approach. However, we can also see that these improvements are most prominent when the ML model is trained with the same type of noise present during evaluation, the case emphasized in Figure 7(a) in red. In Figure 7(b), we consider a more complete evaluation of the case where noisy data is imaged using a sequence trained with the same noise model. There is the same degree of improvement as the noiseless case of Figure 3, featuring dramatic improvements in the gCNR and cystic resolution.

Finally, we verify these results by directly implementing the optimized encoding sequence discussed in Section 3.1 in physical hardware. In the following experiments, we use the Verasonics Vantage 256 research scanner (Verasonics, Inc., Kirkland, WA) with the P4-2v phased array transducer (3 MHz transmit frequency, 64 elements, 0.3 mm pitch, 10V transmit voltage). Received channel data were stored for processing offline, where they were decoded and used to produce an image in the same manner as in simulation. For comparison, we use sequences of 15 transmits that generate planewaves of varying span (15, 60, and 150 degrees). To reproduce the effects of arbitrarily scaled weights (rather than restricting to $\pm 1$ through pulse inversion) we use the Verasonics apodization function to scale the duty cycle of the transmit excitation so that it matches the output amplitude.

The first experiment utilizes an ATS 539 multi-purpose phantom (CIRS, Inc., Norfolk, VA) to confirm the ability of optimized encoding sequences to improve imaging on tissue-like material. Within the phantom we image anechoic lesions of varying position and radius in background speckle, and compute the gCNR of each with respect to surrounding background speckle. We show these gCNR values below each lesion in Figure 8(a). For a fixed lateral position, we acquire data using each encoding sequence at each of seven elevation positions of the phantom. In this way, we image anechoic lesions in the same positions with different realizations of background speckle, and compute the average gCNR across all realizations.

We can see that these results match closely with our numerical simulations, in that our optimized sequence is able to produce the depth of focus and resolution of a very narrow span of planewaves, while still maintaining a full FOV throughout the range. Importantly, these improvements persist beyond the fixed viewing window of the training data, which extends only to a depth of 55 mm. This further emphasizes that our encoding sequence has not simply specialized to the specific characteristics of the training data.

We can also recreate our simulated point target images, which we can in turn use to measure the cystic resolution throughout the imaging domain (albeit on a much coarser grid). To accomplish this, we image a custom, single target wire phantom (0.03 mm tungsten wire) in a water tank. Using the AIMS III Hydrophone Scanning System (Onda Corporation, Sunnyvale, CA), we mechanically translate the suspended P4-2v transducer around the fixed wire so that the target can be imaged from a 6$\times$5 grid of target positions ($\sim$10 mm lateral spacing, $\sim$8 mm axial spacing).

Although we only ever image one point target at a time with this experimental setup, in Figure 8(b), we present a montage of these images distributed appropriately in space. To ensure a fair visual comparison of each point target, we apply a uniform level of gain to each image in the domain. This allows us to very clearly see the effect of our optimized sequence on the FOV of the imaging system, where points in the upper corners of the viewing window are either undetectable or suffer from very poor resolution. In contrast, the optimized encoding sequence provides a uniform degree of resolution at all points in the domain.

Furthermore, we can see that the off-axis scattering artifacts present are greatly reduced when imaged with the optimized sequence. As before, we quantify this effect using the cystic resolution, which we can compute for each point target in Figure 8(b). Although our ability to perfectly capture this metric is slightly hampered by irregularities in the experimental setup, we still observe near uniform improvement over each conventional sequence, analogous to results observed in simulation.

Because the minimization problem presented in Equation 10 is not convex, there is no expectation that any optimization procedure can find the global minimum, should one even exist. Instead, the particular local minimum approached by our optimization procedure is highly dependent on the initial condition. Because locally optimal encoding sequences vary in quality, it is important to select an initial condition that is at or near a good local minimum. The interpretability of the parameters of our ML model offers a distinct advantage over conventional deep learning applications in accomplishing this, as we can intelligently explore the effects of different initial conditions.

For example, an instance of our ML model initialized with a planewave encoding will converge to transmissions that are steered in the same way, albeit with less well-defined wavefronts. On the other hand, initializing with a truncated Hadamard code causes the parameters to converge to a sequence with no visible steering. We can see these qualitative behaviors in Figure 9(a). Despite this, we have found experimentally that a truncated Hadamard sequence is the superior initial condition, despite being a generally lower quality transmission sequence. This can be seen through the $\ell^{2}$ loss plotted over the course of training epochs for each initial condition in Figure 9(b). This can be explained by the conventional Hadamard sequence allowing the ML model to begin with near-optimal FOV, as well as both positive and negative apodization weights. Except for cases in which delays or weights must be fixed as in Section 3.2, each optimized encoding sequence used in this paper was generated by an ML model with this configuration as the initial condition.

In our exploration of this topic, we have found that recovering a known, conventional sequence is exceptionally difficult except in the simplest toy problems. For example, although planewave transmissions are generally performant, the ML model cannot identifiably recover the same geometric properties that make them desirable. On the other hand, the optimal characteristics of the resulting sequence are quite opaque: while optimized sequences appear to be a simple perturbation of the initial condition, performing this perturbation randomly fails to achieve the same improvements. This suggests that it would not be possible to discover these sequences without the use of the proposed ML model.

In Figure 9(b), we also observe that for either initialization, our ML model performs nearly identically between our training and testing data. This is a highly desired property in this context, as ultrasound imaging can be theoretically formulated as a convolution of a spatially varying point spread function across each scatterer in the domain, the results of which are then summed to form an image. Therefore, an effective manner of optimizing image quality for arbitrary data would be to optimize this PSF itself, concentrating it at zero offset (Goudarzi and Rivaz, 2022).

Although our loss function does not reference features of the PSF directly, our choice to use underdeveloped speckle as training data implicitly encourages the same type of improvement at a lesser computational burden. This is in part because the scatterers within each sample are spread out enough that their off-axis scattering artifacts do not necessarily overlap one another, as would be the case with standard background speckle. This means the $\ell^{2}$ loss function can only be meaningfully decreased by acoustically concentrating the energy of each scatterer inwards, thereby improving the PSF. At the same time, scatterers across all samples in the training data are densely distributed throughout the imaging domain. This means that learning improvements that span the entire FOV can be obtained with relatively few samples of type of training data.

Although we exclusively use a loss function that compares produced B-mode images with some ground-truth imaging target, there is still much flexibility within that choice. As stated in Section 2.4, a natural choice for the imaging target is an image derived from ground-truth multistatic data. While this produces target images without any encoding artifacts, we have found that an ML model trained with these targets produces poor image quality relative to alternatives. This is because DAS beamforming with even unencoded multistatic data results in a PSF with undesirable features such as the visible side lobe artifacts present throughout the domain. When the ML model is trained on such target images, the transmit sequence is implicitly encouraged to accurately reconstruct these artifacts. Instead, we see that training against the exact ground-truth contrast produces images with a higher quality PSF, as seen in Figure 10.

Additionally, we see the most uniform improvement in image quality when we “pad” images of underdeveloped speckle so that the final viewing window is itself located in an anechoic field. This ensures that the encoding sequence is not unnecessarily dependent on the specific domain geometry. The consequences of this are best seen through the plots of cystic resolution in Figure 11, where the first encoding sequence is trained on images whose extent is exactly that of the final viewing window. In contrast, the second encoding sequence is trained on images located in a slightly wider anechoic void. As we can see, adding this empty space causes the PSF to more smoothly vary throughout the domain, and failing to do so results in additional artifacts for point targets near the center of the region. In effect, these streaks depict positions for which scattering artifacts extend beyond the imaging domain, and the introduced error is not captured by the loss function.

By using our custom implementation of a differentiable beamformer, we are able to train our ML model according to improvements in an imaging metric, which is a novel departure from existing work that only considers the recovery of the multistatic data set. Importantly, we can demonstrate the existence of transmit sequences that produce images with higher resolution and contrast, even though the encoding itself produces a worse reconstruction of $\boldsymbol{\mathcal{U}}$ from $\boldsymbol{\mathcal{S}}$.

This is primarily because a perfect reconstruction of $\boldsymbol{\mathcal{U}}$ will only produce the original unencoded image, and this type of image lacks many desirable properties, as discussed in Section 4.3. To create a concrete example of this phenomenon, we compare to an ML model within our framework that is instead configured to directly optimize the recovery of the multistatic data set by minimizing

Consider Figure 12. We see that this new configuration successfully generates an encoding sequence that produces a lower STA loss than the previously considered encoding sequence optimized for image recovery. Yet despite this numerical improvement, there is a clear superiority in image quality when the optimization process considers image formation, leading to improvements both visually and in terms of the gCNR averaged over each lesion.

We take this as a counterexample to the common convention that a better conditioned encoding necessarily corresponds to higher quality images. To account for the effects of Tikhonov regularization, we measure the condition number of the collection of encoding matrices $\boldsymbol{\mathcal{H}}$ as $\kappa=\|\boldsymbol{\mathcal{H}}_{0}\|\cdot\|\boldsymbol{\mathcal{H}}^{\dagger}_{0}\|$, doing so because the encoding matrix for the lowest frequency $\boldsymbol{\mathcal{H}}_{0}$ has the poorest conditioning among the collection $\boldsymbol{\mathcal{H}}$. We see that both encoding sequences bring about reasonably well-conditioned matrices following Tikhonov regularization, but there are still clear qualitative differences between the two images.