Octree Grid Topology Preserving Geometric Deformable Model ... - IACL
Octree Grid Topology Preserving Geometric Deformable Model for Three-Dimensional Medical Image Segmentation Ying Bai∗ , Xiao Han† , and Jerry L. Prince∗ ∗
Johns Hopkins University, Baltimore MD 21218 † CMS, Inc., St. Louis, MO 63132
Abstract. Topology-preserving geometric deformable models (TGDMs) are used to segment objects that have a known topology. Their accuracy is inherently limited, however, by the resolution of the underlying computational grid. Although this can be overcome by using fine-resolution grids, both the computational cost and the size of the resulting surface increase dramatically. In order to maintain computational efficiency and to keep the surface mesh size manageable, we have developed a new framework, termed OTGDMs, for topology-preserving geometric deformable models on balanced octree grids (BOGs). In order to do this, definitions and concepts from digital topology on regular grids were extended to BOGs so that characterization of simple points could be made. Other issues critical to the implementation of OTGDMs are also addressed. We demonstrate the performance of the proposed method using both mathematical phantoms and real medical images.
1
Introduction
Front propagation using level set methods [1] and their application in deformable models – geometric deformable models (GDMs) [2–4]are well established and extensively used in medical image segmentation. Topology preserving geometric deformable models (TGDMs) [5–7] were recently introduced in order to provide the ability to maintain topology of segmented objects while preserving the other benefits of GDMs. For example, in medical imaging many organs to be segmented have boundary topologies equivalent to that of a sphere. While many applications such as visualization and quantification may not require topologically correct segmentations, there are some applications — e.g., surface mapping and flattening and shape atlas generation — that cannot be achieved without correct topology of the segmented objects. GDMs represent the evolving surface implicitly as a level set of a higher dimensional function. The resolution of the implicit surface is therefore restricted by the resolution of the sampling grid that defines the level set function, as demonstrated in Figs. 1(a)–(c). Accurate solution and representation of shapes with fine anatomical details (e.g., the folded sulci and gyri on the cortex) requires the use of a fine resolution grid. This dramatically increases the computation time of GDMs and produces surface meshes with prohibitive size, however, especially
Fig. 1. Implicit surface resolution: the dotted contour is the truth contour; the solid contour is the implicit contour embedded in each sampling grid. (a) A coarse resolution grid cannot resolve contour details. (b) A refined grid represents the truth contour better. (c) A more refined grid provides a more accurate representation. (d) An adaptive grid with local refinement provides an accurate and efficient multiresolution shape representation.
on highly resolved 3D medical images. Adaptive grid techniques [8–13] address the resolution problem of GDMs by locally refining the sampling grid in order to resolve details and concentrate computational efforts where more accuracy is needed (as shown in Fig. 1(d)). Although numerical schemes to implement level set methods on adaptive grids are well developed, there is little literature on the issue of defining digital connectivity rules for adaptive grids. Without such rules, it is difficult to guarantee homeomorphisms between the implicit surfaces and the corresponding boundaries of segmented objects on an adaptive grid. Digital connectivity rules for adaptive grids are also necessary in order to design a topology preserving level set method on adaptive grids. The method introduced by Han et al. [5] for regular grids maintains the topology of the implicit surface by controlling the topology of the corresponding binary object segmentation. This is achieved by applying the simple point criterion [14] from the theory of digital topology [15], preventing the level set function from changing sign at non-simple points. Until now, this topology preservation mechanism could not be used on adaptive grids because there was no characterization of “simple points” on adaptive grids. In this paper, we propose a new topology-preserving level set method based on the balanced octree grids (BOGs) (i.e., octree grids for which the maximum cell edge length ratio between adjacent grid cells is 2). We first briefly review the digital topology framework for the adaptive grid that we recently proposed [16]. We then present a topology preserving geometric deformable model for adaptive octree grids (OTGDM), which is based on our new characterization of simple points on BOGs that extends the original characterization on the uniform grid in [14]. Several experiments are used to demonstrate the performance of OTGDM on both computational phantoms and real medical images.
2
Digital topology framework on BOGs
In [16], we extended basic digital topology concepts to BOGs, providing a unique and topologically consistent digital embedding scheme for implicit surfaces defined on BOGs. In the following, we summarize the concepts of “neighbor points” and “invalid cases”, which will be used in the characterization of “simple points” on BOGs presented later. The concept of neighbor points is fundamental in classical digital topology theory [15]. In [16], grid points on an octree grid are defined to be edge(E)neighbors, square(S)-neighbors, or cube(C)-neighbors if they share an edge, a
Fig. 2. 3D neighborhoods on balanced octree grids.
face, or a cube, respectively, of leaf cells of the octree. (Leaf cells are cells that have no child cells.) Fig. 2 shows an example of neighborhoods on a BOG. Fig. 2(a) shows a uniform neighborhood and Figs. 2(b)–(d) show examples of non-uniform neighborhoods. The white circle in each figure indicates the root point of the neighborhoods; black squares are the E-neighbors; white squares are the points that are added to the E-neighbors to yield the S-neighbors; and gray squares are the points that are added to S-neighbors to yield the C-neighbors. Analogous definitions of neighborhood, adjacency, path, and connectivity can be found in [16]. Using the above neighborhood and connectivity definitions, inconsistencies can still occur at the interface(s) between grid cells of different resolution (referred as transition face) as illustrated in Fig. 3. Assume E-connectivity for the
Fig. 3. Examples of invalid cases on interface of resolution transition.
black foreground points, and S-connectivity for the white background points. In Fig. 3(a), the black point is shared only by the fine resolution grid cells and forms a single connected component on the transition face. This foreground component does not exist, however, if we look at the coarse side of the transition face. Similarly, in Fig. 3(b), the two white background points are disconnected and form two connected components at the fine resolution side of the transition face, whereas they are connected as one connected component by the coarse cell. As a result, the embedded isosurface will have a discontinuity at the transition face. To eliminate these inconsistencies and thus guarantee unique and valid surface embedding on BOGs, we define such cases as invalid case [16] where the numbers of connected components (for both foreground and background) formed by the grid points on two sides of a transition face are not equal. Invalid cases are explicitly checked and prevented to happen in the OTGDM algorithm, as will be described later.
3
A Topology-Preserving Level Set Method on BOGs
In this section, we present a new topology-preserving level set method on BOGs. The overall algorithm is first summarized and the details about each step are then discussed. We adopt the narrowband framework [17] in the following implementation and we assume a general GDM model as can be summarized by the following equation: ∂Φ(x, t) = [Fprop (x, t) + Fcurv (x, t)]|∇Φ(x, t)| + F adv (x, t) · ∇Φ(x, t) ∂t
(1)
where Fprop , Fcurv , and F adv denote user-designed force (or speed) terms that control the model deformation. In particular, Fcurv , the curvature force, controls the regularity (smoothness) of the implicit surface. Fprop and F adv are two forms of image forces (scalar and vector respectively) that drive the surface to the desired object boundary. Octree-based TGDM algorithm 1. Initialize the adaptive grid according to the initial surface topology and adaptation metric (cf. Section 3.1). Initialize the level set function to be the signed distance function of the initial surface. 2. Build the narrow band by finding all grid points within a distance threshold of the implicit surface (zero level set of the current level set function). 3. Update the level set function at each point in the narrow band iteratively as follows: (a) Compute the new value of Φ(x, t) using Eq. (1). (b) If there is no sign change, accept the new value and move on to the next point. Otherwise, go to Step 3(c). (c) Test whether the sign change at this point yields a valid configuration (cf. Section 2). If yes, go to Step 3(d). Otherwise, move on to the next point. (d) Test whether the current point is a simple point by computing two topological numbers (cf. Section 3.2)). If the point is simple accept the new value. Otherwise, set the level set function to be a small number with the same sign as its original value. 4. If the zero level set is near the boundary of the current narrow band, reinitialize the level set function to be a signed distance function and go to Step 2. 5. Test whether the zero level set has stopped moving (i.e., no sign change happens at any point inside the narrow band in two or three consecutive iterations). If yes, stop; otherwise, go to the next iteration. 6. Extract the zero-value isosurface using an adaptive connectivity-consistent marching cells algorithm (cf. Section 3.3). A few comments about OTGDM. First, the reinitialization step is a straightforward extension of the fast marching method to the non-uniform cartesian grid. Different grid sizes are handled by the modified finite difference operator [17]. Second, the simple point check can be omitted, and the algorithm becomes a standard geometric deformable model on an adaptive octree grid (OSGDM).
3.1
Adaptive grid generation
We generate a BOG following two criteria: first, the BOG should embed the initial implicit surface with correct topology; second, the BOG should adapt its resolution according to the geometrical shape of the final surface. We discuss these two considerations below. Initial surface topology To guarantee that the final surface has the correct topology, we must start with an initial surface that has the correct topology. Assume that we are given an initial implicit surface defined on the uniformly sampled grid; by definition, it has the correct topology, which can be arbitrary and unknown. To initialize an implicit surface on a BOG while preserving topology, we apply the bottom-up cell merging algorithm that we presented in [16]. The algorithm starts from the original uniform grid and treats it as an octree grid that is at its finest possible resolution. The leaf cells of this octree grid are then traversed level-by-level from bottom to up. At each level, the leaf cells are evaluated one-by-one to see if they can be merged without changing the topology of the underlying isosurface and without generating invalid cases. Details of this algorithm can be found in [16]. After applying the cell-merging algorithm and balancing the grid, we obtain a BOG that embeds the initial surface with the correct topology. Image-based adaptivity metric We now refine the initial BOG so that it has finer cells where image details predict the need for a higher resolution surface representation. The concept of an adaptivity metric, which is derived from the image volume, is used. This metric estimates the local geometrical properties of the final surface boundary — if it should come to rest at the given image location — and the BOG is adapted accordingly. A classical adaptivity metric is the magnitude of the image gradient (cf. [8]), wherein computational grid is refined at high gradient regions and coarsened elsewhere. This metric cannot help to reduce the size of the final surface mesh on the adaptive octree grid, however, since the grid will be uniformly refined along the entire object boundary. In this work, we have used the max-curvature — i.e., the larger absolute value of the two principal curvatures — estimated from the volumetric gray-level image using the method in [18]. To achieve more robust estimation (in noise, for example), we apply an anisotropic smoothing to the image before using curvature estimation (cf. [19]). Once curvature κ(x) is estimated, we define a refinement rule to be: l(x) = i,
if ti−1
Johns Hopkins University, Baltimore MD 21218 † CMS, Inc., St. Louis, MO 63132
Abstract. Topology-preserving geometric deformable models (TGDMs) are used to segment objects that have a known topology. Their accuracy is inherently limited, however, by the resolution of the underlying computational grid. Although this can be overcome by using fine-resolution grids, both the computational cost and the size of the resulting surface increase dramatically. In order to maintain computational efficiency and to keep the surface mesh size manageable, we have developed a new framework, termed OTGDMs, for topology-preserving geometric deformable models on balanced octree grids (BOGs). In order to do this, definitions and concepts from digital topology on regular grids were extended to BOGs so that characterization of simple points could be made. Other issues critical to the implementation of OTGDMs are also addressed. We demonstrate the performance of the proposed method using both mathematical phantoms and real medical images.
1
Introduction
Front propagation using level set methods [1] and their application in deformable models – geometric deformable models (GDMs) [2–4]are well established and extensively used in medical image segmentation. Topology preserving geometric deformable models (TGDMs) [5–7] were recently introduced in order to provide the ability to maintain topology of segmented objects while preserving the other benefits of GDMs. For example, in medical imaging many organs to be segmented have boundary topologies equivalent to that of a sphere. While many applications such as visualization and quantification may not require topologically correct segmentations, there are some applications — e.g., surface mapping and flattening and shape atlas generation — that cannot be achieved without correct topology of the segmented objects. GDMs represent the evolving surface implicitly as a level set of a higher dimensional function. The resolution of the implicit surface is therefore restricted by the resolution of the sampling grid that defines the level set function, as demonstrated in Figs. 1(a)–(c). Accurate solution and representation of shapes with fine anatomical details (e.g., the folded sulci and gyri on the cortex) requires the use of a fine resolution grid. This dramatically increases the computation time of GDMs and produces surface meshes with prohibitive size, however, especially
Fig. 1. Implicit surface resolution: the dotted contour is the truth contour; the solid contour is the implicit contour embedded in each sampling grid. (a) A coarse resolution grid cannot resolve contour details. (b) A refined grid represents the truth contour better. (c) A more refined grid provides a more accurate representation. (d) An adaptive grid with local refinement provides an accurate and efficient multiresolution shape representation.
on highly resolved 3D medical images. Adaptive grid techniques [8–13] address the resolution problem of GDMs by locally refining the sampling grid in order to resolve details and concentrate computational efforts where more accuracy is needed (as shown in Fig. 1(d)). Although numerical schemes to implement level set methods on adaptive grids are well developed, there is little literature on the issue of defining digital connectivity rules for adaptive grids. Without such rules, it is difficult to guarantee homeomorphisms between the implicit surfaces and the corresponding boundaries of segmented objects on an adaptive grid. Digital connectivity rules for adaptive grids are also necessary in order to design a topology preserving level set method on adaptive grids. The method introduced by Han et al. [5] for regular grids maintains the topology of the implicit surface by controlling the topology of the corresponding binary object segmentation. This is achieved by applying the simple point criterion [14] from the theory of digital topology [15], preventing the level set function from changing sign at non-simple points. Until now, this topology preservation mechanism could not be used on adaptive grids because there was no characterization of “simple points” on adaptive grids. In this paper, we propose a new topology-preserving level set method based on the balanced octree grids (BOGs) (i.e., octree grids for which the maximum cell edge length ratio between adjacent grid cells is 2). We first briefly review the digital topology framework for the adaptive grid that we recently proposed [16]. We then present a topology preserving geometric deformable model for adaptive octree grids (OTGDM), which is based on our new characterization of simple points on BOGs that extends the original characterization on the uniform grid in [14]. Several experiments are used to demonstrate the performance of OTGDM on both computational phantoms and real medical images.
2
Digital topology framework on BOGs
In [16], we extended basic digital topology concepts to BOGs, providing a unique and topologically consistent digital embedding scheme for implicit surfaces defined on BOGs. In the following, we summarize the concepts of “neighbor points” and “invalid cases”, which will be used in the characterization of “simple points” on BOGs presented later. The concept of neighbor points is fundamental in classical digital topology theory [15]. In [16], grid points on an octree grid are defined to be edge(E)neighbors, square(S)-neighbors, or cube(C)-neighbors if they share an edge, a
Fig. 2. 3D neighborhoods on balanced octree grids.
face, or a cube, respectively, of leaf cells of the octree. (Leaf cells are cells that have no child cells.) Fig. 2 shows an example of neighborhoods on a BOG. Fig. 2(a) shows a uniform neighborhood and Figs. 2(b)–(d) show examples of non-uniform neighborhoods. The white circle in each figure indicates the root point of the neighborhoods; black squares are the E-neighbors; white squares are the points that are added to the E-neighbors to yield the S-neighbors; and gray squares are the points that are added to S-neighbors to yield the C-neighbors. Analogous definitions of neighborhood, adjacency, path, and connectivity can be found in [16]. Using the above neighborhood and connectivity definitions, inconsistencies can still occur at the interface(s) between grid cells of different resolution (referred as transition face) as illustrated in Fig. 3. Assume E-connectivity for the
Fig. 3. Examples of invalid cases on interface of resolution transition.
black foreground points, and S-connectivity for the white background points. In Fig. 3(a), the black point is shared only by the fine resolution grid cells and forms a single connected component on the transition face. This foreground component does not exist, however, if we look at the coarse side of the transition face. Similarly, in Fig. 3(b), the two white background points are disconnected and form two connected components at the fine resolution side of the transition face, whereas they are connected as one connected component by the coarse cell. As a result, the embedded isosurface will have a discontinuity at the transition face. To eliminate these inconsistencies and thus guarantee unique and valid surface embedding on BOGs, we define such cases as invalid case [16] where the numbers of connected components (for both foreground and background) formed by the grid points on two sides of a transition face are not equal. Invalid cases are explicitly checked and prevented to happen in the OTGDM algorithm, as will be described later.
3
A Topology-Preserving Level Set Method on BOGs
In this section, we present a new topology-preserving level set method on BOGs. The overall algorithm is first summarized and the details about each step are then discussed. We adopt the narrowband framework [17] in the following implementation and we assume a general GDM model as can be summarized by the following equation: ∂Φ(x, t) = [Fprop (x, t) + Fcurv (x, t)]|∇Φ(x, t)| + F adv (x, t) · ∇Φ(x, t) ∂t
(1)
where Fprop , Fcurv , and F adv denote user-designed force (or speed) terms that control the model deformation. In particular, Fcurv , the curvature force, controls the regularity (smoothness) of the implicit surface. Fprop and F adv are two forms of image forces (scalar and vector respectively) that drive the surface to the desired object boundary. Octree-based TGDM algorithm 1. Initialize the adaptive grid according to the initial surface topology and adaptation metric (cf. Section 3.1). Initialize the level set function to be the signed distance function of the initial surface. 2. Build the narrow band by finding all grid points within a distance threshold of the implicit surface (zero level set of the current level set function). 3. Update the level set function at each point in the narrow band iteratively as follows: (a) Compute the new value of Φ(x, t) using Eq. (1). (b) If there is no sign change, accept the new value and move on to the next point. Otherwise, go to Step 3(c). (c) Test whether the sign change at this point yields a valid configuration (cf. Section 2). If yes, go to Step 3(d). Otherwise, move on to the next point. (d) Test whether the current point is a simple point by computing two topological numbers (cf. Section 3.2)). If the point is simple accept the new value. Otherwise, set the level set function to be a small number with the same sign as its original value. 4. If the zero level set is near the boundary of the current narrow band, reinitialize the level set function to be a signed distance function and go to Step 2. 5. Test whether the zero level set has stopped moving (i.e., no sign change happens at any point inside the narrow band in two or three consecutive iterations). If yes, stop; otherwise, go to the next iteration. 6. Extract the zero-value isosurface using an adaptive connectivity-consistent marching cells algorithm (cf. Section 3.3). A few comments about OTGDM. First, the reinitialization step is a straightforward extension of the fast marching method to the non-uniform cartesian grid. Different grid sizes are handled by the modified finite difference operator [17]. Second, the simple point check can be omitted, and the algorithm becomes a standard geometric deformable model on an adaptive octree grid (OSGDM).
3.1
Adaptive grid generation
We generate a BOG following two criteria: first, the BOG should embed the initial implicit surface with correct topology; second, the BOG should adapt its resolution according to the geometrical shape of the final surface. We discuss these two considerations below. Initial surface topology To guarantee that the final surface has the correct topology, we must start with an initial surface that has the correct topology. Assume that we are given an initial implicit surface defined on the uniformly sampled grid; by definition, it has the correct topology, which can be arbitrary and unknown. To initialize an implicit surface on a BOG while preserving topology, we apply the bottom-up cell merging algorithm that we presented in [16]. The algorithm starts from the original uniform grid and treats it as an octree grid that is at its finest possible resolution. The leaf cells of this octree grid are then traversed level-by-level from bottom to up. At each level, the leaf cells are evaluated one-by-one to see if they can be merged without changing the topology of the underlying isosurface and without generating invalid cases. Details of this algorithm can be found in [16]. After applying the cell-merging algorithm and balancing the grid, we obtain a BOG that embeds the initial surface with the correct topology. Image-based adaptivity metric We now refine the initial BOG so that it has finer cells where image details predict the need for a higher resolution surface representation. The concept of an adaptivity metric, which is derived from the image volume, is used. This metric estimates the local geometrical properties of the final surface boundary — if it should come to rest at the given image location — and the BOG is adapted accordingly. A classical adaptivity metric is the magnitude of the image gradient (cf. [8]), wherein computational grid is refined at high gradient regions and coarsened elsewhere. This metric cannot help to reduce the size of the final surface mesh on the adaptive octree grid, however, since the grid will be uniformly refined along the entire object boundary. In this work, we have used the max-curvature — i.e., the larger absolute value of the two principal curvatures — estimated from the volumetric gray-level image using the method in [18]. To achieve more robust estimation (in noise, for example), we apply an anisotropic smoothing to the image before using curvature estimation (cf. [19]). Once curvature κ(x) is estimated, we define a refinement rule to be: l(x) = i,
if ti−1
