Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

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

Open Access Open Access

Abstract

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 coefficient 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.

© 2016 Optical Society of America

1. Introduction

Optical coherence angiography (OCA) is a non-invasive vasculature contrast technique based on optical coherence tomography (OCT). OCA has been used in several applications, e.g., ophthalmology [1–5], oncology [6], neurology [7], and gastroenterology [8]. Several angiography methods and algorithms have been developed, including OCA based on Doppler phase-shift detection [9–11], high-pass filtering techniques [12–18], intensity variance [19], and phase-shift variance [20, 21]. Methods based on intensity correlation [22–24] and complex correlation [6, 25–27] have also been presented. Among these, the complex-correlation-based technique uses all of the complex information of OCT, and provides good vasculature contrast [27]. In addition, the correlation-based methods are sensitive to the mobility of scatterers and less sensitive to the Doppler angle [28].

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 signal-intensity-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 polarization-sensitive OCT that has multiple polarization detection channels and enables polarization-sensitive 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.

2. 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.

2.1. 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

ΣG(Δr;r)=[E[|G(r)|2]E[G(r)G*(r+Δr)]E[G(r+Δr)G*(r)]E[|G(r+Δr)|2]],
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
ρG(Δr;r)eiθG(Δr;r)=E[G(r)G*(r+Δr)]E[|G(r)|2]E[|G(r+Δr)|2]

The correlation coefficient ρG is estimated using the sample covariance matrix SG between the measured OCT signal pair g(r) and g(r + Δr) which is obtained as

SG(Δr;r)g(Δr;r)g(Δr;r)w1=[|g(r)|2g(r)g*(r+Δr)g*(r)g(r+Δr)|g(r+Δr)|2]w1,
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
rG(Δr;r)ρ^G(Δr;r)|sG12|sG11sG22,
where sGpq are the entry of the estimated covariance matrix SG 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:

G(r)=S(r)+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

ΣG(Δr;r)=[E[|S(r)|2]+E[|N|2]E[S*(r)S(r+Δr)]E[S*(r+Δr)S(r)]E[|S(r+Δr)|2]+E[|N|2].]

The correlation between measured OCT signals at different spatiotemporal points, ρGr;r), can then be written as the product of the correlation coefficient of the true OCT signals, ρSr;r), and a factor expressing the noise effect, ρSNRr;r), as [29]

ρG(Δr;r)=ρSNR(Δr;r)ρS(Δr;r),
where
ρS(Δr;r)=E[S*(r)S(r+Δr)]E[|S(r)|2]E[|S(r+Δr)|2]
ρSNR(Δr;r)=[1+SNR(r)11+SNR(r+Δr)1]1.

ρSNRr; r) [Eq. (10)] is a decorrelation factor caused by the random noise, where SNR(r) = E[|S(r)|2]/E[|N|2]. Hence, rG [Eq. (5)] does not provide an accurate estimate of the signal correlation ρS.

The bias brG and variance σrG in the correlation estimate rG are expressed as

brG=E[rG]ρS,
σrG=E[rG2]E[rG]2.

The bias decreases with higher SNRs. As SNR ∞, ρSNR 1 and ρG → ρS. Thus, the bias brG decreases. However, the variance of estimation is not always small at high SNRs. The variance σrG 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 σrG will results in the noise of the correlation imaging.

Tables Icon

Table 1. Variance of correlation coefficient estimation σrG.

2.2. SNR-corrected low-bias, low-noise complex correlation estimation

The correlation of signal ρSr;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

rS(Δr;r)|sG12||[sG11q^N][sG22q^N]|,
where sSpq are the entry of the estimated covariance matrix SS at p-th row and q-th column, and q^N is an estimation of the noise power, obtained as
q^N|gn(t)gn(t)t|2t,
where 〈 〉t denotes time-averaging. gn(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 power q^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 w2 as

r¯G(Δr;r)|sG12|w2sG11sG22w2.

Similarly, a low-noise signal correlation estimate is obtained from Eq. (13) as

r¯S(Δr;r)|sG12|w2ε|[sG11q^N][sG22q^N]|w2,
where
ε=sgn(sG11q^N)+sgn(sG22q^N)2.

The sign ε is applied to suppress the bias of the denominator of Eq. (16) after the averaging w2 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.

3. 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 ρSr;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(ri,rj)r¯S(τ;x,z). This can be described as

r¯S(τ;x,z)|g(r)g*(r+(0,0,τ))w1|w2ε|[|g(r)|2w1q^N][|g(r+(0,0,τ))|2w1q^N]|w2.

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.

 figure: Fig. 1

Fig. 1 Schematic block diagrams of cmOCA processing. (a) Overall processing. (b) Covariance matrix estimate with bulk-phase-offset correction.

Download Full Size | PDF

3.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 w1 in Eq. (18) is extended along the temporal direction. There are two approaches to overcome this problem. One is to apply a bulk-phase-offset correction to the estimation and the other is to form an estimation that is insensitive to the bulk-phase-offset.

3.1.1. 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. 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.

Tables Icon

Table 2. Correlation coefficient ρG for each tissue type and SNR.

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:

SG(τ;r)g(τ;r)g(τ;r)w1.
and the data correlation r¯G is then estimated using Eqs. (15) as
r¯G(τ;r)|sG12|w2sG11sG22w2.

To prevent the disruption caused by bulk-phase offset, the estimation w1 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 zmax(τ; x, t), which exhibits a high value of r¯G at x, as

Δϕm(τ;x,t)arg{sG12(τ;r|z=zmax(τ;x,t))w3}.

Using Eq. (21), the bulk-phase-offset is corrected and the covariance matrix is finally estimated as

SG(τ;x,z)[sG11sG12eiΔϕm(τ;x,t){sG12eiΔϕm(τ;x,t)}*sG22]w1.

The signal correlation map r¯S(τ;x,z) is estimated by substituting Eqs. (19) and (21) to (22) into Eq. (16) as

r¯S(τ;x,z)|g(r)g*(r+(0,0,τ))w1eiΔϕm(τ;x,τ)w1|w2ε|[|g(r)|2w1w1q^N][|g(r+(0,0,τ))|2w1w1q^N]|w2.

3.1.2. Phase-noise-immune estimate

In the case of swept-source OCT (SS-OCT) systems, the phase fluctuations caused by asynchronous acquisition [47] and/or random phase noise of coherence revival signals [48] make it difficult to acquire phase-stable signals. In addition, when no k-clock is used [49], the low repeatability of the wavelength sweep of a light source could further degrade the phase stability.

In such cases, another algorithm with phase-noise immunity is applied. The sample covariance matrix SG is estimated as similar to Eq. (3) as:

SG(τ;r)g(τ;r)g(τ;r)w1.

To minimize the effect of the phase fluctuation, the averaging window w1 is selected within time independent region as in Eq. (19). One example is axial averaging [27].

3.2. 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

ρS(τ;r)ρf(τ;r)ρm(τ;x,t),
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
r¯S(τ;x,z)r˜f(τ;x,z)rm(τ;x),
where r˜f(τ;x,z) is the net correlation coefficient corresponding to blood flow. The reduction of r¯S(τ;x,z) due to rm(τ; 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-noise-immune estimate (Section 3.1.2), decorrelation due to bulk-phase-offset is avoided. Our motion-artifact-correction method first estimates rm(τ; 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 rm(τ; 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 SG [Eq. (22) or (24)] in the manner of Eq. (15) as

r¯G(τ;x,z)|sG12|w2sG11sG22w2.

Hence, the signal correlation estimate r¯S at the depth of the maximum of r¯G(τ;x,z) is used as an estimate of rm:

r^m(τ;x)r¯S(τ;x,argmaxZr¯G(τ;x,z)).

The bulk-motion-corrected signal correlation coefficient is finally obtained as

r¯Sf(τ;x,z)r¯S(τ;x,z)r^m(τ;x).

This estimate represents the correlation due to flow, and is used to form bulk-motion-corrected cmOCA images.

4. 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).

4.1. 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 location on the sample is scanned multiple times by this scanning protocol. The correlation is computed at each spatial position among the frames.

 figure: Fig. 2

Fig. 2 Schematic diagram of applied scanning protocol.

Download Full Size | PDF

JM-OCT simultaneously provides four independent OCT signals [49]. These are denoted as g(x, z, f, p), where f is the index of frame and p = 1,2,3,4 is the index of the JM-OCT channels. Note that variables x and z are re-defined as sampling point indices along the lateral and axial directions, respectively. The time correlation between frames f and f + 1 at space point (x, z) of p channel is denoted as ρS(τ; x, z, f, p), where τ is the time separation between successive frames. The correlations, covariance matrices, and their estimates are denoted in the same manner.

4.2. 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 w1, the bulk-phase-offset correction (Section 3.1.1) is preferable to extend the window w1 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

SGM(τ;x,z,f,p)1DdDg(τ;x,z+d,f,p)g(τ;x,z+d,f,p)
r¯GM(τ;x,z,f)pP1q^N(p)|sG12M|pP1q^N(p)sG11MsG22M,
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¯GM in each axial profile to estimate the bulk-phase offset. The axial indices of these pixels are denoted as zkmax(τ;x,f), with the index k identifying these pixels as k = 1,⋯,K. As in Eq. (21), the bulk-phase-offset estimate is obtained from sG12M(τ;x,z,f,p) at depths zkmax(τ;x,f) as
ΔϕmM(τ;x,f)arg{Medx[1PKpPk=1KsG12M(τ;x,zkmax(τ;x,f),f,p)]},
where Medx[] 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 Medx. The bulk-phase-offset correction is applied to SGM [Eq. (30)] with Eq. (32). The covariance matrix estimate is obtained with lateral and frame averaging for w1. The correlation coefficient estimate is then described from Eq. (23) using the noise power of each channel q^N(p) as
r¯SM(τ;x,z)p1q^N(p)|l,d,fL,D,F1g(x+l,z+d,f,p)g*(x+l,z+d,f+1,p)eiΔϕmM(τ;x+l,f)|pε(x,z,p)q^N(p){|[l,d,fL,D,F1|g(x+l,z+d,f,p)|2LD(F1)q^N(p)]×[l,d,fL,D,F1|g(x+l,z+d,f+1,p)|2LD(F1)q^N(p)]|}12,
where L is the size of the averaging kernel along the lateral direction.

Tables Icon

Table 3. Assignment of averaging and variables for each imaging mode.

The bulk-motion-artifact correction (Section 3.2) is then applied to r¯SM(τ;x,z). For this correction, the estimated data correlation r¯GM(τ;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. (28) with r¯GM(τ;x,z) as r^mM(x)r¯SM(τ;x,argmaxzr¯GM(τ;x,z)), where r^mM 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

r¯SfM(τ;x,z)r¯SM(τ;x,z)r^mM(τ;x).

4.2.1. Single-channel mode without variance averaging

To compare with and without the second averaging w2, 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 gC(x,z,f)=1PpPg(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.

rSC(τ;x,z)|l,d,fL,D,F1gC(x+1,z+d,f)gC*(x+1,z+d,f+1)eiΔϕmC(τ;x+1,f)|ε(x,z){|[l,d,fL,D,F1|gC(x+1,z+d,f)|2LD(F1)P2pPq^N(p)]×[l,d,fL,D,F1|gC(x+1,z+d,f+1)|2LD(F1)P2pPq^N(p)]|}12,

The bulk-motion-artifact correction (Section 3.2) is then applied to rSC(τ;x,z). Without the noise power subtraction and the sign ε in the Eq. (35), rGC(τ;x,z) is obtained. The bulk-motion-induced net correlation is estimated using Eq. (28) as r^mC(x)rSC(τ;x,argmaxrGC(τ;x,z)z). The bulk-motion-corrected signal correlation with coherent composite is obtained as

rSfC(τ;x,z)rSC(τ;x,z)r^mC(τ;x).

5. Image formation

5.1. 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¯S1<1 are expressed as black.

5.2. En-face projection

The projection images of cmOCA are calculated as E(x)=i=0N{1M[r¯Sf(τ;x,zi)]}, where

M[x]={1,ifx1<1.x,otherwise.

The dynamic range of the projection image is set to the [1: 99] percentile of each image’s pixel value.

5.3. 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:

E(x)=n=0Njcj(x,zn)ms(j)Δz[ei=0njcj(x,zi)ma(j)ΔzL0],
where cj is the j-th contrast value. ma is an infinitesimally small color translation of rays, ms is the reflection spectrum of rays, and L0 is initial color of rays. Δz is the axial distance per unit data point. Analogous to Lambert–Beer’s law, c, ms, and ma 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:

E(x)=n=0N{1M[r¯Sf(τ;x,zn)]}T[W[ecPVn(x)mΔzL0]],
where cPVn(x)=i=0n{1M[mDOPU(x,zi)]} is cumulative polarization variance from the surface to the depth zn. This means that color of rays changed according to cPVn by e12cPVnmΔz rays are reflected at zn according to cmOCA r¯Sf, color changed again by e12cPVnmΔz, and reflected rays from every depth zn 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

m=[1.1670083.06700.236027.944000.2831.9320000][μm1].

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)V1 and

ecPVn(x)mΔz=Vdiag(ecPVn(x)λ1Δz,ecPVn(x)λ2Δz,ecPVn(x)λ3Δz,ecPVn(x)λ4Δz)V1.

L0=[L0,a0,b0,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 = z0. In this paper, (L0,a0,b0)=(84.17,13.95,82.72). The transformation W[] is wrapping an input color location into the realizable L*a*b* space as:

W[Lin]={WLin(vLin>0)Lin(vLin0),
where v = [0,cosθ, sinθ, −(ao cosθ + bo sinθ)] and
W=[10000cos(2θ)sin(2θ)ao+aocos(2θ)+bosin(2θ)0sin(2θ)cos(2θ)bobocos(2θ)+aosin(2θ)0001].

Here, we set θ = 0.165 rad, ao = 40.3, and bo = 40.9. The color change according to cPVn is represented in L*a*b* color space as shown in Fig. 3.

 figure: Fig. 3

Fig. 3 Color change trajectory in L*a*b* color space according to cumulative polarization variance cPV.

Download Full Size | PDF

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, L0, and W[].

6. 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.

6.1. 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 cross section of the retina of a volunteer’s (male, 35 years old) right eye which had no remarkable abnormalities was obtained [Fig. 4(a)], and two correlation images were calculated. The SNR-independent and -dependent cmOCAs are obtained by mapping the magnitude of the signal correlation coefficient [Eq. (34)] and data correlation coefficient, which is estimated similar to Eq. (34) using r¯GM as r¯GfM(X,Z)r¯GM(X,Z)/r^mM(X). With the estimated data correlation r¯GfM, both the perfused tissue and low-SNR region exhibit a low correlation, as with other correlation-based and phase-shift-based angiography methods [Fig. 4(b)]. With the signal correlation r¯SfM, 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 SNR-corrected correlation distinguishes perfused tissues from low-SNR regions without the use of an arbitrary masking threshold.

 figure: Fig. 4

Fig. 4 Cross sections of the human retina. (a) Coherence composite OCT intensity, (b) SNR-dependent cmOCA, (c) SNR-independent cmOCA. The scale bar indicates 200 μm.

Download Full Size | PDF

6.2. 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 mm2. 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.

 figure: Fig. 5

Fig. 5 Complex correlation mapping image of the human retina. cmOCA (b, d) with multiple polarization channels, as described in Section 4.2, exhibits lower spurious noise than (a, c) single-channel estimation with a coherent composite signal. The cross sections (a, b) and en-face slices (c, d) are taken from the locations indicated by the black arrows. Both cmOCAs used the same kernel size.

Download Full Size | PDF

In contrast, the multi-channel method exhibits high contrast and low noise vessel imaging owing to the averaging process [Fig. 5(b)]. The corresponding en-face slice [Fig. 5(d)] visualizes the fine abnormal vascular network better than the single-channel method [Fig. 5(c)]. This shows that the second-step averaging is preferred for high-contrast vasculature imaging.

6.3. Motion artifact compensation

The performance of bulk-motion-artifact removal and bulk-phase-offset correction was evaluated with in vivo human posterior eye imaging. The left eye of a healthy subject (male, 27 years old) was scanned (6 × 6 mm2). Figure 6 illustrates the effect of bulk-motion-artifact removal (Section 3.2), and compares bulk-phase-offset correction using the new correlation-based estimation method (Section 3.1.1) against that with the complex-mean estimation method [46].

 figure: Fig. 6

Fig. 6 cmOCA cross sections without (a, b) and with (c, d) bulk-motion artifact removal. Bulk-phase-offset correction has been applied using the complex-mean method (a, c) and correlation-based method (b, d). The corresponding en-face projection images are shown in (e)–(h).

Download Full Size | PDF

There is a thick blood vessel along the B-scan direction. With the previous average-based bulk-phase-offset estimation, the tissue above the thick blood vessel also has a low correlation coefficient [Fig. 6(a)]. This is because the erroneous average-based bulk-phase-offset estimation results in an incorrect phase-offset correction, and the phase-offset fluctuation between B-scan pairs reduces the correlation coefficient. Hence, applying the bulk-motion artifact removal with the average-based phase-offset correction over-corrects the correlation coefficient [Fig. 6(c)]. This effect reduces the contrast of thick blood vessels, as shown in Fig. 6(g) (yellow arrows).

With the new correlation-based phase-offset estimation method, the solid tissue exhibits a high correlation coefficient as expected [Fig. 6(b)]. The bulk-motion artifacts over the thick vessels are reduced. On the other hand, a small vessel overlaid on the thick vessel is retained (blue arrows), and the contrast of thick blood vessels in the en-face image is maintained well [Fig. 6(f), yellow arrows].

By applying bulk-motion artifact removal [Figs. 6(g) and (f)], artifacts appear as horizontal lines because the frame-by-frame bulk-motion decorrelation changes are well suppressed (green circles).

6.4. High-stability mode

In some cases, the phase of the swept-source OCT signal cannot be stabilized without a k-clock. The algorithm of Section 4.2 results in significant artifacts. In this case, high phase-noise-immunity mode (Section 3.1.2) can be applied. The details of the implementation is described in Appendix A.

A subject (male, 33 years old, right eye) was scanned (6 × 6 mm2) using JM-OCT, which has an unstable phase. As shown in Figs. 7(a) and 7(d), there are several artifacts that disturb the vasculature image. Using the correlation calculation in Eq. (44), these artifacts are significantly reduced, as shown in Figs. 7(b) and 7(e). Although the range of correlation has decreased to [0.3, 1.0], the contrast in the en-face projection is significantly improved. After motion artifact correction [Eq. (45)], horizontal line artifacts are suppressed (Fig. 7(f)). The high phase-noise-immunity mode is an alternative high-contrast vasculature imaging method for unstable phase systems.

 figure: Fig. 7

Fig. 7 cmOCA images with low-phase-stability SS-OCT system. The images were obtained using (a, d) the multichannel cmOCA algorithm (Section 4.2), (b, e) high phase-noise-immunity mode [Eq. (44)], and (c, f) high phase-noise-immunity mode with bulk-motion-artifact removal [Eq. (45)]. Cross sections (a, b, c) at the location indicated by the black arrow and en-face projections (d, e, f) are shown.

Download Full Size | PDF

6.5. Retinal capillary imaging

To show the performance of the retinal capillary imaging with the cmOCA method, the same data set of Section 6.1 is used. The retinal layers are segmented to visualize the capillary network of different retinal layers. The Iowa Reference Algorithms (Retinal Image Analysis Lab, Iowa Institute for Biomedical Imaging, Iowa City, IA) [58–60] were applied to the coherent composite OCT intensity volume. The segmentation results were overlaid on a cross-sectional cmOCA image [Fig. 8(a)]. The en-face projection images of cmOCA volume visualize the retinal capillary vasculature at the fovea. From the anterior part, the retinal nerve fiber layer (NFL) + ganglion cell layer (GCL) [Fig. 8(b)], inner plexiform and nuclear layers (IPL + INL) [Figs. 8(c)], outer plexiform and nuclear layers (OPL + ONL) [8(d)] are shown.

 figure: Fig. 8

Fig. 8 cmOCAs at the human fovea, (a) cross-sectional cmOCA image with segmented results, en-face projection images at (b) NFL + GCL, (c) IPL + INL, (d) OPL + ONL. The scanning area is 3 × 3 mm2. The scale bar indicates 200 μm.

Download Full Size | PDF

6.6. Eye disease case

A polypoidal choroidal vasculopathy (PCV) (female, 67 years old, left eye) is shown in Fig. 9. The composite projection image [Fig. 9(b)] combines the correlation mapping and DOPU data [53]. Details of the method are described in Section 5.3. This imaging method reveals whether vessels are covered or uncovered by multiple-scattering tissues, i.e., melanin. The orange-yellow color indicates that blood flow signals are unblocked (UB) by mDOPU signals (vasculature is in front of the pigmented tissues), whereas magenta-blue indicates that the flow 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 and cmOCA are assigned as blue and orange, respectively [57].

 figure: Fig. 9

Fig. 9 cmOCA images of polypoidal choroidal vasculopathy. En-face (a) OCT, (b) composition of cmOCA and mDOPU projection images. The scanning area is 6 × 6 mm2. The corresponding area is indicated by the yellow rectangle in fluorescein angiography (c). The cross-sectional (d) OCT + cmOCA and (e) mDOPU images at the center of the scan [yellow line in (b)] are shown. The scale bar indicates 500 μm.

Download Full Size | PDF

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.

7. 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 Xi and Xi+1 in consecutive sets (i-th and i + 1-th sets) was generated by using the method to generate correlated random numbers. Outcomes Zi and Yi of independent complex Gaussian variables with zero mean and standard deviation σY were generated. Then, Yi+1 was generated with the true correlation coefficient ρ as Yi+1 = ρYi + 1ρ2Zi. Yi and Yi+1 will be correlated with the correlation coefficient ρ. Another outcome Ni of independent complex Gaussian variable with zero mean and standard deviation σN was generated and added to provide Xi = Yi + Ni. By changing the ratio σYN, 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 10log10 σY2/σN2 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 w2. 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 w2 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.

 figure: Fig. 10

Fig. 10 The numerical simulation results of correlation coefficient estimates. (a, b, c) Mean and (d, e, f) standard deviation of each correlation estimator, i.e., data correlation (×), signal correlation with previous method using coherent composite (▲), and the new method (●) are plotted with several SNR. The set correlation coefficients are (a, d) 0.99, (b, e) 0.9, and (c, f) 0.8.

Download Full Size | PDF

In our method, we introduced the sign function ε [Eq. (17)] in the denominator of the noise-corrected 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 w2.

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 w2, 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.

 figure: Fig. 11

Fig. 11 Comparison of statistical properties at no signal region. Histograms of correlation estimates (a) without the sign ε, (b) without the sign ε and thresholding at 0 instead of absolute operation, and (c) with the sign ε Eq. (16) are shown. The numerical simulation of correlation estimates was applied at SNR=0.

Download Full Size | PDF

8. 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.

Appendix A. High-stability mode implementation

To implement the phase-noise-immune estimate (Section 3.1.2), no lateral and frame averaging is used to estimate the covariance matrix [Eq. (24)] because both depend on time in the system with scanning beam. Instead, polarization channel and axial averaging are used.

The correlation coefficient in the absence of phase noise is then obtained by substituting the averaging operators and noise power as shown in Table 4 into Eq. (18) as

r¯S(τ;x,z)l,fL,F1|1Dp,dP,Dg(x+l,z+d,f,p)g*(x+l,z+d,f+1,p)|l,fL,F1ε(τ;x+l,z,f){|[1Dp,dP,D|g(x+l,z+d,f,p)|2pPq^N(p)]×[1Dp,dP,D|g(x+l,z+d,f+1,p)|2pPq^N(p)]|}12.

Tables Icon

Table 4. Assignment of averaging and variables of phase-noise-immune mode.

The bulk-motion-artifact correction (Section 3.2) is then applied to r¯S(τ;x,z). For this correction, the phase-noise-immune data-correlation-coefficient r¯G(τ;x,z) is estimated using Eq. (44) without noise subtraction and ε. The net bulk-motion-induced phase-noise-immune data-correlation is estimated using Eq. (28) with r¯G(τ;x,z) as r¯m(τ;x) r¯S(τ;x,argmaxzr¯G(τ;x,z)). The bulk-motion-corrected phase-noise-immune signal correlation coefficient is finally estimated as

r¯Sf(τ;x,z)r¯S(τ;x,z)r^m(τ;x,).

Appendix B. List of notations

SymbolsDescriptions
Gcomplex random variable represents OCT data
Scomplex random variable represents OCT signal
Ncomplex random variable represents additive noise
grealization of OCT measurement
ΣGpopulation covariance matrix of OCT data
SGsample covariance matrix of OCT data
sGpqentries of SG
q^Nestimate of variance of noise N
ρGpopulation correlation coefficient of OCT data
ρSpopulation correlation coefficient of OCT signals
ρSNRpopulation de-correlation factor of OCT signals due to additive noise ρGS
rGsample correlation coefficient of OCT data
rSsample correlation coefficient of OCT signal
r¯Gestimate of data correlation coefficient
r¯Sestimate of signal correlation coefficient
SNRsignal-to-noise ratio of OCT system E [|S|2]/E[|N|2]

Δϕmestimate of bulk-phase-offset
rmdecorrelation factor due to bulk motion
r^mestimate of decorrelation factor rm
r˜fnet decorrelation due to flow
r¯Sfestimate of signal correlation coefficient after bulk-motion artifact removal

pindex of polarization channel (1,2,⋯,P)
findex of frame (1,2,⋯,F)
Daveraging window size along axial direction
Laveraging window size along lateral direction

Acknowledgments

We acknowledge the cooperation of Satoshi Sugiyama from Tomey Corporation (Nagoya, Japan).

References and links

1. Y.-J. Hong, M. Miura, M. J. Ju, S. Makita, T. Iwasaki, and Y. Yasuno, “Simultaneous investigation of vascular and retinal pigment epithelial pathologies of exudative macular diseases by multifunctional optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 55, 5016–5031 (2014). [CrossRef]   [PubMed]  

2. Y. Jia, E. Wei, X. Wang, X. Zhang, J. C. Morrison, M. Parikh, L. H. Lombardi, D. M. Gattey, R. L. Armour, B. Edmunds, M. F. Kraus, J. G. Fujimoto, and D. Huang, “Optical coherence tomography angiography of optic disc perfusion in glaucoma,” Ophthalmology 121, 1322–1332 (2014). [CrossRef]   [PubMed]  

3. Y. Jia, S. T. Bailey, D. J. Wilson, O. Tan, M. L. Klein, C. J. Flaxel, B. Potsaid, J. J. Liu, C. D. Lu, M. F. Kraus, J. G. Fujimoto, and D. Huang, “Quantitative optical coherence tomography angiography of choroidal neovascularization in age-related macular degeneration,” Ophthalmology 121, 1435–1444 (2014). [CrossRef]   [PubMed]  

4. M. Miura, D. Muramatsu, Y.-J. Hong, Y. Yasuno, T. Iwasaki, and H. Goto, “Noninvasive vascular imaging of polypoidal choroidal vasculopathy by Doppler optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 56, 3179–3186 (2015). [CrossRef]   [PubMed]  

5. M. Miura, Y.-J. Hong, Y. Yasuno, D. Muramatsu, T. Iwasaki, and H. Goto, “Three-dimensional vascular imaging of proliferative diabetic retinopathy by Doppler optical coherence tomography,” Am. J. Ophthalmol. 159, 528–538 (2015). [CrossRef]  

6. B. J. Vakoc, R. M. Lanning, J. A. Tyrrell, T. P. Padera, L. A. Bartlett, T. Stylianopoulos, L. L. Munn, G. J. Tearney, D. Fukumura, R. K. Jain, and B. E. Bouma, “Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging,” Nat. Med. 15, 1219–1223 (2009). [CrossRef]   [PubMed]  

7. V. J. Srinivasan, D. N. Atochin, H. Radhakrishnan, J. Y. Jiang, S. Ruvinskaya, W. Wu, S. Barry, A. E. Cable, C. Ayata, P. L. Huang, and D. A. Boas, “Optical coherence tomography for the quantitative study of cerebrovascular physiology,” J. Cereb. Blood Flow Metab. 31, 1339–1345 (2011). [CrossRef]   [PubMed]  

8. T.-H. Tsai, O. O. Ahsen, H.-C. Lee, K. Liang, M. Figueiredo, Y. K. Tao, M. G. Giacomelli, B. M. Potsaid, V. Jayaraman, Q. Huang, A. E. Cable, J. G. Fujimoto, and H. Mashimo, “Endoscopic optical coherence angiography enables 3-dimensional visualization of subsurface microvasculature,” Gastroenterology 147, 1219–1221 (2014). [CrossRef]   [PubMed]  

9. S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, “Optical coherence angiography,” Opt. Express 14, 7821–7840 (2006). [CrossRef]   [PubMed]  

10. B. Braaf, K. A. Vermeer, V. A. D. Sicam, E. van Zeeburg, J. C. van Meurs, and J. F. de Boer, “Phase-stabilized optical frequency domain imaging at 1-μ m for the measurement of blood flow in the human choroid,” Opt. Express 19, 20886–20903 (2011). [CrossRef]   [PubMed]  

11. Y.-J. Hong, S. Makita, F. Jaillon, M. J. Ju, E. J. Min, B. H. Lee, M. Itoh, M. Miura, and Y. Yasuno, “High-penetration swept source Doppler optical coherence angiography by fully numerical phase stabilization,” Opt. Express 20, 2740–2760 (2012). [CrossRef]   [PubMed]  

12. R. K. Wang, “Three-dimensional optical micro-angiography maps directional blood perfusion deep within microcirculation tissue beds in vivo,” Phys. Med. Biol. 52, N531–N537 (2007). [CrossRef]   [PubMed]  

13. Y. K. Tao, A. M. Davis, and J. A. Izatt, “Single-pass volumetric bidirectional blood flow imaging spectral domain optical coherence tomography using a modified Hilbert transform,” Opt. Express 16, 12350–12361 (2008). [CrossRef]   [PubMed]  

14. V. J. Srinivasan, S. Sakadžić, I. Gorczynska, S. Ruvinskaya, W. Wu, J. G. Fujimoto, and D. A. Boas, “Quantitative cerebral blood flow with optical coherence tomography,” Opt. Express 18, 2477–2494 (2010). [CrossRef]   [PubMed]  

15. W. J. Choi, R. Reif, S. Yousefi, and R. K. Wang, “Improved microcirculation imaging of human skin in vivo using optical microangiography with a correlation mapping mask,” J. Biomed. Opt. 19, 036010 (2014). [CrossRef]  

16. C. Blatter, T. Klein, B. Grajciar, T. Schmoll, W. Wieser, R. Andre, R. Huber, and R. A. Leitgeb, “Ultrahigh-speed non-invasive widefield angiography,” J. Biomed. Opt. 17, 070505 (2012). [CrossRef]   [PubMed]  

17. V. J. Srinivasan, J. Y. Jiang, M. A. Yaseen, H. Radhakrishnan, W. Wu, S. Barry, A. E. Cable, and D. A. Boas, “Rapid volumetric angiography of cortical microvasculature with optical coherence tomography,” Opt. Lett. 35, 43–45 (2010). [CrossRef]   [PubMed]  

18. L. An, J. Qin, and R. K. Wang, “Ultrahigh sensitive optical microangiography for in vivo imaging of microcirculations within human skin tissue beds,” Opt. Express 18, 8220–8228 (2010). [CrossRef]   [PubMed]  

19. A. Mariampillai, B. A. Standish, E. H. Moriyama, M. Khurana, N. R. Munce, M. K. Leung, J. Jiang, A. Cable, B. C. Wilson, I. A. Vitkin, and V. X. D. Yang, “Speckle variance detection of microvasculature using swept-source optical coherence tomography,” Opt. Lett. 33, 1530–1532 (2008). [CrossRef]   [PubMed]  

20. J. Fingler, D. Schwartz, C. Yang, and S. E. Fraser, “Mobility and transverse flow visualization using phase variance contrast with spectral domain optical coherence tomography,” Opt. Express 15, 12636–12653 (2007). [CrossRef]   [PubMed]  

21. D. Y. Kim, J. Fingler, R. J. Zawadzki, S. S. Park, L. S. Morse, D. M. Schwartz, S. E. Fraser, and J. S. Werner, “Noninvasive imaging of the foveal avascular zone with high-speed, phase-variance optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 53, 85–92 (2012). [CrossRef]  

22. G. Liu, L. Chou, W. Jia, W. Qi, B. Choi, and Z. Chen, “Intensity-based modified Doppler variance algorithm: application to phase instable and phase stable optical coherence tomography systems,” Opt. Express 19, 11429–11440 (2011). [CrossRef]   [PubMed]  

23. Y. Jia, O. Tan, J. Tokayer, B. Potsaid, Y. Wang, J. J. Liu, M. F. Kraus, H. Subhash, J. G. Fujimoto, J. Hornegger, and D. Huang, “Split-spectrum amplitude-decorrelation angiography with optical coherence tomography,” Opt. Express 20, 4710–4725 (2012). [CrossRef]   [PubMed]  

24. J. Enfield, E. Jonathan, and M. Leahy, “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]   [PubMed]  

25. L. Wang, Y. Wang, S. Guo, J. Zhang, M. Bachman, G. Li, and Z. Chen, “Frequency domain phase-resolved optical Doppler and Deppler variance tomography,” Opt. Commun. 242, 345–350 (2005). [CrossRef]  

26. G. Liu, W. Qi, L. Yu, and Z. Chen, “Real-time bulk-motion-correction free Doppler variance optical coherence tomography for choroidal capillary vasculature imaging,” Opt. Express 19, 3657–3666 (2011). [CrossRef]   [PubMed]  

27. A. S. Nam, I. Chico-Calero, and B. J. Vakoc, “Complex differential variance algorithm for optical coherence tomography angiography,” Biomed. Opt. Express 5, 3822–3832 (2014). [CrossRef]   [PubMed]  

28. G. Liu, A. J. Lin, B. J. Tromberg, and Z. Chen, “A comparison of Doppler optical coherence tomography methods,” Biomed. Opt. Express 3, 2669–2680 (2012). [CrossRef]   [PubMed]  

29. S. Makita, F. Jaillon, I. Jahan, and Y. Yasuno, “Noise statistics of phase-resolved optical coherence tomography imaging: single-and dual-beam-scan Doppler optical coherence tomography,” Opt. Express 22, 4830–4848 (2014). [CrossRef]   [PubMed]  

30. S. Makita, F. Jaillon, M. Yamanari, M. Miura, and Y. Yasuno, “Comprehensive in vivo micro-vascular imaging of the human eye by dual-beam-scan Doppler optical coherence angiography,” Opt. Express 19, 1271–1283 (2011). [CrossRef]   [PubMed]  

31. B. H. Park, M. C. Pierce, B. Cense, and J. F. de Boer, “Jones matrix analysis for a polarization-sensitive optical coherence tomography system using fiber-optic components,” Opt. Lett. 29, 2512–2514 (2004). [CrossRef]   [PubMed]  

32. M. Yamanari, S. Makita, V. Madjarova, T. Yatagai, and Y. Yasuno, “Fiber-based polarization-sensitive Fourier domain optical coherence tomography using B-scan-oriented polarization modulation method,” Opt. Express 14, 6502–6515 (2006). [CrossRef]   [PubMed]  

33. M. Yamanari, S. Makita, and Y. Yasuno, “Polarization-sensitive swept-source optical coherence tomography with continuous source polarization modulation,” Opt. Express 16, 5892–5906 (2008). [CrossRef]   [PubMed]  

34. M. Yamanari, S. Makita, Y. Lim, and Y. Yasuno, “Full-range polarization-sensitive swept-source optical coherence tomography by simultaneous transversal and spectral modulation,” Opt. Express 18, 13964–13980 (2010). [CrossRef]   [PubMed]  

35. B. Baumann, W. Choi, B. Potsaid, D. Huang, J. S. Duker, and J. G. Fujimoto, “Swept source / Fourier domain polarization sensitive optical coherence tomography with a passive polarization delay unit,” Opt. Express 20, 10229–10241 (2012). [CrossRef]   [PubMed]  

36. Z. Wang, H.-C. Lee, O. O. Ahsen, B. Lee, W. Choi, B. Potsaid, J. Liu, V. Jayaraman, A. Cable, M. F. Kraus, K. Liang, J. Hornegger, and J. G. Fujimoto, “Depth-encoded all-fiber swept source polarization sensitive OCT,” Biomed. Opt. Express 5, 2931–2949 (2014). [CrossRef]   [PubMed]  

37. B. Braaf, K. A. Vermeer, M. de Groot, K. V. Vienola, and J. F. de Boer, “Fiber-based polarization-sensitive OCT of the human retina with correction of system polarization distortions,” Biomed. Opt. Express 5, 2736–2758 (2014). [CrossRef]   [PubMed]  

38. R. J. A. Tough, D. Blacknell, and S. Quegan, “A statistical description of polarimetric and interferometric synthetic aperture radar data,” Proc. R. Soc. Lond. A 449, 567–589 (1995). [CrossRef]  

39. R. Touzi, A. Lopes, J. Bruniquel, and P. Vachon, “Coherence estimation for SAR imagery,” IEEE Trans. Geosci. Remote Sens. 37, 135–149 (1999). [CrossRef]  

40. K. Kurokawa, S. Makita, Y.-J. Hong, and Y. Yasuno, “Two-dimensional micro-displacement measurement for laser coagulation using optical coherence tomography,” Biomed. Opt. Express 6, 170–190 (2015). [CrossRef]   [PubMed]  

41. V. X. Yang, M. L. Gordon, A. Mok, Y. Zhao, Z. Chen, R. S. Cobbold, B. C. Wilson, and I. A. Vitkin, “Improved phase-resolved optical Doppler tomography using Kasai velocity estimator and histogram segmentation,” Opt. Commun. 208, 209–214 (2002). [CrossRef]  

42. T. Schmoll, C. Kolbitsch, and R. A. Leitgeb, “Ultra-high-speed volumetric tomography of human retinal blood flow,” Opt. Express 17, 4166–4176 (2009). [CrossRef]   [PubMed]  

43. B. Baumann, B. Potsaid, M. F. Kraus, J. J. Liu, D. Huang, J. Hornegger, A. E. Cable, J. S. Duker, and J. G. Fujimoto, “Total retinal blood flow measurement with ultrahigh speed swept source/Fourier domain OCT,” Biomed. Opt. Express 2, 1539–1552 (2011). [CrossRef]   [PubMed]  

44. B. R. White, M. C. Pierce, N. Nassif, B. Cense, B. H. Park, G. J. Tearney, B. E. Bouma, T. C. Chen, and J. F. d. Boer, “In vivo dynamic human retinal blood flow imaging using ultra-high-speed spectral domain optical Doppler tomography,” Opt. Express 11, 3490–3497 (2003). [CrossRef]   [PubMed]  

45. L. An and R. K. Wang, “In vivo volumetric imaging of vascular perfusion within human retina and choroids with optical micro-angiography,” Opt. Express 16, 11438–11452 (2008). [CrossRef]   [PubMed]  

46. K. Kurokawa, K. Sasaki, S. Makita, Y.-J. Hong, and Y. Yasuno, “Three-dimensional retinal and choroidal capillary imaging by power Doppler optical coherence angiography with adaptive optics,” Opt. Express 20, 22796–22812 (2012). [CrossRef]   [PubMed]  

47. B. Vakoc, S. Yun, J. d. Boer, G. Tearney, and B. Bouma, , “Phase-resolved optical frequency domain imaging,” Opt. Express 13, 5483–5493 (2005). [CrossRef]   [PubMed]  

48. A.-H. Dhalla, D. Nankivil, and J. A. Izatt, “Complex conjugate resolved heterodyne swept source optical coherence tomography using coherence revival,” Biomed. Opt. Express 3, 633 (2012). [CrossRef]   [PubMed]  

49. M. J. Ju, Y.-J. Hong, S. Makita, Y. Lim, K. Kurokawa, L. Duan, M. Miura, S. Tang, and Y. Yasuno, “Advanced multi-contrast Jones matrix optical coherence tomography for Doppler and polarization sensitive imaging,” Opt. Express 21, 19412–19436 (2013). [CrossRef]   [PubMed]  

50. S. Jiao and L. V. Wang, “Jones-matrix imaging of biological tissues with quadruple-channel optical coherence tomography,” J. Biomed. Opt. 7, 350–358 (2002). [CrossRef]   [PubMed]  

51. Y. Lim, Y.-J. Hong, L. Duan, M. Yamanari, and Y. Yasuno, “Passive component based multifunctional Jones matrix swept source optical coherence tomography for Doppler and polarization imaging,” Opt. Lett. 37, 1958–1960 (2012). [CrossRef]   [PubMed]  

52. E. Götzinger, M. Pircher, W. Geitzenauer, C. Ahlers, B. Baumann, S. Michels, U. Schmidt-Erfurth, and C. K. Hitzenberger, “Retinal pigment epithelium segmentation by polarization sensitive optical coherence tomography,” Opt. Express 16, 16410–16422 (2008). [CrossRef]   [PubMed]  

53. S. Makita, Y.-J. Hong, M. Miura, and Y. Yasuno, “Degree of polarization uniformity with high noise immunity using polarization-sensitive optical coherence tomography,” Opt. Lett. 39, 6783–6786 (2014). [CrossRef]   [PubMed]  

54. Commission Internationale de L’Eclairage, Colorimetry, no. 15 in Technical report / CIE (CIE Central Bureau, Vienna, 2004), 3 ed.

55. International Organization for Standardization, Image technology colour management – Architecture, profile format and data structure – Part 1: Based on ICC.1:2010 ISO 15076-1:2010 (International Organization for Standardization, 2010), 2nd eds.

56. International Electrotechnical Commission, Multimedia systems and equipment - Colour measurement and management - Part 2-1: Colour management - Default RGB colour space - sRGB, IEC 61966-2-1 (International Electrotechnical Commission, 1999).

57. M. Geissbuehler and T. Lasser, “How to display data by color schemes compatible with red-green color perception deficiencies,” Opt. Express 21, 9862–9874 (2013). [CrossRef]  

58. M. Abramoff, M. Garvin, and M. Sonka, “Retinal imaging and image analysis,” IEEE Rev. Biomed. Eng. 3, 169–208 (2010). [CrossRef]  

59. K. Li, X. Wu, D. Chen, and M. Sonka, “Optimal surface segmentation in volumetric images-a graph-theoretic approach,” IEEE Trans. Pattern Anal. Mach. Intell. 28, 119–134 (2006). [CrossRef]   [PubMed]  

60. M. Garvin, M. Abramoff, X. Wu, S. Russell, T. Burns, and M. Sonka, “Automated 3-D intraretinal layer segmentation of macular spectral-domain optical coherence tomography images,” IEEE Trans. Med. Imaging 28, 1436–1447 (2009). [CrossRef]   [PubMed]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (11)

Fig. 1
Fig. 1 Schematic block diagrams of cmOCA processing. (a) Overall processing. (b) Covariance matrix estimate with bulk-phase-offset correction.
Fig. 2
Fig. 2 Schematic diagram of applied scanning protocol.
Fig. 3
Fig. 3 Color change trajectory in L*a*b* color space according to cumulative polarization variance cPV.
Fig. 4
Fig. 4 Cross sections of the human retina. (a) Coherence composite OCT intensity, (b) SNR-dependent cmOCA, (c) SNR-independent cmOCA. The scale bar indicates 200 μm.
Fig. 5
Fig. 5 Complex correlation mapping image of the human retina. cmOCA (b, d) with multiple polarization channels, as described in Section 4.2, exhibits lower spurious noise than (a, c) single-channel estimation with a coherent composite signal. The cross sections (a, b) and en-face slices (c, d) are taken from the locations indicated by the black arrows. Both cmOCAs used the same kernel size.
Fig. 6
Fig. 6 cmOCA cross sections without (a, b) and with (c, d) bulk-motion artifact removal. Bulk-phase-offset correction has been applied using the complex-mean method (a, c) and correlation-based method (b, d). The corresponding en-face projection images are shown in (e)–(h).
Fig. 7
Fig. 7 cmOCA images with low-phase-stability SS-OCT system. The images were obtained using (a, d) the multichannel cmOCA algorithm (Section 4.2), (b, e) high phase-noise-immunity mode [Eq. (44)], and (c, f) high phase-noise-immunity mode with bulk-motion-artifact removal [Eq. (45)]. Cross sections (a, b, c) at the location indicated by the black arrow and en-face projections (d, e, f) are shown.
Fig. 8
Fig. 8 cmOCAs at the human fovea, (a) cross-sectional cmOCA image with segmented results, en-face projection images at (b) NFL + GCL, (c) IPL + INL, (d) OPL + ONL. The scanning area is 3 × 3 mm2. The scale bar indicates 200 μm.
Fig. 9
Fig. 9 cmOCA images of polypoidal choroidal vasculopathy. En-face (a) OCT, (b) composition of cmOCA and mDOPU projection images. The scanning area is 6 × 6 mm2. The corresponding area is indicated by the yellow rectangle in fluorescein angiography (c). The cross-sectional (d) OCT + cmOCA and (e) mDOPU images at the center of the scan [yellow line in (b)] are shown. The scale bar indicates 500 μm.
Fig. 10
Fig. 10 The numerical simulation results of correlation coefficient estimates. (a, b, c) Mean and (d, e, f) standard deviation of each correlation estimator, i.e., data correlation (×), signal correlation with previous method using coherent composite (▲), and the new method (●) are plotted with several SNR. The set correlation coefficients are (a, d) 0.99, (b, e) 0.9, and (c, f) 0.8.
Fig. 11
Fig. 11 Comparison of statistical properties at no signal region. Histograms of correlation estimates (a) without the sign ε, (b) without the sign ε and thresholding at 0 instead of absolute operation, and (c) with the sign ε Eq. (16) are shown. The numerical simulation of correlation estimates was applied at SNR=0.

Tables (4)

Tables Icon

Table 1 Variance of correlation coefficient estimation σ r G.

Tables Icon

Table 2 Correlation coefficient ρG for each tissue type and SNR.

Tables Icon

Table 3 Assignment of averaging and variables for each imaging mode.

Tables Icon

Table 4 Assignment of averaging and variables of phase-noise-immune mode.

Equations (45)

Equations on this page are rendered with MathJax. Learn more.

Σ G ( Δ r ; r ) = [ E [ | G ( r ) | 2 ] E [ G ( r ) G * ( r + Δ r ) ] E [ G ( r + Δ r ) G * ( r ) ] E [ | G ( r + Δ r ) | 2 ] ] ,
ρ G ( Δ r ; r ) e i θ G ( Δ r ; r ) = E [ G ( r ) G * ( r + Δ r ) ] E [ | G ( r ) | 2 ] E [ | G ( r + Δ r ) | 2 ]
S G ( Δ r ; r ) g ( Δ r ; r ) g ( Δ r ; r ) w 1 = [ | g ( r ) | 2 g ( r ) g * ( r + Δ r ) g * ( r ) g ( r + Δ r ) | g ( r + Δ r ) | 2 ] w 1 ,
g ( Δ r ; r ) = [ g ( r ) g ( r + Δ r ) ] ,
r G ( Δ r ; r ) ρ ^ G ( Δ r ; r ) | s G 12 | s G 11 s G 22 ,
G ( r ) = S ( r ) + N .
Σ G ( Δ r ; r ) = [ E [ | S ( r ) | 2 ] + E [ | N | 2 ] E [ S * ( r ) S ( r + Δ r ) ] E [ S * ( r + Δ r ) S ( r ) ] E [ | S ( r + Δ r ) | 2 ] + E [ | N | 2 ] . ]
ρ G ( Δ r ; r ) = ρ SNR ( Δ r ; r ) ρ S ( Δ r ; r ) ,
ρ S ( Δ r ; r ) = E [ S * ( r ) S ( r + Δ r ) ] E [ | S ( r ) | 2 ] E [ | S ( r + Δ r ) | 2 ]
ρ SNR ( Δ r ; r ) = [ 1 + SNR ( r ) 1 1 + SNR ( r + Δ r ) 1 ] 1 .
b r G = E [ r G ] ρ S ,
σ r G = E [ r G 2 ] E [ r G ] 2 .
r S ( Δ r ; r ) | s G 12 | | [ s G 11 q ^ N ] [ s G 22 q ^ N ] | ,
q ^ N | g n ( t ) g n ( t ) t | 2 t ,
r ¯ G ( Δ r ; r ) | s G 12 | w 2 s G 11 s G 22 w 2 .
r ¯ S ( Δ r ; r ) | s G 12 | w 2 ε | [ s G 11 q ^ N ] [ s G 22 q ^ N ] | w 2 ,
ε = sgn ( s G 11 q ^ N ) + sgn ( s G 22 q ^ N ) 2 .
r ¯ S ( τ ; x , z ) | g ( r ) g * ( r + ( 0 , 0 , τ ) ) w 1 | w 2 ε | [ | g ( r ) | 2 w 1 q ^ N ] [ | g ( r + ( 0 , 0 , τ ) ) | 2 w 1 q ^ N ] | w 2 .
S G ( τ ; r ) g ( τ ; r ) g ( τ ; r ) w 1 .
r ¯ G ( τ ; r ) | s G 12 | w 2 s G 11 s G 22 w 2 .
Δ ϕ m ( τ ; x , t ) arg { s G 12 ( τ ; r | z = z max ( τ ; x , t ) ) w 3 } .
S G ( τ ; x , z ) [ s G 11 s G 12 e i Δ ϕ m ( τ ; x , t ) { s G 12 e i Δ ϕ m ( τ ; x , t ) } * s G 22 ] w 1 .
r ¯ S ( τ ; x , z ) | g ( r ) g * ( r + ( 0 , 0 , τ ) ) w 1 e i Δ ϕ m ( τ ; x , τ ) w 1 | w 2 ε | [ | g ( r ) | 2 w 1 w 1 q ^ N ] [ | g ( r + ( 0 , 0 , τ ) ) | 2 w 1 w 1 q ^ N ] | w 2 .
S G ( τ ; r ) g ( τ ; r ) g ( τ ; r ) w 1 .
ρ S ( τ ; r ) ρ f ( τ ; r ) ρ m ( τ ; x , t ) ,
r ¯ S ( τ ; x , z ) r ˜ f ( τ ; x , z ) r m ( τ ; x ) ,
r ¯ G ( τ ; x , z ) | s G 12 | w 2 s G 11 s G 22 w 2 .
r ^ m ( τ ; x ) r ¯ S ( τ ; x , arg max Z r ¯ G ( τ ; x , z ) ) .
r ¯ S f ( τ ; x , z ) r ¯ S ( τ ; x , z ) r ^ m ( τ ; x ) .
S G M ( τ ; x , z , f , p ) 1 D d D g ( τ ; x , z + d , f , p ) g ( τ ; x , z + d , f , p )
r ¯ G M ( τ ; x , z , f ) p P 1 q ^ N ( p ) | s G 12 M | p P 1 q ^ N ( p ) s G 11 M s G 22 M ,
Δ ϕ m M ( τ ; x , f ) arg { Med x [ 1 P K p P k = 1 K s G 12 M ( τ ; x , z k max ( τ ; x , f ) , f , p ) ] } ,
r ¯ S M ( τ ; x , z ) p 1 q ^ N ( p ) | l , d , f L , D , F 1 g ( x + l , z + d , f , p ) g * ( x + l , z + d , f + 1 , p ) e i Δ ϕ m M ( τ ; x + l , f ) | p ε ( x , z , p ) q ^ N ( p ) { | [ l , d , f L , D , F 1 | g ( x + l , z + d , f , p ) | 2 L D ( F 1 ) q ^ N ( p ) ] × [ l , d , f L , D , F 1 | g ( x + l , z + d , f + 1 , p ) | 2 L D ( F 1 ) q ^ N ( p ) ] | } 1 2 ,
r ¯ S f M ( τ ; x , z ) r ¯ S M ( τ ; x , z ) r ^ m M ( τ ; x ) .
r S C ( τ ; x , z ) | l , d , f L , D , F 1 g C ( x + 1 , z + d , f ) g C * ( x + 1 , z + d , f + 1 ) e i Δ ϕ m C ( τ ; x + 1 , f ) | ε ( x , z ) { | [ l , d , f L , D , F 1 | g C ( x + 1 , z + d , f ) | 2 L D ( F 1 ) P 2 p P q ^ N ( p ) ] × [ l , d , f L , D , F 1 | g C ( x + 1 , z + d , f + 1 ) | 2 L D ( F 1 ) P 2 p P q ^ N ( p ) ] | } 1 2 ,
r S f C ( τ ; x , z ) r S C ( τ ; x , z ) r ^ m C ( τ ; x ) .
M [ x ] = { 1 , if x 1 < 1. x , otherwise .
E ( x ) = n = 0 N j c j ( x , z n ) m s ( j ) Δ z [ e i = 0 n j c j ( x , z i ) m a ( j ) Δ z L 0 ] ,
E ( x ) = n = 0 N { 1 M [ r ¯ S f ( τ ; x , z n ) ] } T [ W [ e c P V n ( x ) m Δ z L 0 ] ] ,
m = [ 1.167 0 0 83.067 0 0.236 0 27.944 0 0 0.283 1.932 0 0 0 0 ] [ μ m 1 ] .
e c P V n ( x ) m Δ z = V diag ( e c P V n ( x ) λ 1 Δ z , e c P V n ( x ) λ 2 Δ z , e c P V n ( x ) λ 3 Δ z , e c P V n ( x ) λ 4 Δ z ) V 1 .
W [ L in ] = { W L in ( v L in > 0 ) L in ( v L in 0 ) ,
W = [ 1 0 0 0 0 cos ( 2 θ ) sin ( 2 θ ) a o + a o cos ( 2 θ ) + b o sin ( 2 θ ) 0 sin ( 2 θ ) cos ( 2 θ ) b o b o cos ( 2 θ ) + a o sin ( 2 θ ) 0 0 0 1 ] .
r ¯ S ( τ ; x , z ) l , f L , F 1 | 1 D p , d P , D g ( x + l , z + d , f , p ) g * ( x + l , z + d , f + 1 , p ) | l , f L , F 1 ε ( τ ; x + l , z , f ) { | [ 1 D p , d P , D | g ( x + l , z + d , f , p ) | 2 p P q ^ N ( p ) ] × [ 1 D p , d P , D | g ( x + l , z + d , f + 1 , p ) | 2 p P q ^ N ( p ) ] | } 1 2 .
r ¯ S f ( τ ; x , z ) r ¯ S ( τ ; x , z ) r ^ m ( τ ; x , ) .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.