1. Introduction

Optical coherence tomography (OCT) is an attractive tool for probing the structure of biological tissues on a microscopic scale [1]. Ongoing efforts to reveal additional sample information from structural OCT images beyond their resolvable features include polarization-sensitive OCT [2,3], which leverages the polarization properties of the imaged sample but often requires additional imaging system hardware, and OCT attenuation coefficient imaging [4,5], which leverages the scattering properties of the sample based on the backscattered intensity signal alone. The latter is based on the principle that at typical near-infrared OCT wavelengths, light absorption is minimal, thus making the scattering coefficient $\mu _s$ a reasonable approximation for the light attenuation by the sample. $\mu _s$ is a function of the sizes, concentrations, shapes, and the relative refractive indices of scattering particles in a sample and their surrounding medium [6,7]. OCT attenuation coefficient imaging has been investigated as a source of potential contrast for many biological tissues and diagnostic applications, including characterizing atherosclerotic plaques [8–10], delineating brain tumor margins [11], and identifying dysplastic regions of Barrett’s Esophagus [12]. However, whereas the attenuation coefficient is capable of characterizing sub-resolution physical properties of a sample, many of these applications have employed it as a means of qualitative contrast, rather than use it to quantify diagnostically relevant microstructural parameters. Much of the existing body of work has computed attenuation by fitting the OCT A-line signal with exponentially decaying functions, modeled most simply as $I(z) \propto e^{-2 \mu z}$. Thus, the attenuation coefficient over a region that is homogeneous throughout its depth is linearly proportional to the intensity decay rate on a logarithmic scale. This approach is useful for generating en face maps of $\mu$, particularly in samples that have little variation in optical properties with depth, but offers limited utility and accuracy in samples with strong axial variation in composition. Many tissues of clinical imaging interest fall into the latter category, and thus there remains a need for a more generalized attenuation coefficient analysis.

More recently, a model was developed to calculate sample attenuation on a pixel-per-pixel basis in depth based on previous work in ultrasound and depth-compensated OCT [13–15]. This approach is advantageous for the analysis of layered tissues, where the usage of exponential fitting methods degrades spatial context, allows fitting artifacts to dominate [9,16], and renders the resulting attenuation maps difficult to read and interpret. Although this recent modeling is inherently capable of handling variable attenuation coefficients over a sample’s depth, there is a strong assumption of a constant backscattering fraction throughout sample depth. The backscattering fraction, $\alpha$, is defined as the ratio of backscattered light that is detected within the numerical aperture (NA) of the system to the light scattered by the sample in all directions [17]. Like the attenuation coefficient, it is dependent on the scattering cross sections of underlying sample particles, but whereas $\mu$ (and $\mu _s$) scale with both scattering particle concentration (i.e., the number density, or unitless volume fraction divided by the particle volume) and size, $\alpha$ depends only on particle size (with respect to the wavelength of the illuminating light) and refractive index relative to that of the medium. With $\alpha$ to determine particle size, $\mu$ is able to determine particle concentration, enabling a complete analysis of scattering particle composition. However, fixing the backscattering fraction in samples for which it varies along depth introduces errors in attenuation coefficient calculation. Previous work has shown that at a single site, backscattering fraction can vary strongly over different layers of tissue, such as wall and plaque components in coronary arteries [8]. Although it has been identified as a potentially useful parameter diagnostically in intravascular OCT and other applications where scatterer size is of diagnostic interest, accurately and robustly extracting backscattering fraction remains an important challenge.

In this work, we present a hybrid method for depth-resolved attenuation coefficient computation that analyzes the scattering behavior of different layers of a sample independently. Our analysis partitions each A-line of a sample into layers with $\alpha$ and $\mu (z)$ that fluctuate about a representative value. The integral of the backscattered intensity within that layer is proportionate to both the backscattering fraction and attenuation, but the decay rate depends on attenuation alone. This analysis allows us to decouple both contributions in a layer-based approach that, in conjunction with the depth-resolved model, allows for full-resolution computation of attenuation coefficient without influence from the backscattering fraction. Furthermore, unlike previous methods, our approach enables the calculation of a layer-resolved backscattering fraction. We present simulated and experimental results using our hybrid calculation method demonstrating the ability to determine $\mu (z)$ and $\alpha$ values on a layer-resolved basis. We validate our approach in layered phantoms with well-characterized scattering properties and compare the experimentally determined optical properties to theoretical predictions. We also show $\mu (z)$ and $\alpha$ results in porcine intestine ex vivo within the context of the corresponding histology. Notably, our approach increases the accuracy of $\mu (z)$ by reducing the deleterious impact of deep, multiply scattering tissue layers on more superficial, predominantly singly scattering layers. We attribute this benefit to the effective increase in acceptance angle—and thus increase in sample backscattering fraction—caused by multiple scattering. This effective increase in the backscattering fraction renders the conventional depth-resolved method inaccurate for all singly and multiply scattering layers, but using our methods to independently analyze each discrete sample layer, the impact of multiple scattering is confined to multiply scattering layers only. The evident improvement in attenuation coefficient calculation realized with this method stands to make $\mu (z)$ a more consistent and meaningful parameter for quantitative sample characterization, particularly for evaluation of layered tissues. The complementary information provided by backscattering fraction yields even more scattering-based insight into sub-resolution physical properties of tissue, creating an opportunity for additional diagnostic contrast.

2. Methods

2.1 Theory

We define $z$ as the depth coordinate, with $z=0$ defined at the tissue surface. In this Section, for simplicity, we assume that system-dependent effects on the OCT intensity signal, including the confocal function of the focusing lens, sensitivity roll-off, and depth-dependent noise floor have already been accounted for (see 2.6). Thus, the decay of measured backscattered intensity over depth in an arbitrary sample is given by [15]

In the more general case where neither attenuation nor backscattering fraction is constant over a sample’s depth, the same framework can no longer be used: in general, it is not possible to determine two unknowns—$\alpha (z)$ and $\mu (z)$—from a single measured quantity $I(z)$. Instead, however, the sample may be analyzed as a series of separate layers over which $\mu (z)$ fluctuates about a representative slope and $\alpha$ is reasonably constant

In general, we define layers as visibly distinct regions in OCT tomograms resulting from transitions between sample (e.g., tissue) architectures with different underlying compositions of scattering particles (e.g. cell bodies, collagen fibers). These different architectures result in $\alpha$ values that vary across several orders of magnitude (Supplementary Fig. S1). In contrast, $\mu (z)$ values realistically vary within a more limited range. Since $\mu (z)$ is dependent on particle concentration, it is both reasonable and advantageous to allow $\mu (z)$ to vary within tissue layers as a function of the variable density of tissue scatterers. Essentially, we attribute intra-layer fluctuations about the determined $\mu _{\textrm{fit}, x}$ at constant $\alpha$ to variation in $\mu (z)$ alone. We justify this for tissue layers with gradients in scattering particle density (such as epithelial tissues), and thus gradually varying $\mu (z)$, with minimally varying particle size or shape and thus minimal variation in $\alpha$. In general, we expect $\alpha$ to change dramatically only in discrete steps, representing different tissue layers, and any small intra-layer $\alpha$ fluctuation will only have minimal impact on the calculated $\mu (z)$ signal.

For the specific case of three layers $A$, $B$, and $C$ shown in Fig. 1(a), we write Eq. (6) explicitly as

 

Fig. 1. Visual overview of major steps within the layer-resolved $\mu(s)$ algorithm. Layers are identified on an A-line by A-line basis by detecting large intensity changes (a), then a linear slope, $\mu_{\text{fit}, l}$, is fitted to the intensity decay within each layer $l$ (b). The $\mu_{\text{fit}, l}$ determines the slope of extrapolated data that pads each layer’s signal down to the identified noise floor (c). $\mu(z)$ and $\alpha$ are calculated as detailed in 2.1 (d).

Download Full Size | PDF

Additionally, once a representative $\mu _{\textrm{fit}, l}$ has been calculated for each layer, it is possible to calculate a layer-resolved $\alpha _l$ based on the determined $\mu _{\textrm{fit}, l}$ for that layer by integrating over each region with constant $\alpha$. The definite integral of the intensity over a layer l is

The incident light intensity for each Layer $l$, $L_{0, l}$, is able to be described in terms of the incident light intensity on the previous (i.e., immediately more superficial) layer and the total signal attenuation by that layer. To analyze each layer independently, while still taking the overall $L_0$ attenuation into account, we define an axial coordinate system that begins at the start of each layer, where $\tilde {z}_l = z - z_{l-1}$ (and $z_{l-1}$ is the location of the top of the layer). For example, $\tilde {z}_A = z - 0$, $\tilde {z}_B = z - z_A$, and $\tilde {z}_C = z - z_B$. Thus, $L_{0, B} = L_{0, A}e^{-2 \mu _{\textrm{fit}, A} \tilde {z}_B}$ and $L_{0, C} = L_{0, B}e^{-2 \mu _{\textrm{fit}, B} \tilde {z}_C}$, and $L_{0, A}$, the incident light intensity on the first layer which has not undergone any attenuation, is merely the $L_0$ measured for the system. Therefore, $\alpha _A$ for the first layer may be calculated as

2.2 Implications of multiple scattering

The ability to extract accurate particle size and concentration information from $\alpha$ and $\mu$ based on Mie theory is reliant on the assumption that detected light has only experienced a single scattering event. However, the likelihood of detecting multiply scattered light increases with increasing particle size, which favors forward scattering, and increasing particle concentration, which favors additional scattering events [18]. Many biological tissues of clinical imaging interest are known to exhibit multiple scattering, such as lipid rich atheroma, and sites of dysplastic or neoplastic transformation in tissue epithelial layers. The latter example is accompanied by increasing nuclear diameter and hypercellularity (e.g., an effective increase in scattering particle size and concentration) [19]. Although multiple scattering (MS) is known to impact OCT intensity signals despite coherence and confocal gating, it is difficult to definitively determine its incidence using OCT systems with conventional detection optics [20], and few models exist to quantify its impact on the intensity signal [21,22].

In general, multiple scattering couples additional light that would not be detected under single scattering conditions into the OCT system. The effective path length of the multiply scattered signal is longer than the path length of a singly scattered signal from the same point, thus the impact of MS is expected to increase with depth into the tomogram. We hypothesize that the effective scattering angle of multiply scattered light is also larger than that of singly scattered light from the same point; it may incorporate signal at a scattering angle that may exceed the OCT system NA and not be detected under singly scattering conditions.

Avoiding inter-layer $\mathbf {\boldsymbol \mu (z)}$ dependence with our approach. We expect multiple scattering to impact $\mu (z)$ and $\alpha$ measurements predictably due to these effects. With greater coupling of light into the system under MS conditions, we expect $\mu (z)$ to be lower than theoretically predicted, and to decrease further with greater depth into a tomogram. Due to the effective increase in scattering angle, $\theta$, of the detected light, we also expect $\alpha$ measurements to be higher than predicted in a singly-scattering regime (see Eq. (2)). However, although the measured $\mu (z)$ and $\alpha$ values for a multiply scattering layer will still not accurately portray its underlying particle makeup, using our methods, the impact of multiple scattering is confined to just this layer. Using the conventional depth-resolved method, the required integration of the intensity in depth across multiple layers creates an interdependence between them, and thus, when $\alpha$ varies between layers, $\mu (z)$ is under- or over-estimated in other layers with respect to true values, as shown in Fig. 2. We note that this error affects any situation in which $\alpha$ is variable with depth, not only in the event of multiple scattering. Therefore, we emphasize that the advantages of our $\mu (z)$ calculation methods, in the context of multiple scattering, are able to realize accurate scattering parameter calculation in singly-scattering layers in the presence of multiply-scattering layers, whereas in the conventional method, even the $\mu (z)$ values in the singly-scattering layers are rendered inaccurate.

 

Fig. 2. Simulated data demonstrating the effect of variable $\alpha$ on $\mu(z)$ calculation. Changes in $\alpha$ are visible as abrupt changes in intensity in both the representative A-line (a) and B-scan (b). Calculated $\alpha$ values, matching the input $\alpha$ values, are shown (hue) overlaid onto the intensity (luminance) for the B-scan (c). Calculated $\mu(z)$ values are shown for both a single A-line (d) and as a B-scan overlay, using both the layer-based (“Corrected”) method (e) and the conventional (“Uncorrected”) depth-resolved method (f).

Download Full Size | PDF

2.3 Phantoms and tissue

Well-characterized, agarose-based phantoms of homogeneous polystyrene bead suspensions or Intralipid were constructed, sliced thinly, and stacked to achieve the desired layered geometry. Phantoms were composed of 0.35 µm (volume fraction $v_f = 0.002$) and 1 µm ($v_f = 0.004$) diameter polystyrene beads, and 10 % w/v Intralipid ($v_f$ = 0.023 [23]). Their compositions were chosen to result in Mie theory-derived scattering properties representative of different scattering regimes found in biological tissues, ranging from low to moderate (predominantly singly scattering) to high (multiply scattering). All theoretical $\mu$ for relevant samples were calculated using Mie theory. $\mu$ calculations were additionally augmented by concentration-dependent and Extended Huygens-Fresnel analysis [17], but given the low volume fractions used for the bead phantoms, the resulting values were negligibly different from Mie theory alone. For tissue experiments, porcine intestine samples (Lampire Biological Laboratories, US) were imaged ex vivo, fixed, and processed with standard hematoxylin and eosin (H&E) staining.

2.4 OCT system

A frequency domain OCT system, operating at 100 kHz A-line rate and having a spectral full-width-at-half-maximum (FWHM) bandwidth of 90 nm centered at 1300 nm was used for all measurements. The instantaneous coherence length of the laser (Thorlabs, SL131090) was reported to be greater than 1 m, and thus the OCT system exhibits a minimal roll-off due to the particular frequency bandwidth of the electronics. The sensitivity of the system was measured using a mirror and a density filter. A value of 111 dB was calculated as the ratio of the peak mirror intensity to that of the noise floor adjacent to the peak, and the FWHM point spread function of the peak was determined to be 10.8 µm. A $k$-clock generated by the source was used to acquire fringes linear in $k$, and a central-wavelength trigger generated using the source’s internal Mach-Zehnder interferometer was used for synchronization. In the sample arm, a galvanometer scanner (Cambridge Technology/Novanta Corporation, US) scanned the beam focused with a 5X objective lens (LSM03, Thorlabs, NA = 0.05 with 3.4 mm collimated beam diameter) resulting in a $1/e^2$ beam diameter of 21 µm and a FWHM lateral resolution of 10 µm. Following interference with the reference arm, light backscattered from the sample passed through a polarization-diverse receiver (Advanced Fiber Resources, China) to enable balanced polarization diverse detection with high-speed photodiodes (1.6 GHz, Thorlabs). The signal was digitized with a two-channel acquisition card (AlazarTech, Canada; ATS 9373, 4 GS/s) with additional ports for $k$-clock and wavelength-trigger inputs. Data were visualized in real time and recorded using custom-written, in-house LabVIEW software.

2.5 OCT signal simulations

In addition to experimental measurements, OCT intensity data were simulated for validation of $\mu (z)$ and $\alpha$ calculation methods. Layered samples were simulated by setting input $\mu$ and $\alpha$ for each layer to values within ranges expected to be encountered in scattering phantoms or tissue. The input values for $L_0$ and noise floor were set to match experimental conditions. A noise floor of 80 dB magnitude was added to the signal to mimic experimentally acquired data. Speckle was not simulated, as experimental data required speckle mitigation by averaging or other despeckling techniques (as discussed in Sec. 2.6).

2.6 Data analysis

All data were analyzed in MATLAB. Following conventional tomogram reconstruction, the effect of the confocal function, $H(z)$, on the intensity signal was removed from the tomogram using methods based on those published by Dwork et al. [24]. This necessitated acquiring two tomograms of the same sample location at different distances from the objective lens to find the focal depth and Rayleigh range, $z_f$ and $z_R$, without a priori knowledge of $H(z)$. We modified the published methods by adding the distance between the focal planes in the two tomograms as a third fitting parameter, rather than detecting and fixing it. This modification allowed us to account for variations in sample refractive index and improved the robustness of the technique. Prior measurements using the described OCT system confirmed a flat noise floor and negligible sensitivity roll-off over the ranging depth, eliminating the need for these additional pre-processing steps. Attenuation coefficient maps were then computed from the confocal-corrected and out-of-plane averaged data (averaged over $N$ B-scans, where $N=64$ for phantoms, and $N=24$ for tissue) using both the conventional depth-resolved algorithm [Eq. (5)] as well as our approach (detailed in Sec. 2.1), which also was used to generate maps of the sample backscattering fraction. For visualization, attenuation and backscattering maps were overlaid onto the original intensity tomogram using the hue of the colormap to display the magnitude of either $\mu (z)$ or $\alpha$, and the luminance of the colormap to display the corresponding tomogram intensity [25]. Panels for which this overlay is implemented are accompanied by a two-axis colorbar with luminance scaled according to the corresponding intensity tomogram.

Automated layer segmentation. The intent of the layer-based attenuation calculation approach presented here is to avoid inaccuracies in attenuation estimation due to differences in backscattering fractions of individual layers of a sample while maintaining depth-resolution, as well as to compute the backscattering fractions of those layers. Since this is accomplished by analyzing the attenuated intensity in layers independently, accurate detection of layer boundaries is a critical first step in the analysis of the pre-processed tomograms. The built-in MATLAB function findchangepts (statistical toolbox) was used to generate a set of boundary points for each A-line, $I_j(z)$, based on our definition of a "layer" as a sudden dramatic difference in intensity, resulting in a series of sufficiently long, depth-wise regions of relatively constant attenuation. For the data shown in this paper, the minimum length of these depth-wise regions was 20 pixels, or 108 µm. After the detection of each boundary point, the average attenuation ($\mu _{\textrm{fit}, l}$) is computed for the corresponding layer l, which begins at depth $z_{l-1}$ and ends at depth $z_l$. Then, $\mu _{\textrm{fit}, l}$ is used to extrapolate the intensity of that layer down to the noise floor for $z > z_l$. The signal was extrapolated as $Ae^{-2\mu _{\textrm{fit}, l}z}$, where $A$ and $\mu _{{\textrm{fit}}, l}$ were determined by linearly fitting $b + mz$ to log($I(z_{l-1} < z < z_l)$), and by calculating $\mu _{\textrm{fit}, l} = \frac {1}{2}m$ and $A = \exp {(b)}$. As a result of the extrapolation, the effective intensity at the final point of the extrapolated decay, $D$, is 0, allowing the proper implementation of the traditional depth resolved algorithm, $\mu (z) = \frac {I(z)}{2 \int _z^{D} I(u)du}$. The extrapolated points are then discarded, and the relevant region of calculated attenuation (excluding extrapolated points) is placed back in its original position in the tomogram. We note that the approach described above is only one example of automated layer detection, and whichever image segmentation strategy is most appropriate for a given application should be used. For example, analysis in retina tomograms might be carried out more successfully using existing retina segmentation algorithms that already account for the (relative) uniformity of each conventionally-defined layer and the expected contrast between them. These processing steps, as well as those described below, are summarized in Table 1.

To calculate $\alpha$, a 10-pixel (50 µm) saturation artifact at the sample surface was removed. This helped to reduce the impact of specular reflections, which are not accounted for in the theoretical model. Then, the approach described in Sec. 2.1 was applied to calculate $\alpha$ for each layer. This first required calculating the incident light on each layer, $L_{0, l}$. For the first layer, $L_{0, 1}$, an experimentally determined value of 2.7E16 counts-mm was used for the combined, system-specific $L_0 \beta$ term (where, consistent with previous work, $\beta$ describes the overall system detection efficiency [15]). This value was calculated from a sensitivity measurement of the system as the integration of a mirror signal over its full width in microns. It includes the offset of a noise floor of 78 dB magnitude, as no adjustment was made to the intensity values in transforming raw fringes to tomograms. A pixel-to-mm conversion factor was used to account for substantially increased zero padding during peak reconstruction. For each subsequent layer, $L_{0, l}$ was attenuated at the rate of the previous layer, $\mu _{\textrm{fit}, l-1}$, over the length of the previous layer. $\alpha _l$ was then calculable as

Table 1. $\mu (z)$ and $\alpha$ signal processing steps

View Table

A custom segmentation routine (Fig. S2) was used to quantify local nuclear diameter and concentration from high-resolution histology images (40X magnification) of the H&E-stained tissue sections. Nuclei were identified using thresholding in the blue image channel followed by filtering steps. Watershed segmentation was used to isolate nuclei in close proximity to each other and calculate the area of each. Nuclear diameters were calculated assuming an approximately circular particle shape from the particle area. The nuclear-cytoplasmic ratio (NCR) was approximated in epithelial layers by dividing the local area of tissue occupied by nuclei by that occupied by tissue, with any acellular regions (including area within glands, gaps in the tissue, and area beyond edges) excluded. This is a reasonable approach in cell-dense epithelial layers for estimating NCR without needing to segment and measure the cytoplasmic area of individual cells, and is a direct surrogate for our measurement of particle concentration (to which $\mu (z)$ is sensitive).

3. Results

We present a series of experiments to motivate, validate, and evaluate our $\mu (z)$ and $\alpha$ calculation approach. We first applied our approach to simulated OCT intensity data to demonstrate the impact of axially varying $\alpha$ on $\mu (z)$ calculation and our approach’s ability to preserve accurate $\mu (z)$ extraction. We next applied it to experimentally measured OCT data from tissue-mimicking solid phantoms with known scattering properties, validating its usage in real samples. With this data, we confirmed both the greater accuracy of our method in determining $\mu (z)$ and the expected variation in $\alpha$ with particle size with Mie theory. Finally, we applied our methods to samples of excised porcine intestine to demonstrate that our calculated $\mu (z)$ better matches histologically-determined nuclear-cytoplasmic ratio, and better mitigates against the impact of deep, multiply scattering tissue layers on the epithelial $\mu (z)$ signal, than the $\alpha$-uncorrected $\mu (z)$.

To initially validate $\mu (z)$ and $\alpha$ calculation methods, we simulated the intensity signal for a layered sample using physically reasonable attenuation and backscattering fraction values for each layer. Each distinct layer is characterized by a distinct $\alpha$ and a distinct $\mu _{\textrm {fit}}$, representing different combinations of scattering particle size and concentration (see Fig. S2). A representative A-line and the simulated tomogram are shown in [Figs. 2(a) and (b)], respectively. The attenuation coefficient $\mu (z)$ is calculated for the simulated tomogram using both the conventional ($\alpha$-uncorrected) depth-resolved method and our layer-based approach. Directly due to higher backscattering fractions [Fig. 2(c)] of the second and third layers, $\mu (z)$ is greatly underestimated in the first and second layers using the $\alpha$-uncorrected approach. This error results from the difference in $\alpha$ between these two layers, and the inter-layer $\mu (z)$ calculation dependence using the conventional approach, as discussed in 2.1 ("Implications of multiple scattering"). In Fig. 2, $\mu (z)$ is displayed both for a single A-line [Fig. 2(d)], with respect to the input $\mu$ values for each layer, and also for the full tomogram, both corrected [Fig. 2(e)] and uncorrected [Fig. 2(f)] for the variable backscattering fraction with depth. Additionally, the backscattering fraction values calculated from the simulated tomogram using the approach (outlined in Section 2.1) matched their corresponding input values. Both $\alpha$ and $\mu (z)$ are constant throughout each sample layer, but are displayed, here and in Fig. 3, 6, and 7, as overlays on the intensity tomogram (as described in Section 2.6). Previously published simulations using the conventional depth-resolved $\mu (z)$ calculation did not test its validity in the context of axially varying sample $\alpha$ values, thus this effect has not been reported previously. Notably, an intra-layer gradient in $\mu (z)$ is introduced, most prominently in the second sample layer, due to the inter-layer $\alpha$ variation. This results in the suggestion of a strongly varying $\mu (z)$ within this layer, despite a constant input $\mu$ and visibly linear intensity decay. This error highlights that the ramifications of neglecting variable $\alpha$ range from not only constant offsets in value, such as the great underestimation of $\mu (z)$ in the first layer, but also unexpected $\mu (z)$ fluctuation within homogeneous layers.

 

Fig. 3. Quantifying $\mu(z)$ and $\alpha$ in two-layer phantoms. Representative A-lines are shown for Phantom 1 (0.35 $\mu$ m beads / 1 $\mu$m beads) and Phantom 2 (0.35 $\mu$ m beads / Intralipid) (a), as well as their respective B-scans (b, c) and $\alpha$ overlays (d, e). $\mu(z)$ values for a single A-line in each sample are provided in (f), as well as $\mu(z)$ values for matching homogeneous bead phantoms. $\mu(z)$ overlays for Phantoms 1 and 2 are calculated with (g, h) and without (i, j) using the layer-based method. Scale bars in (c) indicate 0.5 mm. HG = homogeneous.

Download Full Size | PDF

Two solid, two-layer phantoms were constructed and used to investigate the applicability of the layer-based methods to experimental data. The top layer of each phantom, consisting of small (0.35 µm) polystyrene beads at a low concentration, was intended to be representative of a well-behaved, predominantly singly and isotropically scattering sample. The bottom layer of the first phantom, consisting of larger (1 µm) beads, should exhibit more forward scattering and a lower backscattering fraction than the smaller beads, though still remain predominantly singly-scattering. The bottom layer of the second phantom, 10% Intralipid, is expected to have a larger backscattering fraction than either of the bead layers based on the distribution of particle sizes [23]. Despite scattering isotropically due to small average particle size, however, the high particle concentration in this layer should give rise to multiple scattering.

Representative A-lines and intensity tomograms for each layered imaging phantom are shown in Fig. 3(a) and Fig. 3(b), respectively. The backscattering fractions [Fig. 3(d, e)] computed for the identical first layer are in good agreement between the two phantoms, and the larger bead (lower $\alpha$, Phantom 1) and Intralipid (higher $\alpha$, Phantom 2) layers behave as expected. We note that although there is a spike in intensity at the interface of the highly attenuating Intralipid layer, $\alpha$ is nonetheless constant throughout this layer. Figure 3(f) plots the corresponding $\mu (z)$ with depth for the A-lines shown in Fig. 3(a) for each phantom using both the layer-based and traditional depth-resolved $\mu$ calculation methods. Additionally, $\mu (z)$ as measured over depth in a homogeneous phantom matching each of the bead layers is plotted at the corresponding layer depths. Here and elsewhere, the $\mu (z)$ for the homogeneous Intralipid phantom is excluded, as the high degree of multiple scattering made it difficult to interpret, which can be appreciated in the $\mu (z)$ plotted for Phantom 2 [Fig. 3(f)]. In both layered phantoms, the $\mu (z)$ calculated using our $\alpha$-corrected method produces a comparable $\mu (z)$ map for the first layers, identical in each phantom, as expected [Fig. 3(g, h)]. However, for Phantom 1, the lower backscattering fraction of the bottom layer with respect to the upper layer results in an overestimation of the $\mu (z)$ in the upper layer [Fig. 3(i)] using the conventional, "uncorrected” depth-resolved method. Conversely, the higher backscattering fraction of the bottom layer in Phantom 2 results in an underestimation of the uncorrected $\mu (z)$ in the upper layer [Fig. 3(j)].

Figure 4 compares experimentally determined attenuation coefficient and backscattering fraction values to Mie theory predictions. The $\mu (z)$, as computed with both the corrected and uncorrected depth-resolved methods, and $\alpha$ values for each layer of the phantoms were averaged over the full region of interest (ROI) shown in Fig. 3. They were compared to both $\mu$ and $\alpha$ values as computed by Mie theory, as well as the $\mu$ and $\alpha$ values calculated from the homogeneous phantoms made to match each individual layer. Again, attenuation coefficient analysis of the Intralipid sample was omitted due to strong influence of concentration-based multiple scattering. The good agreement between the $\mu$ values as calculated by Mie theory, from homogeneous phantoms, and using our layer-resolved method, for both of the bead layers, supports the greater accuracy of the corrected method [Fig. 4(a)]. Similar to simulation results (Fig. 2), the under- and over-estimation of $\mu$ in the top layer of each phantom using the conventional depth-resolved method, and slight underestimation of $\mu$ in the bottom layer of Phantom 1, is apparent [Fig. 4(a)]. In the backscattering fraction analysis, $\alpha$ values were consistently lower than those predicted by Mie theory [Fig. 3(b)]. Possible reasons for this discrepancy are considered in the Discussion. However, the overall trend of increasing $\alpha$ with decreasing particle size was in agreement with expectations, and the consistency of results between the homogeneous and layered phantoms was encouraging.

 

Fig. 4. Comparisons between theoretical and experimental $\mu$ and $\alpha$ values. In the left panel of (a), the Mie theory-determined $\mu$ value for 0.35 $\mu$m beads at their prepared concentration is compared to $\mu$ values experimentally measured in a homogeneous 0.35 $\mu$m bead phantom, in the top layer of Phantoms 1 and 2 using the layer-based $\mu(z)$ method (“Corrected”), and analogously using the conventional (“Uncorrected”) method. Similar comparisons are made for 1.0 $\mu$m beads in a homogeneous phantom and the bottom layer of Phantom 1. In (b), $\alpha$ values calculated theoretically using Mie theory using Eq. (2) and experimentally using Eq. (12) in homogeneous phantoms and in relevant portions of layered phantoms are shown for each bead size and Intralipid.

Download Full Size | PDF

To investigate the impact of accurate consideration of variable $\alpha$ when calculating layer-resolved $\mu (z)$ in tissue specimens, porcine intestine samples were imaged ex vivo analogously to the imaging phantoms, and submitted for histological processing. The resulting high-resolution digital pathology images of Hematoxylin and Eosin (H&E)-stained tissue sections within the imaging ROI were additionally analyzed using the described custom nuclear segmentation scheme to generate maps of regional nuclear diameter and effective nuclear-cytoplasmic ratio (NCR; here, $\frac {\textrm {nuclear area}}{\textrm {tissue area}}$) as proxies to scattering particle size and concentration.

In Fig. 5, the top row consists of the H&E-stained image (a), nuclear diameter map (b), and NCR map (c) for a histologically normal porcine intestine specimen. Throughout the epithelium, which is our primary layer of interest, nuclear diameter is reasonably constant, and the NCR has some fluctuation, mostly related to the presence of glands. The tissue specimen shown in the bottom row (d-f) contained Peyer’s Patches, or pockets of densely-packed immune cells that are commonly found in mammalian lower gastrointestinal tracts. Throughout the superficial epithelium and the Peyer’s Patches, the diameter again varies little and randomly, whereas the NCR is greatly elevated in the immune cell deposits, resulting in a high particle concentration and effective NCR.

 

Fig. 5. Histological analysis of imaged porcine intestinal tissues. The high-resolution digital histology image of H & E-stained fixed tissue (a, d), mapping of regional trends in nucleus diameter (b, e), and mapping of regional trends in effective nuclear-cytoplasmic ratio (NCR) (c, f) are shown for the two tissue samples imaged and analyzed in Fig. 6 (top row) and Fig. 7 (bottom row). The epithelium (Epi), lamina propria (LP), and Peyer’s Patches (PP) are indicated.

Download Full Size | PDF

Figure 6 shows the corresponding intensity (b), attenuation (e, f), and backscattering fraction (c) maps for the tissue sample featured in the top row of Fig. 5(a–c). A representative A-line showing automatic identification of layers and corresponding exponential fits is shown in Fig. 6(a) (and for a more complete view of the automated segmentation results, see Figure S3). The epithelial layer, bordered by the fibrous, highly-backscattering lamina propria layer, has little contrast in the intensity tomogram, $\alpha$ map, and $\mu (z)$ map as calculated using the uncorrected approach. However, using the layer-resolved approach results in not only higher $\mu (z)$ values throughout the epithelium (d, e), due to accounting for the apparent higher $\alpha$ of the lamina propria, but also greater contrast within the epithelium itself that was otherwise absent in the uncorrected $\mu (z)$ map. The higher $\alpha$ of the lamina propria layer explains the apparent underestimation of $\mu (z)$ in the uncorrected map. Overall, the results are supported by the corresponding histology in Fig. 5 showing the natural variation in NCR throughout the epithelium, which should drive the variation in $\mu (z)$ visible in the corrected map.

 

Fig. 6. Attenuation coefficient and backscattering fraction analysis in a sample of porcine intestinal tissue. A representative A-line and B-scan of the sample are shown in panels (a) and (b), with the plotted A-line location provided by the cyan line in (b). Layer distinctions between the epithelium (“Epi”) and lamina propria (“LP”) are made in (b). The calculated $\alpha$ overlay given in (c). $\mu(z)$ plots for a single A-line using both our layer-based (“Corrected”) and conventional (``Uncorrected'') approaches are given in (d), and accompanied by corresponding Corrected (e) and Uncorrected (f) $\mu(z)$ overlays.

Download Full Size | PDF

Figure 7, which corresponds to the tissue containing the Peyer’s Patches, shows a more extreme example of backscattering fraction variation with depth in tissue. Segmentation and fitting is shown for a single A-line in Fig. 7(a) (and for a more complete view of the automated segmentation results, see Figure S3). In the uncorrected attenuation coefficient map [Fig. 7(f)], the high $\alpha$ observable in the lymphoid tissue [Fig. 7(c)] causes such a substantial underestimation of $\mu (z)$ in the superficial epithelium that all attenuation-based contrast is lost and the resulting $\mu (z)$ values are unrealistically low. This contrast is recovered using the layer-based attenuation calculation approach, where more realistic $\mu (z)$ values are retrieved [Fig. 7(d, e)]. Both the $\alpha$ and $\mu (z)$ values in the region superior to the Peyer’s Patches are similar to those in the epithelium of Fig. 6. This is supported by histological results (Fig. 5), which feature similar nuclear diameter and NCR values throughout these regions in each tissue section.

 

Fig. 7. Attenuation and backscattering analysis in a second sample of porcine intestinal tissue containing lymphocyte-rich regions. A representative A-line and B-scan of the sample are shown in panels (a) and (b), with the plotted A-line location provided by the cyan line in (b). Layer distinctions between the epithelium (“Epi”) and Peyer’s Patches (“PP”) are made in (b). The calculated $\alpha$ map given in (c). $\mu(z)$ plots for a single A-line using both the layer-based and traditional approaches are given in (d), and accompanied by corresponding corrected (e) and uncorrected (f) $\mu(z)$ maps.

Download Full Size | PDF

The origin of the lymphoid tissue’s high backscattering fraction may be considered in the context of the scattering particle analysis from the segmented histology images [Fig. 5(d-f)]. From Mie theory, $\alpha$ varies with particle diameter, system NA, refractive index, and medium refractive index (assuming a constant central illumination wavelength). More empirically, $\alpha$, like $\mu (z)$, is sensitive to strong multiple scattering. Under these conditions, detected light may have undergone forward- and side-scattering events before eventual backscattering and coupling into the OCT system. This effectively broadens the angular range of backscattered light collection (or system NA), thus also increasing $\alpha$. Previous optical analyses of lymphocytes put their refractive index parameters within the normal range of epithelial cells [26], suggesting that neither particle nor medium refractive index drove the large increase in $\alpha$ within the Peyer’s Patches. Furthermore, Fig. 5(d-f) demonstrates that the transition into lymphoid tissue was not accompanied by a change in nuclear, or scatterer, diameter, but rather an increase in the effective nuclear cytoplasmic ratio. As discussed in Sec. 2.1, multiple scattering is expected to cause an increase in measured $\alpha$ due to greater effective angular acceptance of light. The high effective NCR is analogous to high particle concentration in the lymphoid tissue, and is likely to have caused a substantial amount of multiple scattering in this region. The alternative explanation for the very high $\alpha$, a marked decrease in particle diameter with respect to the upper layer, was ruled out by the quantitative histological analysis. Therefore, we attribute the high $\alpha$ primarily to multiple scattering, an effect likely also seen in the Intralipid portion of Phantom 2 (Fig. 3). This suggests that multiple scattering will result in an underestimation of nuclear size and an overestimation of nuclear concentration when these parameters are attempted to be extracted from $\alpha$ and $\mu$ directly.

4. Discussion

We present a new OCT-based intensity attenuation coefficient analysis method. Like the previously published depth-resolved method [14,15], our method preserves full-resolution image context of the original tomogram, but also retains the accuracy of the exponential-fitting method, with respect to axial variation of tissue backscattering fraction. Furthermore, it enables the calculation of a layer-resolved backscattering fraction for samples with layers of sufficient homogeneity. To do so, we separate A-lines into layers of presumed similar $\mu$ and $\alpha$ to be analyzed independently on the basis that each parameter contributes uniquely to the amplitude of the intensity signal and its decay. Our approach accounts for the overall light attenuation of more superficial layers in calculating $\alpha$ for a given layer, which does not allow deeper layers, notably those with a substantial contribution of multiple scattering, to influence $\mu (z)$ determination in superficial layers. This latter point makes our approach for $\mu (z)$ calculation particularly advantageous in many tissues of clinical imaging interest over the conventional depth-resolved method.

For proper implementation of this attenuation calculation algorithm and others, correcting for instrument-induced intensity variation along depth is paramount. This includes accounting for the axial point spread function, sensitivity roll-off, and depth-dependent background signal. Although the latter two effects may be mitigated using a one-time calibration, a stricter method to correct for the confocal function, such as the one suggested by Dwork et al., is required [24]. This method has the advantage of greater flexibility in sample placement and instrument geometry, but comes at the cost of acquiring two consecutive measurements with the sample at different heights. Though this is easily achieved in a benchtop setting, it is less desirable for endoscopic or catheter-based in in vivo applications, for which accurate knowledge of system NA and the location of the focal plane in the sample, or a calibration measurement, can be used to eliminate the effect of the axial PSF. Otherwise, given that necessary corrections related to depth-wise intensity variation can be made, the approach presented here can be used entirely in post-processing to evaluate existing data, as it requires no a priori information about the imaged tissue, or hardware supplementation. For an absolute measurement of $\alpha$, it is necessary to measure the incident light irradiance ($L_0$), but even without an $L_0$ measurement, the relative differences in backscattering fraction between layers of a sample can be determined. For instance, when areas of healthy tissue are imaged in a scan, the $\alpha$ signal in all areas can be compared to that of the healthy area to serve as a reference, making the relative information potentially useful without having a rigorously calibrated $L_0$ for precise $\alpha$ values.

Previous works that have used OCT-derived attenuation coefficients, both exponentially fit and depth-resolved, to characterize tissue morphology have been sensitive to variation in scatterer, e.g. nucleus, concentration and size, as $\mu$ increases with both of these parameters. This creates ambiguity over which one of these two principal parameters is the primary driver of the change in $\mu$. Determining both $\mu (z)$ and $\alpha$ throughout a sample could help reduce ambiguity. For example, measuring a change in attenuation coefficient within a region of consistent backscattering fraction would be indicative of a change in particle concentration rather than of size, a conclusion unable to be drawn by measuring either parameter alone.

However, a similar ambiguity exists in the application of the layer-based method. The change in intensity between two adjacent points in depth could be caused by either a change in $\mu (z)$ or in $\alpha$, though the algorithm presented here attributes such changes to $\mu (z)$, under the assumption that $\alpha$, which varies on a much greater scale than $\mu (z)$, can be assumed to be relatively constant over a layer of sufficiently homogeneous composition. Additionally, small changes in $\alpha$ have minimal impact on $\mu (z)$. As shown in Figure S1, the linear variation of $\mu (z)$ and exponential variation of $\alpha$ support the analysis of phantoms and tissues as $\alpha$-differentiated layers. Although it is possible to use $\mu (z)$ to calculate a per-pixel $\alpha (z)$, this violates the decoupling of $\mu$ and $\alpha$ as the slope of the intensity decay and its inter-layer changes, respectively. Without this decoupling, the theory presented in Methods is no longer applicable from the mathematical standpoint that allows the independent treatment of each layer. An interesting case arises when there is a very strong axial gradient in $\mu (z)$ in the absence of any change in $\alpha$. This could result in a layer of initially consistently increasing intensity, which cannot be accounted for using our current methods, which require a net decrease in intensity over each identified layer. It is unclear how often this particular phenomenon occurs in tissue, but efforts to calculate accurate $\mu (z)$ and $\alpha$ under these conditions may merit further exploration.

Quantitative scattering analysis is also complicated by multiple scattering, as evidenced by analysis of the intestinal tissue in Fig. 7. Multiple scattering can be modulated by aspects of the imaging procedure apart from the sample itself, including the coherence gate and the confocal gate of the OCT system, the impacts of which are investigated in Figure S4. These effects may even be combinatorial, including variable dependence on the NA of the imaging system (and thus the confocal gate) based on the sizes of the scattering particles and their scattering anisotropies [27]. However, even for very large scatterers such as red blood cells, which exhibit very strongly forward scattering behavior, the confocal gating still plays an important role in regulating the amount of detected multiply scattered light, as seen in OCT angiography, where higher NAs translate into smaller decorrelation tails [28,29], despite the knowledge that these tails are caused by strong forward scattering [30]. Our results demonstrate that multiple scattering causes an artificial increase in attenuation coefficient, which appears to be accompanied by an increase in backscattering fraction. It is only with histological validation that we can see that this is not due to a decrease in particle size and commensurate increase in particle concentration, but rather multiple scattering resulting from a strong scattering particle concentration increase. Although this effect frustrates the extraction of meaningful optical properties from the second detected tissue layer, it highlights the importance of using the layer-based approach. Without it, the $\mu (z)$ signal in the epithelial layer of both tissue samples is visibly underestimated, preventing meaningful morphological interpretation. By correcting for depth-wise variation in $\alpha$, the integrity of the attenuation signal in the epithelial layer of the tissue samples is preserved. This correction may be critical for the interpretation of $\mu (z)$ in many clinically relevant epithelial sites of interest that are superior to deeper, more predominantly multiply scattering layers. Excepting the Extended Huygens Fresnel model, OCT signal analysis lacks tools for meaningful interpretation of multiply scattering tissue regions, but with the layer-based approach, its effect on singly-scattering regions may be tempered.

Despite good agreement between $\mu$ values calculated by Mie-theory and $\mu (z)$ measured in phantoms, calculated $\alpha$ values were substantially larger than those measured using system $L_0$ as derived from a sensitivity measurement. As only a small fraction of the difference may be attributed to specular reflection loss, further investigation is necessary to discover the source of this discrepancy. However, the consistency of $\alpha$ values across independent phantom measurements and the expected increase in $\alpha$ with a decrease in particle size, are encouraging. Although we do not expect $L_0$ to fluctuate substantially enough across measurements to cause the observed discrepancies in $\alpha$, future measurements will normalize $L_0$ to signal from a calibration mirror (benchtop setup) or a catheter sheath.

Compared to prior approaches for calculating attenuation coefficients, based either on fitting an A-line to a single exponential decay or using Eq. (5) for a depth-resolved approach, the separation of intensity decays into discrete layers could significantly complicate the application of our methods. To avoid developing customized segmentation algorithms for individual applications, in this work, we have implemented a flexible approach for layer determination based on the automated identification of axial regions with similar backscattering fractions and attenuation coefficients. However, we acknowledge that this approach may require adjustment and refinement for general use, and it may be more effective to leverage preexisting segmentation tools for some applications for which these are available, such as retinal OCT. Further image processing challenges remain concerning the definition of a “layer” and the handling of noisy regions that defy linear fitting. Though out-of-plane averaging was employed here to facilitate layer-identification and fitting steps of the algorithm, more sophisticated despeckling tools are advantageous when averaging is not an option [31].

Our approach is attractive for clinical applications examining tissues with discrete layers which are known to, or suspected to, vary strongly in $\alpha$, either due to intrinsic architectural differences or the influence of multiple scattering. This includes the gastrointestinal samples shown in this work, supporting the use of our methods in gastrointestinal endoscopic applications, particularly those investigating dysplasia and neoplasia. Additional relevant clinical targets that could benefit from our methods include coronary arteries, for which $\alpha$ has been previously reported to vary for different plaque and arterial wall components [8], and retinal layers, for which there is interest in analyzing the scattering behavior of superficial layers [32,33], which is likely impacted by the strong scattering of the retinal pigmented epithelium [34]. In general, we hope that the advantages of quantifying $\alpha$, a more accurate $\mu (z)$, and mitigation of multiple scattering that our approach promises, at the expense of additional image processing steps and remaining ambiguity between independent contributions of $\alpha$ and $\mu (z)$ to the overall intensity signal, will extend the accurate application of scattering property-based analysis to these and other tissues for more informative clinical assessment of tissue architectures.

5. Conclusion

In this work, a new approach to calculating depth-resolved attenuation coefficients was explored to avoid errors present in previously proposed approaches related to variation in sample backscattering fraction along depth. This layer-based approach additionally enables the calculation of a layer-resolved backscattering fraction. The methodology was first validated using simulated data and well-characterized imaging phantoms, then found to be practically applicable to tissue analysis. Based on these initial tissue measurements, the layer-based algorithm improves depth-resolved attenuation calculation in predominantly singly-scattering regions such as most epithelia, mitigates against the influence of multiple scattering from inferior layers, and provides backscattering fraction as an additional source of contrast to be explored. We hope that this approach will prove useful in interrogating tissue morphology for the differentiation of healthy and diseased specimens.