SOCR ≫ BPAD1 Website ≫ BPAD GitHub ≫

Chapter 3: Ultrasound and Photoacoustic Imaging

Ultrasound asks what tissue is like mechanically: how stiff, how dense, how fast it moves. Photoacoustics asks what tissue is made of chemically, and then answers in the language of sound so that the answer can escape from centimeters deep. One detector, two questions, one wave equation.

How to Use This Chapter

Ultrasound is the most widely used imaging technology in medicine: real-time, portable, inexpensive, and free of ionizing radiation. Its physics is the physics of mechanical waves, pressure disturbances that propagate through tissue, reflect at boundaries, scatter off small structures, and shift in frequency when they bounce off moving blood. From those few principles follow anatomical imaging (A-, B-, and M-mode), blood-flow imaging (Doppler), and tissue-stiffness imaging (elastography).

Photoacoustic imaging sits at the boundary between this chapter and the previous one. A short laser pulse is absorbed by a molecule, exactly the optical absorption studied in Chapter 2, and the resulting flash of heat makes the tissue expand and launch an ultrasound wave that the same transducers detect. Photoacoustics therefore marries optical contrast (what the tissue is made of) with acoustic resolution (where it is, deep below the surface), escaping the depth limit that scattering imposes on purely optical methods.

The chapter follows a deliberate arc:

\[\underbrace{\text{acoustic waves}}_{\text{propagation, }Z,\ \text{attenuation}} \rightarrow \underbrace{\text{transducers and beamforming}}_{\text{how a beam is made}} \rightarrow \underbrace{\text{imaging modes}}_{\text{A, B, M, Doppler}} \rightarrow \underbrace{\text{artifacts and safety}}_{\text{when assumptions fail}}\]

\[\rightarrow \underbrace{\text{advanced methods}}_{\text{harmonic, CEUS, elastography}} \rightarrow \underbrace{\text{photoacoustics}}_{p_0 = \Gamma\mu_a F} \rightarrow \underbrace{\text{multispectral tomography}}_{sO_2\ \text{unmixing}} \rightarrow \underbrace{\text{disease}}_{\text{clinical decisions}}\]

This is a quantitative, computational chapter that applies the mathematics of Chapter 1 rather than re-deriving it; Table 1 maps every earlier result the chapter depends on to the section that uses it. Every section closes with a Section summary suitable for a lecture slide and a Checkpoint question suitable for discussion. Three end-of-chapter computational laboratories build the core signal-processing skills: pulse-echo ranging, Doppler velocity estimation with aliasing, and photoacoustic spectral unmixing.

Chapter-level clinical puzzle. A cardiologist must watch a heart valve move in real time and measure a jet velocity of 4 m/s at 12 cm depth. A hepatologist must learn whether a liver is scarred without a needle biopsy. An oncologist wants to know whether a tumour is starved of oxygen without injecting a contrast dye. And a radiologist must decide whether the faint echoes inside a small “cyst” are real debris or an artifact of the instrument. One family of methods, sound, and sound generated by light, answers all four. Which physical property does each measurement exploit, and which of the four questions turns out to be impossible with the obvious method?

Acoustic story for the chapter. Follow a single pulse. A piezoelectric crystal vibrates and launches a pressure wave into tissue. At a smooth organ boundary the wave partly reflects, and its echo timing tells us depth (B-mode). Off the tiny red blood cells it scatters weakly in all directions, and because those cells are moving, the scattered echo returns at a shifted frequency (Doppler). Off countless sub-resolution scatterers the echoes interfere coherently and produce speckle. Now replace the crystal’s electrical kick with a nanosecond flash of laser light: a haemoglobin molecule absorbs a photon, warms by a few millikelvin, expands, and itself becomes the source of a new pressure wave (photoacoustics). The same detector listens in every case.

Key idea. Ultrasound images the mechanical properties of tissue (impedance, motion, stiffness). Photoacoustics adds optical-absorption contrast while keeping acoustic resolution, letting us see molecular composition, above all, blood oxygenation, deep inside the body.

Learning Objectives

# Objective Section
1 Relate sound speed and acoustic impedance to tissue mechanical properties, and compute reflection and transmission at boundaries. 3.1
2 Use the impedance relation \(p = Zu\) and the intensity relation \(I = p^2/2\rho c\) to connect pressure, motion, and power. 3.1
3 Convert correctly between nepers and decibels, and between pressure and intensity attenuation coefficients. 3.1
4 Explain attenuation, scattering, and speckle, and justify the frequency–resolution–penetration trade-off quantitatively. 3.1
5 Show that fully developed speckle has a fixed envelope SNR of 1.91 and explain why that proves speckle is not additive noise. 3.1
6 Distinguish axial, lateral, and elevational resolution and predict how each depends on frequency, pulse length, and aperture. 3.1
7 Describe piezoelectric transducers, matching and backing layers, array geometries, and delay-and-sum beamforming. 3.2
8 Explain time-gain compensation and logarithmic compression, and compute the frame-rate budget. 3.2
9 Apply the pulse-echo range equation and the Doppler equation, including angle-error propagation. 3.3
10 Derive the pulsed-Doppler depth–velocity limit \(v_{\max}d_{\max}\le c^2/8f_0\) and use it to decide between PW and CW. 3.3
11 Explain each major ultrasound artifact in terms of the specific scanner assumption it violates. 3.4
12 Explain harmonic imaging, contrast-enhanced ultrasound, and shear-wave elastography, and convert shear-wave speed to Young’s modulus. 3.5
13 Use the Minnaert resonance to explain why microbubble agents work in the diagnostic band. 3.5
14 Compute the mechanical and thermal indices, state the regulatory limits, and apply the ALARA principle. 3.6
15 Derive \(p_0 = \Gamma\mu_a F\), interpret the Grüneisen parameter, and verify the stress- and thermal-confinement conditions numerically. 3.7
16 Relate absorber size to emitted acoustic frequency and explain why detector bandwidth limits photoacoustic resolution. 3.7
17 Explain multispectral unmixing for \(sO_2\), locate the isosbestic point, and identify spectral colouring as the dominant confounder. 3.8
18 Connect each method to its principal disease applications in cardiology, hepatology, and oncology. 3.9
19 Implement pulse-echo ranging, Doppler estimation with aliasing, and photoacoustic unmixing in R. Labs

Notation

Symbol Meaning Units
\(p,\ p_0\) acoustic pressure; photoacoustic initial pressure Pa
\(u\) particle velocity m s\(^{-1}\)
\(c,\ c_s\) longitudinal sound speed; shear-wave speed m s\(^{-1}\)
\(\rho\) mass density kg m\(^{-3}\)
\(K,\ \mu,\ E,\ \nu\) bulk modulus, shear modulus, Young’s modulus, Poisson ratio Pa (— for \(\nu\))
\(Z=\rho c\) acoustic impedance rayl (1 MRayl \(=10^6\) rayl)
\(R,\ R_I,\ T_I\) amplitude reflection coefficient; intensity reflection, transmission ,
\(\alpha,\ \mu_I\) pressure and intensity attenuation coefficients Np m\(^{-1}\)
\(a\) clinical attenuation coefficient dB cm\(^{-1}\) MHz\(^{-1}\)
\(I\) acoustic intensity (\(I_{\mathrm{SPTA}}\) = spatial-peak temporal-average) W m\(^{-2}\)
\(f_0,\ f_D\) transmit frequency; Doppler shift Hz
PRF pulse repetition frequency Hz
\(F_\#\) f-number = focal length / aperture ,
\(N\) number of cycles in a pulse, or number of array elements ,
MI, TI mechanical index, thermal index ,
\(\mu_a,\ \mu_s'\) optical absorption and reduced scattering coefficients (Ch. 2) cm\(^{-1}\)
\(F\) optical fluence J m\(^{-2}\)
\(\Gamma=\beta c^2/C_p\) Grüneisen parameter ,
\(\beta,\ C_p\) thermal expansion coefficient; specific heat K\(^{-1}\), J kg\(^{-1}\)K\(^{-1}\)
\(\varepsilon_i(\lambda),\ c_i\) molar extinction coefficient; concentration cm\(^{-1}\)M\(^{-1}\), M
\(sO_2\) blood oxygen saturation ,
\(\delta\) optical penetration depth (Ch. 2) mm
\(\tau_p,\ \tau_{\text{st}},\ \tau_{\text{th}}\) laser pulse duration; stress and thermal confinement times s

Bridges from Earlier Chapters

Table 1: Results from earlier chapters that this chapter applies.
Earlier result Chapter Role here Section
Wave equation (PDE) 1 acoustic propagation 3.1
Fourier transform 1 beam pattern as the transform of the aperture; spectral Doppler 3.2, 3.3
Sampling and the Nyquist criterion 1 pulsed-Doppler aliasing and the depth–velocity limit 3.3
Exponential decay 1, 2 acoustic attenuation; optical fluence decay 3.1, 3.7
Least squares and linear systems 1 multispectral unmixing 3.8, Lab 3
Error propagation, delta method 1 Doppler angle sensitivity; \(sO_2\) uncertainty 3.3, Lab 3
Rayleigh and Gaussian statistics 1 speckle envelope distribution and its fixed SNR 3.1
Optical absorption \(\mu_a\) 2 photoacoustic contrast 3.7
Penetration depth \(\delta=1/\sqrt{3\mu_a(\mu_a+\mu_s')}\) 2 photoacoustic depth budget 3.7
Isosbestic point; two-wavelength NIRS 2 generalized to \(M\)-wavelength \(sO_2\) unmixing 3.8
Rayleigh scattering (\(\lambda^{-4}\)) 2 acoustic analogue (\(f^{4}\)) for blood 3.1

3.1 Fundamentals of Acoustic Wave Physics

3.1.1 The acoustic wave equation

Sound in tissue is a longitudinal pressure wave: molecules oscillate back and forth along the direction of travel, creating alternating compressions and rarefactions. Three linearized relations produce the wave equation. Conservation of mass (continuity) relates density change to the divergence of particle velocity; Euler’s equation relates particle acceleration to the pressure gradient; and the equation of state relates pressure to density change through the bulk modulus \(K\):

\[\frac{\partial \rho'}{\partial t} = -\rho_0\nabla\!\cdot\!\mathbf{u}, \qquad \rho_0\frac{\partial \mathbf{u}}{\partial t} = -\nabla p, \qquad p = K\frac{\rho'}{\rho_0}.\]

Eliminating \(\rho'\) and \(\mathbf{u}\) gives the wave equation introduced in Chapter 1,

\[\begin{equation} \frac{\partial^2 p}{\partial t^2} = c^2\,\nabla^2 p \tag{1} \end{equation}\]

whose plane-wave solutions \(p(x,t) = p_A\cos(kx - \omega t)\) travel at the speed of sound \(c\), with \(\omega = 2\pi f\), \(k = 2\pi/\lambda\), and

\[\begin{equation} c = f\,\lambda . \tag{2} \end{equation}\]

The elimination also derives the sound speed rather than asserting it:

\[\begin{equation} c = \sqrt{\frac{K}{\rho}} , \tag{3} \end{equation}\]

so stiffer tissue transmits sound faster and denser tissue transmits it slower. Soft tissues vary little in \(c\) (about 1450–1600 m/s, average 1540 m/s), and every scanner assumes this single value when converting echo time into depth. Bone (\(\approx 4080\) m/s) and air (\(\approx 343\) m/s) are the dramatic exceptions, and both cause artifacts (Section 3.4).

Connection to Chapter 1. Equation (1) is the wave equation from the Chapter 1 treatment of partial differential equations. Everything in clinical ultrasound, pulse shapes, beam patterns, Doppler shifts, is a consequence of this one linear PDE plus the boundary conditions imposed by tissue interfaces.

3.1.2 Particle velocity, impedance, and intensity

A pressure wave also moves the medium. The particle velocity \(u\) is the oscillation speed of the tissue itself, of order millimeters per second, not to be confused with the wave speed \(c\) of about 1540 m/s. For a plane wave, Euler’s equation ties the two together through the acoustic impedance:

\[\begin{equation} Z = \rho\,c , \qquad p = Z\,u . \tag{4} \end{equation}\]

Equation (4) is the acoustic analogue of Ohm’s law, with pressure playing the role of voltage, particle velocity that of current, and \(Z\) that of resistance. This is what makes \(Z\) a physical property rather than a bookkeeping device: a high-impedance medium demands a large pressure to produce a given motion. Impedance is measured in rayl (kg m\(^{-2}\) s\(^{-1}\)); biomedical values are conveniently expressed in MRayl.

The intensity, the power carried per unit area, follows directly:

\[\begin{equation} I = \frac{p_A^{2}}{2\rho c} = \frac{p_{\text{rms}}^{2}}{Z} . \tag{5} \end{equation}\]

Equation (5) explains why intensity coefficients are the squares of amplitude coefficients, and it is the bridge to the safety indices of Section 3.6.

Worked Example 3.1 (how much does tissue actually move?). A diagnostic pulse reaches \(p_A = 1\) MPa in soft tissue (\(Z = 1.63\) MRayl, \(\rho = 1050\) kg m\(^{-3}\), \(c = 1540\) m/s). From Eq. (4),

\[u = \frac{p_A}{Z} = \frac{10^{6}}{1.63\times10^{6}} = 0.61\ \text{m s}^{-1},\]

and at 5 MHz the displacement amplitude is \(u/\omega = 0.61/(2\pi\times5\times10^{6}) = 19\) nm, about a fiftieth of a wavelength of visible light. From Eq. (5),

\[I = \frac{(10^{6})^{2}}{2\times1050\times1540} = 3.1\times10^{5}\ \text{W m}^{-2} = 31\ \text{W cm}^{-2}.\]

The tissue barely moves, yet the instantaneous intensity is large. Because diagnostic pulses have duty cycles of order \(10^{-3}\), the time-averaged intensity is far smaller, and it is the time average that governs heating (Section 3.6).

3.1.3 Reflection, transmission, and refraction

When a wave meets a smooth, flat boundary between media of impedance \(Z_1\) and \(Z_2\) at normal incidence, the amplitude reflection coefficient is

\[\begin{equation} R = \frac{Z_2 - Z_1}{Z_2 + Z_1} . \tag{6} \end{equation}\]

Note that \(R\) is signed: when \(Z_2 < Z_1\) it is negative, meaning the reflected pressure wave is inverted in phase. That sign carries real information, and it is why Eq. (6) must be called a coefficient and not a “fraction.” The intensity reflection and transmission coefficients are

\[\begin{equation} R_I = \left(\frac{Z_2 - Z_1}{Z_2 + Z_1}\right)^{2}, \qquad T_I = 1 - R_I = \frac{4 Z_1 Z_2}{(Z_1 + Z_2)^{2}} . \tag{7} \end{equation}\]

The lesson of Eq. (7) is that echoes come from impedance mismatches. Two tissues with similar \(Z\) produce a faint echo and most energy continues deeper; a large mismatch produces a strong echo and leaves little energy to image beyond it.

If the wave strikes a boundary obliquely it also refracts, bending according to Snell’s law \(\sin\theta_t/\sin\theta_i = c_2/c_1\), which displaces structures from their true position (Section 3.4).

From the equation to the clinic. Soft tissue and air differ enormously in impedance, so \(R_I \approx 0.999\): essentially all the sound reflects at a tissue–air interface. This is why the sonographer applies coupling gel to exclude the air layer between probe and skin, and why ultrasound cannot image through lung or bowel gas. A tissue–bone interface reflects about 43%, casting acoustic shadows behind bone.

tissue <- data.frame(
  material = c("Air","Lung","Fat","Water","Blood","Liver","Muscle","Bone","PZT"),
  rho = c(1.2, 400, 950, 1000, 1060, 1060, 1050, 1900, 7500),   # kg/m^3
  c   = c(343, 650, 1450, 1480, 1570, 1550, 1580, 4080, 4000))  # m/s
tissue$Z <- tissue$rho * tissue$c / 1e6                          # MRayl
tissue$material <- factor(tissue$material, levels = tissue$material)

p_imp <- ggplot(tissue, aes(material, Z)) +
  geom_col(fill = bpad_pal[1]) +
  geom_text(aes(label = sprintf("%.2f", Z)), vjust = -0.4, size = 2.6) +
  scale_y_log10() +
  labs(x = NULL, y = "Acoustic impedance Z (MRayl, log scale)",
       title = "Acoustic impedance of tissues") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

## Intensity reflection at selected interfaces, computed from the SAME table
refl <- function(Z1, Z2) ((Z2 - Z1)/(Z2 + Z1))^2
Zof  <- function(m) tissue$Z[tissue$material == m]
interfaces <- data.frame(
  interface = c("Tissue/Air","Tissue/Bone","Tissue/Lung","Fat/Muscle","Liver/Blood"),
  RI = c(refl(Z_soft, Zof("Air")),  refl(Z_soft, Zof("Bone")),
         refl(Z_soft, Zof("Lung")), refl(Zof("Fat"), Zof("Muscle")),
         refl(Zof("Liver"), Zof("Blood"))))
interfaces$interface <- factor(interfaces$interface, levels = interfaces$interface)

p_ref <- ggplot(interfaces, aes(interface, RI)) +
  geom_col(fill = bpad_pal[2]) +
  geom_text(aes(label = ifelse(RI > 0.01, sprintf("%.1f%%", 100*RI),
                               sprintf("%.1e", RI))),
            vjust = -0.5, size = 2.8) +
  scale_y_log10() +
  coord_cartesian(ylim = c(1e-5, 5)) +
  labs(x = NULL, y = "Intensity reflected (log scale)",
       title = "Reflection at interfaces") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

bpad_grid(p_imp, p_ref, ncol = 2)
Acoustic impedance of common materials (left) and the intensity reflection coefficient at clinically important interfaces (right), both on logarithmic axes so that the six orders of magnitude between a soft-tissue boundary and a tissue-air boundary are visible at once. The tissue-air mismatch reflects 99.9 percent of the energy, which is why coupling gel is essential; the liver-blood boundary reflects less than one part in ten thousand, which is why solid organs look nearly uniform apart from speckle.

Figure 1: Acoustic impedance of common materials (left) and the intensity reflection coefficient at clinically important interfaces (right), both on logarithmic axes so that the six orders of magnitude between a soft-tissue boundary and a tissue-air boundary are visible at once. The tissue-air mismatch reflects 99.9 percent of the energy, which is why coupling gel is essential; the liver-blood boundary reflects less than one part in ten thousand, which is why solid organs look nearly uniform apart from speckle.

knitr::kable(data.frame(
  interface = as.character(interfaces$interface),
  R_intensity = signif(interfaces$RI, 3),
  percent = sprintf("%.4f%%", 100*interfaces$RI),
  transmitted = sprintf("%.4f%%", 100*(1 - interfaces$RI))),
  col.names = c("Interface","R_I","Reflected","Transmitted"),
  caption = "Intensity reflection at clinically important boundaries, computed from the impedance table above.")
Table 2: Intensity reflection at clinically important boundaries, computed from the impedance table above.
Interface R_I Reflected Transmitted
Tissue/Air 9.99e-01 99.8990% 0.1010%
Tissue/Bone 4.26e-01 42.5790% 57.4210%
Tissue/Lung 5.25e-01 52.5433% 47.4567%
Fat/Muscle 8.59e-03 0.8594% 99.1406%
Liver/Blood 4.11e-05 0.0041% 99.9959%

3.1.4 Attenuation: nepers, decibels, and the frequency law

As a pulse travels it loses intensity to absorption (acoustic energy converted to heat) and scattering (energy redirected out of the beam). The combined loss is exponential in depth,

\[\begin{equation} I(x) = I_0\, e^{-\mu_I x}, \tag{8} \end{equation}\]

exactly the form of radioactive decay and of the Beer–Lambert optical law in Chapter 2; only the physical mechanism differs. In clinical practice attenuation is quoted in decibels and, to good approximation, grows with frequency:

\[\begin{equation} \text{attenuation (dB)} \approx a\, f\,[\text{MHz}]\, x\,[\text{cm}], \qquad a \approx 0.5\ \text{dB cm}^{-1}\text{MHz}^{-1}\ \text{(soft tissue)} . \tag{9} \end{equation}\]

More precisely, tissue follows a power law \(\alpha = a f^{\,b}\) with \(b \approx 1.0\)–1.5 (about 1.1 in liver, 1.5 in breast; \(b = 2\) in pure water). The linear approximation \(b = 1\) is the clinical working rule and is what Eq. (9) encodes.

Nepers and decibels are not the same unit. Equation (8) is in nepers; Eq. (9) is in decibels. Mixing them produces factor-of-8.686 and factor-of-two errors that are among the most common mistakes in student work. Converting for intensity,

\[\begin{equation} \text{dB} = 10\log_{10}\frac{I_0}{I} = 10\,\mu_I x\log_{10}e = 4.343\,\mu_I x , \tag{10} \end{equation}\]

and for pressure amplitude, \(p = p_0 e^{-\alpha x}\),

\[\begin{equation} \text{dB} = 20\log_{10}\frac{p_0}{p} = 8.686\,\alpha x . \tag{11} \end{equation}\]

Because \(I \propto p^2\) we have \(\mu_I = 2\alpha\), and both expressions return the same number of decibels. The clinical coefficient \(a\) is an amplitude figure.

depth <- seq(0, 20, by = 0.1)
freqs <- c(2.5, 5, 10, 15)
att <- do.call(rbind, lapply(freqs, function(f)
  data.frame(depth = depth, level = -a_dB*f*2*depth,   # round trip
             f = factor(paste0(f, " MHz"), levels = paste0(freqs, " MHz")))))

p_att <- ggplot(att, aes(depth, level, colour = f)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = -60, linetype = "dashed", colour = "grey40") +
  annotate("text", x = 19.5, y = -55, hjust = 1, size = 2.9, colour = "grey35",
           label = "-60 dB round trip") +
  coord_cartesian(ylim = c(-120, 0)) +
  scale_colour_manual(values = bpad_pal[1:4]) +
  labs(x = "Depth (cm)", y = "Round-trip level (dB)",
       title = "Penetration versus frequency")

f_seq <- seq(1, 20, by = 0.2)
p_res <- ggplot(data.frame(f = f_seq, ax = 2*c_soft/(2*f_seq*1e6)*1e3),
                aes(f, ax)) +
  geom_line(linewidth = 0.9, colour = bpad_pal[1]) +
  labs(x = "Frequency (MHz)", y = "Axial resolution (mm)",
       title = "Axial resolution versus frequency",
       subtitle = "2-cycle pulse, c = 1540 m/s")
bpad_grid(p_att, p_res, ncol = 2)
The frequency trade-off, and the unit bookkeeping behind it. Left: higher-frequency pulses attenuate faster, limiting penetration; the dashed line marks a representative usable level of -60 dB round trip. Right: axial resolution improves linearly with frequency. Together they force the depth-versus-detail compromise that determines every clinical probe choice.

Figure 2: The frequency trade-off, and the unit bookkeeping behind it. Left: higher-frequency pulses attenuate faster, limiting penetration; the dashed line marks a representative usable level of -60 dB round trip. Right: axial resolution improves linearly with frequency. Together they force the depth-versus-detail compromise that determines every clinical probe choice.

## Unit bookkeeping at one representative frequency
f_ref    <- 5
a_at_f   <- a_dB*f_ref                       # dB/cm, amplitude
alpha_Np <- a_at_f/NP_TO_DB_P                # Np/cm, pressure
muI_Np   <- 2*alpha_Np                       # Np/cm, intensity
knitr::kable(data.frame(
  quantity = c("clinical coefficient a", "pressure alpha", "intensity mu_I",
               "intensity mean free path 1/mu_I", "depth of -3 dB (one way)",
               "depth of -60 dB (round trip)"),
  value = c(sprintf("%.2f dB/cm", a_at_f), sprintf("%.4f Np/cm", alpha_Np),
            sprintf("%.4f Np/cm", muI_Np), sprintf("%.2f cm", 1/muI_Np),
            sprintf("%.2f cm", 3/a_at_f), sprintf("%.2f cm", 60/(2*a_at_f)))),
  col.names = c("Quantity", "Value at 5 MHz"),
  caption = "Attenuation expressed four ways at 5 MHz. The intensity mean free path is the 1/e depth of intensity, which is neither the 1/e depth of pressure nor the -3 dB depth; distinguishing them is the whole point of Eqs. (np-db-intensity) and (np-db-pressure).")
Table 3: Attenuation expressed four ways at 5 MHz. The intensity mean free path is the 1/e depth of intensity, which is neither the 1/e depth of pressure nor the -3 dB depth; distinguishing them is the whole point of Eqs. (np-db-intensity) and (np-db-pressure).
Quantity Value at 5 MHz
clinical coefficient a 2.50 dB/cm
pressure alpha 0.2878 Np/cm
intensity mu_I 0.5756 Np/cm
intensity mean free path 1/mu_I 1.74 cm
depth of -3 dB (one way) 1.20 cm
depth of -60 dB (round trip) 12.00 cm
cat(sprintf("Consistency check: 1/mu_I = %.2f cm corresponds to %.2f dB of one-way loss,\n",
            1/muI_Np, a_at_f/muI_Np))
## Consistency check: 1/mu_I = 1.74 cm corresponds to 4.34 dB of one-way loss,
cat(sprintf("and 10*log10(e) = %.4f, so %.2f Np = %.2f dB. Both routes agree.\n",
            NP_TO_DB_I, 1.0, NP_TO_DB_I))
## and 10*log10(e) = 4.3429, so 1.00 Np = 4.34 dB. Both routes agree.

3.1.5 Scattering and speckle

Two kinds of redirected sound matter. Specular reflection occurs at large, smooth boundaries (organ capsules, vessel walls) much bigger than a wavelength; these give the bright, orientation-dependent outlines in a B-mode image. Scattering occurs off structures comparable to or smaller than a wavelength. Off very small scatterers such as red blood cells (about 7 µm, far smaller than the ~300 µm wavelength at 5 MHz), scattering is in the Rayleigh regime with cross-section rising as \(f^{4}\), the acoustic analogue of the \(\lambda^{-4}\) optical law of Chapter 2. This weak backscatter from blood is precisely what Doppler imaging detects.

The coherent sum of countless tiny scattered wavelets within a resolution cell produces speckle, the granular texture of every ultrasound image.

Common misconception. “Speckle is noise to be removed.”

Correction. Speckle is a coherent interference pattern, not additive noise. For a fixed probe position it is deterministic and reproducible, and its statistics carry tissue information. The decisive evidence is quantitative: for fully developed speckle the envelope follows a Rayleigh distribution, so its signal-to-noise ratio is

\[\mathrm{SNR} = \frac{\mathbb{E}[A]}{\mathrm{sd}(A)} = \frac{\sqrt{\pi/2}}{\sqrt{2-\pi/2}} = 1.913,\]

a pure number with no dependence on transmit power, receive gain, or scatterer density. Genuine additive noise would fall relative to signal as transmit power increased; speckle does not, because the “noise” is the signal, coherently interfering. Only averaging independent looks of the same tissue, different angles (spatial compounding) or different frequencies (frequency compounding), reduces it, and every such look costs frame rate.

set.seed(11)
n <- 220
smooth2 <- function(M, sx, sy) {
  N <- nrow(M)
  k <- outer(1:N, 1:N, function(i, j)
    exp(-((i - N/2)^2/(2*sy^2) + (j - N/2)^2/(2*sx^2))))
  Re(fft(fft(M)*fft(k), inverse = TRUE))/(N*N)
}
## Coherent sum of many random scatterers, band-limited by the system PSF.
## Anisotropic kernel: finer axially (sy) than laterally (sx).
re  <- smooth2(matrix(rnorm(n*n), n), sx = 10, sy = 3)
im  <- smooth2(matrix(rnorm(n*n), n), sx = 10, sy = 3)
env <- Mod(complex(real = re, imaginary = im))

img_df <- expand.grid(x = 1:n, y = 1:n); img_df$v <- as.vector(env)
p_sp <- ggplot(img_df, aes(x, y, fill = v)) + geom_raster() +
  scale_fill_gradient(low = "black", high = "white", guide = "none") +
  coord_equal() + labs(x = NULL, y = NULL, title = "Simulated speckle field") +
  theme(axis.text = element_blank())

sigma <- sqrt(mean(env^2)/2)
xs <- seq(0, max(env), length.out = 200)
p_hist <- ggplot(data.frame(v = as.vector(env)), aes(v)) +
  geom_histogram(aes(y = after_stat(density)), bins = 60, fill = "grey75", colour = NA) +
  geom_line(data = data.frame(x = xs, d = xs/sigma^2*exp(-xs^2/(2*sigma^2))),
            aes(x, d), colour = bpad_pal[2], linewidth = 0.9, inherit.aes = FALSE) +
  labs(x = "Envelope amplitude", y = "Density", title = "Rayleigh envelope")

gains <- c(0.25, 0.5, 1, 2, 4, 8)
snr_g <- vapply(gains, function(g) mean(g*env)/sd(g*env), numeric(1))
p_snr <- ggplot(data.frame(gain = gains, snr = snr_g), aes(gain, snr)) +
  geom_hline(yintercept = sqrt(pi/2)/sqrt(2 - pi/2), linetype = "dashed",
             colour = bpad_pal[2]) +
  geom_line(colour = bpad_pal[1], linewidth = 0.9) +
  geom_point(size = 2.2, colour = bpad_pal[1]) +
  scale_x_log10() + coord_cartesian(ylim = c(0, 3)) +
  annotate("text", x = 0.3, y = 2.2, hjust = 0, size = 3, colour = bpad_pal[2],
           label = "theory 1.913") +
  labs(x = "Transmit gain (log scale)", y = "Envelope SNR = mean / sd",
       title = "SNR is invariant")
bpad_grid(p_sp, p_hist, p_snr, ncol = 3)
Fully developed speckle. Left: a simulated B-mode speckle field formed by the coherent sum of many sub-resolution scatterers, band-limited by an anisotropic point-spread function; the elongated texture reflects the fact that lateral resolution is coarser than axial. Centre: the envelope histogram follows the Rayleigh distribution (red curve). Right: the envelope SNR is a constant 1.91 across a 32-fold range of transmit gain, which is the quantitative proof that speckle is coherent interference rather than additive noise.

Figure 3: Fully developed speckle. Left: a simulated B-mode speckle field formed by the coherent sum of many sub-resolution scatterers, band-limited by an anisotropic point-spread function; the elongated texture reflects the fact that lateral resolution is coarser than axial. Centre: the envelope histogram follows the Rayleigh distribution (red curve). Right: the envelope SNR is a constant 1.91 across a 32-fold range of transmit gain, which is the quantitative proof that speckle is coherent interference rather than additive noise.

cat(sprintf("measured envelope SNR = %.3f;  theory sqrt(pi/2)/sqrt(2-pi/2) = %.3f\n",
            mean(env)/sd(env), sqrt(pi/2)/sqrt(2 - pi/2)))
## measured envelope SNR = 1.913;  theory sqrt(pi/2)/sqrt(2-pi/2) = 1.913
cat(sprintf("SNR across a %.0f-fold gain range: min %.3f, max %.3f (i.e. constant)\n",
            max(gains)/min(gains), min(snr_g), max(snr_g)))
## SNR across a 32-fold gain range: min 1.913, max 1.913 (i.e. constant)

3.1.6 Resolution in three dimensions

Spatial resolution has three distinct components, not two, and the one most often forgotten causes a clinically important artifact.

Axial resolution, the ability to separate two reflectors along the beam, is set by the spatial pulse length. For a pulse of \(N\) cycles,

\[\begin{equation} \Delta_{\text{axial}} = \frac{N\lambda}{2} = \frac{N c}{2 f} . \tag{12} \end{equation}\]

Lateral resolution, the ability to separate two reflectors side by side in the image plane, is set by the beam width at the depth of interest,

\[\begin{equation} \Delta_{\text{lateral}} \approx F_\#\,\lambda = \frac{\text{focal length}}{\text{aperture}}\,\lambda , \tag{13} \end{equation}\]

so it improves at higher frequency and with a larger, well-focused aperture.

Elevational resolution, perpendicular to the image plane, is set by a fixed mechanical lens on a conventional 1-D array. It is therefore the worst of the three and cannot be adjusted electronically; it degrades away from the lens focus. The consequence is slice-thickness (partial-volume) artifact: echoes from structures just outside the nominal image plane are displayed as if they lay inside it, filling small anechoic structures with spurious low-level echo. A small simple cyst that appears to contain debris is very often showing slice-thickness artifact rather than real internal echo. Matrix and 1.5-D arrays add elevational focusing to reduce it.

The fundamental trade-off. Raising the frequency improves both in-plane resolutions (Eqs. (12)(13)) but worsens penetration (Eq. (9)). Clinical frequency choice is therefore a depth-versus-detail decision: 2–5 MHz for deep abdominal and cardiac imaging, 7–15 MHz for superficial vascular and small-parts imaging, and 20–50+ MHz for skin, eye, and intravascular ultrasound, where only a few millimeters of penetration are needed.

Worked Example 3.2 (the frequency trade in numbers). At 5 MHz in soft tissue, \(\lambda = c/f = 1540/(5\times10^{6}) = 0.31\) mm. A 2-cycle pulse gives \(\Delta_{\text{axial}} = 2(0.31)/2 = 0.31\) mm, and a 60 dB round-trip budget reaches \(60/(2\times0.5\times5) = 12\) cm. Doubling to 10 MHz halves \(\lambda\) and the axial resolution to 0.15 mm, but halves the usable depth to 6 cm. Resolution and depth trade inversely and exactly: the product \(\Delta_{\text{axial}} \times d_{\max}\) is approximately constant at fixed pulse length and dB budget.

Section 3.1 summary.

  • Sound is a longitudinal pressure wave obeying the Chapter 1 wave equation; \(c = \sqrt{K/\rho} \approx 1540\) m/s in soft tissue and \(c = f\lambda\).
  • Impedance \(Z = \rho c\) links pressure to particle motion through \(p = Zu\), and intensity follows as \(I = p_A^2/2\rho c\).
  • Echoes arise from impedance mismatches; \(R\) is signed, \(R_I\) is its square.
  • Attenuation is exponential and rises with frequency; nepers and decibels differ by 4.343 (intensity) or 8.686 (pressure), and \(\mu_I = 2\alpha\).
  • Rayleigh scattering (\(\propto f^4\)) off blood enables Doppler; coherent scattering creates speckle, whose envelope SNR is fixed at 1.91 and therefore is not noise.
  • Three resolutions: axial (pulse length), lateral (aperture), elevational (fixed lens, and the source of slice-thickness artifact).

Checkpoint 3.1. A vascular study needs 0.2 mm axial resolution at 2 cm depth; an abdominal study needs to reach 15 cm. Using Eqs. (12) and (9), choose a frequency for each, and explain quantitatively why no single probe serves both well.

3.2 Instrumentation and Signal Processing

3.2.1 Transducers: the piezoelectric effect

An ultrasound transducer both transmits and receives sound using the piezoelectric effect: certain materials (historically the ceramic PZT, increasingly single crystals and polymer composites) develop a surface charge when mechanically strained and, conversely, deform when a voltage is applied. A voltage pulse makes the element vibrate and emit sound; a returning echo strains the element and generates a measurable voltage. The element is most efficient at its thickness resonance, where the thickness equals half an acoustic wavelength in the piezoelectric material, \(t = \lambda_{\text{PZT}}/2\).

Two construction details determine image quality.

A matching layer on the front face reduces the huge reflection at the transducer–tissue impedance step (PZT is about 30 MRayl, tissue about 1.6 MRayl). A quarter-wave layer of intermediate impedance

\[\begin{equation} Z_m = \sqrt{Z_{\text{PZT}}\,Z_{\text{tissue}}}, \qquad t_m = \frac{\lambda_m}{4}, \tag{14} \end{equation}\]

minimizes that reflection, the acoustic analogue of an anti-reflection optical coating, and the same quarter-wave condition met in Chapter 2.

A backing (damping) layer behind the element absorbs rearward energy and shortens the ringdown, producing a short, broadband pulse. Short pulses mean good axial resolution (Eq. (12)), so there is a deliberate trade between sensitivity (lightly damped, narrowband, long pulse) and resolution (heavily damped, broadband, short pulse).

Z_pzt <- 30; Z_tis <- 1.6
Zm_opt <- sqrt(Z_pzt*Z_tis)

## Reflection from a single layer of impedance Zm, thickness = quarter wave at f0.
## Transmission-line result at the design frequency: R = ((Zm^2 - Z1 Z2)/(Zm^2 + Z1 Z2))^2
Zm <- seq(1, 20, length.out = 500)
R_layer <- ((Zm^2 - Z_pzt*Z_tis)/(Zm^2 + Z_pzt*Z_tis))^2
R_none  <- ((Z_tis - Z_pzt)/(Z_tis + Z_pzt))^2

p_m1 <- ggplot(data.frame(Zm = Zm, R = R_layer), aes(Zm, R)) +
  geom_hline(yintercept = R_none, linetype = "dashed", colour = bpad_pal[2]) +
  geom_vline(xintercept = Zm_opt, linetype = "dotted", colour = "grey40") +
  geom_line(linewidth = 0.95, colour = bpad_pal[1]) +
  annotate("text", x = 12, y = R_none + 0.04, size = 3, colour = bpad_pal[2],
           label = sprintf("no matching layer: %.0f%%", 100*R_none)) +
  annotate("text", x = Zm_opt + 0.4, y = 0.45, hjust = 0, size = 3, colour = "grey30",
           label = sprintf("optimum %.2f MRayl", Zm_opt)) +
  labs(x = "Matching-layer impedance (MRayl)", y = "Intensity reflected",
       title = "Choosing the matching-layer impedance")

## Frequency response of the optimal layer (quarter wave only at f0)
fr <- seq(0.3, 1.7, by = 0.005)          # normalized frequency f/f0
theta <- (pi/2)*fr                        # electrical thickness
Zin <- Zm_opt*(Z_tis + 1i*Zm_opt*tan(theta))/(Zm_opt + 1i*Z_tis*tan(theta))
R_f <- Mod((Zin - Z_pzt)/(Zin + Z_pzt))^2
p_m2 <- ggplot(data.frame(fr = fr, R = R_f), aes(fr, R)) +
  geom_hline(yintercept = R_none, linetype = "dashed", colour = bpad_pal[2]) +
  geom_line(linewidth = 0.95, colour = bpad_pal[3]) +
  labs(x = expression("Normalized frequency "*f/f[0]),
       y = "Intensity reflected",
       title = "Bandwidth of a single matching layer",
       subtitle = "perfect only at f0, but useful across the band")
bpad_grid(p_m1, p_m2, ncol = 2)
The quarter-wave matching layer, Eq. (matching-layer). Without matching, a PZT-tissue interface reflects more than 80 percent of the intensity back into the element. A single quarter-wave layer at the geometric-mean impedance eliminates the reflection at the design frequency, and the improvement remains substantial across the useful bandwidth. The dotted line marks the geometric-mean impedance that achieves the minimum.

Figure 4: The quarter-wave matching layer, Eq. (matching-layer). Without matching, a PZT-tissue interface reflects more than 80 percent of the intensity back into the element. A single quarter-wave layer at the geometric-mean impedance eliminates the reflection at the design frequency, and the improvement remains substantial across the useful bandwidth. The dotted line marks the geometric-mean impedance that achieves the minimum.

cat(sprintf("Unmatched PZT/tissue interface reflects %.1f%% of the intensity.\n",
            100*R_none))
## Unmatched PZT/tissue interface reflects 80.8% of the intensity.
cat(sprintf("Quarter-wave layer at Zm = sqrt(%.0f x %.1f) = %.2f MRayl reflects %.2e.\n",
            Z_pzt, Z_tis, Zm_opt, min(R_layer)))
## Quarter-wave layer at Zm = sqrt(30 x 1.6) = 6.93 MRayl reflects 2.84e-06.
cat(sprintf("At 5 MHz with c_m = 2000 m/s the layer thickness is %.1f um.\n",
            2000/(5e6)/4*1e6))
## At 5 MHz with c_m = 2000 m/s the layer thickness is 100.0 um.

3.2.2 Array geometries and beamforming

Modern probes are arrays of many small elements (often 128–512). Linear arrays fire groups of elements to sweep a rectangular field for vascular and small-parts imaging. Convex (curvilinear) arrays spread the field into a wider sector for abdominal and obstetric imaging. Phased arrays use the whole small aperture at once and steer the beam electronically, ideal for cardiac imaging through the narrow window between ribs.

Beamforming shapes and steers the beam by controlling the relative timing of the elements. On transmit, a curved profile of time delays across the aperture focuses the wavefronts to a chosen depth; adding a linear ramp of delays steers the beam. On receive, echoes from each element are delayed and summed (delay-and-sum), and because echoes from greater depths arrive later, the receive focus can be swept outward in step with depth, dynamic receive focusing, keeping the beam tight over the whole image. Apodization (tapering element weights) trades main-lobe width against side-lobe level.

Element spacing (pitch \(d\)) matters. Grating lobes appear at angles satisfying \(\sin\theta_g = \pm m\lambda/d\). For a beam steered to a maximum angle \(\theta_{\max}\), avoiding them over the whole scan requires

\[\begin{equation} d \le \frac{\lambda}{1 + |\sin\theta_{\max}|}, \tag{15} \end{equation}\]

which reduces to the familiar \(d \le \lambda/2\) for a phased array steered to \(\pm 90^\circ\), and permits a coarser pitch for an unsteered linear array.

Connection to Chapter 1. The far-field beam pattern of an aperture is the Fourier transform of its aperture function, the same transform used for diffraction in Chapter 1. A uniform \(N\)-element array therefore has the Dirichlet (periodic-sinc) beam pattern below, whose main-lobe width scales as \(\lambda/\text{aperture}\) and whose grating lobes appear according to Eq. (15). Delay-and-sum beamforming and Fourier-domain reconstruction are two views of the same linear-systems picture.

array_factor <- function(theta_deg, N, d_over_lambda, w = rep(1, N)) {
  th  <- theta_deg*pi/180
  n   <- 0:(N - 1)
  vapply(th, function(t) {
    ph <- 2*pi*d_over_lambda*n*sin(t)
    20*log10(pmax(Mod(sum(w*exp(1i*ph)))/sum(w), 1e-5))
  }, numeric(1))
}
theta <- seq(-90, 90, by = 0.2); N_el <- 16
w_hann <- 0.5 - 0.5*cos(2*pi*(0:(N_el - 1))/(N_el - 1)); w_hann[w_hann == 0] <- 1e-6

mkbeam <- function(dl, w, ttl, sub, col) {
  ggplot(data.frame(theta = theta, dB = array_factor(theta, N_el, dl, w)),
         aes(theta, dB)) +
    geom_line(colour = col, linewidth = 0.8) +
    coord_cartesian(ylim = c(-50, 2)) +
    labs(x = "Angle (degrees)", y = "Beam response (dB)", title = ttl, subtitle = sub)
}
bpad_grid(
  mkbeam(0.5, rep(1, N_el), "Pitch = lambda/2", "clean main lobe", bpad_pal[1]),
  mkbeam(1.0, rep(1, N_el), "Pitch = lambda",   "grating lobes appear", bpad_pal[2]),
  mkbeam(0.5, w_hann,       "Hann apodization", "side lobes suppressed", bpad_pal[3]),
  ncol = 3)
Array beam patterns, the Fourier transform of the aperture. Left: with half-wavelength pitch only a main lobe and side lobes appear. Centre: with full-wavelength pitch, grating lobes appear near plus and minus 90 degrees and would place echoes in entirely the wrong direction. Right: apodization (Hann weighting of the element amplitudes) suppresses side lobes by more than 20 dB at the cost of a wider main lobe, which is the standard contrast-versus-resolution trade in beamforming.

Figure 5: Array beam patterns, the Fourier transform of the aperture. Left: with half-wavelength pitch only a main lobe and side lobes appear. Centre: with full-wavelength pitch, grating lobes appear near plus and minus 90 degrees and would place echoes in entirely the wrong direction. Right: apodization (Hann weighting of the element amplitudes) suppresses side lobes by more than 20 dB at the cost of a wider main lobe, which is the standard contrast-versus-resolution trade in beamforming.

for (smax in c(0, 30, 45, 90))
  cat(sprintf("steering to +/-%2d deg -> max pitch = %.3f lambda\n",
              smax, 1/(1 + abs(sin(smax*pi/180)))))
## steering to +/- 0 deg -> max pitch = 1.000 lambda
## steering to +/-30 deg -> max pitch = 0.667 lambda
## steering to +/-45 deg -> max pitch = 0.586 lambda
## steering to +/-90 deg -> max pitch = 0.500 lambda

3.2.3 From echo to grey level: TGC and log compression

Two processing steps stand between the raw echo and the displayed pixel, and without them no image would be visible at all: Section 3.1.4 showed that a 5 MHz beam loses 60 dB on a 12 cm round trip.

Time-gain compensation (TGC) applies a depth-increasing amplification that approximately cancels \(e^{-\mu_I x}\), so identical tissue looks equally bright at all depths. The user-adjustable TGC sliders are literally a piecewise-linear reconstruction of the inverse attenuation curve. TGC cannot create information: it amplifies noise equally, so beyond the depth at which echo falls below the noise floor it produces only a bright, uniform, uninformative haze.

Logarithmic compression maps the enormous echo dynamic range, 100–120 dB between a specular bone reflection and weak parenchymal scatter, onto the roughly 8-bit range of a display and of human perception:

\[\begin{equation} G = A\log_{10}\!\left(\frac{E}{E_{\text{ref}}}\right) + B . \tag{16} \end{equation}\]

The operator’s “dynamic range” control sets \(A\): a low value gives a high-contrast image with few grey levels, good for detecting boundaries; a high value gives a soft image that preserves subtle parenchymal texture.

z <- seq(0, 18, by = 0.05)
f_tgc <- 3.5
echo_dB  <- -a_dB*f_tgc*2*z
noise_dB <- -80
raw      <- pmax(echo_dB, noise_dB)
tgc_gain <- a_dB*f_tgc*2*z
after    <- pmin(raw + tgc_gain, 6)
z_floor  <- -noise_dB/(a_dB*f_tgc*2)

df_t <- rbind(
  data.frame(z = z, dB = raw,   trace = "raw echo"),
  data.frame(z = z, dB = after, trace = "after TGC"),
  data.frame(z = z, dB = rep(noise_dB, length(z)), trace = "noise floor"))
p_tgc <- ggplot(df_t, aes(z, dB, colour = trace)) +
  geom_vline(xintercept = z_floor, linetype = "dotted", colour = "grey45") +
  geom_line(linewidth = 0.9) +
  annotate("text", x = z_floor + 0.3, y = -30, hjust = 0, size = 2.9, colour = "grey30",
           label = sprintf("echo reaches\nnoise floor\n%.1f cm", z_floor)) +
  scale_colour_manual(values = c("raw echo" = bpad_pal[1], "after TGC" = bpad_pal[2],
                                 "noise floor" = bpad_pal[7])) +
  labs(x = "Depth (cm)", y = "Level (dB)", title = "Time-gain compensation",
       subtitle = "3.5 MHz round trip")

E_dB <- seq(-100, 0, by = 0.5)
dr <- c(40, 60, 80)
lc <- do.call(rbind, lapply(dr, function(d)
  data.frame(E = E_dB, grey = pmin(pmax(255*(E_dB + d)/d, 0), 255),
             DR = factor(paste0(d, " dB"), levels = paste0(dr, " dB")))))
p_lc <- ggplot(lc, aes(E, grey, colour = DR)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(x = "Echo level (dB relative to maximum)", y = "Displayed grey level (0-255)",
       title = "Logarithmic compression")
bpad_grid(p_tgc, p_lc, ncol = 2)
Time-gain compensation and logarithmic compression. Left: raw echo from uniform tissue falls exponentially and reaches the noise floor near 11 cm; TGC restores uniform displayed brightness up to that depth and thereafter amplifies only noise, which is the physical origin of the bright, featureless haze at the bottom of a deep image. Right: log compression maps a 100 dB echo range onto the displayable grey scale, with the operator's dynamic-range control setting the slope.

Figure 6: Time-gain compensation and logarithmic compression. Left: raw echo from uniform tissue falls exponentially and reaches the noise floor near 11 cm; TGC restores uniform displayed brightness up to that depth and thereafter amplifies only noise, which is the physical origin of the bright, featureless haze at the bottom of a deep image. Right: log compression maps a 100 dB echo range onto the displayable grey scale, with the operator’s dynamic-range control setting the slope.

3.2.4 The frame-rate budget

Each scan line must wait for the deepest echo before the next pulse is sent. With \(N\) lines to depth \(d\), each line costs \(2d/c\) seconds, so

\[\begin{equation} \mathrm{FR} \le \frac{c}{2\,d\,N_{\text{lines}}\,N_{\text{foci}}} . \tag{17} \end{equation}\]

This one inequality governs almost every real-time compromise in ultrasound: sector width, depth setting, line density, the number of transmit focal zones (each multiplies the line count), and colour-box size. Plane-wave / ultrafast imaging breaks the budget by insonating the whole field with a single unfocused transmit and forming all lines in parallel on receive, reaching thousands of frames per second, which is precisely what made shear-wave elastography (Section 3.5.3) and ultrasound localization microscopy possible.

d_cm <- seq(2, 20, by = 0.1)
lines_n <- c(64, 128, 192, 256)
fr_df <- do.call(rbind, lapply(lines_n, function(N)
  data.frame(d = d_cm, FR = c_soft/(2*(d_cm/100)*N),
             lines = factor(paste0(N, " lines"), levels = paste0(lines_n, " lines")))))

ggplot(fr_df, aes(d, FR, colour = lines)) +
  annotate("rect", xmin = 2, xmax = 20, ymin = 25, ymax = 30,
           fill = "grey50", alpha = 0.20) +
  geom_line(linewidth = 0.9) +
  annotate("text", x = 19.5, y = 36, hjust = 1, size = 3, colour = "grey30",
           label = "real-time threshold") +
  scale_y_log10() + scale_colour_manual(values = bpad_pal[1:4]) +
  labs(x = "Imaging depth (cm)", y = "Maximum frame rate (Hz, log scale)",
       title = "Depth times line count sets the frame rate")
The frame-rate budget, Eq. (frame-rate). Achievable frame rate falls as the product of depth and line count; the shaded band marks the 25-30 Hz needed for flicker-free real-time display. Cardiac imaging survives by using relatively few lines over a narrow sector; deep obstetric imaging with high line density is inherently slow, and each additional transmit focal zone divides these numbers again.

Figure 7: The frame-rate budget, Eq. (frame-rate). Achievable frame rate falls as the product of depth and line count; the shaded band marks the 25-30 Hz needed for flicker-free real-time display. Cardiac imaging survives by using relatively few lines over a narrow sector; deep obstetric imaging with high line density is inherently slow, and each additional transmit focal zone divides these numbers again.

knitr::kable(data.frame(
  exam  = c("Cardiac (sector)", "Vascular", "Abdominal", "Obstetric (deep)"),
  depth = c(16, 4, 15, 15), lines = c(96, 128, 192, 256), foci = c(1, 2, 1, 1),
  FR = round(c_soft/(2*c(16,4,15,15)/100*c(96,128,192,256)*c(1,2,1,1)), 1)),
  col.names = c("Examination","Depth (cm)","Lines","Focal zones","Max frame rate (Hz)"),
  caption = "Frame-rate ceilings for representative examinations. The vascular study uses two focal zones, which halves its otherwise very comfortable budget.")
Table 4: Frame-rate ceilings for representative examinations. The vascular study uses two focal zones, which halves its otherwise very comfortable budget.
Examination Depth (cm) Lines Focal zones Max frame rate (Hz)
Cardiac (sector) 16 96 1 50.1
Vascular 4 128 2 75.2
Abdominal 15 192 1 26.7
Obstetric (deep) 15 256 1 20.1

Section 3.2 summary.

  • Piezoelectric elements interconvert voltage and pressure; thickness resonance sets the centre frequency.
  • A quarter-wave matching layer at \(Z_m=\sqrt{Z_1Z_2}\) removes an 82% reflection; backing shortens the pulse for axial resolution at the cost of sensitivity.
  • Beamforming is controlled element delay: focusing and steering on transmit, dynamic delay-and-sum on receive; apodization trades main-lobe width for side-lobe level.
  • Grating lobes are avoided when \(d\le\lambda/(1+|\sin\theta_{\max}|)\).
  • TGC undoes attenuation until the echo reaches the noise floor; log compression fits 100 dB onto an 8-bit display.
  • Frame rate is bounded by \(c/(2dN)\); ultrafast plane-wave imaging escapes the bound.

Checkpoint 3.2. A cardiologist widens the sector from 60° to 90° and adds a second focal zone, keeping line density and depth constant. Using Eq. (17), estimate the factor by which the frame rate falls, and state whether the study will still support real-time wall-motion assessment.

3.3 Conventional Imaging Modes

3.3.1 A-, B-, and M-mode

All grayscale ultrasound rests on the pulse-echo range equation. The scanner emits a pulse, then times the returning echo. Since the pulse travels to the reflector and back,

\[\begin{equation} d = \frac{c\,t}{2}, \tag{18} \end{equation}\]

with \(c = 1540\) m/s assumed. A-mode (amplitude) plots echo amplitude versus depth along a single line, still used in ophthalmology for precise axial distances. B-mode (brightness) maps echo amplitude to pixel brightness and assembles hundreds of adjacent scan lines into the familiar two-dimensional grayscale image. M-mode (motion) repeatedly samples a single line and displays depth versus time; sacrificing the second spatial dimension buys very high temporal resolution, ideal for valve motion in echocardiography.

From the equation to the clinic. Equation (18) is why an echo at \(t = 130\ \mu\)s corresponds to \(d = 1540\times130\times10^{-6}/2 = 10\) cm. The same relation bounds the frame rate (Eq. (17)) and, as Section 3.3.3 shows, bounds the measurable blood velocity.

3.3.2 Doppler ultrasound

When sound reflects off a moving target its frequency shifts. For blood moving with velocity \(v\) at angle \(\theta\) to the beam, the round-trip shift is

\[\begin{equation} f_D = \frac{2 f_0\, v \cos\theta}{c}, \qquad\Longrightarrow\qquad v = \frac{c\, f_D}{2 f_0 \cos\theta} . \tag{19} \end{equation}\]

The factor of two arises because the blood is both a moving receiver and a moving re-emitter. (Equation (19) is the first-order form; the exact expression carries a further factor \(1/(1 - v\cos\theta/c)\), which for \(v/c \sim 10^{-3}\) changes the result by less than 0.1% and is universally neglected.)

The \(\cos\theta\) term makes the measurement angle-dependent, and differentiating Eq. (19) shows exactly how dangerous that is:

\[\begin{equation} \frac{\Delta v}{v} = \tan\theta\;\Delta\theta . \tag{20} \end{equation}\]

At \(\theta = 20^\circ\) a \(5^\circ\) angle error costs 3% in velocity; at \(\theta = 70^\circ\) the same error costs 24%; at \(\theta = 80^\circ\) it costs 50%. This is the quantitative basis of the clinical rule to keep \(\theta \le 60^\circ\).

Three implementations trade range information against velocity range:

  • Continuous-wave (CW) Doppler uses separate transmit and receive elements, measures very high velocities without ambiguity, but cannot tell where along the beam the flow occurred.
  • Pulsed-wave (PW) Doppler sends pulses at a pulse-repetition frequency (PRF) and range-gates the echoes, localizing the flow, but samples the Doppler signal at the PRF and is therefore Nyquist-limited.
  • Colour-flow imaging estimates the mean Doppler shift over a region (typically by autocorrelation) and overlays direction and speed on the B-mode image.
f0_d <- 5e6
v_g <- seq(0, 1.5, by = 0.01)
angles <- c(0, 30, 60, 80)
dop <- do.call(rbind, lapply(angles, function(a)
  data.frame(v = v_g, fD = 2*f0_d*v_g*cos(a*pi/180)/c_soft/1e3,
             angle = factor(paste0(a, " deg"), levels = paste0(angles, " deg")))))
p_d1 <- ggplot(dop, aes(v, fD, colour = angle)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = bpad_pal[1:4]) +
  labs(x = "Blood velocity (m/s)", y = "Doppler shift (kHz)",
       title = "Doppler shift versus velocity and angle",
       subtitle = "f0 = 5 MHz, c = 1540 m/s")

th_g <- seq(0, 85, by = 0.5)
err  <- 100*tan(th_g*pi/180)*(5*pi/180)
p_d2 <- ggplot(data.frame(theta = th_g, err = err), aes(theta, err)) +
  geom_vline(xintercept = 60, linetype = "dashed", colour = bpad_pal[2]) +
  geom_line(linewidth = 0.95, colour = bpad_pal[1]) +
  annotate("text", x = 58, y = 60, hjust = 1, size = 3, colour = bpad_pal[2],
           label = "clinical limit 60 deg") +
  coord_cartesian(ylim = c(0, 100)) +
  labs(x = "Beam-to-flow angle (degrees)", y = "Velocity error (%) from a 5 deg angle error",
       title = "Why the 60-degree rule exists")
bpad_grid(p_d1, p_d2, ncol = 2)
The Doppler equation and its angle sensitivity. Left: Doppler shift grows with blood velocity and collapses toward zero as the beam-to-flow angle approaches 90 degrees. Right: the fractional velocity error produced by a 5 degree angle-estimation error, Eq. (doppler-angle-error); the error is tolerable below 60 degrees and diverges beyond it, which is the entire justification for the clinical 60-degree rule.

Figure 8: The Doppler equation and its angle sensitivity. Left: Doppler shift grows with blood velocity and collapses toward zero as the beam-to-flow angle approaches 90 degrees. Right: the fractional velocity error produced by a 5 degree angle-estimation error, Eq. (doppler-angle-error); the error is tolerable below 60 degrees and diverges beyond it, which is the entire justification for the clinical 60-degree rule.

knitr::kable(data.frame(
  angle = c(0, 20, 45, 60, 70, 80),
  cos_theta = round(cos(c(0,20,45,60,70,80)*pi/180), 3),
  err_5deg = sprintf("%.1f%%", 100*tan(c(0,20,45,60,70,80)*pi/180)*(5*pi/180))),
  col.names = c("Angle (deg)","cos(theta)","Velocity error from 5 deg angle error"),
  caption = "Angle-error propagation in Doppler velocimetry, Eq. (doppler-angle-error).")
Table 5: Angle-error propagation in Doppler velocimetry, Eq. (doppler-angle-error).
Angle (deg) cos(theta) Velocity error from 5 deg angle error
0 1.000 0.0%
20 0.940 3.2%
45 0.707 8.7%
60 0.500 15.1%
70 0.342 24.0%
80 0.174 49.5%

3.3.3 The depth–velocity limit of pulsed Doppler

Pulsed-wave Doppler samples the flow signal once per pulse, at rate PRF. By the Nyquist criterion of Chapter 1 the largest unambiguous Doppler shift is \(\mathrm{PRF}/2\), so

\[\begin{equation} v_{\max} = \frac{c\,\mathrm{PRF}}{4 f_0 \cos\theta} . \tag{21} \end{equation}\]

Exceed it and the spectrum wraps around, aliasing, displaying fast forward flow as reverse flow. But Eq. (21) alone makes the PRF look like a free parameter, and it is not. The scanner must wait for the deepest echo of interest before firing again:

\[\begin{equation} \mathrm{PRF} \le \frac{c}{2\,d_{\max}} . \tag{22} \end{equation}\]

Substituting Eq. (22) into Eq. (21) gives a hard, instrument-independent bound:

\[\begin{equation} v_{\max}\, d_{\max} \le \frac{c^{2}}{8 f_0 \cos\theta} . \tag{23} \end{equation}\]

Equation (23) says that depth and velocity trade against one another on a hyperbola whose position is fixed by the transmit frequency alone. Every remedy for aliasing is a move along or across this curve: raising the PRF is possible only if you accept a shallower gate; lowering \(f_0\) moves the whole hyperbola outward but costs resolution; increasing \(\theta\) helps through \(\cos\theta\) but inflates the velocity error through Eq. (20).

d_grid <- seq(0.01, 0.22, length.out = 400)
f_list <- c(2.0e6, 3.5e6, 5.0e6, 7.5e6)
hyp <- do.call(rbind, lapply(f_list, function(f0)
  data.frame(d_cm = d_grid*100, v = c_soft^2/(8*f0*d_grid),
             probe = factor(sprintf("%.1f MHz", f0/1e6),
                            levels = sprintf("%.1f MHz", f_list/1e6)))))
clinical <- data.frame(
  task = c("Carotid, normal", "Carotid stenosis", "Portal vein", "Aortic stenosis jet"),
  d_cm = c(2.5, 2.5, 8, 12), v = c(0.8, 2.5, 0.2, 4.0))

ggplot(hyp, aes(d_cm, v, colour = probe)) +
  geom_line(linewidth = 0.9) +
  geom_point(data = clinical, aes(d_cm, v), inherit.aes = FALSE,
             size = 2.6, colour = "grey15") +
  geom_text(data = clinical, aes(d_cm, v, label = task), inherit.aes = FALSE,
            vjust = -0.9, size = 3, colour = "grey20") +
  scale_y_log10() + scale_colour_manual(values = bpad_pal[1:4]) +
  coord_cartesian(xlim = c(0, 22), ylim = c(0.05, 8)) +
  labs(x = "Maximum imaging depth (cm)", y = "Nyquist velocity (m/s, log scale)",
       title = "Pulsed-wave Doppler cannot have both depth and velocity",
       subtitle = "curves: v_max * d_max = c^2 / (8 f0 cos(theta)), theta = 0")
The pulsed-Doppler depth-velocity limit, Eq. (depth-velocity). Each curve is the hyperbola for one transmit frequency; the accessible region for a given probe lies below and to the left of its curve. Points mark four clinical measurements. Aortic stenosis at 4 m/s and 12 cm depth sits far outside every pulsed-Doppler curve, which is precisely why peak stenotic jets must be measured with continuous-wave Doppler, which never samples and therefore never aliases.

Figure 9: The pulsed-Doppler depth-velocity limit, Eq. (depth-velocity). Each curve is the hyperbola for one transmit frequency; the accessible region for a given probe lies below and to the left of its curve. Points mark four clinical measurements. Aortic stenosis at 4 m/s and 12 cm depth sits far outside every pulsed-Doppler curve, which is precisely why peak stenotic jets must be measured with continuous-wave Doppler, which never samples and therefore never aliases.

for (f0 in c(2e6, 5e6))
  cat(sprintf("f0 = %.1f MHz: budget %.4f m^2/s -> at 12 cm, v_max = %.2f m/s\n",
              f0/1e6, c_soft^2/(8*f0), c_soft^2/(8*f0*0.12)))
## f0 = 2.0 MHz: budget 0.1482 m^2/s -> at 12 cm, v_max = 1.24 m/s
## f0 = 5.0 MHz: budget 0.0593 m^2/s -> at 12 cm, v_max = 0.49 m/s
cat(sprintf("A 4 m/s jet at 12 cm demands %.2f m^2/s, which is %.1fx the 2 MHz budget.\n",
            4*0.12, 4*0.12/(c_soft^2/(8*2e6))))
## A 4 m/s jet at 12 cm demands 0.48 m^2/s, which is 3.2x the 2 MHz budget.
cat("No PRF setting reaches it: continuous-wave Doppler is physically required.\n")
## No PRF setting reaches it: continuous-wave Doppler is physically required.

From the equation to the clinic. This resolves a puzzle that the Nyquist relation alone leaves open. Peak aortic-stenosis velocity is always measured with CW Doppler “because it does not alias”, but why not simply raise the PRF? Because at 12 cm depth even a 2 MHz probe caps out near 1.2 m/s, and a stenotic jet is three to four times that. The choice of CW is forced by physics, not convention, and the cost paid is the loss of range information: CW reports the highest velocity anywhere along the beam.

Section 3.3 summary.

  • \(d = ct/2\) underlies A-, B-, and M-mode.
  • \(f_D = 2f_0 v\cos\theta/c\) measures velocity; angle error propagates as \(\Delta v/v=\tan\theta\,\Delta\theta\), which is why \(\theta\le60^\circ\).
  • PW Doppler is Nyquist-limited at \(v_{\max}=c\,\mathrm{PRF}/4f_0\cos\theta\), and the PRF is itself capped by depth.
  • Combining them gives \(v_{\max}d_{\max}\le c^2/8f_0\cos\theta\): the fundamental depth-velocity trade of pulsed Doppler, and the reason CW exists.

Checkpoint 3.3. A PW study of a stenotic jet aliases. Using Eq. (23), name three settings you could change to raise the Nyquist velocity, state the physical cost of each, and identify the one clinical situation in which none of them suffices.

3.4 Artifacts: When the Assumptions Fail

Every ultrasound artifact is the violation of one of four assumptions the scanner makes: that sound travels at exactly 1540 m/s, in a straight line, that all echoes come from the main beam, and that each echo returns after a single reflection. Artifacts are therefore not defects to be memorized but applied physics, and two of them are diagnostically useful. Table 6 organizes them by the assumption each violates.

Table 6: Ultrasound artifacts organized by the assumption each violates.
Artifact Assumption violated Physical mechanism Clinical reading
Acoustic shadowing , (energy conservation) Large \(Z\) mismatch or high absorption removes energy Calculus, bone, gas
Posterior enhancement uniform attenuation Beam through low-attenuation fluid retains energy Confirms a lesion is cystic
Reverberation single reflection Repeated bouncing between two strong reflectors, each round trip displayed deeper Equally spaced lines; comet-tail / ring-down
Mirror image single reflection Strong specular reflector (diaphragm) creates a virtual duplicate beyond it False “lesion” above the diaphragm
Refraction straight-line travel Snell’s law bends the beam at an oblique interface Lateral displacement; edge shadowing
Speed-of-sound (range) error \(c = 1540\) m/s Fat is slower, so structures display too deep Measurement error through thick fat
Side-lobe / grating lobe echoes from the main beam only Off-axis energy assigned to the main-beam line Spurious echo inside anechoic structures
Slice thickness zero-thickness image plane Elevational beam width averages out-of-plane echoes Pseudo-debris in small cysts
## Speed-of-sound range error through a fat layer
fat_cm <- seq(0, 6, by = 0.05)
c_fat  <- 1450
true_target <- 10                                  # cm, true depth of the target
t_total <- (fat_cm/100)/c_fat + ((true_target - fat_cm)/100)/c_soft
app_depth <- c_soft*t_total*100
p_a1 <- ggplot(data.frame(fat = fat_cm, err = app_depth - true_target),
               aes(fat, err*10)) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_line(linewidth = 0.95, colour = bpad_pal[1]) +
  labs(x = "Thickness of intervening fat (cm)",
       y = "Displayed depth error (mm)",
       title = "Speed-of-sound range error",
       subtitle = "target truly at 10 cm; fat at 1450 m/s vs assumed 1540 m/s")

## Reverberation ladder
true_d <- 3                                        # cm to the strong reflector
orders <- 1:6
R_strong <- 0.6                                    # amplitude reflection at the pair
rev_df <- data.frame(order = orders, depth = orders*true_d,
                     amp = R_strong^(orders - 1))
p_a2 <- ggplot(rev_df, aes(depth, amp)) +
  geom_segment(aes(xend = depth, yend = 0), linewidth = 1.3, colour = bpad_pal[2]) +
  geom_point(size = 2.4, colour = bpad_pal[2]) +
  geom_text(aes(label = ifelse(order == 1, "true", paste0("x", order))),
            vjust = -0.8, size = 3) +
  coord_cartesian(ylim = c(0, 1.2)) +
  labs(x = "Displayed depth (cm)", y = "Echo amplitude (a.u.)",
       title = "Reverberation ladder",
       subtitle = "reflector truly at 3 cm; copies at integer multiples")
bpad_grid(p_a1, p_a2, ncol = 2)
Two artifacts computed rather than asserted. Left: the speed-of-sound range error, Eq. (pulse-echo), for a structure viewed through a layer of fat. Because fat propagates at 1450 m/s while the scanner assumes 1540 m/s, every structure beyond the fat is displayed too deep, and the error grows linearly with fat thickness. Right: reverberation between the transducer face and a strong specular reflector; each successive round trip is displayed at an integer multiple of the true depth with progressively lower amplitude, producing the characteristic ladder of equally spaced lines.

Figure 10: Two artifacts computed rather than asserted. Left: the speed-of-sound range error, Eq. (pulse-echo), for a structure viewed through a layer of fat. Because fat propagates at 1450 m/s while the scanner assumes 1540 m/s, every structure beyond the fat is displayed too deep, and the error grows linearly with fat thickness. Right: reverberation between the transducer face and a strong specular reflector; each successive round trip is displayed at an integer multiple of the true depth with progressively lower amplitude, producing the characteristic ladder of equally spaced lines.

cat(sprintf("Through 4 cm of fat, a target at 10 cm is displayed at %.2f cm (%.1f mm too deep).\n",
            c_soft*((4/100)/c_fat + ((true_target-4)/100)/c_soft)*100,
            10*(c_soft*((4/100)/c_fat + ((true_target-4)/100)/c_soft)*100 - true_target)))
## Through 4 cm of fat, a target at 10 cm is displayed at 10.25 cm (2.5 mm too deep).

Two artifacts that are measurements. Posterior acoustic enhancement, increased brightness deep to a structure, arises because the beam lost less energy crossing fluid than crossing the surrounding tissue that the TGC curve was set for. It therefore confirms that a lesion is fluid-filled. Acoustic shadowing deep to a small bright focus confirms a gallstone or calcification. An artifact you understand is a measurement; an artifact you do not understand is a misdiagnosis.

Checkpoint 3.4. A round anechoic structure in the liver shows bright echoes posterior to it and faint internal echoes near its upper margin. Name the two artifacts present, state which assumption each violates, and say what they jointly tell you about the lesion.

3.5 Advanced Ultrasound Techniques

3.5.1 Harmonic imaging

At diagnostic pressures, sound propagation through tissue is slightly nonlinear: compressions travel marginally faster than rarefactions, so an initially pure tone at \(f_0\) progressively distorts and generates energy at the harmonic \(2f_0\). Tissue harmonic imaging transmits at \(f_0\) but forms the image from the received \(2f_0\) component. Because the harmonic builds up only along the strong central beam and grows with depth, it suppresses the weak-signal clutter, reverberation, and side-lobe artifacts that plague the fundamental, markedly improving contrast, especially in technically difficult patients.

3.5.2 Contrast-enhanced ultrasound and the Minnaert resonance

Microbubbles, gas cores a few micrometers across, stabilized by a lipid or protein shell, are injected intravenously as a pure intravascular contrast agent. Driven by the ultrasound field, they oscillate; near resonance the oscillation is strongly nonlinear and asymmetric (they resist compression more than expansion), emitting harmonic and subharmonic signals that tissue does not. Contrast-specific pulse sequences isolate those bubble-only signals, allowing real-time imaging of perfusion.

But why do microbubbles resonate in the diagnostic band? A gas bubble in an incompressible liquid is a mass–spring oscillator: the gas provides the stiffness, the entrained liquid the inertia. Its resonance is the Minnaert frequency

\[\begin{equation} f_{\text{res}} = \frac{1}{2\pi R_0}\sqrt{\frac{3\gamma P_0}{\rho}}, \tag{24} \end{equation}\]

with \(R_0\) the equilibrium radius, \(\gamma\) the polytropic exponent of the gas, \(P_0\) the ambient pressure, and \(\rho\) the liquid density.

minnaert <- function(R, gamma = 1.4, P0 = 101325, rho = 1000)
  (1/(2*pi*R))*sqrt(3*gamma*P0/rho)
R <- seq(0.4, 6, length.out = 400)*1e-6

ggplot(data.frame(diam_um = 2*R*1e6, f_MHz = minnaert(R)/1e6), aes(diam_um, f_MHz)) +
  annotate("rect", xmin = 1, xmax = 5, ymin = 0.1, ymax = 20,
           fill = bpad_pal[1], alpha = 0.12) +
  annotate("rect", xmin = 0.5, xmax = 12, ymin = 1, ymax = 7,
           fill = bpad_pal[2], alpha = 0.12) +
  geom_line(linewidth = 1, colour = "grey15") +
  annotate("text", x = 2.6, y = 15, label = "capillary-passable size",
           size = 3, colour = bpad_pal[1]) +
  annotate("text", x = 9, y = 2.6, label = "diagnostic band",
           size = 3, colour = bpad_pal[2]) +
  scale_x_log10() + scale_y_log10() +
  coord_cartesian(xlim = c(0.8, 12), ylim = c(0.2, 20)) +
  labs(x = "Bubble diameter (um, log scale)",
       y = "Resonance frequency (MHz, log scale)",
       title = "Microbubble resonance falls inside the diagnostic band")
The Minnaert resonance, Eq. (minnaert). A microbubble contrast agent must be smaller than a red blood cell to cross the pulmonary capillary bed (vertical shaded band, 1 to 5 micrometers diameter), and bubbles of exactly that size resonate in the 1 to 7 MHz diagnostic band (horizontal shaded band). The overlap of the two bands is a numerical coincidence of gas stiffness, blood density, and capillary calibre, and it is the reason contrast-enhanced ultrasound is possible at all.

Figure 11: The Minnaert resonance, Eq. (minnaert). A microbubble contrast agent must be smaller than a red blood cell to cross the pulmonary capillary bed (vertical shaded band, 1 to 5 micrometers diameter), and bubbles of exactly that size resonate in the 1 to 7 MHz diagnostic band (horizontal shaded band). The overlap of the two bands is a numerical coincidence of gas stiffness, blood density, and capillary calibre, and it is the reason contrast-enhanced ultrasound is possible at all.

knitr::kable(data.frame(
  diameter_um = c(1, 2, 3, 4, 6),
  f_res_MHz = round(minnaert(c(1,2,3,4,6)/2*1e-6)/1e6, 2),
  passes_capillary = c("yes","yes","yes","yes","marginal")),
  col.names = c("Diameter (um)","Resonance (MHz)","Crosses lung capillaries?"),
  caption = "Minnaert resonance for clinically relevant bubble sizes.")
Table 7: Minnaert resonance for clinically relevant bubble sizes.
Diameter (um) Resonance (MHz) Crosses lung capillaries?
1 6.57 yes
2 3.28 yes
3 2.19 yes
4 1.64 yes
6 1.09 marginal

The enabling coincidence. A microbubble must be under about 6 µm to traverse the pulmonary capillaries and reach the systemic circulation. Equation (24) puts the resonance of a bubble that size at 1–7 MHz, exactly the diagnostic band. The agent is therefore driven at resonance by an ordinary clinical probe, where its oscillation is strongly nonlinear and its harmonic emission far exceeds that of tissue. Contrast-enhanced ultrasound exists because of a numerical accident.

Microbubbles are destroyed by the beam. Above MI \(\approx 0.4\) the bubbles undergo inertial collapse. This is exploited deliberately in flash-replenishment perfusion imaging, a high-MI pulse destroys all bubbles in the plane and the refill rate measures flow, but it means routine CEUS must be performed at low MI, typically below 0.2. Section 3.6 defines the index.

3.5.3 Shear-wave elastography

Elastography images tissue stiffness, a property palpation has always exploited but that grayscale ultrasound does not show. In shear-wave elastography (SWE), a focused “push” pulse displaces tissue at depth, launching transverse shear waves that travel outward. The push arises from acoustic radiation force, the momentum transferred by absorption:

\[\begin{equation} F = \frac{2\alpha I}{c}, \tag{25} \end{equation}\]

a body force per unit volume, with \(\alpha\) the pressure absorption coefficient (Np m\(^{-1}\)) and \(I\) the local time-averaged intensity. Ultrafast plane-wave imaging (Section 3.2.4) then tracks the shear-wave speed \(c_s\), which is set by the tissue shear modulus \(\mu\) and density \(\rho\):

\[\begin{equation} c_s = \sqrt{\frac{\mu}{\rho}} . \tag{26} \end{equation}\]

For nearly incompressible soft tissue (Poisson ratio \(\nu \approx 0.5\)), the clinically reported Young’s modulus is

\[\begin{equation} E = 2\mu(1+\nu) \approx 3\mu = 3\rho\, c_s^{2} . \tag{27} \end{equation}\]

Note the enormous separation of scales: \(c_s\) is 1–10 m/s while the longitudinal speed \(c\) is 1540 m/s, because the shear modulus of soft tissue (kPa) is six orders of magnitude below its bulk modulus (GPa). Soft tissue resists compression like water and resists shear like jelly, and it is the shear property that disease changes.

Self-contained mechanics note. Equations (26)(27) are the small piece of solid mechanics this chapter needs. Stiffness is captured by the shear modulus \(\mu\); soft tissue is essentially incompressible, so \(E \approx 3\mu\); and SWE measures \(c_s\) to recover \(E\). This is the mechanical analogue of measuring a wave speed to infer a material property — exactly what Eq. (3) did for the bulk modulus.

Three assumptions behind every stiffness number. Equation (27) assumes the tissue is (i) linearly elastic, (ii) isotropic, and (iii) homogeneous over the measurement region. Real tissue is viscoelastic, so \(c_s\) is frequency-dependent and the reported value is a group speed over the push bandwidth; skeletal muscle is strongly anisotropic, so stiffness depends on probe orientation relative to the fibres; and a measurement box spanning a vessel or a focal lesion averages inhomogeneous material. This is why device-specific thresholds do not transfer between manufacturers, and why liver stiffness must be measured in a standardized intercostal position during a breath-hold.

rho_l <- 1000
cs <- seq(0.5, 4, by = 0.01)
df_e <- data.frame(cs = cs, E = 3*rho_l*cs^2/1000)
bands <- data.frame(
  ymin  = c(0, 7, 10, 14), ymax = c(7, 10, 14, 50),
  stage = factor(c("F0-F1 (none/mild)","F2 (significant)","F3 (severe)","F4 (cirrhosis)"),
                 levels = c("F0-F1 (none/mild)","F2 (significant)",
                            "F3 (severe)","F4 (cirrhosis)")))

ggplot() +
  geom_rect(data = bands, aes(xmin = 0.5, xmax = 4, ymin = ymin, ymax = ymax,
                              fill = stage), alpha = 0.25) +
  geom_line(data = df_e, aes(cs, E), linewidth = 1, colour = "grey15") +
  annotate("segment", x = 2.0, xend = 2.2, y = 3*rho_l*2.0^2/1000,
           yend = 3*rho_l*2.2^2/1000, colour = bpad_pal[2], linewidth = 1.1,
           arrow = arrow(length = unit(0.15, "cm"))) +
  annotate("text", x = 2.3, y = 11, hjust = 0, size = 3, colour = bpad_pal[2],
           label = "+10% in cs\n-> +21% in E") +
  scale_fill_manual(values = c(bpad_pal[3], bpad_pal[5], bpad_pal[2], bpad_pal[4])) +
  coord_cartesian(ylim = c(0, 50)) +
  labs(x = "Shear-wave speed cs (m/s)", y = "Young's modulus E (kPa)",
       title = "From shear-wave speed to tissue stiffness")
Shear-wave elastography converts a measured wave speed into Young's modulus through E = 3 rho cs^2, Eq. (youngs-modulus). Shaded bands show representative liver-stiffness ranges used to stage fibrosis. Because the relation is quadratic, measurement error in cs is doubled in E: the annotated example shows that a 10 percent speed error becomes a 21 percent stiffness error, which is why several acquisitions are averaged in clinical practice. Exact thresholds are device- and aetiology-specific.

Figure 12: Shear-wave elastography converts a measured wave speed into Young’s modulus through E = 3 rho cs^2, Eq. (youngs-modulus). Shaded bands show representative liver-stiffness ranges used to stage fibrosis. Because the relation is quadratic, measurement error in cs is doubled in E: the annotated example shows that a 10 percent speed error becomes a 21 percent stiffness error, which is why several acquisitions are averaged in clinical practice. Exact thresholds are device- and aetiology-specific.

knitr::kable(data.frame(
  cs = c(1.0, 1.5, 2.0, 2.5, 3.0),
  mu_kPa = round(rho_l*c(1,1.5,2,2.5,3)^2/1000, 2),
  E_kPa  = round(3*rho_l*c(1,1.5,2,2.5,3)^2/1000, 2),
  interpretation = c("normal","normal to mild","significant fibrosis",
                     "severe fibrosis","cirrhosis")),
  col.names = c("cs (m/s)","Shear modulus (kPa)","E (kPa)","Typical reading"),
  caption = "Shear-wave speed, shear modulus, and Young's modulus for liver. Note that E scales as the square of speed.")
Table 8: Shear-wave speed, shear modulus, and Young’s modulus for liver. Note that E scales as the square of speed.
cs (m/s) Shear modulus (kPa) E (kPa) Typical reading
1.0 1.00 3.00 normal
1.5 2.25 6.75 normal to mild
2.0 4.00 12.00 significant fibrosis
2.5 6.25 18.75 severe fibrosis
3.0 9.00 27.00 cirrhosis
cat(sprintf("Scale separation: bulk modulus K = rho c^2 = %.2f GPa, ",
            rho_soft*c_soft^2/1e9))
## Scale separation: bulk modulus K = rho c^2 = 2.49 GPa,
cat(sprintf("shear modulus mu = rho cs^2 = %.1f kPa at cs = 2 m/s.\n", rho_l*4/1000))
## shear modulus mu = rho cs^2 = 4.0 kPa at cs = 2 m/s.
cat(sprintf("Ratio K/mu = %.0f: tissue resists compression a million times more\n",
            (rho_soft*c_soft^2)/(rho_l*4)))
## Ratio K/mu = 622545: tissue resists compression a million times more
cat("strongly than it resists shear, and disease changes the shear property.\n")
## strongly than it resists shear, and disease changes the shear property.

Section 3.5 summary.

  • Harmonic imaging forms the image from tissue-generated \(2f_0\), rejecting clutter that never becomes nonlinear.
  • Microbubbles are intravascular tracers whose Minnaert resonance (Eq. (24)) happens to lie in the diagnostic band; they are destroyed above MI \(\approx 0.4\).
  • Shear-wave elastography pushes with radiation force (Eq. (25)), tracks \(c_s\), and reports \(E \approx 3\rho c_s^2\); the quadratic relation doubles relative error, and elasticity, isotropy, and homogeneity are all assumptions.

Checkpoint 3.5. Microbubbles and tissue both generate harmonics. Explain why CEUS still achieves high agent-to-tissue contrast, and state what happens to that contrast if the operator raises the output power to “see better.”

3.6 Acoustic Output, Bioeffects, and Safety

Ultrasound is non-ionizing, but non-ionizing is not the same as without bioeffects. Ultrasound deposits energy, and two distinct mechanisms can damage tissue. Each has its own displayed index, and every clinical scanner shows both.

3.6.1 Thermal bioeffects

Absorbed acoustic energy becomes heat at a volumetric rate \(q = 2\alpha I\), with \(\alpha\) the pressure absorption coefficient and \(I\) the time-averaged intensity. The regulated quantity is the spatial-peak temporal-average intensity \(I_{\mathrm{SPTA}}\), and the displayed Thermal Index is

\[\begin{equation} \mathrm{TI} = \frac{W_p}{W_{\mathrm{deg}}}, \tag{28} \end{equation}\]

the ratio of acoustic power used to the power estimated to raise tissue temperature by 1 °C. Variants exist for soft tissue (TIS), bone (TIB), and transcranial bone (TIC), because bone absorbs far more strongly than soft tissue and heats at its surface.

3.6.2 Mechanical bioeffects: cavitation

The rarefactional half-cycle can pull dissolved gas out of solution and drive bubble oscillation. Stable (non-inertial) cavitation produces sustained oscillation and microstreaming; inertial cavitation produces violent collapse, with local temperatures and pressures high enough to damage cells. Cavitation likelihood scales with peak rarefactional pressure and inversely with frequency, because a longer rarefactional half-cycle gives a bubble more time to grow. That scaling is encoded in the Mechanical Index

\[\begin{equation} \mathrm{MI} = \frac{p_{r.3}\,[\mathrm{MPa}]}{\sqrt{f_c\,[\mathrm{MHz}]}}, \tag{29} \end{equation}\]

where \(p_{r.3}\) is the peak rarefactional pressure derated for attenuation on the way in,

\[\begin{equation} p_{r.3} = p_r \cdot 10^{-0.0345\, f_c\, z / 20} \tag{30} \end{equation}\]

with \(f_c\) in MHz and \(z\) in cm.

pr <- seq(0, 4, by = 0.02)
fc <- c(2, 3.5, 5, 10)
mi <- do.call(rbind, lapply(fc, function(f)
  data.frame(pr = pr, MI = pr/sqrt(f),
             f = factor(sprintf("%.1f MHz", f), levels = sprintf("%.1f MHz", fc)))))

p_mi <- ggplot(mi, aes(pr, MI, colour = f)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = 1.9, linetype = "dashed", colour = bpad_pal[2]) +
  geom_hline(yintercept = 0.4, linetype = "dotted", colour = bpad_pal[5]) +
  annotate("text", x = 0.05, y = 2.05, hjust = 0, size = 3, colour = bpad_pal[2],
           label = "FDA limit MI = 1.9") +
  annotate("text", x = 0.05, y = 0.55, hjust = 0, size = 3, colour = bpad_pal[5],
           label = "microbubble destruction ~0.4") +
  coord_cartesian(ylim = c(0, 3)) +
  scale_colour_manual(values = bpad_pal[1:4]) +
  labs(x = expression("Derated peak rarefactional pressure "*p[r.3]*" (MPa)"),
       y = "Mechanical Index", title = "Mechanical Index and cavitation risk")

lim <- data.frame(
  application = factor(c("Ophthalmic","Fetal / neonatal","Cardiac",
                         "Peripheral vessel","General"),
                       levels = c("Ophthalmic","Fetal / neonatal","Cardiac",
                                  "Peripheral vessel","General")),
  ISPTA = c(50, 94, 430, 720, 720))
p_lim <- ggplot(lim, aes(application, ISPTA)) +
  geom_col(fill = bpad_pal[1]) +
  geom_text(aes(label = ISPTA), vjust = -0.4, size = 3) +
  scale_y_log10() + coord_cartesian(ylim = c(10, 1500)) +
  labs(x = NULL, y = expression(I[SPTA.3]*" limit (mW cm"^-2*", log scale)"),
       title = "FDA acoustic output limits") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))
bpad_grid(p_mi, p_lim, ncol = 2)
Acoustic safety. Left: the Mechanical Index, Eq. (mechanical-index), against derated peak rarefactional pressure for four transmit frequencies, with the FDA limit of 1.9 and the approximately 0.4 threshold above which microbubble contrast agents are destroyed. Because MI falls as one over the square root of frequency, a high-frequency probe may transmit higher pressure at the same index. Right: regulatory intensity limits by application, spanning more than an order of magnitude between ophthalmic and general-purpose use.

Figure 13: Acoustic safety. Left: the Mechanical Index, Eq. (mechanical-index), against derated peak rarefactional pressure for four transmit frequencies, with the FDA limit of 1.9 and the approximately 0.4 threshold above which microbubble contrast agents are destroyed. Because MI falls as one over the square root of frequency, a high-frequency probe may transmit higher pressure at the same index. Right: regulatory intensity limits by application, spanning more than an order of magnitude between ophthalmic and general-purpose use.

knitr::kable(data.frame(
  depth_cm = c(2, 5, 10, 15),
  derate_5MHz = round(10^(-0.0345*5*c(2,5,10,15)/20), 3),
  derate_2MHz = round(10^(-0.0345*2*c(2,5,10,15)/20), 3)),
  col.names = c("Depth (cm)","Derating factor at 5 MHz","at 2 MHz"),
  caption = "Attenuation derating, Eq. (derating). The in-situ pressure is always lower than the pressure measured in water, and the correction grows with both frequency and depth.")
Table 9: Attenuation derating, Eq. (derating). The in-situ pressure is always lower than the pressure measured in water, and the correction grows with both frequency and depth.
Depth (cm) Derating factor at 5 MHz at 2 MHz
2 0.961 0.984
5 0.905 0.961
10 0.820 0.924
15 0.742 0.888

Two indices, two mechanisms, two dependencies. TI tracks time-averaged intensity and therefore rises with duty cycle and dwell time, it is a heating index. MI tracks peak rarefactional pressure within a single cycle and is independent of duty cycle, it is a cavitation index. A low-duty-cycle, high-pressure pulse can have a high MI and a negligible TI; a long, gentle spectral-Doppler acquisition can do the reverse. Reading only one of them gives an incomplete picture of the exposure.

ALARA, and why the operator matters more than the machine. Diagnostic ultrasound has no confirmed adverse effect at diagnostic output levels, but the governing principle is ALARA, As Low As Reasonably Achievable. The operator, not the manufacturer, controls output, mode, and dwell time. Spectral Doppler and M-mode concentrate energy on a single line and deposit far more energy per unit volume than a sweeping B-mode. Obstetric spectral Doppler in the first trimester is specifically discouraged for this reason, and “keepsake” fetal imaging without medical indication is discouraged by every major professional body. The MI and TI are displayed precisely so that this judgment can be made at the bedside.

Section 3.6 summary.

  • Two mechanisms, two indices: TI for heating (time-averaged intensity, duty-cycle dependent) and MI for cavitation (peak rarefactional pressure, duty-cycle independent).
  • \(\mathrm{MI} = p_{r.3}/\sqrt{f_c}\), capped at 1.9; \(I_{\mathrm{SPTA.3}}\) capped at 720 mW cm\(^{-2}\) generally but 50 mW cm\(^{-2}\) for ophthalmic use.
  • Derating accounts for attenuation between transducer and target.
  • ALARA is an operator responsibility; Doppler and M-mode deposit more energy per volume than B-mode.

Checkpoint 3.6. A 2.5 MHz probe produces a derated peak rarefactional pressure of 2.2 MPa. Compute the MI, state whether it is within the FDA limit, and say whether this setting would be acceptable during a contrast-enhanced study.

3.7 The Photoacoustic Effect

3.7.1 From photon to pressure

Photoacoustic imaging turns light into sound inside the body. A nanosecond laser pulse is absorbed by a molecule (Chapter 2), depositing heat in a few billionths of a second. The abrupt temperature rise causes thermoelastic expansion: the heated region pushes outward, launching a broadband ultrasound wave that propagates to the surface, where a conventional transducer detects it. Light provides the contrast (which molecules absorbed the pulse) while sound provides the resolution and depth.

Under stress and thermal confinement, the initial pressure rise is

\[\begin{equation} p_0 = \Gamma\, \mu_a\, F, \tag{31} \end{equation}\]

where \(\mu_a\) is the optical absorption coefficient (the Chapter 2 quantity), \(F\) is the local optical fluence (J m\(^{-2}\)), and \(\Gamma\) is the dimensionless Grüneisen parameter

\[\begin{equation} \Gamma = \frac{\beta\, c^{2}}{C_p}, \tag{32} \end{equation}\]

with \(\beta\) the thermal-expansion coefficient, \(c\) the sound speed, and \(C_p\) the specific heat. \(\Gamma\) measures how efficiently absorbed heat becomes pressure. It is strongly temperature-dependent: for water \(\Gamma = 0.11\) at 20 °C but \(0.20\) at body temperature, and it is larger for lipid-rich tissue. Quoting \(\Gamma\) without a temperature is meaningless.

Connection to Chapter 2. Equation (31) is linear in \(\mu_a\), so the photoacoustic image is, up to the fluence \(F\) and the factor \(\Gamma\), a map of optical absorption, the very property that produced contrast in pulse oximetry, NIRS, and Raman in Chapter 2. Photoacoustics is best understood as Chapter 2’s absorption contrast read out acoustically, which is what lets it reach centimeters instead of millimeters.

3.7.2 When is \(p_0 = \Gamma\mu_a F\) valid?

Equation (31) assumes the laser pulse ends before the absorber can either expand or conduct heat away. For an absorber of characteristic size \(d_c\),

\[\begin{equation} \tau_{\text{st}} = \frac{d_c}{c}, \qquad \tau_{\text{th}} = \frac{d_c^{2}}{4\,\alpha_{\text{th}}}, \tag{33} \end{equation}\]

with \(\alpha_{\text{th}} \approx 1.4\times10^{-7}\) m\(^2\) s\(^{-1}\) the thermal diffusivity of tissue. The pulse duration must satisfy \(\tau_p \ll \tau_{\text{st}} \ll \tau_{\text{th}}\).

The absorber size also fixes the frequency of the emitted wave: a source of size \(d_c\) radiates a bipolar pressure transient of duration \(\sim d_c/c\), hence

\[\begin{equation} f_{\mathrm{PA}} \approx \frac{c}{d_c} . \tag{34} \end{equation}\]

d_c <- 10^seq(-6, -3, length.out = 300)
conf <- rbind(
  data.frame(d_um = d_c*1e6, t_ns = d_c/c_soft*1e9,             condition = "stress confinement"),
  data.frame(d_um = d_c*1e6, t_ns = d_c^2/(4*alpha_th)*1e9,     condition = "thermal confinement"))

p_conf <- ggplot(conf, aes(d_um, t_ns, colour = condition)) +
  annotate("rect", xmin = 1, xmax = 1e3, ymin = 5, ymax = 10,
           fill = "grey50", alpha = 0.30) +
  geom_line(linewidth = 0.95) +
  annotate("text", x = 1.3, y = 1.6, hjust = 0, size = 3, colour = "grey25",
           label = "5-10 ns laser pulse") +
  scale_x_log10() + scale_y_log10() +
  scale_colour_manual(values = bpad_pal[c(1,5)]) +
  labs(x = "Absorber size (um, log scale)", y = "Confinement time (ns, log scale)",
       title = "Both confinement conditions are easily met")

p_freq <- ggplot(data.frame(d_um = d_c*1e6, f_MHz = c_soft/d_c/1e6),
                 aes(d_um, f_MHz)) +
  annotate("rect", xmin = 1, xmax = 1e3, ymin = 1, ymax = 50,
           fill = bpad_pal[4], alpha = 0.15) +
  geom_line(linewidth = 0.95, colour = bpad_pal[4]) +
  scale_x_log10() + scale_y_log10() +
  labs(x = "Absorber size (um, log scale)",
       y = "Emitted frequency ~ c/d (MHz, log scale)",
       title = "Absorber size sets the acoustic frequency",
       subtitle = "shaded: typical clinical array bandwidth")
bpad_grid(p_conf, p_freq, ncol = 2)
Confinement conditions and emitted frequency, Eqs. (confinement) and (pa-frequency). Left: stress confinement scales linearly with absorber size and thermal confinement quadratically, so the two diverge; a typical 5 to 10 nanosecond laser pulse (shaded) satisfies both by orders of magnitude for every structure of biological interest. Right: the acoustic frequency an absorber emits, with the shaded band showing a typical clinical array bandwidth. Small absorbers radiate above the detector passband, which is the physical reason photoacoustic microscopy and tomography require different detectors.

Figure 14: Confinement conditions and emitted frequency, Eqs. (confinement) and (pa-frequency). Left: stress confinement scales linearly with absorber size and thermal confinement quadratically, so the two diverge; a typical 5 to 10 nanosecond laser pulse (shaded) satisfies both by orders of magnitude for every structure of biological interest. Right: the acoustic frequency an absorber emits, with the shaded band showing a typical clinical array bandwidth. Small absorbers radiate above the detector passband, which is the physical reason photoacoustic microscopy and tomography require different detectors.

sizes <- c(8, 10, 50, 500)
knitr::kable(data.frame(
  absorber = c("red blood cell","capillary","arteriole","small vessel"),
  size_um  = sizes,
  stress_ns  = round(sizes*1e-6/c_soft*1e9, 1),
  thermal_ms = signif((sizes*1e-6)^2/(4*alpha_th)*1e3, 3),
  f_PA_MHz   = round(c_soft/(sizes*1e-6)/1e6, 1)),
  col.names = c("Absorber","Size (um)","Stress conf. (ns)","Thermal conf. (ms)",
                "Emitted freq. (MHz)"),
  caption = "Confinement times and emitted frequency. The six-order-of-magnitude gap between the two confinement times is why a single nanosecond pulse satisfies both conditions everywhere.")
Table 10: Confinement times and emitted frequency. The six-order-of-magnitude gap between the two confinement times is why a single nanosecond pulse satisfies both conditions everywhere.
Absorber Size (um) Stress conf. (ns) Thermal conf. (ms) Emitted freq. (MHz)
red blood cell 8 5.2 0.114 192.5
capillary 10 6.5 0.179 154.0
arteriole 50 32.5 4.460 30.8
small vessel 500 324.7 446.000 3.1

Why the detector, not the laser, limits photoacoustic resolution. A 10 µm capillary emits near 150 MHz while a 500 µm vessel emits near 3 MHz. A clinical 5–15 MHz array is blind to the capillary’s signal not because the signal is absent but because it lies outside the passband. This is why photoacoustic microscopy uses high-frequency focused detectors and reaches micrometer resolution at millimeter depth, while photoacoustic tomography uses clinical-band arrays and reaches sub-millimeter resolution at centimeter depth. The two are the same physics with different receivers.

3.7.3 What actually limits the depth

The initial pressure is proportional to the local fluence, and fluence decays with the Chapter 2 diffuse penetration depth \(\delta = 1/\sqrt{3\mu_a(\mu_a + \mu_s')}\).

delta_of <- function(mu_a, mu_sp) 1/sqrt(3*mu_a*(mu_a + mu_sp))
cases <- data.frame(
  label = c("532 nm (visible)", "650 nm", "800 nm (optical window)", "1064 nm"),
  mu_a  = c(4.00, 0.30, 0.10, 0.30),      # cm^-1, representative soft tissue
  mu_sp = c(19.0, 14.3, 10.9, 7.2))       # cm^-1
cases$delta_mm <- 10*delta_of(cases$mu_a, cases$mu_sp)
cases$label <- factor(cases$label, levels = cases$label)

zmm <- seq(0, 60, by = 0.2)
pa <- do.call(rbind, lapply(seq_len(nrow(cases)), function(i)
  data.frame(z = zmm, sig = exp(-zmm/cases$delta_mm[i]), lab = cases$label[i])))

ggplot(pa, aes(z, sig, colour = lab)) +
  geom_hline(yintercept = 1e-3, linetype = "dashed", colour = "grey35") +
  geom_line(linewidth = 0.9) +
  annotate("text", x = 59, y = 1.6e-3, hjust = 1, size = 3, colour = "grey30",
           label = "representative detection floor") +
  scale_y_log10() + coord_cartesian(ylim = c(1e-5, 1)) +
  scale_colour_manual(values = bpad_pal[c(2,5,1,3)]) +
  labs(x = "Depth (mm)", y = "Relative fluence, and hence p0 (log scale)",
       title = "Optical attenuation sets the photoacoustic depth limit")
The photoacoustic depth budget. Initial pressure is proportional to local fluence, which decays exponentially with the Chapter 2 penetration depth; signal falls by an order of magnitude every 2.3 penetration depths. In the optical window near 800 nm, delta is about 5 to 6 mm, giving a usable range of roughly 3 to 5 cm before the signal reaches a representative detection floor. At 532 nm the range collapses to a few millimeters, which is why visible-wavelength photoacoustics is a microscopy technique and near-infrared photoacoustics is a tomography technique.

Figure 15: The photoacoustic depth budget. Initial pressure is proportional to local fluence, which decays exponentially with the Chapter 2 penetration depth; signal falls by an order of magnitude every 2.3 penetration depths. In the optical window near 800 nm, delta is about 5 to 6 mm, giving a usable range of roughly 3 to 5 cm before the signal reaches a representative detection floor. At 532 nm the range collapses to a few millimeters, which is why visible-wavelength photoacoustics is a microscopy technique and near-infrared photoacoustics is a tomography technique.

knitr::kable(data.frame(
  wavelength = as.character(cases$label), mu_a = cases$mu_a, mu_sp = cases$mu_sp,
  delta_mm = round(cases$delta_mm, 2),
  depth_1e3_mm = round(cases$delta_mm*log(1e3), 1)),
  col.names = c("Wavelength","mu_a (1/cm)","mu_s' (1/cm)","delta (mm)",
                "Depth at 10^-3 (mm)"),
  caption = "Diffuse penetration depth from Chapter 2 and the resulting photoacoustic depth budget.")
Table 11: Diffuse penetration depth from Chapter 2 and the resulting photoacoustic depth budget.
Wavelength mu_a (1/cm) mu_s’ (1/cm) delta (mm) Depth at 10^-3 (mm)
532 nm (visible) 4.0 19.0 0.60 4.2
650 nm 0.3 14.3 2.76 19.1
800 nm (optical window) 0.1 10.9 5.50 38.0
1064 nm 0.3 7.2 3.85 26.6

Closing the loop with Chapter 2. The penetration depth used here is exactly the quantity derived in Chapter 2, and the optical window computed there is the same window that governs photoacoustic depth. Photoacoustics does not escape the optical attenuation of Chapter 2, light still has to get in. What it escapes is optical scattering on the way out, because the return trip is acoustic and ultrasound scatters two to three orders of magnitude less than light in tissue. That asymmetry, one-way optical instead of two-way, is the entire depth advantage.

Photon-to-phonon story. A 760 nm photon is absorbed by a deoxyhaemoglobin molecule in a capillary. Within nanoseconds the local temperature rises by a few tens of millikelvin; the blood expands; a pressure transient of order kilopascals radiates outward as a wideband ultrasonic ping at roughly 150 MHz. The same probe that does B-mode listens for it, and, if it is a clinical array, hears only the lower-frequency contributions from larger vessels. Change the wavelength and a different chromophore answers: the basis of the multispectral imaging in Section 3.8.

3.8 Multispectral Photoacoustic Tomography

3.8.1 System architecture

A photoacoustic tomography system has three parts: light delivery (a tunable pulsed laser and fibre optics flooding the region at a chosen wavelength), acoustic detection (an array of transducers surrounding or facing the tissue), and image reconstruction (algorithms that back-propagate the detected pressure waves to their sources). Reconstruction is closely related to the back-projection tomography met for CT in Chapter 5: each detector records a projection of the source distribution, and the image is recovered by reversing the propagation.

3.8.2 Spectral unmixing and oxygen saturation

Multispectral photoacoustic tomography (MSOT) acquires images at several optical wavelengths and exploits the fact that each chromophore has a distinct absorption spectrum. After correcting for fluence, the absorption at wavelength \(\lambda\) is a linear combination of chromophore concentrations:

\[\begin{equation} \mu_a(\lambda) = \sum_i \varepsilon_i(\lambda)\, c_i , \tag{35} \end{equation}\]

where \(\varepsilon_i(\lambda)\) is the known molar absorption spectrum of chromophore \(i\) and \(c_i\) its concentration. Measuring at \(M\) wavelengths gives an over-determined linear system solved by least squares. With oxy- and deoxyhaemoglobin as the unknowns, the clinically central output is blood oxygen saturation

\[\begin{equation} sO_2 = \frac{[\mathrm{HbO_2}]}{[\mathrm{HbO_2}] + [\mathrm{Hb}]}, \tag{36} \end{equation}\]

obtained without any injected dye, purely from endogenous haemoglobin.

Connection to Chapters 1 and 2. Equation (35) is exactly the two-wavelength NIRS inversion of Chapter 2, now generalized: more wavelengths and more chromophores turn the \(2\times2\) system into an \(M\times n\) least-squares problem solved with the linear algebra of Chapter 1. Lab 3 implements it end to end, including the conditioning question of which wavelengths to choose.

3.8.3 The isosbestic point

An isosbestic point is a wavelength at which two forms of the same molecule absorb equally. For haemoglobin the principal near-infrared isosbestic point lies near 798–805 nm; several more exist in the visible (about 390, 422, 452, 500, 529, 545, 570, and 584 nm) and are used in high-resolution photoacoustic microscopy, where penetration is not the constraint.

Isosbestic wavelengths serve three purposes. Measuring at one gives total haemoglobin concentration, hence blood volume, independent of oxygenation. Measuring on either side of one gives maximum sensitivity to saturation, because the ordering of the two spectra reverses there. And an isosbestic measurement provides a saturation-independent reference against which fluence variation can be partially normalized.

## Molar extinction coefficients (cm^-1 M^-1), Prahl compilation -- identical
## source to Chapter 2, guaranteeing cross-chapter numerical consistency.
hb_tab <- data.frame(
  lam  = c(650, 700, 730, 760, 780, 800, 805, 820, 850, 860, 900, 940, 980, 1000),
  HbO2 = c( 320,  290,  390,  586,  710,  816,  800,  888, 1022, 1058, 1198, 1214, 1154, 1132),
  Hb   = c(3227, 1794, 1548, 1402, 1191,  762,  800,  726,  693,  691,  762,  693,  726,  726))

lam_g  <- seq(650, 1000, by = 2)
e_HbO2 <- exp(approx(hb_tab$lam, log(hb_tab$HbO2), lam_g)$y)
e_Hb   <- exp(approx(hb_tab$lam, log(hb_tab$Hb),   lam_g)$y)
iso    <- lam_g[which.min(abs(e_HbO2 - e_Hb))]

## Lipid and water: shapes only, scaled onto the same axis for comparison
lipid <- 120 + 700*exp(-((lam_g - 930)/22)^2)
water <-  40 + 620*exp(-((lam_g - 975)/30)^2)

df_chr <- rbind(
  data.frame(lambda = lam_g, absorb = e_HbO2, chromophore = "HbO2 (molar ext.)"),
  data.frame(lambda = lam_g, absorb = e_Hb,   chromophore = "Hb (molar ext.)"),
  data.frame(lambda = lam_g, absorb = lipid,  chromophore = "Lipid (scaled shape)"),
  data.frame(lambda = lam_g, absorb = water,  chromophore = "Water (scaled shape)"))
df_chr$chromophore <- factor(df_chr$chromophore,
  levels = c("HbO2 (molar ext.)","Hb (molar ext.)",
             "Lipid (scaled shape)","Water (scaled shape)"))

ggplot(df_chr, aes(lambda, absorb, colour = chromophore)) +
  geom_line(linewidth = 0.9) +
  geom_vline(xintercept = iso, linetype = "dotted", colour = "grey30") +
  annotate("text", x = iso + 8, y = 2600, hjust = 0, size = 3.1, colour = "grey25",
           label = sprintf("isosbestic\n%d nm", iso)) +
  scale_y_log10() +
  scale_colour_manual(values = c("HbO2 (molar ext.)" = bpad_pal[2],
                                 "Hb (molar ext.)"   = bpad_pal[1],
                                 "Lipid (scaled shape)" = bpad_pal[5],
                                 "Water (scaled shape)" = bpad_pal[3])) +
  labs(x = "Wavelength (nm)",
       y = expression("Absorption (cm"^-1*" M"^-1*", log scale)"),
       colour = NULL, title = "Endogenous photoacoustic contrast")
Near-infrared absorption of the dominant endogenous chromophores. Haemoglobin curves are molar extinction coefficients from the Prahl compilation, the same source used in Chapter 2, so the two chapters agree numerically. The HbO2 and Hb curves genuinely cross at the isosbestic point near 800 nm, located here from the data rather than asserted: below it deoxyhaemoglobin absorbs about six times more, above it oxyhaemoglobin absorbs more. That reversal is what makes label-free sO2 estimation possible, and it is why the two pulse-oximeter wavelengths of Chapter 2 straddle the same point.

Figure 16: Near-infrared absorption of the dominant endogenous chromophores. Haemoglobin curves are molar extinction coefficients from the Prahl compilation, the same source used in Chapter 2, so the two chapters agree numerically. The HbO2 and Hb curves genuinely cross at the isosbestic point near 800 nm, located here from the data rather than asserted: below it deoxyhaemoglobin absorbs about six times more, above it oxyhaemoglobin absorbs more. That reversal is what makes label-free sO2 estimation possible, and it is why the two pulse-oximeter wavelengths of Chapter 2 straddle the same point.

cat(sprintf("NIR isosbestic located from the data: %d nm (eps = %.0f cm^-1 M^-1)\n",
            iso, e_HbO2[lam_g == iso]))
## NIR isosbestic located from the data: 798 nm (eps = 805 cm^-1 M^-1)
cat(sprintf("Hb/HbO2 absorption ratio: %.1f at 700 nm and %.2f at 900 nm.\n",
            e_Hb[lam_g == 700]/e_HbO2[lam_g == 700],
            e_Hb[lam_g == 900]/e_HbO2[lam_g == 900]))
## Hb/HbO2 absorption ratio: 6.2 at 700 nm and 0.64 at 900 nm.
cat("The reversal across the isosbestic point is what a two-wavelength measurement\n")
## The reversal across the isosbestic point is what a two-wavelength measurement
cat("exploits -- exactly as pulse oximetry does in Chapter 2.\n")
## exploits -- exactly as pulse oximetry does in Chapter 2.

Endogenous versus exogenous contrast. Endogenous absorbers are already present: oxy- and deoxyhaemoglobin (the basis of \(sO_2\) and vascular imaging), melanin (skin, melanoma), lipids (plaque, breast), and water. They require no injection and underlie most clinical photoacoustics. Exogenous agents are introduced to add or target contrast: organic dyes such as indocyanine green (FDA-approved), gold nanoparticles and nanorods tuned to the near-infrared, and genetically encoded reporters. These extend photoacoustics toward molecular and targeted imaging at the cost of an injected agent.

Spectral colouring is the dominant confounder. Equation (31) contains the local fluence \(F(\lambda, z)\), which is itself wavelength-dependent and unknown, because the overlying tissue absorbs preferentially at some wavelengths. Deep in tissue the illuminating spectrum has therefore been reshaped, “coloured”, before it reaches the target. Unmixing raw \(p_0\) as if it were \(\mu_a\) biases \(sO_2\), and the bias grows with depth, typically toward apparently lower saturation. Correcting it requires either a model-based light-transport inversion, an invariant such as an isosbestic reference, or a calibration phantom. Lab 3 quantifies the size of the resulting error.

Sections 3.7–3.8 summary.

  • \(p_0 = \Gamma\mu_a F\) converts optical absorption into detectable ultrasound; \(\Gamma\) is temperature-dependent (0.11 at 20 °C, 0.20 at 37 °C for water).
  • Confinement requires \(\tau_p \ll d_c/c \ll d_c^2/4\alpha_{\text{th}}\), satisfied by a nanosecond pulse for every biological absorber.
  • Absorber size sets emitted frequency \(f \approx c/d_c\), so the detector bandwidth, not the laser, limits achievable resolution.
  • Depth is limited by optical fluence decay with the Chapter 2 penetration depth; the advantage over pure optics is that only the inbound trip is optical.
  • MSOT unmixes chromophore spectra by least squares to yield label-free \(sO_2\); the isosbestic point near 800 nm gives total haemoglobin; spectral colouring is the dominant confounder.

Checkpoint 3.7. Why does photoacoustics reach centimeters while confocal microscopy (Chapter 2) reaches only a few hundred micrometers, even though both rely on optical absorption? Identify precisely which leg of the round trip differs.

3.9 Clinical and Disease Applications

The four diagnostic axes of this chapter, anatomy and flow (Doppler), stiffness (elastography), perfusion (CEUS), and oxygenation (MSOT), map onto three disease domains.

3.9.1 Cardiovascular disease

Doppler echocardiography is the workhorse of cardiology. Spectral and colour Doppler quantify blood velocity to grade valvular stenosis and regurgitation, estimate pressure gradients through the simplified Bernoulli relation \(\Delta P \approx 4v^2\) (with \(v\) in m/s and \(\Delta P\) in mmHg), and map abnormal jets. M-mode resolves rapid valve and wall motion. In vascular medicine, Doppler measures carotid and peripheral-artery velocities to detect stenosis.

Clinical example, and the physics that forces it. In aortic stenosis the peak jet velocity feeds the Bernoulli estimate of the transvalvular gradient, a primary criterion for the timing of valve replacement. A 4 m/s jet implies \(\Delta P \approx 64\) mmHg. That measurement must be made with CW Doppler: Section 3.3.3 showed that at 12 cm depth even a 2 MHz probe has a Nyquist ceiling near 1.2 m/s, so no PW setting can reach 4 m/s. The price is that CW reports the highest velocity anywhere along the beam, so the operator must reason anatomically about where the jet lies.

3.9.2 Liver disease

Chronic liver injury (viral hepatitis, alcohol, metabolic disease) drives progressive fibrosis and ultimately cirrhosis, both of which stiffen the liver. Shear-wave elastography measures this stiffness non-invasively (Eqs. (26)(27)), staging fibrosis along the F0–F4 scale and reducing the need for percutaneous biopsy with its sampling error and complication risk. Serial measurements track progression or treatment response.

Clinical example. A liver stiffness rising past roughly 12–15 kPa is strongly suggestive of cirrhosis and prompts surveillance for portal hypertension and hepatocellular carcinoma. Because \(E \propto c_s^2\), a 10% error in the measured speed becomes a 21% error in the reported stiffness, which is why guidelines require multiple acquisitions in a standardized intercostal position during a breath-hold, and why thresholds do not transfer between devices.

3.9.3 Oncology

Tumours are betrayed by their vasculature and metabolism. CEUS reveals the chaotic neovascular perfusion of malignant lesions and characterizes focal masses by their enhancement kinetics. Multispectral photoacoustics maps tumour microvasculature and, by unmixing oxy- and deoxyhaemoglobin, quantifies hypoxia, a hallmark of aggressive, treatment-resistant tumours, without injected contrast. Falling \(sO_2\) or rising vascularity can signal progression, and their reversal can indicate response, making these candidate imaging biomarkers in the sense developed in Chapter 8.

3.10 Choosing a Method

From the clinical question to the measurement.

Question Method Physical quantity Binding limit
Where is the anatomy? B-mode impedance mismatch frequency versus depth
How fast is flow, and where? PW Doppler Doppler shift, range-gated \(v_{\max}d_{\max}\le c^2/8f_0\)
What is the peak velocity? CW Doppler Doppler shift, ungated no range information
Is flow present across a region? colour / power Doppler mean shift / integrated power angle dependence; wall filter
How fast does this structure move? M-mode echo position versus time one line only
Is the tissue stiff? shear-wave elastography \(c_s \Rightarrow E \approx 3\rho c_s^2\) assumes elastic, isotropic, homogeneous
Is it perfused, and how? CEUS microbubble nonlinear echo intravascular only; MI \(\lesssim 0.2\)
Is it hypoxic? multispectral photoacoustics \(\mu_a(\lambda) \Rightarrow sO_2\) fluence / spectral colouring
Is it cystic or solid? B-mode artifacts posterior enhancement, shadowing requires operator recognition
What is the total blood volume? photoacoustics at the isosbestic point \(\mu_a(800\ \text{nm})\) saturation-independent by construction

Method \(\rightarrow\) property \(\rightarrow\) disease. Doppler \(\rightarrow\) blood velocity \(\rightarrow\) valvular and vascular disease. Elastography \(\rightarrow\) stiffness \(\rightarrow\) liver fibrosis. CEUS \(\rightarrow\) perfusion \(\rightarrow\) lesion characterization. MSOT \(\rightarrow\) oxygenation \(\rightarrow\) tumour hypoxia. One chapter, four tissue properties, four clinical questions.

3.11 Computational Laboratories

The three laboratories build the core signal-processing skills of the chapter. Each is self-contained and runnable.

Lab 1: Pulse-echo timing, impedance, and the speed-of-sound artifact

Goal. Simulate a one-dimensional layered medium, compute the impedance mismatch at each boundary, synthesize the A-mode echo train, and quantify the depth error the scanner introduces by assuming a single speed of sound.

media <- data.frame(
  name  = c("Fat","Muscle","Liver","Bone"),
  thick = c(2, 3, 4, NA),                   # cm; last layer is a half-space
  c     = c(1450, 1580, 1550, 4080),        # m/s
  rho   = c(950, 1050, 1060, 1900))         # kg/m^3
media$Z <- media$rho*media$c/1e6            # MRayl

f_MHz <- 3.5
n_if  <- nrow(media) - 1
res <- data.frame(interface = character(n_if), depth_true = numeric(n_if),
                  R = numeric(n_if), t_us = numeric(n_if),
                  depth_app = numeric(n_if), amp = numeric(n_if),
                  stringsAsFactors = FALSE)

cum_depth <- 0     # true cumulative depth (cm)
cum_time  <- 0     # one-way travel time (s)
trans2way <- 1     # cumulative two-way intensity transmission through prior interfaces

for (i in seq_len(n_if)) {
  cum_depth <- cum_depth + media$thick[i]
  cum_time  <- cum_time  + (media$thick[i]/100)/media$c[i]
  Ri <- ((media$Z[i+1] - media$Z[i])/(media$Z[i+1] + media$Z[i]))^2
  att_dB <- a_dB*f_MHz*2*cum_depth          # round-trip attenuation
  echoI  <- trans2way*Ri*10^(-att_dB/10)
  res$interface[i]  <- paste(media$name[i], media$name[i+1], sep = "/")
  res$depth_true[i] <- cum_depth
  res$R[i]          <- Ri
  res$t_us[i]       <- 2*cum_time*1e6                    # round trip, microseconds
  res$depth_app[i]  <- c_soft*cum_time*100               # cm, assuming c_soft
  res$amp[i]        <- sqrt(echoI)                       # pressure-amplitude proxy
  trans2way <- trans2way*(1 - Ri)^2                      # forward and back
}
res$err_mm <- 10*(res$depth_app - res$depth_true)

knitr::kable(res, digits = c(0, 2, 4, 1, 3, 5, 2),
  col.names = c("Interface","True depth (cm)","R (intensity)","Echo time (us)",
                "Apparent depth (cm)","Echo amplitude","Depth error (mm)"),
  caption = "Pulse-echo interface analysis. The depth error is positive at every interface because the cumulative path speed is below the assumed 1540 m/s.")
Table 12: Pulse-echo interface analysis. The depth error is positive at every interface because the cumulative path speed is below the assumed 1540 m/s.
Interface True depth (cm) R (intensity) Echo time (us) Apparent depth (cm) Echo amplitude Depth error (mm)
Fat/Muscle 2 0.0086 27.6 2.124 0.04141 1.24
Muscle/Liver 5 0.0000 65.6 5.048 0.00064 0.48
Liver/Bone 9 0.4228 117.2 9.022 0.01715 0.22
ggplot(res) +
  geom_segment(aes(x = depth_app, xend = depth_app, y = 0, yend = amp),
               linewidth = 1.3, colour = bpad_pal[1]) +
  geom_point(aes(x = depth_app, y = amp), colour = bpad_pal[1], size = 2.4) +
  geom_point(aes(x = depth_true, y = 0), colour = "grey45", size = 2.2, shape = 17) +
  geom_text(aes(x = depth_app, y = amp, label = interface),
            vjust = -0.8, size = 3) +
  coord_cartesian(ylim = c(0, max(res$amp)*1.3)) +
  labs(x = "Depth (cm); stems = apparent, triangles = true",
       y = "Echo amplitude (a.u.)",
       title = "Lab 1: simulated A-mode echo train")
Simulated A-mode echo train through a fat-muscle-liver-bone stack at 3.5 MHz. Echo height reflects both the impedance mismatch and the cumulative two-way attenuation, so the bone interface dominates despite lying deepest. Grey markers show the true depth of each interface and coloured stems the apparent depth the scanner displays: every apparent depth is an overestimate, because the beam has crossed fat and liver, both slower than the assumed 1540 m/s.

Figure 17: Simulated A-mode echo train through a fat-muscle-liver-bone stack at 3.5 MHz. Echo height reflects both the impedance mismatch and the cumulative two-way attenuation, so the bone interface dominates despite lying deepest. Grey markers show the true depth of each interface and coloured stems the apparent depth the scanner displays: every apparent depth is an overestimate, because the beam has crossed fat and liver, both slower than the assumed 1540 m/s.

cat(sprintf("Largest error: %+.2f mm at the %s interface (%.1f%% of true depth).\n",
            res$err_mm[which.max(abs(res$err_mm))],
            res$interface[which.max(abs(res$err_mm))],
            100*max(abs(res$err_mm))/10/res$depth_true[which.max(abs(res$err_mm))]))
## Largest error: +1.24 mm at the Fat/Muscle interface (6.2% of true depth).

Discussion. The bone interface returns by far the strongest echo (large \(Z\) mismatch) and leaves little energy to image beyond it, the origin of acoustic shadowing. Note also that every apparent depth exceeds the true depth. The scanner assumes \(c_0 = 1540\) m/s, but the beam has crossed fat (1450 m/s) and liver (1550 m/s), so the true average speed along the path is below the assumed value and the round-trip time is longer than the scanner expects. The error is largest in relative terms for the shallowest interface, whose path is pure fat, and shrinks as faster muscle is added.

Sign convention. A layer slower than 1540 m/s pushes structures deeper in the display; a layer faster than 1540 m/s pulls them shallower. Fat is the common clinical case and it is the slow one, which is why measurements through a thick subcutaneous layer are systematically overstated.

Try it yourself. Replace the fat layer with a 4 cm layer of muscle (1580 m/s) and rerun. Does the error change sign? Then set every layer to exactly 1540 m/s and confirm the error vanishes, a useful check that the arithmetic, not the physics, is being tested.

Lab 2: Doppler estimation, aliasing, and the depth budget

Goal. Generate the slow-time signal that a pulsed-wave Doppler system samples at the PRF, estimate the Doppler frequency spectrally, show how velocities above the Nyquist limit alias, and connect the result to the depth–velocity budget of Eq. (23).

f0_l  <- 5e6
theta_l <- 0
PRF   <- 4000
v_nyq <- c_soft*PRF/(4*f0_l*cos(theta_l))

simulate_doppler <- function(v, N = 256, snr = 8, seed = 7) {
  set.seed(seed)
  fD <- 2*f0_l*v*cos(theta_l)/c_soft
  n  <- 0:(N - 1)
  x  <- exp(1i*2*pi*fD/PRF*n)                       # slow-time signal sampled at PRF
  x  <- x + (rnorm(N) + 1i*rnorm(N))/sqrt(2*snr)
  P  <- Mod(fft(x))^2
  freq   <- (0:(N - 1))*PRF/N
  freq_c <- ifelse(freq >= PRF/2, freq - PRF, freq) # centre on [-PRF/2, PRF/2)
  ord    <- order(freq_c)
  fD_meas <- freq_c[which.max(P)]
  list(df = data.frame(freq = freq_c[ord], power = P[ord]/max(P)),
       fD_true = fD, fD_meas = fD_meas,
       v_meas = c_soft*fD_meas/(2*f0_l*cos(theta_l)))
}
mk_plot <- function(r, v_true, ttl) {
  ggplot(r$df, aes(freq, power)) +
    geom_line(colour = bpad_pal[1]) +
    geom_vline(xintercept = r$fD_true, colour = bpad_pal[3], linetype = "dashed") +
    geom_vline(xintercept = c(-PRF/2, PRF/2), colour = "grey55", linetype = "dotted") +
    labs(x = "Doppler frequency (Hz)", y = "Normalized power", title = ttl,
         subtitle = sprintf("true %.2f, measured %.2f m/s (Nyquist %.2f)",
                            v_true, r$v_meas, v_nyq))
}
caseA <- simulate_doppler(0.20)
caseB <- simulate_doppler(0.50)

v_true_g <- seq(0, 1.2, by = 0.01)
v_meas_g <- vapply(v_true_g, function(v) simulate_doppler(v, snr = 40)$v_meas, numeric(1))
p_fold <- ggplot(data.frame(vt = v_true_g, vm = v_meas_g), aes(vt, vm)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey55", linetype = "dashed") +
  geom_hline(yintercept = 0, colour = "grey75") +
  geom_vline(xintercept = v_nyq*c(1, 3), linetype = "dotted", colour = bpad_pal[2]) +
  geom_line(colour = bpad_pal[2], linewidth = 0.9) +
  labs(x = "True velocity (m/s)", y = "Measured velocity (m/s)",
       title = "Aliasing folds the velocity axis")

bpad_grid(mk_plot(caseA, 0.20, "Below Nyquist: correct"),
          mk_plot(caseB, 0.50, "Above Nyquist: aliased"),
          p_fold, ncol = 3)
Pulsed-wave Doppler spectra. Left: a velocity below the Nyquist limit is recovered correctly, the measured peak coinciding with the true Doppler frequency. Centre: a velocity above the limit aliases, its spectral peak wrapping past the PRF/2 boundary so that fast forward flow is displayed as slower reverse flow. Right: measured velocity against true velocity across the range, showing the characteristic sawtooth of aliasing; the diagonal is truth and each fold occurs at an odd multiple of the Nyquist velocity.

Figure 18: Pulsed-wave Doppler spectra. Left: a velocity below the Nyquist limit is recovered correctly, the measured peak coinciding with the true Doppler frequency. Centre: a velocity above the limit aliases, its spectral peak wrapping past the PRF/2 boundary so that fast forward flow is displayed as slower reverse flow. Right: measured velocity against true velocity across the range, showing the characteristic sawtooth of aliasing; the diagonal is truth and each fold occurs at an odd multiple of the Nyquist velocity.

cat(sprintf("PRF = %.0f Hz, f0 = %.1f MHz -> v_nyq = %.3f m/s\n",
            PRF, f0_l/1e6, v_nyq))
## PRF = 4000 Hz, f0 = 5.0 MHz -> v_nyq = 0.308 m/s
cat(sprintf("This PRF is only legal to a depth of c/(2*PRF) = %.1f cm.\n",
            c_soft/(2*PRF)*100))
## This PRF is only legal to a depth of c/(2*PRF) = 19.2 cm.
cat(sprintf("Depth-velocity budget at 5 MHz: %.4f m^2/s, so at %.1f cm the ceiling\n",
            c_soft^2/(8*f0_l), c_soft/(2*PRF)*100))
## Depth-velocity budget at 5 MHz: 0.0593 m^2/s, so at 19.2 cm the ceiling
cat(sprintf("is %.2f m/s -- exactly the Nyquist velocity above, as Eq. (depth-velocity) requires.\n",
            c_soft^2/(8*f0_l*(c_soft/(2*PRF)))))
## is 0.31 m/s -- exactly the Nyquist velocity above, as Eq. (depth-velocity) requires.

Discussion. With PRF = 4 kHz and \(f_0\) = 5 MHz the Nyquist velocity is 0.31 m/s. The 0.20 m/s flow is recovered correctly; the 0.50 m/s flow exceeds \(v_{\max}\) and its spectral peak wraps past \(\pm\mathrm{PRF}/2\), so it is measured as a slower flow in the reverse direction. The right-hand panel shows the resulting sawtooth: measured velocity is a folded version of true velocity, and the folding is entirely deterministic.

The closing calculation makes the point that Section 3.3.3 argued in general: this PRF is only legal down to 19.3 cm, and the depth–velocity budget evaluated at that depth returns exactly the Nyquist velocity computed independently. The two constraints are one constraint.

Try it yourself. Raise the PRF to 8 kHz and confirm that \(v_{\max}\) doubles, then compute the maximum depth that PRF permits and check whether your clinical target still lies inside it. Next set \(\theta = 60^\circ\) and observe that \(v_{\max}\) doubles again, then use Eq. (20) to compute what a \(5^\circ\) angle error now costs.

Lab 3: Photoacoustic spectral unmixing, conditioning, and spectral colouring

Goal. Recover oxy- and deoxyhaemoglobin concentrations from multi-wavelength photoacoustic measurements by least squares, compute \(sO_2\), and then quantify the two things that degrade it in practice: poor wavelength choice (conditioning) and spectral colouring (uncorrected wavelength-dependent fluence).

## Molar extinction coefficients at the chosen wavelengths, interpolated from the
## same Prahl table used in Section 3.8 -- so Lab 3 and the figure above agree.
wl_good <- c(700, 730, 760, 800, 850, 900)     # straddles the isosbestic point
wl_bad  <- c(850, 870, 890, 910, 930, 950)     # all on one side of it

eps_at <- function(wl) cbind(
  HbO2 = exp(approx(hb_tab$lam, log(hb_tab$HbO2), wl, rule = 2)$y),
  Hb   = exp(approx(hb_tab$lam, log(hb_tab$Hb),   wl, rule = 2)$y))

E_good <- eps_at(wl_good)
E_bad  <- eps_at(wl_bad)

unmix <- function(y, E) {
  beta <- solve(t(E) %*% E, t(E) %*% y)
  beta[beta < 0] <- 0
  as.numeric(beta[1]/max(beta[1] + beta[2], .Machine$double.eps))
}

run_experiment <- function(E, wl, colouring = 0, noise = 0.05, seed = 3) {
  set.seed(seed)
  true_sO2 <- rep(seq(0, 1, by = 0.05), each = 12)
  vapply(true_sO2, function(s) {
    c_true <- c(s, 1 - s)                       # [HbO2, Hb], total Hb = 1
    mua <- as.numeric(E %*% c_true)
    ## spectral colouring: fluence falls with wavelength-dependent attenuation
    fluence <- exp(-colouring*(wl - min(wl))/100)
    y <- mua*fluence*(1 + rnorm(length(mua), 0, noise))
    unmix(y, E)
  }, numeric(1)) -> est
  data.frame(true_sO2 = true_sO2, est_sO2 = est)
}

good  <- run_experiment(E_good, wl_good)
bad   <- run_experiment(E_bad,  wl_bad)
colrd <- run_experiment(E_good, wl_good, colouring = 0.45)

mk_sc <- function(df, ttl, sub, col) {
  ggplot(df, aes(true_sO2, est_sO2)) +
    geom_abline(slope = 1, intercept = 0, colour = "grey55", linetype = "dashed") +
    geom_point(alpha = 0.45, colour = col, size = 1.4) +
    coord_equal(xlim = c(0, 1), ylim = c(0, 1)) +
    labs(x = expression("True "*sO[2]), y = expression("Recovered "*sO[2]),
         title = ttl, subtitle = sub)
}
rmse <- function(d) sqrt(mean((d$est_sO2 - d$true_sO2)^2))
bias <- function(d) mean(d$est_sO2 - d$true_sO2)

bpad_grid(
  mk_sc(good,  "Well-conditioned",
        sprintf("kappa = %.1f, RMSE = %.3f", kappa(E_good), rmse(good)), bpad_pal[3]),
  mk_sc(bad,   "Ill-conditioned",
        sprintf("kappa = %.1f, RMSE = %.3f", kappa(E_bad), rmse(bad)), bpad_pal[2]),
  mk_sc(colrd, "Spectral colouring",
        sprintf("bias = %+.3f, RMSE = %.3f", bias(colrd), rmse(colrd)), bpad_pal[1]),
  ncol = 3)
Photoacoustic spectral unmixing. Left: with a well-chosen wavelength set straddling the isosbestic point, the least-squares estimate tracks the true saturation across the whole physiological range despite 5 percent measurement noise. Centre: restricting the wavelengths to one side of the isosbestic point makes the design matrix ill-conditioned and the estimates scatter badly, even though the noise level is unchanged. Right: uncorrected spectral colouring, a wavelength-dependent fluence that the model ignores, produces a systematic bias rather than scatter, and the bias grows with depth. Scatter can be averaged away; bias cannot.

Figure 19: Photoacoustic spectral unmixing. Left: with a well-chosen wavelength set straddling the isosbestic point, the least-squares estimate tracks the true saturation across the whole physiological range despite 5 percent measurement noise. Centre: restricting the wavelengths to one side of the isosbestic point makes the design matrix ill-conditioned and the estimates scatter badly, even though the noise level is unchanged. Right: uncorrected spectral colouring, a wavelength-dependent fluence that the model ignores, produces a systematic bias rather than scatter, and the bias grows with depth. Scatter can be averaged away; bias cannot.

knitr::kable(data.frame(
  scenario = c("straddling isosbestic (700-900 nm)",
               "one side only (850-950 nm)",
               "straddling, with spectral colouring"),
  condition_number = round(c(kappa(E_good), kappa(E_bad), kappa(E_good)), 1),
  RMSE = round(c(rmse(good), rmse(bad), rmse(colrd)), 4),
  bias = round(c(bias(good), bias(bad), bias(colrd)), 4)),
  col.names = c("Wavelength set","Condition number","RMSE","Bias"),
  caption = "Two distinct failure modes in photoacoustic unmixing. Ill-conditioning inflates variance; spectral colouring introduces bias. They require different remedies.")
Table 13: Two distinct failure modes in photoacoustic unmixing. Ill-conditioning inflates variance; spectral colouring introduces bias. They require different remedies.
Wavelength set Condition number RMSE Bias
straddling isosbestic (700-900 nm) 3.6 0.0273 0.0004
one side only (850-950 nm) 48.8 0.3158 -0.0042
straddling, with spectral colouring 3.6 0.3210 -0.2926

Discussion. Each simulated voxel mixes oxy- and deoxyhaemoglobin in proportions set by the true \(sO_2\); the photoacoustic amplitude at each wavelength is their weighted sum (Eq. (35)). Three lessons follow.

Conditioning is a wavelength-selection problem. The well-chosen set straddles the isosbestic point, where the two spectra reverse their ordering, so the columns of \(E\) are well separated. The poor set lies entirely above the isosbestic point, where both spectra are flat and nearly proportional; the columns become close to collinear, the condition number rises, and the same noise produces far larger errors. This is precisely the conditioning argument made for two-wavelength NIRS in Chapter 2, wavelength choice is a numerical-stability decision, not a convenience.

Colouring produces bias, not scatter. Uncorrected wavelength-dependent fluence leaves the estimator precise but wrong, and averaging more measurements does not help. In real MSOT the bias grows with depth and typically drives apparent \(sO_2\) downward, which is dangerous precisely because tumour hypoxia is the finding of interest.

The two failures need different fixes. Ill-conditioning is fixed by choosing wavelengths that straddle an isosbestic point or by adding wavelengths; colouring is fixed only by modelling or measuring the light transport, or by anchoring to an invariant.

Try it yourself. Add a third chromophore (melanin, roughly \(\propto\lambda^{-3.5}\)) to the design matrix and rerun with the same six wavelengths. Does the condition number rise? Then add three more wavelengths and see whether the extra measurements recover the lost accuracy, the practical question every MSOT protocol has to answer.

3.12 Conclusions and Discussion

Ultrasound and photoacoustic imaging show how far a single physical idea, a propagating pressure wave, can be pushed clinically. From the wave equation and a handful of boundary conditions follow the whole grayscale anatomical image (impedance mismatches and echo timing), the whole flow image (Doppler shifts off moving blood), and, with one extra ingredient, the whole stiffness image (shear waves launched by radiation force). Photoacoustics then reuses the very same detectors to capture an entirely different kind of contrast: by letting a laser pulse rather than a piezoelectric crystal be the sound source, it carries the optical-absorption contrast of Chapter 2 to depths that light alone can never reach.

Three themes recur.

Trade-offs, not rules. Frequency buys resolution at the cost of penetration; pulse length trades axial resolution against sensitivity; depth and line count trade against frame rate; pulsed Doppler trades velocity range against range localization; and photoacoustics trades optical resolution for a hundredfold gain in depth. The most important of these is quantitative and often omitted: \(v_{\max}d_{\max}\le c^2/8f_0\) is not a guideline but a bound, and it is why continuous-wave Doppler exists. Recognizing these as instances of shared physics, rather than a list of disconnected clinical rules, is what lets a practitioner choose settings rationally.

Assumptions are visible in the artifacts. The scanner assumes a single sound speed, straight-line propagation, main-beam-only echoes, and single reflection. Every artifact in Section 3.4 is one of those assumptions failing, and two of them, posterior enhancement and shadowing, are so informative that they function as measurements. An artifact you understand is data; an artifact you do not understand is a misdiagnosis.

Mathematical economy. Almost every quantitative result rests on tools already developed: the wave equation and Fourier analysis (Chapter 1) for propagation, beam patterns, and spectral Doppler; the Nyquist criterion (Chapter 1) for aliasing; exponential attenuation (Chapters 1–2) for penetration; least squares and conditioning (Chapter 1) for spectral unmixing; Rayleigh statistics (Chapter 1) for speckle; and optical absorption and penetration depth (Chapter 2) for photoacoustic contrast and its depth limit. The MSOT oxygen-saturation calculation is, quite literally, the two-wavelength NIRS inversion of Chapter 2 scaled up, including its conditioning requirement that the wavelengths straddle the isosbestic point.

Several frontiers are advancing quickly. Ultrafast plane-wave imaging acquires thousands of frames per second, enabling shear-wave elastography, functional ultrasound of brain activity, and ultrasound localization microscopy, which super-resolves microvasculature below the diffraction limit by tracking individual microbubbles. Handheld and point-of-care systems, increasingly with on-device machine learning, are moving diagnostic ultrasound to the bedside and to low-resource settings. And clinical photoacoustic systems are maturing toward label-free oxygenation imaging of breast lesions, skin cancer, and tumour hypoxia. In every case the images produced here become the inputs to the processing pipelines of Chapter 7 and the modelling and evaluation pipelines of Chapter 8.

Key points

  • Sound in tissue is a longitudinal pressure wave; \(c=\sqrt{K/\rho}\approx1540\) m/s, \(c=f\lambda\), and \(Z=\rho c\) links pressure to particle motion via \(p=Zu\).
  • Echoes come from impedance mismatches; \(R\) is signed and \(R_I=R^2\).
  • Attenuation is exponential and frequency-dependent; nepers and decibels differ by 4.343 (intensity) or 8.686 (pressure), and \(\mu_I=2\alpha\).
  • Speckle is coherent interference with a fixed envelope SNR of 1.913, not additive noise.
  • Three resolutions, axial (pulse length), lateral (aperture), elevational (fixed lens) — and the third causes slice-thickness artifact.
  • Beamforming is controlled delay; the beam pattern is the Fourier transform of the aperture; grating lobes need \(d\le\lambda/(1+|\sin\theta_{\max}|)\).
  • TGC undoes attenuation until the noise floor; log compression fits 100 dB onto a display; frame rate is bounded by \(c/2dN\).
  • \(d=ct/2\) for ranging; \(f_D=2f_0v\cos\theta/c\) for flow, with \(\Delta v/v=\tan\theta\,\Delta\theta\).
  • Pulsed Doppler obeys \(v_{\max}d_{\max}\le c^2/8f_0\cos\theta\); CW Doppler is required for deep high-velocity jets.
  • Artifacts are assumption failures; posterior enhancement and shadowing are diagnostic.
  • Harmonic imaging rejects clutter; microbubbles resonate in the diagnostic band by the Minnaert relation; SWE reports \(E\approx3\rho c_s^2\) under elasticity, isotropy, and homogeneity assumptions.
  • MI and TI index cavitation and heating respectively; ALARA is an operator responsibility.
  • \(p_0=\Gamma\mu_aF\) under stress and thermal confinement; absorber size sets emitted frequency \(f\approx c/d_c\), so detector bandwidth limits resolution.
  • MSOT unmixes haemoglobin spectra for label-free \(sO_2\); conditioning demands wavelengths straddling the isosbestic point, and spectral colouring biases the result.

Connections to other BPAD chapters

  • Chapter 1 (Mathematical and Statistical Foundations). The wave equation; Fourier methods for beam patterns and spectral Doppler; the Nyquist criterion for aliasing; exponential decay for attenuation; least squares, conditioning, and the condition number for unmixing; Rayleigh statistics for speckle; and the delta method for Doppler angle-error and \(sO_2\) uncertainty propagation.
  • Chapter 2 (Optical and Thermal Methods). Optical absorption of haemoglobin, melanin, and lipid supplies photoacoustic contrast; the diffuse penetration depth \(\delta=1/\sqrt{3\mu_a(\mu_a+\mu_s')}\) sets the photoacoustic depth budget; the isosbestic point and the two-wavelength NIRS inversion generalize directly to MSOT; and the quarter-wave matching condition is the acoustic analogue of an anti-reflection optical coating.
  • Chapter 4 (MRI). Both modalities encode space through a physical gradient and both are limited by sampling; MRI’s \(k\)-space and ultrasound’s slow-time Doppler sampling are two applications of the same Nyquist reasoning.
  • Chapter 5 (X-ray Imaging and CT). Photoacoustic tomographic reconstruction shares the back-projection framework formalized there; both are inverse problems built on line or surface integrals. The contrast is in dose: ultrasound and photoacoustics deposit energy but are non-ionizing.
  • Chapter 6 (Nuclear Medicine). Both photoacoustics and emission imaging detect signals generated inside the patient rather than transmitted through, and both are therefore limited by how much signal escapes rather than by how much is delivered.
  • Chapter 7 (General Medical Image Processing). B-mode, Doppler, elastography, and MSOT images feed the denoising, registration, segmentation, and quantification pipelines developed there, including speckle-reduction methods that must respect the finding in Section 3.1.5 that speckle carries tissue information.
  • Chapter 8 (Data Modeling, AI, and Machine Learning). Elastography stiffness, CEUS perfusion kinetics, and MSOT \(sO_2\) are candidate imaging biomarkers, and the evaluation, calibration, and clinical-utility machinery of that chapter is what determines whether they are usable.

Glossary

Term Meaning
A-mode Echo amplitude versus depth along a single line
Acoustic impedance \(Z=\rho c\); relates pressure to particle velocity
ALARA As Low As Reasonably Achievable; the governing exposure principle
Aliasing Spectral wraparound when a Doppler shift exceeds PRF/2
Apodization Tapering element weights to trade main-lobe width for side-lobe level
Attenuation coefficient Exponential loss rate; in Np m\(^{-1}\) or dB cm\(^{-1}\) MHz\(^{-1}\)
B-mode Two-dimensional grayscale brightness image
Beamforming Focusing and steering by controlled element delays
Cavitation Bubble activity driven by rarefactional pressure; stable or inertial
CEUS Contrast-enhanced ultrasound using microbubble agents
Derating Correcting measured pressure for attenuation to estimate in-situ exposure
Doppler shift Frequency change from a moving scatterer, \(f_D=2f_0v\cos\theta/c\)
Elastography Imaging of tissue stiffness
Elevational resolution Resolution perpendicular to the image plane; set by a fixed lens
Grating lobe Spurious beam from element pitch exceeding the limit in Eq. (15)
Grüneisen parameter \(\Gamma=\beta c^2/C_p\); efficiency of heat-to-pressure conversion
Harmonic imaging Imaging from the tissue-generated \(2f_0\) component
Isosbestic point Wavelength at which two forms of a molecule absorb equally (~800 nm for Hb)
M-mode Depth versus time along one line; high temporal resolution
Mechanical Index \(\mathrm{MI}=p_{r.3}/\sqrt{f_c}\); cavitation-risk index, FDA limit 1.9
Minnaert resonance Bubble resonance frequency, Eq. (24)
MSOT Multispectral optoacoustic (photoacoustic) tomography
Neper Natural-logarithm attenuation unit; 1 Np = 8.686 dB in pressure
Particle velocity Oscillation velocity of the medium, not the wave speed
Photoacoustic effect Pulsed light absorbed, converted to heat, launching an acoustic wave
PRF Pulse repetition frequency
Posterior enhancement Increased brightness deep to a low-attenuation structure
Rayleigh scattering Scattering off objects much smaller than a wavelength; \(\propto f^4\)
Reverberation Repeated inter-reflection displayed as equally spaced deeper copies
Shear wave Transverse wave; speed \(c_s=\sqrt{\mu/\rho}\), 1–10 m/s in soft tissue
Slice-thickness artifact Out-of-plane echoes displayed as in-plane; from elevational width
Specular reflection Reflection from a smooth boundary larger than a wavelength
Speckle Coherent interference texture with fixed envelope SNR of 1.913
Spectral colouring Wavelength-dependent fluence attenuation that biases unmixing
Stress confinement Pulse shorter than the acoustic transit time of the absorber
TGC Time-gain compensation; depth-dependent amplification
Thermal Index Ratio of applied power to that raising tissue by 1 °C
Thermal confinement Pulse shorter than the thermal diffusion time of the absorber

Review Questions

  1. Starting from continuity, Euler’s equation, and the equation of state, outline how the acoustic wave equation and \(c=\sqrt{K/\rho}\) arise together.
  2. Why does soft tissue have a nearly constant speed of sound, and what artifacts arise when that assumption fails? In which direction does fat shift a displayed depth?
  3. Explain physically why a large acoustic-impedance mismatch produces both a strong echo and a strong shadow.
  4. Distinguish the amplitude reflection coefficient from the intensity reflection coefficient. Why can only one of them be negative, and what does the sign mean?
  5. State the relations \(p=Zu\) and \(I=p_A^2/2\rho c\), and use them to explain why tissue moves only nanometers under a megapascal pulse.
  6. Convert 0.5 dB cm\(^{-1}\) MHz\(^{-1}\) at 4 MHz into a pressure attenuation coefficient in Np cm\(^{-1}\) and an intensity coefficient. Why do the two differ by a factor of two?
  7. Distinguish specular reflection, Rayleigh scattering, and speckle. Which enables Doppler imaging, and why?
  8. Why is the envelope SNR of fully developed speckle a constant, and what does that constancy prove about the nature of speckle?
  9. Separate axial, lateral, and elevational resolution: what determines each, and which can the operator not adjust?
  10. What is the purpose of the matching layer and the backing layer, and what does each cost?
  11. Why does beam steering with a phased array risk grating lobes, and what pitch is required to avoid them at \(\pm45^\circ\)?
  12. Explain time-gain compensation. Why does the bottom of a deep image often appear uniformly bright and featureless?
  13. Derive the frame-rate bound and explain how plane-wave imaging escapes it.
  14. State the Doppler equation and derive the angle-error propagation relation. Why is \(60^\circ\) the clinical limit?
  15. Derive the pulsed-Doppler depth–velocity bound and explain why it forces the use of CW Doppler for aortic stenosis.
  16. For each of shadowing, posterior enhancement, reverberation, mirror image, and slice-thickness artifact, state which scanner assumption fails.
  17. Why do microbubbles resonate in the diagnostic band, and why must CEUS be performed at low MI?
  18. Explain how shear-wave elastography converts a measured wave speed into a stiffness, and list the three mechanical assumptions involved.
  19. Contrast the Mechanical and Thermal Indices: what mechanism does each track, and how does each depend on duty cycle?
  20. In \(p_0=\Gamma\mu_aF\), what does each factor represent, and why is the image essentially a map of optical absorption?
  21. State the stress- and thermal-confinement conditions and explain why a nanosecond pulse satisfies both for every biological absorber.
  22. Why does absorber size determine the emitted acoustic frequency, and why does that make detector bandwidth the limit on photoacoustic resolution?
  23. How does MSOT compute blood oxygen saturation without an injected dye? Name the two distinct ways the estimate can fail and the remedy for each.
  24. What does a measurement at the isosbestic point report, and why is it independent of oxygenation?

Problems

Problem 1. Wavelength and axial resolution

A 7.5 MHz probe images soft tissue (\(c=1540\) m/s). (a) Find the wavelength. (b) For a 3-cycle pulse, find the axial resolution.

Problem 2. Reflection at a soft-tissue interface

At a fat/muscle boundary, \(Z_{\text{fat}}=1.38\) MRayl and \(Z_{\text{muscle}}=1.66\) MRayl. Compute the amplitude reflection coefficient, the intensity reflection coefficient, and the intensity transmission coefficient. State the sign of \(R\) and its physical meaning.

Problem 3. Why gel?

For soft tissue (\(Z=1.63\) MRayl) to air (\(Z=0.0004\) MRayl), compute \(R\) and \(R_I\), and explain in one sentence why coupling gel is required.

Problem 4. Attenuation and penetration

Soft-tissue attenuation is 0.5 dB cm\(^{-1}\) MHz\(^{-1}\). For a 5 MHz beam imaging a structure at 8 cm depth, find the total round-trip attenuation in dB and the fraction of intensity remaining.

Problem 5. Pulse-echo depth

An echo returns 91 µs after the pulse is transmitted. At what depth is the reflector?

Problem 6. Doppler shift and aliasing

A 4 MHz probe insonates blood flowing at 0.6 m/s at \(\theta=60^\circ\), \(c=1540\) m/s. (a) Compute the Doppler shift. (b) What minimum PRF avoids aliasing? (c) What maximum imaging depth does that PRF permit?

Problem 7. Nyquist velocity

A pulsed-wave system uses \(f_0=5\) MHz and PRF = 5 kHz at \(\theta=0\). What is the maximum unambiguous velocity?

Problem 8. Quarter-wave matching layer

A PZT element (\(Z=30\) MRayl) couples to tissue (\(Z=1.6\) MRayl). (a) What matching-layer impedance minimizes reflection? (b) With \(c_m=2000\) m/s and a 5 MHz centre frequency, what thickness is required? (c) What fraction of intensity would reflect without the layer?

Problem 9. Elastography

Shear-wave elastography returns \(c_s=2.5\) m/s (\(\rho=1000\) kg m\(^{-3}\)). Compute the shear modulus and Young’s modulus in kPa and interpret. If \(c_s\) carries a 10% uncertainty, what is the uncertainty in \(E\)?

Problem 10. Photoacoustic initial pressure

A region has \(\Gamma=0.2\), \(\mu_a=4\) cm\(^{-1}\), and receives \(F=20\) mJ cm\(^{-2}\). Compute \(p_0\) in pascals, checking units explicitly.

Problem 11. Why multiple wavelengths?

  1. Why must MSOT use at least two wavelengths to separate oxy- and deoxyhaemoglobin?
  2. What physiological quantity does a measurement at the isosbestic point report, and why?

Problem 12. Harmonic imaging

Explain why imaging at \(2f_0\) improves contrast and reduces artifact relative to imaging at the fundamental.

Problem 13. Nepers and decibels

Soft tissue at 3 MHz has \(a=0.5\) dB cm\(^{-1}\) MHz\(^{-1}\). (a) Convert to the pressure attenuation coefficient \(\alpha\) in Np cm\(^{-1}\) and the intensity coefficient \(\mu_I\). (b) At what depth has the intensity fallen to \(1/e\)? (c) At what depth has the pressure fallen to \(1/e\)? (d) Explain why the two answers differ by a factor of two.

Problem 14. The depth–velocity limit

A 3 MHz cardiac probe images to 14 cm. (a) What is the maximum PRF? (b) The maximum unambiguous velocity at \(\theta=0\)? (c) A mitral regurgitation jet reaches 5 m/s. Show that no PRF permits PW measurement, and state what must be used instead.

Problem 15. Doppler angle-error propagation

Show that \(\Delta v/v=\tan\theta\,\Delta\theta\). Evaluate for \(\Delta\theta=5^\circ\) at \(\theta=20^\circ\) and at \(\theta=70^\circ\), and explain the clinical guidance to keep \(\theta\le60^\circ\).

Problem 16. Mechanical index

A 2.5 MHz probe produces a derated peak rarefactional pressure of 2.2 MPa. (a) Compute the MI. (b) Is it within the FDA limit? (c) Would this setting be acceptable during a contrast-enhanced study, and why?

Problem 17. Photoacoustic safety and signal

A 700 nm pulse delivers the ANSI maximum permissible skin exposure of 20 mJ cm\(^{-2}\) to blood with \(\mu_a=4\) cm\(^{-1}\), \(\rho=1060\) kg m\(^{-3}\), \(C_p=3600\) J kg\(^{-1}\)K\(^{-1}\). (a) Compute the per-pulse temperature rise. (b) Compute \(p_0\) for \(\Gamma=0.2\). (c) Explain why a technique that generates kilopascals of pressure raises tissue temperature by only millikelvin.

Problem 18. Frame rate

A cardiac sector uses 90 lines to 16 cm depth with two transmit focal zones. (a) Compute the maximum frame rate. (b) Does it support real-time wall-motion assessment? (c) What single change would most improve it, and what is the cost?

Solutions

Table 14: Numerical answers, computed from the chapter constants so that solutions and text cannot disagree.
Problem Computed answer
1(a) lambda = 0.205 mm
1(b) 0.308 mm
2 R = +0.0921, R_I = 0.00848, T_I = 0.99152
3 R = -0.9995, R_I = 0.9990
4 40 dB, fraction 1.0e-04
5 7.01 cm
6(a) 1558 Hz
6(b) > 3117 Hz
6(c) 24.7 cm
7 0.385 m/s
8(a) 6.93 MRayl
8(b) 100.0 um
8(c) 80.8%
9 mu = 6.25 kPa, E = 18.75 kPa, dE/E = 20%
10 16.0 kPa
13(a) alpha = 0.1727 Np/cm, mu_I = 0.3454 Np/cm
13(b) 2.90 cm
13(c) 5.79 cm
14(a) 5500 Hz
14(b) 0.706 m/s
14(c) needs 0.70 m^2/s vs budget 0.0988 -> CW required
15 3.2% at 20 deg, 24.0% at 70 deg
16(a) MI = 1.39
17(a) dT = 0.0210 K
17(b) p0 = 16.0 kPa
18(a) 26.7 Hz

Solution 1

  1. \(\lambda = c/f = 1540/(7.5\times10^6) = 2.05\times10^{-4}\) m \(= 0.205\) mm.
  2. \(\Delta_{\text{axial}} = N\lambda/2 = 3(0.205)/2 = 0.308\) mm.

Solution 2

\(R = (1.66-1.38)/(1.66+1.38) = +0.0921\). Because \(Z_2 > Z_1\) the coefficient is positive: the reflected wave is not phase-inverted. Then \(R_I = R^2 = 0.00849\) (0.85%) and \(T_I = 1 - R_I = 0.9915\) (99.15%). Most energy continues deeper, which is why soft-tissue boundaries give faint but useful echoes and why deep structures remain visible.

Solution 3

\(R = (0.0004 - 1.63)/(0.0004 + 1.63) = -0.9995\): the large negative value signals both a near-total reflection and a phase inversion. \(R_I = R^2 = 0.9990\), so 99.9% of the intensity reflects. Coupling gel excludes the air layer between probe and skin so that sound can enter the body at all.

Solution 4

Round-trip path is \(2\times8 = 16\) cm, so attenuation \(= 0.5\times5\times16 = 40\) dB and the fraction remaining is \(10^{-40/10} = 10^{-4}\): 0.01% of the intensity returns. This is why TGC (Section 3.2.3) is indispensable.

Solution 5

\(d = ct/2 = 1540\times91\times10^{-6}/2 = 0.0701\) m \(= 7.01\) cm.

Solution 6

  1. \(f_D = 2(4\times10^6)(0.6)(\cos60^\circ)/1540 = 1558\) Hz.
  2. PRF must exceed \(2f_D = 3117\) Hz.
  3. That PRF permits a depth of \(c/(2\,\mathrm{PRF}) = 1540/(2\times3117) = 24.7\) cm, comfortably more than a carotid study needs, so in this case the depth constraint is not binding. Contrast this with Solution 14, where it is.

Solution 7

\(v_{\max} = c\,\mathrm{PRF}/(4f_0\cos\theta) = 1540\times5000/(4\times5\times10^6) = 0.385\) m/s.

Solution 8

  1. \(Z_m = \sqrt{Z_{\text{PZT}}Z_{\text{tissue}}} = \sqrt{30\times1.6} = 6.93\) MRayl.
  2. \(\lambda_m = c_m/f = 2000/(5\times10^6) = 400\) µm, so \(t = \lambda_m/4 = 100\) µm.
  3. Without the layer, \(R_I = ((1.6-30)/(1.6+30))^2 = 0.808\): 80.8% of the intensity would reflect straight back into the element, which is why matching layers are not optional.

Solution 9

\(\mu = \rho c_s^2 = 1000\times6.25 = 6250\) Pa \(= 6.25\) kPa; \(E = 3\mu = 18.75\) kPa. A stiffness near 19 kPa is well above normal (a few kPa) and indicates advanced fibrosis, warranting evaluation for cirrhosis (exact thresholds are device-specific). Because \(E \propto c_s^2\), a 10% uncertainty in \(c_s\) becomes a 20% uncertainty in \(E\), the relative error doubles, which is why multiple acquisitions are averaged.

Solution 10

Convert: \(\mu_a = 4\) cm\(^{-1}\) \(= 400\) m\(^{-1}\); \(F = 20\) mJ cm\(^{-2}\) \(= 200\) J m\(^{-2}\). Then \(p_0 = \Gamma\mu_a F = 0.2\times400\times200 = 1.6\times10^4\) Pa \(= 16\) kPa. Units check: m\(^{-1}\cdot\)J m\(^{-2} =\) J m\(^{-3} =\) Pa.

Solution 11

  1. There are two unknown concentrations, so at least two independent equations, two wavelengths, are needed. More wavelengths over-determine the system and reduce noise sensitivity, provided they are chosen to keep the design matrix well conditioned (Lab 3). (b) At the isosbestic wavelength (\(\approx800\) nm) the two extinction coefficients are equal, so \(\mu_a \propto [\mathrm{HbO_2}] + [\mathrm{Hb}]\): the measurement reports total haemoglobin, hence blood volume, independent of saturation.

Solution 12

Nonlinear propagation builds a \(2f_0\) component that grows with depth and pressure, so it is strongest along the high-intensity central beam. Weak-signal sources of artifact — reverberation, side lobes, clutter, body-wall aberration, carry little energy and generate negligible harmonic. Filtering to the \(2f_0\) band therefore keeps the clean on-axis signal while rejecting these artifacts.

Solution 13

  1. \(\alpha = a f/8.686 = 0.5\times3/8.686 = 0.1727\) Np cm\(^{-1}\); \(\mu_I = 2\alpha = 0.3454\) Np cm\(^{-1}\).
  2. Intensity falls to \(1/e\) at \(1/\mu_I = 2.90\) cm.
  3. Pressure falls to \(1/e\) at \(1/\alpha = 5.79\) cm.
  4. Because \(I \propto p^2\), intensity decays twice as fast in the exponent as pressure does. The two depths are not alternative answers to the same question; they answer different questions, and stating which quantity is meant is mandatory.

Solution 14

  1. \(\mathrm{PRF}_{\max} = c/(2d_{\max}) = 1540/(2\times0.14) = 5500\) Hz.
  2. \(v_{\max} = c\,\mathrm{PRF}/(4f_0) = 1540\times5500/(4\times3\times10^6) = 0.706\) m/s.
  3. A 5 m/s jet at 14 cm demands \(v_{\max}d_{\max} = 0.70\) m\(^2\) s\(^{-1}\), while Eq. (23) caps the 3 MHz budget at \(c^2/8f_0 = 0.0988\) m\(^2\) s\(^{-1}\) — a factor of seven short. Raising the PRF beyond 5500 Hz would place the range gate shallower than the mitral valve, so no setting works: CW Doppler is required, at the cost of losing range information.

Solution 15

From \(v = cf_D/(2f_0\cos\theta)\), differentiating with respect to \(\theta\) at fixed \(f_D\) gives \(dv/d\theta = v\tan\theta\), hence \(\Delta v/v = \tan\theta\,\Delta\theta\) with \(\Delta\theta\) in radians. For \(\Delta\theta = 5^\circ = 0.0873\) rad: at \(20^\circ\), \(\tan20^\circ\times0.0873 = 3.2\%\); at \(70^\circ\), \(\tan70^\circ\times0.0873 = 24.0\%\). The error diverges as \(\theta\to90^\circ\) because \(\cos\theta\to0\) in the denominator of the Doppler equation. Keeping \(\theta\le60^\circ\) bounds the penalty near 15% for a \(5^\circ\) misestimate, which is the origin of the clinical rule.

Solution 16

  1. \(\mathrm{MI} = 2.2/\sqrt{2.5} = 1.39\).
  2. Yes, comfortably below the FDA limit of 1.9.
  3. No, not for contrast. Microbubbles undergo inertial collapse above MI \(\approx 0.4\), so an MI of 1.39 would destroy the agent within the imaging plane almost immediately. Routine CEUS is performed below MI 0.2; a high-MI pulse is used only deliberately, in flash-replenishment perfusion imaging.

Solution 17

  1. \(\Delta T = \mu_a F/(\rho C_p) = (400)(200)/(1060\times3600) = 0.0210\) K, about 21 mK.
  2. \(p_0 = \Gamma\mu_a F = 0.2\times400\times200 = 16\) kPa.
  3. The two quantities are governed by different material properties. Temperature rise depends on heat capacity, which is large for water-rich tissue; pressure rise depends on \(\Gamma = \beta c^2/C_p\), and the factor \(c^2 \approx 2.4\times10^6\) m\(^2\) s\(^{-2}\) is enormous. A tiny thermal perturbation therefore produces a substantial acoustic perturbation. This is precisely why photoacoustics is a safe, non-thermal technique that nonetheless generates a readily detectable signal.

Solution 18

  1. \(\mathrm{FR} = c/(2dN_{\text{lines}}N_{\text{foci}}) = 1540/(2\times0.16\times90\times2) = 26.7\) Hz.
  2. Marginally, it sits at the lower edge of the 25–30 Hz real-time band, adequate for wall-motion assessment at normal heart rates but not for a tachycardic patient.
  3. Dropping to a single transmit focal zone doubles the frame rate to 53 Hz at the cost of a narrower depth of field, so the beam is tight only near the single focus. Narrowing the sector is the alternative and costs field of view.

References and Further Reading

  1. Szabo, T. L. (2014). Diagnostic Ultrasound Imaging: Inside Out, 2nd ed. Academic Press.
  2. Cobbold, R. S. C. (2007). Foundations of Biomedical Ultrasound. Oxford University Press.
  3. Hill, C. R., Bamber, J. C., and ter Haar, G. R., eds. (2004). Physical Principles of Medical Ultrasonics, 2nd ed. Wiley.
  4. Prince, J. L. and Links, J. M. (2014). Medical Imaging Signals and Systems, 2nd ed. Pearson.
  5. Bushberg, J. T., Seibert, J. A., Leidholdt, E. M., and Boone, J. M. (2020). The Essential Physics of Medical Imaging, 4th ed. Wolters Kluwer.
  6. Duck, F. A. (1990). Physical Properties of Tissue: A Comprehensive Reference Book. Academic Press.
  7. Nyborg, W. L. (2001). Biological effects of ultrasound: development of safety guidelines. Ultrasound in Medicine and Biology, 27(3), 301–333.
  8. Church, C. C. (2005). Frequency, pulse length, and the mechanical index. Acoustics Research Letters Online, 6(3), 162–168.
  9. Sigrist, R. M. S., Liau, J., Kaffas, A. E., Chammas, M. C., and Willmann, J. K. (2017). Ultrasound elastography: review of techniques and clinical applications. Theranostics, 7(5), 1303–1329.
  10. Sarvazyan, A. P., Rudenko, O. V., Swanson, S. D., Fowlkes, J. B., and Emelianov, S. Y. (1998). Shear wave elasticity imaging. Ultrasound in Medicine and Biology, 24(9), 1419–1435.
  11. Tanter, M. and Fink, M. (2014). Ultrafast imaging in biomedical ultrasound. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 61(1), 102–119.
  12. Errico, C. et al. (2015). Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging. Nature, 527, 499–502.
  13. Wang, L. V. and Wu, H.-I. (2007). Biomedical Optics: Principles and Imaging. Wiley.
  14. Wang, L. V. and Hu, S. (2012). Photoacoustic tomography: in vivo imaging from organelles to organs. Science, 335, 1458–1462.
  15. Beard, P. (2011). Biomedical photoacoustic imaging. Interface Focus, 1(4), 602–631.
  16. Cox, B., Laufer, J. G., Arridge, S. R., and Beard, P. C. (2012). Quantitative spectroscopic photoacoustic imaging: a review. Journal of Biomedical Optics, 17(6), 061202.
  17. Prahl, S. A. Optical Absorption of Hemoglobin. Oregon Medical Laser Center compilation. https://omlc.org/spectra/hemoglobin/
  18. American National Standards Institute (2022). ANSI Z136.1: American National Standard for Safe Use of Lasers. Laser Institute of America.
  19. Dinov, I. D. (2021). Data Science: Time Complexity, Inferential Uncertainty, and Spacekime Analytics. De Gruyter.
  20. SOCR Probability and Statistics EBook. See also BPAD Chapters 1, 2, 5, 7, and 8.
SOCR/BPAD Resource Visitor number Web Analytics SOCR Email