Since it was reported on in 1991,1 optical coherence tomography (OCT) has gradually become the tool of choice for noninvasive three-dimensional imaging of biological and industrial samples.2 Particularly important is its unique combination of noninvasive, noncontact, high spatial resolution, high sensitivity, and high speed, the features that have attracted increased interest in a wide range of preclinical and clinical applications, for example, in ophthalmology, cardiology, neurology, dermatology, etc.3 Traditionally, OCT is used for imaging tissue morphology. However, one of the recent developments has been in the use of OCT for high sensitive detection of tissue microscopic motion/deformation in vivo. Because the tissue motion is directly related to its biomechanical properties, the monitoring and imaging of this motion provides important information about tissue mechanical compliance.
The first report that used OCT to detect tissue micro-deformation was demonstrated under external compression load by Schmitt in 1998,4 introducing the term of OCT elastography (OCE). Since then, various techniques were augmented to OCT for the purpose of tracking tissue motion. In the early development of OCE, the method was based on a straightforward approach that employs intensity-based speckle tracking to provide tissue displacement under static or quasi-static mechanical load.5–7 However, the intrinsic limitations of the traditional speckle tracking algorithm make this method difficult to provide accurate quantitative elastography, largely due to speckle intensity and shape variations.8
To mitigate the problems in the use of speckle tracking, the phase information intrinsically available in the complex OCT signals has been explored to accurately measure minute tissue motion with high sensitivity down to a subnanometer scale.9–11 Thanks to the rapid development of Fourier domain OCT that significantly improved the image acquisition speed and sensitivity, it is now possible to measure and/or image the minute dynamic tissue motion in real time, making an important step toward in vivo applications. In this regard, increased attention has been paid to the use of OCT to track and monitor mechanical waves propagating within tissue, induced by either external12 or internal13 mechanical stimuli. The methods of external stimuli include the use of acoustic radiation force14 or simply by a contact mechanical vibrator.12,15 Recently, we introduced the use of phase-sensitive OCT (PhS-OCT) to sensitively detect surface acoustic waves (SAW) propagating on skin16 and cornea17 surfaces, leading to quantification of the elasticity of underlying tissues. Others have used PhS-OCT to detect the shear wave18,19 and longitudinal wave.20 The assumption here is that the phase changes detected by PhS-OCT represent the true localized motion within tissue. While this is true in the detection of blood flow velocity,21,22 we notice that the detected phase changes do not accurately represent the localized tissue motion, particularly when detecting the induced shear waves within tissue, leading to the measurements being severely contaminated by a motion artifact.
Image artifacts have been an important topic of research in nearly all medical imaging modalities because they may degrade image quality and cause inaccurate clinical interpretation of images. Artifacts in motion detection are even more critical in OCE because the error in the detected mechanical responses could lead to inaccurate inversion results of tissue elasticity.
The purpose of this paper is to describe the theoretical and experimental investigations of artifacts that can arise in the motion detection with PhS-OCT when mechanical stimulation is applied on sample surface or surface motion is induced unintentionally. While spectral-domain OCT (SD-OCT)23,24 and swept source OCT (SS-OCT)25,26 are based on similar spectral interferometric principles and can be used to detect motion using phase-resolved approach, the motion signal from these techniques could be contaminated by similar artifacts. We present the results that provide an understanding of the motion artifacts when detecting the tissue motion in the underlying structures. An effective method of compensating for the motion artifacts is proposed and demonstrated. The results emphasize the need for surface motion compensation when measuring the mechanical response for elastography or other biomedical applications related to tissue motion detection.
Phase-Resolved OCT Motion Detection
Figure 1(a) shows the schematic of a basic configuration of an SD-OCT system. A broadband light is used as the illumination source, which is split into reference beam and sample beam. Reference beam is delivered to a stationary mirror, while sample beam is focused onto the specimen. The coupler recombines the backscattered light from both reference and sample arms to form a spectral interferogram, which is subsequently detected by a line-scan camera-based spectrometer. Hence, each pixel on the CCD array detects optical power as a function of optical wavenumber, .
Optical frequency-domain imaging is accomplished with a swept frequency laser source, thus also termed as SS-OCT.26 As shown in Fig. 1(b), SS-OCT differs from SD-OCT in the light source used and the interference detection method, in which the reflected lights from reference mirror and sample pass through two circulators and recombine by a coupler before being detected by a balanced photodetector.
The photocurrent associated with the interference fringes at the detector can be expressed by27
In the case of SD-OCT, is the integration over A-scan exposure time, while in SS-OCT, is instantly sampled by a high-speed photodetector. For either case, a Fourier transformation of (realistically, discrete Fourier transformation) will give the complex-valued domain signal.27
As seen in Eq. (2), the phase of is proportional to . Hence, the difference in optical path length (OPL), i.e., the difference in , can be calculated from the phase difference of the complex signals in the domain. For the phase-resolved method,28 the instantaneous displacement of a local scatterer between two successive A-scans is equivalent to the instantaneous Doppler shift represented by the derivative of the phase with respect to time, which can be written as
Here denotes phase operator, is the central wavelength of the light source, is the Doppler frequency shift, is the A-scan rate, and is the Doppler angle between optical axis and motion direction. In most cases, one assumes the sample displacement is only in optical axial direction, meaning that only displacement vector projection in the axial direction is detected; Eq. (3) can be simplified to
Artifacts Induced by Air Interface Motion
Currently, the use of PhS-OCT often assumes that evaluated between adjacent A-scans is caused by the displacement of the scatterer at the location . However, OCT phase signal is in fact only sensitive to the path-length difference between the OPL of reference mirror and the OPL of the scatterer at the location in the interferometer. In terms of measured phase change, it is accumulated over the optical path along the sample beam. Thus, is the true representation of the displacement of the scatterer only when there is no OPL perturbation above the scatterer of interest. This requirement is met when (1) there is no tissue movement above the scatterer and (2) refractive index along the beam propagating path in the sample arm is uniform over time. While these conditions are generally satisfied in the use of phase-resolved technique to evaluate the blood flow velocity, they become complicated when using PhS-OCT to monitor/trace the propagating waves within sample induced by a mechanical stimulus. Mechanical stimulus may induce different mechanical waves propagating within tissue, for example, shear wave, Rayleigh wave (surface acoustic waves), and longitudinal waves. Even more complicated is the case of heterogeneous properties of the sample because, by definition, elastic heterogeneities produce reflected waves. All these mechanical waves would affect the accuracy of PhS-OCT monitoring of the movement of the scatterer at the location due to the dynamic mechanical wave of interest.
In the current study, where we are interested in using PhS-OCT to track the shear wave propagating within tissue, the phase difference detected is induced by dynamic OPL change. Assuming that the tissue is uniform in refractive index, two main factors can cause a change in OPL: (1) the motion of the scatterer due to the shear wave of interest and (2) the tissue surface motion at air–tissue interface that perturbs the OPL along the detection beam. The latter becomes a severe artifact in the dynamic motion detection. Because the phase change is accumulated along the optical path, the detected value is the summation of phase changes due to these two motions. This situation is shown in Fig. 2.
Note that in our case, we did not consider the change of refractive index due to compression wave for two reasons: (1) The change of refractive index caused by compression is very small. Even in an extreme condition of 1% density change under shock compression,29 the refractive index change is only . This magnitude is too small to cancel out the OPL change in our study. (2) For gels and soft tissues, the assumption of incompressibility is widely used for modeling mechanical waves in tissue,30,31 where the wave amplitude is small (micron scale).
Assume that the scatterer of interest is located at the depth of away from the sample surface, and the distance is between sample surface and the fiber coupler where the reference beam and reflected sample beam interfere; the total equivalent OPL is then
The OPL difference between the times and , which gives rise to , is then
Here, the term is caused by the true motion of the localized scatterer. is caused by the surface motion, e.g., SAW induced by the mechanical stimulus. This latter term would give rise to an artifact to the experimental measurements of the motions beneath the sample surface in the domain. Therefore, within the B-frame image of the displacement field that is determined by the PhS-OCT system, a propagating surface wave will lead to stripe artifacts according to the above analysis.
From Eq. (7), when the scatterer is close to surface, i.e., , the scatterer displacement will approach the surface displacement, i.e., . Then, we have
Thus, if the detected phase near the sample surface is , we have
In this study, a PhS-OCT system based on SD-OCT configuration was used to observe the artifacts and validate the proposed phase compensation approach. The light source was a superluminescent diode (SLD) (DenseLight Semiconductors Ltd., Singapore) emitting 1.31 μm central wavelength laser with spectral bandwidth. The spectrometer consisted of a diffraction grating, focusing lens, and InGaAs line scan camera (Sensors Ltd., USA) with 1024 pixels and 14-bit precision. The imaging was accomplished by a galvo-scanner that scanned the probe beam over the sample. The synchronization among galvo-scanner, mechanical wave generation, and data acquisition was controlled precisely by a house-developed LabVIEW program.
A low-voltage, stacked piezoelectric actuator (Thorlabs, New Jersey) was used to generate mechanical waves propagating in the sample. A power amplifier (AE Techron, Indiana) with a controllable voltage gain was used to optimally drive the actuator’s capacitance. A slender stainless steel rod with a small polished tip was fixed to the actuator as the contact to the sample. The actuator assembly was then mounted to a precision positioning stage, enabling precise advancement of the actuator-tip to gently contact to the sample surface (see Fig. 3).
Similar to our previous studies,32 we used tissue-mimicking phantoms to verify our proposed approach. Agar is an easily accessible material that can be used to produce tissue phantoms with controllable mechanical properties similar to those of human soft tissue.33 To produce these phantoms, a proper amount of agar-agar powder (Fisher Scientific Inc., USA) was stirred into boiling distilled water until completely dissolved. A few drops of milk were added as scattering particles to facilitate OCT detection.
Scan Protocol and Reconstruction of Displacement Field
To track shear waves propagating in tissue, the PhS-OCT was operated in M-B mode, wherein a sequence of 256 A-scans (one M-scan) was captured at every spatial location sequentially within the B-scan (total 128 locations) mode while the actuator repeatedly fired the stimulus. Thus, a complete M-B scan consisted of A-scans. The timing of the M-B scan protocol is shown in Fig. 4. The system A-scan rate was set to 46.992 kHz, while M-scan repetition rate was 140 Hz. Each M-B scan acquisition takes . The 128 M-B scans were rearranged to form 256 B-M scans. In this way, the rearranged B-M scan had an equivalent frame rate the same as the A-line rate (). This high frame rate could provide a direct visualization of the propagating mechanical waves within a range of 1 to 7 kHz, commonly with wavelength of 0.1 to 1 mm in soft tissue, which is appropriate for OCT imaging.
To obtain the displacement field in each frame, the commonly used phase-resolve method was adopted. Before the Fourier transformation, the spectral interferograms data were resampled into the k space evenly distributed,34 and then the depth-resolved phase data were extracted from the complex OCT signals.9 Localized motion was determined by evaluating the phase differences between adjacent A-scans within the M-scan [see Eqs. (4) and (10)]. The operation was repeated for all M-scans, and displacement at a given pixel over time was obtained, i.e., , assuming and are lateral and depth positions, respectively.
Modeling of Surface and Shear Wave Pattern
To verify the displacement field in the phantoms, we performed a simulation based on a finite difference model to provide theoretical waveforms for comparison with the artifact-contaminated images. Solutions to wave propagation problems by finite-difference methods (FDM) have been widely used in seismology (e.g., Ref. 35) and biomedical ultrasound (e.g., Ref. 36) for forward modeling the elastic waves, because of their ability to accurately predict the wave propagation in homogeneous media. In this investigation, we adopted an FDM model with the second-order temporal and fourth-order spatial approximations.37 The setup for the simulation model is shown in Fig. 5.Ref. 37.
On the free surface (), we impose the boundary conditions of
Here, is the inward normal of . Model grid size and time step size are determined by the numerically stable requirement of the FDM. The model is set as a square of with the model grid size at and the fixed time step at 10 μs. The simulation was implemented with MATLAB code, from which the snapshots of wave field as a time sequence were obtained.
Presence of Artifact
To demonstrate the presence of the artifact, a layered glass plate phantom is fabricated to create a circumstance that mechanical wave is propagating within a layer of soft phantom between a free surface and a glass plate, while at the bottom there is an air gap above the substrate rigid body (steel plate), see Fig. 6(a). With this configuration, the steel plate can be safely assumed to be static, i.e., without mechanical waves propagating in it, because the glass plate and a layer of air isolate the acoustic waves from the soft phantom. And also the stimulation is not able to launch detectable motion in steel with given frequency and amplitude. In the experiments, the stimulation tip was driven by a stack piezoelectric actuator, and the tip was emitting 5 kHz sine waves of 5 cycles. All the phantoms were in diameter and at least 15 mm deep, to avoid any reflected wave from the boundaries. The vibration amplitude was . The tip was a wedge-like shape, with a width of 1.7 mm edge. The waves generated by this tip can therefore be considered approximately planar, having negligible out-of-plane displacement, which simplifies the experiments for two-dimensional analysis. When attaching the tip to the sample, the OCT system was used to provide the real-time B-scan of the region of interest.
After processing PhS-OCT data, a video that provides visualization of the propagating mechanical wave was generated. The phase map was encoded with bipolar color map, where red color denotes the motion toward the OCT light beam. The phase image was then overlaid onto the OCT B-scan structural image, with higher transparency for lower amplitude of motion, as shown in Fig. 7. Before mapping the phase, a mask was applied to the phase image. The mask was generated from the structural image so that any pixel with low OCT signal intensity, e.g., from transparent structures, was masked out. This masking method was applied to all the phase images in this paper. Four snapshots of the dynamic wave (displacement field) are shown in Figs. 7(a1) to 7(a4); each snapshot is noted by the time offset from . The same applies to all the snapshot images in this paper. The bright line at the bottom of each image is the surface of steel plate. The phase artifacts are clearly indicated in these images, with smaller but opposite phase values to the surface phase, as indicated by arrows in Fig. 7. This supports the analyses presented in Sec. 2.2 that a positive phase on a surface will cause a negative phase (artifact) at all locations underneath the surface. Hence, on the B-scan image, the stripe artifacts are created by the surface motion.
The phase compensation method was then applied on these wave field dataset. The results are shown in Figs. 7(b1) to 7(b4). After compensation, the phase artifacts on the steel plate and also on the glass plate surface are eliminated. Corrections to the wave pattern propagating within tissue phantom are also observed, with the amplitude diminishing to zero from surface to bottom. Note that in this figure, there still exists phase signal (red or blue) within some region that is of low OCT signal. This is largely due to the imperfection of the masking method used. This, however, does not affect the demonstration of artifact compensation.
To show the strength of artifact, the displacement waveform at different depths are plotted in Fig. 8, where the compensated and uncompensated waveform at three depths/locations [marked as “+” in Fig. 7(a1)] are given. At depth #1, the vibration amplitude is the largest because this location is close to the phantom surface where the artifact effect is less obvious. As the depth goes deeper to depth #2, both the phase and amplitude show a severe distortion. Because in the deeper region the true displacement becomes smaller while the artifact amplitude is uniform at all depths, the error rate increases significantly at the deeper region. As shown at the bottom pair of waveforms, there is no displacement signal on the steel plate after compensation as expected, demonstrating the effectiveness of the compensation method. Note that the artifact could lead to the increase or decrease of the apparent amplitude (see depth 2 and steel plate) depending on whether it impacts on the true displacement positively or negatively. In the case of “Phantom depth 2,” the artifact has negative effect; therefore, the compensated curve is higher in magnitude than the measured curve.
Wave Field Simulation
As the simulation was intended to create a numerical solution to the wave field in a homogeous phantom, the dimensions and physical parameters were set to mimick the soft phantom. The longitudinal wave (P-wave) and shear wave (S-wave) velocity were set to 1050 and . The simulation modeled a homogeous phantom with a stimulation source emitting 3 cycles of 5 kHz sine waves on surface. A snapshot of displacement field ( direction component) is presented in Fig. 9(a). In this figure, to avoid any reflections, only a section is shown, which is relatively far from any fixed boundaries. Figure 9(b) is the measured wave propagating pattern at the same time instance of the stimulation. As the stimulation in the experiment was applied with a small angle, the wave front observed in the stimulation presents some difference from that in the simulation. In general, the wavelength, amplitude, and the wave pattern reached good agreements between the experiment and the simulation.
The result of simulation demonstrates the coexisting shear wave and surface wave (Rayleigh wave) modes. Rayleigh wave propagates on the free surface of media, with an active depth of about one wavelength. There is no strict division boundary between the bulk shear wave and Rayleigh wave. The wave field snapshot (Fig. 9) indicates the shear waves with a wavelength of and Rayleigh wave with a similar velocity. The displacement amplitudes are in good agreement with that from prior literature about the amplidude map of shear waves generated by surface point source.38
Shear Waves in Homogenous Phantom
As the wave pattern (spatial-temporal displacement field) is critical for the reconstruction of wave velocity map, the elimination of motion artifact has always been an important issue for all wave imaging modalities. Here, we use a homogeous phantom that is large enough to avoid reflected waves to show the pure shear wave images. As shown in Fig. 6(b), the stimulation source is attached to phantom surface. Figures 10(a1) to 10(a4) show snapshots at four different time instants when the wave propagates within the phantom. The stripe artifacts are clearly noticeable in all the frames. Also, in comparison with the theoretical wavepattern in Fig. 9, the wavefronts exhibit obvious distortion. Figures 10(b1) to 10(b4) show the wave field frames with surface phase compensation. The stripe artifacts and wave distortion are effectively suppressed. The compensated wave pattern still shows a small difference when compared to Fig. 9. This is because the excitations are launched with an angle in the experiments, and also the attenuation is not modeled in the theoretical model. However, the phase images are much more realistic than the artifact-contaminated ones.
There is a circumstance where the artifacts can be mitigated without compensation. This is shown in Fig. 6(c) where the phantom is placed in a container. We then poured a layer of mineral oil on the surface. To avoid the oil surface motion, we covered the oil layer with a thin glass plate. In this case, the oil layer works as an optical coupling media because the oil has a refractive index that is close to that of phantom, and there is no surface motion on the oil due to the glass plate. The measured wave patterens in this oil-covered phantom are presented in Figs. 10(c1) to 10(c4). It can be seen that the images (without compensation applied) are basically artifact-free. However, as the mineral oil still exhibits a slightly higher refractive index than the phantom ( versus ), there exists an opposite stripe, i.e., a very weak amplitude stripe artifact with same phase direction could be observed [as indicated with arrows in Figs. 10(c1) to 10(c4)]. This also confirmed the validity of the compensation method, because for this oil–phantom interface, () is a negative value, opposite to the case of air–tissue interface.
In this study, the phase-resolved motion detection methods in OCT are briefly reviewed, upon which the formation and characteristics of surface ripple artifact are theoretically and experimentally proved. Furthermore, a method of compensating for these artifacts is proposed accordingly. The artifacts are first observed on a phantom with a moving part and a static part, and then on a homogenous phantom with propagating shear waves and surface waves. For both the cases, the compensation method showed the ability to remove the stripe artifacts and correct the distorted waveforms. As the described circumstances are common in dynamic OCE methods, the compensation algorithm is expected to have wide applications. This investigation also suggests that special attention should be given to the phase errors whenever nonhomogeneous optical medium exists, as the varying refraction expresses a lens effect for phase.
In this work, the refractive index is considered constant and uniform in the specimen. Although a complex sample with large variation of refractive index would complicate the proposed compensation method, it is still reasonable to conclude that, using a rough estimation of refractive index, most of the artifact errors could be eliminated, providing much a higher SNR than leaving this artifact uncorrected. Because the artifact amplitude is () times the surface displacement, the variation of refractive index has a small variation range in the same type of biological tissue.
One difficulty for this compensation technique is in the situation where the sample surface is uneven or rough. In this case, the optical path would become more complex, leading to complications for interpreting the results. As there already exist very effective algorithms to inversely compute the OCT beam path for refraction correction,39 it is feasible to use this compensation method based on the beam path corrections. However, it is more computation-intensive. Nevertheless, the method given above is only one possible approach for the artifact correction; more effective and robust approaches remain to be explored.
It is hoped that the theoretical analysis and methods reported here will help improve the accuracy of OCT motion detection techniques and advance the developments of a wide range of related biomedical applications.