Forward model for quantitative pulse-echo speed-of-sound imaging

TL;DR

This paper introduces an improved forward model for quantitative pulse-echo speed-of-sound imaging that enhances the robustness and accuracy of tissue SoS reconstruction in ultrasound tomography, demonstrated through phantom and volunteer studies.

Contribution

The paper presents a novel forward model that accounts for specific phase shifts and position errors, significantly improving SoS imaging accuracy over previous models.

Findings

Enhanced SoS reconstruction accuracy in phantom studies

Improved robustness of quantitative results with the new model

Preliminary volunteer results support clinical potential

Abstract

Computed ultrasound tomography in echo mode (CUTE) allows determining the spatial distribution of speed-of-sound (SoS) inside tissue using handheld pulse-echo ultrasound (US). This technique is based on measuring the changing phase of beamformed echoes obtained under varying transmit (Tx) and/or receive (Rx) steering angles. The SoS is reconstructed by inverting a forward model describing how the spatial distribution of SoS is related to the spatial distribution of the echo phase shift. CUTE holds promise as a novel diagnostic modality that complements conventional US in a single, real-time handheld system. Here we demonstrate that, in order to obtain robust quantitative results, the forward model must contain two features that were not taken into account so far: a) the phase shift must be detected between pairs of Tx and Rx angles that are centred around a set of common mid-angles, and b) it must account for an additional phase shift induced by the error of the reconstructed position of echoes. In a phantom study mimicking liver imaging, this new model leads to a substantially improved quantitative SoS reconstruction compared to the model that has been used so far. The importance of the new model as a prerequisite for an accurate diagnosis is corroborated in preliminary volunteer results.

Tables

Table 1. Table 1 : Combinations of Tx ( ϕ n ) \phi_{n}) , Rx ( ψ m , n subscript 𝜓 𝑚 𝑛 \psi_{m,n} ) and common mid-angles ( γ 𝛾 \gamma ) used in this study (in ∘ ). Tracking is performed between angle combinations having the same mid-angle γ 𝛾 \gamma , indicated by fields having the same gray value. No tracking was performed for γ 𝛾 \gamma = ± plus-or-minus \pm 25 ∘ .

${\psi }_{m,n}$	${\varphi }_n$	-25	-15	-5	5	15	25
-25		-25	-20	-15	-10	-5	0
-15		-20	-15	-10	-5	0	5
-5		-15	-10	-5	0	5	10
5		-10	-5	0	5	10	15
15		-5	0	5	10	15	20
25	$\gamma$	0	5	10	15	20	25

Equations

$$s(t,n,k)=G(t-t_0({\text{r}}^{\prime },n,k),{\text{r}}^{\prime },k)\cdot exp\left[2\pi i{f}_0\left(t-t_0({\text{r}}^{\prime },n,k)\right)\right]$$

$$\begin{array}{rl}u(\text{r}={\text{r}}^{\prime },n)& =\sum _k\left[s\left({\hat{t}}_0(\text{r},n,k),n,k\right)\right]=\sum _k\left[s\left({\hat{t}}_0({\text{r}}^{\prime },n,k),n,k\right)\right]\\ & =\sum _k\{G\left(({\hat{t}}_0-t_0)({\text{r}}^{\prime },n,k),{\text{r}}^{\prime },k\right)\cdot exp\left[2\pi i{f}_0\left({\hat{t}}_0-t_0\right)({\text{r}}^{\prime },n,k)\right]\}\end{array}$$

$$u(n)=\sum _k\{G(-{\tau }_{tx}(n)-{\tau }_{rx}(k),k)\cdot exp\left[2\pi i{f}_0\left(-{\tau }_{tx}(n)-{\tau }_{rx}(k)\right)\right]\}$$

$$u(n)\simeq exp\left(-2\pi i{f}_0{\tau }_{tx}(n)\right)\sum _k\left[G\left(-\overline{{\tau }_{tx}}-{\tau }_{rx}(k),k\right)\cdot exp\left(-2\pi i{f}_0{\tau }_{rx}(k)\right)\right]$$

$$\Delta \Theta ({\text{r}}^{\prime },n,n^{\prime })\simeq arg\left[{\int }_{-0.5\Delta x}^{0.5\Delta x}{\int }_{-0.5\Delta z}^{0.5\Delta z}d\hat{\mathbf{r}}\left\{u({\text{r}}^{\prime }+\hat{\mathbf{r}},n)\cdot u^{\star }({\text{r}}^{\prime }+\hat{\mathbf{r}},n^{\prime })\right\}\right]$$

$$\Delta \Theta ({\text{r}}^{\prime },n,n^{\prime })\simeq 2\pi f_0\Delta {\tau }_{tx}({\text{r}}^{\prime },n,n^{\prime })=2\pi f_0\left\{{\tau }_{tx}({\text{r}}^{\prime },n^{\prime })-{\tau }_{tx}({\text{r}}^{\prime },n)\right\}$$

$${\tau }_{tx}({\text{r}}^{\prime },n)={\int }_{j_n}^{{\text{r}}^{\prime }}dl\left\{\frac{1}{c(\text{r})}-\frac{1}{\hat{c}(\text{r})}\right\}\equiv {\int }_{j_n}^{{\text{r}}^{\prime }}dl\Delta \sigma (\text{r})$$

$$u(n)=exp\left(-2\pi i{f}_0{\tau }_{tx}(n)\right)\sum _k\left\{G\left(-\overline{{\tau }_{tx}}-\overline{{\tau }_{rx}},n,k\right)\cdot exp\left(-2\pi i{f}_0{\tau }_{rx}(k)\right)\right\}$$

$$\Delta \Theta ({\text{r}}^{\prime },n,n^{\prime },m)=arg\left[{\int }_{-0.5\Delta x}^{0.5\Delta x}{\int }_{-0.5\Delta z}^{0.5\Delta z}d\hat{\mathbf{r}}\left(u({\text{r}}^{\prime }+\hat{\mathbf{r}},n,m)\cdot u^{\star }({\text{r}}^{\prime }+\hat{\mathbf{r}},n^{\prime },m)\right)\right]$$

$$\begin{array}{rl}\Delta \Theta ({\text{r}}^{\prime },n,n^{\prime },m)& =2\pi f_0\left(\Delta {\tau }_{tx}({\text{r}}^{\prime },n,n^{\prime })+\Delta {\tau }_{rx}({\text{r}}^{\prime },m,n,n^{\prime })\right)\\ & =2\pi f_0\left({\tau }_{tx}({\text{r}}^{\prime },n^{\prime })-{\tau }_{tx}({\text{r}}^{\prime },n)+{\tau }_{tr}({\text{r}}^{\prime },m,n)-{\tau }_{tr}({\text{r}}^{\prime },m,n^{\prime })\right)\end{array}$$

$$d(n,m)=\frac{\hat{c}\left[{\tau }_{tx}(n)+{\tau }_{rx}(m,n)\right]}{2\cos \left[\frac{1}{2}\left({\varphi }_n-{\psi }_{m,n}\right)\right]}$$

$$\Lambda (n,m)=\frac{\frac{\hat{c}}{f_0}}{2\cos \left[\frac{1}{2}\left({\varphi }_n-{\psi }_{m,n}\right)\right]}$$

$$\Delta d=\frac{\hat{c}\left[{\tau }_{tx}({\varphi }_n+\Delta \varphi )+{\tau }_{rx}({\psi }_{m,n}-\Delta \varphi )\right]}{2\cos \left[\frac{1}{2}\left({\varphi }_n-{\psi }_{m,n}+2\Delta \varphi \right)\right]}-\frac{\hat{c}\left[{\tau }_{tx}({\varphi }_n)+{\tau }_{rx}({\psi }_{m,n})\right]}{2\cos \left[\frac{1}{2}\left({\varphi }_n-{\psi }_{m,n}\right)\right]}$$

$$\Lambda =\frac{\frac{\hat{c}}{f_0}}{2}$$

$$\begin{array}{rl}\Delta \Theta ({\text{r}}^{\prime },n,n^{\prime },m)& =\frac{2\pi \Delta d}{\Lambda }\\ & =2\pi f_0\left\{\frac{{\tau }_{tx}(n)+{\tau }_{rx}(m,n^{\prime })}{\cos \left[\frac{1}{2}\left({\varphi }_{n^{\prime }}-{\psi }_{m,n^{\prime }}\right)\right]}-\frac{{\tau }_{tx}(n)+{\tau }_{rx}(m,n)}{\cos \left[{\varphi }_n-{\psi }_m\right]}\right\}\\ & =2\pi f_0\left\{\frac{{\tau }_{tx}({\varphi }_n+\Delta \varphi )+{\tau }_{rx}({\psi }_{m,n}-\Delta \varphi )}{\cos \left[\frac{1}{2}\left({\varphi }_n-{\psi }_{m,n}+2\Delta \varphi \right)\right]}-\frac{{\tau }_{tx}({\varphi }_n)+{\tau }_{tx}({\psi }_{m,n})}{\cos \left[\frac{1}{2}\left({\varphi }_n-{\psi }_{m,n}\right)\right]}\right\}\end{array}$$

$$\Delta \sigma ({\text{r}}^{\prime })$$

$$\Delta \Theta ({\text{r}}^{\prime },n,n^{\prime })$$

$$\Delta \Theta ({\text{r}}^{\prime },n,n^{\prime },m)$$

$$\Delta {\Theta }_{j,k}(n,n^{\prime })=$$

$$2\pi f_0\left[\sum _{j^{\prime },k^{\prime }}{w}_{(j,k,j^{\prime },k^{\prime })}^{{\varphi }_{n^{\prime }}}\Delta {\sigma }_{(j^{\prime },k^{\prime })}-\sum _{j^{\prime },k^{\prime }}{w}_{(j,k,j^{\prime },k^{\prime })}^{{\varphi }_n}\Delta {\sigma }_{(j^{\prime },k^{\prime })}\right]$$

$$\Delta {\Theta }_{j,k}(n,n^{\prime },m)=$$

$$2\pi f_0[\left(\sum _{j^{\prime },k^{\prime }}{w}_{(j,k,j^{\prime },k^{\prime })}^{{\varphi }_{n^{\prime }}}\Delta {\sigma }_{(j^{\prime },k^{\prime })}+\sum _{j^{\prime },k^{\prime }}{w}_{(j,k,j^{\prime },k^{\prime })}^{{\psi }_{n^{\prime },m}}\Delta {\sigma }_{(j^{\prime },k^{\prime })}\right)-$$

$$\left(\sum _{j^{\prime },k^{\prime }}{w}_{(j,k,j^{\prime },k^{\prime })}^{{\varphi }_n}{\sigma }_{(j^{\prime },k^{\prime })}+\sum _{j^{\prime },k^{\prime }}{w}_{(j,k,j^{\prime },k^{\prime })}^{{\psi }_{n,m}}{\sigma }_{(j^{\prime },k^{\prime })}\right)]$$

$$\Delta {\Theta }_{j,k}(n,n^{\prime },m)=$$

$$\frac{2\pi f_0}{\cos \left[\frac{1}{2}\left({\varphi }_{n^{\prime }}-{\psi }_{n^{\prime },m}\right)\right]}\left(\sum _{j^{\prime },k^{\prime }}{w}_{(j,k,j^{\prime },k^{\prime })}^{{\varphi }_{n^{\prime }}}\Delta {\sigma }_{(j^{\prime },k^{\prime })}+\sum _{j^{\prime },k^{\prime }}{w}_{(j,k,j^{\prime },k^{\prime })}^{{\psi }_{n^{\prime },m}}\Delta {\sigma }_{(j^{\prime },k^{\prime })}\right)-$$

$$\frac{2\pi f_0}{\cos \left[\frac{1}{2}\left({\varphi }_n-{\psi }_{n,m}\right)\right]}\left(\sum _{j^{\prime },k^{\prime }}{w}_{(j,k,j^{\prime },k^{\prime })}^{{\varphi }_n}\Delta {\sigma }_{(j^{\prime },k^{\prime })}+\sum _{j^{\prime },k^{\prime }}{w}_{(j,k,j^{\prime },k^{\prime })}^{{\psi }_{n,m}}\Delta {\sigma }_{(j^{\prime },k^{\prime })}\right)$$

$$\Delta \mathbf{\theta }=\mathbf{M}\Delta \mathbf{\sigma }$$

$$\Delta \mathbf{\sigma }$$

$${\mathbf{M}}^{\text{inv}}$$

$$C(\Delta \sigma )=\Vert \Delta \mathbf{\theta }-\mathbf{M}\Delta \mathbf{\sigma }{\Vert }_2^2+\Vert {\gamma }_x{\mathbf{D}}_x\Delta \mathbf{\sigma }+{\gamma }_z{\mathbf{D}}_z\Delta \mathbf{\sigma }{\Vert }_2^2$$

$$u(x,z,\varphi )$$

$${\hat{t}}_0(x,z,\varphi ,k)$$

Full text

\marginsize

3.5cm2.5cm1cm1cm

Forward model for quantitative pulse-echo speed-of-sound imaging

Patrick Stähli, Maju Kuriakose, Martin Frenz*, Michael Jaeger

Institute of Applied Physics, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland

June 24, 2019

*Corresponding author: martin.frenz@iap.unibe.ch

1 Introduction

Ultrasound (US) is an integral part of today’s clinical diagnostic practice since it provides us with real time display and flexible free-hand probe guidance while using non-ionizing radiation. The compact, portable and comparably inexpensive US systems make their use favourable for general practitioners, emergency units and bed-side care. On the downside, classical gray-scale B-mode US often suffers from low sensitivity and non-specific contrast, resulting in difficulties in differential diagnosis for certain disease types [1, 2, 3]. To complement B-mode images with additional structural and functional information in a single multi-modal system, much effort has been placed in recent years in developing new ultrasound-based modalities. This includes Doppler flow imaging, ultrasound elastography [4, 5, 6, 7, 8, 9] and optoacoustic imaging [10, 11, 12, 13].

Based on the dependence of the SoS on the tissues mechanical properties, SoS imaging is another promising modality that can help identifying disease-related changes of tissue composition and structure.

Breast ultrasound computed tomography (UCT) showcases the potential of SoS imaging on the example of breast cancer diagnosis [14]. In UCT, US is through-transmitted through the breast from many angles, and the spatial distribution of SoS is reconstructed based on analysing the detected US signals. In the straight ray approximation of ultrasound propagation, the time-of-flight (ToF) of pulsed signals is assigned to line integrals of slowness (inverse of SoS) along straight lines connecting sender and receiver locations, linking the ToF data to a Radon transform of the slowness distribution [15, 16, 17, 18]. This allows a fast linear reconstruction of the SoS, e.g. via the filtered backprojection, but with the disadvantage of a low spatial resolution due to diffraction and refraction. Refraction can be accounted for in an iterative way in a bent-ray ToF approach. Diffraction, however, requires methods that use the full signal as opposed to the ToF. Diffraction tomography based on the 1st order Born approximation is linear and fast, but limited to low contrast SoS variations so that its application to the breast requires a SoS prior with "reasonable low resolution" [19]. The best in-vivo images were so far obtained using non-linear iterative full-wave inversion schemes, but with the disadvantage of a high numerical cost [20, 21].

Whereas UCT is able to achieve high resolution and quantitative SoS images of breast cancer, the big disadvantage is that it is a through-transmission technique, which limits its use to acoustically transparent parts of the human body. To provide SoS imaging with all the flexibility of conventional handheld US and without the need for specialized equipment, SoS imaging must be based on echo US for a one-sided access. This would allow to not only image the SoS in transparent parts of the body, but also in any other part that is routinely examined using echo US, e.g. for the diagnosis of cancer other than in the breast, of fatty liver diseases or for the assessment of plaque composition inside the carotid artery.

Reflection-mode SoS imaging is feasible when taking into account a multiple-scattering process beyond the 1st order Born approximation. This problem can in principle be solved as in UCT in a non-linear full-wave inversion. Even though promising theoretical results were obtained using this approach in 2D digital breast phantoms [22], no clinical results have been presented so far, potentially because it requires low frequencies outside the bandwidth of conventional clinical US probes. Various alternative approaches were investigated that reconstruct the SoS based on pulse-echo signals. Techniques that estimate the average SoS between the transducer and the focal depth reported accurate SoS measurements of uniform tissue [23, 24, 25]. These methods, however, have low accuracy in the presence of SoS inhomogeneities. The crossed-beam tracking method is based on measuring the round-trip time through the intersection of two scanning beams (one for transmit, one for receive), and provides a spatial resolution on the order of 10 mm [26, 27]. Recently, an approach was proposed that reconstructs the local SoS based on axial variations of the average SoS determined by optimizing the transmit and receive focusing. The technique was verified in phantoms without lateral SoS variations, where it provided an axial resolution of about 7 mm [28].

We recently developed a reflection-mode technique, named computed ultrasound tomography in echo-mode (CUTE), which allows a real time determination of the SoS with promising spatial and contrast resolution [29, 30]. CUTE is based on analysing the spatial distribution of the echo phase in beamformed (using e.g. conventional delay-and-sum algorithm) radio-frequency (rf) mode US images. A deviation of the true SoS from the value assumed for beamforming results in a mismatch between the anticipated and the actual round-trip time of US propagation (henceforth termed ‘aberration delay’). A changing value of the aberration delay when detecting echoes under varying angles of ultrasound transmission and/or reception consequently results in a phase shift of these echoes, which is quantified in a spatially resolved way by e.g. Loupas type phase correlation[31]. This concept is closely related to approaches that have been developed for aberration correction of superficial SoS variations that act like a phase-distorting screen directly in front of the transducer array. They are based on analysing the differential echo phase as function of transducer element position in the channel data [32, 33, 34, 35, 36, 37]. CUTE goes beyond these approaches: By determining the phase shift in the beamformed images as opposed to the channel data, lateral resolution of phase shift data is achieved also away from the aperture. Based on a model of how the spatial distribution of SoS relates to the spatially resolved phase shift, the former can thus be reconstructed with a spatial resolution on the order of a couple of mm. Apart from serving as a source of diagnostic information, knowledge of the spatial distribution of SoS allows for aberration correction beyond the phase screen assumption [38].

The core of CUTE is the forward model that relates the SoS to echo phase shift, and this forward model is closely linked to the data acquisition scheme. Previously, we proposed an implementation of CUTE where the rf-mode images were acquired under a variety of transmit (Tx) angles, whereas the echoes were detected with a constant receive (Rx) aperture [29]. The model (henceforth termed ’old model’) was based on following key assumptions:

Because the Rx aperture is constant, the echo phase shift is entirely determined by the changing aberration delay along the changing Tx propagation directions. 2. 2.

The echo phase shift is proportional to the difference in aberration delay along different round-trip paths.

In [29, 30, 38], this forward model was formulated in the frequency-domain (FD). The FD formulation, however, cannot account for the absence of phase shift data in regions of low echo intensity, which leads to artefacts in the SoS image. Therefore we proposed that quantitative imaging requires a space-domain (SD) approach [39]. Sanabria et al. verified in a simulation study the advantage of using a SD instead of a FD approach [40]. However, in spite of promising results when using the SD approach, we realized that the old model leads to inconsistent results, depending on the sample geometry.

Here we present a new, more general model that solves this shortcoming by modifying the two key assumptions of the old model as follows:

Even if the Rx aperture is constant while changing the Tx angle, the echo phase shift contains the influence of a virtually changing Rx angle. To avoid ambiguities, we propose to switch from a pure Tx-steering approach to simultaneously steering both, the Tx- and the Rx- angles, around a variety of common mid-angles, in an approach similar to the common mid-point gather that has been used in seismic imaging [41] as well as in US [35, 36, 37]. 2. 2.

The a priori unknown SoS distribution leads to an error of the reconstructed position of echoes. This error not only depends on the aberration delay, but independently also on the values of the steering angles. As a result, the initially assumed simple proportionality between the aberration delay and phase shift is no longer valid.

In the theory section, we present a theoretical revision of the old model and the development of the new model. In a phantom study, we then compare the old and the new model in layered structures mimicking the abdominal wall and liver in preparation of one of our target clinical applications, the diagnosis of liver disease. The presented results reveal that the new model results in substantially improved quantitative SoS imaging among different layered phantom geometries. The article concludes with a discussion of the implications of the new model for the clinical use of CUTE, corroborated by a volunteer result that demonstrates the promise of the new model for a quantitative diagnosis of liver SoS.

2 Theory

The theory was developed assuming a linear array probe and a 2D geometry of SoS distribution and sound propagation, but can be readily adapted to curved arrays for 2D sector scans, as well as to matrix arrays for 3D imaging. Fig. 1(a) illustrates the measurement geometry. A linear US probe (parallel to coordinate x) is placed on top of the tissue (at coordinate z = 0), and transmits ultrasound pulses using a variety of different transmit settings, indexed by $n$. This can either be a variety of different Tx angles as proposed in earlier studies [29, 30], or a sequential activation of elements or groups of elements. For the sake of visual clarity, a sequential activation of single elements is considered in Fig. 1(a). Two different elements (indexed by $j_{n}$ and $j_{n^{\prime}{}}$) transmit an US pulse, one at a time. The two divergent wave fronts reach a point $\textbf{r}^{\prime}=(x^{\prime}{},z^{\prime}{})$ inside the tissue from two different angles $\phi_{n}$ and $\phi_{n^{\prime}{}}$. Intrinsic ultrasound reflectors located in a vicinity around $\textbf{r}^{\prime}{}$ lead to echoes that propagate back to the surface where they are detected by the different sensor elements (index $k$) of the array probe resulting in radio frequency (rf) signals $s(t,n,k)$ with the time indexed by $t$.

In line with previous studies [29, 30, 38], we propose that a separate rf-mode image $u(\textbf{r},n)$ is reconstructed from the detected signals for each Tx setting $n$. The images are reconstructed from the complex (analytic) rf-signal (crf-signal), where the imaginary part of the amplitude is the Hilbert transform (along $t$) of the initial real-valued rf-signal. The complex signal generated by a point reflector located at $\textbf{r}^{\prime}$ can be modelled as a product of a complex exponential carrier modulated with a complex-valued envelope $G$:

[TABLE]

where $f_{0}$ is the centre frequency, and $t_{0}$ is the actual round-trip time of the echo. The function $G$ characterises the spatio-temporal impulse response of the system, and depends on $\textbf{r}^{\prime}$ as well as on the receiving element $k$. The dependence on $k$ takes into account that the signal may be influenced by the angle-dependent receive response of the elements, but also by spatial variations of SoS. The reconstructed crf-amplitude in a point $\textbf{r}=\textbf{r}^{\prime}$ when using delay-and-sum beamforming (DAS) is:

[TABLE]

where $\hat{t}_{0}$ is the anticipated round-trip time for point r. Deviations of the true SoS $c(\textbf{r}^{\prime})$ from the anticipated value $\hat{c}(\textbf{r}^{\prime})$ lead to a deviation of the actual round-trip time from the anticipated value. This deviation, termed aberration delay, consists of two parts: one is the delay $\tau_{tx}(\textbf{r}^{\prime},n)$ of the transmitted wave front when propagating towards $\textbf{r}^{\prime}$, and the other one is the delay $\tau_{rx}(\textbf{r}^{\prime},k)$ of echoes propagating from $\textbf{r}^{\prime}$ to element $k$, so that $t_{0}=\hat{t}_{0}+\tau_{tx}+\tau_{rx}$. With the aberration delays, Eq. 2 ($\textbf{r}^{\prime}$ and r are omitted for notational simplicity) becomes:

[TABLE]

The old model

The old model was based on the following rationale: With the bandlimited frequency response of typical clinical US probes, $G$ varies ‘slowly’ compared to the oscillations of the exponential carrier. For a sufficiently small difference between Tx angles $\phi_{n}$ and $\phi_{n}^{\prime}$, the change in $\tau_{tx}$ will be small compared to the temporal variation of $G$, so that the value of $G$ can be assumed constant and $\tau_{tx}$ can be replaced by its average $\overline{\tau_{tx}}$ in the envelope. Eq. 3 is then simplified to:

[TABLE]

On the right-hand side of Eq. 4, only the complex pre-factor depends on $n$. In the old model, the first processing step of CUTE is the experimental determination of the change in this pre-factor as a function of $n$. This is achieved by determining a map of local echo phase shift $\Theta(\textbf{r}^{\prime},n,n^{\prime})$, defined as the phase angle of the (locally averaged) point-wise Hermitian product between the crf-mode images obtained with Tx settings $n$ and $n^{\prime}$:

[TABLE]

Averaging over a limited area around $\textbf{r}^{\prime}$ (‘tracking kernel’), with size $\Delta x$ by $\Delta z$, improves robustness of the phase determination. The size of the tracking kernel defines the trade-off between spatial and contrast resolution of CUTE. Provided that the change in $\tau_{tx}$ is smaller than half the oscillation period (to avoid phase aliasing), the measured echo phase shift is - according to Eq. 4 - related to the Tx aberration delay $\tau_{tx}(\textbf{r}^{\prime},n)$, as:

[TABLE]

Similar to UCT, including US diffraction and refraction in the forward model is in principle feasible, but would result in a time-consuming iterative reconstruction. In view of real time SoS imaging, we thus adhere to the straight ray approximation and relate $\tau_{tx}(\textbf{r}^{\prime},n)$ to line integrals of slowness deviation $\Delta\sigma(\textbf{r},n)$, from Tx elements $j_{n}$ to points $\textbf{r}^{\prime}$.

[TABLE]

By inverting the forward model consisting of Eq. 6 and Eq. 7 the SoS $c$ is reconstructed from the measurements $\Delta\Theta$.

New model

As mentioned in the Introduction, two fundamental changes to the CUTE methodology lead to the new model: In a first step, the common mid-angle (CMA) approach is described. Based on that, the new model is developed in a second step, which takes into account the position error of the reconstructed echoes.

**(a) Common mid-angle tracking

**An important prerequisite to the old model expressed in Eq. 6 was the assumption that $G$ did not depend on $n$. This assumption is plausible when considering a point reflector. In a more general situation, however, the relation between the aberration delay and the measured echo phase shift becomes ambiguous. This is illustrated for two extreme (but commonly found) cases:

• •

Specular reflector: The echo from a reflector intersecting with $\textbf{r}^{\prime}$ propagates back to the probe along a main direction $\psi_{n}$ – determined by the mirror law – that varies when varying $\phi_{n}$.

• •

Uniform diffuse scattering: The echoes from a dense (below spatial resolution cell) and random distribution of ‘point’ reflectors interfere at the aperture. The signal detected by an element $k$ decorrelates with changing $n$, because the relative round-trip times of the different echoes change in a different way. The decorrelation rate depends on how many echoes interfere, i.e. on the a priori unknown (due to aberrations) Tx focus diameter at $\textbf{r}^{\prime}$.

In both cases – and also in a more general sense – the envelope function $G$ in Eq. 1 depends on n. In the same way as for the Tx angles, we assume sufficiently small differences between Rx angles so that the change in $\tau_{rx}$ will be small compared to the oscillation period. Then $\tau_{rx}$ can be replaced by its average $\overline{\tau_{rx}}$ in the envelope and Eq. 4 can be re-written:

[TABLE]

Because $G$ varies with changing $n$, the sum in Eq. 8 depends on n. Thus, the phase shift according to Eq. 5 not only depends on the Tx aberration delay but also on the Rx aberration delay in the complex exponential in the right-hand side of Eq. 8. To avoid ambiguities, we introduce a fundamental change to the CUTE methodology: common mid-angle (CMA) tracking. In this technique, crf-mode images are reconstructed for pairs of Tx angles $\phi_{n}$ and Rx angles $\psi_{m}$ that have the same mid-angle. One has to keep in mind that Rx beamsteering goes hand-in-hand with reducing the Rx angular aperture, at the cost of a reduced lateral resolution of the crf-mode image. CUTE is based on evaluating echo phase shifts in the beamformed crf-mode image as opposed to the channel data, specifically because this allows to spatially resolve the influence of the SoS on the echo phase shift of different echoes. When reconstructing images with sharp $\phi_{n}$ and $\psi_{m}$, no lateral resolution would be obtained in the crf-mode images. Therefore, to achieve lateral resolution, images have to be reconstructed with non-zero Tx and Rx angular apertures, which in combination define the size of a spatial resolution cell. Even though, experimentally, angles have to be defined with non-zero angular aperture, we continue assuming that they are sharply defined for the remainder of this chapter. Later on, we will discuss the practical implication of the non-zero aperture on SoS reconstruction. Similar to the common mid-point approach that was previously used in aberration correction [38], the CMA approach makes use of signal redundancy: Signals obtained under angle pairs $(\phi_{n},\psi_{m})$ grouped around the same mid-angle $\gamma=0.5\left(\phi_{n}+\psi_{m}\right)$ are well correlated, independent of the spatial distribution of US reflectors. For illustration, Fig. 1(b) shows a scale-up of a spatial resolution cell centred around point $\textbf{r}^{\prime}$, containing a number of reflectors. Each reflector is the pivot of an isochrone, a line that connects all points that would lead – for a specific combination of $(\phi_{n},\psi_{m})$ - to the detection of an echo at the same time. When changing $\phi_{n}$ to $\phi_{n^{\prime}}$ without changing $\psi_{m}$ (Fig. 1(c)), the distances between the isochrones change, leading to decorrelation of the signal that is detected by elements located in direction $\psi_{m}$. When changing $\psi_{m}$ in opposite direction (Fig 1(d)), to $\psi_{m}^{\prime}$ such that the mid-angle in the pairs $(\phi_{n},\psi_{m})$ and $(\phi_{n^{\prime}},\psi_{m^{\prime}})$ is fixed (so that $\gamma=0.5(\phi_{n}+\psi_{m})=0.5(\phi_{n^{\prime}}+\psi_{m^{\prime}}))$, then the distances between isochrones is not altered.

We thus assume that crf-mode images $u(\textbf{r},n,m)$ are reconstructed where $n$ indexes the Tx setting resulting in a Tx angle $\phi_{n}$, $m$ indexes the mid-angle setting $\gamma_{m}$, and both together result in a receive angle $\psi_{(m,n)}$. Accordingly, Eq. 5 and Eq. 6 are adapted for the CMA tracking to:

[TABLE]

[TABLE]

**(b) Reconstructed echo position error model

**To derive the second fundamental change to the CUTE methodology, we introduce a geometric perspective of how aberration delays relate to the echo phase (Fig. 2).

We assume a reflector (point or plane) that is detected with a specific combination $(\phi_{n},\psi_{m,n})$. The total aberration delay (in Fig. 2(a) larger than zero) leads to a shift of the reconstructed position of the echo, away from the true location of the reflector, to the isochrone for which the anticipated round-trip time agrees with the actual round-trip time. The distance $d$ between the location of the reflector and the isochrone is – according to the geometric considerations depicted in Fig. 2(a) – given by:

[TABLE]

To determine the phase shift, the spatial oscillation period of the reconstructed echo is also required. This follows the same argumentation (see Fig. 2(b)): The spatial period $\Lambda$ is given by the distance between two isochrones that are separated by one period in time:

[TABLE]

An important prerequisite of the old model was the evaluation of the echo phase shift at the true location of the ultrasound reflector (r was chosen equal to $\textbf{r}^{\prime}$).

However, it is not the phase shift at the true location of the reflector that determines the value of the measured phase shift, but the shift of the reconstructed position of the echo. As a second fundamental modification to the CUTE methodology, this is taken into account by the new model. Fig. 2(c) illustrates that the phase shift is given by the ratio between $\Delta d$, the shift of the echo centroid and $\Lambda$. Formally, when changing $\phi_{n}$ by $\Delta\phi$ to $\phi_{n^{\prime}}$ (so that $\Psi_{m,n}$ changes by $-\Delta\phi$ according to $\Psi_{m,n^{\prime}}$), $d$ changes accordingly by:

[TABLE]

To simplify the maths, we assume that the signals are bandpass filtered in a way that the spatial period $\Lambda$ does not change between combinations $(n,m)$ and $(n^{\prime},m^{\prime})$ (this can easily be achieved in practise) and takes, for example, the value:

[TABLE]

With this convention, the new model for the phase shift then is:

[TABLE]

In comparison to the old model, the new model (eq. 15) includes CMA tracking and a division of the aberration delay by the different cosines for the different ($n$, $n^{\prime}$).

3 Materials and Methods

Space domain forward and inverse model

For simplicity of formulation of the numerical models and without loss of generality, we have chosen a plane wave pulse-echo acquisition approach so that each pixel in the image area is insonified with the same discrete set $\{\phi_{n}\}$ of equidistantly spaced Tx angles. The equidistant spacing of the Tx angles allows choosing – for the CMA approach – a set $\{\gamma_{m}\}_{n}$ of common mid-angles that result in a set $\{\psi_{m}\}_{n}$ of Rx angles that is identical to $\{\phi_{n}\}$. This allows to make use of data redundancies as explained below. For clarity, we introduce a new notation for the angle set, $\{\Phi_{l}\}$, with $\phi_{n}\in\{\Phi_{l}\}$ and $\psi_{m,n}\in\{\Phi_{l}\}$. For this study, $\{\Phi_{l}\}$ was chosen to be [-25∘, -15∘, -5∘, 5∘, 15∘, 25∘] for $l$ = [1 .. 6]. Table 1 summarizes the resulting combinations of $\phi$, $\psi$ and $\gamma$. As described in the Theory section, tracking can only be performed between angle combinations having the same mid-angle $\gamma$, indicated in table 1 by fields having the same gray value. Since only one combination of $\phi$, $\psi$ exists for both common mid angles $\gamma=\pm 25^{\circ}$, no tracking was performed for these $\gamma$.

Choosing identical sets of Tx and Rx angles allows making use of data redundancy: In an ideal case, the crf-mode images resulting from interchangeable $\left(\phi_{n}\mid\psi_{m,n}\right)$ angle pairs $\left(\Phi_{l}\mid\Phi_{l^{\prime}}\right)$ and $\left(\Phi_{l^{\prime}}\mid\Phi_{l}\right)$ are identical by time reversal symmetry. The phase correlation between $\left(\Phi_{l}\mid\Phi_{l^{\prime}}\right)$ and $\left(\Phi_{l+1}\mid\Phi_{l^{\prime}-1}\right)$ and the one between $\left(\Phi_{l}\mid\Phi_{l^{\prime}}\right)$ and $\left(\Phi_{l^{\prime}-1}\mid\Phi_{l+1}\right)$ thus provide redundant data. In a more realistic case, slight differences exists between Tx and Rx beamforming. This asymmetry is removed by averaging of the redundant phase correlations. In addition, the averaging results in phase shift values that are more robust against phase aliasing, and the reduced total number of phase shift maps allows a faster SoS reconstruction.

Depending on the Tx angle, the Tx wave front reaches only a part of the imaging plane, leaving areas of missing echoes (’Tx shadows’) towards the edges of the image. As mentioned in the Introduction, the previously proposed FD formulation cannot account for missing data regions, which complicates a quantitative SoS reconstruction. To avoid such drawbacks, the forward models in this study were implemented in space domain (SD). For this purpose, the measured spatial distribution of echo phase shift as well as the distribution of slowness to be reconstructed were discretised on the same 2D Cartesian grid, as

[TABLE]

Thereby, the number of nodes $(j,k)$ in $x-$ and $z-$ direction is $N_{j}$ and $N_{k}$, respectively. The nodes sample the chosen dimensions $(X,Z)$ of a rectangular image area (in our study 38.4 by 50 mm) with a spatial resolution of $(\Delta x,\Delta z)$ (in our study 0.96 by 1 mm).

The goal of the experimental part is to compare the old and the new model in a phantom study. Thus, both the old and the new model were implemented as described by Eq. 17 and 19, respectively, derived from Eq. 6, 7 and 15. However, to demonstrate that the improvement of the new model not only depends on the CMA tracking, we also implemented a forward model where only CMA tracking - but not the echo position error - is taken into account (Eq. 18).

• •

Old model

[TABLE]

• •

CMA model

[TABLE]

• •

New model

[TABLE]

Thereby the $w_{(j,k,j^{\prime},k^{\prime})}^{\Phi_{l}}$ are the integration weights defining the discrete line integrals along the US propagation direction with the angles $\Phi_{l}$.

Any of the forward models in Eq. 17 to 19 can be written in matrix notation:

[TABLE]

To reconstruct the values $\Delta\sigma_{j^{\prime},k^{\prime}}$ from the measurements $\Delta\theta_{j,k}$, a Tikhonov pseudo-inverse of $\boldsymbol{M}$ is used:

[TABLE]

which minimizes the expression

[TABLE]

Thereby, $\boldsymbol{D}_{x}$ and $\boldsymbol{D}_{z}$ are finite difference operators in $x$ and $z$, respectively, and $\gamma_{x}$ and $\gamma_{z}$ are regularisation parameters. Regularisation of the spatial gradient of the slowness deviation enforces a smooth slowness profile while letting the mean SoS vary freely.

Data acquisition system

For the experimental study, we used a Vantage 64 LE (Verasonics Inc., WA, USA). This research ultrasound system allows simultaneous ultrasound transmission on 128 channels and parallel digitisation of signals on 64 elements, with real time data transfer via a PCI Express link to a host computer for further processing. The system was connected with an L7-4 linear vascular probe (ATL Philips, WA, USA) for pulse-echo signal acquisition. This probe features 128 elements at a 0.29 mm pitch (resulting in 38.4 mm aperture length), and a bandwidth from 4 to 7 MHz with 5 MHz centre frequency. We implemented a dedicated scan sequence for the acquisition of plane wave pulse-echo data. As mentioned earlier, the target Tx and Rx steering angle step size was 10∘. Given the SoS contrast of the chosen phantoms, however, this step size resulted in phase aliasing due to aberration delay changes above half the oscillation period at the centre frequency. To avoid phase aliasing, we chose a tracking angle step size of 2∘. Moreover, for synthetic Tx focusing (via coherent plane wave compounding [42, 43]), data had to be acquired with an even smaller step size of 0.5∘. Plane wave data was thus acquired with Tx angles $\phi$ ranging from -27.5∘ to 27.5∘ in 0.5∘ steps. To use the full probe aperture on receive, echoes were recorded twice for each Tx angle, once on each of two 64-element sub-apertures.

CRF-mode image reconstruction

Crf-mode images were reconstructed off-line for each Tx angle from the raw crf data using a delay-and-sum (DAS) algorithm. Crf-mode amplitudes were calculated on a rectangular grid, with dimensions 38.4 mm in $x$ (aperture length) by 50 mm in $z$ direction and with spatial resolution $dx=0.29$ mm in $x$ (pitch) by $dz=0.037$ mm in z. At each grid node with coordinates $(x,z)$, the anticipated round-trip times $\hat{t}_{0}(x,z,\phi,k)$ to each receiving element k were calculated based on the anticipated SoS $\hat{c}$.

[TABLE]

The weights $\alpha(x,z,k)$ were designed to provide (within the limits of the aperture length) a receive angular aperture of $\pm$ 30∘.

For focused Tx steering along each specific target angle $\phi$, crf-mode images were coherently compounded over a Tx angle range $[\phi-2.5^{\circ}:0.5^{\circ}:\phi+2.5^{\circ}]$ centred around the respective target angle, resulting in a synthetic Tx focus with an angular aperture of 5∘. This was done for a target Tx angle range $\phi\in[-25^{\circ}:2^{\circ}:25^{\circ}]$.

For the CMA approach, each crf-mode image was transformed to a sequence of crf-mode images for different receive steering angles $\Psi\in\left[-25^{\circ}:2^{\circ}:25^{\circ}\right]$. This was achieved by spatial filtering of the already reconstructed crf-mode images. For this purpose, the 2D discrete Fourier transform was first calculated. Then the spatial frequency spectrum was multiplied with an angle-dependent weighting function to highlight a specific mid-angle, in that way implicitly determining the corresponding receive angle. The angular radius of the weighting function was chosen $\pm$ 1.25∘, resulting in an implicit Rx angular aperture radius of 2.5∘. This approach was chosen, as opposed to an explicit Rx beamsteering, to avoid abundant memory usage and it was implemented as part of the tracking algorithm.

Phase tracking

The echo phase shift was determined according to Eq. 6 (old model) and Eq. 9 (CMA approach). The tracking kernel size was chosen $\Delta x$ = 2 mm and $\Delta z$ = 2 mm. Echo phase tracking was performed with the 2∘ angle steps and phase shift maps were accumulated over successive angle steps to obtain the echo phase shift over 10∘ steps between Tx angles $\Phi_{l}$ and $\Phi_{l+1}$ (conventional tracking) and Tx/Rx angle pairs $(\Phi_{l}|\Phi_{l^{\prime}})$ and $(\Phi_{l+1}|\Phi_{l^{\prime}-1})$ (CMA approach).

SoS reconstruction implementation details

For SoS reconstruction, the phase shift maps were downsampled onto a Cartesian grid with 40 $(x)$ by 50 $(z)$ pixels covering 38.4 mm $(x)$ by 50 mm $(z)$. The forward models according to Eq. 20 and 21 were correspondingly formulated for a SoS distribution sampled on the same grid, resulting in square system matrices with (40$\cdot$50) by (40$\cdot$50) elements per angle step. The reconstructed SoS images, however, showed artefacts in the top 5 mm due to element cross talk. Further, artefacts occurred towards the lower edge of the image, independent of the initial choice of the depth range (reason yet unclear). To mask out these artefacts, the SoS images were cut in z-direction to a depth range from 5 to 40 mm for the final presentation.

The regularization parameters $\gamma_{x}$ and $\gamma_{z}$ (see. Equ. 22) are subject to a trade-off between reducing artefacts (by enforcing a smooth slowness profile) and maximising spatial resolution. For this study, they were chosen so as to clearly distinguish the different phantom compartments.

Phantom design and materials

For this study, three different phantoms where designed, each representing as close as possible a particular geometry of the anatomical structure of the abdominal wall and the liver (see Fig. 3).

1. (a)

Two layer phantom: This phantom contained two horizontal (parallel to $x$) layers, mimicking a single fat layer (F1) on top of liver tissue (L). The special feature of this type is the complete absence of lateral SoS variations but a pronounced axial variation. This geometry occurs when imaging the liver sagittally through the linea alba. 2. (b)

Four layer phantom: In comparison to the two layer phantom, two additional layers were added, mimicking the rectus abdominis muscle (M) and the post peritoneal fat layer (F2). The purpose of this phantom was to evaluate the ability of the different models in resolving axial SoS variations of thin layers, representative for a sagittal section lateral to the linea alba. 3. (c)

Laterally varying muscle diameter (LVMD) phantom: In contrast to the four layer phantom, the rectus abdominis (M) layer deviates in this phantom from the parallel layer structure and consists of two wedge-shaped areas, as when imaging the liver in a transverse section through the linea alba.

The phantom components were based on oil-in-gelatin emulsions [44, 45, 46, 47]. In a first step, a gelatin base solution was prepared by dissolving 20 wt% porcine gelatin (Geistlich Spezial Gelatine, health and life AG, Switzerland) in 75∘C H2O and adding 2 wt% cellulose (Sigmacell Cellulose Type 20, Sigma Aldrich, Switzerland) to provide uniform diffuse echogenicity. This aqueous gelatin base solution had a SoS of 1555 m/s and was used to mimic muscle (M) and liver (L). To mimic fat, medium-chain triglycerides oil (Ceres-MCT Oil, Puravita, Switzerland) (SoS = 1350 m/s) was slowly blended under continuous stirring into the aqueous gelatine base solution using a Visco Jet cone-stirrer (VISCO JET Agitation Systems, Germany). During this process, small oil droplets were formed and captured via hydrophobic interaction by the lipophilic part of the gelatine strings. After cooling the emulsion, the oil droplets were trapped within the gelatine matrix. The resulting SoS was determined by the emulsion’s relative MCT oil weight content, in this study 0.65, resulting in a SoS of 1420 m/s. The mentioned SoS values were determined using a through-transmission time-of-fight set-up (accuracy $\pm$ 5 m/s) and will serve in the following as reference values for the SoS images. Note that, even though the absolute SoS values may deviate from real tissue depending on literature references, we have taken care to chose a representative SoS contrast.

4 Phantom Results

Besides the old and the new model (Eq. 17 and 19), an intermediate model was implemented where only the CMA tracking was considered (CMA model, Eq. 18). The purpose of the CMA model was to assess the influence of the CMA tracking on the results separately from the consideration of the echo position error.

The spatial distributions of the reconstructed SoS are showed in color-code (Fig. 4 - 6 (rows)). The color scale was chosen in a trade-off between covering a large SoS range and allowing to perceive small SoS variations (e.g artefacts). As a consequence, image regions that deviate $\gtrsim$ 70 m/s from the true upper and lower SoS appear saturated.

As described in Materials and Methods, the reconstruction of the crf-mode images was based on assuming a uniform SoS value $\hat{c}$. To examine the influence of $\hat{c}$ on the final SoS result, crf-mode images were reconstructed with four different $\hat{c}$ covering the range of SoS values between fat-mimicking material (1420 m/s) and muscle/liver tissue mimicking material (1555 m/s) (Fig. 4 - 6 columns). The rightmost columns show the mean and standard deviation of the reconstructed SoS within each individual compartment.

**Old model

**The SoS images reconstructed with the old model (see Fig. 4 to 6 - top rows) have in common that they all show artefacts at the top corners, independent of the phantom as well as the a priori SoS $\hat{c}$.

In the two layer phantom, the SoS of the fat mimicking compartment (F1) strongly depends on the a priori SoS $\hat{c}$. Whereas a reasonable SoS - apart form the artefacts - is reconstructed for $\hat{c}$ of 1420 m/s, higher $\hat{c}$ lead to a distinct underestimation of the SoS in said compartment. In contrast, the mean SoS of the liver mimicking compartment (L) varies only little when changing $\hat{c}$. However, except for $\hat{c}$ of 1465 m/s, the images contain a high level of artefact inside the liver mimicking compartment, suggesting unrealistic SoS inhomogeneities.

In the four layer phantom, the mean SoS of the muscle mimicking compartment (M) is reconstructed best for a $\hat{c}$ of 1420 m/s. Higher $\hat{c}$ lead to strong artefacts at the right side of the image, resulting in a decreased mean SoS. Despite these artefacts, however, an increasing $\hat{c}$ leads to an increase in SoS in the middle region of the muscle mimicking compartment. A similar effect can be observed for the post peritoneal fat mimicking compartment (F2), where an increasing $\hat{c}$ leads to a decrease in SoS in the middle region of the layer, but an increase in artefacts at the right side of the image. Inside the liver mimicking compartment, the mean SoS varies between 1541$\pm$34 m/s for a $\hat{c}$ of 1420 m/s and 1571$\pm$40 m/s for a $\hat{c}$ of 1555 m/s. Similar to the two layer phantom, a $\hat{c}$ of 1420 m/s leads to a decrease in SoS towards the bottom of the image, whereas higher $\hat{c}$ lead to artefacts of high SoS at the bottom of the image.

In case of the LVMD phantom, all $\hat{c}$ result in high level of artefact around the wedge-shaped muscle mimicking compartment, making a clear identification difficult. Moreover, regardless of $\hat{c}$, the old model is unable to reconstruct the true SoS of the liver mimicking compartment.

**CMA model

**SoS images reconstructed with the CMA model (Fig. 4 to 6 - middle rows) show reasonable mean SoS for the fat mimicking compartments (F1 & F2) in all phantoms when $\hat{c}$ is close to 1420 m/s. However, independent of the phantom, an increasing $\hat{c}$ leads to a distinct decrease of the reconstructed SoS in any fat mimicking compartment.

In the two layer phantom, only a $\hat{c}$ of 1465 m/s leads to reasonable mean SoS in the liver mimicking compartment. Similar to the old model, also in the CMA approach, a $\hat{c}$ of 1420 m/s leads to a distortion of SoS towards the bottom of the image.

In the four layer phantom, the muscles mimicking compartments mean SoS is overestimated for low $\hat{c}$, whereas high $\hat{c}$ results in an underestimated mean SoS. In contrast, inside the liver mimicking compartment, the mean SoS is underestimated for any $\hat{c}$. Furthermore, similar artefacts toward the bottom of the image as in the two layer phantom are found also in the four layer phantom.

Contrary to the old model, the wedge-shaped muscles in the LVMD phantom can be clearly perceived when the reconstruction is performed with the CMA approach. The best agreement of the muscles mimicking compartments SoS with the reference value is found for a $\hat{c}$ of 1420 m/s. An increase in $\hat{c}$, however, leads to a decrease of the reconstructed SoS of the wedge-shaped muscles. Analogous to the old model, the CMA model is, independent of $\hat{c}$, unable to reconstruct the liver mimicking compartments SoS in agreement with the reference value.

**New model

**In comparison to the old and the CMA model, the new model leads to substantially improved SoS images among all phantoms, where the best results are achieved when $\hat{c}$ is close to the true SoS of the first fat mimicking layer (F1), i.e. 1420 m/s. For such a $\hat{c}$, the mean SoS of the fat (F1) as well as the liver (L) mimicking compartment deviate only little from the reference value in any phantom. The SoS of the muscle mimicking compartment (M) in the four layer phantom is reconstructed in agreement to the reference value, although only in the middle region of the image. Toward the lateral edges, a decrease in SoS is observed, leading to a decreased mean SoS. In the LVMD phantom, the mean SoS of the wedge-shaped muscle is underestimated because of the limited spatial resolution of the reconstructed images.

Whereas SoS images reconstructed with the old and the CMA model strongly depend on the a priori SoS $\hat{c}$, the new model shows a higher stability against variations of the SoS $\hat{c}$. However, also the new model shows a slight increase in the mean SoS when $\hat{c}$ is increased, most pronounced in the top fat mimicking compartment(F1). Furthermore, all images have in common that higher $\hat{c}$ leads to a higher level of artefacts, leading to an increased standard deviation in the individual compartments.

5 Discussion and Conclusion

Our phantom study revealed two main limitations of the old model in layered samples: First, the absolute SoS values as well as the SoS contrast between different phantom compartments strongly depend on the a priori SoS $\hat{c}$. Secondly, in phantoms with lateral SoS variations, i.e. the LVMD phantom, the old CUTE model completely fails in reconstructing the SoS of the different phantom compartments.

Two fundamental changes to the CUTE methodology presented in this study lead to the new model that solves these shortcomings. CMA tracking, the first change, determines the echo phase shift not only by a pure Tx-steering approach but in a simultaneous steering of both, the Tx- and the Rx- angle, around a variety of common mid angles.

The results showed that when only CMA tracking is implemented, an improvement in the level of artefacts is observed when imaging samples with laterals SoS variations. Thus, whereas in the SoS images determined with the old model the wedge-shaped muscles of the LVMD phantom could not be perceived, they are clearly visible in the SoS images reconstructed with the CMA model. This confirms that there is an implicit Rx-dependence of the aberration delay $\Delta\tau$ in the phase shift, which was not considered in the old model. However, similar to the old model, also in the CMA model the absolute SoS values as well as the SoS contrast between the different phantom compartments depend strongly on $\hat{c}$. This reveals that the strong dependency of the SoS images on $\hat{c}$ does not origin from a dependency of the receive beamforming on the erroneously assumed background SoS. Anyhow, also the CMA model failed completely in reconstructing the liver’s SoS in the LVMD phantom.

In addition to the CMA approach, the new model takes into account the phase shift that originates from the erroneous reconstructed position of the echoes. This leads, compared to the other models, to substantially improved quantitative SoS reconstructions among all phantoms. This is most pronounced in the LVMD phantom, where only the new model is able to reconstruct the SoS of the liver compartment close to the reference value, best when $\hat{c}$ is close to the true SoS of the top layer.

Even though also the new model shows a slight influence of $\hat{c}$ on the mean SoS of the phantom compartments, it is small compared to the old and the CMA model. We hypothesize that this influence originates from experimental conditions that were not yet accounted for, such as e.g. near-field effects or the 2D approximation of the SoS distribution and sound propagation. Further, an increased level of artefacts in the SoS images is observed when $\hat{c}$ strongly deviates from the SoS of the first layer. This observation is assigned to increased phase noise in the echo phase shift maps, caused by the reduced Tx and Rx focus quality.

One important assumption made in this study is the straight ray approximation of US propagation neglecting diffraction and refraction. Given the lateral resolution of the final SoS image (few mm) in relation to the wavelength (0.2 mm at 5 MHz), diffraction will have played a minor role. With the realistically high SoS contrast of the presented phantoms, however, one would expect that refraction cannot be neglected. None-the-less, our results demonstrate that the straight ray approximation was still reasonable in the sense that correct quantitative results could be obtained. Refraction will have a more noteable influence on the reconstructed SoS images in samples with different geometry, e.g spherical structures. Refraction could then be compensated for e.g. using a more accurate forward model, such as an iterative bent-ray approach or the Rytov approximation.

Another fundamental assumption is that of sharply defined Tx/Rx angles, but finite angular apertures are needed and were used in the experiments to enable a laterally resolved echo phase shift measurement (and thus SoS image). One may argue that the round-trip time of echoes thus correspond to average slowness integrals over wedge-shaped areas rather than along thin lines. Our results demonstrate that the simplification to line integrals had only minor influence (if any) on the quantitativeness of the SoS. This can be understood from two perspectives: in the first one, the reconstruction of local SoS is based on the local variation of echo phase shift rather than the echo phase shift itself [29]. Therefore, independent on whether echo phase is influenced by the slowness over a broad or a narrow beam area above the Tx/Rx focus, at the focus it is determined by variations of slowness on the spatial scale of the spatial resolution of the focus. From the second perspective, the angular aperture used for Tx/Rx focusing (5∘) together with the wavelength (0.2 mm at 5 MHz) result in a focal width of about 3 mm. The active probe aperture needed to focus to the lower edge of the image is about 3.5 mm and thus only slightly larger than the focal width, indicating collimated beams rather than focused ones. At the lateral resolution of echo phase shift detection given by the focal width, the rays can therefore be approximated as lines.

The focus of the present study was to validate the new model. So far we did not attempt to investigate the clinical accuracy of CUTE. We would like to point out, however, that the liver phantoms used in this study were specifically chosen to provide a realistic SoS contrast between different compartments as well as realistic spatial distribution of different compartments. Thus, our liver phantom results suggest that CUTE can achieve sufficient accuracy for a diagnostic application: the new model is able to reconstruct the SoS of the liver mimicking tissue with a deviation from the true SoS of only a few m/s with a standard error of $\approx$ 2.5 m/s in the best and $\approx$ 10 m/s in the worst case. Imbault et al. [25] reported that the average SoS of the liver of healthy patients is $\approx$ 1570 m/s. For patients having steatosis of grade 1 (according to the Brunt scale), the average SoS decreases by about 40 m/s. Thus, CUTE (with the proposed new model) has potential for a quantitative diagnosis of fatty liver disease by complementing real-time b-mode US in a single hand-held device. To illustrate that CUTE can be readily applied to an in-vivo case, Fig. 7 shows the B-mode and SoS images reconstructed with the old and new model of a healthy volunteer’s liver. In comparison to the old model, the new model reflects the true anatomical situation and shows a more homogeneous SoS distribution within the liver, as expected in healthy tissue. Moreover, the reconstructed mean SoS of the liver is 1564$\pm$4 m/s and thus reflects healthy liver tissue [48].

Transabdominal hepatic imaging in practise is conventionally performed using curvilinear probes operating at lower frequencies. Whereas this study used a linear probe, the methodology can readily be adapted to curvilinear probes. The advantage of a linear compared to a curvilinear probe is the larger angle range that can be used for Tx/Rx beamsteering, as the same pixel is within the "view" of a smaller number of elements when these elements are pointing into diverging directions. Also, the ratio between the probe aperture length and the image size limits the available angle range in full liver imaging for points that are located far away from the probe. This leads to reduced axial resolution of the SoS image with increasing depth. In situations where one is interested in an average quantitative reading rather than the spatial distribution of SoS, the low axial resolution may not be important. For differential diagnosis of a spatially confined tumor that is located deep inside the liver but well visible on B-mode US, spatial resolution could be enforced in a prior-guided SoS reconstruction in order to regain quantitative results. In situations where the lack of spatial resolution becomes a problem, using a very large aperture linear array probe instead of a curved one could be a viable solution.

In summary, the new model is a key requirement for quantitative handheld reflection-mode SoS imaging. Thanks to the fast angle scanning of US systems and the low computational cost of CUTE, quantitative SoS images can thus be routinely displayed in parallel to conventional B-mode US. CUTE therefore can be used to image any organ that can be examined using echo US, for example the female breast to improve breast cancer diagnostics. In situations where the echoes are not static, e.g. inside the large blood vessels, blood flow and tissue clutter inhibit phase tracking and may cause SoS artefacts in the reconstruction. To overcome this limitation, we proposed a technique [49] where the first order echoes from the moving blood cells are separated from the static echo clutter via a clutter wall filter as in Doppler flow imaging. This technique extends the application of CUTE also to arteries, e.g. for the assessment of plaque inside the carotid artery and may also further improve liver imaging, since it can reduce SoS artefacts which may occur around hepatic veins. Furthermore, knowing the spatial distribution of SoS also allows reconstructing conventional B-mode US images with reduced aberration artefacts [38].

6 Acknowledgement

This research was funded in part by the Swiss National Science Foundation (project number 205320_178038) and the European Union’s Horizon 2020 research and innovation programme under grant agreement No 732411, Photonics Private Public Partnership, and is supported by the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number 16.0162. The opinions expressed and arguments employed herein do not necessarily reflect the official view of the Swiss Government. The authors thank René Nyffenegger for his technical assistance.