Elastography: strain, ARFI and shear-wave imaging
Elastography images a mechanical property that B-mode cannot see: the shear modulus, which spans five orders of magnitude across soft tissues while the bulk modulus barely changes. Every method has two halves: a way to deform the tissue (hand compression, acoustic radiation force, an external vibrator) and a way to measure the resulting displacement (the speckle-tracking estimators of Module 11, run at kHz frame rates). Shear-wave imaging adds a third half that decides the accuracy: the estimator that turns a space–time displacement field into a speed, and the assumptions (elastic, isotropic, bulk, non-dispersive) that turn the speed into kPa.
You will be able to
- Convert between \(E\), \(\mu\) and \(c_s\) with \(\mu = \rho c_s^{2}\) and \(E \approx 3\mu\), and state every assumption that conversion makes
- Predict strain-image SNR, contrast and resolution from kernel length, LSQ span and applied strain (the strain filter), and explain why strain ratios underestimate modulus ratios
- Compute the radiation force \(F = 2\alpha I/c\), the resulting ARFI displacement in µm, its recovery time and the heating per push
- Design and critique shear-wave speed estimators: time-to-peak and cross-correlation regression, k–f phase velocity, local wavenumber, directional filtering, and their resolution–precision trade-offs
- Separate phase from group velocity in viscoelastic tissue (Kelvin–Voigt, power law) and recognize guided-wave (Lamb) dispersion in arteries, cornea and thin layers
Assumes
- Module 02 · The DFT and FFT by hand and Module 04: 2D DFT, windowing, k-space
- Module 11 · Correlation & time-delay estimation: delay estimation, NCC, Kasai/Loupas and kernel-length trade-offs
- Module 15 · Wave physics, nonlinear propagation and harmonic imaging: attenuation and acoustic intensity; Module 22 · Acoustic output, bioeffects and regulatory measurement: output limits
- Module 23 · Beam patterns, k-space and compounding: plane-wave compounding; Module 29 · Doppler and flow imaging: slow-time estimators
30.1Tissue mechanics: three moduli, two wave speeds and one conversion
A linear, isotropic, elastic solid is described by two constants. Elastography uses three names for them. The Young's modulus \(E\) is the ratio of stress to strain in uniaxial loading with free sides, the shear modulus \(\mu\) the ratio of shear stress to shear angle, and the bulk modulus \(K\) the ratio of pressure to volume change. Each wave sees a different combination:
- \(c_s\)
- shear (transverse) wave speed, m/s
- \(c_L\)
- longitudinal (ultrasound) speed, ≈ 1540 m/s
- \(\mu\)
- shear modulus, Pa (soft tissue 0.3–30 kPa)
- \(K\)
- bulk modulus, Pa (≈ 2.2–2.4 GPa, close to water)
- \(E\)
- Young's modulus, Pa
- \(\nu\)
- Poisson's ratio, dimensionless
- \(\rho\)
- mass density, ≈ 1000 kg/m³ (1050–1060 for liver and muscle)
How incompressible? Writing Poisson's ratio in terms of \(K\) and \(\mu\),
- \(\nu\)
- Poisson’s ratio, dimensionless
- \(K,\ \mu\)
- bulk and shear modulus, Pa, as in (30.1)
- \(\approx\)
- valid for \(\mu \ll K\), i.e. \(\nu \to 1/2\) (nearly incompressible tissue)
Three consequences follow. First, shear waves are 300–1500 times slower than ultrasound, so a plane-wave imager at 5–10 kHz can film them: in 100 µs a 2 m/s wave moves 0.2 mm while the ultrasound crosses 15 cm. Second, shear wavelengths are long: \(\lambda_s = c_s/f\) is 4 mm at 500 Hz for 2 m/s, 40 mm at 50 Hz, compared with 0.31 mm for 5 MHz ultrasound. Spatial resolution of any speed map is limited by \(\lambda_s\) and by the estimation kernel, not by the B-mode resolution. Third, the conversion \(E = 3\rho c_s^{2}\) quietly assumes linear elasticity, isotropy, incompressibility, an unbounded medium and no dispersion. Every one of those can fail in tissue, and Sections 30.7–30.8 quantify how. That is why the SRU and QIBA recommend reporting shear-wave speed in m/s as the measured quantity and treating kPa as a derived, model-dependent number.
| Method | Excitation | Measured | Output | Band |
|---|---|---|---|---|
| Strain (quasi-static) elastography | hand or physiological compression | axial displacement → strain | relative strain image, strain ratio | < 10 Hz |
| ARFI imaging | focused radiation-force push, 30–1000 µs | on-axis displacement, recovery | relative displacement image (µm) | impulsive |
| Point SWE (pSWE) | one push | shear-wave arrival times in a small ROI | one speed (m/s or kPa) | ≈ 50–500 Hz |
| 2D SWE (SSI, comb push) | pushes swept in depth (supersonic) or several beams | wavefield at 2–10 kHz frame rate | speed map | ≈ 50–800 Hz |
| Transient elastography | external piston, one 50 Hz cycle | A-line displacement vs depth and time | one stiffness (kPa) | 50 Hz centre |
| Harmonic / reverberant SWE, MRE | continuous vibration | steady-state wavefield phase | speed or complex modulus map | 20–500 Hz |
30.2Quasi-static strain imaging: displacement, its derivative and the strain filter
Ophir and colleagues (1991) compressed tissue by about 1 % with the probe, tracked the axial displacement between a pre- and a post-compression RF frame, and displayed its derivative. The tracking is the delay estimation of Module 11: an RF kernel of duration \(T\) is cross-correlated between frames and the lag of the peak is converted to displacement with the round-trip factor,
- \(\hat\tau(z)\)
- estimated echo delay change of the kernel centred at depth \(z\), s (negative when the tissue moves toward the probe)
- \(u\)
- axial displacement, m; one sample at 40 MHz is 19.25 µm
- \(\varepsilon\)
- axial strain, positive in compression
- \(z\)
- depth, m
- \(c\)
- longitudinal sound speed, m/s
- \(\partial\)
- partial derivative
Differentiation amplifies noise, so strain is never a two-point difference in practice. The standard estimator is a least-squares slope over \(N\) displacement estimates spaced \(\Delta\); for independent displacement errors of standard deviation \(\sigma_u\) its noise is
- \(N\)
- number of displacement estimates in the LSQ span
- \(\Delta\)
- spacing of the displacement estimates (kernel step), m
- \(\sigma_u\)
- displacement jitter, m (Walker–Trahey bound plus decorrelation, Module 11)
- \(\mathrm{SNR}_e\)
- elastographic SNR: mean strain over its standard deviation in a uniform region
- \(z_i,\ \hat u_i\)
- depth (m) and estimated displacement (m) of the \(i\)-th kernel
- \(\bar z,\ \bar u\)
- their means over the \(N\) kernels
- \(\hat\varepsilon,\ \bar\varepsilon\)
- least-squares strain estimate and its mean in the region
- \(\sigma_\varepsilon\)
- standard deviation of the strain estimate, dimensionless
The strain filter. \(\mathrm{SNR}_e\) is not monotonic in the applied strain. At small strain, \(\sigma_u\) is set by electronic noise and the Walker–Trahey floor, so \(\mathrm{SNR}_e\) grows in proportion to \(\varepsilon\). At large strain the kernel is compressed by \(\varepsilon T\) relative to its partner, the echoes decorrelate, \(\sigma_u\) grows and eventually peak-hopping errors appear. Plotting \(\mathrm{SNR}_e\) against \(\varepsilon\) gives a band-pass curve that Varghese and Ophir (1997) called the strain filter: every combination of kernel, bandwidth and SNR has a window of strains, typically 0.1–3 %, inside which the elastogram is usable. Between the two edges there is a ceiling set by speckle itself: a kernel's correlation peak tracks the amplitude-weighted centroid of its scatterers, which sits a random fraction of the kernel away from its centre, so under strain \(\varepsilon\) the displacement error is of order \(\varepsilon\) times that offset and \(\mathrm{SNR}_e\) stops growing. Temporal stretching of the post-compression frame by the mean strain (companding) moves the upper edge of that window up by realigning the echoes before correlation.
Normalization: why strain is not stiffness. The strain image is a stiffness image only if the stress is uniform. In a 1D layered medium under uniform stress the strain ratio equals the modulus ratio exactly. In 2D and 3D, stress concentrates in stiff inclusions and the strain contrast is compressed. For a circular inclusion in an incompressible matrix under plane strain (Eshelby's result for a cylinder),
- \(\varepsilon_{\mathrm{inc}},\ \varepsilon_\infty\)
- uniform strain inside the inclusion and the remote background strain
- \(r\)
- shear-modulus contrast of the inclusion
- \(\mu_{\mathrm{inc}},\ \mu_{\mathrm{bg}}\)
- shear modulus of the inclusion and of the background, Pa
Strain imaging: NCC tracking, LSQ strain, kernel size and the strain filter
Simulation1D: 2000 point scatterers in 40 mm, Gaussian 5 MHz pulse (60 % bandwidth), \(f_s\) = 40 MHz, white noise at the set SNR. A stiff layer (18–24 mm) of modulus ratio \(R\) under uniform stress (series springs) so the true strain ratio is \(R\). NCC with parabolic subsample fit and a guided search (±12 samples around the previous lag), LSQ strain over \(N\) estimates 0.25 mm apart; same model as m30_strain_imaging.m. No lateral or elevational decorrelation.30.3Acoustic radiation force and ARFI imaging
A wave that is absorbed gives up momentum as well as energy. The power absorbed per unit volume from a plane wave of intensity \(I\) is \(2\alpha I\) (intensity decays as \(e^{-2\alpha z}\) with \(\alpha\) the amplitude coefficient), and dividing by \(c\) converts power to momentum flux. The result is a body force along the beam (Nyborg; Torr 1984; the radiation-force balance of Module 22 weighs the same momentum):
- \(F\)
- body force per unit volume along the beam, N/m³
- \(\alpha\)
- amplitude absorption coefficient, Np/m (1 Np = 8.686 dB); usually approximated by the attenuation coefficient, which overestimates \(F\) by the scattering fraction
- \(I\)
- temporal-average intensity during the push, W/m² (1 W/cm² = 10⁴ W/m²)
- \(c\)
- longitudinal sound speed, m/s
The force is enormous per unit volume but acts on a column only about 0.8 mm wide and 5–10 mm long, for 30–1000 µs. What the tissue does with it depends on how the push duration \(T\) compares with the time a shear wave needs to leave the focal column, \(\sigma_r/c_s\) (≈ 0.2 ms for \(\sigma_r\) = 0.35 mm and 1.5 m/s). For short pushes the tissue cannot react elastically during the push and simply accelerates; for long pushes it approaches the static solution, in which the column is held by shear stresses:
- \(T\)
- push duration, s
- \(\sigma_r\)
- lateral standard deviation of the force (beam) profile, m
- \(F(\mathbf r')\)
- force distribution around the focus, N/m³
- \(u_{\mathrm{static}}\)
- steady displacement at the focus; this scalar kernel equals the incompressible Kelvin point-force solution \(F(1+\cos^{2}\theta)/(8\pi\mu r)\) only on the force axis (\(\theta\) = 0), so for a column of finite width it overestimates \(u_{\mathrm{static}}\) (by ≈ 20 % for the column of Demo 30.2)
- \(\rho\)
- mass density, kg/m³
- \(c_s,\ \mu\)
- as in (30.1)
- \(u(T)\)
- displacement at the focus after a push of duration \(T\), m
- \(F\)
- force density at the focus (in the short-push form), N/m³
- \(\mathbf{r}\prime,\ d^{3}r\prime\)
- integration point (m) and volume element (m³)
Momentum and heat come together. Without conduction, the temperature rise during one push is
- \(\rho C\)
- volumetric heat capacity, ≈ 3.7 MJ/(m³·K) for soft tissue
- \(\Delta T\)
- temperature rise, K
- \(\alpha,\ I,\ F,\ c\)
- as in (30.6)
- \(T\)
- push duration, s
The ARFI sequence (Nightingale et al. 2002) is: two or more reference tracking pulses on the push axis, the push (a long, focused, high-voltage burst from the imaging array), then tracking pulses at a PRF of 5–10 kHz for 3–10 ms, each correlated against the reference with the Kasai or Loupas estimator of Module 11. Displacements are 1–20 µm: 3–70 times the jitter floor of ≈ 0.3 µm (≈ 0.1 µm at high SNR), so the measurement itself is easy. The hard parts are elsewhere. The first 0.2–0.5 ms after the push are lost (transducer ring-down, reverberation and, for single-element tracking, the push and track beams share the aperture). The force is unknown in absolute terms, because \(\alpha\) and the in-situ \(I\) are unknown, so ARFI displacement maps are relative, and the depth dependence of \(I\) (focal gain, attenuation) must be normalized out. And the push excites shear waves, which is exactly what Section 30.4 uses.
Choosing the push frequency. For a fixed source intensity and \(\alpha = a f\), the force at depth \(z\) scales as \(f\,e^{-2afz}\), maximal at \(2afz = 1\): with \(a\) = 0.5 dB/cm/MHz (0.0576 Np/cm/MHz) at 4 cm, \(f^{*}\) = 2.2 MHz. A fixed aperture adds focal gain \(\propto f^{2}\), moving the optimum to \(3/(2az)\) = 6.5 MHz. Real systems push at 2.5–5 MHz, where transducer efficiency, heating of the lens and the MI limit also enter.
ARFI displacement: push duration, intensity, absorption and stiffness
SimulationAxisymmetric (r, z) finite-difference solution of the scalar Kelvin–Voigt shear-wave equation \(\rho\ddot u = \mu\nabla^{2}u + \eta\nabla^{2}\dot u + F\), JS, with absorbing layers; Gaussian force column (\(\sigma_r\) = 0.35 mm, \(\sigma_z\) = 2.5 mm). The step response is computed once per stiffness and viscosity; any push duration is the difference of two shifted step responses (linearity), and \(I\), \(\alpha\), \(f_0\) only scale the amplitude. Scalar model: it drops the \(\nabla(\nabla\cdot\mathbf{u})\) coupling of the elastic equations, so its static limit for this column is about 20 % above the incompressible Kelvin value and the off-axis field is approximate. Same model as m30_arfi_displacement.m.30.4Making shear waves: supersonic pushes, external vibrators and harmonic fields
A single push launches a shear wave that spreads cylindrically from the focal column (amplitude \(\propto r^{-1/2}\)) and spherically beyond its ends (\(\propto r^{-1}\)), so it is usable only over a few millimetres of depth and 5–10 mm laterally. Supersonic shear imaging (Bercoff, Tanter and Fink 2004) fires pushes at successive depths faster than the shear wave can leave them. The source moves along the beam axis at speed \(V\), the cylindrical wavelets interfere constructively on a Mach cone, and two quasi-plane shear waves leave it, tilted symmetrically about the lateral axis:
- \(V\)
- source (push) speed along the beam: focal step \(\Delta z\) per push duration \(T_p\), m/s
- \(M\)
- Mach number, must exceed 1
- \(\varphi\)
- half-angle of the cone; also the tilt of the wavefront normal from the lateral axis
- \(c_{x,\mathrm{app}}\)
- apparent speed measured by a purely lateral time-of-flight estimator
- \(c_s\)
- shear speed, as in (30.1)
An ultrafast plane-wave imager (Module 23) then films the field at 2–20 kHz. One frame of compounding costs only a few tilted transmits, so the whole wavefield \(u(x, z, t)\) is captured from one push sequence, and pushes on the left and right of the region are combined to fill the shadow behind each source. Variants trade the one Mach cone for several simultaneous beams (comb push; Song et al. 2012) or for pushes outside the imaging plane.
Transient elastography (Sandrin et al. 2003, FibroScan) does not use radiation force. A piston on the probe tip delivers one cycle at 50 Hz, a single-element transducer coaxial with it tracks displacement along one A-line at kHz rate, and the slope of the displacement front in the depth–time plane gives the speed, reported as \(E = 3\rho c_s^{2}\) in kPa. The 50 Hz excitation means \(\lambda_s\) = 20–60 mm, comparable to the measurement depth window (a few centimetres, 25–65 mm for the standard adult probe): the measurement happens in the near field of a finite source, where compressional and shear components mix and diffraction biases the apparent speed. Calibration against phantoms and fixed depth windows keep it reproducible; it is not a bulk-wave measurement in the textbook sense, and its kPa are not interchangeable with SWE kPa at 200–500 Hz in dispersive tissue.
- \(\Delta z/\Delta t\)
- slope of the shear front in the depth–time displacement image, m/s
- \(E_{\mathrm{TE}}\)
- Young’s modulus from transient elastography, Pa
- \(\rho\)
- mass density, kg/m³
- \(\lambda_s,\ c_s,\ f\)
- shear wavelength (mm), speed (m/s) and frequency (Hz)
Harmonic (continuous) methods drive the tissue at one or a few frequencies until a steady state forms, then measure the phase of the displacement at that frequency. MR elastography (Muthupillai et al. 1995) and ultrasound reverberant-field and crawling-wave methods are in this family. With a steady state the natural estimator is the local wavenumber \(k = \partial\varphi/\partial x\) or a Helmholtz inversion of the full field (Section 30.6), and the frequency is known by construction, which removes the phase–group ambiguity of Section 30.7. The price is standing waves and reflections, which is why directional filtering was first developed for MRE.
Supersonic push: Mach cone, wavefront tilt and the lateral TOF bias
Simulation2D finite-difference (leapfrog) solution of the scalar shear-wave equation in the imaging plane, JS: 121 × 141 cells of 0.2 mm, elastic, absorbing layers, Gaussian force columns (σx = 0.3 mm, σz = 1.2 mm) at successive depths 2 mm apart, each lasting 2 mm / V. The lateral TTP fit is taken on the row where the plane part of the right-hand front crosses x = 3–8 mm. Scalar model: displacement along z, propagation in (x, z), the force is a line source in elevation. The dashed lines are the Mach front predicted by (30.9). Same model as m30_supersonic_mach.m.30.5Speed estimation I: time of flight, time to peak and cross-correlation
The tracked data are a space–time matrix \(u(x_i, t_n)\) at one depth: \(N\) lateral positions spaced \(\Delta x\), sampled at the tracking PRF or plane-wave frame rate. The simplest estimators reduce each trace to one arrival time and fit a line. With the time to peak (TTP) of displacement or particle velocity at each position, and ordinary least squares,
- \(t_i\)
- arrival time (peak, or a correlation delay) at lateral position \(x_i\), s
- \(\sigma_t\)
- standard deviation of each arrival-time estimate, s
- \(N,\ \Delta x\)
- number and spacing of lateral positions in the fit; \(L = (N-1)\Delta x\) is the lateral kernel
- \(\hat c\)
- estimated speed, m/s (the inverse of the fitted slope)
- \(t_0\)
- origin time of the shear front, s
- \(\bar x,\ \bar t\)
- means of positions and arrival times
- \(c,\ \sigma_c\)
- shear speed and its standard deviation, m/s
- \(x_i\)
- lateral position of the \(i\)-th trace, m
TTP is cheap and robust to amplitude changes, but each \(t_i\) uses only the samples near the peak. Cross-correlation of whole traces uses all of them:
- \(v\)
- particle velocity (or displacement) trace, after any filtering
- \(\hat\tau_{ij}\)
- delay between traces at \(x_i\) and \(x_j\), with parabolic subsample interpolation, s
- \(x_i,\ x_j\)
- lateral positions of the two traces, m
- \(t_n,\ n\)
- time samples and their index
- \(\tau\)
- trial delay, s
- \(\arg\max_\tau\)
- delay that maximizes the cross-correlation
- \(\hat c\)
- estimated shear speed, m/s
Which speed? In an elastic, non-dispersive medium every estimator returns \(c_s\). In viscoelastic tissue the wave is dispersive (Section 30.7) and each estimator returns a different weighted average: TTP of a pulse follows the group velocity at the dominant frequency of the trace it is applied to, or more (the peak of particle velocity marks the steepest rise of the displacement, which the faster high-frequency content carries); cross-correlation sits closer to the phase velocity at that frequency; and the dominant frequency of particle velocity is higher than that of displacement (multiplying by \(j\omega\) tilts the spectrum up). So "TTP on displacement" (Palmeri et al. 2008; with RANSAC outlier rejection, Wang et al. 2010) and "TOF on particle velocity" (SSI) disagree in dispersive tissue by a few percent to over ten percent, for reasons that have nothing to do with noise. The QIBA phantom studies found exactly such system-to-system spreads in viscoelastic phantoms. Robust variants (Radon-sum TOF, Rouze et al. 2010; RANSAC line fits; bidirectional fits from left and right pushes) improve precision but do not change which speed is measured.
Practical constraints. The trace must be recorded long enough to contain the whole pulse at the farthest \(x\) (10 mm at 1 m/s is 10 ms), the frame rate must sample the pulse (10 kHz gives 20 samples per period at 500 Hz) and the region next to the push (≈ 1–2 mm) is excluded, because the source is not a line and the near field is not a travelling wave.
One push, four estimators: TTP, cross-correlation and k–f phase velocity on a dispersive wave
Approx. modelCylindrical shear wave from a 200 µs line push (σ = 0.35 mm) in a Kelvin–Voigt medium: far-field 2D Green's function \(U(r,\omega) \propto S(\omega)\,e^{-jk(\omega)r}/(\mu^{*}\sqrt{k r})\) with complex \(k = \omega\sqrt{\rho/\mu^{*}}\), evaluated on an FFT grid and inverted (JS); content below about 30 Hz removed. White noise at the set SNR (re the peak at 3 mm). Estimators run on the noisy traces from 3 mm to the window end. Same model as m30_shear_wave_tof_kf.m.30.6Speed estimation II: k–f phase velocity, local wavenumber, directional filters and maps
Time-domain estimators return one speed per kernel and blend frequencies. A 2D Fourier transform of the space–time data separates them. For a wave travelling in \(+x\), \(v(x,t) = \mathrm{Re}\int A(\omega)\,e^{j(\omega t - k(\omega)x)}\,d\omega\), so its energy lies on the curve \(k(\omega)\) in the \((k, f)\) plane, and reading the ridge row by row gives the phase velocity at every frequency (Bernal et al. 2011; Deffieux et al. 2009 for the related phase-gradient spectroscopy):
- \(V(k,f)\)
- 2D spectrum of the windowed space–time data; a \(+x\) wave peaks at \(k = +k(\omega)\) with this sign convention
- \(w(x)\)
- lateral window (Hann) of length \(L\), m
- \(k_{\mathrm{peak}}(f)\)
- wavenumber of the spectral maximum at frequency \(f\), rad/m, refined by peak interpolation
- \(\lambda\)
- shear wavelength at \(f\), m
- \(x,\ t\)
- lateral position (m) and time (s)
- \(v(x,t)\)
- particle velocity (or displacement) field
- \(k,\ f\)
- lateral wavenumber (rad/m) and temporal frequency (Hz)
- \(j\)
- imaginary unit
- \(\delta c,\ \delta k\)
- resolution in speed (m/s) and in wavenumber (rad/m)
- \(L\)
- window length, m
Local wavenumber. At one frequency \(f_0\) the temporal Fourier coefficient \(A(x; f_0)\) has a phase that advances by \(k\) per metre. Its spatial derivative is a local wavenumber, and a local phase velocity, with a kernel as small as two points:
- \(A(x; f_0)\)
- temporal DFT of the trace at \(x\), evaluated at \(f_0\)
- \(\Delta\)
- half the lateral kernel, m; requires \(2k\Delta < \pi\), i.e. \(4\Delta < \lambda\), or the phase wraps
- \(\hat k(x),\ \hat c(x)\)
- local wavenumber (rad/m) and speed (m/s) at position \(x\)
- \(f_0\)
- analysis frequency, Hz
- \((\cdot)^{*}\)
- complex conjugate
- \(\arg\)
- phase angle, rad
Directional filtering. Every estimator above assumes one wave travelling in one direction. Reflections from an inclusion boundary or an organ capsule, the wave from the opposite push, and guided modes violate that, and TOF estimators see a superposition whose apparent arrival time is biased. In the \((k, f)\) plane, waves travelling in \(+x\) and \(-x\) occupy opposite quadrant pairs, so a mask separates them (Manduca et al. 2003 for MRE; Deffieux et al. 2011 for transient SWE):
- \(H\)
- directional mask that keeps \(+x\) waves (with the sign convention of (30.13)); smooth its edges in practice to limit ringing
- \(v_{+}\)
- the right-going part of the field
- \(V,\ k,\ f\)
- as in (30.13)
- \(\mathcal{F}^{-1}\)
- inverse 2D Fourier transform back to \((x, t)\)
From traces to a map. A 2D SWE map evaluates one of these estimators in a sliding kernel at every pixel: TTP or cross-correlation between \(x - \Delta\) and \(x + \Delta\), or a local phase gradient. The kernel sets resolution (≈ \(2\Delta\) laterally, plus the depth averaging) and precision ((30.11)), the push geometry sets coverage (each push shadows itself, so maps fuse left and right pushes), and the directional filter decides whether boundaries produce artefacts. Demo 30.5 runs that chain on a simulated inclusion.
A 2D SWE map: inclusion, reflections, directional filter and three local estimators
Simulation2D finite-difference scalar shear-wave equation, JS: 121 × 81 cells of 0.25 mm (x 0–30 mm, z 10–30 mm), elastic, harmonic-mean interface moduli, absorbing layers; background 1.5 m/s, circular inclusion of radius 4 mm centred at (17, 20) mm; line push at x = 5 mm for 150 µs; 10 kHz frames for 20 ms; white noise added to particle velocity. Directional filter: 2D FFT per depth row (x × t, zero-padded to 128 × 256), quadrant mask (30.15). Maps: TTP slope over ±Δ, cross-correlation between x ± Δ, or phase gradient (30.14) at 250 Hz.30.7Viscoelasticity and dispersion: phase velocity, group velocity and the model behind the kPa
Tissue stores and dissipates energy. In the frequency domain this replaces \(\mu\) by a complex modulus \(\mu^{*}(\omega)\), the wavenumber becomes complex, and speed and attenuation become functions of frequency. The simplest model puts a spring \(\mu\) in parallel with a dashpot \(\eta\) (Kelvin–Voigt):
- \(\mu^{*}\)
- complex shear modulus, Pa
- \(\eta\)
- shear viscosity, Pa·s
- \(c_{\mathrm{ph}} = \omega/\mathrm{Re}\,k\)
- phase velocity, m/s
- \(\alpha = -\mathrm{Im}\,k\)
- shear-wave attenuation, Np/m
- \(\omega\)
- angular frequency, rad/s
- \(\mu,\ \rho\)
- shear (storage) modulus, Pa, and density, kg/m³
- \(j\)
- imaginary unit
- \(k\)
- complex wavenumber, rad/m
A pulse contains many frequencies, each travelling at its own phase velocity. The envelope travels at the group velocity
- \(c_g\)
- group velocity, m/s
- \(\omega\)
- angular frequency, rad/s
- \(c_{\mathrm{ph}}\)
- phase velocity, m/s, as in (30.16)
- \(\mathrm{Re}\,k\)
- real part of the complex wavenumber, rad/m
Tissue data are better described by a power law (a spring-pot, or fractional derivative model; KVFD adds a spring in parallel). With \(\mu^{*} = \mu_0 (j\omega/\omega_0)^{n}\), \(0 \le n \le 1\),
- \(n\)
- power-law exponent: 0 = elastic, 1 = Newtonian fluid
- \(c_0\)
- \(\sqrt{\mu_0/\rho}\) at the reference frequency \(\omega_0\), m/s
- \(\omega\)
- angular frequency, rad/s
- \(c_{\mathrm{ph}},\ \alpha,\ c_g\)
- phase velocity, shear-wave attenuation and group velocity, as in (30.16), (30.17)
- \(\tan,\ \cos\)
- argument in rad
Three practical consequences. Vendors disagree partly because they push at different frequencies, track at different frame rates and apply different estimators, so they sample different points of the same dispersion curve. Dispersion is a biomarker: shear-wave dispersion slope (m/s per kHz), attenuation and the power-law exponent are reported as viscosity surrogates for inflammation and steatosis, and Demo 30.6 computes all three. kPa hides the model: \(E = 3\rho c^{2}\) applied to a group velocity in a dispersive medium is neither the storage modulus nor the loss modulus at any frequency.
Dispersion: Kelvin–Voigt versus power law, phase versus group velocity, and a pulse that changes shape
Exact formulaCurves from (30.16)–(30.18); the group velocity by numerical differentiation of Re k. Pulses: 1D plane-wave propagation \(U(x,\omega) = U(0,\omega)\,e^{-jk(\omega)x}\) of a Gaussian-modulated pulse (80 % bandwidth), FFT at 20 kHz (exact for the plane wave in that medium, no geometric spreading). "Apparent TTP speed" = 10 mm divided by the difference in time to peak of the displacement envelope between 5 and 15 mm.30.8Guided waves, anisotropy, and measuring stiffness in a clinic
Guided waves. When the tissue layer is thinner than a shear wavelength (arterial wall 0.5–1.5 mm, cornea about 0.55 mm, myocardium, skin, tendon, bladder wall), the wave is not a bulk wave. It is a Lamb mode of a plate (or a tube), and the fundamental antisymmetric (flexural) mode A0 that a push excites is strongly dispersive even in a perfectly elastic plate. For an incompressible plate of thickness \(h\) in vacuum, with \(k = \omega/c\) and \(q' = \omega\sqrt{1/c^{2} - 1/c_t^{2}}\) (subsonic, \(c < c_t\)), the Rayleigh–Lamb antisymmetric equation reduces to
- \(c\)
- A0 phase velocity at angular frequency \(\omega\), m/s
- \(c_t\)
- bulk shear speed of the plate material, m/s
- \(h\)
- plate thickness, m
- \(k\)
- wavenumber \(\omega/c\), rad/m
- \(q\prime\)
- transverse decay constant \(\omega\sqrt{1/c^{2} - 1/c_t^{2}}\), rad/m
- \(\omega\)
- angular frequency, rad/s
- \(\tanh\)
- hyperbolic tangent
- \(c_{A0}\)
- phase velocity of the A0 (flexural) mode, m/s
- \(\omega\)
- angular frequency, rad/s
- \(h,\ c_t\)
- plate thickness (m) and bulk shear speed (m/s), as in (30.19)
So in an artery, the cornea or a thin capsule, the speed depends on thickness and frequency as much as on stiffness, and the inversion must fit a dispersion curve (k–f data against (30.19) or a tube model), not apply (30.1). Fluid loading lowers the curve further, and in a tube the wall pressure (pre-stress) changes \(c_t\) itself through acoustoelasticity, which is why carotid SWE values swing over the cardiac cycle.
Anisotropy. Skeletal muscle, myocardium, tendon and kidney cortex are fibrous. Modelled as transversely isotropic and incompressible, the horizontally polarized (SH) shear wave travelling at angle \(\theta\) to the fibres has
- \(\mu_L\)
- shear modulus in planes containing the fibre axis (wave along the fibres), Pa
- \(\mu_T\)
- shear modulus in the plane transverse to the fibres, Pa
- \(\theta\)
- angle between propagation direction and fibres
- \(c_{SH}(\theta)\)
- shear-wave speed at angle \(\theta\), m/s
- \(\rho\)
- mass density, kg/m³
Clinical measurement. The liver is the main application, and its protocol (SRU consensus update, Barr et al. 2020) encodes everything in this module: measure 1.5–2.0 cm below the capsule (reverberation and the capsule's guided-wave behaviour above, attenuation below), with the push perpendicular to the capsule, in neutral breath-hold after at least 4 hours of fasting (post-prandial portal flow raises stiffness), 10 measurements for point SWE (5 acceptable for 2D SWE with a quality map), and accept the median only if IQR/median ≤ 30 % in kPa (≤ 15 % in m/s: the factor 2 is the quadratic conversion, \(\delta E/E = 2\,\delta c/c\)). The vendor-neutral "rule of four" for ARFI-based systems:
| Liver stiffness | E (kPa) | cs (m/s), E = 3ρc² | Interpretation (SRU 2020) |
|---|---|---|---|
| normal | ≤ 5 | ≤ 1.3 (5.07 kPa) | high probability of being normal |
| rules out cACLD | < 9 | < 1.7 (8.67 kPa) | without other clinical signs |
| grey zone | 9–13 | 1.7–2.1 | suggestive of cACLD, may need confirmation |
| cACLD | > 13 | > 2.1 (13.2 kPa) | highly suggestive of compensated advanced chronic liver disease |
| CSPH | > 17 | > 2.4 (17.3 kPa) | probability of clinically significant portal hypertension |
Reproducibility and QIBA. The RSNA Quantitative Imaging Biomarkers Alliance measured the same elastic and viscoelastic phantoms on many commercial systems. The lessons (Palmeri et al. 2021): the 95 % interval of between-system speed differences was ±9.6 % in elastic phantoms and ±15.3 % in viscoelastic ones (different push and tracking bands sample different points of the dispersion curve, Section 30.7); system and site mattered more than focal depth; and digital phantoms (simulated data run through each vendor's estimator) separate estimator bias from acquisition bias. That is why the profile asks for m/s, a stated depth, and dispersion-aware reporting. For a pipeline engineer the corollary is an error budget: estimator bias in dispersive media, guided-wave and boundary effects, reflection bias without directional filtering, noise ((30.11)), and the conversion itself.
Guided waves in a thin layer: the A0 Lamb mode and the stiffness you would wrongly report
Approx. modelBlue: exact numerical root of the incompressible free-plate Rayleigh–Lamb antisymmetric equation (30.19) (bisection, JS). Pink dashed: Kirchhoff thin-plate limit. Gold dashed: immersed-plate empirical fit of Couade et al. 2010, K = 0.5 (30.20), an approximation (no tube curvature, no pre-stress, no viscosity). Same root finder as m30_lamb_a0.m.30.9Visual summary
Twelve questions an elastography expert answers in five seconds. Each figure is computed from the formula or model named in its badge.
30.10How experts think about it
An elastography expert separates three layers and never lets them blur: the excitation (what force, where, for how long, at what frequency content), the motion estimate (a delay estimator with known bias and jitter, Module 11), and the inversion (which wave model turns arrival times or wavenumbers into a speed, and which material model turns the speed into kPa). Most published disagreements trace to the third layer: dispersion, guided modes, reflections and anisotropy change the speed that a given estimator reports, while the noise budget is usually the smaller problem.
\(c_L = \sqrt{(K + 4\mu/3)/\rho}\) barely notices \(\mu\); \(c_s = \sqrt{\mu/\rho}\) is all \(\mu\) (30.1). That is why B-mode contrast says almost nothing about stiffness, and why a plane-wave imager at a few kHz can film shear waves that move 1540/1.7 ≈ 900 times slower than its own pulses.
Incompressibility gives \(E = 2\mu(1+\nu) = 3\mu\) to six digits (30.2). The square doubles relative errors: a 5 % speed bias is a 10 % stiffness bias, which is why the SRU quality criterion is 30 % in kPa but 15 % in m/s.
In a dispersive medium \(3\rho c^{2}\) depends on which speed you put in. KV with \(c_0\) = 1.5 m/s and η = 1 Pa·s gives \(c_{\mathrm{ph}}\) = 1.66 and \(c_g\) = 2.01 m/s at 200 Hz (30.16)–(30.17): 8.3 versus 12.1 kPa from the same tissue at the same frequency, against 6.75 kPa of elastic modulus.
Stress concentrates in stiff inclusions and avoids soft ones; for a circular inclusion in incompressible plane strain the inside strain is \(2\varepsilon_\infty/(1+r)\) (30.5). Only layers in series reproduce the modulus ratio. Strain ratios are therefore monotonic but compressed and geometry-dependent.
At small strain electronic jitter is fixed, so SNRe ∝ ε; in the middle a speckle term (each kernel tracks its amplitude-weighted scatterer centroid, an error proportional to ε) caps SNRe; at large strain the kernel stretch \(\varepsilon T\) decorrelates the echoes and SNRe collapses (the strain filter). In the 1D simulation of m30_strain_imaging.m: SNRe 3.1, 4.7, 5.2, 3.1 and 0.4 at 0.05, 0.3, 1, 3 and 6 %; companding moves the upper edge out.
Absorbed power per volume is \(2\alpha I\); dividing by \(c\) gives momentum per volume per second (30.6). The force is linear in absorption and intensity, so ARFI images also map attenuation, and the force is unknown in absolute terms in vivo.
For \(T \ll \sigma_r/c_s\) the column accelerates freely (30.7). In the axisymmetric model, quadrupling μ lowers the peak from 3.51 to 1.98 µm for a 200 µs push and from 14.1 to 5.0 µm for 1 ms, against a factor 4 in the static limit. ARFI displacement is a compressed, push-dependent proxy of stiffness.
With ρC ≈ 3.7 MJ/(m³·K), 500 W/cm² for 200 µs at 23 Np/m gives 12.4 mK (30.8). Fifty pushes delivered faster than conduction removes heat (focal time constants of order seconds) add up to 0.6 K at the focus, and the transducer face heats too. That caps the 2D SWE refresh rate near 1 Hz.
The fronts leave the Mach cone along its normal (30.9). In the 2D FD simulation with \(c_s\) = 2 m/s, the lateral TTP fit reads 1.99, 2.19 and 2.80 m/s at M = 7.5, 2.5 and 1.5, against \(c_s/\cos\varphi\) = 2.02, 2.18 and 2.68. High Mach numbers or direction-aware fits remove the bias.
An OLS slope over \(N\) points spaced Δx has variance \(12\sigma_t^{2}/(\Delta x^{2}N(N^{2}-1))\) (30.11). Halving L at fixed spacing multiplies the noise by 2.8. With correlated errors between traces the absolute number is optimistic, the scaling is not.
On one simulated push (c0 = 1.5 m/s, η = 0.4 Pa·s, 30 dB) the three estimators spread by 21 % with standard deviations of 0.05, 0.007 and 0.014 m/s (m30_shear_wave_tof_kf.m). TTP on particle velocity sits above \(c_g\) at the trace centroid; XC tracks \(c_{\mathrm{ph}}\) there; k–f gives the curve.
Differentiation multiplies the spectrum by ω. The velocity trace of Demo 30.4 centres at 237 Hz, the displacement trace near 105 Hz, and a positively dispersive medium returns a higher speed for the higher band: 1.86 versus 1.62 m/s by TTP at η = 0.4 Pa·s.
Fourier resolution is \(2\pi/L\) (30.13); the phase difference \(2k\Delta\) must stay below π (30.14). The two estimators bracket the resolution–spectral-detail trade-off: k–f for dispersion and guided waves, local phase for maps.
Reflections superpose a counter-propagating wave on the incident one; single-frequency phase estimators are hit hardest. In Demo 30.5 the phase-gradient estimate inside a 3 m/s inclusion goes from 2.51 ± 0.99 to 3.02 ± 0.29 m/s with the quadrant mask (30.15); whole-pulse cross-correlation was already robust.
KV's dashpot makes the modulus grow linearly with ω, so speed grows as √ω at high frequency (30.16). A spring-pot gives \(c \propto \omega^{n/2}\) with attenuation tied to \(n\) (30.18): n = 0.2 means ×1.26 per decade. Fit KV only over a narrow band and quote the band.
The flexural A0 mode is dispersive even in an elastic plate (30.19)–(30.20). For h = 1 mm and \(c_t\) = 3 m/s the exact root gives 1.37 m/s at 200 Hz and 2.68 m/s at 2 kHz. Arteries, cornea and capsules need a dispersion fit, ideally on the high-frequency band.
For the SH wave in a transversely isotropic medium (30.21). Gennisson et al. (2010) reported Voigt elasticities of 5.9 kPa along and 1.6 kPa across the brachialis fibres at rest, and 101 versus about 1 kPa under a 4 kg load. A probe rotated by 20° away from the fibres already loses about 4 % of the along-fibre speed with these values.
\(\lambda_s = c_s/f\) (30.10): at 50 Hz, 20–60 mm for 1–3 m/s. Diffraction of the piston and the coupled compressional field bias the apparent speed; fixed depth windows and calibration make it reproducible, and its kPa are not interchangeable with SWE kPa measured at 200–500 Hz in dispersive liver.
30.11Cheat sheet
| Concept | Formula | Ultrasound use |
|---|---|---|
| Shear speed and modulus | \(\mu = \rho c_s^{2}\), \(E \approx 3\rho c_s^{2}\) | 2 m/s = 4 kPa (μ) = 12 kPa (E) |
| Incompressibility | \(\tfrac12 - \nu \approx \mu/(2K)\) | 6×10⁻⁷ for soft tissue; never ν = 0.49 in FEM |
| Displacement from delay | \(u = c\tau/2\) | 19.25 µm per sample at 40 MHz |
| LSQ strain noise | \(\sigma_\varepsilon = (\sigma_u/\Delta)\sqrt{12/(N(N^{2}-1))}\) | kernel step and span design |
| Strain contrast of a cylinder | \(\varepsilon_\infty/\varepsilon_{\mathrm{inc}} = (1+r)/2\) | strain ratio under-reads modulus ratio |
| Radiation force | \(F = 2\alpha I/c\) | 1.5×10⁵ N/m³ at 500 W/cm², 23 Np/m |
| Inertial ARFI displacement | \(u(T) = FT^{2}/(2\rho)\) | 3.0 µm after a 200 µs push |
| Heating per push | \(\Delta T = 2\alpha IT/(\rho C)\) | 12 mK per 200 µs push |
| Mach cone | \(\sin\varphi = c_s/V\), lateral bias \(1/\cos\varphi\) | SSI push sequence design |
| TOF regression precision | \(\sigma_c/c \approx \sqrt{12}\,c\,\sigma_t/(L\sqrt N)\) | kernel length versus precision |
| k–f phase velocity | \(c_{\mathrm{ph}}(f) = 2\pi f/k_{\mathrm{peak}}(f)\), \(\delta k/k = \lambda/L\) | dispersion, guided waves |
| Local wavenumber | \(\hat k = -\arg[A(x+\Delta)A^{*}(x-\Delta)]/(2\Delta)\), \(\Delta < \lambda/4\) | local phase-velocity maps |
| Directional filter | keep \(kf > 0\) quadrants | remove reflections before TOF |
| Kelvin–Voigt | \(\mu^{*} = \mu + j\omega\eta\), \(c_{\mathrm{ph}} \approx c_0(1 + 3x^{2}/8)\) | viscosity fits over a narrow band |
| Group velocity | \(c_g = c_{\mathrm{ph}}/(1 - (\omega/c_{\mathrm{ph}})\,dc_{\mathrm{ph}}/d\omega)\) | what TOF reports in dispersive tissue |
| Power law | \(c_{\mathrm{ph}} \propto \omega^{n/2}\), \(\alpha = (\omega/c_{\mathrm{ph}})\tan(n\pi/4)\), \(c_g = c_{\mathrm{ph}}/(1-n/2)\) | dispersion and attenuation biomarkers |
| A0 Lamb mode, thin plate | \(c \approx \sqrt{\omega h c_t/\sqrt3}\) (free), \(/(2\sqrt3)\) (immersed) | artery, cornea, capsule |
| Transverse isotropy (SH) | \(\rho c^{2} = \mu_T\sin^{2}\theta + \mu_L\cos^{2}\theta\) | muscle, myocardium fibre direction |
| Liver rule of four (SRU 2020) | 5 / 9 / 13 / 17 kPa = 1.3 / 1.7 / 2.1 / 2.4 m/s | vendor-neutral ARFI cut-offs |
30.12Interview and qualifying-exam questions
I1Why can a scanner whose pulses travel at 1540 m/s measure tissue stiffness at all? · screen
Because it does not measure stiffness with its own pulses: it uses them as a motion sensor. Compressional speed depends on \(K + 4\mu/3\) with \(K/\mu \sim 10^{6}\), so it is insensitive to \(\mu\). Elastography deforms the tissue (compression, radiation force, vibration) and tracks the displacement with speckle tracking at kHz rates; the shear wave moves at 1–10 m/s, about 1000 times slower, so each frame freezes it. Follow-up: what frame rate? Several samples per shear period: 5–10 kHz for 200–800 Hz content.
I2Convert 1.8 m/s to kPa and list every assumption you just made. · screen
\(E = 3\rho c^{2}\) = 3 × 1000 × 3.24 = 9.7 kPa (μ = 3.24 kPa). Assumptions: linear elastic (no dispersion, so phase = group), isotropic, incompressible (ν = ½), unbounded (no guided modes), ρ = 1000 kg/m³, and a plane shear wave measured along its propagation direction. Follow-up: which fails first in liver? Dispersion (viscosity), then near-capsule effects; in muscle, anisotropy; in arteries, guided waves.
I3Derive the radiation force and estimate the displacement of a 200 µs push. · onsite
Power absorbed per volume \(2\alpha I\), divided by \(c\) for momentum flux: \(F = 2\alpha I/c\). At 4 MHz, 0.5 dB/cm/MHz (23 Np/m), 500 W/cm²: 1.5×10⁵ N/m³. During a short push the column accelerates at \(F/\rho\) = 150 m/s², so \(u(T) = FT^{2}/2\rho\) = 3.0 µm; it coasts a little after (3.5 µm peak in the axisymmetric model at 1.5 m/s). Follow-up: heating, 2αIT/ρC = 12 mK per push; and why the in-situ force is unknown (α, derating, focal gain).
I4Your strain images look great at 1 % compression and terrible at 5 %. Why, and what do you change? · onsite
The strain filter: at 5 % the kernel is compressed by εT relative to its partner (0.1 mm over a 2 mm kernel, a third of a wavelength at 5 MHz), the echoes decorrelate and peak hopping starts. Fixes: global or adaptive companding (stretch the post frame by the expected strain), shorter kernels, multi-step (accumulate small inter-frame strains), guided search and outlier rejection. Follow-up: and at 0.05 %? Then jitter dominates; use longer kernels, larger LSQ span or more compression.
I5Explain supersonic shear imaging. Why does the Mach number matter for estimation? · onsite
Pushes are fired at successive depths faster than the shear wave can leave, so the source moves at \(V > c_s\) and the wavelets add on a Mach cone, giving two quasi-plane waves of larger amplitude and depth extent than one push. The fronts are tilted by \(\varphi = \arcsin(1/M)\); a lateral-only TOF reads \(c_s/\cos\varphi\): 0.5 % high at M = 10, 15.5 % at M = 2. Follow-up: how to fix? Higher M, 2D direction estimation, or directional filters.
I6TTP regression, cross-correlation or k–f: which would you put in a product, and why? · onsite
Depends on the claim. For a robust single number in liver (pSWE), TTP or cross-correlation regression with outlier rejection (RANSAC, Radon sum): cheap, short kernels, stable at low SNR; XC is more precise than TTP because it uses the whole trace. For maps, local XC between x ± Δ or local phase gradient with directional filtering. For viscosity or thin structures, k–f phase velocity, because it gives \(c(f)\) and lets you fit a model. And whichever it is, specify the frequency band, because in dispersive tissue they disagree (1.85 / 1.55 / 1.53 m/s in m30_shear_wave_tof_kf.m). Follow-up: how do you validate? Viscoelastic phantoms plus digital phantoms through your exact pipeline.
I7Derive the Kelvin–Voigt phase velocity and explain why group velocity exceeds it. · onsite
\(\mu^{*} = \mu + j\omega\eta\), \(k = \omega\sqrt{\rho/\mu^{*}}\); write \(\mu^{*} = \lvert\mu^{*}\rvert e^{j\theta}\), Re k = \(\omega\sqrt{\rho/\lvert\mu^{*}\rvert}\cos(\theta/2)\), and \(\cos(\theta/2) = \sqrt{(\lvert\mu^{*}\rvert+\mu)/2\lvert\mu^{*}\rvert}\) gives (30.16). Since \(c_{\mathrm{ph}}\) rises with ω, \(d\omega/dk > \omega/k\): \(c_g = c_{\mathrm{ph}}/(1 - (\omega/c)\,dc/d\omega)\). Numbers: 1.5 m/s, 1 Pa·s, 200 Hz: 1.66 and 2.01 m/s. Follow-up: why do people prefer power laws? KV disperses as √ω at high frequency, much faster than tissue.
I8Your 2D SWE map shows a bright halo and stripes upstream of a stiff lesion. Diagnose. · research
Reflections from the lesion boundary: the reflected wave travels back through the background, superposes on the incident one and biases arrival times or phases (stripes at half-wavelength spacing upstream), while the halo is the kernel straddling the boundary (partial volume over 2Δ). Fix with directional filtering in k–f (or angular sectors in 2D), smaller kernels near edges, and compounding of left and right pushes; check the residual with a simulated inclusion run through the exact pipeline. Follow-up: cost of the filter? Edge smearing over about a wavelength and distortion of near-field components.
I9How would you measure the stiffness of the carotid wall with SWE? · research
Not with \(3\rho c^{2}\): the 0.6–1 mm wall carries guided (Lamb, tube) modes. Push in the wall, track along it at high frame rate, compute the k–f spectrum, extract \(c_{\mathrm{ph}}(f)\), and fit a plate or tube model: the immersed-plate fit \(c^{2} = \omega h c_t/(2\sqrt3)\) (Couade et al. 2010) in the high-frequency part, or the full leaky-Lamb or tube dispersion. Measure thickness from B-mode, synchronize with the ECG because pressure changes \(c_t\), and report the band. Follow-up: error from thickness? \(c_t \propto c^{2}/h\), so a 10 % thickness error is a 10 % \(c_t\) and 21 % \(E\) error.
I10Two vendors measure the same patient at 1.45 and 1.62 m/s. Is one of them wrong? · research
Not necessarily. 12 % is the size of the dispersion effects in liver: different push durations and tracking bands sample different frequencies, TTP on displacement versus TOF on particle velocity weights different spectra, ROI depth and focal configuration differ, and kPa conversion doubles any gap. The QIBA phantom work found exactly such system dependence growing in viscoelastic phantoms. The right comparison is c(f) over a stated band, ideally from a k–f or spectroscopy method, plus a phantom cross-calibration. Follow-up: what would you put in the 510(k) bench testing? Elastic and viscoelastic phantoms over depth, with the estimator's frequency response stated.
I11Where do elastography, sound-speed estimation and aberration correction meet? · research
All three are delay estimation problems on the same channel data. Sound-speed errors defocus the tracking beams and bias the displacement estimates (and the depth scale, so the measured shear speed scales with the assumed c). Aberration lowers tracking SNR and correlation, raising jitter. Conversely, a measured SoS map (Module 32) and an aberration correction (Module 33) improve tracking, and steatosis raises attenuation and lowers SoS while it changes shear dispersion, so multiparametric QUS + SWE models share one forward model. Follow-up: how much does a 2 % c error change the reported shear speed? Lateral distances are set by the array geometry and frame times by the clock, so to first order little; depth registration of the ROI shifts by 2 %.
30.13MATLAB lab
Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run.
m30_strain_imaging.m1D RF pre/post compression, NCC tracking, LSQ strain, stiff layer, strain filter with and without companding · change W, N, SNR, strains · runtime ≈ 3 s
%% m30_strain_imaging.m
% Purpose : 1D quasi-static strain imaging. Pre- and post-compression RF
% lines are synthesized from 2000 point scatterers (5 MHz Gaussian
% pulse, 60 % bandwidth, fs = 40 MHz). A stiff layer (18-24 mm,
% modulus ratio R) sits in a uniform-stress (series spring) medium.
% Displacement by normalized cross-correlation with a guided search
% and parabolic subsample fit; strain by least squares over N
% estimates. Then the strain filter: SNRe versus applied strain in a
% homogeneous medium, with and without companding (stretching the
% post-compression line by the mean strain).
% Change : epsApp, R, W (kernel, m), N (LSQ span), snrdB, the strain list.
% Expect : at 1 % strain, 2 mm kernel, N = 9, 30 dB: background strain
% 1.111 %, layer 0.370 %, measured strain ratio close to 3, SNRe of
% order 5 (speckle-limited); companding raises SNRe at large strain.
% Eshelby cylinder: strain ratio (1+r)/2 = 2.5 and 5.5 for r = 4, 10.
rng(1);
%% Parameters
c = 1540; fs = 40e6; f0 = 5e6; B = 0.6; Z = 0.04; nsc = 2000;
dzs = c/(2*fs); ns = round(Z/dzs); % 19.25 um per sample
sf = B*f0/2/sqrt(2*log(2)); st = 1/(2*pi*sf); % pulse envelope std (s)
half = ceil(4*st*fs);
zsc = rand(nsc,1)*Z; amp = randn(nsc,1);
n1 = randn(ns,1); n2 = randn(ns,1);
z1 = 0.018; z2 = 0.024; w = z2 - z1;
% RF synthesis (vectorized with accumarray)
%% Run one configuration
W = 2e-3; N = 9; snrdB = 30; epsApp = 0.01; R = 3;
ebg = epsApp*Z/(Z - w + w/R); eL = ebg/R;
ufun = @(z) -( (z<z1).*ebg.*z + (z>=z1 & z<z2).*(ebg*z1 + eL*(z-z1)) + (z>=z2).*(ebg*z1 + eL*w + ebg*(z-z2)) );
fprintf('True strain: background %.3f %%, layer %.3f %% (ratio %.2f)\n', ebg*100, eL*100, ebg/eL);
fprintf('Max displacement %.0f um = %.1f samples\n', -ufun(Z)*1e6, -ufun(Z)/dzs);
idx0 = (-half:half);
synth = @(zk) accumarray(reshape(min(max(round(2*zk/c*fs) + idx0, 1), ns), [], 1), ...
reshape(repmat(amp, 1, numel(idx0)) .* exp(-((min(max(round(2*zk/c*fs) + idx0, 1), ns) - 1)/fs - 2*zk/c).^2/(2*st^2)) ...
.* cos(2*pi*f0*((min(max(round(2*zk/c*fs) + idx0, 1), ns) - 1)/fs - 2*zk/c)), [], 1), [ns 1]);
addn = @(y, nz, snr) y + sqrt(mean(y.^2)/10^(snr/10))*nz;
x = addn(synth(zsc), n1, snrdB);
y = addn(synth(zsc + ufun(zsc)), n2, snrdB);
nw = round(W/dzs); step = 13; S = 12;
zc = []; uu = []; prev = 0;
for c0 = 41:step:(ns - nw - 41)
xs = x(c0:c0+nw-1); ex = sum(xs.^2); lags = (prev-S):(prev+S); cc = -2*ones(size(lags));
for q = 1:numel(lags)
i0 = c0 + lags(q); if i0 < 1 || i0 + nw - 1 > ns, continue; end
ys = y(i0:i0+nw-1); cc(q) = sum(xs.*ys)/sqrt(ex*sum(ys.^2));
end
[~, q] = max(cc); d = 0;
if q > 1 && q < numel(lags), den = cc(q-1) - 2*cc(q) + cc(q+1); if den < 0, d = 0.5*(cc(q-1) - cc(q+1))/den; end, end
prev = lags(q); zc(end+1) = (c0 - 1 + nw/2)*dzs; uu(end+1) = (lags(q) + d)*dzs; %#ok<AGROW>
end
ne = numel(zc) - N + 1; ze = zeros(1,ne); ee = zeros(1,ne);
for i = 1:ne
zz = zc(i:i+N-1); u2 = uu(i:i+N-1); zm = mean(zz);
ee(i) = -sum((zz - zm).*(u2 - mean(u2)))/sum((zz - zm).^2); ze(i) = zm;
end
bg = ze > 0.004 & ze < 0.016; bg2 = ze > 0.028 & ze < 0.036; ly = ze > 0.0195 & ze < 0.0225;
SNRe = mean(ee(bg))/std(ee(bg));
fprintf('Kernel %.1f mm, N = %d, %d dB: SNRe (4-16 mm) = %.1f, strain std = %.3f %%\n', W*1e3, N, snrdB, SNRe, std(ee(bg))*100);
fprintf('Measured strain ratio = %.2f (true %.1f)\n', mean([ee(bg) ee(bg2)])/mean(ee(ly)), R);
fprintf('Eshelby cylinder strain ratio (1+r)/2: r = 4 -> %.1f, r = 10 -> %.1f\n', (1+4)/2, (1+10)/2);
figure(1);
subplot(2,1,1); zt = linspace(0, Z, 400);
plot(zt*1e3, ufun(zt)*1e6, 'k--', zc*1e3, uu*1e6, 'b'); ylabel('u (\mum)'); title('Displacement: true (dashed) and NCC estimate');
subplot(2,1,2);
plot([0 18 18 24 24 40], [ebg ebg eL eL ebg ebg]*100, 'k--', ze*1e3, ee*100, 'b'); xlabel('depth (mm)'); ylabel('strain (%)'); title('LSQ strain');
%% Strain filter (homogeneous medium)
epsList = [0.0005 0.001 0.003 0.01 0.03 0.06 0.1];
SN = zeros(2, numel(epsList));
for ie = 1:numel(epsList)
e = epsList(ie);
yh = addn(synth(zsc*(1 - e)), n2, snrdB);
for comp = 0:1
yy = yh; sh = 0;
if comp, sh = e; yy = interp1((0:ns-1)', yh, (0:ns-1)'*(1 - sh), 'linear', 0); end
zc = []; uu = []; prev = 0;
for c0 = 41:step:(ns - nw - 41)
xs = x(c0:c0+nw-1); ex = sum(xs.^2); lags = (prev-S):(prev+S); cc = -2*ones(size(lags));
for q = 1:numel(lags)
i0 = c0 + lags(q); if i0 < 1 || i0 + nw - 1 > ns, continue; end
ys = yy(i0:i0+nw-1); cc(q) = sum(xs.*ys)/sqrt(ex*sum(ys.^2));
end
[~, q] = max(cc); d = 0;
if q > 1 && q < numel(lags), den = cc(q-1) - 2*cc(q) + cc(q+1); if den < 0, d = 0.5*(cc(q-1) - cc(q+1))/den; end, end
prev = lags(q); zcc = (c0 - 1 + nw/2)*dzs; zc(end+1) = zcc; uu(end+1) = (lags(q) + d)*dzs - sh*zcc; %#ok<AGROW>
end
ne = numel(zc) - N + 1; ee = zeros(1,ne); ze = zeros(1,ne);
for i = 1:ne, zz = zc(i:i+N-1); u2 = uu(i:i+N-1); zm = mean(zz); ee(i) = -sum((zz-zm).*(u2-mean(u2)))/sum((zz-zm).^2); ze(i) = zm; end
sel = ze > 0.004 & ze < 0.036; SN(comp+1, ie) = mean(ee(sel))/std(ee(sel));
end
fprintf('strain %6.2f %%: SNRe = %7.1f (plain), %8.1f (companded)\n', e*100, SN(1,ie), SN(2,ie));
end
figure(2); loglog(epsList*100, SN(1,:), 'r-o', epsList*100, SN(2,:), 'g--s'); grid on;
xlabel('applied strain (%)'); ylabel('SNR_e'); legend('plain', 'companded'); title('Strain filter, 2 mm kernel');
What to look for: true strain 1.111 % (background) and 0.370 % (layer); 400 µm (20.8 samples) total displacement; SNRe = 5.3 and measured strain ratio 2.98 for a true 3; strain filter 3.1, 4.1, 4.7, 5.2, 3.1, 0.4 at 0.05, 0.1, 0.3, 1, 3, 6 % (collapse at 10 %); companding 3.9 → 330 across the same range (uniform 1D strain: an upper bound).
m30_arfi_displacement.mAxisymmetric FD step response, push-duration sweep by superposition, two stiffnesses, TTP speed, heating, optimal push frequency · change cs, eta, I, a, f0 · runtime ≈ 3 s
%% m30_arfi_displacement.m
% Purpose : ARFI displacement from a radiation-force push. Axisymmetric (r,z)
% finite-difference solution of the scalar Kelvin-Voigt shear-wave
% equation rho*u_tt = mu*lap(u) + eta*lap(u_t) + F, Gaussian force
% column (sigma_r = 0.35 mm, sigma_z = 2.5 mm). The step response is
% computed once; a push of duration T is s(t) - s(t-T) (linearity).
% F = 2*alpha*I/c sets the amplitude.
% Change : cs (shear speed), eta, I (W/cm^2), a (dB/cm/MHz), f0, push list.
% Expect : alpha = 23.0 Np/m, F = 1.50e5 N/m^3 at 500 W/cm^2; peak on-axis
% displacement about 3.5 um for a 200 us push at 1.5 m/s (end-of-
% push value F*T^2/(2 rho) = 3.0 um); long pushes approach the static
% limit; time-to-peak at 2 and 4 mm gives back c_s; heating 12 mK.
rng(1);
rho = 1000; c = 1540;
I = 500; a = 0.5; f0 = 4; % W/cm^2, dB/cm/MHz, MHz
alpha = a*f0/8.686*100; % Np/m
F = 2*alpha*I*1e4/c; % N/m^3
fprintf('alpha = %.2f Np/m, F = 2 alpha I / c = %.3g N/m^3\n', alpha, F);
pushes = [50 100 200 400 800 1000]*1e-6;
dts = 1e-5; tEnd = 7e-3; ns = round(tEnd/dts) + 1;
dr = 1e-4; dz = 2e-4; nr = 81; nz = 101; sr = 0.35e-3; sz = 2.5e-3;
r = (0:nr-1)'*dr; z = ((0:nz-1) - (nz-1)/2)*dz;
f = exp(-r.^2/(2*sr^2)) * exp(-z.^2/(2*sz^2)); % nr x nz unit force
% absorbing layers
dd = zeros(nr, nz); spR = 20; spZ = 15;
for i = 1:nr, for j = 1:nz
d = 0; if i > nr - spR, d = d + ((i - (nr - spR))/spR)^2; end
jj = min(j - 1, nz - j); if jj < spZ, d = d + ((spZ - jj)/spZ)^2; end
dd(i,j) = d;
end, end
lap = @(A) [4*(A(2,:) - A(1,:))/dr^2; ...
(A(3:end,:) - 2*A(2:end-1,:) + A(1:end-2,:))/dr^2 + (A(3:end,:) - A(1:end-2,:))./(2*r(2:end-1)*dr); ...
(-2*A(end,:) + A(end-1,:))/dr^2] + ...
([zeros(nr,1) A(:,1:end-1)] - 2*A + [A(:,2:end) zeros(nr,1)])/dz^2;
csList = [1.5 3]; eta = 0.5;
for ic = 1:numel(csList)
cs = csList(ic); mu = rho*cs^2;
lam = 8/dr^2 + 4/dz^2; dt = 0.7*2/(cs*sqrt(lam)); if eta > 0, dt = min(dt, 0.7*2*rho/(eta*lam)); end
damp = exp(-dd*0.6*cs*dt/(spR*dr)*6);
u = zeros(nr, nz); v = u; S = zeros(3, ns); jc = (nz+1)/2; t = 0; nxt = 1;
while nxt <= ns
while nxt <= ns && (nxt-1)*dts <= t + 1e-12, S(:,nxt) = [u(1,jc); u(21,jc); u(41,jc)]; nxt = nxt + 1; end
v = (v + dt/rho*(mu*lap(u) + eta*lap(v) + f)).*damp; u = u + dt*v; t = t + dt;
end
% static limit (30.7): (1/(4 pi mu)) * int f/|r| dV
[RR, ZZ] = meshgrid(((1:60) - 0.5)*0.05e-3, ((1:120) - 60.5)*0.1e-3);
Is = sum(sum(exp(-RR.^2/(2*sr^2) - ZZ.^2/(2*sz^2))./sqrt(RR.^2 + ZZ.^2).*2*pi.*RR))*0.05e-3*0.1e-3;
ust = F*Is/(4*pi*mu);
fprintf('\nc_s = %.1f m/s, eta = %.1f Pa s (dt = %.2f us): static limit %.2f um\n', cs, eta, dt*1e6, ust*1e6);
tt = (0:ns-1)*dts;
for T = pushes
k = round(T/dts); y = S(1,:) - [zeros(1,k) S(1,1:end-k)];
[pk, ip] = max(y); irec = find(y(ip:end) < 0.5*pk, 1) + ip - 1;
fprintf(' push %4.0f us: peak u = %6.2f um at %.2f ms, F T^2/(2 rho) = %6.2f um, recovery to 50%% %.2f ms after peak\n', ...
T*1e6, F*pk*1e6, tt(ip)*1e3, F*T^2/(2*rho)*1e6, (irec - ip)*dts*1e3);
end
k = round(200e-6/dts); Y = S - [zeros(3,k) S(:,1:end-k)];
tp = zeros(1,3); for q = 1:3, [~, ip] = max(Y(q,:)); de = Y(q,ip-1) - 2*Y(q,ip) + Y(q,ip+1); tp(q) = (ip - 1 + 0.5*(Y(q,ip-1) - Y(q,ip+1))/de)*dts; end
fprintf(' TTP (200 us push) at r = 2 and 4 mm: %.2f and %.2f ms -> speed %.2f m/s\n', tp(2)*1e3, tp(3)*1e3, 2e-3/(tp(3) - tp(2)));
if ic == 1, figure(1); plot(tt*1e3, F*Y'*1e6); xlabel('time (ms)'); ylabel('u (\mum)'); legend('r = 0', 'r = 2 mm', 'r = 4 mm'); title('200 \mus push, c_s = 1.5 m/s'); end
end
fprintf('\nHeating per 200 us push (no conduction): %.1f mK\n', 2*alpha*I*1e4*200e-6/3.7e6*1e3);
ap = a/8.686; zf = 4; % Np/cm/MHz, cm
fprintf('Optimum push frequency at %d cm: %.2f MHz (fixed source), %.2f MHz (focal gain f^2)\n', zf, 1/(2*ap*zf), 3/(2*ap*zf));
What to look for: α = 23.03 Np/m, F = 1.5×10⁵ N/m³; at 1.5 m/s peaks of 0.92, 1.79, 3.51, 6.75, 12.1 and 14.1 µm for 50–1000 µs pushes (static limit 21.8 µm), recovery to 50 % in ≈ 0.75 ms; at 3 m/s 0.55 … 5.0 µm (static 5.4 µm), recovery ≈ 0.3 ms; TTP at 2 and 4 mm returns 1.53 and 3.05 m/s; 12.4 mK per push; optimum push frequency 2.17 MHz (6.51 MHz with f² focal gain) at 4 cm.
m30_shear_wave_tof_kf.mGreen's-function shear wave in a KV medium; TTP, cross-correlation and k–f estimators with 40 noise draws · change eta, SNR, window, frame rate · runtime ≈ 2 s
%% m30_shear_wave_tof_kf.m
% Purpose : Shear-wave speed estimators on one simulated push. Cylindrical
% wave from a 200 us line push (sigma = 0.35 mm) in a Kelvin-Voigt
% medium, far-field 2D Green's function
% U(r,w) ~ S(w) exp(-j k r) / (mu* sqrt(k r)), k = w sqrt(rho/mu*)
% on an FFT grid (frame rate 10 kHz, 512 frames), particle velocity.
% Estimators from 3 to 11 mm: time-to-peak (TTP) regression,
% (peak searched within +-2 ms of the previous trace's peak),
% cross-correlation (XC) delay regression, and the k-f phase
% velocity (2D FFT, Hann window in x). Monte Carlo over noise.
% Change : c0, eta, snrdB, x2 (window end), fr (frame rate), nMC.
% Expect : eta = 0: all estimators return 1.50 m/s. eta = 0.4 Pa s:
% c_ph(200 Hz) = 1.528 m/s, c_g(200 Hz) = 1.584 m/s; TTP on
% particle velocity reads above c_ph (about 1.6-1.9 m/s), XC close
% to c_ph at the trace centroid frequency (std over 40 noise draws
% at 30 dB printed for each estimator).
rng(1);
rho = 1000; c0 = 1.5; fr = 10e3; Nt = 512; Tp = 200e-6; sr = 0.35e-3;
dx = 0.2e-3; x = (1e-3:dx:20e-3)'; nx = numel(x);
x1 = 3e-3; x2 = 11e-3; snrdB = 30; nMC = 40;
m = (1:Nt/2)'; f = m*fr/Nt; w = 2*pi*f;
etaList = [0 0.4];
for eta = etaList
mu = rho*c0^2; ms = mu + 1j*w*eta; k = w.*sqrt(rho./ms);
if any(imag(k) > 0), k = conj(k); end % choose the decaying root
S = (1 - exp(-1j*w*Tp))./(1j*w); hp = 1 - exp(-(f/30).^2); sp = exp(-real(k).^2*sr^2/2);
V = zeros(Nt, nx);
for ix = 1:nx
G = exp(-1j*k*x(ix) + 1j*pi/4)./sqrt(k*x(ix))./ms;
Z = 1j*w.*S.*G.*hp.*sp; % particle velocity spectrum
full = zeros(Nt,1); full(2:Nt/2) = Z(1:end-1); full(Nt/2+1) = real(Z(end)); full(Nt/2+2:end) = conj(Z(end-1:-1:1));
V(:,ix) = real(ifft(full));
end
sel = find(x >= x1 - 1e-9 & x <= x2 + 1e-9); xs = x(sel);
cph = @(ff) 2*pi*ff./real(2*pi*ff.*sqrt(rho./(mu + 1j*2*pi*ff*eta)));
hh = 0.5; cg = @(ff) 2*pi*2*hh./(real(2*pi*(ff+hh).*sqrt(rho./(mu + 1j*2*pi*(ff+hh)*eta))) - real(2*pi*(ff-hh).*sqrt(rho./(mu + 1j*2*pi*(ff-hh)*eta))));
P = abs(fft(V(:, sel(1)))).^2; fcen = sum(P(2:Nt/2).*f(1:end-1))/sum(P(2:Nt/2));
pk = max(abs(V(:, sel(1))));
cT = zeros(nMC,1); cX = cT; cK = cT;
for it = 1:nMC
Vn = V(:, sel) + pk*10^(-snrdB/20)*randn(Nt, numel(sel));
% TTP
tt = zeros(numel(sel),1);
for q = 1:numel(sel)
if q == 1, rg = 1:round(0.01*fr); else, rg = max(1, imPrev - round(0.002*fr)):min(Nt, imPrev + round(0.002*fr)); end
[~, im] = max(Vn(rg,q)); im = im + rg(1) - 1; imPrev = im; d = 0;
if im > 1 && im < Nt, a = Vn(im-1,q); b = Vn(im,q); g = Vn(im+1,q); den = a - 2*b + g; if den < 0, d = 0.5*(a - g)/den; end, end
tt(q) = (im - 1 + d)/fr;
end
pT = [xs ones(size(xs))]\tt; cT(it) = 1/pT(1);
% XC versus first trace, FFT with zero padding
NF = 1024; R0 = fft(Vn(:,1), NF); lg = zeros(numel(sel),1);
for q = 1:numel(sel)
cc = real(ifft(fft(Vn(:,q), NF).*conj(R0))); maxL = round(0.03*fr); [~, im] = max(cc(1:maxL+1)); d = 0;
if im > 1, a = cc(im-1); b = cc(im); g = cc(im+1); den = a - 2*b + g; if den < 0, d = 0.5*(a - g)/den; end, end
lg(q) = (im - 1 + d)/fr;
end
pX = [xs ones(size(xs))]\lg; cX(it) = 1/pX(1);
% k-f at 200 Hz bin
nW = numel(sel); win = 0.5 - 0.5*cos(2*pi*((1:nW)' - 0.5)/nW);
TF = fft(Vn); mb = round(200/(fr/Nt)) + 1; row = TF(mb, :).'.*win;
NK = 256; K = abs(fft(row, NK)); Kn = [K(1); K(end:-1:2)]; % index q -> k = q*dk (wave e^{-jkx})
[~, iq] = max(Kn(1:NK/2)); d = 0; if iq > 1, a = Kn(iq-1); b = Kn(iq); g = Kn(iq+1); den = a - 2*b + g; if den < 0, d = 0.5*(a - g)/den; end, end
kk = (iq - 1 + d)*2*pi/(NK*dx); cK(it) = 2*pi*(mb - 1)*fr/Nt/kk;
end
fprintf('\neta = %.1f Pa s (c0 = %.2f m/s), window %.0f-%.0f mm, %d dB, %d draws\n', eta, c0, x1*1e3, x2*1e3, snrdB, nMC);
fprintf(' centroid frequency of the trace at 3 mm: %.0f Hz; c_ph = %.3f, c_g = %.3f m/s there\n', fcen, cph(fcen), cg(fcen));
fprintf(' theory: c_ph(200 Hz) = %.3f m/s, c_g(200 Hz) = %.3f m/s\n', cph(200), cg(200));
fprintf(' TTP : %.3f +- %.3f m/s\n', mean(cT), std(cT));
fprintf(' XC : %.3f +- %.3f m/s\n', mean(cX), std(cX));
fprintf(' k-f : c_ph(200 Hz) = %.3f +- %.3f m/s\n', mean(cK), std(cK));
if eta > 0
figure(1); subplot(1,2,1); imagesc(x*1e3, (0:Nt-1)/fr*1e3, V.*sqrt(x')); axis xy; ylim([0 20]); xlabel('x (mm)'); ylabel('t (ms)'); title('particle velocity, \eta = 0.4 Pa s');
subplot(1,2,2); ff = 20:5:1000; plot(ff, cph(ff), 'k', ff, cg(ff), 'r--'); xlabel('f (Hz)'); ylabel('speed (m/s)'); legend('c_{ph}', 'c_g', 'location', 'northwest'); title('Kelvin-Voigt dispersion');
end
end
% precision scaling of (30.11)
c = 2; sig_t = 50e-6; L = 10e-3; N = 51; dxs = L/(N-1);
fprintf('\n(30.11): c = 2 m/s, sigma_t = 50 us, L = 10 mm, N = 51: sigma_c/c = %.2f %%\n', 100*c*sig_t/dxs*sqrt(12/(N*(N^2-1))));
What to look for: elastic: 1.501 ± 0.005 (TTP), 1.500 ± 0.001 (XC), 1.500 ± 0.033 m/s (k–f at 200 Hz). η = 0.4 Pa·s: theory \(c_{\mathrm{ph}}\)(200 Hz) = 1.528, \(c_g\) = 1.584 m/s; trace centroid 237 Hz; TTP 1.847 ± 0.052, XC 1.546 ± 0.007, k–f 1.527 ± 0.014 m/s. The precision formula (30.11) gives 0.48 % for the worked example.
m30_viscoelastic_dispersion.mKelvin–Voigt and power-law phase/group velocity, attenuation, dispersion slope, pulse propagation · change c0, eta, n, fc · runtime ≈ 1 s
%% m30_viscoelastic_dispersion.m
% Purpose : Shear-wave dispersion in Kelvin-Voigt and power-law (spring-pot)
% media: phase velocity, group velocity (numerical d(omega)/dk),
% attenuation, the dispersion slope (200-600 Hz) and the shape
% change of a pulse propagated as a 1D plane wave,
% U(x,w) = U(0,w) exp(-j k(w) x).
% Change : c0, etaList, n (power-law exponent), fc (pulse centre frequency).
% Expect : KV, c0 = 1.5 m/s, eta = 1 Pa s: c_ph = 1.54, 1.66, 2.21, 3.16 m/s
% at 100, 200, 500, 1000 Hz; alpha(200 Hz) = 197 Np/m;
% c_g(200 Hz) = 2.01 m/s. Low-frequency approximation
% c0 (1 + 3x^2/8) within 0.1 % at x = 0.3. Power law n = 0.2:
% c_g/c_ph = 1/(1 - n/2) = 1.111 and alpha = (w/c_ph) tan(n pi/4).
rng(1);
rho = 1000; c0 = 1.5; mu = rho*c0^2;
kKV = @(f, eta) 2*pi*f.*sqrt(rho./(mu + 1j*2*pi*f*eta));
cphKV = @(f, eta) 2*pi*f./real(kKV(f, eta));
h = 0.5; cgKV = @(f, eta) 2*pi*2*h./(real(kKV(f + h, eta)) - real(kKV(f - h, eta)));
fl = [100 200 500 1000];
for eta = [0.4 1 2]
fprintf('KV eta = %.1f Pa s: c_ph =', eta); fprintf(' %.3f', cphKV(fl, eta));
fprintf(' m/s; c_g = '); fprintf(' %.3f', cgKV(fl, eta)); fprintf(' m/s; alpha = '); fprintf(' %.0f', -imag(kKV(fl, eta)));
fprintf(' Np/m; slope(200-600 Hz) = %.2f m/s/kHz\n', (cphKV(600, eta) - cphKV(200, eta))/0.4);
end
x = 0.3; fprintf('Low-frequency check x = 0.3: exact %.4f, 1 + 3x^2/8 = %.4f\n', sqrt(2*(1 + x^2)/(1 + sqrt(1 + x^2))), 1 + 3*x^2/8);
fprintf('High-frequency KV limit sqrt(2 w eta/rho) at 1000 Hz, eta = 1: %.2f m/s (exact %.2f)\n', sqrt(2*2*pi*1000*1/rho), cphKV(1000, 1));
% power law with c_ph(100 Hz) = c0
n = 0.2; w0 = 2*pi*100; cr = c0*cos(n*pi/4);
kPL = @(f) (2*pi*f/cr).*(2*pi*f/w0).^(-n/2).*(cos(n*pi/4) - 1j*sin(n*pi/4));
cphPL = @(f) 2*pi*f./real(kPL(f)); cgPL = @(f) 2*pi*2*h./(real(kPL(f + h)) - real(kPL(f - h)));
fprintf('Power law n = %.1f: c_ph = ', n); fprintf(' %.3f', cphPL(fl)); fprintf(' m/s; c_g/c_ph = '); fprintf(' %.4f', cgPL(fl)./cphPL(fl));
fprintf(' (theory %.4f)\n', 1/(1 - n/2));
fprintf(' alpha(200 Hz) = %.1f Np/m, (w/c_ph) tan(n pi/4) = %.1f Np/m; factor per decade 10^(n/2) = %.3f\n', -imag(kPL(200)), 2*pi*200/cphPL(200)*tan(n*pi/4), 10^(n/2));
% pulse propagation (KV eta = 1, fc = 250 Hz)
fs = 20e3; Nt = 1024; t = (0:Nt-1)'/fs; fc = 250; sf = 0.8*fc/2/sqrt(2*log(2)); st = 1/(2*pi*sf);
u0 = exp(-(t - 0.004).^2/(2*st^2)).*cos(2*pi*fc*(t - 0.004)); U0 = fft(u0);
fq = (0:Nt-1)'*fs/Nt; fq(fq > fs/2) = fq(fq > fs/2) - fs;
xs = [5 10 15]*1e-3; tp = zeros(1,3); U = zeros(Nt,3);
for q = 1:3
mp = (1:Nt/2-1)'; kp = kKV(mp*fs/Nt, 1); Zp = U0(mp+1).*exp(-1j*kp*xs(q));
Zf = zeros(Nt,1); Zf(mp+1) = Zp; Zf(Nt-mp+1) = conj(Zp); U(:,q) = real(ifft(Zf));
Z = fft(U(:,q)); Z(fq < 0) = 0; Z(2:Nt/2) = 2*Z(2:Nt/2); e = abs(ifft(Z));
[~, im] = max(e); de = e(im-1) - 2*e(im) + e(im+1); tp(q) = (im - 1 + 0.5*(e(im-1) - e(im+1))/de)/fs;
end
fprintf('Pulse fc = %d Hz, KV eta = 1: envelope TTP 5 -> 15 mm gives %.3f m/s (c_ph(fc) = %.3f, c_g(fc) = %.3f)\n', fc, 10e-3/(tp(3) - tp(1)), cphKV(fc, 1), cgKV(fc, 1));
ff = 10:5:1000;
figure(1); subplot(1,2,1); plot(ff, cphKV(ff, 1), 'b', ff, cgKV(ff, 1), 'r--', ff, cphPL(ff), 'k', ff, cgPL(ff), 'k--');
xlabel('f (Hz)'); ylabel('speed (m/s)'); legend('KV c_{ph}', 'KV c_g', 'PL c_{ph}', 'PL c_g', 'location', 'northwest'); title('c_0 = 1.5 m/s, \eta = 1 Pa s, n = 0.2');
subplot(1,2,2); plot(t*1e3, U); xlabel('t (ms)'); ylabel('u'); legend('5 mm', '10 mm', '15 mm'); title('KV pulse, 250 Hz');
What to look for: KV η = 1 Pa·s: \(c_{\mathrm{ph}}\) = 1.543, 1.659, 2.210, 3.159 m/s and \(c_g\) = 1.631, 2.009, 4.034, 7.099 m/s at 100, 200, 500, 1000 Hz; α(200 Hz) = 197 Np/m; slope 1.88 m/s/kHz (0.48 for η = 0.4, 3.72 for η = 2). Power law n = 0.2: \(c_g/c_{\mathrm{ph}}\) = 1.1111 at every frequency, α(200 Hz) = 123.8 Np/m both ways. The 250 Hz pulse gives an apparent envelope speed of 1.703 m/s, below both \(c_{\mathrm{ph}}\) (1.737) and \(c_g\) (2.278) at its nominal centre, because attenuation strips the high frequencies and shifts the surviving spectrum down.
m30_supersonic_mach.m2D FD supersonic source, Mach cone and the lateral TOF bias versus Mach number · change cs, V, n · runtime ≈ 1 s
%% m30_supersonic_mach.m
% Purpose : Supersonic shear source. 2D finite-difference (leapfrog) solution
% of the scalar elastic shear-wave equation in the imaging plane:
% pushes 2 mm apart along the beam axis, each lasting 2 mm / V, so
% the source moves at V. Measures the lateral time-to-peak speed on
% the row crossed by the plane part of the front and compares with
% c_s / cos(phi), sin(phi) = 1/M (30.9).
% Change : cs, Vlist, n (number of foci).
% Expect : c_s = 2 m/s, 5 foci: V = 15 m/s (M = 7.5, phi = 7.7 deg) reads
% about 2.0 m/s; V = 5 m/s (M = 2.5, phi = 23.6 deg) about 2.2 m/s;
% V = 3 m/s (M = 1.5, phi = 41.8 deg) about 2.7-2.8 m/s.
rng(1);
rho = 1000; cs = 2; n = 5; dx = 0.2e-3; nx = 121; nz = 141; x0 = -12e-3; z0 = 4e-3;
X = x0 + (0:nx-1)*dx; Zc = z0 + (0:nz-1)'*dx; [XX, ZZ] = meshgrid(X, Zc);
tRec = 0.1e-3; nRec = 80; sp = 14;
id = min(min(repmat(0:nx-1, nz, 1), repmat(nx-1:-1:0, nz, 1)), min(repmat((0:nz-1)', 1, nx), repmat((nz-1:-1:0)', 1, nx)));
dmp = max(0, sp - id)/sp;
Vlist = [15 5 3];
for V = Vlist
dzp = 2e-3; Tp = dzp/V; zs = 18e-3 + ((0:n-1) - (n-1)/2)*dzp;
nsub = ceil(tRec/(0.6*dx/(cs*sqrt(2)))); dt = tRec/nsub; damp = exp(-dmp.^2*4*cs*dt/(sp*dx)*3);
u = zeros(nz, nx); v = u; fr = zeros(nz, nx, nRec); t = 0; mu = rho*cs^2;
G = zeros(nz, nx, n); for k = 1:n, G(:,:,k) = exp(-XX.^2/(2*(0.3e-3)^2) - (ZZ - zs(k)).^2/(2*(1.2e-3)^2)); end
for r = 1:nRec
for s = 1:nsub
k = floor(t/Tp) + 1; f = zeros(nz, nx); if k <= n, f = G(:,:,k); end
L = ([u(:,2:end) zeros(nz,1)] + [zeros(nz,1) u(:,1:end-1)] + [u(2:end,:); zeros(1,nx)] + [zeros(1,nx); u(1:end-1,:)] - 4*u)/dx^2;
v = (v + dt/rho*(mu*L + 1e4*f)).*damp; u = u + dt*v; t = t + dt;
end
fr(:,:,r) = u;
end
M = V/cs; phi = asin(1/M); zrow = min(30e-3, 18e-3 + 5.5e-3/sqrt(M^2 - 1)); jr = round((zrow - z0)/dx) + 1;
cols = (61 + 15):2:(61 + 40); tp = zeros(size(cols));
for q = 1:numel(cols)
vel = squeeze(diff(fr(jr, cols(q), :))); [~, im] = max(vel); d = 0;
if im > 1 && im < numel(vel), a = vel(im-1); b = vel(im); g = vel(im+1); den = a - 2*b + g; if den < 0, d = 0.5*(a - g)/den; end, end
tp(q) = (im + d - 0.5)*tRec;
end
p = polyfit(X(cols), tp, 1);
fprintf('V = %4.1f m/s: M = %.2f, phi = %.1f deg, row z = %.1f mm, lateral TTP speed = %.2f m/s, c_s/cos(phi) = %.2f m/s\n', V, M, phi*180/pi, zrow*1e3, 1/p(1), cs/cos(phi));
if V == Vlist(1), figure(1); imagesc(X*1e3, Zc*1e3, fr(:,:,31) - fr(:,:,30)); axis image; colormap(gray); xlabel('x (mm)'); ylabel('z (mm)'); title(sprintf('particle velocity at 3 ms, M = %.1f', M)); end
end
What to look for: V = 15, 5, 3 m/s in 2 m/s tissue: φ = 7.7°, 23.6°, 41.8°; lateral TTP speed 1.99, 2.19, 2.80 m/s against \(c_s/\cos\varphi\) = 2.02, 2.18, 2.68 m/s.
m30_lamb_a0.mA0 Lamb mode of an incompressible plate (Rayleigh–Lamb root), thin-plate and immersed approximations · change h, ct · runtime ≈ 1 s
%% m30_lamb_a0.m
% Purpose : A0 (flexural) Lamb mode of an incompressible elastic plate in
% vacuum, from the Rayleigh-Lamb antisymmetric equation with
% c_L -> infinity:
% (k^2 + q'^2)^2 tanh(k h/2) - 4 k^3 q' tanh(q' h/2) = 0,
% k = w/c, q' = w sqrt(1/c^2 - 1/c_t^2), c < c_t.
% Compared with the thin-plate limit sqrt(w h c_t / sqrt(3)), the
% immersed-plate empirical fit sqrt(w h c_t / (2 sqrt(3))) of
% Couade et al. 2010, and the incompressible Rayleigh speed.
% Change : h (thickness), ct (bulk shear speed), frequency list.
% Expect : Rayleigh speed 0.9553 c_t. h = 1 mm, c_t = 3 m/s: A0 = 1.365 m/s
% at 200 Hz (apparent E 21 % of true), 2.675 m/s at 2 kHz.
rng(1);
Fn = @(c, f, h, ct) ((2*pi*f./c).^2 + (2*pi*f).^2.*max(1./c.^2 - 1/ct^2, 0)).^2.*tanh(pi*f*h./c) ...
- 4*(2*pi*f./c).^3.*(2*pi*f.*sqrt(max(1./c.^2 - 1/ct^2, 0))).*tanh(pi*f*h*sqrt(max(1./c.^2 - 1/ct^2, 0)));
r = roots([1 -8 24 -16]); r = r(abs(imag(r)) < 1e-9); fprintf('Incompressible Rayleigh speed: %.4f c_t\n', sqrt(r));
cases = [1e-3 3; 1e-3 5; 0.55e-3 2]; flist = [50 100 200 500 1000 2000];
for ic = 1:size(cases,1)
h = cases(ic,1); ct = cases(ic,2); fprintf('h = %.2f mm, c_t = %.1f m/s\n', h*1e3, ct);
for f = flist
cg = ct*linspace(1e-3, 0.999999, 600); v = Fn(cg, f, h, ct)./(2*pi*f./cg).^4;
s = find(sign(v(1:end-1)) ~= sign(v(2:end)), 1, 'last'); a = cg(s); b = cg(s+1);
g = @(c) Fn(c, f, h, ct)./(2*pi*f/c)^4; fa = g(a);
for it = 1:60, mm = (a + b)/2; fm = g(mm); if sign(fm) == sign(fa), a = mm; fa = fm; else, b = mm; end, end
cA = (a + b)/2; cK = sqrt(2*pi*f*h*ct/sqrt(3)); cF = sqrt(2*pi*f*h*ct/(2*sqrt(3)));
fprintf(' %5d Hz: A0 = %.3f m/s, thin plate = %.3f, immersed fit = %.3f, (c/c_t)^2 = %.3f\n', f, cA, cK, cF, (cA/ct)^2);
end
end
h = 1e-3; ct = 3; ff = 25:25:2000; cA = zeros(size(ff));
for i = 1:numel(ff)
f = ff(i); cg = ct*linspace(1e-3, 0.999999, 600); v = Fn(cg, f, h, ct)./(2*pi*f./cg).^4; s = find(sign(v(1:end-1)) ~= sign(v(2:end)), 1, 'last');
a = cg(s); b = cg(s+1); g = @(c) Fn(c, f, h, ct)./(2*pi*f/c)^4; fa = g(a);
for it = 1:60, mm = (a + b)/2; fm = g(mm); if sign(fm) == sign(fa), a = mm; fa = fm; else, b = mm; end, end; cA(i) = (a + b)/2;
end
figure(1); plot(ff, cA, 'b', ff, sqrt(2*pi*ff*h*ct/sqrt(3)), 'r--', ff, sqrt(2*pi*ff*h*ct/(2*sqrt(3))), 'k:', ff, 0*ff + ct, 'k-', ff, 0*ff + 0.9553*ct, 'k--');
ylim([0 1.15*ct]); xlabel('f (Hz)'); ylabel('phase velocity (m/s)'); legend('A0 exact', 'thin plate', 'immersed (Couade)', 'c_t', 'Rayleigh', 'location', 'southeast'); title('1 mm plate, c_t = 3 m/s');
What to look for: Rayleigh speed 0.9553 \(c_t\); h = 1 mm, \(c_t\) = 3 m/s: 0.723, 1.003, 1.365, 1.945, 2.373, 2.675 m/s from 50 Hz to 2 kHz, apparent \((c/c_t)^{2}\) from 0.058 to 0.795; the thin-plate formula is within 2 % at 50 Hz and 74 % high at 2 kHz.
30.14Self-check
Q1A shear wave travels at 2.5 m/s. What are μ and E?
\(\mu = \rho c^{2}\) = 1000 × 6.25 = 6.25 kPa; \(E = 3\mu\) = 18.75 kPa (elastic, isotropic, incompressible).
Q2Push at 3 MHz, absorption 0.7 dB/cm/MHz, 800 W/cm² in situ, 300 µs. Force density and heating?
α = 2.1 dB/cm = 24.2 Np/m; \(F = 2\alpha I/c\) = 2 × 24.2 × 8×10⁶ / 1540 = 2.5×10⁵ N/m³; \(\Delta T = 2\alpha IT/\rho C\) = 31 mK.
Q3SSI pushes 1.5 mm apart, 100 µs each, in 3 m/s tissue. Mach number, cone angle, lateral TOF bias?
V = 15 m/s, M = 5, φ = arcsin(0.2) = 11.5°, bias \(1/\cos\varphi - 1\) = 2.1 %.
Q4TTP regression over 8 mm with 33 traces, σt = 100 µs, c = 1.5 m/s. Relative precision (independent errors)?
Δx = 0.25 mm; \(\sigma_c/c = (c\sigma_t/\Delta x)\sqrt{12/(N(N^{2}-1))}\) = 0.6 × √(12/35 904) = 1.1 %.
Q5KV tissue with μ = 4 kPa and η = 2 Pa·s is measured by k–f at 300 Hz. What speed, and what E would a naive conversion report?
x = ωη/μ = 0.94; \(c_{\mathrm{ph}} = c_0\sqrt{2(1+x^{2})/(1+\sqrt{1+x^{2}})}\) = 2 × 1.261 = 2.52 m/s; \(3\rho c^{2}\) = 19.1 kPa against 12 kPa elastic. The excess is viscosity, not stiffness.
Q6A 0.8 mm arterial wall in water shows a 3 m/s phase velocity at 500 Hz. Estimate ct and E with the immersed-plate fit.
\(c_t = c^{2}\,2\sqrt3/(\omega h)\) = 9 × 3.464 / (3142 × 0.0008) = 12.4 m/s; \(E = 3\rho c_t^{2}\) ≈ 460 kPa. A bulk conversion of 3 m/s would have said 27 kPa.
Q7Phase-gradient map at 400 Hz in 2 m/s tissue: largest half-kernel before wrapping?
λ = 5 mm; need \(2k\Delta < \pi\), i.e. Δ < λ/4 = 1.25 mm. At Δ = 1 mm the phase step is 2π·2/5 = 2.51 rad, already close to π with noise.
Q8Ten pSWE measurements have IQR/median = 12 % in m/s. Is the exam reliable by the SRU criterion?
Yes: ≤ 15 % in m/s, equivalent to about 24 % (≤ 30 %) in kPa, because \(\delta E/E \approx 2\,\delta c/c\).
30.15Figures from the literature
Measured data behind the models of this module. Figure numbers checked against the author manuscripts or the publisher's page; links only, no images reproduced.
Guided-wave dispersion in an arterial phantom: the k–f ridge and its fit
M. Couade, M. Pernot, C. Prada, E. Messas, J. Emmerich, P. Bruneval, A. Criton, M. Fink and M. Tanter, "Quantitative assessment of arterial wall biomechanical properties using shear wave imaging," Ultrasound Med. Biol. 36(10), 1662–1676 (2010). doi:10.1016/j.ultrasmedbio.2010.07.004
Look at Fig. 4 (space–time particle velocity in a bulk agar sample versus an agar tube in water: the tube is visibly dispersive) and Fig. 5 (2D Fourier energy of the in vitro data with the extracted dispersion curve and fit). Why: it is Demo 30.4's k–f panel on real data, and the dispersion it shows comes from geometry, not viscosity (Demo 30.7).
The carotid wall in vivo: guided waves and stiffness over the cardiac cycle
Same paper. doi:10.1016/j.ultrasmedbio.2010.07.004
Look at Fig. 11 (shear wave along the carotid wall at 0.5, 2 and 4 ms after the push), Fig. 12 (dispersion curves of arterial phantoms of different shear modulus and of the in vivo common carotid wall, each against a numerical plate solution; the empirical fit (30.20), \(c^{2}\) linear in \(f\), is the paper's eq. 10) and Fig. 13 (shear modulus over one cardiac cycle). Why: \(c^{2} \propto \omega h c_t\) is the signature of a flexural wave; Fig. 13 shows that wall stiffness is pressure-dependent.
Muscle anisotropy and dispersion along and across the fibres
J.-L. Gennisson, T. Deffieux, E. Macé, G. Montaldo, M. Fink and M. Tanter, "Viscoelastic and anisotropic mechanical properties of in vivo muscle tissue assessed by supersonic shear imaging," Ultrasound Med. Biol. 36(5), 789–801 (2010). doi:10.1016/j.ultrasmedbio.2010.02.013
Look at Fig. 6 (group velocity versus probe angle to the fibres, at rest and under a 3 kg load) and Fig. 7 (phase velocity versus frequency along and across the fibres for loads of 0–5 kg). Why: Fig. 6 is (30.21) measured; Fig. 7 shows nearly non-dispersive propagation along the fibres and strong dispersion across them, so even the viscosity is anisotropic.
Liver dispersion in patients, and a literature survey of dispersion values
J. Ormachea and K. J. Parker, "Comprehensive viscoelastic characterization of tissues and the inter-relationship of shear wave (group and phase) velocity, attenuation and dispersion," Ultrasound Med. Biol. 46(12), 3448–3459 (2020). doi:10.1016/j.ultrasmedbio.2020.08.023
Look at Fig. 2 (phase velocity versus frequency in liver at steatosis grades S0, S1, S3, with the fitting bands) and Fig. 9 (published human-liver dispersion values across studies). Why: the power-law fits of (30.18) applied to clinical data, and the spread across studies that Section 30.7 attributes to band and estimator.
Clinical quality maps and push-geometry artefacts
R. G. Barr, S. R. Wilson, D. Rubens, G. Garcia-Tsao and G. Ferraioli, "Update to the Society of Radiologists in Ultrasound liver elastography consensus statement," Radiology 296(2), 263–274 (2020). doi:10.1148/radiol.2020192437
Look at Fig. 1 (quality-assessment displays of several pSWE and 2D SWE systems), Fig. 2 (reverberation artefact below the capsule) and Fig. 4 (artefact when the push is not perpendicular to the capsule). Why: these are the boundary, reflection and direction effects of Sections 30.4–30.6 seen through vendor quality maps, and the reason for the 1.5–2 cm depth rule.
30.16Go deeper
- J. Ophir, I. Céspedes, H. Ponnekanti, Y. Yazdi and X. Li, "Elastography: a quantitative method for imaging the elasticity of biological tissues," Ultrasonic Imaging 13, 111–134 (1991). The origin of quasi-static strain imaging.
- T. Varghese and J. Ophir, "A theoretical framework for performance characterization of elastography: the strain filter," IEEE TUFFC 44(1), 164–172 (1997). SNRe versus strain.
- K. R. Nightingale, M. L. Palmeri, R. W. Nightingale and G. E. Trahey, "On the feasibility of remote palpation using acoustic radiation force," J. Acoust. Soc. Am. 110(1), 625–634 (2001); and K. Nightingale, M. S. Soo, R. Nightingale and G. Trahey, "Acoustic radiation force impulse imaging: in vivo demonstration of clinical feasibility," Ultrasound Med. Biol. 28(2), 227–235 (2002).
- A. P. Sarvazyan, O. V. Rudenko, S. D. Swanson, J. B. Fowlkes and S. Y. Emelianov, "Shear wave elasticity imaging: a new ultrasonic technology of medical diagnostics," Ultrasound Med. Biol. 24(9), 1419–1435 (1998).
- J. Bercoff, M. Tanter and M. Fink, "Supersonic shear imaging: a new technique for soft tissue elasticity mapping," IEEE TUFFC 51(4), 396–409 (2004). doi:10.1109/TUFFC.2004.1295425
- L. Sandrin et al., "Transient elastography: a new noninvasive method for assessment of hepatic fibrosis," Ultrasound Med. Biol. 29(12), 1705–1713 (2003). doi:10.1016/j.ultrasmedbio.2003.07.001
- M. L. Palmeri, M. H. Wang, J. J. Dahl, K. D. Frinkley and K. R. Nightingale, "Quantifying hepatic shear modulus in vivo using acoustic radiation force," Ultrasound Med. Biol. 34(4), 546–558 (2008). TTP regression; the RANSAC version is M. H. Wang, M. L. Palmeri, V. M. Rotemberg, N. C. Rouze and K. R. Nightingale, "Improving the robustness of time-of-flight based shear wave speed reconstruction methods using RANSAC in human liver in vivo," Ultrasound Med. Biol. 36(5), 802–813 (2010). \n
- M. L. Palmeri, A. Milkowski, R. Barr et al., "Radiological Society of North America/Quantitative Imaging Biomarker Alliance shear wave speed bias quantification in elastic and viscoelastic phantoms," J. Ultrasound Med. 40(3), 569–581 (2021). doi:10.1002/jum.15609. The multi-system phantom study behind Section 30.8.
- N. C. Rouze, M. H. Wang, M. L. Palmeri and K. R. Nightingale, "Robust estimation of time-of-flight shear wave speed using a Radon sum transformation," IEEE TUFFC 57(12), 2662–2670 (2010).
- T. Deffieux, G. Montaldo, M. Tanter and M. Fink, "Shear wave spectroscopy for in vivo quantification of human soft tissues visco-elasticity," IEEE Trans. Med. Imaging 28(3), 313–322 (2009).
- A. Manduca et al., "Spatio-temporal directional filtering for improved inversion of MR elastography images," Med. Image Anal. 7(4), 465–473 (2003); T. Deffieux, J.-L. Gennisson, J. Bercoff and M. Tanter, "On the effects of reflected waves in transient shear wave elastography," IEEE TUFFC 58(10), 2032–2035 (2011).
- P. Kijanka and M. W. Urban, "Local phase velocity based imaging: a new technique used for ultrasound shear wave elastography," IEEE Trans. Med. Imaging 38(4), 894–908 (2019).
- M. Bernal, I. Nenadic, M. W. Urban and J. F. Greenleaf, "Material property estimation for tubes and arteries using ultrasound radiation force and analysis of propagating modes," J. Acoust. Soc. Am. 129(3), 1344–1354 (2011).
- Lab cross-links: delay estimation and the strain-filter kernel optimum (Module 11); 2D FFTs (Module 04); output limits and radiation-force balances (Module 22); slow-time estimators (Module 29); sound speed and aberration, which bias tracking (Modules 32, 33).