Skip to main content

Scale-adaptive surface modeling of vascular structures



The effective geometric modeling of vascular structures is crucial for diagnosis, therapy planning and medical education. These applications require good balance with respect to surface smoothness, surface accuracy, triangle quality and surface size.


Our method first extracts the vascular boundary voxels from the segmentation result, and utilizes these voxels to build a three-dimensional (3D) point cloud whose normal vectors are estimated via covariance analysis. Then a 3D implicit indicator function is computed from the oriented 3D point cloud by solving a Poisson equation. Finally the vessel surface is generated by a proposed adaptive polygonization algorithm for explicit 3D visualization.


Experiments carried out on several typical vascular structures demonstrate that the presented method yields both a smooth morphologically correct and a topologically preserved two-manifold surface, which is scale-adaptive to the local curvature of the surface. Furthermore, the presented method produces fewer and better-shaped triangles with satisfactory surface quality and accuracy.


Compared to other state-of-the-art approaches, our method reaches good balance in terms of smoothness, accuracy, triangle quality and surface size. The vessel surfaces produced by our method are suitable for applications such as computational fluid dynamics simulations and real-time virtual interventional surgery.


In surgical planning, treatment evaluation, and medical education, the geometric modeling of vascular structures is of vital importance. Three-dimensional (3D) models can help surgeons better understand the branching patterns and complex topology of vascular structures in a short time for better and quick decision making during surgery by providing straightforward information on the morphology of the vessels, the spatial relationships among these vessels and other relevant anatomic structures, and an intuitive depiction of curvature and depth relations [14].

The surface modeling techniques of vascular tree structures can be broadly classified as either model-based or model-free techniques [1]. Generally, the former methods require vessel centerline extraction and vessel diameter determination from segmented vessels [5]. Based on the centerline model (defined by the centerline and radius), geometric primitives such as cylinders [6] and truncated cones [7] are employed to fit the vessel surface for visualization. Unfortunately, the smoothness of the surface produced by these methods is poor, especially where the vessel branches. At these points, transition is unavoidably discontinuous and therefore has significant artifacts, resulting in very low visual quality. To achieve high-quality surface, other advanced surface representations have been investigated, such as, B-spline surfaces [8], simplex meshes [9], convolution surfaces [10], and subdivision surfaces [11, 12]. Most of these methods yield desirable smooth surfaces; however, they suffer from low accuracy, badly shaped triangles or a large number of polygons.

The most common model-free technique of surface reconstruction in medical visualization is Marching Cubes (MC) [13]. Although effective in capturing overall shape, this technique has two major limitations. One limitation is that the generated surface heavily relies on the chosen isovalue and a slight change in value may result in great change in both the topological and geometrical features of the generated surface. The other limitation is that the visual quality is very low because the generated surface contains strong aliasing artifacts. Furthermore, these artifacts may lead to unstable numerical problems when the generated surface is applied for computational fluid dynamics (CFD) simulations. Recently, Schumann et al. [14, 15] presented a model-free technique that could produce smooth surface from vessel segmentation. The technique is based on multi-level partition unity (MPU) implicits [16] originally dedicated to reconstruct the surface from 3D point clouds.

However, the main drawback of model-based methods is that their assumed models are unable to represent the underlying image data, and are therefore inappropriate for vessel diagnosis where high accuracy of surface representation is required. These methods assume that the cross-section of vessels is circular, whereas the pathologic vessels in clinical practice such as aneurysms might generally have a non-circular shape (e.g. ellipse) [2, 5]. In contrast, model-free methods make no model assumptions and represent the underlying data with high fidelity. Therefore the reconstructed surface from model-free methods could be used for vessel diagnosis.

this paper, we present a model-free approach that relies on a prior vessel segmentation result, point extraction, Poisson equations and adaptive polygonization. With the proposed weight function, the triangulation algorithm in the gap-stitching stage can produce a two-manifold triangulation that maximizes the minimal angle of the triangle. Our approach yields a both morphologically correct and topologically preserved two-manifold smooth surface that is scale-adaptive to the local curvature of the surface by increasing/decreasing the size of triangles in regions with low/high curvatures. In addition, our method generates fewer and better-shaped triangles that are suitable for applications such as CFD computations and finite element analysis, and does not require other additional geometry processing techniques to improve triangle quality or to reduce number of triangles.

This paper is organized as follows. Details of our method are described in Sections Method. The results and discussion are presented in Section Results and Discussion. Finally, our conclusions are given in Section Conclusions.



The pipeline of the presented approach to the geometric modeling of vascular tree structures is illustrated in Figure 1. The pipeline begins with a 3D binary volume data from segmented image (1 represents the voxels of the vessel structures and 0 represents the background). The boundary voxels between the segmented vessel and the background are extracted and used to build a 3D point cloud. Then, a 3D implicit indicator function is computed from the oriented 3D point cloud by solving a Poisson equation. Finally the surface model is generated by a proposed polygonization algorithm for explicit 3D visualization.

Figure 1

Pipeline of the geometric modeling of vascular tree structures. Pipeline of the geometric modeling of vascular tree structures

Point extraction

The point extraction step aims to faithfully represent the boundary of the vessel using 3D point clouds. Generally, point generation from volume data is driven by the voxel grid of the segmented result [14, 17]. Braude et al. [17] used the boundary voxels to directly generate the point cloud; unfortunately their work failed to represent very thin vessel branch which might be represented as lines. To reconstruct the very thin vessel structures and to prevent aliasing artifacts, we use a similar adaptive point extraction technique [14, 15] in our pipeline. The technique relies on the constellation of adjacent object voxels and outer boundary voxels that are closest to the given object voxels, e.g. if there is only one object voxel in the 3D-6-neighborhood, one point will be generated in the center of the boundary face, which is defined as the voxel adjacent to the object neighbor voxel, as illustrated in Figure 2. To sufficiently represent the thin vascular branch, voxels representing thin structures are first identified by a top-hat-transformation with a 3x3x3 structuring element, after which all outer boundary voxels adjacent to thin structures are refined into eight subvoxels [14].

Figure 2

Example of point extraction. Example of generating point (blue) based on constellations of object voxels (grey) and outer boundary voxel (blue)

Normal vector estimation

To define a vector field for computing the indicator function, we need to estimate the normal vectors of the extracted point cloud. Generally, a normal vector is computed based on the image gradient of a given voxel in the segmentation result. However, this method may not be effective when the neighborhood of a voxel is symmetrical because the points are placed in the centers of the boundary faces. In this case, according to Schumann et al. [14, 15], the normal of the face is taken as the normal of the generated point. This method is simple; however, it may achieve undesirable results when the voxel size is large. To avoid this problem and robustly estimate the normal vector for a given sample point, we use a method relying on covariance analysis [18]. The 3x3 covariance matrix C for a sample point P is defined as follows:

C = [ p i , 1 p ¯ p i , k p ¯ ] T [ p i , 1 p ¯ p i , k p ¯ ] p ¯ = j = 1 k p i , j

where {p i,1, p i,2,, p i, k }is the k nearest neighbors of the point P. Since C is symmetric and positive semi-definite, all eigenvalues are real-valued and all eigenvectors form an orthogonal frame. The eigenvector corresponding to the smallest eigenvalue is taken as the normal vector of the point P.

Indicator function

After obtaining the point clouds with oriented normals, many techniques are available to reconstruct them for surface visualization. For a comprehensive introduction to these techniques, please refer to a recent survey [19]. One popular surface reconstruction technique is the implicit function technique. This technique first constructs a 3D function that approximates/interpolates the point samples and then polygonizes the reconstructed surface. The main advantage of this type of technique is that it can reconstruct a watertight surface from different 3D models with any topological complexity. We choose the Poisson surface reconstruction [20] technique to model the vessel surface because, unlike the MPU technique [16], it is robust to recover fine details from noisy data and does not need to resort to heuristic partitioning or blending for surface fitting.

The basic idea behind Poisson surface reconstruction is to utilize the vector field V to compute the indicator function ψ (defined as 1 at points inside the surface and 0 at points outside). At points near the surface, the gradient of ψ is a vector field that is equal to the normal vector field. Hence, the problem of computing the ψ turns into finding a function whose gradient best approximates the V , i.e. min ψ | | ψ V | | . After applying the divergence operator, the problem becomes a standard Poisson problem [20]:

Δ ψ ψ = V

To reconstruct fine details, an adaptive octree ξ defined by the position of the sample points is used to represent the implicit function, and each node O ξ of the tree is associated with a function F o when the following conditions are satisfied [20]: 1) the vector field can be expressed as the linear sum of the F o ; 2) the matrix representation of the Poisson equation can be solved efficiently; and 3) the representation of the indicator function can be accurately evaluated near the surface. For every node o ξ, the F o is set to be the unit-integral centered about the node o and scaled by the size of o:

F o ( q ) = F ( q c w ) 1 w 3 , F ( q ) = f ( q 2 d )

where q is the sample point; c and w are the center and width of the node O, respectively; d is the maximum tree depth; and f is a Gaussian filter with unit variance.

To allow for sub-node precision, the gradient field of the indication function is defined as follows:

V ( q ) = s S o Ω ϖ F o ( q ) n

Where S is the input data with a set of samples s S, each consisting of a point and an inward-facing normal n ; Ω is the set of eight depth d nodes closest to the sample point; and ϖ is a trilinear interpolation weight. After defining the vector field, the indicator function is obtained by solving the Poisson equation Eq. (2) using a conjugate gradient solver.


To explicitly visualize the implicit surface, we need to triangulate the implicit surfaces from the computed indicator function. Our polygonization algorithm consists of two stages: mesh-expanding and gap-stitching.

Mesh-expanding stage

We first introduce how to compute the radius of the curvature of a given point and then describe the mesh-expanding procedure. We choose a point x near the surface, and compute its corresponding surface point p, and the surface normal at p, denoted by n p , using an iterative procedure known as Newton step [21]. The radius of curvature at point p is estimated by calculating the radius of curvature of several geodesics that cover p, and taking the minimum one [22, 23]. Assuming a set of surface point q i with surface normal n q i , are close to the point p, let d i be the distance between p and q i , and θ i be the angle between n p and n q i , as illustrated in Figure 3. The radius of the curvature at point p is then calculated by

r ( p ) = min ( d i 2 sin ( θ i 2 ) )
Figure 3

Illustration of calculating the radius of curvature at a given point. The illustration of calculating the radius of curvature at a given point

Starting from point p, we construct six triangles as an initial mesh. First, the point p is surrounded by constructing a regular hexagon q 1, q 2,,q 6 in the plane tangent to p. To adapt to the local curvatures of the surface, we project q 1, q 2,,q 6 onto the surface to estimate r, the radius of curvature at point p. Then, the generated regular hexagon is adjusted, whose edge lengths are ρ·r, where ρ is a user-defined constant. Finally, the vertices of the hexagon are again projected onto the surface, and denoted as p 1, p 2,,p 6. The triangles formed by the hexagon p 1, p 2,,p 6 with p are the first six triangles of the mesh (see Figure 4). This mesh is considered as the seed element and is expanded by progressively growing triangles from its boundary edges.

Figure 4

The generation process of initial mesh. First generating a regular hexagon (the dashed lines in blue) in the plane tangent to p, and projecting the vertices of the hexagon onto implicit surfaces to estimate curvature. Then adjusting the edge length of the hexagon to ρ·r, and finally, projecting the vertices of hexagon onto surface again to generate an initial mesh, consisting of vertices p 1,,p 6 and p.

Once the initial mesh is generated, its boundary edges are placed into a queue. We then take an edge, denoted as (u, v), from the queue. A new point p is created if u, v and p form an equilateral or close-to-equilateral triangle. The new triangle Δuvp is coplanar and opposite to the existing triangle containing the edge (u, v). After point p is placed onto the implicit surfaces to estimate curvature, it is changed on the original plane such that the lengths of the edges (p, u) and (p, v) are equal toρ·r. Finally, the p is once again placed onto the implicit surfaces (Figure 5).

Figure 5

A fragment illustrating the mesh-expanding procedure. The red lines are boundary edges of the mesh and stored in a queue. From the boundary edges, new triangles are progressively generated. Note that the blue lines combined with one boundary edge will construct a new triangle if it satisfies two specified conditions. If two boundary edges make an angle less than 70 degrees, the three vertices on these two edges are used to produce a new triangle, see the triangle consisting of the green dashed line and the other two boundary edges.

The newly constructed triangle Δuvp is added to the mesh if it satisfies two conditions. One is that both edge (p, u) and edge (p, v) should make an angle of at least 50 degrees (the maximal angle is 70 degrees) with the edge (u, v) in the old mesh. The other condition is that the triangle Δuvp should not approach existing triangles too closely. These two conditions guarantee that the resulting triangles are close-to-equilateral and the gap generated by this stage is not too narrow to sew in the subsequent gap-stitching stage. If any one of the boundary triangle in the mesh (denoted as T) is closer to triangle Δuvp than one-third of the length of the longest edge in T and triangle Δuvp, then the triangle Δuvp is not added to the mesh. Otherwise, the triangle Δuvp is added to the mesh and the boundary edges (p, u) and (p, v) are placed into the queue. The mesh-expanding terminates when the queue is empty.

Gap-stitching stage

To stitch the gap produced in the mesh-expanding stage, as shown in Figure 6 (left), it should first be identified. Since the gap is a closed loop of boundary edges in the mesh, and the mesh produced after mesh-expanding stage is a connected two-manifold surface, the gap can be stitched as described below.

Figure 6

Polygonization of a trifurcate model. A long gap is produced upon the termination of mesh expanding stage (left), and is sewed in the subsequent gap-stitching stage (right).

Let Φ(i, j, k) be a weight function defined on the set of all triangles (p i , p j , p k ) that could be possibly generated in triangulating the polygon p i ,, p j , 0 ≤ i < j < n and let w i, j be the minimum total weight that can be achieved during the triangulation process. Then the triangulation algorithm proceeds as follows:

  1. i.

    For i = 0, 1,, n - 2, let w i, i+1= 0, and for i = 0,1,,n - 3, let w i, i+2: = Φ(i, i+1, i+2). Set j: = 2.

  2. ii.

    Put j: = j+1. For i = 0,1,,n - j - 1 and k = i + j. Let w i , k : = min i < m < k [ w i , m + w m , k + Φ ( v i , v m , v k ) ] . Let O i, k be the index m where the minimum is achieved.

  3. iii.

    If j < n - 1, then return to step 2; otherwise the weight of the minimal triangulation is w 0,n-1.

  4. iv.

    Let Θ: = φ and call the recursive function Trace with the parameters (0, n- 1).

Function Trace (i, k):

If (i + 2) = k, then Θ: = ΘΔv i v i+1 v k ;


  1. 1)

    Let o: = O i, k ;

  2. 2)

    If oi + 1, then Trace (i, o);

  3. 3)

    Θ: = ΘΔv i v o v k ;

  4. 4)

    If ok - 1, then Trace (o, k)

End else.

At the end of the algorithm, Θ contains the required triangulation of polygon p 0,,p n-1. The triangulating steps are similar to the schemes described in [24, 25], but different in defining the weight function Φ(i, j, k). Barequet and Sharir [24] suggested the function as the area of triangle (i, j, k), whereas Liepa [25] designed it by combining the dihedral angles between the neighboring triangles with areas of triangles. However, when the holes are highly irregular, the resulting mesh may not be a topologically preserved two-manifold surface. To avoid this, we propose a strategy wherein the angles of potential triangle (p i , p k , p j ), i < k < j, denoted as T, are first taken into account, with the minimal angle of triangle T being maximized. The dihedral angles and area are then considered the same time. Therefore, we define triples as follows:

Φ ( i , k , j ) = ( α , β , A )

where α is the maximized minimal angle of triangle T, β is the maximal dihedral angle between triangle T and its neighborhoods, and A is the area of triangle T. The ordering in Φ is designed to give precedence to α over β, and β over A:

( α 1 , β 1 , A 1 ) < ( α 2 , β 2 , A 2 ) : ( ( α 1 > α 2 ) ( α 1 = α 2 β 1 < β 2 ) ( α 1 = α 2 β 1 = β 2 A 1 < A 2 ) )

The addition operator sums the area but retains the maximized minimal angle and the "worst" (i.e., largest) dihedral angle:

( α 1 , β 1 , A 1 ) + ( α 2 , β 2 , A 2 ) : = ( m 1 , m 2 , A 1 + A 2 )

where m 1 = min(β 1, β 2), and m 2 = max(β 1, β 2).

With our weighting function, the triangulation algorithm can produce a two-manifold triangulation that maximizes the minimal angle of triangle T. After triangulation, the patching triangles are subdivided to make their density similar to the density of surrounding mesh [25, 26]. An example of gap stitching is illustrated in Figure 6 (right).

Results and Discussion

We have applied the presented model-free approach to the geometric modeling of a variety of vascular trees, namely, cerebral (Figure 7), liver (Figure 8) and aorta trees (Figure 9). The binary segmentation results of the liver and cerebral tree are from our manual segmentation whereas the aorta tree is from the public resource ( Table 1 summarizes the properties of the data tested in this paper. The produced surfaces, as shown in Figure 7, are smooth, especially at transition, and do not show aliasing artifacts. A close visualization shows that the morphology of vascular structures and thin structures (even the thin elongated structures) can be reconstructed with the correct topology.

Figure 7

A cerebral vessel surface model produced by our approach. A cerebral vessel surface model produced by our approach

Figure 8

Comparison of the geometric modeling results on a liver tree. Surface model generated by the MC algorithm (a), MPU-based algorithm (b) and our approach (c). Color-coded visualization of the root mean square curvature distribution for the generated surface using MC algorithm (d), MPU-based algorithm (e) and our approach (f).

Figure 9

Comparison of triangle quality for an aorta tree. Surface model generated by the MC (a), MPU-based method (b), SS-based method (c) and our method (d). The bottom row is a zoomed region corresponding to the rectangle region of the top row.

Table 1 Summary of properties of the data sets Sample table title

We evaluated our method in terms of surface smoothness, surface accuracy, triangle quality, surface size and efficiency on the tested dataset. We also compared our approach with the conventional model-free algorithm, i.e. MC algorithm, and state-of-the-art algorithms, i.e. model-free MPU-based algorithms (MPU-based) [14, 15], and model-based subdivision surface algorithm (SS-based) [11, 12]. In our experiments, the implementation of SS-based is slightly different than in [11, 12]. In [11, 12], after obtaining an initial mesh from the centerline model, the Catmull-Clark scheme [27] is applied to generate the vessel surface; thus the surface is a quadrilateral mesh. Since the surfaces produced by MC, MPU-based and our method are all triangular meshes, for the convenience of comparison, we applied the Loop scheme [27] with three iterations to produce vessel surface (for regular meshes, the surfaces yielded by the Loop scheme and Catmull-Clark scheme are both C2-continuity [27]). For the MPU-based algorithm, we attempted to select parameter settings suggested in [14, 15] with the best results. In our approach the user-defined ρ is set to 0.15 in the mesh-expanding stage and the parameters k is set to 10 in the normal vector estimation stage.

Surface smoothness

We compared our approach with MC and MPU-based algorithm when applied to the same segmentation result. As illustrated in Figure 8, the surface produced by the MC (Figure 8a) has a visually low surface quality and contains a great variety of artifacts, which might disturb the visual interpretation of the vessel surface and therefore affect decision making during diagnosis. In contrast, the surfaces generated by both the MPU-based (Figure 8b) and our approach (Figure 8c) are highly smooth.

To validate the smoothness of the surface, we computed the distribution of the curvature value on the surface. The curvature map can directly help us observe the flaws and roughness of the surface that are not easily identified by human eyes. Here, we computed the root mean square (RMS) curvature of both the maximal k max and the minimal k min principal curvatures respectively, and RMS is defined as ( k 2 max + k 2 min ) / 2 . The principal curvatures are computed based on a finite-differences technique [28].

Figure 8d, Figure 8e, and Figure 8f show the curvature distribution of the MC surface, MPU-based surface and the surface generated by our method respectively. It can be seen that due to the substantial staircase artifacts, the curvature distribution of the MC surface is highly inhomogeneous compared to that of both the MPU-based method and our method. For the thin vessel structures, the distributions of other methods are similar to that of the MC surface. The difference of surface smoothness between our approach and MPU-based algorithm is not very apparent. The reason is that these two approaches both utilize an implicit descriptor as the underlying surface representation. Although the smoothness of the MC surface can be improved with additional geometry processing techniques, e.g. smoothing filter, this unfortunately leads to volume shrinkage, collapse of thin vessel structures and unfaithful representation of the underlying data [29].

3.2 Surface accuracy

We analyzed the accuracy of the generated surface on the assumption that the input binary segmentation result of our pipeline had been validated correctly. To provide a quantitative comparison of the surfaces, we measured the distance error between two surfaces using the MESH tool [30]. The tool utilizes Hausdorff distance to calculate the maximum, mean and RMS errors between two specified surfaces. In this experiment, the surface generated by the MC algorithm is taken as reference surface (although it is not the most accurate technique to visualize the segmentation result, it has been the de facto standard in medical surface visualization and has been widely applied in numerous radiological workstations [5, 15]).

Table 2 lists the mean, maximum and RMS errors for the tested dataset. The maximum error of SS-based method is larger than twice a voxel size, and is also larger than that of both the MPU-based and our method, whose maximum error is less than a voxel size. The reason is that, aside from the simplified model assumptions of circular cross-sections, approximating subdivision scheme leads to volume shrinkage for a closed surface during its convergence to the limit surface. Due to the same point extraction strategy, the errors between MPU-based and our method are very similar. However, the maximum error remains larger than half of a voxel size. This occurs in the feature regions of the vessel surface, such as small concave and convex regions, which are not represented by sufficient points, even though an adaptive subsampling technique is applied [14, 15].

Table 2 Accuracy of MPU-based method, subdivision surface-based method and our method for the tested dataset Mean, max and RMS denote mean distance error, maximum distance error and root mean square distance error.

Triangle quality

Like surface accuracy, triangle quality is an important factor in achieving accurate results in many simulations. CFD simulations require the input surface to be free from block and staircase artifacts; therefore, the triangle meshes should have a good quality with regard to edge ratio [31]. Degenerated triangles such as thin and elongated triangles may lead to numerical unstability in CFD simulation, and may even make the simulation impossible. Smooth transition at the points where the vessel surface branches is also a prerequisite. Furthermore, the triangle size should not change abruptly, and surface regions with high curvature are desirably represented by small triangles.

Figure 9 demonstrates the comparison of triangle quality for the aorta tree. It can be seen once again from the zoomed region that the MC surface contains strong staircase artifacts and badly shaped triangles, and is not naturally smooth at the transition of the branches. The surface generated by the MPU-based method has a smooth transition, but also contains degenerated triangles. Although these triangles can be removed by invoking additional mesh quality improvement techniques [15], special care must be taken to preserve vital surface features during the optimization process. Because of the good underlying property of subdivision surface, the surface produced by the SS-based method is composed of well-shaped triangles, with smooth transition at the branches. However, the surface accuracy of the SS-based method is low (see Section Surface size) and unsuitable for CFD simulations [31]. Similar to the SS-based method, the surface of our method also has a smooth transition without badly shaped triangles. Additionally, the triangle size yielded by our method is adaptively scaled to the local differential geometric surface characteristics. The surface areas with high curvature are represented by smaller triangles, whereas triangles become large in the relative low-curvature region. Meanwhile, the triangle size from small to large is changed gradually.

Edge ratio is one measure for triangle quality, and is defined as τ = |t|0/|t|, where |t|0 is the minimum edge length of a given triangle, and |t|is the maximum length [32]. It is straightforward to see that τ ≤ 1, whereas the equality occurs when the triangle is equilateral, and τ converges to zero if the triangle is highly needle-shaped. An ideal triangular mesh has triangles that are all equilateral, and should adapt to local surface properties, such as curvature. Figure 10 shows the distributions of edge ratio for aorta tree. It can be seen that the MC and MPU-based methods yield approximately 5% very badly shaped triangles (τ < 0.1) such as thin elongated triangles, and less than 6% well-shaped triangles (τ ≥ 0.8). Moreover, these two methods suffer from enormous standard deviations of ratio. In contrast, approximately 50% of generated triangles (τ ≥ 0.8) by SS-based and our methods are close-to-equilateral. Unsurprisingly, the two latter methods with small standard deviation yield no thin elongated triangles. However, our method produces some triangles (less than 1%) with relative small ratio (τ = 0.3). These triangles are constructed to stitch the gap in the polygonization.

Figure 10

The distribution of edge ratio of the aorta tree. The distribution of edge ratio of the aorta tree.

Surface size

Here surface size refers to the total number of triangles and vertices approximating a surface. Surface size affects surface accuracy, surface rendering speed, and human interactive response. In the MPU-based method, Bloomental's implicit polygonizer [33] is applied to produce a triangular mesh. In the polygonization, if the grid size is set too large, the size of generated triangle is also very large, resulting in a loss of finer details such as thin vessel structures. If the grid size is too small, the polygonization process will be time-consuming; the surface size will be huge and will slow down interactive rendering frame rates. Mesh simplification techniques are therefore usually invoked as a subsequent step to reduce the surface size. However, this step may lead to a loss of both topological and geometrical features, and may even produce degenerated triangles during the simplifying process. In the SS-based method, the smoothness of the surface is increased with the iteration of subdivision; unfortunately, the number of polygons grows exponentially. Each iteration of subdivision yields a three-time increase in polygon due to the underlying topological refinement rules. In our method, the size of triangle can be adapted to the local curvature of the surface. This feature can save many triangles in representations. In actuality, using many small triangles to represent surface regions with low curvature, such as a flat region, does not significantly improve surface smoothness but increases surface size.

Table 3 reports several statistics for the surface sizes of the three tested data. Compared with the MC and MPU-based method, our method generally saves approximately 20% and 10% in the number of triangles, respectively, whereas for a vessel tree consisting of many highly curved branchings such as cerebral tree, the saving is more than 30% and 20% respectively. The surface size produced by the SS-based method at third iteration is smaller; however, it will be greatly larger at fourth iteration than that of MC and MPU-based method. Taking the cerebral tree as an example, the number of vertices and triangles after four iterations of subdivision are up to 102443 and 204976, resulting in triangles smaller than the original voxel.

Table 3 Surface sizes for three vessel surfaces generated by the MC, MPU-based method, SS-based method and our method

Computational efficiency

Our method utilizes implicit function to describe vascular structures; therefore, it requires an evaluator for the indicator function defined at all extracted points in space. The function is obtained by solving Poisson equations using the efficient linear solvers [34]. The computational cost mainly depends on the complexity and the resolution of input objects. In the mesh-expanding stage, the time is largely spent on calculating curvature radius and movement of points onto surface. However, in the gap-stitching stage, due to the O(n3) performance complexity of the triangulation algorithm, the time for this stage comprises approximately one-fifth of the entire time. However, the proposed approach is much slower when compared to the MC, MPU-based or SS-based method. Taking the cerebral tree, the most complex tested dataset, as an example, the overall time for our method is 127 seconds, whereas the MC only requires 5 seconds. The performances of this data for the MPU-based method and SS-based method are 53 and 66 seconds respectively.


We have presented a model-free method for the geometric modeling of vascular structures. Our method yields both a topologically correct two-manifold and a geometrically smooth vessel surface. An important feature of the presented method is that it produces a surface that is scale-adaptive to the local curvature of the surface. This minimizes the number of triangles in the representation, leading to faster interactive rendering frame rates and saving much computational time in post-processing procedures, such as real-time blood flow simulations and collision detections in virtual interventional surgery.

We validated our method to a variety of vascular structures and compared the results with other state-of-the-art techniques, both model-based techniques and model-free techniques, in terms of surface smoothness, surface accuracy, triangle quality, surface size and efficiency. Compared to the MC and MPU-based methods, the surface generated by our method achieves comparable accuracy; however, it is more suitable for applications that require high-quality triangulations such as CFD computations and finite element analysis, because our method yields smaller surface size, better-shaped triangles and no thin elongated triangles. Therefore, invoking additional geometry processing techniques to improve mesh quality or to reduce surface size is not necessary for the presented method. Mode-based methods, such as the SS-based method, can produce smooth surface and well-shaped triangles. The simple circular model assumption results in a low accuracy that is inappropriate for vessel diagnosis or CFD simulations, but can be used for certain situations where accuracy is not very important, such as in medical educations. Fortunately, very recent work [35] showed that with an elliptical model assumption, the surface accuracy of SS-based method could be improved. The experimental results demonstrate that our method reaches a better balance with regard to surface accuracy, surface smoothness, triangle quality and surface size.

The investigation of computational efficiency has revealed the limitation of our method in its current implementation. Fortunately, implementing the time-consuming steps, such as the triangulation step on CUDA, a parallel computing engine developed by NVIDIA [36], seems to be a promising solution to the limitation and might be part of future work. Although the presented method makes no model assumption and achieves high accuracy, it does not imply that our method can be directly applied to diagnostic tasks, because our pipeline takes the binary segmentation result as input, and supposes that the segmentation is validated correctly. Therefore, combining the validation of input data with the pipeline as a preprocessing step is also planned for future studies.


  1. 1.

    Oeltze SB, Preim B: 3D visualization of vasculature: an overview. Visualization in medicine and life sciences 2007, 19–39.

    Google Scholar 

  2. 2.

    Piccinelli M, Veneziani A, Steinman DA, Remuzzi A, Antiga L: A framework for geometric analysis of vascular structures: application to cerebral aneurysms. IEEE Trans Med Imag 2009, 28(8):1141–1155. 10.1109/TMI.2009.2021652

    Article  Google Scholar 

  3. 3.

    Joshi A, Qian X, Dione DP, Bulsara KR, Breuer CK, Sinusas AJ, Papademetris X: Effective visualization of complex vascular structures using a non-parametric vessel detection method. IEEE Trans Visual Comput Graph 2008, 14(6):1603–1610. 10.1109/TVCG.2008.123

    Article  Google Scholar 

  4. 4.

    Ding S, Ye Y, Tu J, Subic A: Region-based geometric modelling of human airways and arterial vessels. Comput Med Imag Graph 2010, 34(2):114–121. 10.1016/j.compmedimag.2009.07.005

    Article  Google Scholar 

  5. 5.

    Volkau I, Ng TT, Marchenko Y, Nowinski WL: On geometric modeling of the human intracranial venous system. IEEE Trans Med Imag 2008, 27(6):745–751. 10.1109/TMI.2007.911004

    Article  Google Scholar 

  6. 6.

    Masutani Y, Masamune K, Dohi T: Region-growing-based feature extraction algorithm for tree-like objects. Proc of visualization in biomedical computing 1996, 161–171.

    Google Scholar 

  7. 7.

    Hahn H, Preim B, Selle D, Peitgen H: Visualization and interaction techniques for the exploration of vascular structures. IEEE Visualization 01 2001, 395–402.

    Google Scholar 

  8. 8.

    Höhne K, Pflesser B, Pommert A, Riemer M, Schubert R, Schiemann T, Tiede U, Schumacheret U: A realistic model of the inner organs from the visible human data. Proc of the international conference on medical image computing and computer-assisted intervention 2000, 776–785.

    Google Scholar 

  9. 9.

    Bornik A, Reitinger B, Beichel R: Reconstruction and representation of tubular structures using simplex meshes. Proc of winter school of computer graphics 2005, 611–65.

    Google Scholar 

  10. 10.

    Oeltze S, Preim B: Visualization of vascular structures with convolution surfaces: method, validation and evaluation. IEEE Trans Med Imag 2005, 25(4):540–549. 10.1109/TMI.2004.843196

    Article  Google Scholar 

  11. 11.

    Felkel P, Wegenkittl R, Bühler K: Surface models of tube trees. Proc of the computer graphics international 2004, 70–77. full_text

    Google Scholar 

  12. 12.

    Luboz V, Wu X, Krissian K, Westin CF, Kikinis R, Cotin S, Dawson S: segmentation and reconstruction technique for 3D vascular structures. Proc of the international conference on medical image computing and computer-assisted intervention 2005, 43–50.

    Google Scholar 

  13. 13.

    Lorensen W, Cline H: Marching cubes: a high resolution 3D surface construction algorithm. Proc of ACM SIGGRAPH 1987, 163–169. 10.1145/37402.37422

    Google Scholar 

  14. 14.

    Schumann C, Oeltze S, Bade R, Preim B, Peitgen H: Model-free surface visualization of vascular trees. Eurographics/IEEE-VGTC symposium on visualization 2007, 283–290.

    Google Scholar 

  15. 15.

    Schumann C, Neugebauer M, Bade R, Preim B, Peitgen H: Implicit vessel surface reconstruction for visualization and simulation. Int J CAS 2008, 2(5):275–286.

    Google Scholar 

  16. 16.

    Ohtake Y, Belyaev A, Alexa M, Turk G, Seidel H: Multilevel partition of unity implicits. ACM T Graph 2003, 22: 463–470. 10.1145/882262.882293

    Article  Google Scholar 

  17. 17.

    Braude I, Marker J, Museth K, Nissanov J, Breen D: Contour-base surface reconstruction using MPU implicit models. Graph Models 2007, 69: 139–157. 10.1016/j.gmod.2006.09.007

    Article  Google Scholar 

  18. 18.

    Pauly M, Gross M, Kobbelt L: Efficient simplification of point-sampled geometry. IEEE Visualization 02 2002, 163–170. full_text

    Google Scholar 

  19. 19.

    Schall O, Samozino M: Surface from scattered points: a brief survey of recent developments. 1st International workshop on semantic virtural environments 2005, 138–147.

    Google Scholar 

  20. 20.

    Kazhdan M, Bolitho M, Hoppe H: Poisson surface reconstruction. Proc of the fourth Eurographics symposium on geometry processing 2006, 61–70.

    Google Scholar 

  21. 21.

    Hartmann E: A marching method for the triangulation of surfaces. The Visual Comput 1998, 14(3):95–108. 10.1007/s003710050126

    Article  Google Scholar 

  22. 22.

    Araujo B, Jorge J: Curvature dependent polygonization of implicit surfaces. Proc of 17th Brazilian symposium on computer graphics and image processing 2004, 266–273. full_text

    Google Scholar 

  23. 23.

    Karkanis T, Stewart AJ: Curvature-dependent triangulation of implicit surfaces. IEEE Comput Graph 2001, 22(2):60–69. 10.1109/38.909016

    Article  Google Scholar 

  24. 24.

    Barequet G, Sharir M: Filling gaps in the boundary of a polyhedron. Comput Aided Geom Des 1995, 12(2):207–229. 10.1016/0167-8396(94)00011-G

    MathSciNet  Article  Google Scholar 

  25. 25.

    Liepa P: Filling holes in meshes. Eurographics symposium on gemetry processing 2003, 200–205.

    Google Scholar 

  26. 26.

    Shewchuk JR: Triangle: engineering a 2 d quality mesh generator and delaunay triangulor. Proc of the 1st workshop on applied computational geometry 1996, 123–133.

    Google Scholar 

  27. 27.

    Zorin D, Schroder P: Subdivision for modeling and animation, SIGGRAPH 2000 course notes. 2000, 69–77.

    Google Scholar 

  28. 28.

    Rusinkiewicz S: Estimating curvatures and their derivatives on triangle meshes. Symposium on 3D data processing, visualization, and transmission 2004, 486–493. full_text

    Google Scholar 

  29. 29.

    Bade R, Haase J, Preim B: Comparison of fundamental mesh smoothing algorithms for medical surface models. Simulation and visualization 2006 2006, 289–304.

    Google Scholar 

  30. 30.

    Cignoni P, Rocchini C, Scopigno R: Metro: measuring error on simplified surfaces. Comput Graph Forum 1998, 17(2):167–174. 10.1111/1467-8659.00236

    Article  Google Scholar 

  31. 31.

    Cebral JR, Castro MA, Appanaboyina S, Putman CM, Millan D, Frangi AF: Efficient pipeline for image-based patient-specific analysis of cerebral aneurysm hemodynamics: technique and sensitivity. IEEE Trans Med Imag 2005, 24(4):457–467. 10.1109/TMI.2005.844159

    Article  Google Scholar 

  32. 32.

    Pébay PP, Baker TJ: Analysis of triangle quality measures. Math Comp 2003, 72(244):1817–1839. 10.1090/S0025-5718-03-01485-6

    MathSciNet  Article  Google Scholar 

  33. 33.

    Bloomental J: An implicit surface polygonizer. Graphic gemsIV 2004, 324–349.

    Google Scholar 

  34. 34.

    Toledo S, Chen D, Rotkin V: TAUCS: a library of sparse linearsolvers. 2010. []

    Google Scholar 

  35. 35.

    Wu X, Luboz V, Krissian K, Cotin S, Dawson S: Segmentation and reconstruction of vascular structures for 3D real-time simulation. Med Image Anal, in press.

  36. 36.

    CUDA zone 2010. []

Download references


This work was supported by National Natural Science Foundation of China (Grants No.60803108, No.30700165) and China Postdoctoral Science Foundation (Grant No. 20100470560).

Author information



Corresponding author

Correspondence to Jianhuang Wu.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

JW designed the study, performed the experiments and was responsible for the data analysis. MW and YL contributed to the implementation of the algorithm. MW, XM, and FJ contributed to the discussion and the suggestion through this project. QH contributed to the discussion and, along with JW, drafted the manuscript. All authors have read and approved the final manuscript.

Authors’ original submitted files for images

Rights and permissions


Reprints and Permissions

About this article

Cite this article

Wu, J., Wei, M., Li, Y. et al. Scale-adaptive surface modeling of vascular structures. BioMed Eng OnLine 9, 75 (2010).

Download citation


  • Computational Fluid Dynamic Simulation
  • Subdivision Surface
  • Vessel Surface
  • Surface Size
  • Marching Cube