2.1.

Tissue-Phantom Preparation

Tissue phantoms with optical and mechanical properties comparable to that of soft tissues were used in these experiments, similar to those used previously.³⁵ For the homogeneous phantoms, they were fabricated from a mixture of pure non–cross-linking polydimethylsiloxane (PDMS) fluid (50-cSt viscosity, ClearCo, Inc., Bensalem, Pennsylvania), and a cross-linking room-temperature PDMS base and its associated curing agent (Sylgard 184 A and B, respectively, Dow Corning, Inc. Midland, Wisconsin). The Young's moduli of the phantoms were adjusted over the range from 2 to 40 kPa by changing the ratio of the pure PDMS and the PDMS base over the range from 16:1 to 5:1, whereas the ratio between the PDMS base and its curing agent was held constant at 10:1. The relation between the mixture ratio and Young's modulus was estimated based on previous work.³⁶ Titanium dioxide (TiO₂) powder [Sigma-Aldrich Co., LLC (St. Louis, Missouri), No. 224227, mean size 1 μm, <5 μm] was embedded at a concentration of 4 mg/g in the tissue phantoms to function as tissuelike optical scatterers for holographic imaging. The phantom solutions were mixed thoroughly on an ultrasonicator for several tens of minutes at room temperature and then poured into a 126-mm-diam cylindrical container, with a thickness of ∼30 mm. The phantom dimensions were chosen to be large compared to the attenuation length of Rayleigh waves to eliminate significant boundary reflections. To further eliminate the boundary reflections, a few millimeters of the phantom adhering to the sidewall were removed. All phantoms were cured in an oven at 80°C for one day and subsequently at room temperature for another day. To prepare the two-layer phantoms, the same process was repeated for the top layer. The solution mixture was poured onto the cross-linked lower layer, and the whole phantom was cured as for the homogeneous phantoms. To construct a phantom with a local inhomogeneity, we prepared two homogeneous phantoms with the same mixing ratio and stacked them together, between which a thin metal ring with a 6.6 mm o.d., 3 mm i.d., and thickness of 0.5 mm was sandwiched.

2.2.

Digital Holographic Interferometry

Figure 1 shows the in-line, image-plane holographic interferometry system for monitoring SAWs. A 633-nm helium-neon laser of 5 mW power (Thorlabs) was divided into sample and reference arms by a 92:8 beam splitter, where the majority of the optical power was directed to the sample arm. Two lenses in the sample arm acted as a 2.5× beam expander with the lens spacing adjusted so that the output of the last lens was slightly divergent, illuminating the phantom with a beam nominally 30 mm in diameter. A 2–in.-diam mirror (downward reflecting mirror in Fig. 1) directed the beam downward onto the horizontal surface of the phantom and subsequently directed the diffusely reflected light from the phantom through a 200-mm focal length, 2–in.-diam imaging lens, which imaged the light onto the sensor of a fast CMOS camera (Photron SA3). The phantom was mounted on a vertical stage with which the height was adjusted to focus the image. The camera frame rate was adjusted between 1 and 7.5 kHz to sample images at a rate ≥10× the frequency of the SAW. At these frame rates, the corresponding pixel array size varied between 1024 × 1024 and 256 × 256, respectively. The pixel size is 17 μm, and the imaging area for the full 1024 × 1024 pixel array was 16 × 16 mm². Because there are highly spatially modulated speckles in the image, an adjustable aperture was placed close to the imaging lens. The aperture size was adjusted until the speckles could be distinguished by the camera sensor. The reference arm beam was passed through a spatial filter (i.e., the combination of a short focal length lens and a pinhole) and a 150-mm focal-length positive lens (f₄ in Fig. 1), before being directed into the camera. The distance from the pinhole to the lens f₄ was adjusted so that interference fringes between the reference and sample could be identified over the whole measurement area. Because the optical path of an edge ray from the sample to the camera is greater than that of a central ray, the output light from the lens f₄ was chosen to be divergent. Ideally, the optical path difference (OPD) between edge rays and central rays should be balanced for the two arms. A neutral density filter was inserted in the reference arm to optimize the interference signal within the dynamic range of the camera. A retroreflector was used in the reference arm to adjust the OPD between the two arms so that it was within the coherence length of the laser, ∼200 mm.

Fig. 1

Schematic diagram of holographic elastography experimental setup. NDF: neutral density filter, f₁ = 40 mm, f₂ = 100 mm, f₃ = 5 mm; f₄ = 150 mm, and f₅ = 200 mm.

The image of the phantom interacted with the reference beam, forming a local holographic image. As SAWs propagated across the phantom surface, surface displacement along the axis that bisects the propagation directions of the incident and detected scattered fields (approximately the axis normal to the phantom surface) induced a phase shift in the interferogram due to the changing OPD. As such, this system was predominantly sensitive to out-of-plane sample motions. For generation of SAWs, a piezoelectric transducer (Kinetic Ceramics, Inc., Hayward, California) with a travel range of 20 μm was placed into contact with the phantom immediately to one side of the illumination area. The transducer contact area could be regarded as a point because its diameter is much less than the SAW wavelengths in our experiments. The transducer was driven by a high-voltage power supply controlled by a function generator. Typically, the transducer was initiated at a fixed frequency and the system was allowed to reach a steady state (∼1−2 s), after which the camera recorded the holographic images in real time. On-board image data were transferred to the computer for further processing, as described later.

2.3.

SAW Phase Extraction

The SAWs propagate outward from the point of actuation, which cause temporally modulated phase shifts in the holographic image. When the phantom's thickness is much greater than the SAW wavelength, we expect the SAW to be characterized as a Rayleigh wave.³³ In that case, the displacement in the z direction (out of plane) can be expressed as [TeX:] $\Delta z\left({\vec r,t} \right) = A\cos ({\omega t + \vec k \cdot \vec r})$ $\Delta z\left(\vec{r},t\right)=A\cos \left(\omega t+\vec{k}·\vec{r}\right)$ , where A is the SAW amplitude (assumed to be slowly varying in [TeX:] $\vec r$ $\vec{r}$ , roughly as the inverse square root of [TeX:] $\left| {\vec r} \right|$ $\left|\vec{r}\right|$ ), [TeX:] $\vec r$ $\vec{r}$ is the position in the (x,y) surface plane, ω is the angular frequency of the SAW, and [TeX:] $\vec k$ $\vec{k}$ is the spatial frequency of the SAW. The intensity recorded on the camera can then be written as follows:

To extract the phase of the SAW, that is, [TeX:] $\beta _{{\rm SAW}} \left({\vec r,t} \right)\break = \omega t + \vec k \cdot \vec r$ ${\beta }_{\mathrm{SAW}}\left(\vec{r},t\right)=\omega t+\vec{k}·\vec{r}$ , in the limit of small SAW amplitudes (Δz ≪ λ) and sufficiently fast temporal sampling (ωΔt ≪ 2π, where Δt is the time between successive frames), it is possible to use a temporal phase-shifting algorithm similar to Refs. 37 and 38, such that

where a group of three successive frames at times t – Δt, t, and t + Δt are processed to produce the SAW phase image at time t, and [TeX:] $\mathord{\buildrel{\lower3pt \hbox{\scriptscriptstyle\rightharpoonup}}\over r}$ $\stackrel{\rightharpoonup }{r}$ is the (x,y) pixel position in the frame. Briefly, we can understand why this works by rewriting the phase-shifting algorithm in terms of time derivatives of I and using the small angle approximation,

Plugging these back into Eq. 3 reveals how the extracted phase β_SAW using the phase-shifting algorithm in Eq. 2 is equivalent to the phase of the mechanical wave, [TeX:] $\omega t + \vec k \cdot \vec r$ $\omega t+\vec{k}·\vec{r}$ . It should be noted that our method is nontraditional in that the reference beam is not phase shifted and the optical phase is not obtained directly; instead, we employ a priori knowledge of the mechanical wave frequency to obtain the mechanical phase. This has the advantage of skipping the unnecessary evaluation of the optical phase, in favor of directly assessing the mechanical phase. One consequence is that we do not know at any instant in time whether the SAW is in its upswing or downswing and, therefore, we can only measure β_SAW modulo π. In our phase images, therefore, we represent β_SAW = –π/2 to π/2 as a gray-scale value from black to white, respectively.

Figure 2a (Video 1) shows an example interference pattern measured by the camera, [TeX:] $I({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}}\over r},t})$ $I\left(\stackrel{\rightharpoonup }{r},t\right)$ . In the video, the evenly distributed fringes aligned from lower right to upper left are background fringes caused by the protection glass in front of the camera sensor due to multiple reflections (Fabry–Perot). These are stationary and do not affect the measurement result using the temporal phase-shifting algorithm. The SAW phase image, [TeX:] $\beta _{{\rm SAW}} ({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}}\over r},t})$ ${\beta }_{\mathrm{SAW}}\left(\stackrel{\rightharpoonup }{r},t\right)$ , which was extracted from Fig. 2a according to Eq. 3, is shown in Fig. 2b (Video 2). It reveals a SAW propagating away from the actuator position located outside of the image to the right. Note that the phase-wrap position from β_SAW = π/2 to –π/2 provides a highly contrasted way to visualize the SAW wavefront, which we also verified with simulations. We then obtain the SAW phase velocity by tracking a specific wavefront in the image series. Although we found human inspection to be sufficient for SAW velocity measurement within the scope of this work, in future work we plan to develop an automated method for SAW velocity extraction.

Fig. 2

Images of a phantom during mechanical excitation by a piezoelectric transducer with a driving frequency of 180 Hz. The camera frame rate is 2 kHz, and the measurement area is 10 × 10 mm². (a) Holographic images recorded by the camera (Video 1) and (b) the associated extracted SAW phase according to the methods in Sec. 2.3 (Video 2). (Video 1, QuickTime, 356 KB); (Video 2, QuickTime, 472 KB)

It is important to note that modulation of the speckle is the dominant effect giving rise to the SAW phase signal. This is because the phantoms were designed so that the diffuse scattering from the TiO₂ micropowder dominates over the specular reflection from the surface in order to more accurately mimic the optical signal from real tissues. This experimental configuration can actually be thought of as a specific case of digital holography known as electronic speckle pattern interferometry.^{37, 39, 40}

3.1.

Surface Acoustic Waves in Homogeneous Phantoms with Varying Young's Moduli

Our hypothesis is that SAWs on a soft-tissue–like, homogeneous, elastic half-space will behave as Rayleigh waves.³³ consistent with other observations.²³ Here, we study this phenomenon using digital holographic interferometry. We generated tissue phantoms with Young's moduli ranging from 2 to 40 kPa to match the range of elastic moduli encountered in soft tissue such as the human breast.²⁹ According to elastic wave theory, the phase velocity of Rayleigh waves in a homogeneous half-space depends on the Young's modulus (E) and density of the material (ρ) according to²³

To produce SAWs, we excited the phantoms with a 150-Hz sinusoidal signal on the piezoelectric transducer. This choice of frequency was such that the corresponding SAW wavelength (and, thus, penetration depth) was smaller than the phantom depth (30 mm). This produces a surface wave with constant phase velocity due to the elastic homogeneity of the material. Using the SAW phase-extraction method previously described, we quantified the wave speed as a function of Young's modulus and plotted the results in Fig. 3. Each data point was obtained by averaging eight measurements for a single sample. The solid line is the theoretical curve calculated with Eq. 6 using ρ = 960 kg/m³. The dashed line is a least-squares best-fit curve with one fitting parameter, k ₀, such that c _R = k ₀√E. The fitting value of k ₀ = 0.525 is within 6% of the theoretical value of 0.560. The R ² value of the fit was 0.99. For additional validation, we measured the SAW speeds for one homogeneous phantom while sinusoidally driving the transducer at different frequencies (75−300 Hz) and found that the measured velocity remained constant with a standard deviation of ∼7.5% of the mean. This demonstrates that the induced out-of-plane SAW motions measured by holographic interferometry are well described as Rayleigh waves.

Fig. 3

Variation in SAW speed as a function of the Young's modulus in homogeneous tissue phantoms follows the expected relationship of Eq. 6.

3.2.

Surface Acoustic Wave Velocity Dispersion in Two-Layer Phantoms

Because the penetration depth of Rayleigh waves is proportional to the wavelength, it is known that a variation in the material's elasticity versus depth gives rise to dispersion in the SAW velocity.^{17, 18, 19, 20, 21} We hypothesize that SAWs in soft tissues will exhibit a velocity according to Eq. 6 given by the average material elasticity E within the surface depth probed by the wave. In other words, we expect the long-wavelength (low-frequency) SAWs to probe the tissue response over a large depth range, while the short-wavelength (high-frequency) SAWs should probe only the upper tissue layers.

In order to test this hypothesis, we generated two two-layer phantoms, one with a soft upper layer and stiff lower layer, and a second with a stiff upper layer and soft lower layer. For the first phantom, the Young's modulus E of the upper layer was 7 kPa [corresponding to a theoretical velocity of 1.4 m/s according to Eq. 6], and the lower layer was 22 kPa (corresponding to a velocity of 2.6 m/s). The top layer was 4 mm thick, and the bottom layer was 10 mm thick. The measured SAW speed at different driving frequencies is shown in Fig. 4a, where the dashed lines indicate the theoretical velocities in a homogenous medium composed of the upper or lower layer. For the second phantom, the Young's modulus E of the upper layer was 38 kPa (corresponding to the velocity of 3.4 m/s) and the lower layer was 11 kPa (corresponding to the velocity of 1.8 m/s). The top layer was 6 mm thick, and the bottom layer was 30 mm thick. The measured SAW speed at different driving frequencies is shown in Fig. 4b.

Fig. 4

SAW speed variation with driving frequency for (a) a phantom with upper softer layer and (b) a phantom with upper stiffer layer.

From the two plots, it can be seen that the SAW velocity gradually transitions from close to the velocity of the lower layer to that of the upper layer with the increase of the driving frequency. This expected behavior is observed in both cases, with the stiffer layer either being on top or underneath, such that a positive dispersion is indicative of decreasing stiffness versus depth, and a negative dispersion that of increasing stiffness versus depth. In future work, this measured dispersion profile can be used to quantitatively predict the depth-dependent elasticity either in the perturbative regime^{18, 19} or by treating the exact solution as an inverse problem.^{17, 20}

Qualitatively, we can compare the transition wavelength between lower and upper wave velocities to the layer thickness. This is obtained at the point at which the slope dv/df is the maximum and is observed to occur at approximately 100 and 450 Hz in each of the phantoms, respectively. The associated SAW wavelengths, using the velocities measured at each of these frequencies, are 18 and 7 mm, respectively. In the first phantom with a 4-mm-thick soft upper layer, the transition wavelength of 18 mm is much larger than the layer thickness, suggesting significant decay of the wave has reduced its effective probing depth. In the second phantom with a 6-mm-thick stiff upper layer, the transition wavelength was ∼7 mm. Depending on Poisson's ratio, out-of-plane Rayleigh wave motion is reduced to ∼20–33% of its surface value at one wavelength in depth.³⁴ We might therefore expect the transition to occur when the wavelength is on the order of the upper layer thickness. However, comparing the transition wavelength with upper layer thickness in the two phantoms suggests that the softer phantom material exhibits greater decay than the stiffer material. This suggestion is also supported by the observation that the SAW velocity in the phantom in Fig. 4b does not quite reach that of the soft lower layer at low frequency. If the attenuation in the soft layer is too rapid, then the wave will tend to bias toward a velocity given by the stiffer upper layer. Further investigation is needed to determine whether the softer phantom material exhibits more dissipative loss, which might be expected because the softer phantom contains a lower concentration of cross-linking PDMS (which is dominantly elastic) relative to the non–cross-linking PDMS fluid (which is dominantly viscous).

In experiments, we also note that the observable range of SAW frequencies (frequencies where SAWs had sufficient amplitude to provide a detectable phase variation) were limited in practice. No SAWs were observed at <20 Hz, while at >80 Hz, SAWs were observed in every experiment. We found the high-frequency limit to lie between 300 Hz for phantoms with a soft upper layer and 800 Hz for phantoms with a stiff upper layer. Improvements in phase sensitivity or actuator sweep may increase the detectable frequency range, in order to increase the depth penetration. Thus, we expect the spatial resolution for elastographic reconstruction to be fundamentally limited by the SAW bandwidth and signal-to-noise ratio. To speed the data acquisition, a chirped (frequency-swept) SAW waveform can be applied to rapidly measure the dispersion curve; the rate of chirp can be slowed to provide the desired frequency resolution, depending on the needed spatial resolution.

3.3.

Surface Acoustic Wave Propagation in Phantoms with a Local Inhomogeneity

To investigate the potential of this technique for detecting local inhomogeneities, we prepared an otherwise homogeneous soft phantom (7 kPa) with an embedded thin metal ring as described in Sec. 2.1. The embedding depth was 8 mm below the surface. Figure 5 shows the phase reconstruction of the SAW propagation at frequencies of 65 Hz (Video 3) and 89 Hz (Video 4) using a frame rate of 1 kHz, and 200 Hz (Video 5) using a frame rate of 2 kHz. These frequencies nominally correspond to wavelengths of 22, 16, and 7 mm, respectively. The measurement area was 16 × 16 mm² at 1 kHz and 10 × 10 mm² at 2 kHz. The effect of the ring is clearly observed as a deformation of the SAW propagation, and the shape of this deformation appears to change at different frequencies due to the varying penetration depth. In particular, at the frequency of 200 Hz, the impact of the ring is minimized because its wavelength of 7 mm is shorter than the ring depth of 8 mm, while the impact is larger at the shorter wavelengths. The impact appears to be largest at 89 Hz, although significant wavefront deformation is observed in the upper left corner of the images at 65 Hz as well. This demonstrates that this technique is sensitive to local inhomogeneities at least as deep as 8 mm. The capability to identify local inhomogeneities is useful because tumors and other disease states are often characterized by a palpable, local elastic abnormality.

Fig. 5

SAW propagation in a phantom with a local inhomogeneity at frequencies of 65 Hz (a) (Video 3), 89 Hz (b) (Video 4), and 200 Hz (c) (Video 5). The image at 200 Hz has a smaller field of view due to the increased frame rate. The dashed circle indicates the lateral position of the metal ring. (Video 3, QuickTime, 728 KB); (Video 4, QuickTime, 878 KB); (Video 5, QuickTime, 986 KB)