## Abstract

The attenuation coefficient (AC) is an optical property of tissue that can be estimated from optical coherence tomography (OCT) data. In this paper, we aim to estimate the AC accurately by compensating for the shape of the focused beam. For this, we propose a method to estimate the axial PSF model parameters and AC by fitting a model for an OCT signal in a homogenous sample to the recorded OCT signal. In addition, we employ numerical analysis to obtain the theoretical optimal precision of the estimated parameters for different experimental setups. Finally, the method is applied to OCT B-scans obtained from homogeneous samples. The numerical and experimental results show accurate estimations of the AC and the focus location when the focus is located inside the sample.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Optical coherence tomography (OCT) has been widely used to capture structural information of tissues in clinical tasks such as the diagnosis of retinal and vascular diseases. Previously, extracting valuable information embedded in the signal intensity of constituting tissues has been investigated. An optical property such as the attenuation coefficient (AC) offers valuable information that can be estimated from the intensity of the OCT signal. It has the potential to act as a biomarker for the diagnosis and monitoring of chorioretinal diseases [1], breast tumor lesions [2], renal tumors [3,4], oral cancer [5], rectal cancer [6] and several other applications such as atherosclerotic plaque characterization [7–9]. Several methods based on single [10,11,12] and multiple [13,14] scattering of light have been presented for estimating the attenuation coefficient in a homogeneous medium using OCT. Recently, a depth-resolved single-scattering based method has been developed by Vermeer *et al.* [15] for estimating attenuation coefficients in inhomogeneous mediums, e.g. in tissues. For all methods that estimate the attenuation coefficient, the OCT signal must be corrected for: 1) the depth-dependent noise floor [16]; 2) the so-called roll-off, i.e. the depth-dependent signal decay caused by discrete signal detection and resolution limitations of the detection process [15,17]; and 3) the axial point spread function (PSF), which, for a Gaussian-shaped beam, is governed by the effective Rayleigh length around the focus position of the beam [15,18]. Compensation for noise and roll-off is nowadays a standard procedure, which can be done with a function obtained from a fit to reference data [17]. However, in order to correct for the axial PSF, in many cases its model parameters need to be estimated from the acquired data since the effective Rayleigh length and focus depend on the optical system, e.g. in case the cornea and lens. We showed in previous work how the attenuation coefficient is sensitive to an error in the estimated parameters of the axial PSF model [19]. Therefore, accurate and precise estimation of these parameters is required to achieve an unbiased and precise estimation of attenuation coefficient. Various methods have been developed to estimate the attenuation coefficient of the tissue while taking into account the effect of the beam shape that influence the acquired OCT signal. Smith et al. [20] compensate for the effect of focus using an existing model of the shape of the beam. However, in their work the parameters of the shape of the beam need to be known in advance. In many medical applications, such as ophthalmology, the location of the focal point varies and there is a need for a method to automatically estimate the focus location to compensate for the effect of the beam shape in the estimation of the attenuation coefficient. Stefan et al. [21] introduced a method to estimate the attenuation coefficient using two B-scans to first estimate the location of focus and afterwards estimating the attenuation coefficient from a single scattering model of the OCT light after compensating for the effect of beam shape. This method is dependent on having identical A-lines to be able to eliminate the effect of attenuation coefficients. However, this proposed method was only tested with static samples where the identical physical location in both B-scans is feasible and the factors such as beam’s angle of incidence can be controlled to ensure a similar tissue attenuation coefficient. Another limitation of this method is the necessity to have access to two scans from the same position in the tissue. However, in many clinical data, such as retinal scans, only one averaged measurement of the same tissue’s location is available.

In this paper, we aim to achieve an accurate estimate of attenuation coefficient by compensating for all of the aforementioned effects on the recorded OCT signal. To do so, we propose a method to estimate the axial PSF model parameters (focus depth and Rayleigh length in the medium) and attenuation coefficient by fitting a single scattering based model for a homogenous sample OCT signal to the recorded OCT signal after subtraction of the depth-dependent noise floor and compensating for roll-off. In addition, a Cramér-Rao analysis is performed to theoretically determine the attainable precision of the estimated parameters and to investigate the limitations of the proposed procedure for various experimental configurations. Monte Carlo simulations of the estimation method are performed to evaluate the robustness of the method and compare the precision of the theoretical lower bound produced by the Cramér-Rao analysis and to show a possible bias in the estimated parameters. Finally, the method is applied to B-scans obtained with an experimental OCT system from homogeneous samples with various concentrations of TiO_{2} particles dispersed in silicone to assess the precision and accuracy of the method.

## 2. Method

In this section, we introduce a method for accurate estimation of the attenuation coefficient in a homogenous (or single layer) sample by compensating the recorded OCT signal for the noise floor, roll-off, and axial PSF.

#### 2.1 Estimating the model parameters

In a single scattering model of light presented by Faber et al. [22], the Fourier-domain OCT signal at physical depth *z* in a homogeneous sample may be expressed by,

*R*(

*z*)), the axial PSF modeled by a Cauchy function at focus position

*z*

_{0}and scaled by the Rayleigh length

*z*[18], and the signal attenuation modeled by the attenuation coefficient

_{R}*μ*and scaling factor

*C*. The second term,

*N*(

*z*) is the depth-dependent noise floor and can be obtained by averaging over a large number of A-lines without a sample in the sample arm of the OCT system. The intensity of the OCT signal has an exponential distribution caused by speckle noise. However, due the central limit theorem, by averaging over a sufficiently large number of neighboring A-lines (>30 based on rule of thumb) with exponential distributions, the averaged OCT signal at depth

*z*tends toward a normal distribution ${{\cal N}}[{m(z ),m{{(z )}^2}/N} )]$, with

*m*(

*z*) being the expected value of the exponential distribution, and $m{(z )^2}/N$ the variance of the resulting normal distribution. The third term

*ɛ(z)*represents this speckle noise. In addition, the roll-off can be measured and the signal can be corrected for roll-off by performing the operation $A(z) = (S(z) - N(z))/R(z)$ [15].

We estimate the model parameters of the axial PSF and the attenuation coefficient using a maximum likelihood estimator. For this, Eq. (1) was fitted to the measurements. For an averaged A-line $A(z )$ with ${N_D}\; $ data-points as a function of *z*, the independent parameters $C,\; $*µ*, *z _{0}* and ${z_R}$ can be estimated by minimizing the sum of squared residuals, $\chi ,$ given by,

*j*is an index that denotes the data-point number on each averaged A-line.

#### 2.2 Model selection and evaluation

To design a reliable model-based method for estimating the axial PSF from recorded data, we studied the influence of integrating prior information into the model, such as a known or joint model parameter among multiple averaged A-lines, to reduce the degrees of freedom and aiming to thereby improve the estimation precision of the remaining parameters. Moreover, the attainable precision of the estimated parameters $\{{{\theta_1}, \ldots ,{\theta_N}} \}= \{ C,\mu ,{z_0},{z_R}\}$ needs to be calculated for various experimental setups. Exploring the precision of the estimated parameters such as the focus depth into the sample, the Rayleigh length, illumination intensity and the attenuation coefficient of the medium, enables us to optimize the experimental design and to know the limitations of the proposed method. For these purposes, a Cramér-Rao analysis was applied using a derivation of the Fisher information matrix for a Gaussian noise model (see equations (9)-11 from Caan et al. [23]). Cramer-Rao analysis is limited to finding the minimal variance of the model parameters assuming an unbiased estimator. To evaluate the optimal precision of the estimated parameters and to compare different models and experimental setups, we use the relative errors, as provided by the diagonal elements of the relative Cramér-Rao lower bound (rCRLB) matrix [23]. The diagonal elements are the relative theoretical lower bounds on the variance of the unbiased estimators of each parameter. We intuitively considered an estimation error lower than 10% to be acceptable for the purpose of this paper.

Multiple averaged A-lines can be used to estimate the model parameters. For this, let $A(z) = \{ {A_1}(z),{A_2}(z),\ldots ,{A_{{N_A}}}(z)\}$ be a set of ${N_A}$ averaged A-lines with ${N_D}\; $ data-points on each A-line. For a matrix of ${N_A} \times {N_D}$ averaged A-lines, any unknown parameter among *C, µ*, *z _{0}*,

*z*$\; $ can be estimated by minimizing the sum of squared residuals, $\chi ,$ given by,

_{R}## 3. Results

In this section, we first present the statistical analysis and numerical simulations to study the performance of the different models in Table 1 and estimation method using Cramér-Rao analysis and Monte Carlo simulation, respectively. This provides insight into the available information embedded in the data for different experimental setups and models with different degrees of freedom. Finally, we present the experimental results on a homogeneous phantom to assess the real-life performance of the proposed method.

#### 3.1 Model selection by Cramér-Rao analysis

A Cramér-Rao analysis was performed to assess the amount of information present in the data and the impact thereof on the attainable precision for all model parameters. Equation (1) was used to simulate OCT depth profiles. A simulated (thick) homogeneous sample with a refractive index of 1.44, a physical thickness of 1 mm and an attenuation coefficient of 0.72$\; \textrm{m}{\textrm{m}^{ - 1}}$ was located at the zero-delay line. The model parameters in Eq. (1) were set to ${z_R}$ = 42 µm, *C* = 2.5${\times} $10^{4} and ${z_0}$ = 160 µm inside the sample. Each A-line consisted of 788 pixels and the physical axial pixel size $\Delta{z}$ of the system in air was set to 1.27 µm. For a realistic simulation, the OCT signal was distorted by exponential noise with the intensity-dependent mean at each depth, equal to the expected values of *S(z)* in Eq. (1). To reduce the noise, we averaged over single A-lines as explained in section 2.1. An example of a simulated single A-line is presented in Fig. 1(a), together with an averaged (over 500 simulated A-lines) OCT signal. The noise of the averaged A-line resembles an intensity-dependent Gaussian distribution due to the central limit theorem. We obtained the intensity-dependent standard deviations for all averaged A-lines using 1000 observations of the simulated averaged-A-lines. In Fig. 1(b), we show the rCRLB values after averaging *N* single A-lines. The diagonal elements of the rCRLB matrix represent the optimal precision of the estimated model parameters for the aforementioned intensity-dependent standard deviations. As is shown, by averaging over 500 A-lines, the estimation error of *µ* remains below 10%. Averaging of 500 A-lines is therefore used for the simulated and measured OCT signals in the following sections. The calculated rCRLB matrix for the simulated signals shown in Fig. 2. As seen in this figure, the model parameter estimation errors remain below 5%.

In addition, we investigate the precision of the estimated parameters for different degrees of freedom imposed on the model. The diagonal elements of the rCRLB matrix for the seven models shown in Table 1 are depicted in Fig. 3. As it can be observed, incorporating prior knowledge by fixing the parameter values on different combinations of ${z_0}$, *C*, and ${z_R}$ results in a better precision of parameter *µ* for depth-variant noise as indicated by the smaller values for Model 1…Model 6 compared to Model 7 as defined in Table 1.

#### 3.2 Experiment design by Cramér-Rao analysis

To assess the attainable precision under different experimental conditions, we calculated the rCRLB matrices as a function of one of the model parameters while keeping the other ones fixed. The rCRLB values are shown in Fig. 4 for a range of parameter values. In this figure, the horizontal dashed lines indicate the acceptable error (below 10%) and the vertical dashed lines indicate the set parameter values in the simulation and also were considered to be fixed for the other plots in this figure. In Fig. 4(a), it can be observed that when the focus location is above the sample, the estimation error for parameters ${z_o}$ and *µ* is exceeding 10%. In addition, when the focus location is less than 0.08 mm inside the sample, the estimation error for *µ* is within 10% to 40% and it remains below 10% for the deeper focus locations. The estimation errors for *C* and ${z_R}$ remain below 10% when the focus is located inside the sample. Figure 4(b) shows that for a Rayleigh length below 300 µm the estimation errors of *µ* remains below 10%. The estimation error of ${z_0}$ increases to above 10% for Rayleigh lengths larger than 250 µm. The estimation errors for *C* and *z _{R}* remain below 10% for Rayleigh lengths below 500 µm. By varying the attenuation of the sample, Fig. 4(c) shows that the estimation errors of

*µ*as well as ${z_0}\; $ remain below 10% for the attenuation coefficient values between 0.25 mm

^{−1}and 8 mm

^{−1}. Figure 4(d) shows that by increasing the light intensity the precision of all estimated parameters remains the same due to the intensity-dependent noise.

#### 3.3 Estimation accuracy and precision: Monte Carlo simulation

To investigate if the theoretical lower bounds on the precision estimated by CRLB can be attained by our estimation method, Monte Carlo simulations were performed. In the simulated data, the location of the focus was varied between the surface and 0.6 mm inside the sample; the other model parameters were set to ${z_R}$ = 42 µm, *C* = 2.5${\times} $10^{4} and *µ* = 0.72 mm^{−1}. Next, we simulated 500 averaged A-lines distorted by Gaussian noise with an intensity-dependent standard deviation, as explained in section 3.1, for different focus locations. The method in section 2.1 was applied to estimate the model parameters using the *fmincon* function of MATLAB [Curve Fitting Toolbox, MATLAB 2013; The MathWorks, Natick, MA] using interior-point optimization with a termination tolerance set to ${10^{ - 15}}$, and the maximum number of iteration and function evaluations set to ${10^5}$.

Prior knowledge of the sample under investigation in combination with known properties of the optical system are useful to set suitable initial parameter values. The initial value of *C* (2.5${\times} $10^{4} (arb. unit)) was set by choosing an arbitrary A-line and taking the average of the intensity values at all depths within the sample. To investigate the effect of the initial parameter values on the estimation results, the initial parameter values were varied individually, over the following ranges: 0.01 mm ≤$\; \; {z_R}\; $≤ 0.2 mm, 0 mm ≤ ${z_0} \le $ 2 mm, 10^{3} ≤$\; C\; $≤7${\times} $10^{4} and 0 mm^{−1}≤ *µ* ≤ 6 mm^{−1}, while the other parameters were set to the aforementioned initial parameter values. Figure 5 shows the CoV and bias of the estimated parameters for different settings of the initial values. As can be seen, the CoVs remain below 10% and bias error below 1%.

The coefficient of variation (CoV) of the estimated parameters in Monte Carlo simulation, the rCRLB values for different parameters and the estimation bias as a function of focus location, are shown in Fig. 6. For varying ${z_0}$, the initial values for the unknown parameters were set to ${z_R} = \; $50 µm, *C* = 2$\; \times $10^{4} (arb. unit), *µ* = 1 mm^{−1} and the values of $\; {z_0}$ were set to 0.2 mm above the expected focus locations for each averaged A-line. We can observe in Fig. 6(a) that the estimation error of the parameters by Monte Carlo simulation is below 12% when the focus location is inside the sample. In addition, Fig. 6(b) shows an acceptable bias error in the Monte Carlo simulation.

To summarize the results, it has been shown that knowing more model parameters results in a better estimation precision of *µ* (Fig. 3). Additionally, the proposed approach ideally can estimate the model parameters with an acceptable precision (below 10%) when the number of averaged A-lines is larger than 500, the location of the focus is inside the sample, the Rayleigh length is below 0.4 mm and the attenuation coefficient of the sample is more than 0.2 mm^{−1} and less than 6 mm^{−1} (Fig. 4). Therefore, these limitations for Rayleigh length and the attenuation coefficient were considered in designing our experimental setup. In addition, we obtained acceptable results using interior-point solver. This routine was used for estimating the model parameters in the real measurements explained in the next section.

#### 3.4 Experimental setup

In this section, we apply our method to the measurements obtained from different samples with an experimental OCT system. We apply them to estimate the model parameters based on either single or multiple B-scans. The B-scans of three thick or semi-infinite samples with 0.05 wt. %, 0.1 wt. % and 0.25 wt. % of TiO_{2} in silicone, with the zero delay location positioned 0.4 mm above the sample surface, are recorded with various locations of the focal plane from the samples’ surfaces using a Ganymede-II-HR Thorlabs spectral domain OCT system (GAN905HV2-BU) [24]. The system has a centre wavelength of 900 nm and a bandwidth of 195 nm and a Thorlabs scan lens (LSM02-BB) with 18 mm focal length. The system’s axial and lateral resolutions were 3 µm (in air) and 4 µm, respectively, and the axial and lateral physical pixel size in air was 1.2$7\; \times $ 2.9 µm with 1024 pixels on each A-line.

First, the focus position was manually set to the sample’s surface by optimizing the surface structure’s sharpness in the centre of enface image created by the OCT camera. Next, 90 B-scans were obtained at various locations of the focal plane by changing the axial location of the lens in the sample arm with a physical step size of 11.25 µm within a range of ${\pm} $ 0.5 mm around the initial focus location. Figure 7(a) shows a B-scan of 0.05 wt. % TiO_{2} in silicone with adjusted focus location at 0.26 mm from the focus location on the surface, as estimated by optimizing the sharpness. Figure 7(b) shows the concatenation of the averaged A-lines (from each B-scan) as a function of focus position. The location of the sample remained fixed in the B-scan for various focus positions by adjusting the optical path length of the system’s reference arm. As can be seen in Fig. 7(b), the highest intensity on the surface deviates from the centre of the image, which indicates a shift in the aforementioned adjustment of the focus position on the surface.

Several pre-processing steps have been performed on the measured spectra. First, the reference arm intensity, which was measured automatically by the system for every acquisition, was removed. Second, to compensate for the nonlinear spacing, the spectra were interpolated based on the reported wavelength distribution on the detector by the system manufacturer. Finally, the OCT signals were generated by squaring the Fourier transform of the interpolated signal. As mentioned in section 2.1, the averaged noise floor was obtained, by averaging over a large number of A-lines while having no sample in the sample arm, and subtracted from the A-lines. Afterwards, the roll-off of the system was measured and the A-lines were corrected for it. To estimate the roll-off, scans of 0.25 wt. % TiO_{2} for different locations of the surface within the OCT image depth range were obtained by changing the optical path length of the reference arm, and the sensitivity decay was obtained by smoothing and interpolating the average intensities of the corresponding region of interest in the scans. Finally, regions above and inside the sample that only contained noise were removed for each B-scan [7].

The system’s Rayleigh length was estimated by the following procedure. For each B-scan obtained from the sample with 0.05 wt. % of TiO_{2} in silicone, the averaged A-line was calculated and the arbitrary data-point at the physical depth of 63 µm from the sample’s surface were recorded. This physical depth should be close enough to the surface to obtain a sufficiently high SNR and far enough from the surface to avoid the data-points being affected by reflection artefact. The initial values of ${z_o}$ and ${z_R}$ were obtained by fitting the following model to the recorded data-points (Fig. 7(c)),

*n*of the medium [18]. In addition, the shifted focus positions were transformed from physical to optical distance. The optical Rayleigh lengths in air and silicone (with refractive index of ${n_{sample}} = \; $1.44) were estimated to be 29.3 µm and 60.8 µm (${z_R}_{Si} = {n_{Si}}^2.{z_R}_{air}$), respectively, shown in Fig. 7(c). The physical Rayleigh length in silicone was calculated to be 42 µm.

#### 3.5 Estimating the model parameters in a single B-Scan

The unconstrained model of Eq. (1) was fit to the averaged A-lines obtained from 500 single A-lines of each B-scan. The initial values were set to $\mathrm{\mu} = \; $3 mm^{−1}, *C* = 5 ${\times} $ 10^{3}, ${z_R}$ = 56 µm, and the values of $\; {z_0}$ were set to 0.2 mm above the expected focus locations for each Averaged A-line. ${z_0} =$ (the expected focus location $- $ 0.2 mm). The results of the fit to the measurements of 8 B-scans acquired with focus positions of {−0.5, −0.3, −0.1, 0.15, 0.3, 0.45, 0.6, 0.75} mm are shown in Fig. 8. The estimated parameter values for all recorded B-scans are shown in Fig. 9. Since parameters ${z_R}$ and *µ* are constant among the B-scans, we expect them to have similar estimated values. However, as can be seen, when the focus is above the sample, the estimated parameters are far from the expected values. Additionally, when the focus location is within the depth of 0.1 mm inside the sample, the estimated attenuation coefficients seems to be significantly different compare to their estimated values at larger depths. This also confirms the simulation results in section 3.1.

For a location of focus inside the sample, the estimated focus locations closer to the surface are in better agreement with the expected focus locations than for focus positions very deep into the sample.

#### 3.6 Estimating the model parameters using multiple B-Scans

Multiple B-scans of the same sample acquired with different focus positions can also be combined to estimate the model parameters as was explained in section 2.2. In this case, the model parameters $\mathrm{\mu} $ and ${z_R}$ were considered common while the focus position ${z_o}$ and *C* were left free to vary among the B-scans. Only the B-scans in which the focus was inside the sample were considered in this experiment. The initial parameter values of the fitting process were the same as the values mentioned in section 3.3. To investigate if the estimation result depends on the number of the B-scans used, different numbers of the B- scans (2, 4, 8, 16 and 32) acquired from the sample with 0.05 wt. % TiO_{2} in silicone were used. The estimated parameters *µ*$,\; {z_R}$ and ${z_0}$ are shown in Fig. 10 for different numbers of combined B-Scans. As can be seen, the estimated attenuation coefficient and Rayleigh length do not change significantly when more than 8 B-scans were used.

We combined 8 B-scans in estimating the attenuation coefficient $\mathrm{\mu} $ for samples with different concentrations of TiO_{2} in silicone (0.05 wt. %, 0.1 wt. % and 0.25 wt. %). The resuts are shown in Table 2. The standard errors for the estimated attenuation coefficients, Rayleigh lengths and the focus locations were 0.01 mm^{−1}, 0.0001 mm, and less than 0.01 mm, respectively. In theory, the relationship between the particle concentration and the attenuation coefficient should be linear. The estimation result shows a reasonable correlation between the TiO_{2} weight concentration and the estimated attenuation coefficient [with *R ^{2 }*= 0.99 calculated over all B-scans].

## 4. Discussion

We presented a method to estimate the sample attenuation coefficients from OCT measurements compensated for the effects of the axial PSF. The recorded signals were modeled by assuming single-scattering in a homogeneous sample accounting for the system’s roll-off, noise and focused beam shape (axial PSF). The model parameters of the focused axial PSF were estimated from experimental OCT data. Our goal was to achieve accurate estimation of the attenuation coefficient to enable reliable quantitative analysis of a sample under investigation. The numerical study predicted the performance, and hence the limitations, of our model for different experimental conditions. We observed that for a Rayleigh length smaller than 0.3 mm the estimation error of attenuation coefficient is smaller than 10% for a sample with attenuation coefficient of 0.73 mm^{−1} and with the location of the focused beam inside the sample. The signal decay caused by the effect of axial PSF for the location of the focus inside and close to the sample’s surface for larger Rayleigh lengths is not as significant as for smaller Rayleigh lengths. Accurate estimation of the attenuation coefficient for the location of focus above the sample is not feasible since it is nearly impossible to distinguish between the decay of the recorded signal caused by attenuation and by the axial PSF.

In experiments with phantoms of different weight concentrations TiO_{2} in silicone, good fits of the model to the measurements using a single B-scan as well as multiple B-scans acquired for different focus positions were obtained. In the latter case the parameters related to sample and optics were shared, whereas the focus position and signal strength were allowed to vary among the different B-scans. Only for focus positions very deep into the sample, the fitted model parameters started to deviate from the measurements closer to the sample’s surface. This might be due to the background noise subtraction or an incorrect estimation of roll-off.

The estimated attenuation coefficients using single B-scans are varying among the B-scans and tend to decrease when the focus is shifting to larger depths. For the focus locations above the sample and close to the surface inside the sample, the estimated attenuation coefficients are significantly different compared to the estimated values at larger depths. This confirms the numerical results depicted in Fig. 4(a). We observed that the incorrect estimation of model parameter *C* can significantly influence the estimation of the attenuation coefficients. Therefore, finding a method to fix this parameter would significantly improve the results of estimating the attenuation coefficient from a single B-scan.

We also combined multiple B-scans acquired at different focus positions to estimate the model parameters. The estimated Rayleigh length and locations of focus in different B-scans could be estimated with less than 3% error while using 8 B-scans. Using more B-scans does not show a significant improvement in estimating the model parameters. To be able to compare the numerical and experimental results properly, the true attenuation coefficient values of the samples should be known. While these attenuation coefficients are unknown, we do know the concentration of TiO_{2}. We showed that in samples with different concentrations of TiO_{2} dispersed in silicone there is a linear relation between the 0.05 wt. % and 0.1 wt. % TiO_{2} weight concentration and the estimated attenuation coefficients. However, the estimated attenuation coefficient for 0.25 wt. % TiO_{2} is slightly lower than the expected value. It has been shown previously that the measured attenuation coefficient falls below the expected values for increasing particles concentration due to an increase of the amount of multiple scattering [25]. Applying a multiple-scattering model can result in a better correlation between the measurements and the OCT light model.

A limitation of this work is the assumption of isotropic scattering and weak concentrations of scatterers such that a single-scattering model suffices.

In ophthalmology we encounter a shift of the focal plane due to accommodation of the human eye and movements of the eye and the head. The method can also be applied to data in which the focus position remains fixed for all B-scans. In addition, in ophthalmology, the recorded OCT scans obtained from a clinical OCT system are usually averaged over 10-100 B-scans during acquisition to improve the image quality of the scans. To obtain the precision explained in this research (by averaging over 500 A-lines), we only need to average over a small area in the image (5-50 neighboring A-lines). To obtain a higher spatial resolution, it is advised to reduce the number of averaging among neighboring A-lines while increasing the averaging rate during acquisition. In future work, we will extend the proposed method to estimate the attenuation coefficient values of a multi-layer sample from the OCT images. This is especially relevant in applications such as ophthalmology where each retina layer has different optical properties.

## Funding

ZonMw (91212061).

## Acknowledgment

This research was funded by the Netherlands Organization for Health Research and Development (ZonMw) TOP grant (91212061). We would like to acknowledge Dr. Dirk H. J. Poot (Biomedical Imaging Group Rotterdam, departments of Medical Informatics and Radiology, Erasmus Medical Center Rotterdam, Rotterdam, The Netherlands) for his technical support in this work.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **K. A. Vermeer, J. van der Schoot, H. G. Lemij, and J. F. de Boer, “RPE-normalized RNFL attenuation coefficient maps derived from volumetric OCT imaging for glaucoma assessment,” Invest. Ophthalmol. Visual Sci. **53**(10), 6102–6108 (2012). [CrossRef]

**2. **R. A. McLaughlin, L. Scolaro, P. Robbins, C. Saunders, S. L. Jacques, and D. D. Sampson, “Parametric imaging of cancer with optical coherence tomography,” J. Biomed. Opt. **15**(4), 046029 (2010). [CrossRef]

**3. **K. Barwari, D. M. de Bruin, E. C. Cauberg, D. J. Faber, T. G. van Leeuwen, H. Wijkstra, J. de la Rosette, and M. P. Laguna, “Advanced diagnostics in renal mass using optical coherence tomography: a preliminary report,” J. Endourol. **25**(2), 311–315 (2011). [CrossRef]

**4. **K. Barwari, D. M. de Bruin, D. J. Faber, T. G. van Leeuwen, J. J. de la Rosette, and M. P. Laguna, “Differentiation between normal renal tissue and renal tumors using functional optical coherence tomography: a phase I in vivo human study,” BJU Int. **110**(8b), E415–E420 (2012). [CrossRef]

**5. **P. H. Tomlins, O. Adegun, E. Hagi-Pavli, K. Piper, D. Bader, and F. Fortune, “Scattering attenuation microscopy of oral epithelial dysplasia,” J. Biomed. Opt. **15**(6), 066003 (2010). [CrossRef]

**6. **Q. Q. Zhang, X. J. Wu, T. Tang, S. W. Zhu, Q. Yao, B. Z. Gao, and X. C. Yuan, “Quantitative analysis of rectal cancer by spectral domain optical coherence tomography,” Phys. Med. Biol. **57**(16), 5235–5244 (2012). [CrossRef]

**7. **F. J. van der Meer, D. J. Faber, D. M. B. Sassoon, M. C. Aalders, G. Pasterkamp, and T. G. van Leeuwen, “Localized measurement of optical attenuation coefficients of atherosclerotic plaque constituents by quantitative optical coherence tomography,” IEEE Trans. Med. Imag. **24**(10), 1369–1376 (2005). [CrossRef]

**8. **C. Xu, J. M. Schmitt, S. G. Carlier, and R. Virmani, “Characterization of atherosclerosis plaques by measuring both backscattering and attenuation coefficients in optical coherence tomography,” J. Biomed. Opt. **13**(3), 034003 (2008). [CrossRef]

**9. **G. van Soest, T. Goderie, E. Regar, S. Koljenovic, G. L. J. H. van Leenders, N. Gonzalo, S. van Noorden, T. Oka-mura, B. E. Bouma, G. J. Tearney, J. W. Oosterhuis, P. W. Serruys, and A. F. W. van der Steen, “Atherosclerotic tissue characterization in vivo by optical coherence tomography attenuation imaging,” J. Biomed. Opt. **15**(1), 011105 (2010). [CrossRef]

**10. **J. M. Schmitt, A. Knüttel, M. Yadlowsky, and M. A. Eckhaus, “Optical-coherence tomography of a dense tissue: statistics of attenuation and backscattering,” Phys. Med. Biol. **39**(10), 1705–1720 (1994). [CrossRef]

**11. **R. O. Esenaliev, K. V. Larin, I. V. Larina, and M. Motamedi, “Non-invasive monitoring of glucose concentration with optical coherence tomography,” Opt. Lett. **26**(13), 992–994 (2001). [CrossRef]

**12. **A. I. Kholodnykh, I. Y. Petrova, K. V. Larin, M. Motamedi, and R. O. Esenaliev, “Precision of Measurement of Tissue Optical Properties with Optical Coherence Tomography,” Appl. Opt. **42**(16), 3027–3037 (2003). [CrossRef]

**13. **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**(3), 484–490 (2000). [CrossRef]

**14. **V. D. Nguyen, D. J. Faber, E. van der Pol, T. G. van Leeuwen, and J. Kalkman, “Dependent and multiple scattering in transmission and backscattering optical coherence tomography,” Opt. Express **21**(24), 29145–29156 (2013). [CrossRef]

**15. **K. A. Vermeer, J. Mo, J. J. Weda, H. G. Lemij, and J. F. de Boer, “Depth-resolved model-based reconstruction of attenuation coefficients in optical coherence tomography,” Biomed. Opt. Express **5**(1), 322–337 (2014). [CrossRef]

**16. **B. Ghafaryasl, K. A. Vermeer, J. F. de Boer, M. E. J. van Velthoven, and L. J. van Vliet, “Noise-adaptive attenuation coefficient estimation in spectral domain optical coherence tomography data,” in Proceedings of IEEE 13th International Symposium on Biomedical Imaging (ISBI), 706–709 (2016).

**17. **S. H. Yun, G. J. Tearney, B. E. Bouma, B. H. Park, and J. F. de Boer, “High-speed spectral-domain optical coherence tomography at 1.3 mu m wavelength,” Opt. Express **11**(26), 3598–3604 (2003). [CrossRef]

**18. **T. G. van Leeuwen, D. J. Faber, and M. C. Aalders, “Measurement of the axial point spread function in scattering media using single-mode fiber-based optical coherence tomography,” IEEE J. Select. Topics Quantum Electron. **9**(2), 227–233 (2003). [CrossRef]

**19. **B. Ghafaryasl, K. A. Vermeer, J. Kalkman, T. Callewaert, J. F. de Boer, and L. J. van Vliet, “Accurate estimation of the attenuation coefficient from axial point spread function corrected OCT scans of a single layer phantom,” Proc. SPIE **10483**, 104832B (2018). [CrossRef]

**20. **G. T. Smith, N. Dwork, D. O’Connor, U. Sikora, K. L. Lurie, J. M. Pauly, and A. K. Ellerbee, “Automated, depth-resolved estimation of the attenuation coefficient from optical coherence tomography data,” IEEE Trans. Med. Imaging **34**(12), 2592–2602 (2015). [CrossRef]

**21. **S. Stefan, K. S. Jeong, C. Polucha, N. Tapinos, S. A. Toms, and J. Lee, “Determination of confocal profile and curved focal plane for OCT mapping of the attenuation coefficient,” Biomed. Opt. Express **9**(10), 5084–5099 (2018). [CrossRef]

**22. **D. J. Faber, F. J. van der Meer, M. C. G. Aalders, and T. G. van Leeuwen, “Quantitative measurement of attenuation coefficients of weakly scattering media using optical coherence tomography,” Opt. Express **12**(19), 4353–4365 (2004). [CrossRef]

**23. **M. W. A. Caan, H. G. Khedoe, D. H. J. Poot, A. J. den Dekker, S. D. Olabarriaga, K. A. Grimbergen, L. J. van Vliet, and F. M. Vos, “Estimation of diffusion properties in crossing fiber bundles,” IEEE Trans. Med. Imaging **29**(8), 1504–1515 (2010). [CrossRef]

**24. **T. Callewaert, J. Dik, and J. Kalkman, “Segmentation of thin corrugated layers in high-resolution OCT images,” Opt. Express **25**(26), 32816–32828 (2017). [CrossRef]

**25. **J. Kalkman, A. V. Bykov, D. J. Faber, and T. G. van Leeuwen, “Multiple and dependent scattering effects in Doppler optical coherence tomography,” Opt. Express **18**(4), 3883–3892 (2010). [CrossRef]