## Abstract

The attenuation coefficient provides a quantitative parameter for tissue characterization and can be calculated from optical coherence tomography (OCT) data, but accurate determination requires compensation for the confocal function. We present extensive measurement series for extraction of the focal plane and the apparent Rayleigh length from the ratios of OCT images acquired with different focus depths and compare these results with two alternative approaches. By acquiring OCT images for a range of different focus depths the optimal focus plane difference is determined for intralipid and titanium oxide (TiO_{2}) phantoms with different scatterer concentrations, which allows for calculation of the attenuation coefficient corrected for the confocal function. The attenuation coefficient is determined for homogeneous intralipid and TiO_{2} samples over a wide range of concentrations. We further demonstrate very good reproducibility of the determined attenuation coefficient of layers with identical scatter concentrations in a multi-layered phantom. Finally, this method is applied to in vivo retinal data.

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

## 1. Introduction

Optical coherence tomography (OCT) is primarily used to image the morphology of tissue by visualizing different structures in arbitrary units to diagnose pathologies in medical specialties such as ophthalmology, cardiology, dermatology, and urology. However, the state of a tissue can also be described in absolute units by its inherent optical properties, such as the attenuation coefficient, providing a quantitative tissue parameter. Vermeer et al. [1] introduced and experimentally demonstrated a method that allows estimation of the attenuation coefficient with a pixel level resolution from OCT depth profiles. If the focal plane is located close to or within the sample the confocal function for the OCT measurement geometry [2] has to be taken into account to precisely determine the attenuation coefficient. Smith et al. [3] suggested a method for automated, depth-resolved estimation of the attenuation coefficient with the focal plane located within the sample; however, this method requires a priori knowledge of the confocal function. Determining the confocal function constitutes estimating the position of focus and the apparent Rayleigh length.

The confocal function also depends on the sample that is examined. In many applications, especially ophthalmology, the confocal function parameters, i.e., focus location and apparent Rayleigh length vary. It would be beneficial to extract these parameters from data acquired with the matching setup and sample or subject. Ghafaryasl et al. [4] have suggested a method to estimate the confocal function and the attenuation coefficient directly from a single scan. However, their method is currently limited to homogeneous samples. Dwork et al. [5,6] have demonstrated that the ratio of two B-scans can be used to estimate the focus location and the apparent Rayleigh length and that this can be applied to inhomogeneous samples. Stefan et al. [7] extended this technique by accounting for the curvature and tilt of the focal plane.

In this paper, we analyze the applicability and limitations of extracting the focus depth and/or the apparent Rayleigh length from OCT images using three different methods (Section 3.7 focus series, Section 3.8 ratio fit, and Section 3.10 direct fit). We find that the most robust and versatile method to extract the confocal function parameters in the sample is to use the ratio fit of two B-scans acquired with different focal depths, which is based on the work of Dwork et al. [5,6]. We investigate this method by imaging homogeneous intralipid samples and titanium oxide (TiO_{2}) phantoms of different scatterer concentrations with an OCT system through a series of different focal depths from the sample surface. The data series are analyzed for the optimal focal depth offset for confocal function determination. The obtained focus depths and Rayleigh lengths are used to correct for the confocal function and subsequently calculate the attenuation coefficients. For comparison, we also test a technique where we directly fit the confocal function and attenuation coefficient in a single A-scan. Finally, the ratio fit method is demonstrated for a layered phantom and in vivo retinal data by determining the curved focal planes and calculating attenuation coefficient maps.

## 2. Theory

In this section, we summarize the main aspects of the Gaussian beam model needed to understand how raw OCT data can be corrected for the confocal function. Single backscattered light detected within an OCT A-scan can be described by [1]

where*z*is the depth within the A-scan,

*r*(

*z*) is the sensitivity roll-off,

*h*(

*z*) is the confocal function,

*α*is the backscattering fraction, and

*µ*(

_{t}*z*) is the attenuation coefficient. The attenuation coefficient

*µ*is the sum of the scattering coefficient

_{t}*µ*and the absorption coefficient

_{s}*µ*,

_{a}*µ*=

_{t}*µ*+

_{s}*µ*. For a Gaussian beam, the confocal function

_{a}*h*(

*z*) can be described by [2] with the focus position

*z*and the apparent Rayleigh length

_{f}*z*For a scattering sample with refractive index

_{R}.*n*the apparent Rayleigh length can be calculated as [2] with the central wavelength $\lambda $ and the $1/{e^2}$ radius of the lateral intensity distribution of the focused Gaussian beam ${\omega _0}$ which can be estimated as with the focal length

*f*and the beam diameter at the pupil

*d*.

To determine *z _{f}* and

*z*one can exploit the common characteristics of two scans of the same sample location [5]. If two A-scans of the same structure are acquired with different focus depth settings, all the terms in Eq. (1) are the same, except for

_{R}*h*(

*z*). Therefore, the ratio of the intensity of the numerator A-scan

*I*

_{1}(

*z*) and the denominator A-scan

*I*

_{2}(

*z*) can be described by the ratio of the two confocal functions as

*z*=

_{f}*z*is the axial shift between the two focus depths, and

_{f2}- z_{f1}*z*and

_{f1}*z*are the focus positions of A-scan number one and two, respectively. Equation (5) shows that taking the ratio of two A-scans acquired with a focus offset allows for the determination of

_{f2}*z*and

_{f}*z*, thereby estimating the confocal function parameters [6].

_{R}## 3. Methods

#### 3.1 OCT system

OCT measurements were performed with a high-resolution spectral domain OCT device based on the SPECTRALIS platform (Heidelberg Engineering GmbH, Heidelberg, Germany) with a central wavelength of 853 nm and a spectrometer bandwidth of 137 nm. Using the active tracking capability of the device many B-scans of the same location can be averaged together. The focus of the device can be changed by moving the objective relative to the rest of the device using a focusing knob. The focus setting in diopters (1/m) is displayed in the graphical user interface of the acquisition software and saved with each data set.

#### 3.2 Data acquisition and B-scan processing

For each sample, a focus series of B-Scans was acquired at various depths of the focal plane referenced to the sample surface. All the data processing was performed on the raw spectral data stored by the device using custom-written software in Python 3.7 (Python Software Foundation) and MATLAB 2020b (The MathWorks, Inc.).

First, to reduce fixed pattern noise [8], the raw spectra were divided by the mean spectrum of all A-scans and centered around zero by subtracting one. Second, the nonlinear relation between wavenumber k and wavelength λ was compensated for by cubic interpolation to linearize the signal in wavenumber [9]. Third, to suppress sidelobes caused by the Fourier transformation of a finite interval, the spectrum was windowed using the familiar Hamming window (a_{0} = 0.53836 and a_{1} = 0.46164). Fourth, numerical dispersion compensation was performed [10,11]. Fifth, the Fourier transformation of the signal was squared. Sixth, the depth dependent background was subtracted, which was measured by acquiring a B-scan without a sample. Seventh, to account for sensitivity changes with depth, the obtained B-scans were corrected for roll-off which was measured by moving a mirror through the sampling depth of the OCT scan [12]. Finally, to reduce speckle noise and average over imaging artifacts, 120 B-scans were averaged (768 A-scans/B-scan) to produce one B-scan in which each A-scan was averaged 120 times. These A-scans were used for further processing.

#### 3.3 Model eye

For all phantom measurements, a model eye was used featuring an 18 mm focal length planoconvex lens (01LPX015, Melles Griot, anti-reflection coating on the convex side) and a water filled chamber that can hold different types of samples.

#### 3.4 Intralipid samples

To assess the effect of different scattering strengths, a serial dilution ranging from 6.22% volume to 0.062% volume intralipid (Table 1) was produced by diluting an intralipid 20% emulsion (I141-100ML, Sigma-Aldrich, Saint Louis, MO, USA) with distilled water. Cuvettes (Macro cell 100-5-20, Hellma GmbH & Co. KG, Müllheim, Germany) were filled with the dilutions and imaged inside the model eye with a scan angle of 9 degrees. The nominal scattering coefficients for the OCT wavelength of the different dilutions were calculated by linear interpolation of the scattering coefficient for intralipid 20% described by Michels et al. [13]. The refractive index of the samples was assumed to be equal to the refractive index of water and the absorption coefficient of the intralipid samples was assumed to be negligible (*µ _{s}* ≫

*µ*and

_{a,}*µ*=

_{t}*µ*) [13].

_{s}#### 3.5 Titanium oxide phantoms

Silicone-based phantoms were fabricated similar to the method outlined by de Bruin et al. [14], where the scattering coefficient is controlled by varying the % weight of TiO_{2} in the final elastomer. A TiO_{2} (634662, Sigma Aldrich, Saint Louis, MO, USA) and PDMS (481955, Sigma Aldrich) stock dispersion was thoroughly incorporated in an elastomer base (SYLGARD 184 Silicone Elastomer, Dow, Inc., Midland, MI, USA) then cast into the final phantom geometry and cured. Seven bulk phantoms and one layered phantom were fabricated with TiO_{2}% weight concentrations ranging from 0.011 to 0.730% weight, outlined in Table 2. The layered phantom (Fig. 1) was created by casting three thin films of the substrate with different scatterer concentrations on top of each other and afterwards stacking five of these layers. This procedure was chosen to ensure that films with the same nominal scattering coefficient have the exact same physical properties. In our experience, even multiple pours of the same formulation can be unreliable in producing the same physical properties. Although the samples are flat, the apparent curvature of the sample, that can be observed in Fig. 1, is a typical OCT imaging artifact originating from the optical path length differences of the beam depending on the scan angle. The curvature depends on the position of the scan pupil relative to the focusing lens and on the lens itself. All phantoms were imaged inside the model eye with a scan angle of 9 degrees. The refractive index of the phantoms was assumed to be 1.42 [14], and the absorption coefficient was neglected, *µ _{t}* =

*µ*, as no absorbers were added.

_{s}#### 3.6 In vivo data

To demonstrate the applicability of our analysis to in vivo data, scans of the retinas of two subjects were obtained at the locations shown in Fig. 12(C) and (F) using the tracking capability of the device to ensure acquisition at the same location for all focus settings. To avoid accommodation during the measurement an external fixation target was used. No drugs were applied to dilate the pupil or paralyze accommodation. This paper represents a feasibility study with the focus on the evaluation and correction algorithms and not a systematic analysis of clinical data. A written declaration of informed consent was obtained from both subjects.

#### 3.7 Expected focus depth and Rayleigh length estimation using the focus series

To estimate the theoretical expected focus depth and the actual apparent Rayleigh length, versus the theoretical apparent Rayleigh length calculated using Eq. (3), a focus series fit was performed. To suppress high frequency noise, a two-dimensional Gaussian filter with a standard deviation of 3 pixels was applied to the B-scans and to further reduce noise a mean A-scan was calculated for each B-scan by averaging the 61 central A-scans (location shown in Fig. 2(A)). Next, the intensity just below the surface of the sample was determined for every focus position (Fig. 2(B), red dashed line). This location was chosen to avoid the surface reflection and minimize multiple scattering contributions expected to increase with depth. Each determined intensity value corresponds to a focus setting in diopters (1/m) which was saved by the system together with the B-scan data. The focal length *f* can be calculated from the diopter value *D* in 1/m and the focal length in air of the lens used for imaging *f _{lens}* in meters according to

*n*. Because the medium above the sample has a different refractive index than the sample itself, the theoretical focus distance from the sample surface

*Δz*was determined to fit for the Rayleigh length. Therefore, the focus setting at which the sample surface would be perfectly in focus

_{fsample}*D*was estimated by fitting

_{sample}*A*([(

*D*-

*D*)/

_{sample}*R*]

^{2}+1)

^{−1}to the intensity values, with

*D*,

_{sample}*R*and

*A*as fitting parameters. Using Eq. (6)

*f*was calculated from

_{sample}*D*. The theoretical focus distance from the sample surface was calculated as

_{sample}*Δz*, where

_{fsample}= f - f_{sample}*f and f*were calculated with the respective refractive index depending on if the focus was above the sample (

_{sample}*D*>

*D*) or within the sample (

_{sample}*D*<

*D*). To estimate the apparent Rayleigh length the Gaussian beam confocal function (Eq. (2)) multiplied by a constant

_{sample}*A*was fitted to the intensity values over

*Δz*as described by Faber et al. [15] with

_{fsample}*A*,

*z*and

_{f}*z*as fitting parameters (Fig. 2(C)).

_{R}#### 3.8 Determination of confocal function using the ratio fit

To extract the confocal function parameters, Eq. (5) describing the intensity ratio of two A-scans (numerator *I*_{1}(*z*) and denominator *I*_{2}(*z*)) acquired with two different focus depths was used in the decibel form as

*z*and

_{f1}*z*. This expression was fitted to the ratio of two A-scans acquired with a different focus setting using a Levenberg-Marquardt nonlinear least squares algorithm (MATLAB function fitnlm, example shown in Fig. 3). The focus offset

_{R}*Δz*is determined from the respective refraction settings of the OCT system. The signal to noise ratio (SNR) is low above the sample surface and deep within the sample. To limit the fitting range to a region with high SNR the fit was weighed with the sum of the intensity of the two A-scans on a linear scale and limited to the part below the sample surface. Any negative or zero ratio values were omitted. The ratio fit was performed for all possible combinations of different focus settings resulting in matrices as depicted in Fig. 4 for intralipid sample I3 with a nominal scattering coefficient of 2 mm

_{f}^{-1}(see Supplement 1 for the other intralipid concentrations). The mean focus depth and Rayleigh length for the respective numerator A-scan (i.e., for each horizontal row of Fig. 4) was determined by calculating the mean of the fit results each weighted with the standard error of the fit. Using the resulting mean focus depth and Rayleigh length the confocal function was estimated for each focus depth.

#### 3.9 Calculation of attenuation coefficient

To correct for the axial point spread function (PSF), the A-scans were divided by the estimated confocal function described by the mean focus depth and Rayleigh length obtained for each respective focus setting. Subsequently, the depth resolved attenuation coefficient profile *µ _{t}*(

*i*) was calculated using the method as described by Girard et al [16] and Vermeer et al. [1] as

*I*(

*i*) is the intensity at the

*i*’th pixel along the A-scan and Δ the pixel size. We refer to this method as the depth resolved method. It is based on two main assumptions: (1) all the light is attenuated within the range of the A-scan; and (2) the backscattered light detected by the OCT system is a fixed fraction of the attenuation coefficient, i.e., the backscattering fraction is constant throughout the sample depth.

For the homogeneous phantoms, a reference value for the attenuation coefficient was estimated using the commonly applied slope fitting method [15,17] where a single scattering model in the form *I*(*z*) ∝ exp(-2*µz)* is fitted to the A-scans.

#### 3.10 Direct fit for confocal function and attenuation coefficient

Additionally, for the homogeneous phantoms, the method introduced by Ghafaryasl et al. [4] was tested. It has the advantage that focal position and Rayleigh length do not have to be determined beforehand and only a single scan is needed, removing the necessity to acquire two scans from the exact same position. The single scattering model

*z*, apparent Rayleigh length

_{f}*z*, attenuation coefficient

_{R}*µ*and scaling factor

_{t}*C*as fitting parameters using MATLAB’s fitnlm function. However, due to its single scattering coefficient this method can only be applied to homogeneous samples.

## 4. Results

Sections 4.1 to 4.3 describe the results obtained from the homogeneous samples, which were analyzed by calculating a mean A-scan from the 61 central A-scans of an averaged B-scan. This simplified our analysis because within these A-scans we assumed a uniform focus depth. The results for the layered phantom (Section 4.4) and the in vivo data (Section 4.5) were obtained by analyzing the whole averaged B-scans which were filtered as described in Section 4.4.

#### 4.1 Ratio fit data processing

As an example, Fig. 4 illustrates the ratio fit result matrices for the intralipid sample I3 (nominal *µ _{s}* of 2 mm

^{-1}), showing the results for all possible combinations of different focus depths fitted in Eq. (7) (color-coded for the standard errors of each fit). The fit matrices for the other intralipid concentrations are provided in the Supplement 1.

Ratio combinations with both focus depths (numerator and denominator) above the sample show a very large standard error and mostly fitting the model (Eq. (7)) fails which is reflected in fitted apparent Rayleigh lengths of zero. The lowest fit errors are obtained for ratios where at least one focus location is located well inside the sample with a focus difference between the two tested A-scans of more than 100 µm. For focus locations very deep in the sample an increasing fitting error can be observed.

Calculating the error-weighted mean fit result for each horizontal line of the ratio fit result matrices (Fig. 4) results in the values given in Fig. 5. The mean fitted focus depths agree well with the expected focus depths obtained from the focus settings of the device, and the fitted mean Rayleigh lengths are located around the expected Rayleigh length calculated based on the beam parameters (Eq. (3)).

Using the confocal function parameters recorded in Fig. 5, the respective mean A-scans shown in Fig. 6(A) were corrected for the confocal function. The uncorrected mean A-scans (Fig. 6(A)) acquired with different focus depths clearly show the effect of the confocal function, which is eliminated in the corrected mean A-scans (Fig. 6(B)).

#### 4.2 Fitted Rayleigh length

We determined the Rayleigh length using three different methods (as described in Section 3.7 focus series fit, Section 3.8 ratio fit, and Section 3.10 direct fit) for the homogeneous intralipid and TiO_{2} phantoms with several different scatterer concentrations. For each method and scatterer concentration, the average Rayleigh length was determined (Fig. 7) by calculating the weighted mean of the Rayleigh lengths calculated for focus depths from 0 µm to 1200 µm, weighted with the corresponding standard error. For nominal scattering coefficients greater than 2 mm^{−1} the standard deviation of the direct fit exceeded 50 µm, which we considered as a failed fit and values are therefore not shown. Figure 7 shows that the fitted apparent Rayleigh length increases with increasing scatterer concentration for the ratio fit method. The same trend can be seen for the direct fit results. However, for the focus series fit the trend of increasing apparent Rayleigh length with increasing scatterer concentration is not as pronounced. These values were calculated from the intensities just below the surface minimizing the contribution of multiple scattering. For the intralipid samples (Fig. 7(A)) the data series are more consistent and show less fluctuation than for the TiO_{2} phantoms (Fig. 7(B)). Here, it is worth mentioning that the averaged OCT B-scans of the liquid intralipid samples (Fig. 2(A)) appear much more uniform and show less speckle noise compared to the B-scans of the solid TiO_{2} phantoms (Fig. 1) because Brownian motion averages out the speckle.

#### 4.3 Attenuation coefficient

To demonstrate the validity of the ratio fit method for determining the focus position and Rayleigh length in the sample of a particular mean A-scan, the attenuation coefficients calculated from the mean A-scans acquired with different focus settings are depicted in Fig. 8 (corrected for the fitted Rayleigh lengths and focus depths (Fig. 5)). Both, the depth resolved method and the slope fit method result in a calculated attenuation coefficient that is more or less constant as a function of the focus depth. If the confocal function is not corrected for, a reliable estimation of the attenuation coefficient is not possible which can also be seen in Fig. 8, where the values for the uncorrected mean A-scans vary considerably.

The weighted mean attenuation coefficients for focus locations from 0 µm to 1200 µm below the surface, weighted with the standard error of the corresponding fit, are shown in Fig. 9 for intralipid (A) and TiO_{2} (B) samples with different scatterer concentrations. The direct fit results were obtained directly from uncorrected mean A-scans, whereas slope and depth resolved values were calculated from the confocal function corrected mean A-scans using the parameters from the ratio fit. For scattering coefficients lower than 6 mm^{-1} the calculated attenuation coefficients agree well with the theory independent of the algorithm used. For strong scattering, the attenuation coefficients are consistently underestimated by both the depth resolved and the slope method (the direct fit failed in this range).

#### 4.4 Layered phantom

In general, the focal plane cannot always be assumed to be flat and is rather curved and/or tilted within the sample [7]. For the layered TiO_{2} phantom, the focal planes of the different focus settings were determined by applying the ratio fit method to every A-scan extracting any curvature or tilt of the focal plane. Prior to the confocal function determination, a two-dimensional Gaussian filter with a standard deviation of 5 pixels and a moving mean filter in the x direction with a window size of 11 pixels was applied to suppress high frequency noise. The parameters for the two filters were chosen as a compromise between noise reduction and conserving spatial resolution. For the ratio fit we were interested in the confocal function which varies slowly across the B-scan, therefore preserving hard edges was not a priority and some blurring was accepted. To reduce noise in the raw two-dimensional map of the extracted focus plane and Rayleigh lengths, we fitted a second order polynomial to the extracted raw *z _{f}* and

*z*values and used the smoothed values for the further analysis. The resulting smooth curved focal planes can be seen in Fig. 10(A). The resulting confocal functions were used to correct the unfiltered B-scans, from which attenuation coefficient maps were calculated using the depth resolved method. The attenuation coefficient map obtained for the B-scan with the focal plane position indicated by the white arrow in Fig. 10(A) is shown in Fig. 10(B). Layers with the same scatterer concentration exhibit similar absolute attenuation coefficient values. The mean attenuation coefficients were calculated for each layer marked in Fig. 10(C) and every focal plane location. The resulting values are given in Fig. 10(D) demonstrating the reproducibility of the attenuation coefficient calculation and the effectiveness of the confocal function correction. Independent of the depth within the phantom, layers produced from the same pour (A1 = A2, B1 = B2, C1 = C2) exhibit similar attenuation coefficients. The obtained

_{R}*µ*are in the same range as the nominal

_{t}*µ*(Table 2) indicated with the dashed lines but have an offset.

_{s}#### 4.5 In vivo data

The in vivo OCT images from two normal human retinas were calculated by averaging 120 B-scans for every focus setting using the registration information from the SPECTRALIS system. The data was processed in the same way as the layered TiO_{2} phantom images, extracting the curved focal planes and the Rayleigh length, correcting for the confocal function, and calculating attenuation coefficient maps. Figure 11 shows an example of the resulting attenuation coefficient maps for subject 1 and 2 for a focus position in the lateral center of the scans 200 µm and 147 µm below the inner limiting membrane (ILM), respectively. The different retinal layers visually show comparable attenuation coefficients for the two subjects, however, the map of subject 1 shows lower values in the center of the B-scan, where the retinal nerve fiber layer (RNFL) is thinnest. Note, that the scan for subject 1 was acquired closer to the fovea than for subject 2 (see Fig. 12(C) and (F)). A-scans with apparent blood flow present (Fig. 12(A) and (D), marked in red) were excluded from further calculations because the strong absorption of blood clearly violates the assumption of a constant backscattering fraction for the depth resolved method.

The mean attenuation coefficients calculated for the different layers marked in Fig. 12(A) and (D) for different focus locations are shown in Fig. 12(B) and (E) for subject 1 and 2, respectively. The segmentation was computed using open-source software [18,19]. Table 3 lists the mean and standard deviation of the attenuation coefficients acquired with central focus depths of 0 µm to 500 µm below the ILM. The *µ _{t}* for subject 1 are mostly slightly lower and show a higher standard deviation than for subject 2.

## 5. Discussion

We investigated a technique to extract the Rayleigh length and focus depth from a set of OCT B-scans acquired with different focus settings. This technique is based on the method introduced by Dwork et al. [5,6] and extended by Stefan et al. [7]. We applied this ratio fit method to a range of different foci and showed experimentally that the ratio fit works best in a focus depth range from about 0 µm to 500 µm and a focus offset of more than 100 µm by analyzing scattering phantoms of different concentrations. The ratio fit tends to fail if both B-scans used for the ratio calculation are acquired with a focus location above the sample. The obtained focus depth and Rayleigh length were used to correct for the confocal function and subsequently calculate the attenuation coefficients, which agreed well with the nominal scattering coefficients expected from linear extrapolation. The observed deviation of the calculated attenuation coefficients from the nominal scattering coefficients greater than 6 mm^{−1} can be explained with the increasing presence of multiple and dependent scattering [20–22]. As a reference we also used a method where we directly fitted the confocal function and attenuation coefficient in a single A-scan [4], this has the advantage that no change in focus is needed during the measurement. However, the applicability is limited to homogeneous samples with a scattering coefficient lower than 4 mm^{−1}.

We observed a concentration dependent increase of the apparent Rayleigh length (Fig. 7) when using the ratio fit method. This method is based on the single scattering model and might be increasingly inaccurate with a higher frequency of multiple scattering events. This is expected to happen with both an increased scatterer concentration and at larger imaging depths. The experimental determination of the Rayleigh length using the focus series fit also showed a trend of longer Rayleigh lengths with increasing scatterer concentration, however, this trend is less pronounced than with the ratio fit method. In the focus series fit, multiple scattering effects should be less pronounced because only one location at the top of the sample was analyzed.

The applicability of the ratio fit method was also demonstrated for a layered phantom, where the two-dimensional focal planes could consistently be determined. The reproducibility of the determined attenuation coefficients was excellent, as demonstrated by the equal attenuation coefficients found for the same scattering layers at different depths from the sample surface.

Finally, the ratio fit method was used to extract the confocal function from in vivo retinal scans and calculate attenuation coefficient maps. For the two subjects analyzed we observed a small mismatch between the calculated attenuation coefficients of the retinal layers, however, this might be explained by the different locations of the B-scans relative to the fovea and natural intra- and inter-subject variability. The results on just two subjects seem to indicate that the attenuation coefficient of the RNFL becomes smaller closer to the fovea (Fig. 11(A) and Table 3). Additionally, the less consistent focus dependent attenuation coefficients for subject 1 (Fig. 12(B)) might be explained by accidental accommodation during the measurement making the focus reading of the device inaccurate and possibly also leading to averaging B-scans with different focus locations.

Reliably estimating the attenuation coefficient would add an additional parameter to quantify pathological changes in tissue and help to diagnose those changes earlier. In the field of ophthalmology, it has been reported that the attenuation coefficient might help to distinguish glaucomatous eyes from healthy eyes [23,24] and that diagnostic information could be derived for diabetes [25] and pituitary adenoma [26]. Additionally, extracting the confocal function parameters and subsequently correcting for the axial PSF could be of great value for other quantitative OCT analysis.

For the theory of the ratio fit and the depth resolved attenuation coefficient calculation to apply (as needed for the calculations presented here), we are limited by prior conditions that need to be fulfilled. Explicitly, Eq. (1-6) describing the intensity distribution and the confocal function are only valid if Gaussian optics can be applied, i.e., low numerical aperture and undistorted probe beam are given [15]. Both Eq. (5) (ratio of two A-scans) and Eq. (8) (depth resolved attenuation coefficient profile) are for single scattering only and become more inaccurate with increasing multiple scattering. Furthermore, the depth resolved method can only be applied reliably if the detected light is a fixed fraction of the attenuated light, and if the signal is completely attenuated within the depth of the scan [1]. Some of the mentioned requirements are only partly met for retinal imaging and further refinement is needed.

## Funding

Heidelberg Engineering GmbH; Health~Holland, Topsector Life Sciences & Health; Holland High Tech, Topsector High Tech Systems and Materials.

## Acknowledgments

We acknowledge funding by Heidelberg Engineering GmbH. The collaboration project is co-funded by the PPP Allowance made available by Health∼Holland, Topsector Life Sciences & Health, and by Holland High Tech, Topsector High Tech Systems and Materials, to stimulate public-private partnerships.

## Disclosures

JK: Heidelberg Engineering GmbH (F,E), JFdB: Heidelberg Engineering GmbH (F), JF: Heidelberg Engineering GmbH (E).

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## Supplemental document

See Supplement 1 for supporting content.

## References

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

**2. **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. Sel. Top. Quantum Electron. **9**(2), 227–233 (2003). [CrossRef]

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

**4. **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-infinite media,” Biomed. Opt. Express **11**(11), 6093–6107 (2020). [CrossRef]

**5. **N. Dwork, G. T. Smith, J. M. Pauly, and A. K. Ellerbee Bowden, “Automated estimation of OCT confocal function parameters from two B-scans,” in Conference on Lasers and Electro-Optics, OSA Technical Digest (online) (Optical Society of America, 2016), AW1O.4.

**6. **N. Dwork, G. T. Smith, T. Leng, J. M. Pauly, and A. K. Bowden, “Automatically determining the confocal parameters from OCT B-scans for quantification of the attenuation coefficients,” IEEE Trans. Med. Imaging **38**(1), 261–268 (2019). [CrossRef]

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

**8. **N. Nassif, B. Cense, B. H. Park, S. H. Yun, T. C. Chen, B. E. Bouma, G. J. Tearney, and J. F. de Boer, “In vivo human retinal imaging by ultrahigh-speed spectral domain optical coherence tomography,” Opt. Lett. **29**(5), 480–482 (2004). [CrossRef]

**9. **M. Wojtkowski, R. Leitgeb, A. Kowalczyk, T. Bajraszewski, and A. F. Fercher, “In vivo human retinal imaging by Fourier domain optical coherence tomography,” J. Biomed. Opt. **7**(3), 457–463 (2002). [CrossRef]

**10. **J. F. de Boer, C. E. Saxer, and J. S. Nelson, “Stable carrier generation and phase-resolved digital data processing in optical coherence tomography,” Appl. Opt. **40**(31), 5787–5790 (2001). [CrossRef]

**11. **D. Hillmann, T. Bonin, C. Luhrs, G. Franke, M. Hagen-Eggert, P. Koch, and G. Huttmann, “Common approach for compensation of axial motion artifacts in swept-source OCT and dispersion in Fourier-domain OCT,” Opt. Express **20**(6), 6761–6776 (2012). [CrossRef]

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

**13. **R. Michels, F. Foschum, and A. Kienle, “Optical properties of fat emulsions,” Opt. Express **16**(8), 5907 (2008). [CrossRef]

**14. **D. M. de Bruin, R. H. Bremmer, V. M. Kodach, R. de Kinkelder, J. van Marle, T. G. van Leeuwen, and D. J. Faber, “Optical phantoms of varying geometry based on thin building blocks with controlled optical properties,” J. Biomed. Opt. **15**(2), 025001 (2010). [CrossRef]

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

**16. **M. J. Girard, N. G. Strouthidis, C. R. Ethier, and J. M. Mari, “Shadow removal and contrast enhancement in optical coherence tomography images of the human optic nerve head,” Invest. Ophthalmol. Vis. Sci. **52**(10), 7738–7748 (2011). [CrossRef]

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

**18. **S. J. Chiu, X. T. Li, P. Nicholas, C. A. Toth, J. A. Izatt, and S. Farsiu, “Automatic segmentation of seven retinal layers in SDOCT images congruent with expert manual segmentation,” Opt. Express **18**(18), 19413–19428 (2010). [CrossRef]

**19. **P.-y. Teng, “Caserel-an open source software for computer-aided segmentation of retinal layers in optical coherence tomography images,” Zenodo. 10.5281/zenodo.17893 (2013).

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

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

**22. **M. Almasian, N. Bosschaart, T. G. van Leeuwen, and D. J. Faber, “Validation of quantitative attenuation and backscattering coefficient measurements by optical coherence tomography in the concentration-dependent and multiple scattering regime,” J. Biomed. Opt. **20**(12), 121314 (2015). [CrossRef]

**23. **J. van der Schoot, K. A. Vermeer, J. F. de Boer, and H. G. Lemij, “The effect of glaucoma on the optical attenuation coefficient of the retinal nerve fiber layer in spectral domain optical coherence tomography images,” Invest. Ophthalmol. Vis. Sci. **53**(4), 2424–2430 (2012). [CrossRef]

**24. **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. Vis. Sci. **53**(10), 6102–6108 (2012). [CrossRef]

**25. **D. C. DeBuc, W. Gao, E. Tátrai, L. Laurik, B. Varga, V. Ölvedy, W. E. Smiddy, R. Tchitnga, A. Somogyib, and G. M. Somfái, “Extracting diagnostic information from optical coherence tomography images of diabetic retinal tissues using depth-dependent attenuation rate and fractal analysis,” in * Biomedical Optics and 3-D Imaging*, OSA Technical Digest (Optical Society of America, 2012), BTu3A.74.

**26. **M. Sun, Z. Zhang, C. Ma, S. Chen, and X. Chen, “Quantitative analysis of retinal layers on three-dimensional spectral-domain optical coherence tomography for pituitary adenoma,” PLoS One **12**, e0179532 (2017).