1.

Introduction

Optical coherence tomography (OCT) serves as an emerging optical imaging modality that provides high resolution and deep tissue penetration cross-sectional images of biological tissues. One major functional utility of OCT is vascular imaging. Traditional Doppler OCT (DOCT) vascular imaging involves techniques which are phase sensitive. Either phase changes between A-scans¹ or B-frames² are calculated which are directly correlated to flow velocities of optical scatterers. More sophisticated phase-sensitive techniques, such as optical microangiography (OMAG)³ and phase variance optical frequency domain imaging,⁴ are developed to improve the microcirculation imaging down to the arterioles/venules level. However, aliasing and A-scan rates requirements often have complicated straightforward implementation in terms of instrumentation. To circumvent strict hardware requirements, speckle-based flow imaging techniques come into play.

Speckle, caused by the interference of randomly back-scattered photons from biological tissues, is an inherent phenomenon in coherent imaging modalities. Various optical techniques exploit the intrinsic speckle modulation of flowing scatterers versus the relatively stable speckle patterns of static tissues to map out the microvasculature in biological tissues. Speckle variance OCT (svOCT) calculates the variance estimator of the statistics of the time-varying voxel amplitude to differentiate the significant speckle modulation from the surrounding bulk tissues.^5,6 Correlation mapping OCT (cmOCT),⁷ which maps out the decorrelation across B-mode frames or even A-scans⁸ to determine the vascular regions is an alternative way to calculate the flow-induced speckle modulation.

Despite being powerful and simple, one major drawback of svOCT is that the technique is signal strength dependent.⁹ Various techniques have been developed to circumvent this drawback, such as logarithmic speckle transformation,⁹ split-spectrum decorrelation,¹⁰ double-decorrelation mapping,¹¹ and so on. Decorrelation techniques, however, require the application of subjectively determined masks. In this paper, with the emergence of high speed OCT platforms, a new flow mapping method based on the time-domain histogram of individual pixels (not to be confused with histogram equalization) is presented. Results demonstrated that this new method can potentially enhance the visualization of deeper blood vessels while being less susceptible to motion and can potentially be employed in skin OCT angiographic study of complex vascular irregularities.

2.

Theory

Although speckle formation is a random event, for individual pixels the overall resultant temporal speckle modulation in an OCT B-mode series follows specific distributions.^9,12 The speckle amplitude distribution of a voxel completely occupied by fluid flow, as in the case of larger arterioles and venules, follows the Rayleigh distribution in Eq. (1)

To distinguish among these distributions, traditional svOCT calculates the variance estimator to map out the vasculature. Variance for the Rician distribution, however, is dependent on the noncentrality parameter,⁹ which is one of the contributing factors for svOCT to detect the microvessels. In the Rayleigh regime or low signal regions, however, svOCT signals are substantially weaker. Moreover, svOCT is also sensitive to specular reflection at interfaces. These deficiencies result in nonuniform vessel detection, as is demonstrated in the Results section.

To address these problems, consider normalizing the amplitude axis of the time-domain histogram of each pixel by the minimum and maximum amplitudes with an identical number of bins for every pixel. In this scaleless domain, amplitude distributions can be characterized mainly by their shape but not by the signal strength. To distinguish between different distributions, consider the imaginary parts of the characteristic functions of Rayleigh Eq. (4), Rician Eq. (5), and Gaussian Eq. (6) distributions¹⁴

Fig. 1

(a) $\mathrm{Im}(M)$ of Gaussian distribution; (b) $\mathrm{Im}(M)$ of Rayleigh distribution.

3.

Method

To test this hypothesis, a phantom consisting of both solid and liquid phases of Agrose with 2% titanium dioxide was prepared. The liquid phase was stirred to allow scatterer flow. A Thorlabs (Newton) OCT system in the 1300 nm range with 16 kHz A-scan rate and 100 nm bandwidth was used to acquire the B-mode sequence of the phantom over a 1-mm width. For in vivo imaging, human dorsal hand skin regions were imaged by a 1300 nm spinning polygon laser system with a 108 kHz A-scan rate and a 95-nm bandwidth over a $3\times 3{\mathrm{mm}}^2$ area. Patients’ skin was gently pressed against a microscopic slide to minimize the macromotion while real-time speckle variance was observed to make sure there was the least possible disturbance to blood flow. All procedures were approved by the St. Michael’s Hospital (Toronto, ON) Research Ethics Review Board.

After image acquisition, the data were imported into MATLAB® for postprocessing. In each B-frame, a $30\times 30\text{pixels}$ noise area was first selected, which in essence can be fixed at the lower part of the image which always consists of noise. The root-mean-square value of this noise area was subtracted from each B-frame and thresholded above zero. $N=15$ frames were then used to generate a histogram normalized by the minimum and maximum with 15 bins. A 1024 point inverse fast Fourier transform (iFFT) was then performed on the histogram to obtain the characteristic functions $M(\omega )$. The imaginary part $\mathrm{Im}[M(\omega )]$ was then taken out and the following contrast mechanism $S(x,y)$ was computed: $S(x,y)={\int }_0^a\mathrm{Im}[M(x,y,\omega )]\mathrm{d}\omega$, where $a$ was empirically chosen to be 150 to avoid high frequency noises resulting from interpolation performed in 1024 point iFFT and $(x,y)$ represents the spatial location of a B-scan cross section. The resultant map was filtered by a circular disk averaging filter of a radius of 3 pixel and finally thresholded by subtracting the median of the upper half of the map.

4.

Results

Figure 2(a) shows the structural image of the Agarose phantom. The left side was the solid phase and the right side was the randomly flowing fluid phase, indistinguishable from the structural image. Figure 2(b) shows the corresponding svOCT image ($N=4$) and Fig. 2(c) shows the image processed by our new proposed algorithm. Figure 2(d) shows the plot along the yellow line in Fig. 2(b) (svOCT), and Fig. 2(e) shows the plot along the yellow line in Fig. 2(c) (proposed algorithm). The deeper portion of the ${\mathrm{TiO}}_2$ suspension movement was clearly enhanced in Fig. 2(c). Also, the interface of the two phases (green arrow and Rician zone) in Fig. 2(b) is significantly stronger than the central area of the fluid zone, whereas in Fig. 2(c), both zones have equal brightness, demonstrating the intensity insensitivity of the new algorithm. A contrast scaling mechanism was determined from the phantom results and applied to all subsequent in vivo imaging results.

Fig. 2

(a) Structural image of the Agrose phantom (1 mm by 1 mm); (b) corresponding svOCT image, (c) image calculated by our proposed algorithm, (d) plot of the depth profile along the yellow line in (b), (e) plot of the depth profile along the yellow line in (c).

Figure 3(a) shows the structural image of an in vivo human fingernail root of a healthy volunteer, where blood vessels are abundant (9000 B-Scans with 512 A-Scans each). Figure 3(b) shows the corresponding vascular map projection using svOCT. Figure 3(c) shows the corresponding vascular map projection calculated by our proposed algorithm. The shadowing line is likely caused by protrusion of the epidermal tissue. Since current state-of-the-art vascular mapping algorithms are mostly sensitive to motion, and svOCT is sensitive to signal degradation, there are frequent appearances of bend-like artifacts in the projection map and some of the small vessels are obscured, despite high speed swept source laser systems. Figure 3(d) shows the structural image of an in vivo skin image of a vascular lesion (telangiectasia) in a patient with hereditary hemorrhagic telangiectasia (HHT) (5000 B-Scans with 512 A-Scans each), where the colored boxes indicate the different depths of the projection images created in the subsequent subfigures with matching colors (Figs. 3(e), 3(g), and 3(i) for svOCT and Figs. 3(f), 3(h), and 3(j) for our proposed algorithm). Figures 3(k) and 3(l) are the combined projections corresponding to svOCT and our proposed method, respectively. HHT is a genetic disease characterized by the presence of densely packed dilated microvessel networks, including telangiectasia of the skin and mucosa. In Fig. 3(l), there is evidence of an increased size of the blood vessels compared to normal ones [Fig. 3(c)]. Moreover, the motion artifacts and signal degradation are partially reduced in our proposed algorithm, as shown in Figs. 3(c) and 3(l), where the blurry effect of the vessels disappears. Figure 3(g) shows a histopathology of the resected telangiectasia from the same patient, which shows an area of densely packed blood vessels, which are mainly venules together with some larger arterioles as depicted in the OCT projection images. Differences between Figs. 3(l) and 3(m) can be attributed to the collapse of the microvessel channels in the ex vivo state in the absence of arterialized intravascular pressure in these abnormal channels.

Fig. 3

(a) Structural image of in vivo healthy human volunteer fingernail root (3 mm lateral by 1.25 mm axial); (b) corresponding svOCT projection of the same area; (c) corresponding projection image of human fingernail root calculated by our proposed algorithm; (d) structural image of in vivo human HHT skin telangiectasia (3 mm lateral by 1 mm axial); (e), (g), (i) projection by svOCT as indicated by the colors in (d); (f), (h), (j) projection by our proposed method as indicated by the colors in (d); (a–l) scale bar represents $300\mu \mathrm{m}$; (m) histopathology of the HHT lesion, CD34 staining demonstrating the abundant dilated microvessels. The red arrows point toward an example of dilated arterioles.

Figure 4 shows another HHT lesion (5000 B-scans with 512 A-scans each). Figure 4(a) shows the structural image of the HHT lesion and Figs. 4(b) and 4(c) show the projections of the lesion processed by our proposed algorithm with the colors indicating different depths as depicted in Fig. 4(a) at two time points 4 weeks apart. In both of the images, the surrounding normal vessels are shown feeding into the central plexus of the microvessels of the malformation. The plexus of the microvessels is clearly demonstrated at both time points, although the vessel network appears denser and more complex at the later time point.

Fig. 4

(a) Structural image of another HHT lesion (3 mm lateral by 1 mm axial), (b) corresponding projection processed by our proposed algorithm, (c) projection of the lesion 4 weeks after (b). Different colors represent the different depths as depicted in (a). Scale bar represents $300\mu \mathrm{m}$.

5.

Discussion

Our proposed method, as demonstrated in the phantom study, has the potential to map out deeper flow or dynamically moving particles located in regions of low OCT signal strength. Our HHT imaging results also show that our proposed algorithm can potentially map out the morphology of complex diseased microvascular network with higher fidelity and accuracy (Fig. 3) due to more uniform mapping in all areas by the normalized nature of the histograms. This allows us to distinguish the shape of the histograms based on their characteristic functions.

Our method also does not require the application of subjectively determined masks which can vary in a case-by-case basis. However, since background noise also follows the Rayleigh distribution, an additional step to threshold out the background noise needs to be performed.

Moreover, unlike most other vascular imaging methods in OCT, our proposed method is less sensitive to motion, since for small motions the pixels in the flowing and static regions still correspondingly exhibit the characteristic histogram features. To limit motion, dermatologic OCT data are often acquired with limited physical restriction of the imaging area as in our case. However, patients can still exhibit involuntary tremor, especially internal tissue motion, which is severe in older patients. Moreover, over-restriction will cause the patient discomfort and also prevents blood flow in some of the compressed vasculature. Sparse and rapid data acquisition is another way to mitigate motion effects; however, this may limit the data available for structural and other methods of analysis. Our method can also mitigate the dependence on sophisticated adaptive optics systems.

One drawback of the method is that more frames are required than with the traditional svOCT or correlation methods. However, these constraints can most likely be mitigated in the near future with the advent of high speed OCT systems which can reach up to 1.6 MHz A-scan rate.¹⁵ Processing can also be dramatically speeded up by the use of massively parallel GPU computing, ¹⁶ which is readily available in many computer systems.

There are alternative ways of imaging microvasculature such as laser speckle imaging. While similar to OCT in nature, laser speckle imaging is simpler and equally effective in microcirculation imaging and the measurement of decorrelation characteristics can even yield quantitative perfusion information regarding the tissue of interest, although it is limited in depth.¹⁷ Thus, while laser speckle imaging can help clinicians perform detailed diagnostics about surface vessels, our OCT technique, designed to give an accurate depiction of deeper microcirculations, can complement laser speckling imaging to enable a more comprehensive overall diagnosis of complex vascular lesions. In addition, a multimodal imaging approach can potentially be an effective way of achieving simultaneous acquisition of functional and molecular information regarding microvasculatures,¹⁸ paving the way for one-stop optical diagnosis.

There have been very few studies to date describing the microvascular structure of telangiectasia in HHT. Braverman et al. studied the HHT lesions from excised human biopsy tissues.¹⁹ With our technique, an in vivo description of the lesions can hopefully provide noninvasive insight into the structure of the malformations, and allow us to carry out a time-lapse study of the lesion dynamics as well as changes in lesion structure with therapeutic interventions.

A previous study has investigated the possibility of analyzing the speckle contrast of time-varying speckle to quantify the flow of scatterers.²⁰ Further research of our proposed method may enable high fidelity quantitative perfusion microangiography of complex malformed microvessel networks.

6.

Conclusions

In conclusion, we have proposed a vasculature mapping method for imaging skin microvasculature, which uses the time-domain histogram information of each pixel to determine the types of amplitude distributions. Since our approach is analogous to scaleless analysis, visualization enhancement of vessels in low signal regions is observed with less sensitivity to patient movement and without the application of binary masks. With these benefits, the differentiation between normal and abnormal microvessels in healthy and diseased tissues can potentially be more robust in research and clinical settings. Potential future work includes clinical imaging of skin vasculature abnormalities and wide-field skin angiography of complex vasculature, such as HHT or skin carcinoma angiogenesis and monitoring of their therapeutic progress. With further improvements to our technique, our method has the potential to be translated to analyze other areas such as endoscopic imaging.