problem domain, i.e. a subdivision of the domain into non-overlapping subdo- ..... Nonrigid Registration Using Free-Form Deformations: Application to Breast MR.
Quadrilateral Meshes for Finite Element-Based Image Registration Marcelo Siqueira1 , Tessa Sundaram1 , Suneeta Ramaswami2, Jean Gallier1 , and James Gee1 1
University of Pennsylvania, Philadelphia, PA 19104, USA {marcelos,tessa}@seas.upenn.edu 2 Rutgers University, Camden, NJ 08102, USA
Abstract. The use of finite element (FE) analysis in biomedical research has prompted the need to construct meshes from images. Despite the availability of several tools for generating meshes for FE-based applications, most are not directly applicable to image data. Additionally, they may not produce good-quality meshes when presented with geometrically complex domains such as those in the human body. Here we investigate the construction of quadrilateral meshes from 2D images of the human brain. We describe an algorithm to derive quadrangulations from triangulations of polygonal regions, an approach to construct quadrilateral finite element meshes from segmented images using the aforementioned algorithm, and an analysis of the quality and performance of the resulting meshes when applied to a FE-based image registration.
1
Introduction
Image registration is the process of finding a spatial correspondence between two images, [1]. This correspondence is usually achieved by optimizing a function that balances the similarity between the image pair with the deformation caused by the applied registration transformation. Transformations can be rigid (such as translation or rotation) or non-rigid. The former are used to globally align the image pair, while the latter are necessary when regional distortions are required to effectively match the two images. Image registration has found many uses in biomedical research. For example, it can be employed to extract motion information about an organ (such as the heart or lung), [2], or to study the morphological changes caused by tumor growth within a region, [3]. Broit, [4], developed a method for non-rigid image registration in which one image, modeled as an elastic continuum, is warped to match the appearance of another. The variational formulation of the registration problem makes it possible to implement his technique using the finite element method (FEM), [5]. The FEM is a very powerful method for numerically solving differential equations or variational problems that arise during structural modeling in engineering and the applied sciences. One of its key features is the specification of a mesh over the problem domain, i.e. a subdivision of the domain into non-overlapping subdomains called finite elements. If the problem domain is a subset of the Euclidean
plane, triangular or quadrilateral elements are typically utilized. The accuracy of a problem’s solution and the efficiency with which it is obtained using a particular FE implementation are highly dependent on a number of mesh parameters, such as element number, size, shape, connectivity, and physical behavior, [6]. Triangular meshes have been extensively investigated by the meshing community, and their theoretical properties are now well understood. Algorithms for generating provably good triangular meshes have been proposed, [7]. However, the generation of good quadrilateral meshes is not as well understood. A few algorithms exist that generate quadrilateral meshes of bounded size, [8], and bounded maximum polygonal angle, [9], but there are no known algorithms that optimize size and shape for more complex geometric domains. Such qualityguaranteed meshes are more desirable for certain FE-based applications. The use of the FEM in medical applications has led researchers to explore the construction of meshes that attempt to account for the geometric complexity of the anatomy under study. Most of the available tools for automated mesh generation cater to the construction of meshes for structural engineering applications, [10], but these are typically optimized for much simpler shapes than those encountered in biology, [11]. This paper describes an algorithm to generate anatomically adaptive quadrilateral meshes that can be used with the two-dimensional FE-based non-rigid registration method in [5]. The algorithm takes a triangular mesh built from an image segmentation as input, and produces a convex quadrangulation which preserves local variations in element size and imposes an upper bound on the number of elements in the mesh. Registrations are performed using a variety of meshes, and the results demonstrate that anatomically-specific quadrilateral meshes produce equivalent results to much denser grids of quadrilaterals in a fraction of the time. Results using brain MR data are presented here, but the algorithm is in no way restricted to any particular imaging modality or body part.
2
Preliminaries
A polygonal region R is a connected region of the plane whose boundary is one polygon or the union of a finite number of disjoint polygons. When the boundary of R consists of 1 + k polygons with k > 0, we say that R has k holes. Vertices and edges of a polygonal region are the vertices and edges of its bounding polygons. We denote the set of vertices (edges) of R by V (R) (E(R)). In this paper, we are primarily interested in closed and bounded polygonal regions, i.e., polygonal regions that contain their bounding polygons and for which we can define a “bounding rectangle” of finite size. Unless otherwise specified, we refer to these closed and bounded polygonal regions simply as polygonal regions. A triangulation T of R is a subdivision of R into triangles satisfying the following two conditions: (1) if t1 and t2 are two triangles of T then t1 ∩ t2 is either empty, a common vertex of t1 and t2 , or a common edge of t1 and t2 ; and (2) the union
of all triangles of T , the domain |T | of T , is R. We denote the set of vertices of T by V (T ). If R has n vertices and k holes, and T is a triangulation of R with m vertices such that V (R) ⊆ V (T ), we have that T has t = 2m + 2k − 2 − mb triangles, where mb is the number of vertices of T and n ≤ mb ≤ m. The dual graph of T is an undirected graph that has a vertex for each triangle of T and an edge for each pair of adjacent triangles of T . A spanning tree TG of an undirected graph G is a subgraph of G that is a tree and that contains all the vertices of G. A spanning tree of the dual graph of a triangulation is always a binary tree. A quadrangulation Q of R is a subdivision of R into quadrilaterals satisfying the following two conditions: (1) if q1 and q2 are two quadrilaterals of Q then q1 ∩ q2 is either empty, a common vertex of q1 and q2 , or a common edge of q1 and q2 ; and (2) the union of all quadrilaterals of Q, the domain |Q| of Q, is R. We denote the set of vertices of Q by V (Q). We say that Q is a convex quadrangulation of R if every quadrilateral q of Q is convex, i.e., the quadrilateral q has no angle greater than 180 degrees. If each of the four angles of every quadrilateral q in Q is strictly less than 180 degrees, we say that Q is a strictly convex quadrangulation of R.
3
Convex Quadrangulations
Convex and strictly convex quadrangulations are the two classes of quadrilateral meshes that can be used in FE applications. Several algorithms for generating convex quadrangulations have been proposed, [10]. Some of these algorithms adopt an indirect approach: the polygonal region is first triangulated, and then the triangulation is converted into a convex quadrangulation of the polygonal region, [12], [13], [14], [15]. Since there are several algorithms for generating provably good triangulations of polygonal regions, quadrangulation algorithms that employ the indirect approach may take advantage of this fact by attempting to preserve the features of the input triangulation. Johnston et al., [12], introduced an algorithm that uses several heuristics for automatically converting a triangulation of a general polygonal region into a strictly convex quadrangulation. Their algorithm runs in O(t2 ) for a triangulation with t triangles. Everett et al., [13], presented an algorithm that generates at most 8t 3 strictly convex quadrilaterals and runs in linear time in t. Shimada et al., [14], proposed an indirect algorithm that optimizes several mesh quality criteria, but becomes very expensive when quadrangulating polygonal regions with complex geometry. Owen et al., [15], presented another indirect algorithm that considers some of the mesh quality criteria in [14] during the conversion process. Their algorithm has the advantage of being simpler than the latter, and has been implemented in commercial software. However, it does not guarantee a bounded-size output quadrangulation, even though it accounts for mesh quality. It is unclear how the size and quality of the output quadrangulation would be affected by an input triangulation of a complex domain.
Ramaswami et al., [8], presented an algorithm called Q-percolation for converting triangulations of general polygonal regions into bounded-size quadrangulations, where the bounds are considerably better than those in [13]. However, Q-percolation does not guarantee that the output quadrangulation is convex. By modifying the Q-percolation algorithm, we have developed a new algorithm (CQpercolation) that converts triangulations of general polygonal regions with holes into convex quadrangulations. CQ-percolation always generates convex quadrangulations with at most 3t 2 quadrilaterals, which is a better bound than the one in [13]. CQ-percolation also preserves the local element size variation of the input triangulation, which is highly desirable if the input triangulation was built to adaptively subdivide the polygonal region. The underlying premise of both Q-percolation and CQ-percolation is to quadrangulate a small group of triangles at a time until the entire quadrangulation is complete. Each group represents a triangulation of a small, simple polygonal region, such as a non-convex quadrilateral, pentagon, hexagon, or septagon. Both algorithms form a group of triangles at a time and convert it into a partial or complete quadrangulation before considering the next group. By using a spanning tree of the dual graph of the input triangulation, Q-percolation and CQ-percolation systematically group triangles together so that no isolated triangles remain in the resulting decomposition. The two algorithms differ somewhat in their approaches to the grouping and conversion processes.
(a)
(b)
(c)
Fig. 1. Mesh generation from image contours. (a) Contours from a segmented image. (b) Triangulation. (c) Quadrangulation.
Although CQ-percolation does not take into account the shape quality of the quadrilaterals produced, our registration results in Sec. 5 show that the quality of the output quadrangulations is comparable to those of their triangular counterparts. By using well-known methods to improve and smooth mesh topology, we further improve our mesh quality by making our quadrangulations strictly convex.
4
Meshes from Images
In order to use our quadrangulation algorithm to generate a mesh from an image, we must first triangulate the image domain. Since most triangulation algorithms take a polygonal region as input, we must convert the raw, pixel-based representation of an image into a polygonal representation. This is done by segmenting the image and then generating a closed polygonal contour for every distinct anatomic boundary of interest. The contours are generated by adaptively sampling the set of pixels defining each boundary, and then connecting them with straight line segments. The resulting polygonal contours are always disjoint, but may define several nested layers of polygonal contours. Figure 1a shows the contours obtained from a segmented image of a human brain. Once we have the contours of each distinct anatomic feature of the segmented image, we group them together to define polygonal regions. Next, we generate a triangular mesh for each polygonal region. This must be done systematically since the triangulation and quadrangulation algorithms are allowed to insert points in the polygonal contours. We first mesh the innermost polygonal regions, and only allow insertion of vertices in the outer contours of the polygonal regions. This protects any contours that are part of a polygonal region that has already been meshed. Figure 1b shows a triangular mesh built from the polygonal regions defined by the contours in figure 1a, and figure 1c shows the quadrilateral mesh obtained by CQ-percolation of the triangular mesh. Note that there is a rectangular contour enclosing the contours defining the actual image. This rectangular contour is intentionally inserted to support the image registration algorithm.
5
Results
(a)
(b)
(c)
(d)
Fig. 2. Inputs to and outputs from the FE-based registration using a particular quadrangulation. (a) Coronal brain MR image. (b) Warped version of (a) with quadrilateral mesh overlaid. (c) Registration result using the mesh in (b). (d) RMS error between (b) and (c).
In order to evaluate the quadrilateral meshes generated by our algorithm, we performed a quantitative, registration-based analysis using a pair (A, B) of 2D
images of the human brain. A (figure 2a) was a 256x256-pixel coronal MR image, and B (figure 2b) was the result of applying a cubic polynomial warp to A. We used the approach described in Section 4 to generate a variety of triangular and quadrilateral meshes of A and B. Next, we used these meshes to register A to B with the FE-based elastic registration method in the Insight Segmentation and Registration Toolkit (ITK) (http://www.itk.org).
Mesh #Elements #Vertices #Int. Pts. Runtime (s) RMS 1 1 2 2 3 3 4 4 5 6 7 8 9 10 11 12 13 14 15 16
2921 2991 3549 3549 4914 4914 8254 8254 1024 4096 16384 65536 1645 1941 2581 4318 1773 2073 2747 4499
1472 1472 1790 1790 2481 2481 4173 4173 1089 4225 16641 66049 1654 1957 2605 4364 1785 2089 2771 4545
1 3 1 3 1 3 1 3 4 4 4 4 4 4 4 4 4 4 4 4
10 17 12 21 17 30 32 52 12 46 205 1001 21 24 33 59 23 27 35 62
17.95 17.25 17.35 16.98 17.20 16.81 16.78 16.62 18.56 16.99 16.11 15.93 17.67 17.12 16.93 16.68 17.44 17.47 16.92 16.66
Table 1. Summary of registration results using different types of meshes. Meshes 1– 4 are triangular meshes. Meshes 5–8 are uniform grids of quadrilaterals with 8 × 8, 4 × 4, 2 × 2 and 1 × 1 pixel elements, respectively. Meshes 9–12 are quadrilateral meshes generated by the algorithm in Sec. 3 and the approach in Sec. 4 from triangular meshes 1–4, respectively. Meshes 13–16 are the result of topological and geometric improvements over meshes 9–12, respectively.
We created triangulations with minimum angles of 20o , 25o , 30o , and 33o using Triangle, [16], a publicly available mesh generator. These meshes served as inputs to our quadrangulation algorithm. We also used quadrilateral meshes whose elements were arranged in uniform grids of various resolutions; these were generated automatically within the registration application. Each mesh was sampled using different numbers of integration points (required by the FE code for Gaussian integration over the element) to examine the relationship between element size and number of integration points with respect to registration accuracy. All registrations were run for 20 iterations on a PC with an Intel Pentium III processor running Windows 98. We evaluated the results of the registrations by
calculating the RMS (root mean squared) difference between the intensity values of corresponding pixels over the entire image domain of (A, B). Each row in table 1 shows the results obtained for a particular mesh and number of integration points. The meshes generated by our quadrangulation algorithm contain about 60% as many elements as their input triangulations, although the quadrilateral meshes have more vertices. This reduction in element density does not appear to affect the RMS error of the registration, since the triangular and quadrilateral meshes perform comparably. However, the runtimes for the quadrilateral meshes are slightly higher than those of their triangular counterparts. This is most likely due to the higher vertex number in the quadrilateral meshes as well as the added expense of computing solutions using quadrilaterals in the FE code. As seen in figure 2d, the majority of the RMS error is associated with edge mismatch along the sulci of the brain. The remainder of the image is well-matched by the registration. Note that meshes 11 and 15 have less than 75% as many quadrilaterals as mesh 6 but still report a lower RMS. This result supports the observation that anatomically-adapted meshes provide better registration results in less time than uniform grids of similar size. This is mainly due to the fact that adaptive meshes, such as meshes 1–4 and 9–16, can approximate the inherent geometry of an image more precisely and position smaller elements along the boundary (where an accurate warp is more difficult to compute). In addition, by using adaptive meshes in which every element belongs to only one region, we easily support the modeling of different image regions with anatomically appropriate material properties. This is more difficult to achieve using uniform grids without dramatically increasing both the mesh resolution and, consequently, the runtime (see mesh 8).
6
Conclusions and Future Work
We present an algorithm to convert triangulations of polygonal regions with polygonal holes into convex quadrangulations. Our algorithm has a runtime linear in the number of triangles, and offers better bounds than similar algorithms that also produce strictly convex quadrilateral meshes of bounded size. Using this algorithm and traditional techniques to improve mesh quality, we describe an approach for generating strictly convex quadrilateral meshes that approximate the geometry of an image and prevent elements from belonging to more than one distinct image region. In addition, we discuss a quantitative analysis of these quadrilateral meshes, their triangular counterparts and uniform grids, using a particular implementation of a FE-based, image registration method. Our analysis demonstrates that less dense quadrilateral meshes derived from triangulations effectively capture the transformation of one image into another compared to triangular and uniform quadrilateral meshes of similar density. Future work will focus on registering additional pairs of images with this protocol, incorporating anatomically appropriate material properties into the scheme, improving correspondence along image boundaries and rigorously ex-
amining bounds on angles of quadrilateral meshes generated from triangular meshes with quality constraints.
References 1. Maintz, J.B.A., Viergever, M.A.: A survey of medical image registration. Medical Image Analysis 2 (1998) 1–3 2. Gee, J., Sundaram, T., Hasegawa, I., Uematsu, H., Hatabu, H.: Characterization of Regional Pulmonary Mechanics from Serial MR Data. In T. Dohi, R.K., ed.: Medical Image Computing and Computer-Assisted Intervention - MICCAI 2002. (2002) 762–769 3. Rueckert, D., Sonada, L.I., Hayes, C., Hill, D.L.G., Leach, M.O., Hawkes, D.J.: Nonrigid Registration Using Free-Form Deformations: Application to Breast MR Images. IEE Transactions on Medical Imaging 18 (1999) 712–721 4. Broit, C.: Optimal Registration of Deformed Images. PhD thesis, Department of Information and Computer Science, University of Pennsylvania, Philadelphia, PA, USA (1981) 5. Gee, J., Bajcsy, R.: Elastic Matching: Continuum Mechanical and Probabilistic Analysis. In Toga, A., ed.: Brain Warping. Academic Press (1999) 6. Bern, M., Plassmann, P.: Mesh Generation. In Sack, J., Urrutia, J., eds.: Handbook of Computational Geometry. Elsevier Science (1999) 7. Bern, M., Eppstein, D.: Mesh Generation and Optimal Triangulation. In Hwang, F., Du, D.Z., eds.: Computing in Euclidean Geometry. World Scientific (1992) 8. Ramaswami, S., Ramos, P., Toussaint, G.: Converting Triangulations to Quadrangulations. Computational Geometry: Theory and Applications 9 (1998) 257–276 9. Bern, M., Eppstein, D.: Quadrilateral Meshing by Circle Packing. In: Proceedings of the 6th International Meshing Roundtable. (1997) 7–19 10. Owen, S.: A Survey of Unstructured Mesh Generation Technology. In: Proceedings of the 7th International Meshing Roundtable, Dearborn, MI (1998) 239–267 11. Mohamed, A., Davatzikos, C.: A Simple Approach for Three-Dimensional FiniteElement Mesh Generation from Segmented Medical Images. submitted to Computer Vision and Pattern Recognition, 2003 (2003) 12. Johnston, B., Sullivan, J., Kwasnik, A.: Automatic Conversion of Triangular Finite Meshes to Quadrilateral Meshes. International Journal of Numerical Methods in Engineering 31 (1991) 67–84 13. Everett, H., Lenhart, W., Overmars, M., Shermer, T., Urrutia, J.: Strictly Convex Quadrilateralizations of Polygons. In: Proceedings of the 4th Canadian Conference on Computational Geometry. (1992) 77–82 14. Shimada, K., Liao, J.H., Itoh, T.: Quadrilateral Meshing with Directionality Control through the Packing of Square Cells. In: Proceedings of the 7th International Meshing Roundtable, Dearborn, Michigan, USA (1998) 61–76 15. Owen, S., Staten, M., Cannan, S., Saigal, S.: Advancing Front Quadrilateral Meshing Using Triangle Transformations. In: Proceedings of the 7th International Meshing Roundtable, Dearborn, Michigan, USA (1998) 409–428 16. Shewchuk, J.: Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator. In: First Workshop on Applied Computational Geometry, Philadelphia, PA, USA (1996) 124–133
