Short-time series optical coherence tomography angiography and its application to cutaneous microvasculature

We present a new optical coherence tomography (OCT) angiography method for imaging tissue microvasculature in vivo based on the characteristic frequency-domain flow signature in a short time series of a single voxel. The angiography signal is generated by Fourier transforming the OCT signal time series from a given voxel in multiple acquisitions and computing the average magnitude of non-zero (high-pass) frequency components. Larger temporal variations of the OCT signal caused by blood flow result in higher values of the average magnitude in the frequency domain compared to those from static tissue. Weighting of the signal by the inverse of the zero-frequency component (i.e., the sum of the OCT signal time series) improves vessel contrast in flow regions of low OCT signal. The method is demonstrated on a fabricated flow phantom and on human skin in vivo and, at only 5 time points per voxel, shows enhanced vessel contrast in comparison to conventional correlation mapping/speckle decorrelation and speckle variance methods. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Angiography has provided a valuable means to investigate and assess the vasculature in normal and diseased tissue [1].It has been performed using invasive, high-resolution (submicrometer) microscopy with 100 µm-1 mm field of view (FOV) and using non-invasive methods that have worse resolution but broader FOV, such as magnetic resonance imaging, computed tomography and position emission tomography [1][2][3].Bridging these two extremes is an extension of optical coherence tomography (OCT) [4], optical coherence tomography angiography (OCTA), which offers label-free, non-invasive imaging of small blood vessels, including arterioles, capillaries and venules, at intermediate resolution and FOV [5].Whereas OCT image contrast is determined by the level of backscattering in tissue, OCTA images the microvascular network via motion-induced changes in the OCT signal.OCTA achieves a resolution and FOV in the ranges of 2-20 μm and a few mm to ~20 mm, respectively [6], and an imaging depth limited to <1 mm, due to the strong scattering of OCT light in tissue.Thus, in vivo applications in humans focus mainly on transparent or superficial tissues, including the eye (largely the posterior retinal and choroidal microvasculature) [7][8][9], and the dermal skin [10,11].
OCTA has been widely applied in clinical ophthalmology with commercially available imaging systems [12,13].It has also been used to study cutaneous microvasculature and its temporal alteration in normal human subjects and patients with cutaneous conditions, such as wounds and burn scars [11,14,15], psoriasis [16], and skin cancer [17,18].The clinical applicability of OCTA imaging will likely be further improved by extending the FOV to image larger tissue areas [6], extending the imaging depth to visualize deeper vessels [19], improving the sensitivity of detection of the smallest vessels with low flow contrast [20], and improving the processing speeds to achieve real-time visualization.
OCTA identifies blood vessels by differences in the OCT signal versus time between that arising from moving scatterers in blood and that due to the surrounding largely static tissue.Such flow-induced differences are encoded in both the amplitude and phase of the complex OCT signal [5].OCTA variants utilize one of: temporal changes in the OCT amplitude signal (e.g., speckle variance and correlation mapping/speckle decorrelation) [10,11,[21][22][23]; the phase signal (e.g., Doppler OCT and phase variance) [24][25][26]; or the complex OCT signal (e.g., complex differential variance and OCT-based microangiography (OMAG)) [27][28][29].
Related methods that exploit the vasculature signature in the detected signal in the frequency domain, the main topic of interest here, have been explored for both velocimetry and angiography [30][31][32][33].In early work on velocimetry, Barton et al. [30] analyzed the frequency content of the OCT intensity signal using a sliding window.The ratio of the highest to the lowest frequencies was used to determine flow velocities in cross-sectional images.Also for velocimetry, Szkulmowski et al. proposed a joint spectral and time domain method to process multiple OCT A-scan spectra (n = 40) acquired from the same locations, estimating flow velocities via independent Fourier transformations in time and in the optical frequency domain [31].Szkulmowska et al. extended the method to three-dimensional (3D) velocimetry and angiography, and reduced the number of spectra acquired (n = 16) from each lateral location [32].For angiography, Matveev et al. recently utilized the fast-changing OCT signal components to visualize vessel structures by Fourier transform of closely adjacent, but not colocated, complex A-scans (acquired at a density of 16,384 A-scans/mm) to extract the frequency spectrum.This was followed by high-pass filtering to eliminate the low frequency (static) components and inverse Fourier transformation to highlight the vessel structures [33].The exceptionally high number of samples acquired within each B-scan practically limits the area that can be analyzed.
In this paper, we propose an alternative method that analyzes the frequency spectrum of the detected OCT signal taken from multiple acquisitions at a given voxel.The higher average magnitude of non-zero frequency components in the short-time series OCT signal intensity induced by the blood flow is used to identify vessels.Weighting by the OCT signal is incorporated to improve the contrast of blood vessels with low OCT signal and mitigate surface reflection artifacts.We characterize the method by studying the vessel contrast achieved from a flow phantom and from human skin tissue in vivo.We compare the shorttime series OCTA method to commonly used intensity-based OCTA methods, including speckle decorrelation (correlation mapping) and speckle variance [11,21,34].The results demonstrate, for a modest increase in acquisition times for a given A-scan rate, improved vessel contrast and visibility, in particular, for small vessels.Further, its simplicity lends itself to rapid calculation.These advantages suggest its potential for future applications.

Short-time series OCTA algorithm
The basic assumption underlying the method is that blood flow induces stronger non-zero frequency components in the OCT signal than those induced by the surrounding static tissue [31,33].As with other OCTA methods, the method first requires the acquisition of co-located B-scans (i.e., B-scans from the same lateral location) at multiple time points, throughout an acquisition volume.The OCT intensity signal (i.e., the modulus of the complex amplitude of the OCT signal) at the same voxel locations comprises a discrete time series with the nth sample at location (x, y, z) denoted by: 1 ( , , ; ) ( , , ; ( 1) ), where (x, y, z) is the voxel coordinate in the fast scanning, slow scanning and depth axes, respectively; I represents the OCT intensity signal as a function of the voxel coordinate with time point t n for n, an integer ranging from 1 to 2N + 1, where 2N + 1 is the total number of co-located B-scans (i.e., total number of time samples) acquired from the same lateral location; and T is the time interval between co-located B-scans.
The time series at each voxel in Eq. ( 1) is discrete Fourier transformed to obtain the complex frequency signal with the frequency components F denoted by: 0 ( , , ; ) ( , , ; ), where f 0 is the interval between neighboring discrete frequencies, determined by 1/[(2N + 1)T]; and m is the index of the (two-sided) frequency components ranging from -N to N. The average magnitude of the complex frequency signal at non-zero frequencies is then calculated as Alternatively, if many B-scans are acquired for analysis (i.e., 2N + 1 ≥ 29 for the scanning parameters used in this study), instead of a single frequency component, a narrow band centered on the zero-frequency component is excluded (i.e., high-pass filtered).This narrow band should be optimized for a particular tissue and setup, and will depend on the frequency spectrum recorded from static tissue.The optimization for human skin tissue, recorded using our system parameters, is shown in the Results section.However, there, we demonstrate that only a small number of co-located B-scans (~5) is required for practical imaging of the vessel network with our method.Thus, the elimination of only the zero-frequency component, as shown in Eq. (3), applies.After Fourier transformation, voxels with low OCT signal intensity lead to a correspondingly low magnitude of the complex non-zero frequency components, even if there is flow.To enhance the flow detectability at low OCT signal levels, we incorporate weighting by the inverse of the OCT signal intensity (i.e., zero-frequency component scaled by the number of co-located B-scans), given by: , 0 ( , , ; ) ( , , ) 2 1 ( , , ) , where ( , , ) I x y z is the mean OCT signal intensity and ( , , ;0) F x y z is the zero-frequency component of the 2N + 1 time samples at the same location.To avoid division by zero and over-emphasizing the signal in regions with excessive noise, ( , , ) I x y z is first averaged, and thresholded at an empirically chosen signal level of 16 dB above the noise floor to replace the low signal with the threshold.We used an averaging window of 3 × 3 pixels in the crosssectional plane, approximately 1.4 and 1.9 times the lateral and axial resolutions, respectively.It is used both for our method and for the accompanying speckle decorrelation calculation, empirically chosen to improve the signal-to-noise ratio (SNR) of angiography without significantly degrading the imaging resolution.We assume that an odd number of co-located B-scans are acquired for each lateral location for simplicity, but even numbers are applicable as well.The vessel contrast created by Eq. (3) and the further improvement introduced by weighting according to Eq. (4) will be demonstrated and discussed in Section 3.

OCT scanning of flow phantom and human skin
OCT scans were acquired using a commercial spectral-domain scanner (an upgraded TELESTO II, Thorlabs Inc., USA) to demonstrate the short-time series OCTA method on both a flow phantom and in vivo on normal human skin.The system has a center wavelength of 1300 nm and provides an imaging resolution of 5.5 µm (in air) and 13 µm, respectively, axially and laterally (as defined by the vendor).The scanner was operated at 76 kHz (Ascan/s), below its maximum of 146 kHz.Scans were acquired in one of two modes: 2D scanning by acquiring 200 co-located B-scans from a single lateral location with a FOV of 6 × 3.6 mm (1024 × 1024 pixels) in x and z directions, respectively; and 3D scanning with a FOV of 6 × 1.8 × 3.6 mm in x, y and z directions.In 3D scanning mode, 240 lateral (y) locations were scanned with a set of 5 co-located B-scans acquired from each location, using the same pixel sizes in x and z directions as in the 2D mode.It took approximately 4 and 21 s to acquire a scan, respectively, in the 2D and 3D scanning modes.In addition, the time interval between B-scans was 17.8 ms (~56 B-scans/s) for both 2D and 3D modes, leading to a discrete frequency spectrum with components up to 28 Hz.
For comparison to short-time series OCTA, speckle variance in the same 3D scans was calculated over the 5 co-located B-scans using the method presented by Mariampillai et al. [21].Speckle decorrelation was determined for each adjacent pair of co-located B-scans using the formula described in [15] with a window of 3 × 3 pixels in the fast scanning and depth axes.This led to 4 decorrelation B-scans from each lateral location, which were then averaged to generate a single enhanced decorrelation B-scan.In addition, the speckle decorrelation and speckle variance was weighted by the averaged and thresholded OCT signal at the corresponding pixels to reduce the noise [10].The same threshold level used in short-time series OCTA was used here.The same lateral averaging window (3 × 3 pixels) was applied to the short-time series and speckle variance images to ensure a fair comparison.
Blood vessels were mainly compared over a depth range of 300 µm from the skin tissue surface (determined from the OCT depth scan by assuming an average refractive index of 1.4) to ensure sufficiently strong signals from all three methods.For each method, the maximum OCTA signal of each A-scan in this depth range was used to generate a projection image of vessels.For visualization, the same colormap was used in the projection and cross-sectional OCTA images.The lower and upper thresholds were set at, respectively, the 50% and 99.5% points of the cumulative distribution function of the OCTA signal in the image.These thresholds were empirically chosen to maximize the vessel contrast without loss of vessels with low signal.The thresholds were then used to normalize the vessel signal for each method to the same range (0-1).For quantification, each projection image was processed to measure the vessel area density, defined as the ratio of the total vessel area to the total tissue area in the thresholded vessel image.The threshold was set using Otsu's method for each image [35].This method assumes two classes of pixels (i.e., foreground and background pixels) following a bi-modal histogram to calculate the optimal threshold to separate the two classes of pixels.The vessel images were thresholded using the optimal threshold estimated by this method to determine the vessel area and density.
The silicone flow phantom was fabricated in house by mixing Elastosil P7676A and P7676B fluid (Wacker Chemie AG, Germany) with titanium dioxide in a 3D-printed plastic container [36].The container was customized with two holes in the sidewalls to hold a small glass capillary (outer diameter: 80 µm; inner diameter: 50 µm) that mimicked a vessel orthogonal to the incident optical beam.After curing, the result was a capillary embedded in the silicone that mimicked the static tissue.The capillary was then connected to a syringe filled with a polystyrene microsphere suspension (nominally 0.5-µm or 2-µm diameter) to mimic the blood flow.The syringe was connected to a pump (Fusion 200, Chemyx Inc., USA) to introduce and control the flow speed.The scattering properties of the phantom were adjusted by tuning the ratio of titanium dioxide to Elastosil P7676A and P7676B so that the phantom had a signal attenuation that approximately matched the attenuation of normal human skin.
Human subjects (n = 4) were recruited for in vivo OCT scanning with ethics approval from the Human Research Ethics Committee of The University of Western Australia.Written consent was a including one region from t imaging.To r skin surface t center hole (5 marker to che negligible ves imaging prob phantom and

Results
This section f its optimizatio difference bet signal is also compared to t

Vessel c
The blood ve frequency co contrast, obta plots the Four phantom scan respectively, the adjacent magnify the v    The frequ trace) is ~20 frequencies a region (for flo higher (Fig. 1 Similar pl with a single consistent co parameterized In addition, F frequencies h magnitude at static tissue in the phantom from the pul following sec interval is m component, a

Choice o
The frequenc scans, chosen are not practi samples whils Fig. 2, show magnitude fo (dashed blue)

Discussio
The method p location as an content in ord pass) frequen to distinguish OCT scans a network.In c consistently i variance, espe skin tissue, fu will provide f Matveev analysis of ad their final v transformation frequency spe frequency com vessels by de decreased to 3 vessel image.Another diffe scans, instead sampling den method, phas complex OCT intensity sign features of ou for its ready a Mioseev e B-scans along adjacent B-sc domain, using they proposed domain (i.e., in the frequency domain and is not affected by such artifacts, even for as few as 5 time samples (i.e., 5 co-located B-scans).In addition, we chose to use co-located B-scans so as to avoid the degradation of vessel contrast due to the possibly corrupting spatial frequencyinduced frequency content resulting from analyzing adjacent B-scans.
The number of co-located B-scans acquired from the same location is an important parameter for the practical implementation of the short-time series method.We chose 5 in this study so as to minimize the amount of collected data and corresponding total acquisition time, whilst still attaining a high vessel/static tissue contrast in skin.A very large number of colocated B-scans does not significantly improve the vessel contrast, as we showed.This is due to the simultaneous increase of the vessel signal in both the static tissue and the flow regions.In addition, extraneous motion can become a problem as the acquisition time lengthens.A smaller number of co-located B-scans, on the other hand, leads to slightly poorer vessel contrast, but can still visualize most vessels for even 3 or 4 co-located B-scans.When the minimum number (i.e., 2) is used, the formula in Eq. ( 3) is then equivalent to the simplified OMAG method that generates contrast by subtracting the intensity signal in two co-located Bscans [40].This equivalence does not exist as more co-located B-scans are used, and which leads to superior vessel image quality, as compared to only 2 co-located B-scans (i.e., the subtraction method).The consistent optimization of this sampling parameter is as yet unverified and, thus, needs further investigation when applied to imaging microvasculature in other tissues, such as the retina.
One limitation of short-time series used in this study is the small detectable frequency range (i.e., up to 28 Hz), set by the B-scan time interval.Blood flow produces a characteristic peak typically in the kilohertz-frequency range [31], well outside the range accessed by our method.Increasing the OCT scanning speed to shorten the B-scan interval increases the maximum detectable frequency, but is limited by the fast-axis scanning speed of our system.Alternatively, reducing the sample density (i.e., the number of A-scans) in the B-scan can also improve this, but only slightly, as adequate sample density is needed to represent the vessels in the projection images.Thus, short-time series OCTA, as proposed in this study, is not suitable for velocimetry.In future, faster scanning systems may help address this [41].

Conclusions
We have presented a new method, short-time series OCTA, to perform imaging of tissue microvasculature in vivo.Our method uses the flow-induced signature in the frequency domain via Fourier transform of the time series of the OCT signal in five B-scans from the same lateral location.The angiography signal is computed as the average magnitude of the non-zero (high-pass) frequency components, clearly differentiating blood vessels and static tissue, as demonstrated in a flow phantom and in human skin in vivo.Weighting of the angiography signal by the inverse of the mean OCT signal demonstrated improved detection of blood vessels.We determined the practical minimum number of co-located B-scans needed for analysis, as well as confirming that implementation on the OCT intensity signal is practically superior to implementation on the complex signal without motion correction.The imaging performance of short-time series OCTA was assessed by comparison to the commonly used speckle decorrelation and speckle variance methods, showing consistently superior results, evidenced by improved visualization, especially for small vessels, and increased vasculature density of the human cutaneous microvascular network.Given the requirement for only normal sampling parameters, in combination with fast computing, this method should find application as an efficient and straightforward approach to OCTA in various biomedical applications to visualize the microvasculature in vivo.

Funding
Australian Research Council; National Health and Medical Research Council.
Fig. 1 (a) an magni (right) region from static m the elevated .An example own in Fig. 1(b T intensity sign e are determin llary (red arrow The insets in gnal (right).theflow phantom d).Insets show the ngiography images he flow and static spectra originating as, those from the ysical length.
Fig. 2 locate magni region orang axis.I the ratThe avera and for the st does the diffe the ratio betw it reaches a p Fig. 3 and af corres and m3.4 IntensityAs with other OCT signal (i Fig. 4 before projec3.5 CompariTo further ass used intensity variance, appl

Figure 5
Figure 5 s surface to 300 the middle ro (above) and s visualization (Fig. 5(a)) an improved con the enhanced outlined tissu (h), obtained segments are The cross-sec and (i).Such estimated acc the short time decorrelation method result 5. OCTA compari -time series (d) an ified in (b), (e) an urface to 300 µm i the cross-sectiona sponding to (c), (f) shows an exam 0 µm in depth ow to allow r speckle varian of the vessel nd speckle va ntrast of the blo connectivity e regions in Fi using speckle more clearly ctional views o improvement i curacy of appro e-series image and speckle v ts from the imp ison case 1. Proje nd speckle varian nd (h), respectivel in depth.The arro al views in (c), ( ) and (i).White sca mple from forea h.Vessel image ready comparis nce (below).In network that ariance method ood vessels in and visibility ig.5(d), is mag e decorrelation observed, with of these marke is further quan oximately 1%.shown in Fig variance, respe proved vessel c ection of blood ve nce (g).The outlin ly.Vessels from th ows mark the sam (f) and (i), respec ale bars: 500 µm.arm skin, proje es generated b son to images n Fig. 5(d), th is comparable ds (Fig. 5(g)) the short-time of the vessels gnified in Fig. n and speckle v h a representa ed vessels are s ntified by meas .This results i .5(d), in com ctively.The h contrast, as sho essels by speckle ned regions in (a the forearm skin a me vessel segment ctively.(j) OCT Cyan scale bars: 2 ecting the bloo by the short-tim s generated by he short-time s e to that of the ).Further exa e series project s.One such ex 5(e).In comp variance, resp ative example shown, respect suring the vess in a superior a mparison to 21% higher density own in Fig. 5. decorrelation (a), a), (d) and (g) are are projected from in (b), (e) and (h) image (0-41 dB) 200 µm.od vessels from me series meth y speckle deco series method e speckle deco amination indi tion image, ob xample, taken arison to Figs. ectively, sever marked by the tively, in Figs sel area density area density of % and 20% fo in the short-tim

Fig. 7 .
Fig. 7.All thr ed of the micr at zero flow sp