1. Introduction

In linear elastic theory, it has been posed that knowledge of an object's surface deformation in response to mechanical stress is sufficient to uniquely solve for the internal elastic properties [1]. Solutions to this type of boundary value problem have been employed in seismology [2] and non-destructive testing [3] to map interior elastic properties of the earth and of structural materials, respectively. In many cases, elastography is limited by the inability to collect enough surface data to compute an accurate inversion. High frame rate holography is uniquely suited for capturing nanometer amplitude vibrations over a wide area, particularly in soft materials where surface acoustic waves have millimeter wavelengths at kilohertz frequencies. Here we demonstrate high resolution elastography by holographic imaging of surface acoustic (Rayleigh) waves at multiple interrogation frequencies.

The elastic properties of soft materials, such as polymers and human tissues, are of interest for their effects on electronic device performance [4], thermoplastics used in manufacturing [5], and medical diagnostics for breast cancer [6] and hepatic fibrosis [7]. In breast cancer, malignant breast tissue exhibits an elastic modulus 3-10 × greater than its healthy counterpart [6], motivating the development of biomedical imaging modalities that contrast elasticity. Current methods of biomedical elastography [8–13] are based on tracking interior tissue motion during mechanical excitation, which requires (3 + 1)-D data acquisition (3 spatial + 1 temporal dimension). To reduce the dimensionality of the data acquisition, we propose an elastography method based on high frame rate holographic imaging of Rayleigh waves, a (2 + 1)-D imaging system. By computing the solution to the boundary value problem, one may in theory extract all 3 spatial dimensions of elasticity.

As known in seismology [14] and non-destructive testing [15], Rayleigh waves provide depth-resolved elastography because their sensitivity kernel (i.e., the material space over which propagation is affected by elasticity) extends in depth on a scale proportional to the wavelength, as illustrated in Fig. 1 . In effect, lower frequency waves probe more deeply than higher frequency waves, so that the apparent wave dispersion provides depth-resolved elasticity information. The use of Rayleigh waves to quantify the elasticity of layers in soft materials has previously been performed by scanning a single detector away from a point source [16, 17]. Imaging the spatial pattern of surface waves on tissues has been accomplished by pulsed digital holography [18, 19] or dynamic speckle imaging [20], the latter of which was used to generate 2-D surface strain maps. We recently demonstrated imaging of the real-time propagation of Rayleigh waves with high frame rate holography [21], and showed that they qualitatively exhibit the correct dispersive behavior in layered media. It was also recently shown that high frame rate holography can visualize vibrations induced in skin by a bone conduction device [22].

 

Fig. 1 Cartoon diagram of depth-dependent out-of-plane (z) motion of Rayleigh waves as a function of wavelength (not to scale). Scanning the wavelength effectively scans the probing depth of the SAW. The amplitude of this motion is used as the sensitivity kernel in the proposed inverse method.

Download Full Size | PDF

To combine the 1-D depth resolution of Rayleigh waves with 2-D imaging, here we report, for the first time, the application of an inverse method for quantitative elastography by using a sensitivity kernel equal to the Rayleigh wave amplitude in a homogeneous half-space. From multi-frequency holography data we reconstructed depth profiles that represent 1-D axial elastograms, while scattering in the 2-D transverse plane is observed in response to a subsurface elastic inhomogeneity only when it is within the sensitivity kernel. We perform this method in an experimental context of soft biological tissues to demonstrate its validity for biomedical applications, by investigating phantoms matched to the elastic modulus of breast tissue. Our previous work [21] provides a technical basis for the holographic imaging method and the qualitative features of SAW propagation in layered and inhomogeneous phantoms. In this work we demonstrate, for the first time, quantitative elastography in layered phantoms, SAW propagation in a breast phantom with a 5mm tumor-like inclusion, and well-behaved SAW propagation on an ex vivo tissue sample.

2. Methods

2.1 Theoretical model and inverse method

In a homogeneous, elastic half-space, Rayleigh waves exhibit a phase velocity of [23]:

We propose an alternate approach that uses a Rayleigh wave sensitivity kernel equal to the depth-dependent Rayleigh wave amplitude in a homogeneous medium, specifically, the z component of particle displacement, which takes the form of a double inverse exponential [23]. We hypothesize that Rayleigh waves at a given frequency f in a medium with elastic depth profile, E(z), behave as they would in a homogeneous medium with an effective Young’s modulus, E_eff(f), given by the average of E(z) weighted by the sensitivity kernel for the effective homogeneous medium. Thus, our forward model has the form:

For forward computation, Eq. (2) can then be solved iteratively by substituting E_eff(f) into Eq. (1) to obtain the dispersion curve c_R(f) for a given E(z). This model shares similarities with other forward models for an effective c_R in layered media [15, 25, 26], but these require perturbative variations in E(z). We chose to linearize in E instead of $c_R\propto \sqrt{E}$ because the conditions of continuity for elastic waves are linear in E, analogous to permittivity for electromagnetic waves. As will be shown in below, our forward model can predict the dispersion curve c_R(f) even in highly non-perturbative tissue phantoms.

From this model, we can now straightforwardly solve the inverse problem of determining the elastic profile E(z) from an experimentally measured dispersion curve c_R(f). We rewrite Eq. (2) in a discrete vector-matrix form, E_eff = AE, with E being the elastic depth profile to be determined at each discrete depth z, and E_eff a vector representing the measured effective elasticity obtained by inversion of Eq. (1) using the c_R(f) measured at multiple frequencies f. The observation matrix A is the sensitivity kernel k(f,z) integrated over each discrete depth interval and normalized. Due to the Laplace transform nature of Eq. (2), the inverse problem is highly unstable, and regularization is needed. We find that Total Variation Regularization [27, 28], which assumes a piecewise-continuous solution, works best with our tissue phantoms. That is, the vector E is obtained by minimizing:

2.2 Experimental setup

An in-line, image plane, digital holography system, illustrated in Fig. 2(a) , was used to image the SAWs. Briefly, the system was composed of a Mach-Zehnder interferometer and low-power, long coherence length laser. The beam from a 633 nm, 5mW Helium Neon laser (Thorlabs Inc, Newton, NJ) was passed through a 92:8 beam splitter, with 92% of the power directed along the sample path. The sample beam was then passed through a 2.5 × beam expander composed of two lenses, with the lens spacing adjusted so that the beam was slightly divergent. A downward reflecting mirror illuminated an approximately 30 mm diameter area of the sample, and the diffusely reflected light was directed through a 200 mm focal length, 2 in diameter imaging lens, forming an image onto the sensor of a high speed CMOS camera. A 1:1 magnification was used to image a 17 × 17 mm area of the sample onto the 1024 × 1024 pixel array. An aperture in front of the imaging lens was used to adjust the subjective speckle size to be within the spatial resolution of the camera. The beam along the reference arm was passed through a 14 × beam expander, and coupled through a spatial filter (a converging lens combined with a 15 μm pinhole). The reference and sample beams were combined on-axis, forming an interferogram on the camera sensor. A variable neutral density filter was used in the reference arm to adjust the intensity, and optimize the interference pattern within the dynamic range of the camera. The path lengths of the reference and sample arms were carefully adjusted to be well within the coherence length of the laser (~200 mm). As waves propagated across the phantom surface, out-of-plane (z) surface displacement induced an optical path difference producing a phase shift in the interferogram, while in-plane motion was not observed as it was small compared to the subjective speckle size (typically 30-50 μm).

 

Fig. 2 (a) Digital holography system for monitoring Rayleigh waves on tissue phantoms. f₁ = 40mm, f₂ = 100mm, f₃ = 11mm, f₄ = 150mm, f₅ = 200mm. (b) Raw image (interferogram) obtained by the high speed camera (SAW frequency = 96Hz, camera framerate = 3kHz (Media 1). (c) The corresponding phase reconstruction (Media 2).

Download Full Size | PDF

A piezo-electric transducer (Kinetic Ceramics Inc, Hayward, CA) with a maximum travel range of 20 μm was placed just outside the illumination area to excite SAWs on the sample. The contact area between the transducer and sample is considered a point, as its diameter was much smaller than the SAW wavelengths observed. The waves were allowed to reach a steady state (~2-3 seconds after excitation) before the interferograms were recorded. The piezo-electric transducer was sinusoidally driven at various frequencies between 70 Hz and 300 Hz. The framerate of the camera was chosen to be at least 10 × the frequency of the SAW, at values of 1 kHz, 2 kHz or 3 kHz, corresponding to respective array sizes of 1024 × 1024, 640 × 640 and 512 × 512 pixels. This resulted in 17 × 17 mm, 10.6 × 10.6 mm and 8.5 × 8.5 mm imaged areas of the sample, respectively. A typical interferogram obtained by the system is illustrated in Fig. 2(b). Note that the fringes in this image are attributed to a cover glass plate in front of the sensor; these fringes are stationary and do not provide information about the SAWs. Instead, the interference between reference and sample beams can be observed in Media 1 as the slowly changing speckle pattern in response to the SAWs moving across the surface.

The series of images recorded by the camera was then transferred to a computer, where a phase reconstruction program was used to first bandpass filter the time-series of interferograms, and then apply a temporal phase shifting algorithm detailed in Ref. 21. For sinusoidal excitation of the surface, the reconstructed images using this algorithm represent the phase of the surface acoustic wave, modulo π, as illustrated in Fig. 2(c). Note that we do not track SAW amplitude, as only the phase velocity c_R is needed for elastographic mapping. The reconstructed SAW phase maps are highly contrasted along wavefronts where the phase wraps. We measured the shortest distance between successive wavefronts, d = λ/2, and computed the speed of the SAWs according to c_R = 2df. The mean and standard deviation was computed for measurements at 10 locations for each data set. The experiment was repeated at different frequencies to obtain a dispersion curve, c_R(f), for each sample.

2.3 Silicone phantom preparation

To test the accuracy of the proposed inversion technique, silicone elastomer phantoms with optical and mechanical properties that mimic real tissue were prepared according to previously described methods [30]. To perform experiments in a context of human breast imaging, we chose a Young’s modulus 4.15 kPa to mimic that of healthy adipose or fibroglandular human breast tissue, and a Young’s modulus 31.9 kPa to mimic invasive ductal carcinoma, according to mechanical measurements of fresh breast tissues reported in [6]. We will respectively refer to these Young’s moduli as “soft” and “stiff.”

The phantoms were fabricated from a mixture of pure non-cross-linking polydimethylsiloxane (PDMS) fluid (50 cSt viscosity, ClearCo Products, Bensalem, PA), and a cross-linking room-temperature PDMS base and its associated curing agent (Sylgard 184 A and B, respectively, Dow Corning, Inc., Midland, MI). Adjusting the ratio of the pure PDMS fluid to the PDMS base was used to tailor the Young’s modulus. “Stiff” phantoms were obtained using a 5:1 ratio, while “soft” phantoms were obtained using a 12:1 ratio. The ratio of the PDMS base to its curing agent was held constant at the recommended 10:1. Titanium dioxide (TiO₂) powder (Sigma-Aldrich, #224227, particle size <5 μm, St. Louis, MO) was added to the mixture at a concentration of 4 mg/g to function as tissue-like optical scattering agents. The Young’s modulus for each ratio was measured using a TA.XT Texture Analyzer (Texture Technologies Corp., Scarsdale, NY). To prepare phantoms for imaging, mixtures were poured into cylindrical, plexi-glass containers of diameter 126 mm and height 60 mm, and allowed to cure for 24 hours in an isotemp oven at 80°C.

First, two-layer phantoms were prepared to test the ability to extract E(z) by our inverse method. One phantom was composed of a “soft,” 15 mm layer on top of a “stiff,” 40mm layer, while a second phantom was identical to the first but with the layers switched. These were prepared in two steps. First, the 40 mm lower layer was mixed and allowed to cure in the container. The 15 mm upper layer was then poured over the cured lower layer, and the sample was then allowed to cure for another 24 hours.

Next, a phantom was prepared to simulate the mechanical properties of a tumor-like inclusion embedded inside otherwise healthy breast tissue, to test whether SAWs may be sensitive to tumors based on their mechanical properties. The inclusion was approximately cylindrical with a diameter of 5mm and height of 5mm, and embedded at a depth of ~9 mm below the surface of the phantom. This phantom was prepared in four steps. The first 40 mm, “soft” lower layer was cured in the phantom container. Simultaneously, a 35 mm diameter by 5 mm height “stiff” cylindrical sample was mixed and cured. A 5 mm height cylindrical piece was then cut from this sample to act as the tumor inclusion. This piece was then laid on top of the “soft” layer, surrounded by another 5 mm layer of “soft” mixture, and allowed to cure. Finally, the top-most “soft” layer of approximately 9 mm height was poured on top and allowed to cure. The final result was thus an almost homogeneous phantom of Young’s modulus of 4.15 kPa, with an inclusion in the center at a depth spanning from 9 to 14 mm, with Young’s modulus of 31.9 kPa. B-mode ultrasound images were then acquired of the phantom to verify the location of the inclusion before holographic imaging. The ultrasound system was a SonixTOUCH (Ultrasonix Medical Corporation, Richmond, BC, Canada).

3. Results and discussion

3.1 Elastic depth profiling of two-layer phantoms

First, we investigated the ability to extract E(z) using our model in section 2.1 by imaging SAWs on two-layer phantoms with known Young’s moduli. The experimentally obtained dispersion curves c_R(f) for each two-layer phantom are illustrated in Fig. 3 . As described in our earlier work [21], the slope of the dispersion curve depends on the relative stiffness change for increasing depth. From Eq. (1), it can be seen that stiffer materials (i.e. higher Young’s modulus) exhibit a higher phase velocity than softer materials. Thus, the SAW speed at higher frequencies (i.e. lower wavelength, and therefore shallower probing) tends towards what would be expected for a homogeneous phantom composed purely of the material in the upper layer. Conversely, at lower frequencies, the speed of the SAW tends towards that of a phantom composed of the material in the lower layer. Thus, for the phantom with a “soft” upper layer and “stiff” lower layer, the dispersion curve has a negative slope, while the “stiff” upper layer and “soft” lower layer phantom dispersion curve has a positive slope.

 

Fig. 3 Dispersion curves obtained for samples with (a) stiff lower layer and (b) stiff upper layer. The corresponding inversions of the dispersion curve in (a) and (b) are shown in (c) and (d), respectively.

Download Full Size | PDF

The use of phantoms with known Young’s moduli now allows us to test our inverse method for quantitatively extracting E(z). First, we compared the measured dispersion curves with those computed using our forward model (Eq. (2)) and the known distribution of E(z). While the experimental results seem to produce consistently larger values of c_R, the overall shape of the dispersion curves is in good agreement for both phantoms. Next, we performed the inverse method upon the measured dispersion curves to predict E(z), the results of which are illustrated in Figs. 3(c) and 3(d) in comparison with the known E(z) according to the texture analyzer. For both phantoms, the Young’s modulus distribution obtained by inversion matches the true distribution fairly well (< 2.5 kPa on average), and the transition depth between the layers was accurately predicted within one sampling point (Δz = 2.6 mm). For the phantom with the soft upper layer, the inversion shows a transition from the upper layer to lower layer between 15.7 and 18.3 mm (compared to the true value of 15 mm). The inversion result for the phantom with stiff upper layer shows a tighter match to the true distribution, with a transition between the two layers predicted between 13.1 and 15.7 mm.

3.2 SAWs on a homogeneous phantom with a tumor-like inclusion

Digital holographic imaging of the homogeneous phantom containing a tumor-like inclusion is illustrated in Fig. 4 . To verify the location of the inclusion, a B-mode ultrasound image of a region containing the inclusion is shown in Fig. 4(a), in comparison to a region to one side of the inclusion (Fig. 4(e)). The horizontal lines in the ultrasound images correspond to the separately cross-linked layers of the phantom during the preparation process. The corresponding Rayleigh wave phase reconstructions, obtained by our digital holography method, are shown in the x-y plane next to each corresponding x-z ultrasound image in Fig. 4. Note that the frame rate was increased from 1 to 3 kHz as the excitation frequency was increased from 96 to 228 Hz, impacting the noise level on the phase reconstructions.

 

Fig. 4 Images of a breast tissue phantom with a tumor-like inclusion. (a) B-mode ultrasound of a region with the inclusion (marked with dashed circle). SAW phase reconstructions of the surface over the inclusion are shown in (b) for excitation at 96Hz (Media 3) (c) for 168Hz (Media 4), and (d) for 228Hz (Media 5). The dashed circles in (b-d) correspond to the approximate lateral position of the inclusion. SAWs are perturbed by the inclusion only at low frequencies, when the SAW penetration depth is great enough to reach the inclusion. (e) B-mode ultrasound of the same phantom in a control region without the inclusion. SAW phase reconstructions of the surface in this region are shown in (f) for 96Hz (Media 6), (g) for 168Hz (Media 7), and (h) for 228Hz (Media 8).

Download Full Size | PDF

As can be seen, the Rayleigh wavefronts are increasingly perturbed at lower frequencies in the region containing the inclusion. The higher frequency waves and the waves not passing over the inclusion, on the other hand, remain unperturbed, appearing like waves in an elastically homogeneous medium. This is because the Rayleigh wavelengths in the background region decrease from 12 mm at 96 Hz to 5 mm at 228 Hz, so that the sensitivity kernel, which scales as the wavelength, does not extend significantly into the inclusion for the highest frequency waves. Conversely, the lower frequency waves, which have a deeper sensitivity kernel, exhibit scattering and refraction as they pass over the inclusion.

3.3 Imaging of an ex vivo tissue sample

To test the viability of this technique in real tissue, a raw chicken thigh was obtained for holographic imaging. This ex vivo tissue was chosen because of its geometry: the thigh had an ellipsoidal shape, with an approximate length of 105 mm, width of 68 mm and a depth of 40 mm, with a bone running through the middle along the length of tissue at a depth of 20 mm, measured using calipers. The bone, which is much stiffer than the surrounding tissue, is approximately a cylindrical rod of diameter 10 mm. An ultrasound image illustrating the location of the bone and corresponding digital holographic Rayleigh wave phase reconstructions at various frequencies are shown in Fig. 5 .

 

Fig. 5 Images of a raw chicken thigh with a bone through the middle. Ultrasound images are shown in (a), with the bone indicated by the dashed circle. SAW phase reconstructions of SAWs over the region containing the bone are shown in (b) for 84Hz (Media 9) (c) for 168Hz (Media 10), and (d) for 432Hz (Media 11). High frequency SAWs do not penetrate to the depth of the bone and appear more circular.

Download Full Size | PDF

Similar to the results obtained with the phantom containing an inclusion, the higher frequency Rayleigh waves are more circular and regularly spaced, with little apparent scattering, suggesting that the tissue at this depth scale (λ_R ≈11 mm) is elastically homogeneous. The c_R at 432 Hz was approximately 4.78 m/s, corresponding to an E_eff of 72 kPa in the soft tissue region above the bone. The lower frequency Rayleigh waves, on the other hand, have a deeper sensitivity kernel and are highly deformed, consistent with the presence of the much stiffer, subsurface bone. Overall the Rayleigh waves appear to be well-behaved in soft tissue, highlighting the feasibility of our approach for tissue elastography.

4. Conclusion

The results shown here demonstrate the potential for a novel elastography method using digital holography to image surface acoustic waves. In two-layer phantoms we were successfully able to reproduce the 1-D elastic depth profiles, up to a depth of 55 mm, using an ad hoc model of Rayleigh wave dispersion behavior. On average, the breast size in the U.S. population would require an imaging depth between 40 and 80 mm to reach the chest wall without pre-compression [31]. To obtain this depth at average breast elasticity [6] would require sampling SAWs between 12 and 25 Hz, which could be obtained in our experimental setup by more onboard memory or by a lower available frame rate than 1 kHz. Importantly, this method provides greater penetration depth than other optical elastography methods reported to date [11–13], as it is not limited by the penetration depth of light into the tissue, but rather that of the Rayleigh wave.

It was also shown that real-time monitoring of Rayleigh wave propagation provided the ability to identify the approximate transverse position and depth of a tumor-like inclusion based on multi-frequency data. Thus, while the frequency-dependent behavior of SAWs contains information about the elastic properties as a function of depth, the transverse propagation contains information about the lateral profile. In future work, the application of a more complete elastic scattering model of Rayleigh wave behavior may allow one to obtain 3-D elastograms, through a similar inverse problem approach.

In this demonstration, importantly, the SAW propagation exhibited a well-behaved profile in both phantoms and real tissues that should be amenable to quantitative analysis. The ultimate goal is for the approach to be used to study in vivo human tissue. Before this can be achieved, some experimental and theoretical issues need to be addressed. For example, the effect of tissue motions other than that due to SAWs will need to be investigated. In this work, we found that by bandpass filtering at the excitation frequency, the system became more robust to motion artifacts. Chirped (frequency-swept) SAW excitation may also enable more rapid acquisition to minimize motion artifacts. We note that the nanometer displacement sensitivity of digital holography affords measurements using low strains and strain rates ($\epsilon$ < 10⁻⁵ and $\dot{\epsilon }$ < 0.02 s⁻¹, respectively), which may respectively minimize artifacts due to nonlinear and viscoelastic mechanical properties of breast tissue. Another need will be the use of a larger imaging area and multiple transducers for whole breast imaging. The holography system is easily adaptable to larger field-of-view by modifying the imaging optics, as long as sufficient resolution to observe the smallest desired Rayleigh wavelengths (typically 3-4 mm) is maintained. Furthermore, wide field imaging provides a more complete set of data that leads to high resolution elastic reconstructions.

This elastography technique constitutes a novel, non-invasive approach that takes advantage of the flexibility, nanometer-scale displacement sensitivity, and wide area imaging provided by high frame rate digital holography. Furthermore, the use of SAWs provides depth-dependent elasticity information up to several centimeters into soft materials, which is of high relevance for investigating soft tissues and for cancer detection applications.