Acoustic radiation force single track location shearwave spectroscopy for viscoelasticity imaging, methods and systems thereof

Abstract

A method of mapping viscoelasticity using Acoustic Radiation Force (ARF) source Single Track Location (STL) shear wave spectroscopic imaging is described.

Description

DOCKET NO: 1134-247 PCTTITLEACOUSTIC RADIATION FORCE SINGLE TRACK LOCATION SHEARWAVE SPECTROSCOPY FOR VISCOELASTICITY IMAGING, METHODS AND SYSTEMS THEREOF

[0001] This application claims priority from U.S. Provisional Application No. 63 / 696,117, filed September 18, 2024, which is incorporated herein by reference.FIELD

[0002] This application relates to the field of Ultrasound Shearwave Elastography (USWE) for Acoustic Radiation Force based Viscoelasticity Imaging. In particular, technology that generates estimates of Shear Storage Modulus, [ is(a>)\, Shear Loss Modulus, \^L(G)\, Complex Shear Modulus and Loss Tangent derived from sfco ] and , [« / / «>] along with Shear Elastic Modulus [uifyjj] for non destructive evaluation of tissue’s viscoelastic biomechanics using ultrasound.BACKGROUND

[0003] Ultrasound Shearwave Elastography (USWE) provides contrast based on tissue stiffness complementary to conventional US imaging. Although most commercially available USWE systems assume that tissue is elastic, elasticity alone is not a reliable discriminant for many pathologies.

[0004] In spectroscopic estimation, the shear wave signal is averaged over a spatial extent to produce a single estimate of viscoelasticity. However, local reconstruction of viscoelasticity helps image spatial heterogeneity of biomechanical variations in tissues. Some studies have explored the local reconstruction of the shear wave dispersion but did not account for the shear wave attenuation.

[0005] Wave attenuation estimation is challenging due to the loss of amplitude to nearfield diffraction or geometric spreading from the source, which has a finite width and elevational extent in the out-of-scan direction. These effects must be compensated using analytic corrections to the propagation model. The correction factor must be chosen carefully to avoid estimation bias and depends on the acoustic radiation force (ARF) source geometry.

[0006] Frequency shift approaches fit the shearwave amplitude spectrum to a distribution function and estimate the shearwave attenuation based on the smoothing function parameters. However, they are sensitive to the choice of the model and its parameters.

[0007] Cross-spectral techniques estimate the amplitude loss and spectral phase shift from an empirical transfer function (TF) and derive amplitude loss relative to a referencespectrum. In an STL ensemble, the shearwave due to the first ARF push serves as the reference waveform. The differential form of the measurement allows the diffraction and geometric losses to be compensated in the normalized phase and amplitude spectra without additional model parameters or geometry assumptions. The STL acquisition scheme aligns well with the frequency domain TF approach, as the cross-spectra of two identical ARF sources normalized by the auto-spectrum of the reference shearwave capture solely the material-induced propagation loss and can potentially minimize speckle-bias from push locations.

[0008] Conventional rheometers are not designed to assess the viscoelasticity of tissue in-situ. They are based on interrogating bulk properties of tissues, limited only up to several hertz in frequencies.

[0009] The present application provides a novel non-invasive 2D shear wave viscoelasticity imaging (VI) algorithm to reconstruct the viscoelastic properties of tissues and tissue-mimicking materials.SUMMARY

[0010] An aspect of the application is directed to a method of mapping viscoelasticity using Acoustic Radiation Force (ARF) source Single Track Location (STL) shear wave spectroscopic imaging comprising: generating shear waves in a region of interest in a subject from a plurality of spatially offset locations; tracking the shear waves at a common receiver; measuring frequency-dependent shear wave speed; measuring attenuation of the shear waves; estimating by a cross-power spectral density (CPSD) transfer function of phase shift and amplitude decay as a function of propagation distance of the tracked wave; and constructing multifrequency parametric two-dimensional maps capturing local tissue viscoelasticity.

[0011] Another aspect of the present application is directed to a method of estimating quality of viscoelasticity estimates, comprising the step of estimating coherence between tracked shearwaves within an imaging ensemble.

[0012] These and other aspects, objects, features, and advantages of the example embodiments will become apparent to those having ordinary skill in the art upon consideration of the following detailed description of example embodiments.BRIEF DESCRIPTION OF THE DRAWINGS

[0013] An understanding of the features and advantages of the present application will be obtained by reference to the following detailed description that sets forth illustrative embodiments, in which the principles of the application may be utilized, and theaccompanying drawings. The figures herein are for illustrative purposes only and are not necessarily drawn to scale.

[0014] FIG. 1 shows ARFI STL shearwave spectroscopic viscoelasticity imaging reconstruction pipeline. The experimental setup shows a Philips ATL linear transducer mounted on a motion-controlled positioner in scan configuration. The schematic illustrating the STL ARFI algorithm. The ARF “push” is exerted at two spatial points at a particular depth, Xpi and A, which are offset by a specific distance A / ?. The resulting shear waves are then tracked at a position, XT, located at a known distance AT from Xp2. The STL ARFI algorithm is programmed at the transducer face. The “push-track” pair is laterally translated at a distance for each axial depth to generate a two-dimensional elastogram. The In-phase and Quadrature (IQ) modulated tracked dataset is software beamformed using a GPU accelerated DAS beamformer. Beamformed data is passed to a Doppler style displacement estimator, that estimates the displacement between two consecutive echo traces. Estimated displacement is processed one-of two ways (i) Frequency-domain transfer function estimator that estimates “phase delay” and “amplitude decay” over the frequency bandwidth, (ii) Time-domain crosscorrelation “group-delay” estimator that estimates time delay between displacement-peaks. The ratio of A P to the estimated time delay provides a group wave-speed estimate. The phase shift between echoes at a particular frequency provides a phase wave-speed estimate, while the amplitude decay provides an attenuation estimate. The shear wave speed and attenuation estimates are projected to a modulus space by using the plane-wave approximation of shearwave propagation to render 2D multi -frequency uE , uSfc )^ and UI.(OJ) elastograms.

[0015] FIG. 2 shows elastograms in homogeneous gelatin and oil-in-gelatin emulsion phantoms reconstructed at f-econ = 300 Hz. Phantoms PI-P fabricated with Pi: 3% gelatin. Pi: 3% gelatin, and 20 % castor-oil P3: 6% gelatin, and / A: 6% gelatin, and 40% castor-oil. The top-left sub-figure in each panel shows / / £•, top-right sub-figure shows / / s, bottom-left sub-figure shows / / L, and bottom-right sub-figure shows the magnitude of dynamic modulus fi*. All elastograms are rendered up to 25 mm in depth from the transducer face and laterally across a 20 mm field of view. The ARF was focused at a depth of 15 mm, and elasticity values are reconstructed over 8 mm x 12 mm ROI shown by the white box.

[0016] FIGS. 3A-3C shows rheological responses of homogeneous phantom models. In FIG, 3 A, shear modulus responses averaged over a bounding box around the ARF focal zone. The bounding box over which estimates are averaged is shown in white in FIG. 2. FIG. 3B shows the complex modulus \\fi*(frecon)\\ and FIG. 3C shows tan 5 responses for the phantoms.

[0017] FIG. 4 shows imaging configuration for reconstruction with varying ARF focus and reconstruction ROI. Two focal zones are labeled as fi and f2 with two reconstruction ROI’s around the focal zones. For each phantom P1-P4 the upper row shows the shear moduli reconstructed at the ROI around 20 mm focus. Each lower row shows the shear moduli for ROI around 30 mm focus. Each box plot shows scenarios where the ARF focus was set at 20 and 30 mm respectively. A schematic of the imaging scenario is shown on the top left.

[0018] FIG. 5 shows the influence of ARFI source parameters on reconstruction quality in homogeneous phantom 4. Panels (A)-(D) show the s and p.L estimates with varying ARFI F# and its corresponding coefficient-of-variation. Panels (E)-(H) show the s and .L estimates with varying ARFI transmit voltage and its corresponding coefficient-of- variation. Panels (I)-(L) show the ps and p.L estimates with varying ARFI transmit duration and its corresponding coefficient-of-variation.

[0019] FIG. 6 shows influence of ARFI source parameters on reconstruction quality in homogeneous phantom P4. Panels (A)-(D) show the ps and p.L estimates with varying ARFI F# and its corresponding coefficient-of-variation. Panels (E)-(H) show the pS and pL estimates with varying ARFI transmit voltage and its corresponding coefficient-of-variation. Panels (I)-(L) show the ps and p.L estimates with varying ARFI transmit duration and its corresponding coefficient-of-variation.

[0020] FIG. 7 shows variation of Storage (top), and Loss (bottom) modulus across frequencies between 100-800 Hz in phantom [ / fi], The box plot denotes the variance of the modulus data at a particular frequency. Black solid circles denote the mean modulus values. The mean coefficient of variation (CV) across all the frequencies is denoted at the top left.

[0021] FIG. 8 shows Independent Speckle Realization (ISR) average of modulus estimates from STL-elastography scans shown for phantom 4. The first panel shows linear modulus reconstructed across the 8 mm x 12 mm bounding-box ROI shown in FIG. 2 across five ISR averages across trials. The next set of panels show estimates for storage, and loss modulus for various discrete frequency bands of 200, 300, 400, 500, 600, 700, and 800 Hz respectively. The inset shows the coefficient of variation (CV) expressed as a percentage of variation in trial-level means for each ISR trial. The x-axis of each panel plot denotes the trial number, while the y-axis represents the modulus value in [kPa], Each bar plot denotes the range of values reported within the ROI for the particular trial, which is indicated by a circular marker. The sub-title over each panel reflects the modulus reconstructed at a discrete frequency band.

[0022] FIGS. 9A-9C show reconstructions in a viscoelastic phantom for imaging a 10 mm inclusion target in (FIG. 9A) CIRS™ breast phantom, (FIG. 9B) gelatin target in the oilingelatin background, and (FIG. 9C) oil-in-gelatin target in an elastic background. Each figure panel shows a B-mode ultrasound image and / / £•, fis elastogram images, and .L elastogram images at discrete shearwave frequencies. Average modulus values over the background and inclusion ROI, CR, and CNR metrics are used to evaluate the reconstruction quality of the target.

[0023] FIGS. 10A-10C show representative images of the viscoelastic imaging results in porcine liver ex vivo inclusion in gelatin. (FIG. 10 A) B-mode and fiE elastogram images. (FIG. 10B) fis elastogram images andelastogram images at 200 Hz, 400 Hz, 600 Hz, and 1000 Hz. (FIG. 10C) Average modulus values inside the inclusion ROI. (FIG. 10D) CR metrics are used to evaluate the reconstruction quality of the target. The analysis is shown here for a 10 mm (length) 5 mm (depth) bounding box for the background ROI and the inclusion ROI over which the mean modulus estimates are derived.

[0024] FIGS. 11 A-l ID show representative images of the viscoelastic response of a whole porcine liver subject to GA insult. Bounding-box analysis was performed within a bounding box around the landmark indicated by a red star, where GA has introduced with a hypodermic needle (FIG. 11 A) B-mode and elastograms at two time-stamps ”day-0” before and after subjecting the sample at ”day-l”. (FIG. 1 IB) Average complex modulus and elastic modulus values. (FIG. 11C) Average fis and .L values. (FIG. 1 ID) Tan d at ”day-0” and ”day-l”. The analysis is shown here for a 5 mm (length) 5 mm (depth) bounding box near the red landmark over which the mean modulus estimates are derived.

[0025] FIGS. 12A-12B show the impact of glutaraldehyde insult in porcine liver ex- vivo. (FIG. 12 A) Separation of the control group representing “day-0” and the treated group representing “day-1” post-GTA exposure to the sample. The elastogram pixels within the ROI of a single liver lobe in the 3D plane are shown here. (FIG. 12B) The experimental results were repeated on all four major lobes of a liver harvested from a single animal. Comparisons between the treated and control (n=4) are shown by the elastic and dynamic modulus reconstructed at several frequencies. Statistical significance is reported (p < 0.001) for a righttailed t-Test supporting an apparent stiffness increase after GA insult was applied to the liver samples.

[0026] FIG. 13 shows the primary lobe of a rabbit liver placed between gelatin layers over its top and bottom and scanned sagittally.

[0027] FIGS. 14A-14B show a rabbit kidney sample in the water path.

[0028] FIG. 15 shows representative images of the viscoelastic response of the whole rabbit brain embedded in gelatin. The top panels show elastograms of a single echographic scan plane’s storage ( / s), loss ( / L), complex-magnitude ( / / *), and the elastic modulus (JUE). The bottom panels show the average modulus values within the bounding-box ROI reconstructed at multiple discrete shearwave frequencies between 50-800 Hz for all echographic scan planes L1-L5. The analysis is shown here for a 10 mm (length) x 10 mm (depth) bounding box over which the mean modulus estimates are derived.

[0029] FIGS. 16A-16B shows a study of full FOV multi -frequency viscoelastic parametric imaging with ARFI shear waves in human placenta tissue. FIG. 16A shows the elastic and viscoelastic modulus at 250, 400, 650, and 800 Hz. FIG. 16B shows the average viscoelastic estimate is plotted in the ROI for all the analyzed samples.

[0030] FIGS. 17A-17B show variation of the complex shear modulus at various frequencies in a viscoleastic phantom. FIG. 17A shows variation in shear moduli measurements across independent speckle realizations (ISR) in an oil-in-gelatin emulsion. FIG. 17B shows viscoelasticity estimates obtained in a rabbit liver embedded and scanned over a gelatin cast with increasing levels of applied probe pressure.DETAILED DESCRIPTION

[0031] Reference will be made in detail to certain aspects and exemplary embodiments of the application, illustrating examples in the accompanying structures and figures. The aspects of the application will be described in conjunction with the exemplary embodiments, including methods, materials and examples, such description is non-limiting and the scope of the application is intended to encompass all equivalents, alternatives, and modifications, either generally known, or incorporated here. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. One of skill in the art will recognize many techniques and materials similar or equivalent to those described here, which could be used in the practice of the aspects and embodiments of the present application. The described aspects and embodiments of the application are not limited to the methods and materials described.

[0032] As used in this specification and the appended claims, the singular forms "a," "an" and "the" include plural referents unless the content clearly dictates otherwise.

[0033] Ranges may be expressed herein as from "about" one particular value, and / or to "about" another particular value. When such a range is expressed, another embodimentincludes from the one particular value and / or to the other particular value. Similarly, when values are expressed as approximations, by use of the antecedent "about," it will be understood that the particular value forms another embodiment. It will be further understood that the endpoints of each of the ranges are significant both in relation to the other endpoint, and independently of the other endpoint. It is also understood that there are a number of values disclosed herein, and that each value is also herein disclosed as "about" that particular value in addition to the value itself. For example, if the value "10" is disclosed, then "about 10" is also disclosed. It is also understood that when a value is disclosed that "less than or equal to" the value, "greater than or equal to the value" are also disclosed, as appropriately understood by the skilled artisan. For example, if the value " 10" is disclosed, the "less than or equal to 10" and “greater or equal to 10” is also disclosed. When two or more value are disclosed, all possible ranges between any two values are disclosed.

[0034] As used herein, the term “phase delay” refers to the time delay experienced by a propagating wave at a specific frequency component as it passes through a system.

[0035] As used herein, the term “group-delay” refers to the average transit time of a wave through a device under test as a function of frequency.

[0036] As used herein, the term “amplitude decay” refers to the gradual reduction in the maximum displacement or peak value of an oscillating or vibrating system over time.

[0037] As used herein, the term “push-track” refers to spatial locations within an imaging ensemble that illuminates a field-of-view of the imaging where the Acoustic Radiation Force is exerted and the receiver locations where the induced wave is tracked. Overview

[0038] Viscoelasticity quantification with spectral markers of wave propagation can improve pathological characterization. Acoustic Radiation Force (ARF) source Single Track Location (STL) shear wave spectroscopic imaging described herein enables reconstruction of viscoelastic signatures of tissues derived from frequency-dependent shear wave speed and attenuation. A focused ARF pushes at two spatially offset locations and acts as impulse excitations, introducing a body force into the tissue, launching shear waves from the focal push region that are tracked at a single common receiver. The cross-power spectral density (CPSD) transfer function estimates phase shift and amplitude decay as a function of propagation distance and also as a result of wave dispersion and attenuation.

[0039] Exploiting the spatio-temporal relationships of the wave propagation to the medium properties, multifrequency parametric 2D maps capturing local tissue viscoelasticity resolved over 50-1000 Hz are reconstructed. The comprehensive evaluation of this method inin-silico breast pathology models and several ex-vivo porcine livers is demonstrated. The sensitivity of image quality to acquisition parameters is studied. Viscoelastic reconstruction in STL is sensitive to tracking and excitation parameters, which deteriorate at higher frequency values. Reconstruction bandwidths are also sensitive to STL sequencing configurations and significant enough to bias reconstructed modulus values.

[0040] Elastograms of the viscoelastic complex shear modulus achieve 10-15 dB contrast enhancements at higher frequency values than the linear shear modulus image. Images can capture viscoelastic contrast in synthetic lesions and inclusions of porcine liver embedded on a gelatin background. Complex modulus elastograms also captured the spatial onset of tissue viscoelasticity alteration in an ex-vivo liver specimen after glutaraldehyde exposure. Free-hand scanning reveals stability in reconstructed modulus estimates across independent speckle realizations.

[0041] The present disclosure is based on remote palpation with ARF and measuring generated shear waves with the STL scheme. This approach interrogates viscoelastic properties in tissues up to several hundred hertz and performs imaging using standard diagnostic linear array transducers for mapping spatial heterogeneity of tissue rheology. Imaging performance is tested and validated on custom in-silico pathology models. The disclosed viscoelasticity measurement more accurately characterizes tissue rheology than existing medical ultrasound-based techniques that usually consider a tissue elastic.

[0042] The application introduces single-track location shearwave spectroscopic imaging that reconstructs complex modulus elastograms of tissues. Traditional elasticity imaging is limited to reconstructing linear shear modulus elastograms, while the proposed approach can capture both the shear storage and loss modulus elasticity images of tissues.

[0043] The application further provides study of the variation of image quality with changes in transmit and receive beamforming parameters and study the sensitivity of these on modulus reconstruction quality.

[0044] The application exemplifies the feasibility of the disclosed method by imaging elastic and viscoelastic homogeneous and heterogeneous inclusion targets and several ex-vivo small animal organs imaging (kidney, liver, and brain) and human placental specimens. The application further exemplifies the repeatability of the novel disclosed imaging method by reconstructing the viscoelastic parameters by free-hand scans on a viscoelastic phantom and on rabbit liver with increasing probe-induced pressure that mimics motion-induced artifacts in clinical and or pre-clinical imaging scenarios.Theory — Signal model for Viscoelasticity Estimation from Shear-waves

[0045] Consider two noise-corrupted shearwave time-series signals si and S2 generated by two spatially offset radiation force perturbations at Xpi and Xp2 respectively, that are probed at a common receiver at Xp. The scenario is represented in FIG 1. In order to derive a spectral estimator for spectroscopy, the study considers the two shear waves in the frequency domain represented by 5i(m) and 82(02), respectively. The study defines the crosspower transfer function as the ratio of cross-power spectral density (CPSD) of the two tracked wavefields as H (co), and the auto-spectra of the first shear wavefield Si(co). co is the angular frequency in the frequency bandwidth, represents the complex conjugate operation and c is a tunable parameter which controls the numerical stability of the denominator and is chosen as 1% of the auto-power spectrum term of the first shear wavefield.

[0046] Where the amplitude and phase of the cross-spectra transfer function is represented by |Hn| and 12 respectively. The phase spectrum and the amplitude spectrum are calculated from the real and imaginary parts of the FFT, z is the index of the spectral frequency components as shown in Equation 2.

[0047] Assuming a linear dependence of 12 with frequency, we can estimate the phase velocity (c) of the wave at a spatial location by Equation 3.> to X AP cfi>(w) = — -^12 (3)

[0048] AP is the spatial separation between the two excitation locations and is the additional distance traveled by the wave excited at Xpi compared to Xp2, when tracked at Ar. The attenuation coefficient a can be estimated from the amplitude of the normalized CPSD by Equation 4 assuming plane-wave propagation of the waves.

[0049] According to the Weiner-Khinchin theorem, the power spectral density (PSD) of a stationary random process is the Fourier transform of its autocorrelation function (ACF) [A. V. Oppenheim et al., Class Notes for 6.011 : Introduction to Communication, Control and Signal Processing, Spring 2010, Massachusetts Institute of Technology, 2010.]. In the context of two spatially offset shearwave signals, the theorem can be extended to relate the CPSDtransfer function [J / 12] and the normalized cross-correlation function (NCC) [A12] applied to si(f) and 82(f).

[0050] Where R is expressed by Equation 6

[0051] The peak time delay rpeak is the value of T in Equation 6 that maximizes 7?I2(T).The group velocity CG is related to rpeak by Equation 7.

[0052] Assuming minimal dispersion such that the phase delay spectrum 12(co) is approximated linearly around a single dominant frequency mode coo as shown in Equation 8.

[0053] Differentiating both sides of Equation 8 with respect to co we obtain a closed form solution of rpeak.Tpeak

[0054] At a single-frequency moment coo, CG ~ c<j>. Therefore, the group velocity estimate from the peak time delay of the NCC is a special case of the spectral velocity obtained using the normalized CPSD transfer function when the dispersion effect is negligible. Inserting the values of c and a to the relationship between shearwave parameters to linear viscoelastic moduli [J. Vappou, et al. Physics in Medicine Biology, vol. 54, no. 11, p. 3579, May 2009], we can convert the acoustic measurements to the storage ( / / s) and loss (HL) components of the shear viscoelastic modulus as shown in Equation 10. Both c and a are frequency-dependent estimates explicit dependence is omitted below for brevity.

[0055] The s and / / L is related to the absolute value of complex modulus [| / / *|(m)] and the loss-tangent [(tan 5] by Equation 11 and Equation 12 respectively. The complex modulus and loss-tangent are landmarks used to trace the viscoelastic behavior of a material.

[0056] The estimated group velocity from NCC is related to calculate the linear shear modulus (JUE) by Equation 13 assuming quasi-incompressibility with Poisson’s ratio 0.5.P 3pCg 3)

[0057] Frequency bandwidth of the reconstruction scenario is dictated by the sampling period of the shearwave temporal signal. Frequency resolution is controlled by the number of points used to compute the Fast-Fourier transform (FFT) of the temporal shearwave signals. A higher frequency resolution improves the accuracy of H which in turn provides more precise spectral estimates but at the cost of increased sensitivity to noise. Modulus estimates are averaged over a few frequency frames to mitigate the noise. The stability factor e controls the numerical stability of H at the cost of a slight bias. It can be shown that the estimate variance of the wavefield phase shift O12 and thus the c n is related to the recorded shearwave signal SNRs, reconstruction frequency (frecon), center frequency ( ), bandwidth ft , and ultrasound acquisition parameters like Tacqand frame-rate.Jitter of the Cross-Power-Spectral Estimator

[0058] Shearwave decorrelation between the tracked waveforms is formed primarily by the losses in the medium in which propagation occurs. The expected variance (c2), also known as the jitter magnitude of the estimated shear modulus, is based on the variance of delay estimation between the pairs of shearwave fields. In the linear shear modulus estimation, the modulus value is estimated from the peak time delay of the wavefields. A lower bound to the variance of this estimated time delay can be derived, particularly for thecross-correlation estimator assuming that an SNR for tissue velocity field estimation from ultrasonic frames and no diffraction or dispersion losses leading to shearwave decorrelation. However, this estimated variance relation holds when the material is assumed to be linear elastic.

[0059] Different parts of the wavefront travel at different speeds in a dispersive medium. The study extends the analytical form to estimate the jitter in phase delay between the two recorded wavefields that is valid for the propagation of the waves where materials are viscoelastic. Based on Equation 14 (Deffieux et al. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 59, no. 11, pp. 2390-2410, Nov. 2012), the study notes that the variance in the estimated time-delay <T2[T ] of two noise-corrupted shearwave fields recorded at two spatially offset locations but at the same depth depends upon shearwave center frequency Jo, shearwave bandwidthand acquisition duration Tacqand is given by Equation 20 a CV = -

[0060] The shear wave particle velocity SNR is defined as the bandlimited SNR and depends on the signal strength, bandwidth, and acquisition frame rate. However, for ARFI impulse excitation, the shear waves generated are sufficiently broadband. In a dispersive medium, the phase difference n must be calculated using the correct delay for each frequency component T(m). Since the phase difference is frequency dependent, then follows Equation 21®12 = to ' T(to)(21)

[0061] where m=27tf, f is the shearwave reconstruction frequency then follows Equation 22

[0062] Substituting Equation 19 into Equation 21 then obtains Equation 23

[0063] After a few algebraic manipulations and putting m=27tf, then obtains Equation 24 as a function of reconstruction frequency f.

[0064] Equation 23 shows that the variance of the phase difference estimationincreases with the square of the reconstruction frequency (f), and it decreases with increasing SNRs of either trace, center frequency ( / >), bandwidth / bw, and the acquisition time. The / o and fbw are the center frequency and -6dB bandwidth calculated from the CPSD phase spectrum, respectively. In practice, the SNR of the traces i and s are similar under the assumption that they originate from similar ARF sources and local homogeneity in the region between push-beam spacings.

[0065] Assuming co and AP are constants (i.e., their variances are zero), the variance of A012 primarily depends on the variance of 12. Using the propagation of variance for a function of a random variable, shows that Var( / dgx2VarW(25)

[0066] where Y=g(X), g() is a differentiable function and Xis a random variable. In this case,

[0067] Using equation 25 and taking the partial differentiation of c<z>(m) w.r.t 12 then follows Equation 27.

[0068] Putting back Equation 27 in Equation 25 then follows Equation 28.

[0069] Attenuation also broadens the waist of the shearwave trace, thereby creating a phase offset with respect to the reference beam also reducing its amplitude which will appear as a bias in the phase estimation. Equation 28 represents the lower bound of jitter in phasevelocity estimation from phase shift jitter between two shearwave traces generated by laterally offset pushes.

[0070] An aspect of the application is directed to a method of mapping viscoelasticity using Acoustic Radiation Force (ARF) source Single Track Location (STL) shear wave spectroscopic imaging comprising: generating shear waves in a region of interest in a subject from a plurality of spatially offset locations; tracking the shear waves at a common receiver; measuring frequency-dependent shear wave speed; measuring attenuation of the shear waves; estimating by a cross-power spectral density (CPSD) transfer function of phase shift and amplitude decay as a function of propagation distance of the tracked wave; and constructing multifrequency parametric two-dimensional maps capturing local tissue viscoelasticity.

[0071] In some embodiments, viscoelasticity is resolved over a bandwidth of shearwave frequencies.

[0072] In some embodiments, the body region is in a mammal.

[0073] In some embodiments, the spatially offset locations are in a pair that is laterally shifted sequentially along the body region.

[0074] In some embodiments, the method comprises the further step of: generating shear waves by an acoustic radiation force push that is exerted at two spatial points (a pushtrack pair) at a particular depth, Xpi and Xp2, which are offset by a specific distance Ap.

[0075] In some embodiments, the method further comprises the still further step of: tracking the generated shear waves at a position, Xp, located at a known distance AT from Xp2.

[0076] In some embodiments, the method further comprises the step of: laterally translating the push-track pair at a distance for each axial depth to generate a two- dimensional elastogram.

[0077] In some embodiments, the method further comprises the step of: beamforming using a beamformer a modulated tracked dataset that is software.

[0078] In some embodiments, the method further comprises the step of: passing the beamformed data to a Doppler style displacement estimator, that estimates the displacement between two consecutive echo traces.

[0079] In some embodiments, the method further comprises the step of: processing the estimated displacement by frequency-domain transfer function estimator that estimates phase delay and amplitude decay over the frequency bandwidth.

[0080] In some embodiments, the method further comprises the step of: processing the estimated displacement by time-domain cross-correlation group-delay estimator that estimates time delay between displacement-peaks.

[0081] In some embodiments, the method further comprises the step of: projecting the shear wave speed and attenuation estimates to a modulus space by using the plane-wave approximation of shearwave propagation to render 2D multi-frequency / / g,and L(CO) elastograms.

[0082] In some embodiments, the ratio of AP to the estimated time delay provides a group wave-speed estimate, and the phase shift between echoes at a particular frequency provides a phase wave-speed estimate, while the amplitude decay provides an attenuation estimate.

[0083] Another aspect of the present application is directed to a method of estimating quality of viscoelasticity estimates, comprising the step of estimating coherence between tracked shearwaves within an imaging ensemble.

[0084] The present application is further illustrated by the following examples that should not be construed as limiting. The contents of all references, patents, and published patent applications cited throughout this application, as well as the Figures and Tables, are incorporated herein by reference.EXAMPLESExample 1: Elastography Data Acquisition and Processing

[0085] Single Track Location (STL) elasticity imaging sequence (S. Goswami, et al. in IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 66, no. 8, pp. 12921303, Aug. 2019.) and processing was implemented in MATLAB-2020a (Mathworks, Natick, MA) and data was acquired with a L7-64LE research ultrasound scanner (Verasonics Inc., Kirkland, WA, USA) schematically shown in FIG. 1. Within each imaging ensemble, over the aperture at transducer face, the two radiation force push beams were laterally offset from each other at a fixed distance (AP ), and a common focused track line at a distance AZ was used to track resulting displacements of ultrasonic scatterers induced by the shearwaves due to the ARF ’’push”. The study experimented with several AP and AZ values to study their impact on our viscoelastic modulus estimates. The study synthesized 30 push-pairs amounting to 60 push beams that covered a reconstruction zone of 20 mm in lateral distance and 30 mm in depth. The imaging object was placed over an acoustic absorber (natural rubber pad). The homogeneous phantoms were scanned under a polystyrene compression plate mounted over the transducer face with its aperture shape cut out. The plateensures uniform compression over the rectangular test objects. The breast elasticity phantom was scanned in air, and ultrasound transmission gel (Aquasonic 100, Parker Labs) was applied to transducer face to ensure coupling. Whole ex-vivo tissue samples were scanned under a Phosphate Buffered Saline solution (Sigma Aldrich). The STL ARFI sequence used for imaging with the following parameters are listed in Table I.TABLE IBE M PARAMETERS FOR STI. ELASTOGRAPHY ACQUISITIONS.Parameters; Push TrackTransmit frequency (MHz) 3.8 5Pulse Duration ( U 50-450 0.2F # 1.5-6 2-2.5Transmit voltage (V) 10-40 10-40Focus (mm) 20-30 20-30PRE (KHz - 7Beam spacing 0.6-7.4 0.6-7.4

[0086] Based on the sequencing parameters specified above, the study acquired raw channel RF that were beamformed using a GPU accelerated Delay-and-Sum (DAS) beamforming (Perrot. Ultrasonics, vol. I l l, p. 106309, Mar. 2021). The beamformed data was processed into analytic signal using Hilbert transform, envelope-detected and subsequently displayed as B-mode image on a logarithmic scale. Additionally, In-phase and Quadrature (IQ) dataset beamformed using the Verasonics software beamformer was used to compute the shear wave displacement field (FIG. 1) by the Loupas 2D Autocorrelation algorithm (Loupas et al., IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 42, no. 4, pp. 672-688, Jul. 1995) with a kernel size of [3,2] samples in the axial and slow-time directions respectively and using the progressive referencing scheme. A bandpass filter composed of a 2nd order Butterworth architecture with an upper and lower cut-off of 1200 and 50 Hz respectively was applied to the 2D spatio-temporal particle displacement dataset.

[0087] Subsequently, elastograms were processed by the frequency domain transferfunction estimator for viscoelastogram estimation and the time-domain cross-correlation estimator for elastogram estimation respectively as shown in FIG 1. Finally, multifrequency viscoelastic modulus maps are generated at discrete ARFI excitation frequencies at a resolution determined by the number of discrete frequency samples selected while taking the FFT of the shearwave pulses. The phase and amplitude spectra of the cross-spectrum transfer function are interpolated using the P-CHIP routine in MATLAB-2023a. Based on oursampling frequency of 7 KHz, we choose 2048 sample FFT resolution and average up to 5-10 samples around a discrete frequency bin to render viscoelastic modulus images. While averaging at the higher and lower edges of our shearwave bandwidth of 50-1000 Hz, each multi -frequency rendering is weighted by inverse noise variance over the entire rendered sample to reduce the random noise in the measurements. All the above signal processing steps for viscoelastogram processing and rendering were performed offline on a dedicated node (with one NVIDIA Tesla A40 GPU with 48GB of memory, two CPU Xeon(R) Gold 6338 with 64 GB RAM) on the institutional Linux cluster (BlueHive, CIRC, University of Rochester).Phantoms & Tissue Preparation

[0088] The homogeneous viscoelastic phantoms were fabricated with oil-in-gelatin emulsion ( and TV), while the elastic phantom models ( i and P3) were created with a water-gelatin mixture. 3 and 6 percent type A, 300 Bloom gelatin (SuperClear Gelatin, Custom Collagen, Addison, IL) was mixed with distilled water to create Pi and P3 respectively, while 3 percent gelatin was mixed to 20 percent by weight castor oil (JT Baker) to create TV and 6 percent gelatin was mixed to 40 percent castor oil to create P . The molten gelatin water solution was first heated until its temperature reached 60°C. The solution was frequently mixed over a hot plate stirrer until the gelatin powder had been completely dissolved in water. A fixed percent weight of oil based on the earlier recipe was added in parts and 1 % regular dish detergent (Wegmans, NY) to emulsify the oil droplets in the oilgelatin mixture. 0.7-1 % by weight of cornstarch (Wegmans, NY) was added to enhance ultrasonic scatter, 5-6 % of 2-Propanol (Sigma-Aldrich) was added to match the speed-of- sound to the soft tissue average of 1540 ms For scanning viscoelastic inclusion targets, we first inserted a cylindrical rod of diameter 10 mm diameter in an emulsion of 4 % pure gelatin. After the background had solidified the cylindrical rod was removed and the resulting cavity was filled with 8 % oil in gelatin emulsion mixture and scanned in sagittal orientation to yield a circular inclusion on the image plane. A 10 mm lesion (composed of stiff Zerdine) that appears hypoechoic in Bmode ultrasound was also scanned in the breast elastography phantom (Computerized Imaging Reference Systems, Norfolk, VA). A fresh porcine liver was procured from a local abattoir (Joe’s Meat Market, Ontario, NY) and its largest lobe was scanned ex-vivo under PBS and also after injecting 50% aqueous glutaraldehyde (GT A) (BBC, Biochemicals) and scanned after 24 Hrs post-exposure at which the hardening had onset. Cylindrical sections of another lobe of the porcine liver of diameter of approximately10 mm was carved out using a biopsy punch and fixed in a gelatin-water base for viscoelastic scans after the gelatin had been solidified.Evaluation Metrics

[0089] The study assessed the linear and viscoelastic elastograms qualitatively by visual inspection and quantitatively using the elastographic signal-to-noise ratio (SNRE), the contrast ratio (CRE), and contrast-to-noise ratio (CNRE) in dB metrics that are defined as follows:

[0090] Where, 'LIH and OH represent the mean and S.D. of elastogram estimate within a homogenous ROI. "p and GI / B represent the mean and S.D. of elasticity estimates of the inclusion and background respectively. The percentage bias between the estimated overall viscoelastic modulus ||p*|| at a particular ARFI shearwave reconstruction frequency (frecon) and the elastic modulus (JUE) is expressed by Equation 18. 100(18)Analysis

[0091] This study developed STL-spectroscopic viscoelasticity imaging, a new technique for mapping local viscoelastic biomechanical properties of tissues. This technique captures pixel-wise shearwave phase velocity and attenuation over the ARF excitation bandwidth from cross-power spectral density (CPSD) of two laterally offset ARF pushes observed at a single spatial observation point. Estimating the phase delay and the amplitude decay of the wave at the observed spatial location the study can estimate the dispersion and attenuation profiles over a frequency bandwidth. Shearwave dispersion and attenuation measurements are then converted into complex shear modulus estimates using the planewave approximation of shear wave propagation at a fixed beam spacing. The study validates its imaging results in phantom models of tissue viscoelasticity and ex vivo animal organs. The primary findings are as follows.

[0092] In homogeneous oil-gelatin emulsions that model tissue viscoelasticity, we observe that fiE underestimates the stiffness over a significant region of shearwave bandwidth. The loss-tangent increases at higher frequencies since .L> jus.

[0093] Modulus estimates vary with ARF focal depth; the shear wave amplitude degrades with the depth due to poor ARF push intensity, and hence the variance of fis and .L are higher around a reconstruction zone located out of focus. The variance of fiE derived using cross-correlation is less than that of the spectral estimator.

[0094] Track AT and push AP beam spacings influence reconstructions of / / £•, fis and .L. The STL measurement scheme is more immune towards "speckle-bias" a systematic error that has been shown previously to occur due to preferential tracking of the "brighter" speckle. In STL the push-beam uncertainty leads to errors epi and ep2 (Equation 18). Therefore, the AP variations are shown to more influence viscoelastic modulus quantities over AT variations. Key impacts of the spacing variations leading to extension of the temporal waist of the shearwave trace that further impacts the bandwidth in frequency domain, which leads to a dependence of the parameters on the frecon.

[0095] ARFI source parameters such as the F / #, Tx voltage and Tx excitation duration impacts viscoelasticity modulus estimates. Specifically, a low F / # and higher transmission voltages and excitation durations lead to more robust / / £•, fis and .L estimations. The .L estimation is more impacted due to variations in source parameters than others.

[0096] Imaging results in 10 mm synthetic in-vitro breast lesion in FIG. 9 A and in a viscoelastic lesion in FIG. 9B show a frequency dependent increase in moduli within the inclusion target. The study observed higher contrast of the lesion from the background at higher shearwave frequencies. At higher frequencies although there is a trade-off to be made between the contrast achieved and the contrast-to-noise.

[0097] The present application provides the potential of viscoelasticity imagining in several porcine, rabbit, and human tissue specimens for reconstructing multi -frequency elasticity estimates. Imaging results in the ex-vivo porcine liver, shows that higher contrasts in the elastograms can be achieved when the moduli are reconstructed at higher frequencies. Whole livers pre- and post-GTA insult reveals an enhanced modulus after exposure. GTA has been known to alter mechanical properties in tissues. The Tan 5 shows an increase at day-1 after GTA insult is injected. The scatter plot in FIG. 12A shows that a high separation is achieved when the pixel values of the elastogram around the evaluation ROI is plotted. A 3D scatter plane of / / L, / / S at 1000 Hz and fiE between the control (day-0) shown in blue and after treatment at day-1 shown in red. The study also observes that modulus values for the controlgroup is concentrated shown by the small sphere, while for the treated group, the values are more spread.

[0098] FIG. 17A shows the viscoelastcity estimates for an oil-in-gelatin emulsion reconstructed at 100, 250, 400 and 1000 Hz acquired under acquisition events with free hand. The estimates show good agreement across different speckle realizations. FIG. 17B shows probe-induced pre-loading on viscoelastic modulus estimation in a liver. With higher initial compression, the dynamic modulus increases from 4-6 kPa at initial compression levels and gradually approaches a modulus value of 12-14 kPa.

[0099] This study introduces a novel method for simultaneously imaging the local elastic and viscoelastic properties of tissue, to improve understanding of complex tissue biomechanics. Linear elasticity imaging alone is insufficient for fully characterizing biological tissues, as they typically exhibit time-dependent stress-strain characteristics resulting in viscoelastic behavior. The challenge in probing viscoelasticity using radiation force imaging is the manifestation of time-dependent stress-strain behavior at multiple time scales. Firstly, the probe contact during a shearwave acquisition can result in viscoelastic behavior over a few hundredths of seconds. A second concurrent time-scale effect is observed due to the displacement induced by the ARF at much smaller timescales in milliseconds, during which the particle displacement caused by the push exhibits viscoelastic behavior. The study addressed imaging the ARF viscoelastic effect by assuming quasi-static contact of the probe to the imaged medium.

[0100] Imaging in homogeneous phantom models with a linear elastic assumption using the cross-correlation time-of-flight estimate reveals high percentage differences between the estimated linear elastic modulus (JUE) and the overall viscoelastic modulus (| / *(co)|) evaluated at a specific frequency co. This is especially true for phantoms with Oil (Pi and PP) that are seen in Table II to underestimate the stiffness at higher frequencies. fiE in pure Gelatin phantoms ( i and P3) overestimate the stiffness at lower frequencies that asymptotically approaches the |p*| value at higher frequencies. Hence the absolute percentage bias for i and P3 are higher at lower frequencies in Table II.Table II: Percentage bias between pE and |p*|.Phantom 200 Hz 400 Hz 600 Hz 800 HzPi 38.97 8.07 3.55 8.61P> 23.41 31.46 43.21 66.04P 39.44 26.57 12.99 11.96 4 22.73 13.07 19.84 109.63

[0101] While imaging heterogeneous imaging targets as shown in FIGS. 9A- 9C, the study observed that spectral estimates of elasticity are significantly higher than the linear elastic modulus. This is more dominant in materials where oil was used to fabricate the phantoms and in ex-vivo livers. The frequency cutoffs that limit the feco for the elastogram images are dictated by the loss of the shearwave signal SNR at higher frequencies.

[0102] The study observed that the reconstruction range of 100-1000 Hz to be most reliable for elastogram reconstructions. At the low-frequency cutoff, minimum wavelength detection considerations limit the frecon of the elastograms fis and fiL. Specifically, phantoms prepared with higher percentage weight of gelatin P3 and P have lower modulus SNR compared to Pi and 2 that are less stiff (FIG. 3 A), hence possess shorter wavelengths at lower frequencies than phantoms with higher stiffness at lower frequencies where wavelength is larger. The high-frequency cut-off is due to onset of attenuation which limits the number of valid pixel estimate in the elastograms. Additionally, shearwave attenuation is estimated from the amplitude of the CPSD of the shearwave signals, which is prone to enhanced levels of measurement noise than the CPSD phase that is used to estimate phase velocity, hence fiL estimates are corrupted more than fis and fiE. Thus, the effect of noise is more evident in the shearwave amplitude estimates.

[0103] The present example of studying the viscoelastic signature of porcine liver inclusions in gelatin shows that higher contrast between the liver and gelatin can be achieved with fis and elastograms at higher reconstruction frequencies. The artifacts in these elastograms may be attributed to the attenuation effects of the wave when propagating through the porcine liver. The liver tissue’s fis and .L are both increased by exposure to GTA at day-1, as shown in FIGS. 11 A-l ID. The region around the site where GTA was injected (indicated by the red landmark) shows enhanced modulus values. Through a comparison of the elastograms prior to and following GTA treatment, the study was able to identify the treated area using the red region indicating higher modulus values at various discrete frequencies. The mean elastic modulus fiE increased from 2 KPa to 15 KPa post-exposure, while | pi* | increases from 30 KPa at ”day-0” to 50 KPa at ”day-l” at the highest fecon of 1000 Hz. When these values are plotted in a scatter plot (FIG. 12 A) a clear separation in the shear modulus plane [« / ., fis~ 1000, / / L-1000] can be visualized.

[0104] The study also performed the spectral viscoelastic imaging method with free-hand scans, the estimates are conserved across all the independent speckle realization acquisitions. Highest variation in acquisition event is observed in .L at 1000 Hz.

[0105] Estimation of shearwave attenuation a(m) is difficult due to noisy measurements of the amplitude decay of the CPSD. This leads to a high coefficient of variations of .L as shown in the measurements. This is improved in the present application using a physics-driven filter to correct non-physical estimates of shear loss modulus. fiE and .L are not simply standalone estimations but, in fact, are intertwined by Hilbert-transform pairs of each other. Further embodiments include a framework to impute non-physical noisy estimates using a Hilbert transform filter to improve elastogram reconstructions over a larger FOV.

[0106] Comparisons of the complex modulus reported by STL-SWEI (between 50-800 Hz) with the bulk viscoelastic modulus interrogated at sub-hertz frequencies (0.01-15.9 Hz) by the DMA revealed a similar trend in increase in modulus measurements with frequency.

[0107] Accordingly, the present application provides STL spectroscopic viscoelastic imaging, a novel spectral imaging technique for mapping tissues’ local viscoelastic biomechanical properties using shearwave dispersion and attenuation parameters. Prior work only reconstructed the frequency dependent shear modulus from the dispersion slope of the phase velocity, thereby underestimating the viscoelastic modulus due to the omission of the wave attenuation phenomenon. Concurrent wave dispersion and attenuation estimation have been limited to spectroscopy, which reports an average estimate of tissue modulus over an ROI. Local analysis of viscoelastic properties over a complete field of view had not been studied extensively. The present method extracted the shearwave dispersion and attenuation measurements from the spectral shift in phase and amplitude of the cross-spectrum of wave generated by ARF at two offset locations. Due to a common probe location, particle displacement estimation suffers less “speckle-bias” effects. The present validations in synthetic viscoelastic phantom models and ex-vivo porcine liver tissues demonstrate the feasibility of this technique, surprisingly resulting in 10-15 dB enhancements in viscoelastographic contrast while imaging the inclusions in loss and storage modulus elastograms at higher frequency ranges of the ARFI bandwidth. STL spectroscopic imaging reveals the contrast between treated and untreated samples of tissues when subjected to chemical glutaraldehyde insults. When visualized in the viscoelastic modulus plane, a clear separation is observed between co-located pixels at “day-0” before exposure and at “day-1” post-exposure. The present application provides a useful method to map complex biomechanical interactions in viscoelastic tissue pathologies that can be applied in preclinical and clinically relevant imaging scenarios.Example 2: Homogeneous Viscoelastic Phantoms

[0108] The study independently quantified the effects of varying the gelatin and oil concentrations on the estimated stiffness. The study also investigated the impact of STL imaging parameters on modulus estimates and image reconstruction quality.Reconstructing Complex Modulus Elastograms in Gelatin Emulsions

[0109] FIG. 2 shows the elasticity maps for homogeneous phantom models. The cross correlation estimator calculated the elastic shear modulus uE while S, ul., and \ *| are estimated from the phase and amplitude of the cross-power spectral density described previously. Qualitatively, the study observed that while imaging at 300 Hz, phantoms fabricated with Oil in Gelatin emulsion reveal enhanced contributions of the loss modulus JUL to the overall stiffness. A reconstruction ROI of 8 mm x 12 mm (delimited by a white bounding box) was selected around the focal zone of the ARF source. Quantitatively, the contributions of both / / s and JUL were pronounced when the samples were imaged at higher frequencies, especially for the Gelatin-Oil emulsion phantoms as shown in top-row of FIG. 3 A. Between, phantom Pi and P3, that differed in gelatin concentration (going from 3% gelatin to 6 % gelatin), / / s increased from 1.5 kPa to 4.2 kPa at 200 Hz, while / / L increased from 1.2 kPa to 2.2 kPa at 200 Hz, while the / / L increased from 0.6 kPa to 3.2 kPa between phantoms 2 and P going from 20% castor oil to 40 % castor oil respectively. The study observed a frequency-dependent increase in the stiffness of the estimated complex modulus | / / *(co)|. fiE (indicated in blue), is derived from the group estimate of the propagating wave and hence the fiE is the group shearwave estimate approximately at the mid-band frequency of the shearwave. The study observed that fiE underestimates the stiffness at higher frequencies. The cross-over between the fiE and the fis modulus appears earlier at lower frequencies (150-200 Hz) in gelatin-oil emulsions as observed in phantoms 2 and P as opposed to gelatin-water emulsions Pi and P3 where this cross-over is delayed (550 Hz). The second row of FIG. 3 A, shows the SNR analysis of the elasticity reconstructions in homogeneous phantoms. The SNReof / / s and / / L deteriorate at both higher and lower frequencies, except for the phantom Pi where we observed an enhancement of SNR for the / / s elastogram at higher frequency. SNRein phantoms PI-PA peaks between 200-500 Hz and falls- off from their respective ranges in individual phantoms beyond those frequencies. As an example for PA fis and .L SNReachieves a peak value of 6 and 4 respectively at 300 Hz and beyond the frequency bandwidth of 200-500 Hz, quickly degrades to less than half its value. This frequency bandwidth is then the usable reconstruction bandwidth for the corresponding viscoelastograms of phantom PA.

[0110] The rheological visualization also shows a second landmark indicating crossover between fis and the fiL. This is known as the “gelation-point” beyond which fis > .L. The cross-over frequency between fis -gL appears to be delayed in Gelatin-Oil emulsion samples (beyond 800 Hz), while for Gelatin-Water emulsion this cross-over appears earlier at 100-300 Hz. The loss-tangent (tan 5 = fiiJ JUS) is considered as a measure of energy dissipation in the medium, tan 5 as shown in Fig. 3 (C) is high in pure-gelatin phantoms (red and blue) since the ns is lower at lower frequencies, while at the higher frequencies around 700-800 Hz pL < Us so the study observes a decrease in tan 5. In P2 and P , however, fiL is greater than fis at higher frequencies; hence an increase in tan 5 is observed around 700-800 Hz.

[0111] All phantoms exhibit an increasing trend in dynamic modulus with frequency which is referred to as one of the two characteristic viscoelastic landmarks of a material as shown in FIG. 3B. P shows the steepest rise, reaching up to 20 kPa at 800 Hz since its composed of the highest composition of both gelatin and oil. P2 and P3 demonstrate intermediate stiffness around 10 kPa. Pi is the least stiff since its fabricated with the lowest gelatin and oil concentration, with the highest dynamic modulus reported for i at 2 kPa.

[0112] The loss tangent (tan b) quantifies energy loss in the medium and its frequency signature is the second viscoelastic landmark of a material, tan b, as shown in FIG. 3C, is high in pure-gelatin phantoms (red and blue) since the fis is lower at lower frequencies. In contrast, at the higher frequencies around 700-800 Hz .L < is, so we observe a decrease in TanA At 800 Hz, the values of Tam) are 0.2 for i and 0.3 for P3 respectively. In P2 and P4, however, .L is greater than fis at higher frequencies; hence an increase in Tarn) is observed around 700- 800 Hz. At 800 Hz, the values of Tan 4 are 0.8 forP2 and 1.2 forP4 respectively, suggesting viscoelastic domination for the oil-in-gelatin phantoms.Focal depth and reconstruction ROI

[0113] The elastic and viscoelastic stiffness estimates reconstructed at 200 Hz for phantoms Pi -P4 imaged at two focal configurations are shown in FIG. 4. Given that the study examines homogeneous phantoms, differences between the two groups are not anticipated, except for imaging artifacts. The two groups labeled in the box plots as 20 mm and 30 mm represent the location of the ARF focus. The measurement scheme on the upper left consisted of two 8 mm x 12 mm reconstruction zones for the upper and lower rows in each set. Each upper row is the elasticity measurement with the ROI centered at a depth of 20 mm from the scan surface, while each lower row describes the measurement with the ROI centered around 30 mm from the transducer face.

[0114] The study observed four focal-zone and ROI pair configurations for each estimated modulus parameter for every phantom model. Overall, the study observed that when the reconstruction ROI is located within the focal zone, the modulus estimates are higher; also, the variance increases with the concentration of oil and gelatin. Higher frequency estimates of the modulus report lower usable bandwidth of the CPSD amplitude spectrum hence the modulus estimates report higher variance.

[0115] The wave quality degrades beyond the focal zone hence for reconstructions outside the focus degrade severely. For phantoms PI-P graded in Gelatin and Oil concentrations, the study observed that spectral estimate of modulus at 200 Hz is noisy around the ROI, since the shear wave intensity beyond the focal zone is low, the spectral SNR degrades, and hence modulus estimates show a high variance.

[0116] Mean values of fiE remain consistent between ROI sizes and focus configurations. We observe that when the reconstruction ROI is located within the focal zone, the modulus estimates are higher; additionally, the variance increases with the concentration of oil and gelatin, evident by the spread of the data around the mean. Higher frequency estimates of the modulus obtained beyond the usable bandwidth of the CPSD amplitude spectrum show higher variance.

[0117] Shear wave amplitude diminishes with distance from the focal zone and depth. Elasticity reconstructions outside the focus and at the 30 mm depth are noisier than at the 20 mm depth. Across all phantoms P1-P4, we observe that the variance of the viscoelastic modulus at 200 Hz is higher compared to fiE. We match the ARF push and track focus in this set of experiments; hence, when the reconstruction ROI is not set at the push or track focus, we observe an increase in estimate variance shown by increased spread in data around mean for the blue data points in upper row panels compared to red data points and vice-versa for lower row panels for all the phantoms.“Push” & “Track” geometry

[0118] This study examined the impact of beam spacings in an STL ensemble and the extent to which AP AT influence the variance of arrival time estimations of the shear wave at the tracking locations. This study examined such variations of the viscoelastic phantom / A as presented in FIG. 5. In the STL scheme, the shear wave speed CSTL is estimated between two push locations Pi and P2. The extra propagation path traveled by the shearwave from push beam 1 compared to push beam 2 is AP (not shown). The wave speed under a linear elasticity assumption is estimated as a ratio of the path length difference over the arrival time difference between pushes, as shown below in Equation 18.

[0119] The estimated arrival times induced by the two pushes TPI and TP2 with associated estimation errors epi and ep2. These errors are due to the misalignment of the push beam, while er is the error due to the misalignment of the track beam. Since, the study was tracking ultrasonic displacements at a common location, the STL scheme is resilient towards preferential tracking of the bright speckle therefore the eT is common to both expressions in the denominator of the expression in Equation 18. Panel (A) of FIG. 5 studies the variation in reconstruction quality with tracking beam separations (AT). The push-beam spacing, in this case, AP , is fixed at a large value of 24 (7.4 mm), and AT is varied from 2-24 (0.6-7.4 mm). A large AP ensures that we provide enough propagation distance to the wave generated at 1st location, so that near-field diffraction or geometric spreading effects do not bias the estimates and that we just study the impact of AT . The viscoelastic modulus parameters are then plotted as a function of the ARF excitation frequency and the beam separations. Shown here is the plot of the shearwave frequency on the x-axis and AT on the y-axis. An enlarged track beam separation increases the effective propagation distance of the shearwaves generated by the ’’pushes” Pi and . In the time domain, an increase in AT shifts the wave’s leading edge at a further in time for a fixed acquisition event time, ensuring that it is adequately sampled, thereby avoiding leading-edge truncation. Due to the shift, the spectral content of the wave is changed.

[0120] The study presents results of reconstructed fis andfor an ROI of height 6 mm x and width 12 mm around the focal zone situated at a depth of 20 mm at each such shear wave frequency fecon and AT combinations. The study observes that gs and gL estimates increase with AT and frecon. The CV of estimates increases with f recon and decreased with increase in AT. At 200 Hz, the study observes that fis and .L estimates increase with AT and frecon. The coefficient of variation (CV) of the estimates increases with frecon and decreases with an increase in AZ. From Panels (B) and (C) of FIG. 5 it was observed that the viscoelastic modulus ( / s, fif) reconstructed at 200 Hz and the elastic modulus fis) increases with AT up to 2.4 mm after which stable estimates (decreased sensitivity) of the modulus parameter at higher beam spacings were observed. At narrow AT spacing (0.6 mm), the leading edge of the wave from the later push P2 is truncated while tracking. In contrast, for the widest spacing (7.4 mm), the tracked wave undergoes widening, which shifts the spectral content, effectivelyreducing the usable bandwidth of the CPSD and, therefore, the spectral range over which the modulus is reconstructed.

[0121] As shown in Panel (D) of FIG. 5, the CV is lowest for / / £• (< 0.2) for track spacings between 2.4-7.4 mm, and the .L at 200 Hz is more impacted by the choice of spacings reflected by a higher CV (0.3) across the parameter space. Higher reconstruction frequencies are more severely impacted with CV of .L > 2 for the largest AT.

[0122] The “speckle-bias” albeit lesser in magnitude associated with the STL scheme is due to variance induced by changes in the spacing between push beams (AP ). In this experiment the study varied push beam spacings AP between 0.6-7.4 mm while AT remained fixed at 4.9 mm. Panels (E-H) of FIG. 5 show that the apparent moduli decrease with AZ except for 200 which shows a slight increase, while the decrease in fiE and « / .2oo is steeper. The sensitivity of fiE to AP is comparable to AT, both showing a CV below 0.2 across the parameter space. The CV for fis andis around 0.3 for AP, which is slightly lower than the corresponding CV for AZ, which is 0.1. This is evident when comparing Panel (H) of FIG. 5 to Panel (D) of FIG. 5. Additionally, CV decreases with increasing AP hence elastogram image uniformity is achieved at higher beam spacings. Higher CV of the estimates > 0.5 is observed at lower AP that increases at higher freCon (>1.5). An increase in AP leads to a time shift of the tracked displacement from Pi and broadens its temporal extent, which results in a decrease in bandwidth. Although a smaller AP is desirable to maintain a good spatial resolution, yet a much lower AP leads to an increase in noise and to avoid the shear wave trace from Pi to be truncated.ARFI source parameters

[0123] The study varies the parameters of the ARF source controlled by the F / # of the aperture, the transmit (Tx) voltage and its duration driving the ARF "pushes." As shown in panels A-D (top-left) of FIG. 6, this study varied the ARFI F / # on the viscoelastic phantom P4 and observe that when the focus zone is set at 20 mm, track beam spacing is set at 3.1 mm, and the push beam spacing is set at 7.4 mm, the modulus estimates remains consistent across the range of F / # studied here between 1.5 to 6. However, the CV increases with reconstructions with larger F / #. A larger push F / # means a smaller push aperture for a focal depth, resulting in a weaker force and lower shear wave amplitudes, yielding more variance in the estimate of the viscoelastic modulus. We report that with at 1.5 F / #, we can reconstruct the for « / .OOwith a lower CV of 0.25 compared to F / # 6 where CV is 0.65. F / # also impacts the frequency content of the wave generated by the push. A tightly focused beam generates amore localized deposition of acoustic energy, producing broad spectral bandwidth, thereby supporting multi -frequency stiffness estimation.

[0124] Next, the study investigated the impact of ARF transmit (Tx) excitation voltage and its duration on viscoelastic reconstructions; the Tx voltage studied here is changed between 10-40 Volts. Higher Tx voltages produced a more stable estimate of both the moduli as seen in panels E-H of FIG 6. For voltage studied between 20-40 Volts, the CV is below 0.5. A lower Tx voltage yields a weaker push, leading to lower shear waves in amplitude and reduced bandwidth. Therefore, we observe that the viscoelastic moduli are better reconstructed at higher Tx voltages for higher shearwave frequencies. Similarly, increasing Tx duration (Panels I-L) to 440 / / s push yields a highest CV of 0.21 for « / .QlQlcompared to a maximum CV of 0.8 for a 40 / s push.Variations of viscoelastic estimates in homogeneous phantoms across frequencies and independent speckle realizations

[0125] FIG. 7 shows the variation of viscoelastic modulus components, Storage (top) and Loss (bottom) modulus reconstructed across the frequency range of 100-800 Hz in phantom P . We observe that the overall trend that the moduli increase over the frequency bandwidth. The variance of the modulus increases with frequencies and is lowest between 100-300 Hz, thereby supporting the findings in FIGS. 3A-3C, where the SNR of a similar analysis was found to be the highest. The reported mean coefficient of variation (CV) for the storage modulus is slightly lower at 97.2% compared to 98.9% for the loss modulus.

[0126] The study also reports viscoelasticity estimations across five Independent Speckle Realization (ISR) trials, in order to demonstrate the repeatability of measurements and the viscoelastic estimation algorithm. For each such trial the transducer was physically lifted from previous scan position and located at a different scan area within the FOV. This resulted in a completely different speckle realization over the scanned image. FIG. 8 shows the acquisitions across ISR trials on viscoelastic phantom P . The study observes that the trial-to- trial variation in the linear modulus estimate measured by the coefficient-of-variation (CV) of the trial-level means is the least at 3.7%, while it increases progressively with frequency. The CV for [is and [IL at 200 Hz is 0.5 %, and 16.3 % respectively, which increases up to 2.4 % for [is at 500 Hz and 25.4 % at 700 Hz. Generally, for all frequency bands studied here we have observed that the CV of [IL is higher than that of reported CV of [is.

[0127] Both these observation, that higher reported CV with frequency and higher CV for [IL compared to [is is in line with our previous observations of lower SNR of elastogramsreconstructed at higher frequencies and that of loss elastograms reporting lower SNR compared to storage elastograms.Imaging In-vitro inclusion targets in synthetic phantoms

[0128] FIGS. 9A-9C show the results of imaging inclusion targets. Imaging was performed with B-mode guidance to locate an optimal echographic view. The linear elastogram is shown by fiE. The corresponding viscoelasticity images of / / s, / / L, and | / / *| are shown for frequencies between 50 Hz and 1000 Hz. Background and inclusion ROI are appropriately demarcated for computing the average modulus estimates over frequency and evaluating Contrast Ratio (CR) and contrast-to-Noise Ratio (CNR) of reconstructions without additional image smoothing filters.

[0129] FIG. 9A shows representative images of a 10 mm diameter stiff lesion in a breast phantom model (CIRS Model 059). We observed that the background is predominantly elastic; therefore, fiE represents the characteristic stiffness behavior. The inclusion fis shows enhancement over the frequency range relative to the background. The CR reveals a contrast enhancement relative to the background with theover the interrogated frequency range reaching up to 20 dB, while the CR for fis and the overall modulusis lesser and remains constant around 10 dB. The spectroscopy-processed images appear noisy, and hence, the fiE shows the highest CNR of -40 dB compared to the fis and

[0130] FIG. 9B demonstrates imaging of an elastic gelatin target in a viscoelastic oilin-gelatin background. The scenario mimics a situation where the target and background have similar fiE or, in other terms, matched in conventional elasticity measures. However, their respective viscoelastic moduli are different. The mean fiE of the inclusion target is 11 kPa, while that of the background is 15 kPa. The inclusion is visible in the .L elastogram at 700 Hz and 1000 Hz. The | / / *| and .L CR is -2 dB, compared to fiE which is -3.75 dB and therefore higher than the fiE elastogram between 400 800 Hz, offering better target detectability by the viscoelasticity reconstruction.

[0131] FIG. 9C shows the imaging results of a stiff viscoelastic oil-in-gelatin target embedded in an elastic gelatin background. The mean fiE of the inclusion target is 15 kPa, while that of the background is 4 kPa, thus approximately 3.5 times stiffer than its background. The inclusion’s mean fis reaches up to 20 kPa, while the background is 6 kPa at 1000 Hz. The .L is 5 KPa for the inclusion, while 2 KPa for the background is between 400- 100 Hz. The .L best reconstructs the inclusion target between 200-400 Hz. The CR of the reconstructions is low for lower reconstruction frequencies and increases as fis and .L arereconstructed at higher frequencies. However, the CNR decreases when the elastograms are reconstructed at higher frequencies.Ex-vivo animal imaging

[0132] FIGS. 10A-10D show a 10 mm x 5 mm section of a freshly excised porcine liver embedded in a solid gelatin-water emulsion. The 1E elastogram of the liver section is shown in FIG. 10 A. FIG. 10B shows the s (top) and / L (bottom) images of the liver. The boundary box analysis of the biomechanical character of the liver is shown in FIG. 10C. The average fiE within the black ROI that demarcates the liver is 16 kPa (pink). The fis modulus (red) and the .L modulus (green) increase with reported shear wave frequencies up to 1000 Hz. The apparent fis andat 700 Hz are 23 kPa and 7.8 kPa, respectively. The CR of the liver inclusion target evaluated against its background is shown in FIG. 10D. We note that fiE offers a contrast of 10 dB while both fis and .L contrast up to 15 dB and 25 dB, respectively, at 1000 Hz.

[0133] Next, the study imaged the viscoelastic response of a whole porcine liver in a saline water bath, before and after subjecting the sample to the GA insult, as shown in FIGS. 11A-11D. FIG. 11A shows the results in a single porcine liver lobe that was imaged on “day-0,” which means the day it was secured fresh after harvest. The mean fiE at baseline evaluated around a 5 mm x 5 mm ROI at a depth of 10 mm is 2 kPa as shown in the right y- axis ofFIG. 11B. The mean complex modulusshown on the right y-axis of fig 11B at “day-0” is around 10 kPa over a bandwidth of 200-600 Hz and rises to 30 kPa at 1000 Hz. At “day 1”, meaning that the GA was injected and imaged after 24 hrs, fiE is around 15 kPa. The is 30 kPa over 200-600 Hz and increases up to 50 kPa at 1000 Hz post-GA exposure. 5 ml of GA was injected in the tissue at the indicated landmark (red star) using a hypodermic needle. As shown in FIG. 11C, we observe that the native porcine liver at “day-0” the jus evolves between 5-10 KPa, while UL evolves from 0.5-10 KPa between 50-1000 Hz.

[0134] The gelation point of the native liver at “day-0” is at 1000 Hz, while at “day- 1” it drops to 775 Hz in the frequency bandwidth. At “day-1”, the / / s increases up to 25 kPa between 200-400 Hz and then decreases, while / / L initially increases 8.5 kPa and then decreases 400-600 Hz and subsequently increases up to 1000 Hz. FIG. 1 ID shows the tarn) pre and post-GA insult; the loss tangent increases post-exposure at higher frequencies, suggesting modulation of viscolasticity of the sample induced by the GA.

[0135] FIG. 12A represents the separation achieved on the viscoelastic hyper-plane when examined pixels are plotted over the analyzed ROI pre and post-GA insult painted inblue and red, respectively. The three axes represent « / .), us1(mHzand LWMHZevaluated at 1000 Hz respectively. FIG. 12B shows the statistical analysis performed on the mean modulus reported within the ROI for viscoelastic analysis of all four lobes of the sample. A paired t-test showed that elastic and dynamic modulus at multiple frequencies showed that the treated group had significantly larger mean values than the control group.

[0136] FIG. 13 shows representative images of the viscoelastic response of the primary lobe of two rabbit livers. Panels show elastograms of the storage ( / s), loss ( / L), complexmagnitude ( / / *), and the elastic modulus (JUE). Modulus values averaged within the white bounding box ROI and reported as mean over standard deviation over multiple shearwave frequencies. For Rabbit: 1, the mean fiE is 4.34 kPa. For Rabbit:2, the mean fiE is 3.41 kPa.

[0137] FIGS. 14A-14B show representative images of the viscoelastic response of rabbit kidneys. FIG. 14A shows elastograms of the storage fis (top) and loss .L (bottom), as well as the elastic modulus fiE. FIG. 14B shows the viscoelastic modulus values averaged within the white bounding box ROI and reported as mean over standard deviation over multiple shearwave frequencies for the samples tested. The analysis is shown in Panel (B) for a 10 mm (length) 5 mm (depth) bounding box over which the mean modulus estimates are derived. The apparent mean fiE is as follows. Rabbit: 1, Kidney 1 : 19.4 kPa. Rabbit: 1, Kidney 2: 11.1 KPa. Rabbit:2, Kidney 1 : 20 kPa. Rabbit:2, Kidney 2: 8.24 kPa respectively.

[0138] FIG. 15 shows a whole rabbit brain scanned while embedded in a gelatin construct that is scanned across various echographic planes referred to as Z1-Z5. For all the ex-vivo scans, the study reconstructs the storage, loss, complex-magnitude, and elastic modulus elastograms shown in the panels at a few discrete frequency bins, namely 250 Hz, 650 Hz, and 1000 Hz, respectively. The white bounding box in all the elastogram images was used to perform aggregate viscoelasticity analysis over the ROI shown on most right panels, which reported the mean over the standard deviation of reported modulus values. The apparent viscoelastic moduli within the ROI increases with the shearwave reconstruction frequency analyzed here, which is between 50-800 Hz for the bounding-box report, revealing their viscoelastic behavior in the frequency range while interrogated by ARFI.Imaging human placental samples

[0139] FIGS. 16A-16B show viscoelastic characterization results in human placental specimens. The placenta is a key organ that supports the full pregnancy term by acting as a medium for maternal-fetal transport. MR and Doppler ultrasound are established modalities, while elastography is also an emerging tool for placental functional imaging. Several priorstudies have reported the presence of shearwave dispersion in placental tissues ex-vivo. The present inventors have also demonstrated that spatial mapping of the placental tissue by imaging may aid disease state discrimination and analyzing its structural characteristics in health and disease.

[0140] Hence, full FOV multi -frequency viscoelastic parametric imaging with ARFI shear waves is useful to study localized landmarks over the tissue guided by the B-mode and elastic modulus elastogram fiE for viscoelasticity analysis. Viscoelastic modulus quantification in previous studies was limited to Kelvin-Voigt-based shearwave dispersion analysis. However, the local frequency reconstructions of the tissue’s viscoelastic modulus have not been reported. The present example reports results from four human placental tissues imaged in a temperature-controlled saline bath. The images are captured with the fetal surface facing the transducer. The ARF is focused midplane between the fetal and maternal surfaces of the tissue at a depth of 10-15 mm from the surface. FIG. 16A shows the elastic and viscoelastic modulus at 250, 400, 650, and 800 Hz. The average viscoelastic estimate is plotted in the ROI for all the analyzed samples (FIG. 16B).Impact of probe induced pre-compression and freehand scanning

[0141] FIG. 17A shows the variation in shear moduli measurements across independent speckle realizations (ISR) in the oil-in-gelatin emulsion P with each free-hand acquisition event by a single operator of the probe. Each scan event involved physically lifting the probe and scanning a different region of the phantom, including varying scan planes in every acquisition episode. The estimates are averaged around a 10 x 10 mm ROI laterally and at a focal depth of 20 mm. As shown, for each trial we report the dynamic modulus at 250-750 Hz is between 10-12 kPa, while for at 1000 Hz the dynamic modulus is greater than 20 kPa for all trials. The variation in the complex modulus | / / *| is higher than that in the elastic modulus / / £•, and the estimate variance is higher for the complex modulus as the reconstruction frequency increases as indicated by the error bar over the bar plots. Nevertheless, we observe that stable estimates of viscoelasticity can be obtained and potentially clinically realizable under freehand operator deformation.

[0142] FIG. 17B shows viscoelasticity estimates obtained in a rabbit liver embedded and scanned over a gelatin cast with increasing levels of applied probe pressure. In this scenario, we acquired ultrasound data at the exact location without lifting and compression exerted by a mechanical 6-degree-of-freedom positioner. The induced compression was restricted to the small strain regime. Strain stiffening behavior in soft tissue mechanics of the liver and brain has been documented under compression.

[0143] Also, elastography results have reported an increase in mean elastic modulus value evaluated via the linear approximation. Similar to such a trend, a general increase can be observed in the magnitude of the complex viscoelastic modulus with an increase in pressure. In particular, the largest enhancement of the estimated viscoelastic stiffness is reported at higher shearwave frequencies and the highest probe induced pressure exerted on the ex-vivo tissue. At 800 Hz, the study reports an 81% increment in the reported dynamic modulus increasing from 7.5 kPa up to 14 kPa at the highest probe-induced pressure.Increment in viscoelastic modulus with probe compression suggests that the particular tissue demonstrates strain stiffening, thereby we may potentially capture greater contrast in viscoelastographic elastograms when the tissue-of-interest is strain stiffening and exists in-situ in a background region that does not demonstrate such behavior. The apparent nonlinear enhancement in reported viscoelastic modulus with induced pressure can be clinically relevant in scenarios such as breast lesions where malignancy has been found to be correlated to elastic nonlinearity of the elastography reported stiffness. Also notable is that probe-induced compression may also create anomalous results when reporting absolute stiffness values without controlling for the induced pressure.

[0144] While various embodiments have been described above, it should be understood that such disclosures have been presented by way of example only and are not limiting. Thus, the breadth and scope of the subject compositions and methods should not be limited by any of the above-described exemplary embodiments but should be defined only in accordance with the following claims and their equivalents.

[0145] The above description is for the purpose of teaching the person of ordinary skill in the art how to practice the present invention, and it is not intended to detail all those obvious modifications and variations of it which will become apparent to the skilled worker upon reading the description. It is intended, however, that all such obvious modifications and variations be included within the scope of the present invention, which is defined by the following claims. The claims are intended to cover the components and steps in any sequence which is effective to meet the objectives there intended, unless the context specifically indicates the contrary.

Claims

WHAT IS CLAIMED IS:

1. A method of mapping viscoelasticity using Acoustic Radiation Force (ARF) source Single Track Location (STL) shear wave spectroscopic imaging comprising: generating shear waves in a region of interest in a subject from a plurality of spatially offset locations; tracking the shear waves at a common receiver; measuring frequency-dependent shear wave speed; measuring attenuation of the shear waves; estimating by a crosspower spectral density (CPSD) transfer function of phase shift and amplitude decay as a function of propagation distance of the tracked wave; and constructing multifrequency parametric two-dimensional maps capturing local tissue viscoelasticity.

2. The method of Claim 1, where viscoelasticity is resolved over a bandwidth of shearwave frequencies.

3. The method of Claim 1, where the body region is in a mammal.

5. The method of Claim 1, where the spatially offset locations are in a pair that is laterally shifted sequentially along the body region.

6. The method of Claim 1-5, comprising the further step of: generating shear waves by an acoustic radiation force push that is exerted at two spatial points (a push-track pair) at a particular depth, Xpi and Xp2, which are offset by a specific distance Ap.

7. The method of Claim 6, comprising the further step of: tracking the generated shear waves at a position, Xp, located at a known distance AT from Xp2.

8. The method of Claim 7, comprising the further step of: laterally translating the push-track pair at a distance for each axial depth to generate a two-dimensional elastogram.

9. The method of Claim 8, comprising the further step of: beamforming using a beamformer a modulated tracked dataset that is software.

10. The method of Claim 9, comprising the further step of: passing the beamformed data to a Doppler style displacement estimator, that estimates the displacement between two consecutive echo traces.

11. The method of Claim 10, comprising the further step of: processing the estimated displacement by frequency-domain transfer function estimator that estimates phase delay and amplitude decay over the frequency bandwidth.

12. The method of Claim 10, comprising the further step of: processing the estimated displacement by time-domain cross-correlation group-delay estimator that estimates time delay between displacement-peaks.

13. The method of Claim 1-12, comprising the further step of: projecting the shear wave speed and attenuation estimates to a modulus space by using the plane-wave approximation of shearwave propagation to render 2D multi -frequency E, ps(co), and PL(CO) elastograms.

14. The method of 6-13, wherein the ratio of AP to the estimated time delay provides a group wave-speed estimate, and the phase shift between echoes at a particular frequency provides a phase wave-speed estimate, while the amplitude decay provides an attenuation estimate.

15. A method of estimating quality of viscoelasticity estimates, comprising the step of estimating coherence between tracked shearwaves within an imaging ensemble.

Images