A Multiresolution PDE-Based Deformable Surface for Medical Imaging Applications

We recently developed a multiresolution PDE-based deformable surface whose deformation behavior is governed by partial differential equations (PDEs) such as the weighted minimal surface flow. Comparing with the level-set approach, our new model has better control of the mesh quality and model resolution, and is much simpler to implement since all the computations are local. The new deformable model is very useful for a variety of medical imaging applications including boundary reconstruction, surface visualization, data segmentation, and topology discovery. In this paper, we demonstrate both the accuracy and robustness of our model on areas such as medical image segmentation through a number of experiments on both real (MRI/CT) and synthetic volumetric datasets.


INTRODUCTION
Since the seminal work of the snake model proposed by Kass et al. [1] in 1988, deformable models have significantly gained popularity in medical imaging analysis and related fields. Conventional deformable models are essentially parametric models [2][3][4]. In general, they have concise representation and it is easy to incorporate user constraints in them. However, it is very difficult for parameterized deformable models to represent shapes of arbitrary topology. Recently, implicit snakes were proposed by Malladi et al. [5] and Caselles et al. [6] with the ability to handle topology changes. Their schemes are based on the modeling of propagating fronts which are the level set of certain scalar functions. Nonetheless, the topological flexibility is accompanied by a time-consuming integration method that must be applied to a higher-dimensional space, and the desirable shape must be explicitly evaluated using the marching-cube-like [3] techniques in a separate postprocessing stage.
In order to bridge the divide between explicit parametric models and implicit level-set models [7,8], most recently, researchers have proposed a new kind of deformable model-topology-adaptive explicit deformable models such as the topology-adaptive snake proposed by McInerney and Terzopoulos [9,10] and the discrete triangle model of Lachaud and Montanvert [11]. Following the similar research direction, we had also developed a new topologyadaptive explicit deformable model. Our new model is a multiresolution PDE-based deformable surface whose deformation behavior is governed by PDEs such as the weighted minimal surface flow. The new model can either grow from the inside or shrink from the outside. Furthermore, our new model supports different LODs through the use of both global subdivision and local/adaptive subdivision. Comparing with existing topology-adaptive explicit models [9][10][11], our new model has better control of the mesh quality and model resolution, and is much simpler to implement since all the computations are local. The new deformable model will be very useful for a variety of medical image applications including boundary reconstruction, surface visualization, data segmentation, and topology discovery. In this paper, we illustrate some of the experimental results of our new model on medical image segmentation for both real (MRI/CT) and synthetic volumetric datasets.
Since the proposed deformation framework is quite general, it can also be applied to other application domains that are outside the scope of biomedical imaging. For example, it can be applied for 3D shape reconstruction from 2D 2 International Journal of Biomedical Imaging multiple-view images and unorganized point clouds [12]. It can also be applied for interactive mesh editing, sketching, morphing, as well as shape warping [13].

Model initialization
The model can be initialized either interactively by the user as a small seed model that can grow from the inside, or automatically by the system as a bounding box that will shrink from the outside. Any closed polyhedron can be used as a seed model. If the model is initialized as a simple bounding box shrink-wrapping from the outside, several iterations of global refinement are usually needed to increase the model resolution so as to correctly capture the geometry and topology of the shape. In this paper, we use Loop's subdivision scheme [14] for the purpose of global refinement.

Model growing
After the model is initialized, the model will grow (or shrink) and deform. The deformation of the model is governed by a PDE called weighted minimal surface flow proposed by Caselles et al. [6]: (1) is the 3D deformable surface, t is the time parameter, and S 0 is the initial shape of the surface. H is the mean curvature value of the surface, ρ N is the unit normal vector of the surface, and v is a constant speed that will enable the convex initial shape to capture nonconvex, arbitrary complicated shapes. The nonzero constant velocity is also useful to avoid the model getting stuck into the local minimum during the evolution process. g is the monotonic, nonincreasing, nonnegative weight function that enables the model to interact with the datasets, and will stop the deformation of the model when it reaches the object boundary. In this paper, g is defined as the commonly used 3D edge detector: Here, I is the volumetric density function, and G σ * I is the smoothed density function by convoluting with a Gaussian filter with variance σ. The surface evolution process is approximated using an explicit iterative equation: where F(p, t) is the evolution speed of the surface at the current position p and at the current time t, and is calculated by the right-hand side of (1). To calculate the mean curvature of the surface, we employ the discrete curvature estimator proposed by Desbrun et al. [15]: Here x j is one of the vertices at the one-neighborhood of x i . α j and β j are the two angles opposite to the edge connecting the two vertices x i and x j , and H is the mean curvature vector at vertex x i . In essence, (1) controls how each point in the deformable surface should move in order to minimize the weighted surface area. The detected object is then reconstructed by the steady-state solution of the equation S t = 0 (i.e., the velocity is zero).
In order to control both the smoothness of the model and the size of each triangle during the model-deformation phase, we must allow the model to be able to dynamically change its degrees of freedom during the deformation process. This is achieved using the local subdivision. If the area of an active face (on the model) is larger than a certain userdefined threshold, then this face will be subdivided into four smaller triangles by splitting the middle positions of its three edges.

Model relaxation
To ensure that the numerical simulation of the deformation process proceeds smoothly, we must maintain the regularity of the mesh such that the mesh has a good node distribution, a proper node density, and a good aspect ratio of the triangles. This is achieved by the incorporation of the tangential Laplacian operator, and three mesh operations: edge split, edge collapse, and edge swap. Laplacian operator, in its simplest form, moves repeatedly each mesh vertex by a displacement equal to a positive scale factor times the average of the neighboring vertices. Consider a mesh vertex P and its neighbors Q 1 , Λ, Q n , the Laplacian operator U is The tangential Laplacian operator is used to maintain a good node distribution and is defined as where n is the mesh normal at vertex P and C is a positive constant. Edge split and edge collapse are used to keep an appropriate node density. An edge split is triggered if the edge length Ye Duan 3 is bigger than the maximum edge length threshold. Similarly, an edge will be collapsed if its length is smaller than the minimum edge length threshold. Edge swapping is used to ensure a good aspect ratio of the triangles. As suggested by Kobbelt et al. [16], this can be achieved by forcing the average valence to be as close to 6 as possible. An edge is swapped if and only if the quantity p∈Δ (valence(p) − 6) 2 is minimized after the swapping.

Topology modification
In order to recover a shape of arbitrary, unknown topology, the model must be able to modify its topology properly whenever necessary. In general, there are two kinds of topology operations: (1) topology merging, and (2) topology splitting. We will explain these two operations in the following two subsections, respectively.

Topology merging
We propose a novel method called "lazy merging" to handle topology merging. The basic idea is that whenever two nonneighboring vertices are too close to each other, they will be deactivated (i.e., not allowed to move). Topology merging will happen only after the deformation of the model stops and all the vertices become nonactive. There are three steps in the topology merging operation: (1) collision detection, (2) merging-vertices clustering, and (3) multiple-contours stitching.

Collision detection
Collision detection is done hierarchically in two different levels: coarser-level and finer-level. Coarser-level collision detection is mainly for the purpose of collision exclusion. For each active vertex V , we will calculate its distance to all other nonneighboring active vertices. Vertices whose distance to the current vertex V is bigger enough so that no collision will happen between them will be excluded from collision detection in the next level. Otherwise, they will be passed to the finer-level collision detection. For each face f with three corner points (u, v, w) that is adjacent to one of the vertices being passed into the finer level of collision detection, we will calculate the distance between a number of sample points αu + βv + γw of the face f with barycentric coordinates α + β + γ = 1 and the current vertex V . If at least one of these distances is smaller than the collision threshold, the two corresponding vertices will be marked as merging vertices and will be deactivated. To further speed the performance, a uniform occupancy grid can be superimposed on the domain space for faster collision detection. Each vertex of the object will belong to a grid cell, and each grid cell will store the index/pointer of the vertices that belong to the current grid cell. At the beginning of each deformation step, the occupancy grid needs to update its vertices information. This can be done locally since only a few vertices will move at each deformation step, and it usually takes constant time, or at most O(n).

Merging-vertices clustering
After all the merging vertices have been deactivated, we need to divide them into several connected clusters. We randomly pick any merging vertex and find all of its connected merging vertices by a breadth-first search. We recursively do this until all the merging vertices belong to certain merging vertex clusters. Then for each cluster, we will remove all its interior vertices, and put all its boundary vertices into a linked list. We called this algorithm "lazy merging," and is based on the following observation: when two or more propagating fronts are merging with each other, only the boundary regions will remain, all the interior regions will disappear (i.e., removed).

Multiple-contours stitching
After the merging-vertex linked lists have been created, we need to stitch them together. We propose a "multiplecontours stitching" algorithm that is based on the proximity information between merging vertices. The algorithm iteratively connects each vertex in the linked list with its corresponding merging vertex in another linked list, and creates triangle strips that connect the two contours. If there are more than two contours that need to be stitched together, then there may be holes generated in the center of multiple contours. A hole-filling operation is conducted by inserting a new vertex in the center of the gap and connecting it with all the vertices in the gap. Figure 1 illustrates the three steps involved in the "multiple-contours stitching" algorithm.

Topology splitting
Topology splitting occurs when a part of the surface becomes too narrow and intends to shrink itself to a single point. In this scenario, the surface has to split up into two parts precisely at that location. In this paper, we use a method similar to the one proposed by Kass et al. [1]. In particular, if there exist three neighboring vertices which are interconnected to each other, but the face consisting of these three vertices does not belong to the model (i.e., a virtual face), then if the length of any of the three edges of the virtual face is smaller than the minimum edge length threshold and thus needs to be collapsed, a split operation is triggered. For example, in Figure 2, face ABC represents a virtual face that needs to be split because the length of the edge BC is smaller than the threshold. To properly handle the split operation, we divide the surface exactly at this location by cutting it into two open sub-surfaces. Then we close the two split-in-two surfaces using two faces A1B1C1 and A2C2B2 whose orientations are opposite to each other. Finally, we reorganize the neighborhood around the new faces A1B1C1 and A2B2C2, while removing the old vertices A, B, and C from the model.

Model refinement
Once an initial shape of the object is recovered, the model can be further refined several times to improve the fitting accuracy. In this paper, we have implemented two kinds of model refinement: global refinement and local/adaptive refinement. The decision of which method to be employed can be made either interactively by the user (whether he/she prefers a more uniformed mesh or an adaptively sampled mesh), or automatically by the system. The system can make a technically sound decision by calculating the variance of the fitting accuracy of the current model. If the variance of the fitting accuracy is very low, then the underlying object must be relatively smooth and global refinement will be a good choice. Otherwise, adaptive refinement will be used to recover the fine details embedded in the underlying Figure 2: Topology splitting by splitting the virtual face ABC into two faces A1B1C1 and A2C2B2 whose orientations are opposite to each other.
object. Global refinement is conducted by Loop's subdivision scheme [14]. Adaptive refinement is guided by the fitting accuracy. Various kinds of metrics can be used to evaluate the fitting accuracy over each triangle. For example, if the object boundary is well defined, then the fitting accuracy can be defined as the maximum or average distance of the triangle to the object boundary. However, for most of the medical image datasets including those obtained from image modalities such as CT or MRI, the object boundary is unknown; in this case, the curvature of the triangle can serve as the fitting accuracy. The curvature of the triangle can be defined as the summation of the mean curvature of its three vertices, which can be calculated by the aforementioned discrete curvature estimator of (4) in Section 2.2. At each level of adaptive refinement, all the triangles whose fitting accuracy is below the user-specified threshold will be quadrisected. The deformation of the model will resume only among those newly refined regions. Several levels of adaptive refinement can be applied until a user-specified fitting accuracy has been met.

RESULTS
In this section, we will show some experimental results obtained from our new model on both real and synthetic volumetric datasets in medical imaging applications. Note that in all the following figures, red regions represent parts of the model that are still active and deforming, blue region represent parts of the model that have already reached the boundary of the object and are not active anymore. The input dataset of Figure 3 is an MRI brain dataset of size 91 by 109 by 91 voxels. Figure 3(a) is the seed model (shown as green lines) enclosing the volumetric dataset (shown as cyan points). Figure 3(b) to Figure 3(d) are the three snapshots during the deformation process. Figure 3(e) is the wireframe view of the initial recovered shape. Figure 3(f) is the shape after one level of global refinement. Figure 3(g) is the shape after another level of adaptive refinement. Figure 3(h) is the rendered view of the same model of Figure 3(g). The input dataset of Figure 4 is obtained from CT scanning of a phantom of the vertebra. The data size is 128 by 120 by 52 voxels. Figure 4 of the model has been correctly modified. Figure 4(f) is the wireframe view of initial recovered shape of the model. Figure 4(g) is the refined shape after one level of global refinement. Figure 4(h) is the rendered view of the same shape in Figure 4(g). The input dataset of Figure 5 is a lobster dataset generated by CT scan with 128 by 128 by 128 voxels. Figure 5(a) shows the initial seed model (shown as tiny red dot) inside the volumetric dataset. This illustrates that the new model cannot only shrink from the outside of the dataset, but can also grow from the inside of the dataset. The middle two figures Figure 5(b) and 5(c) are the two snapshots of the model before and after the topology modification. Figure 5(d) is the wireframe view of initially recovered shape. Figure 5(f) is the wireframe view of the shape after one level of global refinement. Figure 5(g) is the shape after another level of adaptive refinement. Highly detailed features such as the tails are well represented. Figure 5(h) is the rendered view of the same model of Figure 5(g). The input of Figure 6 is a synthetic volumetric dataset generated by the union of two disjoint tori with 100 by 100 by 100 voxels. Just like previous examples, Figure 6(a) is the model initialization stage. Figure 6(b) to Figure 6(c) are the two snapshots of the model deformation process. The topology modification process is illustrated in Figure 6(d) to Figure 6(e). Figure 6(f) is the wireframe view of the initial recovered shape. Figure 6(g) is the refined shape after one level of global refinement. Figure 6(h) is the rendered view of the same model in Figure 6(g). Table 1 summarizes the information of the recovered shape, including the number of vertices, edges, and faces for each model, and the running time. The running time is measured on an Intel Pentium 4 M 1.6 GHz Notebook PC with 384 MB internal memory.

CONCLUSION
We have proposed a topology-adaptive multiresolution deformable surface that is very useful for medical image applications such as boundary reconstruction, surface visualization, data segmentation, and topology discovery. The deformation of the new model is governed by PDEs such as the weighted minimal surface flow. In comparison with existing topology-adaptive explicit models, our model is natural and Ye Duan 7 intuitive, easy to use and implement, hierarchically flexible, and has much better control of the model quality and resolution. Significantly different from the level-set approach, our new model always maintains the explicit representation of the geometry and topology, making it a powerful tool for topology interrogation and geometric recovery. More importantly, our model is uniquely suited for modeling timevarying datasets (because of its PDE relevance) such as 3D motion tracking for medical imaging analysis. We are currently exploring research topics along this direction.