## Abstract

Optical properties, such as the attenuation coefficients of multi-layer tissue samples, could be used as a biomarker for diagnosis and disease progression in clinical practice. In this paper, we present a method to estimate the attenuation coefficients in a multi-layer sample by fitting a single scattering model for the OCT signal to the recorded OCT signal. In addition, we employ numerical simulations to obtain the theoretically achievable precision and accuracy of the estimated parameters under various experimental conditions. Finally, the method is applied to two sets of measurements obtained from a multi-layer phantom by two experimental OCT systems: one with a large and one with a small Rayleigh length. Numerical and experimental results show an accurate estimation of the attenuation coefficients when using multiple B-scans.

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

## 1. Introduction

Optical coherence tomography (OCT) is an interferometric imaging technique that can generate high-resolution three-dimensional images of biological tissues. Many tissues, such as the retina or the blood vessel wall are composed of a number of distinct tissues, each having its own optical tissue properties. OCT has been widely used to capture structural information of tissues for clinical tasks such as the diagnosis of retinal and vascular diseases. Recently, there has been a growing interest in the field of ophthalmology for utilizing optical tissue properties, such as the attenuation coefficient, for diagnosis and disease progression. The attenuation coefficient can be estimated from the intensity or amplitude of the OCT signal and can be used as a biomarker for the diagnosis and monitoring of diseases [1–6] as well as for tissue characterization [7–9]. Several methods based on single [10–12] and multiple [13,14] scattering of light have been presented for estimating the attenuation coefficient in a homogeneous medium using OCT. Only a few methods have been developed to estimate the attenuation coefficient of the tissue while taking into account the effect of the beam shape on the acquired OCT signal. Smith et al. [15] compensated for the effect of beam shape by correcting the OCT scan for the confocal detection efficiency using an existing model [8]. However, in their work the parameters of the shape of the beam need to be known before the estimation of the attenuation confident. In many medical applications, such as ophthalmology, the location of the focus varies between B-scans 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. [16] introduced a method to estimate the attenuation coefficient using two B-scans with different focus locations to first estimate the location of focus and subsequently estimate the attenuation coefficient from a single scattering model of the OCT light after compensating for the effect of beam shape. However, this method is dependent on having identical lateral beam locations of the sample to unambiguously determine the effect of the beam. In addition, this method was only tested with static samples where an identical physical location in both B-scans is well feasible and factors such as the beam angle of incidence can be controlled to ensure a similar intensity of the measured OCT signal. The limitation of having access to two scans from exactly the same position in the tissue is a problem in many clinical data, such as retinal scans, where only one averaged measurement of the same tissue’s location is available. Recently, we presented a method to achieve an accurate estimate of attenuation coefficient for semi-infinite samples by compensating for the effects of the beam shape [17]. The proposed method estimates the axial point spread function (PSF) model parameters (Rayleigh length and focus depth) and attenuation coefficient by fitting a single-scattering model to the measured OCT signal of a homogenous sample. Monte Carlo simulations quantified the maximum expected accuracy and precision of our proposed method. However, while the method could estimate the attenuation coefficients of the materials from measurements of uniform samples, most biological tissues such as the retina are layered, and hence, the method cannot be applied straightforwardly. Therefore, there is a need, e.g. in ophthalmology, for a method to reliably estimate the attenuation coefficient properties in multi-layer samples.

In this paper, we demonstrate simultaneous attenuation coefficient and beam focus position estimation in multi-layer samples. We investigate whether using multiple B-scans, acquired with different focus positions, improves the estimation of the attenuation coefficient. In addition, simulations of the OCT signal with speckles were used to numerically evaluate the accuracy and precision of the estimated parameters by the proposed method. The numerical analysis is needed for the evaluation of the proposed method since the true attenuation coefficient values of the materials are not known in real experiments. Finally, actual measurements of multi-layer phantoms composed of layers with different concentrations of TiO_{2} dispersed in silicone is used to investigate the accuracy of our method. For this, the method is applied to B-scans of this phantom obtained with two different experimental OCT systems: one with a small and one with a large Rayleigh length.

## 2. Method and experiments

In this section, the OCT signals for multi-layer samples including the effect of the axial point spread function is described. In addition, we introduce a method for accurate estimation of the attenuation coefficient and focus position in a multi-layer sample after compensating the recorded signal for the average noise level and roll-off. We simulate a realistic OCT signal including speckle formation and intensity noise. A Monte Carlo simulation is used to investigate the accuracy and precision of the AC estimation method. Finally, we present the experimental results on a multi-layer phantom to assess the performance of the proposed methods in practice.

#### 2.1 Estimating the model parameters

We simultaneously estimate the model parameters of the axial PSF and the attenuation coefficient per layer using a maximum likelihood estimator. For an inhomogeneous sample, a single-scattering physical model of averaged intensity of the OCT light at sampled physical depth *z* can be expressed as,

*R*(

*z*) is the signal intensity decay caused by roll-off as the ratio of spectral resolution over pixel resolution [18]. In the first term the signal variation is caused by three factors (from left to right): a scaling factor

*L,*which comprises several factors such as input power, detector efficiency, coherence length and integrated phase function; the signal decay caused by a depth-resolved exponential signal attenuation modelled by ${\mu _t}(z )$ (the total attenuation coefficient), which is the sum of ${\mu _a}(z )$ (absorption coefficient) and ${\mu _s}(z )$ (the scattering coefficient). The backscattering coefficient is $\; {\mu _{b,NA}}(z )$ = ${\mu _s}(z ).{p_{NA}}(z )$, with ${p_{NA}}(z )$ being the fraction of scattered photons detected by the OCT system [11,19]. Finally, the axial PSF is modelled by a Lorentzian function at focus position ${z_0}$ and scaled by the Rayleigh length in the medium, ${z_R}(z )$ [20]. The second term$\; D(z )$ is the noise level offset. The intensity of the OCT signal inside a (homogenous) layer is distorted by speckle and has an exponential distribution [18]. However, based on the central limit theorem, the uncorrelated intensity noise in the average of a large enough number (

*M*) of neighbouring A-lines (

*M*> 30 based on rule of thumb) converges to a normal distribution ${{\cal N}}[{m(z ),M} )]$, with

*m*(

*z*) being the expected value of the exponential distribution, and $m{(z )^2}/M$ the variance of the resulting normal distribution. For an averaged A-line, the third term

*ɛ*(

*z*) is a Gaussian noise which represents the effect of additive background noise and speckles with normal distributions. The acquired OCT signal can be corrected for the effect of noise floor [15] and roll-off [15,16].

To reduce the complexity of the fit model, we consider ${\mu _t}$ = ${\mu _s}$ and call it *µ*, by assuming ${\mu _a}$ ${\ll} \; {\mu _s}$, which is a good approximation for the wavelengths typically used for ophthalmic OCT. In addition, we assumed ${p_{NA}}(z )$ and refractive indices to be constant for all depths. A constant refractive index results in a constant Rayleigh length for all depths. With the aforementioned assumptions, a model of the OCT signal with ${N_D}\; $data-points ${z_j}$ on each A-line can be expressed as,

*C =*$L\cdot {p_{NA}}$ is a scaling factor. To estimate the model parameters, we use multiple averaged A-lines with different focus locations ${z_{{o_i}}}$. Therefore, for a set of ${N_A}$ averaged A-lines ${I_1}({{z_j}} ),\; {I_2}({{z_j}} ),\; \ldots ,\; {I_N}_A({{z_j}} )$, the unknown parameters

*C*and$\; {z_{{o_i}}}$, and the attenuation coefficient at each depth, can be estimated by minimizing the log-weighted sum of squared residuals, $\chi ,$ given by,

_{i}*j*is an index that denotes the data-point number on each A-line. The proposed method requires knowledge of the thickness of each layer per A-line and therefore segmentation of the multi-layer sample is a prerequisite.

#### 2.2. Simulation of OCT signal and Monte Carlo simulation

To study the performance of the proposed method numerically, we applied it to simulated OCT signals. For this, we integrated the effect of axial PSF to an existing single-scattering based simulation of the OCT signals, which are distorted by speckle and signal intensity noise of a multi-layer sample [21].

OCT signals were simulated for a system with a Gaussian-shaped spectrum with a center wavelength of 1000 nm and a full width at half maximum of 73 nm. The front surface is located at 0.2 mm from the zero-delay. The simulated sample has four homogeneous layers each having thickness (*d*) and attenuation coefficients (*µ*) of: ${d_1}$ = 170 µm, ${d_2}$ = ${d_3}$ = ${d_4}$ = 100 µm and *µ _{1}*$\; = $ 4.4 mm

^{-1},

*µ*=

_{2 }*µ*/2.5 = 1.76 mm

_{1}^{-1},

*µ*=

_{3 }*µ*/5 = 0.88 mm

_{1}^{-1}and

*µ*= ${\mu _1}$/2.5 = 1.76 mm

_{4}^{-1}, respectively. To simulate these attenuation coefficient values, the fraction of the scattered intensity was assumed to be 0.5 and the averaged particle size was set to 700 nm. The concentration of particles in each layer was

*P*= 7%,

_{1}*P*= 2.8%,

_{2}*P*= 1.4%,

_{3}*P*= 2.8%. The model is based on the single scattering approximation. In addition, to calculate the scattering properties of the particles using Mie-scattering, the refractive indices of the sample and particles were set to$\; {n_{med}}$ = 1.44 and$\; {n_{part}}$ = 1.48. The Rayleigh length of the axial PSF’s model was set to either ${z_R}$ = 432 µm or ${z_R}$ = 57.6 µm in the medium, and the location of the focus varied between ${\pm} $1 mm from the surface of the sample with a step size of 20 µm. We simulated 500 averaged A-lines, with 1024 data-points per A-line, distorted by additive Gaussian noise with a standard deviation of 5% of the signal intensity. All the simulations were implemented in Matlab 2017b on a Dell Latitude with a dual core CPU (2.60 GHz) and 8 GB of RAM.

_{4}For each Rayleigh length, the method in section 2.1 was applied to estimate the model parameters *C, z _{0}* and

*µ*(

_{n}*n*= 1,2,3,4), while fixing the

*z*values to the known Rayleigh lengths, using the interior–point technique of Matlab (Curve Fitting Toolbox, Matlab 2017b; The MathWorks, Natick, MA) with a termination tolerance set to ${10^{ - 15}}$, and the maximum number of iteration and function evaluations set to ${10^6}$. The initial values for the unknown parameters were set to

_{R}*C*= 5$.$10

^{4}, ${z_0} = \; $0.6 mm and ${\mu _1}$ $= $ 3 mm

^{-1}, ${\mu _2} = $ 3 mm

^{-1}, ${\mu _3} = \; $2 mm

^{-1}and ${\mu _4} = $ 1 mm

^{-1}. In the process of fitting the model of Eq. (2) for homogeneous layers, a lower and upper bound was set for all parameters, i.e. 0 mm

^{-1}$< \mu < $ 6 mm

^{-1}for all layers, 0 $< C < $ 10

^{10}(arbitrary units), 0 mm $< {z_0} < \; $1 mm for focus inside and -1 mm$\; < {z_0} < \; $0 mm for focus above the sample.

To evaluate the accuracy, precision and feasibility of the proposed method numerically, the coefficient of variations (CoV) of the estimated parameters and the bias were calculated for the two selected Rayleigh lengths and different locations of focus around the surface of the sample.

#### 2.3. Phantom experiments

To investigate the practical feasibility of the proposed method, we applied it to measurements obtained from a multi-layer phantom by two experimental OCT system with different Rayleigh lengths. The model parameters were estimated based on either a single or multiple B-scans. The sample consists of four layers with 0.25 wt%, 0.1 wt%, 0.05 wt% and 0.1 wt% of TiO_{2} in silicone. The B-scans of the phantom were recorded for various locations of the focal plane from the sample’s surface using two experimental systems. The first one is a Ganymede-II-HR Thorlabs spectral domain OCT system (SD-OCT) (GAN905HV2-BU) with an estimated Rayleigh length of 36 µm in air. The Rayleigh length of the system was estimated by fitting the axial PSF to a set of measurements obtained by changing the focus location from a sample’s surface [19]. The system has a center wavelength of 900 nm and a bandwidth of 195 nm and a scan lens with 18 mm focal length (LSM02-BB). The axial and lateral physical pixel size in air was 1.27$\; \times $ 2.9 µm^{2} with 1024 pixels per A-line. The second system is a swept-source OCT (SS-OCT) system [22] with an estimated Rayleigh length of 288 µm in air. The system has a center wavelength of 1 µm. The axial and lateral pixel size in air was 3.3$\; \times $ 1.45 µm^{2} with 1024 pixels per A-line. The refractive index of silicone was considered to be 1.44 [11] and assumed to be constant through all layers.

First, the focus position was set inside the sample, but close to sample’s surface by probing the highest intensity in the area of interest during the acquisition. Next, 70 B-scans were obtained with various locations of the focal plane by changing the location of the lens in the sample arm with a physical step size of 20 µm and 15 µm for the systems with large and small Rayleigh length, respectively. Then, the B-scans with the location of focus within a range of ${\pm} $ 0.5 mm around the initial focus location were selected.

Several post-processing steps were performed on the measured A-lines. The averaged noise level was obtained by averaging over a large number of A-lines without any sample in the sample arm, and subtracted from all A-lines. In addition, the signal was corrected for the measured roll-off values of 0.81 and 1.7 for the two systems with large and small Rayleigh length, respectively [17]. Next, the surfaces of the samples were segmented in each B-scan using a minimum cost path search applied to individual B-scans [23]; the locations of the other interfaces were derived from the known thicknesses of the layers. Finally, the averaged A-lines were created by averaging over the central 200 A-lines in each B-scan.

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

The method in section 2.1 was applied to estimate the model parameters as explained in section 2.2. The model was fit to the averaged A-lines from each B-scan, with different locations of focus, where the model parameters *µ _{n}* (

*n =*1, 2, 3, 4),

*C*, and ${z_0}$ were unknown and Rayleigh lengths were set to the values of the experimental system. The initial values were set to ${\mu _1}$ $= \; $3.5 mm

^{-1},$\; {\mu _2}$ $= $ 3 mm

^{-1}, ${\mu _3} = $ 2 mm

^{-1}, ${\mu _4}$ = 3 mm

^{-1}, ${z_o}$ = 0.6 mm for both set of measurements, and

*C*= 10

^{3}(arbitrary units) for the set of measurements with large and

*C*= 10

^{7}for the set of measurements with small Rayleigh lengths. We introduced an upper and lower bound for each model parameters to restrict the optimizer to a reasonable domain. The intervals were set to 0 $\textrm{m}{\textrm{m}^{ - 1}} < \mu < $ 7$\; \textrm{m}{\textrm{m}^{ - 1}}$; and $0 \le {z_0} < \; $1 mm for the location of focus inside and -1 mm$\; < {z_0} < 0\; $ above the sample for both sets of measurements, and 0 $< C < $ 10

^{4}for the set of measurements with large and 0 $< C < $ 10

^{10}for the set with small Rayleigh lengths. Since the real values of the attenuation coefficients of the sample’s layers are unknown, the correlation between the estimated attenuation coefficients and the concentrations of TiO

_{2}in silicone was investigated for each set of measurements.

#### 2.5 Estimating the model parameters using multiple B-Scan

Multiple B-scans with different locations of focus were used to estimate the model parameters. The free running parameters *µ*_{1} through *µ*_{4} were considered to be common for all B-scans, while ${z_{{o_i}}}$ and ${C_{i\; }}$ in Eq. (3) were estimated for each B-scan individually. The Rayleigh length ${z_R}$ was fixed to the known values of the experimental setups. We used the same initialization values and boundary conditions for the optimization process as reported earlier for the single B-scan experiments. In addition, different combinations of the averaged A-lines were used to investigate the possibility to improve the results.

Previously, we showed that in samples with different concentrations of TiO_{2} dispersed in silicone there is a linear relation between the TiO_{2} weight concentration and the estimated attenuation coefficients [17]. The linear relation was calculated by fitting a line to the estimated attenuation coefficients in the averaged A-lines using Matlab’s *linear regression* (Statistics Toolbox, MATLAB 2017; The MathWorks, Natick, MA).

## 3. Results

#### 3.1 Simulation results

The proposed method in section 2.1 was tested with the simulation of the OCT signal as described in section 2.2. The OCT images of a multi-layer sample and the averaged central 200 central A-lines are shown in Fig. 1 for $\; {z_R}$ = 40 µm and ${z_R}$ = 300 µm in air and an initial location of focus ${z_0}$ = 0.1 mm inside the sample. The proposed method was applied to 200 averaged simulated A-lines and the unknown model parameters *C _{,} z_{0},* and ${\mu _n}$ (

*n =*1, 2, 3, 4) were estimated for different locations of the focus and the two selected Rayleigh lengths.

As can be observed in Fig. 2(a), the focus is estimated accurately for the small Rayleigh length when it is located inside the sample. However, Fig. 2(b) shows that for the large Rayleigh length, the estimation of the focus location was inaccurate because the effect of the focus on the intensity along an A-line cannot be distinguished from the signal decay caused by light attenuation. The CoV of the estimated attenuation coefficients of four layers and the corresponding estimation bias are shown in Fig. 2(c-f) as a function of focus location for the two selected Rayleigh lengths. As can be seen, the CoV for both systems remains below 10% when focussing inside the sample. For ${z_R}$ = 40 µm (small Rayleigh length), the bias of the estimated attenuation coefficients remains below 10% when the focus location is less than 0.5 mm inside the sample. For ${z_R}$ = 300 µm (large Rayleigh length), the bias of the estimated attenuation coefficient is above 10% for ${\mu _1}$ and ${\mu _2}$ and below 10% for ${\mu _3}$ and ${\mu _4}\; $if the location of the focus lies inside the sample. This bias is the result of the incorrect estimation of the focus.

The effect of the beam shape on the acquired OCT signal and therefore on the estimation of the attenuation coefficients is more significant for the smaller Rayleigh length. Therefore, for small Rayleigh lengths, an inaccurate estimation of the location of focus results in a less accurate estimation of the attenuation coefficients. This effect is shown in Fig. 2(a-d) where erroneous estimations of focus locations increase the CoVs for the smaller Rayleigh length, while there is no significant change of the CoVs for the larger Rayleigh length.

The linear relation between the estimated attenuation coefficients and the particle concentration of the respective layer as a function of focus location is shown in Fig. 3 for both the small and the large Rayleigh length. This figure shows *R ^{2}*-values indicating the goodness of the linear fit and the corresponding p-values. The

*R*-values larger than 0.95 and 0.98 for small and large Rayleigh lengths, respectively, show a good correlation between the estimated attenuation coefficients with

^{2}*p*< 0.01 and the TiO

_{2}concentrations for all focus positions.

In addition, Fig. 4 shows this linear relation between the estimated attenuation coefficients and the particle concentration of the respective layer in the averaged A-lines acquired with the focus set to 0.1 mm inside the sample as a function of the system’s Rayleigh length from 0.005 mm to 1 mm. This figure shows *R ^{2}*-values indicating the goodness of the linear fit and the corresponding p-values. The

*R*-values larger than 0.92 show a good correlation between the estimated attenuation coefficients with

^{2}*p*< 0.02 and the TiO

_{2}concentrations for all focus positions.

#### 3.2 Experimental results

As mentioned in the previous section, the experimental data of the multi-layer phantom were acquired with an experimental OCT systems using two different Rayleigh lengths. Figure 5 shows two typical examples of recorded B-scans obtained by the systems with the focus set to 100 μm inside the sample. The result of fitting the model to the averaged OCT signals (over 200 central A-lines) for the two systems and a series of focus positions, both above and inside the sample, are shown in Fig. 6 and Fig. 7. As can be seen visually, an acceptable fit to the measurements was obtained when the focus location is within 0.25 mm of the phantom’s surface, in the set of measurement with small Rayleigh lengths and for all positions of focus for measurements with large Rayleigh lengths.

The estimated parameters for each averaged A-line are shown in Fig. 8 for both sets of measurement. We expect to have similar attenuation coefficient values for each layer of the sample irrespective of the focus position. However, it can be seen in Fig. 8 (top row) that the estimated attenuation coefficients vary significantly for different focus positions. The estimated focus positions are shown in Fig. 8 (middle row). For small Rayleigh length, there is a good correlation between the estimated and expected focus position. However, the focus position could not be estimated reliably in the data set obtained with the large Rayleigh length.

In addition, we expected to observe a linear relation between the estimated attenuation coefficients and the particle concentration of the respective layer. Figure 9 (top) shows this linear relation to the estimated attenuation coefficients in the averaged A-lines acquired with the focus set to 0.6 mm inside the sample for both set of measurements. Figure 9 (middle) shows *R ^{2}*-values indicating the goodness of the linear fit for all the averaged A-lines and Fig. 9 (bottom) shows the corresponding p-values. The

*R*-values larger than 0.94 show a good correlation between the estimated attenuation coefficients with

^{2}*p*< 0.02 and the TiO

_{2}concentrations for all focus positions, except for focus positions higher than 0.3 mm above the sample, in the measurements obtained with the small Rayleigh length.

In the next step, we investigated if using multiple averaged A-lines with different focus locations improves the estimation of the attenuation coefficients. For this, we considered the averaged A-lines recorded with a different number of focus positions, i.e. (2, 4, 8 and 16). The attenuation coefficients of the identical layers among the averaged A-lines were considered to be common in the estimation process. The selected focus positions at different depths and their estimated values are shown in Fig. 10 and Fig. 11 for the measurement obtained using the OCT system with the small or the large Rayleigh length, respectively. The imposed (expected) and estimated focus position correlate well in measurements obtained using the OCT system with the small Rayleigh length. However, as was expected, an inaccurate estimation of z_{0} was obtained in measurements obtained using the OCT system with the large Rayleigh length. The estimated attenuation coefficients for both set of measurements (with small and large Rayleigh lengths) are shown in Tables 1 and 2. In addition, the *R*** ^{2}**-values of the linear fits to the estimated attenuation coefficients as a function of the particle concentration are also shown in these tables. For the measurements obtained with the small Rayleigh length, the

*R*

^{2}-values of the fits were higher than 0.99 (

*p*= 0.003) using 8 B-scans, which is slightly better than the results obtained using a single B-scan. For the measurements obtained with large Rayleigh length, the

*R*

**-values of the fits is the highest 0.99 (p-value < 0.006) when using 8 B-scans.**

^{2}## 4. Discussion

In this research, we introduced a method to estimate the attenuation coefficients of a multi-layer sample. A single-scattering model of the OCT signal was assumed while accounting for the system’s roll-off, noise and focused beam shape. To reduce the complexity of the proposed model, we assumed ${p_{NA}}(z ),$ refractive index and *C* parameter to be constant for all depths. The model parameters of the focused axial PSF and the attenuation coefficients of each layer were simultaneously estimated from experimental OCT data.

The numerical study predicted the optimal performance, and hence the inherent limitations, of our model for two experimental systems with different Rayleigh lengths. The simulation results indicate that the proposed method can estimate the selected attenuation coefficients with an acceptable precision and accuracy (CoV $< \; $10%) from the OCT signal obtained from a system with small Rayleigh length when focussing inside the sample. For the system with large Rayleigh length, while the precision of the estimation method is acceptable (CoV $< \; $10%), the bias is large due to erroneous estimation of the focus depth and the effect thereof on the estimated attenuation coefficients.

It has been shown that the average thickness of the retina and choroid of healthy eyes is about 250 µm [24] and 270 µm [25], respectively. To the best of our knowledge, in conventional ophthalmic OCT systems, the Rayleigh length is larger than 250 µm to have chorioretinal structures within the depth of focus. For large Rayleigh lengths, the intensity and consecutively the attenuation coefficient values are less affected by the shape of the beam and therefore it is less critical for accurate estimation of the attenuation coefficient. In high resolution imaging, which in ophthalmology can be done with adaptive optics, correcting for the effect of focus to estimate the attenuation coefficient is crucial. We observed a strong linear correlation (*R ^{2}* > 0.92, p-value < 0.02) between the estimated attenuation coefficients and the particle concentration of the respective layer for different Rayleigh lengths.

Experimental results obtained from a single B-scan of a multi-layer phantom composed of layers with different weight concentrations TiO_{2} in silicone, show an acceptable fit to the measurements when the focus location is within 0.3 mm of the phantom’s surface in the set of measurement with small Rayleigh lengths and for all positions of focus for measurements with large Rayleigh lengths. We observed a good correlation between the estimated and expected focus position for the measurements obtained with the small Rayleigh length where for larger Rayleigh length this correlation was not observed. This is mainly because the larger the Rayleigh length, less changes in the signal intensity is caused by the shape of beam and therefore it is more difficult to obtain the Rayleigh length from intensity data using the model in Eq. (2). However, with increasing Rayleigh lengths, the intensity and consecutively the attenuation coefficient values are less affected by the shape of the beam and therefore it is less essential to take into account this effect. In ophthalmic imaging the Rayleigh length is usually large and therefore the effect of focus location is reduced. However, in finer microscopic scales using adaptive optics, correcting for the effect of focus to estimate the attenuation coefficient is crucial.

Previously, we showed that in samples with different concentrations of TiO_{2} dispersed in silicone there is a linear relation between the TiO_{2} weight concentration and the estimated attenuation coefficients [17]. Using a single B-scan, we could observe that the estimated attenuation coefficients vary significantly for different focus positions while expecting to have similar attenuation coefficient values for each layer of the sample. Despite of this large variation, we observed a strong linear correlation (*R*^{2} > 0.94, p < 0.02) between the estimated attenuation coefficients and the particle concentration of the respective layer, in the measurements obtained with both Rayleigh lengths except for focus positions higher than 0.3 mm above the sample, in the measurements obtained with the small Rayleigh length. Using multiple B-scans with different focus locations for the measurements obtained with the small Rayleigh length system, the *R*^{2}-values of the fits were higher than 0.97 (p-value < 0.008) using 2, 3, 4 or 8 B-scans with different focus locations. The best result was obtained using the combination of 8 B-scans.

The clinical application of our proposed method for a multi-layer sample such as retinal tissue should be further investigated in future research. One of the limitations of the proposed method is that it requires the knowledge of the thickness of each layer per A-line and therefore accurate segmentation of the multi-layer sample such as retina [26,27] is a prerequisite. One of the limitations of our proposed method is that it requires accurate segmentation of the multi-layer sample. Other limitations of our technique, i.e. setting the location of focus inside the sample and the averaging rate should be considered while acquiring the OCT scans. In clinical practice, the operator of the OCT system aims to focus on the surface of the retina using a SLO camera, integrated into the OCT system. To increase the depth of focus in chorioretinal structures, the location of focus should be set to be inside the region of the retina. However due to the head movements and accommodation of the eye lens, the focus location can vary during acquisition. Further investigation is required to study the variation of focus location due to head movements and the eye lens accommodation. In addition, to study the effect of the model assumptions, further investigation in retinal tissue measurements is required.

## 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 Jelmer Weda, BSc, (Department of Physics and Astronomy, Vrije Universiteit Amsterdam, 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. **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]

**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. **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]

**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. Duc 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. **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]

**16. **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]

**17. **B. Ghafaryasl, K. A. Vermeer, J. Kalkman, T. Callewaert, J. F. de Boer, and L. J. van Vliet, “Analysis of attenuation coefficient estimation in Fourier-domain OCT of semi-definite media,” Biomed. Opt. Express **11**(11), 6093–6107 (2020). [CrossRef]

**18. **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]

**19. **V. M. Kodach, D. J. Faber, J. van Marle, T. G. van Leeuwen, and J. Kalkman, “Determination of the scattering anisotropy with optical coherence tomography,” Opt. Express **19**(7), 6131–6140 (2011). [CrossRef]

**20. **T. G. van Leeuwen, D. J. Faber, and M. C. Aalders, “Measurement of the axial point spread function in scattering media using single-mode fibre-based optical coherence tomography,” IEEE J. Sel. Top. Quant. Electr. **9**(2), 227–233 (2003). [CrossRef]

**21. **J. Kalkman, “Fourier-domain optical coherence tomography signal analysis and numerical modeling,” Int. J. Opt. **2017**, 1–16 (2017). [CrossRef]

**22. **B. Braaf, K. A. Vermeer, K. V. Vienola, and J. F. de Boer, “Angiography of the retina and the choroid with phase-resolved OCT using interval-optimized backstitched B-scans,” Opt. Express **20**(18), 20516–20534 (2012). [CrossRef]

**23. **E. W. Dijkstra, “A note on two problems in connexion with graphs,” Numer. Math. **1**(1), 269–271 (1959). [CrossRef]

**24. **B. Alamouti and J. Funk, “Retinal thickness decreases with age: an OCT study,” Br. J. Ophthalmol. **87**(7), 899–901 (2003). [CrossRef]

**25. **V. Manjunath, M. Taha, J. G. Fujimoto, and J. S. Duker, “Choroidal thickness in normal eyes measured using Cirrus HD optical coherence tomography,” Am. J. Ophthalmol. **150**(3), 325–329.e1 (2010). [CrossRef]

**26. **J. Novosel, G. Thepass, H. G. Lemij, J. F. de Boer, K. A. Vermeer, and L. J. van Vliet, “Loosely coupled level sets for simultaneous 3D retinal layer segmentation in optical coherence tomography,” Med. Image Anal. **26**(1), 146–158 (2015). [CrossRef]

**27. **J. Novosel, K. A. Vermeer, J. H. de Jong, Z. Wang, and L. J. van Vliet, “Joint segmentation of retinal layers and focal lesions in 3-D OCT data of topologically disrupted retinas,” IEEE Trans. Med. Imaging **36**(6), 1276–1286 (2017). [CrossRef]