Research Papers: Imaging

Statistical analysis of motion contrast in optical coherence tomography angiography

[+] Author Affiliations
Yuxuan Cheng, Li Guo, Cong Pan, Tongtong Lu, Zhihua Ding, Peng Li

Zhejiang University, College of Optical Science and Engineering, State Key Laboratory of Modern Optical Instrumentation and the Collaborative Innovation Center for Brain Science, Hangzhou, Zhejiang 310027, China

Tianyu Hong

Zhejiang University, College of Biomedical Engineering and Instrument Science, Hangzhou, Zhejiang 310027, China

J. Biomed. Opt. 20(11), 116004 (Nov 02, 2015). doi:10.1117/1.JBO.20.11.116004
History: Received June 22, 2015; Accepted October 2, 2015
Text Size: A A A

Open Access Open Access

Abstract.  Optical coherence tomography angiography (Angio-OCT), mainly based on the temporal dynamics of OCT scattering signals, has found a range of potential applications in clinical and scientific research. Based on the model of random phasor sums, temporal statistics of the complex-valued OCT signals are mathematically described. Statistical distributions of the amplitude differential and complex differential Angio-OCT signals are derived. The theories are validated through the flow phantom and live animal experiments. Using the model developed, the origin of the motion contrast in Angio-OCT is mathematically explained, and the implications in the improvement of motion contrast are further discussed, including threshold determination and its residual classification error, averaging method, and scanning protocol. The proposed mathematical model of Angio-OCT signals can aid in the optimal design of the system and associated algorithms.

Figures in this Article

Optical coherence tomography angiography (Angio-OCT) is capable of contrasting the dynamic blood flow against the static tissue bed with high spatial resolution and high motion sensitivity (down to capillary level) in a depth-resolved manner. Up until now, several Angio-OCT algorithms have been developed for generating the motion contrast.117 Generally, each spatial position should be sampled/imaged several times with a certain time interval by using repeated10,13,16 or dense14,18 scanning protocols in the OCT system. Then the temporal changes in amplitude (or intensity),2,4,10,11 phase,5,9,19 or complex-value7,13,17,20 of OCT signals over such a time interval are analyzed with different processing algorithms, such as speckle variance,24 Doppler variance,7,8,14 phase variance,5,6 differential calculation,18 and correlation mapping.11,21 By circumventing the exogenous contrast injection, such motion-contrast Angio-OCT provides great advantages over conventional flourescence-based angiography.

Knowledge of the statistical properties of Angio-OCT signals would be helpful for further understanding the origin of the motion-contrast and guiding the optimization of the system and associated algorithms. It is well known that the motion-contrast Angio-OCT is mainly based on the temporal dynamics of OCT scattering signals, and the algorithms of differential calculation are widely used for dynamics analysis, including the amplitude differential (AD) and complex differential (CD) algorithms. The temporal statistics of the OCT amplitude signals have been well documented in the literature.22,23 In a similar way, the temporal statistics of the complex-valued OCT signals can be mathematically described based on the knowledge of statistical optics.24 In this study, the statistical properties of AD- and CD-Angio-OCT signals are derived in theory, and the implications of the developed statistical model are briefly illustrated.

This work is organized as follows. (1) The temporal statistics of the complex-valued OCT signals were first described and the mathematical statistics of AD- and CD-Angio-OCT were further deduced in Sec. 3. (2) The statistical properties derived in theory were validated through the flow phantom and live animal experiments in Sec. 4. (3) The potential implications of the statistical properties were discussed in Sec. 5.

Flow Phantom and Animal Preparation

The flow phantom was made of an agarose gel mixed with 5% milk to mimic the static scattering tissue background. A capillary tube with an inner diameter of 0.5 mm was embedded in this tissue-like phantom. A 3% milk solution was pumped into the tube at a constant rate with a syringe pump (KDS 100 series, Stoelting Co., Wood Dale, Illinois) to simulate the flowing blood.

C57BL/6 mice of 8 to 10 weeks old were used in animal experiments. A mouse was anesthetized by intraperitoneal injection of 10% chloral hydrate (4ml/kg). Its head was fixed in a stereotaxic frame (Stoelting), and the scalp was retracted. The skull was thinned by using a saline-cooled dental drill to generate a window of 3mm×3mm area and to facilitate the optical penetration within the cortex at an 850 nm wavelength. All animals were provided by the Experimental Animal Center and treated with the guidelines of the Institutional Animal Care and Use Committee of Zhejiang University.

System Setup and Scanning Protocol

The imaging system in this study was built based on a typical configuration of spectral domain OCT (SDOCT). Briefly, the light source is a broadband superluminescent diode with a central wavelength of 850 nm and a full width at half maximum bandwidth of 100 nm, theoretically offering a high axial resolution of 3.2μm in air. The measured lateral resolution is 15μm. A high-speed spectrometer equipped with a fast line scan CMOS camera was used as the detection unit in our system, providing a 120-kHz line scan rate.

In this study, MB-mode scanning protocol was used for dynamic analysis. Each B-scan was formed by 512 A-lines, determining a rate of 190 fps. Therefore, 1000 repeated B-scans were sequentially acquired in the same cross-section within 5.3s, generating a three-dimensional (3-D) OCT data cube (z,x,n), where n is the B-frame index, equivalent to the time dimension. z and x represent the depth and transverse position, respectively, as shown in Fig. 1(a).

Graphic Jump Location
Fig. 1
F1 :

Flow chart of the statistical analysis in optical coherence tomography angiography (Angio-OCT). Amplitude differential (AD) method is used as an example for illustration. (a) Three-dimensional (3-D) OCT data cube (z,x,n). (b) Temporal distribution of the complex-valued OCT signals. (c) Statistical distribution of the AD-Angio-OCT signals.

Processing Algorithm

The depth-resolved complex reflectivity of a scattering sample is reconstructed by performing Fourier transform of the spectral interference fringe signals in SDOCT. The complex-valued OCT signals of the n’th repeated B-frame is denoted as A˜(z,x,n). A map of the amplitude signals A(z,x,n) is used to generate the structural image. The differences of the amplitude A(z,x,n)18 and complex-valued A˜(z,x,n)17 OCT signals between adjacent B-frames are computed for AD- and CD-Angio-OCT, respectively, as follows: Display Formula

AngioOCTAD=aAD(z,x,n)=A(z,x,n+1)A(z,x,n),(1)
Display Formula
AngioOCTCD=aCD(z,x,n)=|A˜(z,x,n+1)A˜(z,x,n)|,(2)
where aAD and aCD represent the amplitude of AD-Angio-OCT and CD-Angio-OCT signals, respectively. Typically, the absolute value of aAD is used in the final angiograms in AD-Angio-OCT. The adjacent B-frames are acquired at the same cross section with a certain time interval (t). Then thresholds are used to identify the dynamic areas. Due to the bulk motion, prior to the subtraction operation in Eq. (2), the global phase fluctuations are determined and compensated by a histogram-based phase selecting process.9,16,2527

In this study, the measured data were statistically analyzed with the histogram, and then compared with the statistical model proposed in theory. R-square (R2) statistic was measured to evaluate how well the experimental outcome fit the theoretical prediction.28

The flow chart of the mathematical analysis is depicted in Fig. 1. Because the complex-valued OCT signal is used in the CD-Angio-OCT, the temporal distribution of the complex-valued OCT signals was briefly described [Fig. 1(b)]. Then the mathematical statistics of Angio-OCT signals were deduced in theory [Fig. 1(c)]. The following derivations are mainly based on the knowledge of random phasor sums and transformations of random variables which have been well elaborated in Ref. 24. We focus on the physical plausibility and the associated assumptions used for solving OCT problems.

Temporal Statistics in Optical Coherence Tomography

For brevity, at a given position, the complex-valued OCT signals A˜(z,x) are denoted by a phasor ad/sexp(jθd/s)Display Formula

A˜(z,x)=ad/sexp(jθd/s),(3)
where the subscripts d and s represent the dynamic and static signals, respectively. As a result of the coherence gating used in OCT, the phasor ad/sexp(jθd/s) is a complex addition of many small phasors, arising from a collection of small scatters that are distributed within the OCT voxel of interest.29

In the dynamic regions, the small scatters are mainly contributed by the moving red blood cells (RBCs). The corresponding light scattering signals are time variant. Thus, the resultant phasor adexp(jθd) can be regarded as a complex sum of a large number of small, random phasors, usually referred to as random phasor sum. Display Formula

adexp(jθd)=l=1Mβlexp(jϕl)=rd+j·id,(4)
where β, ϕ, and M are the amplitude, phase, and total number of the small independent phasor in the dynamic regions, respectively. rd and id are the real and imaginary parts of the resultant phasor, respectively. The properties of random phasor sum have been detailed in Ref. 24. Briefly, in the limit of a very large M, the resultant phasor adexp(jθd) is a circular complex Gaussian random variable.24 And the joint probability density function (PDF) of rd and id is Display Formula
fRdId(rd,id)=12πσd2exp(rd2+id22σd2),(5)
where σd2 represents the variance of the real and imaginary parts. Using the transformation rd=adcosθd, id=adsinθd, the temporal PDF of amplitude ad is derived.24Display Formula
fAd(ad)={adσd2exp(ad22σd2)ad00otherwise.(6)
Therefore, in the dynamic regions, the amplitude of the OCT signals obeys a Rayleigh distribution with mean π/2σd and variance (2π/2)σd2, which is in good agreement with the literature.30

In the static regions, the light scattering signals are time invariant, and can be regarded as a strong, constant phasor. In this case, the system noise becomes the primary random contribution, and the phasor anoiseexp(jθnoise) of the noise can be regarded as a weak random phasor sum Display Formula

anoiseexp(jθnoise)=l=1Nαlexp(jφl)=rnoise+j·inoise,(7)
where α, ϕ, and N are the amplitude, phase, and total number of the small independent phasor in the static regions, respectively. rnoise and inoise are the real and imaginary parts of the random phasor sum, respectively. The weak random phasor sum is a circular complex Gaussian random variable with zero mean and standard deviation σs. In most situations of interest, the OCT signal C is much stronger than the system noise. Thus, the resultant phasor asexp(jθs) equals a strong constant phasor C plus a weak random phasor sum, as follows: Display Formula
asexp(jθs)=C+anoiseexp(jθnoise)=(C+rnoise)+j·inoise=rs+j·is,(8)
where rs and is are the real and imaginary parts of the resultant phasor, respectively. The joint PDF of rs and is follows a two-dimensional Gaussian distribution24Display Formula
fRsIs(rs,is)=12πσs2exp[(rsC)2+is22σs2],(9)
with Cσs, the approximation is made that Display Formula
asC+rnoise.(10)
Variations in the amplitude as are primarily caused by the real part rnoise of the weak random phasor sum, thus we have Display Formula
fAs(as){12πσsexp[(asC)22σs2]as00otherwise.(11)
In the static regions, the amplitude of the OCT signals obeys a Gaussian distribution with mean C and variance σs2, which is also in agreement with the literature.30

Angio-Optical Coherence Tomography Statistics
Amplitude differential-angio-optical coherence tomography

In AD-Angio-OCT, substituting Eq. (3) into Eq. (1) yields Display Formula

AngioOCTAD=aADd/s(n)=ad/s(n+1)ad/s(n).(12)
In the dynamic regions, we have Display Formula
aADd(n)=ad(n+1)ad(n).(13)
Due to the moving of RBC, the variables ad(n+1) and ad(n) which represent the amplitude of the scattering light can be regarded to be independent and random, and they follow the same but independent Rayleigh distribution, as described by Eq. (6). The statistics of the random variable aADd can be derived as follows:24Display Formula
fAADd(aADd)=+fAd(waADd)fAd(w)dw=π4σd(1aADd22σd2)[1erf(aADd2σd)]exp(aADd24σd2)+aADd4σd2exp(aADd22σd2),(14)
where the w is an intermediate variable. The erf is the Gauss error function. Unfortunately, we do not have a simplified analytic expression of Eq. (14). According to numerical simulation in MATLAB® as shown in Fig. 2, Eq. (14) is extremely close to a Gaussian distribution with zero mean and variance σd2: Display Formula
fAADd(aADd)12πσdexp(aADd22σd2).(15)
And the absolute value |aADd| follows a truncated Gaussian distribution with mean σd/2π and variance (12/π)σd2: Display Formula
f|AADd|(|aADd|){22πσdexp(|aADd|22σd2)aADd00aADd<0.(16)
In the static regions, the OCT scattering signal C is time invariant, and remains constant in all the B-frames for a given spatial point. Substituting Eq. (10) into Eq. (12), we have Display Formula
aADs(n)=as(n+1)as(n)rnoise(n+1)rnoise(n).(17)
The rnoise(n+1) and rnoise(n) are random variables, which obey the same but independent Gaussian distribution [Eq. (9)]: Display Formula
fRnoise(rnoise)=12πσsexp(rnoise22σs2).(18)
The statistics of the random variable aADs can be deduced by a transformation of random variables: Display Formula
fAADs(aADs)=+fRnoise(vaADs)fRnoise(v)dv=12π(2σs)exp[aADs22(2σs)2],(19)
where v is an intermediate variable. It is a Gaussian distribution with zero mean and variance 2σs2. Similar to Eq. (16), the absolute value |aADs| obeys a truncated Gaussian distribution with mean 2σs/π and variance (24/π)σs2: Display Formula
f|AADs|(|aADs|)={22π(2σs)exp[|aADs|22(2σs)2]aADs00aADs<0.(20)

Graphic Jump Location
Fig. 2
F2 :

Numerical simulation of AD-Angio-OCT signals in the dynamic region. The solid lines are the numerical simulation of Eq. (14) with different variances. The dash lines are the corresponding Gaussian distributions of Eq. (15). The R-square value (R2) is 0.99.

Complex differential-angio-optical coherence tomography

In the CD-Angio-OCT, Eq. (2) can be rewritten as Display Formula

AngioOCTCD=aCDd/s(n)=|ad/s(n+1)exp[jθd/s(n+1)]ad/s(n)exp[jθd/s(n)]|=[rd/s(n+1)rd/s(n)]2+[id/s(n+1)id/s(n)]2.(21)
Here, we define four intermediate random variables ξd/s and ηd/s: Display Formula
ξd/s=rd/s(n+1)rd/s(n)={rd(n+1)rd(n)dynamic regionrnoise(n+1)rnoise(n)static region,(22)
and Display Formula
ηd/s=id/s(n+1)id/s(n)={id(n+1)id(n)dynamic regioninoise(n+1)inoise(n)static region.(23)
Referring to Eq. (5), rd(n+1), rd(n), id(n+1), and id(n) follow the same but independent Gaussian distribution. According to Eqs. (17)–(19), the subtraction of two independent Gaussian distributions remains Gaussian with a modified variance. Thus, the random variables ξd and ηd have the same but independent Gaussian distribution with zero mean and variance 2σd2. Referring to the derivation from Eqs. (5) to (6), in the dynamic regions, the amplitude aCDd obeys a Rayleigh distribution with mean πσd and variance (4π)σd2, as follows: Display Formula
fACDd(aCDd)={aCDd2σd2exp(aCDd24σd2)aCDd00aCDd<0.(24)
Similarly, in the static regions, the statistics of the variable aCDs obey a Rayleigh distribution with mean πσs and variance (4π)σs2, as follows: Display Formula
fACDs(aCDs)={aCDs2σs2exp(aCDs24σs2)aCDs00aCDs<0.(25)
In this section, the mathematical statistics of the AD-Angio-OCT and CD-Angio-OCT were further deduced. According to Eqs. (16), (20), (24), and (25), the Angio-OCT statistics depend on the variances of the OCT statistics, i.e., σd2 and σs2. As a summary, the OCT and Angio-OCT statistics of the dynamic and static signals were tabulated, as shown in Table 1.

Table Grahic Jump Location
Table 1Statistics in optical coherence tomography and Angio-OCT.
Flow Phantom Imaging

Figure 3(a) is a representative OCT structural cross section of the flow phantom. The transparent tube is clearly visualized. The regions inside and outside the tube correspond to the dynamic fluid and the static solid gel, respectively. Figures 3(b) and 3(c) show the corresponding cross-sectional angiograms of AD- and CD-Angio-OCT, respectively. As indicated in Fig. 3(a), in the recorded cross section, two spatial points were randomly selected from the regions of static solid-gel (zs,xs) and dynamic fluid (zd,xd), respectively, and used for the statistical analysis of the OCT amplitude signals and the Angio-OCT signals.

Graphic Jump Location
Fig. 3
F3 :

(a) Representative structural cross section of flow phantom displayed in log scale. The corresponding cross-sectional angiograms (b) AD-Angio-OCT and (c) CD-Angio-OCT. The asterisk and cross indicate the selected two points from the static and dynamic regions, respectively. (b, c) The yellow arrows indicate the different performance between AD and CD algorithms. Four adjacent frames are averaged for presentation.

The temporal statistics of the OCT amplitude signals at the selected points are reported in Fig. 4. The histograms present the statistical distributions of the measured data. The variances σd2 and σs2 are computed from the measured data. Substituting the variances σd2 and σs2 into Eqs. (6) and (11) yields the theoretical predictions of OCT statistics, as plotted by the dashed curves in Fig. 4. The high R-square statistics (R2>0.95) indicate the well-matched agreement between the experimental outcome and the theoretical model.

Graphic Jump Location
Fig. 4
F4 :

Temporal statistics of the OCT amplitude signals at the (a) dynamic and (b) static positions, which are indicated by the cross and asterisk in Fig. 3(a), respectively. The histograms are the statistical distributions of the measured data. The dashed curves correspond to the theoretical predictions. The high R-square statistics (R2>0.95) indicate that the experimental outcome matches well with the theoretical model. The variances σd2 and σs2 are computed from the measured OCT data.

Figure 5 reports the statistical analysis of the Angio-OCT signals. Figures 5(a) and 5(b) correspond to the dynamic and static regions in AD-Angio-OCT, respectively. Figures 5(c) and 5(d) correspond to the dynamic and static regions in CD-Angio-OCT, respectively. The histograms and solid curves show the experimentally measured and theoretically predicted distributions of the Angio-OCT signals, respectively. The theoretical predictions are generated by substituting the variances σd2 and σs2 into Eqs. (16), (20), (24), and (25). The R-square statistics (R2) are higher than 0.95, indicating good agreement between the experimental outcome and the theoretical prediction.

Graphic Jump Location
Fig. 5
F5 :

Statistics of the Angio-OCT signals in flow phantom imaging. (a, b)The statistics of the dynamic and static signals in AD-Angio-OCT, respectively. (c, d) The statistics of the dynamic and static signals in CD-Angio-OCT, respectively. The R-square statistics (R2) are higher than 0.95.

In Vivo Brain Imaging

Figure 6(a) is a representative structural cross section of the mouse brain in vivo, in which it is challenging to discriminate the dynamic blood flow from the static tissue bed. Figures 6(b) and 6(c) show the corresponding cross-sectional angiograms of AD-Angio-OCT and CD-Angio-OCT, respectively. Similar to the phantom experiment, in the recorded cross section, two spatial positions were randomly selected from the regions of static tissue bed (zs,xs) and dynamic blood flow (zd,xd), respectively, for the following statistical analysis. The variances σd2 and σs2 are computed from the measured data, and then substituted into Eqs. (6), (11), (16), (20), (24), and (25) for predicting the OCT and Angio-OCT statistics.

Graphic Jump Location
Fig. 6
F6 :

(a) Representative structural cross section of mouse brain in vivo displayed in log scale. (a, b) The corresponding cross-sectional angiograms are produced by AD-Angio-OCT and CD-Angio-OCT algorithms, respectively. The asterisk and cross indicate the selected two points from the static tissue bed and the dynamic blood flow, respectively. (b, c) The yellow arrows indicate the different performance between AD and CD algorithms. Four adjacent frames are averaged for presentation.

The temporal statistics of the OCT amplitude signals at the selected static (zs,xs) and dynamic (zd,xd) positions are reported in Fig. 7. The histograms show the experimentally measured distributions. The dashed curves in Fig. 7 correspond to the theoretical predictions. The R-square statistics (R2) are higher than 0.95.

Graphic Jump Location
Fig. 7
F7 :

Temporal statistics of the OCT amplitude signals from (a) dynamic blood flow and (b) the static tissue bed in mouse brain in vivo, which are marked by the cross and asterisk in Fig. 6(a), respectively. The histograms are the distributions of the experimental data. The dashed curves in (a) and (b) correspond to the theoretical fitting using Eqs. (6) and (11). The R-square statistics (R2) are higher than 0.95. The variances σd2 and σs2 are computed from the measured OCT data.

Figure 8 reports the temporal statistics of the Angio-OCT signals of the in vivo experiment. Figures 8(a) and 8(b) correspond to the blood flow and static tissue bed in AD-Angio-OCT, respectively. Figures 8(c) and 8(d) correspond to the blood flow and static tissue bed in CD-Angio-OCT, respectively. The histograms and the solid curves show the experimentally measured and theoretically predicted distributions of the Angio-OCT signals, respectively. The R-square statistics (R2) are higher than 0.95. The experimental data match well with the theoretical predictions of the proposed statistical models.

Graphic Jump Location
Fig. 8
F8 :

Statistics of the Angio-OCT signals in mouse brain in vivo imaging. (a, b) The statistics of the dynamic and static signals in AD-Angio-OCT, respectively. (c, d) The statistics of the dynamic and static signals in CD-Angio-OCT, respectively. The R-square statistic (R2) is higher than 0.95.

Based on the model of random phasor sums, the temporal statistics of the complex-valued OCT signals were mathematically described, as expressed in Eqs. (5), (6), (9), and (11). The dynamic and static signals exhibit intrinsic differences in the statistics of OCT amplitude. Despite the differences, the distributions of the dynamic and static signals still have a large overlap (referring to the error areas in Fig. 9), making it challenging to directly separate the dynamic blood flow from the static tissue bed in the OCT structural image [Figs. 3(a) and 6(a)]. In Angio-OCT, extra processing algorithms, such as the widely used AD and CD algorithms, are applied to the original OCT signals. Using mathematical transformations and reasonable approximations, the statistics of the AD-Angio-OCT and CD-Angio-OCT signals were further derived. In Angio-OCT, the overlap between the distributions of dynamic and static signals is greatly reduced [referring to the error areas in Figs. 9(b) and 9(c)]. The origin of motion–contrast in Angio-OCT has been mathematically explained in this work.

Graphic Jump Location
Fig. 9
F9 :

Normalized statistical distributions in OCT: (a) AD-Angio-OCT, (b) CD-Angio-OCT, and (c) σd/σs=4.

The proposed statistical model can be used for guiding the threshold determination in Angio-OCT. Currently, the threshold is set empirically, and the signals above the threshold are classified as the dynamic regions. The ratio of the misclassified signals can be defined as the classification error rate (CER), i.e., dynamic signals below the threshold plus the static signals above it. In Angio-OCT, the minimal CER is determined by the residual overlap between the distributions of dynamic and static signals. Accordingly, the cross point of the two distribution curves can be considered as the optimal threshold. Any offset from the optimal value would lead to an increased CER. The optimal thresholds in the AD-Angio-OCT (TAD) and the CD-Angio-OCT (TCD) should meet the following conditions: Display Formula

{f|AADs|(TAD)=f|AADd|(TAD)inADfACDs(TCD)=fACDd(TCD)inCD.(26)
Substituting Eqs. (16), (20), (24), and (25) into the expressions above, we obtain: Display Formula
TAD=2σdσs1σd22σs2ln(σd2σs),(27)
Display Formula
TCD=2σdσs2σd2σs2ln(σdσs).(28)
Accordingly, the thresholds of the AD- and CD-Angio-OCT are determined as indicated by the dashed lines in Fig. 9. Then the minimal CER of the AD-Angio-OCT (CERAD) and the CD-Angio-OCT (CERCD) can be calculated: Display Formula
2·CERAD=TADf|AADs|(aADs)daADs+0TADf|AADd|(aADd)daADd=1+erf[2k22ln(k2)]erf[k2ln(k)k21],(29)
Display Formula
2·CERCD=TCDfACDs(aCDs)daCDs+0TCDfACDd(aCDd)daCDd=1exp[2ln(k)k21]+exp[2k2ln(k)k21],(30)
where we define k=σd/σs. In practice, a set of training data can be collected using MB-mode scanning protocol in the region of interest. The empirical threshold is first used to separate the dynamic and static regions. Based on the dynamic and static data, the parameters of σd2 and σs2 are learned, and the theoretical thresholds are determined. Initial proof of concept of the threshold determination was validated in this study. As shown in Figs. 3 and 6, the theoretical thresholds work well in the homogeneous tissues. Tissues of different scattering properties correspond to different parameters σd2 and σs2. In Fig. 6, the cortex was used for parameter learning and threshold determination, and consequently there exist apparent classification errors in the cranium using the threshold of cortex.

Although it has been recognized that the method using the complex-valued signals offers higher motion-contrast by combining both the amplitude and phase information, the performance of the CD- and AD-Angio-OCT can be further understood from the theoretical model. As reported in Fig. 10, the CD-Angio-OCT shows a lower CER in most situations, i.e., a superior motion-contrast, which can be confirmed in Figs. 3(b), 3(c), 6(b), and 6(c) as indicated by the yellow arrows. In Fig. 6, the averaged CERs of the AD and CD methods are 0.26 and 0.12 in the region of cortex, respectively. However, it should be noted that the CD algorithm is extremely sensitive to the phase fluctuation, and consequently poses a high requirement for the system phase stability and a large computational load for the phase compensation. The phase compensation works well on the situations, such as flow phantoms and stable animal models, but it is challenging in the clinical circumstances. Fortunately, several motion–tracking techniques have been developed in ophthalmic OCT systems for motion correction.

Graphic Jump Location
Fig. 10
F10 :

Numerical comparison of classification error rate (CER) between AD-Angio-OCT and CD-Angio-OCT. Complex differential (CD) method has a lower CER than AD.

Averaging is widely used in Angio-OCT for high contrast. The developed model is helpful for guiding the design of the averaging approaches. According to the model, averaging of independent angiograms offers reduced CER and improved motion-contrast. Assuming two independent angiograms with the PDF fAngioOCT expressed by Eqs. (16), (20), (24), and (25), the PDF of the averaged Angio-OCT signals fAngioOCT¯ has a simple relation with the PDF of the original signals: Display Formula

fAngioOCT¯=fAngioOCT*fAngioOCT,(31)
where * represents the convolution computation. Figure 11 reports the normalized statistical distributions of the averaged Angio-OCT signals. Compared with the original distribution in Figs. 9(b) and 9(c), the averaged signals show a lower CER. If the angiograms are totally dependent, no contrast improvement can be obtained. Because the Angio-OCT signals of either the dynamic or static regions are random in the time dimension, repeated angiograms with a large time interval (T) can be considered as totally independent and can be used for averaging. However, the repeated imagings are performed in B-scans in typical Angio-OCT,10,13,16 and consequently, the large time interval leads to an increased imaging time which is not desired due to the influence of bulk motion. In contrast, similar to the averaging approaches used for speckle reduction,31 independent angiograms can also be achieved by methods such as wavelength diversity, angular diversity, and polarization diversity, and it can be explained in theory that the split-spectrum algorithm offers improved contrast.10

Graphic Jump Location
Fig. 11
F11 :

Normalized statistical distributions of the averaged (a) AD-Angio-OCT and (b) CD-Angio-OCT signals. σd/σs=4.

In Eqs. (13), (22), and (23), it is assumed that the dynamic signals in the n’th and (n+1)th B-frames are totally uncorrelated and independent under the condition of a sufficient time interval (t). Typically, the time interval (t) is determined by the B-frame rate in the interframe Angio-OCT,10,13,16 which is 5.3 ms in our system. Such a time interval is sufficient for the fast blood flow, but not for the slow one. Taking the approximation where the lateral resolution of OCT is 15μm, and the velocity of RBC in the capillaries is 1mm/s,32 the required t is around 15 ms. According to the proposed model, in spite of being sensitive to the slow motion, the Angio-OCT with a short time interval would result in an increased CER and limited motion-contrast. Thus, there exists a tradeoff between the imaging speed and motion-contrast in Angio-OCT.

Angio-OCT suffers from shadow artifacts extending below the vessels, and the artifacts frustrate the automated 3-D analysis of vascular networks.1 Due to the forward scattering of RBC, the static signals below the vessels are influenced by the dynamic multiple-scattered signals from the blood flow, and present a Rician distribution.30 The statistical differences between the shadow and flow areas can be analyzed, which may be helpful for suppressing the shadow artifacts.

Although the mathematical derivation is focused on the AD and CD Angio-OCT in this study, it can be transferred to the phase-based methods. According to Eqs. (5) and (9), the temporal PDFs of phase θd and θs follows uniform and Gaussian distributions, respectively.24 Then the statistics of the phase-based Angio-OCT can be deduced.

There are several limitations in the current model. First, in the dynamic regions, it is assumed that a large number of scattering RBCs are randomly distributed within the OCT voxel of interest, and the interframe OCT signals contributed by RBC are totally independent. The assumption can be well satisfied in the large blood vessels, but not in the capillaries. In the capillaries, RBC flow one by one at a slow speed (<1mm/s).32 The temporal statistics of the capillary signals in Angio-OCT will be investigated in future study. Second, in the static regions, we assume that the OCT scattering signal is a strong phasor and remains constant in all the B-frame for a given spatial point, referring to Eqs. (17), (22), and (23). When the bulk motion happens, there exists a relative change of the strong phasor. In particular, the boundaries of the layered tissues have great changes in both amplitude and phase. As shown in Figs. 6(b) and 6(c), obvious artifacts can be observed in the boundaries of the cranium, and the CD method shows more artifacts in the boundaries due to the considerable phase changes.

Based on the model of random phasor sums, the temporal statistics of the complex-valued OCT signals were mathematically described. Using mathematical transformations and reasonable approximations, the temporal statistics of AD- and CD-Angio-OCT signals were derived and were found to obey different statistical distributions. The theories were further validated through both the flow phantom and live animal experiments. Using the model developed in this work, the origin of the motion–contrast in Angio-OCT is mathematically explained, and the possible implications in the improvement of motion–contrast are further discussed, including the threshold determination and its residual classification error, averaging method, and scanning protocol. The CD-Angio-OCT shows a lower CER than the AD method when the phase compensation works well. The proposed mathematical model of Angio-OCT signals can aid in the optimal design of the system and associated algorithms.

We acknowledge financial supports from National Natural Science Foundation of China (Nos. 61475143, 11404285, 61335003, 61327007, and 61275196), Zhejiang Provincial Natural Science Foundation of China (No. LY14F050007), National Hi-Tech Research and Development Program of China (No. 2015AA020515), Zhejiang Province Science and Technology Grant (No. 2015C33108), Fundamental Research Funds for the Central Universities (No. 2014QNA5017), and Scientific Research Foundation for Returned Scholars, Ministry of Education of China.

Vakoc  B. J.  et al., “Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging,” Nat. Med.. 15, , 1219 –1223 (2009).CrossRef
Mariampillai  A.  et al., “Speckle variance detection of microvasculature using swept-source optical coherence tomography,” Opt. Lett.. 33, , 1530 –1532 (2008).CrossRef
Mariampillai  A.  et al., “Optimized speckle variance OCT imaging of microvasculature,” Opt. Lett.. 35, , 1257 –1259 (2010).CrossRef
Barton  J., and Stromski  S., “Flow measurement without phase information in optical coherence tomography images,” Opt. Express. 13, , 5234 –5239 (2005).CrossRef
Motaghiannezam  S. M. R., , Koos  D., and Fraser  S. E., “Differential phase-contrast, swept-source optical coherence tomography at 1060 nm for in vivo human retinal and choroidal vasculature visualization,” J. Biomed. Opt.. 17, (2 ), 026011  (2012).CrossRef
Fingler  J.  et al., “Mobility and transverse flow visualization using phase variance contrast with spectral domain optical coherence tomography,” Opt. Express. 15, , 12636 –12653 (2007).CrossRef
Yu  L., and Chen  Z., “Doppler variance imaging for three-dimensional retina and choroid angiography,” J. Biomed. Opt.. 15, (1 ), 016029  (2010).CrossRef
Liu  G.  et al., “Intensity-based modified Doppler variance algorithm: application to phase instable and phase stable optical coherence tomography systems,” Opt. Express. 19, , 11429 –11440 (2011).CrossRef
Makita  S.  et al., “Optical coherence angiography,” Opt. Express. 14, , 7821 –7840 (2006).CrossRef
Jia  Y.  et al., “Split-spectrum amplitude-decorrelation angiography with optical coherence tomography,” Opt. Express. 20, , 4710 –4725 (2012).CrossRef
Enfield  J., , Jonathan  E., and Leahy  M., “In vivo imaging of the microcirculation of the volar forearm using correlation mapping optical coherence tomography (cmOCT),” Biomed. Opt. Express. 2, , 1184 –1193 (2011).CrossRef
Wang  R. K.  et al., “Three dimensional optical angiography,” Opt. Express. 15, , 4083 –4097 (2007).CrossRef
Wang  R. K.  et al., “Depth-resolved imaging of capillary networks in retina and choroid using ultrahigh sensitive optical microangiography,” Opt. Lett.. 35, , 1467 –1469 (2010).CrossRef
Liu  G.  et al., “A comparison of Doppler optical coherence tomography methods,” Biomed. Opt. Express. 3, , 2669 –2680 (2012).CrossRef
Leitgeb  R. A.  et al., “Doppler optical coherence tomography,” Prog. Retinal Eye Res.. 41, , 26 –43 (2014).CrossRef
An  L., , Shen  T. T., and Wang  R. K., “Using ultrahigh sensitive optical microangiography to achieve comprehensive depth resolved microvasculature mapping for human retina,” J. Biomed. Opt.. 16, (10 ), 106013  (2011).CrossRef
Lee  J.  et al., “Motion correction for phase-resolved dynamic optical coherence tomography imaging of rodent cerebral cortex,” Opt. Express. 19, , 21258 –21270 (2011).CrossRef
Choi  W. J.  et al., “Improved microcirculation imaging of human skin in vivo using optical microangiography with a correlation mapping mask,” J. Biomed. Opt.. 19, (3 ), 036010  (2014).CrossRef
Fingler  J.  et al., “Volumetric microvascular imaging of human retina using optical coherence tomography with a novel motion contrast technique,” Opt. Express. 17, , 22190 –22200 (2009).CrossRef
Liu  G.  et al., “Real-time bulk-motion-correction free Doppler variance optical coherence tomography for choroidal capillary vasculature imaging,” Opt. Express. 19, , 3657 –3666 (2011).CrossRef
Liew  Y. M.  et al., “In vivo assessment of human burn scars through automated quantification of vascularity using optical coherence tomography,” J. Biomed. Opt.. 18, (6 ), 061213  (2013).CrossRef
Mahmud  M. S.  et al., “Review of speckle and phase variance optical coherence tomography to visualize microvascular networks,” J. Biomed. Opt.. 18, , 050901  (2013).CrossRef
Lozzi  A.  et al., “Image quality metrics for optical coherence angiography,” Biomed. Opt. Express. 6, (7 ), 2435 –2447 (2015).CrossRef
Goodman  J. W., Statistical Optics. , p. 567 ,  Wiley-Interscience ,  New York  (1985).
Li  P.  et al., “Pulsatile motion of the trabecular meshwork in healthy human subjects quantified by phase-sensitive optical coherence tomography,” Biomed. Opt. Express. 4, , 2051 –2065 (2013).CrossRef
Zotter  S.  et al., “Visualization of microvasculature by dual-beam phase-resolved Doppler optical coherence tomography,” Opt. Express. 19, , 1217 –1227 (2011).CrossRef
Yang  V. X. D.  et al., “Improved phase-resolved optical Doppler tomography using the Kasai velocity estimator and histogram segmentation,” Opt. Commun.. 208, , 209 –214 (2002).CrossRef
Steel  R. G. D., and Torrie  J. H., Principles and Procedures of Statistics with Special Reference to the Biological Sciences. ,  McGraw-Hill Book Co. ,  New York, London  (1960).
Motaghiannezam  R., and Fraser  S., “Logarithmic intensity and speckle-based motion contrast methods for human retinal vasculature visualization using swept source optical coherence tomography,” Biomed. Opt. Express. 3, (3 ), 503 –521 (2012).CrossRef
Cheng  K. H.  et al., “Histogram flow mapping with optical coherence tomography for in vivo skin angiography of hereditary hemorrhagic telangiectasia,” J. Biomed. Opt.. 19, (8 ), 086015  (2014).CrossRef
Drexler  W., and Fujimoto  J. G., Optical Coherence Tomography: Technology and Applications. ,  Springer ,  Berlin, London  (2008).
Kleinfeld  D.  et al., “Fluctuations and stimulus-induced changes in blood flow observed in individual capillaries in layers 2 through 4 of rat neocortex,” Proc. Natl. Acad. Sci. U. S. A.. 95, , 15741 –15746 (1998).CrossRef
© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation

Yuxuan Cheng ; Li Guo ; Cong Pan ; Tongtong Lu ; Tianyu Hong, et al.
"Statistical analysis of motion contrast in optical coherence tomography angiography", J. Biomed. Opt. 20(11), 116004 (Nov 02, 2015). ; http://dx.doi.org/10.1117/1.JBO.20.11.116004


Figures

Graphic Jump Location
Fig. 1
F1 :

Flow chart of the statistical analysis in optical coherence tomography angiography (Angio-OCT). Amplitude differential (AD) method is used as an example for illustration. (a) Three-dimensional (3-D) OCT data cube (z,x,n). (b) Temporal distribution of the complex-valued OCT signals. (c) Statistical distribution of the AD-Angio-OCT signals.

Graphic Jump Location
Fig. 2
F2 :

Numerical simulation of AD-Angio-OCT signals in the dynamic region. The solid lines are the numerical simulation of Eq. (14) with different variances. The dash lines are the corresponding Gaussian distributions of Eq. (15). The R-square value (R2) is 0.99.

Graphic Jump Location
Fig. 3
F3 :

(a) Representative structural cross section of flow phantom displayed in log scale. The corresponding cross-sectional angiograms (b) AD-Angio-OCT and (c) CD-Angio-OCT. The asterisk and cross indicate the selected two points from the static and dynamic regions, respectively. (b, c) The yellow arrows indicate the different performance between AD and CD algorithms. Four adjacent frames are averaged for presentation.

Graphic Jump Location
Fig. 4
F4 :

Temporal statistics of the OCT amplitude signals at the (a) dynamic and (b) static positions, which are indicated by the cross and asterisk in Fig. 3(a), respectively. The histograms are the statistical distributions of the measured data. The dashed curves correspond to the theoretical predictions. The high R-square statistics (R2>0.95) indicate that the experimental outcome matches well with the theoretical model. The variances σd2 and σs2 are computed from the measured OCT data.

Graphic Jump Location
Fig. 5
F5 :

Statistics of the Angio-OCT signals in flow phantom imaging. (a, b)The statistics of the dynamic and static signals in AD-Angio-OCT, respectively. (c, d) The statistics of the dynamic and static signals in CD-Angio-OCT, respectively. The R-square statistics (R2) are higher than 0.95.

Graphic Jump Location
Fig. 6
F6 :

(a) Representative structural cross section of mouse brain in vivo displayed in log scale. (a, b) The corresponding cross-sectional angiograms are produced by AD-Angio-OCT and CD-Angio-OCT algorithms, respectively. The asterisk and cross indicate the selected two points from the static tissue bed and the dynamic blood flow, respectively. (b, c) The yellow arrows indicate the different performance between AD and CD algorithms. Four adjacent frames are averaged for presentation.

Graphic Jump Location
Fig. 7
F7 :

Temporal statistics of the OCT amplitude signals from (a) dynamic blood flow and (b) the static tissue bed in mouse brain in vivo, which are marked by the cross and asterisk in Fig. 6(a), respectively. The histograms are the distributions of the experimental data. The dashed curves in (a) and (b) correspond to the theoretical fitting using Eqs. (6) and (11). The R-square statistics (R2) are higher than 0.95. The variances σd2 and σs2 are computed from the measured OCT data.

Graphic Jump Location
Fig. 8
F8 :

Statistics of the Angio-OCT signals in mouse brain in vivo imaging. (a, b) The statistics of the dynamic and static signals in AD-Angio-OCT, respectively. (c, d) The statistics of the dynamic and static signals in CD-Angio-OCT, respectively. The R-square statistic (R2) is higher than 0.95.

Graphic Jump Location
Fig. 9
F9 :

Normalized statistical distributions in OCT: (a) AD-Angio-OCT, (b) CD-Angio-OCT, and (c) σd/σs=4.

Graphic Jump Location
Fig. 10
F10 :

Numerical comparison of classification error rate (CER) between AD-Angio-OCT and CD-Angio-OCT. Complex differential (CD) method has a lower CER than AD.

Graphic Jump Location
Fig. 11
F11 :

Normalized statistical distributions of the averaged (a) AD-Angio-OCT and (b) CD-Angio-OCT signals. σd/σs=4.

Tables

Table Grahic Jump Location
Table 1Statistics in optical coherence tomography and Angio-OCT.

References

Vakoc  B. J.  et al., “Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging,” Nat. Med.. 15, , 1219 –1223 (2009).CrossRef
Mariampillai  A.  et al., “Speckle variance detection of microvasculature using swept-source optical coherence tomography,” Opt. Lett.. 33, , 1530 –1532 (2008).CrossRef
Mariampillai  A.  et al., “Optimized speckle variance OCT imaging of microvasculature,” Opt. Lett.. 35, , 1257 –1259 (2010).CrossRef
Barton  J., and Stromski  S., “Flow measurement without phase information in optical coherence tomography images,” Opt. Express. 13, , 5234 –5239 (2005).CrossRef
Motaghiannezam  S. M. R., , Koos  D., and Fraser  S. E., “Differential phase-contrast, swept-source optical coherence tomography at 1060 nm for in vivo human retinal and choroidal vasculature visualization,” J. Biomed. Opt.. 17, (2 ), 026011  (2012).CrossRef
Fingler  J.  et al., “Mobility and transverse flow visualization using phase variance contrast with spectral domain optical coherence tomography,” Opt. Express. 15, , 12636 –12653 (2007).CrossRef
Yu  L., and Chen  Z., “Doppler variance imaging for three-dimensional retina and choroid angiography,” J. Biomed. Opt.. 15, (1 ), 016029  (2010).CrossRef
Liu  G.  et al., “Intensity-based modified Doppler variance algorithm: application to phase instable and phase stable optical coherence tomography systems,” Opt. Express. 19, , 11429 –11440 (2011).CrossRef
Makita  S.  et al., “Optical coherence angiography,” Opt. Express. 14, , 7821 –7840 (2006).CrossRef
Jia  Y.  et al., “Split-spectrum amplitude-decorrelation angiography with optical coherence tomography,” Opt. Express. 20, , 4710 –4725 (2012).CrossRef
Enfield  J., , Jonathan  E., and Leahy  M., “In vivo imaging of the microcirculation of the volar forearm using correlation mapping optical coherence tomography (cmOCT),” Biomed. Opt. Express. 2, , 1184 –1193 (2011).CrossRef
Wang  R. K.  et al., “Three dimensional optical angiography,” Opt. Express. 15, , 4083 –4097 (2007).CrossRef
Wang  R. K.  et al., “Depth-resolved imaging of capillary networks in retina and choroid using ultrahigh sensitive optical microangiography,” Opt. Lett.. 35, , 1467 –1469 (2010).CrossRef
Liu  G.  et al., “A comparison of Doppler optical coherence tomography methods,” Biomed. Opt. Express. 3, , 2669 –2680 (2012).CrossRef
Leitgeb  R. A.  et al., “Doppler optical coherence tomography,” Prog. Retinal Eye Res.. 41, , 26 –43 (2014).CrossRef
An  L., , Shen  T. T., and Wang  R. K., “Using ultrahigh sensitive optical microangiography to achieve comprehensive depth resolved microvasculature mapping for human retina,” J. Biomed. Opt.. 16, (10 ), 106013  (2011).CrossRef
Lee  J.  et al., “Motion correction for phase-resolved dynamic optical coherence tomography imaging of rodent cerebral cortex,” Opt. Express. 19, , 21258 –21270 (2011).CrossRef
Choi  W. J.  et al., “Improved microcirculation imaging of human skin in vivo using optical microangiography with a correlation mapping mask,” J. Biomed. Opt.. 19, (3 ), 036010  (2014).CrossRef
Fingler  J.  et al., “Volumetric microvascular imaging of human retina using optical coherence tomography with a novel motion contrast technique,” Opt. Express. 17, , 22190 –22200 (2009).CrossRef
Liu  G.  et al., “Real-time bulk-motion-correction free Doppler variance optical coherence tomography for choroidal capillary vasculature imaging,” Opt. Express. 19, , 3657 –3666 (2011).CrossRef
Liew  Y. M.  et al., “In vivo assessment of human burn scars through automated quantification of vascularity using optical coherence tomography,” J. Biomed. Opt.. 18, (6 ), 061213  (2013).CrossRef
Mahmud  M. S.  et al., “Review of speckle and phase variance optical coherence tomography to visualize microvascular networks,” J. Biomed. Opt.. 18, , 050901  (2013).CrossRef
Lozzi  A.  et al., “Image quality metrics for optical coherence angiography,” Biomed. Opt. Express. 6, (7 ), 2435 –2447 (2015).CrossRef
Goodman  J. W., Statistical Optics. , p. 567 ,  Wiley-Interscience ,  New York  (1985).
Li  P.  et al., “Pulsatile motion of the trabecular meshwork in healthy human subjects quantified by phase-sensitive optical coherence tomography,” Biomed. Opt. Express. 4, , 2051 –2065 (2013).CrossRef
Zotter  S.  et al., “Visualization of microvasculature by dual-beam phase-resolved Doppler optical coherence tomography,” Opt. Express. 19, , 1217 –1227 (2011).CrossRef
Yang  V. X. D.  et al., “Improved phase-resolved optical Doppler tomography using the Kasai velocity estimator and histogram segmentation,” Opt. Commun.. 208, , 209 –214 (2002).CrossRef
Steel  R. G. D., and Torrie  J. H., Principles and Procedures of Statistics with Special Reference to the Biological Sciences. ,  McGraw-Hill Book Co. ,  New York, London  (1960).
Motaghiannezam  R., and Fraser  S., “Logarithmic intensity and speckle-based motion contrast methods for human retinal vasculature visualization using swept source optical coherence tomography,” Biomed. Opt. Express. 3, (3 ), 503 –521 (2012).CrossRef
Cheng  K. H.  et al., “Histogram flow mapping with optical coherence tomography for in vivo skin angiography of hereditary hemorrhagic telangiectasia,” J. Biomed. Opt.. 19, (8 ), 086015  (2014).CrossRef
Drexler  W., and Fujimoto  J. G., Optical Coherence Tomography: Technology and Applications. ,  Springer ,  Berlin, London  (2008).
Kleinfeld  D.  et al., “Fluctuations and stimulus-induced changes in blood flow observed in individual capillaries in layers 2 through 4 of rat neocortex,” Proc. Natl. Acad. Sci. U. S. A.. 95, , 15741 –15746 (1998).CrossRef

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging & repositioning the boxes below.

Related Book Chapters

Topic Collections

PubMed Articles
Advertisement
  • Don't have an account?
  • Subscribe to the SPIE Digital Library
  • Create a FREE account to sign up for Digital Library content alerts and gain access to institutional subscriptions remotely.
Access This Article
Sign in or Create a personal account to Buy this article ($20 for members, $25 for non-members).
Access This Proceeding
Sign in or Create a personal account to Buy this article ($15 for members, $18 for non-members).
Access This Chapter

Access to SPIE eBooks is limited to subscribing institutions and is not available as part of a personal subscription. Print or electronic versions of individual SPIE books may be purchased via SPIE.org.