Deep learning for sound speed imaging using pulse-echo ultrasound

Abstract

Provided are methods of creating a trained deep neural network, in embodiments the methods comprising: providing a number of natural images; converting each natural image into a speed of sound (SoS) map; for each of the SoS maps, deriving a theoretical time-shift map; and training a deep neural network based on time-shift and SoS map pairs, each of the pairs comprising one of the theoretical time-shift maps and its corresponding SoS map; wherein the deep neural network is capable of reconstructing a speed of sound (SoS) image from a measured time-shift map generated from ultrasound.

Description

BACKGROUND OF THE INVENTIONField of the Invention

[0001] The present invention is directed to the field of deep learning networks and, in particular, deep neural networks for sound speed imaging. Embodiments provide methods of creating a trained deep neural network by providing a number of natural images; converting each natural image into a speed of sound (SoS) map; for each of the SoS maps, deriving a theoretical time-shift map; and training a deep neural network based on time-shift and SoS map pairs, each of the pairs comprising one of the theoretical time-shift maps and its corresponding SoS map; wherein the deep neural network is capable of reconstructing a speed of sound (SoS) image from a measured time-shift map generated from ultrasound.Description of Related Art

[0002] Ultrasound is a widely available, safe, and portable medical imaging modality, but its usefulness in differential diagnosis (e.g., differentiating between benign and malignant tumors) is limited by unspecific contrast for certain disease types. Traditional B-mode ultrasound imaging is based on echoes backscattered from interfaces of two media that have different acoustic impedance values. Novel contrast mechanisms can be achieved by imaging biomechanical properties that change with disease processes, with sound speed being a promising candidate. Many of the current clinical scanners provide capabilities to image the shear wave speed and have achieved tremendous success in liver fibrosis assessment and tumor classification, among other applications.

[0003] Despite these successes, shear-wave imaging has many known limitations. Longitudinal sound speed (typically ranging between 1400 and 1650 meters per second in soft tissue), a different type of sound speed, has much higher values than shear wave speed (typically a few meters per second) and represents a different biomechanical property. The longitudinal sound speed, hereinafter called the sound speed, is shown to differ between different tissue types, differ between normal versus diseased tissues, and correlated with the amount of fat in the liver. Therefore, the sound speed has promising applications for tissue characterization and differential diagnosis and may provide capabilities that are currently unavailable from B-mode images. Moreover, knowing the speed of sound distribution can improve the quality of B-mode images.B-mode image quality is negatively affected by the so-called aberration correction. The aberration can be more accurately corrected if the speed sound distribution in tissue is known.

[0004] Several methods exist in the literature to image the speed of sound in the human body. Each method has its own limitations, summarized as follows.

[0005] Ultrasound Computed Tomography is based on the through-transmission technique, that is, transmitting the ultrasound waves from one side of the body / organ, and detecting the sound waves from the other side. Therefore, this method requires specially designed ring-shaped ultrasound transducers, and requires that the organ being scanned fit into the ultrasound. The method cannot be used on widely available clinical ultrasound scanners that work in the pulse- echo mode and cannot be used to image most parts of the body (e.g., abdominal image).

[0006] Physics-based pulse-echo sound speed imaging works in pulse-echo mode, without requiring specifically designed ring-shaped ultrasound transducers and can be used to image most parts of the body. The method consists of two major steps. The first step includes transmitting sound waves from slightly different angles, receiving the sound waves with the same transducers, and computing the time-shift map (that is, difference between the echo receive time) between different transmission angles. The second step is to use the time-shift map to reconstruct the sound speed map by solving an inverse problem. This second step is very time-consuming and cannot be performed in real time, however, real-time imaging is essential for ultrasound. Also, the accuracy of the method is limited by the simplified assumptions of the physical models.

[0007] Supported by the universal approximation theorem, deep learning (DL) has the potential to accurately reconstruct the SoS, given 1) substantial training data and 2) consistency between training and testing domains. Moreover, the DL method is fast in the inference stage and fulfills the real-time imaging requirement. Driven by these promises, researchers developed frameworks to reconstruct SoS from the received channel data.

[0008] While DL methods demonstrate high accuracy in controlled simulation environments, their performance significantly declines when applied to real-world data. This discrepancy underscores a generalizability challenge when transitioning from simulations to experimental data. A root cause of the challenge is the complexity of simulating imaged media, where acoustic properties vary across tissues and remain hard to accurately anticipate. Prevalent as they are, unrealistic assumptions on acoustic properties, e.g., attenuation and acoustic non-linearity, by theory, can introduce considerable differences in amplitude and reverse local phases of pulse echoes. These differences, consequently, negatively impact neural networks that directly use pulse echoes as input.

[0009] Instead of using raw channel data, this invention explores the idea of using domain- specific time-shift measurements as network input. This idea is in light of the physical principle that the spatial variation of SoS causes shifts in travel time. This principle has recently been used to develop non-DL approaches for SoS imaging. Herein the inventor proposes a DL framework that reconstructs SoS distribution from time-shift maps. Specifically, the problem is formulated as an image-to-image translation task, whose solution is learned using a convolutional neural network. Compared with physics-based non-DL reconstruction approaches that formulate a linear forward model, this formulation is expected to achieve more accurate reconstruction results (as the network may be able to implicitly model the physical factors that are not considered in an explicit non-DL inversion model) and a faster reconstruction speed (due to the low computational burden during model inference). From the DL perspective, the use of time-shift maps as the input compared with using the raw channel data as the input has three perceived benefits: 1) improving the generalizability across simulation settings; 2) enabling the use of the straight-ray forward model for efficient data synthesis; and eventually 3) improving the generalizability from simulation to experimental data.

[00010] Under this overall idea, this invention proposes a nonlinear SoS inversion approach that employs deep neural network to learn from theoretical and simulation data without experimental data for training, wherein time-shift maps are used as the input to the DL model. To maximize the information provided to the DL model, the inventor designs an approach that uses multiple unique time-shift maps as the model input. Specifically, common mid-angle beamforming was performed at multiple mid-angles, each yielding a unique time-shift map. Also, the use of common mid-angle beamforming enhances the quality of time-shift maps.

[00011] Some techniques include those disclosed in US Patent Nos. 10,799,208; 10,271,819; and 11,397,167; US Published Patent Application No. 2021 / 0177380; International Patent Application Publication No. WO 2020 / 097298, and in journal articles (such as M. Bernhardt, V. Vishnevskiy, R. Rau, and O. Goksel, “Training Variational Networks with Multidomain Simulations: Speed-of-Sound Image Reconstruction”, IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol 67, no 12, 2020, doi: 10.1109 / TUFFC.2020.3010186; Ruby et al. "Breast cancer assessment with pulse-echo speed of sound ultrasound from intrinsic tissue reflections: Proof-of-concept." Investigative radiology 54.7 (2019): 419-427), which references are hereby incorporated by reference herein in their entireties.SUMMARY OF THE INVENTION

[00012] Aspects of embodiments of the invention include Aspect 1, which is a method of creating a trained deep neural network, the method comprising: providing one or more natural image(s); converting each of the one or more natural image(s) into a speed of sound (SoS) map; for each of the SoS maps, deriving a theoretical time-shift map; training a deep neural network using time-shift and SoS map pairs, wherein each of the pairs comprises one of the SoS maps and the theoretical time-shift map derived therefrom; and optionally, fine-tuning the deep neural network using simulated time-shift and SoS map pairs, each of the pairs comprising a SoS map and its corresponding time-shift maps obtained through acoustic simulations, e.g., full-wave simulations; wherein the deep neural network is capable of reconstructing a speed of sound (SoS) image from a measured time-shift map generated from ultrasound.

[00013] Aspect 2 is an imaging method comprising: imaging an object using ultrasound; converting, e.g., in real-time, ultrasound signals obtained from the imaging into one or more measured time-shift map(s); inputting the measured time-shift map(s) into the deep neural network of Aspect 1; and obtaining as an output a speed of sound (SoS) image reconstructed from the measured time-shift map(s).

[00014] Aspect 3 is an imaging method, e.g., a real-time imaging method, comprising: imaging an object using ultrasound; converting ultrasound signals obtained from the imaging into one or more measured time-shift map(s); providing a speed of sound (SoS) image by analyzing the measured time-shift map with a deep neural network.

[00015] Aspect 4 is a device comprising: memory storing instructions of a program; and a processor that generates a speed of sound (SoS) image; wherein the processor is capable of executing the instructions to reconstruct the SoS image from a measured time-shift map derived from ultrasound signals; and wherein the reconstruction of the SoS image is provided in response to input of the measured time-shift map(s) into a trained deep neural network.

[00016] Aspect 5 is the method or device of any of Aspects 1-4, wherein the deep neural network is not trained on: actual / experimental ultrasound data / signals; and / or geometrical information of the object to be imaged.

[00017] Aspect 6 is the method or device of any of Aspects 1-5, wherein one or more or each of the theoretical and / or simulated time-shift maps account for receive (Rx) delay.

[00018] Aspect 7 is the method or device of any of Aspects 1-6, wherein the measured time-shift map accounts for Rx delay.

[00019] Aspect 8 is the method or device of any of Aspects 1-7, wherein one or more or each of the theoretical and / or simulated time-shift maps is derived using dynamic programming and / or vectorization.

[00020] Aspect 9 is the method or device of any of Aspects 1-8, wherein the measured time-shift map is derived using dynamic programming and / or vectorization.

[00021] Aspect 10 is the method or device of any of Aspects 1-9, wherein the deep neural network is trained by: converting a representative time-shift map into an SoS map; inputting the SoS map into a theoretical time-shift map model; obtaining a predicted time-shift map for the SoS map; and comparing and minimizing any difference between the predicted time-shift map and the representative time-shift map.

[00022] Aspect 11 is the method or device of any of Aspects 1-10, wherein the deep neural network is trained by data parallelism involving multiple GPUs.

[00023] Aspect 12 is the method or device of any of Aspects 1-11, wherein the measured time- shift map is obtained in real time from ultrasound imaging of a patient.

[00024] Aspect 13 is the method or device of any of Aspects 1-12, wherein the training is performed in a manner, such as to a degree, to obtain an ultrasound-based system-independent, e.g., generalizable, trained deep neural network.

[00025] Aspect 14 is the method or device of any of Aspects 1-13, wherein the reconstructed SoS image is based on the inputting into the trained deep neural network of: only a single measured time-shift map or multiple measured time-shift maps; and no geometric information of the object to be imaged.

[00026] Aspect 15 is the method or device of any of Aspects 1-14, wherein no B-mode image is used as guidance and / or the deep neural network is not trained on experimental / actual ultrasound data / signals.

[00027] Aspect 16 is the method or device of any of Aspects 1-15, wherein: no ultrasound data / signals are used in the training of the deep neural network; i.e., only the theoretical ultrasound data are used.

[00028] Aspect 17 is the method or device of any of Aspects 1-16, wherein the deep neural network has a U-net architecture.

[00029] Aspect 18 is the method or device of any of Aspects 1-17, wherein the deep neural network has a GAN architecture.

[00030] Aspect 19 is the method or device of any of Aspects 1-18, wherein data provided for the deep neural network is provided in a radiofrequency (RF) format.

[00031] Aspect 20 is the method or device of any of Aspects 1-18, wherein data provided for the deep neural network is provided in an in-phase and quadrature (IQ) format.

[00032] Aspect 21 is the method or device of any of Aspects 1-20, wherein the object or patient comprises a body, or part thereof, such as a human or animal body, including tissue and / or one or more organs thereof.

[00033] Aspect 22 is the method or device of any of Aspects 1-21, wherein the theoretical and / or the measured time-shift maps between adjacent steering angles are computed by using RF or IQ data from adjacent steering angles.

[00034] Aspect 23 is the method or device of any of Aspects 1-22, wherein the time-shift maps are concatenated to make a single robust time-shift map.

[00035] Aspect 24 is the method or device of any of Aspects 1-23, comprising: determining if an SoS image needs to be learned by estimating the theoretical time-shift map, instead of performing a full-wave simulation, using a physically-derived forward model; and identifying and selecting therefrom additional training data.

[00036] Aspect 25 is the method or device of any of Aspects 1-24, comprising selectively running simulations using the additional training data selected.

[00037] Aspect 26 is the method or device of any of Aspects 1-25, wherein the measured time-shift maps are derived from data obtained from transmitting at multiple steering angles.

[00038] Aspect 27 is the method or device of any of Aspects 1-26, wherein the measured time-shift maps are derived from data obtained from different sized apertures.

[00039] Aspect 28 is the method or device of any of Aspects 1-27, wherein the Rx path is incorporated into the forward model and / or the beamforming process.

[00040] Aspect 29 is the method or device of any of Aspects 1-28, comprising using mirrored receive aperture, wherein the transmit (Tx) path is at Tx angle @Tx, the Rx path is at an angle QRX, and the Rx angle Rx is kept the same as the Tx angle QTx.

[00041] Aspect 30 is the method or device of any of Aspects 1-29, wherein in the forward model, an integration term is added representing Rx delay along the Rx path:√1 + tan² QRX Z ∫ 0 dz'σ(x + ztanqRx – z'tanqRx, Z')

[00042] Aspect 31 is the method or device of any of Aspects 1-30, wherein the deep neural network comprises a physics-based deep learning training framework, such as is shown in FIG. 2.

[00043] Aspect 32 is the method or device of any of aspects 1-31, further comprising determining a total time of flight as t(r) = ∫p 1 c(rp) dl .

[00044] Aspect 33 is the method or device of any of claims 1-32, further comprising determining a time delay introduced by beamforming as At(r) = ∫p ( 1 c(rp) 1 Cbf. ) dl.

[00045] Aspect 34 is the method or device of any of claims 1-33, further comprising determining a phase shift as (r) = 2πfc [ ∫P1 ( 1 cP1(r) 1 cbf ) dl cos (PTX1-PRX1) 2 ∫P2 ( 1 cp2(r1) 1 cbf ) dl cos (@TX2-PRX2) 2 ]

[00046] Aspect 35 is the method or device of any of claims 1-34, further comprising determining a time shift as (r) = 1 2πfc (r).

[00047] Aspect 36 is the method or device of any of claims 1-35, further comprising expressing time shift as t = L.s.

[00048] Aspect 37 is the method or device of any of claims 1-36, further comprising performing beamforming using a common mid-angle approach with a pre-specified mid-angle defined as Omid = PTX+PRX 2 .

[00049] Aspect 38 is the method or device of any of claims 1-37, further comprising determining a phase shift based on cross-correlation for a pair of beamformed RF frames according to:R(PTX1, PRX1RX2, PRX2. PRX2) = LPF(bRF(QTX₁, PRX₁) ⊗ bRF* (QTX2, PRX2)).BRIEF DESCRIPTION OF THE DRAWINGS

[00050] The accompanying drawings illustrate certain aspects of implementations of the present disclosure and should not be construed as limiting. Together with the written description the drawings serve to explain certain principles of the disclosure.

[00051] FIG. 1 is a graph showing transmission with different steering angles.

[00052] FIGS. 2A-D are a theoretical time-shift map (FIG. 2A), a measured time-shift map (FIG. 2B), the reconstructed SoS map (FIG. 2C), and the ground-truth SoS map (FIG. 2D) generated in a DL model without considering Rx delay.

[00053] FIGS. 3A-B are images showing a benefit of using the mirrored receive aperture technique in the models according to embodiments of the invention.

[00054] FIG. 4 is an illustration of Tx and Rx.

[00055] FIGS. 5A-C are a theoretical time-shift map (FIG. 5A), a measured time-shift map (FIG. 5B), and the reconstructed SoS map (FIG. 5C) using a DL model with Rx delay taken into account.

[00056] FIG. 6 is a schematic diagram of the deep learning pipeline according to an embodiment of the invention.

[00057] FIG. 7 compares estimated and ground-truth SoS maps for a training set and testing set according to an embodiment of the invention.

[00058] FIG. 8 is a schematic diagram comparing image reconstruction using the traditional analytical approach versus deep learning according to embodiments of the invention.

[00059] FIG. 9 is a graph showing the estimation accuracy among different SoS groups.

[00060] FIG. 10 is a series of ground-truth and reconstructed SoS maps according to embodiments of the invention.

[00061] FIG. 11 shows conversion from a natural image to a theoretical time-shift map.

[00062] FIG. 12 is a training set including theoretical time-shift maps, reconstructed SoS maps, and ground-truth SoS maps according to an embodiment of the invention.

[00063] FIG. 13 is a training set including theoretical time-shift maps, reconstructed SoS maps, and ground-truth SoS maps according to an embodiment of the invention.

[00064] FIG. 14 is an SoS imaging approach according to an embodiment of the invention.

[00065] FIG. 15 is an illustration showing the calculation of an infinitesimal time delay associated with the round-trip path of a transmission according to an embodiment of the invention.

[00066] FIG. 16 is an illustration of a method to obtain time-shift maps to be used as inputs to the DL model according to an embodiment of the invention.

[00067] FIGS. 17A-B are illustrations showing SoS maps containing ellipses and layered structures, mimicking lesions, and the layered tissue structure of the abdominal wall according to an embodiment of the invention.

[00068] FIG. 18 shows a natural image converted to grayscale and proportionally mapped to the SoS values according to an embodiment of the invention.

[00069] FIGS. 19A-B are drawings showing digital phantoms mimicking tissue structures according to an embodiment of the invention.

[00070] FIG. 20 shows a set of ground-truth and reconstructed SoS maps according to an embodiment of the invention.

[00071] FIG. 21 shows a series of ground-truth and reconstructed SoS maps prepared with three different simulation conditions according to an embodiment of the invention.

[00072] FIG. 22 is a graph showing the estimated error for each of the simulation conditions shown in FIG. 21.

[00073] FIG. 23 is a series of cross-section profiles of reconstructed SoS maps according to embodiments of the invention.

[00074] FIG. 24 provides experimental results for a tissue-mimicking abdominal phantom (top row) and breast phantom (middle and bottom rows).

[00075] FIGS. 25A-D provide exemplary time-shift maps obtained for a breast phantom (FIGS. 25B-D), along with B-mode image for reference (FIG. 25A).DETAILED DESCRIPTION OF VARIOUSEMBODIMENTS OF THE INVENTION

[00076] Embodiments of the present invention include methods and devices for achieving real-time pulse-echo sound speed imaging through deep learning while addressing the limitations of existing deep learning approaches. Exemplary methods are explained as follows:

[00077] 1. Transmit ultrasound waves using an array transducer with different steering angles (i.e., imaging from different directions), and record the echo signals backscattered from the region to be imaged (FIG. 1). The echo signals can be either in the radiofrequency (RF) format or the in phase and quadrature (IQ) format. Transmit waves can be performed using traditional focused beams, plane waves, or divergence waves. Using IQ data is more advantageous than using RF data, because IQ data is more widely available than RF data in clinical scanners, and IQ data requires lower sampling frequencies and less storage than RF data.

[00078] 2. Apply a physical model that considers both the transmit and receive delays for beamforming. Beamforming can be performed using either RF data or IQ data. To ensure generalizability, beamforming should be performed using a prespecified f-number. Note that the f-number has an impact on the obtained time-shift map. (The f-number is tuned by changing the number of receive [Rx] elements.) In theory, the f-number impacts the speckle size and Rx path. The lateral size of speckle cells is determined by 0.87λz / D', where z / D' is effectively the f-number. It indicates that a small f-number leads to a smaller speckle size, which may impact the performance of the speckle-tracking algorithm. In addition, the f-number impacts the Rx path (whether broad or narrow) and therefore the Rx time delay. A smaller f-number empirically leads to more accurate time-shift maps (closer to the theoretical time-shift maps).

[00079] 3. Compute the time-shift map between adjacent steering angles by using the RF or IQ data from adjacent steering angles. Concatenate the individual time-shift maps to make a single robust time-shift map. One way of computing the time-shift map is to use the method described in the following paper (which is incorporated by reference herein in its entirety) to reduce the sampling frequency requirement for detecting small time shifts. M. O'Donnell, A. R. Skovoroda, B. M. Shapo, S. Y. Emelianov, “Internal Displacement and Strain Imaging Using Ultrasonic Speckle Tracking”, IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol 41, no 3, 1994, doi: 10.1109 / 58.285465.

[00080] 4. Train a deep neural network to reconstruct the speed of sound image, without requiring real experimental data:

[00081] a) Generate training data using a large number of natural images to fully cover the wide range of potential geometry of SoS distributions;

[00082] b) Pretrain the model using theoretical time-shift maps and use more accurate theory to compute the theoretical time-shift map (see Example 1 below);

[00083] c) Fine-tune the model using acoustic simulations (such as k-Wave) to yield simulated time-shift maps;

[00084] d) As acoustic simulations are slow, apply the method described in Example 2 to improve efficiency;

[00085] e) To be able to train a large amount of data for high-resolution sound image reconstruction, use the methods described in Examples 3-5.

[00086] 5. Input the single robust time-shift map into the trained deep neural network to reconstruct the speed of sound image, without requiring a B-mode image as guidance or experimental training data.

[00087] EXAMPLE 1: Incorporating Rx path into the model. If the model (M. Jaeger, G. Held, S. Peeters, S. Preisser, M. Grünig, and M. Frenz, "Computed Ultrasound Tomography in Echo Mode for Imaging Speed of Sound Using Pulse-Echo Sonography: Proof of Principle", Ultrasound Med. Biol., vol 41, no 1, 2015, doi: 10.1016 / j.ultrasmedbio.2014.05.019, incorporated by reference herein in its entirety) does not calculate the Rx delay, there can be a mismatch between the theoretical time-shift map and the measured map. To fill the gap, the Rx path can be incorporated into the forward model and the beamforming process.

[00088] Results without considering Rx delay are shown in FIGS. 2A-C.

[00089] The ground-truth SoS can be found in FIG. 2D.

[00090] Compared to the theoretical time-shift map (FIG. 2A, Tx =15° relative to Tx=0°), the measured time-shift map (FIG. 2B) has an additional beam on the left (indicated with an arrow). As a result, the reconstructed SoS map has a hallucination on the left of the actual inclusion (FIG. 2C).

[00091] Mirrored receive aperture can be implemented to keep the SPF approximately invariant to transmit angles and therefore results in better correlation estimation. Its effect can be demonstrated using the strong-scatterers case, as shown in FIGS. 3A-B. FIG. 3A shows that without the mirrored receive aperture, tilted SPFs due to steering are observed. FIG. 3B shows the SPFs lie horizontally, regardless of the Tx steering.

[00092] When steering Tx waves with different angles, the SPFs rotate accordingly, which disrupts the phase shift estimation.

[00093] Considering Rx in modeling and beamforming: FIG. 4 provides an illustration of Tx and Rx.

[00094] The Rx path is specified in the Rx beamforming by using a small aperture at the Rx angle QRx. The Rx angle Prxwas kept the same as the Tx angle qrx (i.e., mirrored receive aperture (A. S. Podkowa and M. L. Oelze, “The convolutional interpretation of registration-based plane wave steered pulse-echo local sound speed estimators", Phys. Med. Biol., vol 65, no 2, 2020, doi: 10.1088 / 1361-6560 / ab6071, incorporated by reference herein in its entirety)). In the forward model, an integration term is added representing Rx delay along the Rx path:√1 + tan² QRX Z ∫ 0 dz'σ(x + ztanqRx – z'tanqRx, Z')Equation 1

[00095] Results obtained by taking into consideration the Rx path are shown in FIGS. 5A-C.

[00096] The theoretical time-shift map taking into consideration Rx delay is shown in FIG. 5A, using the same pair of Tx angles (15° to 0°). The measured time-shift map (FIG. 5B) now has a better alignment with the theoretical one. As a result, the reconstructed SoS map (FIG. 5C) successfully reveals the inclusion without previous hallucination (FIG. 2C).

[00097] By incorporating Rx path into the DL model, there is better alignment between the theoretical and measured time-shift maps and the reconstructed SoS map aligns better with the ground truth (FIG. 2D).

[00098] EXAMPLE 2: Using time-shift maps as input for deep learning-based speed of sound estimation. A deep learning (DL) pipeline was developed, which uses the time-shift map as an input for speed of sound (SoS) estimation. Results demonstrated the feasibility of such an approach.

[00099] 1. Method[000100]1.1 The deep learning pipeline (FIG. 6)[000101]The input to the neural network is a measured time-shift map, when changing Tx angle from 0° to 15°. The neural network follows an auto-encoder architecture, containing a down-sampling path for feature extraction and an up-sampling path for image reconstruction. The output is an SoS map.[000102]1.2 Generation of training data[000103]First, a number of SoS maps are generated with randomness. Each SoS map includes a homogeneous (SoS = 1540m / s) background and an inclusion. The location, size, and SoS value (selected from 1480-1600m / s range) of the inclusion are randomized.[000104]For each SoS map, the transmitted plane waves of a variety of angles were simulated. The echoes were collected, beamformed, and the time shift was calculated using different Tx angles. The time shifts between adjacent Tx angles are subtle, so they were concatenated together and the time shift map was obtained between 0° Tx angle and 15° Tx angle.[000105]Each pair, comprising a time-shift map and a ground-truth SoS map, serves as an input- output pair for training the neural network.[000106]1.3 Network training and testing[000107]A total of 62 pairs of time-shift and ground-truth SoS maps were collected, with 20 pairs for testing, while reserving the remaining 42 pairs for training.[000108]To evaluate the performance of the training and testing set, the mean absolute error (MAE) and root mean square error (RMSE) metrics were used to quantify the discrepancy between an estimated SoS map and the ground-truth SoS map.[000109]2. Results[000110]FIG. 7 presents samples of estimated SoS maps and ground-truth SoS maps. Overall, the neural network is able to reveal the inclusion in both the training set and the testing set. To quantify the performance, Table 1 presents MAE and RMSE of the estimation. Results demonstrate that the neural network maintains predictive power in the testing sets for SoS estimation.Table 1. Performance of SoS estimationTraining setTesting setMAE (m / s)1.351.49RMSE (m / s)3.003.50[000111]In this work, a deep learning pipeline enabling the SoS estimation using time-shift maps was developed. The performance of the approach shall also be benchmarked with 1) non-DL reconstruction approaches, and 2) a conventional DL approach using raw channel data as input.[000112]EXAMPLE 3: Comparing Deep Learning with the Analytical Method for Speed-of- Sound Estimation[000113]1. Merits of DL over the analytical method[000114]FIG. 8 puts the DL-based approach and the analytical method in context. Both methods take a time shift map as input, and both output an SoS map. While solving the same problem, the DL-based approach appears to have the merits of i) a reduced computational time, which enables real-time imaging and ii) improved estimation accuracy, potentially due to eased model assumption.[000115]1.1 Reduced computational time to enable real-time imaging[000116]The analytical method formulates SoS estimation as an inverse problem (Sanabria, 2018). The optimization process involved in solving it is computationally expensive. To give a perspective, it took 91 seconds to estimate an SoS map with 73 × 89 pixels. (The computational expense is in part due to the large matrix size and the resultant large inverse problem. As described by Sanabria et al., the SoS is first flattened into a vector of 73 × 89 = 6497 elements. Then a 6497 × 6497 matrix is needed to describe the linear transformation between an SoS map and a time-shift map.) Long computational time limits the frame rate and therefore prohibits the real-time application of SoS imaging.Table 2. Computational time of deep learning vs. analytical methodComputational time per SoS map (s)Equivalent frame rate (Hz)Analytical method910.011Deep Learning0.002500[000117]By contrast, the estimation done by deep learning is much more time-efficient, since only a forward pass is required. Table 2 compares the computational time required by deep learning vs. the analytical method. When estimating an SoS map with 73 × 89 pixels, deep learning demonstrates 45454 times efficiency improvement compared with the analytical method, and the resultant 500 Hz frame rates appear sufficient for real-time SoS imaging.[000118]1.2 Improved accuracy[000119]The analytical method holds a straight-ray assumption, which ignores wave propagation in the forms of refraction and diffraction. In this regard, deep learning does not pose such a straight-ray assumption and thus has the potential to approximate a more realistic propagation process.[000120]A neural network was trained with 70 pairs of time shift maps and ground-truth SoS maps. The network was then tested on 30 pairs of testing data. These testing data were also used to evaluate the performance of the analytical method. Table 3 compares the performance of deep learning and the analytical method, where deep learning demonstrates a slight gain in estimation accuracy.Table 3. Accuracy of deep learning vs. analytical methodRMSE (m / s)Analytical method5.99Deep Learning5.67RMSE: Root mean square error[000121]To probe the distribution of RMSE, the RMSE among different SoS groups is visualized in FIG. 9. Such a plot can also be found in Stähli et al. (P. Stähli, M. Kuriakose, M. Frenz, M. Jaeger, “Improved forward model for quantitative pulse-echo speed-of-sound imaging", Ultrasonics, vol 108, 2020, doi: 10.1016 / j.ultras.2020.106168, incorporated by reference herein in its entirety). Pixels are grouped by the ground-truth SoS value, and the RMSES are calculated for these groups. The background SoS is 1540 m / s. It is seen that for the SoS largely deviating from the 1540 m / s, the analytical method has larger errors. These larger errors may be attributed to the regularization used in the analytical method, which penalizes large deviations from the 1540 m / s and from adjacent pixels.[000122]EXAMPLE 4: SoS estimator trained with theoretical time-shift maps[000123]Simulated time-shift maps can be used to estimate SoS distribution, however, running each simulation (1 hr per SoS map) was time-consuming, making it practically prohibitive to collect sufficient samples for model training. The inventor has discovered that theoretical time-shift maps derived, e.g., using a forward model are much faster to obtain.[000124]The inventor demonstrates below the feasibility of using theoretical time-shift maps to train an estimator. Particularly, the trained estimator is interrogated with 5 patterns that are unseen in the training phase.[000125]1. Model training. The estimator adopts a U-net architecture. It takes theoretical time- shift maps as input and then outputs an SoS map. To enhance the diversity of training data, the theoretical time-shift maps are derived from a massive ImageNet dataset (O. Russakovsky et al., "ImageNet Large Scale Visual Recognition Challenge", Int. J. Comput. Vis., vol 115, no 3, 2015, doi: 10.1007 / s11263-015-0816-y, incorporated by reference herein in its entirety). With details provided in section 3, below, the estimator is trained with 10,000 time-shift maps and performs well in the training and testing sets. In the sections below, its performance is further validated in 5 tissue-mimicking SoS patterns that are unseen in the training phase.[000126]2. Reconstruct from the theoretical time shift maps[000127]Various regularization techniques were attempted: (i) mode filter applied on natural-image-converted SoS maps (aggregating fine details) and (ii) weight decay applied on U-net architecture. The parameters of regularizations are specified in each column of FIG. 10. Particularly, column 2 (window size of mode filter = 1, weight decay = 0) is equivalent to having no regularization.[000128]The inventor has demonstrated the use of theoretical time-shift maps for model training. After training with time-shift maps, the model is generalizable to unseen patterns. With desired generalizability, the model can serve as the backbone for further endeavors. For example, simulated data (or experimentally obtained data) can be used to fine-tune this model, either through incremental-learning techniques or through an auxiliary task.[000129]A key bottleneck is the dissimilarity between the theoretical time-shift map and the measured (simulated) time-shift map. A fundamental difference is that when measuring time shifts the phase shift is detected at distant spots, instead of universally spread in an area. This discrepancy cascades even when smoothing is performed.[000130]3. Data Processing and Training Performance[000131]Data processing - From ImageNet image to time-shift map[000132]FIG. 11 shows conversion from a natural image to a theoretical time-shift map.[000133]Performance for the training set is shown in FIG. 12 and for the testing set in FIG. 13.[000134]EXAMPLE 5: The overview of an exemplary SoS imaging approach, according to an embodiment of the invention, is illustrated in FIG. 14. The ultrasound transducer transmits a series of plane waves to insonify the medium from various angles. Upon receiving channel data, a customized beamforming algorithm is performed to obtain frames of beamformed RF data, where the phase shifts in the RF signals between different transmit angles are tracked and translated to time shifts. The time-shift maps contain full information on the SoS distribution and serve as the input to a neural network that reconstructs the SoS image. The use of the time-shift maps as the input not only enhances the model's generalizability and robustness but also enables a set of tailored data augmentation techniques critical to the successful training of the network.[000135]In principle, SoS distribution determines the travel time of sound waves and is therefore reflected in certain time-domain measurements. Consider a scatterer at location r and an acoustic pulse that takes a round-trip path P to reach the scatterer and return. The total time of flight can be expressed as an integration:t(r) = ∫p 1 c(rp) dlEquation 2[000136]where r'p denotes a location on the path P, and c(rp) denotes the SoS value at location r'p. dl is an infinitesimal displacement along a path. Equation 2 shows that t(r) is location dependent. However, the received channel data are a mixture of backscattered signals from various locations. To relate signals to particular locations, delay-and-sum (DAS) beamforming is often performed with an assumed SoS value cbf. When the actual SoS deviates from the assumed cbf, a time delay is introduced by beamforming:At(r) = ∫p ( 1 c(rp) 1 Cbf. ) dlEquation 3[000137]This time difference contains rich information on the SoS distribution along the propagation path. However, this difference cannot be measured on a single beamformed RF signal because there is no reference to compare with. Therefore, a pair of RF signals beamformed with different transmit (TX) angles (i.e., steering angles defined during transmission) or different receive (RX) angles (defined during beamforming) was considered. Because the time difference is small, measurement can be performed by tracking local phase shift (r) between the beamformed RF pair. The phase shift (r) is related to the time difference At(r) and the SoS c(r) by (Stähli, M. et al., “Improved forward model for quantitative pulse-echo speed-of-sound imaging," Ultrasonics, vol. 108, p. 106168, 2020):Φ(r) = 2πε [ At₁ (r) cos (PTX₁ - PRX₁) 2 At₂(r) cos (ΦTX2 - PRX2) 2 ] = 2πf [ ∫p₁ ( 1 cp₁ (r) 1 cbf ) dl cos (PTX, PRX) 2 ∫ Sp₂ ( 1 cp₂ (r) 1 cbf ) dl cos (PTX-9RX2) 2 ]Equation 4[000138]where fc is the center frequency of the pulse, qrx₁ and qrx₁ are the TX and RX angles of the i-th propagation path Pi, and i ∈ {1,2}, as shown in FIG. 15. The denominator accounts for the angle between the propagation path and the oscillation direction of RF, assuming two RF signals oscillate along the same direction.[000139]Equation 4 shows that the measured phase shift (r) is dependent on the transducer center frequency fc. Under the same SoS distribution, a higher center frequency leads to faster signal oscillation, causing a larger phase shift. To cancel the effect of center frequency, Φ(r) is divided by 2nf to obtain time shift τ(r):τ(r) = Φ(r) ·· 1 2πfcEquation 5[000140]The task of SoS imaging involves reconstructing SoS values from these time-shift values. This reconstruction is typically achieved by inverting the forward model (Equation 4) and assuming a straight-ray propagation path for P₁ and P2. However, the actual propagation of the acoustic pulse does not follow a straight ray, due to the presence of SoS variation. Even when proceeding with the idealized straight-way forward model, the inversion process is computationally expensive. To achieve higher realism and computational efficiency, we propose to utilize DL techniques for the reconstruction.[000141]These two categories of methods are described as follows.[000142]Ray-tracing modeling: Under the straight-ray assumption, time shift (r) can be modeled as a linear summation of slowness values (reciprocal of the SoS, i.e., 1 c(rp) ) along a straight ray path. With this linearity, the relationship between time shift and slowness can be expressed as a matrix multiplication:τ = L.SEquation 6[000143]where vector τ∈ RNqNxNz contains a total of No vectorized time shift maps, with each map having NxN₂ pixels. Matrix LE RNQNxNzXNxNz describes a forward model that relates measured τ to the underlying s. Vector s ∈ RNxNz denotes the vectorized slowness map and can be obtained by inverting the forward model with regularizations. Finally, the SoS map C∈ RNxXNz is obtained by reshaping s and taking element-wise reciprocal.[000144]DL: While the ray-tracing model offers a simplified solution, in reality, the relationship between time shift and SoS is more complicated because: 1) propagation of sound waves does not necessarily follow a straight path, and involves refraction, diffraction, and interference; and 2) processing of echoes involves beamforming and phase-shift tracking that have combined effects on the resulting time-shift measurement. Effects from these two sources contribute to a more realistic forward model Fprop+proc(C), describing the relationship between measured time shift maps TE IRNXNxXNz and the underlying SoS map C:T = Fprop+proc(C)Equation 7[000145]Accurate formulation of Fprop+proc(C) remains a challenging task. Inverting Fprop+proc(C) is also computationally expensive due to nonlinear components, e.g., the SoS- dependent refraction, hampering the potential for real-time application. Owing to these difficulties in forward modeling and inversion, in embodiments of the invention, the reconstruction function is learned directly using DL:C = Fprop+proc(T) = G(T)Equation 8[000146]where G(T) denotes a trained neural network.[000147]Upon receiving channel data, beamforming and phase shift tracking were performed to obtain T. In the context of SoS imaging, it is challenging to accurately track phase shift, due to an artifact introduced by changing TX angles (Stähli, 2020; Jaeger, M. and Frenz, M., “Towards clinical computed ultrasound tomography in echo-mode: Dynamic range artefact reduction," Ultrasonics, vol. 62, pp. 299–304, 2015; Podkowa, A.S. and Oelze, M. L., “The convolutional interpretation of registration-based plane wave steered pulse-echo local sound speed estimators," Phys. Med. Biol., vol. 65, no. 2, p. 025003, 2020). The artifact arises from the rotated oscillation direction of beamformed RF when TX angles are changed. To address this issue, in embodiments of the present invention, the method utilizes a common mid-angle tracking approach (Stähli, 2020; Jaeger, 2015; Haun, M. A. et al., “Overdetermined least-squares aberration estimates using common-midpoint signals," IEEE Trans. Med. Imaging, vol. 23, no. 10, pp. 1205–1220, 2004) that aligns the oscillation direction of the RFs. This approach results in high tracking accuracy (Stähli, 2020).[000148]Beamforming and phase-shift tracking[000149]The detailed process for how to obtain time-shift maps to be used as inputs to the DL model is illustrated in FIG. 16 and described as follows.[000150]Wave transmission. A sequence of 71 steered plane waves were transmitted. The sequence covered TX angles -17.5° to 17.5° with a 0.5° interval. Following each wave transmission, echoes were recorded in all transducer elements. The 71 transmissions led to 71 frames of channel data.[000151]Beamforming. Each frame of channel data was beamformed using the DAS algorithm with dynamic receive focusing with a constant f-number (=3). Importantly, beamforming was performed using a common mid-angle approach, that is, for a pre-specified mid-angle @mid defined as:Tx + 4 RX mid , 2Equation 9[000152]The RX angle QRX of the central element of the RX aperture used for beamforming was jointly determined by Фтx and Omid through Equation 9. Therefore, the width of the RX aperture was determined by the f-number and the center of the RX aperture was determined by the angle RX[000153]To enrich the information available in the final time-shift maps to be used for DL, beamforming was performed using multiple mid-angles: Omid = 0°, -7.5°, and 7.5°. Specifically, each of the 71 channel data frames was beamformed using Omid = 0°. The channel data frames obtained with QTX ≤ 2.5° were additionally beamformed using Omid = -7.5°, and the channel data frames obtained with Фтх ≥ −2.5° were additionally beamformed using mid = 7.5°. This process resulted in three groups of beamformed RF frames, each obtained from a separate mid-angle.[000154]After beamforming, coherent plane wave compounding was performed to reduce clutter. Compounding was performed by summing the neighboring frames (TX angle difference ≤2.5°) within the same mid-angle Omid. In the end, a beamformed RF image bRF(PTX, PRX) was uniquely obtained for a combination of qrx and PRX.[000155]Phase shift tracking. The phase shift was calculated based on the cross-correlation for a pair of beamformed RF frames.R(PTX₁) PRX1, PRX2, PRX2) = LPF(bRF(QTX₁ PRX₁) ⊗ bRF*(QTX2, PRX2))Equation 10[000156]The operator ⊗ represents the element-wise product. A low-pass filter LPF() was applied to filter out speckle noise (Podkowa, 2020). The bRF* denotes the complex conjugate of the beamformed RF frame. The phase shift (r) was calculated by taking the arctangent of R at position r.[000157]The phase shift was tracked between the beamformed images that shared a common mid-angle mia. In total, three common mid-angles (i.e., -7.5°, 0, and 7.5°) were used, leading to three time-shift maps, as shown in FIG. 16. Time shifts were converted from the measured phase shift using Equation 5. To prevent aliasing, phase tracking was performed in small angular steps that were later summed up (Stähli, 2020).[000158]Neural network[000159]Several theoretical studies suggest the suitability of convolutional operations in characterizing SoS-related time delays (Bezek, C. D. et al., “Learning the imaging model of speed- of-sound reconstruction via a convolutional formulation," arXiv Prepr. arXiv2309.00453, 2023; Beuret, S. et al., "Refraction-aware integral operator for speed-of-sound pulse-echo imaging," in IEEE International Ultrasonics Symposium, IUS, 2020, pp. 1-4). The convolutional neural network is a type of network architecture designed upon convolution operations and has proven successful in ultrasound imaging (Tehrani, A. K. Z. and Rivaz, H., "Displacement estimation in ultrasound elastography using pyramidal convolutional neural network," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 67, no. 12, pp. 2629-2639, 2020; Huang, X. et al., "Deep learning for ultrasound beamforming in flexible array transducer,” IEEE Trans. Med. Imaging, vol. 40, no. 11, pp. 3178–3189, 2021; Nair, A. A. et al., "A fully convolutional neural network for beamforming ultrasound images," in IEEE International Ultrasonics Symposium, IUS, 2018, pp. 1-4). In embodiments of the invention, a particular convolutional network, such as U-Net (Ronneberger, O. et al., “U-net: Convolutional networks for biomedical image segmentation," in Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 2015, vol. 9351, pp. 234-241) can be used for SoS reconstruction. In the encoding path, each encoding block contains two 3×3 convolutional layers with batch normalization (BN) and Tanh activations, followed by a 2×2 max pooling layer. A dropout layer with a 0.3 dropout rate was added after each encoding block to reduce overfitting. In the decoding path, each decoding block included a 2×2 transpose convolutional layer with a stride of 2 and two 3×3 convolutional layers with BN and Tanh. Finally, 3×3 convolutional layers were added as the output layer. Hyperparameters of the U-Net are determined on a validation set including 200 full-wave simulation samples.[000160]The loss function was composed of 1) the mean squared error (MSE) loss that quantifies pixel-wise difference, and 2) the structural similarity index measure (SSIM) loss that quantifies structural differences. The network was trained with the Adam optimizer, a learning rate of 1 x 10-3, and a batch size of 32.[000161]The model training was performed entirely using synthesized and simulation data.[000162]The use of time-shift maps as the input enables the efficient generation of a large number of synthesized time-shift and SoS map pairs to facilitate model training because computing the time-shift map from a given SoS map using the straight-ray theory is orders of magnitude faster than simulating the channel data from a given SoS map if the channel data were used as the input. However, the straight-ray computation is less accurate and less realistic than full-wave simulations. To combine the advantages of both computational methods (i.e., straight-ray computation and full-wave simulation), model training was divided into two stages:[000163](i) The pre-training stage: the model was trained with 1.2 million sample pairs computed by using the straight-ray forward model without requiring full-wave simulations; and[000164](ii) The fine-tuning stage: the model was fine-tuned with 1,000 sample pairs obtained from full-wave simulations where the time-shift maps were obtained from the simulated channel RF data for the given SoS maps.[000165]Full-wave simulation[000166]The k-Wave toolbox (Treeby, B. E. and Cox, B. T., “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields," J. Biomed. Opt., vol. 15, no. 2, pp. 021314-021314, 2010) was used for full-wave acoustic simulations to generate the training data required during the fine-tuning stage. A GE 9L-D transducer (GE Healthcare, Chicago, Illinois, USA) was simulated, with 192 transducer elements, a 0.23-mm pitch, and a 5.2-MHz center frequency. To reduce the cross-talking effect, a highly absorbing medium (30 dB / cm-MHz) was stuffed between transducer elements. The simulation medium was defined in two dimensions. A total of 1024×1024 grid points were simulated, including perfectly matching layers with a thickness of 20 grid points on all four sides to absorb waves. The spatial resolution of the simulation was 0.046 mm, and the temporal resolution was 5.97 ns with a Courant-Friedrichs- Lewy number of 0.2.[000167]In part of the simulation medium, the SoS maps were generated in a randomized way. The SoS maps contained ellipses and layered structures, mimicking lesions, and the layered tissue structure of the abdominal wall, as illustrated in FIGS. 17A-B. The SoS values varied from 1400 m / s to 1600 m / s. The background mass density was set to be 1000 kg / m³. Scatterers were simulated in the density domain, where the mass density of the scatterers varied between -3% and +3%.[000168]Domain-specific data augmentation[000169]Data synthesis with straight-ray computation. The physics model assuming straight rays is widely utilized in SoS imaging literature, resulting in a linear forward model applicable through matrix multiplication. This forward model was employed to efficiently generate time-shift maps. This process involves a linear transformation applied to SoS maps to produce corresponding time-shift maps, requiring less than 0.01 seconds for each SoS map, in contrast to the over 15 minutes needed for equivalent simulation data. Due to its time efficiency, a ray-tracing approach can effectively synthesize a substantial amount of training data.[000170]To diversify the image patterns, SoS maps were synthesized from natural images in the ImageNet dataset (2017) (Deng, J. et al., “ImageNet: A large-scale hierarchical image database,” in 2009 IEEE conference on computer vision and pattern recognition, 2010, pp. 248–255). The natural images were resized to 73 × 89 pixels to be consistent with the neural network input. Each natural image is converted to grayscale and proportionally mapped to the SoS values ranging from 1400 m / s to 1600 m / s. Finally, the images were converted to time-shift maps (FIG. 18). As a result, a total of 1.28 million input-output pairs were synthesized. These synthesized input-output pairs were used to extensively pretrain the neural network. After this pretraining stage, the more accurate but computationally-expensive full-wave simulated sample pairs were then used to fine-tune the model.[000171]Symmetry. Due to the symmetry of the imaging system, when an SoS map is horizontally flipped, the corresponding time-shift map shall also flip. In this way, a valid new pair of input-output maps can be generated. Not only does this technique double the amount of training data, but it also aids the neural network in learning the symmetry of the system.[000172]Induced noise. A difference between reality and simulation is that the signal tends to be noisier and less coherent in reality. As a result, the actual time shifts tend to be less salient, leading to reduced amplitude of time shifts. Such a reduction can happen at random locations, so a randomized two-dimensional Gaussian mask was used for noise generation.[000173]Proof-of-concept study on straight-ray synthesized time-shift maps[000174]In embodiments of the invention, a DL-based approach estimates SoS distribution by performing an image-to-image translation, rather than inversing a linear forward model. This choice represents an exploration beyond the linear forward modeling, relaxing the underlying straight-ray assumption. With advantages to be further demonstrated, the first goal was to establish that the proposed time-shift DL at least possesses the quality of a linear problem solver.[000175]Given a ground truth SoS map, the corresponding ray-tracing synthesized time-shift maps are calculated using the state-of-the-art linear forward model (Stähli, 2020). By inverting the synthesized time-shift maps to an SoS map, the DL is essentially solving the linear forward problem.[000176]To intensively test the learning capabilities of DL, time-shift maps were generated from the SoS maps that had complicated patterns. These SoS maps were derived from the ImageNet dataset. A total of 1.2 million samples, including the SoS maps and the corresponding synthesized time-shift maps, were used to train the DL model. The trained DL model was then tested on 80,000 hold-out samples.[000177]Robustness test on simulation data[000178]The acoustic signals captured by ultrasound transducers are influenced by a variety of acoustic properties and pulse transmit settings. A robust DL model is supposed to accurately discern the SoS distribution from other coinciding acoustic properties and imaging system settings. To assess the model's robustness, a controlled study was derived wherein the SoS maps remained constant while other parameters were modified. A robust model is expected to produce relatively consistent outcomes. This study was conducted in simulation due to the convenience of adjusting medium properties and pulse settings in simulations.[000179]Two types of fluctuations were simulated to test model robustness: density distribution and pules settings.[000180]Density distribution is an important medium property that impacts sound wave propagation. When generating training data, the density was set to be uncorrelated with the SoS, preventing information leakage from the SoS map. In reality, however, soft tissues may have density values positively correlated with the SoS (Mast, T. D. "Empirical relationships between acoustic parameters in human soft tissues," Cit. Acoust. Res. Lett. Online, vol. 1, no. 2, pp. 37– 42, 2000). This scenario was also simulated for robustness testing.[000181]Pulse settings may differ across transducer probes. Even in the same transducer, parameters such as the center frequency are often tuned to optimize image quality. A typical example is to decrease the center frequency for better penetration through the abdominal wall. To test such a scenario, we decreased the center frequency from the 5.2 MHz used in the training data to 4 MHz for testing and performed the simulation.[000182]Full-wave simulations were performed on the digital phantoms mimicking tissue structures as illustrated in FIGS. 19A-B.[000183]Phantom experiments[000184]The experimental data were acquired using the GE 9L-D linear transducer probe connected to a Vantage 256 system (Verasonics, Kirkland, WA, USA). The transducer probe features a 5.2-MHz center frequency, 192 elements, and a 0.23-mm pitch. Two types of tissue- mimicking phantoms were imaged: 1) Multi-modality breast biopsy and sonographic trainer phantom (Model 073) and 2) Triple modality 3D abdominal phantom (Model 057A) (CIRS Inc., Norfolk, VA, USA). The reference SoS value of each phantom compartment was provided by the phantom manufacturer.[000185]The performances of three reconstruction methods were compared as follows: 1) a state-of-the-art physics-based inversion that applies linear programming on the time-shift maps to yield the SoS map, 2) DL using channel data as the input, and 3) the proposed physics-guided DL that uses the time-shift map as the input, according to embodiments of the invention. The characteristics of these three methods are summarized in Table 4.[000186]Table 4. Summary of the Reconstruction MethodsMethodInputStraight-ray assumptionModel linearityPhysics-based inversionTime-shift mapsYesLinearChannel-data DLChannel dataNoNon-linearTime-shift DLTime-shift mapsNoNon-linear[000187]The physics-based inversion involved beamforming and phase tracking. The SoS map was obtained by solving a linear optimization problem (Sanabria, S. J. et al., “Spatial domain reconstruction for imaging speed-of-sound with pulse-echo ultrasound: Simulation and in vivo study," Phys. Med. Biol., vol. 63, no. 21, p. 215015, 2018; Stähli, 2020). Following Sanabria, 2018, the problem was solved using MATLAB CVX toolbox, a disciplined convex programming framework (Grant, M. and Boyd, S., "CVX: Matlab software for disciplined convex programming, ver 2.1,", 2017).[000188]The channel-data DL method was implemented using the U-net as the backbone, and the network architecture design followed (Feigin, M. et al., “A deep learning framework for single- sided sound speed inversion in medical ultrasound," IEEE Trans. Biomed. Eng., vol. 67, no. 4, pp. 1142-1151, 2020; F. Khun Jush et al., "Deep learning for ultrasound speed-of-sound reconstruction: Impacts of training data diversity on stability and robustness," Mach. Learn. Biomed. Imaging, vol. 2, pp. 202-236, 2023).[000189]For a fair comparison, all three methods were set to have the same output size of 73 x 50 pixels. This reconstruction view was located in the center of the input time-shift maps. This region was located at the intersection of the plane waves with different Tx angles.[000190]FIG. 20 provides the exemplary SoS maps obtained from the testing set. Visually, the reconstructed SoS maps share high similarities with the ground truths. A quantitative evaluation throughout the 80,000 test samples yielded a root mean square error (RMSE) of 10.5 m / s. Relative to the ground truth SoS values, the reconstruction had the root mean squared percent error (RMSPE) of 0.68%. The low errors indicate the capability of DL in solving the classic inverse problem for SoS imaging.[000191]By closely examining the reconstructed SoS map, it is observed that the reconstruction error concentrates on the bottom region. This error at the bottom region is expected, given that the time-shift signal occurs in the posterior region, whereas the bottom pixels lack the posterior region beneath them. A study using physics-based inversion also reported an artifact at the bottom (Jaeger, M. et al., "Computed ultrasound tomography in echo mode for imaging speed of sound using pulse-echo sonography: Proof of principle," Ultrasound Med. Biol., vol. 41, no. 1, pp. 235- 250, 2015). To avoid the artifact, it is suggested to set the time-shift maps deeper than the SoS maps, allowing adequate incorporation of time-shift information at the bottom.[000192]Simulation study[000193]FIG. 21 displays the ground-truth SoS maps and the reconstructed SoS maps, where each row corresponds to a different digital phantom. Initially, an idealized condition was considered, whose density assumption and pulse settings were identical to those used for generating the training data. However, in reality, properties of the medium (such as density distribution) vary across tissues, making these properties challenging to anticipate. Additionally, the center frequency of the pulse can deviate from its nominal value. As a result, a pre-specified simulation condition might not accurately reflect these real-world variations. To assess the model's robustness under these conditions, the reconstruction performance was evaluated at two alternated simulation conditions. Results were reported for a total of three simulation conditions:[000194]Original condition was the condition for model training, with a 5.2-MHz center frequency and the mass density being uncorrelated with SoS.[000195]Altered condition 1 had an altered density assumption, with mass density linearly correlated with SoS.[000196]Altered condition 2 had an altered center frequency, with the center frequency changed from 5.2 MHz to 4 MHz.[000197]As shown in FIG. 22, in the original condition, all three methods demonstrated the ability to correctly register the layered boundaries and inclusion location. Among these reconstruction methods, the physics-based inversion showed lower SoS contrast than DL-based methods and the ground truth. This low contrast could be attributed to the spatial-gradient regularization used by the physics-based inversion. The spatial-gradient regularization, sometimes known as total-variation regularization, was shown necessary for solving the ill-posed inverse problem in SoS imaging (Sanabria, 2018; Stähli, 2020). On the flip side, the regularization penalized the spatial variation of SoS values, resulting in reduced SoS contrast in both layered medium and circular inclusions.[000198]Results in altered conditions highlight the susceptibility of channel-data DL. When the density assumption is altered, the channel-data DL shows reduced SoS contrast and a faint inclusion boundary. In the case of altered center frequency (i.e., altered condition 2), the SoS contrast was further decreased on layered media but remained strong on the inclusion. The perseverance of inclusion was hypothetically related to the distinctive hyperbolas in channel data. Under the change of center frequency, the hyperbolas could still be perceived by channel-data DL to identify the inclusion. The reconstruction errors, measured by RMSE, are presented in FIG. 22. All three methods achieved their best performance in the original simulation condition and then declined in the altered conditions. Among the three methods, the channel-data DL showed the largest increase in RMSE, from 26.2 m / s to 37.5 m / s. The two other methods, i.e., physics-based inversion and time-shift DL, experienced significantly less performance decline when applied to altered simulation conditions.[000199]Cross-section profiles of the reconstructed SoS maps are provided in FIG. 23. These profiles were obtained along the vertical line at the center of the four-layer phantom. The ground-truth SoS values exhibited sharp changes at the layered boundary. These edges were effectively captured by the physics-based inversion. However, the SoS variation in the physics-based inversion was consistently smaller than the ground truth. In terms of the channel-data DL, the cross section of SoS map has a good alignment with the ground truth from z = 14 mm to z = 19 mm in the original condition but appeared flatter at other location locations and in altered simulation conditions. This bias was reduced in the time-shift DL, with reconstructed SoS values fluctuating around the ground truth.[000200]Phantom experiments[000201]FIG. 24 presents the experimental results conducted on a tissue-mimicking abdominal phantom (top row) and breast phantom (middle and bottom rows). The abdominal phantom, from top to bottom, had SoS values of 1430 m / s for the fat-mimicking layer, 1555 m / s for the muscle- mimicking layer, and 1540 m / s for the liver-mimicking region as provided by the vendor. The breast phantom had SoS values of 1460 m / s for the top layer, 1520 m / s for the background, and 1580 m / s for the hypoechoic lesions as provided by the vendor. Presented in FIG. 24 are (from the left to right): B-mode images, reference SoS maps created using vendor-provided SoS values, and SoS maps reconstructed using the physics-based inversion, channel-data DL, an ablated version of time-shift DL (without the proposed domain-specific data augmentation described previously), and the full proposed time-shift DL.[000202]In the B-mode images, bright reflections resembling the noise present in real tissue can be observed on the breast phantom.[000203]It is observed that the physics-based inversion effectively differentiated the layered structure at reconstructed SoS maps, especially for the abdominal phantom. For the breast phantom, although the inclusions are not as salient, there are indeed contrasts reconstructed on the inclusion regions.[000204]In contrast, the channel-data DL approach had challenges in distinguishing tissue- mimicking compartments.[000205]Overall, the proposed time-shift DL exhibited the best reconstruction accuracy and the highest contrast-to-noise ratio (CNR). Compared with the physics-based inversion, the inventive approach reduced RMSE by 26% and the mean absolute error (MAE) by 32%. The CNR is improved by 52%, allowing a better differentiation of the breast lesions.[000206]As an ablation study, the time-shift DL without domain-specific data augmentations performed worse than the full model proposed time-shift DL method, showing the effectiveness of the proposed domain-specific data augmentations.[000207]Quantitative results comparing these methods are presented in Table 5, which shows consistent results with visual observations made in FIG. 24.[000208]Table 5. Reconstruction Performances in Phantom ExperimentsMethodRMSE (m / s)MAE (m / s)CNRPhysics-based inversion21.616.71.71Channel-data DL80.172.90.48Time-shift DL (Ablation, without data augmentations)16.713.12.09Time-shift DL (proposed)16.011.32.60[000209]In embodiments, the present invention comprises a robust DL approach that is capable of reconstructing SoS in tissue-mimicking phantoms, without requiring experimental training data.[000210]Pulse-echo SoS imaging has been a longstanding challenging task. Despite the availability of equipment (i.e., multi-element linear electronic arrays) for over half a century, the proof-of-concept study of SoS imaging has only emerged within a decade (Jaeger, 2015, Ultrasound Med. Biol.). In recent years, advancements in DL and simulation software provided an opportunity for pursuing SoS imaging in a data-driven manner. Several DL approaches were proposed to reconstruct SoS maps directly from raw channel data. While these approaches showed admirable results in simulation data, they encountered a drastic performance decline in the experiment data. The experimental results described herein demonstrated that the trained channel- data DL is susceptible to the alternation of medium density assumption and the fluctuation in the transmitted waveform. In contrast, the time-shift DL shows better robustness, thanks to the domain-specific phase-tracking technique that effectively extracts SoS-related information from the raw data. FIGS. 25A-D provide exemplary time-shift maps obtained on a breast phantom (FIGS. 25B-D), along with B-mode image for reference (FIG. 25A). A pattern resembling a rocket's tail flame is observed posterior to the inclusion, revealing a high spatial affinity.[000211]The present invention has been described with reference to particular embodiments having various features. In light of the disclosure provided above, it will be apparent to those skilled in the art that various modifications and variations can be made in the practice of the present invention without departing from the scope or spirit of the invention. One skilled in the art will recognize that the disclosed features may be used singularly, in any combination, or omitted based on the requirements and specifications of a given application or design. When an embodiment refers to "comprising" certain features, it is to be understood that the embodiments can alternatively "consist of" or "consist essentially of" any one or more of the features. Any of the methods disclosed herein can be used with any of the systems disclosed herein or with any other system, and vice versa. Other embodiments of the invention will be apparent to those skilled in the art from consideration of the specification and practice of the invention.[000212]It is noted in particular that where a range of values is provided in this specification, each value between the upper and lower limits of that range is also specifically disclosed. The upper and lower limits of these smaller ranges may independently be included or excluded in the range as well. The singular forms "a," "an," and "the" include plural referents unless the context clearly dictates otherwise. It is intended that the specification and examples be considered as exemplary in nature and that variations that do not depart from the essence of the invention fall within the scope of the invention. Further, all of the references cited in this disclosure are each individually incorporated by reference herein in their entireties and as such are intended to provide an efficient way of supplementing the enabling disclosure of this invention as well as provide background detailing the level of ordinary skill in the art.

Claims

1.A method of creating a trained deep neural network, the method comprising:providing one or more natural image(s);converting each of the one or more natural image(s) into a speed of sound (SOS) map;for each of the SoS maps, deriving a theoretical time-shift map;training a deep neural network using time-shift and SoS map pairs, wherein each of the pairs comprises one of the SoS maps and the theoretical time-shift map derived therefrom; andoptionally, fine-tuning the deep neural network using simulated time-shift and SoS map pairs, each of the pairs comprising a SoS map and its corresponding time-shift maps obtained through acoustic simulations, e.g., full-wave simulations;wherein the deep neural network is capable of reconstructing a speed of sound (SoS) image from a measured time-shift map generated from ultrasound.2.An imaging method comprising:imaging an object using ultrasound;converting, e.g., in real-time, ultrasound signals obtained from the imaging into one or more measured time-shift map(s);inputting the measured time-shift map(s) into the deep neural network of claim 1; andobtaining as an output a speed of sound (SoS) image reconstructed from the measured time-shift map(s).3.An imaging method, e.g., a real-time imaging method, comprising:imaging an object using ultrasound;converting ultrasound signals obtained from the imaging into one or more measured time-shift map(s);providing a speed of sound (SoS) image by analyzing the measured time-shift map(s) with a trained deep neural network.4.A device comprising:memory storing instructions of a program; anda processor that generates a speed of sound (SoS) image;wherein the processor is capable of executing the instructions to reconstruct the SoS image from measured time-shift map(s) derived from ultrasound signals; andwherein the reconstruction of the SoS image is provided in response to input of the measured time-shift map(s) into a trained deep neural network.5.The method or device of any of claims 1-4, wherein the deep neural network is not trained on:actual / experimental ultrasound data / signals; and / orgeometrical information of the object to be imaged.6.The method or device of any of claims 1-4, wherein one or more or each of the theoretical and / or simulated time-shift maps account for receive (Rx) delay.7.The method or device of any of claims 1-4, wherein the measured time-shift map accounts for Rx delay.8.The method or device of any of claims 1-4, wherein one or more or each of the theoretical and / or simulated time-shift maps is derived using dynamic programming and / or vectorization.9.The method or device of any of claims 1-4, wherein the measured time-shift map is derived using dynamic programming and / or vectorization.10.The method or device of any of claims 1-4, wherein the deep neural network is trained by:converting a representative time-shift map into an SoS map;inputting the SoS map into a theoretical time-shift map model;obtaining a predicted time-shift map for the SoS map; andcomparing and minimizing any difference between the predicted time-shift map and the representative time-shift map.11.The method or device of any of claims 1-4, wherein the deep neural network is trained by data parallelism involving multiple GPUs.12.The method or device of any of claims 1-4, wherein the measured time-shift map is obtained in real time from ultrasound imaging of a patient.13.The method or device of any of claims 1-4, wherein the training is performed in a manner, such as to a degree, to obtain an ultrasound-based system-independent, e.g., generalizable, trained deep neural network.14.The method or device of any of claims 1-4, wherein the reconstructed SoS image is based on the inputting into the trained deep neural network of:only a single measured time-shift map or multiple measured time-shift maps; and no geometric information of the object to be imaged.15.The method or device of any of claims 1-4, wherein no B-mode image is used as guidance and / or the deep neural network is not trained on experimental / actual ultrasound data / signals.16.The method or device of any of claims 1-4, wherein:no measured ultrasound data / signals are used in the training of the deep neural network;i.e., only the theoretical and / or simulated ultrasound data are used.17.The method or device of any of claims 1-4, wherein the deep neural network has a U-net architecture.18.The method or device of any of claims 1-4, wherein the deep neural network has a GAN architecture.19.The method or device of any of claims 1-4, wherein data provided for the deep neural network is provided in a radiofrequency (RF) format.20.The method or device of any of claims 1-4, wherein data provided for the deep neural network is provided in an in-phase and quadrature (IQ) format.21.The method or device of any of claims 1-4, wherein the object or patient comprises a body, or part thereof, such as a human or animal body, including tissue and / or one or more organs thereof.22.The method or device of any of claims 1-4, wherein the theoretical and / or the measured time-shift maps between adjacent steering angles are computed by using RF or IQ data from adjacent steering angles.23.The method or device of any of claims 1-4, wherein the time-shift maps are concatenated to make a single robust time-shift map.24.The method or device of any of claims 1-4, comprising:determining if an SoS image needs to be learned by estimating the theoretical time-shift map, instead of performing an acoustic simulation, using a physically-derived forward model; andidentifying and selecting therefrom additional training data.25.The method or device of any of claims 1-4, comprising selectively running simulations using the additional training data selected.26.The method or device of any of claims 1-4, wherein the measured time-shift maps are derived from data obtained from transmitting at multiple steering angles.27.The method or device of any of claims 1-4, wherein the measured time-shift maps are derived from data obtained from different sized apertures.28.The method or device of any of claims 1-4, wherein the Rx path is incorporated into the forward model and / or the beamforming process.29.The method or device of any of claims 1-4, comprising using mirrored receive aperture, wherein the transmit (Tx) path is at Tx angle QTx, the Rx path is at an angle Rx, and the Rx angle QRx is kept the same as the Tx angle QTx.30.The method or device of any of claims 1-4, wherein in the forward model, an integration term is added representing Rx delay along the Rx path:√1 + tan² RX Z ∫ 0 dz'σ(x + ztanqRx z'tanorx, z')31.The method or device of any of claims 1-4, wherein the deep neural network comprises a physics-based deep learning training framework.32.The method or device of any of claims 1-4, further comprising determining a total time of flight as t(r) = ∫p 1 c(rp)' dl.33.The method or device of any of claims 1-4, further comprising determining a time delay introduced by beamforming as At (r) = ∫p ( 1 c(rp) 1 Cbf ) dl.34.The method or device of any of claims 1-4, further comprising determining a phase shift as Φ(r) = 2πfc [ ∫P1 ( 1 cP1(r) 1 cbf ) dl cos ( TX1-RX1 2 ) ∫P2 ( 1 cP2(1) 1 cbf ) dl cos ( TX2-RX2 2 ) ]35.The method or device of any of claims 1-4, further comprising determining a time shift as τ(r) = Φ(r) 1 2πfc .36.The method or device of any of claims 1-4, further comprising expressing time shift as t = L.S.37.The method or device of any of claims 1-4, further comprising performing beamforming using a common mid-angle approach with a pre-specified mid-angle defined as @mid = QTX+PRX 2 .38.The method or device of any of claims 1-4, further comprising determining a phase shift based on cross-correlation for a pair of beamformed RF frames according to:R(PTX1, PRX1» PRX2» PRX2) = LPF(bRF(QTX₁, PRX₁) ⊗ bRF* (QTX2, PRX2)).

Images