Improved motion contrast and processing efficiency in OCT angiography using complex-correlation algorithm

The imaging system used in this study was built based on the configuration of spectral domain OCT (SDOCT) as shown in figure 1. A broadband superluminescent diode (Superlum, Carrigtwohill, Ireland, Broadlighters D855-HP2) was used as the light source 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 output light was delivered into a 2 × 2 fiber coupler and then split into the reference arm and the sample arm. In the sample arm, an X–Y galvanometer was adopted for scanning a 3D volume, and an objective lens with the focal length of 40 mm was used to focus the probing light beam on the region of interest, yielding a measured lateral resolution of ∼15 μm. The OCT detection unit in our system was a high-speed spectrometer, equipped with a fast line-scan CMOS camera (Basler, Ahrensburg, Germany, Sprint spL4096-140k) providing a 120 kHz line-scan rate and 2048 active pixels. The spectrometer had a designed spectral resolution of 0.062 nm, providing an imaging range of ∼2.9 mm on each side of the zero delay line in air. The data processing unit in our system was a Dell Precision T5610 workstation installed with a single Intel Xeon Processor E5-2609 v2 and 16G RAM.

Zoom In Zoom Out Reset image size

As reported in [26], a general expression for the sensitivity of SDOCT is ${\rm{SNR}}=\eta {P}_{{\rm{s}}}\tau /{E}_{\upsilon },$ where $\eta$ is the quantum efficiency (∼38% in our system) of the detector, ${P}_{{\rm{s}}}$ is the sample arm power returning to the detection unit, $\tau$ is the integration time of detector (∼$8.4\;\mu {\rm{s}}$ in our system), and ${E}_{\nu }=hc/{\lambda }_{{\rm{c}}}$ (${\lambda }_{{\rm{c}}}$ is central wavelength of the light source, i.e. 850 nm in our system) is the photon energy. With an incident optical power of 2 mW on the perfect reflector sample surface, a theoretic signal-to-noise ratio (SNR) of ∼104 dB is predicted. The theoretic predication is in well agreement with the experiment measurement. As shown in figure 2, a system sensitivity of 100 dB was measured for the peak position at 0.5 mm with a calibrated sample arm attenuation of −52 dB. The experimental SNR in dB was calculated as twenty times the base-10 logarithm of the ratio of the A-scan peak amplitude to the standard deviation of the noise floor which was taken at the location of the A-scan peak by blocking the sample arm.

Zoom In Zoom Out Reset image size

In this study, a step-scanning protocol with repeated B-scans was used to acquire 3D datasets. The volumetric scanning was composed of 512 A-scans in the fast-scan (X) direction and 200 steps in the slow-scan (Y) direction. Each 512 A-scans form one B-frame and five repeated B-frames were acquired at each step, resulting in a total of 1000 B-frames for each 3D dataset. The scanning range covers a 2 × 2 mm spatial distance in the X and Y directions. The B-frame rate was 190 frames s⁻¹, determining a total acquisition time of ∼5.3 s for each volumetric dataset. The repeated scanning is designed to suppress the decorrelation noise in blood flow imaging [27].

The depth-resolved complex analytic OCT spatial signal is reconstructed by performing Fourier transform of the raw spectral interference fringe signal. At each step of scanning, the complex analytic signal of the nth repeated B-frame is denoted as ${\tilde{A}}_{n}(z,x)={A}_{n}(z,x)\mathrm{exp}[{\rm{i}}{\phi }_{n}(z,x)],$ where, n is the B-frame index, and the coordinates ($z,x)$ are indices along the Z and X directions.

2.2.1. Parallel processing of BIS correction and flow signal extraction

Large BTMs, particularly the respiratory motion, may result in BISs between adjacent B-frames, which consequently cause a global decorrelation and hinder the identification of the dynamic flow regions. In this study, we follow the assumptions that the BISs are rigid with no deformation and rotation, and out-of-plane motion is also limited for well pre-process and fixation of the sample.

The processing flowchart of CC-Angio-OCT algorithm is depicted in figure 3. At each scanning step, adjacent complex-valued B-frames ${\tilde{A}}_{n}(z,x)\;$ and ${\tilde{A}}_{n+1}(z,x)$ are paired. The paired B-frames are up-sampled using cubic spline interpolation in both Z (axial) and X (lateral) directions, and then analyzed by the normalized cross-correlation algorithm

$$\begin{eqnarray}&&{{\rm{CC}}}_{n}(z,x)\\ &&=1-\displaystyle \frac{|\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{\tilde{A}}_{n}(z+p,x+q){\tilde{A}}_{n+1}^{*}(z+p,x+q)|}{\sqrt{\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{[{A}_{n}(z+p,x+q)]}^{2}}\sqrt{\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{[{A}_{n+1}(z+p,x+q)]}^{2}}},\end{eqnarray} \tag{ 1 }$$

where, * means the complex conjugate. The frame index n ranges from 1 to 4 in this work. P and Q are the window sizes in the Z and X directions, respectively, and p and q are the corresponding indices of the pixels within the window. Moving the 2D window across the entire frame, a cross-sectional CC map is generated. In the CC map, the value of each pixel, as an index of decorrelation, lies in the range [0, 1].

Zoom In Zoom Out Reset image size

BISs lead to a global decorrelation between adjacent B-frames. In order to eliminate the effect of BISs, the mean value $\bar{{{\rm{CC}}}_{n}(z,x)}$ of all the pixels is treated as an index of the global decorrelation and sub-pixel alignment can be achieved by minimizing the mean value $\bar{{{\rm{CC}}}_{n}(z,x)}.$ In the aligned CC map, the areas with high decorrelation can be classified as the dynamic flow regions. Then the aligned CC map is down-sampled to its original sampling to generate a cross-sectional CC angiogram. The up-sampling performed here is to achieve an alignment precision of sub-pixel and the followed down-sampling is to restore the original size of dataset. In the proposed CC algorithm, the flow signal extraction is completed in parallel with the BIS correction.

Additionally, in cross-sectional CC angiograms, the static tissue background with weak OCT signals also exhibit high decorrelation due to the large influence of noise, which may be misclassified as dynamic blood flow. To remove the noise-induced decorrelation artifacts, a 'structural mask', which is generated by performing a binary threshold on the OCT amplitude signal ${A}_{n}(z,x),$ is applied to the CC map. The threshold value is empirically chosen to be above the mean background value.

With the trade-off among SNR, resolution and computation time, the up-sampling factor in this study takes the value four and the window size P and Q in equation (1) take the value 12. In each scanning step, five repeated B-scans are performed, producing a total of four cross-sectional angiograms. The four angiograms are averaged to improve SNR. A total of 200 cross-sectional angiograms are obtained to reconstruct a 3D map of vascular network.

2.2.2. Immunity against GPFs

Equation (1) can be further expressed as

$$\begin{eqnarray}&&{{\rm{CC}}}_{n}(z,x)\\ &&=1-\displaystyle \frac{|\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{A}_{n}(z+p,x+q){A}_{n+1}(z+p,x+q)\mathrm{exp}[j{\rm{\Delta }}{\phi }_{n}(z+p,x+q)]|}{\sqrt{\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{[{A}_{n}(z+p,x+q)]}^{2}}\sqrt{\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{[{A}_{n+1}(z+p,x+q)]}^{2}}},\end{eqnarray} \tag{ 2 }$$

where, ${\rm{\Delta }}{\phi }_{n}(z,x)={\phi }_{n+1}(z,x)-{\phi }_{n}(z,x)$ is the phase difference between the paired B-frames. Small BTMs and random environmental vibrations result in GPFs. Thus, ${\rm{\Delta }}{\phi }_{n}(z,x)$ can be written as

$$\begin{eqnarray}&&{\rm{\Delta }}{\phi }_{n}(z,x)={\rm{\Delta }}{\phi }_{n}^{{\rm{F}}}(z,x)+{\rm{\Delta }}{\phi }_{n}^{{\rm{GPF}}}(z,x),\end{eqnarray} \tag{ 3 }$$

where, ${\rm{\Delta }}{\phi }_{n}^{{\rm{F}}}(z,x)={\phi }_{n+1}^{{\rm{F}}}(z,x)-{\phi }_{n}^{{\rm{F}}}(z,x)$ is the flow-induced phase change. ${\rm{\Delta }}{\phi }_{n}^{{\rm{GPF}}}(z,x)$ is the GPF-induced phase change. Following the assumptions used in [21], ${\rm{\Delta }}{\phi }_{n}^{{\rm{GPF}}}(z,x)$ is considered as a constant within a limited 2D window, of which the size is $P\times Q$ pixels in equation (2). Substituting equation (3) into equation (2), we have

$$\begin{eqnarray}&&{{\rm{CC}}}_{n}(z,x)\\ &&=1-\displaystyle \frac{|\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{A}_{n}(z+p,x+q){A}_{n+1}(z+p,x+q)\mathrm{exp}[j{\rm{\Delta }}{\phi }_{n}^{{\rm{F}}}(z+p,x+q)]|}{\sqrt{\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{[{A}_{n}(z+p,x+q)]}^{2}}\sqrt{\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{[{A}_{n+1}(z+p,x+q)]}^{2}}}.\end{eqnarray} \tag{ 4 }$$

From the expression above, we can know that the CC mapping is insensitive to the GPFs.

In this study, CD-Angio-OCT algorithm is also conducted for comparison,

$$\begin{eqnarray}&&{{\rm{CD}}}_{n}(z,x)=|{\tilde{A}}_{n+1}(z,x)\mathrm{exp}[-{\rm{i}}{\rm{\Delta }}{\phi }_{n}^{{\rm{GPF}}}(z,x)]-{\tilde{A}}_{n}(z,\;x)|.\end{eqnarray} \tag{ 5 }$$

Different from the CC method, the CD is sensitive to GPFs. The phase term ${\rm{\Delta }}{\phi }_{n}^{{\rm{GPF}}}(z,x)$ is used to compensate the GPF, which is determined by a histogram-based phase selecting algorithm [3, 17, 25, 28, 29]. The phase compensation is conducted A-line by A-line, and consequently the procedure is extremely time-consuming. Exactly, in our system for each 3D volumetric dataset (i.e. 2048 × 512 × 1000 pixels), the processing time of the CD algorithm is ∼4455 s, while that of the CC algorithm is ∼401 s which is about 10 times faster.

2.2.3. Superior motion contrast

The dynamic scattering of the moving blood cells causes changes in both the amplitude and phase of OCT signals. In the conventional IC method, only the amplitude information ${A}_{n}(z,x)$ is utilized for extracting the flow signal, which can be expressed by

$$\begin{eqnarray}&&{{\rm{IC}}}_{n}(z,x)\\ &&=1-\displaystyle \frac{\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{A}_{n}(z+p,\;\;x+q){A}_{n+1}(z+p,x+q)}{\begin{array}{l}\sqrt{\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{[{A}_{n}(z+p,x+q)]}^{2}}\\ \sqrt{\displaystyle {\sum }_{p=0}^{P-1}\displaystyle {\sum }_{q=0}^{Q-1}{[{A}_{n+1}(z+p,x+q)]}^{2}}\end{array}}.\end{eqnarray} \tag{ 6 }$$

Comparing the CC and IC algorithms (i.e. equations (1) and (6)), the only difference is the numerator of the equations. For clarity, ${A}_{n}(z+p,x+q){A}_{n+1}(z+p,x+q)$ $\mathrm{exp}[j{\rm{\Delta }}{\phi }_{n}^{{\rm{F}}}(z+p,x+q)]$ is denoted by a vector $M(p,q).$ As illustrated in figure 4, due to the influence of phase, the length of the vector sum $|{\sum }^{}M|$ is not larger than the sum of the vector length ${\sum }^{}|M|.$ Then, we have

$$\begin{eqnarray}{{\rm{CC}}}_{n}(z,x)={{\rm{IC}}}_{n}(z,x) & {\rm{static}}\;\mathrm{region},\\ {{\rm{CC}}}_{n}(z,x)\gt {{\rm{IC}}}_{n}(z,x) & {\rm{dynamic}}\;{\rm{region}}{\rm{.}}\end{eqnarray} \tag{ 7 }$$

In the static regions, the flow-induced phase change ${\rm{\Delta }}{\phi }_{n}^{{\rm{F}}}$ approximates to zero in theory, while in the dynamic regions, the phase change ${\rm{\Delta }}{\phi }_{n}^{{\rm{F}}}$ is not constant and exhibits spatial fluctuations (e.g. parabolic distribution) within the lumen in reality. According to equation (7), it is reasonable to conclude that, by combining both the amplitude and phase information, the CC-Angio-OCT offers a higher motion contrast than the IC method does.

To validate the feasibility and performance of the proposed CC-Angio-OCT algorithm, both phantom and in vivo animal experiment were conducted.

Zoom In Zoom Out Reset image size

The flow phantom was made of an agarose gel mixed with ∼5% milk to mimic the static scattering tissue background and a capillary tube with an inner diameter of 0.5 mm was embedded in this tissue-like phantom. Three percent milk solution was pumped into the tube at a constant rate by syringe pump (KDS 100 series, Stoelting Co., Wood Dale, Illinois) to simulate the flowing blood. The phantom was placed on a translation stage, and BISs in the X and Z directions were manually introduced during OCT scanning.

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

Figure 5 shows the results of the BIS correction and flow signal extraction in phantom experiments. Figure 5(a) is a representative cross-sectional OCT structural image of the phantom. The transparent tube can be visualized clearly. The regions inside and outside the tube represent the dynamic fluid and the static solid gel, respectively. Figure 5(b) shows the corresponding cross-sectional CC angiogram without BIS correction. Due to the BISs, both the static and dynamic regions exhibit high decorrection, and cannot be distinguished from each other. Using the sub-pixel alignment, sub-pixel shifts (3/4 pixel in the Z direction and 1/2 pixel in the X direction) are determined. The motion contrast appears in the aligned CC map without any further processing. Figure 5(c) presents the CC angiogram with the BIS correction conducted only in the Z direction, and figure 5(d) reports the result with the correction conducted in both Z and X directions.

Zoom In Zoom Out Reset image size

With the BISs corrected, the immunity against GPFs of CC method was then tested in comparison with CD method. The left column in figure 6 presents the CD (a) and CC (b) angiograms without phase compensation, while the right one is the corresponding flow images after phase compensation. Comparing figure 6(b) with 6(d), it can be observed that the CD method is highly sensitive to GPFs. Due to the influence of GPFs, no motion contrast can be visualized before phase compensation in the CD method. In comparison, the CC method exhibits a high immunity against GPFs, referring to figures 6(a) and (c).

Zoom In Zoom Out Reset image size

Figure 7 presents a comparison of the motion contrast between the IC and CC methods. Figures 7(a) and (b) are the IC and CC angiograms, respectively. They are displayed in the same dynamic range and a motion-contrast improvement can be visualized in the CC flow image. The quantitative comparison is provided in figures 7(c) and (d). Figure 7(c) shows the representative signal profiles of the IC and CC methods, which are extracted from the same depth position indicated by the yellow dashed lines in figures 7(a) and (b), respectively. By comparison, although the CC profile exhibits a little higher background in the static part, the signals in flow part seem much more powerful. By calculating the mean value of the flow ($\bar{{\rm{Flow}}})$ and solid gel ($\bar{{\rm{Tissue}}})$ regions respectively, the motion contrast, which is defined as the ratio of $\bar{{\rm{Flow}}}$ and $\bar{{\rm{Tissue}}},$ can be quantified. Based on the statistical analysis of 200 flow images, the motion contrast of the CC method is much higher than that of the IC method (10.41 versus 6.64, ∼1.5 times) as reported in figure 7(d).

Zoom In Zoom Out Reset image size

Figure 8(a) is a representative cross-section of the OCT structural image, showing the depth-resolved morphological features of the mouse cerebral cortex, such as cranium and cortex. Figure 8(b) is the corresponding cross-sectional CC angiogram without phase compensation, in which the depth-resolved vasculature can be clearly visualized, as indicated by the arrows. Figures 8(c) and (d) are the cross-sectional CD angiograms before and after GPF compensation, respectively. It can be perceived that the CD angiogram is highly sensitive to the GPFs, and the residual static background is apparent in the cranium region.

Zoom In Zoom Out Reset image size

Figure 9 presents a comparison of the motion-contrast between the IC and CC-Angio-OCT of the mouse cortex in vivo. Figures 9(a) and (b) are the representative cross-sectional angiograms using the IC and CC algorithms, respectively. They are displayed in the same dynamic range. Higher contrast of flow signals can be clearly visualized in the CC method. The signal profiles of the same flow area (indicated by the white arrows in figures 9(a) and (b)) are extracted along the same position (marked by the yellow dashed lines in figures 9(a) and (b)), as plotted in figure 9(c). The blood flow and static tissue regions can be distinguished. The motion contrast, i.e. the ratio of $\bar{{\rm{Flow}}}$ and $\bar{{\rm{Tissue}}},$ can be quantified. Totally, 100 flow areas are likewise selected randomly from 200 cross-sectional angiograms for statistical analysis, as plotted in figure 9(d). The mean values (mean ± standard deviation) of motion contrast are 48.78 ± 5.82 and 73.29 ± 5.54 in IC and CC methods, respectively, demonstrating ∼1.5 times improvement in motion contrast of the proposed CC method.

Zoom In Zoom Out Reset image size

Figure 10 presents a 3D CC angiogram to demonstrate the feasibility of the CC algorithm for in vivo applications. Figure 10(a) is the 3D rendering of vasculature. The corresponding en face maximum intensity projection view is displayed in figure 10(b). Figure 10(c) shows a representative cross-sectional structural image overlaid with the corresponding blood flow signals, which was extracted from the position indicated by the yellow dashed line in figure 10(b).

Zoom In Zoom Out Reset image size