1.

Introduction

Most brain imaging techniques, such as optical imaging of intrinsic signal (OIS) and functional magnetic resonance imaging (fMRI), rely on blood oxygen content for imaging neuronal activity.^{1, 2, 3} Thus, hemodynamic changes play an important role in investigating the coupling relationship between neurovascular response and functional brain mapping.⁴ Cerebral hemodynamics is considered to be closely related to several compartments, including arterial, venous, arteriolar, venular, and other capillary dynamics.^{5, 6}

To collect information on the detailed activity of cerebral vessels, researchers have constructed many vascular models by investigating the dynamics of the vascular compartments using optical imaging of intrinsic signal. Friston developed a nonlinear dynamic model of hemodynamic response to neural activity by linking the normalized deoxyhemoglobin changes to stimulation in the venous compartment.⁷ This model was extended by Zheng via construction of a three-compartment hemodynamic model, including the arterial, venous, and capillary compartments, for further study of the relationship between hemodynamic changes and neural activity.^{8, 9} Vanzetta and Grinvald developed the dual-wavelength imaging method to identify the microvascular compartments for analyzing the activity-dependent dynamics of the cortical blood volume and oximetry.⁵ Recently, vascular compartment dynamics has been examined with 3-D laminar optical tomography (LOT) and in vivo two-photon microscopy during somatosensory stimulation.¹⁰ More detailed information about hemodynamics has been revealed using these newly developed techniques. Schiessl employed independent component analysis (ICA) to study hemodynamic responses by separating the vascular compartments to improve the signal-to-noise ratio (SNR) in functional brain mapping.¹¹ They studied the distinct response dynamics of the arterial, venous, and capillary compartments by identifying each region using ex vivo vascular casts following the LOT data and two-photon image acquisition. Luo investigated the different responses between venules and arterioles by laser speckle imaging of the cerebral blood flow in microcirculation.^{12, 13} By separating the arteries and veins, it would be feasible to investigate in greater detail the spatiotemporal characteristics of arterial, venous, and capillary bed regions independently to determine their coupling relationships with neural activity using 2-D optical imaging.

Automatic separation of arteries and veins is also important for computer-aided diagnosis of vascular diseases, e.g., aneurysm, stenosis, and retinal arteriosclerosis.^{14, 15, 16} To enhance the visualization of the vascular structure, it is critical to segment the vascular network and to distinguish arteries and veins. However, due to the complexity of the vascular network, it is time consuming and impractical to do this manually, suggesting a need for robust and quick techniques to identify and separate arteries from the vascular network with a minimal amount of manual intervention.

According to our known literature, little work has been published to automatically identify arteries and veins in the cerebral cortex. A few papers have reported on artery-vein separation using retinal fundus image and magnetic resonance angiography (MRA). Yu and Akita identified the vessel types by utilizing the intensity difference between arteries and veins in a retinal fundus image.^{17, 18} A piecewise Gaussian model for profiling and differentiating retinal arteries and veins was proposed to describe the intensity distribution by considering the central reflex characteristic of the vessels.¹⁹ Grisan developed a so-called divide et impera strategy with color information for automatic discrimination of retinal vessels into arteries and veins.²⁰ However, these described methods utilize only the intensity profile and may be prone to errors in the optical image of the cerebral cortex. More recently, Narasimha-Iyer combined structural and functional features that are the ratio of optical densities obtained through dual-wavelength illumination to separate retinal arteries and veins with an improved result.²¹ In MRA and computed tomography (CT), several methods to separate arteries and veins have also been proposed for many potential applications in the head, neck, lungs, heart, abdomen, and lower extremities.^{22, 23} These separation strategies include acquisition methods, post-processing techniques, and image-processing methods.²² A software package called 3DVIEWNIX-AVS for artery-vein separation was developed to improve the vascular visualization of contrast-enhanced MRA images.²⁴ Although the above strategies have been successfully applied in 3-D MRA images, most of them, such as the phase-contrast and time-resolved acquisition approaches, may be unsuitable to apply to on 2-D optical imaging data because of their specificity related to MRA.

In this paper, we present an automated method for artery-vein separation in cerebral cortical images recorded with optical imaging of intrinsic signal. Differing from prior works, this artery-vein separation method utilizes the vessel dynamic information based on the fact that the power spectra of the intrinsic physiological oscillations distribute differently on arteries than on veins. First, we extracted the vascular network from one cortical image, then determined the spectral features of the intrinsic signals by the multi-taper method (MTM) and statistical F-test. Finally, we separated the arteries and veins using the fuzzy c-means (FCM) method and region-growing approach. These procedures are described in this paper followed by a discussion, in which we analyze the limitations of our method and characterize the phase distributions of the intrinsic signals.

2.

Materials and Methods

2.1.

Optical Imaging of Intrinsic Signal

Optical imaging of intrinsic signal (OIS) is a functional brain-mapping technique that is increasing in popularity. OIS measures neural activities by imaging the reflected and scattered light from the active cerebral cortex. OIS offers both a higher spatial and temporal resolution than that achieved by most other alternative imaging techniques, such as fMRI and positron emission tomography (PET), for imaging cortical functional architecture.^{1, 25, 26, 27}

In this study, all the image data were recorded using the IMAGER 3001 system (Optical Imaging Inc., New York). A schematic of the OIS system is shown in Fig. 1 . The camera system consisted of a charge-coupled device (CCD) array with high precision ( $12\text{-bit}$ depth, $\sim 60\text{-}\mathrm{dB}$ dynamic range) and a set of front-to-front $f1.2$ camera lenses with a diameter of $50\mathrm{mm}$ . The camera was mounted to an XYZ regulation holder that was fixed on a vibration-isolated table for steady focusing on the region of interest (ROI) of the cerebral cortex. The cortical surface was epiilluminated with green light $(546\pm 10\mathrm{nm})$ from a stabilized tungsten halogen lamp via light guides. This wavelength is one of the isosbestic points for the absorption of oxyhemoglobin and deoxyhemoglobin and can reflect cerebral blood volume changes more significantly than red light $(600\sim 630\mathrm{nm})$ in the visible spectrum. Also, imaging with this green light can maximize the contrast of the vasculature and provide morphological information for segmentation of the blood vessels. A personal computer was used to capture and store the image data. The size of one image could reach $1024\times 1024\text{pixels}$ in a factual imaging area of $12\times 12{\mathrm{mm}}^2$ (approximate $12\times 12\mu {\mathrm{m}}^2$ per pixel), which offered high-resolution images for the vessel-type separation procedure.

Fig. 1

Intrinsic optical imaging system. A stabilized halogen lamp provides the illumination source of a specific wavelength using a light filter. The high-resolution CCD camera has a $12\text{-bit}$ depth, $\sim 60\text{-}\mathrm{dB}$ dynamic range, and two reverse-connected camera lenses. The personal computer stores the collected image data.

2.3.

Extraction of the Vessel Network and its Centerline

The first step of the artery-vein separation was to segment the vasculature from one cortical gray image chosen arbitrarily from the image sequence. Though few vessel segmentation methods are proposed for the cortical images, numerous approaches for retinal vessel extraction have been presented in the literature.^{30, 31, 32, 33, 34} We extracted the vessel network using an extended method proposed by Li ³⁰ with a Gaussian-Laplacian filter and automated threshold selection.^{31, 35} The boundary of the extracted vessel was refined using the Gaussian matched filter constructed by the nonlinear Levenberg-Marquardt method.³⁰ Figures 2, 2, 2 show a retinal image and the segmented vessels by Li’s method and our initial segmentation method, respectively. Among the four leading algorithms in the retinal DRIVE database (see Table 1 ), this initial vessel segmentation method ensures that false vessels produce the lowest acceptable level of true positives, which is practicable for our artery-vein separation method.

Fig. 2

Extraction of the vasculature and its centerline network. (a) A retinal fundus image used in Ref. 30. (b) Vessels detected by method used in Ref. 30. (c) Vessels extracted by our initial vessel segmentation method. (d) Raw cortical gray image captured with the optical imaging system. (e) Vasculature segmented from the raw gray image. (f) Centerline of the vascular network obtained by morphological operation of thinning. The branching points are marked by a larger gray value.

Table 1

Performance comparisons of the initial vessel segmentation method and four other leading algorithms on the retinal DRIVE database. The performance comparison shows that our initial vessel segmentation method is practicable for the separation of arteries and veins, ensuring false vessels produce the lowest acceptable level of true positives.

Method	TPR	FPR	Ac
Al-Diri (2009)	0.7282	0.0449	0.9258
Martinez-Perez (2007)	0.7246	0.0345	0.9344
Staal (2004)	0.7194	0.0227	0.9441
Jiang (2003)	0.6478	0.0374	0.9212
Our initial method	0.6651	0.0189	0.9430

TPR: truth positive rate; FPR: false positive rate; Ac: accuracy. Here, the TPR, FPR, and Ac are computed depending on the pixel number of each class.

An extracted cortical vessel network denoted by VB is shown in Fig. 2, and the raw gray image can be seen in Fig. 2. To extract the vessel centerline from the vessel network, VB was first converted into a binary image. Then a morphological thinning operator³⁵ was applied to the binary image to obtain the vascular centerlines, as shown in Fig. 2. The branching and crossover points were detected and labeled on the centerline network. For each intersection point, the number of its neighboring pixels with logical 1 should be no less than three under the definition of eight-connectivity. All of the detected branching and crossover points were assigned with a value differing from the other pixel values in the centerline network. These intersection points segmented the vessels into pieces on which the artery-vein separation procedure was implemented. Any two neighboring vessel segments were judged whether they belonged to the same vessel based on the similarity of their tangential directions in the neighborhood of the noted intersection point.

2.4.

Feature Extraction in Frequency Domain

The time course for each pixel was transformed into the frequency domain to extract the spectral features. A mean time course for all image pixels is shown in Fig. 3 . If $x\left(n\right)$ denotes the average gray level of the $n$ ’th image, $n=1,\dots ,N$ , where $N$ is the number of the sampling image in each trial, then $N=T∕\Delta t$ , where $T$ is the duration of one trial. Using the MTM,^{36, 37} $x\left(n\right)$ was transformed into an expression of power spectral density (PSD) denoted by $P\left(k\right)$ , $0⩽k⩽N$ .

Fig. 3

(a) Mean time course for all pixels of the captured image sequence. (b) Power spectral density of the mean time course as estimated by the MTM algorithm in Matlab. RF and HF, as well as eigenfrequencies F1 and F2, are labeled. The bandwidth between the two black dashed lines indicates the $1–7.5\text{-}\mathrm{Hz}$ frequency range for computing the spectral power ratio values. The power summation in this spectral range is used as the denominator. (c) Results of statistical F-test spectral components showing. HF had a $>95\mathrm{\%}$ significance level (horizontal dashed line at the bottom), while RF, F1, and F2 had a significance level $>99\mathrm{\%}$ (dashed line at the top).

The MTM utilizes multiorthogonal data tapers known as a discrete prolate spheroidal sequence (DPSS) to estimate the spectrum in the adaptive weighted forms. Briefly, for a time sequence $x\left(n\right)$ , $n=0,\dots ,N-1$ , the frequency eigencoefficients can be computed by^{36, 37}

Eq. 1

$$y_k\left(f\right)=\sum _{n=0}^{N-1}x\left(n\right)v_n^{\left(k\right)}{e}^{-i2\pi fn},k=0,\cdots ,K-1,$$

where $\left\{x\left(n\right)\right\}$ is the sequence of data samples, $f$ is the frequency in standard units, and $v_n^{\left(k\right)}$ is the $k$ ’th DPSS data window or taper defined as the real, unit-energy sequences on $(0,N-1)$ . The number of the data tapers is determined by $K\approx \lfloor 2NW\rfloor$ , where $W$ is the bandwidth in standardized units. When the range of the spectrum is not too large, the multitaper estimate can be approximated as

Eq. 2

$$\widehat{P}\left(f\right)=\frac{1}{K}\sum _{k=0}^{K-1}{\left|y_k\left(f\right)\right|}^2.$$

For a periodic signal at frequency $f_0$ with amplitude $\mu$ , the least-squares estimate of $\mu$ is given by

Eq. 3

$$\widehat{\mu }\left(f\right)=\frac{\sum _{k=0}^{K-1}{y}_k\left(f\right)V_k(f-f_0)}{\sum _{k=0}^{K-1}{V}_k^2(f-f_0)}.$$

The statistical F-test is presented as follows:

Eq. 4

$$F\left(f\right)=\frac{\frac{1}{2}{\left|\widehat{\mu }\left(f\right)\right|}^2\sum _{k=0}^{K-1}{\left|V_k(f-f_0)\right|}^2}{\frac{1}{2K-2}{r}^2(f,\widehat{\mu }\left(f\right))},$$

where $\widehat{\mu }\left(f\right)$ is the complex amplitude at frequency $f$ estimated by ordinary least-squares regression; $V_k(f-f_0)$ , the $k$ ’th Slepian function, is the Fourier transform of the corresponding Slepian sequence; and $r^2[f,\widehat{\mu }\left(f\right)]$ is the minimized residual sum of squares, so the ratio of the energy in the F-test can be explained by assuming a line component to the residual energy. The standardized F variance-ratio test has 2 and $2K-2$ degrees of freedom.³⁷ If the $F$ value at frequency $f_0$ is larger than the given threshold, $f_0$ is thought to be a periodic component rather than the background noise.

Due to the limited frequency resolution and the different physiological sources for respiration and heartbeat,^{38, 39} signals below $1\mathrm{Hz}$ , known as low-frequency oscillations with a center component of $0.1\mathrm{Hz}$ ,³⁸ were not prepared for discussion in the following vessel separation procedure. In the frequency range of greater than $1\mathrm{Hz}$ , four significant linear spectral components were detected by spectrum estimation and statistical F-test at significance level $\alpha =0.05$ . These four components were HF, RF, the first significant frequency denoted as F1, and the second significant frequency denoted as F2. For an SD rat, RF is approximately $2\mathrm{Hz}$ and HF is approximately $5\mathrm{Hz}$ ; the actual values of RF and HF are dependent on experimental individuals and their physiological state.³⁹ For example, in the shown data, RF was approximately $1.6\mathrm{Hz}$ and HF approximately $5.7\mathrm{Hz}$ [see Fig. 3]. Based on the F-test results on each spectral component, the $F$ values indicated that there were significant linear components at the four frequency points [for HF, $p<0.05$ ; for RF, F1, and F2, $p<0.01$ ; see Fig. 3]. Additionally, the spectral power of these four linear components occupied over 10% of the gross power ranging from $1\text{to}7.5\mathrm{Hz}$ . Therefore, RF, F1, F2, and HF were all selected as the feature frequencies for artery-vein separation.

A spectral power-ratio calculation in the range of $1\sim 7.5\mathrm{Hz}$ was implemented at each pixel to reduce the effect of uneven illumination and to enhance the contrast between vessels and cortical parenchyma. Four power-ratio images corresponded to the four eigenfrequencies. Each power-ratio image was obtained by computing the ratio of the spectral power at the corresponding eigenfrequency to the summation of the spectral power ranging from $1\text{to}7.5\mathrm{Hz}$ at each pixel. The spectral amplitudes and power-ratio values constructed four couples of new gray intensity images, which were termed the spectral amplitude images and power-ratio images of RF, F1, F2, and HF, respectively (see Fig. 4 ).

Fig. 4

(a) Normalized spectral amplitude 1, 2, 3, and 4, and (b) normalized power-ratio distributions 1, 2, 3, and 4 of RF, F1, F2, and HF, respectively $(\mathrm{AMP}=\text{amplitude})$ .

Both sides of most veins had higher RF amplitude than in venous central region. The RE is considered to reflect the periodical changes of the cerebral blood volume in the veins. The statistical significance of RF was also an important factor for choosing RF as one of the eigenfrequencies. The frequency of F2 behaved much stronger in the capillaries and thin vessels than in large vessels. The amplitude of F2 was relatively small in the main veins, and its power ratio was relatively small in both arteries and veins. The F2 frequency was considered to mainly reflect the red blood cells (RBCs) that traveled discretely in the vessels observed under OIS. This phenomenon was in accordance with observations by Vanzetta,⁴⁰ and can be seen clearly in the recorded video after being bandpass filtered. Due to the thick walls of the main blood vessels, RBC movement may not be seen in the large arteries and veins using optical imaging. The amplitudes of F1 and HF were significantly different in the arteries than the veins. Both F1 and HF (especially HF) were observed to be stronger in the arteries than in the veins. The power-ratio image of HF revealed more differences in amplitude intensity distribution between the major arterial and venous regions, which provided more detailed information for the separation of arteries and veins. Therefore, both the amplitude and power ratio of HF were chosen as the elements of the feature vector for artery-vein separation. The arteries moved and oscillated in a different pattern with both sides of the arteries exhibiting higher HF amplitudes than the central regions, which presumably was caused by the vasomotion and transverse oscillation section by section at the frequency of the heartbeat. The eigenfrequency of F1 presented partly the same physiological response as HF. F1 had larger amplitude intensities in the arteries than RF and F2. The optical intrinsic signals would be a mixture in the artery-vein overlapping place, which must be taken into account in the artery-vein separation procedure.

The cerebral cortex consists of three compartments—arteries, veins, and the capillary bed—according to the spectral distribution of the four eigenfrequencies [see Figs. 4 and 4]. For each compartment, four typical rectangular regions were manually selected with square regions, and the mean pixel sequence of each area was then transformed into the spectral domain by MTM to compare their spectral differences (see Fig. 5 ).

Fig. 5

Comparison of the power spectra in arteries, veins, and capillary bed. (a) Illustration of the manually selected typical feature regions with three class regions located in the arteries, veins, and capillary bed. Each of the four regions in the same compartment is marked with the same color. (b) Power spectra estimation of the pixel sequences in the arteries, veins, and capillary bed. All the spectra were estimated by the MTM algorithm.

Figure 5 shows that the spectral magnitudes of HF and F1 were much higher than that of the other frequencies in the $>1\text{-}\mathrm{Hz}$ range in the arterial regions. In the venous regions, the spectral magnitudes of RF were higher in the $>1\text{-}\mathrm{Hz}$ frequency range, and the absolute spectral amplitudes of HF were comparable to RF. In the capillary bed, HF and F2 were the two dominant frequencies. Therefore, RF and HF were considered to be the main factors for identifying arteries and veins. In the F-test values, HF was less significant than RF, F1, and F2, which was probably because the arterial compartment occupied a small area in the ROI of the cerebral cortex.

2.5.

Artery-Vein Separation Using the Fuzzy C-Means Method

To separate the arteries and veins, the FCM clustering method⁴¹ was first used to determine the arterial seeds. Considering the uneven spectral power distribution of the intrinsic physiological signals, the seeded region-growing method⁴² was then applied to form integral vessel pieces from the classified seeds on the centerline of the vascular network. FCM clustering is an elementary but popular approximate clustering approach. It can be used to classify a data set with $N$ items into $c$ classes. The principle of this clustering approach is described as follows: Mathematically, the objective function of FCM for partitioning a dataset ${\left\{x_j\right\}}_{j=1}^N\subset R^d$ into $c$ clusters can be given by

Eq. 5

$$J_m(U,\Theta )=\sum _{j=1}^N\sum _{i=1}^c{u}_{ij}^m{\left|x_j-{\theta }_i\right|}^2,$$

where $U={\left(u_{ij}\right)}_{c\times N}$ is the fuzzy membership degree matrix that satisfies

$$U∊\{u_{ij}∊[0,1]|\forall i,0<\sum _{j=1}^N{u}_{ij}<N|\forall j,\sum _{i=1}^c{u}_{ij}=1\},$$

and $\Theta =\{{\theta }_1,{\theta }_2,\cdots ,{\theta }_c\}$ stands for the centroids of the clusters. The notation $|\mathbf{\cdot }|$ stands for the Euclidean norm, and $N$ is the number of elements in the dataset. The parameter $m$ represents the weighting exponent on each fuzzy membership and can be arbitrarity selected from the domain of $(1,\infty )$ .

The FCM clustering approach was used to classify the arterial seeds from the four couples of spectral amplitude and power-ratio images. A morphologically clean filter was then performed on the clustered image to remove the point noise. Subsequently, an intersection operator was conducted on the vascular centerline to obtain the seeds of each compartment. The seeded region-growing method was applied to the vessel segments to reconstruct the arteries and veins based on the classified seeds. Two criteria were defined in the growing region:

• 1. The homogeneity criterion: If the intensity of the candidate pixel was in the tolerance region of the average intensity of the premerged region, the candidate pixel was merged. In our procedure, the vessel network was converted to a binary image. So if the candidate pixel was equal to logical 1 and had a neighboring seeded pixel with eight-connectivity, it was merged into the region.

• 2. The local stop criterion: If the neighboring candidate pixel was the crossover or branching point, the growing process was stopped.

The five steps of the separation procedure, with the arteries, veins, and vascular network denoted by $A$ , $V$ , and $VN$ , respectively, are described below:

Step 1: Extract the vascular network from the grayscale cortical image and obtain the centerline of the vascular network with the method described in Sec. 2.3.

Step 2: Determine the arterial seed sets by the FCM clustering method. After the clustering returns, pixels in the cluster with maximum mean spectral power of HF are judged to be the arterial parts. The pseudo-code of the FCM method is given below.

Fuzzy c-means clustering
Begin initialize the cluster number $c$ , the parameter $\mathit{m}$ , here $\mathit{m}=2$ , and an appropriate small positive number $\epsilon >0$ , determine the initialized fuzzy membership matrix ${\mathit{U}}^{\left(0\right)}$ , set $s=0.{\mathit{U}}^{\left(0\right)}$ is initialized by randomly generation. The sum of each $c$ columns of ${\mathit{U}}^{\left(0\right)}$ is equal to unity, as required by FCM.
Do (1) Compute ${\Theta }^s$ depend on ${\mathit{U}}^{\left(s\right)}$ :
$${\theta }_i^{\left(s\right)}=\raisebox{1ex}{$\sum _{j=1}^N{u}_{ij}^m{x}_j$}\!\left/ \!\raisebox{-1ex}{$\sum _{j=1}^N{u}_{ij}^m$}\right.,i=1,2,\dots ,c$$
(2) Renew the matrix value ${\mathit{U}}^{\left(s\right)}$ to ${\mathit{U}}^{(s+1)}$
$$u_{ij}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sum _{k=1}^c{\left(\frac{d_{ij}}{d_{kj}}\right)}^{2∕(m-1)}$}\right.$$
Until $\parallel U^{\left(s\right)}-U^{(s-1)}\parallel <\epsilon$ , or $s=s+1$ .
Return $\mathit{U}$
End

Step 3: Find the arterial seeds on the centerline of the vascular network by applying the morphological operator of intersection ‘∩’.

Step 4: Run the region-growing procedure to obtain the integral arterial centerline segments. The pseudo-code of region growing for arteries is given below.

Region growing
Begin Define $\mathit{S}$ as the seed set, $\mathit{A}={\mathit{0}}_{R\times C}$ as the arterial centerline segments, $R$ and $C$ are the row and column of the gray image respectively.
Do start from each element in $\mathit{S}$ , if the homogeneity criterion is satisfied for a candidate pixel, set its value to 1 in $\mathit{A}$ and merge this pixel into $\mathit{S}$ .
Until no new pixel is added
End

Step 5: Execute the morphological operation of dilation on $A_{R\times C}$ . The size of the dilation kernel was set to $7\times 7\text{pixels}$ depending on the estimation of the maximum vessel diameter. The arteries and veins were then obtained by $A=VN\cap A_{R\times C}$ and $V=VN\backslash A$ , respectively. Here, the notation ‘\’ represents the difference between the two sets.

The FCM clustering method was implemented on the feature vectors consisting of the four couples of eigenspectra images on the vasculature. After considering the artery-vein overlap and their relative location, i.e., whether the artery was located on the top of the vein or vice versa, we set the number of clusters to four in the FCM procedure. Implementing the clustering on the vasculature has the advantage of reducing the capillary disturbance to the arteries and veins, but it has the disadvantage of losing some vessel segments, which were absent from the extracted vascular network. We also tried to cluster the spectral distributions in the whole image. However, six or more classes must be defined to separate arteries and veins efficiently, taking account of the additional components such as arterioles, venules, and capillary regions for interchanging oxygen and nutrients. If an active area is stimulated, more predefined clusters are required. The FCM algorithm can converge to the same final clustering result when given different initial membership matrices in different runs. Figure 6 illustrates the convergence of FCM with different initial membership matrices on one image data point.

Fig. 6

Convergence validation of FCM method with different initial membership matrices on one image data point. (a) Initial membership matrices in six runs. (b) Corresponding final clustering membership for each cluster.

3.

Results

A sampling cortical image $(80\times 80\text{pixels})$ sequence was first adopted to validate the separation method. This sequence contained three main pieces of vessels, including both arteries and veins. Arterial seeds were obtained using the FCM clustering procedure from the amplitude and power-ratio mappings of RF, F1, F2, and HF [see Figs. 7 and 7 ], and then the arteries and veins were identified from the vascular network using the region-growing method. The separation results are shown in Fig. 7, where the red parts are arteries and the blue parts are veins. Figure 7 shows the raw gray image.

Fig. 7

(a) Spectral amplitude distributions of RF, F1, F2, and HF. (b) Power-ratio distributions of the four eigenfrequencies. (c) Raw cortical image. (d) Artery-vein separation result. The red and blue regions represent arteries and veins, respectively. (Color online only.)

The separation results from the sampling cortical image validated the five steps of the separation procedure described above. We applied this procedure to the other cortical images with more arteries and veins. Figure 8 shows one example result of the separated arteries and veins. Four veins labeled V1 through V4 and five arteries labeled A1 through A5 [see Figs. 8 and 8] were extracted by the FCM and region-growing methods using the spectral amplitude and power ratio images of RF, F1, F2, and HF, as well as the corresponding vascular network. All the separated arterial and venous segments, plotted in red and blue, respectively, are shown in Fig. 8. The separation result with the arterial and venous centerline superimposed on the raw cortical image is shown in Fig. 8. We marked two typical misjudgments of the vessel types as C1 and C2 [see Fig. 8] and will discuss them below. Ground-truth data were obtained from the 12 cortical gray images by the physiological expert, who manually labeled significant vessels and indicated their types. The vessels whose type exhibited ambiguity were excluded in the ground-truth data. The separation results were compared against the ground-truth data to estimate the performance. The true positive rate (TPR) was computed by the total length of the ground truth arterial/venous segments divided by the total length of the correctly classified arterial/venous segments. The sum of the true positive rate and misclassification rate was 1 for one vessel type.²¹ The TPR of our separation method reached 97.7% for the arteries and 98.4% for the veins. The misclassification of arteries as veins was 2.3%, and the misclassification of veins as arteries was 1.6%.

Fig. 8

Example result of artery-vein separation using the FCM clustering approach and region growing method. (a) Venous centerline segments labeled as V1, V2, V3, and V4. (b) Extracted arterial centerline segments labeled A1 through A5. (c) Artery-vein contrast image: red parts represent arteries, blue parts represent veins, and cyan pixels show the overlapping parts of arteries and veins. (d) Superimposition of arteries and veins on the raw gray image. The two typical limitations in the artery-vein separation results are denoted by C1 and C2. The points marked by C1 represent the first limitation where arterial fragments belonging to the same artery were divided into pieces by larger veins. The points marked by C2 represent the second limitation, a misjudgment of arterial segments. (e) Estimated power spectrum in arterial, venous, and capillary bed regions. The black line represents power spectra in the capillary bed, the blue line represents power spectra in venous regions, and the red line represents power spectra in arterial regions. (Color online only.)

Following the artery-vein separation, the power spectral density of the pixel sequences in the arteries, veins, and capillary bed were estimated by MTM to compare their oscillational distribution characteristics in the frequency domain. In the arterial regions, the power spectra of F1 and HF were significantly larger than those in the venous regions, while the power spectra of RF and F2 had greater values in the veins [as shown in Fig. 8]. The three mean time courses were normalized before the spectrum estimation. The spectrum in the capillary bed displayed a median value compared with those in the arterial and venous regions, which may be explained by the contribution of the remaining arterioles and venules. Overall, these three spectral densities suggest that the different oscillational characteristics may provide useful information for the detailed hemodynamic analysis by separating arteries and veins.

The veins separated by using the venous maximal RF power in the region-growing procedure are also shown in Fig. 9 for comparison. Some small veins were missed with this strategy, which presented a relatively weak performance in contrast with the approach above. Hence, we chose the subtraction strategy to obtain the veins.

Fig. 9

Comparison of the separated veins with two strategies. (a) Venous centerline extracted by the region-growing strategy. (b) Venous centerline extracted by the subtraction strategy.

4.

Discussion

The artery-vein separation method used in this study had two limitations, as marked in Fig. 8. First, the area enclosed by C1 contained typical artery and vein crossover points, termed break points. In C1, the arterial fragments, which factually belonged to the same artery, were divided into two more parts by the thick veins lying over the arteries. In some situations, it is important to connect the adjacent arterial fragments that are divided by the veins for functional brain mapping analysis or hemodynamic modeling. In this case, because the directions of the two arterial parts were very close, it was feasible to link the two parts together using an image-processing method. However, if an artery travels below the vein for a distance and its direction alters significantly on the other side of this vein, it would be a challenging problem to link the arterial fragments together unless certain heuristic information has been given. The second limitation of the vessel-type separation is illustrated in the area marked by C2. C2 was located at the crossover regions of the arteries and veins with abundant arterioles. One vessel part at the inner top of C2 was judged to be an artery and the other part at the inner bottom of C2 to be a vein. In C2, there were many small arterioles and venules located nearby, and their inevitable interference with the neighboring major arteries and veins likely resulted in misjudgments of the vessel types.

We attempted to separate the arteries and veins using the phase distribution of the four eigenfrequencies, as shown in Fig. 10 . In the phase distribution of RF, both the arteries and veins exhibited a bilateral antiphase, which provided little contribution to the vessel-type separation but was useful for vessel extraction. Fortunately, phases of F1 and HF exhibited significant differences between the regional arteries and veins. The phase maps of F1 and HF depicted a bilateral antiphase distributed along all the vessels. However, the phenomenon of bilateral antiphase in the arteries was altered section by section, i.e., the bilateral phases in one section were the reverse of the neighboring section, which was likely caused by the arterial wiggle. In contrast, the phases in the veins were maintained almost constant throughout all sections along the two sides of the vessels. This was mainly caused by the RF vasomotion. Therefore, based on these features, we could determine whether the vessel was an artery or not. However, the phase distribution of F2 provided information mainly in small vessels; due to the limited spatial resolution, the very small vessels were difficult to correctly identify.

Fig. 10

Spectral phase maps of (a) RF, (b) F1, (c) F2, and (d) HF.

The optical intrinsic signals captured from rats anesthetized by pentobarbital sodium also had four significant frequencies in the domain of $1\sim 7.5\mathrm{Hz}$ (data not shown), which is in accordance with the spectral distribution of signals from rats anesthetized by urethane. Thus, overall the results from the present study suggest that the four eigenfrequencies rely on only the physiological state of the animal, and are intrinsic physiological oscillations. The RF and HF were relatively stronger in veins and arteries than in other small vessels, whereas F1, and especially F2, exhibited a relatively high amplitude in small vessels, including some in the deeper layers of the cerebral cortex. In this sense, RF and HF were the dominant elements for the separation of arteries and veins. The duration of image acquisition was set to $40\mathrm{s}$ in this work, which ensured an adequate frequency resolution for RF, F1, F2, and HF. However, the $40\text{-}\mathrm{s}$ duration was not critical for artery-vein separation. A shorter duration can be feasible as long as it provides the proper frequency resolution. The minimum duration would vary for individuals and their physiological state during the experiments.

The proposed method is distinct from prior work applied to retinal imaging and MRA. Almost all the prior work identified arteries and veins by utilizing the color or intensity information of the vessels. The performance of the dual-wavelength method²¹ was reported to be 97% and 90% accuracy for the retinal arteries and veins, respectively. The performance of our method reached 97.7% and 98.4% for the cerebral cortical arteries and veins, respectively, by utilizing the signal spectral features. We believe that the performance of our artery-vein separation method could be further improved under multiple-wavelength illumination. Nevertheless, applying our proposed method to retinal images would be a challenging task due to the problem of eye motion and blinking during the acquisition of a series of images.

Though our method works on vessels with considerable width, its success is limited on very thin arterioles and venules. Very thin vessels with a width of one pixel $(\sim 12\mu \mathrm{m})$ or so are difficult to extract from the cortical gray image. For thin vessels, a successful classification method depends on the reliable extraction of the vessel hierarchy in vasculature and accurate labeling with their parent vessels.

5.

Conclusion

The artery-vein separation method proposed in the present study was achieved using the spectral amplitudes and power ratios of the intrinsic physiological oscillations of RF, F1, F2, and HF. Each of the four feature frequencies was significant (F-test, $p<0.05$ ). The amplitudes and power ratios of the four frequencies in all the pixels constituted the feature vector, which was clustered using the FCM method to find the arterial seeds. By combining the extracted vessel network, this method was able to separate arteries and veins by 2-D optical imaging of intrinsic signal. The separation method was useful for hemodynamic modeling in different vessel types and may be valuable for cerebral surgery or diagnosis of vascular diseases by arterial morphological analysis. This method might also be extended to the identification of vessel types in depth-resolved optical imaging and focusing deeper into the cerebral cortex. A potential future study is to investigate the distributions of other intrinsic physiological signals, such as low-frequency oscillations, and the analysis of their power spectral differences in arteries and veins with optical imaging of intrinsic signal.