Noise-bias and polarization-artifact corrected optical coherence tomography by maximum a-posteriori intensity estimation

: We propose using maximum a-posteriori (MAP) estimation to improve the image signal-to-noise ratio (SNR) in polarization diversity (PD) optical coherence tomography. PD-detection removes polarization artifacts, which are common when imaging highly birefringent tissue or when using a ﬂexible ﬁber catheter. However, dividing the probe power to two polarization detection channels inevitably reduces the SNR. Applying MAP estimation to PD-OCT allows for the removal of polarization artifacts while maintaining and improving image SNR. The e ﬀ ectiveness of the MAP-PD method is evaluated by comparing it with MAP-non-PD, intensity averaged PD, and intensity averaged non-PD methods. Evaluation was conducted in vivo with human eyes. The MAP-PD method is found to be optimal, demonstrating high SNR and artifact suppression, especially for highly birefringent tissue, such as the peripapillary sclera. The MAP-PD based attenuation coe ﬃ cient image also shows better di ﬀ erentiation of attenuation levels than non-MAP attenuation images.

Although the contrast properties of these functional OCTs are well investigated, the property of the most basic contrast, i.e., OCT intensity, had not been extensively investigated. In a recent publication, we have highlighted a commonly encountered problem due to low SNR in standard OCT B-scans [23]. For example, in posterior eye imaging, the vitreous body appears with a weak random signal, where one expects a much lower signal because its purpose is to be transparent.
A conventional method to improve image quality involves repeatedly acquiring multiple Bscans at the same location, and then averaging the signal intensity [24,25]. Although this intensity averaging is effective for reducing speckle contrast and/or improving image SNR, it does not reduce bias [23], which we denote as "noise-offset". To compensate for this noise-offset, histogram equalization [26] or a simple intensity threshold may be applied, but this would sacrifice the quantitative nature of the signal.
Another commonly encountered problem is the appearance of polarization artifacts in standard OCT images, especially if the OCT system is equipped with a flexible fiber probe [27] or if it images highly birefringent tissue, such as the peripapillary sclera [24,28]. Polarization diversity (PD) detection can remove polarization artifacts by summing the OCT signal intensities of the vertical and horizontal polarization detection channels [28,29]. However, dividing the signal power into two detection channels inevitably reduces the sensitivity. This is because the noise power is doubled when using two detectors instead of one. Therefore, PD-detection still results in a reduction of image SNR, compared with standard OCT.
Hence, we have identified two problems in OCT imaging which we wish to solve, i.e, noiseoffset and the corruption of images by polarization artifacts.
To solve these two problems, we present an image composition method based on PDdetection to suppress polarization artifacts along with a maximum a-posteriori (MAP) intensity estimation to reduce noise-offset while preserving the quantitative nature of the signal intensity.
MAP estimation has been applied to OCT intensity [23] and phase estimation [30], and to birefringence estimation in polarization sensitive OCT [31][32][33]. Among these, MAP OCT intensity estimation allows for estimation of signal intensity by utilizing the signal and noise statistics. If used on the image composition of repeatedly obtained B-scans, it provides an image with less noise-offset than images obtained from conventional intensity averaging methods [23]. In addition, because the estimated intensity has less noise-induced offset, it is suitable for further statistical analysis. The additional information provided by the precision and reliability of the MAP intensity and attenuation coefficient estimates also enhances the quantitative accuracy of further analysis of OCT signals. For example, we can rationally reject untrustable signals.
Here, we combine PD-OCT with MAP estimation. This combination will provide a more quantitative estimation of the total light energy. This method is compared with combinations of non-PD-OCT and intensity averaging, PD-OCT and intensity averaging, and non-PD-OCT and MAP estimation. Non-PD-OCT is emulated by the coherent composition of two PD-detection channels. These comparisons show the superiority of the MAP estimation methods.
In addition, we demonstrate depth-localized attenuation coefficient imaging [34] based on the OCT intensity estimated by MAP estimation. We show that the attenuation coefficient estimation from MAP intensity images results in a broader dynamic range and better differentiation of estimated attenuation levels compared with attenuation images derived from averaged intensity images.

MAP estimation of OCT amplitude and intensity
In this section, we describe the theory of MAP estimation of OCT intensity. This theory was previously described in [23], but repeated for completeness. More details are also provided here on the implementation of the algorithm. We first describe a MAP estimation of the OCT amplitude, and then show that the square of the MAP estimate of the amplitude is equivalent to the MAP estimate of the OCT intensity.

MAP estimation of OCT amplitude
To reduce noise-offset and quantitative nature of the PD-OCT method, we utilize MAP estimation of the OCT signal intensity [23]. The MAP estimation method uses the statistics of the OCT signal and noise. It is assumed that both the real and imaginary parts of an OCT signal are affected by independent and identically distributed (i.i.d.) additive white Gaussian noise. Hence, the OCT signal amplitude is modeled by a Rice distribution. The probability density function of the observed OCT signal amplitude, a, given the "true" signal amplitude, α, is given by [35] p a | α, σ 2 = a σ 2 exp where σ 2 is the variance of the real or imaginary parts of the additive white Gaussian noise, which are equal in value by assumption, and I 0 is the 0-th order modified Bessel function of the first kind. In practice, the depth-dependent noise variance σ 2 is estimated from the noise data, which is obtained by taking A-scans while obstructing the probe beam. The noise variance, σ 2 , is defined as the average of the variances of the real and imaginary parts. Note that in some of the existing literature, the symbol of σ 2 denotes a noise energy which is the sum of the variances of the real and imaginary parts. In contrast, we follow the notation introduced by Goodman [35]. By treating the underlying true amplitude α as a variable, and the observed values a and σ 2 as parameters, the likelihood of the true signal amplitude under specific observations of signal amplitude a and noise variance σ 2 can be expressed as: The combined likelihood function for a set of independent and identically distributed (i.i.d) measured amplitudes a = {a 1 , · · · , a n , · · · , a N }, which were obtained from repeated B-scans at the same location, is given by Therefore, MAP estimation of the true signal amplitude from this set of measurements is given as the value of α which maximizes the posterior distribution, l α; a, σ 2 π(α), where π(α) is the prior distribution of the true amplitude. In our case, the prior distribution is assumed to be uniform (non-informative).

MAP estimation of OCT intensity
By defining the true value of the OCT intensity as υ = α 2 , the MAP estimation of the intensity is expressed asυ = arg max υ l υ υ; i, σ 2 π(υ) where l υ υ; i, σ 2 is the likelihood function of the intensity and π(υ) is the prior distribution of the true intensity. The set of measured intensity is square of amplitudes i = {i 0 , · · · , i n , · · · , i N } = {a 2 1 , · · · , a 2 n , · · · , a 2 N }. If we assume a uniform prior π(υ), it can be shown that A proof of this is given in Appendix A. For convenience of implementation, we first compute the MAP estimate of the amplitude, and then square it to obtain the MAP estimate of the intensity.

Reliability of MAP estimation
From the likelihood ratio statistics of the combined likelihood [Eq. (4)], one may also obtain the 68% credible interval of the amplitude estimation [36]. The likelihood ratio statistic T (α) is given by Theory states that this test statistic has a χ 2 1 -distribution (i.e., a χ 2 -distribution with one degree of freedom) [36], hence the 68% credible interval is given as the region of α where T (α) ≤ 0.99. Half of this interval is an approximate estimate for the standard deviation of the amplitude estimation σ α , and the amplitude estimation error is defined as σ 2 α . The reciprocal of the estimation error 1/σ 2 α is taken as the precision of the amplitude estimates. Note that this precision measure only accounts for the estimation variance, not the estimation bias.
A similar method of credible interval calculation was used in our earlier publication [23], but it was defined by a different threshold T (α) ≤ 3.84, which provides the 95% credible interval. We have changed the 95% credible interval to the 68% credible interval because it is a good approximation to the two standard deviation interval of the MAP amplitude estimate [37,38].
The precision of the MAP intensity estimation is then computed by error propagation, in which σ 2 α is propagated from the amplitude to the intensity. In particular, the uncertainty of the MAP intensity estimate (σ υ ) is defined as σ υ ≡ 2ασ α . The precision of the MAP intensity estimate is then defined as 1/σ 2 υ . Higher intensity regions are expected to have higher fluctuations, and hence lower precision. However, this precision is mainly dominated by the intensity itself, and is not a good measure of the estimate reliability. Therefore, it is informative to also define an estimate reliability measure as a squared-intensity-to-error ratioυ 2 /σ 2 υ . In decibel scale, it is expressed as 20 log 10 (υ/σ υ ). It is also noteworthy that the intensity reliability is proportional to the squared-amplitude-to-error ratio according to:υ 2 /σ 2 υ =α 4 / 4α 2 σ 2 α =α 2 / 4σ 2 α .

Numerical implementation of the probability density function
Note that because the zeroth order Bessel function of the first kind, I 0 (z), is of order O exp(z 2 ) , it is numerically divergent and cannot be used in a numerical implementation.
To overcome this problem, we use the exponentially-scaled modified Bessel function of the first kind. Therefore, in our numerical implementation, the probability density function, Eq. (1), is given by p a n | α, σ 2 = a n σ 2 exp −(a 2 n +α 2 ) 2σ 2 I 0 a n α σ 2 = a n σ 2 exp −(a n −α) 2 −2a n α 2σ 2 I 0 a n α σ 2 = a n σ 2 exp −(a n −α) 2 2σ 2 I 0 a n α σ 2 exp − a n α where the final square-bracket component comes from the exponentially scaled modified Bessel function of the first kind. We numerically compute the first exponential part and the part in square brackets independently, then multiply them. The same algorithms of lookup table generation and peak searching were used in this manuscript compared with our previously published method [23]. However, the algorithm was newly implemented in Python 2.7.11, while it was implemented in Matlab 2014b in [23].
The computation time for a 500-amplitude-level by 200-noise-level lookup table is 14.1 s for lookup table generation and 27.3 s for the MAP estimation using a Windows 10 PC, with an Intel Core i7-5930K processor and 32GB of RAM.

OCT image compositions
We assume that multiple, N, B-scans are obtained at the same location of the sample by a PD-detector, resulting in 2N frames. That is, N B-scans × 2 PD-detection channels. In this section, we describe four image composition methods used to create a single composite image from the 2N frames. The main purpose of this section is to present a polarization-artifact-free, high-contrast OCT image composition method by using PD-detection and MAP intensity estimation (MAP PD-OCT or MPD). This composition method is described in Section 2.2.2. In addition, the standard composition methods also presented are: first, a polarization-artifact-free OCT based on intensity averaging (standard PD-OCT or SPD, Section 2.2.1); second, intensityaveraging combined with coherently combined PD-detected signals (standard non-PD image or SnPD, section 2.2.3); third, MAP estimation combined with coherently combined PD-detected signals (MAP non-PD image or MnPD, Section 2.2.3).

Standard PD-OCT
The PD-detection method uses two complex OCT signals from two orthogonal detectionpolarizations E h (z) and E v (z), where the subscripts h and v denote horizontal and vertical polarization, respectively. A standard PD-OCT (SPD) image is obtained by averaging N frames for each PD-detection channel and summing the averaged frames as: where the over-line indicates averaging over the N frames, and the subscript SPD is for standard PD. Because the optical energies of the two detection polarizations are summed in this image, I SPD is free from polarization-artifacts. On the other hand, intensity averaging along the frames results in significant signal bias in low-signal-intensity regions, which is denoted by noise-offset in this manuscript.

High-contrast PD-OCT by MAP intensity estimation (MAP PD-OCT)
MAP intensity estimation can remove the noise-offset found in standard PD-OCT, while retaining the polarization-artifact-free state of PD-OCT. High-contrast polarization-artifact-free PD-OCT can be obtained by: where the hat represents the MAP intensity estimate over frames [Eq. (5)], and the subscript MPD is for MAP-PD. The estimation error of I MPD is the summation of the estimation errors of |E h (z)| 2 and |E v (z)| 2 , so the precision of I MPD is defined as its reciprocal, 1/σ 2 υ,h and σ 2 υ,v are the intensity estimation errors of the horizontal and vertical detection polarizations (see the last paragraph of Section 2.1.3). The reliability is then defined as I 2 MPD /σ 2 MPD .

Standard and MAP non-PD-OCT
Non-PD OCT can be emulated using OCT signals obtained by PD-detection. A single frame of a pseudo-non-PD-OCT image is obtained by coherent composition [17] of OCT signals from the two orthogonal polarization channels according to: where θ is a depth-independent relative phase offset, defined as where z is the pixel depth. As is evident in this equation, the coherent composition is a complex average with adaptive phase correction. Because coherent composition is effectively complex averaging, it suppresses noise and improves the SNR.
The N frames from the non-PD-OCT can be combined either by intensity averaging (standard contrast) or MAP intensity estimation. The standard, non-PD-OCT contrast is obtained by averaging the non-PD-OCT frames in their intensity as: where the subscript SnPD is for standard non-PD. This image suffers from noise-offset and polarization artifacts, but would have a higher SNR than would standard PD-OCT. High-contrast non-PD-OCT is obtained by combining the non-PD-OCT frames using MAP intensity estimation: where the subscript MnPD is for MAP non-PD. This image has a low noise-offset because of the MAP intensity estimation. Although it is affected by polarization artifacts, the noise-suppression effect of complex averaging provides a higher SNR in comparison to MAP PD-OCT. The properties of the four composition methods are summarized in Table 1.

Attenuation coefficient calculation
Attenuation coefficient images are generated for each of the four types of composite images by applying a method previously presented by Vermeer et al. [34]. Here, the depth-dependent attenuation coefficient μ is computed as where z i is the depth of i-th depth pixel, I (z i ) is the intensity of the composite OCT image, Δ is the inter-pixel distance, and M is the number of pixels per A-line.

Signal-roll-off correction
For attenuation coefficient estimation, the OCT intensity is corrected to account for the depthdependent sensitivity roll-off. The depth-dependent SNR (SNR(z i )) was measured using a mirror sample and a neutral density (ND) filter at approximately each 275-μm depth interval from 0 to 3 mm in air (or each 200-μm depth interval in tissue). The depth-dependent signal decay curve C(z i ) is then computed from the SNR and the depth-dependent noise energy (σ 2 (z i )) as: This signal decay curve is transformed to logarithmic scale, then fit by a quadratic function, and is used as a correction factor. Two correction factors are obtained independently for the two PD-detection channels, referred to as C h (z i ) and C v (z i ) for the horizontal and vertical channels, respectively. Another correction factor, C nPD (z i ), is obtained from the coherent composite (non-PD) OCT signal. The PD-OCT signals are corrected by using C h (z i ) and C v (z i ) in: where |E h (z i )| 2 and |E v (z i )| 2 represents intensity averaging or MAP estimation of |E h (z i )| 2 . From Eq. (16), we may then obtain I PD (z i ), which are the corrected standard or MAP PD-OCT intensities, and are substituted into I (z i ) of Eq. (14). Note that this signal-roll-off correction was performed only for attenuation coefficient imaging, but not for standard OCT imaging. It should be noted that an estimator with noise floor subtraction is applied in Ref [34]. However, the subtraction is not applied to the intensity averaging estimator here. To make consistent with intensity imaging comparison, attenuation calculations based on average-only and MAP estimates are used.

Precision and reliability of the attenuation coefficient estimation
As we obtain the error from the MAP intensity estimate, the error in intensity can be used to calculate the attenuation coefficient precision. This calculation is performed by the method of error propagation.
The estimation error of the attenuation coefficient is computed by using the error propagation method based on Eq. (14), which relates the MAP intensity estimation error, σ 2 υ , to the attenuation coefficient error, σ 2 μ , by where σ 2 υ (z i ) is the estimation error of the OCT intensity defined in Section 2.1.3. According to Eq. (14), the partial derivatives in this equation are evaluated as ∂μ(z i )/∂I (z i ) = μ(z i )/I (z i ) and ∂μ(z i )/∂I (z k ) = −μ(z i )/ M k=i+1 I (z k ) for i j. Therefore, the estimation error of the attenuation coefficient at depth z i can be expressed by The first term in the equation can be interpreted as the error contribution from the pixel of interest, and the second term is the error contribution due to all of the pixels below it. Using this relation, one may calculate the attenuation coefficient precision, 1/σ 2 μ(z ) . Because the precision is mainly dominated by the attenuation coefficient, it is also informative to calculate the squared-attenuation-to-error ratio μ 2 (z i )/σ 2 μ (z i ). This may also be expressed according to decibel scale as 20 log 10 μ(z i )/σ μ (z i ) . This ratio can help determine the regions in which the attenuation image is not conveying any meaningful information, so it is a measure of the reliability of the attenuation estimation.

Jones-matrix OCT system
We used a 1.06-μm multifunctional Jones-matrix OCT for the experimental study. This system uses a wavelength swept laser light source (AXSUN 1060, Axsun Technology Inc., MA, USA) with a center wavelength of 1,060 nm, a scanning bandwidth of 123 nm, full-width-athalf-maximum bandwidth of 111 nm, and a scanning rate of 100 kHz. Two incident (probe) polarization states are multiplexed by a polarization-dependent delayer, in which one polarization travels a longer optical path than the other. As a result, the OCT signal of the delayed input polarization appears farther from the zero-delay position than does the input polarization which was not delayed. The interference signal detection is performed by a PD-detection module. In this way, the orthogonal output polarizations are measured by two independent dual-balanced photodetectors, and we obtain four OCT images which correspond to each entry of the Jones matrix.
The axial resolution and axial pixel separation in tissue are 6.2 μm and 4.0 μm, respectively. The system sensitivity was measured to be 91 dB. 490 A-lines are taken per B-scan and the image is truncated into 480 depth-pixels per line.
More details of the hardware and software of this system are described in Refs.

PD and non-PD image formation from Jones-matrix OCT
The PD image composition methods described in Section 2.2 are based on PD-detection, which provides two OCT signals with two orthogonal output (detection) polarizations. On the other hand, the JM-OCT system provides four OCT signals, because it also multiplexes the two input polarization states. To apply the image composition methods, we need to create two OCT signals, that emulate the PD-detection signals from the four OCT signals.
To emulate non-Jones matrix PD-detection, we combine the two multiplexed incident polarization signals. The mutual phase difference between two incident polarizations is first estimated by: where the summation is over all pixels in an A-line. This equation represents two independent equations for θ h and θ v as identified by the subscript. For θ h , the subscripts "h, v" should be read as h, while for θ v it should be read as v. We use the same convention in the subsequent equations. E (1) h,v and E (2) h,v represent OCT signals obtained from the first and the second incident polarizations, as indicated by the superscripts (1) and (2), respectively. The reconstructed polarization diversity signals are then given by: In Section 2.2.3, we describe that the pseudo-non-PD-OCT is obtained as a complex composition of two PD-OCT signals [Eq. (11)]. However, in this particular study, the pseudo-non-PD-OCT is obtained by directly applying the complex composition of the four OCT signals of JM-OCT, as described in Section 3.6 of Ref.
[17], rather than Eq. (11). The reason is that these two methods are theoretically equivalent, and the latter is computationally less intensive.

Measurement and validation protocol
To evaluate and compare the performances of the composition methods, the (right) optic nerve head (ONH) and (right) macula of a 29-year-old male subject were imaged. This subject was without any marked disorders except non-pathological myopia (-7.45D spherical equivalent). Four repeated B-scans, with lateral widths of 6.0 mm, were taken at a single position, and the four types of compositions were created. We have previously shown that a four-fold scan is a good compromise between image quality and acquisition speed for systems of this sensitivity range and acquisition rate [23].
The compositions were generated by a custom-made program written in Python 2.7.11 with numerical computation packages Numpy 1.9.2-1 and Scipy 0.15.1-2. The signal intensity ratio (SIR) between the retinal pigment epithelium and the vitreous were computed as a metric for the performance evaluation of the composition methods.
The data acquisition protocol adhered to the tenets of the Declaration of Helsinki, and was approved by the institutional review board of the University of Tsukuba. Figure 1 shows the intensity images from the macula and ONH. The images are composed using intensity averaging in Fig. 1(a)-1(d) (first row, standard non-PD-and PD-OCT in Table 1) and MAP estimation in Figs. 1(e)-1(h) (second row, MAP non-PD-and PD-OCT in Table 1 By comparing the images made by intensity averaging (first row) and the MAP images (second row), it can be seen that the intensity averaged images show a higher noise-offset in the low intensity regions, and have lower contrast compared with the corresponding MAP images.

Intensity imaging
The SIR is measured to be 6.7 dB higher in the MAP PD-OCT image [ Fig. 1(h)] than in the corresponding averaged image [ Fig. 1(d)]. This indicates that MAP estimation over repeated frames is more effective in improving image contrast than intensity averaging over the same number of repeated frames. This is also evident qualitatively, as the MAP-estimated images in the second row appear with better contrast than the averaged images in the first row. Hence, we conclude that MAP estimation provides better contrast than intensity averaging.
By comparing non-PD (left four) and PD (right four) images, it is evident that PD-detection and image composition suppresses polarization artifacts significantly. For example, non-PD images show polarization artifacts in the peripapillary sclera of the ONH, denoted by red arrows in Figs. 1(b) and 1(f), while they are strongly suppressed in the PD images [ Figs. 1(d) and 1(h)].    Fig. 1(h). Black means the corresponding pixel has intensity larger than or equal to -50 dB and white means less than -50 dB. There are many pixels in the MPD ONH image with intensities less than -50 dB. They are located in both the vitreous and deep regions.
On the other hand, the non-PD-OCT images show slightly better contrast than the corresponding PD-OCT images. The non-PD MAP macula image [ Fig. 1(e)] has a SIR around 2.7 dB higher than the PD-OCT MAP image [ Fig. 1(g)]. The greater SIR of the non-PD image is due to the complex averaging of four Jones elements, rather than just two Jones elements as in PD-OCT.
For PD-OCT, the probe beam power is split into two polarization detection channels, while in standard (non-PD) OCT, it is not. Dividing the power evenly into two detection channels results in a 3-dB sensitivity loss, because the noise power is doubled when using two detectors instead of one.
Hence, PD-OCT suffers from a SNR penalty compared with non-PD-OCT, although it can effectively suppress polarization artifacts. Figure 2 shows contour plots of 2D histograms of MAP and average intensities of corresponding pixels in non-PD [ Fig. 2(a)] and PD [ Fig. 2(b)] images of the ONH. The red line is the line of equal intensities for MAP estimation and intensity averaging. MAP estimation intensity is broadly spread at low averaged intensity. It indicates that the MAP composition method can estimate far lower intensities than the intensity averaging. The cluster of pixels at the MAP intensity of -85.8 dB in Fig. 2(a) and of -87.5 dB in Fig. 2 Figure 3 shows the image histograms, in which each histogram corresponds to the images at the same location in Fig. 1 Fig.  3 also confirm the observations in Fig. 1. That is, significantly lower signal intensity is found with MAP estimation (second row), compared with averaging (first row). As shown in Fig. 4, when back-projecting the low intensity pixels (pixels with intensities less than -50 dB) to their spatial locations, we find that the low intensity pixels in the MAP intensity images are broadly distributed in the vitreous and deep regions. According to the histograms of the averaged images (Fig. 3, first row), the low intensities appear to be shifted up. This suggests that a large estimation noise-offset exists in the low intensity regions in the averaged images. Among the averaged images (first row), the upward shift was slightly higher in PD-OCT (right two) than in non-PD-OCT (left two). This is because the number of intensity-averaged frames used to form a PD-OCT image is twice the number for non-PD-OCT. Namely, each of the two frames of the non-PD-OCT image is formed by the complex averaging (coherent composition) of two frames, rather than intensity averaging. Complex averaging does not result in a signal shift in low intensity regions.

. The histograms show that applying MAP [Figs. 3(e)-3(h)] results in a lower noise-offset than intensity averaging [Figs. 3(a)-3(d)]. The histograms in
The single peaks seen in the histograms of the MAP intensity images at the low intensity values are situated at the lowest possible estimable value in that particular numerical estimation scheme. For our purposes, those corresponding pixels can be considered to have zero intensity. The same issue has also been discussed in our previous publication [23].
When examining the maps of intensity estimation precision [Figs. 5(a)-5(d)], it can be seen that high intensity regions, such as the retinal pigment epithelium (RPE) have low precision (with large error), while the low intensity regions, such as the vitreous, have a high precision (and small errors). This is logical, because the higher intensity regions are expected to have higher intensity fluctuations. On the other hand, the reliability of estimation, as measured by the squared-intensity-to-error ratio, [Figs. 5(e)-5(h)] shows that the estimation reliability is usually higher in the higher intensity regions. Figure 6 shows contour plots of 2D histograms of MPD intensity and reliability [ Fig. 6(a)], and MPD intensity and precision [ Fig. 6(b)] of the macula. In general, there is positive correlation between intensity and reliability, and negative correlation between intensity and precision. There is a large cluster of pixels with the intensity of -86.7 dB (red arrows). -86.7 dB is the predefined minimum value of the estimation, i.e., all pixels having an intensity equal or lower than this value are regarded as having this minimum value. Because of the low SNR, these pixels also have a low reliability. However, as the intensity values are low, the fluctuation of the original OCT intensity is low, so the estimated precision is high. There are also 98 pixels that have high intensity (around 0 dB of MAP intensity), high reliability (greater than 70 dB), and intermediate precision (50 dB) as indicated by the black arrows in Fig. 6. By back-projecting these pixels in the original image, it was found that these are isolated high intensity pixels located at the retinal surface and RPE.

Attenuation imaging
Attenuation coefficient imaging provides information of the light-scattering properties of the tissue, rather than the back-scattered intensity information obtained from standard OCT intensity images. We computed the attenuation coefficient images from the ONH and the macula (Fig.  7) from the composite intensity images described in Section 2.2 by the method described in Section 3 [34].
It can be seen that, for the images of the vitreous shown, the MAP attenuation images (Fig.  7 second row) have lower estimated attenuation coefficients than the intensity averaged attenuation images (Fig. 7, first row). Hence, the MAP attenuation images show higher contrast and dynamic range than do averaged attenuation images.
Among the MAP attenuation coefficient images (Fig. 7, second row), the non-PD images non-PD-OCT can be a good choice for imaging regions with no or little birefringence.
There appears to be vertical line artifacts in the averaged images (first row). These line artifacts are less apparent in the MAP images. Figure 8 shows contour plots of 2D histograms of attenuation coefficient values from the MAP estimation and intensity averaging at corresponding pixels in the ONH images. The red line indicates the pixels of equal attenuation between the MAP and averaging estimation methods. It can be seen that the MAP attenuation images have a wider dynamic range than the attenuation images calculated from intensity averaged images.
In both Non-PD [ Fig. 8(a)] and PD [ Fig. 8(b)] image composition methods, there is a high number of pixels in agreement, as shown by the large counts on the equi-attenuation line. However, there is also a large cluster of pixels that the MAP estimation method estimates lower attenuations, indicated by the region below the equi-attenuation line. For a large number of pixels in the vitreous and deep region, the MAP estimator estimated values less than 100 times lower than the intensity averaged method.
In Fig. 8(a), there is a small cluster of 274 pixels that have a MAP estimated attenuation coefficient values greater than 10 4 mm −1 and intensity averaged attenuation coefficient values greater than 1 mm −1 . These pixels come from the deepest location in the attenuation images [Figs. 7(b) and 7(f)]. The estimated attenuation coefficients at this location are larger than those of above region. This overestimation is pointed out by Vermeer et al. [34] as violation of the assumption of the attenuation reconstruction theory. Because MAP estimation estimates lower intensity compared to averaging as presented in Section 5.1, this probably makes the overestimation of attenuation more prominent. As shown in Figs. 10(b) and 10(f), MAP estimation precisions and reliability are low at the deepest location. Hence, the MAP estimation precision and reliability may be able to treat this estimation error. There is a similar number of artifactual pixels (160 pixels) as shown in Fig. 8(b), which also come from the deepest location in the PD attenuation images, [Figs. 7(d) and 7(h)].
The histograms, Fig. 9, do not show a significant difference between the non-PD (left four)  and PD-OCT images (right four). However, it is evident that the averaged images (first row) have a reduced dynamic range and poor discrimination between attenuation coefficient levels, compared with the MAP images (second row). The MAP image histograms show a broader dynamic range of attenuation coefficients, and more numerous and better defined peaks. The precision of the images [Figs. 10(a)-10(d)] suggests that the non-PD images [Figs. 10(a) and 10(b)] have a slightly higher precision than the PD images [Figs. 10(c) and 10(d)], especially in the low intensity regions, such as the vitreous. The precision decreases with depth due to the reduced number of pixels used for estimation, and the increasing error contribution from the second term in Eq. (18). The reliability (squared-attenuation-coefficient-to-error ratio) [Figs. 10(e)-10(h)] maps show that the reliability is higher where the signal strength and signal SNR are higher. By using these reliability maps, we can conclude that the attenuation coefficients in the deep regions are not reliable.

Conclusion
PD-detection and image composition removed polarization artifacts that were apparent in OCT images of the peripapillary sclera. On the other hand, non-PD-OCT images show a slightly higher SIR than PD-OCT images, but also contained polarization artifacts. The images composed by MAP estimation always show better image contrast than the corresponding intensity averaged images. The combination of MAP composition and PD-detection is successful because it can compensate for the reduction in SNR caused by the division of probe power during PD detection, while still suppressing polarization artifacts. In light of these results, we conclude that the combination of MAP and PD-OCT is the best choice for birefringent samples, such as the ONH. In contrast, the combination of MAP and non-PD-OCT is a good option for less birefringent regions, such as the macula. One of the important purposes of this study is to obtain accurate attenuation coefficient values, which is a quantitative measure of the optical property of tissue. As the attenuation coefficient is based on the backscattered light intensity, quantitative light intensity information is required. Two problems then arise that hinder the acquisition of quantitative light intensity: noise-offset in low OCT intensity regions and polarization artifacts.
In this paper we described an image composition method that combines polarization diversity detection and MAP estimation, to remove noise-offset and polarization artifacts. By applying model-based attenuation coefficient reconstruction to quantitative light intensity [34], one may obtain fully quantitative attenuation.
The resulting quantitative intensity images, with noise-offset and polarization artifact correction, provide superior contrast for subjective observation compared with conventional OCT. Moreover, the quantitative attenuation images computed from the quantitative light intensity provide a more accurate estimation of the tissue optical properties. This is especially important for quantitative or automated diagnosis.