1. Introduction

In direct analogy to ultrasound imaging, optical coherence tomography (OCT) acquires cross-sectional images of biological tissues by measuring back-scattered light from different depths [1,2]. The contrast of OCT originates from the in-sample spatial fluctuation of refractive index, which ranges from ∼1.3 to ∼1.5 for most biological tissues [3]. Such back-scattering intensity-based contrast information can be insufficient for tissue characterization and disease diagnosis, especially when a longer center-wavelength (e.g. 1300 nm) light source is used [4,5]. On the other hand, optical attenuation properties of biological tissues, including the absorption and scattering coefficients, can be extracted from OCT data as valuable complementary parameters for tissue characterization, abnormality detection, and disease diagnosis [4,6–9]. Conventional approach to quantifying tissue attenuation involves typically modeling the tissue as a homogeneous slab and fitting part of OCT A-line signal to a single-exponential decay curve [10–14]. Such simplified treatment smears the depth-resolved ability of OCT and disregards the depth-dependent attenuation coefficient of practical heterogeneous biological tissues. Recently, a depth-resolved attenuation estimation algorithm has gained popularity [15–17]. Under the single-scattering framework and assuming that almost all light is attenuated within the recorded imaging depth, this algorithm estimates local attenuation coefficients with axial resolution and much improved numerical robustness than the straightforward piecewise fitting method [15]. The model in discrete domain is transcribed below:

When applying this algorithm to practical OCT images of other biological tissues, however, extra care needs to be taken to avoid potential pitfalls. First, the algorithm doesn’t take into account noise floor in practical OCT A-lines, which becomes dominant at deeper imaging depth. Such noise floor, stemming from multiple-scattered photons [18,19] and system noise [2], violates the single-scattering assumption and compromises the overall attenuation estimation accuracy [20]. Additionally, we find that such excessive noise floor leads to considerable underestimation bias (as presented shortly). One intuitive solution is to completely exclude the noise floor-dominant portion of an A-line (from some specific depth beyond) from analysis; the resultant truncated A-line profile, although complying better with the single-scattering model, breaks the second assumption that almost all light gets attenuated within the imaging depth. This in turn leads to substantial overestimation bias, which grows significantly near the end of the depth range of interest [15,20,21]

In light of these practical pitfalls, Liu et al. [20] has recently proposed an improved algorithm which mitigates the growing error near the depth limit of analysis by adding the discarded sum of undetected signal beyond the boundary back to the denominator of Eq. (1), thereby yielding

Herein we systematically analyze the underling roots and corresponding consequences of under- and over- estimation issues that afflict the existing method, and then propose an alternative, robust depth-resolved attenuation coefficient calculation algorithm which rectifies the underestimation at shallower depths and overestimation at deeper locations for an incompletely decayed OCT A-line. The efficacy of this new algorithm is first validated via simulated numerical phantoms, and then further verified using practical data from both single- and multi-layer silicone-TiO₂ phantoms and ventral tongue leukoplakia tissues. Comparison with existing methods proposed by Vermeer et al. [15] and Liu et al. [20] is also presented.

2. Attenuation coefficient estimation model

2.1 Estimation bias of the contemporary method

The continuous form of $\tilde{\mu }$ in Eq. (1) can be written as:

Equation (4) above indicates that the presence of noise floor in practical OCT A-line profiles can cause substantial underestimation of attenuation coefficients in the useful imaging range (i.e., for depth $z < F$). The longer the noise floor portion, the more significant the underestimation. To verify the phenomenon, we synthesized a hybrid A-line profile $A$(z) by splicing a noise-free homogenous decay portion with attenuation coefficient of 3 mm⁻¹ (ranging from 0.5 mm to 1.7 mm) with an about 2.2-mm-long noise floor extracted from real-world OCT data acquired from a scattering phantom using a swept source OCT system (Fig. 1(a)). Indeed, the $\tilde{\mu }(z )$ profile estimated using Eq. (1) (blue curve, Fig. 1(b)) falls significantly below the theoretical value (black curve, Fig. 1(b)) in the exponential decay range of interest (from around 1.0 mm to 1.7 mm here). With the noise floor cut to half of the previous length (i.e., ∼1.1 mm long, Fig. 1(c)), as expected, the overall underestimation error gets less severe (blue curve, Fig. 1(d)), as manifested by the elevated attenuation coefficient value at transition depth $\tilde{\mu }({z = 1.7\, \textrm{mm}} )$ = ∼0.75 mm⁻¹ and the delayed onset of 3dB underestimation depth ${z_{3\textrm{dB}}}$ = ∼1.50 mm in Fig. 1(d), compared with corresponding values in Fig. 1(b) ($\tilde{\mu }({z = 1.7\, \textrm{mm}} )$ = ∼0.4 mm⁻¹ and ${z_{3\textrm{dB}}}$ = ∼1.35 mm, respectively).

 

Fig. 1. Estimation bias of Vermeer’s method on OCT data with noise floor. a, c, e. Synthesized OCT A-line data with the same length of effective signals (red dashed line), but different noise floor lengths, corresponding to the cases of long (a), short (c), and no (e) noise floor, respectively. g. Part of e resulting from overcutting of the noise floor. b, d, f, h. Corresponding depth-dependent attenuation profiles $\tilde{\mu }(z )$ obtained by Vermeer’s algorithm (blue curve) compared against ground truth profiles ${\mu _{\textrm{Ideal}}}(z)$ of matching lengths (black curve).

Download Full Size | PDF

To avoid such underestimation bias, one straightforward option is to completely exclude the noise floor portion from attenuation calculation (Fig. 1(e)), but in this way the resultant attenuation profile suffers in turn from the overestimation bias which grows with depth (Fig. 1(f)). Especially, if the noise floor gets overly truncated (Fig. 1(g)), the overestimation bias becomes more prominent (Fig. 1(h)). Such growing overestimation error near the end of the analyzing depth is also prominent in Fig. 1(b) and Fig. 1(d), as well as noticed in the original paper [15]. To gain further insight into this universal overestimation bias, we take a closer look at the approximation made in Eq. (3) and rewrite the equation as

As the underestimation bias can be eliminated by entirely excluding the noise floor portion from calculation (as demonstrated in Figs. 1(e)–1(h)), correcting the overestimation bias for general incompletely decayed A-lines is therefore critical for improving the overall estimation accuracy throughout the remaining single scattering-dominant depth range. In light of this, we propose a new estimation strategy to eliminate the progressive overestimation error near the calculation depth limit, as detailed in the following section.

2.2 Overestimation-free depth-resolved attenuation estimation

A spatially coherent beam propagating through a turbid medium gets attenuated due to both scattering and absorption. Considering a depth-resolved attenuation coefficient profile $\mu (z )$, the differential attenuation equation, according to Lambert-Beer’s law, can be described as:

3. Methods

3.1 Phantom fabrication and ex vivo imaging procedure

To experimentally verify the accuracy of the proposed depth-resolved algorithm, silicone phantoms with different concentrations of TiO₂ particles were fabricated following a published protocol [24] to mimic single- and multi-layer tissues. First, a ∼10-mm-thick single-layer homogeneous phantom was formed from a mixture of silicone and 0.2 w% concentration of TiO₂. Then to build a multi-layer phantom, four thin homogenous phantoms with 0.05 w%, 0.1 w%, 0.2 w%, and 0.3 w% TiO₂ concentrations, respectively, were first made separately; the thickness of each thin phantom was controlled to be around 300-350 µm by slowly pouring the viscous, uncured silicone-TiO₂ mixture into a petri dish under the guidance of real-time OCT imaging (with the refractive index of silicon-TiO₂ assumed to be 1.40). The surface of the phantom became smooth automatically owning to its fluidity. After leaving the petri dishes still on a level surface to cure for 24 hours, these thin phantom layers were then carefully harvested and stacked together; care was taken to make sure no inter-layer air gap exists.

Tissue imaging study was conducted in accordance with the protocol approved by the Ethics Committee of Tianjin Stomatological Hospital. A healthy volunteer and a patient with leukoplakia were enrolled for in vivo and ex vivo imaging studies. They were informed of experiment contents and consents were obtained prior to each experiment. OCT images of leukoplakia were collected ex vivo by imaging tissues harvested from the patient’s ventral tongue, while OCT images of normal ventral tongue were acquired in vivo from the healthy volunteer at a location corresponding to the leukoplakia position of the patient.

3.2 System setup and axial PSF calibration

All imaging experiments were conducted on a home-made swept-source OCT system with a measured axial resolution of ∼15 µm and a lateral resolution of ∼17 µm (both in air). The system consists of a 100-kHz sweep source laser (HSL-20-100-B, Santec) centering at 1310 nm with a tuning range of ∼87 nm, a Mach-Zehnder interferometer and a balanced photodetector (1817-FC, New Focus). Sensitivity roll off is an important factor in attenuation estimation. Thanks to the long coherence length of 20.3 mm, roll-off within 2 mm depth is almost unchanged, and is measured to be −2 dB over the entire imaging depth of 5.7 mm. Based on the fact that tissue imaging depth seldom exceeds 2 mm, the effect of roll-off is neglected here. A well-developed sensitivity roll-off expression can be used if necessity [25]. In practice, the scaling factor h in Eq. (7) is depth-dependent due to focusing effects, which necessitates experimental calibration. A classical and simplified axial PSF model is adopted here [26]:

4. Results

4.1 Digital phantom simulation

To validate the accuracy of our algorithm, we started with five single-layer, homogenous numerical phantoms with distinct attenuation coefficients. The imaging depth was set to 2.5 mm, with the phantom surface situated 0.2 mm deep, and the axial pixel size (i.e., resolution cell) set to 2.5 µm. For demonstration, the numerical simulations are free from noise floor and confocal effect. Thus, underestimation is not an issue here. The simulated OCT A-line profiles (magnitude squared) constructed according to Eq. (7) is plotted in Fig. 2(a). The attenuation coefficient profiles $\tilde{\mu }[n ]$ calculated using Eq. (1) and $\mu [n ]$ given by our algorithm (i.e., using Eq. (10)) are compared in Fig. 2(b). While matching the ground truth in the beginning, all $\tilde{\mu }[z ]$ profiles (solid curves, Fig. 2(b)) soar up monotonically with depth; also noticeable is that, for less scattering digital phantoms, the overestimation bias is generally more pronounced and sets in from a shallower depth, which conforms our previous analysis following Eq. (5). In contrast, the new method herein proposed compensates this error source by including a rectification factor $1 - {e^{ - 2\delta \cdot \mathop \sum \nolimits_{n + 1}^N \mu [m ]}}$ in Eq. (10). The resultant attenuation profiles (dashed lines, Fig. 2(b)), as expected, match the theoretical values well over the entire depth range. We also applied the method proposed in [20] on the same noiseless data, and obtained very similar attenuation profiles as our algorithm (data not shown).

 

Fig. 2. Numerical simulation results of single- and multiple-layer digital phantoms. a. Simulated OCT signals (in logarithmic scale) of five different single-layer phantoms, with the corresponding attenuation coefficients indicated in matching color. b. Attenuation profiles recovered from OCT signals in a with both Vermeer’s method (solid lines) and the new algorithm herein proposed (dashed lines). c. Simulated OCT signal (in logarithmic scale) from a noisy 5-layer heterogeneous phantom with attenuation coefficients of 0.5 mm⁻¹, 1 mm⁻¹, 3 mm⁻¹, 2 mm⁻¹, and 0.5 mm⁻¹ (from top to bottom), respectively. d. Attenuation profiles recovered from the OCT signal in c with Vermeer’s method (blue curve) and the new algorithm herein proposed (yellow curve).

Download Full Size | PDF

To further validate the proposed method, we also simulated a 5-layer numerical phantom. Staring from 0.2 mm deep, five 0.516-mm-thick layers of different attenuation coefficients were concatenated (Fig. 2(c)). Gaussian noise was then added with a signal-to-noise ratio (SNR) of ∼26 dB, which was deliberately set poorer than normal OCT, in order to examine the robustness of both estimation algorithms. As shown in Fig. 2(d), both methods are robust against noise. However, the attenuation profile $\tilde{\mu }(z )$ given by the previous method [15] is afflicted with prominent soaring tails in the bottom two layers (blue curve in Fig. 2(d), while the $\mu (z )$ profile given by the proposed method (yellow curve in Fig. 2(d)) matches well the ground truth (orange curve in Fig. 2(d)). The attenuation profile calculated by using the Ref. [20] method also matches the ground truth well (data not shown), indicating that that method achieves similar performance as ours on synthetic data with no noise floor.

4.2 Phantom experiments

A single-layer quasi-infinite phantom (∼10 mm thick) was imaged to experimentally assess the performance of the algorithm. After correcting the axial PSF as described in Section 3.2, and realigning A-lines to get a flat phantom surface for the convenience of later comparison, we first determined the magnitude of noise-floor based on the bottom ∼2.5 mm portion of the averaged A-line (purple double arrow, Fig. 3(b)), and then determined the truncation depth as the location where the decaying OCT magnitude is 5dB higher than that of noise-floor (∼2.3 mm here; blue arrow, Fig. 3(b)) [20]. Then within the remaining single scattering-dominant, high-fidelity shallower depth range, a best estimation of the ${\mu _{\textrm{Ideal}}}$ was determined by fitting the logarithm of the A-line magnitude to a linear model log(A(z)) = -µz + ɛ (red curve in Fig. 3(b)) yielding ${\mu _{\textrm{Ideal}}}$ ∼ 2.13 mm⁻¹ with a coefficient of determination of 0.98. This value serves as the “ground-truth” value against which other estimation results will be compared.

 

Fig. 3. OCT intensity and attenuation coefficient images of a single-layer silicone-TiO₂ phantom. a. Normalized OCT image. b. Averaged OCT magnitude of all A-lines (black curve) for noise floor determination. The noise floor portion is demarcated by purple double arrow, and the transition depth is indicated by a blue arrow. The effective signal portion used for global fitting of ${\mu _{\textrm{Ideal}}}$ is overlaid with the resultant single-exponential decay curve (red curve). c. $\tilde{\mu }$-image calculated with noise floor (w/ NF) using Eq. (1) after axial PSF correction and surface alignment. d. Mean attenuation profiles averaged over all A-lines of c, e, f and g. The black curve is the globally fitted value. Attenuation profile of c from 2.3 mm to 5.7 mm is plotted separately in the right panel. e, f, g. $\tilde{\mu }$-, $\hat{\mu }$- and $\mu $-images calculated from truncated A-lines without noise floor. Bright (dark) striped artefacts representing over- (under-) estimation bias are indicated by red (yellow) arrowheads in $\hat{\mu }$- and $\mu $-images. All OCT images shown are 5 mm width in the lateral direction.

Download Full Size | PDF

First, without removing the noise floor, $\tilde{\mu }$-image was calculated based on Eq. (1) using A-line data throughout the entire imaging depth of 5.7 mm (Fig. 3(c)). The resultant attenuation image exhibits the characteristic underestimation bias in the relatively shallow depth range (about 1.7 to 2.3 mm) as we have seen in Fig. 1. The averaged profile over all such-calculated A-lines is also shown in Fig. 3(d) (orange curve), where its deviation from the globally fitted ${\mu _{\textrm{Ideal}}}$ value (black curve, Fig. 3(d)) is prominent. The $\tilde{\mu }$-image also exhibits rapidly growing overestimation error near the depth limit (orange profile in the right extension panel of Fig. 3(d)). Nevertheless, the strongly scattering (and thus more attenuating), randomly distributed clumps of particles are more distinguishable in the $\tilde{\mu }$-image than in the intensity image, demonstrating the enhanced imaging contrast afforded by attenuation images.

Then, with data in this truncated depth range, the $\tilde{\mu }$-image given by Ref. [15], $\hat{\mu }$-image given by Ref. [20] and $\mu $-image given by our method are juxtaposed in Figs. 3(e)–3(g), and the respective average attenuation profiles across all A-lines are compared in Fig. 3(d). First, we notice that there’s no underestimation for any attenuation image or averaged attenuation profile, underpinning the effectiveness and necessity of removing the noise floor from depth-resolved attenuation calculation. Second, the growing overshooting bias of $\tilde{\mu }(z )$ estimation near the end of depth range is manifest in both $\tilde{\mu }$-image (Fig. 3(e)) and the average attenuation profile (blue curve, Fig. 3(d)). This again corroborates our previous analysis about the intrinsic overestimation tendency of the $\tilde{\mu }(z )$ estimation algorithm as revealed in Fig. 1(f), where cutting off the noise floor introduces overestimation to the informative signal range. In comparison, mean $\mu (z )$ profile estimated with our algorithm matches the fitted value (${\mu _{\textrm{Ideal}}}$ ∼ 2.13 mm⁻¹) across the entire selected depth range of analysis without systematic under- or over-estimation bias. Mean $\hat{\mu }(\textrm{z} )$ profile agrees with the ground truth and mean $\mu (\textrm{z} )$ profile well over most depth range except near the end, where it exhibits considerable overshooting bias (beyond 2.25 mm, Fig. 3(d)). Corresponding to this depth range, the bottom part of the $\hat{\mu }$-image, however, contains both extra bright overestimated regions (red arrowheads, Fig. 3(f)) and dark striped underestimated shadows (yellow arrowheads, Fig. 3(f)), implying that the overall overshooting error observed in mean $\hat{\mu }(\textrm{z} )$ profile is simply an ensemble effect. Down to each individual A-line, the attenuation profile could suffer from either over- or under-estimation bias in this depth range. In contrast, the $\mu $-image obtained with our algorithm contains almost no bright and fewer dark stripes (yellow arrowhead, Fig. 3(g)). Given that both algorithms start with the same fitted boundary attenuation value $\mu [N ]$ for each A-line, and that the fitted $\mu [N ]$ value is error-prone due to the noisy nature of OCT data (especially at such deep locations), this observation suggests that our algorithm is more robust to initial estimation error (see Discussion for details).

To evaluate the performance of the proposed algorithm in heterogeneous samples, a four-layer scattering phantom was fabricated. From top to bottom, the TiO₂ concentration in each layer is 0.05 w%, 0.1 w%, 0.3 w% and 0.2 w%, respectively, with the intensity image shown in Fig. 4(a). Again, the $\tilde{\mu }$-image calculated using the entire A-line data without cutting off the noise floor part (Fig. 4(b)) suffers eye-catching overestimation error near the lower boundary. The corresponding mean attenuation profile across all A-lines throughout the imaging field is shown in Fig. 4(f) (orange curve) and compared with the layer-wise fitted values (black line) used as “ground truth”. One can observe that $\tilde{\mu }$ values in layers III and IV are severely underestimated. Note also that theoretically, the mean attenuation values are supposed to be proportional to the respective TiO₂ particle concentrations, and stay mostly uniform within each layer. In practice, the fitted values fall below expectations in layers III-IV with higher TiO₂ concentrations, which results probably from multiple scattering [27] and the inherent difficulty of fully mixing and evenly distributing larger amount of TiO₂ particles within silicone solution.

 

Fig. 4. OCT image and the corresponding attenuation images of a multi-layer silicone-TiO₂ phantom. a. Normalized OCT image. b. $\tilde{\mu }$-image calculated over the entire imaging depth of 5.7 mm with noise floor (w/ NF). c, d, e. $\tilde{\mu }$-, $\hat{\mu }$- and $\mu $-images calculated without noise floor (w/o NF). Bright (dark) striped artefacts representing over- (under-) estimation bias are indicated by red (yellow) arrowheads in the $\hat{\mu }$-image. f. Comparison of mean attenuation profiles averaged over all constituent A-lines of subfigures b, c, d and e, respectively. g. Quantitative comparison of the average value and standard error of mean attenuation profiles shown in f over all depths within each layer. h. A 3D volumetric $\mu $-image demonstrating the consist estimation accuracy of our algorithm. All OCT images shown are 5 mm width in the lateral direction.

Download Full Size | PDF

As the lower boundary of this heterogeneous phantom (prepared in a Petri dish) is obvious in the intensity image, we simply treated all depth beyond the boundary as noise floor and excluded it from further calculation. As expected, attenuation images calculated with noise floor removed (Figs. 4(c)–4(e)) are free of such systematic underestimation bias. The respective mean attenuation profiles of all A-lines are also compared in Fig. 4(f). Two large peaks at the top and bottom of the phantom reveal strong reflection (thereby manifested as strong attenuation) due to the large variation of refractive index between air and phantom. The mean $\tilde{\mu }$(z) profile calculated without noise floor (green curve in Fig. 4(f)) is again polluted by significant overshooting bias near the bottom of layer IV (i.e., the depth limit of analysis here), while the mean $\hat{\mu }(z )$ and $\mu (z )$ profiles (blue and red curves, Fig. 4(f)) generally follow the fitted value much closer expect that $\hat{\mu }(z )$ values in layer IV are slightly overestimated compared with the fitted value. Correspondingly, scattering particles in layer IV of $\hat{\mu }$- and $\mu $-images are much cleaner and more discernible (Figs. 4(d)–4(e)) than those of $\tilde{\mu }$-image (Fig. 4(c)).

We then turn our attention to the uniformity of these attenuation profiles within each layer. The average and standard error of mean attenuation coefficients within each layer obtained by different calculation methods are compared in Fig. 4(g). All methods yield more uniform attenuation profiles within the first two, relatively less scattering layers, as one can also tell from attenuation images and profiles (Figs. 4(c)–4(f)). Within the bottom two layers: 1) both mean $\tilde{\mu }(z )$ profiles (with or without noise floor) demonstrate significant intra-layer variation, besides deviating greatly from the respective fitted values; 2) mean $\hat{\mu }(z )$ profile demonstrates larger variability than mean $\mu (z )$ profile obtained by our algorithm. The increased non-uniformity of $\hat{\mu }$(z) values is also manifested by conspicuous, randomly distributed bright or dark striped artefacts in layer IV of the $\hat{\mu }$-image (red and yellow arrowheads, Fig. 4(d)). The less occurrence of such over- and under-shooting bias in the $\mu $-image (Fig. 4(e)) suggests again that our algorithm is more robust to initial estimation bias of $\mu [N ]$ (see Discussion for details).

Finally, a three-dimensional volumetric visualization of sequentially stacked 500 $\mu $-images obtained by the proposed algorithm is shown in Fig. 4(h). The layered cake appearance demonstrates the consistent accuracy of the proposed algorithm over the entire field of view.

4.3 Human ventral tongue leukoplakia imaging

With the efficacy of our algorithm validated via both numerical simulation and phantom experiments, we applied it to clinical OCT data. Figure 5(a) is an OCT image of normal ventral tongue mucosa collected from of a 24-year-old healthy volunteer in vivo. Figure 5(b) is an ex vivo OCT image of a leukoplakia tissue section harvested from the ventral tongue mucosa of a 63-year-old male patient. Noise floor of both images has already been truncated according to the same aforementioned 5dB criterion. Corresponding attenuation images obtained by using the method developed in this article are shown in Figs. 5(c) and 5(d), respecively, and corresponding $\hat{\mu }$-images constructed by the method proposed in Ref. [20] in Fig. 5(e) and 5(f).

 

Fig. 5. OCT images and the corresponding attenuation images of normal ventral tongue mucosa and a resected leukoplakia tissue sample. a. Normalized OCT image of normal ventral tongue mucosa. b. Normalized OCT image of resected leukoplakia tissue from ventral tongue mucosa. c, d. Corresponding $\mu $-images calculated with our algorithm of a and b respectively. The upper and lower boundaries of the epithelium layer (EP) are marked by yellow dash curves, while the dysplasia-invaded portion of EP is delineated with a red dash curve in d. e, f. $\hat{\mu }$-images calculated with the algorithm by Liu et al. [20] of a and b respectively. Bright (dark) striped artefacts representing over- (under-) estimation bias are indicated by red (yellow) arrowheads. All OCT images shown are 5 mm width in the lateral direction.

Download Full Size | PDF

We first delineate the epithelium (EP) layer based on both intensity and attenuation images. While the lower boundary of EP layer is more clear in the attenuation images (Figs. 5(c) and 5(d)), its upper boundary is hardly discernible there, but has to be inferred from the corresponding intensity images (Figs. 5(a) and 5(b)), reflecting the benefits gained from leveraging these two complementary contrast mechanisms simultaneously. Compared with normal epithelium of the healthy volunteer (Figs. 5(a) and 5(c)), the lower portion of EP layer that is invaded by dysplasia exhibits stronger scattering (Fig. 5(d)), while the more superficial, normal-looking part of EP layer remains low scattering as the normal EP layer. The difference between the dysplasia region and the remaining uninvaded EP layer is more distinct in the attenuation image (delinated by a red dash curve, Fig. 5(d)) than in the intensity image (Fig. 5(b)), underscoring the resolving power of our estimation algorithm. Apart from increasing the attenuation coefficients, dysplasia invasion also thickens the EP layer here, resulting in a peak thickness of ∼650 µm (double arrow, Fig. 5(d)), compared with a maximal thickness of ∼400 µm in the normal EP layer (double arrow, Fig. 5(c)).

While comparable structural and attenuation information can also be extracted from corresponding $\hat{\mu }$-images (Figs. 5(e)–5(f)), as we have seen in the cases of TiO₂ phantoms, the lower part of both $\hat{\mu }$-images are afflicted with irregularly distributed extra dark (yellow arrowhead, Figs. 5(e)–5(f)) and extra bright (red arrowheads, Figs. 5(e)–5(f)) striped artifacts near the lower boundary, suggesting that the accuracy of $\hat{\mu }(z )$ estimation is sensitive to the fitted boundary attenuation coefficients µ[N], while the method proposed in this article is more robust (see Discussion for details).

5. Discussion

Both the algorithm by Liu et al. [20] and our algorithm require the lower boundary attenuation coefficient $\mu $[N] as an input parameter (or boundary conditions). With the incident beam decays monotonically inside a scattering medium like biological tissue, the signal-to-noise ratio of OCT signal from deeper locations is generally inferior, and therefore the estimation of boundary $\mu [N ]$ is subject to over- or under-shooting error. How such initial error at the boundary affects the overall accuracy of attenuation estimation determines the robustness and subsequently the applicability of the algorithm to practical biological imaging. From estimation results above, we have seen that our algorithm turns out more robust than the method by Liu et al. [20], as manifested by the better uniformity (especially at the bottom portion of phantom images) of $\mu $-images. To gain further insight into this observation, we attempt a theoretical analysis on how initial error propagates for both algorithms.

Assume a noise-free OCT A-line profile complying with the single-scattering model and the initial estimation of $\mu [N ]$ deviates from its true value, denoted by ${\mu^{\ast}}[N ]$, by a ratio (or percentage) of p, i.e. $\mu [N ]= ({1 + p} ){\mu^{\ast}}[N ]$. Note that the two algorithms converge under this noise-free setup, so the true attenuation values ${\mu^{\ast}}[n ]$ for $n = 1,\, \ldots ,\, N$ are the same for both algorithms; in the following we focus on how $\hat{\mu }[n ]- {\mu^{\ast}}[n ]$ given by Liu et al.’s algorithm and $\mu [n ]- {\mu^{\ast}}[n ]$ given by the reported algorithm evolve over depth.

First, following Eq. (2), the error in Liu et al.’s algorithm ends up

Then for our algorithm, the error evolves as

The reported attenuation estimation procedure can be further optimized in several aspects. For example, while the model for confocal parameters estimation we adopted here is widely used and yields satisfactory results, new automatic methods that determine confocal parameters from two OCT B-frames at different focal planes or incident angles can further improve the estimation accuracy [28,29]. Second, similar to existing methods, our algorithm is essentially derived from the single-backscattering model of OCT signal [19]. Extending depth-resolved attenuation calculation to account for multi-scattered (incoherent) component in general OCT signal warrants further investigation [12,27]. Further, as in most OCT-based attenuation profiling algorithms, the physical thickness of scattering tissues or phantoms was converted from optical path length by assuming a constant refractive index, therefore the attenuation coefficients were evaluated essentially against the optical path length. Integrating our algorithm with novel OCT imaging configurations that can retrieve refractive index information [30–32] and thus extract physical depth-based attenuation information can potentially yield new insights into quantitative tissue scattering characterization for biomedical research and clinical applications. Finally, it is noteworthy that attenuation images obtained with our algorithm can help highlight boundaries between tissue (sub)layers, thus offering a valuable feature dimension for automatic tissue layer segmentation in OCT images of various tissues [33–37].

6. Conclusion

In summary, we proposed an accurate and robust depth-resolved attenuation estimation algorithm that remedies the inherent under- and over-estimation bias of existing methods, thereby enabling depth-dependent attenuation calculation for general OCT A-line data containing incomplete decay and noise floor, and allowing consistent attenuation estimation accuracy for various scattering samples over the effective OCT imaging depth. Compared with other algorithms, our method demonstrates improved numerical robustness, as experimentally observed and theoretically confirmed. Results from both numerical simulations, phantom experiments and human tissues verified the excellent performance of this new algorithm and its universal applicability.