## Abstract

An automatic segmentation technique has been developed and applied to two renal micro-computer tomography (CT) images. With the use of a 20-μm voxel resolution image, the arterial and venous trees were segmented for the rat renal vasculature, distinguishing resolving vessels down to 30 μm in radius. A higher resolution 4-μm voxel image of a renal vascular subtree, with vessel radial values down to 10 μm, was segmented. Strahler ordering was applied to each subtree using an iterative scheme developed to integrate information from the two segmented models to reconstruct the complete topology of the entire vascular tree. An error analysis of the assigned orders quantified the robustness of the ordering process for the full model. Radial, length, and connectivity data of the complete arterial and venous trees are reported by order. Substantial parallelism is observed between individual arteries and veins, and the ratio of parallel vessel radii is quantified via a power law. A strong correlation with Murray's Law was established, providing convincing evidence of the “minimum work” hypothesis. Results were compared with theoretical branch angle formulations, based on the principles of “minimum shear force,” were inconclusive. Three-dimensional reconstructions of renal vascular trees collected are made freely available^{1} for further investigation into renal physiology and modeling studies.

- kidney
- Strahler ordering
- vascular statistics
- renal modeling
- vascular reconstruction

renal vascular structure has been the focus of a number of studies analyzing the affects of conditions such as hypertension and diabetes (12, 17, 23, 24, 28, 29, 38) as well as studies attempting to characterize renal function mathematically (10, 26, 27). The kidney has been shown to play a pivotal role in blood volume, vascular tone, electrolyte control, as well as blood filtration. Strong structural heterogeneity and correlates to functional physiology (3, 18) make the study of renal vasculature essential to understanding renal function and consequently, blood toxicity, volume, and pressure.

The current anatomical ontology characterizes renal vasculature based on its functional role in the circulation. Blood flow in the rat kidney enters through the renal artery splitting into two arterial trees that enter via the pelvic region dividing into 6–10 (8 in the data set we segment below) interlobar arteries (IA). These vessels lead along the corticomedullary boundary (3) eventually feeding into the arcuate arteries (AA). Derivatives of these arteries are the cortical radial arteries (CRA), which penetrate through the cortical labyrinth (34). Renal arterial anatomy is described in detail in Reference 3 and the microvasculature in Reference 18. For standard nomenclature, refer to Kriz and Bankir (25).

Microcomputed tomography (micro-CT) provides the ability to produce high resolution detailed images from which information on the vascular structures can be collected (4, 5, 6). The development of automatic segmentation methods to extract relevant information from these data sets and the use of organization techniques provide a novel set of tools for analyzing and quantifying the structure of vascular trees.

Significant work has gone into measuring and enumerating the topology of vascular trees (7, 37, 43, 45, 46), particularly recent work in the heart (20–22). Because of the amount and the form of data in these trees, organizational techniques are essential to quantifying structural anatomy and gaining insight into physiological processes. A number of ordering methods have been proposed to tabulate and organize vascular data (19, 44); however, the most commonly applied is the so-called Strahler Ordering method proposed by Strahler (40) and Horsfield (15, 16). Strahler ordering begins by labeling all terminal arterioles as *order 0*. Following vessel segments upstream, if two *order 0* vessels adjoin at a bifurcation, the order of the parent vessel is labeled 1, whereas if an *order 0* vessel merges with a vessel of dissimilar order, the parent assumes the value of the highest ordered vessel (40). This scheme, based entirely on the vascular structure and connectivity, has been shown to correlate well with vessel radius (7, 19, 37, 45, 46). As blood flow through the systemic circulation is, structurally, influenced most significantly by vessel diameters, Strahler ordering provides insight into blood distribution and transport.

The focus in this study was the vasculature in the unilobar cortex of the rat kidney. Renal rat models are commonly used for developing a scientific understanding of renal function (2–4, 12, 17, 23, 24, 28, 29, 38). Raw data from micro-CT images of vascular casts (11) were segmented by using a method outlined below. The extracted three-dimensional arterial and venous trees were then analyzed for an entire renal cortex at a voxel resolution of 20 μm and for a cortical subtree at a voxel resolution of 4 μm. Despite working with two partial images, we developed an iterative technique to correlate ordering and statistical data from both to assess the complete vascular structure. Finally, a quantitative error analysis technique is proposed to assess the validity of the applied methods. The three-dimensional data sets of these segmented vascular trees are made available to both the biological modeling and physiological communities^{1} for further investigation.

## METHODS

#### Images and segmentation.

Garcia-Sanz et al. (11), using micro-CT imaging, acquired scanned voxel images of rat renal vasculature. Renal circulation was infused at physiological pressures using the protocol described by Bentley et al. (4). Briefly, the sample was perfused using a radiopaque silicone polymer, effectively delineating the vasculature from the organ tissue via a uniform coloration specific to the vascular cast. The resulting voxel grayscale images were broken into two-dimensional images using the National Institutes of Health sponsored program ImageJ.^{2} Renal images of 20-μm (showing an entire kidney) and 4-μm (showing a vascular subtree) voxel resolution were used for this analysis (4).

Further image processing was conducted in CMGUI, a three-dimensional graphical environment developed at the Bioengineering Institute.^{3} The extracted images were uploaded, and isosurface rendering was used to locate vessel walls. Isosurfaces were stored as a point group, consisting of up to a million points, dispersed on arterial and venous walls. With the use of this point group, the vascular trees were skeletonized while approximate vessel radii were calculated, a process referred to in this paper as vessel emaciation, resulting in vessel diameters, lengths, and connectivity calculations.

With the use of the point group generated from the kidney images, emaciation of vessels was conducted based on the assumption that individual vessel segments form elliptical cross sections, and on average the length between two nearby vessels is greater than that of their respective diameters. With these assumptions, peripheral points on vessel walls were projected inward by surface projections, emaciating all vessels. This process is demonstrated in Fig. 1 on the 20-μm kidney image. Figure 2 provides a general schematic of the segmentation process where the detailed vector algebra for this process is outlined in detail in appendix a. The final result is a stream of points at the centerline of all the vessels as seen in Fig. 1*F*. Tracking, based on a snake algorithm adapted to handle the three-dimensional point group, was utilized to connect points and successfully recreate the original vasculature. With the use of the original point group to define the walls of the vasculature, the fully emaciated set of points were progressively linked and centered within the vessel cross section. This process was repeated until all points had been tracked and associated to a vessel segment, resulting in a list of spatial positions, radii, and local connectivity indicative of the original vascular image.

#### Organizational technique.

The scheme adopted for this work was Strahler ordering (outlined in the introduction). Because of the branching morphology of renal vasculature, the Strahler ordering technique required additional rules to handle higher-order junctions such as trifurcations (which occurred between 3% and 5.5% for all segmented branch points in this study across the different network subtrees) that were observed more often than reported in other ordering studies (20). These rules are illustrated in Fig. 3 (dagger denotes location of trifurcations). The convention adopted was to apply Strahler ordering to those branching vessels with multiple daughters by increasing the parent vessel segment an order if two or more vessels of the same and highest order merge (†_{2}), otherwise labeling the parent vessel as that of the highest order daughter vessel (†_{1}).

The Strahler ordering scheme works logically from capillaries or venules upstream and is, consequently, sensitive to pruned subtrees. Each image in this study was incomplete due to resolution constraints and partial casting of vessels. To accommodate this, an iterative technique was adopted to define the proper Strahler order of each vessel terminal. Under this scheme a conservative estimate of true *order 0* terminal vessels was made while all other presumably pruned terminal vessels were left undefined.^{4} Ordering proceeded through the subtrees until undeclared terminal orders inhibited further progress. Based on the statistics generated from this initial ordering, undeclared terminal vessels were progressively defined according to their terminal bounds until the ordering was complete.^{5} This process was iterated by reapplying statistical ranges of the initial ordering to all terminals and reordering until convergence was obtained (i.e., no statistical deviation from one iteration to the next). Statistical ranges from the high-resolution image (4 μm) were then used as an estimate for the bounds of terminal vessels segmented in the low-resolution image (20 μm).

#### Ordering error analysis.

Central to effective ordering of any vascular tree is the proper assignment of vessel terminals. However, the assumption of a terminal vessel ordering profile, as the previous section suggests, imparts error associated with commuting subtrees with the entire vascular tree. A method was designed to assess how these errors affect results of the Strahler classification. For example, Fig. 3 illustrates common artifacts seen in medical images of vasculature (denoted with ‡). The gray marker infers a bifurcation where incorrect ordering has no effect on the upstream ordering, as changing the first-order vessel to a terminal zero-order vessel results in the same upstream ordering. Consequently, the probability that vessels were ordered incorrectly due to this terminal is zero. The dark markers, however, infer bifurcations where not applying the iterative technique entirely changes the ordering upstream as well as the calculated statistics (i.e., changing these terminal vessels to zero order causes a restructuring of the upstream ordering).

To assess these errors one might change terminal orders of the vascular network and compute the effect. This process, however, would be computationally intensive, especially for large-scale vascular trees. The framework developed provides a means for quantifying the sensitivity of a vascular tree to likely perturbations of terminal vessels, assuming the given ordering scheme is a close match to the correct ordering.

We define *ρ*_{i} = [*ρ**ρ**ρ*]^{T} to be a vector of probable outcomes for a vessel segment *S*_{i}, where the components are the probabilities of that segment increasing (*ρ*), remaining (*ρ*), or decreasing (*ρ*) an order. The transference of probability through vessel segments (and consequently the overall effect) is mediated by the series of bifurcations, or operators, and Strahler ordering rules at each. These operators summarize all possible outcomes given *ρ* for the daughter segments and modify the upstream segment probabilities accordingly. Operators, denoted as *ζ*, have subscripts that label a specific branch point and the superscript *T* referring to the matrix transpose (for complete details on operators, refer to appendix b).

When calculating the probability profile of the *j*th terminal vessel segment *S*_{j}, with radius *r*_{j} that is of order *i,* we define the upper (r_{max}) bound by the same diameter defined criterion proposed by Jiang et al.^{6} (19) shown in *Eq. 1*, where *r*_{i} and *σ*_{i} are the average radius and standard deviation of order *i* vessels. (1) Note that the lower bound for order *i*+1 is set as the upper bound for order *i*. Using a normal distribution to describe error in ordering terminal vessels (with mean *r*_{j} and standard deviation *α*_{r}*r*_{j}, where *α*_{r} is a constant), we have assumed that terminal probabilities can be assigned by the area under the curve delineated by the bounds in *Eq. 1* as seen in Fig. 4.

To calculate *ρ*_{i} or the probability distribution of the *i*th network segment (*S*_{i}), the method requires the definition of a pathway *t*_{i}. Considering the subtree beneath *S*_{i}, *t*_{i} is an index of segments that defines a single path down the *i*th subtree to a terminal vessel. For example, if one such path from *S*_{i} to a terminal was given to go through segments *S*_{1}, *S*_{2,} and *S*_{3}, *t*_{i} = [*S*_{1} *S*_{2} *S*_{3}]. Choice for defining *t*_{i} is arbitrary but requires only movement downstream. A different path for the same segment (*S*_{i}) results in a modification of operators (which are defined by a chosen branch or branch pair) but no modification in the probability *ρ*_{i}. If *ρ* is the probability of the terminal vessel by *t*_{i}, the overall probability profile of *S*_{i} can be written as seen in *Eq. 2*. (2) The result is a directed product of operators acting on the terminal values of *ρ*. A sample calculation using the vascular tree from Fig. 3 is presented in appendix a. This calculation is a concise and efficient means for predicting how incorrectly ordering pruned terminals will alter ordering and, consequently, the statistics of the tree classified by order.

#### Calculating lost networks.

Kassab et al. (20, 21) presented the idea of the connectivity matrix (CM). This matrix, *C*, is order dependent with component *C*_{mn} representing the number of elements of order *m*, which arise from an order *n* vessel divided by the total number of order *n* vessels. The CM gives insight into the connectivity of a vascular tree, but perhaps more importantly it provides a means to approximate the number of vessels in a subtree by order when the appropriate feeding vessel is known. This has implications on blood flow and transport to downstream tissues. Unfortunately, in our study the CM is based on incomplete images and is consequently not accurate. However, by statistically defining terminal vessels, many pruned vascular trees can be identified. Let *n*_{v,O(n)} refer to the number of vessels of order *n* and *n*_{v,O(n)←O(m)} refer to the number of vessels of order *m* arising from a vessel of order *n*. If the prime superscript uses the entire network as a sample, and the *P* superscript is used to identify the subnetworks that consist of a pruned vessel (i.e., a terminal vessel with an order greater than the lowest vascular tree order) and all subtrees that branch from the pruned vessel at nonterminal positions. The CM is defined in *Eq. 3*. (3) It is important to note that the lowest order vessels of a network are not tabulated as part of the pruned vessels (as for within the image of interest, these vessels will not bifurcate). This formulation allows complete subtrees that arise from a pruned vessel to contribute to the global statistics of the network while eliminating the potential of underestimating the subtrees arising from that pruned vessel.

A complication of ordering terminals based on radial profile is that specific sets may lack low or high order terminal vessels. In this case, the 20-μm image lacked *order 0* and *1* vessels, whereas the 4-μm image lacked orders over 5 (arterial) and 7 (venous) under the Strahler scheme. To properly describe the CM, a weighted average (shown in *Eq. 4*) was used to amalgamate the results of both image sets based on the estimated number of total vessels in each sample. In *Eq. 4*, the weight *w*= *n*− *n*+ *C**n*, and the superscripts refer to either of the two image resolutions analyzed. (4) It should be noted that this formulation was used when the values from both data sets for *C*_{mn} were not zero.

The number of vessels in each given order of a subtree defined as the network downstream from an inlet (root vessel) can be determined by using the connectivity matrix. This is done by calculating the number of daughter vessels branching form the inlet and then iteratively progressing down the subtree. This process can be reduced to solving the set of simultaneous equations given below in *Eq. 5*. (5) Where *C* is the connectivity matrix calculated in *Eq. 4*, *I* is the identity matrix, and *A* is the *r*th row of *C*, where *r* is the order of the root vessel. Solving *Eq. 5* for *N* provides as its components *N _{i}*, the number of vessels of a given order

*i*.

## RESULTS

From the segmentation methods, both a full-scale rat kidney at 20 μm voxel resolution (Fig. 5) and a portion of kidney at 4 μm voxel resolution (Fig. 6) were reconstructed. The combination of these two image reconstructions were used in the ensuing statistical studies.

Each set was ordered using the Strahler method outlined. The average luminal radius of each order as well as the approximate number of vessels is given for the venous tree (Table 1) and arterial tree (Table 2). Because of vessel pruning and significant vessel tapering in the kidney, giving equal weight to complete vessels and proximal ends of pruned vessels can result in incorrect calculations. To account for this, the calculation of means and standard deviations was done by using a weighted average based on vessel length. These data are also given graphically for the venous tree in Fig. 7.

General vessel connectivity for the arterial tree is listed in Table 3, where the low-level vessel branching was included from the high resolution set as prescribed in the methods. The same information is given for the venous tree (Table 4).

Vessel length data are given in Table 5 for all trees organized according to Strahler ordering. Though vessel length often has a poor correlation with order, this information can be useful in calculating global vascular properties such as subtree volume.

Because of the parallel architecture of the arterial and venous trees, an empirical relation was derived to correlate the two. Figure 8 shows a relation between the parallel venous and arterial vascular trees. A power law relation is fitted to these data showing the ratio of arterial and venous radii for varying venous radii.

Error analysis of the ordering scheme is given in Fig. 9 for each tree. Assignment of terminal probabilities are given in detail in methods and varied according to radial factor (*α*_{r}). The results are shown graphically in Fig. 10.

Shown in Fig. 11 is the total downstream subtree volume at a specific bifurcation by the radius of the feeding vessel for the venous tree. With the use of the CM (Table 4) to estimate number of subtree vessels and the vessels statistics (Table 5), an expected subtree volume was plotted by average radius per order.

Cross-sectional area by order is plotted in Fig. 12 for the arterial and venous trees. The cross-sectional area was calculated as using information about the average radius and number of vessels from Tables 1 and 2.

A benefit over traditional cast studies is the quantification of spatial aspects of vasculature, such as branching angles that we have calculated from the two three-dimensional vectors, which define the branch angle geometry at each bifurcation. This allows for spatial heterogeneity, which has a physiological basis, but also allows for the experimental testing of hypotheses on vascular morphology. One such theory by Zamir (47) and Zamir and Bigelow (48) proposed that arterial branching angles are designed to minimize total shear force based on the ratio of daughter vessel radii (either linear or quadratic). Figure 13 shows experimental data compared with these the proposed optimums.

Another optimization theory is that proposed by Murray (31) where the fundamental structure of vascular trees is such that it minimizes work. Under ideal conditions where flow rate is conserved, Murray's Law states that the sum of each daughter vessel (*r*_{d}) cubed is equivalent to the parent vessel (*r*_{p}) cubed, or *r*= ∑*r*. Figure 14 shows the results for the arterial tree compared with those predicted by Murray's Law.

## DISCUSSION

The applied mathematical reduction of micro-CT images is a useful and relatively straightforward means for converting an image into a quantified reconstruction of vascular structure. Though data were framed and handled discretely by using point set representations, the same technique can be applied to the voxel image itself, providing a major boost in computational performance. The limitation of this process is that it requires sufficiently clear imaging of the vascular walls and boundaries. This proposed method is not intended as an image-processing technique (because it does nothing to account or accommodate for unclear imaging) but as a tool for automated segmentation.

A number of studies have looked at the effects of structure on function, primarily in the analysis of disease states. Recent studies have assessed vascular changes in kidneys following chronic bile duct ligation (33), hypercholesterolemia (5, 6), chronic nitric oxide inhibition (9), renal artery stenosis (50), and chronic potassium depletion (8). A significant setback to research in the morphological affects of physiological conditions is the lack of a concrete basis for comparison. Studies have reported vascular diameters that vary significantly (24, 38, 44, 11) for the same vessel structure. To integrate information from these sorts of studies, it is necessary to look at specific regions of measurement to provide a context for evaluation of total renal physiological function.

Tables 1 and 2 give a quantitative correlation between radius and order. The class of afferent arterioles is represented as *order 0* to *1* vessels. Published radial values for afferents have spanned from 7 μm (24) to 14 μm (38) with a number of studies reporting in that range (17, 28, 36). This is consistent with vessels measured in the 4-μm data set. Cortical radial arteries span *orders 2* to *6*. These vessels bifurcate resulting in a significant span of orders as well as change in vessel radii. This is reflected in the literature, where deviations in radius range from 46.4 μm (42) to 68.2 μm (41). Because this particular class of vessels is so broad, the terminology only provides a loose way to correlate quantitative information. *Orders 6* to *7* comprise the arcuate arteries running along the corticomedullary boarder. As is the case with the cortical radial arteries, the arcuate arteries experience repeated bifurcations and radial reduction. Because the kidney experiences no anastomoses between arterial branches, arcuate arteries eventually reduce to cortical radial arteries. This could be a possible explanation of deviation seen in the literature where radial averages have been given as low as 96.5 μm (23) and as high as 153.5 μm (11). The interlobar arteries comprise *orders 8* to *9*, whereas the feeding renal arteries are of *orders 9* to *10*. Few other works have focused on the venous tree within the kidney (interlobular, *orders 2–6*; arcuate, *orders 6–7*; interlobar, *order 8–9*).

The estimated numbers of vessels by order are given in Tables 1 and 2. Though the actual number of vessels in the kidney have yet been tabulated, the number of glomeruli in a normotensive rat kidney has been estimated at ∼30,000 (2, 29, 32). With the assumption of an approximate 1:1 ratio between the afferents and glomeruli, the predicted number of glomeruli would be 29,566 ± 5,965 according to our study. However, because of resolution, incomplete perfusion, and the visualization methods used, the majority of tracked afferents (4-μm voxel image) and tracked cortical radial arteries (20-μm voxel image) were lost and not linked to the arterial tree. Though accounting for pruning in the connectivity matrix significantly improved results, to obtain a more accurate representation the arterial tree was correlated tightly to the venous tree.

A striking feature of renal vasculature is its parallel nature. In all images, arteries and veins could easily be distinguished by their path duplicity and disparity in vessel radii. Parallelism was maintained through the corticomedullary boundary between arcuate arteries and veins. Orders are often seen running in parallel (though the arterial orders tend to pair with veins typically two orders higher). At higher resolution, deviations could be seen for single segments where there appear to be two arterial vessels following a vein. However, on bifurcation of the vein, proper ratios are restored. More often a direct parallel architecture was seen. Figure 8 shows the ratio of *r*_{a}-to-*r*_{v} for various values of *r*_{v}. The fit, with constants *α* and *β*, is of the form *r*_{a}/*r*_{v} = exp(α *r*), chosen based on the observation that ∂*r*_{a}/∂*r*_{v} ≤ 1 in all parallel vessels. That is, local change in arterial radii is small relative to the change in parallel venous radii.

At a fundamental level, parallelism in vascular trees has been suggested as a means for prevention of vessel intersection during angiogenesis. However, functional hypotheses to vascular parallelism have been suggested. Thermoregulation has been seen as an explanation for parallelism, occurring primarily in arteries and veins in the same sheath (35). “Cross talk” or the so-called “counter-current exchange” hypothesis has been demonstrated in renal vasa recta and tubules (13) as well as in villi of the small intestine (14). Bassingthwaighte et al. (1) suggest a mechanical advantage for distributing blood within the microvasculature of the myocardium as a possible explanation for parallelism. The reason for vascular parallelism in the kidney has not been explained but could provide insight into overall renal function.

From a technical standpoint, the parallel vessel architecture had negative implications for both the imaging modality and segmentation method. As a consequence of the renal arterial and venous interactions, often the spatial scale between these parallel vessels was small as Fig. 2*C* demonstrates. Despite the fine voxel resolution of 4 μm, the distinction between arteries and veins was lost. Furthermore, one may notice from Fig. 2*C* that the influence of an ambiguous boundary between the vessels has a greater impact on the total diameter of the artery than the vein due to the shear size and shape of the two vessels. Thus extrapolating the arterial tree further became virtually impossible. However, because of this parallelism, it was assumed that the arterial tree should continue to parallel the venous tree. This assumption was implemented by way of cross correlating the terminal ordering profile of the venous tree with that of the arterial tree. This provides a more complete description of the renal arterial tree than could otherwise be given by using the parallelism inherent to the vasculature.

Another complication of renal vascular parallelism is that a major assumption of the segmentation technique was the spatial scale between vessels was greater than that of the vessel lumen. This contradiction resulted in a number of seed points that failed under tracking, requiring more computation time to obtain the entire network.

Figure 8 illustrates the average probability a vessel will remain its given order (*ρ*_{s}) for varying allowed error modulated by the radial factor (see methods). This process assumes there is some error due to pruning or in the calculated radius from segmentation and consequently the iterative terminal assignment. The probability of *ρ*_{s} for a specific order is not of general concern; however, it is the interplay of all orders that elude to errors and/or significant changes in ordering given errors image acquisition or methods applied. Because of the incomplete data sets, low *ρ*_{s} suggests areas in the network requiring further investigation.

Figure 9, *A* and *B*, shows the profile for the arterial and venous trees, respectively. This analysis predicts reasonable stability at high order even for large error due to a buffering affect of downstream operators, (i.e., the gray ‡ in Fig. 3). Thus we have confidence that the method is robust for the orders assigned via Strahler ordering. Specifically, these plots show that the relative change in order of each tree is not prone to drastic changes in the ordering due to a different terminal ordering profile.

Figure 11 illustrates the usefulness of statistically representing structural data and the efficacy of the techniques applied. The plot shows subtree volume for both the image reconstructions as well as a prediction based on the CM and vessel statistics. The prediction shows agreement and integration between both reconstructions.

The optimality comparison of Fig. 13 shows branch angle data from the three-dimensional anatomical data compared with the theories proposed by Zamir (47) and Zamir and Bigelow (48). The significance of such data is in understanding the means for which vessels branch within vascular systems. These results show a poor correlation with the “minimum shear force” hypothesis set forth. In fact, the use of a power of zero for the cost function resulted in a curve with a significantly lower mean standard error. These results suggest the need for more applicable cost functions, studies that are achievable with full three-dimensional reconstructions.

Murray's Law, which suggests that the relation between daughter and parent vessels is governed by minimizing work, was compared with data from both arterial trees. The deviation between the line of best fit and Murray's Law was approximately 1%. This provides convincing evidence that Murray's Law effectively characterizes the radii change across bifurcations over a range of spatial scales. Furthermore, this adds additional support of the minimum work hypothesis (31).

Blood flow to glomeruli and the resulting glomerular filtration rates remain difficult to calculate in vivo. However, estimates of renal perfusion are attainable through blood flow simulations. Models such as Poiseuille flow (49), lumped parameter models (30), and more complex models such as one-dimensional fluid flow through passive distensible vessel networks (39) have been proposed. These models require information on structure due to the strong dependence of vascular resistance on radius. Determining the blood flow through a vascular tree is dependent on connectivity, vessel length, and most significantl, radius. Correlations to total cross-sectional area (determined entirely by vascular radius) determine the affect of flow through the vascular network, where increasing total area is proportional to decreasing flow and visa versa (Fig. 12). The data from this study provide both a specific three-dimensional data set and the statistical information required for reconstruction and morphologically based blood flow analysis.

The renal vascular reconstructions developed through the course of this study will be freely available online for download and use, providing the academic community with a structurally accurate anatomical model. Not only useful for general structural analysis and basis for comparisons, it provides a modeling framework for further investigations into renal hemodynamics.

#### Applied methods.

Cast studies such as those by Kassab et al. (20, 21) have provided invaluable knowledge to the scientific community on vascular structure. Casting studies are manually intensive; however, they provide a complete description of vascular trees. The benefit in automatic segmentation techniques is the promise of a quick and efficient means for extracting relevant spatial and structural data from these images. At present, automatic segmentation of vascular images can only hope to extract portions of morphological data due a lack of established techniques and problems in processing medical images robustly. This paper attempts to address these issues that present significant challenges to scientific research in vascular morphology and cardiovascular disease states.

Strahler ordering has been used in a number of studies to organize structural information about vascular networks (7, 37, 45, 46). Kassab et al. (20, 21) and Jiang et al. (19) proposed a modification to the Strahler scheme by modifying orders according to the bounds as those written in *Eq. 1*. This analysis was applied to the renal arterial tree, with successful convergence; however, the 20-μm voxel venous tree did not achieve convergence after 20 iterations. This was due to instances where the absolute difference in neighboring mean radii was less than or equal to the absolute difference in their standard deviations. Thus the modification resulted in poor solutions or no solution at all. For these reasons, though the scheme worked well for the arterial tree, it was not used as it was not well suited for both arterial and venous trees.

The ordering error analysis applied to the renal vasculature provides a means for assessing the accuracy of the assigned orders by their probability of change. From a qualitative standpoint, the ordering profile seemed to agree well with known terminology and vessel morphology. However, because much of the process relies on developed statistics for terminal ordering, it is important to know the confidence in ordering results. The drawback of this approach is the assumption that the network does not contain operators intermittently that were not accounted for (i.e., vessel pruning occurred at terminals, not midway along network segments). Such added operators could alter the behavior of the upstream network but remain difficult to account for.

The unique approach given to commute networks via declaring terminal orders using statistical profiles has its limitations. Taking data on a specific subtree and applying it to an entire vascular network assumes a limited sample size at higher orders within the subtree as well as insignificant regional variation. Because of the structure of renal vasculature and the agreement of collocated data with both sets (Fig. 11), sufficient overlap seems to exist to corroborate this approach.

In conclusion, from micro-CT images of rat renal vasculature, high- and low-resolution data sets have been automatically segmented, and an anatomically accurate model of the arterial and venous trees was constructed. The segmented models have been combined and classified using Strahler ordering and the statistics of radius, length, and connectivity reported by order. Further development in techniques for assessing vascular images and data were described. These morphological data provide a statistical basis from which renal vascular topology can be generated. Furthermore, the three-dimensional geometry (available online) enables further analysis of global spatial properties and local branch angle relations, as well as a basis for simulations of renal perfusion.

## APPENDIX A

The notation used in the method below defines *a*→*b* as a length and direction going from *a* to *b*, “ ^ ” refers to a unit direction, and the subscripts *i*, *j*, and *k* refer to specific points within a defined group. From the group of points (Figs. 1*A* and 2*A*), each point *p*_{i} was selected and a local neighborhood *ω* of points near *p*_{i} defined. Figure 2*B* illustrates the point *p*_{i} and its neighborhood *ω* seen in dark grey. If the transmural direction *u* could be determined using the set *ω*, the inward direction could be obtained for a particular point. The approximate transmural direction *u* was determined as the average sum of cross products (formally stated in *Eq. A1*) between all combinations of points in *ω*. In *Eq. A1*, *p̂*_{i→j} and *p̂*_{i→k} refer to the direction from *p*_{i} to the points in *ω*, *p*_{j}, and *p*_{k}, *n*_{ω} is the number of points in *ω*, and AP_{r} is an additive projection function ensuring all directions point inward or outward from the vessel. (A1) Based on the assumptions outlined above, the closest point along ±*u* was assumed to reside on the opposing vessel wall as seen in Fig. 2*C*. There are a number of cases where this assumption is invalid (i.e., at a bifurcation or near parallel vessels) that are dealt with below. With the knowledge of the direction and distance to the opposing wall, the *i*th point is then projected to the center of the vessel as shown in Figs. 1*B* and 2*D*.

The new partially emaciated group (best illustrated in Fig. 1*B*) consists of points that lie along the centerlines of all vessels and are surrounded radially by the original point group. The further emaciation of the new point group is required to filter out errors in the initial projection. It is important to note that large errors in the initial projection lead to isolated or grouped points with little to no longitudinal neighbors, providing a quantitative means for error identification.

In the new partially emaciated point group, a neighborhood of points is established based on the information stored about the initial projection of *p*_{i} to its new position p̃_{i}. This group *ω′* is the collection of points within the assumed “radius” of the vessel, which is simply calculated as the distance of *p*_{i}→*p̃*_{i}. Applying *Eq. A2*, the resulting sum points transluminally. Although the group *ω′* contains scattered radial noise due to error in the initial projection of all points, the random occurrence of this error is relatively homogeneous and its effects on the calculation of *Eq. A2* approximately cancel. (A2) The estimate of the axial direction at *p̃*_{i} is useful because it not only orients the vessel direction, but defines the vessel cross section (i.e., Fig. 2*D* is oriented so the reader is viewing along the transluminal direction). Figure 2*D* illustrates the expected result when comparing the new point group to the original, where the point *p̃*_{i} is surrounded by the original set of vessel wall points. However, if the initial projection was incorrect, the point *p̃*_{i} would not have this relation to the original set. To utilize this relationship, a section of the vessel (chosen as twice the assumed “radius”) was used. The axial normal from *Eq. A2* was used to project points into a plane (as seen in Fig. 2*D*) calculated using *Eq. A3*, where *j* subscript refers to the *j*th point of the cross section from the original point group. (A3) Common to poorly projected points are holes in the wall cross section. Whereas the vessel bifurcations (seen in Fig. 2*A*) also produce holes in the wall, projecting points into a plane through a section of the bifurcation eliminates this problem. To identify holes numerically, the planar cross section was divided into angles as indicated by Fig. 2*D*. The angle φ_{j} is defined by the direction *p̃*_{i→j} and *p̃*_{i→φ} based on their dot product as detailed in *Eq. A4* and seen in Fig. 2*D*. However, the dot product only uniquely defines angles from 0 to π, thus it is also necessary to compare the vector from *p̃*_{i} to an in-plane point with the resultant vector of *p̃*_{i→φ} and the axial vector *u*_{z} (*Eq. A4*). (A4) In each section, the closest set of points to *p̃*_{i} were used to reposition *p̃*_{i} to the centerline of the vessel, and the radius was recalculated as the average distance of *p̃*_{i} → *p̃*_{j,2D}. Points that had holes or empty sections were eliminated from the new point group. By repeating this process, points converged toward the centerline of the vessels (shown in Fig. 1, *C*–*E*), and erroneous points were removed.

## APPENDIX B

Bifurcation operators are listed below for all possible combinations of input orders. It is assumed that the vector of probable outcomes from *i* is being modified by the operator formed from *j*.

To determine the probability *ρ*^{i} that segment *S*_{i} from Fig. 3 will increase, remain, or decrease its order, we can express it in terms of downstream probabilities. Here the superscript indicates the vessel segment from Fig. 3 seen as *Eq. B1*. (B1)

Labeling the operator *ζ* according to the upstream segment, the probability can be rewritten (*Eq. B2*). (B2) If we identify the pathway for this calculation *t*_{i}, which is identical to the sequence *S*_{i},*S*_{9},*S*_{6}, the total calculation is seen in *Eq. B3*. (B3) This calculation also relies on the calculation of *ρ*^{3}. *Equation B4* allows us to calculate the proper parameters for *ζ*_{i}. (B4) The combination of all these equations allows us to calculate how *S*_{i} will change in order due to any combination of terminal probabilities.

## GRANTS

D. A. Nordsletten is supported by scholarships from the New Zealand Mathematics Institute and its Applications and the University of Auckland International Doctoral Scholarship. N. P. Smith acknowledges support by the Marsden Fund of the Royal Society of New Zealand through Grant 04-UOA-177 and the National Institutes of Health (NIH) through the NIH/NIBB Multi-Scale Grant RO1-EB-00582501. J. Carlos Romero's lab, under the sponsorship of National Heart, Lung, and Blood Institute Grant HL-16496, performed the kidney preparation. High-resolution micro-CT data were obtained at the Brookhaven National Synchroton Light Source, Brookhaven National Laboratory and was supported by the United States Department of Energy, Division of Materials Science and Division of Chemical Sciences by contract DE-AC02-76CH00016. NIH Grant EB000305 sponsored lower-resolution micro-CT data collection.

## Acknowledgments

The authors thank the reviewers for insightful and constructive comments in regard to this work. Authors offer special thanks to P. E. Beighley and S. M Jorgensen for contributions to data collection, J. Lam for background study, and A. S. Beatty for helpful comments.

## Footnotes

↵1 http://www.physiome.org.nz/publications/nordsletten_blackett_ritman_bentley_smith_2005.

↵2 ImageJ is an image processing tool sponsored by National Institutes of Health and downloadable at http://rsb.info.nih.gov/ij/

↵3 CMGUI or Continuum Mechanics Graphical User Interface is a open source and available online at http://www.bioeng.auckland.ac.nz.

↵4 Estimates were determined based on fitting a histogram of terminal vessel radii to a normal distribution and ordering 0 for all radii within a standard deviation (SD).

↵5 The high and low radial bounds for an order are defined as shown in

*Eq. 1*, where*i*indicates the order of interest.↵6 Jiang et al. (19) applied

*Eq. 1*to reduce the radii overlap between successive orders and modify the process of assigning Strahler orders. It is important to note that*Eq. 1*is applied, in this study, in a different context.

- Copyright © 2006 by the American Physiological Society