Photoacoustic microscopy (PAM)1, 2, 3 is a newly developed photoacoustic (PA) imaging4, 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,8 and skin melanoma tumors in mice9 with a axial resolution, a lateral resolution, and a maximum imaging depth. Functional properties such as total hemoglobin concentration and hemoglobin oxygen saturation 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.7
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 axis) converted from the time-resolved PA signals recorded at each horizontal (along the plane) position. The recording duration of the PA signal is determined by the desired depth range along the axis and the sound velocity in soft tissues . A 1-D lateral scanning of the ultrasonic detector along the axis produces a two-dimensional (2-D) cross-sectional image (referred to as a B-scan image), and a complete raster-scanning along the plane produces the final volumetric image.
Besides direct three-dimensional (3-D) visualization of the volumetric image,7 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 axis, according to the design of the animal holding.3 Thus, we only calculated MAP images along the axis to avoid motion artifacts (Fig. 1).
By employing a dark-field illumination with a large illumination area,3 an overly strong PA signal from the skin surface is avoided, which potentially overshadows the weaker PA signals generated from the subcutaneous vessels.1 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.
The skin profile was detected manually in the early works of PAM.1 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 and a large numerical aperture (0.44),2, 3 which results in a focal zone of only . Lateral resolution is degraded at the out-of-focus region. Since the raster scanning is usually conducted along the 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;16 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.
Materials and Methods
In vivo images of the normal subcutaneous microvasculature were acquired from healthy Sprague-Dawley rats ( , Charles River Breeding Laboratories, Wilmington, Mass.) and immunocomprised nude mice ( , Harlan Co., Indianapolis, Ind.).An in vivo image of the tumor angiogenesis was acquired from a BALB/c mouse ( , Harlan Co., Indianapolis, Ind.). A subcutaneous tumor was produced by inoculation of murine colon carcinoma tumor cells (CT26.WT, American Type Culture Collection, Manassas, Va.) with injection volume under the scalp. Imaging was conducted after inoculation when a bump at the inoculation site was visually observable.
Before imaging, the region of interest was depilated with a commercial human hair removing lotion (Surgi Cream, Ardell International, Los Angeles, CA). A dose of of Ketamine plus of Xylasine was administered intramuscularly to anesthetize the animals. During experiments, the animal motion was minimized by a breathing anesthesia system (E-Z Anesthesia, Euthanex Corporation, Palmer, Pa.). The flow rate of the mixture of medical grade oxygen with 1% (volume to volume ratio) vaporized isoflurane was . The arterial blood oxygenation and the heart rate of the animals were monitored by a pulse oximeter ( , Nonin Medical, Plymouth, Minn.). In all the in vivo experiments, the scanning step size was along both the and axes. Each A-line was recorded for , which corresponded to a depth range of in biological tissues.
All experimental animal procedures were carried out in conformity with the guidelines of the National Institutes of Health.17 The laboratory animal protocol for this work was approved by the University Laboratory Animal Care Committee of Texas A&M University, where all experiments were carried out.
Skin Profile Detection
In PAM, modern edge-detection algorithms18, 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 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 smoothing20, 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 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 ) and the amplitude of PA signal generated from a subcutaneous blood vessel (denoted as ) can be classified into three major categorizes: (1) and are comparably strong [Fig. 2b]; (2) dominates [Fig. 2c]; and (3) dominates [Fig. 2d]. Figure 2b demonstrates the most ideal case as we can easily select the location of 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 and , 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 peaks in the first-step rough estimation for the 1-D skin profile in each B-scan image. 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 strongest PA peaks are identified and the distance between each pair of consecutive peaks ( , ) are further calculated. Once the largest pair distance is found, 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, was 2. In Fig. 2b, where and are comparable, the rough estimation gives the correct skin location. In Fig. 2c, where 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 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.
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 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,22 the nonparametric, robust, locally-weighted-regression smoothing20, 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”20 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 A-lines) was much shorter than the data acquisition time , 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 axis. The final 2-D skin profile was further acquired after a Gaussian low-pass spatial filtering.
After the 2-D skin profile was acquired from a flat scan, an optimal distance between the ultrasonic detector and the skin surface was chosen based on the ultrasonic focal length and the imaged mean vessel depth. One horizontal position of the ultrasonic detector, where the best ultrasonic focusing of a subcutaneous vessel was observed, was set as the reference position. The observation of the ultrasonic focusing is based on the imaged cross-section of the vessels in a B-scan image. Small objects (such as microvessels) demonstrate similar “butterfly” shape with symmetries along both the and axes when imaged within the ultrasonic focal region. For a more detailed description about ultrasonic focusing in PAM, please refer to Ref. 3. Distances from all other scanning positions to the reference position were then calculated along the axis. The distances were further converted to motor rotation steps to vertically translate the ultrasonic detector. In the end, the auto-fit scan was conducted in a way that the coordinate was adjusted at each horizontal position to maintain the preset optimal distance within the whole region of interest. Since the auto-fit scan requires the information from an initial flat scan, the data acquisition time is doubled.
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 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, 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.
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 . 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.
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.
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 , which was much larger than the ultrasonic focal zone . 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 density23 and vessel branching24 for further studies of cancer biology.
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, 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.
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.
This project is sponsored by National Institutes of Health grants R01 EB000712 and R01 NS46214 to L. V. Wang. We thank Dr. George Stoica and Dr. Meng-lin Li for their help in animal experiments. We also thank Mrs. Michelle Schoenecker for proofreading the manuscript.