Bayesian analysis of depth resolved OCT attenuation coefficients

Optical coherence tomography (OCT) is an optical technique which allows for volumetric visualization of the internal structures of translucent materials. Additional information can be gained by measuring the rate of signal attenuation in depth. Techniques have been developed to estimate the rate of attenuation on a voxel by voxel basis. This depth resolved attenuation analysis gives insight into tissue structure and organization in a spatially resolved way. However, the presence of speckle in the OCT measurement causes the attenuation coefficient image to contain unrealistic fluctuations and makes the reliability of these images at the voxel level poor. While the distribution of speckle in OCT images has appeared in literature, the resulting voxelwise corruption of the attenuation analysis has not. In this work, the estimated depth resolved attenuation coefficient from OCT data with speckle is shown to be approximately exponentially distributed. After this, a prior distribution for the depth resolved attenuation coefficient is derived for a simple system using statistical mechanics. Finally, given a set of depth resolved estimates which were made from OCT data in the presence of speckle, a posterior probability distribution for the true voxelwise attenuation coefficient is derived and a Bayesian voxelwise estimator for the coefficient is given. These results are demonstrated in simulation and validated experimentally.

Optical coherence tomography (OCT) is an optical technique which allows for volumetric visualization of the internal structures of translucent materials. Additional information can be gained by measuring the rate of signal attenuation in depth. Techniques have been developed to estimate the rate of attenuation on a voxel by voxel basis. This depth resolved attenuation analysis gives insight into tissue structure and organization in a spatially resolved way. However, the presence of speckle in the OCT measurement causes the attenuation coefficient image to contain unrealistic fluctuations and makes the reliability of these images at the voxel level poor. While the distribution of speckle in OCT images has appeared in literature, the resulting voxelwise corruption of the attenuation analysis has not. In this work, the estimated depth resolved attenuation coefficient from OCT data with speckle is shown to be approximately exponentially distributed. After this, a prior distribution for the depth resolved attenuation coefficient is derived for a simple system using statistical mechanics. Finally, given a set of depth resolved estimates which were made from OCT data in the presence of speckle, a posterior probability distribution for the true voxelwise attenuation coefficient is derived and a Bayesian voxelwise estimator for the coefficient is given. These results are demonstrated in simulation and validated experimentally.
Optical coherence tomography (OCT) is an imaging modality which allows for the visualization of internal structures of tissues and other translucent materials volumetrically. OCT images give a practitioner insight into qualitative structural information such as layer structure and morphology. However, the extraction of reliable quantitative information from these tissue volumes is an area of current research. One quantitative measure of interest is the rate of signal decay in depth known as the attenuation coefficient 1 . The attenuation coefficient compounds effects of absorption and scattering losses in depth which can be related to physiological properties such as blood content and tissue organization [1][2][3] . Currently, methods to extract the attenuation coefficient fall into one of two categories: layerwise extraction through curve fitting 4 and depth resolved or voxelwise extraction 5 .
In the layerwise approach, the layers of media are segmented, and then an exponentially decaying model is fit to each A-scan of the OCT signal in the least squares sense [3][4][5] . From this perspective, the attenuation coefficient is a bulk measure which assigns a single, deterministic number to each segment of an A-scan. However, a measured A-scan will contain fluctuations due to speckle 6,7 . OCT speckle is the voxel-to-voxel variation of OCT amplitude, due to random variations in the spatial position of scattering particles within the imaging voxel. Randomly placed scatterers within the voxels will thus return scattered fields with random amplitude and phase-leading to intensity fluctuations at the detector. While the origin of speckle is deterministic at the microscopic level, in practice the measured signal is well modeled as a realization of a random process equivalent to randomly varying the exact microscopic position of the scattering particles in the bulk of the media 8 . One common technique to overcome the speckle variations is lateral averaging 2,3,9 , where neighboring A-scans are averaged together prior to fitting. Lateral averaging can be an effective technique at reducing speckle variations but at a severe cost to lateral resolution. If the sample is not perfectly static, as is the case in liquid samples with particles undergoing Brownian motion or sufficiently dynamic living samples, consecutive A-scans taken at the same location can be averaged together to reduce speckle variations at the cost of effective acquisition time 10 . In either case, the layerwise fitting assumes complete uniformity in the composition and statistics of the layer segment in depth and lateral averaging makes the same assumption over a volumetric region. www.nature.com/scientificreports/ A depth resolved (DR) approach, initially developed for ultrasound image quantification 11 , was adapted by Vermeer for use in OCT and has become popular in recent years [12][13][14] . This approach removes the assumption of material uniformity in depth and allows variations in the attenuation coefficient in three dimensions. The DR approach assumes the material is weakly absorbing, making this technique related to voxelwise OCT scattering parameter inference methods [15][16][17] which have a long history in OCT signal processing. This method has been further refined by Liu 18 to better handle boundary effects caused by finite imaging depth. In either formulation, reconstructions of the attenuation coefficient will be highly variable due to the influence of speckle 14 . Thus, as before, lateral averaging is often still employed to get a more consistent result 12 . Conceptually, the DR approach allows one to recover some amount of the natural variability of optical properties within the tissues. While the advantages of the DR approach are manifest, the result of this approach in the presence of intensity variations due to speckle leads to reconstructions in which the recovered attenuation coefficient itself has large variations.
The propagation of speckle variation into the recovery of an otherwise deterministic coefficient has clear implications for the accuracy of the attenuation parameter inferred at a single voxel. Since the exact measured intensity is effectively random, one can in general expect the inferred coefficient to be effectively random as well. One way to handle the inference of parameters in these circumstances would be to adopt a Bayesian perspective. In this paradigm, instead of simply seeking an estimate for the value, one seeks the posterior distribution, which quantifies how probable each attenuation value is 19 . In these methods, accurate physical models about measurement uncertainty are combined with prior information about the objects which are being measured. Utilizing the posterior distribution allows for the identification of estimation biases and the quantification of uncertainty by giving access to statistics about the inferred attenuation coefficient. A better understanding of uncertainty can have direct clinical implications by helping to inform practitioners of how much they can trust a given inference. Furthermore, this approach opens the door to probabilistic tissue classification tasks such as tumor grading where the likelihood of various outcomes must be compared.
In this manuscript we model the effect of speckle on the inference of OCT attenuation coefficients using a Bayesian approach. The interaction between the DR reconstruction technique and the speckle variation is considered, and a probability distribution for the measurements made under physically realistic speckle variations is derived. Following this, we derive a prior distribution for a simple system using statistical mechanic principles. Finally, we combine these to derive a probability distribution for the attenuation coefficient itself and define a Bayesian voxelwise estimator for the mean attenuation coefficient. These results are then demonstrated in simulation and experimentally.
Paper structure. The goal of this work is to construct the posterior distribution for the voxelwise attenuation coefficient and to validate it using numerical experiments and tissue phantom measurements. The posterior distribution assigns a meaningful probability to every possible value of the true attenuation coefficient. Here, the true attenuation coefficient is defined as the attenuation coefficient of the mean OCT signal without speckle fluctuations. Using the existing DR method 12 , the attenuation coefficient at each voxel can be estimated from the measured OCT signal. These depth resolved estimates are denoted as μ . These estimates depend on the intensity at each voxel which fluctuates due to speckle. Because of these voxelwise fluctuations, the estimated value of the attenuation coefficient at that point will likely differ from the true coefficient. The posterior probability distribution gives the probability that the true value of the attenuation coefficient is equal to µ oct given that our depth resolved estimate was equal to μ.
Mathematically, the posterior distribution can be written as the conditional probability distribution P µ oct μ . Conditional probabilities can be rewritten as product of two easier to derive probability distributions using Bayes' theorem. This theorem states that the posterior distribution is given by, In this expression, P μ µ oct is called the likelihood function which represents the probability of estimating μ given that the true attenuation coefficient is equal to µ oct . The distribution denoted by P(µ oct ) is called the prior distribution for the unknown µ oct . The prior probability allows the incorporation of additional information into the statistical model and is often used as a way to establish bounds or to bias solutions towards realistic values. The marginal probability P μ is a normalizing factor and can be computed via integration. Using this relation, we can find the posterior distribution by solving two easier problems: finding the likelihood function and finding the prior distribution.
Before these two distributions can be derived we must first have a mathematical model for the measured OCT signal so we can make depth resolved attenuation coefficient estimates. In "Modeling intensity decay", a model which describes the mean signal decay is given. This model assumes that the measurements are made on a weakly absorbing medium and that the majority of measured light is single scattered. Next, in "A statistical model of the OCT amplitude and intensity", the effect of speckle on this OCT signal is considered and the probability distribution for the measurement is given. The likelihood function, P μ µ oct , is derived in "Analyzing the DR reconstruction distribution" by analyzing the speckle variations and is verified experimentally in "Experimental verification and results" by measuring the distribution of depth resolved attenuation coefficient estimates for a very homogeneous phantom. In "Constructing a prior distribution", the prior probability, P(µ oct ) , is derived using basic physical principles. This prior gives the background probability for finding a particular value of µ oct at any point in the sample without any additional measurement information. Following this we define a Bayesian estimator for the attenuation coefficient in "Bayesian parameter estimator". In "Simulation results" we simulate OCT signals with realistic variations to test our assumptions and statistical model.

Methods
Modeling intensity decay. In many practical OCT systems, the decay of the OCT intensity with depth can be adequately described using a single exponential decay model 5,20,21 , understanding the form of the OCT signal is necessary prior to understanding the attenuation itself. The attenuation coefficient is a material property, which depends on the absorption and scattering properties of tissue and is not a function of the measurement system. However, several system dependent factors can also contribute to measured signal attenuation such as the confocal point spread function and the sensitivity roll off function for OCT systems based on detection in the Fourier domain 4,22 . A model which takes all of these effects into account was described in detail in earlier work 20,22 . Typically for an OCT system, the signal decay due to the confocal PSF and the sensitivity roll off function can be independently measured, and subsequently, the resulting OCT data can be corrected for these effects.
For the sake of analysis, we will assume that the measured signal has already been calibrated for these system dependent effects. A more thorough discussion of this can be found in Supplemental Information S3. We denote the corrected OCT signal at depth z as I(z; µ b,NA (z), µ oct (z)) where µ b,NA (z) is the depth dependent back-scattering coefficient (the probability per unit length that light is back-scattered into the detection numerical aperture). The depth dependent attenuation coefficient, µ oct (z) , and the back-scattering coefficients depend on both the scattering coefficient µ s and the absorption coefficient µ a . These coefficients describe the probabilities of scattering and absorption per unit length, respectively. For weakly scattering samples, with negligible contributions from multiple scattered light, µ oct = µ s + µ a .
Following Vermeer 12 , we further assume that the tissue is very weakly absorbing ( µ a ≈ 0 ), and, a constant fraction of the attenuated light is back-scattered at every point in the tissue. We denote this fraction as β NA and define µ b,NA = β NA µ oct . Physically, this implies that the system is highly scattering dominant, i.e., there is very little absorbed light in the system when compared to the total attenuated light. Lastly, we assume the measurements are made with a fixed axial resolution denoted by z . Combining these assumptions the discretized quantity is defined which describes the mean value of the OCT signal with depth in a certain region at depth z = N z where N is the pixel index and given an incident intensity I inc . We use the shorthand µ oct (N) = µ oct (N�z) . Provided that the inverse of the attenuation coefficient is relatively small compared with the pixel size, its value is given by 12 As recently noted by Liu 18 , the tail of the series in the denominator in Eq. (3), meaning all of the terms in the sum after some large term K, can be computed when an estimate for an attenuation coefficient at that point in the sample is available. This is given by A statistical model of the OCT amplitude and intensity. The measured OCT signal is the amplitude of the backscattered field, which contains contributions from scatterers within the measurement volume, each contributing to the resulting field with their respective random amplitude and phase. This scattering results in OCT signal amplitude fluctuations called speckle. When there are sufficiently many scattering events within a single voxel, the speckle is called fully developed 23 and the measured signal becomes effectively random. In this case, the statistics for the signal amplitude, A, are well described by the Rayleigh distribution 8 given by where I denotes the mean intensity value. This formula gives the probability of measuring amplitude A when the mean signal is given by √ �I� . When OCT measurements are made, typically intensity is measured and not amplitude. Given a Rayleigh distributed amplitude of the form given in Eq. (5) it can be shown that the intensity 24 , which is the square of the amplitude, follows which is an exponential distribution with parameter I .
Analyzing the DR reconstruction distribution. This section considers the estimation of µ oct (N) from intensity measurements in the presence of speckle modeled by Eq. (6). In this case, instead of measuring the mean intensity I N directly we can only measure I N which is exponentially distributed with parameter I N . Because the constituent parts of Eq. (3) are now random, the estimate will be itself a random variable. The estimated random variable is denoted as . www.nature.com/scientificreports/ Following Vermeer 12 we consider the attenuation coefficient at the N th point and truncate the series in the denominator at M which in practice corresponds to the maximum imaging depth Z max with M > N giving Consider the denominator, and let, The variable D N is the sum of M − (N + 1) independent exponentially distributed random variables I i , taken from distributions parameterized only with average I i . Thus, D N will be distributed as a hypoexponential distribution and has mean because the I i 's are independent. If M is sufficiently larger than N, Eq. (4) implies that It is known that reconstruction artifacts 12,18 make the inferred coefficient unreliable near the deepest point of an A-scan. In practice, the reconstructed attenuation coefficient made from this approach must be discarded near the bottom of a scan and estimated using a different method 18 .
One useful measure of how much a random variable deviates from the mean called the coefficient of variation, and is denoted C v . This quantity is defined as the standard deviation divided by the mean. It can be shown that for a hypoexponential variable the coefficient of variation is always less than 1 as shown in the Supplemental Information S1. In practice, we find that C v ≪ 1 as demonstrated in Fig. 4 and described in detail in "Results".
Next, letting allows formula (8) to be rewritten as Because η has zero mean with a very small C v one can expect η N D N to be small. Using this as justification, consider the Taylor approximation At leading order, the reconstruction of the attenuation coefficient is given by Intuitively, this means the denominator of Eq. (8) is approximately constant at the scale set by the mean. Therefore, the probability distribution of μ will be given by rescaling the distribution of I N . Rescaling Eq. (6) yields Next, using the approximation for the tail of D given in Eq. (11) with K = N and substituting I N D N with µ oct (N) yields the probability distribution Therefore, the reconstructed coefficient at leading order will be exponentially distributed around the mean attenuation parameter. The accuracy of this estimate is demonstrated in Fig. 4.
This approach can be extended to the time-averaged case, where k independent co-registered measurements have been made. To do this, first the k estimates for the attenuation coefficient, denoted by . www.nature.com/scientificreports/ µ i (N), i = 1, 2m, . . . , k , should be constructed using Eq. (8). Then, assuming the measurements are independent, the likelihood is given by Constructing a prior distribution. In this section, a prior distribution for the variation in attenuation coefficient in a layer is derived based on physical principles. As an initial theoretical step we consider a simplified media of dispersed scattering particles with negligible absorption. Following Chandrasekhar 25 it is assumed that the system is a single layer, with N p dispersed particles throughout. Let be the ratio of the volume of a single voxel to the volume of the entire scanned layer. Provided that where k is the number of co-registered scans and µ oct is the layer mean of the DR estimates. The specific proportionality constant is given by integrating the numerator of Eq. (1) over all possible values of µ oct . Considering the case where only a single independent scan can be made the posterior distribution for the attenuation coefficient at depth N is given by This distribution describes the probability of the mean coefficient at voxel N. Assuming that each voxel is independent, a joint posterior distribution for the attenuation coefficient map for the entire A, B or C scan can be written as where R is the total number of voxels in the scan, µ µ µ oct is an R × 1 vector of true coefficients and μ μ µ is the R × 1 vector of voxelwise estimates for the attenuation coefficient. Figure 1 shows two posterior distributions plotted using Eq. (24) which use two different values for the DR estimate. These examples demonstrate the impact that the initial DR estimate has on the shape and position of the posterior distribution for the attenuation coefficient.
Bayesian parameter estimator. In Bayesian formulations of parameter estimation problems, when a single number prediction for the coefficient must be made, a Maximum a Posteriori (MaP) approach is often employed 14,26 . This approach gives the attenuation coefficient which maximizes the posterior distribution. However, as can be seen in Fig. 1 for sufficiently small DR estimates, the posterior distribution becomes bimodal and . www.nature.com/scientificreports/ the MaP estimate will nearly coincide with the low DR estimate for the attenuation coefficient. As demonstrated in Fig. 1b this peak is relatively narrow and contains little probability mass. Because of this, the maximum a posteriori is a bad representation of the entire probability distribution. The mean of the posterior distribution is agnostic to the bimodality of the distribution and provides a more stable and realistic estimate for the attenuation parameter. Therefore, when a single value estimate is desired, the quantity can be computed.

Experimental verification and results.
To verify the likelihood model from Eq. (17), the DR attenuation formula is applied to phantom data and a histogram is computed to compare against theory. The data was collected with a Santec IVS 2000 swept source OCT system with a central wavelength of 1309 nm, axial resolution of 12 micron in air and lateral resolution of 25.5 micron . The phantom was made by suspending silica beads manufactured by BaseClear with mean diameter of 0.47 micron and a refractive index of 1.425 in water at a volume fraction of 0.08. Water is assumed to have a phase refractive index of 1.32 and a group refractive index of 1.34 28 . Using Mie theory, the scattering cross section is given by 1.9 × 10 −9 mm 2 and the total attenuation coefficient is 3.2 mm −120 . This value is realistic for tissue 21,29 . An OCT B-scan of the phantom is shown in Fig. 2a. Using these values and Eq. (22) we can see that the expected variance for the attenuation coefficient is �µ oct � · ζ = 0.0020 mm −2 which is very small when compared with the variance of the exponential distribution which is �µ oct � 2 = 11.5 mm −2 . Since the speckle variance dominates the distribution of attenuation coefficients the reconstruction should look like Eq. (17). This is demonstrated in Fig. 2c. Figure 3 demonstrates the effect of the posterior mean estimator defined in Eq. (26) when compared with lateral averaging. Fig. 3a,b show the OCT attenuation coefficient B and A-scans respectively generated from the same OCT B-scan used in Fig. 2. This phantom is very homogeneous so we expect that the variation is almost entirely generated from speckle, thus it is reasonable to assume if sufficiently many A-scans are averaged together then the resulting attenuation coefficient should look constant. Figure 3d shows the resulting OCT attenuation coefficient after laterally averaging 1000 A-scans together. Figure 3c shows the result of the mean estimator defined in Eq. (26) applied to the A-scan from panel (b). There is little remaining variation in the signal when compared with standard lateral averaging.
Simulation results. To validate and better understand the statistical model from "Methods", a series of simulations were preformed. In Fig. 4a, a B  (b) This panel shows a constructed posterior distribution which is Bi-modal and has two local maxima. For a given layer mean, the constructed distribution develops a second peak if the DR estimate used to construct the posterior is sufficiently small. This second peak can make the Maximum a Posteriori difficult due to non-convexity. In many cases, the maximum value of the Posterior distribution may sit very near the origin on this second peak. As demonstrated in this panel, often the total amount of probability mass under the addition peak is relatively small, meaning that while the initial peak is overwhelmingly the maximum likelihood. Thus, the Maximum of the posterior distribution is a poor representation for the distribution itself. In these cases an estimate for the mean is a better choice. This posterior was constructed with a layer mean of �µ oct � = 0.4 mm −1 , ζ = 6.87 × 10 −2 mm −1 , and a DR estimate of μ = 0.015 mm −1 . www.nature.com/scientificreports/ coefficient of µ oct = 2.00 mm −1 . Once the deterministic signal is modeled we generate the OCT signal per voxel as a realization of an exponential random variable with parameter given by the true coefficient as in Eq. (6). This random realization can be seen in Fig. 4a. The attenuation coefficient was estimated using the DR method given in Eq. (8) and is shown in Fig. 4b. The reconstruction equation becomes inaccurate near the bottom of the measurement volume, preventing accurate estimation. To avoid these inaccuracies the deepest 30% of the reconstructed attenuation coefficients were truncated. The 30% value was arrived at by inspection. In Fig. 4c we fit an exponential model to the histogram of the reconstruction and see that the best fit parameter agrees with  www.nature.com/scientificreports/ our model to the 2nd decimal point. In Supplemental Information S2 we show that the truncation error from Eq. (14) leads to an error in the variance of μ on the same order as our fit error. To avoid artifacts the bottom 30% of the predicted attenuation coefficient is discarded. Figure 5c shows a posterior mean estimate for the attenuation coefficient which was computed with Eq. (26) voxelwise. In general, the mean attenuation coefficient for the layer, µ oct , would not be known ahead of time to compute the prior distribution. To account for this, we used the mean of the truncated DR attenuation estimate for the whole scan in Eq. (26). The estimate given by the mean of the posterior distribution for the attenuation coefficient can give much more accurate estimates for the true coefficient than using the standard DR technique, as demonstrated in Fig. 5.

Discussion
In this paper the impact of speckle fluctuations on the depth resolved recovery of the OCT attenuation coefficient has been addressed. When making an OCT measurement, effectively random voxelwise intensity fluctuations are present in the signal due to speckle, and as a result, the voxelwise mean attenuation coefficient can not be exactly determined. Utilizing a statistical understanding of speckle fluctuations and prior physical knowledge, the posterior distribution for the attenuation coefficient was derived from first principles. This probability distribution better characterizes the voxelwise attenuation coefficient because it allows for the weighing of relative likelihoods and the quantification of uncertainty by measuring the variance of the attenuation posterior distribution.
While the statistical framework derived in this paper is general, the applicability is limited by the assumptions made for the underlying depth resolved reconstruction technique. The DR reconstruction technique, given in Eq. (8), requires that the absorption of light be negligible when compared the total amount of attenuated light. This assumption is restrictive in the materials and wavelengths of light the DR technique can be applied to. However, for the materials and wavelengths used in most common biomedical applications of OCT this assumption is valid. Furthermore, when the probability distribution for the reconstructed coefficient in Eq. (17) was derived, it was assumed that the coefficient of variation of the denominator in Eq. (8) is sufficiently small such that the denominator can be treated as constant. This does seem to be valid in numerical simulations and experiments, however, it is not clear if this is generally true. www.nature.com/scientificreports/ Additional physical assumptions are made during the derivation of the prior distribution for the attenuation coefficient given in "Constructing a prior distribution". The prior distribution allows for the use of physical knowledge about the attenuation coefficient to introduce bounds and bias the probabilities towards realistic values. The derivation given in "Constructing a prior distribution" was made assuming the measured object contained uniform idealized scattering particles with no absorption. While this assumption may not hold for most tissue systems, a normally distributed prior is still a safe choice due to the fact that superpositions of random fluctuations tend to look normally distributed. In real tissue, the parameter ζ in Eq. (24) is difficult to define, as the meaning of the effective scattering cross section is ambiguous. However, it is still reasonable to assume that the true attenuation coefficient is normally distributed around the mean. The variance of the prior must be provided or inferred by other methods. There are techniques to estimate this parameter from the data such as empirical Bayesian methods 30 , however, the implementation of these techniques can be nontrivial and a robust verification must be performed before the method could be used clinically. While this is outside of the scope of this paper, the Bayesian model presented here serves an an initial step towards the goal of estimating these parameters more robustly in tissue, and elucidates the impact of speckle on the recovered coefficients.
The use of physically accurate statistical models for the attenuation coefficient has several potential advantages. The variance of the posterior distribution provides a way to quantify uncertainty in reconstructions. Furthermore, estimation bias from higher order moments of the posterior can be quantified as well. The likelihood ratio statistic 26 can be computed using the physically accurate likelihood function given in Eq. (17). This statistical test gives a practitioner a sense of how likely a parameter is to fall within a specified range. In situations where a practitioner may want to have a single number to understand the attenuation in a system, the mean of the posterior can be computed as demonstrated in Fig. 5. In Fig. 6 we measure the error in the estimates for both the DR and mean of posterior estimators as the scattering cross section and attenuation coefficient is varied.
Another potential application domain is in OCT image segmentation where attenuation analysis is used to correct for signal decay and as a contrast enhancement tool 13,31 . As we have discussed in this manuscript, the resulting attenuation image can be very highly variable due to the speckle fluctuations in the original signal. If the attenuation image is to be segmented, these fluctuations may lead to segmentation inaccuracies. Denoising algorithms could combine our exponential likelihood with a spatial priors, such as total variation 14 which would increase the likelihood of the piecewise constant attenuation coefficients. This could be used to improve segmentation accuracy by removing speckle fluctuations from the attenuation image. This approach may be applicable even in the case of absorbing media because image segmentation does not require extraction of accurate attenuation values, only sufficient contrast between layers.
This work is an initial theoretical step towards fully quantifying and characterizing uncertainties in voxelwise OCT attenuation coefficient recovery in order to better understand the resulting estimates. The likelihood function from Eq. (17) accurately models the voxelwise measurement uncertainty of the attenuation coefficient due to speckle. This likelihood function gives insight into the voxelwise statistics of the DR attenuation images. The posterior distribution for the mean value of the attenuation coefficient, given in Eq. (24), allows parameter estimation to be performed in a consistent and reliable manner by using the posterior mean estimator given in Eq. (26). Furthermore, the posterior distribution derived in this paper can be used to quantify the variance Figure 5. This figure shows estimates of the attenuation coefficient for simulated OCT data using the standard DR and the Bayesian estimator given in (26). The OCT data was simulated with parameters I inc = 1 × 10 7 , µ oct = 2.00 β NA = 0.3 , σ scat = 1 × 10 −6 mm 2 , a lateral resolution of x = 0.02 mm and z = 0.0068 mm in a domain which is 13.6 mm deep. After the attenuation coefficient was inferred using the DR method the bottom 30% of pixels are discarded to avoid reconstruction artifacts. Both the simulations and figure creation were done in Matlab 2019a 27 , https ://www.mathw orks.com/. (a) This panel shows the ground truth attenuation coefficient for the simulation. This ground truth is a realization of the prior distribution given in Eq. (22). (b) This image shows the reconstructed attenuation coefficient using the DR method given in Eq. www.nature.com/scientificreports/ in estimates, which gives insight into uncertainty. While this is a promising approach, further research is still needed to find the best way to apply these techniques to clinical practice. Clearly, the incorporation of of this value into the prior impacts our uncertainty in our Bayesian estimate.