CC BY
0Contents lists available at ScienceDirect
NeuroImage
journal homepage: www.elsevier.com/locate/ynimg
Automatic anatomical labeling of the complete cerebral vasculature in mouse models
Sahar Ghanavati *, Jason P. Lerch, John G. Sled
Department of Medical Biophysics, University of Toronto, Toronto, Ontario M5G 2M9, Canada Mouse Imaging Centre, The Hospital for Sick Children, 25 Orde St., Toronto, Ontario M5T3H7, Canada
ARTICLE INFO
ABSTRACT
Article history: Accepted 15 March 2014 Available online 28 March 2014
Keywords:
Cerebrovascular anatomy Automatic segmentation Anatomical labeling Attributed relational graph Maximum a posteriori Stochastic relaxation Classification Micro-CT Mouse
Study of cerebral vascular structure broadens our understanding of underlying variations, such as pathologies that can lead to cerebrovascular disorders. The development of high resolution 3D imaging modalities has provided us with the raw material to study the blood vessels in small animals such as mice. However, the high complexity and 3D nature of the cerebral vasculature make comparison and analysis of the vessels difficult, time-consuming and laborious. Here we present a framework for automated segmentation and recognition of the cerebral vessels in high resolution 3D images that addresses this need. The vasculature is segmented by following vessel center lines starting from automatically generated seeds and the vascular structure is represented as a graph. Each vessel segment is represented as an edge in the graph and has local features such as length, diameter, and direction, and relational features representing the connectivity of the vessel segments. Using these features, each edge in the graph is automatically labeled with its anatomical name using a stochastic relaxation algorithm. We have validated our method on micro-CT images of C57Bl/6J mice. A leave-one-out test performed on the labeled data set demonstrated the recognition rate for all vessels including major named vessels and their minor branches to be > 75%. This automatic segmentation and recognition methods facilitate the comparison of blood vessels in large populations of subjects and allow us to study cerebrovascular variations.
© 2014 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license
(http://creativecommons.org/licenses/by/3.0/).
Introduction
The cerebrovasculature is a network of blood vessels including arteries, veins, sinuses, arterioles, venules and capillary beds, which supplies blood to the brain. The malfunction of this system can lead to morbidity and death. Many of the degenerative diseases such as Alzheimer's disease have been proposed to have an underlying vascular etiology (Bell and Zlokovic, 2009; de la Torre, 2004; de la Torre and Stefano, 2000). The characteristics of the cerebrovasculature may affect the probability of occurrence of such disorders. However, our understanding of the development, structuring and anatomical variations of the cerebrovascular system is limited. Phenotyping the cerebral vasculature and analyzing the vascular variations are essential for understanding genetic factors in the development, characteristics and diseases of the cerebrovasculature.
Due to the high variability of the cerebrovasculature, inter-subject comparison in a large number of subjects is required in order to characterize variations of the blood vessels. A method for studying and comparing cerebral anatomical structures is to generate an atlas and registering new subjects to it, in order to identify the corresponding structures in
* Corresponding author at: Department of Medical Biophysics, University of Toronto, Toronto, Ontario M5G 2M9, Canada.
E-mail address: sahar.ghanavati@mail.utoronto.ca (S. Ghanavati).
each individual. This approach has enabled a comparative study of cerebral structures across multiple subjects (see Dorr et al. (2008) for example). However, unlike neuroanatomical regions, where there is homology between different subjects, the differences between cerebral vasculature of any two subjects are so large that common registration methods fail to correct for such differences. These differences include the number of branches, size, curvature, and branching locations in the vascular network and can include the absence of a major artery (Beckmann, 2000; Kitagawa et al., 1998).
There is a rich literature on automatic methods of vascular segmentation (Babin et al., 2013; Caselles et al., 1997; Forkert et al., 2013; Frangi etal., 1998; Fridman etal.,2004; Krissian et al., 2000; Lorigoetal., 2001; Piccinelli et al., 2009; Qian et al., 2009; Shang et al., 2011; Sofka and Stewart, 2006). A review of different approaches of vessel segmentation can be found in Suri et al. (2002), Lesage et al. (2009). However, only a few automated methods for labeling vasculature have been investigated in the past. Graph matching algorithms have been employed to identify and label the vascular systems such as the airway trees (Tschirren et al., 2005) and the bronchial branches (Mori et al., 2000, 2009). In most of the early works such as the labeling of the bronchial branches, a knowledge-based algorithm based on fixed topological constraints on the bifurcations was used, which is likely to fail in the case of large variations from the reference such as missing vessels (Mori et al., 2000,
http: //dx.doi.org/10.1016/j.neuroimage.2014.03.044
1053-8119/© 2014 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/3.0/).
2009). Tschirren et al. proposed an automated labeling of the airway trees where the branches are first segmented and skeletonized from the CT images and are represented by a directed acyclic graph (DAG) (Tschirren et al., 2005). A branch point matching algorithm is then performed in order to find corresponding major branches in different data sets. In order to perform automatic labeling, an average labeled tree is generated from manually labeled trees which represent the geometrical and topological features of the bifurcations. Labels are assigned to the target tree by matching it to the average labeled hierarchical tree using association graphs. In another approach, supervised machine learning algorithms were used to automatically identify and label the abdominal arteries (Mori et al., 2010). For each vessel branch, features including the direction, length and diameter are extracted and the main abdominal aorta is found by the largest diameter. A pair-wise labeling is then performed based on the local features of the branches. Despite the high recognition rate of these tree matching algorithms, their applications are limited by their dependence on the topology of hierarchical structures.
There has been only a few works investigating the labeling of the cerebral arteries. Maximum a posteriori probability (MAP) estimation was used to automatically label the circle ofWillis (Bogunovic et al., 2013). In this approach, the arteries in the circle ofWillis are segmented from the magnetic resonance angiography and are modeled as a rooted attributed graph. Local features, such as the radius and direction of the vessels, are computed for each bifurcation. Topological information is defined based on comparison with a reference labeled graph. The optimal labeling is formulated as a MAP. They labeled 11 bifurcations of the circle ofWillis with 95% accuracy. In another approach, Bilgel et al. used a combination of random forests on local features of vessels with a belief propagation of vascular tree connectivity to label the anterior cerebral arteries in 3D time of flight MR angiography (TOF MRA) (Bilgel et al., 2013). The application of these methods has been demonstrated for a few cerebral arteries and may not generalize to the more complex graph of the complete cerebrovasculature. Specifically, the labeling method used by Bogunovic et al. (2013) is based on finding maximal cliques between the graph and a reference graph. Constructing a reference graph such as that constructed for the circle ofWillis is an elegant approach for labeling sub-trees in the cerebrovasculature. However, the high variations in the branching of the cerebrovasculature make it very complicated to construct an encyclopedic reference graph for the complete cerebrovascular graph. To the best of our knowledge, no previous work has been published on the labeling of the complete cerebrovasculature.
The difficulty of high resolution imaging of cerebral vessels of laboratory animals such as mice, the high complexity of the cerebrovasculature and the density of vessels in the brain makes the study of this interesting
anatomy challenging. As a result, there has been little detailed description of the whole cerebrovasculature of mice. A few reference works are the mouse cerebrovascular atlas by Dorr et al. (2007), and the description of the mouse cerebrovascular system detailed by Scremin and Holschneider (2011). Developing a framework to automatically recognize and label all the vessels in the mouse cerebrovascular system will enable the study of large populations in greater detail and can lead to a better understanding of the development, phenotypes and diseases of vasculature in the brain. In this paper, we detail the development of such a framework.
Methods
Our automatic segmentation and recognition of the vasculature consist of sub-processes as shown in the pipeline in Fig. 1. First, an ex vivo high resolution micro-CT image is acquired from the complete cerebrovascular system of the mouse brain filled with contrast agent. The details of preparation and imaging of the specimens can be found in Ghanavati et al. (2014). The images are then registered to a common space by manually selecting landmarks on the major vascular branches as described by Chugh et al. (2009). The common space used is the vascular atlas created by Dorr et al. (2007). The registered images are normalized so that the background and foreground intensity of all the subjects will be in a similar range. The original CT scans contain cerebral vessels as well as facial arteries and veins and upper jaw. Therefore, the brain in each scan is delineated manually using the CT-Analyser software (Bruker micro-CT, Belgium), so that final images only contain the cerebral blood vessels. Where an MRI scan has been acquired in addition to micro-CT, this step can be replaced by an automatic method using a brain mask constructed from the anatomical MRI and registered to the CT scan. Next, the cerebrovasculature is automatically segmented and is represented as a connected graph. Finally, each graph is automatically labeled using a reference set of labeled cerebrovasculatures (the training set) through a stochastic optimization scheme. In the following sections, we will describe the steps of the pipeline in detail.
Vascular segmentation
The first step in the automatic labeling of vasculature is to segment the vessels and to represent them as a graph. The vessel segmentation was fully automated. The method closely follows the tubular object extraction method of Fridman et al. (2004) and is detailed in Rennie et al. (2011). We have modified the method to eliminate the need for manual definition of seeds required for the initialization of the segmentation algorithm. The distance transform of a binarized version of the image is calculated and used to place an initial seed at the center of the largest
Classification Phase Fig. 1. The pipeline for the automatic segmentation and labeling of the cerebrovascular system.
vessel. After the completion of the first iteration of vessel segmentation, a binary image is created from the segmented tree and its negation is combined by logical-and with the binarized original image obtained by thresholding. The threshold value is chosen manually, half way between the average intensity of Microfil and the average background intensity. The distance transform of this new binary image is computed and the seed selection process is repeated once for each connected component of the binary image. The vessel segmentation is repeated with new seeds until no new seed can be added. It is important to note that the tracing of vessels is performed on the original gray level image once the seeds have been identified. The algorithm returns one or more connected graphs whose vertices (nodes) are located at the center line of vessels (Fig. 2). In the present work 97 ± 1% of vessel segments belonged to a single connected component. As the remaining disconnected components consisted of a small number (< 8) of short segments (< 10 voxels in length), these were neglected in subsequent analysis. Each vertex has features including the radius, the Cartesian spatial location, and the direction of the normal vector of the vessel cross-section which pass through the vertex. Using the vertex attributes, a tubular model of the segmented vasculature can be constructed in which the mean diameter of the vessels can be demonstrated, as shown in Fig. 2. These features are used to convert the segmented vasculature into an attributed graph G = (V, E), in which vertices V represent the vessel branching points and edges E represent the vessel segments contained between two branching points or a branching point and a terminal node (where the graph segment ends without further bifurcation) and have local features including the length and diameter.
Manual labeling
In order to perform the automatic labeling using supervised machine learning algorithms, we needed a set of training data, which also provides a ground truth for the validation of the method using a
leave-one-out (LOO) test. As a result, we have developed an interactive visualization software "Brain-view" using Coin-3D Graphics Development kit (Kongsberg Oil & Gas Technologies, Sandvika, Norway) and Qt Develop Cross-platform application and UI framework (Nokia, Oslo, Norway), in order to visualize the segmented vascular graph in 3D and to enable the labeling of each vessel segment by a click of the mouse. The open source software can be downloaded from https://github.com/sghanavati/brain-view2. The software allows the labeling of each vessel segment by one click on the segment and assigns a unique numerical and color code for each chosen label to the labeled vessel segment. As a result, the manual labeling can be performed in a rapid and intuitive fashion. The labels from the reference labeled micro-CT by Dorr et al. (2007) can be propagated onto the graph segments visualized in Brain-view, to be used as an initial estimate to guide the manual labeling. An example of a manually labeled sample is shown in Fig. 3. Some of the major arteries and veins are identified in the 3 views of the sample.
Features
Similar to Bilgel et al. (2013) and Bogunovic et al. (2013), we calculated two types of features on the graph: local and relational. For each of the edges (vessel segments) in the segmented graph, a set of local features was extracted. The local features consisted of the mean diameter, length, curvature, tortuosity, orientation, midpoint position in 3D Cartesian space, and the anatomical regions containing the vessel. Geometrically, the curvature of a line is defined as R, where R is the radius of the osculating circle of the curve at the midpoint. For simplicity, we calculated R based on the distance between the midpoint of the edge and the median point, where the median point is the point with equal distance from each end of the edge located on a straight line connecting the two ends of the edge. Tortuosity (t) is defined by the ratio of the length of the curve (L) to the distance between the ends of it (C): t = £.
Fig. 2. Overview of the automatic vessel tracking results in sagittal view. (a) The maximum intensity projection of a micro-CT image, (b) the center line of the segmented vasculature, (c) the 3D iso-surface model constructed from the micro-CT image using marching cubes, and (d) the tubular object generated by the automatic segmentation color coded by the vessel diameter attribute.
Fig. 3. The cylindrical representation of the segmented vasculature and the color-coded labeling. Top left is the sagittal view. Top right is the dorsal view and the bottom is the ventral view of the sample. Some of the major arteries and veins are identified on the image. Specifically, the arteries forming the circle of Willis can be seen at the bottom image. The arteries include internal carotid arteries (ICA), middle cerebral arteries (MCA), anterior cerebral arteries (ACA), olfactory artery, superior cerebellar arteries (SCA), basilar artery, and vertebral arteries. The veins depicted include transverse sinus and superior sagittal sinus.
Orientation of each vessel segment was defined by a unit vector connecting the first and last vertices of the vessel segment, in the +Z semi-hemisphere in the Cartesian space. All the micro-CT vascular images are registered by manually selected landmarks to the mouse brain anatomical atlas developed by Dorr et al. (2008). The proximity of the midpoint of the vessel to each anatomical region was calculated. Relational features take the connectivity of the edges into account. Assuming a Markovian property on the graph, we only use the immediate adjacent edges to define the connectivity feature. For each edge, four immediate neighbors are assumed, two at each end. For consistency, in the case of a terminal edge, we consider it to be adjacent to two null edges. The relative diameter of the current edge with the adjacent edges and the angle it makes with each of its adjacent edges were retained as relational features. For the terminal edges, the relative diameter is set to 1 and the angle is set to 0. Finally, an adjacency matrix is formed using the training set, which defines the probability of two adjacent vessel segments having a given pair of labels. For a set of size N for possible labels, the adjacency matrix is N x (N + 1), whose elements give the probability of each label being adjacent to all N possible labels or to the null edge. For example, the (m, n)-th element of the adjacency matrix is p(lm\ln), indicating the probability of m-th label being adjacent to the n-th label in the whole training data set.
Stochastic relaxation
The problem of automatic recognition and labeling of the cerebro-vasculature has many similarities to the automatic labeling of cortical folding patterns (Mangin et al., 1995; Riviere et al., 2002). Both the sulci and cerebral vascular systems can be represented by an attributed relational graph (ARG), where each edge is a random variable and local
features are presented as attributes of each corresponding edge (Di Ruberto et al., 2002). The labeling problem can be represented as a Markov random field (MRF), assuming that the posterior probability of the label of each edge depends only on a closed neighborhood of that edge. Fig. 4 shows the Markovian property assumption on some arteries. As it was shown for the labeling of the circle of Willis (Bogunovic et al., 2013), a probabilistic Bayesian framework can be formed on the graph of the cerebral vasculature and the labeling problem can be formulated as a MAP estimator on the local and relational features of the graph edges, and the optimum configuration is achieved by finding
Azygos artery
Fig. 4. The Markovian assumption illustrated on a sub-tree of the cerebral arteries. The MCA segment highlighted in red has exactly four immediate neighbor segments, ICA, ACA and two MCA branches. Under the Markovian property assumption, the posterior probability of labeling the indicated segment as MCA depends only on the labels of its immediate neighbors and not on the labeling of a non-adjacent neighbor such as the azygos anterior cerebral artery.
argmaxLp(L\G), where L is the set of assigned labels for the given graph G, and p(L\G) is an aggregation of posterior probabilities of all edges in the graph given the set of labels L.
In accordance to the definition of ARG, we assume that each feature is an independent and identically-distributed (iid) random variable, and is represented by a univariate Gaussian distribution. The direction feature consisting of 3 components of the direction unit vector is estimated by the Kent distribution (5-parameter Fisher-Bingham distribution), which is equivalent to the Gaussian distribution on a unit sphere (Kent, 1982; Mardia, 1972; Mardia and Jupp, 2009). A Bayesian framework is employed to assign labels to the edges of a test data set using the estimated distributions from the training data set. For each edge ej in the training data set with the feature vector <&ej, a posterior probability for each of the possible labels is calculated:
pttljp^ |li)p(li) (1)
where U is the i-th possible label, 4y is the feature vector of edge j. The first term is the likelihood of the label lt for the given set of features 4>e. and the second term is the prior probability of the given label.
The likelihood of the label U for the given set of features <&ej is calculated as follows:
pK lk) = n Pd (2)
The labeling of the edges is also dependent on their adjacent edges, based on the Markovian property of the graph. As a result, the posterior probability of a given label for a given edge is updated as follows:
p je,)~p(% lli)p(li) n pl) (3)
' V ' ' e^N(ej) k
where, N(e,) is the set of adjacent edges to e,.
In order to find the optimum labeling configuration for the whole graph with consistent labels for the vessel segments in relation with each other, we employ the simulated annealing stochastic relaxation method detailed in Geman and Geman (1984). Simulated annealing has previously been used to label cortical folding patterns in the adult brain (Mangin et al., 1995; Riviere et al., 2002). We define an energy function, U as an aggregation of edge label-dependent posterior probabilities:
U = -X log(p(l|))■ (4)
Given this definition for the energy function U on the whole graph, finding the optimum labeling configuration, arg maxLp(L\G) by MAP estimator can be translated into minimizing the energy function. The minimum energy state can be estimated using simulated annealing stochastic relaxation method (Gibbs sampler) (Geman and Geman, 1984). Adopting from the statistical mechanics, the probability distribution of
Fig. 5. The manual labeling result of 3 samples in the study. Top row is the dorsal view and the bottom row is the ventral view of the samples. All labeled vessels are given a unique label and color code.
the labeling configurations (states of the system) of the graph edges can be modeled by Gibbs distribution:
defined as
P(L) = 1 e-U(L)/T
where Z is the normalizing constant, also known as the partition function and is defined as an aggregation over the domain ofall possible labeling configurations:
Z = £ e
-U(()/T
T is the temperature parameter of the system which is controlled by an annealing factor n, and U(L) is the energy of the system at the given configuration L.
The stochastic relaxation algorithm for minimizing an energy function is detailed in Duda et al. (2001). In brief, the algorithm consists of first initializing the edge labels randomly, initializing the temperature of the system to a high temperature T0, and select the annealing parameter n. We chose annealing parameter of 0.99 by experiment. At each iteration, an edge e is picked randomly and the transition energy AU(le) = U(e/lenew)-U(e/lecurrent) for replacing its current label with all other possible labels is calculated. The actual label transition is done by drawing a label from the label transition probability density function
£ e-AU(le)=T
■. At every global iteration, when the energy tran-
sition was calculated for all the edges, the temperature is multiplied by n. The global iterations are repeated until the global energy of the system U(l) converges, meaning that the changes in the global energy of the system are smaller than a threshold in at least 10 consecutive global iterations.
The temperature parameter controls the probability of the edge label replacement e~AU(le)/T. At high temperatures, the changes of the system's state are more random and there is a higher probability to replace a label with one with a higher energy. This provides a wider search in the possible state space and prevents the system from becoming confined to a local optimum far from the global optimum by providing the possibility of jumping out of the local optimum through going to a higher energy state. As the temperature decreases, the transitions are forced towards the labels with a lower energy and the system converges to a low energy configuration (Geman and Geman, 1984).
Experiments and results
We have prepared 15 brain samples of female C57Bl/6J adult mice of age 16-20 weeks old by the Microfil CT contrast agent using the perfusion protocol detailed in Ghanavati et al. (2014). The Microfil contrast agent has minimal shrinkage and distortion effects (Chugh et al., 2009;
Fig. 6. The comparison of manual and automatic labeling and the automatic labeling error for a brain sample. The manual labeling (left), the result of the automatic labeling (middle), and the automatic labeling errors highlighted in red (right) are shown for each vasculature. Top row is the dorsal view and the bottom row is the ventral view of each sample.
Fig. 7. The mean recognition rate for each label propagated on one of the samples. The mean recognition rate of 50% is mapped in color red and mean recognition of 100% is mapped in color white. (a) The dorsal view, (b) ventral view, and (c) sagittal view of the sample.
Fig. 8. The correct recognition rate for the 7 samples in the 5 tests of LOO for 3 methods: MAP on local features (local MAP), gradient descent optimization (GD), and stochastic relaxation optimization (SR). (a) The correct recognition rate by the percentage of vessel segments correctly labeled compared to the ground truth manual labeling, and (b) The volume weighted recognition rate. The error bars show the standard deviation of the corresponding recognition rates.
Fig. 9. The distribution of the energy state of the system for the 7 test samples relative to the recognition rate in the three automatic labeling schemes: the local MAP (square), the gradient descent (point), and stochastic relaxation (triangle).
Cortell, 1969). The measurements performed by Rennie et al. showed the Microfil shrinkage to be less than 2.6% in diameter (Rennie et al., 2007). 4 of the samples were rejected due to a complication in the surgical steps during the perfusion process, according to the troubleshooting flowchart in Ghanavati et al. (2014) or damage to the surface vessels in the brain dissection. The 11 successfully perfused samples have been CT scanned (SkyScan 1172 high-resolution micro-CT, Bruker micro-CT, Belgium) with a10W X-ray tube at 50 kV and 200 |jA, and an isotropic resolution of 7 |am. The field of view covers the whole brain approximately 11 mm x 16 mm x 8 mm or 1500 x 2300 x 1100 voxels. The mean scan time of each sample was about 7 h and the total absorbed radiation dose of each specimen was approximately 20 Gy. Out of the successful perfusions, 4 have been rejected after reviewing the micro-CT images because of gaps in the filling of major vessels, especially in the large sinuses.
Segmentation and labeling
The vasculature in the 7 micro-CT images was segmented. The automatic segmentation of each image took on average about 100±18hon
a Linux x 86-64 Intel Core i7 machine. The mean number of edges in the segmented graphs was 2826 ± 415. The smallest diameter of vessels that were detecting was 20 |jm. The point spread function of the scanner, the signal to noise ratio, and the inherent limitations of the algorithm all contribute to determining the smallest detectable vessel (Marxen et al., 2004). We used the comparison of tracked vasculature with the binarized image as a means to assess the completeness of the tracked vessels by our algorithm. Our vessel tracking algorithm may fail at branching points and miss a whole tree, in which case, the comparison with the binarized image acquired by thresholding (although its accuracy is poorer than the vessel tracking algorithm) can highlight the missed trees. The comparison of the binarized original micro-CT image with the binarized image constructed from the segmented vascu-lature graph showed that the automatic vascular segmentation covered 91.5 ± 1.04% of the total volume of vessels in the micro-CT images in our study.
In total, 54 labels were used to label the cerebrovascular graph manually in the Brain-view software. The manual labeling of each sample took on average about 5 days of work. The labels were given a unique color code in order to facilitate the manual labeling procedure. The major vessels that are consistent across the population have established names. An example of some of these named vessels can be seen in Fig. 3. The smaller branches of the major named vessels were given a separate unique label and color code. In Fig. 5, the dorsal and ventral views of 3 of the manually labeled samples are shown. The color coding of the labels are the same as in Fig. 3.
Automatic labeling
Due to the small number of samples in the study, we have run LOO test on the 7 C57Bl/6J. We have tested alternate optimization procedures for the automatic labeling of the cerebrovasculature. The first was the MAP estimator on the local features. The second was to optimize the defined energy function using a gradient descent optimizer, and finally, the minimization of the energy function using the stochastic relaxation scheme as explained above. In order to examine the repro-ducibility of the results, each labeling scheme was repeated 5 times for each sample. The automatic labeling algorithm of each specimen takes on average 10.5 ± 2.5 h with a single threaded implementation on a Linux x 86-64 Intel Core i7 machine.
Fig. 6 shows a comparison of the ground truth labeling and the result of the automatic labeling of one of the data sets labeled in the LOO test. The vessel segments that were mislabeled in the automatic labeling
Table 1
The mean overall recognition rate and the mean recognition rate of some of the major arteries, veins and sinuses in the 7 samples in the LOO test.
Correct recognition rate (%)
Vessel numbers
Volumetric
Mean ± STD Max Min Mean ± STD Max Min
Overall 79.02 ± 2.99 84.54 73.14 78.06 ± 3.47 83.36 70.69
Internal carotid 93.69 ± 8.36 100 75 93.17 ± 11.03 100 69.93
Middle cerebral 80.21 ± 17.58 100 57.5 79.61 ± 13.42 100 55.1
Anterior cerebral 86.01 ± 12.09 100 66.67 92.24 ± 7.89 100 79.73
Azygos anterior cerebral 86.91 ± 11.86 94.74 60 86.45 ± 9.41 97.27 73.93
Olfactory 87.33 ± 5.25 93.1 74.68 79.62 ± 2.15 84 76.73
Superior cerebellar 90 ± 16.09 100 66.67 96.03 ± 10.36 100 66.67
Basilar 72.26 ± 14.17 100 57.14 83.96 ± 10.98 100 67.77
Vertebral 81.48 ± 24.75 100 51.11 92.87 ± 14.2 100 65.76
Thalamoperforating 92.31 ± 8.78 100 66.67 90.05 ± 7.39 100 69.95
Transverse sinus 78.77 ± 14.96 100 66.67 95.94 ± 4.69 100 87.18
Superior sagittal sinus 76.47 ± 25.72 100 50 99.74 ± 0.28 100 99.45
Great cerebral vein of Galen 89.59 ± 13.99 100 66.67 89.8 ± 15.31 100 68.41
Caudal rhinal vein 87.19 ± 9.4 96.03 64.32 88.88 ± 9.21 98.26 63.16
Rostral rhinal vein 76.69 ± 10.67 100 61.67 80.76 ± 12.75 100 61.9
Superior olfactory vein 98.81 ± 2.71 100 92.86 99.34 ± 1.49 100 96.06
Longitudinal hippocampal vein 85.38 ± 9.72 100 66.67 91.68 ± 8.01 100 60.91
Thalamostriate vein 93.99 ± 9.93 100 66.67 94.57 ± 10.25 100 71.01
were highlighted in red. As it can be seen, the major vessels that are more consistent across the population have a higher accuracy in the automatic labeling, as a result of the presence of similar samples in the training data set. In order to quantify the performance of the automatic labeling algorithm, the correct recognition rate was calculated by measuring the proportion of vessel segments correctly labeled with respect to the ground truth to the total number of vessels in the graph. The correct recognition rate was averaged over the 5 repetitions of the algorithm for each of the samples in the LOO test and is presented in Fig. 8. Volume weighted recognition rate was also calculated by the percentage of the vascular volumes that were correctly labeled with respect to the total cerebrovascular volume in the graph. The total volume of vessels is approximated by the sum of the volume of cylindrical objects representing the vessel tubular objects in the segmented graph. The volume weighted recognition rate is a better indicator of the automatic labeling accuracy over major vessels, since they occupy the larger portion of the total cerebrovascular volume. The segment number recognition rate is a better indicator of the automatic labeling accuracy over smaller branches, which are larger in number but are shorter and thinner vessels, and thus occupy a smaller volume. In order to present an overview of the automatic labeling accuracy by vessel label, we propagated the mean recognition rate of each vessel label onto the vasculature of one of the samples and mapped the mean recognition rate onto a red color map, which is shown in Fig. 7. As it can be seen, the major arteries and veins have a lighter shade of red compared to their smaller branches, indicating a higher recognition rate. The minimum recognition rate of labels over the 7 specimens is 50%.
The distribution of the energy state relative to the recognition rate for the 5 repetitions of LOO test on each of the 7 samples is shown in Fig. 9, for the MAP with local features, gradient descent and stochastic relaxation methods.
Despite the high variations among the vascular systems of different brain samples, there exist stereotyped main branches that are present in all samples (Larrivee et al., 2009). It is important to evaluate the
performance of the automatic labeling algorithm on these main vessels that are more consistent among the population. As a result, the mean recognition rate was estimated on the major arteries, veins and sinuses over the 7 test cases in 5 trials of the LOO test with the stochastic relaxation method and the result is shown in Table 1.
Discussion
We have perfused the cerebrovasculature of 7 C57Bl/6J mice with Microfil and imaged them with micro-CT. The vasculature was segmented and labeled automatically using our proposed method based on a stochastic relaxation technique. We tested our method by 5 trials of the LOO test. The results showed the correct recognition rate of the method to be > 75%.
It is important to note that for the mislabeled vessel segments with MAP estimator on the local features, the correct label is usually the second or third best candidate (with second or third highest posterior probability). On the other hand, despite the high inter-subject variabilities, the neighborhood structure of the vessels in the cerebral vasculature is fairly consistent among individuals (e.g. the internal carotid arteries (ICA) always bifurcates into MCA and ACA and the basilar artery always bifurcates into posterior cerebral artery (PCA)). As a result, the neighborhood information of the vessels can be a powerful feature to increase the accuracy of the labeling process. As seen from Fig. 9, the stochastic relaxation converges to a lower energy state and in most of the test cases to a higher recognition rate compared to the gradient descent minimization scheme. In our 7 test cases, except for one outlier, the labeling result of the stochastic relaxation technique was more accurate than the gradient descent method. However, the gradient descent minimization is much faster than the stochastic relaxation method since it converges very quickly to the nearest local minimum. In some cases the result from the gradient descent minimization method can be sufficient, especially if the method's run time is an import factor.
Fig. 10. The confusion matrix of the automatic labeling. The probability of mislabeling vessel segments of a given label with other labels is calculated and depicted. The entries on the main diagonal show the probability of labeling vessel segments correctly. The labels are numbered on the axes. The complete list of label names are provided in the Supplementary materials.
Table 2
Mean, standard deviation and coefficient of variation of some arterial features calculated in 7 labeled C57B1/6J specimens.
Mean ± SD Min Max CV(%)
ICA diameter (|jm) 230 ± 15 203 258 6.5
Mean ACA diameter (|m) 172 ± 34 134 214 19.7
Azygos arteries length (mm) 14.75 ± 3.7 11.44 19.56 25.1
MCA branches volume (|m3) 316 ± 69 209 440 21.8
ACA perfusion territory volume (cm3) 97.72 ± 17.85 60.34 113.77 18.3
MCA perfusion territory volume (cm3) 226.79 ± 32.63 166.57 279.25 14.4
SCA perfusion territory volume (cm3) 77.28 ± 19.66 52.11 103.65 25.4
In order to examine the mislabeled vessel segments, we have constructed the confusion matrix of the automatic labeling by calculating the probability of labeling the vessel segments of a given label with each of the possible labels over the 7 specimens (Fig. 10). The main diagonal of the matrix shows the probability of labeling vessel segments with their ground truth labels. The other entries of the confusion matrix show the probability of mislabeling vessel segments of a given label. As we expected, the probabilities of the main diagonal are the highest for each label. Most of the off-diagonal non-zero probabilities in the confusion matrix are related to labeling errors such as mislabeling the daughter branches of ACA to daughter branches of MCA, the daughter branches of PCA to PCA, the daughter branches of SCA to SCA, and the longitudinal hippocampal vein to great cerebral vein of Galen. As seen from the confusion matrix and Fig. 7, most of the errors in the automatic labeling method happen in the small branches of a label that are adjacent or near the small branches of another label. These branches often have similar length, diameter and approximate position, and as a result, the posterior probability of the 2 labels is close. This can lead to the mislabeling of a terminal branch with the label of a neighbor terminal
branch. The main labeling errors in large vessels happen either at the bifurcation of a label to two other different labels, or at the branching point of a major vessel and its daughter branches, which have been assigned a separate label. In the latter case, the manual selection of the exact point where a major vessel with many branches, such as MCA and ACA, is assigned a new label is to some extent arbitrary and not always consistent among all cases. We have tried to have the manual labeling of the test cases consistent as much as possible, however, some variations are inevitable given the high variability of the branching points and daughter branches of major vessels such as MCA.
The automatic labeling of over 3000 vessel segments in the mouse cerebrovasculature is a complex machine learning problem and requires a large training data set and descriptive features. In this study, we have developed a database of labeled cerebrovasculature samples to be used as the training data set for the automatic labeling algorithm and have achieved a correct recognition rate of more than 75%. Correcting the mislabeled vessel segments on the outcome of the automatic labeling is much faster than manual labeling of a complete cerebrovasculature. The labeled cerebrovasculature database and automatic labeling algorithm developed and presented in this paper will be used to expand the training data set. We expect that the expansion of the training data set improves the accuracy and generalization of the automatic labeling algorithm proposed here, similar to the work on cortical sulci recognition, where the expansion of the database that could represent the high inter-subjects variability of sulcal patterns improved the accuracy and generalization of the method (Mangin et al., 2004; Perrot et al., 2011; Riviere et al., 2002). In addition, the training data set needs to be expanded to include other mouse strains with significantly different cerebrovasculature. The feasibility of using a strain mixed training data set versus single strain training set for the automatic labeling algorithm should be examined as well.
Fig. 11. Average perfusion territory map by arterial labels (left), and the accuracy map of the average territories (right) of 7 C57Bl/6J specimens; in coronal, axial and sagittal planes. The colors of the perfusion territories match the color codes of their perfusing arterial labels. The perfusion territories marked on the figure are perfused by the terminal branches of left and right middle cerebral arteries (MCA), olfactory artery, cortical surface anterior cerebral artery (ACA), azygos anterior cerebral artery (Azygos), thalamoperforating arteries, left and right superior cerebellar arteries (SCA), and anterior inferior cerebellar arteries (AICA).
Given the labeled segmented vasculatures, it is possible to compare and quantify the features of different vessels across the specimens. Moreover, arterial and venous vessels can be analyzed separately based on the labeled data. Examples of features that one can compute for specific arteries are mean vessel diameter, total sub-tree vessel length, and total intravascular volume. An illustrative selection of these parameters is reported in Table 2, along with standard deviation, minimum and maximum for the 7 labeled C57Bl/6J specimens. The mean diameter of arteries calculated by our method is comparable to published values measured by in vivo angiography methods (Kidoguchi et al., 2006), indicating the minimal shrinkage and distortion effect of Microfil perfusion and sample preparation. The coefficient of variation for these different metrics varied between 6.5% and 25.4%. The analysis of the labeled vasculature can also provide qualitative information such as the existence of collateral circulation in the circle of Willis. For example, three of the C57Bl/6J specimens had no posterior communicating artery. The other four specimens had a unilateral posterior communicating artery with a mean diameter of 125 ± 11 |jm. For each specimen, a coarse perfusion territory map was constructed by assigning the label of the nearest arterial terminal edge to each voxel. An average perfusion territory map and the accuracy map of the average territories were created by voxel voting over the 7 C57Bl/6J specimens (Fig. 11).
Conclusion
Comparison of cerebrovascular variations can broaden our understanding of this anatomical structure. The vascular morphology is highly variable from one individual to another. The variations include different number of branches, variation in the location of bifurcation, absence of certain vessel segments, bifurcation or trifurcation, etc. The high complexity and variability of the vascular structure make it very difficult to compare vessels among individuals which are essential in studies of the vasculature development and pathologies. The automatic recognition of the vessels through registration to an atlas would also fail due to the large inter-subject variations. The variation ofvascular structures, such as the location and angle of branches, missing or additional branches, and smoothness and diameter of vessel segments from one subject to another make an accurate co-registration of the vasculature almost impossible. We have developed a framework to segment and label all the vessels in the cerebral vasculature automatically. We have shown briefly that this automatic method facilitates and expedites the comparison and analysis of many cerebrovascular data in a reasonable time. In the future we will utilize this method to compare anatomical variations of cerebro-vasculature caused by different mouse strains and to study the geno-type-phenotype relationships in the development of cerebral vessels. In addition, the processing methods introduced in this paper can be adapted with minimal changes to in vivo images for applications such as longitudinal studies. Currently, the typical in vivo angiography resolution is limited between 50 and 100 |jm, which is not high enough to visualize the fine arterioles and venules. We anticipate that as in vivo angiography methods continue to improve this will be a good application for the presented methodology.
Acknowledgments
This research was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) 418539-2012 and the Canadian Institutes of Health Research (CIHR) 231389.
Appendix A. Supplementary data
Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.neuroimage.2014.03.044.
References
Babin, D., Pizurica, A., Vylder, J.D., Vansteenkiste, E., Philips, W., 2013. Brain blood vessel segmentation using line-shaped profiles. Phys. Med. Biol. 58,8041.
Beckmann, N., 2000. High resolution magnetic resonance angiography non-invasively reveals mouse strain differences in the cerebrovascular anatomy in vivo. Magn. Reson. Med. 44, 252-258.
Bell, R.D., Zlokovic, B.V., 2009. Neurovascular mechanisms and blood-brain barrier disorder in Alzheimer's disease. Acta Neuropathol. 118,103-113.
Bilgel, M., Roy, S., Carass, A., Nyquist, P.A., Prince, J.L., 2013. Automated Anatomical Labeling of the Cerebral Arteries Using Belief Propagation. SPIE Medical Imaging, International Society for Optics and Photonics (866918-866918).
Bogunovic, H., Pozo,J., Cardenes, R., San Roman, L., Frangi, A., 2013. Anatomical labeling of the circle of Willis using maximum a posteriori probability estimation. IEEE Trans. Med. Imaging 32,1587-1599.
Caselles, V., Kimmel, R., Sapiro, G., 1997. Geodesic active contours. Int. J. Comput. Vis. 22, 61-79.
Chugh, B.P., Lerch, J.P., Yu, L.X., Pienkowski, M., Harrison, R.V., Henkelman, R.M., Sled, J.G., 2009. Measurement of cerebral blood volume in mouse brain regions using micro-computed tomography. NeuroImage 47,1312-1318.
Cortell, S., 1969. Silicone rubber for renal tubular injection. J. Appl. Physiol. 26,158-159.
de la Torre, J.C., 2004. Is Alzheimer's disease a neurodegenerative or a vascular disorder? Data, dogma, and dialectics. Lancet Neurol. 3,184-190.
de la Torre, J., Stefano, G., 2000. Evidence that Alzheimer's disease is a microvascular disorder: the role of constitutive nitric oxide. Brain Res. Rev. 34,119-136.
Di Ruberto, C., Rodriguez, G., Casta, L., 2002. Recognition of Shapes by Morphological Attributed Relational Graphs. Proceedings of Workshop sulla Percezione e Visione delle Macchine, Siena, Italy.
Dorr, A.E., Sled, J.G., Kabani, N., 2007. Three-dimensional cerebral vasculature of the CBA mouse brain: a magnetic resonance imaging and micro computed tomography study. NeuroImage 35,1409-1423.
Dorr, A.E., Lerch, J.P., Spring, S., Kabani, N., Henkelman, R.M., 2008. High resolution three-dimensional brain atlas using an average magnetic resonance image of 40 adult C57Bl/6J mice. NeuroImage 42, 60-69.
Duda, R.O., Hart, P.E., Stork, D.G. (Eds.), 2001. Pattern Classification, 2nd ed. John Wiley & Sons Inc., Hoboken, NJ., USA.
Forkert, N.D., Schmidt-Richberg, A., Fiehler, J., Illies, T., Moller, D., Saring, D., Handels, H., Ehrhardt, J., 2013.3D cerebrovascular segmentation combining fuzzy vessel enhancement and level-sets with anisotropic energy weights. Magn. Reson. Imaging 31, 262-271.
Frangi, A.F., Niessen, W.J., Vincken, K.L., Viergever, M.A., 1998. Multiscale Vessel Enhancement Filtering. Medical Image Computing and Computer-Assisted Interventation. Springer, pp. 130-137.
Fridman, Y., Pizer, S.M., Aylward, S., Bullitt, E., 2004. Extracting branching tubular object geometry via cores. Med. Image Anal. 8,169-176.
Geman, S., Geman, D., 1984. Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE Trans. Pattern Anal. Mach. Intell. 6, 721-741.
Ghanavati, S., Yu, L.X., Lerch, J.P., Sled, J.G., 2014. A perfusion procedure for imaging of the mouse cerebral vasculature by X-ray micro-CT. J. Neurosci. Methods 221, 70-77.
Kent, J., 1982. The Fisher-Bingham distribution on the sphere. J. R. Stat. Soc. 44, 71-80.
Kidoguchi, K., Tamaki, M., Mizobe, T., Koyama, J., Kondoh, T., Kohmura, E., Sakurai, T., Yokono, K., Umetani, K., 2006. In vivo X-ray angiography in the mouse brain using synchrotron radiation. Stroke 37,1856-1861.
Kitagawa, K., Matsumoto, M., Yang, G., Mabuchi, T., Yagita, Y., Hori, M., Yanagihara, T., 1998. Cerebral ischemia after bilateral carotid artery occlusion and intraluminal suture occlusion in mice: evaluation of the patency of the posterior communicating artery. J. Cereb. Blood Flow Metab. 18,570-579.
Krissian, K., Malandain, G., Ayache, N., Vaillant, R., Trousset, Y., 2000. Model-based detection of tubular structures in 3D images. Comput Vis. Image Underst 80,130-171.
Larrivee, B., Freitas, C., Suchting, S., Brunei, I., Eichmann, A., 2009. Guidance of vascular development. Circ. Res. 104,428-441.
Lesage, D., Angelini, E.D., Bloch, I., Funka-Lea, G., 2009. A review of 3D vessel lumen segmentation techniques: models, features and extraction schemes. Med. Image Anal. 13, 819-845.
Lorigo, L., Faugeras, O., Grimson, W., Keriven, R., Kikinis, R., Nabavi, A., Westin, C.F., 2001. Curves: curve evolution for vessel segmentation. Med. Image Anal. 5,195-206.
Mangin, J.F., Regis, J., Bloch, I., Frouin, V., Samson, Y., Lopez-Krahe, J., 1995. A Markovian Random Field Based Random Graph Modelling the Human Cortical Topography. First Int Conf on Computer Vision, Virtual Reality and Robotics in Medicine. Springer, New York, pp. 177-183.
Mangin, J.F., Riviere, D., Cachia, A., Duchesnay, E., Cointepas, Y., Papadopoulos-Orfanos, D., Scifo, P., Ochiai, T., Brunelle, F., Regis, J., 2004. A framework to study the cortical folding patterns. NeuroImage 23 (Suppl. 1), S129-S138.
Mardia, K.V., 1972. Statistics of Directional Data. Academic Press, San Diego, CA., USA.
Mardia, K.V., Jupp, P.E., 2009. Directional Statistics. John Wiley & Sons Ltd., Hoboken, NJ., USA.
Marxen, M., Thornton, M.M., Chiarot, C.B., Klement, G., Koprivnikar, J., Sled, J.G., Henkelman, R.M., 2004. MicroCT scanner performance and considerations for vascular specimen imaging. Med. Phys. 31, 305.
Mori, K., Hasegawa, J., Suenaga, Y., Toriwaki, J., 2000. Automated anatomical labeling of the bronchial branch and its application to the virtual bronchoscopy system. IEEE Trans. Med. Imaging 103-114.
Mori, K., Ota, S., Deguchi, D., Kitasaka, T., Suenaga, Y., Iwano, S., Hasegawa, Y., Takabatake, H., Mori, M., Natori, H., 2009. Automated Anatomical Labeling of Bronchial Branches Extracted From CT Datasets Based on Machine Learning and Combination Optimization and its Application to Bronchoscope Guidance. Proceedings of the 12th
International Conference on Medical Image Computing and Computer-Assisted Intervention: Part II. Springer-Verlag, Berlin, Heidelberg, pp. 707-714.
Mori, K., Oda, M., Egusa, T., Jiang, Z., Kitasaka, T., Fujiwara, M., Misawa, K., 2010. Automated Nomenclature of Upper Abdominal Arteries for Displaying Anatomical Names on Virtual Laparoscopic Images. Medical Imaging and Augmented Reality, pp.353-362.
Perrot, M., Riviere, D., Mangin, J.F., 2011. Cortical sulci recognition and spatial normalization. Med. Image Anal. 15, 529-550.
Piccinelli, M., Veneziani, A., Steinman, DA., Remuzzi, A., Antiga, L., 2009. A framework for geometric analysis of vascular structures: application to cerebral aneurysms. IEEE Trans. Med. Imaging 28,1141-1155.
Qian, X., Brennan, M.P., Dione, D.P., Dobrucki, W.L., Jackowski, M.P., Breuer, C.K., Sinusas, A.J., Papademetris, X., 2009. A non-parametric vessel detection method for complex vascular structures. Med. Image Anal. 13,49-61.
Rennie, M.Y., Whiteley, K.J., Kulandavelu, S., Adamson, S.L., Sled, J.G., 2007. 3D visualisation and quantification by microcomputed tomography of late gestational changes in the arterial and venous feto-placental vasculature of the mouse. Placenta 28, 833-840.
Rennie, M.Y., Detmar, J., Whiteley, K.J., Yang, J., Jurisicova, A., Adamson, S.L., Sled, J.G., 2011. Vessel tortuosity and reduced vascularization in the fetoplacental arterial tree after maternal exposure to polycyclic aromatic hydrocarbons. Am. J. Physiol. Heart Circ. Physiol. 300, H675-H684.
Riviere, D., Mangin, J.F., Papadopoulos-Orfanos, D., Martinez, J.M., Frouin, V., Regis, J., 2002. Automatic recognition of cortical sulci of the human brain using a congregation of neural networks. Med. Image Anal. 6, 77-92.
Scremin, O.U., Holschneider, D.P., 2011. The Mouse Nervous System. In: Watson, C., Paxinos, G., Puelles, L. (Eds.), Vascular Supply. Academic Press, San Diego, USA.
Shang, Y., Deklerck R., Nyssen, E., Markova, A., de Mey, J., Yang, X., Sun, K., 2011. Vascular active contour for vessel tree segmentation. IEEE Trans. Biomed. Eng. 58,1023-1032.
Sofka, M., Stewart, C., 2006. Retinal vessel centerline extraction using multiscale matched filters, confidence and edge measures. IEEE Trans. Med. Imaging 25,1531-1546.
Suri, J.S., Liu, K., Reden, L., Laxminarayan, S., 2002. A review on MR vascular image processing: skeleton versus nonskeleton approaches: part II. IEEE Trans. Inf. Technol. Biomed. 6, 338-350.
Tschirren, J., Mclennan, G., Palagyi, K., Hoffman, E.A., Sonka, M., 2005. Matching and anatomical labeling of human airway tree. IEEE Trans. Med. Imaging 24,1540-1547.
