A novel spectral analysis technique of OCT images is demonstrated in this paper for classification and scatterer size estimation. It is based on SOCT autoregressive spectral estimation techniques and statistical analysis. Two different statistical analysis methods were applied to OCT images acquired from tissue phantoms, the first method required prior information on the sample for variance analysis of the spectral content. The second method used k-means clustering without prior information for the sample. The results are very encouraging and indicate that the spectral content of OCT signals can be used to estimate scatterer size and to classify dissimilar areas in phantoms and tissues with sensitivity and specificity of more than 90%.
© 2010 Optical Society of America
A great challenge in the fight against cancer is the timely diagnosis of the disease which greatly enhances the chance of survival. Optical imaging can play an important role in the process providing powerful, non-invasive, diagnostic tools. However, in order to reach their full potential, optical techniques must be able to delineate the changes associated with the early stages of cancer (known as dysplasia and carcinoma in situ.) These include sub-micron variations such as disparities in nuclear number and size, sub-nuclear modifications, and variations in other cellular organelles. Currently, this level of detail is only available by histopathology or, non-invasively, by confocal or multi-photon microscopy . Unfortunately, neither of the two modalities has been extensively implemented due their complexity and their limited penetration in tissue. Optical Coherence Tomography (OCT) is a non-invasive tissue imaging technique that generates in vivo cross-sectional images of tissue microstructure. OCT provides a resolution of 1-10 μm with a penetration of a few mm. However, some dysplastic changes are not clearly discernible even at that resolution. One possible remedy for this limitation may be the detection of such changes indirectly, based on the size-dependent spectral variations they induce on the backscattered light.
Measuring spectral modulation to determine the scatterer size and distribution has been the subject of Light Scattering Spectroscopy (LSS) . Cellular nuclei, organelles and other scatterers in epithelial tissues, considered as spherical particles whose interactions with light are described by Mie theory [3,4], produce characteristic spectral modifications. Significant research effort has been focused on measuring the size of scatterers in tissue [5,6], detecting malignancies in vivo  and in situ , and characterizing particle sizes in tissue phantoms [9,10]. Some limitations of LSS stem from the use of a low numerical aperture (NA) and limited penetration due to multiple scattering. Low Coherence Interferometry (LCI) has been employed to resolve the multiple scattering issues .
Size-dependent, spectral information, similar to LSS, can also be extracted from OCT signals using spectroscopic techniques. Spectroscopic OCT (SOCT) can enhance the contrast of the OCT images by differentiating tissue based on scatterer properties rather than just the intensity of the backscattered light. The first demonstration of SOCT used the spectral centroid shift as an indicator of spectral modification to enhance the contrast of OCT images . SOCT gradually became a powerful extension of OCT with a variety of applications [12–15]. It has also been demonstrated that SOCT can be combined with LSS [16,17]. Recently, advances have also been reported in the algorithms and computational methods for differentiating tissue pathologies, based on small-angle approximation of the radiative transfer equation (RTE) , and on texture analysis of the OCT images [19,20]. Additionally, near-real time classification of different types of vegetables based on multivariate analysis of the OCT images was proposed . In a recent paper, AM-FM techniques were employed to classify tissue phantoms based on scatterer size. Although such techniques can achieve high sensitivity and specificity, they do so at the expense of spatial resolution .
In this paper, a novel technique is described for the analysis of the spectral contents of OCT signals which, in combination with advanced statistical methods, can extract information about scatterer size, a diagnostically relevant tissue feature which yet remains below the resolution of OCT. In the imaging volume of an OCT voxel many individual scattering sides, such as organelles, nuclei, nucleoli, etc, can exist simultaneously. These multiple scatterers, which are spaced less than a coherence length apart, interfere coherently in a stochastic manner. The statistical properties of this signal depend on scatterer size and concentration. Using SOCT, this information can be retrieved by analyzing the statistical properties of the backscattered spectra. The spectral content of the backscattered light can be obtained by time-frequency analysis of the interferometric OCT signal. Previous research relied on the Short-Time Fourier transform (STFT) and the continuous wavelet transform (Morlet transform) to obtain the spectrum. This paper introduces a novel technique for calculating and processing the spectral content of OCT images based on autoregressive techniques and Principal Component Analysis (PCA) with their application demonstrated on appropriate phantoms. The particle size can be calculated from a set of linear equations and the samples can be classified into categories using linear discriminant analysis or k-means clustering. The results are very encouraging and indicate that the spectral content of OCT signals can be used to estimate scatterer size and to classify tissue with high accuracy.
2. Materials and methods
2.1. Mie theory predictions
Cellular organelles in epithelial tissue can be modeled as spheroidal scatterers whose interactions with light are described by Mie theory. Based on this approach several methods for scatterer size estimation have been developed such as curve fitting the experimentally measured spectrum to the theoretical predictions or recognizing oscillation patterns from different scatterer spectra. These methods use light sources with a much broader wavelength range as compared to the broadband light sources used in OCT. However, changes can be detected even within a relatively narrow spectral range. Figure 1 shows the different Mie spectra for each scatterer size that correspond to the wavelength range of the OCT light source used during the experiments described below. These calculations were based on the Fortran BHMIE subroutine  which has been modified and converted to Matlab code. The medium refractive index was set to 1.47 (measured experimentally for the gel medium described in the experimental section) and the sphere refractive index was set to 1.59 (provided by the manufacturer.) The spectral differences induced by different scatterer sizes are more than obvious even between 1μm and 2μm size particles. These calculations indicate that such changes are relevant and could be detected in tissues as well since the scattering cross-sections of cells and organelles are in the same range as those of the phantoms (Table 1). The refractive indices of the relevant structures were taken from the literature .
2.2. Spectral estimation
In SOCT, the spectral content of backscattered light can be estimated by time-frequency analysis of the OCT signal. In the literature, the Short-Time Fourier transform (STFT) and the continuous wavelet transform (Morlet transform) have been used for obtaining the backscattering spectrum. In both of these techniques, performance is constrained by the time-frequency uncertainty principle, which states that there exists an inherent trade-off between the spectral resolution and the time resolution . Improvement in one implies degradation in the other. STFT techniques are also particularly prone to noise degradation of the signal.
In general, Time-Frequency Distributions (TFD) can be separated into three categories: (1) Linear TFDs, (2) Cohen’s class TFDs and (3) model-based TFDs. A comparative performance study of different methods from all three categories has been performed but did not result in a recommendation regarding the optimal method for obtaining SOCT spectra . Autoregressive (AR) spectral estimation was used in this study to calculate the spectral density of each section of the OCT image. AR spectral estimation belongs to the model-based TFDs. Burg’s method of AR spectral estimation was used in this paper. This technique has been applied to a variety of applications such as data forecasting, speech coding and recognition, model-based spectral analysis, model-based interpolation, signal restoration. It has also been used for ultrasonic tissue characterization and scatterer size estimation [25,26]. Model-based methods of spectral estimation, such as Burg’s method, converge faster for shorter signals and are computationally more efficient than classical, non-parametric, techniques [27,28]. Analysis of the same OCT images with classical methods (FFT and autocorrelation) was performed for comparison purposes.
The SOCT signals are stochastic process and have finite energy characterized by a power spectral density (PSD). Lets define the autocorrelation function of a stationary stochastic process x(t)
where E is defined as the mean operator. The autocorrelation estimation for the sampled sequence x(n) of the signal is
Based on the Wiener-Khintchine theorem the PSD is
Inserting 2 into Eq. (3) results in
This form of the PSD is defined as the Periodogram and it is a popular spectral estimation technique because its computation is especially convenient. It can be computed by the use of the Discrete Fourier Transform (DFT) using the efficient Fast Fourier Transform (FFT) algorithm.
An AR model for a sampled signal x(n) is given by the difference equation
where w(n) is a zero mean white noise process with variance σp, αk are the AR coefficients, and p is the order od the AR model. The PSD for AR models is
For the AR model the constant |b 0|2 can be replaced by the error variance σ P 2 in the PSD giving
The AR coefficients can be obtained using the Levinson-Durbin recursion  given by
The reflection coefficients, k̂p , are found recursively by minimizing the least-square error between the original signal, x(n), and the AR predicted estimate. Each recursion is based on arbitrary p,k estimates and does not require calculating the autocorrelation function. The PSD is then computed from the frequency response of the prediction filter and has the form:
Successful spectral analysis of complex samples, independent of the employed time-frequency technique, is hindered by a significant amount of spectral modulation due to the random positioning of the scatterers inside the imaging volume. This modulation makes distinguishing the true scattering pattern very difficult. Theoretical studies have showed that spectral analysis is possible given the assumption that the scattering from different locations has the same general wavelength-dependent scattering profile . Biologically this is valid. Epithelial cell architecture demonstrates a relatively uniform distribution of scatterers in axial and transverse directions for micron-scale distances. Experimental observations agree that one SOCT measurement is not enough to distinguish the wavelength-dependent scattering pattern. One simple method to resolve this problem is to average the OCT signal around the scan of interest before obtaining the backscattered spectrum .
2.3. Spectroscopic OCT system and phantoms
For the purposes of this study, a time domain OCT system has been used. The light source was a superluminscent diode (Superlum, Broadlighter) at 1285 nm center wavelength and a 55nm bandwidth providing ~13 μm axial resolution. A fiber optic-based Michelson interferometer was implemented with the reference arm scanned by translating a retro-reflector mounted on a galvanometer arm. The reference arm was scanned slowly, at 20 Hz, to allow for adequate over sampling given the maximum acquisition frequency of the analog-to-digital converter. The sample arm consisted of collimator (OZ Optics) and a 30 mm achromatic lens (Thorlabs) which provided ~18 μm transverse resolution. The signal from the detector was then digitized using a 16-bit 1.25 MHz data acquisition board (National Instruments, PCI-6251).
The system used for these experiments does not incorporate dynamic focusing. The confocal parameter of the system was ~0.4 mm. The depth of focus was relatively large compared to the area under investigation (1-1.5 mm), so no PSF correction was applied. In each imaging session, the focus was set to ~0.4 mm below the top surface of the sample. A more fundamental issue, however, is the existence of depth-dependent spectral differences, as predicted by Thrane et al, which would not be straightforwardly corrected by PSF correction . As a consequence, all the calculations, whether attempting diameter estimation or classification, were performed separately for each individual depth, at steps of 50 μm. When this approach was applied, the estimate uncertainty was not found to vary for the different depths.
The novel analysis technique was first applied on samples consisting of polystyrene microspheres embedded in an acrylamide gel. These samples contained 1, 2, and 4 μm diameter spheres (Polysciences, Polybead) at such concentrations to achieve similar backreflected signal strength for the three gels. Images were acquired at ~12x the Doppler frequency to assure adequate sampling. Following the acquisition of the images, the power spectral density of each section of the data was calculated and multivariate and clustering analysis was performed on the images, as described below.
A small ex vivo study was also performed, in order to examine the application of the method on more complex biological samples such as neuronal tissue. The tissue sample consisted of a rabbit sciatic nerve ligated at the proximal and distal ends and harvested immediately post mortem. The samples were preserved in cold phosphate buffered saline solution and imaged within less than one hour after excision.
2.4. Spectral analysis and classification
The algorithms were first tested on the microsphere samples. For each phantom an image consisting 200 A-scans was acquired. The image size was 1 mm axial × 2 mm transverse. The first 100 A-scans were used as a training set and the remaining 100 as a sample set. After the images were acquired, they were divided into 2400 (axial) × 25 (transverse) pixel neighborhoods using a sliding window. This window size was found to offer a good balance between an increase in sensitivity and specificity and loss of spatial resolution. For the phantom samples, the sliding window overlap was set to 50% axial × 96% transverse (1200 × 24 pixels.) Given the window size and overlap, 1800 spectra were calculated per phantom image which required a total of ~7 minutes per image (using Matlab R2006b, on a 2 GHz personal computer.) The overlap was set to 90% × 96% (2160 × 24 pixels) for the tissue images to assure better spatial resolution.
The average spectrum for each neighborhood was calculated using both Burg’s method (p = 100) and FFT (212 points), for comparison purposes. Each spectrum was normalized to the system spectrum. Although not necessary for the analysis, the normalization accentuates the differences between the spectra which makes studying the spectral variations easier by reducing the intensity variations from the exponential attenuation of the signal. Furthermore, for the Fourier transform based method, median filtering (15 × 15 pixels) was applied to the spectrum in an effort to reduce the effects of the high level of noise present. In the case of the FFT, the autocorrelation of the power spectrum was used, since the full width magnitude around 75% is an indicator of the scatterer size, a feature already exploited by others [31,32].
The statistical analysis of the spectra began with a reduction of the number of variables per sample using PCA. PCA can identify patterns in the data, and express the data in such a way as to highlight their similarities and differences . This is achieved by transforming to a new set of variables, the principal components (PCs), which are mutually-orthogonal linear combinations and are ordered based on the variance of the original variables. PCA was performed on the training data set and then a so-called Feature Vector (FV) was constructed by taking the first n PCs and forming a matrix with these PCs as its columns
It was found that a FV consisting of the first 35 PCs describes more than 99.99% of the training data group variance. To obtain a new orthogonal sample (sa) and training (tr) sets with no redundant information, the FV can be used by left multiplying both group data.
2.4.1. Discriminant based classification
The scatterer-based analysis began with a classification scheme. Using the reduced data, processed as described above, the samples were classified in one of three categories: 1, 2, or 4 μm diameter spheres. First, one-way Multivariate Analysis of Variance (MANOVA) was applied to the training group P̂ xx PCAtr to obtain new linear combinations of variables with maximum separation between categories. For the training phase, a vector dtr was created which contained the diameter of the spheres for each training sample. Size matrix dtr in MANOVA is defined as a categorical variable for grouping the training data. MANOVA resulted in eigenvectors on which the data was projected.
For this classification, a discriminant quadratic type function was used which fits a multivariate normal density to each group, with covariance estimates stratified by group. The classification scheme was applied to both microsphere phantoms and ex vivo tissue samples. The sensitivity and the specificity of the classification were computed for samples of phantoms of spheres with various scatterers diameter.
2.4.2. Scatterer diameter estimate
In addition to the classification scheme, an estimate of the diameter of the microspheres was calculated. A linear dependence was assumed between particle diameter and frequency. The training data group P̂ xx PCAtr was used to find the solution to a set of linear equations
Then the coefficient matrix A can be used to estimate the scatterer size of the sample group
The success of the algorithm was assessed by evaluating the mean error and its standard deviation for all estimates.
2.4.3. k-means clustering based classification
The above methods require some prior knowledge regarding the sample under examination. An unsupervised technique with no a-priori information regarding the spectral content or scatterer size could be a very useful tool. One of the simplest unsupervised algorithms is k-means clustering  which can classify or group a given data set into a certain number of clusters k. The procedure begins by defining k centroids, one for each cluster, the next step is to take each point belonging to a given data set and associate it to the nearest centroid. Then the centroids are recalculated to minimize an objective function, in this case a squared error function
where ∥x i (j) - cj∥ is a distance measure between a data point x i (j) and the cluster centroid cj.
The results of the proposed algorithms are greatly affected by the treatment of the spectra as well as the construction of the feature vector (FV). Different attributes were investigated such as the number of the PCs, the number of the parameters for the AR model, and the number of points for the FFT. In addition, an approximate derivative of the spectrum as well as the square power have been included in the PCA to improve the results by better distinguishing spectral differences. The final choice of parameters was based on optimizing the sensitivity and specificity of the classification results to more than 90% while preserving, as much as possible, the spatial resolution and minimizing the computational effort. The algorithms were optimized for both the Burg and the FFT methods so that the results would be comparable. Figure 2 shows the autocorrelation and average Burg spectra of the three different regions with 1, 2, and 4 μm spheres.
Tables 2 and 3 summarize the results of the analysis. Results from the Fourier based method, as expected were worse due to high level of noise. The mean error for the size estimation was more than 65% and the classification sensitivity was as low ~77%. Burg’s method successfully classified the scatterers by category with sensitivity and specificity more than 90% and the mean error for the size estimation was less than 17%. K-means clustering, even though it assumes no a-priori information regarding the sample, was even more successful in classifying the three regions with sensitivity and specificity of more than 93%.
In order to visualize the results in two dimensions, a Hue-Saturation-Luminance (HSL) scheme was used to overlay the results of the classification on the intensity image. HSL allows the assignment of up to three independent values for every pixel, the cluster number was set for hue (H) and the OCT intensity image was set for luminance (L) with saturation (S) set to 1. Thus, pseudo color images were created in order to display the results more effectively. Figure 3 and 4 illustrates the results.
The clustering technique described above was also applied to the images collected from tissues. In the case of the nerve images, several areas not discernible in the usual intensity OCT image were clearly identified in the cluster image. Those areas probably correspond to fascicles of different neural types, although that claim must be verified by examination of histological cross-sections before any definitive claims can be made.
One issue not addressed in the above experiments is the validity of the technique in the presence of a strong absorber. Such an absorber would modify the spectral content of the backscattered signal, not only of the containing layer but of all inferior layers as well. If the absorber is known, this effect can be estimated and corrected. However, if not corrected, it would result in spectral modifications which would affect the results of the algorithm. Fortunately, most endogenous molecules in tissue exhibit a flat absorption profile, in the 1250-1350 nm range, which does not significantly affect the spectral variations indicative of the scattering properties of tissue .
In this manuscript, a novel technique for the analysis of OCT signals is presented. It is based on autoregressive spectral estimation, classification, and clustering analysis and permits the segmentation of OCT images based on differences in backscattered spectral content. The proposed method is shown to be very effective in differentiating spectral differences which depend on scatterer diameter. Such information can be very useful diagnostically since malignant tissue is known to exhibit characteristic scatterer changes. When fully developed, this technique can be used to differentiate dissimilar areas of the tissue under investigation. One of the major advantages of this method is the fact that it does not require knowledge of the variations in spectral content of the tissue under investigation. Such an unsupervised technique could make available a tool for differentiating tissues with different size scatterers with no a-priori information regarding spectral content or scatterer size. Once verified with histology and pre-clinical studies, such a tool can prove extremely valuable for the investigation of disease tissue features which yet remain below the resolution of OCT.
This research was supported in part by the Research Promotion Foundation of Cyprus.
References and links
2. H. C. V. de Hulst, Light Scattering by Small Particles (Dover Publications, New York, 1981).
3. C. Bohren and D. Huffman, Absorption and Scattering of Light by Small Particles, Wiley Science Paperback Series (John Wiley and Sons, New York, 1998).
5. L. Perelman, V. Backman, M. Wallace, G. Zonios, R. Manoharan, A. Nusrat, S. Shields, M. Seiler, C. Lima, T. Hamano, I. Itzkan, J. Van Dam, J. Crawford, and M. Feld, “Observation of periodic fine structure in reflectance from biological tissue a new technique for measuring nuclear size distribution,” Phys. Rev. Lett. 80(3), 627–630 (1998). [CrossRef]
6. V. Backman, V. Gopal, M. Kalashnikov, K. Badizadegan, R. Gurjar, A. Wax, I. Georgakoudi, M. Mueller, C. Boone, R. Dasari, and M. Feld, “Measuring cellular structure at submicrometer scale with light scattering spectroscopy,” IEEE Sel. Top Quantum Electron. 7(6), 887–893 (2001). [CrossRef]
7. M. B. Wallace, L. T. Perelman, V. Backman, J. M. Crawford, M. Fitzmaurice, M. Seiler, K. Badizadegan, S. J. Shields, I. Itzkan, R. R. Dasari, J. Van Dam, and M. S. Feld, “Endoscopic detection of dysplasia in patients with Barrett’s esophagus using light-scattering spectroscopy,” Gastroenterology 119(3), 677–682 (2000). [CrossRef]
8. M. S. Feld, V. Backman, M. B. Wallace, L. T. Perelman, J. T. Arendt, R. Gurjar, M. G. Müller, Q. Zhang, G. Zonios, E. Kline, T. McGillican, S. Shapshay, T. Valdez, K. Badizadegan, J. M. Crawford, M. Fitzmaurice, S. Kabani, H. S. Levin, M. Seiler, R. R. Dasari, I. Itzkan, and J. Van Dam, “Detection of preinvasive cancer cells,” Nature 406(6791), 35–36 (2000). [CrossRef]
10. A. Wax, C. Yang, V. Backman, M. Kalashnikov, R. R. Dasari, and M. S. Feld, “Determination of particle size by using the angular distribution of backscattered light as measured with low-coherence interferometry,” J. Opt. Soc. Am. A 19(4), 737–744 (2002). [CrossRef]
11. U. Morgner, W. Drexler, F. X. Kärtner, X. D. Li, C. Pitris, E. P. Ippen, and J. G. Fujimoto, “Spectroscopic optical coherence tomography,” Opt. Lett. 25(2), 111–113 (2000). [CrossRef]
12. D. J. Faber, E. G. Mik, M. C. G. Aalders, and T. G. van Leeuwen, “Light absorption of (oxy-)hemoglobin assessed by spectroscopic optical coherence tomography,” Opt. Lett. 28(16), 1436–1438 (2003). [CrossRef]
13. C. Xu, J. Ye, D. L. Marks, and S. A. Boppart, “Near-infrared dyes as contrast-enhancing agents for spectroscopic optical coherence tomography,” Opt. Lett. 29(14), 1647–1649 (2004). [CrossRef]
14. B. Hermann, K. Bizheva, A. Unterhuber, B. Povazay, H. Sattmann, L. Schmetterer, A. Fercher, and W. Drexler, “Precision of extracting absorption profiles from weakly scattering media with spectroscopic time-domain optical coherence tomography,” Opt. Express 12(8), 1677–1688 (2004). [CrossRef]
15. D. Adler, T. Ko, P. Herz, and J. Fujimoto, “Optical coherence tomography contrast enhancement using spectroscopic analysis with spectral autocorrelation,” Opt. Express 12(22), 5487–5501 (2004). [CrossRef]
17. C. Xu, P. S. Carney, W. Tan, and S. A. Boppart, “Light-scattering spectroscopic optical coherence tomography for differentiating cells in 3D cell culture,” Proc. SPIE 6088(608804), 1–10 (2006).
18. I. V. Turchin, E. A. Sergeeva, L. S. Dolin, V. A. Kamensky, N. M. Shakhova, and R. Richards-Kortum, “Novel algorithm of processing optical coherence tomography images for differentiation of biological tissue pathologies,” J. Biomed. Opt. 10(6), 064024 (2005). [CrossRef]
19. X. Qi, M. V. Sivak, G. Isenberg, J. E. Willis, and A. M. Rollins, “Computer-aided diagnosis of dysplasia in Barrett’s esophagus using endoscopic optical coherence tomography,” J. Biomed. Opt. 11(4), 044010 (2006). [CrossRef]
20. C. A. Lingley-Papadopoulos, M. H. Loew, M. J. Manyak, and J. M. Zara, “Computer recognition of cancer in the urinary bladder using optical coherence tomography and texture analysis,” J. Biomed. Opt. 13(2), 024003 (2008). [CrossRef]
21. A. Dunn and R. Richards-Kortum, “Three-dimensional computation of light scattering from cells,” IEEE J. Sel. Top. Quantum Electron. 2(4), 898–905 (1996). [CrossRef]
22. F. Bazant-Hegemark and N. Stone, “Near real-time classification of optical coherence tomography data using principal components fed linear discriminant analysis,” J. Biomed. Opt. 13(3), 034002 (2008). [CrossRef]
23. C. Pitris, A. Kartakoullis, and E. Bousi, “AM-FM techniques in the analysis of optical coherence tomography signals,” J. Biophoton. 2(6–7), 364–369 (2009). [CrossRef]
24. C. Xu, F. Kamalabadi, and S. A. Boppart, “Comparative performance analysis of time-frequency distributions for spectroscopic optical coherence tomography,” Appl. Opt. 44(10), 1813–1822 (2005). [CrossRef]
26. K. Wear, R. Wagner, and B. Garra, “A comparison of autoregressive spectral estimation algorithms and order determination methods in ultrasonic tissue characterization,” Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on 42(4), 709–716 (1995). [CrossRef]
27. J. Proakis and D. Manolakis, Digital Signal Processing, Principles, Algorithms and Applications (Prentice Hall, 2006).
28. P. Stoica and R. L. Moses, Introduction to Spectral Analysis (Prentice Hall, 1997).
29. R. N. Graf and A. Wax, “Temporal coherence and time-frequency distributions in spectroscopic optical coherence tomography,” J. Opt. Soc. Am. A 24(8), 2186–2195 (2007). [CrossRef]
30. I. T. Jolliffe, Principal Component Analysis, 2nd ed. (Springer-Verlag, 2002).
31. J. B. MacQueen, “Some Methods for classification and Analysis of Multivariate Observations,” Proceedings of 5-th Berkeley Symposium on Mathematical Statistics and Probability, University of California Press1, 281–297 (1967).
33. C. L. Tsai, J. C. Chen, and W. J. Wang, “Near-infrared Absorption Property of Biological Soft Tissue Constituents,” J. Med. Bio. Eng. 21(1), 7–14 (2001).