Abstract
Advanced signal reconstruction in polarization-sensitive optical coherence tomography (OCT) frequently relies on an accurate determination of the signal noise floor. However, current methods for evaluating the noise floor are often impractical and subjective. Here we present a method using the degree of polarization uniformity and known speckle intensity statistics to model and estimate the OCT noise floor automatically. We establish the working principle of our method with a series of phantom experiments and demonstrate the robustness of our noise estimation method across different imaging systems and applications in vivo.
© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
Optical coherence tomography (OCT) reveals the microstructure of biological tissue by measuring the path length difference of backscattered light [1]. Polarization-sensitive OCT further determines the polarization state of the detected light [2,3]. Evaluating the depth-dependent evolution of the polarization state permits reconstruction of tissue birefringence. Similarly, the degree of polarization uniformity (DOPU) is a quantitative metric expressing the randomness of the polarization states within a local region around each pixel in an OCT tomogram [4–6]. Variations in DOPU can be attributed to specific optical sample properties and have been found to provide useful information in a variety of applications, such as retinal, wound repair, and coronary plaque imaging [5,7,8]. However, measurement noise also significantly affects the detected polarization states by causing a sample-independent randomization and reduction of DOPU [9]. To avoid this noise-induced loss of DOPU and to improve DOPU interpretation, methods for correcting DOPU have been developed, but they rely on the signal-to-noise ratio and, hence, knowledge of the noise floor [10]. Signal-to-noise dependent corrections have also served to improve birefringence imaging [11] and the reconstruction of the attenuation coefficient [12].
In OCT, the signal-to-noise ratio is defined as the amplitude of the squared norm of the complex-valued tomogram divided by the mean of the same quantity in the absence of any meaningful signal. Ideally, this noise floor is dominated by the reference light shot noise. Determining the noise floor in an OCT system is typically accomplished in one of two ways. First, the sample arm can be blocked while recording data, so only the noise of the reference light and electronic system noise are captured. Alternatively, after data acquisition, tomograms can be manually segmented to find a region of interest (ROI) that “appears devoid” of a signal to estimate the noise floor. Each of these methods has considerable drawbacks. First, blocking the sample arm removes potential noise contributions arising from reflections in the optical train of the sample arm. Additionally, in the case of fiber-based, or enclosed, OCT systems, blocking the sample arm may not be practical. Manual segmentation of an ROI devoid of a signal in an OCT tomogram is subjective, and it is challenging to robustly determine a truly “signal-less” region. Segmentation can also be impractical in large datasets where the sample location may vary, and the noise floor needs to be evaluated over time. A fully automated method for computing noise post-data acquisition without the need for modifying optical pathways or manual segmentation could therefore benefit PS-OCT data analysis.
To address this challenge, we developed a fully automated method that utilizes the link between noise-induced DOPU loss and speckle statistics to estimate noise in an OCT system. This fully automated method can be performed on any OCT dataset provided both intensity and DOPU signals are available. This method does not require any mechanical blocking of optical pathways, separate B-scans of noise, or manual segmentation.
Data analysis was conducted using MATLAB (The MathWorks Inc., Natick, Massachusetts). Unless noted, data was collected using a polarization-sensitive optical frequency domain imaging (OFDI) system described previously in detail by Lippok et al. [13]. The system has a central wavelength of 1310 nm with a 110 nm sweep range and 50 kHz repetition rate. Using polarization diverse detection ($h, \nu$), the signals corresponding to two orthogonal illuminating polarization states (1, 2) were simultaneously acquired through multiplexing with acousto-optic modulation. A MATLAB reconstruction algorithm was then used to build intensity tomograms for each of the four polarization channels ($h 1 , h 2 , \nu 1 , \nu 2$), or sums thereof, and compute maps of DOPU independently for both illuminating polarization states ($h 1 , \nu 1$ and $h 2 , \nu 2$). DOPU maps were computed by local averaging of Stokes vector elements [6]. DOPU values ranged from 0 (completely depolarized) to 1 (completely polarized).
In OCT, the coherent superposition of light backscattered by sample structure results in overlapping interference patterns that create speckle, which follow well-defined statistical properties related to the sample properties. The more general case of the incoherent sum of several speckle intensity patterns can be described by a gamma probability distribution function (PDF) [14]. Because of their large dynamic range, tomograms are displayed in logarithmic scale. After log-transformation, the resulting intensity distributions correspond to the log-gamma distribution
Utilizing these properties, we devised a method to automatically estimate the noise floor of an OCT intensity tomogram. Pixels from the intensity tomogram were ignored for estimating noise if their corresponding DOPU value was above a DOPU threshold (${{\rm DOPU}_t}$). To find the optimal threshold, ${{\rm DOPU}_t}$ was incrementally decreased from 1 (includes all pixels) to ignore more and more correlated pixels. For each ${{\rm DOPU}_t}$, PDFs were fit to the left-hand side of the distribution of the remaining pixels. The left side was used to avoid fitting the very undefined shape of pixel distributions of higher intensity values. The left side was comprised of the lowest pixel values and therefore dominated by characteristic noise. To ensure that enough pixels remain to produce a reliable PDF, as ${{\rm DOPU}_t}$ was lowered, we utilized the Dvoretzky–Kiefer–Wolfowitz inequality for a 1% difference at a 95% confidence interval and ignored ${{\rm DOPU}_t}$ values, which resulted in too few pixels to reliably fit the PDF [15]. Afterwards, the ${{\rm DOPU}_t}$ with the corresponding highest $\alpha$ was selected for each polarization channel, identifying the optimum trade-off between the inclusion of a sufficient number of pixels and rejecting signal-dominated pixels. This balance resulted in a PDF matching closely to the known statistics of pure noise. To avoid the required assumption that each polarization channel must have equal mean intensity when summing multiple channels, an optimal ${{\rm DOPU}_t}$ was found for each channel separately by optimizing its shape parameter $\alpha$ to approach 1. After determining the optimum ${{\rm DOPU}_t}$, the mean of the log-gamma PDF was used to directly determine the noise floor
where $I$ is the mean (noise floor), $\alpha$ is the shape parameter, and $\beta$ is the scale of the log-gamma PDF.In Fig. 2(a), a histogram of pixel intensities from the noise tomogram in Fig. 1 is shown. PDFs normalized to histogram counts compare well to the histograms for all pixels and pixels remaining after ${{\rm DOPU}_t}$ filtering. In Fig. 2(b), though nearly identical, the PDFs of all channels after ${{\rm DOPU}_t}$ filtering compared to the PDF of all pixels is shown. This reference PDF (dashed red) is the mean of the PDFs of all four channels found using all pixels in the tomogram. In this ideal case of pure noise, all PDFs show strong agreement.
To validate our method in a more representative sample, we imaged a 300 µL well containing intralipid. In Fig. 3(a), the upper bright line corresponds to the bottom plastic layer of the 96-well plate. Signal from a mirror placed in the sample arm for phase calibration is seen near the bottom of the tomogram. The red ROI is an area devoid of meaningful signal used for reference. The noise pixels identified by thresholding with the optimal ${{\rm DOPU}_t}$ are displayed in a binary image [Fig. 3(b)]. Identified noise pixels are seen outside of the sample area, plastic layer, and phase calibration line, indicating our method is robust to artifacts and other saturation-level signals. Next, the PDFs of each polarization channel were compared to the average PDF of the red reference ROI [Fig. 3(c)]. PDFs of the corresponding polarization states and noise ROI show strong agreement. Histograms comparing all, noise ROI (red square), and ${{\rm DOPU}_t}$ pixels with associated normalized PDFs are shown in Fig. 3(d). The impact of varying the ${{\rm DOPU}_t}$ on the shape parameter $\alpha$ and the number of pixels remaining is shown in Fig. 3(e). As ${{\rm DOPU}_t}$ increases, more pixels are included in the PDF model, causing the distribution to diverge from the ideal $\alpha = 1$ log-gamma PDF. In the case of ${{\rm DOP}_{t}}=1 $, all pixels pass the DOPU filter, resulting in an arbitrarily shaped total pixel distribution [Fig. 3(d)]. This arbitrary shape is the reason the left-hand side of the distribution should be used for fitting the PDF. Left-side fitting also further aides in filtering out pixels with a potentially meaningful signal. The mean noise [Eq. (2)] from the red ROI and DOPU thresholded pixels was found to be $63 . 0 \pm 0 . 6 \;{\rm dB}$ and $62 . 4 \pm 0 . 2 \;{\rm dB}$, respectively, when averaged over the four polarization channels.
Finally, to demonstrate our method is robust, data from three different tissue types acquired using three different systems was analyzed: first, human cardiovascular data acquired by Villiger et al. [Figs. 4(a) and 4(b)] [8]; second, human retinal data acquired in vivo with depth-multiplexed polarization-sensitive OFDI system described by Lippok et al. [Figs. 4(c) and 4(d)] [16]; and, lastly, swine airway data acquired in vivo using a modified commercial NvisionVLE Imaging System (NinePoint Medical, Bedford, MA) [Figs. 4(e) and 4(f)]. The noise floor–adjusted tomograms in absolute SNR scale are shown in the left column [Figs. 4(a), 4(c), and 4(e)]. The PDFs of all polarization channels after DOPU thresholding were compared to the PDF averaged across the four channels within a noise ROI (red) [Figs. 4(b), 4(d), and 4(f)]. Good agreement was observed. The estimated noise in the red ROIs for each sample in raw dB was $56 . 0 \pm 1 . 1 \;{\rm dB}$, $11 . 0 \pm 0 . 6 \;{\rm dB}$, and $67 . 4 \pm 0 . 4 \;{\rm dB}$, respectively. The noise determined using DOPU thresholding was $57 . 7 \pm 0 . 1 \;{\rm dB}$, $11 . 0 \pm 0 . 3 \;{\rm dB}$, and $68 . 59 \pm 0 . 06 \;{\rm dB}$, respectively. The estimated noise was computed by averaging the noise floor found for each channel. The lower noise floor in the retinal data [Figs. 4(c) and 4(d)] is caused by scaling of the data at recording. Crucially, our DOPU filtering method works independently of such scaling.
One potential limitation of this method may arise when analyzing a highly depolarizing sample. We validated this method in tissue samples, which included regions of sample-induced depolarization, such as lipid-rich plaque and retinal pigment epithelium, without encountering any obvious impairment. Potentially, a sample consisting of highly concentrated gold nanorods could cause a noticeable effect. However, even in the case of gold nanorods, the optimal ${{\rm DOPU}_t}$ thresholds utilized by this method are well below the measured DOPU values of all but the most extremely concentrated nanorod samples [13]. Additionally, our fitting of the left-hand side of the intensity distribution serves to minimize the impact of areas with a high signal and low DOPU. In practical terms, we do not expect tissue-induced depolarization to impact our noise floor evaluation. A further limitation is the assumption of a flat noise floor, which is not necessarily true. A depth-dependent noise model could similarly be fit to the noise-dominated pixels that our method finds throughout the tomogram.
In summary, we have presented an automated method for computing the noise floor in PS-OCT tomograms, utilizing the link between intensity, DOPU, and speckle statistics, across a variety of systems and samples. Unlike common strategies for computing the noise floor, this method does not require modification of any optical pathways (blocking), manual ROI segmentation, or separate noise B-scans. Lastly, given the noise floor is computed one polarization channel at a time, this method should be compatible with data acquired using single input state PS-OCT systems regardless if it is a time-domain, frequency-domain, or spectral-domain system.
Funding
National Institute of Biomedical Imaging and Bioengineering (P41EB015903); National Cancer Institute (5F99CA212217).
Acknowledgment
The authors would like to thank Dr. Norman Lippok, Dr. Nestor Uribe-Patarroyo, Dr. Boy Braaf, and NinePoint Medical, Inc., for providing raw data to test our algorithm.
Disclosures
B.E.B.: NinePoint Medical, Inc. (I, P).
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, Science 254, 1178 (1991). [CrossRef]
2. M. R. Hee, E. A. Swanson, J. G. Fujimoto, and D. Huang, J. Opt. Soc. Am. B 9, 903 (1992). [CrossRef]
3. J. F. de Boer, T. E. Milner, M. J. C. van Gemert, and J. S. Nelson, Opt. Lett. 22, 934 (1997). [CrossRef]
4. S. G. Adie, T. R. Hillman, and D. D. Sampson, Opt. Express 15, 18033 (2007). [CrossRef]
5. E. Götzinger, M. Pircher, W. Geitzenauer, C. Ahlers, B. Baumann, S. Michels, U. Schmidt-Erfurth, and C. K. Hitzenberger, Opt. Express 16, 16410 (2008). [CrossRef]
6. N. Lippok, M. Villiger, and B. E. Bouma, Opt. Lett. 40, 3954 (2015). [CrossRef]
7. W. C. Y. Lo, M. Villiger, A. Golberg, G. F. Broelsch, S. Khan, C. G. Lian, W. G. Austen, M. Yarmush, and B. E. Bouma, J. Invest. Dermatol. 136, 84 (2016). [CrossRef]
8. M. Villiger, K. Otsuka, A. Karanasos, P. Doradla, J. Ren, N. Lippok, M. Shishkov, J. Daemen, R. Diletti, R. J. van Geuns, F. Zijlstra, G. van Soest, P. Libby, E. Regar, S. K. Nadkarni, and B. E. Bouma, JACC Cardiovasc. Imaging 11, 1666 (2018). [CrossRef]
9. S.-W. Lee, J.-H. Kang, J.-Y. Yoo, M.-S. Kang, J.-T. Oh, and B.-M. Kim, J. Biomed. Opt. 13, 054032 (2008). [CrossRef]
10. S. Makita, Y. Hong, M. Miura, and Y. Yasuno, Opt. Lett. 39, 6783 (2014). [CrossRef]
11. D. Kasaragod, S. Makita, Y.-J. Hong, and Y. Yasuno, Biomed. Opt. Express 8, 653 (2017). [CrossRef]
12. K. A. Vermeer, J. Mo, J. J. A. Weda, H. G. Lemij, and J. F. de Boer, Biomed. Opt. Express 5, 322 (2014). [CrossRef]
13. N. Lippok, M. Villiger, A. Albanese, E. F. J. Meijer, K. Chung, T. P. Padera, S. N. Bhatia, and B. E. Bouma, Nat. Photonics 11, 583 (2017). [CrossRef]
14. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (W. H. Freeman, 2006).
15. A. Dvoretzky, J. Kiefer, and J. Wolfowitz, Ann. Math. Stat. 27, 642 (1956). [CrossRef]
16. N. Lippok, B. Braaf, M. Villiger, W. Y. Oh, B. J. Vakoc, and B. E. Bouma, J. Biophoton. 12, e201800156 (2019). [CrossRef]