1. Introduction

Optical coherence tomography (OCT) is a micron-scale, non-contact three-dimensional imaging technique that employs low coherence light to derive the structural information within optically scattering media, such as biological tissues [1]. Since its advent, it has been widely used in ophthalmology, dermatology, dentistry, lung, gynecology, urology, kidney transplantation, oncology, artwork examination, and nondestructive material testing [2]. Its popularity is attributed to its high imaging resolution on a micrometer scale and non-destructive nature.

In the sub-field of the anterior segment (AS-) OCT, various research groups have explored and reported iridocorneal angle imaging [3–6]. These systems include the time-domain AS-OCT [3], spectrum domain AS-OCT [4], and the swept-source AS-OCT system [5] in the research field, as well as some current clinic systems in the industries [6]. All these works were developed for diagnostic purposes but not directly related to the monitoring of the treatment region, for example, in glaucoma surgery. Visualization of the iridocorneal angle details provides ophthalmologists/researchers a direct view of the ocular outflow structures, such as the trabecular meshwork (TM), Schlemm’s canal (SC), collector channels (CCs), and vessels in the corneoscleral limbus region [7]. The TM and the nearby SC are both ring structures and are sometimes hard to find on the OCT images given their deep-seated nature and shadowing effect caused by the overlying blood vessels [8]. CCs [9] are smaller opening structures typically 30 μm in diameter with large variation and occasionally appear connecting the outer wall of the SC and episcleral veins. These structures, affected by some physiological or pathological changes, are directly or indirectly responsible for reduced aqueous outflow in human eyes, leading to elevated eye pressure and glaucoma development. To restore the natural ocular aqueous outflow, one can drain into some part of the TM tissue [7] by different methods, e.g, medical, surgical, and laser treatments [10,11]. We attempted to develop a numerical dispersion compensation to obtain a high-resolution OCT imaging system, providing researchers some insights into details of the treatment areas for glaucoma surgeries. However, the axial resolution of a typical AS-OCT is typically insufficient to resolve the angle details. Since axial resolution relates to the Fourier transform of the spectra, one needs to exploit an ultra-broadband light source. However, the dispersion, which is due to the frequency-dependent refractive index of a dispersive material, increases with the increased optical source bandwidth, causing fringe signal chirping and point spread function (PSF) broadening [2].

Various strategies have been used to achieve an optimal axial resolution for high-resolution OCT applications by dispersion compensation either via hardware or software-based methods. Hardware dispersion compensation includes inserting optical materials [12] or using a grating-based phase control delay line [13], both of which are capable of correcting some degree of dispersions. To further improve the axial resolution, numerical dispersion correction is usually needed. It is low-cost, provides continuous tunability, and can be easily applied to other OCT platforms. Some numerical methods that directly measure the dispersion component, have been reported by Singh et al. [14] and Attendu et al. [15], where two mirror measurements at conjugate positions of the zero-delay line are used. Other numerical methods are to find the dispersion mismatch to be compensated, a phase term ${\textrm{e}^{ - \textrm{i}\Delta \mathrm{\Phi }}}$, multiplied by the dispersed cross-spectral density function resulting in dispersion-removed interferograms. Cense et al. [16] proposed a method to obtain the phase function $\Delta \mathrm{\Phi }$ by polynomial fitting and obtaining the non-linear orders of the unwrapped phase of the coherence function of a well-reflecting reference point. Using mirror measurements on the sample arm at a range of different optical path differences, Uribe-Patarroyo et al. [17] performed the simultaneous k-linearization and dispersion compensation based on the fact that the correct k-mapping function is achievable by obtaining the depth-independent chirp of the dispersive fringe. Wojtkowski et al. [18] proposed an iterative procedure to attain the phase function by optimizing the image sharpness as the objective function, which compensates for second and higher-order dispersions. Furthermore, Yasuno et al. [19] used the information entropy of the image as the optimization metric. Most recently, Hillmann et al. [20] reported that short-time Fourier transform (STFT) with small enough windows can be used to extract signal peak shifts by cross-correlating the sub-bandwidth reconstructed images, after which the linear part of the phase function was approximated by integrating the peak shifts $\Delta z(k )$ over k. Ni et al. [21] further attempted to minimize the information entropy of the spectral centroid image as the sharpness metric in the iterative optimization, where centroids of the depth-resolved spectra obtained by STFTs of fringes were calculated for each pixel corresponding to the original image. Note that some of the above-mentioned software methods are only systematic dispersion compensation (a static method correcting only device dispersion) or both the systematic and sample dispersion compensations (a live method dynamically correcting sample dispersion). In this manuscript, since the sample dispersion was found to be negligible, we are solely targeting the systematic dispersion compensation, which can be pre-computed and calibrated only once for an SD-OCT system.

In this paper, we proposed an alternative numerical method using time-frequency analysis (TFA) technique with iterative optimization for the dispersion compensation in an SD-OCT by using only one single mirror fringe. Specifically, STFT (one of the TFA techniques) was performed on the resampled spectral interferogram of a mirror reflection, which resulted in a 2-D depth-wavenumber plot (spatial-spectral domain). The purpose of the automated optimization process is to minimize the wavenumber-dependent ridge variance in the spatial-spectral domain, including the variation of mirror peak position shifts in z-space after STFT and PSFs broadening, both of which are caused by the dispersion. The proposed method is simple and robust, using only one spectral interferogram, and can reduce the fringe chirping and the PSF broadening. Image resolution improvement has been demonstrated when the dispersion was compensated for measurements of the mirror and detailed biological micro-structures in the iridocorneal angle of human cadaver eyes.

2. Theory

2.1 Dispersion in SD-OCT

The interference fringes, after background subtraction, are given by Eq. (1) [2], where k is the angular wavenumber, I(k) is the spectral interferogram, $\rho $ is the responsivity of the detector, S(k) is the power spectrum of the light source, ${R_R}$ and ${R_{Sn}}$ are reflectivities in reference and sample arms, $2\Delta \it{z_n} \cdot {n_g}$ is the optical path difference at sample depth z_n between the two arms which is scaled by the group refractive index of the material ${n_g}$ at depth z_n, $\mathrm{\Delta \Phi }(\textrm{k} )$ is the wavenumber-dependent dispersion correction term.

The phase function $\Delta \mathrm{\Phi }(k )$ can be expanded as a Taylor series around the central wavenumber ${k_0}:$

Since ${a_1}$ is a positional shift term, typically it is sufficient to compensate the second and third-order dispersion compensations with the phase correction term rewritten in Eq. (3) [18], where ${a_2}\; $ and ${a_3}\; $ are the second and third-order dispersion compensation coefficients, ${k_0}$ is the central angular wavenumber. Once $\Delta \mathrm{\Phi }(k )$ is obtained, the Fourier transform of the dispersive fringe multiplied with ${e^{ - i\Delta \mathrm{\Phi }(\textrm{k} )}}$ produces a depth profile with improved axial resolution.

2.2 Time-frequency analysis (TFA)

Time-frequency analyses (TFAs) are a group of techniques that represent both the time and frequency domains simultaneously, which is particularly useful for detecting non-stationary signals [22]. The most common form is the STFT. Briefly, STFT performs a size-fixed, localized segmental windowed Fourier transform, repeatedly sliding through a signal over time, after which it takes the summation of the series of windowed Fourier transforms. In practice, STFT of a resampled spectral interferogram produces a two-dimensional graphic representation of the energy distribution in the wavenumber and depth directions. While it is simple, window size selection and overlap of windows are of extreme importance when STFT is applied to the spectral interferogram because there is a trade-off between spatial and spectral resolutions in those two directions. A narrower spectral window size results in good spectral resolution but poor spatial resolution, while a longer spectral window size will improve the spatial resolution in z with poor spectral resolution. Since we aimed to favor a better spatial resolution to ensure the accuracy in peak position determination in z-space, we chose a large spectral window length with 1024 pixel points. Details of evaluation in window size and overlaps chosen will be discussed later. Figure 1 illustrates how STFT works and can be potentially used for targeting the dispersion mismatch using a mirror measurement. The dispersion is exhibited by the sub-bandwidth reconstructed PSFs’ shift and broadening, although the main cause of the latter comes from the reduced spectral bandwidth by windowing.

 

Fig. 1. Illustration of how STFT works on a typical resampled spectral interferogram of a mirror positioned at the sample arm. For illustration, (A) a complete fringe was partially captured into 5 overlapping segments by a sliding window centered at k(i), where i = 1,2,…,5. (B) Each segmented fringe was then processed via conventional reconstructions into depth profiles, highlighted with different colors corresponding to the window color in (A). With unmatched dispersion in the system, the windowed bandwidth reconstructed PSFs are shifted and broadened accordingly. (C) The central depths of these sub-bandwidth reconstructed PSFs against different central k values of windows are plotted and linearly fitted, which show a direct display of the dispersion. (D) To correct the dispersion mismatch, a method is needed to align all of these shifted peaks to the correct depth position as well as to sharpen the PSFs, as will be described in detail in Sec. 2.3.

Download Full Size | PDF

2.3 Automatic dispersion compensation algorithm using TFA

To start the dispersion compensation term calculation, a resampled spectral interferogram, which has an equal-spacing distribution of the wavenumbers among P = 2048 pixels, is input into an iterative loop that optimizes a user-defined objective function $Var({{a_2}} )$ as defined in Eq. (4). First, STFT was performed on the spectral interferogram under optimization, which resulted in a 2-D depth-wavenumber plot. The detailed MATLAB document for STFT usage can be found at https://www.mathworks.com/help/signal/ref/stft.html. The ridge was then calculated by finding the maximum depth value at each wavenumber being evaluated. The total number of wavenumbers being evaluated (or the total number of local sliding windows) ${K_{eval}}$ in an STFT process is determined by Eq. (5), where M is the length of the sliding window and L is the overlap length between any two consecutive sliding windows, and the ⌊ ⌋ symbols denote the floor function. According to Eq. (5), ${K_{eval}}$ is proportional to the overlap length L but inversely proportional to the sliding window size M. Therefore, one can have more wavenumbers being evaluated by having much denser sliding windows or narrower windows, by consuming more computing time. On the other hand, narrow windows have broadened PSFs, which could lead to difficulty in ridge detections. The selection of window size and overlaps has an impact on the algorithm's output, and details will be discussed in Sec. 4.1. The k-dependent ridge variance, caused by dispersion, was then minimized at the dispersion compensation coefficients ${a_2}$, where $z_{{a_2}}^i$ is the calculated depth at i^th wavenumber, and ${\bar{z}_{{a_2}}}$ is the average of depths across ${K_{eval}}$ equally spaced wavenumbers.

The automated system outputs the optimized dispersion compensation coefficients $a_2^{OPT}$ and dispersion compensation term $\Delta {\mathrm{\Phi }^{OPT}}(k )$, as illustrated in Fig. 2. $\Delta {\mathrm{\Phi }^{OPT}}(k )$ is a P-pixel vector stored in local and is used for future dispersion compensation during routine OCT image reconstruction. Like the method described in Ref. [18], this procedure can be done for any higher-order dispersion compensations. However, since our system manifests unnoticeable third or higher-order dispersions, we calibrated the dispersion mismatch only at the second order. Note that to obtain the desired dispersion compensation vector using this algorithm, only one A-line resampled spectral interferogram is needed.

 

Fig. 2. Flowchart of our proposed SD-OCT dispersion compensation algorithm. The objective function $\textrm{Var}({{\textrm{a}_2}} )$ denotes ridge variance, where ${\textrm{a}_2}$ is a variable of second-order dispersion compensation coefficient under optimization. STFT: short-time Fourier transform; OPT: optimized.

Download Full Size | PDF

3. Materials and methods

3.1 OCT system

A schematic diagram of the custom-built SD-OCT imaging system is shown in Fig. 3. The broadband laser light source (SuperLum Ltd., Ireland) has a bandwidth of $\mathrm{\Delta \lambda } = 165{\; }nm$ at 3 dB spectrum width centered at 850 nm. The source light passing into a single-mode 50:50 fiber coupler (TW850R5A2; Thorlabs, Newton, NJ) is equally separated into two parts, directed to the reference arm and sample arm, respectively. A pair of galvanometric scanning mirrors (Cambridge Technology Inc., Bedford, MA) was implemented which enables simultaneously horizontal (X-axis) and vertical (Y-axis) scanning of the iridocorneal angle of the eye. To match the optical path lengths in the two arms, a motor-controlled mirror was used in the reference arm. A dispersion compensator (LSM54DC1; Thorlabs, Newton, NJ) was used in the reference arm to compensate for the objective lens which has an effective focal length of 54 mm (LSM54-850; Thorlabs, Newton, NJ). The two lights interfere with each other at the coupler and are then detected by a spectrometer (Cobra-S 800, Wasatch Photonics, NC), which comprises a diffraction grating and a high-speed 12-bit complementary metal-oxide-semiconductor (CMOS) line-scan camera.

 

Fig. 3. Schematic diagram of the SD-OCT system. An optical path length mismatch is corrected by using a movable mirror motor-controlled in the reference arm. The mirror or anterior segment of the cadaver eye placed on an eye holder is in the sample arm, as shown in the close-up for the cadaver eye setup. SLD: superluminescent diode; S: spectrometer; FC: fiber coupler; PC: polarization controller; C: collimator; DC: dispersion compensator; MCM: motor-controlled mirror; M: mirror; YGM: Y-axis galvanometer mirror; XGM: X-axis galvanometer mirror; OL: objective lens.

Download Full Size | PDF

The line-scan camera has P = 2048 pixels and runs at a 20 kHz rate. Processed images contain $1024 \times 1000$ pixels per frame, and display images at a rate of 10 frames per second. The measured average axial resolution across 2-mm depths and theoretical lateral resolution are 2.7 $\mu m$ and 8.2 $\mu m$ in air, respectively. The OCT imaging system has a sensitivity of ∼110 dB and is capable of imaging the eye to a depth of approximately 1.5 $mm$ with ∼4.0 $mW$ incident light power on the sample.

3.2 Sample preparation

The human cadaver eyes for OCT imaging were obtained from San Diego Eye Bank within 24-hour postmortem. The study was in compliance with the Declaration of Helsinki. Since we also use these eyes for experiments of femtosecond laser trabeculoctomy where the laser is targeted on the TM region, unrelated structures are removed (living eyes do not encounter this issue because eyes are vibrant and clear) [23,24]. Specifically, the eyes were carefully dissected to keep the TM intact, while removing other components, such as the posterior segment, iris, lens, ciliary body, and uveal tissues. The prepared eyes were then stored in the Optisol corneal storage medium (Bausch and Lomb, Rochester, NY) and refrigerated at 4 °C. Before imaging, the enucleated eye was mounted on a customized eye holder as illustrated in Fig. 3, perfused with Dulbecco's modified Eagle's medium (DMEM) containing 5 µg/mL amphotericin, and 100 µg/mL streptomycin (MilliporeSigma, MA), and placed in an incubator (Sheldon Manufacturing, Inc., Cornelius, OR) at 37 °C, 5% CO₂, and 90% humidity for 30 min. The eye was then placed on a movable mechanical stage with 5 degrees of freedom, i.e., X, Y, Z translations, rotation, and tilting under the OCT imaging system. The OCT beam was focused on the limbus to image the iridocorneal angle. The cornea surface was kept moist to prevent dehydration during imaging. Note that cadaver tissue swelling will make the relationship between ex vivo and in vivo dimensions uncertain.

3.3 Methods

To test the feasibility of the proposed dispersion compensation algorithm, a 19.6 mm, N-BAK1 glass manufactured block (LSM04DC; Thorlabs, Newton, NJ) was deliberately added into the OCT beam, which introduces dispersion mismatch between the two arms. A mirror was placed at the focus on the sample arm while the delay line was adjusted in the reference arm such that the PSF was located at 200 $\mu m$ from zero-delay. Ten spectral interferograms were taken at this location and put into the algorithm independently to obtain the dispersion compensation vectors. We tested the effectiveness and robustness of the algorithm by shifting PSFs at different locations as well as applying the dispersion compensation term to human cadaver eye data. In addition, to verify of the algorithm works with more dispersion, we added two same glass blocks as did previously described, with the PSF located at 800 $\mu m$ from zero-delay.

In performing STFT, a 1024-pixel-sized Hann window with 99% overlapping was applied to the 2048-pixel sized spectral interferogram totaling 94 sliding spectral windows in this work. The Simplex search method (“fminsearch” function with default stopping criteria in MATLAB R2021b, MathWorks, Inc., Natick, MA) was used to minimize the ridge variance. Our dispersion compensation algorithm as detailed in Sec. 2.3 takes only around 0.5 seconds in a MacBook Air 13-inch, 2017 (1.8 GHz Dual-Core Intel Core i5, 8 GB 1600 MHz DDR3).

4. Results

4.1 Selecting the overlap length and window size for STFT processing

Figure 4(A) and 4(B) show the effects of overlap ratio, which is defined by the ratio of overlap length to the window size, on the dispersion compensated axial resolution and the averaged compute time at an array of window sizes. Figure 4(C) demonstrates that a strong linear relationship exists between the averaged compute time and the total number of windows. To evaluate the influence of window size on the dispersion compensated axial resolution and compute time, we first fix the overlap ratio to 0.9, and the results are displayed in Fig. 4(D). More details can be found in Appendix Table 2.

 

Fig. 4. (A) Axial resolution and (B) averaged compute time $\Delta t$ in milliseconds to generate a dispersion compensation vector against a range of the overlap ratios, i.e., 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, and 0.99. Note that data are not available at some overlap ratios for window size M=1024 (in red). (C) A coefficient of determination of 0.99888 indicates that there is a strong correlation between the total number of windows K_eval and $\Delta t$. (D) Obtained axial resolution (y-axis left) and $\Delta t$ (y-axis right) versus a series of window sizes.

Download Full Size | PDF

Table 1. Statistics of the extracted a₂ values from 10 mirror fringes with induced dispersion block(s)

View Table | View all tables in this article

In Fig. 5, we show a few examples to represent the impact of the window size on the spectral and spatial distribution, a pair of competing factors in an STFT. The k-axis is the wavenumbers being evaluated, which are consist of a series of the central k values of each sliding window. For different combinations of window size and overlaps, the total number of the sliding windows and the interval in k will change, resulting in a difference on the k-axis. For example, a large window has a larger starting k value on the k-axis than a smaller window does, and a larger overlap has a finer k increment on the k-axis. Even so, either a large or small window will cover the spectrum from one end but might miss some information from the other end of the spectrum. As expected, significant loss of spectrum will only occur if large window sizes are used with small overlaps.

 

Fig. 5. Spectral and spatial distribution using STFT with various window sizes (A) M=1024, (B) M=256, (C) M=64, and (D) M=16. Overlap ratios = 0.9. Red is the detected ridge of the distribution. Failure of edge detection in D.

Download Full Size | PDF

4.2 Dispersion compensation using mirror data

We plotted the energy maps before and after dispersion compensation in Fig. 6(A) and 6(B) for comparison. Using the proposed method discussed in Sec. 2.3, we obtained the dispersion compensation term in Fig. 6(C). It showed a quadratic curve with a larger amount in the lower and higher frequency bands relative to the central working frequency, because we consider only a₂ optimization in our study. Higher-order dispersion can also be performed if the dispersion is considerable. To avoid false ridge detections from the strong direct current component, spectral interferograms of the mirror at the focus of the sample arm were collected while shifting the PSF at 200 $\mu m$. Figure 6(D) shows the PSFs comparison before and after dispersion compensation.

 

Fig. 6. Spectral and spatial distributions (A) before and (B) after dispersion compensations. The red lines centered on the energy distributions indicate the detected ridges. (C) is the extracted dispersion using the proposed method, and (D) shows PSFs before and after dispersion compensations.

Download Full Size | PDF

To test the robustness of the proposed method, we repeated 10 independent spectral interferograms on the algorithm as summarized in Table 1.

The dispersion compensation vector was found to be depth-independent and, once obtained, can compensate for any other depths’ measurements. Figure 7 summarizes the sensitivity roll-off performance of the imaging system over a range of 2 mm in air, corresponding to 1.5 mm in eye tissues.

 

Fig. 7. Sensitivity roll-off performance of the calibrated SD-OCT system by the proposed method, in which the reflectivity in dB is calculated by taking the 20·log of the amplitude of Fourier transform of the spectral interferograms.

Download Full Size | PDF

4.3 Comparison to polynomial fitting method using Hilbert transform

The polynomial fitting method is done by performing a Fourier transform of the spectral interferogram, symmetrically gating the signal, inverse Fourier transformation, Hilbert transformation to obtain the phase angle, unwrapping the phase, and finally fitting an n-order polynomial to smooth the curve [16]. Figure 8(A) and 8(B) show the corrected PSFs using the polynomial fitting method at a range of fitting orders and using our proposed method in two different amounts of dispersion introduced into the OCT system. When one dispersion glass block was introduced, the axial resolution using the polynomial fitting method can achieve as sharp as 2.4 $\mu m$. However, the resolution was raised to 3.0 $\mu m$ or larger when 2 dispersion glass blocks were added. Larger sidelobes were observed using the polynomial fitting method when more dispersion existed. The proposed method achieved approximately 2.6–2.7 $\mu m$ in both cases. The extracted dispersion curves between different methods are compared in one dispersion block (Fig. 8(C)) and two dispersion blocks (Fig. 8(D)). Our quadratic, convex dispersion compensation can “linearize” the raw, concave unwrapped phases. The dispersion curves extracted by the polynomial fitting method are also convex, thus capable of correcting some dispersion. In addition, the extracted dispersions using the polynomial fitting method are contingent on the fitting orders, especially when a large amount of dispersion is induced.

 

Fig. 8. PSF comparisons between the proposed method and the polynomial fitting method with various fitting orders in an SD-OCT by introducing (A) one dispersion glass block, and (B) two dispersion glass blocks. Note that the “proposed” PSF is superimposed on each of the “polyfit” PSFs. All of them should locate at the same depth, but they are intentionally shifted to avoid crowd. (C, D) Measured dispersion curves of both techniques (solid orange: polynomial fitting methods; dash orange: proposed method) with the raw unwrapped phase (blue, left y-axis) from the Hilbert transform are plotted for the two cases of dispersion, respectively.

Download Full Size | PDF

4.4 Imaging for iridocorneal angle in human cadaver eyes

The current gold standard for angle imaging associated with glaucoma is gonioscopy [25], which can view the TM surface through the cornea but it’s extremely subjective and operator-dependent. For the potential glaucoma treatment application, we implemented the OCT imaging system specifically for the iridocorneal angle in human eyes. A 2mm×1.5 mm area in human cadaver eyes was shown in Fig. 9. To remove the speckle noise, 10 sequential images at the same location were taken to obtain an averaged image. The dispersion compensated images were sharper and showed more clearly the anatomical details in the iridocorneal angle. For instance, SC, CCs, and adjacent vessels were clearly visualized on the dispersion compensated images while barely visible on uncorrected images. Note that the polynomial fitting method also represents quite similar results to our proposed method therefore not repeated here.

 

Fig. 9. Typical averaged (A, B) horizontal and (C, D) vertical scanning OCT images of the iridocorneal angle of human cadaver eyes with a field of view $2mm \times 1.5mm$, compared between (A, C) without and (B, D) with the proposed dispersion compensation algorithm. The beam probe was arranged approximately perpendicular to the limbus to ensure the best visibility of TM, SC, and surrounding microstructures. SC, CCs, and adjacent vessels were found on the dispersion compensated images while some of them were not visible or blurred on the uncorrected images. Scale bar is 200 $\mu m$. AC: anterior chamber; TM: trabecular meshwork; SC: Schlemm’s canal; CC: collector channel.

Download Full Size | PDF

5. Discussion

Optical dispersion remains an unneglectable issue in the development of an OCT system, especially when using light sources with broad bandwidth, such as the 165 nm spectral bandwidth used in our study. In addition, the dispersion is difficult to detect in the spectral interferogram domain. To address these challenges, we developed a new TFA automated iterative framework that extracted the dispersion compensation by optimizing the energy distribution of a dispersive spectrum from a reflection of a mirror. The mirror data showed that dispersion induced discontinuity or disturbances and k-dependent fluctuation of ridges in the spatial-spectral domain. This variance was minimized via an iterative procedure together with STFT. An intuitive observation is that, if we align the central depths of PSFs to eliminate the k-dependence in the spatial-spectral domain, the FWHM of the resultant 1-D PSF can be greatly reduced after ridge variance minimization (Fig. 1, Fig. 6(A) and 6(B)).

Selecting appropriate parameters (overlaps and window size) in an STFT affects the results of dispersion compensation on both achieved axial resolution and compute time. Although the overlap ratio has little effect on the axial resolution, it does matter when a large window is used. This is because, for example, if an M=1024 window was used, some overlap was needed to ensure a sufficient number of windows K_eval was used. Otherwise, it will lead to the failure of dispersion extraction, such as in the cases of overlap ratios $\le 0.6$ (Fig. 4(A) and 4(B), Appendix Table 2). Non-overlap can be considered if median or smaller sizes of windows (M=512, 256, 128, 64, 32) are used, considering less compute time is required (Fig. 4(B), Appendix Table 2). In addition, window size selection is closely related to the axial resolution and compute time. The window size does not influence the axial resolution for most window sizes; but to some degree when the windows are too small (M=16 and 8), the axial resolution dramatically worsens (Fig. 4(D)). This is possibly due to the use of extremely narrow windows, which means higher spectral resolution in the k-axis, causing the spatial resolution to greatly degrade, leading to the failure of ridge detection (Fig. 5(D)). A smaller window also has a problem with a longer compute time (Fig. 4(B), Appendix Table 2). To summarize, when using an STFT for dispersion compensation in our algorithm, one could either select a larger window size (M=1024) with some overlaps (0.7 or higher in our case), or select a median or median small size window (M=512, 256, 128, 64, 32) with no or small overlaps. The shaded area in Fig. 4(B) incorporates a broad range of optional choices if 1 second of averaged compute time is considered to be acceptable.

Previous studies have attempted to use STFT and/or iterative optimization methods for correcting the dispersion mismatch between the two arms [20,21]. By simply scanning a mirror at the sample arm in this work, we presented an alternative method to compensate for the dispersion in an SD-OCT. Unlike using cross-correlation to find the peak shifts [20], we used a simple maximum search method to detect the locations of ridges, followed by an iterative optimization procedure for ridge variation minimization. In addition, Ni et al.’s work [21] focused on the optimization of the sample centroid image’s information entropy, as a further refinement of the technique demonstrated by Wojtkowski et al. [18] and Yasuno et al. [19] towards optimizing an image’s contrast and entropy, respectively. A simple and widely used polynomial fitting method proposed by Cense et al. used Hilbert transform to extract the phase function $\Delta {\Phi }$. This method worked pretty well when small dispersion was introduced in our experiment (Fig. 8(A)). However, axial resolution degraded as the dispersion mismatch grew (Fig. 8(B)). Our method showed consistent corrected PSFs in either situation (Fig. 8(A) and 8(B), Table 1).

In addition, our algorithm is simple, automatic, and robust for repeated measurements. The dispersion compensation terms obtained from the 10 spectral interferograms of a mirror positioned in the sample arm have only a −0.46% (−0.37%) coefficient of variation (Table 1). Another advantage of our algorithm is that it provides a broad range of parameter selection for an STFT, as discussed before. Lastly, it is convenient to use STFT to directly and dynamically visualize how the spectral non-linearity and dispersion mismatch represents at each step through the entire OCT processing (in Appendix C).

 

Fig. 10. STFT of a mirror’s spectral interferogram with (A) background subtraction, (B) normalization to reference spectrum, (C) resampling, (D) windowing, and (E) the proposed dispersion compensation. The color bar indicates the amplitude of the STFT.

Download Full Size | PDF

The resolution and imaging capability of the angle details were improved with the proposed method using the STFT iterative technique. The dispersion compensation algorithm enabled high-fidelity visualization of CCs adjacent to the SC, which was lost on the uncompensated images. The imaging capability of these micro-structures is comparable to Kagemann et al.’s work [8]. Furthermore, our horizontal and vertical imaging can assist to access angle details more efficiently, which is crucial before and during glaucoma surgeries. For instance, in Fig. 9(B), the horizontal scanning allows ophthalmologists to have a different view of the continuous structures i.e., TM/SC/CCs connected with each other rather than only the cross-sectional images in the vertical scanning in Fig. 9(D). This will be beneficial for the treatment planning of the surgical volumes in our ongoing laser studies in human cadaver eyes [23,24] because laser drilling through TM that is close to the opened SC/CCs might have a positive effect on intraocular pressure reduction. Given the importance of better image visualization and characterization of the ocular outflow structures in the iridocorneal angle for glaucoma surgery, future studies include further image quality improvement of the TM/SC/CC areas, as well as studies of OCT image-guided femtosecond laser trabeculoctomy in cadaver eyes for the treatment of glaucoma [23,24].

6. Conclusion

In conclusion, we presented a new numerical dispersion compensation algorithm for an SD-OCT, for imaging the iridocorneal angle of human cadaver eyes. The dispersion compensation term can be calculated with an automatic iterative procedure that minimizes the k-dependent ridge variance via optimization of energy redistribution of a mirror’s spectral interferogram in the spatial-spectral domain using STFT. Lastly, we demonstrate the feasibility of the proposed method for dispersion compensation in an SD-OCT by evaluating both the mirror and angle imaging of human cadaver eyes.

Appendix

A. Parameters related to STFT processing

Table 2. ^aEffects of overlap ratio and window size on the dispersion compensated axial resolution and compute time.

View Table | View all tables in this article

B. Deriving the measured and theoretical a₂ values of the glass dispersion block

Assuming that each dispersion block includes the measured dispersion of ${a_{2,\; \; meas}}$ in the unit of (${\times} {10^{ - 11}}{m^2})$ and the intrinsic dispersion of ${a_{2,\; inc}}$, though small, exists within the system before the introduction of any dispersion block. From Table 1, we have

(6)$${a_{2,\; \; meas}} + {a_{2,\; inc}} ={-} 4.118\; $$

(7)$$2 \cdot {a_{2,\; \; meas}} + {a_{2,\; inc}} ={-} 8.5078$$

By solving the above equation set, we get

(8)$${a_{2,\; \; meas}} ={-} 4.3898 \times {10^{ - 11}}{m^2}$$

(9)$${a_{2,\; inc}} = 0.2718 \times {10^{ - 11}}{m^2}$$

The ‘+’ and ‘-’ signs indicate that the dispersion tends to reside on the side of the sample or reference arm, respectively. The above calculation assumes that the algorithm is capable of extracting the dispersion perfectly in both, one, or two dispersion blocks. The other extreme possibility is that there is no system intrinsic dispersion but just the algorithm’s limitations in small dispersion detection. The averaged ratio of the second-order dispersion coefficients of the two dispersion blocks to that of one dispersion block is measured to be $\frac{{ - 8.5078}}{{ - 4.118}} = $ 2.07 using the proposed method, comparable to the theoretical value of 2. A possibility could be both the minor system intrinsic dispersion and the algorithm’s imperfectability exist. A relatively low dispersion extraction capability of our algorithm in both one and two dispersion blocks indeed seems possible but it is excluded by the following discussions.

We attempt to analytically derive the theoretical a₂ value for the completeness of this paper. Since there is no theoretical formula for characterizing the wavelength-dependent refractive index, derivation of a perfect theoretical a₂ value is not achievable. Even though, an empirical Sellmeier equation can be used to describe the relationship between the refractive index n and wavelength $\lambda $ for a particular transparent medium, where the Sellmeier coefficients for the N-BAK1 glass are ${C_1} = 1.12365662,\; {C_2} = 0.309276848,\; {C_3} = 0.881511957,\; {B_1} = 0.00644742752\; \mu {m^2},\; \; {B_2} = 0.0222284402\; \mu {m^2},\; {B_3} = 107.297751\; \mu {m^2}$ (https://refractiveindex.info/).

(10)$${n^2} - 1 = \frac{{{C_1} \cdot {\lambda ^2}}}{{{\lambda ^2} - {B_1}}} + \frac{{{C_2} \cdot {\lambda ^2}}}{{{\lambda ^2} - {B_2}}} + \frac{{{C_3} \cdot {\lambda ^2}}}{{{\lambda ^2} - {B_3}}}\; $$

From Eq. (1), we have

(11)$$\mathrm{\Phi } = k \cdot 2d = \frac{w}{c} \cdot 2d \cdot {n_g} + \Delta \mathrm{\Phi }$$

Below also defines the wavenumber k and the group refractive index ${n_g}$ in a dispersive medium:

(12)$$k = \frac{{wn}}{c}$$

(13)$${n_g} = n + w\frac{{\partial n}}{{\partial w}}$$

Therefore, the phase function or dispersion compensation term is calculated by:

(14)$$\Delta \mathrm{\Phi } = k \cdot 2d - \frac{{2dw}}{c} \cdot \left( {n + w\frac{{\partial n}}{{\partial w}}} \right) ={-} \frac{{2d{w^2}}}{c}\frac{{\partial n}}{{\partial w}}\; \; $$

where the speed of light in vacuum $c = 299792458\; \frac{m}{s}$ and the thickness of the glass dispersion block $d = 19.6\; mm$_. Note that the wavenumber k is angular frequency dependent, or the angular frequency w is wavenumber dependent, the refractive index n is wavelength dependent, and the wavenumber and the wavelength are related by $k = \frac{{2\pi }}{\lambda }$_.

Comparing the second-order terms in Eq. (2) and Eq. (3) yields:

(15)$${a_{2,theory}} = \; {\left. { - \frac{1}{2} \cdot \frac{{{\partial^2}\Delta \mathrm{\Phi }(k )}}{{\partial {k^2}}}} \right|_{{k_0}}}$$

Therefore, it is straightforward that ${a_{2,theory}}$ is achievable by calculating the $\frac{{{\partial ^2}\Delta \mathrm{\Phi }(k )}}{{\partial {k^2}}}$ term in theory, considering the definition of the partial derivative as well as the chain rule, the product and/or quotient rule. The full expression of $\frac{{{\partial ^2}\Delta \mathrm{\Phi }(k )}}{{\partial {k^2}}}$ can be eventually evaluated as a function of n and $\lambda $, and these are directly accessible from the Sellmeier equation.

However, it is unpractical to manually derive $\frac{{{\partial ^2}\Delta \mathrm{\Phi }(k )}}{{\partial {k^2}}}$, and even undoable to get its full expression in Mathematica online (Wolfram Mathematica, Wolfram Research, Inc., Champaign, IL) due to its high compute complexity. The Mathematica can indeed compute the value of ${\left. {\frac{{{\partial^2}\Delta \mathrm{\Phi }(k )}}{{\partial {k^2}}}} \right|_{{k_0}}}$ and related Mathematica code is provided in italics below, where ${k_0} = 7581092.7\; {m^{ - 1}}\; $ (central k value after k-linearization during the OCT processing).

c1 = 1.12365662;

c2 = 0.309276848;

c3 = 0.881511957;

b1 = 0.00000000000000644742752;

b2 = 0.0000000000000222284402;

b3 = 0.000000000107297751;

c = 299792458;

d = 0.0196;

l[k_]:=2*Pi/k

n[k_]:=(1 + c1*l[k]^2/(l[k]^2-b1)+c2*l[k]^2/(l[k]^2-b2)+c3*l[k]^2/(l[k]^2-b3))^(1/2)

w[k_]:=c*k/n[k]

Phi[k_]:=−2*d*(w[k])^2*n'[k]/(c*w'[k])

D[Phi[k],{k,2}]/.k->7581092.7

By hitting “shift + enter”, the above code directly gives us an output of $- 3.28396 \times {10^{ - 10}}$. Therefore, the theoretical a₂ value of the dispersion block is ${a_{2,theory}} = 1.64198 \times {10^{ - 10}}\; {m^2}$, which is approximately 4 times compared to our absolute measured ${a_2}$ value. This discrepancy is not attributed to our algorithm’s limitation in extracting the dispersion. As we test, the ${a_{2,theory}}$ value (or other values close to the ${a_2}$ measured value) for the dispersion compensation broadens PSFs and blurs the OCT images compared to that when the ${a_2}$ measured value is applied. We speculate that the assumption using the empirical fitting Sellmeier equation and the compute precision of small numbers in the commercial software might be responsible for this discrepancy.

C. STFT analysis for an entire OCT processing workflow

The spatial-spectral plot provides direct access to the two-dimensional representation of the step-by-step signal processing of a spectral interferogram of a single reflection of a mirror, which illustrates how the spectral non-linearity and mismatch between the two arms changes over the entire image processing. The energy distribution is gradually changing as shown in Fig. 10 A-E. The energy distribution band is downward sloping and locally concentrated at some k numbers after background subtraction (Fig. 10(A)), downward sloping and more uniformly distributed after normalization (Fig. 10(B)). The chirping issue was largely solved by the resampling step, but it still remained observable ridge variances (Fig. 10(C)). Therefore, our algorithm is applied to the resampled spectral interferograms to further extract the underlying dispersion that causes these ridge variances along the k-axis. To avoid unnecessary spectral leakage effects, a Hann window was used to suppress the fringe signals at two ends (Fig. 10(D)). Little ridge variance was observed after the proposed dispersion compensation (Fig. 10(E)).

Funding

National Institutes of Health (5R01 EY030304-02); Research to Prevent Blindness (RPB-203478).

Acknowledgments

We thank Dr. Taiji Chen for insightful discussions for the derivation of the theoretical a₂ value in this work.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical Coherence Tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef]  

2. J. A. Izatt, M. A. Choma, and A.-H. Dhalla, “Theory of Optical Coherence Tomography,” in Optical Coherence Tomography: Technology and Applications, 2nd ed. (Springer, 2015).

3. J. A. Izatt, M. R. Hee, E. A. Swanson, C. P. Lin, D. Huang, J. S. Schuman, C. A. Puliafito, and J. G. Fujimoto, “Micrometer-Scale Resolution Imaging of the Anterior Eye In Vivo With Optical Coherence Tomography,” Arch Ophthalmol 112(12), 1584–1589 (1994). [CrossRef]  

4. P. Li, L. An, G. Lan, M. Johnstone, D. Malchow, and R. K. Wang, “Extended imaging depth to 12 mm for 1050-nm spectral domain optical coherence tomography for imaging the whole anterior segment of the human eye at 120-kHz A-scan rate,” J. Biomed. Opt. 18(1), 016012 (2013). [CrossRef]  

5. B. Potsaid, B. Baumann, D. Huang, S. Barry, A. E. Cable, J. S. Schuman, J. S. Duker, and J. G. Fujimoto, “Ultrahigh speed 1050 nm swept source / Fourier domain OCT retinal and anterior segment imaging at 100,000 to 400,000 axial scans per second,” Opt. Express 18(19), 20029–20048 (2010). [CrossRef]  

6. M. Ayres, R. Smallwood, A. M. V. Brooks, E. Chan, and X. Fagan, “Anterior segment optical coherence tomography angiography,” J. Vis. Commun. Med. 42(4), 153–157 (2019). [CrossRef]  

7. E. R. Tamm, “The trabecular meshwork outflow pathways: Structural and functional aspects,” Exp. Eye Res. 88(4), 648–655 (2009). [CrossRef]  

8. L. Kagemann, G. Wollstein, H. Ishikawa, R. A. Bilonick, P. M. Brennen, L. S. Folio, M. L. Gabriele, and J. S. Schuman, “Identification and Assessment of Schlemm’s Canal by Spectral-Domain Optical Coherence Tomography,” Invest Ophthalmol Vis Sci 51(8), 4054–4059 (2010). [CrossRef]  

9. C. R. Hann, A. J. Vercnocke, M. D. Bentley, S. M. Jorgensen, and M. P. Fautsch, “Anatomic Changes in Schlemm’s Canal and Collector Channels in Normal and Primary Open-Angle Glaucoma Eyes Using Low and High Perfusion Pressures,” Invest Ophthalmol Vis Sci55, 5834–5841 (2014). [CrossRef]  

10. R. N. Weinreb, T. Aung, and F. A. Medeiros, “The pathophysiology and treatment of glaucoma: A review,” JAMA 311(18), 1901–1911 (2014). [CrossRef]  

11. R. Conlon, H. Saheb, and I. I. K. Ahmed, “Glaucoma treatment trends: a review,” Can. J. Ophthalmol. 52(1), 114–124 (2017). [CrossRef]  

12. W. Drexler, U. Morgner, F. X. Kaertner, C. Pitris, S. A. Boppart, X. D. Li, E. P. Ippen, and J. G. Fujimoto, “In vivo ultrahigh-resolution, functional optical coherence tomography,” Opt. Lett. 24(17), 1221–1223 (1999). [CrossRef]  

13. G. J. Tearney, B. E. Bouma, and J. G. Fujimoto, “High-speed phase- and group-delay scanning with a grating-based phase control delay line,” Opt. Lett. 22(23), 1811–1813 (1997). [CrossRef]  

14. K. Singh, G. Sharma, and G. J. Tearney, “Estimation and compensation of dispersion for a high-resolution optical coherence tomography system,” J. Opt. 20(2), 025301 (2018). [CrossRef]  

15. X. Attendu and R. M. Ruis, “Simple and robust calibration procedure for k-linearization and dispersion compensation in optical coherence tomography,” J. Biomed. Opt. 24(5), 056001 (2019). [CrossRef]  

16. B. Cense, N. A. Nassif, T. C. Chen, M. C. Pierce, S.-H. Yun, B. H. Park, B. E. Bouma, G. J. Tearney, and J. F. de Boer, “Ultrahigh-resolution high-speed retinal imaging using spectral-domain optical coherence tomography,” Opt. Express 12(11), 2435–2447 (2004). [CrossRef]  

17. N. Uribe-Patarroyo, S. H. Kassani, M. Villiger, and B. E. Bouma, “Robust wavenumber and dispersion calibration for Fourier-domain optical coherence tomography,” Opt. Express 26(7), 9081–9094 (2018). [CrossRef]  

18. M. Wojtkowski, V. J. Srinivasan, T. H. Ko, J. G. Fujimoto, A. Kowalczyk, and J. S. Duker, “Ultrahigh-resolution, high-speed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12(11), 2404–2422 (2004). [CrossRef]  

19. Y. Yasuno, Y. Hong, S. Makita, M. Yamanari, M. Akiba, M. Miura, and T. Yatagai, “In vivo high-contrast imaging of deep posterior eye by 1-µm swept source optical coherence tomography and scattering optical coherence angiography,” Opt. Express 15(10), 6121–6139 (2007). [CrossRef]  

20. D. Hillmann, T. Bonin, C. Lührs, G. Franke, M. Hagen-Eggert, P. Koch, and G. Hüttmann, “Common approach for compensation of axial motion artifacts in swept-source OCT and dispersion in Fourier-domain OCT,” Opt. Express 20(6), 6761–6776 (2012). [CrossRef]  

21. G. Ni, J. Zhang, L. Liu, X. Wang, X. Du, J. Liu, and Y. Liu, “Detection and compensation of dispersion mismatch for frequency-domain optical coherence tomography based on A-scan’s spectrogram,” Opt. Express 28(13), 19229–19241(2020). [CrossRef]  

22. L. Cohen, TIME-FREQUENCY ANALYSIS, 1st ed. (Prentice Hall, 1995).

23. E. Mikula, G. Holland, H. Srass, C. Suarez, J. V. Jester, and T. Juhasz, “Intraocular Pressure Reduction by Femtosecond Laser Created Trabecular Channels in Perfused Human Anterior Segments,” Transl. Vis. Sci. Technol. 10(9), 22 (2021). [CrossRef]  

24. E. R. Mikula, F. Raksi, I. I. Ahmed, M. Sharma, G. Holland, R. Khazaeinezhad, S. Bradford, J. V Jester, and T. Juhasz, “Femtosecond Laser Trabeculotomy in Perfused Human Cadaver Anterior Segments : A Novel, Noninvasive Approach to Glaucoma Treatment,”Transl. Vis. Sci. Technol.11(3), 28 (2022). [CrossRef]  

25. H. G. Scheie, “Width and pigmentation of the angle of the anterior chamber: a system of grading by gonioscopy,” AMA Arch Ophthalmol. 58(4), 510–512 (1957). [CrossRef]