|Home | About | Journals | Submit | Contact Us | Français|
Decorrelation noise limits the ability of phase-resolved Doppler optical coherence tomography systems to detect smaller vessels exhibiting slower flow velocities, which limits the utility of the technique in many clinical and biological settings. An understanding of the statistical properties of decorrelation noise can aid in the optimal design of these systems and guide the development of noise mitigating strategies. In this work, the statistical properties of decorrelation noise are derived from the underlying statistics of the coherent imaging system and validated through comparison with empirical results and Monte Carlo modeling. Expressions for the noise distribution and the noise variance as a function of relevant imaging system parameters are given, and the implications of these results on both system and algorithm design are discussed.
Doppler optical coherence tomography (DOCT) techniques ,  allow microvascular imaging within the superficial 1–2 mm of highly scattering tissue. Applications of the technology extend from diagnostic medicine , where changes in vascular function are hallmarks of a broad range of diseases, to biological studies of vascular physiology in disease models –. While pilot studies of DOCT in these areas have generally confirmed its promise, they have fallen short of establishing clear clinical or biological utility. In large part, this is due to limitations in the flow sensitivity of DOCT imaging that prevent the reliable detection of small vessels exhibiting flow velocities below approximately 500 μm/s.
Flow contrast in phase-resolved DOCT (PR-DOCT) , which is an optical analog of Doppler ultrasound imaging, is derived from modulations in phase of light induced by its reflection from moving scatterers within the bloodstream. These motion-induced phase shifts are revealed by differencing two measurements of the same location acquired at two time points. The accuracy of these measurements, specifically the phase noise level of these measurements, determines the sensitivity of the system. Initial work on PR-DOCT imaging utilized time-domain OCT techniques that require a mechanical rapid scanning optical delay which due to poor repeatability introduced large errors. Through the adoption of Fourier-domain techniques which eliminate the need for scanning optical delays ,  these instrumentation-based noises have been eliminated , . Within this regime of negligible instrumentation noise, if we assume that the two measurements to be differenced can be made to sense exactly the same physical structure of the tissue, then any detected phase shift between them in excess of that predicted by the signal-to-noise ratio (SNR) can be attributed to motion of scattering sites within the sample. However, a number of mechanisms prevent repeated measurement of identical areas in the tissue. The resulting displacements between these measurements decorrelate the detected phase signal, inducing a phase difference, termed decorrelation noise, that degrades the flow sensitivity of the system. The significance of this noise to the overall sensitivity of PR-DOCT systems has been demonstrated empirically in a number of recent works , , and mitigating approaches have been proposed . However, despite its prevalence, a description of the statistical properties of decorrelation noise has not been reported. Knowledge of these properties can help to guide the optimal design of imaging systems and can influence the design of algorithms deriving blood flow and associated vascular parameters from the phase difference data. The aim of this work is first to derive these properties, and second to briefly illustrate their implications in two aspects of the design of a PR-DOCT imaging system.
The organization of this work is as follows: Section II describes the origin of decorrelation noise in PR-DOCT; Section III presents the derivation of its statistical properties; Section IV validates these derivations through comparison with both empirical results and Monte Carlo modeling; and finally Section V discusses the noise properties and its implications on imaging performance.
The operating principles of PR-DOCT have been described in detail in a number of prior works , , , . The principles relevant to this work are briefly summarized in this section, and the reader is referred to the aforementioned references for greater detail. OCT techniques allow the measurement of the complex reflectivity of a scattering sample as a function of depth (z) coordinate. Scanning the imaging beam along a transverse (x) dimension enables two-dimensional (x,z) mapping of the complex reflectivity. Conventional OCT images are generated from the magnitude of the scattered signal and allow imaging of tissue microarchitecture based on endogenous variation in scattering cross-sections. In PR-DOCT, a second image highlighting the vasculature is derived from the phase of the scattered signal. To generate this image, a series of measurements are acquired, and the phase difference between these measurements is displayed . The phase difference between two complex reflectivity measurements is therefore the fundamental signal parameter in PR-DOCT imaging.
Decorrelation noise results from the comparison of two reflectivity measurements acquired at slightly displaced transverse locations. Constant velocity beam scanning along the x-dimension, for example, induces a deterministic displacement between phase measurements, and the resulting decorrelation noise has been termed beam-scanning noise , . Meaningful displacements can additionally result from non-deterministic sources such as an instability of the beam scanning mechanism. We have observed that, even with an unscanned beam, the dithering of a stationary but actively controlled galvanometric mirror induces decorrelation noise at limiting levels. In this work, we adopt the more general term decorrelation noise to encompass all noise resulting from a comparison of measurements acquired at displaced locations regardless of the specific mechanism inducing that displacement. Through a number of empirical studies, it has been shown that the magnitude of decorrelation noise exceeds all other noise sources over a range of operating parameters that encompasses typical system configurations , .
Fig. 1 defines generally the parameters of a pair of beams yielding a single phase-difference measurement. Although we assume Gaussian beam profiles in this work, the methodology supports arbitrary beam profiles. The two beams are labeled beam A and beam B and result in complex reflectivity measurements at a given depth denoted by the phasors ΓA and ΓB, respectively. As a result of the coherent imaging properties of OCT, these measurements are known to follow a complex Gaussian distribution , . The phase difference between them, θ, is the parameter of interest and is illustrated in Fig. 2a. The aim of this work is to derive the statistical distribution of this phase difference as a function of the fractional displacement given by the ratio of the beam displacement, δx, to the beam (1/e) radius, ωo.
A few notes are appropriate before proceeding. First, in this work we are concerned with measurements made within tissue and thus we assume that the reflected signals originate from a large number of randomly positioned scatterers within a particular voxel. The analysis is therefore not applicable to the imaging of highly organized samples with large-scale (relative to the imaging resolution) spatial correlations in complex reflectivity, e.g., a mirror. Second, because of phase wrapping ambiguities in the PR-DOCT measurement process, the phase difference is a circular parameter and is meaningfully defined only over a range of 2π necessitating the use of circular statistics . Third, the statistical distributions of the phase difference at limiting values of fractional displacement are known a priori: as (δx/ωo) becomes increasingly large, the measurements ΓA and ΓB become fully uncorrelated and the distribution of the phase difference should asymptote to a uniform distribution from −π to π; as (δx/ωo) approaches zero, the measurements ΓA and ΓB become identical and the distribution should asymptote to a delta function centered at zero.
We begin the derivation by reformulating the problem definition illustrated in Fig. 2a to a form that is free of statistically dependent parameters. The beam B is decomposed into a linear combination of beam A and a third beam C with profile EC(x,y) that is orthogonal to beam A,
Because the OCT system is linear, the measured complex reflectivity ΓB resulting from beam B can be likewise decomposed into two terms, each being associated with one of the terms of the modal decomposition in Eq. 1,
The phasor ΓC is defined as the measured reflectivity associated with beam C. The squaring of the coefficients cBA and cBC in this expression results from the modal coupling that occurs both on illumination (i.e., the fraction of the illumination beam B that mode-matches to beam A) and detection (i.e., the fraction of the beam A that mode-matches back into the illumination/collection fiber that defines beam B). This is in agreement with a prior model of OCT detection .
Expressions relating the coefficients cBA and cBC to the fractional displacement (δx/ωo) can be found using mode-overlap integrals. Using the definition of the field profiles of beams A and B implied by Fig. 1, the coefficient cBA is given as
To find a similar expression for the coefficient cBC, we first note that the ensemble averaged squared magnitudes of ΓA and ΓB are necessarily equivalent. Further, we are free to define arbitrarily the relative scaling of the beam C and the coefficient cBC, i.e., we can reduce the magnitude of EC(x,y) by ½ if the coefficient cBC is increased by a factor of 2. According to our formalism, the measured reflectivity ΓC is associated with the beam C defined by EC(x,y), and therefore scales with this beam profile. Thus, we must first define the scaling of EC(x,y) to obtain a specific expression for cBC in Eq. 2. For simplicity, we assume that EC(x,y) is defined such that all three measurements ΓA, ΓB, and ΓC have equivalent average squared magnitudes. Under this condition, Eq. 2 implies that the coefficients satisfy , where we have used the statistical independence of ΓA and ΓC that results from the orthogonality of the profiles of their associated beams. In combination with Eq. 3, this expression defines the coefficient cBC, and allows the measurement ΓB to be expressed as
Next we introduce a rotated coordinate system (X′,Y′) such that the X′-axis is defined by the angle of the phasor ΓA. The phase difference, θ, which is invariant upon such a rotational transformation, is given by the phase angle of the phasor ΓB in this rotated coordinate system, as illustrated in Fig. 2c. We can express the phasor ΓB in the rotated coordinate system as
Solving for the distribution of the angle of this phasor from the known underlying distributions of |ΓA| and ΓC reveals the distribution of the phase difference, θ.
As described above, both ΓA and ΓC follow circular Gaussian statistics defined by the two-dimensional probability distribution function
where the x and y coordinates represent the real and imaginary components of Γ respectively. The average squared magnitude of Γ is equal to the parameter . Noting that the average squared magnitudes of ΓA and ΓC are equal to one by definition, the distribution of the second term of Eq. 5 is therefore described by Eq. 6 with the parameter . The coordinate rotation used to arrive at the expression of Eq. 5 introduces a functional dependence on the magnitude of the measurement ΓA. The distribution of the magnitude of a circular Gaussian variable described by Eq. 6 is given by a Rayleigh distribution ,  defined by the distribution function
The distribution of the first term of Eq. 5 is therefore described by Eq. 7 with the parameter μ2 = α2. The distribution of the phasor ΓB in the rotated coordinate frame can be derived from these distributions of its constituent terms. Specifically, the distribution of ΓB is given by the convolution of the distributions of its two terms ,
where the convolution integral is evaluated only for x′>0 as a result of the one-sided nature of the Rayleigh distribution (Eq. 7).
The following substitutions x = rcos(θ) and y = rsin(θ) recast this distribution into polar coordinates,
Finally, integration over the radial dimension leaves the distribution as a function of only the angle, θ,
Equation 10 defines the decorrelation noise distribution over the range from −π to π as a function of the coefficient α defined by the fractional displacement, . Note that in the limit of an infinite fractional displacement, the coefficient α approaches zero and the distribution converges to a uniform distribution. Likewise, in the limit of no fractional displacement, the coefficient α approaches one, and it can be shown that Eq. 10 approaches a delta function centered at zero. Fig. 3 plots this distribution for three values of the fractional displacement. It is immediately apparent from Fig. 3 that the distribution possess significant power in its tails. The excess kurtosis ,  of the distribution is a measure of this property. Values near zero indicate similarity with normal distributions (or their circular statistical analogs such as the wrapped normal or Von Mises  distributions) wherein relatively little power exists in the distribution tails. Values greater than zero indicate increasing power in the tails with corresponding narrowly peaked distributions about the center. The excess kurtosis of the decorrelation noise distribution for fractional displacements below 0.1 exceeds ~35, indicating that the tails play a more significant role (relative to Gaussian noise) in determining the overall variance of the distribution. This property can be used to inform the choice of data processing and filtering algorithms (see Section V).
Validation of the derived decorrelation noise distribution presented in Eq. 10 is achieved through comparison with both measured noise statistics and numerically simulated noise statistics. These validating methods are described here.
Noise levels were acquired using a phase-resolved optical frequency domain imaging (OFDI) system. The system design and its ability to make phase-resolved measurements has been reported previously , . The OFDI system utilized a wavelength-swept narrowband laser source with a 10 kHz repetition rate. The laser swept from 1215 to 1360 nm with a linewidth of ~0.17 nm. The system featured both balanced and polarization diverse detection. A 1% tap coupler in the sample arm was used to generate a continuous phase recalibration signal that allowed phase variations due to timing jitter to be removed, resulting in SNR-limited phase measurements absent decorrelation noise  Beam scanning was achieved with a high-speed galvanometer (Cambridge Technology 6220H, 5 mm mirror) driven by the analog output of the acquisition board. To prevent digitization errors from affecting beam scanning, the analog voltage waveform used to drive the galvanometer was generated over the board’s full range (12 bits) and subsequently attenuated to the desired amplitude by a variable voltage divider. A focused beam diameter (intensity 1/e2) of 24 μm was used. Doppler processing of the images followed established techniques , , . Bulk motion artifacts were removed by subtracting from each depth line in the Doppler image its measured median phase shift . The phase differences at each depth from each polarization channel were recorded. Noise measurements were obtained by imaging a homogeneous rubber scattering phantom. Imaging was performed over a transverse field of view (FOV) of 3.8 mm, and phase differences were calculated between adjacent Alines. The number of A-lines per frame was varied to modify the fractional displacement. Noise was measured in a region of interest extending from 170 μm to 500 μm in depth and encompassing the center 85% percent of the image to exclude artifacts resulting from either galvanometer flyback or from specular reflections at the phantom surface. Measurements from a single polarization channel were used. Fig. 4(a–b) presents reflectivity and phase difference images over a subset of the ROI over which statistics were measured.
A Monte Carlo model of decorrelation noise was developed as a second validating method. This model follows the problem definition of Fig. 1 and prior models of OCT detection . The scattering of the sample as a function of x-coordinate is represented by an array of random complex variables, Ri, with real and imaginary components uniformly distributed over [−1,1]. The beam (field) profile array, Ei, is Gaussian in shape, and generated over n points falling to 1/e5 at the edges. The complex reflectivities, ΓA and ΓB, were calculated according to
where Δ is an integer that defines the shift of beam B relative to beam A, and the array R is generated with (n+Δ) elements. Fig. 4c illustrates this procedure. Each randomly generated array R yields a pair of complex reflectivities ΓA and ΓB and a single phase difference, θ = angle(ΓAΓB*). Repeating this process for N arrays generates a distribution of phase differences. Because of the high Kurtosis in the form of the decorrelation noise distribution, accurate measures of the distribution required large values of N. In this work, a value of 1 × 106 was used for both parameters n and N.
To validate the distribution of Eq. 10, decorrelation noise was measured and simulated as described in sections IV.A and IV.B for fractional displacements of 0.078 and 0.156. Histograms of the resulting noise were normalized by the total sample count and histogram bin width to express the probability density in units of (1/radians) matching that given by Eq. 10. Fig. 5 plots the analytic, measured, and Monte Carlo noise distributions for these fractional displacements showing excellent agreement between all methods.
As discussed earlier, knowledge of the specific properties of a dominate noise source can guide the optimal design of an imaging system and influence the selection of signal processing and noise filtering algorithms. In this section, we present two such examples along these lines. First, we use these statistics to guide the selection of the fractional displacement, which generally involves a compromise between the need for (i) rapid imaging enabled by larger fractional displacements and (ii) high sensitivity that requires smaller fractional displacements. Using Eq. 10, we can quantify the magnitude of decorrelation noise as a function of the fractional displacement,
This integral can be evaluated straightforwardly using numerical integration. We note that the noise magnitude predicted by Eq. 12 conflicts with a previously reported expression  but agrees both with measured and simulated noise magnitudes. To illustrate the utility of Eq. 12 in the design of an imaging system, we follow the method of  which examines the relative magnitude of decorrelation noise to that resulting from the SNR limit. In Fig. 6, the fractional displacement (lower y-axis) is plotted as a function of its resultant noise magnitude. Note that, for large displacements, the noise magnitude asymptotes to radians - the standard deviation of a uniform distribution over the range [−π,π). Below fractional displacements of approximately 10%, the noise follows an approximately linear dependence on displacement. Also plotted in Fig. 6 is the SNR (upper y-axis) required to achieve a given noise magnitude in the absence of decorrelation noise. This relationship ,  is approximated by Δθ2 ≈ SNR−1. Although this expression is inexact because it neglects the circular nature of the phase difference, simulations indicate the error in its predicted RMS noise magnitude is less than 3% for SNRs exceeding 10 dB. If we assume an SNR of light scattered from vessels of interest ranging from 15 dB to 30 dB, typical for vessels located at a depth of ~500 μm, one can see from Fig. 6 that fractional displacements less than 0.064 and 0.009 respectively are required to maintain nearly SNR-limited performance. These fractional displacements along with the focused beam size and system A-line rate determine maximum imaging speeds such that decorrelation noise remains at or below SNR-limited performance levels.
As a second example, we use the specific distribution of decorrelation noise described by Eq. 10 to quantify the performance of three simple filtering approaches commonly applied in Doppler imaging to reduce noise. An examination of the distribution presented in Fig. 3 highlights the high kurtosis of this noise, suggesting that relatively few samples of large deviation from the mean disproportionately affect the overall noise performance, and that filtering approaches with outlier rejection properties are well matched to decorrelation noise. We analyzed the performance of straightforward mean and median filters (applied directly to the phase difference data), and an amplitude weighted mean filter that weights each phase difference by a power of the amplitude of one of the component signals, |ΓA|k. Using the numerical model, filters were applied to decorrelation noise samples of size n = [3, 10, 20, 30, 40], and the residual error of each was recorded for 250,000 independent sets for each set size. Fig. 6 plots the standard deviation of this error for each of the three filters as a function of the sample size. In Fig. 6, the amplitude weighted filter used the squared magnitude (k = 2) as the weighting parameter. The dependence of the estimator error as a function of the weighting parameter k for set size n = 20 is shown in the inset and indicates that the error is minimized near this value of k = 2. As anticipated, the mean filter performs poorly relative to the median filter due to the poor outlier rejection properties of the former. Further, the amplitude-weighted filter outperforms both the mean and median filter. The origin of outlier rejection for the filter in this application results from its ability to de-emphasize those phase differences for which the magnitude of the first measurement, |ΓA|, is small relative to its ensemble average. A small magnitude |ΓA| increases the susceptibility of the measurement to phase perturbation by the phasor |ΓC| (see Fig. 2c), and therefore lower values of |ΓA| are characteristic of the outlier points (data not shown). We note that, although we have used the magnitude of the phasor ΓA as the weighting, the complimentary phasor ΓB could be used with the same results by symmetry, or the average of the two phasor magnitudes could be used. Further, it is interesting to note that amplitude weighting by the squared magnitude (k = 2) is additionally optimal if the system is SNR limited assuming Gaussian noise statistics, and the weighted mean filter can thus be applied efficiently across both of these noise limiting regimes.
The examples above serve both to describe specifically the implications of the derived statistical properties of decorrelation noise on the design of PR-DOCT systems as well as to illustrate generally the practical utility of these properties in the design of phase-sensitive OCT imaging systems.
The statistical distribution of decorrelation noise in PR-DOCT imaging was derived from the underlying statistical properties of OCT signals and using optical mode overlaps to quantify correlations. Expressions for the noise distribution and variance as a function of the fractional displacement between beams are given and validated through comparison with both measured and numerically modeled noise distributions. The utility of these noise properties in the design of imaging systems was illustrated in two examples. The improved performance that results can directly enhance the ability of PR-DOCT systems to visualize smaller vessels or to image more rapidly, both of which can lead to enhancements in the clinical or biological utility of the technology.
This work was supported in part by the U.S. National Institutes of Health under grants R21 CA125560 and K25 CA127465.