1.

Introduction

Genesis and progression of disease are accompanied by morphological changes in tissues on length scales ranging from intracellular organelles to macroscopic tissue structures. These lead to changes in the spatial distribution of the complex refractive index, which in turn leads to changes in the absorption and scattering properties that can be measured using optical techniques. The main hypothesis underlying many applications of biophotonics is that, by measuring the optical properties, diagnosis or monitoring of tissue disease state or treatment is possible. Consequently, knowledge of the accuracy and precision of the methods to assess these optical properties is paramount.

The optical property accessible with optical coherence tomography (OCT) measurements is the attenuation coefficient, which describes the decay rate of the OCT signal with depth.^1,2 It is commonly extracted by non-linear least squares curve fitting (CF) of a single scattering-based model to the OCT signal.³ The main cause of imprecision in the determination of the attenuation coefficient is the inherent random fluctuation of the OCT signal, due to speckle and (shot) noise. Speckle is the voxel-to-voxel variation of OCT amplitude, caused by the spatial variation of the refractive index in tissue.^4,5 Randomly placed scatterers within the voxels will return scattered fields with random amplitude and phase, leading to intensity fluctuations at the detector.⁶ We have recently derived a simple expression for the minimal attainable precision with which the attenuation coefficient can be determined using CF based on the so-called Cramér–Rao (CR) lower bound:⁷

CF suffers from the drawback, as can be inferred from Eq. (1), that a finite sized AFR is necessary to achieve sufficient precision, which may preclude measurement of the attenuation coefficient of thin layers, such as in the retina or the arterial wall, or regions near the basal membrane.⁸ In recent years, the depth-resolved estimation (DRE) method has grown popular as an alternative to CF. Introduced to the OCT field by Vermeer et al.,⁹ it was inspired by earlier work on shadow removal in OCT¹⁰ and on ultrasound attenuation compensation,¹¹ once again demonstrating one of many conceptual similarities between both modalities. Practical improvements of the method were introduced by Liu et al.,¹² Smith et al.,¹³ and Dwork et al.¹⁴ The principal allure of the method is the (apparent) pixel-wise determination of the attenuation coefficient, which may circumvent the need for an AFR that extends far into depth.¹⁵ However, assessment of the accuracy and precision of DRE estimation is scarce. The aim of this study therefore is to determine the accuracy and precision of the DRE of the attenuation coefficient.

2.

Theory

Under the assumption of single backscattering from a homogeneous medium with stationary optical properties, the OCT signal versus depth $z$ is modeled as a single exponential decay combined with the confocal point spread function and sensitivity roll-off:²

We proceed to compute the definite integral ${\mathbb{I}}_{\mathcal{E}}(z)={\int }_z^{\mathcal{E}}I(z^{\prime })\mathrm{d}{z}^{\prime }$ of Eq. (3), which runs from the depth $z$ up to the end of the available (or used) data range $\mathcal{E}$, to estimate the attenuation coefficient as (Appendix A):

Compared to the original formulation by Vermeer, Eq. (4) contains a regularization term $I(\mathcal{E})/{\widehat{\mu }}_{\mathcal{E}}$ in the denominator to compensate for the finite data range.¹⁶ Here, ${\widehat{\mu }}_{\mathcal{E}}={\widehat{\mu }}_{\mathrm{OCT}}(\mathcal{E})$ is an independently obtained estimate for the attenuation coefficient at the end of the data range $\mathcal{E}$, which may be found, e.g., by CF or from transmission measurements.¹⁷ Due to speckle and noise fluctuations, the OCT signal is itself an inherently fluctuating quantity of which Eq. (2) represents the average. Inspection of Eq. (4) suggests that these fluctuations will be largely averaged out only in the denominator $\mathbb{I}(z)$ due to the integration. Indeed, as shown by Fiske et al.,¹⁸ the attenuation coefficient retrieved by the DRE follows the same statistical distribution as the OCT intensity $I(z)$, which is a Rayleigh distribution when the OCT signal is represented on amplitude basis, or an exponential distribution when the OCT signal is represented on an intensity basis as in this article:

The mean value $⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩$ can be obtained from a large set of estimations of ${\widehat{\mu }}_{\mathrm{OCT}}(z)$, in practice over some spatial range around $z$ and/or from several A-scans at the same position. Combining Eqs. (3) and (4), we theoretically obtain

Thus, $⟨{\widehat{\mu }}_{\mathrm{OCT}}(z)⟩\to {\mu }_{\mathrm{OCT}}$ at a location sufficiently far from $\mathcal{E}$, whereas $⟨{\widehat{\mu }}_{\mathrm{OCT}}(z)\to {\widehat{\mu }}_{\mathrm{OCT}}(\mathcal{E})⟩$ as $z$ approaches $\mathcal{E}$. Vermeer considered the effect of discretization of $I(z)$, i.e., each datapoint $I[i]$ corresponds to the integration of Eq. (3) over a finite pixel size $\mathrm{\Delta }$ around $z$. The exact, discretized version of Eq. (4) reads (Appendix B)

Often, a simplified version of Eq. (7) is used by linearization of the logarithmic and exponential terms [perhaps inspired by the closer visual resemblance to Eq. (4)]

Moreover, some authors further omit the (linearized) regularization term ${\mathcal{C}}_L=I[i_{\max }]/{\widehat{\mu }}_{\mathrm{OCT}}[i_{\max }]$ from Eq. (8). Use of these approximations is discouraged as they come with the penalty of reduced accuracy. The analysis in Appendix B reveals that Eq. (8) systematically overestimates the attenuation coefficient in the order of ${\mu }_{\mathrm{OCT}}^2\times \mathrm{\Delta }$.

We now seek the precision with which $⟨{\widehat{\mu }}_{\mathrm{OCT}}[i]⟩$ can be estimated with maximum likelihood (ML) from the Fisher information associated with the data using a CR analysis. Conceptually, Fisher information measures the amount of information that a dataset provides about the parameters of a model for the data. The CR lower bound, the inverse of Fisher information, measures the highest precision with which the parameters can be estimated using ML methods. In the case of CF of OCT intensity (or amplitude) values, the parameter of interest is the attenuation coefficient. In the present case, somewhat trivially, the parameter of interest is the mean of the set of ${\widehat{\mu }}_{\mathrm{OCT}}$ -estimations obtained by Eq. (7), which are distributed according to Eq. (5). The ML estimator of the mean of an exponential distribution is simply the arithmetic mean of the estimations. The Fisher information for $M$ independent estimations from an exponential distribution is ${Ϝ}_{\exp }=M/{⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩}^2$ so the CR lower bound, expressed as standard deviation becomes ${\sigma }_{{\mu }_{\mathrm{OCT},\mathrm{DRE}}}=⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩/\sqrt{M}$. Commonly, A-scans are pre-averaged prior to the application of the depth-resolved estimation. This changes the distribution of recovered attenuation coefficients to an approximately normal distributions when the number of over averages is $N\gtrsim 30$, with mean equal to $⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩$ and variance equal to ${\sigma }_{\mathcal{N}}^2={⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩}^2/N$. The corresponding Fisher information for the estimated attenuation coefficient based on $M$ independent measurements then becomes $Ϝ=MN/{⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩}^2$ and the CR lower bound, expressed as standard deviation is thus ${\sigma }_{{\mu }_{\mathrm{OCT},\mathrm{DRE}}}=⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩/\sqrt{MN}$. From a frequentist statistician’s point of view, this quantity represents the standard deviation of the normal distribution of attenuation coefficient values, which would be obtained if the DRE analysis were repeated many times. We can use this interpretation to calculate the decrease in precision under the influence of noise. The result of the lower bound of this analysis, which can be found in Appendix C, is

Comparing Eq. (9) to Eq. (1), we see that the precision in the DRE estimation method is directly proportional to the pixel-wise estimate of the attenuation coefficient, whereas it is independent of $⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩$ for CF. It also shows that a higher precision can be obtained using DRE compared to CF, when the AFR becomes smaller than $\sqrt{3c_{\mathcal{R}}}/{\widehat{\mu }}_{\mathrm{OCT}}$ or roughly two mean free paths.

3.

Methods

To validate the accuracy and precision derived in Eqs. (6) and (9), we performed numerical simulations based on OCT scans from a homogeneous medium. Details of our simulation procedure can be found in Ref. 7. Briefly, single A-scans with randomly varying amplitude were generated based on $A_{\mathrm{sim}}(z_i)=\sqrt{-\frac{4}{4-\pi }·{\sigma }_A^2(z_i)·\ln ({\xi }_i)}$, where ${\sigma }_A^2(z)=\frac{c_{\mathcal{R}}}{4}(I(z)+{⟨\zeta ⟩}^2)$ is the amplitude variance, $I(z)$ is given by Eq. (3), and ${\xi }_i$ is a uniformly distributed number between 0 and 1 drawn for each pixel. This procedure assures that the amplitudes $A_{\mathrm{sim}}(z)$ follow a Rayleigh distribution corresponding to fully developed speckle and the contribution of shot noise. $N=100$ of these A-scans are first squared, then averaged so the resulting averaged intensities at each depth position are normally distributed. Thereafter, we subtract the mean noise floor $⟨\zeta ⟩$.

Each squared, averaged A-scan is processed using Eq. (7) to estimate the attenuation coefficient. To demonstrate the spurious effect of linearization, we also analyzed the data using Eq. (8) with the regularization term ${\mathcal{C}}_L$ omitted. This procedure was repeated ${10}^4$ times to obtain a distribution of ${\widehat{\mu }}_{\mathrm{OCT}}[i]$ estimations at each depth position. Comparing the mean of this distribution to the input attenuation coefficient yields the accuracy of the method, and the precision is given by the distribution’s standard deviation.

We used comparable simulation parameters as previously⁷ reported for a direct comparison between the precisions of least squared fitting and the DRE method used in this article. In the simulations, $p_{\mathrm{NA}}$ was set to unity and an arbitrary scaling factor of ${2500}^2$ was included. Simulations were performed both with and without shot noise included. In the latter case, signal fluctuations are caused only by speckle. In the former case, the mean noise level was fixed at $\zeta =13.5$. Values of 2 and $5{\mathrm{mm}}^{-1}$ were used for the attenuation coefficient, which leads to a maximum SNR expressed in decibels of 60 and 64 dB, respectively. In all simulations, we used a value of ${\hat{\mu }}_{\mathrm{OCT}}(\mathcal{E})=5{\mathrm{mm}}^{-1}$ for the estimation of the attenuation coefficient at the end of the data range $\mathcal{E}$.

4.

Results

Figure 1(a) shows an example of simulated $N=100$ times averaged A-scans, including shot noise, obtained by the procedure outlined in the previous section. The arrow indicates the position $z_c=1.18\mathrm{mm}$, at which the intensity signal hits the noise floor with the condition $I[i]=⟨\zeta ⟩$ for the simulation with ${\mu }_{\mathrm{OCT}}=5{\mathrm{mm}}^{-1}$ as input parameter. For ${\mu }_{\mathrm{OCT}}=2{\mathrm{mm}}^{-1}$ the intensity signal does not reach the noise floor.

Fig. 1

(a) Example of simulated, $N=100$ times averaged A-scans including shot noise ($\zeta =13.5$) with ${\mu }_{\mathrm{OCT}}=2$ and $5{\mathrm{mm}}^{-1}$ using the model in Eq. (3). The backscatter efficiency within the NA $p_{\mathrm{NA}}$ was set to unity. The arrow indicates the position $z_c=1.18\mathrm{mm}$, where the OCT signal hits the noise floor. (b) The distribution of intensity values at $z=0.42\mathrm{mm}$ obtained from ${10}^4$ independent simulations for ${\mu }_{\mathrm{OCT}}=2{\mathrm{mm}}^{-1}$.

Figure 1(b) shows a histogram of intensity values at a depth of $z=0.42\mathrm{mm}$ using ${\mu }_{\mathrm{OCT}}=2{\mathrm{mm}}^{-1}$ obtained from of ${10}^4$ independent simulations. It shows that the averaged intensities, obtained by pre-averaging $N=100$ A-scans, are indeed to good approximation normally distributed.

The assessment of the accuracy of the DRE method is shown in Fig. 2, in the absence of shot noise, and in Fig. 3, in the presence of shot noise for attenuation coefficients of 2 and $5{\mathrm{mm}}^{-1}$. Both figures show the estimated attenuation coefficients versus depth using Eqs. (7) and (8) with ${\mathcal{C}}_L=0$ and the theoretical prediction of Eq. (6). Figure. 3(a) shows the DRE algorithm applied to a single averaged A-scan and demonstrates the remaining fluctuation in the estimations of ${\mu }_{\mathrm{OCT}}$; even after pre-averaging $N=100$ A-scans. The data shown in Fig. 3(b) is averaged over ${10}^4$ independent simulations and therefore permits closer comparison of theory and simulations. For the data shown in both figures, a value of ${\hat{\mu }}_{\mathcal{E}}=5{\mathrm{mm}}^{-1}$ is used in the regularization term $\mathcal{C}$ in Eq. (7). The results in Fig. 2(b) demonstrate that the estimated attenuation coefficient differs $<1\%$ from the true value up to a depth of 1.98 mm for ${\mu }_{\mathrm{OCT}}=2{\mathrm{mm}}^{-1}$ and, for ${\mu }_{\mathrm{OCT}}=5{\mathrm{mm}}^{-1}$, over the entire depth range. It approaches, in both cases, ${\widehat{\mu }}_{\mathcal{E}}$ near the end of the available data range. If, on the other hand, the linearized approximation Eq. (8) with ${\mathcal{C}}_L$ omitted is used, the attenuation coefficient has a fixed offset with respect to the true value and does therefore not stay within the 1% mark and, furthermore, tends to infinity at the end of the data range.

Fig. 2

Accuracy of ${\widehat{\mu }}_{\mathrm{OCT}}$ determination using the DRE method in absence of shot noise ($\zeta =0$). Attenuation coefficients of ${\mu }_{\mathrm{OCT}}=2$ and $5{\mathrm{mm}}^{-1}$ were used to simulate $N=100$ times averaged A-scans according to the model of Eq. (3) with $p_{\mathrm{NA}}$ set to unity. The number of sample points is $M=1$, as the reconstruction is done per pixel. (a) Data points, show ${\hat{\mu }}_{\mathrm{OCT}}$ determined using the exact, discretized Eq. (7) and using its linearized approximation Eq. (8), ${\mathcal{C}}_{\mathrm{L}}=0$. Solid lines show the theoretical prediction of $⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩$ using Eq. (6). The mean per pixel of ${10}^4$ independent, averaged A-scans is compared to the theoretical value in (b).

Fig. 3

Accuracy of ${\widehat{\mu }}_{\mathrm{OCT}}$ determination using the DRE method in the presence of shot noise ($\zeta =13.5$). Attenuation coefficients of ${\mu }_{\mathrm{OCT}}=2$ and $5{\mathrm{mm}}^{-1}$ were used to simulate $N=100$ times averaged A-scans according to the model of Eq. (3) with $p_{\mathrm{NA}}$ set to unity. The number of sample points is $M=1$, as the reconstruction is done per pixel. The arrow indicates the depth position $z_c=1.18\mathrm{mm}$ where the OCT signal with ${\mu }_{\mathrm{OCT}}=5{\mathrm{mm}}^{-1}$ hits the noise floor. (a) Data points, shows ${\widehat{\mu }}_{\mathrm{OCT}}$, calculated using the exact, discretized Eq. (7), and using its linearized approximation Eq. (8), ${\mathcal{C}}_{\mathrm{L}}=0$. Solid lines show the theoretical prediction of $⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩$ using Eq. (6). The mean per pixel of ${10}^4$ independent, averaged A-scans is compared with the theoretical value in (b).

The more realistic case, when noise is included, is depicted in Fig. 3. Analogously to Fig. 2, the attenuation coefficients, calculated from one single, averaged A-scan (a) and their means per pixel (b) from ${10}^4$ independent estimations is shown. Including noise into the calculation results into a strong fluctuation of the estimated attenuation coefficient as soon as the signal hits the noise floor at $z_C$. This trend is expected as a result of using the full depth range for the attenuation coefficient estimation. The sum term in the denominator of Eqs. (7) and (8) is padded by random noise values fluctuating around zero after ${\mathrm{z}}_{\mathrm{c}}$ (assuming the average noise floor is correctly subtracted). The contribution of noise to the sum term and therefore the effect of including the noise region into the calculation increases with depth, as is clearly seen more clearly in Fig. 3(b) for depths larger than $z_c$. However, including the noise area into the calculation does not show a significant effect the attenuation estimation in the depth region before the signal hits the noise floor and differs $<1\%$ of true attenuation value up to a depth of 1.18 mm for ${\mu }_{\mathrm{OCT}}=5{\mathrm{mm}}^{-1}$ and, for ${\mu }_{\mathrm{OCT}}=2{\mathrm{mm}}^{-1}$, over the entire depth range similar to the results in Fig. 2.

We proceed to compare the standard deviation ${\sigma }_{{\hat{\mu }}_{\mathrm{OCT}}}$ of the distribution of estimated attenuation coefficients at each depth with calculations based on Eq. (9) in Fig. 4(a) (without added shot nose) and Fig. 4(b) (with shot noise). It is shown in Fig. 4(a) that the results obtained using both the exact Eq. (7) and the approximation Eq. (8) with ${\mathcal{C}}_L=0$ are in good agreement with the predictions of Eq. (9) except very near to the end of the available data range. Since by Eq. (9), the precision is proportional to the mean estimated attenuation coefficient at each depth [e.g., the results of Figs. 3(b) and 4(b)], it is found that the precision is slightly higher when Eq. (8) is used instead of Eq. (7). When shot noise is present [Fig. 4(b)], there is good qualitative agreement between the simulations and the predictions of Eq. (9) with the largest deviations occurring at depths where the signal is close to the noise floor.

Fig. 4

Precision of ${\widehat{\mu }}_{\mathrm{OCT}}$ determination using the DRE method (a) in the absence ($\zeta =0$). and (b) in the presence of shot noise ($\zeta =13.5$). Attenuation coefficients of ${\mu }_{\mathrm{OCT}}=2$ and $5{\mathrm{mm}}^{-1}$ were used to simulate $N=100$ times averaged A-scans according to the model of Eq. (3) with $p_{\mathrm{NA}}$ set to unity. The number of sample points is $M=1$, as the reconstruction is done per pixel. The lower bound on the precision ${\sigma }_{{\hat{\mu }}_{\mathrm{OCT},\mathrm{DRE}}}$ (solid lines) was obtained from Eq. (9) and compared to the standard deviation ${\sigma }_{{\widehat{\mu }}_{\mathrm{OCT}}}$, (data points) obtained using DRE estimation based on the exact, discretized Eq. (7) and its linearized approximation [Eq. (8), ${\mathcal{C}}_{\mathrm{L}}=0$]. The arrow in (b) indicates the depth position $z_c=1.18\mathrm{mm}$ where the OCT signal with ${\mu }_{\mathrm{OCT}}=5{\mathrm{mm}}^{-1}$ hits the noise floor.

Finally, we compared our previously reported lower bounds for the CF method⁷ with the precision we derived in this article for the DRE method in the presence of shot noise ($\zeta =13.5$). Figure 5 shows the numerically obtained CR lower bound (the minimal precision for the CF approach) in dependence of $N$ for an AFR of $328\mu \mathrm{m}$ ($M=41$ points, $\mathrm{\Delta }=8\mu \mathrm{m}$) located well before $z_c$ such that the SNR in the AFR is $>20\mathrm{dB}$. The black dashed line represents the analytical CRLB calculated using Eq. (1) and overlaps with the numerically obtained curves, thus demonstrating the validity of Eq. (1) for low noise levels as well as the independence of the lower bound on the value of the attenuation coefficient itself. The precision of the DRE method calculated by Eq. (9) (also using $M=41$ points, with end-of-range values ${\hat{\mu }}_{\mathcal{E}}$ set to the true value of ${\mu }_{\mathrm{OCT}}$.) does show a dependence on the attenuation coefficient. All curves follow a $1/\sqrt{N}$ trend, whereas ${\sigma }_{{\widehat{\mu }}_{\mathrm{OCT}}}$ for the DRE method is smaller than for the CF method. Inspection of Eqs. (1) and (9) in the low-noise limit quickly reveals that the DRE method outperforms the CF method in precision, when ${\mu }_{\mathrm{OCT}}·|\mathrm{AFR}|\le \sqrt{3·c_R}\approx 1.8$.

Fig. 5

Comparison between the lower bound of the precision of the attenuation coefficient based on the DRE [Eq. (9)] and CF⁷ methods. The precision is shown in dependency on the number of averages $N$. Shot noise was included in the calculations ($\zeta =13.5$). A region of interest from $z_{\min }=0.04\mathrm{mm}$ to $z_{\max }=0.36\mathrm{mm}$, with a length of $328\mu \mathrm{m}$ and $M=41$ sample points (step size $\mathrm{\Delta }=8\mu \mathrm{m}$), was used in both methods. The end-of-range values ${\widehat{\mu }}_{\mathcal{E}}$ used in the DRE method were set to the true value of the attenuation coefficient ${\mu }_{\mathrm{OCT}}$. Note that the curves for the CF method overlap in this high-SNR limit.

5.

Discussion

Quantification of the attenuation coefficient requires thorough assessment of the accuracy and precision with which it can be estimated from OCT data. In recent years, the DRE method has emerged as an attractive alternative to the conventional approach of CF. We have derived expressions for the accuracy and precision of attenuation coefficient determined by the DRE method and validated those with numerical simulations.

The accuracy of the DRE method is given by Eq. (6). This equation includes a regularization term that sets the attenuation coefficient ${\hat{\mu }}_{\mathcal{E}}$ at the end of the available data range (either at the end of the A-scan, or the part of the A-scan that is included in the analysis). Omitting the regularization term essentially sets it value to infinity. This will result in an inaccurate attenuation estimation at the end of the data range, which can be overcome by choosing a region for attenuation reconstruction, which is far from the end of the data range. However, an accurate estimation can only be achieved when the attenuation coefficient is estimated using Eq. (7), which properly takes into account discretization of the OCT signal, as shown in Figs. 2 and 3. The approximate expression Eq. (8), which is often found in the literature, leads to a consistent overestimation of the attenuation coefficient in the order ${\mu }_{\mathrm{OCT}}^2\times \mathrm{\Delta }$, where $\mathrm{\Delta }$ is the pixel increment, and for that reason its use is highly discouraged. Next to inappropriate use of Eq. (8) [rather than Eq. (7)] that leads to a systemic offset, other factors that may contribute to a loss of accuracy are inadequate noise subtraction¹⁶ or incomplete compensation of the point spread function and roll-off function.¹⁹ Incorrect estimation of ${\hat{\mu }}_{\mathcal{E}}$ clearly leads to loss of accuracy at the end of the data range.

The expression for precision [Eq. (9), Appendix C] was derived under the premise that a number of $N$ A-scans are averaged prior to application of the DRE method. Averaging results in approximately normally distributed averaged intensity values, which is the underlying assumption of the derivation in Appendix C. In the absence of noise, this leads to a CR lower bound on the precision for the estimation of the mean attenuation coefficient as ${\sigma }_{\widehat{{\mu }_{\mathrm{OCT}}}}=⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩/\sqrt{MN}$. Contrary to our own advice in the previous paragraph, we continued to employ the linearized Eq. (8) with the purpose of arriving at a compact expression for the precision in the presence of noise. As evidenced from the results in Fig. 4(b), the resulting Eq. (9) is in good qualitative agreement with the simulation data and captures the effect of increasing SNR on the precision of the estimated attenuation coefficient. The key feature of Eq. (9) is that the precision depends on the mean estimated attenuation coefficient itself. Therefore, any loss in accuracy will directly lead to a loss of precision.

The dependence of precision on $⟨{\widehat{\mu }}_{\mathrm{OCT}}⟩$ is also the main difference between the precision obtained through NLLS CF with two free running parameters, as expressed in Eq. (1), which depends on the extent of the AFR.⁷ Comparing both methods (Fig. 5), we see that the DRE method can obtain a better precision when ${\mu }_{\mathrm{OCT}}|\mathrm{AFR}|\lesssim 1.8$ when the same number of A-scans ($N$) is pre-averaged, and the same number of independent data points ($M$) is included in the analysis. Both methods thus require spatial support to achieve sufficient precision (the required level of precision may well depend on the application). In this sense, the term “depth resolved estimation” is somewhat misleading because in practice information from some spatial region must be included for the analysis.

6.

Conclusion

In this article, we derived and validated the accuracy and precision of the depth resolved estimation method of the attenuation coefficient ${\mu }_{\mathrm{OCT}}$ in optical coherence tomography. We showed that a commonly used simplification of the method results in loss of accuracy in the order of $\mathrm{\Delta }\times {\mu }_{\mathrm{OCT}}^2$ where $\mathrm{\Delta }$ is the sampling resolution and is therefore not recommended as being used for OCT-attenuation reconstruction. Furthermore, we derived an analytical expression for the precision of ${\hat{\mu }}_{\mathrm{OCT}}$, which proportionally scales with its expectation value and inversely with the square root of the number of independent sample points included in the analysis. Lastly, we compared our outcome with the precision obtained using a CF procedure and provided an easy applicable rule of thumb to determine which method will have a better precision. Our theoretical framework gives valuable insight regarding accuracy and precision of parametric imaging based on a depth-resolved reconstruction of the attenuation coefficient and is, given its wide and easy-to-use applicability, an important advance toward design and improvement of standardized OCT-experiments, which are, e.g., used for tissue characterization in the clinic.

7.

Appendix A

We model the mean OCT signal intensity as function of depth in Eq. (3) using a single exponential decay function, assuming that the confocal point spread function and the sensitivity roll-off are fully compensated and that a constant mean noise floor is subtracted. We first compute the integral ${\mathbb{I}}_{\infty }(z)=-{\int }_z^{\infty }I(z^{\prime })\mathrm{d}{z}^{\prime }$ of Eq. (3), which yields

Taken together with Eq. (3), we can solve for the attenuation coefficient as

In practice, data are only available over a finite range, up to $z=\mathcal{E}$. We compute the definite integral ${\mathbb{I}}_{\mathcal{E}}(z)=-{\int }_z^{\mathcal{E}}I(z^{\prime })\mathrm{d}{z}^{\prime }$, which yields

Solving Eq. (12) for ${\mathbb{I}}_{\infty }(z)$, and substituting the result in Eq. (11) gives

Finally, Eq. (13) can be rearranged to solve for the attenuation coefficient giving the finite-range equivalent of Eq. (11)

When $z$ approaches the end of range $\mathcal{E}$, the finite integral in the denominator term goes to zero; therefore $2{\mathbb{I}}_{\mathcal{E}}(z)\to 0$ and the estimation of the attenuation coefficient in Eq. (14) will tend to infinity. A better strategy is to regularize the depth resolved estimation using an independently obtained estimate for the value of the attenuation coefficient at the end of the data range,¹² ${\hat{\mu }}_{\mathcal{E}}={\mu }_{\mathrm{OCT}}(\mathcal{E})$. This allows us to rewrite Eq. (13) as an estimator of ${\mu }_{\mathrm{OCT}}$:

Equation (15) reveals that the estimate ${\widehat{\mu }}_{\mathrm{OCT}}(z)$ approaches the true value ${\mu }_{\mathrm{OCT}}$ at a $z$-position sufficiently far from the end of range $\mathcal{E}$, whereas ${\widehat{\mu }}_{\mathrm{OCT}}(z)$ approaches the estimate ${\widehat{\mu }}_{\mathcal{E}}$ as $z$ approaches $\mathcal{E}$.

8.

Appendix B

Vermeer et al.⁹ considered the effect of discretization of $I(z)$. Each data point $I[i]$ corresponds to the integration of Eq. (3) over a finite pixel size $\mathrm{\Delta }$ around $z$. They show that the discretized version of Eq. (11) reads

We can use Eq. (16) to obtain an expression for $\mathcal{C}$. We have ${\widehat{\mu }}_{\mathrm{OCT}}[i_{\max }]=\frac{1}{2\mathrm{\Delta }}\ln (1+\frac{I[i_{\max }]}{{\sum }_{j=i_{\max }+1}^{\infty }I[j]})=\frac{1}{2\mathrm{\Delta }}\ln (1+\frac{I[i_{\max }]}{\mathcal{C}})$ so that $\mathcal{C}=I[i_{\max }]/(\exp (2{\mu }_{\mathcal{E}}\mathrm{\Delta })-1)$. Here, ${\mu }_{\mathcal{E}}={\mu }_{\mathrm{OCT}}[i_{\max }]$ is an independent estimate of the attenuation coefficient at the end of the range, as before (Appendix A).

Quite often, approximate forms of Eq. (7) / Eq. (17) are found in literature, which are obtained by linearization of the logarithmic and exponential terms. Then, factor $\mathcal{C}\approx I[i_{\max }]/2{\mu }_{\mathcal{E}}\mathrm{\Delta }$ and upon expanding the logarithmic term $\ln (1+x)=x-\frac{1}{2}{x}^2+\dots$:

9.

Appendix C

We seek the precision of the DRE method in the presence of noise. We make use of the fact that the attenuation coefficient will be approximately normally distributed, and that the precision is given by the standard deviation ${\sigma }_{{\widehat{\mu }}_{\mathrm{OCT}}}$ of that distribution. We use the “simplified” form of the depth resolved method

After pre-averaging $N>30$ times, the intensity values are normally distributed. The term $D(i)=\sum _{j=i+1}^{M}I(j)+\mathcal{C}$ in the denominator is then the sum of normally distributed random variables plus a constant, which yields a new normal variable with mean ${\mathcal{m}}_D(i)=\sum _{j=i+1}^M{\mathcal{m}}_I(j)+\mathcal{C}$ and variance ${\sigma }_D^2(i)=\sum _{j=i+1}^M[{\sigma }_I^2(j)+{\sigma }_{\varsigma }^2]$; that is, the means and variances simply add up. When the coefficients of variation of either the nominator ${\delta }_I={\sigma }_I/{\mathcal{m}}_I$ or the denominator ${\delta }_D={\sigma }_D/{\mathcal{m}}_D$ in the ratio is $<1$ (in fact, both are), the result is also normally distributed with mean

The variance is given as

Fiske et al.¹⁸ showed that the coefficient of variation ${\delta }_D\ll 1$ already, even without pre-averaging so we neglect that term in Eq. (21). See also Appendix D for further justification.

Thus, upon averaging $N>30$ of individual A-scans followed by subtraction of the mean noise level $⟨{\varsigma }_I⟩$ [Eq. (2)], the signal $I(i)$ is obtained, mean ${\mathcal{m}}_I\equiv ⟨I(i)⟩$ and variance ${\sigma }_I^2=(\frac{{⟨I(i)⟩}^2}{N}+\frac{{⟨{\varsigma }_I⟩}^2}{N})=(\frac{{\mathcal{m}}_I^2}{N}+\frac{{⟨{\varsigma }_I⟩}^2}{N})$. These relations indicate that subtracting the mean noise floor does not remove the fluctuations caused by noise. The required expression for ${\delta }_I^2={\sigma }_I^2/{\mathcal{m}}_I^2$ (the square of the coefficient of variation CV) is

10.

Appendix D

Let $D(i)=\sum _{j=i+1}^{M}I(j)+\mathcal{C}$ be the sum term in the denominator of Eq. (19). It is the sum of normally distributed random variables plus a constant, which yields a new normal variable with mean ${\mathcal{m}}_D(i)=\sum _{j=i+1}^M{\mathcal{m}}_I(j)+\mathcal{C}$ and variance ${\sigma }_D^2(i)=\sum _{j=i+1}^M[{\sigma }_I^2(j)+{\sigma }_{\varsigma }^2]$; that is, the means and variances simply add up. The square of the coefficient of variation is

We will first work out the term without noise, then the term with noise, making use of the fact that the analysis is based on the exponential decay model so we can write inside the summations: ${\mathcal{m}}_I(j)={\mathcal{m}}_I(i)\exp (-2{\mu }_{\mathrm{OCT}}\mathrm{\Delta }·(j-i))$.

Let $k=j-i$, $x=e^{-2{\mu }_{\mathrm{OCT}}\mathrm{\Delta }}$ and use the identity:

Then, the mean of $D(i)$ follows as

And the variance as

Inserting the expressions for the SNR and regularization term $\mathcal{C}$:

To arrive at a more compact, albeit approximate expression we first linearize the exponentials

Then assume ${\mu }_{\mathrm{OCT}}\mathrm{\Delta }\ll 1$ and rearrange to

Comparing to Eq. (22), we see that ${\delta }_D^2(i)={\delta }_I^2(i)\times {(2{\widehat{\mu }}_{\mathcal{E}}\mathrm{\Delta })}^2/(M-i)$. Since ${({\hat{\mu }}_{\mathcal{E}})}^{-1}$ is in the order of $mm$ and $\mathrm{\Delta }$ is in the order of $\mu \mathrm{m}$, neglecting ${\delta }_D^2(i)$ as is done in Appendix C is justified.