|Home | About | Journals | Submit | Contact Us | Français|
We describe an approach to construct hexahedral solid NURBS (Non-Uniform Rational B-Splines) meshes for patient-specific vascular geometric models from imaging data for use in isogeometric analysis. First, image processing techniques, such as contrast enhancement, filtering, classification, and segmentation, are used to improve the quality of the input imaging data. Then, lumenal surfaces are extracted by isocontouring the preprocessed data, followed by the extraction of vascular skeleton via Voronoi and Delaunay diagrams. Next, the skeleton-based sweeping method is used to construct hexahedral control meshes. Templates are designed for various branching configurations to decompose the geometry into mapped meshable patches. Each patch is then meshed using one-to-one sweeping techniques, and boundary vertices are projected to the lumenal surface. Finally, hexahedral solid NURBS are constructed and used in isogeometric analysis of blood flow. Piecewise linear hexahedral meshes can also be obtained using this approach. Examples of patient-specific arterial models are presented.
Recently, patient-specific modeling was proposed as a new paradigm in simulation-based medical planning. Physicians, using computational tools, construct and evaluate combined anatomical/physiological models to predict the outcome of alternative treatment plans for an individual patient. A comprehensive framework has been developed to enable the conduct of computational vascular research [1, 2]. Blood flow simulations provide physicians with physical data to help them devise treatment plans.
Isogeometric analysis is a new computational technique that improves on and generalizes the standard finite element method. It was first introduced in , and expanded on in . In an effort to instantiate the concept of isogeometric analysis, an analysis framework based on NURBS was built. Mathematical theory of this NURBS-based approach was put forth in . NURBS is not the only possible basis for isogeometric analysis but it is certainly the most highly developed and widely utilized. For an introductory text on NURBS, see Rogers . A more advanced treatment of the subject is given in Piegl and Tiller . Other geometric modeling techniques that have potential as a basis for isogeometric analysis include: A-patches , T-splines , and subdivision . These warrant further investigation.
Figure 1 shows one such model, obtained from patient-specific imaging data. We have designed a set of procedures which allows us to create solid NURBS vascular models directly from patient-specific data. We have named this process the vascular modeling pipeline, which can be divided into four main steps:
The remainder of this paper is organized as follows: Section 2 reviews related previous work. Section 3 describes the meshing pipeline and preprocessing for our geometric modeling approach. Section 4 talks about solid NURBS construction and isogeometric analysis. Section 5 explains the skeleton-based sweeping method and decomposition templates for various branching configurations. Section 6 presents three numerical examples. Section 7 draws conclusions and outlines planned future work.
Sweeping, or Inline meshing, is one of the most robust techniques to generate semi-structured hexahedral meshes. One-to-one sweeping requires that the source and target surfaces have similar topology. The source surface is meshed with quadrilaterals , which are swept through the volume using linking surfaces as a guide .
However, few geometries satisfy the topological constraints required by one-to-one sweeping. In the CUBIT project  at Sandia National Labs, a lot of research has been done to automatically recognize features and decompose geometry into mapped meshable areas or volumes. Various many-to-one and many-to-many sweeping methods have been developed [14, 15, 16, 17]. Care should also be taken in locating internal nodes during the sweeping process [18, 19].
Medial axis is the locus of points that are minimally equidistant from at least two points on the geometry's boundary. The medial axis transform provides an alternative representation of geometric models that has many useful properties for analysis modeling . Applications include decomposition of general solids into subregions for mapped meshing, identification of slender regions for dimension reduction and recognition of small features for suppression. The medial surface subdivision technique  decomposes the volume into map-meshable regions, which are then filled with hex elements using templates.
Medial axis has been used to construct hexahedral meshes for CAD objects. The skeleton-based modeling methods were developed for solid models . Quadros et al. used a skeleton technique to control finite element mesh size . Besides other unstructured mesh generation methods [24, 25, 26], a skeleton-based subdivision method has also been used in biomedical applications, such as a below-knee residual limb and external prosthetic socket , and bifurcation geometry in vascular flow simulation . However, trifurcations and more complex branchings also exist in the human artery tree. Therefore, decomposition templates for arbitrary branching configurations are desirable and are constructed in this paper.
As the most highly developed and widely utilized technique, NURBS [6, 7, 29] has evolved into an essential tool for a semi-analytical representation of geometric entities. Sometimes NURBS solid models are taken as input for finite element mesh generation . Anderson et al. proposed a fast generation of NURBS surfaces from polygonal mesh models of human anatomy . An enhanced algorithm was developed for NURBS evaluation and utilization in grid generation . In isogeometric analysis , NURBS basis functions are used to construct the exact geometry, as well as the corresponding solution space.
The input images are often of poor quality which makes it difficult to generate quality meshes for regions of interest. To circumvent this problem we pass the raw imaging data through a preprocessing pipeline where the image quality is improved by enhancing the contrast, filtering noise, classifying, and segmenting regions of various materials. The surface model is then extracted from the processed imaging data, and the vessel path is obtained after skeletonizing the volume bounded by the surface. First we modify the geometry by removing unnecessary components, then extract the skeleton. The generated path can also be edited according to simulations, e.g., adding a path for a Left Ventricle Assist Device (LVAD) in the thoracic aorta model (Figure 12). A skeleton-based sweeping method is then used to generate hexahedral control meshes for solid NURBS construction and isogeometric analysis. Figure 2 shows the meshing pipeline. The preprocessing step of our skeleton-based meshing approach is described below, including image processing, isocontouring and geometry editing, and path extraction.
We choose a fast localized method for image contrast enhancement . The basic idea is to design an adaptive transfer function for each individual voxel based on the intensities in a suitable local neighborhood. A bilateral pre-filtering coupled with an evolution driven anisotropic geometric diffusion PDE (partial differential equation)  is utilized to remove noise. Sometimes we need to classify the voxels into several groups, each of which corresponds to a different type of material. We choose an approach which relies on identification of the contours by membership of seed points which are located by the gradient vector diffusion . A variant of the fast marching method is adopted  to segment the imaging data to find the clear boundary of each voxel group belonging to a certain category.
There are two main isocontouring methods from imaging data: Primal Contouring (or Marching Cubes ) and Dual Contouring . In this application we choose Dual Contouring to extract the isosurface, because it tends to generate meshes with better aspect ratios. We then modify the model to suit our particular application. This can be done in various ways, for example, by removing unnecessary components, adding necessary components which are not constructed from imaging data, denoising the surface, etc. After getting the vessel path, we can edit it according to simulation requirements. For example, we can add a path for the left ventricle assist device (LVAD) in the thoracic aorta model (Figure 12).
The vertex set of the extracted and possibly repaired geometry is then used to create an interior path lying in the middle of the blood vessels. We define a squared distance function which assigns to any point x 3, the minimum square distance to the vertex set. We further compute the index 1 and index 2 saddle points of this distance function and compute the unstable manifold of these two types of critical points. The identification of the critical points along with their indices and the computation of the unstable manifold are done efficiently via the Voronoi and its dual Delaunay diagram of the point set. The details of this method can be found in . We adopt this method of path generation because it has several advantages which are useful for the patient specific modeling of blood vessels. One advantage is that it can handle noisy input gracefully. Often the noise present in the data is not fully eliminated after the preprocessing stage. In the path generation step we employ another stage of filtering which helps to construct a clean skeletal path for the extracted geometry. Secondly, the extracted geometry may have flat regions where it is not straight forward to obtain a linear path. Fortunately our starring scheme, as described in , eliminates these spurious features and create the one-dimensional path. The results of this path generation step are shown on various datasets (Figures 1, ,11,11, ,1212).
In a NURBS-based isogeometric analysis a physical domain in 3 is defined as a union of patches. A patch, denoted by Ω, is an image under a NURBS mapping of a parametric domain (0,1)3
In the above, (ξ, η, ζ,)'s are the rational basis functions, and Ci,j,k's 3 are the control points. In the definition of the rational basis, Ni,p(ξ)'s, Mj,q(η)'s, and Lk,r(ζ)'s, are the univariate B-spline basis functions of polynomial degree p, q, and r; wi, j,k's, strictly positive, are the weights.
In isogeometric analysis the geometry generation step involves construction of a control mesh, which is a piecewise multi-linear interpolation of control points, and the corresponding rational basis functions. The initial mesh encapsulates the ‘exact geometry’ and, in fact, defines it parametrically.
For the purposes of analysis, the isoparametric concept is invoked (see Hughes ). The basis for the solution space in the physical domain is defined through a push forward of the rational basis functions defined in (2) (see  for details). Coefficients of the basis functions, defining the solution fields in question (e.g., displacement, velocity, etc.), are called control variables.
As a consequence of the parametric definition of the ‘exact’ geometry at the coarsest level of discretization, mesh refinement can be performed automatically without further communication with the original description. This is an enormous benefit. There are NURBS analogues of finite element h- and p-refinement, and there is also a variant of p-refinement, which is termed k-refinement, in which the continuity of functions is systematically increased. This seems to have no analogue in traditional finite element analysis but is a feature shared by some meshless methods. For the details of the refinement algorithms see .
The isogeometric approach is fundamentally higher-order. For example, in order to represent circles, cylinders and spheres, rational polynomials of at least quadratic order are necessary. The generation of refined NURBS bases of all orders is facilitated by simple recursion relationships. The versatility and power of recursive NURBS basis representations are truly remarkable. Equation systems generated by NURBS tend to be more homogeneous than those generated by higher-order finite elements and this may have some benefit in equation solving strategies. NURBS satisfy a ‘variation diminishing’ property. For example, they give monotone fits to discontinuous control data and become smoother as order is increased, unlike Lagrange interpolation polynomials which oscillate more violently as order is increased. NURBS of all orders are non-negative pointwise. This means that every entry of the NURBS mass matrix is non-negative. These properties are not attained in finite element analysis. On the other hand, NURBS are not interpolatory. They are fit to nets of control points and control variables. This aspect is less transparent to deal with than the corresponding finite element concepts of interpolated nodal points and nodal variables but somewhat similar to the situation for meshless methods. There are many robust algorithms to create very complex geometries with NURBS.
Blood vessels are tubular objects, therefore we choose the sweeping method to construct hexahedral control meshes for NURBS-based isogeometric analysis.
In the sweeping method, a templated quadrilateral mesh of a circle is projected onto each cross-section of the tube, then corresponding vertices in adjacent cross-sections are connected to form a hexahedral mesh. A hexahedral NURBS control mesh should satisfy the following four requirements:
We choose to parameterize the template cross-section as follows. One parametric direction is associated with a closed circular curve, while another parametric direction is associated with a radial coordinate. Rational quadratic basis is used to define the circular curve with a control polygon given by the linear interpolation of the green and blue points shown in Figure 3. For the template shown in Figure 3a, the control polygon is a square consisting of 8 control nodes, while in Figure 3b, it is an octagon. Note that the circular cross-section is unchanged geometrically and parametrically as more control points are chosen for its representation. The green control points lie on the circle, while the blue control points do not. This is due to the fact that the rational basis is interpolatory at the green points and is not interpolatory at the blue points. Also note that each interpolatory control point has two neighboring non-interpolatory points that are collinear with it. This construction guarantees the resultant circular curve to be G1-continuous. Later, when we discuss data fitting, it is only the interpolatory points that get projected onto the true surface. The non-interpolatory points are adjusted to preserve the collinearity in order to obtain a G1-continuous cross-section.
In the process of sweeping, we translate the cross-section template to the selected locations on the path, and rotate it to make its normal vector pointing in the direction tangent to the path as shown in Figure 4. This gives the third parametric direction for the solid NURBS representation. The hexahedral control mesh is constructed by connecting the corresponding control nodes in adjacent cross-sections. Piecewise linear hexahedral meshes can also be generated at the same time by projecting all boundary vertices to the vessel surface, or by interpolating the elements of the solid NURBS geometry.
One-to-one sweeping requires that the source and target surfaces have similar topology. Generally, arterial models do not satisfy this requirement, therefore we need to decompose arterial networks into mapped meshable regions. In this section, we will discuss various decomposition templates for different branching configurations. An n-branching is formed when n branches join together, where n ≥ 3. When n = 3, it is a bifurcation; when n = 4, it is a trifurcation; when n > 4, we call this situation higher order branching.
In the human vascular system, most branchings are bifurcations. However, trifurcations or higher order branchings also exist. For example, there are several trifurcations in the coronary arteries (Figure 11) and the abdominal aorta (Figure 1). In the following, we will discuss decomposition templates for all possible branching configurations.
For every intersection, a so-called master arterial branch is chosen. Typically, it is an artery with the largest diameter. Suppose the master branch consists of two sub-branches (Branch 1 and Branch 2), and the slave branch is Branch 3, as shown in Figure 5a. The axes of Branch 1, 2 and 3 are Axis 1, 2 and 3 respectively (Axis 1 and Axis 2 may not be collinear). There is one basic case, shown in Figure 5, and all bifurcations can be decomposed into three map-meshable regions by a variant of this basic template.
Figure 5 shows the path, the constructed hexahedral control mesh, the solid NURBS mesh, and the piecewise linear hexahedral meshes of the bifurcation template. The bifurcation geometry is decomposed into three patches: the master branch contains two patches (red and green), and the slave branch has one patch (yellow). Here we choose Level-1-template (Figure 3a) for each cross-section, as the master and slave branches have similar diameters. The bifurcation template also works for finer cross-sections.
When the master branch and the slave branch have different diameters, the control nodes of some cross-sections are distributed unevenly in order to generate better intersection regions. Figure 6 shows two control meshes and their corresponding solid NURBS meshes. The master branch control polygon is deformed from a square to a trapezoid so as to accommodate a slave branch with a smaller diameter. Note that the NURBS basis changes accordingly so as to preserve the circular cross-section, and the quality of the intersection geometry is improved as can be seen in Figure 6 and Figure 7, where the axes of the master and the slave branches are non-orthogonal, or non-coplanar. Although deforming the control polygon of the master branch gives better results as compared to the non-deformed case, for the intersection of branches with high diameter ratios we advocate the use of a finer template for the master branch, such as a Level-2-template or a Level-3-template.
Trifurcation has one master branch and two slave branches. According to the position of slave branches relative to the master branch, we classify all possible trifurcations to fall into five irreducible cases. All other trifurcations can be decomposed into map-meshable regions by extending the five basic decomposition templates.
Case 1: The two slave branches are distributed along the peripheral direction of the master branch, and they are in opposite relative to the master branch (the angle between them is around 180°). The same cross-section template can be used for the master and slave branches.
Case 2: The two slave branches are distributed along the peripheral direction, and the angle between them is arbitrary. Finer cross-section template is chosen for the master branch.
Case 3: The two slave branches are distributed along the axial direction of the master branch, and they intersect with each other.
Case 4: The two slave branches are distributed along the axial direction of the master branch, and they do not intersect with each other. This situation degenerates into two bifurcations.
Case 5: The two slave branches do not intersect with the master branch at the same point, but they intersect with each other. In this situation, two bifurcations merge into one trifurcation.
If a Level-1-template is selected as the cross-section of the master branch, then there are at most two slave branches along the peripheral direction as shown in Figure 8 (Case 1). If the two slave branches are not opposite relative to the master branch, or the two slave branches have different diameters from the master branch, then a Level-1-template is not suitable, and we need to choose finer cross-section templates, such as a Level-2-template or a Level-3-template. Similarly, if a Level-2-template is selected as the cross-section of the master branch, then there can be at most four slave branches along the peripheral direction. A Level-3-template allows at most eight slave branches along the peripheral direction. In Case 2 of Figure 8, the two slave branches are distributed along the peripheral direction and they are not opposite, therefore we choose the finer cross-section template for the master branch (Level-2-template), while the slave branch may have coarser cross-sections (Level-1-template).
Case 3 and Case 4 have the same path, but Case 4 degenerates into two bifurcations because its two slave branches do not intersect with each other even though their axes intersect. There is another special situation (Case 5) where two slave branches do not intersect with the master branch at the same intersection point in the skeleton, but the two intersection points are very close and the two slave branches intersect with each other. This situation contains two bifurcations in the skeleton, but it should be considered as one trifurcation. Therefore, when we choose branching configurations, both the path and the vessel size should be considered.
Remark: In n-branching, n should be decided not only by the path, but also by the diameter of each slave branch. In other words, if neighboring slave branches intersect with each other, then it is n-branching. Otherwise, it degenerates into several m-branchings, where m < n. On the other hand, several m-branchings may merge into one n-branching if its slave branchings intersect with each other.
Here we discuss three basic templates for n-branching when n > 4. Relative to the master branch, there are only two directions to arrange slave branches, the peripheral and axial directions of the master branch. All other n-branching configurations can be obtained by combining the three basic ones.
Case 1: There are three or more slave branches distributed along the peripheral direction of the master branch. Figure 9 shows one example of four slave branches along the peripheral direction. Level-2-template is selected for the master branch. If there are more than four slave branches, the master branch needs to have a finer cross-section. The cross-section template of slave branches can be coarser.
Case 2: There are three or more slave branches distributed along the axial direction of the master branch. Neighboring slave branches intersect with each other.
Case 3: There are three or more slave branches distributed along the axial direction of the master branch. Slave branches do not intersect with each other. n-branching degenerates into several m-branchings (m < n).
Several lower order branchings may merge into a higher order one, for example, one bifurcation and one trifurcation can merge into a 5-branching. Case 1, Case 2 and Case 3 can be combined together to form more complex configurations. Figure 9 shows one example of 7-branching. It has four slave branches along the peripheral direction and two slave branches along the axis of the master branch.
After the sweeping step, each circular cross-section needs to be projected onto the vessel surface as shown in Figure 10. First the interpolatory control points (green points) are moved in the radial direction to the true surface. Then, the non-interpolatory points (blue points) are placed at the intersection of the lines tangent to the true surface passing through the two neighboring interpolatory points.
There are situations when the tangent lines do not intersect inside the fan region defined in Figure 10b, or do not even intersect with each other when they are parallel. This may occur when the cross-section template is not sufficiently fine to capture features of the true surface, or when the true surface is noisy. This situation will also result in an overlap in the geometry. In order to avoid overlap, we force the non-interpolatory point to stay inside the fan region (the sector between two radial rays) by placing it at the midpoint (indicated by the magenta color) of the segment connecting the two interpolatory points. Finally, the location of the interpolatory points is changed so as to preserve G1-continuity of the surface.
After projecting each cross-section to the vessel surface, we construct hexahedral control meshes and generate solid NURBS for patient-specific vascular models. The geometric error can be reduced by choosing a finer template for each cross-section.
The proposed construction of a NURBS solid mesh for isogeometric analysis has the following implications:
In this section we present applications of the meshing pipeline to three patient-specific vascular models: a model of a portion of the coronary tree, a model of the thoracic aorta, and a model of the abdominal aorta. Isogeometric analysis is then used to compute blood flow in the models. In all cases, time-dependent, viscous, incompressible Navier-Stokes equations were used as the blood model. The fluid density and dynamic viscosity were chosen to be representative of blood flow. The first example makes use of the Casson model for the dynamic viscosity while in other examples viscosity was set to a constant value. All models are subjected to a time-periodic inflow boundary condition, which simulates the input from a beating heart. The arterial wall is assumed rigid in the first example. Examples two and three present fluid-structure interaction calculations in which the wall is assumed to be elastic (see Bazilevs et al.  for the details of the mathematical formulation). The rigid wall simulation was performed on a single processor, while the elastic wall simulations were done in parallel.
Data for this model was obtained from CT Angiography imaging data of a healthy male, over 55 years of age. Large motions of the heart, as it supplies blood to the circulatory system, decreases the quality of the imaging data, and makes construction of patient-specific coronary models a challenging task. Nevertheless, we managed to extract a portion of the coronary tree for the purposes of creating an analysis suitable model. Results of the isocontouring algorithm are shown in Figure 11a. Figures 11b-11d show the path, the control mesh, and the solid NURBS model of the arterial segment. The model was used to study drug delivery processes in arteries. The drug concentration in the blood is modeled as a passive scalar governed by an unsteady advection-diffusion equation. Figure 11e shows the isosurface of the drug concentration at 50% colored by the blood velocity magnitude, revealing that the flow is unsteady, and has many complex features.
Data for this model was obtained from CT Angiography imaging data of a healthy male over-30 volunteer. A patient-specific model of the thoracic aorta was constructed by running through the meshing pipeline. An extra branch, representing a left ventricular assist device (LVAD), was added to the arterial model. Evaluation of LVADs, as well as other electromechanical devices used to support proper blood circulation, is of great interest to the cardiovascular community. The path, the control mesh, and the solid NURBS model are shown in Figures12a-12c. Figure 12d shows a result of the fluid-structure interaction simulation. Note that the inlet and the three smaller outlet branches were extended for the purposes of analysis.
Data for this model was obtained from 64-slice CT angiography of a healthy male over 55 years of age. Various stages of the meshing pipeline are illustrated in Figure 1a--1f.1f. Figure 1g shows a result of the fluid-structure simulation. A computational study using a truncated geometrical model of this aorta was performed in . We used 85 seconds for path extraction and 8 seconds for control mesh construction on a 64-bit dual-AMD 2GHz linux system, and 20 seconds for solid NURBS generation on an Intel 3GHz linux system.
We have developed a four-stage process to construct analysis suitable geometric models from patient-specific vascular imaging data with a goal of using them in isogeometric analysis of blood flow in arteries. We have focused on the NURBS modeling, and did not treat other geometrical modeling technologies, such as A-patches, T-splines, and subdivision. We would like to investigate these techniques in the future.
We have successfully applied our method to three patient-specific examples, which involve a model of a part of the coronary arterial tree, a thoracic aorta model, and an abdominal aorta model. As part of the future work, we would like to apply the techniques described here to modeling and analysis of the human heart.
Y. Zhang was partially supported by the J. T. Oden ICES Postdoctoral Fellowship at the Institute for Computational Engineering and Sciences. This research of Y. Zhang, S. Goswami, and C. Bajaj was supported in part by NSF grants EIA-0325550, CNS-0540033, and NIH grants P20-RR020647, R01-GM074258, R01-GM073087. This support is gratefully acknowledged. We would also like to thank Fred Nugen, Bob Moser, and Jeff Gohean for providing us with the data for the thoracic aorta model.