1.

Introduction

Photoacoustic microscopy (PAM)^{1, 2, 3} is a newly developed photoacoustic (PA) imaging^{4, 5, 6} technology that has significant potential for various biomedical research and clinical applications. It is safe to animals and humans and has been used to image the subcutaneous microvasculature in rats and humans,^{2, 7} acute skin burns in swine,⁸ and skin melanoma tumors in mice⁹ with a $15\text{-}\mu \mathrm{m}$ axial resolution, a $45\text{-}\mu \mathrm{m}$ lateral resolution, and a $3\text{-}\mathrm{mm}$ maximum imaging depth. Functional properties such as total hemoglobin concentration and hemoglobin oxygen saturation $\left(\mathrm{s}{\mathrm{O}}_2\right)$ also have been imaged in vivo in single blood vessels in rats.^{2, 10} Based on its high optical absorption contrast and large maximum imaging depth, PAM has great advantages over existing high-resolution optical imaging modalities for imaging both the anatomy and functions of subcutaneous microvasculature.⁷

The system description and the image formation procedure of PAM have been reported in detail elsewhere.^{2, 3} Here, we only reiterate the definitions of the key elements in image formation and the calculation of the maximum-amplitude-projection (MAP) image.^{11, 12, 13} As shown in Fig. 1, the A-line is a one-dimensional (1-D), depth-resolved image (along the $z$ axis) converted from the time-resolved PA signals recorded at each horizontal (along the $x\text{-}y$ plane) position. The recording duration of the PA signal is determined by the desired depth range along the $z$ axis and the sound velocity in soft tissues $(1.5\mathrm{mm}∕\mu \mathrm{s})$ . A 1-D lateral scanning of the ultrasonic detector along the $x$ axis produces a two-dimensional (2-D) cross-sectional image (referred to as a B-scan image), and a complete raster-scanning along the $x\text{-}y$ plane produces the final volumetric image.

Fig. 1

Image formation in PAM.

Besides direct three-dimensional (3-D) visualization of the volumetric image,⁷ a 2-D MAP visualization, which has been widely used in radiological angiography, is usually adequate if the depth information of vessels is less interesting. MAP exploits the fact that within the volumetric data of vasculature, the PA amplitudes generated from the vessels are much higher than those from the surrounding tissues; thus, the structure of the vasculature can be well captured.^{2, 9, 10} MAP images can be calculated by performing ray casting and searching for the maximum amplitude on the ray line along an arbitrary direction. In in vivo PAM imaging on small animals, the major motion artifact is caused by animal breathing, which exists primarily along the $z$ axis, according to the design of the animal holding.³ Thus, we only calculated MAP images along the $z$ axis to avoid motion artifacts (Fig. 1).

By employing a dark-field illumination with a large illumination area,³ an overly strong PA signal from the skin surface is avoided, which potentially overshadows the weaker PA signals generated from the subcutaneous vessels.¹ However, due to light scattering in the tissue, PA signals from the skin surface are still detected, as shown in Fig. 2a . Moreover, because of the uneven distribution of pigments (such as melanin) in the skin, the PA amplitudes generated from the skin surface vary with horizontal positions. In some A-lines, the PA amplitude from the skin signal can be much higher than that from a vessel signal [Fig. 2c]. If the MAP image is calculated directly, PA amplitudes from the skin surface rather than those from the vessels are projected onto the final 2D image, resulting in errors in the visualization of the true microvasculature. Therefore, detecting the skin profile and subsequently removing the skin signals are necessary before MAP images are calculated.

Fig. 2

Relationship between the $A_S$ and $A_V$ in a typical in vivo B-scan image: (a) B-scan image, (b) an A-line where $A_S$ is comparable to $A_V$ , (c) an A-line where $A_S$ is stronger than $A_V$ , and (d) an A-line where $A_S$ is weaker than $A_V$ (S: skin surface; V: vessel; 2nd: the second strongest peak).

The skin profile was detected manually in the early works of PAM.¹ First, several points on the skin surface were visually selected in multiple B-scan images, then a 2-D interpolation through the selected points was applied to estimate the skin profile. This manual operation is time-consuming, however, since the number of B-scan images often exceeds 200. If spectral measurement is conducted, the same manual operation needs to be repeated for all the optical wavelengths. As a result, an automatic algorithm for detection is necessary.

Automatic skin profile detection is also important for imaging samples with large variations in skin contours. To achieve high lateral resolution, PAM employs an ultrasonic detector with a high center frequency $\left(50\mathrm{MHz}\right)$ and a large numerical aperture (0.44),^{2, 3} which results in a focal zone of only $0.3\mathrm{mm}$ . Lateral resolution is degraded at the out-of-focus region. Since the raster scanning is usually conducted along the $x\text{-}y$ plane at a fixed vertical coordinate (flat scan) when the fluctuation of the skin contour is much larger than the focal zone, vessels out of the focal region cannot be resolved correctly. Thus, further studies, such as the validation of Murray’s law,^{14, 15} that rely on the imaged vessel structure can be significantly affected. A virtual-point-detector-based synthetic aperture focusing technique was developed to restore the degraded out-of-focus lateral resolution numerically;¹⁶ however, only 1-D restoration is applicable to in vivo experiments since motion artifacts between B-scan images are too strong. Moreover, the change of acoustic velocity in water and tissue is not considered. Based on the anatomic observation of the subcutaneous microvasculature, the major vessels are within a layer of limited thickness beneath the skin surface; therefore, good ultrasonic focusing of such a vessel layer can be maintained if the ultrasonic detector can follow the skin contour during the raster scanning (auto-fit scan). Thus, detecting the skin profile is a prerequisite for the auto-fit scan.

We have developed an automatic algorithm to detect the skin profile by analyzing A-lines and nonparametric smoothing in each B-scan image. PA signals generated from the skin surface were then removed from the volumetric data to calculate the MAP images. Also, by applying auto-fit scan, vessels that were previously blurred in the flat scan became resolved. In this paper, we describe such an automatic algorithm and the procedure of the auto-fit scan. Comparison of in vivo MAP images before and after the skin signal removal and comparison of in vivo MAP images from the flat scan and the auto-fit scan, respectively, are provided. The application of the auto-fit scan is further demonstrated by in vivo imaging of a tumor vasculature.

2.

Materials and Methods

2.2.

Skin Profile Detection

In PAM, modern edge-detection algorithms^{18, 19} cannot be applied directly because the imaged skin surface shows no continuity [Fig. 2a]. Here, skin profile detection is facilitated by analyzing the features of A-lines and is carried out on a B-scan basis (along the $x$ axis). The detection of a 1-D skin profile in a B-scan image consists of two steps. The first step is a rough estimation of the skin location based on analyzing the relationship between the amplitudes of PA signals generated from the skin surface and from subcutaneous blood vessels in each A-line. The second step is a nonparametric robust local regression smoothing^{20, 21} of the rough-estimation result. After the skin profiles are detected in all the B-scan images, the robust local regression smoothing is performed again along the $y$ axis. Finally, the complete 2-D skin profile is acquired by applying a 2-D Gaussian low-pass spatial filtering.

Within the depth range, most A-lines have only two dominating PA peaks, which are generated from the skin surface and a subcutaneous vessel, respectively. The relationship between the amplitude of PA signal generated from the skin surface (denoted as $A_S$ ) and the amplitude of PA signal generated from a subcutaneous blood vessel (denoted as $A_V$ ) can be classified into three major categorizes: (1) $A_S$ and $A_V$ are comparably strong [Fig. 2b]; (2) $A_S$ dominates [Fig. 2c]; and (3) $A_V$ dominates [Fig. 2d]. Figure 2b demonstrates the most ideal case as we can easily select the location of $A_S$ to be the position of the skin surface. However, such an easy method is not valid for other situations. In addition to the major relationships between $A_S$ and $A_V$ , other rare situations also exist. In some cases, more than two dominating PA peaks exist within an A-line when, for example, multiple PA peaks can be generated from multiple blood vessels at different depths at the same horizontal position. In other cases, an A-line could show only one or no dominating PA peak, i.e., it only has a PA peak generated from the skin surface or it consists only of random noise and background signals without distinctive PA peaks from either the skin surface or blood vessels. All these situations are generalized to have $N$ peaks in the first-step rough estimation for the 1-D skin profile in each B-scan image. $N$ is a varying number that can be adjusted according to each special situation.

In B-scan images, A-lines are analyzed in an order that follows the data acquisition sequence. For each A-line, the depth locations of the first $N$ strongest PA peaks $(L{P}_i,i∊[1,N])$ are identified and the distance between each pair of consecutive peaks ( $D_j=L{P}_{j+1}-L{P}_j$ , $j∊[1,N-1]$ ) are further calculated. Once the largest pair distance $(D_{j\max }=\max \left\{D_j\right\})$ is found, $L{P}_{j\max }$ is considered the location of the skin surface in that A-line. Such an operation is based on the observation that the distance between the PA peak generated from the skin surface and the PA peak generated from a shallower blood vessel is generally larger than the distances between any other consecutive pairs within the same A-line.

After all A-lines within a B-scan image as shown in Fig. 2a are processed as described above, a 1-D rough estimation of the skin profile is acquired as shown in Fig. 3a, where the detected skin locations at the positions of the three A-lines shown in Fig. 2 are labeled. In this special case, $N$ was 2. In Fig. 2b, where $A_S$ and $A_V$ are comparable, the rough estimation gives the correct skin location. In Fig. 2c, where $A_S$ dominates, the second strongest PA peak selected by the algorithm is near the real skin surface. Therefore, the estimated skin location is close to the true location; however, in Fig. 2d, where $A_V$ dominates, the second strongest PA peak is behind the PA peak from the vessel. Hence, the rough estimation gave an error. Several similar erroneous estimations also can be observed in Fig. 3a.

Fig. 3

Detected 1-D skin profiles: (a) comparison of the rough and smoothed estimations of the skin profile and (b) overlay of the smoothed skin profile onto its original B-scan image.

After the rough estimation, a better estimation of the 1-D skin profile was achieved by removing the error points (outliers) through smoothing. The assumption behind the smoothing operation was that the majority of the estimated skin locations in the first-round rough estimation were either true or were close to their true skin locations, and the outliers were randomly distributed along the $x$ axis. Hence, a method is needed that can extract the backbone based on the estimated major skin locations and that will not be affected by the outliers. Moreover, such a method should not require human intervention. Among the several smoothing technologies that have been reported,²² the nonparametric, robust, locally-weighted-regression smoothing^{20, 21} was found to fit our applications well based on the trade-off between computation time and capability to reject outliers from the first-round rough estimation. Compared with parametric smoothing, a nonparametric smoothing showed a greater advantage on the degree of freedom since it does not require a single expression of the smoothing result. Hence, nonparametric smoothing offers better approximations than parametric smoothing when applied to the irregularly shaped skin profile. Moreover, robust regression “guards against deviant points distorting the smoothed points”²⁰ and the local weighting allows for the variation in the skin profile. The robust regression requires longer computational time than the direct parametric and nonparametric regression smoothing. Since this skin profile detection was conducted offline and the time required for regression (around two minutes for a complete volumetric data containing $120\times 100$ A-lines) was much shorter than the data acquisition time $\left(40\min \right)$ , computational cost was not a concern here. The smoothing procedure was carried out in a piecewise manner for each point span. The length of a point span is defined by the ratio of the number of points within the span to the total number of A-lines in a B-scan image.

The smoothed skin profile, whose span was 0.4, is presented in Fig. 3a. It can be seen that the outliners are rejected and the underlying backbone is extracted. Figure 3b shows the overlay of the smoothed result onto the original B-scan image, where a good agreement is demonstrated. After the 1-D smoothed skin detection is completed for each B-scan image, the same smoothing procedure was performed again along the $y$ axis. The final 2-D skin profile was further acquired after a Gaussian low-pass spatial filtering.

3.

Result and Discussion

Both the skin profile detection algorithm and the auto-fit scan were tested on a known phantom. A brass sphere with a diameter of $8\mathrm{mm}$ was imaged as shown in Fig. 4 . Due to the nature of this phantom, PA signals were only generated from the surface of the brass sphere. In Fig. 4a, the detected sphere surface (red dashed line; color online only) is overlaid onto the B-scan image. It can be clearly observed that the detected profile overlaps the imaged sphere surface, which indicates that our algorithm found the surface profile correctly. In this special situation, $N$ was 2. Once the skin profile was detected, an auto-fit scan was conducted. The middle location in the B-scan image from the flat scan was taken as the reference position. The B-scan image from the auto-fit scan and its correspondingly detected surface profile are shown in Fig. 4b. Since the distance between the sphere surface and the ultrasonic detector was maintained, both the newly imaged surface and detected surface profile appear to be flat.

Fig. 4

Detecting the surfaces of a brass ball using the automatic algorithm in both (a) flat scan and (b) auto-fit scan.

Once the skin profile was detected in a volumetric image, the skin signals were removed by setting PA values to zero within a range that starts from the beginning of the A-line to a certain depth beneath the corresponding detected skin position. Such a depth represented the skin thickness, which was an empirical value varying with animals.

The comparison between in vivo MAP images before and after skin signal removal is shown in Fig. 5 . The span here was 0.6 and the skin thickness was $0.2\mathrm{mm}$ . The detailed vessel structure that is overshadowed by stronger signals from the skin surface, as shown in Fig. 5a, is clearly seen in Fig. 5b. However, since the image was acquired from a flat scan, degraded lateral resolution (resulting in vessel blurring) due to the skin contour variation can be observed at the top and bottom regions in Fig. 5b, as indicated by the arrows.

Fig. 5

In vivo MAP images of subcutaneous microvasculature in a Sprague-Dawley rat (a) before and (b) after the skin signal removal.

The auto-fit scan overcomes the blurring problem because the ultrasonic focus is maintained on the subcutaneous vessel layer within the entire field of view. The MAP images acquired from a flat scan and from an auto-fit scan, respectively, are compared in Fig. 6 . In Fig. 6a, the 1-D detected skin profile from the flat scan scan was overlaid onto the B-scan image, where large skin contour variations were observed. Figure 6b shows the B-scan image acquired by auto-fit scan from the same cross-section as in Fig. 6a. Although the skin contour in the volumetric image acquired from the auto-fit scan is flat, it can be recovered based on the detected skin profile from the flat scan. In Fig. 6c, good ultrasonic focusing was achieved within the central region, whereas the top left region (within the dashed box) was out of the ultrasonic focus due to the changing skin contour. As a result, no vessel structure could be observed within the dashed box due to degraded lateral resolution. In comparison, good ultrasonic focusing was achieved within the area that was previously out of focus in the auto-scan, and a clear microvascular structure was evident [Fig. 6d]. Such an image quality improvement warrants the accuracy of further study based on the anatomical structure of the imaged vessel network.

Fig. 6

Comparison of images acquired from flat scan and auto-fit scan: (a) B-scan image from the flat scan overlaid with the smoothed estimation of its skin profile at the location marked by the dashed line in (c), (b) B-scan image from the auto-fit scan overlaid with the smoothed estimation of its skin profile at the same location as marked in (c), (c) MAP image from the flat scan, and (d) MAP image from the auto-fit scan.

A direct application of the auto-fit scan can be found in the imaging of subcutaneous tumor angiogenesis where large variations in skin contour exist. The vasculature of a subcutaneously inoculated tumor is shown in Fig. 7 . The elevation of the skin bump around the tumor region was $1.9\mathrm{mm}$ , which was much larger than the ultrasonic focal zone $\left(0.3\mathrm{mm}\right)$ . Using the auto-fit scan, the ultrasonic transducer successfully followed the skin contour and, as a result, the vessel structure within both the tumor and the neighboring regions were both clearly resolved regardless of the large variations in vertical elevations. Such a capability to resolve tumor vasculature in vivo provides the foundation to extract information such as microvessel density²³ and vessel branching²⁴ for further studies of cancer biology.

Fig. 7

MAP image of vascular network in a subcutaneous tumor acquired from auto-fit scan. The clouds labeled as 1 and 2 in the image were found to be subcutaneous hemorrhages.

Since PAM is a new technology with unique features in its B-scan image as discussed above, established image processing technologies cannot be directly applied to PAM images. Although many edge detection algorithms have been developed, they rely on the change of grayscale values (or their higher orders of derivatives) and generally require a closed object boundary. Moreover, most of the algorithms need prior knowledge of the image content and human inputs either to start or stop the computation. In our algorithm, once the only two parameters, $N$ and the span, are specified, no knowledge of the skin and vasculature and no human interventions are necessary. Moreover, it works well with open object boundary in cross-sectional images from PAM and, presumably, images from optical coherence tomography.

We plan to further improve the skin profile detection by either achieving real-time skin profile estimation or reducing the time required for the first-round flat scan. First, because detection of the skin profile used in our experiments was conducted offline, a complete flat scan was required before the auto-fit scan. Hence, the second-round auto-fit scan might face a different skin contour due to animal motion in in vivo experiments. To solve this, an ultrasonic pulse echo measurement can be employed at each horizontal location to detect the location of the skin surface and adjust the vertical position of the ultrasonic transducer in real time. The ultrasonic measurement can be performed between two laser pulses so that the total data acquisition time is not prolonged. However, this method requires modifying the imaging system to switch between ultrasonic and photoacoustic detections although these two detections share the same ultrasonic detector.

Second, the data acquisition time in the first-round flat scan can be reduced by minimizing the number of scanning steps and enlarging the scanning step size. The same skin profile detection algorithm can then be applied to acquire a sparse skin profile that can be further interpolated to match the desired number of scanning steps. Moreover, the computation time for the application of our algorithm will also be reduced as there are fewer A-lines to be processed and smoothed.

4.

Summary

In summary, we have developed an algorithm to automatically detect the skin profile in volumetric PAM data. Such an algorithm does not require human involvement and has been demonstrated to be essential in calculating MAP images. Furthermore, an auto-fit scan is developed for PAM to achieve good ultrasonic focusing of the subcutaneous vasculature when the variation of the skin contour is beyond the ultrasonic focal zone. This algorithm may also find applications in other imaging techniques such as ultrasonic imaging and optical coherence tomography, where depth-resolved A-lines are acquired.