Full correction for spatially distributed speed-of-sound in echo ultrasound based on measuring aberration delays via transmit beam steering

Aberrations of the acoustic wave front, caused by spatial variations of the speed-of-sound, are a main limiting factor to the diagnostic power of medical ultrasound imaging. If not accounted for, aberrations result in low resolution and increased side lobe level, over all reducing contrast in deep tissue imaging. Various techniques have been proposed for quantifying aberrations by analysing the arrival time of coherent echoes from so-called guide stars or beacons. In situations where a guide star is missing, aperture-based techniques may give ambiguous results. Moreover, they are conceptually focused on aberrators that can be approximated as a phase screen in front of the probe. We propose a novel technique, where the effect of aberration is detected in the reconstructed image as opposed to the aperture data. The varying local echo phase when changing the transmit beam steering angle directly reflects the varying arrival time of the transmit wave front. This allows sensing the angle-dependent aberration delay in a spatially resolved way, and thus aberration correction for a spatially distributed volume aberrator. In phantoms containing a cylindrical aberrator, we achieved location-independent diffraction-limited resolution as well as accurate display of echo location based on reconstructing the speed-of-sound spatially resolved. First successful volunteer results confirm the clinical potential of the proposed technique.


Introduction
Clinical ultrasound (US) imaging is becoming an inherently multimodal diagnostic tool, as the traditional dual display of echogenicity (B-mode) and blood flow (Doppler) is gradually complemented by novel modalities that provide diagnostically valuable information of tissue composition. Examples are strain imaging (Thomas et al 2006, Kato et al 2008, Hong et al 2009 or shear wave elastography (Nightingale et al 2002, Athanasiou et al 2010, Bavu et al 2011 of tissue stiffness, or optoacoustic/photoacoustic imaging (Manohar et al 2007, Heijblom et al 2012, Nuster et al 2013, Alles et al 2014, Jansen et al 2014, Wang et al 2014, Zalev et al 2014 of haemoglobin concentration, blood oxygenation, and lipid content. An important limiting factor to US-based modalities is the distortion of US waves when propagating through tissue with an unknown spatial distribution of speed-of-sound. Speed-of-sound varies by several percent depending on tissue composition  and condition , leading to acoustic wavefront distortions of both transmitted and reflected ultrasound , Freiburger et al 1992, Hinkelman et al 1994. In absence of any other prior information, US image reconstruction is typically based on assuming a homogeneous speed-of-sound equal to the average of soft tissue of e.g. 1540 m s −1 . The mismatch between assumed and actual distribution of speed-of-sound results in degradation of transmit (Tx) and receive (Rx) focusing, limiting resolution and contrast, and acoustic refraction leads to geometric distortions (Austin et al 2004, Spieker et al 2004 if not accounted for.
Great effort has therefore been dedicated to the development of techniques to mitigate the image degradation caused by acoustic aberrations (aberration correction, AC). A significant part of such techniques were built around the assumption that the inhomogeneous speed-ofsound acts like a phase screen in front of the US probe. Measurement of the different arrival phases of strong coherent echoes at the individual transducer elements along the receive aperture allows determination of this phase screen under the assumption that such echoes are generated by point-like ultrasound reflectors, so-called guide stars or beacons (Flax and O'Donnel 1988, Nock et al 1989, Rachlin 1990). This allows modifying the timing of transmitted and received ultrasound in a way that diffraction-limited focusing can be obtained behind the phase screen. However, perfect guide stars are rare in a clinically realistic situation, where the image region is dominated by diffuse scattering and specular reflections. In the first case, the resulting incoherent echoes prohibit the determination of the aberration profile at the aperture; in the second case, the arrival time of echoes is mainly determined by the geometry of the reflector and does not accurately represent the aberration delays. In both cases, a narrowly focused transmit beam allows the generation of a virtual guide star by restricting potential sources of echoes to a small volume in the range of the resolution cell in the focus of the beam (Flax and O'Donnel 1988). On the downside, degradation of the transmit focus owing to aberrations again reduces the robustness of the technique (Walker and Trahey 1997) and requires iterative approaches to eliminate residual ambiguities O'Donnel 1988, O'Donnell and. To solve this problem, techniques have been proposed that allow determination of the inter-element timing errors based on pulse-echo signal redundancies when acquiring pulseecho signals with a changing Tx aperture (Rachlin 1990, Li et al 1997. Aperture data based AC techniques are well suited in situations where the aberrations can be modelled as a phase screen directly in front of, or close to, the US probe. This model is appropriate to describe, for example, the aberrating effect of the abdominal wall. The irregular thickness and composition (adipose and muscle tissue) imposes a phase delay per transducer element that can be regarded as roughly independent of the position of the Tx/Rx focus in a region that is located far away from the probe. If the speed-of-sound inside the liver is homogeneous, i.e. no additional aberrations are introduced by the liver tissue itself, estimation of the aberrating effect of the abdominal wall based on the arrival time of a statistical ensemble of echoes allows a successful global correction of both the Tx and the Rx focus. This has been shown to lead to reduced clutter in the interior of blood vessels, and better delineation of vessel boundaries, renal capsule and spleen (Rigby et al 2000). However, the assumption of a homogeneous speed-of-sound may become invalid with e.g. an irregular fatty infiltration. In the presence of a spatially distributed volume aberrator, the Rx aberration delay per probe element depends on the location of echoes and can thus not be simplified by a phase screen. The spatial distribution of speed-of-sound of a volume aberrator must then be taken into account for AC in order to achieve spatially invariant diffraction-limited resolution. Various techniques to correct for volume aberrators have been proposed (Liu and Waag 1994, Mallart and Fink 1994, Krishnan et al 1996, Ng et al 1997, among them a signal-redundancy based approach that is suited for AC of diffuse echoes (Ng et al 1997). These approaches are still focused on analysing the relative arrival phase variations across the Rx aperture, with the goal to improve spatial resolution. To our knowledge, no attempt has been made so far to also correct for the residual beam-steering errors that exist in lateral as well as in axial direction owing to an unknown absolute value of the speed-of-sound. However, such a correction would be highly desirable because beam steering errors lead to geometrical distortions of the resulting image and thus limit the accuracy of US guided interventions.
In this study, we present a novel technique which facilitates AC in situations with spatially distributed aberrators, and which resolves several limitations of previous techniques. Our technique relies on directly detecting the effect of acoustic aberration in the reconstructed pulse-echo data, by measuring the changing phase of reconstructed echoes as a function of a changing Tx steering angle at the location where these echoes are generated. The detected spatial distribution of the aberration phase allows a quantitative reconstruction of the distribution of the local speed-of-sound (Jaeger et al 2015), facilitating full AC in presence of a spatially distributed aberrator. As a result, diffraction-limited resolution is obtained independent of location, and beam-steering errors are eliminated leading to an accurate geometric display of echo position. Because echoes from different locations are spatially decoupled in the reconstructed pulse-echo data, the ambiguities linked to specular reflections are avoided, and our technique works equally well with coherent and diffuse echoes. This paper is dedicated to a proof-of-principle phantom study of AC in pulse-echo US using a linear array probe. In addition we provide first volunteer results where we demonstrate strongly improved resolution after AC inside the forearm of a volunteer.

Model of acoustic aberration
In the following we assume that a linear array probe is used for US imaging (see figure 1(a)). The probe transmits a plane ultrasound transient at an angle φ (Tx angle) relative to the axial direction z, via the transmit timing of the individual array elements. The transmitted ultrasound is scattered inside the tissue, and the echoes are detected by the same probe. An inhomogeneous speed-of-sound results in wavefront distortions of both the Tx and the echo wavefront (only Tx wavefront is shown in figure 1(a)). For this study, we assume straight ray propagation and define the aberration delay τ(r, α) as the actual propagation time of ultrasound along a line with angle α that connects r to the probe surface, minus the anticipated time assuming a homogeneous speed-of-sound c 0 . Thereby α is either identical to φ when analysing aberrations along the Tx angle (shown in figure), or to the receive angle ψ (see later on) when analysing receive aberrations. The value of τ(r, α) is given by the line integral of the deviation Δσ(r) of the true slowness σ(r) = 1/c(r) from the assumed value σ 0 = 1/c 0 . For moderate speed-of-sound variations, effects of ultrasound refraction can be described within the concept of straight ray propagation (see section 5). Even though it seems contradictory, this simplification was therefore not a strong limitation to the validity of this study.
Conventionally, acoustic aberrations are investigated in the aperture data and accommodated as part of image reconstruction by adapting the per-element delays in the delay-and-sum beamforming. For the sake of simplicity, we prefer to adopt a convenient model of pulse-echo image generation, which describes the effect of aberration on local echoes directly in the reconstructed image. This model will, as we see later on, also allow efficient a posteriori correction of the already reconstructed image. In a first step we define the reconstructed RF data u(r, φ) in the non-aberrated case (when the true speed-of-sound is homogeneous and equal to c 0 ) as a convolution of an echogenicity function ε(r) with a translation-invariant point-spread function (PSF) h(r, φ) that depends on the Tx angle φ: Thereby the translation invariance is not a requirement but simplifies further notation. In preparation for analysing the influence of aberrations, we decompose u(r, φ) into directional components U(r, φ, ϕ), where ϕ is the centroid angle of the spatial frequency spectrum of U (see figure 1(b)). The decomposition is achieved by convolution of u with a finite set of directional filtering kernels K(r, ϕ) (where ϕ ∈ [ϕ 1 … ϕ Nϕ ]) that are assumed to be orthogonal and that cover the whole part of the spatial frequency spectrum that can possibly be sensed by the probe: Equation (2) suggests that U(r, φ, ϕ) is a convolution of ε(r) with the component H(r, ϕ) of the directional decomposition of the initial PSF h(r, φ). Thereby w(ψ) is the corresponding weight of H in the decomposition of h. The rationale behind this model is following: U(r, φ, ϕ) is the part of u that is sensed by the combination of the Tx angle φ and the Rx angle ψ, where ψ is determined according to the mirror law (see figure 1(b)). The product εH can be envisioned as the back-projected echo of a point scatterer situated in point r′, beam-formed from a collection of transducer elements located at a view angle ψ. The weights w(ψ) contain the angular sensitivity of the probe, as well as angular apodisation weights that are potentially used for image reconstruction. Because clinical ultrasound probes typically have a 'narrow' bandwidth, H is modelled as a product of a slowly varying complex envelope A(r, ϕ) and a complex exponential carrier: Thereby k is a wave vector at the transducer centre frequency f 0 , pointing into either direction φ or ψ.
In the next step aberrations are incorporated into our model, continuing with the concept of H being the beam-formed echo that is detected with the combination of the Tx angle φ and the Rx angle ψ. The effect of the inhomogeneous speed-of-sound is a shift of this echo along angle ϕ, perpendicular to the lines of constant phase of H. The magnitude of this shift is given by the sum of the aberration delays τ(r, α) along angles α = φ and α = ψ. Assuming that the shift is significantly smaller than the spatial scale of variations of the envelope A, the shift can be expressed as a phase factor, and the aberrated directional component U′ becomes: The full aberrated image u′ is then the superposition of the directional components of all angles ϕ: We would like to draw attention to the fact that the aberrated PSF h′ defined in equation (5) can no longer be assumed translation invariant: Because the aberration delays of different Rx angles ψ do not necessarily depend on the spatial location in the same way, their combined influence on h′ is spatially dependent. Assuming that the spatial variations of τ(r′, φ) are slow compared to the size of the PSF h′, r′ can be replaced by r and the corresponding exponential can be factored out of the integral: The assumption on the extent of spatial variations of τ is reasonable as long as the size of the beamformed PSF h′ is close to the diffraction-limited spatial resolution of the system. To what extent this assumption may be valid will be discussed in section 5. The didactic conclusion from equations (5) and (6) is following: Rx aberrations alone determine the shape of the PSF h′ and thus the spatial resolution. As a part of their influence on the shape of h′, Rx aberrations may lead to a lateral shift of echoes owing to refraction. Tx aberrations, on the other hand, have no influence on the shape of the PSF, even though h′ depends on φ. This is illustrated by rewriting the definition of h′ as a summation over Rx angles ψ: Equation (7) illustrates that the relative weight and shift of directional components H (r, ϕ), at angles ϕ relative to a centre angle ϕ 0 = ½φ, are not influenced by φ and thus their superposition will have the same shape independent of φ. The only influence of φ on h′ is via the selection of a changing centre angle ϕ 0 , which leads to a rotation of the PSF with changing φ. Apart from this rotation of the PSF, the only influence of φ on the image is via the phase factor in front of the integral in equation (6). This generates a local shift of echoes perpendicular to the lines of constant phase and thus mainly in axial direction. In contrast to Rx aberrations, Tx aberrations have therefore an insignificant influence on the lateral position of echoes. Any strong change in shape and lateral position is therefore caused by aberration and refraction in the Rx beamforming.

Aberration delay determination
In addition to giving intuitive insight into how aberrations influence the reconstructed image, equation (6) demonstrates the possibility of sensing Tx aberrations in the local echo phase of the image: When scanning the Tx angle φ, the phase factor in front of the integral probes the aberration delays τ for positions r as a function of φ. Unfortunately, τ cannot be determined absolutely because the value of the integral in equation (6) is not known a priori. However, the relative change Δτ of τ when changing the Tx angle by an increment Δφ can be detected. Provided that the Rx aperture is left unchanged, the influence of Rx aberrations is independent of the Tx angle change, and Δτ can be detected as the phase angle of the correlation of the images before and after the angle change: ] 2 i 0 (8) Thereby (*) indicates complex conjugation. The accuracy of equation (8) is limited by the rate of decorrelation owing to the rotation of the PSF with Δφ, as well as by phase aliasing when Δφ results in a too large change Δτ. For this reason, Δφ needs to be small enough to allow a robust determination of Δτ. Within this limitation, accumulative tracking of the echo phase over a sequence of angles φ ∈ [φ 1 , φ 2 , … φ Nφ ] allows determining the aberration delays Δτ for any angle φ n relative to a reference angle φ ref :

Aberration correction
Echo tracking provides relative aberration delays. For full aberration correction (AC), however, the absolute value of the aberration delays is required. For this purpose, the spatial distribution of the absolute speed-of-sound c(r) is reconstructed in a first step, from the relative aberration delay maps Δτ(r, φ ref , φ n ), employing our recently developed CUTE (computed ultrasound tomography in echo-mode) (Jaeger et al 2015). In a second step, the absolute aberration delays τ(r, φ n ) can be simulated for any angle φ using any of the various established numerical techniques. Once the absolute aberration delays are known for any propagation direction (where 'any' refers to just as many as required, which will be specified later), this information can be used for AC. A common way of implementing AC is via the adaptation of the per-element delays in the delay-and-sum beamforming. For this study, however, we continue along our model of image generation to develop a technique of AC that can be applied a posteriori to the already reconstructed image. We start over with equation (6), which we rewrite by expanding the aberrated PSF h′ into its directional components: Assuming that the spatial variations of τ(r′, ψ) are slow compared to the spatial extent of H, the exponential of τ(r′, ψ) can be factored out of the integral: Analogue to the decomposition of u (equation (2)), the aberrated image u′ is decomposed into directional components U′ using the same filtering kernels K. Using the assumption of a slowly varying τ(r′, ψ) from above, U′ becomes: For the last step in equation (12) we used the orthogonality of U and K for ϕ′ ≠ ϕ. Equation (12) states that by directional filtering of u′ using the same filtering kernels K as for u, one obtains components U′ (.,.,ϕ′) of u′ that are identical to the components U(.,.,ϕ′) of u, apart from phase factors containing the aberration delays along both the Tx and the Rx angles. Consequently, AC can be performed on the aberrated directional components U′, to obtain the non-aberrated directional components U: The non-aberrated full image u can be obtained by superposition of the non-aberrated directional components U. At this point we remind the reader that the representation of the echo shift as a mere phase shift was based on assuming that the aberration delays are small compared to the spatial scale of variations of the envelope. However, the goal of this study is the demonstration that full AC leads to the accurate geometric display of echo location, in situations where aberrations are sufficiently strong to introduce a significant shift of the envelope. Therefore AC must include a correction of the envelope in addition to a correction of the echo phase. For this study, correction of the local shift of the envelope was implemented via interpolation of U′ such that the estimated non-aberrated directional components Û were: where In summary, the proposed algorithm for a posteriori AC is: • Directionally filter the aberrated frame u′(r, φ) for angles ϕ n to obtain aberrated directional components U′(r, φ, ϕ n ). The directional filtering can be performed either in the spatial frequency domain, or by applying filtering kernels in space domain. • Calculate for each direction ϕ n the corrected Û(r, φ, ϕ n ) according to equation (14).
Thereby, the ϕ n constitute a finite set of angles as described earlier. The number N ϕ of required angles is on one hand determined by the angular aperture of the probe, on the other hand by the resolution of angles that can be independently sampled by the probe. Because the number of independently sampled angles cannot be larger than the number of transducer elements N e (e.g. 128), a simple way of defining the angle resolution is by dividing the magnitude of the full angle range by N e , resulting in an angle set with N ϕ = N e . In the experimental part of this study we will use even a smaller number, i.e. N ϕ = 121.

Materials and methods
A phantom was built from agar (A7921-100G, Sigma-Aldrich, Germany) with cellulose (Sigmacell Type 20, Sigma-Aldrich, Germany) for ultrasound scattering, and ethanol for speed-of-sound tuning. An aqueous solution with 1% agar and 4% cellulose (by mass) was prepared and dissolved at 80 °C. In addition, 10% ethanol (by volume) was added after the solution had cooled down to a temperature of 60 °C. This solution ('solution 1') was cast into a pre-cooled cubic mould to obtain the phantom bulk with dimensions 8 cm (x-axis) by 5 cm (z-axis) by 4 cm (y-axis). The mould contained prefixed sewing threads that served as distinct acoustic scatterers, as well as a solid cylinder that could be removed after solidification of the bulk material, to provide a cylindrical hole with 10 mm diameter. In parallel, a part of the same solution was cast into a cylindrical mould to provide a matching cylindrical inclusion having identical speed-of-sound as the bulk. A second solution ('solution 2') was prepared, without ethanol to provide a lower speed-of-sound, and only 2% cellulose. This solution was cast into a second cylindrical mould to provide a second inclusion having a lower speed-of-sound and a slightly lower echogenicity than the bulk. The lower echogenicity was needed to identify the inclusion location on B-mode US. Figure 2(a) shows a sketch of the phantom. The linear US probe was placed onto the phantom facet on the side where the cylindrical inclusion was located, with the imaging plane aligned perpendicular to both the cylindrical inclusion and the threads. For pulse-echo data acquisition, we used a Verasonics ® V1-64 ultrasound system (Verasonics Inc., WA, USA). This research system allows ultrasound transmission on 128 channels (two-board version) and signal acquisition on 64 channels simultaneously, as well as data transfer to a host computer via a PCI Express link. We implemented a dedicated script for the acquisition of pulse-echo data from a HDI L7-4 broadband linear vascular probe (ATL Philips, WA, USA). This probe features 128 elements at 0.29 mm pitch, and a bandwidth from 4 to 7 MHz with 5 MHz centre frequency. The script steered the acquisition of plane-wave pulse-echo data using a total of N φ = 121 different Tx steering angles φ ranging from −30° to 30° in steps of 0.5°. From this data, US images were reconstructed off-line using a Fourier-domain reconstruction algorithm. For each Tx angle, the final US image was obtained by coherent averaging of plane-wave images from a collection of adjacent Tx angles covering an angular aperture of 5°. This synthetic aperture technique resulted in images with significantly reduced clutter, compared to planewave images, and thus more robust echo phase tracking.
In a first experiment, the cylindrical hole of the phantom bulk was filled with the cylindrical inclusion made of solution 1, resulting in a phantom with homogeneous speed-of-sound. The resulting images served as a reference for quantifying the maximum achievable spatial resolution and for determining the accurate position of echoes in a situation without aberrations. For a second experiment, the cylindrical inclusion was then replaced by the one providing a negative speed-of-sound contrast relative to the bulk. The setup was designed in a way that the inclusions could be exchanged without moving the phantom, thus allowing a direct comparison of the images taken without and with aberration in terms of spatial resolution and location of echoes.
Aberration correction (AC) was performed according to the methodology described in section 2. For tracking the local echo phase according to equation (8) we employed base-band correlation (O'Donnell et al 1994). The point-wise Hermitian product of the complex (after Hilbert transform) reconstructed RF data of successive frames was calculated, and then lowpass filtered by convolution with a tracking kernel of size Δx (laterally) and Δz (axially). The local phase shift was then estimated as the argument of the resulting complex map. With φ = −30° as a starting angle, the echo phase shift in-between successive frames was gradually accumulated for Tx angles φ up to +30° with a step size of 2.5°, resulting in relative aberration delay maps Δτ(r, φ 0 = −30°, φ). From this data, the relative delay maps were recalculated for a reference angle of φ 0 = 0 as Δτ′(r, 0 , φ) = Δτ(r, 0 , φ) − Δτ(r,−30°, φ) according to equation (9). The 2D distribution of speed-of-sound was then reconstructed from the relative delay maps that were obtained at TX angles φ n ∈ [−30°,−20°,−10°, 10°, 20°, 30°] relative to a reference φ ref = 0°, using the reconstruction algorithm described elsewhere (Jaeger et al 2015). Then the absolute aberration delay maps τ(r, φ) were calculated assuming straight ray propagation, for a total number of 120 propagation angles ranging from −30° to 30° in 0.5° steps and used for AC according to equation (14). Directional filtering of the complex reconstructed RF data was performed by decomposition of the spatial Fourier transform into spectral sectors with 0.25° aperture, corresponding to Rx angle sectors with 0.5° aperture in agreement with the resolution of the simulated propagation angles.
The sewing threads provided distinct point-like targets for demonstrating the effect of aberration on the spatial resolution and position of coherent echoes. On the other hand, it was one of our experimental findings (see section 4) that the apparent spatial resolution of the B-mode speckle pattern in regions of random echoes did not significantly depend on aberration. Blurring of random echoes could only be observed after spatial compounding, i.e. the incoherent averaging of the envelope of images obtained under different Tx angles. For spatial compounding we used a Tx angle range of −20° to 20° in 2.5° steps. The compound B-mode US image of the phantom with the aberrating inclusion in figure 2(b) clearly shows the blurring of the coherently scattering threads as well as of the random background speckle pattern, owing to the aberrations below the cylindrical inclusion. The inclusion is visible at around 15 mm depth and 20 mm lateral position, as a circular region having slightly lower echogenicity than the bulk medium.

Results
All B-mode images are shown in a logarithmic scale covering a range of 60 dB after normalisation to the respective maximum intensity level. Figures 3(a) and (b) are the images of the phantom with the non-aberrating inclusion (same speed-of-sound as background), without (TX angle zero) and with spatial compounding, respectively. The threads are seen as bright point-like structures in the lower part of the image. Even though spatial compounding is often described as a technique for speckle reduction, comparison of the two images demonstrates a nearly identical speckle pattern. The speckle-reducing effect of spatial compounding is based on speckle decorrelation owing to the rotating PSF (see section 2). In the absence of acoustic aberrations, speckle decorrelation is slow because echoes are reconstructed with diffractionlimited resolution and thus minimum spatial coherence length. Figures 3(a) and (b) constitute the 'ground truth' showing the optimum achievable resolution and the accurate location of echoes when the speed-of-sound is homogeneous and accurately known. Figures 3(c) and (d) are the corresponding B-mode images after inserting the inclusion having a lower speed-of-sound than the background. As a result of aberrations, significant blurring of the threads below the inclusion is observed. The dependence of the blurring on the lateral position of echoes is the result of the 2D distribution of the aberrator. In addition to defocusing (broadening of the main receive lobe), the inhomogeneous speed-of-sound resulted in an increased sidelobe level. In its extreme this lead to a double dot appearance of threads 3 to 5. On the other hand, no blurring is observed of the speckle pattern in the region of random echoes below the inclusion. The explanation is that the reconstructed RF amplitude of random echoes is always a superposition of harmonic components with random phase and random amplitude, independent of aberration. The lateral correlation length of diffuse echoes in the B-mode image depends only on the angular aperture of the Rx focus. Whereas the speckle statistics was not visibly influenced, a close look reveals that the speckle pattern itself was significantly changed by aberrations compared to figure 3(a). Note also the two dark streaks extending downwards from the left and right edges of the inclusion. These shadows result from the phase cancellation of the Tx wave-front when propagating along boundaries in-between regions having different speed-of-sound. Figure 3(d) is the compounded image as already shown in figure 2(b). Note that phase cancellation shadows point into different directions in the different single Tx angle images and thus disappear after averaging in the spatial compound image. With spatial compounding, significant blurring is observed in the speckle pattern below the inclusion. It is an interesting observation that the speckle pattern itself was not significantly changed by a changing Tx angle, but the blurring resulted from a lateral shift of this pattern. This shift could be observed even though, according to equation (7), the lateral position of reconstructed echoes is mainly determined by the Rx beamforming independent from the Tx angle. In agreement with equation (7), no lateral shift was observed for the point targets. The reason for the apparent lateral shift in regions of random echoes will be discussed in the next section. Figure 4(a) is the CUTE image that shows the reconstructed spatial distribution of the speed-of-sound. The inclusion is visible as a well-defined region of darker appearance visualising a lower speed-of-sound than the background. The bright area below the inclusion suggesting a higher speed-of-sound has been shown by simulations to be an artefact caused by acoustic refraction. In addition, strong spatial variations of the reconstructed speed-of-sound at 35 mm depth are also artefacts that were caused by the echoes from the point targets. A lateral profile through the CUTE image at the depth of the inclusion ( figure 4(b)) illustrates the inclusion contrast of around −60 m s −1 corresponding to −4% of a background speedof-sound of 1520 m s −1 . Reduction of artefacts is subject to ongoing research, and they were circumvented in this study by considering only the top 20 mm of the CUTE image containing the inclusion for calculation of the absolute aberration delays.
The absolute aberration delays were subsequently used for aberration correction (AC) according to equation (14). Figure 4(c) is the result of AC when applied to the image obtained with a single Tx angle. After AC, all threads appear with uniform spatial resolution independent of their lateral position. In terms of resolution, the image looks similar to when the speed-of-sound was homogeneous (compare to figure 3(a)) suggesting that diffraction-limited spatial resolution was obtained. A close look also reveals that, apart from the phase cancellation shadows, the speckle pattern below the inclusion looks identical to the one that was found in the homogeneous case (compare to figure 3(a)). Figure 4(d) is the result of compounding after correcting the images from the different Tx angles. The compounding image shows a uniform speckle pattern with no blurring suggesting that the lateral shift of speckles with changing Tx angle was eliminated.
To precisely evaluate the resolution and geometrical display of echo position after AC, lateral and axial profiles of the point targets 3 to 5 are given in figures 5 and 6. The profiles show the squared envelope (echo intensity) in linear scale, normalised to the respective maximum intensity. Figure 5 shows the lateral profiles of the point targets, in the image of the homogeneous phantom (reference, grey line), and in the images of the inhomogeneous phantom without AC (solid black line) and with AC (dotted line). Whereas the lateral profiles of the aberrated echoes confirm the pronounced double-dot appearance and differ strongly from the reference profiles, the profiles after AC are nearly identical to the reference for all three point-targets. This comparison confirms that diffraction-limited resolution as well as accurate lateral position of echoes was obtained with AC. In comparison to the reference, the corrected profiles show a slightly increased sidelobe level indicating residual beam-forming errors. A potential explanation of such residual beamforming errors will be given in the next section. Figure 6 allows a comparison of the axial profiles. According to this figure, aberration led to an axial shift of the echoes relative to their true position. This shift is explained by the unanticipated propagation delay caused by the negative inclusion contrast, and it was fully corrected by AC as the corrected profiles look identical to the reference.

Discussion and conclusion
Our results demonstrate that aberrations caused by a spatially distributed aberrator can be fully compensated, based on measuring spatially resolved the TX angle-dependent aberration delays in the reconstructed image. Unlike a flat aberrator that acts like a phase screen in front of the probe, the cylindrical aberrator in this study was located at a distance away from the probe and therefore the amount to which it interfered with the Rx beamforming varied with the relative position of ultrasound reflectors inside the phantom. Nonetheless, the technique proposed in this paper was fully capable of correcting for the introduced aberrations, such that diffraction-limited resolution was obtained of the thread targets as well as diffuse echoes independent of the lateral position. In addition, full compensation for beam-steering errors and thus an accurate geometric display of echo position was achieved. To keep the mathematical fundament of the technique simple, a number of assumptions was introduced that will be discussed in the following: (a) The aberration delay was assumed to vary on a spatial scale larger than the size of the aberrated PSF. As a consequence it was possible to decouple the effects of Tx and Rx aberrations on the image formation in equation (6). According to this equation, the relative aberration delay when changing the Tx angle results in an axial shift of echoes, which qualified phase tracking as a method of quantifying the relative aberration delay in equation (8). The experimental results, however, indicate that the underlying assumption was not fully satisfied. The apparent lateral shift of diffuse echoes with changing Tx angle (see description of figure 3(d)) that was not anticipated by equations (6) and (7) can be  explained with a lateral extent of the aberrated PSF that was larger than the lateral variations of the Tx delay profile. This made it possible that interfering echoes from randomly distributed scatterers were axially displaced by a different amount when changing the Tx angle, leading to an apparent lateral shift of the interference pattern (see figure 7 for illustration). The perfect outcome of the AC results indicates that assumption (i) was not a stringent requirement for the applicability of the technique. (b) The applicability of the tracking approach was based on assuming that the tracking angle step Δφ is sufficiently small to avoid decorrelation of the PSF and phase aliasing. Decorrelation results from the changing interference of adjacent echoes with the rotation of the PSF. Correspondingly, the decorrelation rate is expected to increase with the lateral extent of the aberrated PSF and thus with the grade of aberrations. A detailed investigation of the influence of aberrations on the accuracy of echo phase tracking is beyond the scope of the present study. Nevertheless we would like to mention that an angle step as large as 5° allowed robust phase tracking in the phantom study. Since a step of 2.5° was finally chosen, PSF decorrelation is unlikely to have been a significant source of tracking errors. Phase aliasing can occur when the relative aberration delay change Δτ is larger than half the oscillation period of the probe's centre frequency. In the phantom study, the speed-of-sound contrast and the diameter of the inclusion were 4% and 10 mm, respectively, leading to an absolute aberration delay of 0.27 µs for ultrasound propagating through the middle of the inclusion. The relative delay change when changing the angle by only 2.5° was much smaller and consequently far below the 0.2 µs period of the 5 MHz probe. Phase aliasing could, however, become a limitation in presence of a higher speedof-sound contrast. If increasing the number of angles is not an option due to limitations to the acquisition time, a more robust tracking approach will then be required. Possible choices are the combined autocorrelation technique (Shiina et al 2002), a multi-frequency phase tracking technique , or a multi-step approach where the speedof-sound in a superficial tissue layer is used for pre-compensation of the echo phase in the underlying medium prior to phase tracking in the underlying medium. (c) The applicability of the proposed compensation for aberrations a posteriori in the already reconstructed image was based on assuming that the spatial variations of τ are slow compared to the spatial extent of the PSF H of the directional components (see equations (10) and (11)). This condition was clearly not fulfilled in the phantom study: The width of H was determined by the 0.25° aperture of the directional filtering sectors, corresponding to 0.5° aperture of Rx angle sectors. Employing Abbé's equation for an ultrasound ray with Θ = 0.5° aperture and an acoustic wavelength λ = 0.3 mm, the width of H was around λ/Θ = 34 mm. Therefore the width of H was much larger than the scale of spatial variations of τ. This explains the increased sidelobe level of the corrected images in figure 4 compared to the non-aberrated images in figure 3, and it is remarkable how well AC still worked out. The 4% speed-of-sound contrast used for the phantom study is realistic and explains why the proposed compensation technique works also clinically (see below). The speed-of-sound in tissue can however vary on an even larger scale of up to 10%. At such extreme contrast levels we expect that the technique will finally fail: If the relative displacement of interfering directional echoes is larger than half the spatial oscillation period, the interference pattern will be altered. The correction approach according to equation (11) can only correct for the local shift but not for a changing interference pattern of directional echoes. Note that, whereas conditions (i) and (ii) referred to the core of the study, i.e. the demonstration of the possibility of determining aberrations spatially resolved in the reconstructed image, condition (iii) is not fundamental to the key message of the paper. Any technique of incorporating the relative or absolute aberration delays for generation of a corrected image is equally valid, among them the incorporation of the aberration delays in the delay-and-sum (DAS) beamforming. With a large speed-of-sound contrast, we expect that a corrected DAS approach will yield superior performance than the proposed a posteriori correction technique. (d) A further assumption was the straight-ray approximation of ultrasound propagation. We would like to draw attention to the fact that this assumption is not as stringent as it might seem at a first glance: The straight ray approximation is rather a conceptual trick than an actual assumption on how sound propagates. The differential propagation delay along adjacent 'rays' from a point reflector to the probe leads to a relative shift of the corresponding directional components of the echo. This leads to a deviation of the intersection of the directional components from the true position of the reflector, which intrinsically contains the effect of refraction. Above a certain speed-of-sound contrast, however, the actual bending of the 'straight rays' needs to be taken into account, because the Tx angle for which the aberration delay is detected will differ from the angle for which the aberration needs to be compensated. In addition, strong refraction/diffraction will lead to caustics in the wave field, i.e. the intersection of initially parallel or even divergent 'rays', and thus to multiple arrivals of different portions of the Tx wave front at the same reflector from different angles. This makes the phase tracking ambiguous, and the assignment of phase delays to propagation angles erroneous. With the moderate speed-of-sound contrast in the present phantom study, the influence of caustics was apparently insignificant.
To accommodate for refraction and caustics in a situation with higher speed-of-sound contrast, we propose an iterative approach where, in a first step, the speed-of-sound is reconstructed inside a superficial tissue layer where refraction effects are still weak. Based on the result, the images can be pre-corrected such that more robust echo tracking will then be possible inside the region of the underlying tissue. We hypothesize that, with the use of an accurate model of sound refraction, AC might be possible also with significant refraction and caustics.
As a preliminary demonstration that the proposed AC for a spatially distributed aberrator works clinically, figure 8 shows images of the forearm of a volunteer. The conventional B-mode US (figure 8(a)) revealed significant lateral blurring of the point-like echoes of muscle fibre fascicles inside a distinct region denoted by a dashed circle. Left to this region, the apparent lateral resolution was similar to the resolution in the region of the superficial muscle layer. The different spatial resolution of echoes depending on their lateral position indicates the presence of a spatially distributed aberrator that was located well away from the probe. Based on the evidence, we assume that the aberrations were caused by the tissue layer denoted by an arrowhead, in-between the superficial and the deeper muscle layers. With AC ( figure  8(b)), the spatial resolution became independent of the lateral position of reflectors, demonstrating that AC for a spatially distributed aberrator was achieved.
The present study focused on pulse-echo US using a linear (1D) array probe, and thus on AC for a two-dimensionally varying speed-of-sound in the imaging plane of the probe. However, the methodology can be readily adapted to 3D imaging using 1.5D or 2D array probes, where it can improve the elevational in addition to the lateral resolution. It can also be extended to improve performance of other US-based techniques that are used in combination with pulse-echo US, such as optoacoustic (photoacoustic) imaging (improve Rx focus), acoustic radiation force imaging (improve radiation force focus and thus tissue displacement), or even therapeutic applications of high intensity focused ultrasound (improve dose delivery). The technique could even be promising when translated to optical coherence tomography (OCT) which is conceptually close to echo US. Interferometric synthetic aperture microscopy (ISAM) (Ralston et al 2007) is an OCT pendant to synthetic aperture US. In a medium with constant refractive index, this approach gives a uniform optical resolution independent of the position relative to the focal plane of the microscope. An inhomogeneous index of refraction, however, is likely to cause a notable degradation of this resolution. Scanning of the angle of the transmitted optical radiation could be achieved by partially blocking the collimated beam before entering the microscope objective, and the resulting phase shift of the reconstructed image could then be used for aberration correction of this image in the same way as is was done in this study. This could yield improved resolution, contrast and thus larger imaging depth of ISAM images of biological tissue.