Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow

doi:10.1016/j.cma.2007.02.009

Computer Methods in Applied Mechanics and Engineering

Volume 196, Issues 29–30, 15 May 2007, Pages 2943-2959

https://doi.org/10.1016/j.cma.2007.02.009 Get rights and content

Abstract

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, luminal 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 luminal 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.

Introduction

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 [49], [6]. 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 [32], and expanded on in [20]. 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 [7]. 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 [41]. A more advanced treatment of the subject is given in Piegl and Tiller [38]. Other geometric modeling techniques that have potential as a basis for isogeometric analysis include: A-patches [4], T-splines [44], and subdivision [14]. These warrant further investigation.

Fig. 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:

(1)
Preprocessing – in scanned computed tomography (CT) or magnetic resonance imaging (MRI) data, the intensity contrast may not be clear enough, noise exists, and sometimes the blood vessel boundary is blurred. Therefore, we use image processing techniques to improve the quality of CT/MRI data, such as contrast enhancement, filtering, classification, and segmentation.
(2)
Path extraction – the goal is to find arterial paths. Vascular surface models can be constructed from the preprocessed imaging data via isocontouring. The skeleton is then extracted from the surface model using Voronoi and Delaunay diagrams. This skeletonization scheme is suitable for noisy input and creates one-dimensional clean skeletons for blood vessels.
(3)
Control mesh construction – a skeleton-based sweeping method is developed to construct hexahedral NURBS control meshes by sweeping a templated quad mesh of a circle along the arterial path. Templates for various branching configurations are presented which decompose the geometry into mapped meshable patches using the extracted skeleton. Each patch can be meshed using one-to-one sweeping techniques. Some nodes in the control mesh lie on the surface, and some do not. We project nodes lying on the surface to the vascular surface. The blood vessel wall can be built by radially extending the surface outside 10–15% of the distance to the center line (see, e.g. [26]). It should be noted that the arterial wall is comprised of multiple layers, and may have significant regional variations in the thickness. Estimating the wall thickness correctly presents a significant challenge to the state-of-the-art imaging technology. In this work we make use of the MRI and CT data that contains no information on the wall thickness, and, as a result, we select the values that are consistent with what we found in the literature as well as private communication with vascular surgeons.
(4)
NURBS construction and isogeometric analysis – after generating hexahedral control meshes, we construct solid NURBS geometric models and employ isogeometric analysis to simulate blood flow. Piecewise linear hex meshes can also be obtained. Three numerical examples, coronary, thoracic and abdominal arteries, are presented.

The remainder of this paper is organized as follows: Section 2 reviews related previous work. Section 3 describes the meshing pipeline and preprocessing. Section 4 discusses a skeletonization scheme. Section 5 talks about solid NURBS construction and isogeometric analysis. Section 6 explains the skeleton-based sweeping method and decomposition templates. Section 7 presents three numerical examples. Section 8 draws conclusions and outlines planned future work.

Section snippets

Sweeping method

Sweeping, or $2 \frac{1}{2}$ -D 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 [9], which are swept through the volume using linking surfaces as a guide [16].

However, few geometries satisfy the topological constraints required by one-to-one sweeping. In the CUBIT project [1] at Sandia National Labs, a lot of research has been

Meshing pipeline and preprocessing

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

Skeletonization

The subsequent meshing process relies on generating a skeleton for the object to be meshed. Extracting the skeleton of a three dimensional object is a research problem that has drawn much attention for its wide applicability in graphics, solid modeling and in other diverse areas of science and engineering. Some of the techniques include, to name a few, topological thinning [11], distance field based methods [8], [12], [29], [62], potential field based methods [18], thinning via medial geodesic

Solid NURBS construction and isogeometric analysis

In a NURBS-based isogeometric analysis a physical domain in $R^{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)³ $Ω = {x = (x, y, z) \in R^{3} | x = F (ξ, η, ζ), 0 < ξ, η, ζ < 1},$ where $F (ξ, η, ζ) = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \sum_{k = 1}^{l} R_{i, j, k}^{p, q, r} (ξ, η, ζ) C_{i, j, k},$ $R_{i, j, k}^{p, q, r} = \frac{N_{i, p} (ξ) M_{j, q} (η) L_{k, r} (ζ) w_{i, j, k}}{\sum_{\hat{i} = 1}^{n} \sum_{\hat{j} = 1}^{m} \sum_{\hat{k} = 1}^{l} N_{\hat{i}, p} (ξ) M_{\hat{j}, q} (η) L_{\hat{k}, r} (ζ) w_{\hat{i}, \hat{j}, \hat{k}}} .$ In the above, $R_{i, j, k}^{p, q, r} (ξ, η, ζ)$ ’s are the rational basis functions, and $C_{i, j, k}$ ’ $s \in R^{3}$ are the control points. In the definition of the

The skeleton-based sweeping method

Blood vessels are tubular objects, therefore we choose the sweeping method to construct hexahedral control meshes for NURBS-based isogeometric analysis.

Numerical examples

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

Conclusions and future work

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 hexahedral solid 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 use the sweeping method to construct hexahedral control meshes,

Acknowledgements

An early version of this paper appeared in 15th International Meshing Roundtable conference [61]. 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

References (62)

F. Cirak et al.
Integrated modeling, finite-element analysis, and engineering design for thin-shell structures using subdivision
Comput.-Aided Design
(2002)
J.A. Cottrell et al.
Isogeometric analysis of structural vibrations
Comput. Methods Appl. Mech. Engrg.
(2006)
T.J.R. Hughes et al.
Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement
Comput. Methods Appl. Mech. Engrg.
(2005)
R.L. Ogniewicz et al.
Hierachic Voronoi skeletons
Pattern Recognit.
(1995)
O. Sahni et al.
Efficient anisotropic adaptive discretization of the cardiovascular system
Comput. Methods Appl. Mech. Engrg.
(2006)
C.A. Taylor et al.
Finite element modeling of blood flow in arteries
Comput. Methods Appl. Mech. Engrg.
(1998)
Y. Zhang et al.
Adaptive and quality quadrilateral/hexahedral meshing from volumetric data
Comput. Methods Appl. Mech. Engrg. (CMAME)
(2006)
Cubit mesh generation toolkit, web site:...
C.W. Anderson, S. Crawford-Hines, Fast generation of nurbs surfaces from polygonal mesh models of human anatomy, in:...
C. Armstrong, D. Robinson, R. McKeag, T. Li, S. Bridgett, R. Donaghy, C. McGleenan, Medials for meshing and more, in:...

C. Bajaj et al.

Modeling with cubic A-patches

ACM Transactions on Graphics

(1995)

C. Bajaj, Q. Wu, G. Xu, Level Set Based Volumetric Anisotropic Diffusion, in: ICES Technical Report 301, the University...

Y. Bazilevs et al.

Isogeometric fluid–structure interaction analysis with applications to arterial blood flow

Comput. Mech.

(2006)

Y. Bazilevs et al.

Isogeometric analysis: approximation, stability and error estimates for h-refined meshes

Math. Models Methods Appl. Sci.

(2006)

I. Bitter et al.

Penalized distance volumetric skeleton algorithm

IEEE TVCG

(2001)

T. Blacker

A new approach to automated quadrilateral mesh generation

Int. J. Numer. Meth. Engrg.

(1991)

T. Blacker, The Cooper Tool, in: 5th Int. Meshing Roundtable, 1996, pp....

G. Borgefors, I. Nystrom, G.D. Baja, Computing skeletons in three dimensions, Pattern Recognit. 32 (7)...

S. Bouix et al.

Flux driven fly throughs

IEEE Conf. Comput. Vision Pattern Recognit.

(2003)

CGAL Consortium, CGAL: Computational Geometry Algorithms Library,...

Cocone, Tight Cocone Software for surface reconstruction and medial axis approximation,...

W.A. Cook et al.

Mapping methods for generating three-dimensional meshes

Comput. Mech. Engrg.

(1982)

N. Cornea et al.

Curve skeleton applications

IEEE Visualiz.

(2005)

N. Cornea et al.

Computing hierarchical curve skeletons of 3D objects

Visual Comput.

(2005)

L. Costa, Multidimensional scale space shape analysis, in: IWSNHC3DI, 1999, pp....

M. de Berg et al.

Comput. Geometry: Algorithms Appl.

(1997)

T.K. Dey, J. Giesen, S. Goswami, Shape segmentation and matching with flow discretization, in: F. Dehne, J.-R. Sack, M....

T.K. Dey, S. Goswami, Tight cocone: a water-tight surface reconstructor, in: Proc. 8th ACM Sympos. Solid Modeling...

T.K. Dey, J. Sun, Defining and computing curve-skeletons with medial geodesic functions, in: Sympos. Geom. Proces.,...

T.K. Dey et al.

Approximating the medial axis from the Voronoi diagram with convergence guarantee

Algorithmica

(2004)

A. Figueroa et al.

A coupled momentum method for modeling blood flow in three-dimensional deformable arteries

Comput. Methods Appl. Mech. Engrg.

(2006)

Cited by (369)

Echocardiogram-based ventricular isogeometric cardiac analysis using multi-patch fitted NURBS
2024, Computer Methods in Applied Mechanics and Engineering
Monitoring the cardiac function often relies on non-invasive vital measurements, electrocardiograms, and low-resolution echocardiograms. The monitoring data can be augmented with numerical modeling results to support treatment-risk assessment. Often, medical images are not suitable for high-fidelity modeling due to their spatial sparsity and low resolution. In this work, we present a workflow that converts a limited number of two-dimensional (2D) echocardiogram images into analysis-suitable left-ventricle geometries for isogeometric analysis (IGA). The image data is fitted by reshaping a multi-patch NURBS template. This fitting procedure results in a smooth interpolation with a limited number of degrees of freedom. In addition, a coupled 3D-0D cardiac model is extended with a pre-stressing method to perform ventricular-cardiac-mechanic analyses based on in vivo images. The workflow is benchmarked using population-averaged imaging data, which demonstrates that both the shape of the left ventricle and its mechanical response are well predicted in sparse-data scenarios. The case of a clinical patient is considered to demonstrate the suitability of the workflow in a real-data setting.
Simultaneous Boundary and Interior Parameterization of Planar Domains Via Deep Learning
2024, CAD Computer Aided Design
In isogeometric analysis (IGA), domain parameterization plays a crucial role as it significantly impacts the results of subsequent numerical analyses. Previous literature has often treated domain parameterization as two separate processes: boundary correspondence and interior parameterization. However, this division tends to decrease the quality of the final parameterization. To address this issue, this article proposes a learning-based method that simultaneously parameterizes the boundary and interior of a simply connected planar computational domain using a deep neural network. In our approach, we train a neural network to represent the inverse parameterization, which maps the computational domain to the parametric domain, based on the given boundary curve of the planar computational domain. The key component of our neural network is a Multilayer Perceptron (MLP), which takes the coordinates of a point in the computational domain as input and outputs the corresponding parametric values. The loss function of our method comprises three terms. First, it considers the area distortion and angular distortion of the parameterization, which are determined by the Jacobian matrix of the parameterization. Second, it incorporates the boundary difference, quantified by the Hausdorff distance between the boundary of the parametric domain and the image of the inverse mapping of the computational boundary. By optimizing this loss function, we achieve a robust parameterization. Once the network is trained, we employ a tensor-product B-spline function to fit the mapping from the parametric domain to the computational domain. To evaluate the effectiveness of our method, we conduct experiments on the MPEG-7 dataset. The experimental results demonstrate that our approach yields parameterization results with lower distortion and higher bijectivity ratio compared to state-of-the-art approaches.
Subdivision-based IGA Coupled EIEQ Method for the Cahn–Hilliard Phase-Field Model of Homopolymer Blends on Complex Surfaces
2023, CAD Computer Aided Design
In this paper, we construct an IGA-EIEQ coupling scheme to solve the phase-field model of homopolymer blends on complex subdivision surfaces, in which the total free energy contains a gradient entropy with a concentration-dependent de-Gennes type coefficient and a non-linear logarithmic Flory–Huggins type potential. Based on the EIEQ method, we develop a fully-discrete numerical scheme with the superior properties of linearity, unconditional energy stability, and second-order time accuracy. All we need to do with this fourth-order system is to solve some constant-coefficient elliptic equations by applying a new nonlocal splitting techniqueWe then provide detailed proof of the unconditional energy stability and the practical implementation process. Subdivision approaches show a robust and elegant description of the models with arbitrary topology. Subdivision basis functions serve to define the geometry of the models and represent the numerical solutions. Subdivision-based IGA approach provides us with a good candidate for solving the phase-field model on complex surfaces. We successfully demonstrate the unity of employing subdivision basis functions to describe the geometry and simulate the dynamical behaviors of the phase-field models on surfaces with arbitrary topology. This coupling strategy combining the subdivision-based IGA method and the EIEQ method could be extended to a lot of gradient flow models with complex nonlinearities on complex surfaces.
Modeling and hexahedral meshing of cerebral arterial networks from centerlines
2023, Medical Image Analysis
Computational fluid dynamics (CFD) simulation provides valuable information on blood flow from the vascular geometry. However, it requires extracting precise models of arteries from low-resolution medical images, which remains challenging. Centerline-based representation is widely used to model large vascular networks with small vessels, as it encodes both the geometric and topological information and facilitates manual editing. In this work, we propose an automatic method to generate a structured hexahedral mesh suitable for CFD directly from centerlines.
We addressed both the modeling and meshing tasks. We proposed a vessel model based on penalized splines to overcome the limitations inherent to the centerline representation, such as noise and sparsity. The bifurcations are reconstructed using a parametric model based on the anatomy that we extended to planar n-furcations. Finally, we developed a method to produce a volume mesh with structured, hexahedral, and flow-oriented cells from the proposed vascular network model.
The proposed method offers better robustness to the common defects of centerlines and increases the mesh quality compared to state-of-the-art methods. As it relies on centerlines alone, it can be applied to edit the vascular model effortlessly to study the impact of vascular geometry and topology on hemodynamics. We demonstrate the efficiency of our method by entirely meshing a dataset of 60 cerebral vascular networks. 92% of the vessels and 83% of the bifurcations were meshed without defects needing manual intervention, despite the challenging aspect of the input data. The source code is released publicly.
Stability of non-linear flow in heterogeneous porous media simulations using higher order and continuity basis functions
2023, Journal of Computational and Applied Mathematics
Using isogeometric analysis, we present a highly optimized FORTRAN code for simulations of three-dimensional non-linear flow problem in heterogeneous media. This includes deriving the weak formulations, discretizing with the forward Euler scheme, and implementing computation of the right-hand side. Some properties of exact solutions of heat transport and non-linear flow in heterogeneous media problems have been derived to validate the numerical solution. Stability analysis of the employed explicit time-stepping scheme has been performed. First, the properties of the discrete form of the heat transfer problem have been systematically investigated, leading to the conjecture concerning the sufficient condition of stability of the scheme, which numerical computations have subsequently confirmed. A similar analysis has been carried out for the non-linear flow in heterogeneous media problem. For this problem, the bound on the time step was conjectured to decrease exponentially, which has been verified experimentally. The simulation’s behavior close to the stability limit has been investigated and found to be considerably more complex than for the heat transport problem. Detailed experimental analysis of the time complexity of the particular parts of the simulator algorithm has been described. In addition, dependence on both mesh size and order of continuity of B-spline basis has been investigated.
A connection element method: Both a new computational method and a physical data-driven framework—-Take subsurface two-phase flow as an example
2023, Engineering Analysis with Boundary Elements
In traditional reservoir simulation methods, the computational domain is often discretized into a grid system (FDM, FVM, FEM, etc.) or a mesh reduction system (point source solution, BEM, meshless methods, etc.). This paper develops a connection-based method, which uses connection elements to discretize the computational domain from the flow perspective (named as connection element method, CEM). Based on a generalized finite difference approximation of the pressure diffusion term, each connection element can be characterized by two parameters, namely connection transmissibility and connection volume. The connection transmissibility measures the diffusion capacity of the connection element, and the connection volume measures the volume of the connection element as a geometric entity. Then, the flow equation can be solved directly on the connection elements, by either a fully implicit scheme or a sequential coupled scheme. The fully implicit scheme solves pressures and saturations simultaneously with high accuracy, but the computational cost is very high when the number of connections is large. The sequential coupled scheme calculates the pressure at the nodes of connection elements at the coarse-scale, then solves the fine-scale saturation along the connection elements. In this way, high computational efficiency can be achieved with accuracy. Furthermore, CEM can efficiently and intuitively obtain the possible flow paths between connection elements and reveal the interaction between the nodes by the path-tracking algorithm in graph theory. In addition, CEM with a sequential coupled scheme can be regarded as a physical network framework with small degrees of freedom, as well as a good data-driven framework. In the end, two numerical examples and one field case are presented to demonstrate the superiority of the developed CEM.

View all citing articles on Scopus

^☆: Visit http://www.ices.utexas.edu/~jessica/paper/vascularmodel.

View full text

Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow☆

Abstract

Introduction

Section snippets

Sweeping method

Meshing pipeline and preprocessing

Skeletonization

Solid NURBS construction and isogeometric analysis

The skeleton-based sweeping method

Numerical examples

Conclusions and future work

Acknowledgements

Comput.-Aided Design

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Pattern Recognit.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg. (CMAME)

Modeling with cubic A-patches

ACM Transactions on Graphics

Isogeometric fluid–structure interaction analysis with applications to arterial blood flow

Comput. Mech.

Isogeometric analysis: approximation, stability and error estimates for h-refined meshes

Math. Models Methods Appl. Sci.

Penalized distance volumetric skeleton algorithm

IEEE TVCG

A new approach to automated quadrilateral mesh generation

Int. J. Numer. Meth. Engrg.

Flux driven fly throughs

IEEE Conf. Comput. Vision Pattern Recognit.

Mapping methods for generating three-dimensional meshes

Comput. Mech. Engrg.

Curve skeleton applications

IEEE Visualiz.

Computing hierarchical curve skeletons of 3D objects

Visual Comput.

Comput. Geometry: Algorithms Appl.

Approximating the medial axis from the Voronoi diagram with convergence guarantee

Algorithmica

A coupled momentum method for modeling blood flow in three-dimensional deformable arteries

Comput. Methods Appl. Mech. Engrg.