A Discrete Spring Model for Generating Fair Curves and Surfaces
Atsushi Yamada1, Kenji Shimada2, Tomotake Furuhata1, and Ko-Hsiu Hou2
1
Tokyo Research Laboratory, IBM Japan Ltd., LAB-S73 623-14, Shimotsuruma, Yamato, Kanagawa, 242-8502, JAPAN [e-mail], [e-mail]
2
Mechanical Engineering, Carnegie Mellon University Pittsburgh, PA 15213 [e-mail], [e-mail]
The ability to generate fair curves and surfaces is important in computer graphics (CG), computer-aided design (CAD), and other geometric modeling applications. In this paper, we present an iteration-based algorithm for generating fair polygonal curves and surfaces that is based on a new discrete spring model. In the spring model, a linear spring, whose length approximately represents a curvature radius, is attached along the normal line of each polygon node. Energy is assigned to the difference of the lengths, that is, the difference in curvature radius, of neighboring springs. Our algorithm then minimizes the total energy by an iterative approach. It accepts as inputs (1) an initial polygonal curve (surface), which consists of a set of polygonal segments (faces) and a set of nodes as polygon-vertices, and (2) constraints for controlling the shape. The outputs are polygonal curves (surfaces) with smooth shapes. We also describe a method for improving the performance of our iterative process to obtain a linear execution time. Our algorithm provides a tool for the fair curve and surface design in an interactive environment. Keywords: geometric modeling, fair surface design, polygonal models, energy minimization
1. Introduction
The purpose of this paper is to provide an algorithm for generating fair curves and surfaces for use in the fields of computer graphics (CG) and computer-aided design (CAD). Generation of fair shapes is a major topic in shape design [5][11][19][21][25][26][27][29]. It is also required in applications such as fitting of smooth shapes to scattered points [8][23], texture mapping [16], and so on. The curves and surfaces treated in this paper are represented in polygonal form. The inputs of the algorithm are (1) an initial polygonal curve (surface) consisting of a set of nodes and a set of polygonal segments (faces), and (2) constraints for controlling the shape. The algorithm moves the nodes to suitable positions while
minimizing the curvature variation under the given constraints. The outputs are smoothly shaped polygonal curves (surfaces). The polygonal surfaces treated in this paper are not limited to triangular meshes; they also include quadrilateral meshes. Theoretically, nsided faces such as pentagons or hexagons may also be included in the meshes. Our algorithm follows an iterative approach of the Gauss-Seidel type: node positions are iteratively updated under the given constraints. To update the node positions, two types of spring force are applied to each node: (1) a force acting in the normal direction, to optimize the curvature variation and (2) a force in the direction perpendicular to the normal, to optimize the node distribution. Then main idea of this paper lies in the discrete spring model that produces the former force minimizing curvature variation. Since curvature is a natural measurement of fairness, our spring model produces fair curves and surfaces. As Taubin [26] points out, most energy-minimization approaches are expensive in terms of time and space. This paper also provides a method for achieving a linear execution time. The remainder of the paper is organized as follows. In Section 2, we summarize previous work. After defining a discrete spring model in Section 3, we present a curve (surface) modeling algorithm based on the spring model in Section 4. Section 5 describes some results obtained by using the algorithm, and Section 6 summarizes the paper.
2. Previous Work
A considerable amount of work has been done on fair surface modeling. We consider that most of the previous work related to polygonal surface modeling can be classified into two types: finitedifference approaches and finite-element approaches. Our approach belongs to the former category. In both types of approach, node positions are either updated iteratively or are calculated by solving a large sparse linear system. In finitedifference approaches, values are evaluated only at nodes, while in
finite-element approaches, values are evaluated over faces by using numerical integration. One key factor that distinguishes several approaches within each category is the spring model that is used to move the node positions. Irrespective of the classification, when the displacement of a deformation is large; that is, when the shape of the initial polygonal surface is very different from the final shape, the deformation reduces to a non-linear problem. Therefore, if solving a sparse linear system does not yield a suitable answer, the system must be solved iteratively. Large deformations appear often in the process of surface design. To construct a robust algorithm for such cases, we select an iterative approach. The following two subsections survey both types of approach to surface modeling. In the third subsection, some other applications that use fairing as a part of their algorithm are also surveyed.
2.1 Finite-Difference Approaches to Surface Modeling
Laplacian smoothing is the easiest way to generate a fair surface. This approach iteratively moves the position of a node to the barycenter of its neighboring nodes. It has the characteristic that the area of the surface is minimized under given constraints, and is often used to improve the geometrical irregularity of a mesh in the field of finite-element meshing [12]. But when it is applied to surface modeling, the problem arises that a sharp tip is generated in the neighborhood of a node fixed by a constraint. Moreover, it is not possible to control normals by giving normal constraints. For surface modeling, Szeliski and Tonnesen [25] used a particle system, in which each particle is defined by its position and normal, and the population of particles is controlled automatically. To obtain a fairly triangulated surface, they minimized the weighted sum of three energy factors: (1) coplanarity, (2) co-normality, and (3) co-circularity. The co-planarity condition causes neighboring nodes to lie on each other' tangent s plane, the co-normality condition modifies irregular twisting of two neighboring nodes, and the co-circularity condition preserves constant curvature on edges connecting neighboring nodes. Whereas their approach uses three factors, our spring model achieves a similar effect by using just one factor. Mallet [18][19] provided a general formulation for discrete surface interpolation with many adjustable parameters. When harmonic weighting [19] is selected among his proposed parameters - as in his results -, his approach minimizes the sum of the distances, each of which is from a node to the barycenter of its neighboring nodes. Taubin [26] proposed a fair surface design approach in which
high-frequency terms such as noise are removed by using a technique based on Fourier analysis. Due to its linear execution time, his approach is powerful for polygonal surfaces with millions of nodes, such as one captured by using a range scanner. The approach involves iterative repetition of a Gaussian operator, which moves a node to a suitable position on a line connecting a node and the barycenter of its neighboring nodes. This position is determined in such a way as to avoid shrinkage of the entire shape. Basically, Taubin's approach is suitable for removing noise from a given initial shape without causing shrinkage. The main difference between our approach and those of Mallet and Taubin is that in our approach curvature variation is used as a measurement of minimization. Curvature is a natural measurement of fairness; consequently, our approach produces a fairer shape than theirs. Welch and Witkin [29] also proposed a mesh-based modeling method for surfaces with arbitrary topology. After defining a local surface equation in the neighborhood of each node, which fits the latter's neighboring nodes in the least-square sense, they minimized a fairness norm based on curvature, which is calculated from the local surface. The local surface is temporarily used to evaluate the curvature at a node, and the final outputs are only node positions, not surfaces of faces. In the sense that curvature is adopted as a fairing measurement, our approach is similar to theirs. However, their approach has to solve a 5 m least-square linear system (where m is the number of neighboring nodes) for each node at each iteration. This is a time-consuming process if the number of nodes is large. Kuriyama and Tachibana [16] applied Gaussian operators not only to node positions but also to second-order derivative vectors at nodes. The magnitude of a second-order derivative vector is equal to the curvature; therefore, their approach can be interpreted as using Gaussian operators to optimize curvature. Kobbelt and his colleagues [14] also proposed a similar approach. Eck and Jaspert [7] proposed a fairing algorithm of planar and spatial curves based on the measurement of discrete curvature and torsion. Their approach timidly moves only one node in each iteration, where the new position of the node is determined on a line connecting the node and the weighted barycenter of its neighboring nodes. The suitable position on the line is also calculated in a heuristic way.
2.2 Finite-Element Approaches to Surface Modeling
In comparison with finite-difference approaches, finite-element approaches generally yield a higher-quality surface and generate an explicit surface equation for each face, while they are
computationally more expensive. Celniker and Gossard's approach [5] generates a C1continuous surface by using a finite-element technique. They applied Zienkiewicz' shape function [30], which was originally s proposed for finite element analysis, to attain C1 continuity in surface modeling. Their approach minimizes a weighted combination of the energy factors of a membrane and a thin plate. Welch and Witkin [28] also proposed an approach very similar to that of Celniker and Gossard, although the former has more general formulations for both energy factors and shape control constraints. Moreton and Sequin [21][22] presented a different type of finite-element approach in which they use biquintic Bezier patches and a fairness norm based on measures of curvature variation. Their approach produces the most impressive surface, but it is too expensive for use in an interactive environment.
the spring lengths |Pi - H|| and |Pj - H|| of a V-shape formed by Pi, H, and Pj. The spring length approximately represents the curvature radius; therefore, keeping the spring length equal is equivalent to minimizing variation in the curvature radius. Suppose node Pj is fixed by a constraint. If |Pi - H|| is smaller than |Pj - H||, as shown in Figure 2(A), node Pi moves to a new position along the normal Ni in the direction that enlarges |Pi - H|| to the size of |Pj - H||. In contrast, if |Pi - H|| is larger, node Pi moves to a new position along the normal Ni in the direction that shortens the |Pi - H|| to the size of |Pj - H||. In a stable configuration, Pi and Pj are considered to be on a circular arc whose center is at H and whose curvature radius is |Pj - H|| (=||Pi + dPi - H|). | Because of this V-shaped configuration of virtual springs, our spring model is named "V-Spring."
Nj Ni Pj Pi
2.3 Other Applications Using Fairing
Halstead and his colleagues [11] provided a solution to the problem that Catmull-Clark's subdivision surface scheme [4] gives a shrunken surface without interpolating given control points. In their process, a fairness factor proposed by Celniker and Gossard is used in conjunction with the Catmull-Clark scheme. The aim of Eck and Hoppe' work [8] is to generate tensor product surfaces s from scattered points. They use a thin plate as a fairness factor to remove undulations in surfaces. Levy and Mallet [16] applied their discrete smooth interpolation approach [18][19] to the nondistorted texture-mapping problem. Koch and his colleagues [15] applied the thin-plate approach to a system for simulating facial surgery, and DeCarlo and his colleagues [6] applied the thin-plate approach to geometric modeling of human faces. As can be seen in the above literature in this subsection, fairness factors are widely used in the area of CG and CAD alongside factors unique to each study, which depend on their applications. Our approach has the potential to be applied for such applications.
H
!
"
#
# #
. Here we formally define a displacement dPij, from a current position to a new position, of a node Pi by the force of the spring model on the supposition that Pj is fixed. Let Pi and Pj be two nodes, and let Ni and Nj, be unit normal vectors associated with the nodes, respectively. We assume the inner product dot(Ni, Nj) to be positive; therefore, if two normal vectors with their negative inner product are given, let either Ni or Nj be the direction-reversed vector. An intersection of the two normal lines is denoted by H (see Figure 2(A)). Let ti and tj be real values satisfying
Pi - H = ti Ni, Pj - H = tj Nj.
(1)
3. Definition of the Discrete Spring Model: V-Spring
Figure 1 shows a planar curve and its normal lines at some sampling points on the curve. Consider the normal lines at two neighboring sampling points Pi and Pj. For a curvature-continuous curve, if Pi approaches Pj, the intersection point H of the normal lines converges to a center of curvature at Pj [3]. Therefore, our idea is to attach a linear spring, as shown in Figure 2(A), to each normal line of a node consisting in a polygonal curve (surface). The linear spring works to keep equal
|ti| and |tj| correspond to the distances |Pi - H|| and |Pj - H||, respectively, because Ni and Nj are unit vectors. We define the displacement dPij of node Pi by our spring model as
dPij = (tj - ti ) Ni.
(2)
(tj - ti) in the equation is calculated as follows. For the equation in which H is deleted from Equation (1), by calculating the inner products with Ni and Nj, we obtain the following equations:
dot (( Pj - Pi + tiNi - tjNj ), ( Ni )) = 0,
dot (( Pj - Pi + tiNi - tjNj ), ( Nj )) = 0.
(3)
Nj Pj Ni dPij Pi Pi
By solving Equation (3) in terms of ti and tj, we obtain
(1 - cos (a)), t = dot ((P - P ), (- cos( a ) N + N )) (1 - cos (a)),
ti = dot ((Pj - Pi ), (- Ni + cos(a ) Nj ))
j j i i j 2 2
where cos(a) = dot( Ni, Nj ). Consequently, (tj - ti) in Equation (2) is determined as follows:
tj - ti = dot ((Pj - Pi ), ( Ni + Nj )) dot ((Pj - Pi ), ( Ni + Nj )) . = 1 + dot ( Ni, Nj )
Nj Pj
H Ni Pi
(1 + cos( a) )
(4)
Pi
dPij
Therefore, Equation (2) is written as
dot (( Pj - Pi ), ( Ni + Nj )) dPij = 1 + dot ( Ni, Nj ) Ni.
Hj Hi
(5)
The denominator of Equation (4) always has a non-zero value, because cos(a) (=dot(Ni, Nj,)) is assumed to be positive. One may consider from Figure 2(A) that a numerical error occurs when two normal lines are parallel because of the absence of H. However, Equation (5) does not use H directly; therefore, our spring model is stable even in the case of parallel normal lines. In this case, Equation (5) is deduced to the following equation: dPij = dot ((Pj - Pi ), N )N , ( N = Ni = Nj ). In the case of a planar curve, normal lines along Ni and Nj always have an intersection. However, if we consider a non-planar curve or a surface, normal lines do not always intersect at a point. Therefore, we cannot use Equation (1) as it is for non-planar cases. Instead of an intersection point H in Equation (1), we use Hi and Hj, which are the feet of the shortest line segment connecting two normal lines (see Figure 2(B)). Then, for non-planar cases, we modify Equation (1) to obtain the following equations:
Pi - H i = t i N i , Pj - Hj = tj Nj.
Figure 2: Definition of the spring model. (A) Displacement dPij for a planar curve. (B) Displacement dPij for a non-planar curve or a surface.
4. Curve/Surface Modeling Using a Discrete Spring Model
4.1 Overview of the Algorithm
A polygonal surface is defined as a pair of a set of nodes Pi (i = 1, ..., n), and a set of polygonal faces. We define neighboring nodes Pj (j = 1, ..., m) of Pi to be a set of nodes connected to Pi by polygonal edges. Our algorithm can be classified as an iterative method of GaussSeidel type. The positions and normals of nodes are updated in each iteration, and the iterations are continued until the termination condition is satisfied. In one iteration, the updated latest positions are always used to calculate the positions of other nodes. Figure 3 shows an overview of our algorithm.
Let Pi be the position of the i-th node Let n be the number of nodes While( termination condition is not satisfied ){ For( all nodes Pi (i = 1, ..., n) ){ Step 1: Calculate pseudo-normal Ni } For( all nodes Pi (i = 1, ..., n) ){ Step 2: Calculate displacement dPi caused by the force exerted by V-Spring Step 3: Calculate displacement dPu, i caused by the force for regularizing node distribution Step 4: Pi = Pi + ( dPi + dPu, i ) } }
(6)
For non-planar cases, the equations corresponding to Equation (3) are
dot (( Hj - Hi ), ( Nj )) = 0. dot (( Hj - Hi ), ( Ni )) = 0,
(7)
By solving Equation (7), we can derive an equation for (tj - ti) in Equation (2) that is the same as Equation (4).
Figure 3: Overview of the algorithm.
$
The following subsections 4.2, 4.3, and 4.4 describe Steps 1, 2, and 3 in Figure 3, respectively. In Section 4.5 we describe several constraints. After describing the termination condition in Section 4.6, we give a method for reducing the execution time in Section 4.7. Up to Section 4.7, our discussion concerns surface modeling. Section 4.8 extends the discussion to curve modeling.
Ni dPi Ni dPi P2 Pi-1 Pi +1 Pi P1 Pm
Pi
4.2 Pseudo-normal Calculation
The first step of the algorithm is to calculate a unit normal vector Ni for each node Pi of a polygonal surface (Step 1 in Figure 3). The polygonal surface is a discrete model; therefore, the unit normal vector must be calculated only approximately. In our implementation, we calculate the unit normal Ni by averaging the normals of polygonal faces adjoining the node. There are more sophisticated ways to calculate the normal. For example, one way is to calculate the normal from a sphere fitted to the target node and its neighboring nodes in the least-square sense. However, this is more time-consuming, and fails when the target node and its neighboring nodes are on the same plane or the neighboring nodes are placed so that the target node is a saddle point. Another way, used by Taubin [26], is to calculate the normal as the average of the vectors from target node to each of its neighboring nodes. This is faster, but it fails when the target node and its neighboring nodes are on the same plane. Our choice is the more basic but the more robust way. In the early phase of iterations, the normal vector is unreliable; however, as the iterations proceed, the normal vector converges to the reliable normal of the fair surface.
Figure 4: Node displacement by forces exerted by neighboring nodes. (A) Planar curve case. (B) Surface case.
4.4 Node Displacement for Regularizing Node Distribution
In finite-difference approaches, it is important to maintain the regular node distribution during the iteration process, because uneven distribution of nodes results in incorrect estimation of curvature. To obtain the regular node distribution, we use a variation of a Laplacian smoothing operator [12], which is a popular and effective way for removing the irregularity. The Laplacian operator moves a node to the barycenter of its neighboring nodes. However, applying the regular Laplacian operator offsets the displacement dPi in Equation (8). Therefore, our idea is to use only a component dPu,i, that is perpendicular to the normal Ni, of the displacement of the Laplacian operator. only Using this component creates two displacements dPi and dPu,i perpendicular to each other. Therefore, the two displacements do not offset each other. The displacement dPu,i is written as follows:
dPu , i = dPu 0, i - dPu1, i, m dPu 0, i = Pj m - Pi, dPu1, i = dot (dPu 0, i, Ni ) Ni. j =1
4.3 Node Displacement by the Force of V-Spring
Node Pi obtains forces from its neighboring nodes Pj (j = 1, ..., m). Each force works to keep the edge from Pi to Pj in a circular arc. Each of the displacements dPij (j = 1, ..., m) is calculated from Equation (5). Node Pi moves to a new position along the normal Ni by the weighted average dPi of the displacements dPij (j = 1, ..., m) as follows:
dPi =
4.5 Constraints
Constraints are considered to be external forces for controlling the shape of a surface. Various kinds of constraints can be considered, depending on the requirements of applications. In this section, we describe several important constraints used in surface modeling. We classify the constraints into two types: direct and indirect constraints.
w dP w ,
j ij j j =1 j =1
m
m
(8)
where wj (j = 1, ..., m) are weights for the averaging. In our implementation, the weight wj is determined by the inverse of the length of the edge connecting Pi and Pj.
wj = 1 Pi - Pj , ( j = 1, 2, ..., m) .
Direct Constraints
Direct constraints are given directly for a certain node. Major %
constraints are: - Positional constraints, - Normal constraints. A positional constraint fixes a node to a certain position during the iterations; therefore, Steps 2 and 3 in Figure 3 are skipped for the fixed node. A normal constraint fixes a normal of a node to a certain direction. For the normal fixed node, the pseudo-normal calculation (Step 1 in Figure 3) is skipped and the given normal is assigned. In our approach, positional constraints must be given at least to end nodes for the case of an open polygonal curve and to boundary nodes for the case of an open polygonal surface. That is why the use of a Laplacian operator causes shrinkage of the shape. If it is necessary to modify the shape of the boundary curves of a surface, start from the modeling of boundary curves and go on to the modeling of the surface bounded by the boundary curves.
dPc , i =
j =1
m
wj dPc , ij
w ,
j j =1
m
where wj (j = 1, ..., m) are weights for averaging. In our implementation, the weight wj is determined by the inverse of the distance from node Pi to its foot Qj to the tangent plane at Pi (see Figure 6(B)).
wj = 1 Pi - Qj , ( j = 1, 2, ..., m) .
Indirect Constraints
Indirect constraints are not directly connected to certain nodes. During the iterations, connections between the constraints and the nodes are updated dynamically. One major constraint we introduce here is scattered points. This is an important application in CAD and CG for generating a smooth surface fitted to scattered points in the least-square sense. In the application, the constraint by a scattered point affects the point on a surface closest to it. Therefore, in our discrete model, we update the connection between the scattered point and its closest node during the iterations. We have to modify the algorithm slightly to introduce the indirect constraints; the modified algorithm is shown in Figure 5. In Figure 5, all the steps except 2 and 5 are the same as in Figure 3. In Step 2, each scattered point is connected to its nearest node. In practice, it is not necessary to perform Step 2 at every iteration; once in every several iterations is enough. To find the nearest node, it is enough to search neighbor of the previous nearest node in first and second orders, except when performing a first search. In Step 5, external forces from connected scattered points Vj (j = 1, ..., m) are applied to a node Pi. The displacement of Pi by Vj is calculated as follows:
dPc , ij = kj dot (Vj - Pi, Ni ) Ni,
The weighting means that the point closer to Pi contributes more. When Qj is identical with or very close to Pi, a sufficiently large value is assigned to its wj to avoid division by zero. In Equation (9), weight kj is used to make a trade-off between fairing and keeping proximity to the scattered point Vj. A larger kj approximates Vj more closely. In the case of a planar curve, least-square fitting is geometrically interpreted in such a way that, when a spring is attached to each line from a scattered point to the nearest point on a curve, the sum of the internal energies of the springs are minimized (see Figure 6(A)). Our approach is geometrically interpreted as being to attach a spring to a line from a scattered point to its foot in the tangent plane at the nearest node (see Figure 6(B)). We therefore consider that our approach is approximately equivalent to least-square fitting. The advantage of our approach is that it can be applied not only to surfaces with regular topology, such as tensor product surfaces, but also to surfaces with arbitrary topology. In addition, theoretically, n-sided polygons such as pentagons or hexagons may be included in the polygonal surface.
Let Pi be the position of the i-th node Let n be the number of nodes Let Vi be the i-th scattered point Let l be the number of scattered points While( termination condition is not satisfied ){ For( all nodes Pi (i = 1, ..., n) ){ Step 1: Calculate pseudo-normal Ni } For( all scattered points Vi (i = 1, ..., l) ){ Step 2: Make connection of each scattered point to its nearest node } For( all nodes Pi (i = 1, ..., n) ){ Step 3: Calculate displacement dPi caused by the force exerted by V-Spring Step 4: Calculate displacement dPu, i caused by the force for regularizing node distribution Step 5: Calculate displacement dPc, i caused by the force exerted by scattered points Step 6: Pi = Pi + (dPi + dPu, i + dPc, i) } }
(9)
where kj denotes a weight assigned to Vj. The displacement dPc,ij is interpreted as a kj -times component of vector Vj - Pi along the normal Ni. In the same manner, we calculate dPc,i1, ..., dPc,im for all connected scattered points. The final displacement of node Pi by the forces of its connected scattered points is determined by weighted-averaging of dPc, ij (j = 1, ..., m), as follows:
Figure 5: Overview of the fitting algorithm. &
Vm Qm
Ni P i Q2 V2 V1 Q1
As the post-smoothing proceeds, the rough shape approaches the precise shape. In the post-smoothing, iterations are performed at each level until the termination condition described in Section 4.6 is satisfied. Our termination condition is normalized by the resolution; this provides an efficient way of determining the time at which to move on.
Figure 6: Least-square fitting of a planar curve to scattered points. (A) Geometric interpretation of general least-square fitting. (B) Geometric interpretation of our fitting approach.
M1 Pre-smoothing
M1 Post-smoothing
4.6 Termination Condition
To determine when to terminate iterations, the maximum among the norms of all node displacements is compared with a given threshold e. If the maximum norm is less than the threshold e, the iterations are terminated. The size of the displacement depends on the resolution of the polygonal surface; therefore the displacement should be normalized by the sizes of polygonal faces. In our implementation, we normalize the displacement of each node by the average length of its neighboring edges. The maximum norm of the normalized displacements is then compared with the threshold e.
M2
M2
M3
Figure 7:
'
#
#
4.8 Extension to Curve Modeling
Basically, the algorithm for curve fairing is the same as the one for surface fairing. In each iteration, the node Pi is moved to a new position by the forces exerted by two neighboring nodes Pi-1 and Pi+1 (see Figure 4(A)). The case of a planar curve does not involve any extension of the surface case; however, in the case of a non-planar curve, the calculation of the pseudo-normal Ni is more difficult than in the surface case. From a sequence of nodes Pi-1, Pi, and Pi+1, we calculate the unit tangent Ti and the unit binormal Bi as follows:
Ti = ( Pi +1 - Pi -1) Pi +1 - Pi -1 , Bi = ( Pi - Pi -1) ( Pi +1 - Pi )
4.7 Performance Improvement
Iterative approaches of the Jacobi or Gauss-Seidel types tend to remove high frequencies quickly, while they take many iterations to remove low frequencies. Consequently, if the number of nodes is extremely large, it takes many iterations to achieve convergence. A promising way to reduce the execution time is to employ multigrid methods [2][10][20]. Kobbelt and his colleagues [14] are positively using the multi-grid method to model dense meshes. By using the multi-grid method, a linear execution time can be achieved. The multi-grid method requires polygonal surfaces with several different levels of resolutions. Such surfaces can be obtained by the use of mesh simplification algorithms [1][9][13][24][27]. For example, let M1, M2, and M3 be three levels of polygonal surfaces, where M1 is the finest and M3 is the coarsest (see Figure 7). The V-cycle multi-grid method applies iterations for polygonal surfaces with different levels, in the sequence {M1, M2, M3, M2, M1}. The first half of the process, going down from M1 to M2, is called presmoothing, and the second half of the process, coming up from M3 to M1, is called post-smoothing. In the pre-smoothing, some iterations are performed at each level in order to remove noise. On the coarsest level M3, a rough shape is predicted by the solution.
(Pi - Pi -1) (Pi +1 - Pi )
,
where denotes the outer product. As the outer product of Bi and Ti, we calculate the unit principal normal Ni as follows:
Ni = Bi Ti Bi Ti .
If Pi-1, Pi, and Pi+1 are collinear, Bi and Ni are zero vectors; then the displacement dPi is a zero vector according to Equation (5). The best way to obtain a fair curve is to apply V-Spring forces in the directions of both Bi and Ni; however, in practice, applying a force only in the direction of Ni gives a fair curve, even if the problem is a non-planar case.
5. Results
Figure 8 shows a result of fair surface generation obtained by using our algorithm. Figure 8(A) shows the initial mesh with sharp corners, that is, sudden changes in curvature, while figures (B) and
(C) show the faired results with two different sets of direct constraints. As can be seen in the figures, the algorithm produced fairer surfaces. The surface of Figure 8(B) resulted when the positions of nodes on the inner and outer boundaries were kept unchanged by using the positional constraints. The surface of Figure 8(C) was generated by constraining the normals of the nodes on the inner and outer boundaries, in addition to the positions of these nodes. Figure 10 shows another result of fairing in a shaded image. Figure 9 shows the stability of our algorithm by giving it an extremely noisy mesh as initial data. We added random noise of large amplitude to the mesh in Figure 8(A), and then processed the noisy mesh by using our algorithm. Five iterations of presmoothing generated the smoother surface shown in Figure 9(B). With further processing, the mesh converged to the fair surface shown in Figure 9(C). Table 1 shows the execution time of our algorithm measured for meshes of various resolutions. It indicates that the time grows approximately linearly to the complexity of the given meshes. Table 1: Execution time for fair surface generation. The execution time is measured for the surface shape with the constraints shown in Figure 8(C). Column (a) in the table contains data for the resolution shown in Figure 8(C). Columns (b), (c), and (d) contain data for different resolutions with the same surface shape. These data are measured under the following conditions:
CPU: Pentium II 450 MHz, System: Windows NT 4.0, Threshold e for termination condition: 0.001, Number of levels for multi-grid: 6, Number of iterations at each level of pre-smoothing: 5. (a) (b) (c) Number of nodes 553 1610 3644 Number of faces (triangles) 985 3018 6983 Execution time (sec) 2.38 13.8 36.4
and fitted to the scattered points.
...