Frequency-Dependent Attenuation Reconstruction with an Acoustic Reflector

Refraction is the change in wavefront direction at interfaces with acoustic impedance variation, and diffraction is the bending of waves, e.g., around high impedance barriers [41]. We herein assume that the speed-of-sound and density, and hence acoustic impedance, in soft tissues vary marginally such that refraction and diffraction effects can be neglected. Accordingly, band-limited ultrasound pulses at a center frequency $f$ are assumed to propagate along rays. Note that, due to constructive and destructive interferences, acoustic intensify from a transducer element has a response as a function of wavefront (ray) propagation direction $\theta$, which causes natural focusing of energy in the main beam axis as well as side lobe effects. The amplitude of the reflector for transmitting with channel $t$ and receiving with channel $r$ can then be described by

$$A_{t,r}=A_{t}(\theta,f)\cdot R(\theta)\cdot S_{r}(\theta,f)\cdot\exp\bigg{(}\!\!-\!\!\!\int_{\text{ray}_{t,r}}\!\!\!\!\alpha(x,y,f)\ \mathrm{d}l\bigg{)},$$

(1)

where $R(\theta)$ is the incident-angle dependent reflection coefficient at the reflector interface, and $A_{t}(\theta,f)$ and $S_{r}(\theta,f)$ are, respectively, the band-limited amplitude of Tx element $t$ and the sensitivity of Rx element $r$, at center frequency $f$ and transducer directivity $\theta$. The exponent describes frequency-dependent amplitude decay based on the line integral of attenuation $\alpha$ along ray_t,r from element $t$ to $r$.

To estimate attenuation $\alpha$ from $A_{t,r}$, any confounding factors in (1), such as transducer impulse-response influence in both $A_{t}(\theta,f)$ and $S_{r}(\theta,f)$ as well as reflection characteristics in $R(\theta)$, must first be compensated for. Considering the impulse response, the band-limited signals are normalized with a calibration experiment in water, where the speed-of-sound $c_{\mathrm{water}}$ and the attenuation coefficient $\alpha_{\mathrm{water}}$ are known from the literature, given the water temperature and imaging frequency. By using the same predetermined frequency bandwidths and reflector distances in the water calibration as in the subsequent tissue experiments, $A_{t}(\theta,f)$ and $S_{r}(\theta,f)$ components can be effectively estimated and factored out. However, $R(\theta)$ is still likely to differ between any calibration and application, because of the speed-of-sound mismatch between water and tissue. Such reflection coefficient changes at the acoustic reflector interface can nevertheless be analytically estimated using Snell’s law, as the wavelength ($\approx 300\,\mu\mathrm{m}$) is small compared to the reflector dimensions. One can then write the reflection coefficient at the interface of the reflector and some medium $k$={water,tissue} as follows:

$$R_{k}(\theta)=\frac{m_{k}\cos(\theta)-n_{k}\sqrt{1-\frac{\sin^{2}(\theta)}{n_{k}^{2}}}}{m_{k}\cos(\theta)+n_{k}\sqrt{1-\frac{\sin^{2}(\theta)}{n_{k}^{2}}}},$$

(2)

where speed-of-sound ratio $n_{k}$=$c_{\mathrm{reflector}}/c_{k}$ and density ratio $m_{k}$=$\rho_{\mathrm{reflector}}/\rho_{k}$. We herein assume $\rho_{\mathrm{tissue}}$ $\approx$ $\rho_{\mathrm{water}}$ $\approx$ 1000 kg/m³ and that the speed-of-sound is dispersion-free (i.e., not frequency dependent) within the range of the utilized imaging frequency bandwidth. A normalized reflection amplitude matrix $B^{\prime}$ shown in Fig. 1e can then be computed with its elements corresponding to Tx $t$ and Rx $r$ being

$$b^{\prime}_{t,r}=\frac{A_{t,r,\mathrm{tissue}}R_{\mathrm{water}}(\theta)}{A_{t,r,\mathrm{water}}R_{\mathrm{tissue}}(\theta)}\,.$$

(3)

Substituting (1) in (3), we can then relate calibrated and band-limited measurements at frequency $f$ to attenuation $\alpha(x,y,f)$ at image location ($x$,$y$) by

$\displaystyle b_{t,r}=\ln b^{\prime}_{t,r}=-\int_{\text{ray}_{t,r}}\!\!\!\!\alpha(x,y,f)\ \mathrm{d}l\approx-\!\!\!\sum_{i\in\text{ray}_{t,r}}\!\!l_{i}\alpha_{i}(f)\,,$

(4)

where the ray integral is discretized as a summation over a reconstruction grid. For this discretization, we use local attenuation value $\alpha(x,y)$ weighted by its path length within the grid element ($x$,$y$).

Given $M$ logarithms of normalized reflection amplitudes $\mathbf{b}\in\mathbb{R}^{M}$ at a specified frequency band, we perform a tomographic reconstruction of the spatial UA distribution on a $N_{1}$$\times$$N_{2}$ spatial grid by formulating the following convex optimization problem:

$$\boldsymbol{\hat{\alpha}}(f)=\underset{\boldsymbol{\alpha}(f)}{\arg\min}\|\textbf{L}\boldsymbol{\alpha}(f)+\mathbf{b}(f)\|_{1}+\lambda\|\textbf{D}\boldsymbol{\alpha}(f)\|_{1},$$

(5)

where $\textbf{L}\in\mathbb{R}^{M\times N_{1}N_{2}}$ is the sparse ray path matrix (cf. Fig. 2) that implements (4) and $\boldsymbol{\alpha}(f)\in\mathbb{R}^{N_{1}N_{2}}$ is the reconstructed image at the center frequency $f$ of the band-pass filtered signal (cf. Fig. 1b). Due to the ill-conditioning of $\mathbf{L}$, a regularization term controlled by weight $\lambda$ encourages spatial smoothness. Similarly to [30, 29, 33, 36], we herein use $\ell_{1}$-norm for both the data and regularization terms for robustness to outliers in, respectively, the measurements and the reconstructed image (edges) [16]. Due to a lack of full angular coverage of measurements, regularization matrix $\mathbf{D}$ implements LA-CT specific image filtering to suppress streaking artifacts orthogonal to missing projections via anisotropic weighting of directional gradients [33]. In axis-aligned directions a Sobel kernel, and in diagonal directions a Roberts kernel is used [29] (see Fig. 2d/e). Similarly to [35], we herein utilize a $\kappa$$=$$0.9$ anisotropic weighting for a stronger regularization in the lateral direction, taking into account the uneven distribution of measurement paths. In this paper we empirically set $\lambda$=0.6 for all experiments , following [31]. For the numerical solution of the optimization problem (5), a limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [4, 15, 18, 39] from the unconstrained optimization package minFunc¹¹1https://www.cs.ubc.ca/~schmidtm/Software/minFunc.html is used.

The true center frequencies after bandpass filtering are calculated as the spectral weighted average of the filtered signal, and these used for further analysis. Note that the true center frequencies (e.g., $[8.1,6.1,3.9]$ MHz as indicated in the northwest corner of each reconstruction ) are slightly different than their nominal values (respectively $[8.5,6.5,3.5]$ MHz) used as bandpass centers. This is because the emitted wide-band pulse is approximately Gaussian with a peak at 5 Mhz, i.e. more energy is concentrated around this center. Consequently, even after a uniform filter, the resulting signals have weighted averages $f$ that are shifted towards the original peak at 5 MHz.

In Fig. 3b, the frequency-dependent attenuation values are shown for two regions, within the inclusion and outside. These plots already yield a spectral fingerprint of the attenuation characteristics in each tissue. To parametrize such spectroscopic attenuation reconstructions in a standardized way, one can assume the conventional power-law relation:

$$\alpha(f)=\alpha_{0}f^{y}\,,$$

(6)

where attenuation is a function of frequency with medium-specific attenuation coefficient $\alpha_{0}$ and medium-specific exponent $y$. Accordingly, separately at each spatial image location, $\alpha_{0}$ and $y$ are estimated by least-squares fitting the power-law function to the band-limited reconstructions, yielding the parameter maps shown in Fig. 3c and d.

Note that the attenuation term in (5) could have alternatively been replaced by the power-law equation (6), for solving all frequency bandwidths simultaneously. This could provide further robustness by regularizing reconstructions across frequencies by their expected power-law relation as well as by allowing for the direct spatial regularization of the two final parameters. However, this would also increase the optimization problem size in several folds, impeding an efficient numerical solution.

Three simulations (#1–#3) shown in Fig. 4 were conducted (case $\#$1 is the one illustrated earlier in Fig. 3), as representative cases with inclusions separated laterally and axially, as well with both positive and negative contrast to background. Given the ground truth maps in Fig. 4 it can be seen that the attenuation coefficient $\alpha_{0}$ is successfully reconstructed in all cases. It is observed that with increasing UA distribution complexity (i.e., from simulation #1 to #3), the reconstruction accuracy decreases. In particular, axially-separated inclusions are reconstructed poorer compared to those laterally separated. This is due to the transmit-receive paths traversing the tissue in a mostly axial direction and hence lateral projections being missing in the resulting LA-CT, as was also reported and discussed in the context of speed-of-sound reconstruction in [31, 33]. In general, within the inclusions a shift of $\alpha_{0}$ towards background values is observed. This is due to the lack of lateral projections combined with the regularization term in (5) axially elongating the inclusions, hence spreading them over a larger area. Since the cumulative UA effect shall stay the same, thus per-pixel inclusion values are effectively averaged with the background, inversely proportional to such artifactual area increase.

In Fig. 4 the frequency exponent $y$-maps should not present any contrast, since the k-Wave simulated ground-truth data does not allow to spatially vary $y$. Reconstructed $y$-maps are seen to exhibit a deviation from such constant ground-truth value, with a fixed background offset and some artifacts along the axial edges of inclusions. The latter artifact originates from the artificial axial elongation of inclusions due to above-described LA-CT nature. This herein manifests itself stronger with increased frequency band, i.e. the higher the frequency band, the longer an inclusion axially is in the reconstruction. This effect can be observed in Fig. 3a, where the inclusion shape at $8.1$ MHz is axially slightly longer compared to that at $3.9$ MHz. This artifactually higher UA values around the axial edges at higher frequencies are then erroneously attributed by the model to increased $y$ values around these edges. Note that an inverse problem with multiple-bands formulated in the same linear system would allow to constrain and regularize the reconstructions at the parametric level (i.e. edges enforced consistently at all frequencies) and hence may prevent such axial edge artifacts in $y$-maps.

where an improvement is observed in every simulated case for every metric. The average improvement in UA estimation accuracy (RMSE) is 17%, with a substantial enhancement in reconstruction contrast: an average CNR improvement of 176% and an average CRF improvement of 80%. In an actual imaging scenario such improved contrast would enable easier identification and characterization of abnormality, especially in challenging cases where more heterogeneous compositions are imaged (as in simulations $\#2$ and $\#3$). In comparison to a single value of frequency-averaged reconstructions, frequency-dependent UA imaging yields two biomarkers, which can either be used in independently to compensate confounding effects or in combination as a multi-parametric imaging biomarker, both of which can thus be potentially more expressive for tissue characterization.

with their corresponding B-Mode images and attenuation reconstructions for the frequency-dependent ($\alpha_{0}$, $y$) parameters as well as the fullband frequency-averaged ($\alpha$) attenuation (as in [31]). For the frequency-resolved approach, an equivalent $\alpha$ is also computed (for comparison to [31]) based on (6) using the reconstructed coefficient and exponent maps. For the phantoms #1 & #2 with muscle inclusions in gelatin, a clear contrast is visible for all UA reconstructions (Fig.6c-f). With our proposed frequency-dependent method, it can be seen that the observed attenuation contrast originates from a variation not only in attenuation coefficient but also in frequency-dependence. For phantom #3, the reconstructed values with our proposed frequency-resolved method are consistent with phantoms #1 and #2, where the inclusion and background are inverted. Compared to the fullband method, $\alpha_{0}$- and $y$-maps are seen to better preserve inclusion geometry with a reduced axial elongation, corroborating the earlier findings in the simulations. This geometric accuracy is especially important in cases, where the inclusions are not directly observable in the B-Mode images, e.g. in phantom #3 where the scattering gelatin and the surrounding muscle being isoechoic.

To suppress axial edge artifacts, especially in $y$-maps, a possible solution would be to first delineate target regions for UA characterization and then to solve the reconstruction problem with UA constrained as piece-wise constant within each such predefined region. Such regions could be delineated manually or automatically from B-Mode images or from $\alpha_{0}$-maps of initial unconstrained UA reconstructions. It is also possible to reconstruct all frequency bands in a single linear system for consistent regularization and an implicit model-fitting across frequencies (e.g., to better leverage higher SNR frequencies), although this would yield much larger inverse-problems for reconstruction.

In general, the background reconstructions are seen to deteriorate towards the two lower corners, likely due to fewer tomographic projections of these regions, as also illustrated in Fig. 2c. Such errors may also propagate to the rest of the image, e.g. deteriorating how well the inclusions, especially those nearby the lower corners, are reconstructed. Consequently, phantom #4 with two inclusions on the sides yields a poorer reconstruction. Indeed, in contrast to the superior performance of the frequency-resolved method for phantoms #1–#3, for phantom #4 it yields poorer results compared to the fullband, with localization and geometry estimation errors. Surprisingly, the two inclusions with the same material are not even characterized similarly in the reconstructions of neither $\alpha_{0}$ nor $y$. Nevertheless, similar to the fullband representation, an attenuation map $\alpha$ (Fig.6e) can be derived from the frequency-dependent $\alpha_{0}$- and $y$-maps using (6), which is then seen to recover any lost contrast, even for phantom #4. This observation is supported quantitatively by average CNR of attenuation across ex-vivo samples being comparable as $(3.3\pm 2.6)$ for the fullband method and $(3.1\pm 2.2)$ using (6) from frequency-resolved parameter maps. These results also demonstrate that both our frequency-resolved and fullband methods can find applicability in different clinical scenarios.

For the scattering gelatin with cellulose (referred as gelatin+), a low standard-deviation in Table LABEL:tab:exphant indicates a high reproducibility considering all three UA parameters (also can be observed in Fig.7b/b’/b”). Despite an expected natural tissue heterogeneity, all reconstructions of the muscular tissue also follow relatively tight distributions with somewhat similar values as seen in Figs.7a/a’/a”. Muscle tissue appears to have a high $\alpha$ and $\alpha_{0}$ whereas a low exponent $y$, and accordingly presents a stark contrast to gelatin+. For the pure gelatin case, we expect an overall attenuation lower than gleatin+ due the lack of scattering that contributes to attenuation. Indeed, for phantom #2 compared to #1 & #3, both $\alpha_{0}$ and $\alpha$ are significantly lower, as expected, cf. Fig.7b/b”/c/c”. Interestingly, $y$ is similar with and without scattering, cf. Fig.7b’/c’. This can be explained by the fact that frequency-dependent attenuation is mostly due to relaxation and viscous absorption, and not due to scattering. Also the cellulose scatterers of size 50 $\mu$m are smaller than the wavelength of the highest frequency band, thus interacting with all involved frequency bands similarly as a scattering source. Such observations unfortunately cannot be made from phantom #4, the reconstructions of which as mentioned above were affected severely by artifacts.

With the separate characterization of the power-law coefficients of frequency-dependent attenuation, observed attenuations can now be extrapolated to a wider hypothetical frequency range using the power-law relation (6). This is presented in Fig. 8 for the muscle and gelatin+ materials using the average values reported in Tab. LABEL:tab:exphant. From Fig. 8 we can conclude that within typical medical ultrasound bandwidth, i.e. under 20 MHz, contrast between the attenuation $\alpha$ of these two materials would be maximized at $5.6$ MHz as $\Delta\alpha=2.8$ dB/cm. Conversely, at $\approx 18$ MHz as seen in Fig. 8, the two materials will likely not exhibit any contrast at all. This demonstrates that using a frequency-averaged UA imaging, the achieved contrast may highly depend on the imaging frequency. Similarly, an imaging frequency unfit for a practical scenario may potentially lead to suboptimal tissue characterization and even potential misclassification, e.g., of a tumor. Our proposed frequency resolved attenuation imaging method redeems such limitation of having to select an optimal imaging frequency, which is likely unknown a-priori, by allowing for a characterization independent of utilized frequency range (and thus imaging transducer and frequency setting).

Herein we only take the speed-of-sound difference between water calibration and target tissue into account for compensating for the reflection coefficient differences, cf. (3). Tissue speed-of-sound difference between the water calibration and target tissue may, however, also affect the angular beam profile of the transducer elements, hence confounding the assumption of amplitude matrix constancy between water and tissue in (3). Nevertheless, we believe such discrepancy of the beam profiles to be minimal, given the successful reconstruction results demonstrated with both our frequency-resolved and fullband attenuation imaging methods.

In our experiments we utilized a transmit pulse of five cycles. This could be shortened for the full spectral bandwidth of the transmit pulse to be extended up to the full transducer bandwidth, hence allowing for a wider range of reconstruction frequencies. For our purposes, within a given transducer bandwidth, either the echos from several narrowband Tx’s can be used to separately analyze the data or a single wideband Tx can be analyzed in subbands with post-processing as presented herein. Our method is thus applicable using existing US transducers and systems, provided access to raw RF data.

Similarly to most earlier works in US-CT, we herein assume a straight ray wavefront propagation, which is however only true for a constant speed-of-sound distribution. Local speed-of-sound inhomogeneities may indeed refract and aberrate wavefronts, resulting in inaccurate amplitude readings as well as the assumed wavepaths (e.g. Fig.2b) being different than actual wavepaths – both possibly leading to inaccuracies in reconstruction. A potential remedy would be to correct any ray refractions based on an a-priori speed-of-sound reconstruction, e.g. using time-of-flight measurements as in [36].

With the proposed frequency dependent UA imaging method, the opposite side of an imaged tissue needs to be accessible with a passive reflector, which limits potential anatomical regions to the breast, the extremities, male genitelia, etc. These are similar to the locations applicable for passive reflector based speed-of-sound imaging as well as for full US-CT methods. It may further be possible to use a passive reflector during open or laparoscopic surgeries, such as using a reflector with a so-called “pick-up” US transducer [38].

For a multi-parametric tissue characterization, our proposed method add s two additional imaging biomarkers complementary to robust speed-of-sound imaging methods [29, 30], which is potentially superior to elastography [17]. Hence, compared to existing ultrasound imaging modalities, a more comprehensive assessment of tissue composition and potential pathology is enabled using the proposed frequency-dependent attenuation imaging.

In addition to using the proposed attenuation imaging for diagnostic purposes, the derived UA distribution can also be utilized to improve other ultrasound imaging modalities. For instance, similarly to [30], where speed-of-sound imaging is used to correct for aberration effects in beamforming, UA maps can be used to compensate signal gain individually for each Rx-element and beamformed pixel combination. This can homogenize B-mode image appearance by removing artificial contrasts, can enhance optoacoustic [25] and functional imaging [28], and can further improve displacement estimations, e.g. in shear wave elastography.

A practical limitation of the herein proposed frequency-resolved attenuation imaging is the computational time required for the reconstructions due to the iterative numerical optimization employed. A deep variational network solution such as in [44, 45] with reconstruction times in the range of milliseconds could allow to overcome this restriction towards real-time imaging.