1. Introduction

Various OCT techniques have been developed for blood vessel imaging, such as optical Doppler tomography (ODT), optical angiography (OMAG), speckle variance OCT (svOCT) and etc [1–23]. All these techniques detect blood flow by comparing OCT signals (phase or amplitude) acquired at different time points (between Ascans or between Bscans) to sense particles in motion (flowing blood cells). OCT based blood vessel imaging does not require exogenous contrast agent, which makes it favorable for in vivo imaging. In particular, OCT angiography can be used to image subcutaneous microvasculature in situ and non-invasively; and is considered to have great potential for the diagnosis of skin diseases [13]. Blood flows within the dermis of skin are usually very slow (in the range of 0.1-0.9μm/s) [14], therefore have to be imaged by high sensitivity angiography that compares OCT signal acquired with larger time interval. For example, subcutaneous microvasculature image can be obtained by calculating interframe instead of inter-Ascan variation of OCT signal.

One of the major challenges for high sensitivity in vivo microvasculature imaging is bulk motion (BM) that comes from sample. BM inevitably exists during the imaging of living animal or human and can significantly degrade image quantity of OCT angiography. This is because motion signal due to BM contributes to every pixel in an angiography image, while motion signal caused by blood flow contributes to a limited number of pixels corresponding to vessels. In this study, we investigate high sensitivity, interframe signal variation based OCT angiography (isvOCA) for subcutaneous microvasculature imaging. isvOCA signal is obtained by calculating pixelwise interframe signal variation between Bscans. According to our theoretical derivation, isvOCA signal has an explicit functional dependency on the magnitude of motion in x, y and z directions. For the visualization of microvasculature, signal of isvOCA is predominantly derived from individual particle’s random motion rather than ensemble directional motion that is quantified by blood flow. isvOCA is particularly vulnerable to BM, because the time interval between the acquisition of subsequent Bscans is long enough for BM to induce detectable changes in OCT signal. Therefore, it is essential to reduce motion artifact to obtain high sensitivity isvOCA.

To remove motion artifact in Doppler OCT angiography, it was usually assumed that the sample movement induced a constant phase term for each Ascan and motion artifact could be removed by subtracting this phase term [15–19]. However, this phase correction method for Doppler OCT cannot be directly applied to isvOCA, because of different contrast mechanisms for blood vessel imaging in these two techniques. isvOCA is sensitive to motion in three dimensions while Doppler OCT signal is only sensitive to motion in the axial direction. The artifact due to BM in speckle variance OCT could be corrected through image registration [23]. Similar to phase correction in Doppler OCT, this method tracked bulk motion by analyzing OCT signal and could only be valid when motion was constrained to one or two dimensions. However, BM is generally a 3D vector with x, y and z components to be determined by three independent variables at any time instant. The quantitative tracking of BM with high accuracy is extremely challenging due to its 3D nature. Nevertheless, when motion gradient over time is small within the interval between Bscans, BM induces the same displacement for every scatterer within the sample. Therefore, BM results in the same change in OCT signal for a collection of pixels. As a result, it is possible to estimate and remove image artifact caused by BM via appreciating its global nature.

For isvOCA, BM results in motion artifact because the non-blood pixel has non-zero signal. On the other hand, angiography signal corresponding to blood flow does not change significantly with BM. It was demonstrated that the high sensitivity of OCT angiography based on interframe signal variation was due to the random motion of scatterers rather than directional flow [7]. Therefore, if the background level associated with BM can be estimated and subtracted, BM induced image artifact should be removed from OCT angiography.

In this study, we analyzed isvOCA signal and conducted systematic motion analysis for isvOCA. We further proposed to quantify the BM induced motion artifact through 2D spatial averaging on isvOCA signal. As the magnitude of BM inevitably varies as time, an adaptive threshold corresponding to instantaneous magnitude of BM was applied to reduce motion artifact. To the best of our knowledge, systematic motion analysis and motion artifact removal through adaptive thresholding have not been investigated before for intensity variation based OCT angiography.

To perform motion analysis, we considered OCT signal at each pixel as a random variable and demonstrated theoretically and experimentally that the expected value of isvOCA signal had an explicit functional dependency on motion. Our results indicated that isvOCA signal at non-vessel pixels had the same expected value <v_BM> due to global BM. Assuming most pixels in an OCT image do not correspond to blood vessels, <v_BM> could be simply assessed through spatial averaging of the entire 2D isvOCA image. Using <v_BM> as an adaptive threshold value, motion artifact reduction could be achieved by thresholding the original cross-sectional isvOCA images. We demonstrated the principle of the proposed motion artifact removal method and evaluated its effectiveness on phantoms. We implemented the algorithm for angiography reconstruction with the proposed motion artifact reduction in real-time using graphic processing unit (GPU) [8, 23, 24] and conducted in vivo imaging of subcutaneous blood vessel network on human subject. It is worth mentioning that the isvOCA is used to provide high sensitivity microscopic visualization of blood vessel network. The quantification of blood flow is beyond the scope of this work.

2. Principle

In the following analysis, a Cartesian coordinate system (x,y,z) is used (x: Bscan direction; y: Cscan direction; z: Ascan direction). I(x,y,z,t) indicates OCT signal obtained at spatial location (x,y,z) and time t. I(x,y,z,t) can be considered as a random variable with expected value I₀. v(x,y,z,t), the temporal variation of signal, is calculated to generate motion image for OCT angiography (isvOCA), as shown in Eq. (1). It is essential to normalize the right hand side of Eq. (1) using I₀². Otherwise, the magnitude of motion image also depends on the strength of scattering.

The expected value of v(x,y,z,t) is calculated in Eq. (2), where < > indicates taking the expected value of a random variable.

To find the functional relationship between motion and <v(x,y,z,t))>, we consider I(x,y,z,t-δt) = I(x-δx,y-δy,z-δz,t), assuming that scatterers at (x,y,z) displace for (δx, δy, δz) within the time interval from t-δt to t. In other words, the following conditions are equivalent: (1) OCT measurement performed at different time points from a moving sample with the beam at the same spatial coordinate; (2) OCT measurement performed with the beam at different spatial locations from a static sample “frozen” at a time point. Such equivalency is shown in Eq. (3)

Equivalence shown in Eq. (3) allows us to apply speckle tracking developed for motion analysis in ultrasound as well as OCT to investigate BM in isvOCA [25–29]. Based on previous studies, we can further writhe Eq. (3) as Eq. (4) where S(x,y,z) is the 3D convolution of scattering distribution function a(x,y,z) with OCT system's 3D point spread function (PSF) h(x,y,z).

Plugging the expression of S(x,y,z) and S(x-δx,y-δy,z-δz) into Eq. (4) and taking the fact that h(x,y,z) is a deterministic rather than random function, we can re-write Eq. (4) as:

Assuming that the speckle is fully developed and thus scatterers in different spatial location are described by identical but independent random variables, we have Eq. (7) in which a₀ is a constant representing the scattering strength.

Taking the sifting property of delta function, Eq. (6) can be simplified.

A simplified model for OCT system allows us to express h(x,y,z) as a Gaussian function (Eq. (9): ω_x, ω_y and ω_z are the 1/e width of the PSF in x, y and z directions; B is a constant factor taking account of system efficiency and responsivity) [30].

Using Eq. (8) and (9), we can express <v(x,y,z,t)> in Eq. (10) where A = (a₀²B²/I₀).

As shown in Eq. (10), the expected value of v(x,y,z,t) has an explicit functional dependency on motion. When δx = δy = δz = 0, <v(x,y,z,t)> = 0. Otherwise, <v(x,y,z,t)> takes a positive value. δx, δy and δz generally depend on spatial and temporal coordinates, x, y, z and t; therefore, <v(x,y,z,t)> varies as spatial and temporal coordinates. According to Eq. (10), the contrast mechanism of isvOCA is significantly different from Doppler OCT. First, Doppler OCT is only sensitive to motion in axial direction (δz) while isvOCA is sensitive to motion with arbitrary projection in δx, δy and δz directions. Second, Doppler phase shift is proportional to δz, therefore motion in opposite directions cancels each other out when contributing to Doppler signal. This does not happen in isvOCA, as image contrast essentially derives from the square of motion in x, y and z directions.

If the only source of motion within a sample is the motion of particles in liquid blood, δx, δy and δz are non-zero only for scatterers within liquid blood. According to Eq. (10), <v(x,y,z,t)> has a small value at pixels corresponding to static solid tissue but has a large value at pixels corresponding to flowing blood. Therefore, temporal variation of structural OCT signal can provide contrast for microvasculature visualization. It is worth mentioning that within flowing blood, particles have random motion as well as directional flow. The magnitude of random motion is significantly different from the ensemble directional motion that is quantified by blood flow. The random motion can be quantitatively evaluated using mean square displacement (MSD) between subsequent Bscans denoted as ${\overline{\delta x}}_{MSD}$, ${\overline{\delta y}}_{MSD}$, and ${\overline{\delta z}}_{MSD}$. When signal decorrelation is primarily due to random motion within liquid blood and any one of ${\overline{\delta x}}_{MSD}$, ${\overline{\delta y}}_{MSD}$, and ${\overline{\delta z}}_{MSD}$ is significantly larger than the spatial resolution in respective direction, <v(x,y,z,t)> approximately equals 2A² that is independent of flow speed. This was demonstrated experimentally in previous work by H. Ren et al [7].

On the other hand, when BM induces identical spatial displacement δx = δx_BM(t), δy = δy_BM(t), δz = δz_BM(t) to all pixels, <v(x,y,z,t)> is no longer zero for solid tissue and the angiography image suffers from motion artifact. If δx_BM(t), δy_BM(t), δz_BM(t) can be quantitatively extracted to restore OCT image, motion artifact in angiography image can be removed. However, simultaneous tracking of δx_BM(t), δy_BM(t), δz_BM(t) is extremely challenging due to the 3D nature of BM. Here we propose a method to reduce motion artifact without quantitative tracking of 3D motion.

It is noted that v(x,y,z,t) has non-zero expected value due to BM. However, this expected value denoted as <v_BM(t)> is the same for solid tissue at different spatial locations, because δx_BM(t), δy_BM(t), δz_BM(t) do not depend on spatial coordinate. Therefore, BM effectively results in an increased background noise level for non-vessel tissue. On the other hand, signal corresponding to blood derives from random motion of particles within flowing blood and would remain the same magnitude when the sample is interrupted by bulk motion. Therefore, we can reduce motion artifact in isvOCA, by assessing the increased background level due to BM (<v_BM(t)>) and apply <v_BM(t)> as an adaptive threshold value to isvOCA image. As the value of <v_BM(t)> is independent of spatial location, we can estimate <v_BM(t)> through spatial averaging if the sample volume does not contain blood vessels. For sample volume containing blood vessels, spatial averaging of v(x,y,z,t) also approximately equals <v_BM(t)>, because only a few pixels within the volume are corresponding to flowing blood and signal variation from blood flow plays a much less significant role compared to bulk motion in spatial averaged signal.

For angiography imaging, we acquire multiple structural 2D images (N frames of Bscan) at the same plane in the elevation direction (the same y coordinate: y = y₀) and calculate the normalized signal variation for the same pixel in different frames, as shown in Eq. (11).

Here discrete signal is related to continuous signal as: I_ij,n = I(x₀ + id_x,y₀,z₀ + jd_z, t₀ + nΔt) with d_x and d_z indicating spatial domain sampling interval.

To quantitatively evaluate BM induced artifact, we calculate the spatial average of each frame of signal variation obtained from Eq. (12) (N_x: number of Ascans in a Bscan; N_z: number of pixels in an Ascan).

Equation (12) indicates taking spatial average over the entire 2D isvOCA image. Such extensive average can provide a reliable approximation of BM induced increase in image background. If spatial average is performed using a smaller amount of pixels, the estimation of BM can have a higher spatial resolution. However, with fewer pixels involved in the averaging, motion signal from blood vessel might contribute substantially in $\overline{\upsilon }$_n and $\overline{\upsilon }$_n cannot be used to quantify artifact due to BM.

According to our previous discussion, $\overline{\upsilon }$_n that is the spatial average of v_ij,n can be used to approximate <v_BM>, the expected value of v(x,y,z,t) induced by BM in non-blood tissue, assuming most pixels in OCT image do not correspond to blood vessels. Therefore, we can use $\overline{\upsilon }$_n as an adaptive threshold value to threshold v_ij,n and reduce motion artifact, as shown in Eq. (13).

Afterwards, we perform temporal average on thresholded OCT signal variation ${\widehat{\upsilon }}_{ij,n}$ to generate cross-sectional microvasculature image at the given y coordinate, as shown in Eq. (14).

The reconstruction of cross-sectional OCT angiography with motion artifact reduction is summarized in Fig. 1. 3D microvasculature visualization can thus be obtained from a series of cross-sectional angiography images acquired at different y locations.

 

Fig. 1 Reconstruction of motion artifact reduced OCT angiography.

Download Full Size | PDF

3. OCT system configuration and signal processing

We used a spectral domain OCT (SD OCT) system at 1.3μm wavelength range for this study. The broadband light source (superluminescent diode, SLD1325 Thorlabs, 100nm bandwidth, corresponding to a 7.4 μm axial resolution) illuminates the reference and sample arm of a fiber-optic Michelson interferometer through a fiber-optic coupler. In the sample arm, a 3X scanning lens (SLM04, Thorlabs) is used in front of the specimen to focus the probing beam and collect back scattered photons. Light returned from sample and reference arms is routed by the fiber-optic coupler to the spectrometer where interference signal is detected by a CMOS InGaAs camera (SUI1024LDH2, Goodrich). A frame grabber (PCIe-1433, National Instrument) takes the interferometric signal from the camera and the signal is further processed. The system synchronization is achieved by triggering the camera exposure and controlling the galvanometer using the same data acquisition card (NI DAQ USB-6212 BNC, National Instrument) with the same time base. Two galvanometers scan the light beam in x and y directions to cover the transverse plane for 3D imaging. All the device controls and signal processing were performed using a host computer (Dell Precision T7600) with software developed in C + + (Microsoft Visual Studio, 2012).

Signal processing for motion artifact reduced angiography is illustrated in Fig. 2(a). In brief, N frames of OCT Bscans are obtained at the same y plane; afterwards, interframe signal variation is calculated to obtain motion image according to Eq. (11); motion artifact is then assessed through spatial averaging according to Eq. (12); adaptive thresholding is applied to each frame of motion image to reduce BM induced artifact as in Eq. (13); multiple frames of artifact reduced motion image is averaged to generate cross-sectional angiography image using Eq. (14). Signal processing was implemented in real-time using computer unified device architecture (CUDA) on a graphic processing unit (GPU, NVIDIA Geforce GTX 780, 2304 CUDA cores at 0.9GHz, 3 GB graphics memory). The speed of signal processing was analyzed using NVIDIA Nsight that is embedded in Microsoft Visual Studio. The average time to process an Ascan on GPU for the angiography image reconstruction is shown in Fig. 2(b). With the cubic spline interpolation (wavelength to wavenumber, or λ to k) and interframe signal variation calculation taking the longest computation time, the total average time to process one Ascan for angiography imaging is about 10μs and is less than the minimum time interval for data acquisition determined by the highest line rate (92kHz) of the CMOS camera. Therefore, the imaging speed of our system is limited by the camera data acquisition rather than signal processing. With the maximum line scan rate of our CMOS camera equal to 92kHz, it takes about 3.6 seconds to generate a 256 × 256 en face angiography image by acquiring 5 Bscans (N = 5) at each y coordinate.

 

Fig. 2 (a) Data processing flowchart of motion artifact reduced angiography; (b) Time expenditure of each processing step preformed on the GPU.

Download Full Size | PDF

4. Results

To show that Eq. (10) is valid for BM in a solid sample, we imaged a solid scattering phantom and introduced bulk motion with known magnitude during imaging. The phantom was fabricated by curing the mixture of epoxy and titanium dioxide that served as scatterer. Bscans were obtained by steering the light beam in x direction using a galvanometer. In addition, we used another galvanometer to translate the beam in y direction to introduce a known beam displacement δy_BM between the acquisitions of subsequent Bscans. As a result, every pixel in a Bscan was displaced for the same magnitude in y direction between adjacent frames. For each value of δy_BM, we acquired multiple Bscans I_ij,n (n = 1,2,…,25), used Eq. (11) to obtain interframe signal variation v_ij,n, used Eq. (12) to obtain the spatial average of v_ij,n denoted as $\overline{\upsilon }$_n, and calculated ${\overline{\upsilon }}_{\delta y_{BM}}$: ${\overline{\upsilon }}_{\delta y_{BM}}=({\displaystyle {\sum }_{n=2}^N{\overline{\upsilon }}_n})/(N-1)$. N = 25 is chosen because this value can provide sufficient temporal averaging to reduce random noise in the obtained values of ${\overline{\upsilon }}_{\delta y_{BM}}$. However, when N is significantly larger than 25, the denoise effect of temporal averaging saturates. Results obtained at different magnitudes of bulk motion in y direction are shown in Fig. 3 as red circles. Error bar at each data point corresponds to standard deviation between $\overline{\upsilon }$_n for n = 2, 3, …, 25, at each δy_BM. Clearly, ${\overline{\upsilon }}_{\delta y_{BM}}$ increases as δy_BM. To further demonstrate the functional dependency of ${\overline{\upsilon }}_{\delta y_{BM}}$ on δy_BM, we used the known values of δy_BM and ${\overline{\upsilon }}_{\delta y_{BM}}$ to fit the following model derived from Eq. (10) with δx = δz = 0 and δy = δy_BM: ${\overline{\upsilon }}_{\delta y_{BM}}=2A^2[1-\exp (-\delta y_{BM}^2/{\omega }_y^2)]+\beta$. ω_y took a value of 17μm according to the optical property of scanning lens used in the sample arm. The constant $\beta$ was included in this model to take account for signal decorrelation due to random perturbation from ambient environment. Least square fitting of this model is shown as the black curve that is highly consistent with measured data points. The fitting lead to $\beta$ = 0.2075 and A = 1.1. The R² statistics of the fitting is larger than 0.99, indicating Eq. (10) appropriately models the functionally dependency of isvOCA signal on the magnitude of bulk motion.

 

Fig. 3 Measured spatial average of intensity variation (red circles with errorbar) and fitted curve (black).

Download Full Size | PDF

To further validate Eq. (10) for BM in different directions, we performed Bscan on the same phantom as described above and introduced BM by applying additional voltage to x and y galvanometers respectively so that BM had both x and y components: δx_BM = δdcosθ and δy_BM = δdsinθ. The direction of motion was determined by θ and the magnitude of motion was determined by δd. With different values of θ and δd, we obtained multiple sets of structural OCT images. Similarly, we calculated v_ij,n, $\overline{\upsilon }$_n and the average of $\overline{\upsilon }$_n denoted as ${\overline{\upsilon }}_{\delta d,\theta }$ for different BM. Results are shown in Fig. 4 with δd as horizontal axis. Data points in different colors correspond to different θ, or different direction of motion. Figure 4 indicates ${\overline{\upsilon }}_{\delta d,\theta }$ remains the same for transverse motion in different direction but with the same magnitude. This is because of the circular geometry of imaging optics that leads to ω_x = ω_y = ω_t. When transverse motion is considered and δz_BM = 0, ${\overline{\upsilon }}_{\delta d,\theta }$ can be expressed as: ${\overline{\upsilon }}_{\delta d,\theta }=2A^2\{1-\exp [-(\delta y_B^2+\delta x_B^2)/{\omega }_t^2]\}+\beta =2A^2[1-\exp (-\delta d^2/{\omega }_t^2)]+\beta$. Clearly, ${\overline{\upsilon }}_{\delta d,\theta }$ does not depends on θ when δx_BM = δdcosθ and δy_BM = δdsinθ and this is consistent with results in Fig. 4. In Fig. 4, larger inconsistency can be observed between results obtained with motion in different directions at larger displacements, because ${\overline{\upsilon }}_{\delta d,\theta }$ suffers from noise that is larger when δd is larger. This is also illustrated in Fig. 3 where error bars are larger for larger displacement.

 

Fig. 4 Spatial average of motion image obtained with BM that has different magnitude and direction.

Download Full Size | PDF

To demonstrate the principle of adaptive thresholding in motion artifact removal, we established a phantom by embedding a transparent polyimide tube in solid scattering medium (cured mixture of epoxy and titanium dioxide) and filling the polyimide tube with bovine milk. We used milk as a scattering fluid to mimic the blood flow because the fat particles in milk scatter light similarly to red blood cells (RBC’s) in blood [31]. We performed structural and isvOCA imaging on the phantom. Figure 5(a) shows the Bscan OCT image (log scale) of the phantom and Fig. 5(b) shows the isvOCA image (linear scale) obtained without bulk motion. In Fig. 5(b), the tube filled with liquid has large signal amplitude due to random motion of protein particles within the liquid milk. In comparison, the solid part of the phantom has small signal amplitude in isvOCA image because the scatterers were static. Notably, the capillary filled with scattering liquid was placed in horizontal plane and we did not use any external mechanism to generate direction flow. Therefore, motion signal was completely derived from random motion of scatterers in liquid medium.

 

Fig. 5 Structural (a) and flow (b) image of a phantom with a polyimide tube filled with bovine milk embedded in solid scattering medium

Download Full Size | PDF

To illustrate how bulk motion degrades the quality of isvOCA image, we introduced BM in x and y directions respectively in addition to raster scanning pattern by applying pre-defined voltage to the galvanometer. This was done in a way similar to our previous experiments on solid phantom. We calculated interframe signal variation v_ij,n. Afterwards, we performed spatial averaging on v_ij,n for solid and liquid parts of the phantom as enclosed by rectangles in Fig. 5(a) to obtain $\overline{\upsilon }$_solid,n and $\overline{\upsilon }$_liquid,n. $\overline{\upsilon }$_solid and $\overline{\upsilon }$_liquid were thereafter obtained through averaging $\overline{\upsilon }$_solid,n and $\overline{\upsilon }$_liquid,n for different values of n. $\overline{\upsilon }$_solid and $\overline{\upsilon }$_liquid (normalized) obtained with different magnitudes of x (red circles) and y (green circles) displacement are shown in Fig. 6(a) and 6(b), respectively. Clearly, $\overline{\upsilon }$_solid increases as increased bulk motion, which is consistent with Eq. (10) and our experimental results in Fig. 3 and 4. On the other hand, $\overline{\upsilon }$_liquid remains approximately the same at different magnitude of bulk motion, because it was the random motion of particles within the tube that resulted in complete decorrelation between signals in different frames. This is consistent with our previous analysis and also consistent with experimental results in a previous study that stated the high sensitivity of microvasculature visualization based on interframe OCT signal variation was due to random Brownian motion of scatterers [7]. As signal at solid part of the sample increases with BM while signal remains the same at liquid part of the sample, suppressing background signal due to BM by thresholding can effectively enhance the visualization of flow in isvOCA.

 

Fig. 6 (a) $\overline{\upsilon }$_solid increase as motion in both x and y directions; (b) $\overline{\upsilon }$_liquid remains almost constant with bulk motion

Download Full Size | PDF

To demonstrate the effectiveness of adaptive thresholding in cross-sectional flow imaging, we obtained multiple Bscans from the same phantom as shown in Fig. 5(a) at the same y coordinate. We calculated interframe signal variation using Eq. (11), and performed temporal averaging on the signal variation obtained to generate flow image. Figure 7(a) and 7(b) were obtained without adaptive thresholding when bulk motion was in x and y directions respectively. The magnitude of displacement between adjacent Bscans is 6μm. Due to bulk motion, flow signal at the solid part of the phantom does not diminish. For comparison, we applied our adaptive thresholding method to the same OCT data set used to obtain Fig. 7(a) and 7(b). The results are shown in Fig. 7(c) and 7(d) where solid part of the sample appears to be dark. To quantitatively evaluate the enhancement in image quality by using adaptive thresholding, we calculated the image contrast defined as C = 20log₁₀[($\overline{S}$_liquid-$\overline{S}$_solid)/$\overline{S}$_solid] where $\overline{S}$_liquid and $\overline{S}$_solid were obtained by calculating the averaged isvOCA signal within rectangles shown in Fig. 5(a). Figure 7(e) and 7(f) show contrast of isvOCA image obtained without (blue) and with (green) adaptive thresholding when BM was in x and y directions. Figure 7(g) and 7(h) show adaptive thresholding could enhance the contrast of isvOCA. Clearly, a 5 to 10 dB enhancement in contrast was achieved through adaptive thresholding.

 

Fig. 7 Cross-sectional flow image obtained without adaptive thresholding when bulk motion was in x dimension (a) and y dimension (b); Cross-sectional flow image obtained with adaptive thresholding when bulk motion was in x dimension (c) and y dimension (d); contrast of flow image, with (green) and without (blue) adaptive thresholding at different magnitude of motion in x direction (e) and y direction (f); contrast enhancement through adaptive thresholding when motion was in x direction (g) and y direction (h).

Download Full Size | PDF

Results in Fig. 7 suggest that the proposed adaptive thresholding method can reduce BM induced artifact more effectively when the magnitude of angiography signal derived from global BM is smaller. This is because, v(x,y,z,t), the interframe signal variation that provides contrast for motion imaging, is a random variable and can have different values in the vicinity of its expected value. When large BM results in large signal decorrelation, pixels corresponding to solid tissue or flow cannot be differentiated by their magnitudes and thus adaptive thresholding would fail to reduce motion artifact. As shown in Fig. 7(c) and 7(d), although background due to BM is significantly suppressed through adaptive thresholding, signal due to particle motion within liquid at the bottom of the polyimide tube also diminishes. Moreover, Fig. 7(g) and 7(h) indicate contrast enhancement obtained from adaptive thresholding decreases when BM has a larger magnitude.

Using real-time software based on GPU, we conducted in vivo subcutaneous microvasculature imaging on healthy human volunteer and obtained motion artifact reduced isvOCA. A cross-sectional structural OCT image and the corresponding cross-sectional isvOCA image of the palm skin are shown in Fig. 8(a) and 8(b). 5 Bscans were obtained at each y coordinate to generate cross-sectional isvOCA. To obtain en face visualization of angiography image, raster scanning was performed in both x (fast axis) and y (slow axis) directions to obtain 3D OCT data volume. To reduce the 3D data volume for en face visualization, we averaged angiography signal in z direction at small (500μm - 650μm), medium (650μm - 800μm) and large (800μm - 1100μm) depths (referred to zero-delay plane) as indicated by blue, green and red arrows in Fig. 8(a). En face blood vessel maps at small, medium and large depths are shown in Fig. 8(c), 8(d) and 8(e). Figure 8(c) shows a smaller amount of superficial capillaries that follow the arrangement of friction ridges of human hand. More blood vessels are seen in Fig. 8(d). Figure 8(e) suggests more and larger blood vessels are located further from the skin surface. We assigned pixel values in Fig. 8(c), 8(d) and 8(e) to blue, green and red channels in a RGB color image as shown in Fig. 8(f). Scale bars in Fig. 8 represent 500μm. isvOCA imaging was also conducted on human fingertip and a skin lesion with results shown in Fig. 9(a) and 9(b). Figure 9(a) shows that microvasculature at fingertip follows the fingerprint, which is consistent with previous studies using OCT and scanning electron microscopy (SEM) [32]. Figure 9(b) shows larger blood vessels connected to the lesion to supply its increased metabolism activity.

 

Fig. 8 (a) structural OCT image of human palm skin (E: epidermis; D: dermis); (b) cross-sectional angiography image that highlights blood vessel; (c) – (e) en face microvasculature image in small, medium and large depths; (f) microvasculature image that encodes signal depth with color.

Download Full Size | PDF

 

Fig. 9 isvOCA image obtained from human fingertip (a) and a small skin lesion (b).

Download Full Size | PDF

As discussed before, adaptive thresholding can effectively remove artifact when BM is sufficiently small in magnitude compared to the random motion of scatterers within liquid blood. To demonstrate this, we conducted in vivo imaging from the palm of our subject. Figure 10(a) and 10(b) were obtained from the palm of subject using the same set of raw structural OCT signal, without and with adaptive thresholding. During data acquisition for Fig. 10(a) and 10(b), the magnitude of bulk motion was minimized with a hand rest to constrain the bulk motion. However, motion artifacts still exist in Fig. 10(a), shown as the horizontal lines. In comparison, Fig. 10(b) does not have horizontal lines and has a significantly improved image contrast, indicating the effectiveness of adaptive thresholding in motion artifact removal. The hand rest was removed in acquiring data for Fig. 10(c) and 10(d). With larger magnitude of motion, adaptive thresholding cannot completely remove motion artifact. Although blood vessels are more visible in Fig. 10(d) compared to Fig. 10(c), the image quality of Fig. 10(d) is inferior to Fig. 10(b) and horizontal line artifacts still exist in Fig. 10(d). This is consistent with results from controlled experiment shown in Fig. 7.

 

Fig. 10 Angiography image obtained (a) with small magnitude of BM without adaptive thresholding; (b) with small magnitude of BM with adaptive thresholding; (c) with large magnitude of BM without adaptive thresholding; (b) with large magnitude of BM with adaptive thresholding.

Download Full Size | PDF

It is worth mentioning that the CMOS camera was operated at its highest speed (92 kHz, line scan rate) during our experiments. For 2D imaging on phantom with results shown in Fig. 3-7, each Bscan contained 1024 Ascans and the time interval between each Bscan was 0.011s. For in vivo imaging of human skin, we acquired 512 Ascans in each Bscan, so that large enough field-of-view could be achieved at a reasonably high frame rate to illustrate morphological feature of subcutaneous microvasculature. Therefore, in vivo imaging was conducted at a 180Hz frame rate or a 0.0056s time interval between adjacent Bscans.

5. Discussion

Adaptive thresholding can reduce motion artifact in isvOCA because signal level of flowing blood does not depend on BM while noise floor does. Physical origin of image contrast in OCT angiography based on interframe signal variation was attributed to random Brownian motion by H. Ren et al [7]. The magnitude of Brownian motion of blood cells can be estimated according to Einstein’s theory (〈[Δd(t)]²〉 = 2Dt) and measured experimentally [33–35]. However, the square root of 〈[Δd(t)]²〉is smaller than 1μm for the experimental setting in this study. In other words, the magnitude of Brownian motion is not large enough to completely decorrelate OCT signal between Bscans. Therefore, Brownian motion of individual blood cells might not be the most appropriate description of random motion that decorrelates OCT signal. A more appropriate modeling is needed.

It is assumed throughout the manuscript that BM resulted in identical displacement for all pixels corresponding to solid tissue. However, biological tissue is not rigid and the assumption of identical displacement is not valid when tissue deforms. However, if the mechanical property of tissue can be considered as homogeneous within sub-regions and we can quantify the impact of BM by averaging signal from a smaller number of neighboring pixels, adaptive thresholding can still be used to reduce motion artifact.

We conducted experiments with BM in x and y directions by simply applying pre-defined voltage to the galvanometers. However, as shown in Eq. (10), <v> depends on motion that can have non-zero components in x, y and z directions. Therefore, our method for motion artifact reduction is also valid for motion with axial component. In fact, during in vivo imaging of microvasculature, BM could be in arbitrary directions and our method successfully removed motion artifact, indicating adaptive thresholding is effective for BM with both lateral and axial components.

Large motion artifact as shown in Fig. 10(c) and 10(d) can be easily removed using a larger value to threshold the interframe signal variation signal. This can be achieved by replacing${\overline{\upsilon }}_n$ in Eq. (13) with$A{\overline{\upsilon }}_n$ where A is a constant larger than 1. In real application, the maximal bulk motion amplitude that current method can remove effectively depends on the speed of blood flow of interest. Artifact due to large motion can be suppressed with a larger threshold which however also excludes signal due to blood flow from being detected, as shown in Fig. 7. As a result, the trade-off between the motion artifact removal and blood flow detection is inherent to the thresholding based motion artifact removal.

Our theoretical analysis in this study is based on fully developed speckle. With the assumption of fully developed speckle, <v(x,y,z,t)> only depends on motion but does not depend on local optical property of tissue, as shown in Eq. (10). This implies that contrast of flow image based on normalized signal variation only derives from motion. However, speckle does not always develop fully in practice. As a result, the correlation between spatially shifted scattering distribution functions is not a Dirac Delta function and the expected value of normalized signal variation depends on both motion and optical property of sample. Nevertheless, signal due to bulk motion in motion image can have the same value if tissue has homogeneous optical property, when fully developed speckle assumption does not hold.