Imaging of microcirculatory in vivo permits the direct visualization of long-term microhemodynamic, which is closely associate with disease process^[1,2], neural dynamics^[3], and functional recovery from pathological states^[4]. Photoacoustic microangiography^[5,6] is an innovation technique used to obtain blood vessels from biological tissues for microcirculatory disease detection^[7]. It employs a modulated light pulse to irradiate biologic tissues. Due to the contrast in the absorption coefficient between blood vessels and the background medium, a high quality image of tissue absorption properties is reconstructed by the received ultrasound signal^[8].

Photoacoustic microangiography provides the direct visualization of blood vessels^[9,10]. Usually these images are interpreted qualitatively. Therefore, a quantitative method for characterizing blood vessels would have several clinical applications. Previous methods include measuring blood vessel diameters^[11], the flow velocity^[12], and the maximum distance to the nearest blood vessel^[13]. It is worth noting that all of these parameters are intensity information. A parameter that can describe the vessel tortuosity would be more beneficial. For example, a change in retinal vessels is an early indicator of coronary heart disease^[14] and stroke^[15]. Vascular remodeling in which vessel tortuosity plays an important role has also been of interest in several fields^[16]. Thus, four measurement parameters are considered to give a more complete description.

In order to quantify the morphology of blood vessels, a segmentation algorithm is needed first. The segmentation algorithm returns a binary map of locations of vessels. Numerous segmentation methods have been proposed and the structure or intensity information is explored. The simplest one is the adaptive threshold algorithm that is based on the intensity. However, the intensity-based techniques are very sensitive to the threshold parameter selection while lack sensitivity to the morphology of the blood vessels, which often leads to over segmentation.

By exploring both the shape and direction of vessels, the Hessian filter is a good fit for multisize blood vessels^[17–19]. Typically, a known range of scales is required to maximize the vesselness function. However, if the maximum scale is chosen large enough for the large vessels, blurriness and enlargement of the smaller vessels will be clearly observed.

In this Letter, a multiscale Hessian filter-based algorithm is proposed. The limitation of the Hessian filter has been analyzed and corrected by a local adaptive threshold algorithm. Additionally, the segmented binary image is utilized to acquire four measurement parameters to further quantify the blood vessels. Finally, the segmentation and quantification method is valid on the whole and a small area of photoacoustic microangiography to identify different tissue characteristics.

The first step of the proposed method is the segmentation progress. The algorithm we considered is based on the multiscale Hessian filter. The multiscale Hessian filter proposed by Frangi^[20] has been widely used in many applications^[21–23]. For a Hessian filter, the local behavior of the second-order gradient image, which is called Hessian matrix, is utilized to identify the boundaries of the blood vessels from the background and other nonvessel tissues. The Hessian matrix H(I)s can be computed as H(I)s=(∂∂x)(∂∂x)I(x)=s2γI(x)*∂2∂x2G(x,s),(1)where I(x) is the original image intensity at location x, γ is the regularization parameter, and G(x,s) is the Gaussian kernel at scale s, that is G(x,s)=12πs2e−‖x‖22s2.(2)By using the convolution of the image with the Gaussian kernels, the multiscale nature that allows for the identification of blood vessels of different sizes has been enforced. In other words, with a reasonable range of scales, all the blood vessels in different sizes could be efficiently segmented.

Figure 1 shows the principle of the Hessian filter. By analyzing the eigenvalues and the eigenvectors of the Hessian matrix, the vessel information can be extracted. For a 3D structure, the three eigenvalues could be acquired and placed in the order |λ1|≥|λ2|≥|λ3|. For an ideal vessel structure, the eigenvalues should satisfy the following condition: |λ3|≈0,|λ3|≪|λ2|,λ2≈λ1.(3)In order to identify the blood vessel from the other nonvessel tissues by the above condition, two geometric ratios are defined as RA=|λ2||λ3|,(4)RB=|λ1|λ2λ3.(5)The first ratio, RA, can distinguish the blob-like structure from other tissues. Meanwhile, the second one, RB, can distinguish between plate-like and line-like structures.

Moreover, in order to identify the blood vessels from the background, a third ratio is introduced that is defined as RC=∑j≤Dλj2,(6)where D is the dimension of the image. The ratio RC will be very low in the background because there is no tissue except noise.

Therefore, the vessel function at scale s is defined as vs(x)={0ifλ2<0orλ3<0(1−e−RA22α2)*e−RB22β2(1−e−RC22θ2)others,(7)where α, β, θ control the sensitivity to the three ratios, respectively. The vessel function shows the probability of the location belongs to a vessel at a certain scale. The function is computed within multiscales and the value of the vessel function is maximized at a scale that approximately matches the vessel size, V(x)=argmaxs{vs(x)},s∈[smin,smax],(8)where smin and smax are the minimum and maximum scale, respectively.

By using the multiscale Hessian filter, different size of blood vessels could be extracted. However, the method is very sensitive to the maximum scale. Figures 2(a)–2(c) show the segmentation results obtained by the multiscale Hessian filter using different maximum scales 1, 7, and 10. The blur and enlargement effect can be clearly observed, especially in Fig. 2(c) for an overly large maximum scale is used. Figures 2(d)–2(f) are the corresponding vessel diameter quantification maps. The results demonstrate that the quantified diameter is very sensitive to the segmented results.

In order to minimize the sensitivity to the maximum scale parameter, a local adaptive threshold method is incorporated. It is utilized in parallel with the multiscale Hessian filter. The final result is acquired by compounding these two results with a weighted average scheme, Iout=α×IH+(1−α)×IT,(9)where Iout is the final segmentation result, IH and IT are the results obtained by the Hessian filter and local adaptive threshold method, respectively. The weight scale α is a scalar value between [0, 1]. As long as the compounding result is not over segmentation, α should be as small as possible.

The performance of the multiscale Hessian-based segmentation algorithm is shown in Fig. 3. Figure 3(a) shows an original photoacoustic microangiography of microcirculation tissue beds in vivo. It is generated by an optical-resolution photoacoustic microscopy system, which has been reported in Ref. [24]. Figure 3(b) shows the segmentation results using the local adaptive threshold method. Here, a window size was defined to determine the local adaptive threshold. The threshold was determined by the mean value within the window and window size is 3 pixels. Because this method just uses the intensity information, it will produce false segmentation. All pixels in the image are treated equally, thus the boundary of the vessels cannot be clearly identified. Moreover, over segmentation will be caused by an inappropriate threshold. Figures 3(c)–3(e) show the segmented results by the Hessian filter and the proposed method. The maximum scale used for Figs. 3(c) and 3(e) is 8, to make sure all of the vessels in the image could be extracted, while the value for Fig. 3(d) is 10, which is larger than the most appropriate one. By visual comparison, it can be observed that the results of our method are much better than the multiscale Hessian filter. The blurring artifact and unwanted enlargement due to the inaccuracy choice of the maximum scale could be appropriately solved even when the value was selected too large.

Based on the better-segmented binary map, four measurement parameters (fractal dimension, vessel length fraction, vessel density, and vessel diameter) will be quantified.

The vessel diameter is the most commonly used parameter. The distance transform of a blood vessel is the minimum number of pixels between each foreground pixel to the boundary of the vessel^[25]. The results should have a maximum value on the vessel centerline. The exact vessel diameter can be measured after correcting by the spatial size of each pixel. Here, the distance transform result is used to represent the quantified vessel diameter. Both Figs. 2(d)–2(f) and Fig. 3(e) are quantified vessel diameter maps.

Vessel density and vessel length fraction are parameters that represent a relative value of the total area occupied by the vessels and the total length of the vessels, respectively^[26]. The vessel density can be acquired directly on the segmented image. It is calculated as the number of the white pixels, which represent the total area covered by the blood vessels, divided by the total number of images, which represents the total area of the imaging area. The vessel density value quantified for Fig. 3(a) is 0.59896.

For the vessel length fraction, skeletal images that represent the total length of the blood vessels are needed. The skeletonization process consists of iteratively deleting the pixels in the outer boundary of the segments until a single pixel width line is obtained^[27]. In other words, the skeletal image can represent the centerline of the vessel. Figure 4 shows a skeletal image of the Fig. 3(a). Therefore, the vessel length fraction can be calculated as the number of white pixels divided by the total number of skeletal images. The vessel length fraction quantified for Fig. 3(a) is 0.13663.

Fractal dimension is a parameter to characterize a self-similar image^[28]. A fractal dimension is a value that gives an indication of how an image fills space into smaller scales. Here, the parameter is utilized to quantify the vessel turtuosity. It has been applied in diverse areas of medicine to describe complex biological structures such as branching patterns of the retina, coronary, and pulmonary arterioles^[29]. It has also been used to quantify the fractal distribution of scatters in tissues, the parafoveal capillary network, and optical coherence tomography images of arteries^[30].

Although the fractal dimension can be computed both on the segmented image and the skeleton map, the result acquired from the skeleton map is more sensitive to changes of the vessels^[31]. Thus, in this Letter, the fractal dimension are all calculated on the skeletal image. A box counting method is utilized to calculate the fractal dimension. It is a method of estimating the fractal dimension from structures that are not perfectly self-similar^[32]. More importantly, the box counting method can also be used in the quantification of the fractal dimension on a small area. The fractal dimension value quantified for Fig. 3(a) is 1.8734.

In the study of several microvascular phenomena, such as angiogenesis (growth of new blood vessels), it is important to quantify small areas of tissue. It can help us find the location of the diseased tissues in the region of interest (ROI). For example, regions close to tumors may present angiogenic blood vessels (higher tortuosity and fractal dimension) compared to the healthy surrounding blood vessels.

For the small area quantification, the image is cropped to create smaller ones. Therefore, the vessel length fraction, vessel density, and fractal dimension would be calculated over the smaller image. Additionally, there is no extra change to quantify the vessel diameter on small areas.

For Fig. 3(a), in order to quantify the small area, a window of a given size is used and the calculated measurement parameter values are stored in the center pixels of the window. By sliding the window across the whole image, a color map can be obtained. The size of the window was appropriately reduced at the borders of the image. In practice, the window size should be large enough to include one or more vessels. The window size we used here is 8×8 pixels.

The vessel length fraction, vessel density, and fractal dimension from the photoacoustic microangiography [Fig. 3(a)] are presented in Figs. 5(a)–5(c), respectively. For ease of visualization, the images were multiplied by the segmented images. The edges of the image present an artifact due to the window having smaller size. The vessel diameter map is the same with the Fig. 3(e).

Two ROIs have been selected as shown in Fig. 6. Region 1 is covering large blood vessels, while Region 2 is covering small vessels and capillaries. The mean and standard deviations of the parameters within the two ROIs are presented in Fig. 7.

In general, the large blood vessels regions has larger diameter values than the smaller ones. Figure 7(a) shows the relative diameter values of the two regions. The vessel length fraction and fractal dimension show smaller values in the regions within large blood vessels compared to the capillaries in Figs. 7(b) and 7(d). This is expected since large vessels are fairly smooth compared to the capillaries, which are more tortuous. Also, a few large vessels have a smaller length compared to several small vessels covering the same area. On the other hand, the vessel density is slightly higher when it covers large blood vessels because the area of the vessels covering the image is larger (Fig. 7(c)). This agrees with the real situation.

In conclusion, a multiscale Hessian-based segmentation and quantification method is proposed for photoacoustic microangiography. In the proposed segmentation algorithm, the blurring and enlargement that are limitations of the Hessian filter are corrected. The results of the algorithm can be utilized to get more effective measurement parameters. The vessel diameter, vessel density, vessel length density, and fractal dimension are quantified to give both intensity and tortuous information of the blood vessels. Moreover, the segmentation and quantification method is applied on a small area within the photoacoustic microangiography to give a quantified color map. This is very important to properly characterize different ROIs within an image. In the future, the proposed method for a small area could be used to monitor the morphological changes in regions close or far away from a diseased region, such as a burn or cancer.