Surface mesh denoising with normal tensor framework

Graphical Models 74 (2012) 130–139

Contents lists available at SciVerse ScienceDirect

Graphical Models journal homepage: www.elsevier.com/locate/gmod

Surface mesh denoising with normal tensor framework Shoichi Tsuchie a,⇑, Masatake Higashi b a b

Nihon Unisys, Ltd., Japan Toyota Technological Institute, Japan

a r t i c l e

i n f o

Article history: Received 2 March 2012 Accepted 30 March 2012 Available online 12 April 2012 Keywords: Mesh denoising Normal tensor Laplacian of curvature Feature-preserving

a b s t r a c t In this paper, we propose a novel method for feature-preserving mesh denoising based on the normal tensor framework. We utilize the normal tensor voting directly for the mesh denoising whose eigenvalues and eigenvectors are used for detecting saliency, and introduce an algorithm that updates a vertex by the Laplacian of curvature which minimizes a difference of the curvature in one neighborhood. By connecting the feature saliency with a distance metric in the normal tensor space, our algorithm preserves sharp features more robustly and clearly for noisy mesh data. Comparing our method with the existing ones, we demonstrate the effectiveness of our algorithm against some synthetic noisy data and realworld scanned data. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Surface mesh denoising has been an important research area and its goal is to eliminate noise or spurious information on the mesh while preserving original features. Therefore, a good denoising algorithm should be able to remove the noise while retaining original features as clearly as possible. In particular, to preserve sharp features such as crease edges and corners indicated in Fig. 1, is our great concern, since they are often blurred if no special care is taken. Various approaches have still been proposed to tackle this problems on end, and are classified into three categories, depending on whether the process is performed under a deterministic, image filtering [7,17,34] or probabilistic point of view [6,23]. In this paper, we focus on a deterministic framework which is governed by an anisotropic diffusion equation. Most of their approaches use a weighed Laplacian smoothing. Representative works are Taubin’s signal processing technique called by kjl method [30] and Desbrun’s mean curvature flow method [5]. But the Laplacian-based methods have two main drawbacks:

⇑ Corresponding author. E-mail address: [email protected] (S. Tsuchie). 1524-0703/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.gmod.2012.03.010

(i) smoothing side-effects which blur sharp features, (ii) shrinkage of the shape as a whole. In order to avoid (i), anisotropic diffusion is discussed by many researchers [1,3,4,8,35]. The basic mechanism is the directional smoothing that at the smaller principal curvature value errors are reduced and at the larger value less. With respect to (ii), Kobbelt et al. [12] applied a surface fairing in CAGD [9], and proposed a discrete fairing method by the second order Laplacian in order to prevent from shrinking in denoising process. But their approach is not enough for mesh denoising with sharp features. In addition, we add the following problem to the above two: (iii) cumbersome setting of tuning parameters. For example, in the algorithms using a Gaussian weighting function, the band width becomes an important parameter for getting better results, but the way how to decide the value is hardly discussed in detail and treated as a thing that should be properly set by users suited for the level of noise. In this paper, we propose a novel method for the feature-preserving mesh denoising based on the normal tensor which is constructed by a sum of covariance matrices of unit facet normal vectors in one-ring neighborhood. The main idea behind our approach is followed by two things. At first, we recognize the characteristic features

131

S. Tsuchie, M. Higashi / Graphical Models 74 (2012) 130–139

KV(i) and the other is facet indices KF(i). @F denotes a set of edges that constitute a facet boundary. An illustration is shown in Fig. 2. 2. Related work

Fig. 1. Smoothing and denoising a noisy octahedron (jVj = 1026,jFj = 2048). From left to right, Noisy data, laplacian smoothing and our denoising.

The governing equation in mesh denoising where we focus on as a deterministic framework, is described by the following diffusion equation:



@uðx; tÞ ¼ div ðDðxÞruðx; tÞÞ

ð1Þ

uðx; 0Þ ¼ f ðxÞ; via surface normal not via curvature or higher order differential invariants [14], since they are sensitive to noise and cannot be directly computed on noisy data. Furthermore, sharp features cannot be primarily expressed by a differentiable manifold, so we cannot represent them precisely as a smooth surface. As reported in [2,28], the surface normal describes the feature saliencies well though it is lower order derivative than curvature. Secondly, in order to prevent from the shrinking problem, we also use a surface fairing technique as Kobbelt et al. [12], and apply it to the feature-preserving mesh denoising, though the existing methods does not often use it. In general, different from Laplacian based smoothing, fairing is related to aesthetics and is often processed with higher order derivative such as the second order Laplacian. Therefore, by the same reason as stated above, we develop a robust algorithm for noise to realize the fairing technique. Our contributions are as follows: by connecting the feature saliency based on the normal tensor framework, which is proposed by Medioni et al. [19] and Page et al. [24], with a distance metric in the normal tensor space, we have realized the following items:  a novel algorithm that updates a vertex by the Laplacian of curvature which minimizes a difference of the curvature in one neighborhood, preserves sharp features more robustly and clearly for noisy mesh data,  the proposed method prevents from the shrinkage problem due to the property of Laplacian of curvature,  our weighting function is simply designed, and the denoising algorithm contains only one tuning parameter that defines a crease angle. The rest of the paper is organized as follows: Section 2 presents a brief overview of related works. The working principle of normal tensor is explained in Section 3. In Section 4, we introduce our method which simultaneously smoothes the geometric flow and noisy mesh. Results and some comparison with other schemes are done in Section 5 and finally we conclude the paper in Section 6. Before proceeding to Section 2, we define some notations used in this paper. Notations: We treat a triangle mesh denoted by M(V, F), where V and F indicate mesh vertices V = {Vij i = 1, 2, . . . , jVj} and triangle facets F = {Fiji = 1, 2, . . . , jFj}, respectively. Here j  j denotes the cardinality of a set. Vertices and facets have each unit normal vector Nv and Nf respectively. And 1-ring neighborhood of a vertex Vi has two cases; one is vertex neighbor indices denoted by

where D(x) is diffusion coefficient, u(x, t) is a certain function such as concentration of a matter and x 2 Rd, t > 0. We show some mathematical foundations for the functions related to our proposed approach below. 2.1. Isotropic diffusion Here we consider d = 1, and D(x) is a constant D0. Eq. (1) is reduced to @ tu(x, t) = D0@ xxu(x, t), which results in the following equation with a simple finite difference scheme:

uðx; t þ DtÞ ¼ uðx; tÞ þ

D0 dt X ðDxÞ2 ~x2Kx

ðuð~x; tÞ  uðx; tÞÞ;

where Kx  {x + Dx,x  Dx}. We apply the above finite difference equation to an unstructured mesh, that is, replace Kx and u(x,t) with KV(i) and mesh vertex Vi respectively. Thus we obtain the following vertex update equation:

V nþ1 i

k V ni þ X

X wij j2KV ðiÞ

wij ðV nj  V ni Þ;

ð2Þ

j2KV ðiÞ

where k > 0 is the iteration step size, which should be a small value for stable calculation. Eq. (2) is well known as Laplacian smoothing, and the weight function wij plays an important role, since simple functions such as angularor area-weight have lost and blurred the feature characteristics due to a diffusion process. On the other hand, the counterpart of the Euclidian Laplacian D on a smooth surface is the Laplace–Beltrami operator DM. Thus we obtain a geometric diffusion equation [27]

@ t x ¼ DM x;

ð3Þ

where x is a point on the manifold M. From differential geometry [33], we know that the mean curvature vector KHn equals the Laplace–Beltrami operator on a surface manifold:

K H n ¼ DM x: Thus the geometric diffusion Eq. (3) is equivalent to the mean curvature flow (MCF) given by

@ t x ¼ K H n: In mesh area, Desbrun et al. [5] proposed the MCF smoothing method using a discrete mean curvature vector K Hi ni given by cotangent formula:

132

S. Tsuchie, M. Higashi / Graphical Models 74 (2012) 130–139

Since in general the above system of equation has no nontrivial solutions, Taubin [31] proposed to solve it in a least squares sense, i.e. to minimize the following error function:

WðVÞ ¼

X X

e f  ðV i  V j ÞÞ2 : ðN k

k2F ði;jÞ2@F k

Using the gradient descent method accompanied by a certain weighting function, the vertex update equation can be written as follows:

Fig. 2. Illustration of variables used in this paper.

K Hi ni ¼

1 X ðcot aj þ cotbj ÞðV nj  V ni Þ; 4S j2K ðiÞ

ð4Þ

k V ni þ X

V nþ1 i

V

2.2. Anisotropic diffusion Anisotropic diffusion process is a powerful denoising approach that performs the directional smoothing and enhancing of the manifold features while preserving adequately curvatures, edges and corners, and has led a great success to a 2D image denoising, originally in Perona and Malik [25]. They proposed to use a suitable, monotone decreasing diffusion coefficient D(x) = g(jru(x)j) instead of a constant in Eq. (1), so that the diffusion decreases across the edge, which is characterized by large jru(x)j. Afterwards, Weickert [32] and others have researched a type of Perona–Malik (P–M) diffusion equation. In mesh area, Clarentz et al. [3,4] and Bajaj and Xu [1] discussed the diffusion process from the viewpoint of differential geometry in detail. Hildebrandt and Polthier [8] proposed an anisotropic MCF and restored the sharp feature more clearly, which is given by the following scheme:

3k þ H A ðV ni Þ; S 1 X wðHe ÞHe N e ; H A ðV i Þ  2 e¼ij;j2K ðiÞ V ni

wk j2KV ðiÞði;jÞ2@F k

  ef N e f  ðV n  V n Þ : wk N j i k k

k2KF ðiÞ

where S is the sum of the areas of the triangles having Vi as a common vertex. The vertex position is updated along the surface normal, and the important features are preserved by a weighting parameter which depends on the principal curvature values.

V nþ1 i

X X

ð5Þ

V

where He and Ne are the mean curvature at the edge and the edge normal. Recently, Zhang et al. [35] proposed the scheme based on the P–M diffusion equation and reported a great success compared with the existing methods.

ð6Þ Taubin [31] lets the weight be constant. As indicated in [28], Ohtake et al. [22] proposed the different error function, but their result becomes to the same as Eq. (6) with the triangle area weighting function. Here we note that using the covariance matrix ef N ef T ~k ¼ N a k k , Eq. (6) can be rewritten in the following:

k V ni þ X

V nþ1 i

X X wk j2KV ðiÞði;jÞ2@F k

~k ðV nj  V ni Þ: wk a

k2KF ðiÞ

The result can be interpreted as the Laplacian smoothing with the weight defined by wk e a k , which is the essence of this paper discussed in Section 4 in detail. 3. Basic framework 3.1. Normal tensor We construct a normal tensor Ai as a weighted sum of a covariance matrix ak of the facet normals N fk ðk 2 KF ðiÞÞ connecting to a mesh vertex Vi:

X

Ai 

k2KF ðiÞ

wk ak ¼

X

wk N fk N fk T ;

ð7Þ

k2KF ðiÞ

where we use the Nelson Max’s weighting method [18], that is, the weight is an area of facet divided by the squares of the lengths of the two edges that connect vertex as follows:

wk 

kðV j  V i Þ  ðV jþ1  V i Þk 2kV j  V i k2 kV jþ1  V i k2

¼

sinh : kV j  V i kkV jþ1  V i k

Since the matrix Ai is a 3  3 and symmetric semi-definite, eigen decomposition generates real eigenvalues

2.3. Face normal field

r1 P r2 P r3 P 0

We assume that noisy facet normals are filtered beforehand by some methods, such as bilateral filtering [16,36], robust statistics [11], and diffusion process [22]. Once the e f (k = 1, 2, . . . , jFj) are obtained, the corfiltered normals N k responding triangle face vertices, V ki ði ¼ 1; 2; 3Þ, are updated by the family of simultaneous linear equations:

with the corresponding eigenvectors E1, E2, and E3 which are mutually orthogonal, Ei  Ej = di,j. As a result, we can rewrite Eq. (7) in the following:

8 ef > > < N k  ðV k2  V k1 Þ ¼ 0; e f  ðV k  V k Þ ¼ 0; N k 3 2 > > : ef N k  ðV k1  V k3 Þ ¼ 0:

0

r1

1

0

0

B Ai ¼ RRRT ¼ ðE1 E2 E3 ÞB @ 0

r2

C T 0 C AðE1 E2 E3 Þ :

0

0

ð8Þ

r3 :

Eq. (8) means that normal tensor Ai can be factored as the production of a rotation matrix R and a diagonal scaling matrix R.

133

S. Tsuchie, M. Higashi / Graphical Models 74 (2012) 130–139

3.2. Feature extraction Medioni et al. [19] defined saliency maps from their normal tensors called by normal tensor voting. The saliency maps use the following relationships for their tensor eigenvalues:

8 > < SSu ¼ r1  r2 ; surface saliency; SCr ¼ r2  r3 ; crease saliency; > : SCo ¼ r3 ; corner saliency:

ð9Þ

Fig. 3 shows the characteristics of eigenvalues of our normal tensor defined by Eq. (7), and we can verify that the eigenvalues satisfy the relationships in Eq. (9). Thereafter, Page et al. [24] suggested the maximum of Eq. (9) determines how we classify the feature saliency at a mesh vertex Vi as follows:

8 v > < SSu : surface N ¼ E1 v maxfSSu ; SCr ; SCo g ¼ SCr : crease T ¼ E3 > : SCo : no orientation;

on a surface mesh (Section 4.3) and apply it to our mesh denoising scheme (Section 4.4). 4.1. Normal tensor metric We define the following vectors

ð10Þ

Eei ¼

where  P 0 is a constant parameter which controls the relative significance of the above feature saliencies as shown in Fig. 4. If we set  ? 0, it always classifies a mesh vertex as a surface. (Conversely  ? 1 will never recognize a vertex as a surface.) This parameter should be fixed considering a certain level of noise. As a rule of thumb, Page et al. [24] proposed the following formula to design for a specific crease angle / (see Fig. 4b):

rffiffiffiffiffiffiffiffiffiffiffiffi! 1 : / ¼ 2 arctan þ1

Fig. 3. Example of eigenvalues for different features.

ð11Þ

For example,  ’ 100 is correspondent to having set a crease angle to 10° and this is the default setting in our several experiments.

pffiffiffiffi

ri Ei

e1; E e 2 , and E e 3, and let B be the matrix whose columns are E respectively. Then the matrix B can be written as

0 pffiffiffiffiffiffi

r1

B B ¼ RR ¼ R@ 0 0 1 2

0 pffiffiffiffiffiffi

r2

0

0

1

C 0 A pffiffiffiffiffiffi

r3

and

 1  1 T 1 1 BBT ¼ RR2 RR2 ¼ RR2 R2 RT ¼ RRRT ¼ A: According to Labelle and Shewchuk [13], B is referred to as deformation tensor, which is used for mesh anisotropy. There is another way to look at this coordinate frame. Given a vector v, we can see that

v T Av ¼ v T BBT v ¼ ðBT v ÞT ðBT v Þ ¼ qT q ¼ kqk2 : 3.3. Normal tensor reconfiguration

ð12Þ

The eigenvectors indicate the principal directions, but the absolute sign is not unknown, since

Eq. (12) shows that the quadric form defines a norm in the tensor space. So we define the following metric distance between two mesh vertex Vi and Vj:

N fk0 ¼ N fk ) N fk N fk

liji 

T

¼ N fk0 N fk0 T :

A

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðV j  V i ÞT Ai ðV j  V i Þ ¼ kBTi ðV j  V i Þk: A liji

ð13Þ

A ljij

We solve this problem in the following steps. Calculate the inner product between E1 and the vertex normal N vi . If Ei  N vi < 0, then we change the sign of E1. Next we consider crease saliency as shown in Fig. 4b. In this case, we only use eigenvector E3 as the crease direction, since the other vectors are geometrically less meaningful. So we replace E1 and E2 with the mesh vertex normal N vi and E3  N vi respectively. Examples are shown in Figs. 5 and 6.

Here we note several things: (a) – in general, (b) this distance indicates the one in the normal direction, therefore if Vj exits on the tangent plane at Vi, for instance, Vi and Vj are on a same plane or on the crease line, then A liji ’ 0, since there is no difference between N Vi and N Vj , and (c) if A is identity matrix, then Eq. (13) equals to the Euclid distance in R3 denoted by lij.

4. Our proposed method

The normal curvature is the curvature of the curve projected onto the plane containing the curve’s tangent and the surface normal direction, and is given by the following formula [21,29]:

Using the normal tensor framework described in the previous section, we prepare the metric distance in the tensor space (Section 4.1) and mean curvature approximation (Section 4.2) for designing a weighting function of Laplacian operator, and construct a moderate tensor field

4.2. Mean curvature approximation

jij ¼

2kN vi  ðV j  V i Þk kðV j  V i Þk2

:

ð14Þ

134

S. Tsuchie, M. Higashi / Graphical Models 74 (2012) 130–139

4.3. Normal tensor smoothing

Fig. 4. Classification of feature saliencies, surface, crease and corner. Each saliency has the particular direction given by the eigenvector.

In the case where the local shape of a mesh is almost isotropic such as a plane and sphere, we cannot decide the appropriate directions of E2 and E3 along to the surrounding shape flow. Fig. 6a shows an example in such case. And in another case that data contains a large level of noise, it might be effective to smooth the tensor before the main denoising process. So we perform Laplacian smoothing for our normal tensor A:

1 Ani þ X

Anþ1 i

Fig. 5. Visualization of the eigenvectors at each vertex on a torus. Left figure shows E2 (left side) and E3 (right side) separately. Right figure shows three eigenvectors.

X

  wAij Anj  Ani ;

wAij j2KV ðiÞ j2KV ðiÞ

ð17Þ

where the weight function wAij plays an important role, since simple functions such as angular- or area-weight have lost and blurred the feature characteristics due to a diffusion process. Therefore, our strategy is as follows. If Ai is the tensor at the crease vertex, then we do nothing for it. Otherwise we perform Eq. (17) with the following weight

wAij ¼

1 : ~ ij k kj

ð18Þ

We fix the number of iteration to 3 in our experiments, since a large number of it has influence on the efficiency of the whole denoising process though it should be decided depending on the convergence of the tensor smoothing. An example is shown in Fig. 6b. Before proceeding to the next section, we note that wAij ~ ij k ¼ 0. In order to avoid this probgoes to infinity when kj ~ ij k if kj ~ ij k < d. In our experlem, we set a small value d to kj iment, we use d = DBL_EPSILON. 4.4. Mesh denoising Fig. 6. Geometric flow on the FanDisk model. For visual convenience, only Ejbf3 (=minimum direction of curvature) is displayed.

Langer et al. [15] proposed a formula to approximate the mean curvature KH instead of Eq. (4). Their method is given by a weighted sum of normal curvatures, which can be seen as an extension of a result from mathematical folklore P that states K H ¼ 1n i ji that is a generalization of KH = (jmax + jmin)/2. Using this folklore, we will present another formula below. (Proof is shown in Appendix A.) Proposition. Mean curvature vector K H i is given by the ~ ij , which are defined as average of a normal curvature vector j follows:

In order to introduce our mesh denoising scheme, we use a surface fairing technique known as the minimization  R of bending energy of a plate, min S j2min þ j2max dS, in CAGD [9]. From [20], the average value of the square of normal curvature for all directions at Vi is given by

X  2 1 3 ~ ij k2 ¼ j2min þ j2max þ jmin jmax : kj jKV ðiÞj j2K ðiÞ 8 8 V

For a closed surface, furthermore, the integral of the GaussR ian curvature, S jmin jmax dS, is a topological constant [10], thus we omit it in our energy minimization formulation. Hence, we obtain the following objective function:

min W ¼ min

X

wi ;

ð19Þ

V i 2MðV;FÞ

where

X 1  K Hi ¼ j~ ij ; jKV ðiÞj j2K ðiÞ V pffiffiffi T 2Bi ðV j  V i Þ j~ ij  ; ? lij lij ? lij

ð15Þ

wi ¼

X

~ ij k2 : kj

j2KV ðiÞ

ð16Þ

where is the projected distance of the edge ViVj on the tangent plane at Vi.

From Eqs. (16) and (13), (20) is written as

wi ¼

X 2ðV j  V i ÞT Ai ðV j  V i Þ :  2 ? j2KV ðiÞ lij lij

ð20Þ

135

S. Tsuchie, M. Higashi / Graphical Models 74 (2012) 130–139

Fixing the denominator as a constant parameter, the gradient of W with respect to Vi is given by ! X 2ðV j  V i ÞT Ai ðV j  V i Þ X 2ðV i  V j ÞT Aj ðV i  V j Þ @W @ ¼ þ ? ? @V i @V i j2K ðiÞ ðlij lij Þ2 ðlij lij Þ2 j2KV ðiÞ V X ðAi þ Aj ÞðV j  V i Þ ¼ 4 : ? 2 ðlij lij Þ j2KV ðiÞ Therefore, using the gradient decent method, the vertex position can be updated as follows:

V nþ1 i

V ni

   Ani þ Anj V nj  V ni X 4ki þ X wij ;  2 ? wij j2KV ðiÞ l l j2KV ðiÞ

ð21Þ

SumAi þ ¼



2 1 ; ~ ij k kj  4 ki ¼ min flij g :

ð22Þ

?

!2 ð24Þ

;

and finally substituting Eq. (24) into Eq. (21), we obtain the following vertex update scheme:

 X Ani þ Anj  n 2ki n V ni þ X V  V : j i wij j2KV ðiÞ ðlAiji Þ2

V nþ1 i

ð25Þ

j2KV ðiÞ

Algorithm 1 sums up the main steps performed during our mesh denoising process. Eq. (22) has very important meanings in our algorithm. In the case where Vi is on the crease edges, the distance A lij ’ 0 due to the feature saliency from the eigenvalue ri. Hence, the crease direction is only denoised, since the A other directions are less weighted due to larger lij , and therefore the crease edges preserve their sharpness. And in the other case, a usual anisotropic denoising is performed by a weighting parameter which corresponds to the square of the inverse normal curvature values. Furthermore, Eq. (25) is equivalent to the Laplacian of normal curvature vector. In order to show this, we rewrite the second term in the right hand side of Eq. (25) as follows: X Ani þ Anj  A

j2KV ðiÞ

ðliji Þ2

V nj  V ni



 1 0 n n n X lij l?ij Anj ðV nj  V ni Þ Ai V i  V j @ A ¼   2 ? ? A lij lij lij lij j2KV ðiÞ l i ij   1 0 T n n n BTi n V ni  V nj X lij l?ij Bn Bj V j  V i i @ A; ¼  ? ? Ai Ai lji lji lij lij j2KV ðiÞ lij lij ?

?

where we assume that Bnj ’ Bni and lij ’ lji in the last expression. Using Eq. (16), we obtain the following result:

V nþ1 i

 X 2ki Bni  n V ni þ X j~ ji  j~ nij : Ai ~ ij klij wij j2KV ðiÞ kj j2KV ðiÞ

A

ðlij Þ2

ðV nj  V ni Þ

end for V nþ1 i

SumAi V ni þ 2ki  SumW i

end for until ++ iter > max_iterjjconvergent

5. Results

From Eqs. (16) and (13), (22) becomes

lij lij pffiffiffi A 2liji

ðAni þAnj Þ

ð23Þ

j2KV ðiÞ

wij ¼

Require Mesh M(V, F) and Normal Tensor A iter = 0 repeat Smooth Normal Tensor A with Eq. (17) for each Vi 2 M(V,F) do SumWi = 0 SumAi = 0 for each j 2KV(i) do SumWi + = wij

ij ij

where we define our weight function wij and iteration step size ki as follows:

wij ¼

Algorithm 1. Mesh Denoising

We have implemented the mesh denoising algorithm described in Section 4 using MS VC++2010 on Intel/XeonÒ X5690 (3.46 GHz)  2CPU computer. All meshes are rendered with flat shading to show faceting. To demonstrate a shrinkage property, feature-preserving capability and its error analysis, we compare our results to two methods: the anisotropic MCF (A-MCF) method defined by Eq. (5) and the face normal field (FNF) method by Eq. (6) based on [22]. Fig. 7 shows the results for noisy torus data. The bottom row in the figure indicates the results where the number of iterations is three times of the top row, respectively. We can see that only our algorithm suppresses the data shrinkage or expansion (see Fig. 7g), since it is based on Laplacian of curvature. Other algorithms make the data continue to deform during the mesh denoising. In general, the A-MCF algorithm as itself still remains the shrinkage problem, Hildebrandt and Polthier [8] proposed an additional process called PMC flow in order to tackle it. In the FNF algorithm, it is able to stop  the deformation in a planar region, since e f  V n  V n ’ 0 in Eq. (6); otherwise it has the data conN j i k tinue to deform in the normal direction. Fig. 8 shows the capability of sharp feature preserving. Stopping criteria of the denoising process, i.e. the total number of iterations is due to our own judgment whether planar parts are visually well smoothed or not. We can see that all methods maintain the sharpness of edges. Particularly among them, our scheme, which is based on the tensor smoothing and the anisotropic denoising of the normal curvature vector, can clearly and robustly detect the crease edges which almost disappear due to noises and small angles at creases (see Fig. 8d). This can be seen more clearly in Fig. 9 which is colored by the maps of root mean square (RMS) curvature [26]. With respect to the error analysis, our result has the smallest error and is nearest to the original data (see Fig. 8h). Here, the error between two meshes is defined by the minimum Euclid distance from

136

S. Tsuchie, M. Higashi / Graphical Models 74 (2012) 130–139

Fig. 7. Deformation during the denoising process. jVj = 1343, jFj = 2686. Top row shows the results of each method. Figures from left to right show the AMCF (anisotropic mean curvature flow: Eq. (5)), the FNF (face normal field: Eq. (6)) and our method (Eq. (25)). Bottom row is the result where the number of iteration is three times of the top for each case.

Fig. 8. Comparison of the feature-preserving capability and its accuracy. jVj = 6475, jFj = 12,946. Figures from (b) to (d) show the results of each denoising method. Those from (f) to (h) show the error maps defined by the distance between (a). The distance within ±0.005 mm is colored by blue.

Fig. 9. Comparison of the curvature maps. Root mean square (RMS) curvature [26] is used.

S. Tsuchie, M. Higashi / Graphical Models 74 (2012) 130–139

137

Fig. 10. Bimba model scanned by a laser scanner. 3.7 million triangles.

Fig. 12. Convergence analysis in our method. The vertical axis indicates P the average difference: i2V kV nþ1  V ni k=jVj. i

Fig. 11. Isidore horse model scanned with a laser scanner. Two million triangles.

Table 1 Model size, parameters and run-time in our method. / is defined by Eq. (11). Run-time is measured on the Intel/XeonÒ X5690(3.46 GHz)  2CPU (=12 core) computer and the prototype program that implements our algorithm is parallelized with openMP multi-threading technique. Model

Torus FanDisk Horse Bimba

jFj

2688 12,946 2.0 M 3.7 M

/ (deg)

10 10 10 10

Iteration Ani

V ni

3 3 3 3

20/60 40 40 40

Run-time Fig. 13. Illustration of variables in the proof. 46/140 ms 421 ms 82 s 148 s

each vertex of the original FanDisk model to the mesh triangles after denoising. The bounding box size of the data used in our test is 4.8 mm (width)  5.2 mm (height)  2.7 mm (depth). So regarding the representative scale of model as 5 mm, 86% vertices have been recovered within the maximum error, 0.1% of the model scale. Fig. 10 shows the results on the raw data measured by a laser scanner. In the upper eyelid sulcus as well as the lower one, our method reproduce the original shapes more sharply than A-MCF, since our weighting function is defined by the square of the inverse of curvature and therefore, the larger curvature value tends to keep the shape

more strongly than a usual curvature weight. Fig. 11 shows another one. In this case, our method not only preserves edge features but also retains some intrinsic texture patterns well. Table 1 lists the model size, parameters used in our method and run-time for each data. Fig. 12 shows the convergence of our algorithm, and 40 iterations in Table 1 are to be appropriate for the convergence of our scheme from a numerical point of view.

6. Conclusion We have proposed a novel mesh denoising method based on the normal tensor framework, which is well known as a powerful technique for detecting feature saliencies robustly. In this paper, we utilize the framework

138

S. Tsuchie, M. Higashi / Graphical Models 74 (2012) 130–139

directly for the mesh denoising. Accordingly, we have shown the effectiveness of our algorithm comparing with the existing ones: our method (i) preserves sharp features more robustly and clearly, (ii) does not make data shrink, and (iii) is free from a cumbersome setting of parameters. Future research includes applying our method to treat the mesh decimation, reconstruction, and so on. Acknowledgment The model of Figs. 10 and 11 are courtesy of the AIM@SHAPE Shape Repository. Appendix A. Proof of Eq. (15)

Proof. We assume the simple situation, that is given by a curve, as shown in Fig. 13. What we want to show is that the following relation is satisfied:



1 X 1 X jij N vi ¼ j~ ij : 2 j2fright;leftg 2 j2fright;leftg

ð26Þ

At first, from Eq. (14) the left hand side of Eq. (26) becomes   1 X 1 2 sin h 2 sin h v 2 sin h 0  jij N vi ¼  þ Ni ¼  : 2 j2fright;leftg 2 l l l 1 since l ¼ kV j  V i k; d ¼ kN vi  ðV j  V i Þk ¼ l sin h. On the other hand, in this case, the coordinates of Vi,Vright and Vleft are given by (0,r), (lcosh, r  lsinh) and (lcosh, r  lsinh) respectively. And the unit normal vectors of the edges ViVright and ViVleft are Nright = (sinh,cosh) and Nleft = (sinh,cosh). Using Eq. (7), the normal tensor Ai results in ! 2

Ai ¼ 2

sin h

0

0

cos2 h

) Bi ¼

pffiffiffi sin h 0 ; 2 0 cos h

?

and using the relation lij ¼ kV j  V i k? ¼ l cos h, the right hand side of Eq. (26) is as follows:

pffiffiffi T 1 X 1 X 2Bi ðV j  V i Þ j~ ij ¼ ? 2 j2fright;leftg 2 j2fright;leftg lij lij pffiffiffi    pffiffiffi sin h 0 l cos h 2 ¼ 2 2 0 cos h l sin h 2l cos h  

pffiffiffi sin h 0 l cos h 2 sin h 0 :  ¼ þ 2 l 1 0 cos h l sin h References [1] C.L. Bajaj, G. Xu, Anisotropic diffusion of surface and functions on surfaces, ACM Transactions on Graphics (Proc. of ACM SIGGRAPH 2003) 22 (1) (2003) 4–32. [2] Z. Bian, R. Tong, Feature-preserving mesh denoising based on vertices classification, Computer Aided Geometric Design 28 (2011) 50–64. [3] U. Clarenz, U. Diewald, M. Rumpf, Anisotropic geometric diffusion in surface processing, in: Proceedings of IEEE Visualization, 2000, pp. 397–405. [4] U. Clarenz, U. Diewald, M. Rumpf, Processing textured surfaces via anisotropic geometric diffusion, IEEE Transactions on Image Processing 13 (2) (2004) 248–261.

[5] M. Desbrun, M. Meyer, P. Schroder, A.H. Barr, Implicit fairing of irregular meshes using diffusion and curvature flow, ACM Transactions on Graphics (Proc. of ACM SIGGRAPH 1999) (1999) 317–324. [6] J.R. Diebel, S. Thrun, A bayesian method for probable surface reconstruction and decimation, ACM Transactions on Graphics 25 (1) (2006) 39–59. [7] S. Fleishman, I. Drori, D. Cohen-Or, Bilateral mesh denoising, ACM Transactions on Graphics (Proc. of ACM SIGGRAPH 2003) 22 (3) (2003) 950–953. [8] K. Hildebrandt, K. Polthier, Anisotropic filtering of non-linear surface features, Computer Graphics Forum 23 (3) (2004) 391–400. [9] J. Hoschek, D. Lasser, Fundamentals of Computer Aided Geometric Design, Wellesley, MA, 1993. [10] P. Joshi, C. Séquin, Energy minimizers for curvature-based surface functionals, Computer-Aided Design & Applications 4 (5) (2007) 607–617. [11] E. Kalogerakis, P. Simari, D. Nowrouzezahrai, K. Singh, Robust statistical estimation of curvature on discretized surfaces, in: SGP ’07 Proceedings of the 5th Eurographics Symposium on Geometry Processing, 2007, pp. 13–22. [12] L. Kobbelt, S. Campagna, J. Vorsatz, H.-P. Seidel, Interactive multiresolution modeling on arbitrary meshes, ACM Transactions on Graphics (Proc. of ACM SIGGRAPH 1998) (1998) 105–114. [13] F. Labelle, J.R. Shewchuk, Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation, in: SCG’03: Proceeding of the 19th Annual Symposium on Computational Geometry, 2003, pp. 191–200. [14] Y.-K. Lai, Q.-Y. Zhou, S.-M. Hu, J. Wallner, H. Pottmann, Robust feature classification and editing, IEEE Transactions on Visualization and Computer Graphics 13 (1) (2007) 34–45. [15] T. Langer, A.G. Belyaev, H.-P. Seidel, Analysis and design of discrete normals and curvatures, Technical Report, Max-Plank-Institute für Infomatik, 2005. [16] K.-W. Lee, W.-P. Wang, Feature-preserving mesh denoising via bilateral normal filtering, in: 9th International Conference on Computer Aided Design and Computer Graphics, 2005, pp. 275–280. [17] Z. Mao, L. Ma, M.Z.X. Xiao, Susan structure preserving filtering for mesh denoising, The Visual Computer: International Journal of Computer Graphics 22 (4) (2006) 276–284. [18] N. Max, Weighs for computing vertex normals from facet normals, Journal of Graphics Tools 4 (2) (1999) 1–6. [19] G. Medioni, M.-S. Lee, C.-K. Tang, A Computatinal Framework for Segmentaion and Grouping, Elsevier, Amsterdam, 2000. [20] E. Mehlum, C. Tarrou, Invariant smoothness measures for surfaces, Advances in Computational Mathmatics (1) (1998) 49–63. [21] M. Meyer, M. Desbrun, P. Schröder, A.H. Barr, Discrete differentialgeometry operators for triangulated 2-manifolds, in: Visualization and Mathematics, 2003. [22] Y. Ohtake, A. Belyaev, I. Bogaevski, Mesh regularization and adaptive smoothing, Computer-Aided Design 33 (11) (2001) 789–800. [23] A.F.E. Ouafdi, D. Ziou, H. Krim, A smart stochastic approach for manifolds smoothing, in: Eurographics Symposium on Geometry Processing 2008, vol. 27(5), 2008, pp. 1357–1364. [24] D.L. Page, A.F. Koschan, Y. Sun, J.K. Paik, M.A. Abidi, Robust crease detection and curvature estimation of piecewise smooth surfaces from triangle mesh approximations using normal voting, in: Proceedings of the International Conference on Computer Vision and Pattern Recognition, vol. 1, 2001, pp. 162–167. [25] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Transactions on Pattern and Machine Intelligence 12 (7) (1990) 629–639. [26] S. Pulla, A. Razdan, G. Farin, Improved curvature estimation for watershed segmentation of 3-dimensional meshes, IEEE Transactions on Visualization and Computer Graphics (2001). [27] S. Rosenberg, The Laplacian on a Riemannian Manifold, Cambridge University Press, 1997. [28] X. Sun, P.L. Rosin, R.R. Martin, F.C. Langbein, Fast and effective feature-preserving mesh denoising, IEEE Transactions on Visualization and Computer Graphics 13 (2007) 925–938. [29] G. Taubin, Estimating the tensor of curvature of a surface from a polyhedral approximation, in: ICCV ’95 Proceedings of the Fifth International Conference on Computer Vision, 1995, pp. 902–907. [30] G. Taubin, A signal processing approach for fair surface design, ACM Transactions on Graphics (Proc. of ACM SIGGRAPH 1995) (1995) 351–358. [31] G. Taubin, Linear anisotropic mesh filtering, IBM Research Report RC22213(W0110-051), 2001.

S. Tsuchie, M. Higashi / Graphical Models 74 (2012) 130–139 [32] J. Weickert, Theoretical foundations of anisotropic diffusion in image processing, Computing, Supplement 11 (1996) 211–236. [33] T.J. Willmore, Riemannian Geometry, Oxford University Press, 1993. [34] S. Yoshizawa, A.G. Belyaev, H.-P. Seidel, Smoothing by example: mesh denoising by averaging with similarity-based weight, IEEE International Conference on Shape Modeling and Applications 23(3) (2006) 38–44.

139

[35] Y. Zhang, A.B. Hamza, Vertex-based diffusion for 3-d mesh denoising, IEEE Transactions on Image Processing 16 (4) (2007) 1036–1045. [36] Y. Zheng, H. Fu, O. Au, C.-L. Tai, Bilateral normal filtering for mesh denoising, IEEE Transactions on Visualization and Computer Graphics 17 (10) (2011) 1521–1530.