Complex Correlated Phase Gradient Variance Based Optical Coherence Tomography Angiography

In this study, we reported a complex-signal-based optical coherence tomography angiography (OCTA) method, called complex correlated phase gradient variance (CCPGV), for mapping high-quality microvascular images. The performance of the newly proposed algorithm is benchmarked against the previously reported phase-resolved Doppler variance, complex differential variance (CDV), and split-spectrum amplitude and phase-gradient angiography (SSAPGA), by both tissue phantoms experiments and in vivo human skin measurements. Compared to the phase-resolved Doppler variance, CCPGV can intrinsically reject undesirable phase shifts caused by bulk motion and trigger jitter. In contrast to CDV and SSAPGA, which are insensitive to phase instability, CCPGV can provide superior motion contrast, as demonstrated by ∼1.4 and 3 times higher contrast in phantom experiments respectively. An increase of 27.3% and 106.3% were found in vivo experiments for CCPGV, making it easier to distinguish vessels from the background static. Benefiting from the advantages, more vessels and better connectivity can be visualized in the en-face angiogram processed by the CCPGV method. We believe that our method will benefit the biomedical community to some extent in disease diagnosis and monitoring.


I. INTRODUCTION
M ICROVESSELS are the thinnest and most widely distributed blood vessels within the human body and are responsible for many critical physiological functions, including regulation of blood pressure, body temperature, blood flow within tissues, delivery of nutrients, and removal of metabolic waste. Structural and morphological changes of the microvessels can be symptomatic of various pathological conditions, such as diabetes [1], Alzheimer's [2], cancer [3], psoriasis [4], and wound healing [5]. Therefore, detection and imaging of the microvascular map can serve as important evidence for disease diagnosis, risk stratification, disease monitoring, and assessment of treatment effectiveness. Accordingly, this boosts many techniques capable of imaging microvessels. Among them, optical coherence tomography angiography (OCTA) has been rapidly gaining attention as a powerful imaging modality to visualize vascular networks within the tissue (e.g., the retina and skin) due to the advantage of high imaging speed, high resolution, and non-invasiveness in clinical applications [6]. Different from the traditional dye-based angiography techniques (i.e., fluorescein angiography and indocyanine green angiography), OCTA eliminates the risk of injecting contrast agents and minimizes the time consumption. The technique leverages blood flow-induced signal variation as an intrinsic contrast mechanism to distinguish vasculature from static tissues. Intending to contrast blood flow within microcirculatory tissue beds, many data processing methods are available to acquire the angiographic signals [6], [7], [8]. According to the components explored, these methods can be broadly classified into three categories: phase-signal-based OCTA techniques, intensity-signalbased OCTA techniques, and complex-signal-based OCTA techniques. Among them, the complex-signal-based methods have been reported to provide better performance than the pure phaseor amplitude-based methods in vascular imaging [9], [10]. However, as with phase-based methods, complex algorithms are highly vulnerable to bulk motion (including the sample or/and machine movement) [11]. Therefore, complex phase correction or compensation methods, such as the median/mean value searching method [12], [13], or histogram-based maximum bin searching method [14], [15], [16], [17], are typically required for high-quality angiography images. Whereas effective, in regions with a relatively large vessel or strong vessel shadow, these methods will fail to generate the correct bulk phase and may require additional correction steps [18]. An additional main issue is trigger jitter, a source of noise unique to swept-source devices in OCTA technology [19]. It is caused by minor sweep triggering time differences. The slight difference would significantly deteriorate the performance of complex-based OCTA algorithms for high-quality angiography. Therefore, apart from a recent study that used a phase-stabilized akinetic laser source [20], several methods based on software or hardware, for example utilizing an additional Mach-Zehnder interferometer [21], external phase [22], [23], and an extra fiber Bragg grating as a fixed wavenumber [24] have been developed to compensate the trigger-jitter induced phase error. Although feasible, these methods either increase system complexity and cost or increase computational time. Recently, two complex-signal-based angiographic algorithms have been reported to be inherently insensitive to these phase errors induced by bulk motion and trigger jitter. One approach is complex differential variance (CDV) [25], [26], which can minimize sub-pixel bulk phase shift caused by axial motion and phase artifact resulting from trigger jitter. Another approach is split-spectrum amplitude and phase-gradient angiography (SS-APGA) [27], which can eliminate the depth-independent bulk phase shift and minimize the trigger-jitter-induced phase error. Benefiting from the advantages, the two methods have drawn tremendous attention in biomedical research [28], [29], [30], [31], [32], [33], [34]. However, they both provide a relatively low motion contrast. Low motion contrast would result in significant misclassification of dynamic and static regions, severely degrading vascular visibility and hindering the interpretation of imaging outcomes such as hemodynamic quantification [35]. Clinically, an angiographic algorithm that can provide robust phase artifacts rejection and superior motion contrast is preferred.
In this paper, we developed a complex correlated phase gradient variance (CCPGV) method for mapping microvasculature. It can not only strongly reject the phase artifacts but also provide an appealing high motion contrast. The performance is verified through the flow phantom experiment by comparison with the three published algorithms. The feasibility of the proposed angiographic method is also demonstrated for human skin imaging with an advanced Jones-matrix swept-source OCT imaging system. It has been concluded that the use of polarization-based angiography has advantages over conventional OCTA by improving vessel contrast and increasing image signal [36], [37].

A. System
The experimental setup of the Jones-matrix swept-source OCT system has been described previously [38] as shown in Fig. 1. Briefly, light from a MEMS-based swept-source (AXP50125-6, AXSUN Technologies, center wavelength = 1310 nm, sweeping range = 100 nm, and scanning rate = 200 kHz) was split into a reference arm and a sample arm. The light entering the sample arm was multiplexed into two orthogonal polarization states with a passive polarization-delay unit and then sent into the probe unit. In the probe unit, the collimated beam is focused on the sample by an objective lens (AC254-060-C-ML, Thorlabs, the effective focal length = 60 mm, working distance = 46 mm) after a two-axis galvanometer mirror scanner (GVSM002-EC/M, Thorlabs). The incident power to the sample is approximately 10.8 mW. After traveling through the sample and backscattering in the tissue, the light was recollected by the same fiber and interfered with the reference light. The interference is then split into horizontally and vertically polarized light and detected by two separate balanced detectors (BPD, PDB480C-AC, Thorlabs). The combination of the depth multiplexing of the incident beam and the polarization diversity detection can simultaneously provide four OCT signals, which correspond to the four entries of the Jones matrix. After detection, the signal from the BPDs is digitized by a dual-channel digitizer (ATS9371, AlazarTech) with a 12-bit resolution and 1.0 GHz bandwidth. The system utilizes K-trigger mode so that no re-calibration is needed.
The system sensitivity was measured to be 98 dB from a mirror close to the zero-delay location by coherently combining the four polarization channels. The imaging depth range is around 2.5 mm in air, which corresponds to 1.8 mm in skin tissue assuming a refractive index of 1.38. The axial resolution and the transverse resolution are measured to be 10 μm and 19.4 μm in the air, respectively.

B. Method
Fig . 2 shows the data processing steps of the CCPGV. The Jones matrix, describing the polarization characteristics of the sample, can be generated directly by performing a Fourier transform of the raw spectral interference fringe signals. A sensitivityenhanced scattering frame is composited by coherently summing the four entries of the Jones matrix [39] and used to generate the most appealing quality of angiograms. The registration between repeated frames at the same location is performed using the subpixel image registration algorithm. To validate the feasibility and performance of the proposed algorithm, no correction or compensation methods in this paper are adopted to eliminate phase artifacts. In the following analysis, x and z represent Cartesian coordinates (x: fast scanning direction, and z: depth direction).
Typically, the time changes can be calculated from the conjugate multiplication of the A-line pair as a complex correlation, defined as: where C x,z,t 1 indicates the complex value at the time t 1 and C * x,z,t 2 is the conjugate of C x,z,t 2 at the time t 2 . A x,z,t 1 and A x,z,t 2 are the amplitude at the time t 1 and t 2 respectively. Here, the phase difference Δφ x,z,t between two repeated B-scans can be calculated by: Theoretically, (2) only obtains the phase shift induced by red blood cell (RBC) movement inside the vessel. However, it also includes phase noise and phase artifact originating from bulk motion and trigger jitter represented as: where Δφ v is the phase induced by RBC movement, Δφ b is the bulk phase which is typically independent of depth, andΔφ n is the random phase noise, which is typically produced by galvanometer scanners [40]. 2πδ(t) z N , a linear function along the axial direction, indicates the phase artifact caused by the trigger jitter. δ is the relative k shift between acquired interferograms and N is the number of data points in an A-scan. According to the nature of the phase artifacts, the gradient along the depth direction can be employed for the phase difference to separate the phase shift induced by RBC from these phase artifacts [27]. Since the derivative of the phase noise term can be neglected [6], [27], (3) could be simplified as: This is the axial gradient of the phase shift induced by RBC in the blood vessels. Hence, the time changes of signals can be calculated by: whereρ is the weight parameter (in this study, ρ = 4). Further, the variance is calculated along the depth direction to improve the signal-to-noise ratio (SNR) and motion contrast. Thus, the flow signal can be acquired using the following equation: where R indicates the repetition number of B-scans at the same location and w(k) is the depth window function of length 2L + 1 (in this study, L = 5). The denominator is the arithmetic average of the intensities to normalize. We note this angiographic algorithm includes both phase and intensity but the angiographic signal is largely immune to the phase shift caused by bulk motion and trigger jitter. Furthermore, in the static regions, the phase in (6) is close to zero, suggesting that the numerator and denominator are approximately equal in magnitude. As a result, f static → 0. Whereas inflow areas, since the complex vectors are at random, the sum of the vector length in the numerator term is always smaller than the sum of the arithmetic average in the denominator. Hence, f flow → 1. Accordingly, this method can offer a superior motion contrast, making it easy to interpret and analyze the angiography signature.

C. Scanning Protocol
A raster scan with repeated B-scans was conducted to obtain the 3D dataset. The volumetric scanning was composed of 512 A-scans in the fast-scan direction, and 400 steps in the slow-scan direction. To suppress the noise in blood flow imaging, four repeated B-scans were acquired at the same step with a time interval of 3.2 ms. The imaging range covered 3.0 mm × 3.0 mm in the lateral and 2.5 mm in depth. The total data acquisition time was 5.12 seconds.

A. Immunity of Phase Artifacts
As mentioned above, the bulk motion and trigger jitter would significantly deteriorate angiographic image quality, especially images processed by phase-sensitive algorithms. Therefore, complex phase correction or compensation is typically a critical requirement to eliminate phase artifacts. Fortunately, our proposed complex-signal-based CCPGV method is intrinsically strongly immune to these phase artifacts. To verify this performance, we acquired OCT data from human skin on the back of the hand and compared the processed angiograms using CCPGV and phase-resolved Doppler variance [41], [42], another angiographic algorithm based on complex signals but sensitive to phase artifacts. Fig. 3 compares the cross-sectional angiograms generated from the same data acquired with the OCT system. The phase noise mainly due to bulk motion and trigger jitter manifests as the characteristic vertical banding, as shown in Fig. 3(b). Fig. 3(c) demonstrates the Doppler variance image with the axial subpixel phase correction [43]. It shows a dramatic improvement in image quality. However, a small amount of minor artifact still remains in the image. This is because many phase compensation algorithms assume a constant phase across imaging depth, but the random electronic noise introduced by the trigger jitter violates this assumption [19]. In contrast, these artifacts are absent in the angiogram processed by the CCPGV method [ Fig. 3(d)]. The white spots with tailed artifacts are blood vessels in the dermis in Fig. 3(c) and (d), as pointed out by the red arrows.

B. Phantom Experiments
To validate the performance of the proposed CCPGV algorithm in enhancing motion contrast, a two-region flow phantom experiment was performed. The static tissue was mimicked by a piece of Teflon. The flowing blood was simulated by pumping a 3% milk solution into the plastic flow channel at a constant rate using a syringe (KDS 100 series, Stoelting Co., Wood Dale, Illinois). The flow speed was ∼1cm/s, similar to moderate blood flow in real tissue [10]. We compared the processed angiograms using CCPGV and two previously reported angiographic algorithms based on complex signals but insensitive to phase artifacts, namely CDV and SSAPGA. The CDV algorithm is defined as [25]: Here, w(k) is the depth window function of length 2L + 1, and C r (x, z) indicates the complex value in r-th B-scans at the lateral location x and depth position z, and C * r is the conjugate Fig. 4. The relationship between the weight parameter and ISNR for CCPGV.
of C r . The SSAPGA algorithm is given as [27]: where M is the number of narrow split spectrum bands, A m r (x, z) indicates the intensity value in r-th B-scans of m-th split spectrum at the lateral location x and depth position z, ρ is the weight parameter and P G is the gradient of the phase difference between the consecutive B-scans.
We note that the depth summation in the CCPGV and CDV may lead to blurring along the depth axis within the range defined by the kernel. We have explored angiographic performance as a function of the width of w(k) and found that optimal imaging is achieved for the kernel of limited widths (full-width half max of 25 μm in the air). This is closely consistent with the results reported in Ref. [25]. Therefore, the depth window size for CDV and CCPGV methods were both set to 11 pixels. Additionally, to select the best split number and weight parameters for highcontrast OCTA images, the image signal-to-noise (ISNR) [44] is used to evaluate OCTA imaging quality. The ISNR of the flow phantom is defined as follows: where flow and static are the average signal values within the flow and the static regions based on the statistical analysis of 400 phantom images, respectively, and σ static is the intensity standard deviation in the static region, representing the noise from the static background. Fig. 4 shows the relationship between the ρ and ISNR for CCPGV. There is a maximum ISNR when ρ is equal to 4, as indicated by the red circle. Therefore, ρ is set to 4 in this study. Fig. 5 shows the relationship between M , ρ , and the ISNR for SSPAGA. For any given number of spectral splits, the maximum ISNR is always positioned at the point where the weight parameter is equal to 2. Therefore, ρ is set to 2 in this work for SSAPGA. Additionally, the ISNR increases as the number of spectral splits increase when ρ is 2. However, the ISNR growth slowed down significantly when the sub-spectral number is greater than 9. Given the burden of data processing, the number of sub-spectra in this paper is set to 9, as the white circle indicates. Fig. 6 shows the representative cross-sectional OCT structural image and corresponding cross-sectional angiograms. The left and right halves of each image represent the signal from static and moving scatterers, respectively. To avoid specular reflection, the tube is artificially tilted at a certain angle so that it is not perpendicular to the probe beam. In addition, the scanning direction of the fast axis is also not perpendicular to the direction of flow, so the imaging area of the moving scatterers is elliptical. The selected plastic channel had an inner diameter of 2 mm, resulting in weakly backscattered signals from deeper regions inside the tube. To eliminate the effect of the noise on the result, the log OCT reflectance and angiograms are thresholded. The threshold level was determined by three standard deviations above the mean of the estimated noise. Fig. 6(b), (c) shows the corresponding angiograms processed by CDV, SSAPGA, and CCPGV methods on the same dataset, respectively. For a fair comparison, they are all displayed in the same dynamic range. The regions of interest (ROIs) for the static and mobile areas are defined by the rectangular boxes of the red solid line and blue dashed line. By qualitatively comparing, the CCPGV method can provide the most powerful OCTA signals for the moving scatterers and can effectively suppress the OCTA signals for the static scatterers.
A quantitative comparison is also conducted using the index of motion contrast (i.e., MC). It is defined as: where the flow and static are the mean value of the flow and static within the ROIs respectively. The higher the motion contrast, the easier it is to distinguish between dynamic and static areas. The low motion contrast would result in significant misclassification of dynamic and static regions, severely degrading the vascular visibility and hindering the interpretation of imaging outcomes such as hemodynamic quantification. Additionally, we also assessed the algorithms' complexity in terms of computational time to demonstrate their practicability for the dynamic monitoring of blood flow. For a fair comparison,  all coding algorithms were run using MATLAB processing language run on a laptop with an Intel Core i7 processor and 32 GB Memory. Table I tabulates the quantitative comparison among different angiography methods. CCPGV was found to provide the most attractive motion contrast, approximately 3 and 1.4 times that of SSAPGA and CDV, respectively. For SSAPGA, the reduction in the axial resolution after spectral segmentation results in lower sensitivity to detect flow, which is manifested as low motion contrast, as shown in Fig. 6(c). For CDV, the implementation of the square root is intended to increase the angiographic signal, which in turn leads to a higher signal in static areas. In contrast, using the proposed CCPGV algorithm, the static and flow area could be more significantly classified. Using the split-spectrum method, the SSAPGA does provide higher ISNR and more uniform background intensity (as demonstrated in Fig. 6(c)). However, the multiple angiography calculations as a result of the split spectrum caused a longer processing time. From Table I, we noted that the CDV and CCPGV methods have the lowest computational cost and were about 3 times faster than the SSAPGA method.
If a two-dimensional averaging window is implemented, the SNR and motion contrast of the angiogram may be also improved [45]. However, it is worth noting that the processing time can also inevitably increase [9], [38], which is clinically undesirable. It has been reported that the 3% milk solution used in the phantom experiment has an enhanced vascular imaging signal due to more backscattering than red blood cells [46]. However, in our study, the same data set from the phantom experiment was processed by three different post-processing algorithms. As a result, the differences in the compared results are mainly due to the different characteristics of the algorithms themselves. Therefore, this conclusion is still relevant. To further validate the clinical feasibility and performance of the proposed method, it is necessary to perform in vivo experiments.

C. In Vivo Experiments
To demonstrate the applicability and performance of the newly proposed method for in vivo mapping of the microcirculation, the palmar region of a 26-year-old healthy male with informed consent was imaged. Before imaging, the subject was seated in an upright position and allowed to acclimatize in the laboratory for 10 minutes. After acclimatization, the palm was placed under the OCT probe and a soft cushion was placed underneath the palm. To flatten the skin surface, an optical window was attached in front of the probe, and the space between the window and the skin was filled with ultrasound gel for refractive index matching. The cross-sectional results are illustrated in Fig. 7. The scattering map, as shown in Fig. 7(a), can reveal typical skin-layered structures, including the weakly scattering epidermis and strongly scattering dermis. Small hyper-scattering spots pointed by the blue arrows are the cross-sections of unevenly distributed ultrasound gel. Fig. 7(b), (c), and (d) are the angiograms processed by the CDV, SSAPGA, and CCPGV methods, respectively. The white spots with tailed artifacts are blood vessels in the dermis, as pointed out by the red arrows, while the other low-intensity data (shown in black appearance) are static backgrounds. Overall, all three algorithms successfully map the vascular, but subtle differences remain in revealing the capillaries. By qualitatively comparing, we can find that Fig. 7(d) shows the highest contrast, indicating an appealing motion contrast can be provided by   Fig. 7(b) to (d)). By using (10), the calculated metrics of motion contrast are 1.9679, 1.2148, and 2.5058, respectively for CDV, SSAPGA, and CCPGV. Compared to CDV (Fig. 7(b)) and SSAPGA (Fig. 7(c)), an increase of 27.3% and 106.3% were found for CCPGV. Therefore, it is concluded that the vessels can be distinguished more easily from the static tissues using the proposed CCPGV method.
It is worth noting that several random dots are displayed in the epidermis, as shown in Fig. 7(b) to (d). Such dots are the artifacts induced by the random noise, as the weakly scattering epidermis is clinically considered to be a non-vascular area of the skin. Therefore, like all angiographic methods, the CCPGV method is prone to noise-induced artifacts in the low-SNR area. As the SNR decrease, the random noise would domain in both the amplitude and phase of the complex-valued OCT signal. In particular, the effect of noise on vascular imaging is more Fig. 9. En-face images of human palm skin. CDV images before (a) and after (d) Hessian-based frangi vesselness filtering; SSAPGA images before (b) and after (e) Hessian-based frangi vesselness filtering; CCPGV images before (c) and after (f) Hessian-based frangi vesselness filtering. Scale bars = 0.5 mm.
pronounced in the air and deep tissue regions, as shown in Fig. 7(b) to (d). Fig. 8 presents the dependence of the calculated angiographic signal on the SNR in the static region. For easier comparison and analysis, the angiographic signals obtained by each of the three methods were normalized.
As depicted in Fig. 8, the CCPGV and SSAPGA values from the static region are approximately zero when the SNR is over 50 dB, whereas the CDV values are approximately zero when the SNR is over 65 dB. Due to the implementation of the square root in the CDV method, the calculated angiographic signals are slightly higher, as shown in Fig. 6(b). However, as the SNR decreases, the SSAPGA value increases the fastest. This indicates that in the low SNR regions, narrower spectrum bands might be more susceptible to noise, resulting in noise-induced artifacts. The CCPGV algorithm, to robustly eliminate the artifacts induced by bulk tissue motion and trigger jitter, adopts the gradient of the phase difference in (8) and combines it with depth directionality summation, rather than the split-spectrum method. The split-spectrum method is an averaging method, which can reduce the spectral-dependent speckle noise in the angiogram [47], [48]. Specifically, the full spectrum is split firstly into several narrower bands by multiplying with Gaussian windows. After splitting the spectrum, the data in each narrower band are processed separately to generate several flow images for the same B-scan location. Finally, the flow signal is acquired by averaging the narrower-band flow images. Although the signalto-noise can be improved to some extent, the sensitivity to flow inevitably is reduced [6] and the processing time is increased [47], as demonstrated in Fig. 6(c) and Table I. In contrast, the depth directional summation is essentially a calculation of the directional variance of the complex vectors [25]. Additionally, the CCPGV method can provide a more acceptable angiographic signal for a given SNR (when SNR exceeds 15 dB) as shown in Fig. 8. Thus, compared to the CDV and SSAPGA methods, the CCPGV method not only provides a significant improvement in motion contrast but also slightly better suppression for random noise without increasing processing costs. To remove the low-SNR regions, a signal-intensity-based threshold has been widely used, as shown in Fig. 6. Alternatively, the thought of the noise-bias correction method [26], SNR-corrected low-noise complex correlation [49], or inverse SNR-decorrelation method [50] can also be borrowed in the future.
To comply with the clinician's reading custom of images with which they are familiar, color-code depth projection is generated from the 3D data. As the superficial epidermis is a non-vascular region of the skin, Fig. 9 only shows the projections over a depth range of 318-818 μm into the skin. The first and second rows respectively demonstrate the results before and after the Hessianbased Frangi vesselness filtering [45]. We can find that this filter can significantly improve the contrast of the skin vessels while suppressing noises. Fig. 9(c) and (f) show the en-face images processed with the proposed CCPGV method. As pointed out by the blue arrows, more vessels and better connectivity can be visualized using CCPGV. However, part of the blood vessel details in the CDV images ( Fig. 9(a) and (d)) and SSAPGA images ( Fig. 9(b) and (e)) cannot be clearly observed due to the lower motion contrast. The low-motion contrast makes them more difficult to distinguish vessels from the background.
As mentioned above, the depth summation may lead to blurring along the depth axis, but this impact is relatively small [25]. Consistent with this conclusion, the en-face angiograms ( Fig. 9(c) and (f)) for human skin demonstrate that the axial blurring is hidden within the out-of-plane direction. In this study, human skin was taken as the sample to evaluate the performance of the proposed OCTA method and the corresponding conclusions were also drawn. Despite the growing interest in dermatological research, it has to be acknowledged that the OCTA technology has been adopted in many other applications, such as retina [52], neuro [53], brain [54], and colon [55]. In the future, we will image a wide range of biological tissues to obtain a more comprehensive understanding and analysis of the blood flow imaging performance, including but not limited to motion contrast.

IV. CONCLUSION
In summary, we present an angiographic algorithm called CCPGV. Although this method is phase-sensitive, it is highly immune to phase artifacts caused by bulk motion and trigger jitter. The phantom experiments and in vivo experiments were then implemented to compare the performance of CCPGV with phase-resolved Doppler variance, CDV, and SSAPGA. In contrast to the phase-resolved Doppler variance, CCPGV can intrinsically reject undesirable phase shifts caused by bulk motion and trigger jitter. Additionally, the results demonstrate that the CCPGV method could present approximately 3 and 1.4 times greater motion contrast in phantom angiograms than the SS-APGA and CDV, respectively. In vivo experiments, the motion contrast can be improved by 27.3% and 106.3%, respectively. Due to these advantages, more vessels and superior vascular connectivity can be visualized. It is also found that this method is less affected by the noise-induced artifacts in the low-SNR area than the SSAPGA and CDV methods. Therefore, we believe our method will benefit the biomedical community in disease diagnosis and monitoring to some degree.