2.1.

Texture and Co-occurrence Matrix

In computer vision, a texture is characterized by gray-level variations in an image that form small similar and repetitive regions.¹³ Thus, in order to describe a texture, it is necessary to analyze the distribution of its gray levels. While a histogram provides information about the distributions of gray levels, it does not take into account the relationship between a central pixel $p$ and its neighbors.¹⁴ This information is important to provide a rich description of the texture. The gray level CM is one of the standard methods to analyze texture in digital image processing¹³ and can be used for describing the information of RSI, as we report in an initial approach.¹⁵ Through second-order statistics, this approach studies the way in which the gray levels are distributed in an image or region and their spatial relationships.^16,17 A CM indicates the frequency at which a pair of pixels with gray values $i$ and $j$, located spatially at a distance of $d$ pixels in the direction $\theta$, occurs in an image $I$. A distance value must be chosen carefully because a large distance may result in texture overlapping; oftentimes, $d=1$ is used. The analysis direction takes values of $\theta =\{0\mathrm{deg},45\mathrm{deg},90\mathrm{deg},135\mathrm{deg}\}$; the analysis in all the directions is obtained by the addition of the four matrices.

Let $p$ be a pixel in the image $I^{H\times W}$ of $H$ rows and $W$ columns, where $p$ takes gray values in the range $[0,G-1]$, $G$ being the range of possible values (256 for an 8-bit image). For instance, if $I$ is an image with five possible gray values from 0 to 4 ($G=5$), as shown in Fig. 1(a), and the parameters $d=1$ and $\theta =45\mathrm{deg}$ are considered for its CM, then the procedure is as follows. A matrix ${\mathrm{CM}}^{G\times G}$ is generated and set to zeros; CM is always square and its dimensions depend on $G$. To assign values to the CM, each pair of pixels in $I$ is analyzed considering the parameters $d$ and $\theta$. Figure 1(b) shows that for the central pixel $p$, the parameters indicate that $q_1$ is the neighbor pixel to analyze. Often, the neighborhood analysis is made symmetrically; hence, the cooccurrence is also extended to $q_2$. Suppose we want to know the cooccurrences of a gray level $i=1$ (reference level) with a gray level $j$ (comparison level), denoted as $\mathrm{CM}(i,j)$. For $j=0$, the cooccurrences with $i$ happen when the coordinates of $p$ are $\{I(\mathrm{1,0}),I(\mathrm{1,1})\}$, therefore $\mathrm{CM}(\mathrm{1,0})=2$. For $j=1$, the cooccurrences with $i$ happen in coordinates $\{I(\mathrm{2,0}),I(\mathrm{2,1})\}$ by considering a neighborhood in $q_1$. Also, in this case, the cooccurrences happen in a symmetrical way, when $p$ has coordinates $\{I(\mathrm{1,1}),I(\mathrm{1,2})\}$ by considering a neighborhood in $q_2$; therefore, $\mathrm{CM}(\mathrm{1,1})=4$. The final CM is obtained by analyzing the cooccurrences between all the possible values of $i$ and $j$, as shown in Fig. 1(c) for this example.

Fig. 1

(a) Example of a small image $I$ with $G=5$, (b) the neighborhood for cooccurrence analysis, and (c) its resulting CM.

From the CM calculus, two important things should be noted: (a) a new representation of the gray levels in $I$ is provided by the CM, which is considered a symmetric two-dimensional histogram, but it does not characterize a texture by itself;^13,18 (b) although the amount of information in $I$ has been reduced in the CM, a large amount of data remain for texture characterization, i.e., $G\times G$ characteristics.

2.2.

Feature Extraction

In order to use the information provided by the CM, some features must be calculated. Haralick et al.¹⁷ proposed 14 textural features that can be calculated from the CM to identify difference among textures. One of the most used features is the inverse difference moment, a measure of the local homogeneity, calculated as

Fig. 2

Example of textures with (a) wide homogeneous regions (T1) and (b) fine granulated structures (T2) and a plot of their respective cooccurrence matrices.

2.3.

Skin Phantoms

These skin phantoms consist of two parts: dermis and epidermis with specific characteristics, as suggested by Saager et al.²¹ For the dermis phantom [Fig. 3(a)], we used a transparent resin with the necessary amount of titanium dioxide (${\mathrm{TiO}}_2$ $1.45\mathrm{mg}/\mathrm{mL}$) to simulate the scattering coefficient of the skin (${\mu }_{\mathrm{s}}=2{\mathrm{mm}}^{-1}$). The epidermis phantom [Fig. 3(b)] was made of silicone (polydimethylsiloxane) mixed with ${\mathrm{TiO}}_2$ powder ($2\mathrm{mg}/\mathrm{mL}$) to simulate the corresponding scattering properties (${\mu }_{\mathrm{s}}=3.0{\mathrm{mm}}^{-1}$) and dried coffee ($10\mathrm{mg}/\mathrm{mL}$) to simulate the absorption coefficient (${\mu }_{\mathrm{a}}=0.2{\mathrm{mm}}^{-1}$).²¹ Different thicknesses of phantom epidermis layers, namely 200, 400, 500, 600, 700, and $900\mu \mathrm{m}$, were placed on top of the phantom dermis to simulate different depths of the vessels [Fig. 3(c)].

Fig. 3

Silicon phantoms for (a) dermis, (b) epidermis, and (c) the scheme of the skin phantom.

2.4.

Experimental Setup

The LSI setup (Fig. 4) consists of laser He–Ne (632.8 nm) that impinges on a static engineered diffuser (Model ED1-C20, Thorlabs Inc.), which in turn illuminates homogeneously a circle with a 1-cm diameter on the phantom’s surface. If the phantom is not homogeneously illuminated, the intensity distribution can influence the corresponding spatial contrast. A lens coupled to a CCD camera (Model Retiga 2000R, QImaging, Canada) images the surface of the phantom on the CCD. The magnification of the optical system is $M=0.285$ and the corresponding speckle size is $\approx 7.9\mu \mathrm{m}$. Alinear polarizer is placed in front of the lens with its transmission axis crossed with the He–Ne polarization to minimize the specular reflection. The CCD camera is connected to a computer to save the captured speckle images. The camera exposure time is set to $T=10\mathrm{ms}$ because this value is quite common for in vivo LSI applications.²² This exposure time is approximately three-orders of magnitude larger than the correlation time reported²³ for a similar LSI experiment at a flow speed of $8\mathrm{mm}/\mathrm{s}$. An infusion pump (Model NE-500, New Era Pump System Inc.) was used to inject the intralipid at 3% in water as blood substitute, whose scattering coefficient is ${\mu }_{\mathrm{s}}\approx 150{\mathrm{cm}}^{-1}$ at 630 nm.²⁴ The pump allows control of the speed of the intralipid through the glass capillary. In this work, we use two different flow speeds: 10 and $20\mathrm{mm}/\mathrm{s}$. A similar range of speeds has been used in similar experiments.²⁵

Fig. 4

Optical setup where the diffuser is used to produce a homogeneous illumination on the skin phantom.

3.1.

Contrast Representation

Briers and Fercher²⁶ demonstrated that it is possible to estimate the relative blood flow speed from the blurring (contrast) of a spatial and temporal integrated speckle pattern. The contrast measure ($C$) is defined as the ratio between the standard deviation ($\sigma$) and the mean intensity gray value ($\overline{\mu }$) over a sliding window $J$ in the image $I$ [Eq. (2)]. Low values of $C$ indicate a reduced range of gray levels (smoothness), whereas high values of $C$ indicate a high dispersion among the gray levels (variability).

The contrast highlights the difference between static (skin) and dynamic (blood vessel) regions of an RSI (Fig. 5), improving their visualization (contrast representation). As shown in the two views of the LSI [Figs. 5(b) and 5(c)], the contrast changes abruptly in any given region, which may hinder vessel boundaries location, and this task becomes harder as the depth increases. In a blood vessel location task, this is a drawback because it implies that vessel limits become harder to define as the depth increases. Therefore, in this work, we propose the use of a homogeneity descriptor based on a CM computing that not only takes into account the spatial relations among gray values but also provides a representation that is less sensitive to the variability of the speckle pattern. This is an important characteristic as it provides a better definition of the vessel location, as described next.

Fig. 5

(a) An RSI of a superficial vessel (i.e., $\rho =0\mu \mathrm{m}$) and (b, c) its contrast image computed with an analysis window of $5\times 5$. The homogeneity computed images using analysis windows of (d, e) $5\times 5$ and (f, g) $31\times 31$, respectively. Left and right columns show the transversal and frontal views of the images, respectively.

3.2.

Homogeneity Representation

It must be taken into account that the vessel region usually occupies a small area in the image; therefore, the homogeneity value must be computed for the small regions. Hence, the image $I^{H\times W}$ is divided into subimages or subregions through an analysis window $J^{h\times h}$ [Fig. 5(a)], and a CM is calculated for this area, just as the contrast image is calculated. With the aim of increasing the precision in the location of the vessel, a certain level of information redundancy is necessary. There must be an overlap between each pair of subregions, i.e., $J_{m+1}$ must have a displacement $\mathrm{\Delta }$ with regard to $J_m$ with $1\le \mathrm{\Delta }\le (h-1)/2$ where $h$ is the size of the analysis window. In this way, we generate a new representation of the RSI through the homogeneity image, shown in Figs. 5(d) and 5(e), in which it is possible to improve the discrimination between dynamic and static regions. As shown, high and low values of homogeneity will be associated to dynamic and static regions, respectively. Homogeneity representation enables smoother transitions among consecutive regions and defines the vessel better.

Often, the analysis window has a small size, such as $5\times 5$ or $7\times 7$, but if a large area is analyzed, it is possible to find statistics that better describe the characteristics of that area. This is particularly interesting because the speckle pattern is generated randomly. If the analysis is made using a small window, a higher dispersion is found [as shown in Figs. 5(b)–5(c) and Figs. 5(d)–5(e)]; but if a larger window is used, that dispersion is attenuated because the occurrence of gray levels becomes higher, i.e., some of them are no longer seen as isolated occurrences and the relationship among pixels is higher, although they still reflect local values. As the use of a larger window improves the homogeneity representation by avoiding strong variations, the use of a larger window with $h=31$ and $\mathrm{\Delta }=5$ will be explored. As shown in Figs. 5(f) and 5(g), the use of a larger window improves the homogeneity representation by avoiding strong variations. The transversal view of $I_{H_g}$ ($P_{H_g}$, hereinafter named as homogeneity profile) shows with more detail that low homogeneity values are related to static region, whereas high homogeneity values are related to dynamic regions, i.e., regions of high contrast. We can also observe that the homogeneity rapidly increases in the regions closest to the vessel until it reaches its maximum level at the center of the vessel. In this new representation, it is clearer to distinguish the static and the dynamic regions, thereby improving the vessel location. Although this homogeneity profile corresponds to an RSI at $\rho =0\mu \mathrm{m}$ in depth (which we call the reference case), this analysis can be extended to deeper vessels.

3.3.

Homogeneity Profile Characterization

In LSI, it is known that the vessel region is associated with a blurring effect due to the blood flow (dynamic region);¹ however, as the depth increases, this effect is reduced, making it difficult to identify the vessel region. Although a certain level of blurring remains in the image, the dynamic region becomes very similar to the static region. Therefore, we need to find other ways to identify the vessel region from an image that appears to be completely noisy. The homogeneity representation will be used for this purpose.

For the proposed methodology, an initial set of RSI (S1) was tested. The set contains a package of images obtained using the experimental setup described in Secs. 2.3 and 2.4. Table 1 describes the characteristics of this package, consisting of 30 RSI images at $\rho =0\mu \mathrm{m}$ in depth. This is done in order to study whether there is a characteristic pattern in the homogeneity behavior. The starting point (reference case) is $\rho =0\mu \mathrm{m}$. All the homogeneity profiles were found to present a distribution with a tall narrow peak and long tails due to the extreme deviation from the central tendency of the values, this feature suggest the use of kurtosis. This deviation was caused by the rapid increase of homogeneity values in the dynamic region [similar to Figs. 5(e)–5(g)].

Table 1

Characteristics of the reference set (S1) of RSI.

The kurtosis is a measure that describes the shape of a distribution based on the ratio between the central tendency and the deviation of the data.²⁷ When values are greatly deviated from the central tendency, the kurtosis takes values $>3$, and the distribution is called leptokurtic. Given that a superficial vessel at $\rho =0\mu \mathrm{m}$ has a characteristic leptokurtic homogeneity distribution, kurtosis can be used as a parameter of reference for vessel location. In order to prove this statement, the kurtosis value ($k$) was calculated for the 30 profiles, according to

Once the reference homogeneity profile has been established, it is necessary to locate the vessel and later to estimate its diameter, as shown in Sec. 4.1. A transversal cut of the homogeneity profile [Fig. 6(a)] shows that it is highly correlated with the actual dimension of the blood vessel [Fig. 6(b)]. Also, it notes that the dotted blue lines correspond to the actual vessel location, whereas the region comprised between a continuous red line and a dotted blue line is associated to the area around the vessel due to the migration of photons from the vessel to the surrounding regions [Figs. 6(a) and 6(b)]. The actual vessel region in the homogeneity representation was found to be best defined in the interval comprised between the maximum value (vessel center) and a standard deviation of the homogeneity profile, i.e., $\max (P_{H_g})\pm \sigma$ (dotted blue lines).

Fig. 6

(a) Peak analysis for vessel location in a homogeneity profile at $\rho =0\mu \mathrm{m}$ and (b) its corresponding RSI.

Nonetheless, the supposed location depends on a distribution with high kurtosis since the standard deviation may comprise a wide or narrow region, depending on the peak width of the homogeneity distribution. As observed in Fig. 7(a), the deeper the vessel, the wider the homogeneity profile (and therefore, the lower the kurtosis). Figure 7(b) shows a close up around the peak and the vertical dotted lines indicate the actual vessel location for reference purposes only. In the same figure, the circular markers correspond to $\max (P_{H_g})-\sigma$ for the different depths of the vessels. Note that the criteria for vessel location established for the reference case is no longer valid for deeper vessels. Hence, the one standard deviation is not a good parameter for locating deep vessel.

Fig. 7

(a) Homogeneity profiles of S2 at different depths. Close up around the peak vessel location (b) by taking $\max (P_{H_g})-\sigma$ and (c) by taking $\max (P_{H_g})-{\sigma }_w$.

We propose a weighted standard deviation (${\sigma }_w$) parameter [Eq. (5)] that compensates the flattening of the distribution caused by the increasing depth. The weight ($w$) depends on the ratio between the kurtosis of the homogeneity profile at a given depth $d$ ($k_d$) and the kurtosis of reference $k_0=4.42$. The deeper the vessel, the lower the ${\sigma }_w$. Thus, $w$ controls the proportion that must be considered from the standard deviation at a given depth (${\sigma }_d$); therefore, it restricts the peak width associated to the vessel region. Figure 7(c) shows the vessel location obtained with ${\sigma }_w$, as we observe it fits best with the actual vessel location,

Fig. 8

General diagram of the proposed methodology.

4.1.

Linear Vessel

In order to assess the proposed methodology under different conditions, we tested it using another set of RSI obtained at different depths as described in Table 2.

Table 2

Characteristics of the set S2 of RSI used for testing the proposed methodology.

For the set S2, the vessel was manually located among the rows [141, 168] for $\rho =0\mu \mathrm{m}$, and this location is used as a reference for future analysis, whereas for other depths, the location of the straight vessels was done automatically by searching for the peak of the homogeneity profiles and performing the kurtosis analysis. The proposed methodology was first tested in straight vessels by processing homogeneity profiles for each column and computing the kurtosis analysis. Figure 9 shows the vessel extracted from S2 at different depths and with a velocity of $10\mathrm{mm}/\mathrm{s}$. As one can see, the located vessel is better defined at lower depths. At higher depths, the image is noisier, and therefore, the homogeneity values in static and dynamic regions become similar and harder to define.

Fig. 9

Automatic vessel location at different depths by the proposed methodology: (a) $\rho =0\mu \mathrm{m}$, (b) $\rho =400\mu \mathrm{m}$, (c) $\rho =700\mu \mathrm{m}$, and (d) $\rho =900\mu \mathrm{m}$; colors are related to the homogeneity values.

The goal of performing in vitro tests at different depths is to demonstrate that ${\sigma }_w$ allows the estimated vessel location to match with the real vessel location, and therefore, it is possible to measure the vessel diameter. For our skin phantoms, it is a known value ($700\mu \mathrm{m}$). Thus, based on the obtained location at each depth, we were able to estimate the vessel diameter and its deviation from its actual value. Figure 10 shows the plot of this estimation for set S2 by comparing the results at two velocities. We can observe that the speed does not seem to affect the estimation significantly of the vessel diameter as the majority of the values lies within 10% of the actual vessel diameter (dashed lines). The plot shows how close the automatic location is to the reference diameter (solid line).

Fig. 10

Vessel diameter estimation for S2 at two different velocities.

4.2.

Bifurcated Vessel

In in vivo situations, the vessels are not straight but takes complicated forms and even bifurcations. In order to apply our proposed kurtosis analysis, a bifurcated vessel with a diameter of $300\mu \mathrm{m}$ and a depth of $\rho =400\mu \mathrm{m}$ was tested [Fig. 11(b)]. In this vessel morphology, the analysis of distributions is performed along each column and each row.

Fig. 11

(a) Homogeneity profile for $c=10$, (b) bifurcated vessel at $\rho =400\mu \mathrm{m}$ in the RSI, (c) homogeneity profile for $r=165$, and the segmentation of the bifurcated vessel according to the estimated location using kurtosis analysis over (d) a contrast and (e) a homogeneity representation.

Figure 11(a) shows the unidimensional profile corresponding to column $c=10$, where the two main peaks indicate a high homogeneity area related to the left branches of the vessel [Fig. 11(b)]. The same analysis is performed but along the row $r=165$ [Fig. 11(c)]. It is important to remark that only the most prominent peaks are considered, and each one is analyzed as an individual distribution. Evidently, these peaks match with the center of the vessel branches but they cover a wider area than the reference as expected. For instance, the upper peak in Fig. 11(a) comprises rows [210, 270]. However, using the proposed parameter ${\sigma }_w$, a better location of the vessel is achieved once again now between rows [233, 251]. The final location of the vessel is given by adding both results, and it is used as a marker. When this marker overlaps with the homogeneity representation, it generates the segmented image of the vessel. Figures 11(d) and 11(e) show the result of automatic vessel location performed in the contrast and homogeneity representations. Although the vessel is well located in both representations, it presents a poorer definition in the contrast representation unlike the homogeneity representation, which retains better the vessel shape. The white lines in Fig. 11(b) indicate the ideal definition and location of the vessel.

4.3.

Discussion

The full width half maximum (FWHM) is a common measure often used to estimate the width of a function or signal.²⁹ In this way, the FWHM and ${\sigma }_w$ accomplish a similar function. Here, we compare the performance for vessel location of the FWHM and $w$ measurements for each depth (set S2). Figure 12 shows the differences in the vessel diameter when the vessel location is estimated with ${\sigma }_w$ and with FWHM. We can observe that the FWHM has a marked deviation from the actual diameter, mainly for higher depths since the homogeneity profiles tend toward a normal distribution. As shown in Fig. 10, the deviation of the automatic location proposed here is smaller than the 10% (dashed lines) of the diameter of the reference (black continuous line).

Fig. 12

Comparison of vessel diameter estimation using ${\sigma }_w$ and FHWM in set S2.

Table 3 shows a comparison among different methods based on LSI analysis to improve the visualization of deep (hundreds of microns) blood vessel. Some of them use an external agent, for example, photothermal-LSI uses a 595-nm laser pulse to thermally excite the blood flow, resulting in a local increase of temperature and thus, a transient decrease in speckle contrast;⁷ this technique should not be applied on temperature-sensitive areas, such as brain. PCTOC uses glycerol to reduce the scattering of in vivo and ex vivo skin samples; 60 min after the application of PCTOC, a sequence of RSI is acquired and processed with the LSI algorithm.⁹ LSI-SFDI uses structured illumination and LSI to separate a 2-mm-diameter shallow vessel (2 mm depth) from a deeper one (4 mm depth) but the vessel diameter is not estimated.²² Postnov et al.¹⁰ proposed an image-processing algorithm based on second-order gradient analysis of LSI images to estimate the vessel diameter and blood flow from an exposed rat mesenterium and a mouse cortex. Our results compare quite well with the previous results but, in addition, can improve the visualization of deep vessels and estimate its diameter without the use of external agents and in a noninvasive way.

Table 3

Reported works for visualizing deep blood vessels.

Although in all image-processing methods there is a loss of information, due to the required processes for generating the homogeneity representation and the kurtosis analysis, there is also a trade-off in the location of the blood vessel that has demonstrated to be accurate. Moreover, the homogeneity representation visually improves the blood vessel definition, and whereas the loss of information is inevitable, new information is also gained (vessel location). It is important to remark that several versions of the information are generated as the information is being processed and then, at the end of our process, the segmented image can be used as a marker. The resulting marker in our proposal can be used just for indicating the vessel location, as the white lines in Fig. 11(b), over the homogeneity representation or over the raw speckle. In this way, original and new information can be complemented. Even more, if a contrast representation is also computed, the marker can be used over it.