Ultra-fast line-field low coherence holographic elastography using spatial phase shifting

: Optical coherence elastography (OCE) is an emerging technique for quantifying tissue biomechanical properties. Generally, OCE relies on point-by-point scanning. However, long acquisition times make point-by-point scanning unfeasible for clinical use. Here we demonstrate a noncontact single shot line-field low coherence holography system utilizing an automatic Hilbert transform analysis based on a spatial phase shifting technique. Spatio-temporal maps of elastic wave propagation were acquired with only one air-pulse excitation and used to quantify wave velocity and sample mechanical properties at a line rate of 200 kHz. Results obtained on phantoms were correlated with data from mechanical testing. Finally, the stiffness of porcine cornea at different intraocular pressures was also quantified in situ .


Introduction
Changes in tissue mechanical properties are often associated with disease etiology and can provide valuable information about tissue health [1]. Tissue elasticity is often estimated by analyzing parameters of transversely propagating elastic waves, such as in supersonic shear imaging (SSI) [2] and magnetic resonance elastography (MRE) [3]. However, the low spatial resolution of the parent imaging techniques limits applications of SSI and MRE for small and thin samples, such as the cornea [4]. Brillouin microscopy can measure the biomechanical properties of tissues and cells [5,6]. However, long acquisition times may be a concern for rapid elastographic imaging, and translating the Brillouin shift to quantitative mechanical parameters is still a challenge. Optical coherence tomography (OCT) is a noninvasive imaging technique that provides cross-sectional images with spatial resolution at the micrometer scale [7]. OCT based elastography, optical coherence elastography (OCE), is a rapidly emerging technique that can quantify local tissue biomechanical properties with micrometer-scale spatial resolution [8][9][10]. Generally, excitation methods fall into two categories: static and dynamic. While static loading is powerful for micro-scale mapping of tissue biomechanical properties such as breast cancer biopsy samples [11], there are challenges using such techniques in vivo. Dynamic OCE relies on monitoring and analyzing transient displacements in tissue. For example, Adie et al. have demonstrated differences in mechanical properties of tumor and surrounding normal tissues by analyzing the spectrum of vibrational frequencies induced by dynamic compression [12]. In addition, Kennedy  Although this method has proven useful for research purposes, long acquisition times (over tens of seconds) make this approach impractical for clinical use due to motion artifacts, repeated excitations, and laser safety concerns.
To reduce the acquisition time for elastic wave imaging we have demonstrated an OCE technique at 1.5 million A-lines/s where successive B-scans were acquired with a single excitation, reducing the total acquisition time to milliseconds [24]. However, this ultra-fast OCE technique had an inherent trade-off between the temporal and transverse spatial resolutions. Nahas et al. demonstrated a full-field OCE system to detect acoustically generated elastic waves in an ex vivo rat brain [25], and Chao et al. showed full-field shear wave imaging using speckle contrast imaging method [26]. However, the maximum detectable wave speed in both systems was limited due to the relatively low camera frame rate. Fourier domain line-field OCT is a high speed imaging modality that provides twodimensional cross-sectional maps without mechanical scanning [27,28]. However, frame rates of a few kilohertz restrict imaging of rapid dynamic motion such as an elastic wave in stiff samples.
Recently, we introduced a line-field low coherence holography system (LF-LCH) based on temporal phase shifting using piezo-electric transducer (PZT) [29]. The application was performed on ex vivo chicken breast at a line rate of 16 kHz. However, the line rate was limited due to the need of multiple interferograms in order to retrieve the phase information. Moreover, phase errors were introduced from the PZT hysteresis, nonlinear PZT movement, and vibrations from surroundings during phase shifts [30]. An alternative to PZT modulation, spatial phase shifting, can detect the optical phase of dynamic subjects by introducing a carrier frequency [31]. This technique has the benefit of no moving parts and a rapid line or frame rate by utilizing only one interferogram for phase retrieval [31].
Traditionally, spatial phase shifting has been based on the Fourier transform (FT) [32]. However, the conventional FT method is a global computational process in which noise from different spatial pixels influence one another in spectral domain [33]. Moreover, FT spatial phase shifting requires clean separation of the interferogram spectrum in the Fourier domain [32]. This requirement means that the conventional FT method cannot accurately retrieve phase information from a curved structure that would form a "closed fringe" pattern [34]. In addition, the FT method is a semi-automatic process that requires appropriate window selection in order to extract the fundamental spectrum [35], which may cause a loss of detail at higher spatial frequencies in the interferogram [36]. In contrast, spatial phase shifting using the Hilbert transform (HT) is a high resolution and fully automatic phase retrieval technique that requires only a single interferogram [37] and retains the optical phase information from all spatial frequencies because of the absence of windowed filtering [36]. Furthermore, this method is simple and has less computational cost for phase computation [38]. However, this spatial phase shifting technique has not been explored in optical elastography methods.
In this work, we demonstrate a spatial phase shifting technique using a HT algorithm in combination with a noncontact line-field low coherence holography system to assess the biomechanical properties of the porcine cornea. Using only one air-pulse excitation, the system accurately estimated the elasticity of homogenous and transversely heterogeneous tissue-mimicking phantoms, which were validated with compressive mechanical testing. By using the HT phase retrieval analysis, the line-field system achieved a line rate of 200 kHz with an acquisition time of a few milliseconds. The ultra-fast line rate of 200 kHz is over ten times faster than in our previous work [29] and among the highest speeds for lateral imaging of elastic wave propagation with optical elastography methods. Furthermore, the elasticity changes of an in situ porcine cornea at various artificially controlled intraocular pressures (IOP) were quantified with the presented system, demonstrating the feasibility of using the proposed technique for tissue biomechanical properties assessment.

Tissue mimicking phantoms and biological samples
Tissue mimicking agar phantoms of various concentrations (1%, 1.5%, and 2% w/w, N = 3 for each concentration) were used to validate the proposed technique. Three transversely heterogeneous phantoms, which were composed of 1% and 2% agar concentration components, were created to test the system sensitivity to spatial variations in stiffness. The boundary between the agar concentrations was aligned orthogonal to the line beam. Uniaxial mechanical compressional testing (Model 5943, Instron Corp., MA, USA) was performed on the homogenous phantoms immediately after the OCE measurements. To test the feasibility of the proposed LF-LCH system, the stiffness of in situ corneas from ex vivo porcine eyes (N = 3, Sioux-Preme Packing Corp., IA, USA) was assessed along nasal/temporal axis (in order to avoid corneal mechanical anisotropy ambiguity [39,40]) at IOPs of 10, 15, and 20 mmHg.

Line-field low coherence holography (LF-LCH)
A schematic of the LF-LCH setup is shown in Fig. 1. In brief, the LF-LCH utilized a low coherence light source with a wavelength range of 840 ± 20 nm. A 10 nm band pass filter was used to increase the coherence gate to ~31 µm. The light was expanded to a 4 mm diameter beam using a fiber collimation package (F810APC-842, Thorlabs, Inc., NJ, USA) and directed to a 50:50 non-polarizing beam splitter cube. A cylindrical lens (CL1) was used to create a line focus. The light from the sample and reference mirror was respectively focused by the 50 mm focus lenses, L1 and L2. The illumination power on sample was 165 µW. The interference pattern was focused on the CCD (spL4096-140km, Basler, Inc, German) by another 50 mm focus lens (L3) and created a line focus of ~20 µm in height and 4 mm wide, illuminating 440 CCD pixels. The line rate of the line scan camera was 200 kHz during the phantom experiments, and the CCD counts represent the digitized quantity of the collected electrons, which is proportional to light intensity. For the corneal measurements, the 50 mm focal length lenses (L1 and L2) were replaced by 13.86 mm focal length aspheric lenses (NA 0.18, C560TME-B, Thorlabs. Inc., NJ, USA), to create a line focus of 5.5 µm x 1.1 mm. The pixel resolution was 2.6 µm with the camera line rate of 200 kHz. The air-pulse port had an incident angle of 45 degrees with sample surface due to the physical space of the sample arm optics. Whole porcine eyes were cannulated with two needles for artificial IOP control by our previously published system [41]. One needle was connected to a micro-infusion pump, and the other needle was connected to a pressure transducer. Measurements were made at IOPs of 10, 15, and 20 mmHg.

Reference tilt
The reference arm mirror was slightly tilted to introduce a spatial carrier frequency that shifted the spatial frequency band to fall below the Nyquist frequency, which was determined by the lateral pixel resolution. In this work, the reference mirror was tilted to ensure that all fringes across the line focus were sampled sufficiently as per Shannon-Nyquist theory. In addition, an appropriate tilt angle can control the fringe period such that there are fewer phase reconstruction errors while de-trending using an averaging filter [32]. Once the abovementioned conditions are satisfied, the tilt angle is empirically optimized to obtain the maximum distance of useful fringe formation. However, if the reference mirror tilt angle is too large, the fringe visibility will decrease due to an unbalanced light intensity ratio between the both arms, and the higher frequencies would not allow for sufficient fringe sampling.

Phase extraction by HT-based analytic method
The phase retrieval work flow and example data from a 1% agar phantom are illustrated in Fig. 2. The raw interferogram along time is plotted in Fig. 2(b) and the signal at 4 ms is illustrated in Fig. 2(c). The intensity distribution along spatial domain can be written as: where I 0 , I m , and φ represent bias intensity, fringe modulation, and target phase, respectively, and i is the spatial sample index along the illumination line beam. Here, the target phase, φ, includes the elastic wave displacement, φ e , and the spatial frequency, s, as written: The instantaneous spatial frequency s can be shown as: where D(x i ,t) is the optical path difference between the reference tilt and surface structure. The fringe contrast can be poor due to the intensity bias and low fringe modulation, as shown in Fig. 2(c). In Fig. 2(c), if the optical path difference between both arms is increased by the reference tilt, there would likely be an increase in the fringe spatial frequency. Generally, the intensity bias was not uniform and varied slower than the spatial carrier component, which easily met the criterion of this phase retrieval method [42]. In addition, the intensity bias can be removed by a local averaging filter, which generally includes at least 3 periods of a fringe cycle [37]. The interferogram, U, after de-trending is: (4) According to signal processing theory, the real function u(x) can be associated with an analytic form of ψ(x) = u(x) + iv(x), where the imaginary part is the Hilbert transform (HT) of the real function u(x) [43]: where * is the convolution operator, and the HT is equivalent to filtering that keeps the amplitude of the spectrum the same but shifts the phase of all frequencies by π/2. Thus, the Hilbert transform of the bias-removed signal U can be written as: Due to the reduced fringe visibility from the variation in the intensity bias, the I m term can be obtained by computing the envelope of the analytic signal ψ(x) [44]: Using Eqs. (4) and (7), the contrast enhanced fringe, as shown in Fig. 2(d), can be obtained as: By applying the phase shifting property of the HT, the target phase φ(x i ,t) can be computed by: The contrast enhanced fringe shown in Fig. 2(d) was linearly interpolated to 1024 pixels, and the resulting wrapped phase computed by Eq. (6) is plotted in Fig. 2(e). The phase was unwrapped in time and the structural phase was removed by subtraction of the mean phase of the initial 1 ms of the measurement. The temporal phase profile φ at a certain spatial pixels k was firstly unwrapped and subsequently subtracted with the mean value of the first 1 ms to obtain the elastic displacement φ e(x = k,t) . This was repeated for each spatial position to remove the structural phase. The resulting spatio-temporal displacements were smoothed by a 15 x 15 median filter. The elastic wave propagation on a homogenous phantom can be clearly seen as the inclined blue strip in Fig. 3(a). Selected temporal displacement profiles spaced 0.45 mm apart are plotted in Fig. 3(b). To assess the system performance, we characterized the system by imaging homogeneous agar phantoms. The RMS phase stability was evaluated based on the HT phase retrieval process without the use of the median filter, and the result was 5.5 nm measured from the surface of a 1% agar sample.

Quantification of elastic wave group velocity
The elastic wave group velocities were obtained by linearly fitting the elastic wave propagation delays to the corresponding propagation distances [21]. The velocity, c g , was then translated to Young's modulus, E, using the surface wave equation [45], 3 2 2 where ν = 0.49 was Poisson's ratio [46]. The density, ρ, was assumed as 1000 kg/m 3 for the agar phantoms, and porcine cornea [24].

Phantoms
The system performance was assessed with homogeneous 1% and 2% agar phantoms, and an example spatial-temporal displacement map and corresponding elastic wave temporal displacement profiles from a 1% sample are shown in Fig. 3. The LF-LCH and spatial-phase shifting technique were further evaluated with side-by-side heterogeneous phantoms composed of 2% and 1% agar. The contrast enhanced fringe map from a heterogeneous phantom is shown in Fig. 4(a), and the corresponding spatio-temporal displacement map is illustrated in Fig. 4(b). The change in the wave velocity can be seen by the shift in the slope, as indicated by the dashed black arrows. The stiffness of the 1% homogenous phantoms (N = 3) was 7.48 ± 0.65 kPa as assessed by LF-LCH and 14.6 ± 1.52 kPa as measured by mechanical testing. The elasticity of 2% homogenous agar phantoms (N = 3) was 50.92 ± 1.12 kPa as estimated by LF-LCH and 55.66 ± 8.08 kPa as measured by mechanical testing. The Young's modulus of the 1% and 2% agar in the heterogeneous phantoms (N = 3) as assessed by LF-LCH was 6.96 ± 0.96 kPa and 48.6 ± 7.02 kPa, respectively. The LF-LCH results showed good agreement between both types of phantoms and with mechanical testing as plotted in Fig. 4(c). These results demonstrate that the LF-LCH system can accurately and rapidly perform elasticity assessment at 200 kHz line rate and is sensitive to the stiffness changes along the propagation path of elastic wave.

Elasticity measurement on porcine cornea
After the preliminary phantom experiments, the stiffness of an in situ porcine cornea at artificially controlled IOPs of 10, 15, and 20 mmHg was assessed with the presented technique. Due to the limited illumination power, the line rate was reduced to 100 kHz to increase the exposure time. Figure 5(a) is the spatial-temporal map of the interference pattern from the cornea at 10 mmHg IOP. Figure 5(b) shows the corresponding intensity distribution at 8 ms, and a closed fringe appears from 0.1 to 0.4 mm due to the curvature of the cornea and background tilt. The spatio-temporal phase map, spatio-temporal displacement map, and sample elastic wave displacement profiles are plotted in Figs. 5(c)-5(e), respectively. The propagation of the elastic wave can be seen across the imaging area of 0.56 mm. The change in elasticity of the cornea as a function of IOP is presented in Fig. 5(f). The Young's modulus of the corneas, as quantified by Eq. (7), was estimated as 14.52 ± 1.21, 30.94 ± 2.46, and 55.33 ± 5.08 kPa respectively at 10, 15, 20 mmHg IOP. These results demonstrate the feasibility of the presented technique to quantify the stiffness of soft tissues.

Discussions and conclusions
The performance of the proposed ultra-fast LF-LCH system in combination with a robust Hilbert transform based phase retrieval algorithm was validated with tissue-mimicking agar phantoms. The results in Fig. 4 show small disagreements between the stiffness as estimated by the proposed elastographic technique and mechanical testing, but this discrepancy corroborates with our previous work showing that the surface wave equation underestimates stiffness as compared to mechanical testing [47]. After the preliminary phantom experiments, the ability of the system to evaluate soft tissue stiffness was assessed on an in situ porcine cornea from ex vivo porcine eyes at various IOPs. Due to the thin-plate like structure of the cornea, the detected elastic wave is likely a Rayleigh-Lamb wave [45]. In this case, we have developed a modified Rayleigh-Lamb frequency equation (mRLFE), which accounts for the fluid-structure interface at the sample posterior surface in addition to the finite thickness and wave dispersion, and may be more suitable for the cornea biomechanical assessment rather than the surface wave equation. Nevertheless, the proposed elasticity assessments were used for quick estimation of corneal elastic properties which agree well with our previous results using EWI-OCE [40], demonstrating that the proposed system provides a simple and effective way to detect changes in tissue biomechanical properties.
The results in Fig. 3 present an idealized case in which the interferogram is continuous across the line focus and homogenous in terms of spatial frequency due to the flat surface of sample, which results in accurate phase retrieval. In real biological samples such as the cornea, curvature can generate a "closed fringe" which has multiple fringe spatial frequencies, as shown in Figs. 5(a) and Fig. 5(b). In this case, this HT algorithm may encounter phase ambiguity (π shift), which in Fig. 5(c) would be from 0.26 mm to the end of the line beam. Since the elasticiy assessment relies on elastic wave displacement amplitude rather than the sample topography, the phase ambiguity can be automatically corrected from the reconstructed displacement map by tracking the displacement peak. Additionally, extremely high carrier frequencies can cause sampling errors, and low carrier frequencies can reduce the resolution and cause phase retrieval errors [42,48]. Both of these scenarios occur in the case of the cornea, as seen at the periphery and apex, respectively. To resolve these problems, we reduced the line focus width to increase the fringe spatial sampling. In addition, a 10 nm band-pass filter was used to increase the coherence length to preserve continuous interference across the line beam. However, imperfections from de-trending and insufficient fringe cycles as shown in Fig. 5(a) may still lead to inaccurate phase computation. Nevertheless, results shown in Figs. [3][4][5] show that the retrieved phase is accurate enough for the wave velocity calculation.
In order to provide accurate velocity calculation, the line rate and interference distance are critical. Considering an elastic wave velocity across ~0.5 mm with a line rate of 100 kHz, at least 5 temporally delayed pixels are generally needed for accurate velocity calculation by cross-correlation. As a result, the LF-LCH system can provide a maximum detectable Young's modulus of ~300 kPa, which covers the vast majority of soft tissues and allows for the measurement of small samples. With higher illumination power and longer continuous interference across the line focus, stiffer samples can be assessed.
Due to employment of the band-pass filter, the power on the sample was only 155 µW. Although the incident power is well below the exposure limit for the cornea and dilated pupil [49], the low power requires a longer camera exposure time. Furthermore, the lower power limited in-depth imaging of scattering samples. However, increasing the coherence gate (~31 µm in the present setup) may be acceptable for air-pulse based elastographic applications since the wavelength of the air-pulse induced elastic wave is on the order of millimeters [45]. Thus, incorporating a source with longer coherence and higher power may be a feasible solution for characterizing depth-wise biomechanical properties, as we have shown in the cornea using EWI-OCE [19,20].
However, the imaging field-of-view is limited by the surface curvature and roughness of biological tissues. The curved surface and sample roughness heavily influence the fringe spacing, which may introduce high spatial frequencies that cannot be captured based on Nyquist-Shannon theory. This may be overcome by either enhancing the lateral pixel resolution by adjusting the system magnification at the expense of the imaging field-of-view or by utilizing a camera sensor with higher pixel resolution. In addition, the curved and rough surface may easily exceed the coherence length that directly determines the maximum imaging field-of-view. Thus, the imaging distance of the proposed system could potentially be improved by optimizing the system magnification and increasing the coherence length to a desired range.
To accurately assess the biomechanical properties of the cornea, developing a robust phase velocity algorithm is critical [50]. The bulk motion evident in Fig. 5(d) (after 24 ms) distorts the phase velocity calculation due to the global property of the Fourier transform [51]. Secondly, the de-trending process using the averaging filter with constant window size may not be suitable for closed fringes in the cornea. In Fig. 5(a), the fringe period gradually increased from positions at 0 mm to 0.08 mm, with a maximum period at 0.22 mm. This modulation may cause errors while de-trending the fringe, and an adaptive window size technique may overcome this. Another option to address the bulk motion removal and closed fringe analysis may be a short-time Fourier transform [51]. Thirdly, phase errors caused by the Hilbert transform may appear on the both sides of the imaging distance due to the finite data length, which can be compensated by simulating a look-up table based on the obtained fringe pattern [52]. To provide a more accurate and robust elasticity assessment on the cornea, all of the aforementioned solutions will be included in our future to incorporate the modified Rayleigh-Lamb wave model with this optical system [45,47,50].
Recently, we have demonstrated a line-field LCH system [29] for a single shot imaging of elastic wave propagation in ex vivo chicken breast using four phase shifting method. However, the phase stability was constrained by the hysteresis of the PZT resulted in a phase stability of 12.6 nm. Due to the larger FOV (~4 mm) and poor pixel resolution (~10 μm per pixel), the application was limited to measuring flat samples. In contrast, this work improved the phase stability to 5.5 nm by removing the PZT, and the line rate was drastically increased from 16 kHz in the chicken breast to 100 kHz for the corneal sample and 200 kHz in a phantom by introducing a reference tilt and reducing the line focus beam width. In addition, the superior pixel resolution of the current setup allowed for fringe sampling from the curved corneal surface. Previously, Song et al. [53] and Singh et al.
[24] demonstrated single shot phase-sensitive OCE using Fourier-domain mode-locked lasers with frame rates up to 16 kHz. In contrast, the presented system has a line rate of 200 kHz with the added benefit of being drastically more cost-effective due to the use of a one-dimensional camera and superluminescent diode source. Moreover, the enhanced system line rate may enable additional applications of LF-LCH when integrated with the high frequency external loading methods such as piezo-vibrators for assessing depth-resolved biomechanical properties [16].
In summary, we have successfully demonstrated a novel ultra-fast LF-LCH system with a robust spatial phase shifting algorithm for noncontact dynamic elastography at a line rate up to 200 kHz. The system successfully imaged the propagation of air-pulse induced elastic waves in homogenous and heterogeneous tissue-mimicking agar phantoms and on an in situ porcine cornea from ex vivo porcine eye at various IOPs. Due to the ultra-fast line rate, a single line scan elastographic assessment required only one air-pulse excitation and only one spatio-temporal map of a few megabytes, suggesting the potential for real-time in vivo optical elastography.

Funding
This work was supported, in part, by the U.S. National Institutes of Health (NIH) grants 2R01EY022362, 1R01HL120140, and U54HG006348 and U.S Department of Defense (DOD) Congressionally Directed Medical Research Programs (CDMRP) grant PR150338.