20 August 2014 Histogram flow mapping with optical coherence tomography for <italic<in vivo</italic< skin angiography of hereditary hemorrhagic telangiectasia
Author Affiliations +
J. of Biomedical Optics, 19(8), 086015 (2014). doi:10.1117/1.JBO.19.8.086015
Speckle statistics of flowing scatterers have been well documented in the literature. Speckle variance optical coherence tomography exploits the large variance values of intensity changes in time caused mainly by the random backscattering of light resulting from translational activity of red blood cells to map out the microvascular networks. A method to map out the microvasculature malformation of skin based on the time-domain histograms of individual pixels is presented with results obtained from both normal skin and skin containing vascular malformation. Results demonstrated that this method can potentially map out deeper blood vessels and enhance the visualization of microvasculature in low signal regions, while being resistant against motion (e.g., patient tremor or internal reflex movements). The overall results are manifested as more uniform en face projection maps of microvessels. Potential applications include clinical imaging of skin vascular abnormalities and wide-field skin angiography for the study of complex vascular networks.
Cheng, Mariampillai, Lee, Vuong, Luk, Ramjist, Curtis, Jakubovic, Kertes, Letarte, Faughnan, HHT Investigator Group, and Yang: Histogram flow mapping with optical coherence tomography for in vivo skin angiography of hereditary hemorrhagic telangiectasia



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-scans1 or B-frames2 are calculated which are directly correlated to flow velocities of optical scatterers. More sophisticated phase-sensitive techniques, such as optical microangiography (OMAG)3 and phase variance optical frequency domain imaging,4 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),7 which maps out the decorrelation across B-mode frames or even A-scans8 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.9 Various techniques have been developed to circumvent this drawback, such as logarithmic speckle transformation,9 split-spectrum decorrelation,10 double-decorrelation mapping,11 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.



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)


where A is the speckle amplitude and σ is the standard deviation of the underlying quadrature. When a pixel is partially populated by static tissues and flowing scatterers, as in the case of smaller vessels such as capillaries, the amplitude distribution follows the Rician distribution, which is the generalization of the Rayleigh distribution by the introduction of the noncentrality parameter in Eq. (2)


where ν is the noncentrality parameter and I0 is the zeroth-order modified Bessel function of the first kind. In general, static tissues are corrupted by complex additive Gaussian noise, with a high specular reflection component and thus large ν, and the corresponding amplitude statistics would converge from Rician to Gaussian distributions (as νσ) in Eq. (3):13



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,9 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) distributions14










where Ψ2 is one of the Horn’s confluent hypergeometric functions, μ is the mean of the Gaussian distribution, and Γ is the Gamma function. Notice that Eqs. (4) and (5) have an exponential decay nature while Eq. (6) also has an oscillatory nature. In addition, from Eqs. (4), (5), and (7), the effect of the exp[(ν2/2σ2)] term in Eq. (5) can be seen to be partially cancelled by the hypergeometric function. Figure 1(a) shows a typical normalized Im[M(ω)Rayleigh] and Fig. 1(b) shows a typical normalized Im[M(ω)Gaussian] in OCT data, whereas the normalized Im[M(ω)Rician] takes an intermediate shape compared when to them. By integrating the histogram as described in the Method section, Rayleigh and Rician distributions can be distinguished from Gaussian distributions.

Fig. 1

(a) Im(M) of Gaussian distribution; (b) Im(M) of Rayleigh distribution.




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×3mm2 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×30pixels 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(ω). The imaginary part Im[M(ω)] was then taken out and the following contrast mechanism S(x,y) was computed: S(x,y)=0aIm[M(x,y,ω)]dω, 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.



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 TiO2 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μ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μm.




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.15 Processing can also be dramatically speeded up by the use of massively parallel GPU computing, 16 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.17 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,18 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.19 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.20 Further research of our proposed method may enable high fidelity quantitative perfusion microangiography of complex malformed microvessel networks.



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.


This work is supported by the National Science and Engineering Council of Canada (NSERC), Canadian Institute of Health Research, and the Canada Foundation for Innovation (CFI). The HHT study is part of a collaboration with the Brain Vascular Malformation Consortium (BVMC) HHT Investigator Group and with support from a grant provided by the National Institutes of Health (NIH) Rare Disease Clinical Research Network (RDCRN). The Brain Vascular Malformation Consortium (BVMC; U54NS065705) is a part of the National Institutes of Health (NIH) Rare Disease Clinical Research Network (RDCRN), supported through collaboration between the NIH Office of Rare Diseases Research (ORDR) at the National Center for Advancing Translational Science (NCATS), and the National Institute of Neurological Disorders and Stroke (NINDS). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

The Brain Vascular Malformation Consortium (BVMC) is an integrated group of academic medical centers, patient support organizations, and clinical research resources focused on clinical research on brain vascular malformations and improving care for patients with hereditary hemorrhagic telangiectasia (HHT), Sturge-Weber syndrome, and familial cerebral cavernous malformations. The BVMC HHT Investigator Group consists of: Murali Chakinala, Marie E. Faughnan, James R. Gossage, Katharine Henderson, Vivek Iyer, Raj Kasthuri, Helen Kim, Timo Krings, Michael Lawton, Doris Lin, Johannes Jurgen Mager, Justin McWilliams, Jamie McDonald, Ludmila Pawlikowska, Jeffrey Pollak, Felix Ratjen, Karen Swanson, Karel terBrugge, Dilini Vethanayagam, Andrew White, Robert I. White Jr., Pearce Wilcox, and William L. Young.



V. X. D. Yanget al., “High speed, wide velocity dynamic range Doppler optical coherence tomography (Part III): in vivo endoscopic imaging of blood flow in the rat and human gastrointestinal tracts,” Opt. Express 11(19), 2416–2424 (2003).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.11.002416Google Scholar


B. Braafet al., “Angiography of the retina and the choroid with phase-resolved OCT using interval-optimized backstitched B-scans,” Opt. Express 20(18), 20516–20534 (2012).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.20.020516Google Scholar


L. AnJ. QinR. K. Wang, “Ultrahigh sensitive optical microangiography for in vivo imaging of microcirculations within human skin tissue beds,” Opt. Express 18(8), 8220–8228 (2010).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.18.008220Google Scholar


B. J. Vakocet al., “Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging,” Nat. Med. 15(10), 1219–1223 (2009).1078-8956http://dx.doi.org/10.1038/nm.1971Google Scholar


A. Mariampillaiet al., “Speckle variance detection of microvasculature using swept-source optical coherence tomography,” Opt. Lett. 33(13), 1530–1532 (2008).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.33.001530Google Scholar


A. Mariampillaiet al., “Optimized speckle variance OCT imaging of microvasculature,” Opt. Lett. 35(8), 1257–1259 (2010).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.35.001257Google Scholar


E. JonathanJ. EnfieldM. J. Leahy, “Correlation mapping method for generating microcirculation morphology from optical coherence tomography (OCT) intensity images,” J. Biophoton. 4(9), 583–587 (2011).JBOIBX1864-063Xhttp://dx.doi.org/10.1002/jbio.v4.5Google Scholar


X. Liuet al., “Quantitative transverse flow measurement using optical coherence tomography speckle decorrelation analysis,” Opt. Lett. 38(5), 805–807 (2013).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.38.000805Google Scholar


R. MotaghiannezamS. Fraser, “Logarithmic intensity and speckle-based motion contrast methods for human retinal vasculature visualization using swept source optical coherence tomography,” Biomed. Opt. Express 3(3), 503–521 (2012).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.3.000503Google Scholar


Y. Jiaet al., “Split-spectrum amplitude-decorrelation angiography with optical coherence tomography,” Opt. Express 20(4), 4710–4725 (2012).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.20.004710Google Scholar


A. DoroninI. Meglinski, “Imaging of subcutaneous microcirculation vascular network by double correlation optical coherence tomography,” Laser Photon. Rev. 7(5), 797–800, (2013).1863-8880http://dx.doi.org/10.1002/lpor.2013.7.issue-5Google Scholar


J. W. Goodman, Speckle Phenomena in Optics, Roberts and Company, Greenwood Village, Colorado (2007).Google Scholar


H. GudbjartssonS. Patz, “The Rician distribution of noisy MRI data,” Magn. Reson. Med. 34(6), 910–914 (1995).MRMEEN0740-3194http://dx.doi.org/10.1002/(ISSN)1522-2594Google Scholar


X. LiuL. Hanzo, “A unified exact BER performance analysis of asynchronous DS-CDMA systems using BPSK modulation over fading channels,” IEEE Trans. Wireless Commun. 6(10), 3504–3509 (2007).WCDMA1536-1276http://dx.doi.org/10.1109/TWC.2007.060017Google Scholar


W. Wieseret al., “Extended coherence length megahertz FDML and its application for anterior segment imaging,” Biomed. Opt. Express 3(10), 2647–2657 (2012).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.3.002647Google Scholar


K. K. C. Leeet al., “Real-time speckle variance swept-source optical coherence tomography using a graphics processing unit,” Biomed. Opt. Express 3(7), 1557–1564 (2012).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.3.001557Google Scholar


V. Kalchenkoet al., “Label free in vivo laser speckle imaging of blood and lymph vessels,” J. Biomed. Opt. 17(5), 050502, (2012).JBOPFO1083-3668http://dx.doi.org/10.1117/1.JBO.17.5.050502Google Scholar


V. Kalchenkoet al., “In vivo characterization of tumor and tumor vascular network using multi-modal imaging approach,” J. Biophoton. 4(9), 645–649 (2011).JBOIBX1864-063Xhttp://dx.doi.org/10.1002/jbio.201100033Google Scholar


I. M. BravermanA. KehB. S. Jacobson, “Ultrastructure and three-dimensional organization of the Telangiectases of hereditary hemorrhagic Telangiectasia,” J. Investig. Dermatol. 95(4), 422–427 (1990).JIDEAE0022-202Xhttp://dx.doi.org/10.1111/jid.1990.95.issue-4Google Scholar


J. BartonS. Stromski, “Flow measurement without phase information in optical coherence tomography images,” Opt. Express 13(14), 5234–5239 (2005).OPEXFF1094-4087http://dx.doi.org/10.1364/OPEX.13.005234Google Scholar

Biographies of the authors not available.

Kyle H. Cheng, Adrian Mariampillai, Kenneth K. Lee, Barry Vuong, Timothy W. Luk, Joel Ramjist, M. Anne Curtis, Henry Jakubovic, Peter Kertes, Michelle Letarte, Marie E. Faughnan, Victor X. Yang, "Histogram flow mapping with optical coherence tomography for <italic<in vivo</italic< skin angiography of hereditary hemorrhagic telangiectasia," Journal of Biomedical Optics 19(8), 086015 (20 August 2014). http://dx.doi.org/10.1117/1.JBO.19.8.086015

Optical coherence tomography




In vivo imaging


Image processing


Back to Top