A Simple and Stable Feature-Preserving Smoothing Method for Contours-Based Reconstructed Meshes Kai Xu*

Yueshan Xiong†

Yanzhen Wang‡

School of Computer Science, National University of Defense Technology, Changsha, 410073, China Ke Tan§ Guangyou Guo¶ Center for Technique Education, PLA General Hospital, Beijing, 100853, China Abstract In this paper, we develop a new feature preserving smoothing method for the irregular and coarse meshes reconstructed from 2D contours. To make the feature detecting robust, a new detecting algorithm using the continuity among adjacent contours is proposed. Then the original mesh is subdivided adaptively according to the detected geometric features and smoothed with isotropic method. With the help of that, our algorithm obtains not only the feature-preserving result of anisotropic methods but also the simple and stable process of isotropic ones. After smoothing, the resolution of mesh can be resumed by moving the vertices in the original mesh to the positions of their corresponding vertices in the smoothed adaptive mesh. For smoothing of contours-based reconstructed meshes, experimental results show that our method can smooth the meshes satisfactorily and preserve their geometric features much better than other implemented methods. CR Categories: I.3.0 [COMPUTER GRAPHICS]: General Keywords: mesh smoothing, feature preserving, feature detecting, adaptive subdivision

1

Introduction

Nowadays, various mesh models are obtained from different measurement techniques or reconstructed from 2D or 3D data sets. The derived models are often characterized by interesting features which represent the geometric features of the real objects. Unfortunately, these meshes are typically disturbed by some noises, which are due to the inevitable physical noises, human errors or some limitations of reconstruction algorithm. Mesh smoothing is a process dedicated to the removal of the This work is supported by the National Natural Science Foundation of China under Grant No.60371036. * e-mail: [email protected] † e-mail: [email protected] ‡ e-mail: [email protected] § e-mail: [email protected] ¶ e-mail: [email protected] Copyright © 2006 by the Association for Computing Machinery, Inc. Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from Permissions Dept, ACM Inc., fax +1 (212) 869-0481 or e-mail [email protected]. GRAPHITE 2006, Kuala Lumpur, Malaysia, November 29 – December 02, 2006. © 2006 ACM 1-59593-564-9/06/0011 $5.00

391

undesirable rough features. Its ultimate goal is to produce highly smooth mesh, for rendering, visualization and other post processing, while still preserve the overall shape and important geometric features of the original mesh. Mesh smoothing is mostly viewed as a diffusion process minimizing the constrained energy. From this point of view, the mesh smoothing algorithms can be classified into two classes: isotropic algorithms and anisotropic ones. Isotropic smoothing adopts isotropic diffusion process to minimize the constrained energy of a surface mesh. Laplacian smoothing is a well-known isotropic method, which can result in serious shrinking and degrading of meshes. Taubin [1995] proposed a low-pass filter by reducing the problem of mesh smoothing to low-pass filtering. Desbrun et al. [1999] extended this approach to irregular meshes using mean curvature flow which can overcome the tangent floating of vertices during smoothing, and introduced the implicit fairing to remove the stability condition, allowing for large time step to speed up the smoothing process. Other improvements followed, such as methods of combining geometry smoothing and parameterization regularization [Ohtake et al. 2000], and curvature averaging [Karbacher et al. 1998]. However, isotropic methods indiscriminately smooth noises and sharp features. Therefore the preservation of sharp features has to resort to a different approach. Feature-preserving mesh smoothing has been proposed more recently [Desbrun et al. 2000; Clarenz et al. 2000; Meyer et al. 2002; Zhang and Fiume 2002; Bajaj and Xu 2003], mostly inspired by image processing work on scale-space and anisotropic diffusion. The idea behind these approaches is to modify the diffusion equation to make it non-linear and/or anisotropic. The curvature tensor determines the local diffusion, thus preserving or even enhancing sharp features. Desbrun et al. [2000] applied the anisotropic curvature diffusion to smoothing of height fields. Clarenz et al. [2000] combined the nonlinear diffusion equations and the theory of geometric evolution problems and proposed an anisotropic smoothing algorithm for general 3D meshes. In both cases, geometric features are detected explicitly and locally at each time step, using either an edge contrast parameter or a principal curvature threshold. Zhang and Fiume [2002] presented MIN-DIST Flow algorithm through minimizing the vertex-to-plane distances to preserve features without explicit edge detecting. They still defined the bridges and corners implicitly by dihedral angle of adjacent triangles. Jones et al. [2003] proposed a non-iterative smoothing algorithm based on robust statistics and local first-order predictors of the surface. Fleishman et al. [2003] presented another anisotropic algorithm called bilateral mesh denoising motivated by the impressive results of bilateral filtering for image denoising. Feature detecting is critical to most anisotropic algorithms. It

often depends on the estimation of geometric properties of meshes, such as principal curvatures, principal directions etc., and consequently depends strictly on the topology connectivity and resolution of meshes. Contours-based 3D reconstruction is an important mesh generation method in visualization of medical data (CT, MRI, etc.). In this method, contours of interesting regions are extracted in every cross-section image with the help of some image processing techniques. 3D meshes are then reconstructed by connecting adjacent 2D contours. Since the adjusting of mesh density is hard for this method, it often provides more irregular and coarse meshes compared to other mesh generation methods. For these meshes, previous methods which simply depend on the estimation of geometric properties can hardly detect geometric features robustly. And more importantly, the previous anisotropic smoothing methods rely on shock formation to preserve details, which affects the numerical conditioning of the diffusion equations [Jones et al. 2003] and adds more instability to the smoothing of those irregular and coarse meshes. In this paper, we propose a new feature detecting approach aiming at detecting geometric features of meshes reconstructed from 2D contours more robustly, and develop a new smoothing method based on adaptive subdivision to break the limitation of numerical stability of most anisotropic methods. The contributions of our algorithm are as follows: Firstly, a robust feature detecting algorithm for contours-based reconstructed meshes is proposed by considering the continuity among the adjacent contours. Secondly, unlike other methods, we build our anisotropic method on subdividing a mesh adaptively according to the information of geometric features and smoothing the adaptive mesh with isotropic method. The isotropic diffusion on the adaptive mesh is equivalent to an anisotropic diffusion on the original mesh, which enables us to obtain the feature preserving result at the cost of linear isotropic diffusion process. And what’s more, the isotropic diffusion makes the problem of numerical stability, which is often a main problem in smoothing an irregular and coarse mesh, easy to deal with. Finally, the resolution of the mesh can be resumed by moving the vertices in the original mesh to the positions of their corresponding vertices in the smoothed adaptive mesh. The paper is organized as follows. In section 2 we briefly describe the general mesh smoothing and present our algorithm. Then we discuss the feature detecting and adaptive subdivision of our algorithm respectively in section 3 and 4. Some problems of implementation are discussed in section 5. Finally, we present some comparisons with the smoothing experimental results in section 6 to show the effectiveness of our algorithm.

2 2.1

Algorithm Notations

For

convenience of discussion, let M denote a mesh, vi = ( xi , yi , zi ) a vertex of mesh M , X the set of vertices of mesh M , and eij the edge (if existing) connecting vi to v j . We will call N1 (i ) the “1-ring neighbors” of vi , i.e., all vertices v j such that there exists an edge eij between vi and v j . Let # N1 (i ) denote the cardinal number of set N1 (i ) . Assume all the contours are organized in order and numbered sequentially, c(i) denotes the serial number of the contour the vertex vi is on. We have 0 < c(i) < # C , where #C is the total number of contours. We use N1c (i ) to denote 1-ring neighbors of vi on contour c(i) , and N1c −1 (i ) and N1c +1 (i ) 1-ring neighbors on c(i) −1 and c(i) + 1 respectively. We have N1 (i) = N1c−1 (i) ∪ N1c (i) ∪ N1c +1 (i) . In

addition, n denotes normal vector of a vertex, and m1 and m 2 are the two principal directions in which the normal curvature of a vertex takes on a minimum and maximum value, denoted by the scalars κ 1 and κ 2 respectively.

2.2

Mesh Smoothing

Most of mesh smoothing techniques use constrained energy minimization. The most used energy functions are membrane function and thin plate function:

E memb ( S ) =

1 ( X u2 + X v2 ) dS 2 ∫S

(1)

1 (2) ( X uu2 + 2 X uv2 + X vv2 ) dS 2 ∫S Their respective derivatives correspond to the Laplacian and the second Laplacian: E thin ( S ) =

L( X ) = X uu + X vv

(3)

L ( X ) = L o L( X ) = X uuuu + 2 X uuvv + X vvvv

(4)

2

One common way to attenuate noises in a mesh is through a diffusion process [Taubin 1995]:

∂X = λ L( X ) ∂t

(5)

In his λ|µ algorithm, Taubin[1995] used a combination of the two derivatives of the form:

∂X = (λ + µ ) L ( X ) − λµ L2 ( X ) ∂t

(6)

In the diffusion equation (5) and (6), the Laplacian operator determines the diffusion tensor, which represents the moving speed of vertices during the diffusion process. Both of the equations are isotropic and cannot preserve sharp features.

2.3

Smoothing Contours

Meshes

Reconstructed

from

The density of the mesh reconstructed from contours can only be adjusted through adjusting the space between cross sections (rate of vertical sampling) or number of vertices in one cross section (rate of horizontal sampling). Because of the difficulty of mesh density adjusting, the reconstructed meshes are often undesirably irregular and coarse, which leads to the inaccurate estimation of the geometric properties and hence the inaccurate feature detecting of those approaches which simply depend on the estimation of the geometric properties on discrete surfaces. By examining the continuity among the adjacent cross sections, we know that the true features often exhibit some continuity among the adjacent cross sections. A robust feature detecting approach can be derived from determining whether a vertex represents a noise (called noise vertex in this paper) or a geometric feature (called feature vertex) by considering the continuity among the adjacent contours. We establish our anisotropic smoothing algorithm based on the detected geometric features and the diffusion equation (5) or (6). We discretize the Laplacian operator to the well-known umbrella operator: 1 (7) U ( vi ) = ∑ v j − vi # N1 (i) j∈N1 (i ) As we use the umbrella operator, the moving speed of a vertex

392

As with some other geometric feature detecting methods [Clarenz et al. 2000; Zhang and Fiume 2002], we classify the geometric features of meshes into edges (or ridges) and corners. Edges often appear as a feature line on the mesh. On a feature line, the curvatures of adjacent contours are close (including magnitude and direction). Also, the general direction of a feature line is close to one of the two principal directions of a vertex on the feature line. In our feature detecting algorithm, the principal directions of vertices are estimated, with which we compare the curvatures of adjacent contours along the principal direction to detect geometric features. As shown in Figure 2(a), along the direction m 2 (one of the two principal directions of vertex vi ), curvatures of contour c(i) ( vi is on it) and its two adjacent contours c(i) −1 and c(i) +1 are close (including magnitude and direction), i.e. the vertex vi exhibits the local continuity and is unlikely a result of random noise. Instead, we expect the vertex to reflect a true geometric feature of the original mesh. We call it an edge vertex in this paper. On the contrary, as shown in Figure 2(b), the above-mentioned curvatures are significantly discontinuous. We expect vi to reflect a random, isolated perturbation, or in other word, a noise.

is proportional to the distance from the vertex to the barycentre of its 1-ring neighbors. The noise vertices often have a big value of U (vi ) and move rapidly during smoothing, which exactly meets our desire for smoothing. But the vertices representing sharp features also have a big U (vi ) , and move too rapidly to attenuate undesirably during smoothing. As a result, we reduce the moving speed of a feature vertex by reducing its U (vi ) value. Specifically, we subdivide the 1-ring neighboring triangles of a feature vertex if its U (vi ) value exceeds a threshold. By subdividing the mesh adaptively and smoothing it with isotropic diffusion, the feature vertices may move more slowly, which contributes to the feature preserving, while the noise vertices still move rapidly, which is beneficial to a fast smoothing. The isotropic diffusion on the adaptive mesh is equivalent to an anisotropic diffusion on the original mesh, so we obtain the feature-preserving result of anisotropic method at the linear cost of isotropic one. Furthermore, the stability conditions of the isotropic diffusion equation (5) and (6) can be satisfied easily (see [Desbrun et al. 1999]). The resolution of the smoothed mesh can be resumed by moving the vertices in the original mesh to the positions of their corresponding vertices in the smoothed adaptive mesh. Figure 1 shows the process of our algorithm through a 2D curve smoothing.

Figure 2: The vertex vi in (a) is likely to represent a geometric feature. The vertex vi in (b) is likely to represent a noise. We need to estimate the curvatures of contour c(i) −1, c(i) , and c(i) +1 along the principal direction m 2 of vi . Let kic denote the three curvatures mentioned-above, where c = −1, 0,1 . As for ki0 , we can simply use the curvature of contour c(i) at vertex vi estimated with vi and the vertices in N1c (i ) . We can estimate ki− 1 and k i1 as follows. Here we take the ki−1 for example. As shown in Figure 3, we use a conic c to fit the vertices in N1c −1 (i) . S is a plane passing through vertex vi with m1 as its normal vector. Suppose curve c intersects plane S at point P. The ki− 1 can be approximated by the curvature of c at P.

Figure 1: In original mesh (a), vertex 2 is detected as a feature vertex while 5 a noise one. The mesh is adaptively subdivided (b) based on the information of geometric features and smoothed (c) by isotropic diffusion process. At last, the resolution of the mesh was resumed in (d). Our algorithm is described in Algorithm 1. Algorithm 1: Mesh smoothing based on feature detecting and adaptively subdivision 1. Detect geometric features on the original mesh; 2. Subdivide the original mesh adaptively with the information of detected geometric features; 3. Smooth the adaptive mesh with the isotropic smoothing method; 4. Resume the resolution of the mesh;

3

Geometric feature Detecting

Noises in contours-based reconstructed meshes are mainly from the image noises, which is hard to be filtered, and the inaccuracy of image segmentation. Thus, a noise often appears as a perturbation on a single contour while a geometric feature often exhibits some continuity on several adjacent contours.

Figure 3: Estimate the curvature of contours c(i) −1 along the principal direction m 2 of vertex vi

393

The inaccuracy caused by local perturbations of a contour can be minimized by curve fitting. Additionally, if all the vertices in N1c (i ) lie on the same side of plane S and hence the intersection point P doses not exist, we would consider the vertex as a noise. If vi is an umbilicus vertex, which means its principal directions are zero vectors, the vertex may be either a feature vertex (see Figure 4(a)) (called a corner vertex in this paper.) or a noise vertex (see Figure 4(b)). In practice, a vertex can be treated as an umbilicus one if its two principal curvatures are close enough.

Figure 5: Local reference frame at vertex vi We transform the coordinates of 1-ring neighbors v j ∈ N 1 (i ) of vertex vi in world reference frame (denoted by ( xi , yi , zi ) ) to the coordinates in the local reference frame (denoted by (ui , vi , hi ) ). Similar to Milroy et al., we used a parabolic quadric surface S ( u , v ) = ( u , v , h ( u , v )) (where h(u , v) = a 0 u 2 + a1uv + a 2 v 2 ) to fit the 1-ring neighbors by weighted least-squares surface fitting, in which we minimize the weighted sum defined in (9):

Figure 4: The vertex vi in (a) is likely to represent a geometric feature. The vertex vi in (b) is likely to represent a noise. We can determine if an umbilicus vertex vi is a corner or a noise vertex with its 1-ring neighbors v j through equation (8):

Ψi =

∑w

j

⋅τ j

H (a 0 , a1 , a 2 ) =

∑ j =1

(8)

  u j  uj     ρj  vj  − v j  a u2 + a u v + a v2  h  1 j j 2 j   0 j  j

2

(9)

To fit a polynomial that best approximates the surface at vertex vi , closer 1-ring neighboring vertices are weighted more heavily, as they best represent the local shape. We choose −d / d ρ j = 10 i , j max , where d i , j and d max have the same meaning as before. We can find the solution of a i ( i = 0 ,1, 2 ) which minimize the H by solving a 3×3 linear system. After that, the principal directions can be estimated with (10) and (11), and the principal curvatures with (12) and (13):

j∈N1 ( i )

In (8), the possible value of τ j may be -1, 0, or 1. If the neighboring vertex v j is detected as a noise vertex, the τ j is set to -1. And contrarily, if v j is detected as a feature vertex, the τ j is set to 1. Otherwise, if v j is also an umbilicus vertex and hasn’t been detected, we set the τ j to 0. The w j is a weight by which closer neighboring vertices are weighted more heavily as they have more influences on vertex vi . In our algorithm, we use −d /d w j = e i , j m ax , where d i , j is the Euclidean distance from vertex vi to v j and dmax = max{di , j , i = 1, 2,..., # N1 (i )} . If the Ψ i of an umbilicus vertex vi exceeds 0, we consider it as a corner vertex, otherwise, a noise vertex. We discuss the estimation of principal directions and principal curvatures of a vertex on a discrete surface in section 3.1, and then describe our geometric feature detecting algorithm in section 3.2.

3.1

# N1 ( i )

(a − a + (a − a ) 2 + a 2 , −a ) 2 0 0 2 1 1 m1 =  2 2  (a1 , a2 − a0 − (a0 − a2 ) + a1 )

Estimation of Principal Directions and Principal Curvatures

Milroy et al. [1997] used weighted least-squares surface fitting to estimate principal directions and principal curvatures of a point in OCS model (Orthogonal cross-section model) in segmentation of OCS model. We adopt his method in estimating of principal directions and principal curvatures of vertices in a triangulated mesh. We first build a local reference frame vi − huv at vertex vi (see Figure 5), where vi is at the origin, and the h axis is aligned with the vertex normal n . The u and v axes lie in the tangent plane at vertex vi . The principal directions m1 and m 2 are orthogonal mutually and lie in the uv plane.

a0 < a2 a0 ≥ a2

(10)

 ( a , a − a + (a − a ) 2 + a 2 ) a < a 1 2 0 0 2 1 0 2 m2 =  2 2 a ≥ a 2 (a2 − a0 − (a0 − a2 ) + a1 , − a1 ) 0

(11)

κ 1 = a 0 + a 2 − (a 0 − a 2 ) 2 + a12

(12)

κ 2 = a0 + a 2 + (a0 − a 2 ) 2 + a12

(13)

If a 0 = a 2 , or | κ 1 − κ 2 |< ε , we set the vertex vi an umbilicus vertex, whose principal directions are both zero vectors.

3.2

Feature Detecting Algorithm

As described above, the ultimate goal of feature detecting is to reduce the moving speed of feature vertices. Because of the umbrella operator used in our isotropic diffusion, which determines the moving speed, we only need to examine those vertices whose value of umbrella operator exceeds the threshold γ . The detecting algorithm is described in Algorithm 2.

394

Algorithm 2: Detecting the geometric feature for all vertices vi in the mesh do begin if U (vi ) > γ then if vi is a umbilicus vertex then calculate Ψ i for vertex vi ; if Ψ i > 0 then mark vi as a corner vertex; else mark vi as a noise vertex; else calculate kic ( c = − 1, 0,1) for vertex vi ; if ki0 ⋅ k i− 1 > 0 && | k i0 − ki− 1 |< δ && ki0 ⋅ k i1 > 0 && | k i0 − k i1 |< δ then mark vi as an edge vertex; else mark vi as noise vertex; end

Figure 6: Feature detecting on a noisy skull model (left top) and a noisy pelvis model (left bottom). Feature vertices are designated by green points and noise vertices by red points (right).

Both the γ and the δ should be chosen for different meshes. The γ can be chosen by calculating the average value of the umbrella operator of all the vertices, which is denoted by U = mean {U (vi ), vi ∈ X } . Our choice is 2 ⋅ U . The larger δ is set, the more vertices will be accepted as feature vertices. In practice, we can find out the best δ for individual mesh by experiments. Figure 6 shows the feature detecting results on noisy skull model and pelvis model by our algorithm. Both of the models are reconstructed from contours extracted from CT data.

4

Adaptive

Subdivision

Based

on

Feature

Information Figure 7: (a) is original mesh, assuming that only vertex 1 is detected as a feature vertex. In order to preserve this feature, we directly modify the moving speed of the feature vertex with equation (12) in (b), which results in serious tangential drifting of the neighboring vertices 2 and 3 after 10 iterations. In (c), instead, we subdivide the neighboring segments of vertex 1 three times. The drifting of vertices 2 and 3 is obviously reduced after the same times of iteration and the shape is preserved better.

As discussed above, our anisotropic diffusion is based on the adaptive subdivision of the original mesh. The algorithm reduces the moving speed of a feature vertex simply by subdividing its 1-ring triangles continually until its value of umbrella operator is reduced below the threshold γ . The times of subdivision needed in practice are generally not more than 3. We find out that if we reduce the moving speed of feature vertices through modifying the diffusion equation (5) or (6) directly, more tangential drifting will be introduced, which is harm to keeping the shape of a mesh [Desbrun et al. 1999]. We demonstrate this phenomenon through a smoothing experiment on a 2D close curve, shown in Figure 7. In the experiment, we use both subdivision scheme and direct modifying scheme to smooth the noisy 2D close curve. In subdivision scheme, we subdivide the neighboring segments of the detected feature vertices by 3 times, and use equation (6) to smooth the subdivided adaptive mesh. Accordingly, we use (12) for the feature vertices and equation (6) for other vertices in direct modifying scheme. Smoothing is performed the same times of iteration in the two schemes.

∂X = 2 − 3 ( λ + µ ) L ( X ) − 2 − 2×3 λµ L2 ( X ) ∂t

As a result, we choose the subdivision scheme in our anisotropic algorithm.

5

Implementation

Serious noises scattered in mesh often make the feature detecting more difficult, in which case low-pass pre-filtering described in [Clarenz et al. 2000] is necessary in order to separate the true features from noises effectively. For the reason of computational accuracy, some applications, such as finite element analysis and deformable modeling based on Finite Element Methods (FEMs), usually expect the models to have denser mesh in regions with more details than others. In other words, a mesh should be adaptive according to the local curvature of its surface. The adaptive subdivision of our algorithm

(12)

395

of the model surface are very important. As a result, we do our smoothing experiments on bone models reconstructed from 2D contours extracted from CT data. We implemented isotropic methods including Laplacian smoothing, Taubin’s λ|µ smoothing [Taubin 1995], and the mean curvature flow (with the explicit integrating scheme) [Desbrun et al. 1995], anisotropic method the MIN-DIST flow [Zhang and Fiume 2002], and the algorithm presented in this paper. From Figure 9 to 12, we uniformly arrange the results as follows. (a) is original mesh, (b) result by Laplacian smoothing (λ = 0.3, 20 iterations), (c) result by Taubin’s λ|µ smoothing (λ = 0.6307, µ = -0.6732, 20 iterations), (e) result by the mean curvature flow (explicit integrating, λ = 0.3, 20 iterations), (e) result by the MIN-DIST flow (λ = 0.5, β = 0.03, 20 iterations), and, and (f) result by our method (use equation (6) in the isotropic smoothing, λ = 0.6307, µ = -0.6732, 40 iterations). All models are flat-shaded in rendering to enhance the faceting effect.

exactly fills the gap between the needs of those applications and the difficulty of mesh density controlling of contour-based mesh reconstruction methods. If the user hopes that the smoothed mesh has the same resolution as the original mesh, our algorithm can easily resume the resolution as follows. Let X 0 denote the vertex set of the original mesh and X a( s ) the vertex set of the subdivided and ( s) , where X 0( s ) is the smoothed mesh. We have X a( s ) = X 0( s ) + X new (s) set of original vertices while X new is the set of new vertices generated in the subdividing. The result mesh with original resolution is composed of the vertices in set X 0( s ) and their connectivity in the original mesh.

6

Experimental results

In visualization of human bones, both the smoothness and details

Figure 9: Smoothing the skull (super maxilla part) model

Figure 10: Smoothing the pelvis model

396

Figure 11: Smoothing the femur (knee joint part) model

Figure 12: Smoothing the femur (hip joint part) model Figure 9 shows the smoothing results of skull (super maxilla part) model. As we can see, Laplacian smoothing shrinks the model seriously. Taubin’s λ|µ smoothing alleviates the shrinking remarkably. Mean curvature flow can preserves some features by fixing up the tangential drifting, but it makes some features too sharp at the same time. Generally speaking, isotropic methods can not preserve the geometric features. MIN-DIST flow preserves most of the features of original mesh, but it enlarges the eye socket and widens the nose bridge simultaneously because of the inaccurate feature detecting. Our algorithm, instead, preserves the overall shape and geometric features more effectively. Figure 10 compares the five smoothing methods on pelvis model. As with former results, isotropic methods attenuate most

geometric features, such as the three features designated in Figure 10(a). Compared with other algorithms including the anisotropic MIN-DIST flow, our algorithm can preserve those important features availably thanks to the robust feature detecting and stable smoothing process. The similar results can be observed in Figure 11 and Figure 12, in which the knee joint part of femur model and hip joint part of femur model are smoothed respectively. As for smoothing of meshes reconstructed from 2D contours, the experimental results shows the advantage of our algorithm in feature preserving compared to the other 4 algorithms. Table 1 compares the execution time of the algorithms on PC (1.7GHz CPU and 512MB memory). Our algorithm

Laplacian

λ|µ

Mean curvature flow

MIN-DIST flow

detecting

subdividing

smoothing

total

skull

0.031

0.063

0.250

0.738

0.009

0.478

0.141

0.628

pelvis

0.031

0.079

0.282

0.871

0.012

0.603

0.159

0.774

femur (knee joint)

0.000

0.015

0.047

0.098

0.000

0.048

0.045

0.093

femur (hip joint)

0.016

0.047

0.157

0.359

0.000

0.051

0.047

0.098

Table 1: Comparison of execution time of the 5 algorithms As the Table 1 reported, the anisotropic algorithms generally spend more time than the isotropic ones. Our algorithm shows computational advantage over the MIN-DIST flow algorithm. Especially for models with fewer features, such as femur (hip joint part) model, lots of execution time can be saved by our algorithm by avoiding a great deal of computations of implicit feature detecting in every iteration needed by MIN-DIST flow

algorithm in virtue of the explicit and robust feature detecting before the smoothing iteration. In Figure 13, we plot the graphs of volume shrinking rate as a function of times of iteration. As shown in the graphs, our algorithm has the least volume changing among the 5 algorithms. (The negative volume shrinking means swelling of a mesh.)

397

Figure 13: Times of Iteration/Volume Shrinking Rate curves of the 5 smoothing algorithms

7

Press / ACM SIGGRAPH, New York. A. Rockwood, Ed., Computer Graphics Proceedings, Annual Conference Series, ACM, 317-324.

Conclusion

Since the contour-based 3D reconstruction is an important mesh generation method in visualization of medical data and geometric modeling for surgery simulation, smoothing of meshes reconstructed from 2D contours becomes an important work. Meanwhile, it often presents special difficulty. Accurate feature detecting on meshes reconstructed from 2D contours is usually difficult for previous anisotropic smoothing algorithms. We propose a new robust feature detecting approach for these meshes by considering the continuity among adjacent contours. With the detected geometric features, we develop a new smoothing algorithm based on adaptively subdivision. The numerical stability of our algorithm can be guaranteed since we use a linear and isotropic diffusion to smooth the adaptive mesh. Smoothing of the up-sampled mesh definitely needs more iterations than that of the original mesh to achieve the same smoothness. This is the main deficiency comes from our method. However, this will not introduce too much computation since the cost of adaptively subdividing a coarse mesh and smoothing the adaptive mesh with the isotropic method is tolerable. As for future work, we plan to extend our adaptive subdivision scheme to arbitrary coarse and irregular meshes because of its advantages in computational cost and numerical stability by combining a robust feature detecting algorithm for arbitrary coarse and irregular meshes.

DESBRUN, M., MEYER, M., SCHRÖDER, P., AND BARR, A. H. 2000. Anisotropic Feature-preserving Denoising of Height Fields and Bivariate Data. In Proc. Graphics Interface 2000, 145-152. FLEISHMAN, S., DRORI, I., AND COHEN-OR, D. 2003. Bilateral Mesh Denoising. ACM Transactions on Graphics, 22, 3, 950-953. JONES, T. R., DURAND, F., AND DESBRUN, M. 2003. Non-iterative, Feature-preserving Mesh Smoothing. ACM Transactions on Graphics, 22, 3, 943-949. KARBACHER, S., AND HAUSLER, G. 1998. A New Approach for Modeling and Smoothing of Scattered 3D Data. In Three-Dimensional Image Capture and Applications, Ellson R. N., Nurre J. H. (eds.), Washington: SPIE Press, 168-177. MEYER, M., DESBRUN, M., SCHRÖDER, P., AND BARR, A. H. 2002. Discrete Differential-geometry Operators for Triangulated 2-manifolds. In Proceeding of Visualization and Mathematics III, Hege H C, Polthier K (eds.), 35-57. MILROY, M. J., BRADLEY, C., AND VICKERS, G. W. 1997. Segmentation of a Wrap-around Model Using an Active Contour. Computer Aided Design, 29, 4, 299-320. OHTAKE, Y., BELYAEV, A., AND BOGAESKI, I. 2000. Polyhedral Surface Smoothing with Simultaneous Mesh Regularization. In Proc. Geometric Modeling and Processing 2000, 229-237.

References BAJAJ, C, AND XU, G. 2003. Anisotropic Diffusion on Surfaces and Functions on Surfaces. ACM Transactions on Graphics, 22, 1, 4-32. CLARENZ, U., DIEWALD, U., AND RUMPF, M. 2000. Anisotropic Geometric Diffusion in Surface Processing. In Proc. IEEE Visualization 2000, 397-405.

TAUBIN, G. 1995. A Signal Processing Approach to Fair Surface Design. In Proceedings of ACM SIGGRAPH 95, ACM Press / ACM SIGGRAPH, Los Angeles. S. Mair, R. Cook, Ed., Computer Graphics Proceedings, Annual Conference Series, ACM, 351-358.

DESBRUN, M., MEYER M., SCHRÖDER P., AND BARR A. H. 1999. Implicit Fairing of Irregular Meshes Using Diffusion and Curvature Flow. In Proceedings of ACM SIGGRAPH 99, ACM

ZHANG, H., AND FIUME, E. L. 2002. Mesh Smoothing with Shape or Feature Preservation. In Advances in Modeling, Animation, and Rendering, Vince J, Earnshaw R (eds.), Springer, 167-181.

398

A Mass-Spring Model for Surface Mesh Deformation Based on Shape Matching Yanzhen Wang, Yueshan Xiong, Kai Xu, Ke Tan, Guangyou Guo

Figure 1: Scenario of human chest simulated by the mass-spring model proposed in this paper. (a) Undeformed chest organs. (b) Deformation of heart model. (c) Deformation of left lung model.

A Simple and Stable Feature-Preserving Smoothing Method for Contours-based Reconstructed Meshes Kai Xu, Yueshan Xiong, Yanzhen Wang, Ke Tan, Guangyou Guo

Figure 6: Feature detecting on a noisy skull model (left top) and a noisy pelvis model (left bottom). Feature vertices are designated by green points and noise vertices by red points (right).

502