Robust, accurate depth-resolved attenuation characterization in optical coherence tomography

: Depth-resolved optical attenuation coeﬃcient is a valuable tissue parameter that complements the intensity-based structural information in optical coherent tomography (OCT) imaging. Herein we systematically analyzed the under- and over-estimation bias of existing depth-resolved methods when applied to real biological tissues, and then proposed a new algorithm that remedies these issues and accommodates general OCT data that contain incomplete decay and noise ﬂoor, thereby aﬀording consistent estimation accuracy for practical biological samples of diﬀerent scattering properties. Compared with other algorithms, our method demonstrates remarkably improved estimation accuracy and numerical robustness, as validated via numerical simulations and on experimental OCT data obtained from both silicone-TiO 2 phantoms and human ventral tongue leukoplakia samples.


Introduction
In direct analogy to ultrasound imaging, optical coherence tomography (OCT) acquires crosssectional 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][7][8][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][11][12][13][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][16][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: whereμ denotes the depth-resolved attenuation coefficient profile given by Ref. [15], A[n] the real-valued magnitude of a complex OCT A-line profile, δ the coherent length of the low-coherence beam (also the axial pixel size here), and N the number of resolvable pixels within each A-line. Note that the A-line magnitude A [n] in OCT imaging is proportional to the E-field magnitude of light returned from the sample arm, and therefore A 2 [n] denotes the depth-dependent light intensity in the scattering sample. 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 yieldinĝ Here the boundary attenuation efficient µ[N] is obtained by fitting a short piece of OCT A-line profile near the effective depth limit to the classical exponential decay model. This essentially substitutes the unwanted noise floor by a synthetic single-scattering tail, thereby reinstating the assumption of almost complete attenuation of Vermeer's algorithm [15]. Nevertheless, we found this method relies heavily on the fitting accuracy of the boundary µ[N], which in practice could lead to conspicuous and randomly distributed bright or dark striped artifacts in the resultant attenuation image. More recently, Amaral et al. [21] proposed a proof-of-concept general model that generalizes Vermeer's model and promises improved attenuation estimation accuracy; however, this method requires physically measuring the transmittance through the selected tissue slab, which is very challenging in real-world in vivo imaging circumstances. Herein we systematically analyze the underling roots and corresponding consequences of underand 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 2 phantoms and ventral tongue leukoplakia tissues. Comparison with existing methods proposed by Vermeer et al. [15] and Liu et al. [20] is also presented.

Estimation bias of the contemporary method
The continuous form ofμ in Eq. (1) can be written as: where D is the maximum imaging depth. The approximation made in Eq. (3) is rooted in the almost complete attenuation assumption which, however, hardly holds in real-world OCT imaging practice, even for highly scattering biological tissues (thus large µ values). In a typical OCT A-line profile, the informative signal decays rapidly with depth and finally transits into a trailing noise floor that extends until the depth limit. As a simplified model, denoting the depth at which such transition occurs by F (e.g., 5 dB above the noise floor intensity), we find that the attenuation coefficient given by Eq. (3) for shallower depths z<F turns out: Here we have introduced a hypothetical continuation of the A-line profile A(z) beyond the transit depth F, denoted by B(z) for z>F, which complies with the single scattering-based signal attenuation model, and therefore has smaller magnitude than the noise floor dominated by multiple-scattered photons or system noise [2,18,19]. 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 −1 (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μ(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μ(z = 1.7 mm) = ∼0.75 mm −1 and the delayed onset of 3dB underestimation depth z 3dB = ∼1.50 mm in Fig. 1(d), compared with corresponding values in Fig. 1(b) (μ(z = 1.7 mm) = ∼0.4 mm −1 and z 3dB = ∼1.35 mm, respectively).
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 At shallower depths (i.e. smaller z) where ∫ D z A 2 (y)dy ∫ ∞ D A 2 (y)dy, the estimatedμ(z) approximates the ground truth µ Ideal (z) well with negligible error. As z grows and approaches the calculation depth limit D, the integral ∫ D z A 2 (y)dy in the denominator decreases monotonically, and therefore the overestimation bias increases accordingly and peaks at depth D. The depth at which such rising-tail overestimation bias becomes prominent is controlled by the magnitude of ∫ ∞ D A 2 (y)dy. The less the light gets attenuated up to the depth limit D, the larger the integral ∫ ∞ D A 2 (y)dy becomes, and the earlier (and stronger) the overestimation bias emerges (and grows), as confirmed by simulation results in Fig. 1. Especially, even when the incident light attenuates almost completely within the imaging range, thereby the magnitude of ∫ ∞ D A 2 (y)dy being small, such rising-tail overestimation bias still emerges as z approaches the limiting depth D and ∫ D z A 2 (y)dy approaches zero. Therefore, such overestimation bias near the depth limit is inherent to the original algorithm [15] both physically and numerically. 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.

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 µ(z), the differential attenuation equation, according to Lambert-Beer's law, can be described as: where L(z) denotes the optical irradiance (W/cm 2 ) of the beam after traveling a distance z in the medium. For OCT imaging of most biological tissues with near-infrared (NIR) light source, the contribution of tissue absorption to µ(z) can be ignored [3], and therefore the attenuation of OCT signal results almost completely from tissue scattering. Assuming that a fixed fraction β out of the attenuated portion of incident irradiance near depth z is back-reflected towards the OCT detector, the back-scattered irradiance, denoted by the square of A-line profile A 2 (z), should comply with: where a scaling factor C denotes the conversion scale from beam irradiance to OCT signal, and L R the irradiance of the constant reference-arm beam. The exponential decay term in the first line of Eq. (7) represents light attenuation in the return path [19]. Note that the factor h is regarded depth-independent for now so that the model focuses purely on signal decay induced by sample attenuation; this corresponds to practical OCT A-line data with the focusing effect calibrated (see Methods). Integrating A-line magnitude squared from z to the imaging depth limit D, one gets: We thus obtain an implicit equation of the attenuation coefficient profile µ(z) as: The corresponding discrete version of Eq. (9) reduces to This equation can be solved recursively starting from the boundary value µ[N], which is in turn fitted from the last typically ∼100-200 µm (sample-dependent) segment of the informative depth range, after locating the tissue surface and cutting off the noise floor. To determine the fitting length, we average 10 neighboring A-lines, and then fit the last 100 µm of the average A-line to a log-linear model log A(z) = −µz + . If the resultant correlation coefficient R (defined as denotes the logarithm A-line values to fit, f i the corresponding fitted values, andȳ . = 1 M M i y i the mean of y i dataset) is smaller than 0.95, the fitting interval will be elongated in 20 µm increment until R achieves 0.95 or until the fitting interval reaches 200 µm (in which case we will accept the fitted value anyway to avoid deviating too far from the boundary). Note that in the context of simple linear regression, R is the square root of coefficient of determination R 2 which explains the extent to which the total variance of raw data is explained by the model [22,23]. In the following sections, the efficacy of the proposed algorithm will be validated on both numerical simulations and actual OCT data from scattering phantoms and clinical tissue experiments; comparison with existing algorithms reported in Ref. [15] and Ref. [20]will also be presented and discussed.

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 2 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 2 . 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 2 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 2 mixture into a petri dish under the guidance of real-time OCT imaging (with the refractive index of silicon-TiO 2 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.

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]: where z f denotes the position of beam waist and z R the Rayleigh length. Their values were estimated by sweeping a flat mirror placed inside a water chamber through the imaging range of the sample arm and fitting the back-reflection signal to Eq. (11) with α = 1 (i.e. the coherent case [4]). Then the axial PSF for scattering phantom imaging was estimated by setting α = 2 in Eq. (11) due to the loss of coherence. For proper calibration, all phantoms were imaged under the same sample-and reference-arm configuration as used in PSF measurement.

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μ[n] calculated using Eq. (1) and µ[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μ[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 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). 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μ(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 µ(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.

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 µ 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 µ Ideal ∼ 2.13 mm −1 with a coefficient of determination of 0.98. This value serves as the "ground-truth" value against which other estimation results will be compared.
First, without removing the noise floor,μ-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 µ Ideal value (black curve, Fig. 3(d)) is prominent. Theμ-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μ-image than in the intensity image, demonstrating the enhanced imaging contrast afforded by attenuation images.
Then, with data in this truncated depth range, theμ-image given by Ref. [15],μ-image given by Ref. [20] and µ-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μ(z) estimation near the end of depth range is manifest in bothμ-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μ(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 µ(z) profile estimated with our algorithm matches the fitted value (µ Ideal ∼ 2.13 mm −1 ) across the entire selected depth range of analysis without systematic under-or over-estimation bias. Meanμ(z) profile agrees with the ground truth and mean µ(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μ-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μ(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 µ-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 µ[N] for each A-line, and that the fitted µ[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 2 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μ-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μ 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 2 particle concentrations, and stay mostly uniform within each layer. In practice, the fitted values fall below expectations in layers III-IV with higher TiO 2 concentrations, which results probably from multiple scattering [27] and the inherent difficulty of fully mixing and evenly distributing larger amount of TiO 2 particles within silicone solution. 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μ(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μ(z) and µ(z) profiles (blue and red curves, Fig. 4(f)) generally follow the fitted value much closer expect that µ(z) values in layer IV are slightly overestimated compared with the fitted value. Correspondingly, scattering particles in layer IV ofμ-and µ-images are much cleaner and more discernible (Figs. 4(d)-4(e)) than those ofμ-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μ(z) profiles (with or without noise floor) demonstrate significant intra-layer variation, besides deviating greatly from the respective fitted values; 2) meanμ(z) profile demonstrates larger variability than mean µ(z) profile obtained by our algorithm. The increased non-uniformity ofμ(z) values is also manifested by conspicuous, randomly distributed bright or dark striped artefacts in layer IV of theμ-image (red and yellow arrowheads, Fig. 4(d)). The less occurrence of such over-and under-shooting bias in the µ-image (Fig. 4(e)) suggests again that our algorithm is more robust to initial estimation bias of µ[N] (see Discussion for details).
Finally, a three-dimensional volumetric visualization of sequentially stacked 500 µ-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.

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 correspondinĝ µ-images constructed by the method proposed in Ref. [20] in Fig. 5(e) and 5(f).
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μ-images (Figs. 5(e)-5(f)), as we have seen in the cases of TiO 2 phantoms, the lower part of bothμ-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μ(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).

Discussion
Both the algorithm by Liu et al. [20] and our algorithm require the lower boundary attenuation coefficient µ[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 µ[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 µ-images. To gain further insight into this observation, we attempt a theoretical analysis on how initial error propagates for both algorithms.
In reality the magnitude of 1/µ * [N] is on the same order of scattering mean free path, which is typically several hundred microns for the near-infrared wavelength region adopted in OCT imaging, while the pixel size of OCT image δ is typically several microns.
This implies that the initial error ratio p, be it over-or under-estimation, will persist over a considerable depth range immediately above the boundary. Since whether the initial boundary attenuation coefficient µ[N] is under-or over-estimated, i.e. the sign of p, is random, this explains the irregularly distributed bright or dark stripe artifacts appearing at deep locations (or layers) ofμ-images (Figs. [3][4][5]. Furthermore, as the calculation proceeds toward upper depths, with n decreasing and N n+1 A 2 [m] growing monotonically, the effect of initial error ratio p gradually diminishes, which matches the improved performance of this algorithm at shallower depths.
Note that the further approximation made in Eq. (16) above follows from applying Eq. (15) iteratively and noticing that the exponential index α[n] = 2δ · N n+1 µ * [m] grows slowly as n decrements. From Eq. (16) one can see that, while the initial error percentage p in µ[N] estimation also propagates to upper depths, its effect is suppressed progressively by an exponential decaying term. This explains the improved robustness of our algorithm as manifested by the improved uniformity of µ-images and mean µ(z) profiles shown above.
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][31][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][34][35][36][37].

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.

Funding
National Natural Science Foundation of China (61875092); Science and Technology Support Program of Tianjin (17YFZCSY00740); Scientific Research Foundation of Graduate School of Southeast University (3307031808D4).