Open Access
Add to BibSonomy
Abstract
Inverse spectroscopic optical coherence tomography (IS-OCT) methods apply inverse problem formulations to acquired spectra to estimate depth-resolved sample properties. In the current study, we modelled the time-frequency-distributions using Lambert-Beer’s law and implemented IS-OCT using backscattering spectra calculated from Mie theory, then demonstrated the algorithm on polystyrene microspheres under idealized conditions. The results are significant because the method generates depth dependent estimates of both the concentration and diameter of scattering particles.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
In medicine and life sciences, the characterization of light scattering, absorption and polarization in biological tissues reveals information on their underlying structure and composition, allowing for diagnosis and monitoring in diverse applications such as arterial plaque in atherosclerosis, the size and nature of tumors, blood perfusion and oxygen saturation levels, and drug and nanoparticle delivery [1]. In particles, at the cellular and subcellular size-scale, light scattering and absorption are wavelength dependent (spectral) properties that can be described using Mie theory [2], a well-established stochastic model that forms the basis for the Monte Carlo model of multiple light-particle interactions [3–6]; alternatively, more recent scattering models also exist that include electromagnetic wave descriptions of light [7,8].
In optical coherence tomography (OCT), low coherence spectroscopy (LCS) and low coherence interferometry (LCI), light from a low coherence source interact with a sample and a Michelson interferometer is used to extract depth-dependent information about the sample from the backscattered light [9]. In spectroscopic OCT (S-OCT) and LCS, broadband illumination and spectrally sensitive signal detection allows for the characterization of spectrally resolved absorption and scattering properties [6] which have been used to quantify media composition from absorption [10–14] and, when combined with Mie scattering models, to quantify particle size [4,5,14–21]. The use of inverse problem approaches to analyze particle size distribution in a sample is prevalent in systems using forward-scattering of light (see for example: [22,23]) but is less common in backscattered light techniques such as OCT where it is known as inverse S-OCT (IS-OCT) [19,20,24–26].
In this paper, we briefly review the theory of S-OCT to show where each of the existing IS-OCT techniques are derived from [19,20,24–26] before extending the technique presented in [26], that used Lambert-Beer’s law to describe the time frequency distributions, to backscattering spectra from Mie theory. Where these earlier techniques have demonstrated capability in estimating particle size [19,20,24–26], the technique presented here provides a quantitative estimate of particle size and concentration. Results are presented from data acquired on solutions of polystyrene microspheres with mean diameters of 0.17, 1.04 and 8.91 µm at concentrations of 0.156, 0.625 and 2.5 volume %. The results demonstrate the capability of IS-OCT to separate experimentally acquired spectra into several Mie scattering components.
1.1 Theory
1.1.1 S-OCT
In OCT, the fields originating from the sample arm $E_S$ and the reference arm $E_R$ interfere with one another to produce the interferogram detected by the photodetector, ${\textrm{I} _{\textrm{D}} }({k,z} )$:
where $k = 2\pi /\lambda $ is the wavenumber, $\lambda $ the wavelength, $2z$ the difference in optical path length between the two arms, and angled brackets indicate an average over the detector acquisition time which is assumed to be much longer than the time of an optical cycle. In a swept source OCT system, a laser with narrow linewidth, $\delta k$, spectrally sweeps across the light source bandwidth while the photodetector acquires the interferogram ${\textrm{I} _{\textrm D} }(k )$, which is then Fourier transformed, ${\mathbb{F}}$, to produce the sample depth dependent OCT signal ${\textrm{I} _{\textrm D} }(z )$[6], Fig. 1(a):
Fig. 1. (a) In OCT, an A-scan is produced by taking the Fourier transform of the entire bandwidth of the interferogram, producing high depth resolution but no spectral resolution. (b) Lambert-Beer’s law can be used to extract backscattering and attenuation coefficients from an A-scan. (c) In S-OCT, a window slides across the interferogram and the STFT produces an A-scan for a narrow spectral band, shown here as different colored lines. The depth-resolved spectral profiles are then generated by extracting depth indexed values from this set of A-scans. (d) The spectrally resolved A-scans can also be modelled using Lambert-Beer’s law, Eq. (4). (e) In IS-OOCT, discrete sets of a priori data from Mie scattering theory [19,20], controlled experiments [21], and literature [26], as well as continuous (non-discrete) parameter functions [24–26] have been used to estimate sample properties through linear inverse problem algorithms.
Assuming that only light from single scattering events are present in the interferogram (first order Born approximation), that light attenuation is due to scattering and absorption (Lambert-Beer’s law), and plane wave illumination from the light source, produces the general form for ${\textrm{I} _{\textrm D} }(z )$ of [27–30], Fig. 1(b):
where $\kappa $ is a system constant which describes the influence of system components such as the photodetector, $\textrm{h} (z )$ is a function which describes the depth-dependent influence of system components such as beam divergence, ${\mu _b}$ the backscattering coefficient, and ${\mu _t} = {\mu _s} + {\mu _a}$ is the attenuation coefficient which is the sum of the scattering and absorption coefficients, ${\mu _s}$ and ${\mu _a}$.Extraction and analysis of the depth-resolved spectral variation from the spectral interferogram is known as spectroscopic OCT (S-OCT), and is performed using a time-frequency-transformation (TFT) method such as a short-time Fourier transform (STFT), double window (DW), or the wavelet transform [6,31]. In the case of STFT, a window, ${\rm w} (k)$, is scanned across the interferogram and Fourier transformed to produce a set of intensity profiles, also called time-frequency distributions (TFD), which are specific to the spectral components encompassed by the window, ${\textrm{I} _{\textrm D} }[{k,z;{\rm w} (k )} ]$. By extracting depth indexed values from this set of TFDs, one constructs depth dependent spectra, ${{\textrm S} _{\textrm{TFD}}}[{k,z;{\rm w} (k )} ]$, Fig. 1(c):
1.1.2 Parameters
Researchers have employed a range of experimental configurations using OCT to extract ${\mu _t}$, ${\mu _a},\,{\mu _s}$, and ${\mu _b}$ from Eq. (3) for tissues and tissue phantoms [12,13,15,29,31], Figs. 1(d) and 1(e). Once extracted, the spectral variation of ${\mu _s}$ and ${\mu _b}$ allows characterization of particle diameter using Mie theory [5,16,32], while the spectral variation of ${\mu _a}$ aids in tissue classification [31]. The parameters ${\mu _a},\,{\mu _s}$, and ${\mu _b}$ are local material properties which, for identical scattering particles, can be calculated from the product of the number density of particles (with units: number of particles per unit volume), $\eta $, the cross-sectional area of the particle, C, and either the backscattering or scattering efficiency, $Q$[2]:
where the subscript $x$ indicates either backscattering or scattering, and $\sigma $ is the effective cross section. ${{\textrm Q} _{\textrm x}}({\phi ,\lambda ,{n_p},{n_m}} )$ are functions of the particle diameter, $\phi $, light wavelength, $\lambda $, and the refractive indices of the particle and the medium, ${n_p}$ and ${n_m}$, and can be found using Mie theory [2,33].While it is not the focus of the current study, other methods have resolved the depth dependent spectra by TFT of the A-scan [31], ${{\textrm I} _{\textrm D} }(z )$, and then estimated the particle size from autocorrelation of the resultant spectra ${{\textrm S} _{\textrm{TFD}}}[{k,z;w (k )} ]$ [4,19,34]. A further, alternative method for characterizing ${\mu _s}$ and ${\mu _b}$ is to describe spatial fluctuations in refractive index in the sample using the refractive index (RI) spatial auto-correlation function. In Refs. [24–26], the authors model this RI function using the Whittle-Matern family of functions that includes a shape factor determining the type of function, and then derive expressions for ${\mu _s}$ and ${\mu _b}$ using the shape factor among other parameters, Fig. 1(e).
1.1.3 IS-OCT to estimate particle size distribution
The linear inverse problem ${\mathbf{Ax}} = {\mathbf d}$ can be used to estimate a sample’s properties, ${\mathbf x}$, from experiment or modelling data, ${\mathbf d}$, via inversion of the coefficient matrix, ${\mathbf A}$, which relates the two. The Moore-Penrose pseudoinverse method for inversion of the coefficient matrix ${\mathbf A}$ is equivalent to the least squares method that minimizes the expression ${||{{\mathbf A}{{\mathbf x}_{\textrm{MP}}} - {\mathbf d}} ||_2}$:
Researchers have applied this linear inverse problem formulation to IS-OCT using several different parameters and inversion algorithms: in Ref. [20], an iterative thresholding algorithm which minimizes the expression ${||{{\mathbf A}{{\mathbf x}_\gamma } - {\mathbf d}} ||_2} - \gamma {||{{{\mathbf x}_\gamma }} ||_1}$, with the L1 norm implemented on the penalty term, was used to estimate diameters of microspheres from spectral interferogram data, Eq. (2), and a coefficient matrix that was populated with Mie spectra for different diameter particles each evaluated across different depths; in Ref. [21], k-means clustering was used to estimate diameters of microspheres from frequencies in the TFT of the depth interferogram ${{\textrm I} _{\textrm D} }(z )$ and a coefficient matrix that was populated using principal component analysis on a set of training data generated in controlled experiments; and, lastly, in Ref. [26], the Nelder-Mead algorithm, which is another iterative method, was used to estimate ${\mu _b}$ and the contributions of ${\mu _a}$ and ${\mu _s}$ to ${\mu _t}$ from the depth dependent spectra, ${{\textrm S} _{\textrm{TFD}}}[{k,z;{\rm w} (k )} ]$, Eq. (4), using the absorption spectra of chlorophyll obtained from literature and the Whittle-Matérn function derivations of scattering and backscattering.
2. Methods
2.1 Linear inverse problem formulation
For the current study, the linear inverse problem was formulated to estimate particle size and concentration from depth dependent spectra, ${{\textrm S} _{\textrm{TFD}}}[{k,z;{\rm w} (k )} ]$, Eq. (4), using a coefficient matrix that was populated with Mie backscattering efficiencies. Using the depth dependent spectra, as in Ref. [26], permitted incoherent averaging of multiple data sets to reduce speckle [37] without loss of signal amplitude – which can occur in the spectral interferogram method [20] if phase instability is present in the spectra, as is common in some wavelength-swept lasers.
In the current study, a set of TFDs were generated, and then each TFD was modelled using Lambert-Beer’s law, Eq. (3), Figs. 1(c) and 1(d). Depth dependent spectra were then constructed from depth indexed values within the set of modelled TFDs, such that:
2.2 OCT system
The OCT system used is polarisation sensitive with separate reference arms for horizontal and vertical polarisation and a quarter waveplate in the sample arm to circularly polarize the incident beam, and has been described in detail in an earlier study [38]. In the current study, this polarisation sensitivity was not utilised because the samples were not birefringent; the total back-reflected intensity was calculated as the sum of horizontal and vertical polarisation intensities. The OCT system used a swept source laser with 50 kHz sweep rate, that sweeps through approximately 1100 wavelengths spaced linearly in k-space (wavenumber space) across the 1259 to 1370 nm bandwidth (Axsun Technologies, AXP50125-6). At each wavelength, the interferograms for each polarisation were measured using balanced photodetectors (Thorlabs, PDB425C) and then digitized at 14 bits per sample and 125 MS/s (National Instruments, NI 5761). A galvo-mirror on the sample arm scanned the incident beam by 0.8 mm laterally across each sample using a 6.25 Hz sawtooth wave while 8,000 A-scans were acquired at the laser sweep rate. This scanning approach produced uncorrelated speckle across the A-scans that was used to reduce speckle during data processing by averaging the resultant TFDs from a total of 80,000 A-scans [37]. The axial resolution was 10 μm in air and signal to noise ratio was 104 dB.
2.3 Microsphere samples
Microsphere samples were prepared from three stock solutions of polystyrene particles in diluted water and 0.02% Sodium Azide, with mean particle diameters of 0.17, 1.04 and 8.91 µm (Spherotec, Inc., SPHERO Polystyrene Particles: PP-015-10, PP-10-10, PP-100-10), Fig. 2. The stock solutions were each vortex stirred and sonicated for 30 minutes, and then a total of nine samples were prepared by diluting each of the three stock solutions (0.17, 1.04 and 8.91 µm mean diameter) with deionised water to concentrations of 0.156, 0.625, and 2.5 volume %. The concentrations were selected to be representative of organelles in biological cells, which for example can vary from 0.01 to 20% in yeast [40], so that the results could aid future development of the method towards biological applications.
Fig. 2. The diameter distribution of the stock solutions of polystyrene microspheres, reproduced from data supplied by the manufacturer.
2.4 Data processing
2.4.1 Time frequency distributions
A set of TFDs (spectrally resolved intensity profiles) were generated for each of the 80,000 interferograms by stepping a 20 nm Gaussian window at 0.1 nm increments across each interferogram and, at each step, zero-padding to 4096 points then applying the Fast Fourier Transform algorithm. The TFDs for each central wavelength were averaged to produce a single set of TFDs with reduced speckle. Depth dependent analysis of the TFDs was realized by extracting depth indexed values from the set of TFDs to construct depth dependent spectra, ${{\textrm S} _{\textrm{TFD}}}[{k,{z_0};w (k )} ]$, Fig. 3(a).
Fig. 3. Variables within Eq. (10). (a) A set of TFDs (spectrally resolved intensity profiles) were generated for each sample, and then depth indexed values from the set of TFDs were used to construct depth dependent spectra. (b) The attenuation coefficient was estimated by fitting a straight line to the natural logarithm of each of the TFDs. (c) The system constant was calculated as the product of the normalized source spectral variation, the photodetector sensitivity (1.0 A/W), and the photodetector gain (2.5 × 105 V/A). (d) The confocal function was calculated for a Gaussian beam in water with waist radius of 46 µm and focal depth of 1.263 mm. (e) and (f) The coefficient matrix was populated with normalized Mie backscattering efficiencies for microsphere diameters: 0.10, 0.14, 0.18, 0.22; 0.9, 1.0, 1.1, 1.2; and, 8, 8.4, 8.8, 9.2, 9.6, 10.0 µm. (g) Backscattering efficiencies with a Euclidean distance less than 1 were removed from the coefficient matrix to improve the matrix condition number to a value comparable to the SNR of the TFDs.
2.4.2 Attenuation coefficient
The attenuation coefficient ${\mu _t}(k)$ was estimated by fitting a straight line to the logarithm of each of the TFDs using the MATLAB function polyval, Fig. 3(b); an approach made possible because the samples were considered homogeneous at the macro-scale. In a heterogeneous sample, the local attenuation coefficient can be estimated by dividing the intensity in the pixel of interest by the sum of the intensities of the pixels located at deeper depths [39].
2.4.3 System constant
The power variation in the laser source spectra was characterised on a mirror covered with a diffusive filter. The TFDs from the mirror, ${{\textrm S} _{\textrm{TFD}}}[{k,{z_{mir}};{\rm w} (k )} ]$, where ${z_{mir}}$ is the depth of the mirror, were then normalised and used to correct for source spectral variation, ${{\textrm H} _{\textrm{mir}}}(k )$, Fig. 3(c). The system constant, $\kappa $, was then calculated as the product of the normalized source spectral variation, the photodetector sensitivity (1.0 A/W), and the photodetector gain (2.5 × 105 V/A).
As a secondary measure of source spectral variation, a calibration function, ${{\textrm H} _{\textrm{cal}}}(k )$, was defined using the desired spectral profile, ${{\textrm S} _{\textrm{MIE}}}[{k,{z_{mir}};{\rm w} (k )} ]$, and the average spectra across the depth range of interest on the 0.17 µm/0.156% microsphere sample, ${{\textrm S} _{\textrm{TFD}}}[{k,{z_{R1:R2}};{\rm w} (k )} ]$, such that:
This calibration function was compared to the source spectral variation but not used to correct data because of the large presence of speckle, Fig. 3(c).
2.4.4 Confocal function
The confocal function, ${\textrm h} (k,z),$ was calculated for a Gaussian beam in water with waist radius of ${r_w}$ = 46 µm and focal depth ${z_f} = $ 1.263 mm, Fig. 3(d), based on measurements of the OCT system. After multiplying data by ${{\textrm H} _{\textrm{mir}}}(k )$ and dividing by ${\textrm h} (k,z)$, the region of interest for analysis was visually selected to extend from the sample surface to the end of the linear region of signal decay on plots of several TFDs with logarithm intensity axis.
2.4.5 Population of coefficient matrix with a priori data
The coefficient matrix ${\mathbf A}$ was populated with backscattering efficiencies that were generated using Mie theory equations implemented in MATLAB [33], Fig. 3(e). The backscattering efficiencies were calculated across the wavelength range of the source, of 1.259 to 1.370 µm (from product test sheet supplied by manufacturer), for microsphere diameters in the range of 0.1 to 10 µm and refractive indices of water and polystyrene microspheres of 1.3225 and 1.59, respectively. The wavelength range of each backscattering efficiency array was then clipped to 1.2771 to 1.3594 µm to account for the window size of the TFT and the data acquisition method that acquired the first 1024 points of the total 1120 points generated by the laser source.
2.4.6 Conditioning of the coefficient matrix
In linear inverse problems, a small matrix condition number indicates good solution stability in the presence of noise in the experiment data. The backscattering efficiencies in the coefficient matrix spanned two orders of magnitude because of the large range of microsphere diameters, which produced a large condition number. Therefore, the backscattering efficiencies for each diameter microsphere were normalized against its maximum value, Fig. 3(f), to reduce the matrix condition number from 8×1015 to 8×1012. Furthermore, the distance matrix, which calculates the Euclidean distance between each column, revealed that several of the backscattering efficiency arrays were similar to one another when observed across the narrow bandwidth of the laser source, particularly across 0.1 to 2 µm range that encompassed the 0.17 and 1.04 µm diameter samples, Fig. 3(g). This minimal variation, which is observed as a low Euclidean distance values, acts to increase the matrix condition number and destabilize the solution. Therefore, the diameter range in the coefficient matrix was reduced to encompass the diameter ranges of the samples: 0.10, 0.14, 0.18, 0.22; 0.9, 1.0, 1.1, 1.2; and, 8, 8.4, 8.8, 9.2, 9.6, 10.0 µm. Columns with a Euclidean distance less than 1 were removed to reduce the matrix condition number to 2×105, which is comparable to the SNR of the TFDs, Fig. 3(e). This targeted and restricted diameter range essentially limits the reconstruction to a categorization problem, where data are fit to different categories, though this may be improved with broader bandwidth sources.
2.4.7 Inverse problem solution
The Lawson-Hanson algorithm iteratively finds a solution to a least squares expression with a non-negative constraint using an active-set method [41] and is implemented in the MATLAB function lsqnonneg. The Lawson-Hanson algorithm was used to fit backscattering spectra from experiment data ${\boldsymbol{\mathrm{\mu}}_{\mathbf b}}$ to the backscattering spectra from Mie theory ${\mathbf{Ax}}$ by minimizing the least squares expression $\mathop {\min }\limits_\textrm{x} {||{{\mathbf{Ax}} - {\boldsymbol{\mathrm{\mu}}_{\mathbf b}}} ||_2}$, Fig. 4(a). The output of the Lawson-Hanson algorithm was an estimate of the particle diameter distribution ${\mathbf x}(z)$, Fig. 4(b). Depth-resolved estimates of the particle diameter distribution were obtained by calculating the backscattering spectra using the appropriate depth-dependent values for each input, where ${\boldsymbol{\mathrm{\mu}}_{\mathbf b}} = {{\textrm S} _{\textrm{TFD}}} \odot {{\textrm h} ^{ - 1}} \odot {[{\kappa \exp ({ - 2{z_0}{\mu_t}} )} ]^{ - 1}}$, Fig. 4(c).
Fig. 4. (a) The Lawson-Hanson algorithm iteratively changes the values in ${\mathbf x}$ to find a solution to the least squares expression $\mathop {\min }\limits_{\mathbf x} \vert\vert{\mathbf{Ax}} - {\boldsymbol{\mathrm{\mu}}_{\mathbf b}}{\vert\vert}_2$ with a non-negative constraint, where ${\boldsymbol{\mathrm{\mu}}_{\mathbf b}}$ is the backscattering spectra calculated from the experiment data, and ${\mathbf{Ax}}$ is the backscattering spectra from Mie theory. (b) The result is an estimate of the diameter distribution, ${\mathbf x}$, with the units of number density (number of particles per unit volume), which have been converted here to volume % for easier comprehension. (c) Depth-resolved diameter distribution estimates were realized by inputting depth-resolved backscattering spectra.
3. Results
The diameter distribution estimates from the surface to a depth of 5 mm are presented in Fig. 5(a) for each of the nine microsphere sample preparations as described in section 2.3. The diameter distribution is expressed as volume concentration (vol %) of each of microsphere diameters included in the coefficient matrix of backscattering Mie spectra: 0.18, 0.9, 1.0, 1.1, 8.0, 8.4, 8.8, 9.2, 9.6, and 10 µm. To aid quantitative analysis, the same data are also presented in box plots in Fig. 5(b) to provide the median, 25th and 75th percentiles of the diameter distribution estimates across depth range. In the box plots, data from each of the diameters have been grouped into three ranges (0.18 µm, 0.9 to 1.1 µm, and 8 to 10 µm) representing the three microsphere samples (mean diameters of 0.17 µm, 1.04 µm, and 8.9 µm).
Fig. 5. Diameter distribution estimates for the nine microsphere sample preparations. (a) Color maps of the diameter distribution estimates with sample depth indicate good consistency in the 0.17 and 1.04 µm samples and poor consistency in the 8.49 µm samples across the 0 to 5 mm depth range. (b) Box plots showing the median, 25th and 75th percentiles of diameter distribution estimates across all depths and a green line indicating the expected value for each sample show the concentrations were underestimated in all samples. In the box plots, data have been grouped into diameter ranges corresponding to the three microsphere samples.
In the three 0.17 µm sample preparations, of 0.156, 0.625 and 2.5 vol.%, the estimated diameter distributions contained a significant concentration of 0.18 µm microspheres across the 0 to 5 mm depth range but consistently estimated the concentrations to be lower than the actual sample preparation by a factor of approximately 4, Fig. 5(a) (left column). A significant concentration of 1.1 µm microspheres and low concentration of large, 8–10 µm, microspheres was also common across the three 0.17 µm sample preparations, Fig. 5(b) (left column). This presence of larger microspheres in the diameter distribution is a possible cause of the underestimated concentrations of 0.18 µm microspheres. In the three 1.04 µm sample preparations, of 0.156, 0.625 and 2.5 vol.%, the estimated diameter distributions contained a significant concentration of 0.9, 1.0, and 1.1 µm microspheres across the 0 to 5 mm depth range, Fig. 5(a) (center column).
The median concentrations in the 0.15 vol.%, 0.625 vol.% and 2.5 vol.% samples were, respectively, 0.03 vol.%, 0.3 vol.%, and 2 vol.%, Fig. 5(b) (center column), indicating relatively poor accuracy albeit an improvement on the 0.17 µm samples. In the three 8.49 µm sample preparations, of 0.156, 0.625 and 2.5 vol.%, the estimated diameter distributions contained significant concentrations of smaller microspheres and relatively low consistency across the 0 to 5 mm depth range depths when compared to the 0.17 and 1.04 µm samples, Fig. 5(a) (right column).
4. Discussion
Mie scattering in OCT [15] and LCS [16] has previously been used to estimate particle size and concentration in samples without depth dependence, and in S-OCT to estimate the size of individual particles with depth dependence [5,19,20]. The current study has demonstrated the potential application of IS-OCT to quantify the depth dependent diameter and concentration of scattering particles using backscattering spectra from Mie theory. The results indicated reasonable accuracy in particle size estimation, good consistency of estimates with depth, and poor concentration accuracy for small particles (0.17 and 1.04 µm), but poor performance in these three areas for large particles (8.49 µm). In the nine sample preparations, the microsphere sizes and concentrations made it highly likely that multiple scatterers were present within each imaging volume (pixel). Incoherent averaging and data filtering approaches were essential to remove artefacts from system instability and speckle generated by multiple scatterers in the volume. Artefacts in the spectra are interpreted as scattering particles and are thought to be a significant contributor to errors in the results.
The limited variation in the Mie backscattering spectra for different particle diameters when observed across a narrow bandwidth hindered the inverse solution stability and was a major limitation of the method in the current study. Shorter wavelength and broader bandwidth light sources can reduce this instability by capturing more of the oscillations in the Mie backscattering spectra, Fig. 6(a), and increasing the Euclidean distance, Fig. 6(b). Such improvements may be achieved using visible wavelength and dual-band (visible and near infrared) OCT systems [42,43].
Fig. 6. (a) Backscattering spectra for polystyrene microspheres diameters 0.1 to 10 µm analyzed across a wavelength range of 0.4 to 1.4 µm and (b) the corresponding distance matrix showing that shorter wavelength and broader bandwidth light sources can improve solution stability by increasing the Euclidean distance and, in doing so, decreasing the coefficient matrix condition number.
With some prior knowledge of the sample, another method to improve solution stability is to remove some of the spectra from the coefficient matrix, as was performed in the current study for the microsphere experiments. Prior knowledge of the sample can also be built into the algorithm in much the same way the ${{\textrm H} _{\textrm{mir}}}(k )$ was implemented; for example, a reference spectral profile might be generated from modelling or measurements across a population, and then spectral deviations from this reference used to indicate changes in size or concentration of scatterers in the sample, where high frequency deviations would indicate a change in large scatterers and low frequency changes a change in small scatterers. In Ref. [26], prior information was used to great effect to split the attenuation coefficient into absorption from chlorophyll and scattering from the sample structure. The frequency of oscillations in the Mie backscattering spectra can also be used to populate the coefficient matrix and is the basis for the interferogram IS-OCT implementation used in Ref. [20]. There, the authors accurately estimated the diameters of 1.0 and 4.5 µm polystyrene particles in air attached to glass slides though with deteriorating accuracy as the spheres occupied the same imaging volume – highlighting the problem of speckle in quantitative measurements as observed in the current study. The IS-OCT studies that fitted Whittle-Matern functions to the TFDs [24–26] did not describe speckle as a problem; possibly because the Whittle-Matern functions characterise the backscattering coefficient using a power law of the form ${k^{4 - D}}$, where $D$ is the shape factor and is greater than zero, which limits the number of oscillations in the fitted curve. While this limitation is not a problem for small diameter microspheres [24], and biological tissues in mammal [24,25] and coral [26], it may not extend well to large diameter scatterers or broader bandwidth OCT systems.
A significant limitation to the discrete particle model approach presented in the current study as well as others [19–21] is that the properties of the scatterers that can be output by the model are inherently limited to the size, shape and refractive index properties of the a priori data used to generate spectra for the coefficient matrix. Any discrepancy between the spectra of the sample and the a priori data are output from the discrete particle model as either a residual error or an erroneous estimate of the sample properties. With this limitation in mind, the discrete particle model is most suited to applications with good prior knowledge of the range of properties of scatterer within the sample, such as quantifying number of virions within a fluid. In contrast, continuous random model approaches such as those that fit Whittle-Matern functions to the sample spectra [24–26] do not suffer limitations from discretized a priori data; an attribute that makes these approaches more suited to biological tissue and other samples with highly varying scatterer sizes, shapes and refractive indices.
The IS-OCT problem becomes non-linear when ${\mu _t}$ cannot be extracted from the data or if non-linear effects such as multiple scattering events are included. Multiple scattering increases with the particle size and concentration. In these scenarios, non-linear formulation and inversion techniques are required, for which there are several techniques that are commonly applied to forward scattering problems [23].
Limitations of the current study were the limited number of diameters used in the coefficient matrix because of the long central wavelength and the narrow bandwidth laser source; the use of homogeneous microsphere samples meant that the depth dependent capabilities of the method were not demonstrated; and, multiple scattering events were not considered when extracting the attenuation coefficient. To improve the method in future, models of multiple scattering [15] might be used to improve estimation of the attenuation coefficient, and non-linear inverse problem methods [23] used to find a solution which extracts the backscattering, scattering and attenuation coefficients from the TFDs. In addition, samples with higher particle concentrations and comprising mixtures of different particle diameters would better represent the light scattering properties of biological tissues.
The approach to particle estimation presented in the current study, Eqs. (8), (9), and (10) can also be applied to the attenuation coefficient to estimate the composition of the media in the imaging volume. This application is described briefly here to highlight a potential basis of future studies. Assuming the contribution to the attenuation of light is the coherent sum of individual attenuation events, the spectral modification by sample attenuation can be expressed as the sum of the individual attenuation coefficients, which may be separated into contributions from scattering by particles and absorption by media. While Mie theory can be used to generate particle diameter dependent spectral absorption profiles, the diameter dependence is minor and so this is omitted in the current formulation and instead the imaging volume as a set of volume concentrations of different media using the relationship: ${\mu _a}(k) = 2.303\tau \alpha {}_a(k)$, where $\tau $ is molar concentration, and $\alpha {}_a(k)$ is the molar extinction coefficient. Analysing for a set of $L$ different scattering particle diameters and J different media absorption profiles, yields:
5. Conclusion
A S-OCT method was presented that estimates depth-dependent particle size and concentration from TFDs and Mie backscattering spectra and is most suited to occurrences where multiple scatterers in the imaging volume prohibit IS-OCT analysis of individual particles. The method was demonstrated on polystyrene microspheres in solution using a coefficient matrix containing a limited set of backscattering spectra to compensate for system errors and narrow source bandwidth.
Funding
Marsden Fund; Royal Society Te Apārangi (UoA1509).
Disclosures
The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.
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. A. N. Bashkatov, K. V. Berezin, K. N. Dvoretskiy, M. L. Chernavina, E. A. Genina, V. D. Genin, V. I. Kochubey, E. N. Lazareva, A. B. Pravdin, M. E. Shvachkina, P. A. Timoshina, D. K. Tuchina, D. D. Yakovlev, D. A. Yakovlev, I. U. Yanina, O. S. Zhernovaya, and V. V. Tuchin, “Measurement of tissue optical properties in the context of tissue optical clearing,” J. Biomed. Opt. 23(09), 1 (2018). [CrossRef]
2. C.F. Bohren and D.R. Huffman, “Absorption and scattering of light by small particles,” 2008: John Wiley & Sons.
3. I. Krasnikov, A. Seteikin, and B. Roth, “Advances in the simulation of light–tissue interactions in biomedical engineering,” Biomedical Engineering Letters 9(3), 327–337 (2019). [CrossRef]
4. F. E. Robles and A. Wax, “Measuring morphological features using light-scattering spectroscopy and Fourier-domain low-coherence interferometry,” Opt. Lett. 35(3), 360–362 (2010). [CrossRef]
5. K.C. Zhou, R. Qian, S. Farsiu, and J. A. Izatt, “Spectroscopic optical coherence refraction tomography,” Opt. Lett. 45(7), 2091–2094 (2020). [CrossRef]
6. H. S. Nam and H. Yoo, “Spectroscopic optical coherence tomography: A review of concepts and biomedical applications,” Appl. Spectrosc. Rev. 53(2-4), 91–111 (2018). [CrossRef]
7. T. Brenner, P. R. T. Munro, B. Krüger, and A. Kienle, “Two-dimensional simulation of optical coherence tomography images,” Sci. Rep. 9(1), 12189 (2019). [CrossRef]
8. P. R. Munro, “Three-dimensional full wave model of image formation in optical coherence tomography,” Opt. Express 24(23), 27016–27031 (2016). [CrossRef]
9. N. Bosschaart, T. van Leeuwen, M.C.G. Aalders, B. Hermann, W. Drexler, and D.J. Faber, “Spectroscopic Low Coherence Interferometry, in Optical Coherence Tomography: Technology and Applications,” W. Drexler and J.G. Fujimoto, eds., 2015, Springer International Publishing: Cham. p. 1163–1187.
10. M. Tanaka, M. Hirano, K. Murashima, H. Obi, R. Yamaguchi, and T. Hasegawa, “1.7-μm spectroscopic spectral-domain optical coherence tomography for imaging lipid distribution within blood vessel,” Opt. Express 23(5), 6645–6655 (2015). [CrossRef]
11. L. Yu, J. Kang, C. Jinata, X. Wang, X. Wei, K. T. Chan, N. P. Lee, and K. K. Y. Wong, “Tri-band spectroscopic optical coherence tomography based on optical parametric amplification for lipid and vessel visualization,” J. Biomed. Opt. 20(12), 126006 (2015). [CrossRef]
12. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Differential absorption imaging with optical coherence tomography,” J. Opt. Soc. Am. A 15(9), 2288–2296 (1998). [CrossRef]
13. N. Bosschaart, M. C. G. Aalders, D. J. Faber, J. J. A. Weda, M. J. C. van Gemert, and T. G. van Leeuwen, “Quantitative measurements of absorption spectra in scattering media by low-coherence spectroscopy,” Opt. Lett. 34(23), 3746–3748 (2009). [CrossRef]
14. J. Yi, J. Gong, and X. Li, “Analyzing absorption and scattering spectra of micro-scale structures with spectroscopic optical coherence tomography,” Opt. Express 17(15), 13157–13167 (2009). [CrossRef]
15. M. Almasian, N. Bosschaart, T. G. van Leeuwen, and D. J. Faber, “Validation of quantitative attenuation and backscattering coefficient measurements by optical coherence tomography in the concentration-dependent and multiple scattering regime,” J. Biomed. Opt. 20(12), 121314 (2015). [CrossRef]
16. N. Bosschaart, D. J. Faber, T. G. van Leeuwen, and M. C. G. Aalders, “Measurements of wavelength dependent scattering and backscattering coefficients by low-coherence spectroscopy,” J. Biomed. Opt. 16(3), 030503 (2011). [CrossRef]
17. D. C. 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]
18. V. Jaedicke, J. Agcaer, F. E. Robles, M. Steinert, D. Jones, S. Goebel, N. C. Gerhardt, H. Welp, and M. R. Hofmann, “Comparison of different metrics for analysis and visualization in spectroscopic optical coherence tomography,” Biomed. Opt. Express 4(12), 2945–2961 (2013). [CrossRef]
19. M. Kassinopoulos, E. Bousi, I. Zouvani, and C. Pitris, “Correlation of the derivative as a robust estimator of scatterer size in optical coherence tomography (OCT),” Biomed. Opt. Express 8(3), 1598–1606 (2017). [CrossRef]
20. H. L. Seck, Y. Zhang, and Y. C. Soh, “Contrast enhancement using Mie spectra representation for spectroscopic optical coherence tomography,” Opt. Commun. 338, 328–335 (2015). [CrossRef]
21. A. Kartakoullis, E. Bousi, and C. Pitris, “Scatterer size-based analysis of optical coherence tomography images using spectral estimation techniques,” Opt. Express 18(9), 9181–9191 (2010). [CrossRef]
22. C. Lin, J. Shen, and T. e. Wang, “Use of non-negative constraint in Tikhonov regularization for particle sizing based on forward light scattering,” J. Mod. Opt. 63(16), 1630–1637 (2016). [CrossRef]
23. H. Levinson, “A novel iterative algorithm for solving nonlinear inverse scattering problems,” 2016. Publicly Accessible Penn Dissertations, 1843.
24. J. Yi and V. Backman, “Imaging a full set of optical scattering properties of biological tissue by inverse spectroscopic optical coherence tomography,” Opt. Lett. 37(21), 4443–4445 (2012). [CrossRef]
25. J. Yi, A. J. Radosevich, Y. Stypula-Cyrus, N. N. Mutyal, S. M. Azarin, E. Horcher, M. J. Goldberg, L. K. Bianchi, S. Bajaj, H. K. Roy, and V. Backman, “Spatially resolved optical and ultrastructural properties of colorectal and pancreatic field carcinogenesis observed by inverse spectroscopic optical coherence tomography,” J. Biomed. Opt. 19(3), 036013 (2014). [CrossRef]
26. G. L. C. Spicer, A. Eid, D. Wangpraseurt, T. D. Swain, J. A. Winkelmann, J. Yi, J. Kühl, L. A. Marcelino, and V. Backman, “Measuring light scattering and absorption in corals with Inverse Spectroscopic Optical Coherence Tomography (ISOCT): a new tool for non-invasive monitoring,” Sci. Rep. 9(1), 14148 (2019). [CrossRef]
27. 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]
28. P. Gong, M. Almasian, G. van Soest, D. M. de Bruin, T. G. van Leeuwen, D. D. Sampson, and D. J. Faber, “Parametric imaging of attenuation by optical coherence tomography: review of models, methods, and clinical translation,” J. Biomed. Opt. 25(04), 1–34 (2020). [CrossRef]
29. J. Xi, Y. Chen, and X. Li, “Characterizing optical properties of nano contrast agents by using cross-referencing OCT imaging,” Biomed. Opt. Express 4(6), 842–851 (2013). [CrossRef]
30. A. Guan, S. Richardson, and S. Hinckley, “Optical coherence tomography modeling incorporating scattering, absorption, and multiple reflections,” J. Opt. Soc. Am. A 37(3), 391–398 (2020). [CrossRef]
31. N. Bosschaart, T. G. van Leeuwen, M. C. G. Aalders, and D. J. Faber, “Quantitative comparison of analysis methods for spectroscopic optical coherence tomography,” Biomed. Opt. Express 4(11), 2570–2584 (2013). [CrossRef]
32. C. Xu, P. S. Carney, and S. A. Boppart, “Wavelength-dependent scattering in spectroscopic optical coherence tomography,” Opt. Express 13(14), 5450–5462 (2005). [CrossRef]
33. C. Mätzler, “MATLAB functions for Mie scattering and absorption, version 2,” IAP Res. Rep 8, 1–24 (2002).
34. A. Wax, C. Yang, and J. A. Izatt, “Fourier-domain low-coherence interferometry for light-scattering spectroscopy,” Opt. Lett. 28(14), 1230–1232 (2003). [CrossRef]
35. R.C. Aster, B. Borchers, and C.H. Thurber, “Chapter Three - Rank Deficiency and Ill-Conditioning, in Parameter Estimation and Inverse Problems,” (Second Edition). 2013, Academic Press: Boston. p. 55–91.
36. I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Comm. Pure Appl. Math. 57(11), 1413–1457 (2004). [CrossRef]
37. O. Liba, M. D. Lew, E. D. SoRelle, R. Dutta, D. Sen, D. M. Moshfeghi, S. Chu, and A. de la Zerda, “Speckle-modulating optical coherence tomography in living mice and humans,” Nat. Commun. 8(1), 15845 (2017). [CrossRef]
38. M. Goodwin, B. Bräuer, S. Lewis, A. Thambyah, and F. Vanholsbeeck, “Quantifying birefringence in the bovine model of early osteoarthritis using polarisation-sensitive optical coherence tomography and mechanical indentation,” Sci. Rep. 8(1), 8568 (2018). [CrossRef]
39. K. A. Vermeer, J. Mo, J. J. A. Weda, H. G. Lemij, and J. F. de Boer, “Depth-resolved model-based reconstruction of attenuation coefficients in optical coherence tomography,” Biomed. Opt. Express 5(1), 322–337 (2014). [CrossRef]
40. M. Uchida, Y. Sun, G. McDermott, C. Knoechel, M. A. Le Gros, D. Parkinson, D. G. Drubin, and C. A. Larabell, “Quantitative analysis of yeast internal architecture using soft X-ray tomography,” Yeast 28(3), 227–236 (2011). [CrossRef]
41. C.L. Lawson and R.J. Hanson, Chapter 23. “Linear Least Squares with Linear Inequality Constraints, in Solving Least Squares Problems,” 1974, Englewood Cliffs, N. J., Prentice-Hall p. 158–173.
42. X. Shu, L. Beckmann, and H. Zhang, “Visible-light optical coherence tomography: a review,” J. Biomed. Opt. 22(12), 1 (2017). [CrossRef]
43. N. Lippok, S. G. Murdoch, K.-L. Wu, and F. Vanholsbeeck, “Dispersion mapping at the micrometer scale using tri-band optical frequency domain imaging,” Opt. Lett. 38(16), 3028–3031 (2013). [CrossRef]
