Visualization of blood vessels is a fundamental task in the evaluation of the health and biological integrity of tissue. Laser speckle contrast imaging (LSCI) is a non-invasive technique to determine the blood flow in superficial or exposed vasculature. However, the high scattering of biological tissue hinders the visualization of those structures. In this paper, we propose the use of principal component analysis (PCA) in combination with LSCI to improve the visualization of deep blood vessels by selecting the most significant principal components. This analysis was applied to in vitro samples, and our results demonstrate that this approach allows for the visualization and localization of blood vessels as deep as 1000 μm.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
LSCI is a technique in which time-integrated speckle patterns generated by low-power laser irradiation (usually < 30 mW) are recorded by a CCD camera [1–3]. Time-integrated speckle pattern analysis enables the visualization of superficial blood flow in exposed or superficial vasculature. LSCI is typically used to map and quantify relative changes in blood flow in response to an intervention and/or an external stimulus. LSCI does not enable visualization of subsurface or small vasculature because of the optical scattering by static structures, such as the skull or epidermis. In order to address this problem, a variety of methods have been proposed, for example, magnetomotive LSCI  use an external alternating magnetic field to induce movement of superparamagnetic iron oxide nanoparticles introduced into the vasculature. The induced motion of the particles causes an increase in the local blood motion and consequently, a decrease in the local speckle contrasts. Another method involves blood heating due to light absorption from a short pulse laser impinging on the skin [5,6]. Usually mid-infrared detectors are used to collect the spatially-dependent infrared emission from skin. In a closely related technique, photothermal LSCI [7,8], a short pulse laser (λ = 595 nm) is used to heat up subsurface blood vessels, so the rapid heating of the molecules causes an increment in the particles motion, which traduces as a decrease in the local speckle contrast. However, in these invasive methods a pulsed laser or external agents are needed to modify the contrast to improve the visualization of blood vessels. In Refs [9,10] the authors use a non-invasive image processing method based on PCA to analyze the dynamic speckle patterns in maize seed, on apple and the drying process of a painted coin. In these works, PCA is used as a filtering process to study spatially and temporally the dynamic of the speckle patterns. It is worth to notice that the dynamic of those three examples samples is much slower than the dynamic of a blood vessel. Recently, Li et al  propose an eigen-decomposition filtering approach to observe in vivo, an angiography of mouse ear pinna, a very similar approach is reported by Wang et al , where, a real-time full-field optical angiography utilizing principal component analysis is proposed, however, the effect of the exposure time (T) and the thickness of the overlaying layer on the blood vessels was not studied.
In this work, we propose the use of PCA in combination with LSCI to improve the visualization of subsurface blood vessels by separating out the static and dynamic components of the backscattering light. The first principal component accounts for as much of the variability in the data as possible, and each succeeding component accounts for as much of the remaining variability as possible [13,14]. By using the Guttan-Kaiser Criterion we define the grouping threshold between the principal components associated with the dynamic and static speckle coming from an in-vitro blood vessel phantom and study the effect of the exposure time and the thickness of the scattering overlaying layer on the blood vessel phantom. Our results show that by using PCA in combination with the traditional LSCI, it is possible to improve the visualization and localization of blood vessels at depths as large as 1000 µm.
2.1 Laser speckle imaging
The intensity image (raw speckle image) is mapped to the contrast one by using a spatial contrast algorithm . In short, the local contrast K is computed, typically using a sliding window of 5x5 pixels, with the equationEq. (2) depend on the ratio of the exposure time and the correlation time, they can be associated with light scattered by dynamic structures, such as the blood running through the blood vessels,
2.2 Principal components analysis
The procedure to obtain the principal components (PCs) begins with the organization of the data in a matrix of dimension M x N, where M represents the number of observations and N represents the number of each set or group of observations that are correlated with each other.
The final step involves the calculation of the uncorrelated principal components scores by
3. Materials and methods
3.1 LSCI optical system
The optical setup of a LSCI system is quite simple: a coherent light source (He-Ne laser, λ = 632.8 nm, 30 mW, Thorlabs) is employed to irradiate the sample under study. A plano convex lens and an optical diffuser are employed to uniformly irradiate an area of approximately 2 cm of diameter. The resulting speckle pattern is imaged into a monochromatic CCD camera (Retiga 2000R, Qimaging, Canada) equipped with a macro lens. The field of view was set to an area of approximately 1 cm of diameter. The image integration time was varied from 0.06 to 40 ms. Video images acquired at 30 frames per second were transferred from the camera to a PC. Custom software written in MATLAB (R2017b, The MathWorks, Inc., Natick, MA) was developed to acquire and process the images. To perform PCA, 30 successive raw speckle images (RSI) were processed. We use 30 images because it is typical number employed in other LSCI algorithms .
3.2 Skin-simulating phantom
The skin-simulating phantom consists of a silicone block containing appropriate concentrations (2 mg/mL) of titanium dioxide (TiO2) microparticles to mimic the scattering coefficient of biological tissue at visible and near-infrared wavelengths . We embedded a glass capillary tube, with an inner diameter of 550 µm, at the surface of the block. A second thin silicone skin-simulating phantom of varying thicknesses (190, 311, 510, and 1000 µm) is placed on top of the glass capillary tube to mimic the overlying tissue. The thin phantoms were fabricated using TiO2 powder (2 mg/mL) and freeze dried coffee powder (10 mg/mL) to simulate the optical scattering and absorption coefficients, respectively, of the epidermis . By using different thicknesses of the overlying layer, we effectively model the static component of speckle on the CCD sensor. We collected data from the phantom with the following configurations: no top layer present (i.e., capillary tube at surface of phantom) and with overlying top-layer of thicknesses δ = 190, 311, 510, and 1000 µm. We used a syringe pump (Harvard Apparatus, Holliston, Massachusetts) to infuse intralipid at 3% in water as a blood substitute into the glass capillary at speed of 4 mm/s, representative of flow in arterioles and venules .
3.3 PCA study
Following the procedure described in Section 2.2, the first step consists in converting each image into a column vector. This can be done by taking each column of the image and arranging one below each other. Each image (now a column vector) is arranged one next to each other to form the matrix (Eq. (5)) as shown in Fig. 1. In our case, each collected images have a size of 640 x 480 pixels, so the column vector dimension is of 1 x 307200 pixels. Finally, the matrix dimension is 30 x 307200 pixels/elements.
4. Results and discussion
4.1 PCA study
In order to analyze the information encoded in the PCs, they were grouped into two (A and B) according to the Guttan-Kaiser Criterion . This criterion is based on the weight of each PC. The group A consist of all PCs whose eigenvalues satisfy , where is the mean value of all the eigenvalues (Eq. (12)). Group B consists of all the PCs whose eigenvalues satisfy that . By definition, the sum of all eigenvalues is 1.
Figure 2 shows the eigenvalues vs the number of PCs, the thick red horizontal line establishes the separation among components according to the Guttan-Kaiser Criterion for three different thickness of the overlaying layer δ at three different exposures times T. We can observe that most of the times only the first PC belngs to group A. Only in 4% of the experiments reported here, the group A contains the first two PC. For a given δ and T values, the first eigenvalue has a much higher value than the rest of them. This fact is more evident for long exposure times and thicker top layers. It means that for long exposure times (i.e. 10.929 ms) and/or thick top layer (i.e. 1000 µm) the set of images are quite similar and therefore the variance of the raw speckle images is concentrated in the first PC and the filtering process is not efficient. On the other hand, for short exposure times and/or thin top layers, the variance is distributed more homogeneously among all PCs and therefore the filtering process is more efficient.
This fact is more evident in Fig. 3, were the sum of eigenvalues in the group A vs the exposure time for different overlaying layer thickness is shown. For long exposure times the sum of eigenvalues approximates to 1, and therefore most of the information is concentrated in the group A, as mentioned above the filtering process is not efficient in this case. A similar trend, but less dramatic, is found as the top layer becomes thicker.
Once the groups A, B and C were generated, their corresponding speckle images (Fig. 4) were recovered using Eq. (11). The contrast images were obtained from each group using a sliding window of 5x5 pixels. It can be observed that contrast values for group A are significantly higher in the area outside the capillary (static region), while in group B the higher values are located in the capillary zone (dynamic region). It means that the static structure is highlighted in group A and the dynamic one in group B. Group C allows us to compare PCA + LSCI technique vs LSCI alone. The same behavior shown in Fig. 4 was observed for short (< 1 ms) exposure times and thickness of the overlying scattering layers employed in this work.
In order to study contrast as function of the exposure time and different overlaying layers, we selected a region of interest (ROI) of 10 x 100 pixels centered on the dynamic region shown in Fig. 5. The contrast values in this region were averaged and a representative single value was obtained. This procedure was repeated for groups A, B, and C.
Figure 6(a) shows the experimental curve of (symbols) as function of the exposure time, obtained from the ROI of group C, for different For a given value, is maximum for short exposure time (T) and decays for larger T, reaching a constant value () for longest exposure times, as predicted by Eq. (3). The observed offset is scaled up with the thickness of the scattering layer. The continuous lines correspond to theoretical fits to Eq. (2). Note the good agreement between theory and experiment.
Figure 6(b) shows a similar analysis for group A. Interestingly, for a given value, remains constant for all the exposure times, i.e., the contribution from the static region of the ROI to the global contrast is independent of T and increases with the thickness of the layer, as expected. Moreover, the asymptotic behavior of seems to be a continuation from those of Fig. 6(a) and the offsets in both figures match almost perfectly. Thus, our analysis suggest that the first PCs contains the static component of and more importantly, we can single out the static component from the dynamic one.
The dynamic scattering contrast is obtained by subtracting contrast of the group A from the one of group C. Figure 7, shows as function of T for different . The symbols correspond to the experimental results and the dashed line corresponds to fittings to Eq. (4). Interestingly the contrast values obtained are almost independent of the thickness of the scattering layers for all the exposure times employed in this study. This result is fundamental to calculate the speckle flow index () in a more accurate way since the influence of the static speckle has been filtered out.
Figure 8 shows a comparison between the traditional LSCI (group C) and PCA-based dynamic contrast (group C-A). From left to right, the first column shows images of for the standard LSCI (group C) for different top layers. As the thickness δ of scattering layer increases, it becomes increasingly harder to identify the blood vessel, being almost impossible to distinguish at δ = 1000 µm. In the second column, we used the segmentation technique known as K-mean, commonly used in digital image processing  to segment the images of group C (labeled as Seg(C)). Basically, in the segmentation technique the pixels of an image are clustered according to its values. Even in this case, it is hard to distinguish the blood vessel phantom at δ = 1000 µm. The third column corresponds to the dynamic contrast obtained from the difference between the traditional LSCI images (Group C) and the images recovered with the first PC (Group A), labeled as C-A images. Its corresponding segmentation (Seg(C-A)) is shown in the last column. An obvious improvement on the visualization of the blood vessel phantom is obtained.
4.2 Kurtosis analysis
Although segmentation helps to visualize deep blood vessels, the neighboring region corresponding to the transition between static and dynamic regions is ill defined. In order to solve this issue, we applied kurtosis analysis to the segmented images. Kurtosis has been successfully used as a parameter to estimate the limits of the blood vessel in LSCI images . Kurtosis analysis is used as a post-processing step to better define the dynamic region (blood vessels) in both C and C-A images (Fig. 8). Figure 9 show the automatic blood vessel location obtained from the kurtosis analysis for δ = 0, 510 and 1000 µm. For more details on the kurtosis analysis see . When the analysis is applied to Seg(C) and δ = 0, the automatic location of the blood vessel (round colored markers) approach to its actual position (parallel vertical black lines), as shown in Fig. 9(a). Nevertheless, as the depth increases (δ>300 µm) location moves away from the reference because the noise around the vessel increases, causing that the vessel profile to have a normal distribution and therefore a lower kurtosis. Thus, for an accurate blood vessel location, it is necessary that the profile has high kurtosis and the standard LSCI algorithm (Seg(C) image) at high depth does not provide such condition. Figure 8(b) shows the profiles obtained from the Seg(C-A) images at the same depths than the previous plot. These profiles keep a high kurtosis level even when depth increases, and the automatic location of the blood vessel is closer to the actual location. This fact demonstrates the effectiveness of PCA as a tool to separate static and dynamic regions from LSCI images.
These results are summarized in Table 1 show that even segmenting image C, the estimated the error for the location of the deepest vessel is around 140%. On the other hand, when we combine PCA and LSCI, the vessel location improves, and the estimated width of the vessel is more accurate, the error is one order of magnitude lower compared with the standard LSCI algorithm.
In this work, we demonstrate that PCA applied to LSCI allow us to separate out the static component from the dynamic component in the raw speckle images achieving visualization of blood vessels as deep as 1000 µm. In addition, employing kurtosis on the dynamic region (such as flow in the capillary) we demonstrated a more accurate estimation of the actual vessel width. It is important to mention that our proposal works well when the variance between the raw speckle images is enough to be distributed among all the PCA, for example when the dynamic of the sample is relatively low or for short exposure times. Otherwise the first component attracts most of the information compared with the rest of the component and the filtering process fails.
NSF (1545852); Mexican Council of Science and Technology (CONACyT) (OISE:PIRE-SOMBRERO)/CONACyT 251992, grant 261148).
T. Spezzia-Mazzocco would like to thank the program Cátedras-CONACyT. R. Chiu would like to thank CONACyT for their support for sabbatical leave. The authors would like to thank Red de Biofotónica CONACYT number F0003-2018-03-294910.
The authors declare that there are no conflicts of interest related to this article.
3. K. R. Forrester, C. Stewart, J. Tulip, C. Leonard, and R. C. Bray, “Comparison of laser speckle and laser Doppler perfusion imaging: Measurement in human skin and rabbit articular tissue,” Med. Biol. Eng. Comput. 40(6), 687–697 (2002). [CrossRef] [PubMed]
5. B. Li, B. Majaron, J. A. Viator, T. E. Milner, Z. Chen, Y. Zhao, H. Ren, and J. S. Nelson, “Accurate measurement of blood vessel depth in port wine stained human skin in vivo using pulsed photothermal radiometry,” J. Biomed. Opt. 9(5), 961–966 (2004). [CrossRef] [PubMed]
6. S. A. Prahl, I. A. Vitkin, U. Bruggemann, B. C. Wilson, and R. R. Anderson, “Determination of optical properties of turbid media using pulsed photothermal radiometry,” Phys. Med. Biol. 37(6), 1203–1217 (1992). [CrossRef] [PubMed]
9. K. M. Ribeiro, R. A. Braga, G. W. Horgan, D. D. Ferreira, and T. Safadi, “Principal component analysis in the spectral analysis of the dynamic laser speckle patterns,” J. Eur. Opt. Soc. 9, 14009 (2014). [CrossRef]
10. J. M. López-Alonso, E. Grumel, N. L. Cap, M. Trivi, H. Rabal, and J. Alda, “Characterization of spatial–temporal patterns in dynamic speckle sequences using principal component analysis,” Opt. Eng. 55(12), 121705 (2016). [CrossRef]
12. M. Wang, C. Guan, W. Mao, H. Xiong, H. Tan, D. Hang, and Y. Zeng, “Real-time full-field optical angiography utilizing principal component analysis,” Opt. Lett. 43(11), 2559–2562 (2018). [CrossRef] [PubMed]
13. S. Jung, A. Sen, and J. S. Marron, “Boundary behavior in High Dimension, Low Sample Size asymptotics of PCA,” J. Multivariate Anal. 109, 190–203 (2012). [CrossRef]
14. H. Abdi and L. J. Williams, “Principal component analysis,” Wiley Interdiscip. Rev. Comput. Stat. 2(4), 433–459 (2010). [CrossRef]
15. A. F. Fercher and J. D. Briers, “Flow visualization by means of single-exposure speckle photography,” Opt. Commun. 37(5), 326–330 (1981). [CrossRef]
17. J. C. Ramirez-San-Juan, R. Ramos-Garcia, G. Martinez-Niconoff, and B. Choi, “Simple correction factor for laser speckle imaging of flow dynamics,” Opt. Lett. 39(3), 678–681 (2014). [CrossRef] [PubMed]
18. Y. Yang, A. Lü, W. Li, and Z. Qian, “Microfluidic-based laser speckle contrast imaging of erythrocyte flow and magnetic nanoparticle retention in blood,” AIP Adv. 9(1), 150031 (2019). [CrossRef]
19. F. Ayers, A. Grant, D. Kuo, D. J. Cuccia, and A. J. Durkin, “Fabrication and characterization of silicone-based tissue phantoms with tunable optical properties in the visible and near infrared domain,” Proc. SPIE Design and, 687007 (2008). [CrossRef]
20. V. S. Kalambur, H. Mahaseth, J. C. Bischof, M. C. Kielbik, T. E. Welch, A. Vilbäck, D. J. Swanlund, R. P. Hebbel, J. D. Belcher, and G. M. Vercellotti, “Microvascular blood flow and stasis in transgenic sickle mice: utility of a dorsal skin fold chamber for intravital microscopy,” Am. J. Hematol. 77(2), 117–125 (2004). [CrossRef] [PubMed]
21. H. F. Kaiser, “a Note on Guttman’S Lower Bound for the Number of Common Factors,” Br. J. Stat. Psychol. 14(1), 1–2 (1961). [CrossRef]
22. R. C. Gonzalez, R. E. Woods, and B. R. Masters, “Digital Image Processing, Third Edition,” J. Biomed. Opt. 14, 029901 (2009).
23. H. Peregrina-Barreto, E. Perez-Corona, J. Rangel-Magdaleno, R. Ramos-Garcia, R. Chiu, and J. C. Ramirez-San-Juan, “Use of kurtosis for locating deep blood vessels in raw speckle imaging using a homogeneity representation,” J. Biomed. Opt. 22(6), 066004 (2017). [CrossRef] [PubMed]