Visualization of in vitro deep blood vessels using principal component analysis based laser speckle imaging

: 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.


Introduction
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][2][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 [4] 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 noninvasive 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 [11] 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 [12], 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.

Laser speckle imaging
The intensity image (raw speckle image) is mapped to the contrast one by using a spatial contrast algorithm [15]. In short, the local contrast K is computed, typically using a sliding window of 5x5 pixels, with the equation where σ is the standard deviation and I < > is the mean intensity of the pixels within the sliding window. The local contrast value is assigned to the central pixel. Then, the window is displaced by one pixel and the contrast is calculated again and then the process is repeated until the whole image is processed. The resulting image is known as the contrast image.
The empirical contrast Eq. (1) can be expressed [16] as function of the correlation time c τ of the backscattered light from the sample and the exposure time T of the CCD camera as [17]: , ρ is the fraction of the dynamically scattered light and β is a correction factor that depends on the ratio of speckle and pixel size. When 1 x  , the contrast reaches an asymptotic value ( S K ) given by: i.e. a time-independent contrast S K value which can be associated to light scattered by static structures, such as the tissue that surrounds the blood vessels. Since the first two terms on Eq.
(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, so, in general, the scattered light contains both static and dynamic components. It may be possible that one prevails over the other like in blood vessels. Thus, the objective of this work is to use PCA in combination with LSCI to separate out 2 S K from 2 D K and therefore improve the visualization of deep blood vessels.

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. 11 12 1 The second step consists in extracting the mean of each set N. This step avoids the statistical dispersion of the data set; this process called centralization is represented by where i Φ is the data column vector centralized around the mean, i γ represents the N-th observation, and ( ) i μ γ is the mean of the observations column vector, and they can be calculated by The third step comprises the calculation of the covariance matrix by mean of the following equation where C represents the covariance matrix, Φ and T Φ the data matrix with zero mean and its transpose, respectively. The resulting covariance matrix is real and symmetric. This allows us the factorization of the matrix C in the following way where Λ is an M x M matrix of orthonormal eigenvectors and λ is the diagonal matrix of eigenvalues. The final step involves the calculation of the uncorrelated principal components scores by As for any other transformation, it is possible to compute the inverse transformation of the PC into the original data using the equation The original data can be recovered, only if the overall set of eigenvectors is used in the inverse transformation. On the other hand, if some eigenvectors are used in the inverse transformation, it is still possible to construct a new data set (approximation) where some undesirable characteristic can be discriminated or removed from the original data.

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 [18].

Skin-simulating phantom
The skin-simulating phantom consists of a silicone block containing appropriate concentrations (2 mg/mL) of titanium dioxide (TiO 2 ) microparticles to mimic the scattering coefficient of biological tissue at visible and near-infrared wavelengths [19]. 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 TiO 2 powder (2 mg/mL) and freeze dried coffee powder (10 mg/mL) to simulate the optical scattering and absorption coefficients, respectively, of the epidermis [20]. 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 [20].

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.   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 hi compare PCA observed for employed in t     Figure 6(b) shows a similar analysis for group A. Interestingly, for a given δ value, 2 K 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 2 K 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 2 K and more importantly, we can single out the static component from the dynamic one.  shows a comp trast (group C-LSCI (group C becomes increa h at δ = 1000 mean, common eled as Seg(C)) according to om at δ = 100 m the differenc ith the first (Seg (C-A)

Conclusio
In this work, component fro analysis gmentation he g to the transit ue, we applied used as a para s analysis is us s) in both C an ined from the k nalysis see [21 he blood vessel k lines), as sho ves away from the vessel prof accurate blood dard LSCI alg gure 8(b) show the previous p d the automatic rates the effect mages. 9. a) Location of sults are summ error for the l mbine PCA and more accurate CI algorithm. 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.