Noise-immune complex correlation for optical coherence angiography based on standard and Jones matrix optical coherence tomography

: This paper describes a complex correlation mapping algorithm for optical coherence angiography (cmOCA). The proposed algorithm avoids the signal-to-noise ratio dependence and exhibits low noise in vasculature imaging. The complex correlation coefﬁcient of the signals, rather than that of the measured data are estimated, and two-step averaging is introduced. Algorithms of motion artifact removal based on non perfusing tissue detection using correlation are developed. The algorithms are implemented with Jones-matrix OCT. Simultaneous imaging of pigmented tissue and vasculature is also achieved using degree of polarization uniformity imaging with cmOCA. An application of cmOCA to in vivo posterior human eyes is presented to demonstrate that high-contrast images of patients’ eyes can be obtained. indicates that the ﬂow signals are blocked (B). The cross-sectional image of the combination of OCT intensity (coherent composite + sensitivity-enhancement with global phase correction) and cmOCA reveals the relationship between morphology and vasculature. The log-scaled OCT intensity

Despite these advantages, complex-correlation-based OCA still suffers from certain problems. First, the correlation coefficient is not only affected by motion, but also by the signal-to-noise ratio (SNR) [29]; i.e., low SNRs create artificially high flow signals. A signalintensity-based threshold has been widely used to remove the low-SNR regions [20, 24, 27, 30], but the correlation signals above the threshold are still affected by the SNR, and this approach is not fully quantitative. The second point is the variance of the estimated complex correlation coefficient. This is limiting the contrast of vasculature imaging. For in vivo imaging, the involuntary motion of samples causes decorrelation and artifacts.
In this paper, a new OCA method, which we call correlation mapping OCA (cmOCA), is derived using the SNR-corrected low-noise complex correlation. We present SNR-correction and variance-reduction theories for the complex correlation coefficient estimation (Section 2). The treatment of bulk tissue motion is described in Section 3.1 and Section 3.2. The cmOCA method is applied with Jones matrix OCT (JM-OCT) [31-37], which is a type of polarizationsensitive OCT that has multiple polarization detection channels and enables polarizationsensitive contrast. Details of the implementation of cmOCA with JM-OCT are presented in Section 4. The combination of degree-of-polarization uniformity (DOPU) and cmOCA for simultaneous visualization of vasculature and the retinal pigment epithelium (RPE) in the eye (Section 5.3) is also described. The cmOCA method is applied to normal and pathologic eyes.

Theory of SNR-corrected low-noise complex correlation coefficient estimates
In this section, the general complex correlation coefficient of OCT signal and its estimate are derived (Section 2.1). The SNR-corrected low-noise estimator of the complex correlation coefficient is then described (Section 2.2). The notations used here and the following sections are listed in Appendix B.

Complex correlation of OCT signals
By assuming that the OCT signal at the spatiotemporal point r = (x, z,t) (where x is the transversal position, z is the axial position, and t is time) is a random variable G(r), the covariance matrix between two OCT signals at r and r + Δr can be expressed as where E [ ] represents the expectation, and E [G(r)] is assumed to be zero. The population correlation coefficient ρ G and phase offset θ G between G(r) and G(r + Δr) are The correlation coefficient ρ G is estimated using the sample covariance matrix S G between the measured OCT signal pair g(r) and g(r + Δr) which is obtained as where g(Δr; r) = g(r) g(r + Δr) , and · w is an averaging operation across a window w. The superscript † denotes the conjugate transpose operation. The correlation coefficient estimate of ρ G is then obtained as where s Gpq are the entry of the estimated covariance matrix S G at p-th row and q-th column.
In practical OCT measurements, the correlation coefficients obtained from measured OCT signals are not only affected by the true OCT signal, but also by noise. The OCT signal to be measured at the spatiotemporal point r can be considered as the summation of the true OCT signal S(r) and some random additive noise N: By assuming S(r) to be a zero-mean process and N to be a zero-mean stationary random process, the covariance matrix [Eq. (1)] can be written as The correlation between measured OCT signals at different spatiotemporal points, ρ G (Δr; r), can then be written as the product of the correlation coefficient of the true OCT signals, ρ S (Δr; r), and a factor expressing the noise effect, ρ SNR (Δr; r), as [29] ρ G (Δr; r) = ρ SNR (Δr; r)ρ S (Δr; r), The bias decreases with higher SNRs. As SNR → ∞, ρ SNR → 1 and ρ G → ρ S . Thus, the bias b r G decreases. However, the variance of estimation is not always small at high SNRs. The variance σ r G depends on ρ G , and rapidly increases as ρ G decreases [38,39]. Hence, the variance is large with small ρ S even when the SNR is high ( Table 1). The large variance σ r G will results in the noise of the correlation imaging.

SNR-corrected low-bias, low-noise complex correlation estimation
The correlation of signal ρ S (Δr; r) can be estimated using the noise power estimation technique presented in our previous studies [29,40]. In the present paper, we further modify this method to achieve a lower estimation noise than the previous approach.
The signal correlation was estimated as where s Spq are the entry of the estimated covariance matrix S S at p-th row and q-th column, and q N is an estimation of the noise power, obtained aŝ where t denotes time-averaging. g n (t) is the measured OCT signal obtained without a sample. This can be pre-determined by blocking the probe beam, or can be estimated from tissue measurement data using OCT signals from non-tissue regions. Here, we assume the noise N is identical in space. However, some OCT systems have a depth-dependent noise power. In this case, the noise powerq N should be treated as a function of depth z. The SNR-correction, however, does not reduce the variance of the estimation.
To reduce the estimation noise, a second averaging is applied to estimate the magnitude of the covariance and the geometric mean of variances. A correlation estimate with lower noise than Eq. (5) can be obtained using an averaging window w 2 as Similarly, a low-noise signal correlation estimate is obtained from Eq. (13) as The sign ε is applied to suppress the bias of the denominator of Eq. (16) after the averaging w 2 when the signal is zero (See Section 7). Note that the estimate of Eq. (16) can be negative or exceed 1. For imaging, the range of estimate r S will be wrapped as described in Section 5.1.

Correlation mapping optical coherence angiography
The present paper aims to demonstrate high-contrast OCA using the correlation coefficient of the true OCT signal, as shown in Section 2.2. We refer to this as correlation mapping OCA (cmOCA).
Typically, a single location on the sample is scanned multiple times, and the temporal correlation is calculated to contrast the blood flow. The temporal correlation between time points t and t + τ at space point (x, z) is denoted as ρ S (Δr; r)| Δr=(0,0,τ) in the notation of Section 2. Hereafter, this is denoted as ρ S (τ; r) for simplicity.
During the estimation of the temporal correlation coefficient, temporal averaging will be applied, and then the spatial distribution of the signal correlation coefficient is obtained, r S (r i , r j ) → r S (τ; x, z). This can be described as Two specific issues for in vivo cmOCA measurements are also discussed. These are associated with the bulk sample motion. The first issue is the bulk phase offset caused by the axial shift of the sample between the two time points for which the correlation is computed. Two approaches to avoid this bulk phase offset are described in Section 3.1. The second issue is the reduction in the correlation coefficient due to the general bulk motion of the sample. This creates artifacts, which affect our primary interest of visualizing only the region with real blood flow. Section 3.2 presents a method to remove such artifacts. The processing flow is described as a block diagram in Fig. 1.

Bulk-motion-induced phase-offset problem
The bulk motion of tissues induces a substantial phase shift, and this significantly decreases the correlation coefficient when the averaging window w 1 in Eq. (18) is extended along the temporal direction. There are two approaches to overcome this problem. One is to apply a bulkphase-offset correction to the estimation and the other is to form an estimation that is insensitive to the bulk-phase-offset.

Bulk-phase-offset correction
To eliminate the bulk-motion effect, we must detect the magnitude of the bulk-phase offset. For this purpose, previous studies have introduced sophisticated methods [9, 11, 41-43]. Although these methods exhibit high detection accuracy, they also suffer from high computational costs. Averaging after bulk-phase-offset Some simpler methods have been used for vasculature imaging [20, 44-46], but their use is limited to specific cases, such as when the non perfused region is dominant in each axial profile. Here, we describe a simple and robust bulk-phase-offset correction method. The key of this method is the estimation of locations of non-perfused tissue. The non perfused region is selected for each axial profile using the characteristics of data correlation coefficient ρ G . According to the relationship between correlation and SNR [Eq. (8)], locations with high correlation coefficient ρ G , are very likely to be tissue without perfusion and simultaneously to have high SNRs. Such high SNRs are preferable for the reliable estimation of the bulk-phase offset. The relationship is summarized in Table 2.
To implement this method, the covariance matrix estimation process [Eq. (3)] is resolved in two steps, i.e., bulk-phase-offset estimation and averaging after the phase correction. A block diagram of this processing framework is shown in Fig. 1(b). In the first step, a sample covariance matrix similar to Eq. (3) is estimated as: and the data correlation r G is then estimated using Eqs. (15) as To prevent the disruption caused by bulk-phase offset, the estimation w 1 of Eq. (19) should be computed by time-independent averaging, e.g., axial averaging. The bulk-phase offset between time points t and t + τ is then determined at depth z max (τ; x,t), which exhibits a high value of Using Eq. (21), the bulk-phase-offset is corrected and the covariance matrix is finally estimated as The signal correlation map r S (τ; x, z) is estimated by substituting Eqs. (19) and (21) to (22) into Eq. (16) as In such cases, another algorithm with phase-noise immunity is applied. The sample covariance matrix S G is estimated as similar to Eq. (3) as: To minimize the effect of the phase fluctuation, the averaging window w 1 is selected within time independent region as in Eq. (19). One example is axial averaging [27].

Removal of bulk-motion artifacts
Motion artifacts can occur in cmOCA images as a result of the de-correlation due to involuntary eye motion. This section describes a method to remove these artifacts from correlation coefficient estimates which is obtained without bulk-phase-offset problem (Section 3.1).
If the bulk tissue motion is rigid translation, and its velocity is quite smaller than the net blood velocity, the correlation ρ S (τ; r) can be decomposed as where ρ f (τ; r) is the correlation coefficient corresponding to blood flow. This is mainly caused by the motion of blood cells. ρ m (τ; x,t) is the correlation coefficient corresponding to bulk motion, which is uniform along the axial direction because single axial profiles are acquired simultaneously by Fourier-domain OCT. By assuming the bulk tissue motion velocity is constant during the temporal estimate window of correlation, the correlation estimate r S (τ; x, z) can be expressed as wherer f (τ; x, z) is the net correlation coefficient corresponding to blood flow. The reduction of r S (τ; x, z) due to r m (τ; x) is the source of the bulk motion artifacts. Because the correlation coefficient is estimated with bulk-phase-offset correction (Section 3.1.1) or phase-noiseimmune estimate (Section 3.1.2), decorrelation due to bulk-phase-offset is avoided. Our motion-artifact-correction method first estimates r m (τ; x), and then removes this quantity from r S (τ; x, z). As discussed in Section 3.1.1, tissue with a higher data correlation is more likely to be non perfused tissue with a high SNR. For such non perfused tissue, the correlation coefficient estimate r S (τ; x, z) will be identical to r m (τ; x), so r S (τ; x, z) in the non perfused tissue can be used to cancel out the bulk motion artifacts. The correlation coefficient r G (τ; x, z) is obtained from the covariance matrix estimate S G [Eq. (22) or (24)] in the manner of Eq. (15) as Hence, the signal correlation estimate r S at the depth of the maximum of r G (τ; x, z) is used as an estimate of r m :r m (τ; x) ≡ r S (τ; x, arg max z r G (τ; x, z)).
The bulk-motion-corrected signal correlation coefficient is finally obtained as This estimate represents the correlation due to flow, and is used to form bulk-motion-corrected cmOCA images.

Method
In this paper, we use Jones matrix OCT (JM-OCT), which is a multi-channel extension of OCT. JM-OCT was originally developed for polarization-sensitive imaging [31-33,50] and has recently been applied to multifunctional imaging [1,49,51]. First, we describe the scanning protocol for cmOCA measurements, and re-define variables to indicate sampling points (Section 4.1). The details of cmOCA implementations with JM-OCT is then described (Section 4.2).

Measurement protocol and correlation for cmOCA
For cmOCA measurements, multiple frames (B-scans) are obtained from the same position on the sample. The schematic diagram of these sampling points is shown in Fig. 2. A single

cmOCA implementation with JM-OCT
To obtain lower bias in correlation estimate, a large number of samples is required to estimate the covariance matrix [Eq. (3)]. To increase the size of the average w 1 , the bulk-phase-offset correction (Section 3.1.1) is preferable to extend the window w 1 along temporal axis. Each averaging operation in the pre-processing of the bulk-phase-offset estimate [Eqs. (19) and (20)] is replaced according to Table 3 as where D is the size of the averaging kernel along the axial direction. Here, we have used the K pixels with the highest correlation coefficients r M G in each axial profile to estimate the bulkphase offset. The axial indices of these pixels are denoted as z max k (τ; x, f ), with the index k identifying these pixels as k = 1, · · · , K. As in Eq. (21), the bulk-phase-offset estimate is where Med represents the median operation along the transverse (x-) direction, which is independently applied to the real and imaginary parts. In this paper, we consider K = 5, and select ±15 lines to compute the median Med . The covariance matrix estimate is obtained with lateral and frame averaging for w 1 . The correlation coefficient estimate is then described from Eq. (23) using the noise power of each channelq N (p) as where L is the size of the averaging kernel along the lateral direction.
The bulk-motion-artifact correction (Section 3.2) is then applied to r M S (τ; x, z). For this correction, the estimated data correlation r M G (τ; x, z) is first obtained by avoiding the noise power subtraction and the sign ε in the Eq. (33). The bulk-motion-induced net correlation is estimated using Eq.
, wherer M m is the estimated bulk-motion-induced correlation during the acquisition of F frames at a lateral point with index x. Finally, the bulk-motion-corrected signal correlation is obtained as To compare with and without the second averaging w 2 , the coherent composite signal was calculated and the same covariance estimation was applied. This mode is corresponding to the previous complex correlation estimation method [29] in standard single-channel OCT. The polarization channels are coherently averaged after the correction of phase offset [49] to obtain the coherent composite signal g C (x, z, f ) = 1 P ∑ P p g(x, z, f , p)e −Δθ (p;1) , where Δθ (p;1) is the phase offset between the channel p and channel 1. Because of the coherent averaging, the noise power will be reduced by the number of polarization channels P. Then, the correlation coefficient estimation was applied as the same as Section 4.2 except the polarization channel averaging.
The bulk-motion-artifact correction (Section 3.2) is then applied to r C S (τ; x, z). Without the noise power subtraction and the sign ε in the Eq. (35), r C G (τ; x, z) is obtained. The bulk-motioninduced net correlation is estimated using Eq. (28) asr C m (x) ≡ r C S (τ; x, arg max z r C G (τ; x, z)). The bulk-motion-corrected signal correlation with coherent composite is obtained as . (36)

Cross sections
Cross-sectional cmOCA images are displayed in grayscale with the estimated correlation coefficient r S . Low correlation corresponds to high perfusion. Hence, a correlation coefficient of 0 is assigned to white regions, and a value of 1 is assigned to black. The locations for which r −1 S < 1 are expressed as black.

En-face projection
The projection images of cmOCA are calculated as The dynamic range of the projection image is set to the [1: 99] percentile of each image's pixel value.

Retinal angiography with degree of polarization uniformity
Using JM-OCT for angiography not only increases the detection channels, but also allows the tissue discrimination ability to be used to enhance the angiographic imaging. Several retinal OCAs use a segmentation algorithm to separate some retinal and/or choroidal layers, thus enhancing each vasculature and/or increasing the power of diagnosis of eye diseases. Using DOPU measurement [52], it is suggested that the RPE can be distinguished. Here, we implement an advanced angiographic projection image construction using DOPU.
First, we generalize the encoding method of 3D multi-contrast information into a 2D projection. This is a kind of color volume ray casting, where each contrast is considered as a material that have reflection, absorption, and/or emission spectra of ray. Rays are introduced from top of the 3D image, passing through the voxels, and back-reflected. According to assigned properties of voxels, rays' colors changed. The composite en-face projection image E with multiple contrasts is calculated by averaging all reflected rays. This is defined as: where c j is the j-th contrast value. m a is an infinitesimally small color translation of rays, m s is the reflection spectrum of rays, and L 0 is initial color of rays. Δz is the axial distance per unit data point. Analogous to Lambert-Beer's law, c, m s , and m a are a pseudo-concentration, pseudo-spectrum of scattering cross section, and pseudo-spectrum of attenuation/emission cross section assigned for each contrast, respectively. The combination of absorption and emission can be considered as single translation in a color space. Note that e () denotes the matrix exponential operation.
In this paper, we utilize cmOCA and mDOPU [53] to create the color en-face image as: where cPV n (x) = ∑ n i=0 {1 − M[mDOPU(x, z i )]} is cumulative polarization variance from the surface to the depth z n . This means that color of rays changed according to cPV n by e 1 2 cPV n mΔz , rays are reflected at z n according to cmOCA r S f , color changed again by e 1 2 cPV n mΔz , and reflected rays from every depth z n are summarized. The projection image will discriminate a vessel by color whether it is located in front of or back of pigmented tissues, such as RPE.
Because the Euclidean distance between two colors in L * a * b * color space is roughly proportional to perceptual color difference between them, m is defined as the infinitesimally small color change in the ICC L * a * b * color space as To treat translation in color space, m is defined as 4-dimensional matrix while L * a * b * color space is 3 dimension. In this case, m is diagonalizable, m = Vdiag(λ 1 , λ 2 , λ 3 , λ 4 )V −1 and e cPV n (x)mΔz = Vdiag(e cPV n (x)λ 1 Δz , e cPV n (x)λ 2 Δz , e cPV n (x)λ 3 Δz , e cPV n (x)λ 4 Δz )V −1 .
T is the initial location in the color space when n=0, i.e., the assigned color at the top of cross-sectional image z = z 0 . In this paper, (L * 0 , a * 0 , b * 0 ) = (84.17, -13.95, 82.72). The transformation W [] is wrapping an input color location into the realizable L * a * b * space as: where Here, we set θ = 0.165 rad, a o = 40.3, and b o = 40.9. The color change according to cPV n is represented in L * a * b * color space as shown in Fig. 3. The transformation T [] translate color in L * a * b * color space to RGB color space by (1) converting the L * a * b * color space into the XYZ color space [54], (2) applying a chromatic adaptation (using Bradford transform) from the D50 to D65 white points [55], and (3) converting XYZ color space to the linear sRGB space [56]. The resulting E is scaled by fitting the range [5, 99.9] percentile of the average RGB value, (R + G + B)/3, to within the range of [0, 1]. As this is in the linear sRGB space, the gamma correction is applied to fit with the sRGB standard [56]. We choose the sRGB color space because it is commonly used in Windows® operating system. Note that we refer the "isolum" color map [57] to design m, L 0 , and W [].

Results
JM-OCT described in elsewhere [49] is used to scan human eyes. The sampling configuration of JM-OCT is the same as previous, with 512 axial lines per frame acquired and the scan repeated four times at the same location (F = 4). The set of frames is scanned at 256 locations to construct a 3D volume set. The scanning speed is 100,000 lines/s, and the time lag between consecutive frames is 6.4 ms. JM-OCT simultaneously measures four complex OCT signals that correspond to four polarization channels (P = 4). For every cmOCA calculation, spatial kernel sizes are fixed to be (L = 3, D = 3) along the x and z directions, respectively. For scattering OCT image, global-phase-corrected sensitivity enhanced imaging with 4 frames followed by the coherently composition of 4 polarization-channel signals are applied [49]. All protocols for measurement were approved by the Institutional Review Board of the University of Tsukuba. Written, informed consent was obtained prior to measurement.

Complex correlation with/without SNR correction
The SNR-independent vasculature images can be obtained using our signal correlation estimator. Figure 4 shows a comparison of SNR-independent and -dependent cmOCAs. A  Fig. 4(b)]. With the signal correlation r M S f , SNR-independent imaging is achieved [ Fig. 4(c)]. Even the SNR of OCT signal at some non-perfused tissue, such as nuclear layers and vitreous humor, is low, its corresponding correlation coefficient exhibits close to 1. The presented cmOCA with SNRcorrected correlation distinguishes perfused tissues from low-SNR regions without the use of an arbitrary masking threshold.

Comparison between single-channel and multi-channel implementations
For comparison, the correlation mapping method for the single-channel case [Eq. (29)] was applied to the complex summed OCT signal using the coherent composition [49]. The coherent composite signal clearly has a higher SNR than that of each channel, because the four polarization channels are coherently averaged. The correlation image is created with the coherent composite signal (Section 4.2.1).
Although the coherent composite signal has a higher SNR, the single-channel implementation exhibits significant spurious noise [ Fig. 5(a)]. The images shows the left eye of a age-related macular degeneration patient (male, 79 years old) with the scanning range of 4.5 × 4.5 mm 2 . When the SNR is low, the correlation between measured data ρ G is low and the correlation estimation exhibits high variance. In addition, the large bias correction for estimating the signal correlation ρ S enhances its variance.
In contrast, the multi-channel method exhibits high contrast and low noise vessel imaging owing to the averaging process [ Fig. 5 and cmOCA are assigned as blue and orange, respectively [57].
An abnormal blood flow signal can be observed in the cross section [ Fig. 9(d), white circle]. The DOPU image [53] reveals some RPE detachment and its discontinuity [ Fig. 9(e), white circle]. The corresponding signal appears in orange in the en-face composite projection image. This shows that the blood flow signal is not blocked by pigmented tissue, i.e., the RPE. The corresponding location exhibits hyperfluorescence in late-phase fluorescein angiography (FA) [ Fig. 9(c)]. The composite imaging may be a good non-invasive alternative to FA for detecting the collocation of abnormal vessels with RPE disruption, such as choroidal neovascularization. Note that no segmentation is applied to generate the composite image, and so no segmentation errors are present.

Discussion
We have presented a method for estimating the SNR-corrected complex correlation coefficient with low noise. To compare the statistical properties of estimators, a numerical simulation was applied. A pair of random numbers X i and X i+1 in consecutive sets (i-th and i + 1-th sets) was generated by using the method to generate correlated random numbers. Outcomes Z i and Y i of independent complex Gaussian variables with zero mean and standard deviation σ Y were generated. Then, Y i+1 was generated with the true correlation coefficient ρ as Y i+1 = ρY i + 1 − ρ 2 Z i . Y i and Y i+1 will be correlated with the correlation coefficient ρ. Another outcome N i of independent complex Gaussian variable with zero mean and standard deviation σ N was generated and added to provide X i = Y i + N i . By changing the ratio σ Y /σ N , arbitrary true SNR can be simulated. These generated signals X were processed with correlation estimators. For each trial, 4 sets of complex random Gaussian variable which have 9 sampling points were generated. The random number generation iterated four times to simulate 4 channels of JM-OCT. The data size for correlation estimate was 3 pairs from 4 sets of 9 data points times 4 channels, which is similar to the experiments. The mean and the standard deviation of estimates are obtained with 10,000 trials. The true correlation coefficient ρ was set to be 0.99, 0.9, and 0.8. The SNR 10 log 10 σ 2 Y /σ 2 N was set from 0 to 20 dB in 1 dB step. Figure 10 shows the comparison of correlation estimates. The previous method with coherent composite and the new estimators give lower bias estimates than the normal data correlation estimate for every true correlation coefficient [ Fig. 10(a), (b), and (c)]. At low SNR, the estimate with coherent composite exhibits lower bias than that of the new estimator. This is reasonable because the coherent composite decreases the noise. However, the standard deviation at high SNR is larger than that of the new estimator. This is because of the absence of the second averaging w 2 . The advantage of the new estimator becomes greater as signal correlation decreases. It means that the imaging noise of blood vessels (where signal correlation is decreased) is lower in the case of the new estimator. Hence, the new estimation scheme is probably a good compromise in the sense of both bias and noise. This fact suggests that although SNR will be decreased, applying second averaging w 2 by splitting the OCT signal, such as spectrally splitting method [23], may provide better vasculature imaging with standard single channel OCT system at the same averaging kernel size.
In our method, we introduced the sign function ε [Eq. (17)] in the denominator of the noisecorrected correlation estimate [Eq. (16)]. Although the absolute operation is applied inside the square root to avoid non-real-valued results in Eq. (16), the proposed estimator retains the sign ε after the noise correction. This modification from the previous method [29, 40] reduces the spurious appearance in the no signal region after second-step averaging w 2 .
In the SNR range used in the abovementioned simulation, the sign ε in Eq. (16) does not affect to the estimator's performance. The difference of the existence of ε becomes apparent at low SNR. The distributions of estimates with 0 SNR (-∞ dB) are shown in Fig. 11. The estimation is regarded as invalid if the estimated correlation is smaller than 0 or larger than 1. Because of the sign ε, the correlation estimates at 0 SNR spread and are less likely to be fall in the range of valid correlation coefficient estimate [0, 1]. With the sign ε, that probability of correlation estimate in the range [0, 1] occurs roughly 2.7 % [ Fig. 11(c)], while this is 32.0 % in the case without ε [ Fig. 11(a)] and 23.9 % in the case with thresholding [ Fig. 11(b)]. Without the second averaging w 2 , the thresholding method and the method using the sign ε provide the same probability. However, the second averaging generates a bias in the denominator of correlation coefficient estimate. This results in the large discrepancy after the second averaging.

Conclusion
Non-invasive vasculature imaging has been applied using complex-correlation-mapping OCA. Using low-bias and low-noise complex correlation estimation and motion artifact corrections, we have obtained better discrimination of perfusing tissue. The combination of perfusion contrast and polarization contrast demonstrates the potential of clinical utility in ophthalmology.