Accuracy and precision of depth-resolved estimation of attenuation coefficients in optical coherence tomography

Abstract

Significance

Parametric imaging of the attenuation coefficient μ_OCT using optical coherence tomography (OCT) is a promising approach for evaluating abnormalities in tissue. To date, a standardized measure of accuracy and precision of μ_OCT by the depth-resolved estimation (DRE) method, as an alternative to least squares fitting, is missing.

Aim

We present a robust theoretical framework to determine accuracy and precision of the DRE of μ_OCT.

Approach

We derive and validate analytical expressions for the accuracy and precision of μ_OCT determination by the DRE using simulated OCT signals in absence and presence of noise. We compare the theoretically achievable precisions of the DRE method and the least-squares fitting approach.

Results

Our analytical expressions agree with the numerical simulations for high signal-to-noise ratios and qualitatively describe the dependence on noise otherwise. A commonly used simplification of the DRE method results in a systematic overestimation of the attenuation coefficient in the order of μOCT2×Δ , where Δ is the pixel stepsize. When μ_OCT · | AFR | ≲ 1.8, μ_OCT is reconstructed with higher precision by the depth-resolved method compared to fitting over the length of an axial fitting range | AFR | .

Conclusions

We derived and validated expressions for the accuracy and precision of DRE of μ_OCT. A commonly used simplification of this method is not recommended as being used for OCT-attenuation reconstruction. We give a rule of thumb providing guidance in the choice of estimation method.

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:⁷

[Formula omitted. See PDF] where |AFR| is the length of the axial fitting range (AFR), [Formula omitted. See PDF] is the number of independent data points in the AFR, and [Formula omitted. See PDF] is the number of A-scans averaged prior to fitting. The constant [Formula omitted. See PDF] originates from the Rayleigh distribution of amplitude values corresponding to fully developed speckle. The lower bound given by Eq. (1) is expressed in the same units as the attenuation coefficient, e.g., [Formula omitted. See PDF] and corresponds to the standard deviation of the normal distribution of attenuation coefficients that would be obtained by repeating the fitting procedure a large number of times, each time with a different, random, realization of the speckle pattern. The main feature of Eq. (1) is that the precision is independent of the attenuation coefficient itself but only depends on the parameters used in the fitting procedure. In the derivation of Eq. (1), we assumed that shot noise was negligible. Our results validated this assumption when signals-to-noise ratios (SNRs) within the AFR exceed 20 dB.⁷

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 [Formula omitted. See PDF] is modeled as a single exponential decay combined with the confocal point spread function and sensitivity roll-off:²

[Formula omitted. See PDF] where [Formula omitted. See PDF] is a conversion factor that includes the detector response, [Formula omitted. See PDF] is the confocal point spread function, [Formula omitted. See PDF] describes the sensitivity roll-off in depth for non-time domain OCT, [Formula omitted. See PDF] is the backscattering coefficient within the numerical aperture (NA) of the detection optics; [Formula omitted. See PDF] is the OCT attenuation coefficient that contains contributions from both scattering and absorption. In the absence of multiple forward scattering, which we will assume henceforth, [Formula omitted. See PDF] is the sum of scattering and absorption coefficients. The backscattering coefficient is proportional to the scattering coefficient [Formula omitted. See PDF] through a phase function and NA dependent factor [Formula omitted. See PDF]. The mean squared noise background is given by [Formula omitted. See PDF]. Upon noise subtraction and following correction for point spread function and roll-off, we arrive at [Formula omitted. See PDF]

We proceed to compute the definite integral [Formula omitted. See PDF] of Eq. (3), which runs from the depth [Formula omitted. See PDF] up to the end of the available (or used) data range [Formula omitted. See PDF], to estimate the attenuation coefficient as (Appendix A):

[Formula omitted. See PDF]

Compared to the original formulation by Vermeer, Eq. (4) contains a regularization term [Formula omitted. See PDF] in the denominator to compensate for the finite data range.¹⁶ Here, [Formula omitted. See PDF] is an independently obtained estimate for the attenuation coefficient at the end of the data range [Formula omitted. See PDF], 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 [Formula omitted. See PDF] 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 [Formula omitted. See PDF], 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:

[Formula omitted. See PDF]

The mean value [Formula omitted. See PDF] can be obtained from a large set of estimations of [Formula omitted. See PDF], in practice over some spatial range around [Formula omitted. See PDF] and/or from several A-scans at the same position. Combining Eqs. (3) and (4), we theoretically obtain

[Formula omitted. See PDF]

Thus, [Formula omitted. See PDF] at a location sufficiently far from [Formula omitted. See PDF], whereas [Formula omitted. See PDF] as [Formula omitted. See PDF] approaches [Formula omitted. See PDF]. Vermeer considered the effect of discretization of [Formula omitted. See PDF], i.e., each datapoint [Formula omitted. See PDF] corresponds to the integration of Eq. (3) over a finite pixel size [Formula omitted. See PDF] around [Formula omitted. See PDF]. The exact, discretized version of Eq. (4) reads (Appendix B)

[Formula omitted. See PDF] where [Formula omitted. See PDF] is the pixel index corresponding to the end of the data range. The factor [Formula omitted. See PDF] is the discretized equivalent of the term [Formula omitted. See PDF] in Eq. (4).

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)]

[Formula omitted. See PDF]

Moreover, some authors further omit the (linearized) regularization term [Formula omitted. See PDF] 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 [Formula omitted. See PDF].

We now seek the precision with which [Formula omitted. See PDF] 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 [Formula omitted. See PDF] -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 [Formula omitted. See PDF] independent estimations from an exponential distribution is [Formula omitted. See PDF] so the CR lower bound, expressed as standard deviation becomes [Formula omitted. See PDF]. 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 [Formula omitted. See PDF], with mean equal to [Formula omitted. See PDF] and variance equal to [Formula omitted. See PDF]. The corresponding Fisher information for the estimated attenuation coefficient based on [Formula omitted. See PDF] independent measurements then becomes [Formula omitted. See PDF] and the CR lower bound, expressed as standard deviation is thus [Formula omitted. See PDF]. 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

[Formula omitted. See PDF] where the SNR is defined per pixel as [Formula omitted. See PDF].

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 [Formula omitted. See PDF] for CF. It also shows that a higher precision can be obtained using DRE compared to CF, when the AFR becomes smaller than [Formula omitted. See PDF] 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 [Formula omitted. See PDF], where [Formula omitted. See PDF] is the amplitude variance, [Formula omitted. See PDF] is given by Eq. (3), and [Formula omitted. See PDF] is a uniformly distributed number between 0 and 1 drawn for each pixel. This procedure assures that the amplitudes [Formula omitted. See PDF] follow a Rayleigh distribution corresponding to fully developed speckle and the contribution of shot noise. [Formula omitted. See PDF] 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 [Formula omitted. See PDF].

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 [Formula omitted. See PDF] omitted. This procedure was repeated [Formula omitted. See PDF] times to obtain a distribution of [Formula omitted. See PDF] 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, [Formula omitted. See PDF] was set to unity and an arbitrary scaling factor of [Formula omitted. See PDF] 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 [Formula omitted. See PDF]. Values of 2 and [Formula omitted. See PDF] 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 [Formula omitted. See PDF] for the estimation of the attenuation coefficient at the end of the data range [Formula omitted. See PDF].

4. Results

Figure 1(a) shows an example of simulated [Formula omitted. See PDF] times averaged A-scans, including shot noise, obtained by the procedure outlined in the previous section. The arrow indicates the position [Formula omitted. See PDF], at which the intensity signal hits the noise floor with the condition [Formula omitted. See PDF] for the simulation with [Formula omitted. See PDF] as input parameter. For [Formula omitted. See PDF] the intensity signal does not reach the noise floor.

Figure 1(b) shows a histogram of intensity values at a depth of [Formula omitted. See PDF] using [Formula omitted. See PDF] obtained from of [Formula omitted. See PDF] independent simulations. It shows that the averaged intensities, obtained by pre-averaging [Formula omitted. See PDF] 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 [Formula omitted. See PDF]. Both figures show the estimated attenuation coefficients versus depth using Eqs. (7) and (8) with [Formula omitted. See PDF] 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 [Formula omitted. See PDF]; even after pre-averaging [Formula omitted. See PDF] A-scans. The data shown in Fig. 3(b) is averaged over [Formula omitted. See PDF] independent simulations and therefore permits closer comparison of theory and simulations. For the data shown in both figures, a value of [Formula omitted. See PDF] is used in the regularization term [Formula omitted. See PDF] in Eq. (7). The results in Fig. 2(b) demonstrate that the estimated attenuation coefficient differs [Formula omitted. See PDF] from the true value up to a depth of 1.98 mm for [Formula omitted. See PDF] and, for [Formula omitted. See PDF], over the entire depth range. It approaches, in both cases, [Formula omitted. See PDF] near the end of the available data range. If, on the other hand, the linearized approximation Eq. (8) with [Formula omitted. See PDF] 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.

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 [Formula omitted. See PDF] 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 [Formula omitted. See PDF]. 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 [Formula omitted. See PDF] (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 [Formula omitted. See PDF]. 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 [Formula omitted. See PDF] of true attenuation value up to a depth of 1.18 mm for [Formula omitted. See PDF] and, for [Formula omitted. See PDF], over the entire depth range similar to the results in Fig. 2.

We proceed to compare the standard deviation [Formula omitted. See PDF] 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 [Formula omitted. See PDF] 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.

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 ([Formula omitted. See PDF]). Figure 5 shows the numerically obtained CR lower bound (the minimal precision for the CF approach) in dependence of [Formula omitted. See PDF] for an AFR of [Formula omitted. See PDF] ([Formula omitted. See PDF] points, [Formula omitted. See PDF]) located well before [Formula omitted. See PDF] such that the SNR in the AFR is [Formula omitted. See PDF]. 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 [Formula omitted. See PDF] points, with end-of-range values [Formula omitted. See PDF] set to the true value of [Formula omitted. See PDF].) does show a dependence on the attenuation coefficient. All curves follow a [Formula omitted. See PDF] trend, whereas [Formula omitted. See PDF] 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 [Formula omitted. See PDF].

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 [Formula omitted. See PDF] 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 [Formula omitted. See PDF], where [Formula omitted. See PDF] 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 [Formula omitted. See PDF] 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 [Formula omitted. See PDF] 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 [Formula omitted. See PDF]. 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 [Formula omitted. See PDF] 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 [Formula omitted. See PDF] when the same number of A-scans ([Formula omitted. See PDF]) is pre-averaged, and the same number of independent data points ([Formula omitted. See PDF]) 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.

5.1. Limitations

Simulations can be performed quickly and at low cost compared to the time and resources required for phantom experiments. Manufacturing of phantoms with precise control of the scattering properties can be challenging whereas simulations can explore a much wider range of parameter space beyond what is feasible in experiments. Simulations allow for precise control and manipulation of individual parameters (or tuning correlation among them), leading to deeper understanding of the underlying mechanisms. Although we have included only a limited number of [Formula omitted. See PDF] values and SNRs in the present article, our simulations are straightforwardly extended to include a wider range of scattering and absorption coefficients; to add point spread function and sensitivity roll-off [Eq. (2)] and introduce layers with their specific optical properties.

Performance of both the CF and DRE methods depends on the appropriateness of the underlying single-exponential decay model of Eq. (3) to describe the light–tissue interaction (assuming instrumental factors are corrected appropriately). In the this study, both the simulation generating OCT data and analysis were based on the same Eq. (3). Therefore, the accuracy and precision derived in this article represent the best values that can be obtained. This best-case scenario may not be true for experimental data, for instance when a small fraction of multiple scattering occurs. Multiple scattering models are available²⁰ and can be adapted for CF, leading to the inclusion of one or more fit parameters describing tissue scattering (e.g., the root mean square scattering angle or scattering anisotropy). However, adaptation of these models for use in DRE seems challenging. In practice, even in the presence of multiple scattering, the part of the signal decay caused by absorption and scattering is often adequately modeled as a single exponential decay, albeit with a decay constant [Formula omitted. See PDF] (because multiple scattering causes more light to be detected, than expected based on the single-scattering model). Thus, the adoption of [Formula omitted. See PDF] allows us to describe tissue attenuation as measured by OCT as an effective parameter that does not require an estimate of the relative weight of single and multiple scattering contributions. When applying the DRE method to multiple layers with varying optical properties, the DRE method generally fails to extract the correct optical properties, unless for each layer [Formula omitted. See PDF] only (e.g., no absorption) and [Formula omitted. See PDF] is a constant throughout the sample.²¹ Whether or not these conditions are met in practice should ultimately be verified by experiments while the level of inaccuracy and imprecision may be estimated using simulations.

5.2. Clinical Implications

Measurements of the attenuation coefficient complement the structural images that OCT provides. The main premise is that quantification of [Formula omitted. See PDF] can be used to distinguish different tissue types (e.g., benign versus malignant). Ideally, with perfect accuracy and precision, the sensitivity and specificity of such an approach are determined by the degree of biological variation within, and between the different tissue types. In practice, however, the accuracy and precision will be finite. The results laid out in this article, as well is in the previous publication⁷ allow us to determine to which degree the observed variation in attenuation coefficients is due to the employed method, and which part can be attributed to biological variation.

6. Conclusion

In this article, we derived and validated the accuracy and precision of the depth resolved estimation method of the attenuation coefficient [Formula omitted. See PDF] in optical coherence tomography. We showed that a commonly used simplification of the method results in loss of accuracy in the order of [Formula omitted. See PDF] where [Formula omitted. See PDF] 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 [Formula omitted. See PDF], 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 [Formula omitted. See PDF] of Eq. (3), which yields

[Formula omitted. See PDF]

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

[Formula omitted. See PDF]

In practice, data are only available over a finite range, up to [Formula omitted. See PDF]. We compute the definite integral [Formula omitted. See PDF], which yields

[Formula omitted. See PDF]

Solving Eq. (12) for [Formula omitted. See PDF], and substituting the result in Eq. (11) gives

[Formula omitted. See PDF]

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

[Formula omitted. See PDF]

When [Formula omitted. See PDF] approaches the end of range [Formula omitted. See PDF], the finite integral in the denominator term goes to zero; therefore [Formula omitted. See PDF] 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,¹² [Formula omitted. See PDF]. This allows us to rewrite Eq. (13) as an estimator of [Formula omitted. See PDF]:

[Formula omitted. See PDF]

Equation (15) reveals that the estimate [Formula omitted. See PDF] approaches the true value [Formula omitted. See PDF] at a [Formula omitted. See PDF]-position sufficiently far from the end of range [Formula omitted. See PDF], whereas [Formula omitted. See PDF] approaches the estimate [Formula omitted. See PDF] as [Formula omitted. See PDF] approaches [Formula omitted. See PDF].

8. Appendix B

Vermeer et al.⁹ considered the effect of discretization of [Formula omitted. See PDF]. Each data point [Formula omitted. See PDF] corresponds to the integration of Eq. (3) over a finite pixel size [Formula omitted. See PDF] around [Formula omitted. See PDF]. They show that the discretized version of Eq. (11) reads

[Formula omitted. See PDF] where [Formula omitted. See PDF] is now the estimate of the average attenuation coefficient in the [Formula omitted. See PDF]’th pixel. Considering that data are only available over a finite data range we write [Formula omitted. See PDF] where [Formula omitted. See PDF] is the pixel index corresponding to the end of range [Formula omitted. See PDF] and [Formula omitted. See PDF].

We can use Eq. (16) to obtain an expression for [Formula omitted. See PDF]. We have [Formula omitted. See PDF] so that [Formula omitted. See PDF]. Here, [Formula omitted. See PDF] 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 [Formula omitted. See PDF] and upon expanding the logarithmic term [Formula omitted. See PDF]:

[Formula omitted. See PDF] after which only the first term is retained. Under that same approximation, the second term in Eq. (18) is approximately equal to [Formula omitted. See PDF] and we conclude that linearization of Eq. (7)/(17) leads to a systematic overestimation of the attenuation coefficient in the order of [Formula omitted. See PDF].

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 [Formula omitted. See PDF] of that distribution. We use the “simplified” form of the depth resolved method

[Formula omitted. See PDF]

After pre-averaging [Formula omitted. See PDF] times, the intensity values are normally distributed. The term [Formula omitted. See PDF] in the denominator is then the sum of normally distributed random variables plus a constant, which yields a new normal variable with mean [Formula omitted. See PDF] and variance [Formula omitted. See PDF]; that is, the means and variances simply add up. When the coefficients of variation of either the nominator [Formula omitted. See PDF] or the denominator [Formula omitted. See PDF] in the ratio is [Formula omitted. See PDF] (in fact, both are), the result is also normally distributed with mean

[Formula omitted. See PDF]

The variance is given as

[Formula omitted. See PDF]

Fiske et al.¹⁸ showed that the coefficient of variation [Formula omitted. See PDF] already, even without pre-averaging so we neglect that term in Eq. (21). See also Appendix D for further justification.

Thus, upon averaging [Formula omitted. See PDF] of individual A-scans followed by subtraction of the mean noise level [Formula omitted. See PDF] [Eq. (2)], the signal [Formula omitted. See PDF] is obtained, mean [Formula omitted. See PDF] and variance [Formula omitted. See PDF]. These relations indicate that subtracting the mean noise floor does not remove the fluctuations caused by noise. The required expression for [Formula omitted. See PDF] (the square of the coefficient of variation CV) is

[Formula omitted. See PDF] where the SNR is now defined as [Formula omitted. See PDF]. Combining with Eq. (21), assuming [Formula omitted. See PDF] independent estimations, the precision expressed as standard deviation becomes [Formula omitted. See PDF]

10. Appendix D

Let [Formula omitted. See PDF] 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 [Formula omitted. See PDF] and variance [Formula omitted. See PDF]; that is, the means and variances simply add up. The square of the coefficient of variation is

[Formula omitted. See PDF]

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: [Formula omitted. See PDF].

Let [Formula omitted. See PDF], [Formula omitted. See PDF] and use the identity:

[Formula omitted. See PDF]

Then, the mean of [Formula omitted. See PDF] follows as

[Formula omitted. See PDF]

And the variance as

[Formula omitted. See PDF]

Inserting the expressions for the SNR and regularization term [Formula omitted. See PDF]:

[Formula omitted. See PDF]

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

[Formula omitted. See PDF] [Formula omitted. See PDF]

Then assume [Formula omitted. See PDF] and rearrange to

[Formula omitted. See PDF]

Comparing to Eq. (22), we see that [Formula omitted. See PDF]. Since [Formula omitted. See PDF] is in the order of [Formula omitted. See PDF] and [Formula omitted. See PDF] is in the order of [Formula omitted. See PDF], neglecting [Formula omitted. See PDF] as is done in Appendix C is justified.

Disclosure

The authors have no relevant financial interest in this article and no potential conflicts of interest to disclose.

Code, Data, and Materials Availability

The code and data are available from the authors upon request.

Acknowledgments

This publication is part of the project “An integrated Optical Coherence Tomography system for medical imaging at 1300 nm” (with project number 16251) of the research program HTSM2017, which is partly financed by the Dutch Research Council (NWO).

Linda B. Neubrand received her BSc and MSc degrees in physics from the Institute of Technology (KIT), Karlsruhe, Germany, in 2017 and 2019. She performed her master’s project in the University of Amsterdam (UvA), Amsterdam, The Netherlands. She is currently pursuing her PhD at the Department of Biomedical Engineering and Physics, Amsterdam, UMC, University of Amsterdam, Amsterdam, The Netherlands.

Ton G. van Leeuwen is the head of the BME and Physics Department at Amsterdam UMC. His research focuses on the physics of the interaction of light with tissue and to use that knowledge for the development, introduction, and clinical evaluation of (newly developed) optical imaging and analysis techniques. Key in this research is to use this knowledge and new devices to gather quantitative functional information of tissue or tissue sample.

Dirk J. Faber received his MSc degree in applied physics from the University of Twente, Enschede, The Netherlands, in 1999, and his PhD from the University of Amsterdam, Amsterdam, The Netherlands, in 2005, based on his work on OCT. He is currently an assistant professor with Department of Biomedical Engineering and Physics, at the Amsterdam University Medical Centers, Academic Medical Center. His current research focuses on the physics of light–tissue interaction and the development of OCT and single fiber reflectance spectroscopy. He is a senior member of SPIE.

References

1. J. M. Schmitt, A. Knüttel and R. F. Bonner, “Measurement of optical properties of biological tissues by low-coherence reflectometry,” Appl. Opt., 32 6032 –6042 https://doi.org/10.1364/AO.32.006032 APOPAI 0003-6935 (1993).

2. P. Gong et al., “Parametric imaging of attenuation by optical coherence tomography: review of models, methods, and clinical translation,” J. Biomed. Opt., 25 (4), 040901 https://doi.org/10.1117/1.JBO.25.4.040901 JBOPFO 1083-3668 (2020).

3. D. J. Faber et al., “Quantitative measurement of attenuation coefficients of weakly scattering media using optical coherence tomography,” Opt. Express, 12 4353 –4365 https://doi.org/10.1364/OPEX.12.004353 OPEXFF 1094-4087 (2004).

4. J. M. Schmitt, S. H. Xiang and K. M. Yung, “Speckle in optical coherence tomography: an overview,” J. Biomed. Opt., 4 (1), 95 –105 https://doi.org/10.1117/1.429925 (1999).

5. M. Almasian, T. G. van Leeuwen and D. J. Faber, “OCT amplitude and speckle statistics of discrete random media,” Sci. Rep., 7 14873 https://doi.org/10.1038/s41598-017-14115-3 SRCEC3 2045-2322 (2017).

6. A. Weatherbee et al., “Probability density function formalism for optical coherence tomography signal analysis: a controlled phantom study,” Opt. Lett., 41 2727 –2730 https://doi.org/10.1364/OL.41.002727 OPLEDP 0146-9592 (2016).

7. L. B. Neubrand, T. G. van Leeuwen and D. J. Faber, “Precision of attenuation coefficient measurements by optical coherence tomography,” J. Biomed. Opt., 27 (8), 085001 https://doi.org/10.1117/1.JBO.27.8.085001 JBOPFO 1083-3668 (2022).

8. R. Wessels et al., “Optical biopsy of epithelial cancers by optical coherence tomography (OCT),” Lasers Med. Sci., 29 1297 –1305 https://doi.org/10.1007/s10103-013-1291-8 (2014).

9. K. A. Vermeer et al., “Depth-resolved model-based reconstruction of attenuation coefficients in optical coherence tomography,” Biomed. Opt. Express, 5 322 –337 https://doi.org/10.1364/BOE.5.000322 BOEICL 2156-7085 (2014).

10. M. J. A. Girard et al., “Shadow removal and contrast enhancement in optical coherence tomography images of the human optic nerve head,” Invest. Ophthalmol. Vis. Sci., 52 (10), 7738 –7748 https://doi.org/10.1167/iovs.10-6925 IOVSDA 0146-0404 (2011).

11. D. I. Hughes and F. A. Duck, “Automatic attenuation compensation for ultrasonic imaging,” Ultrasound Med. Biol., 23 (5), 651 –664 https://doi.org/10.1016/S0301-5629(97)00002-1 USMBA3 0301-5629 (1997).

12. J. Liu et al., “Optimized depth-resolved estimation to measure optical attenuation coefficients from optical coherence tomography and its application in cerebral damage determination,” J. Biomed. Opt., 24 (3), 035002 https://doi.org/10.1117/1.JBO.24.3.035002 JBOPFO 1083-3668 (2019).

13. G. T. Smith et al., “Automated, depth-resolved estimation of the attenuation coefficient from optical coherence tomography data,” IEEE Trans. Med. Imaging, 34 (12), 2592 –2602 https://doi.org/10.1109/TMI.2015.2450197 ITMID4 0278-0062 (2015).

14. N. Dwork et al., “Automatically determining the confocal parameters from OCT B-scans for quantification of the attenuation coefficients,” IEEE Trans. Med. Imaging, 38 (1), 261 –268 https://doi.org/10.1109/TMI.2018.2861570 ITMID4 0278-0062 (2019).

15. S. Liu et al., “Tissue characterization with depth-resolved attenuation coefficient and backscatter term in intravascular optical coherence tomography images,” J. Biomed. Opt., 22 (9), 096004 https://doi.org/10.1117/1.JBO.22.9.096004 JBOPFO 1083-3668 (2017).

16. K. Li et al., “Robust, accurate depth-resolved attenuation characterization in optical coherence tomography,” Biomed. Opt. Express, 11 672 –687 https://doi.org/10.1364/BOE.382493 BOEICL 2156-7085 (2020).

17. M. M. Amaral et al., “General model for depth-resolved estimation of the optical attenuation coefficients in optical coherence tomography,” J. Biophotonics, 12 (10), e201800402 https://doi.org/10.1002/jbio.201800402 (2019).

18. L. D. Fiske et al., “Bayesian analysis of depth resolved OCT attenuation coefficients,” Sci. Rep., 11 2263 https://doi.org/10.1038/s41598-021-81713-7 SRCEC3 2045-2322 (2021).

19. J. Kübler et al., “Investigation of methods to extract confocal function parameters for the depth resolved determination of attenuation coefficients using OCT in intralipid samples, titanium oxide phantoms, and in vivo human retinas,” Biomed. Opt. Express, 12 6814 –6830 https://doi.org/10.1364/BOE.440574 BOEICL 2156-7085 (2021).

20. L. Thrane, H. T. Yura and P. E. Andersen, “Analysis of optical coherence tomography systems based on the extended Huygens–Fresnel principle,” J. Opt. Soc. Am. A, 17 484 –490 https://doi.org/10.1364/JOSAA.17.000484 JOAOD6 0740-3232 (2000).

21. T. M. Cannon, B. E. Bouma and N. Uribe-Patarroyo, “Layer-based, depth-resolved computation of attenuation coefficients and backscattering fractions in tissue using optical coherence tomography,” Biomed. Opt. Express, 12 5037 –5056 https://doi.org/10.1364/BOE.427833 BOEICL 2156-7085 (2021).