Multiresolution Surface Modeling Based on Hierarchical Triangulation

COMPUTER VISION AND IMAGE UNDERSTANDING

Vol. 63, No. 1, January, pp. 1–14, 1996 ARTICLE NO. 0001

Multiresolution Surface Modeling Based on Hierarchical Triangulation1 MARC SOUCY

AND

DENIS LAURENDEAU2

Computer Vision and Digital Systems Laboratory, Laval University, Ste-Fey, Quebec G1K 7P4, Canada Received January 6, 1993; accepted November 22, 1994

There are a set of properties which are highly desirable for a multiresolution surface modeling algorithm. For example, the resulting model should be defined in an objectcentered reference frame and non-uniform sampling of the object surface should be allowed. For triangle-based representations, the surface triangulations should be equiangular in 3-D space in order to increase the upper bound on the curvature of sampled surfaces. Finally, the algorithm should preserve the topology of the object as well as the surface orientation discontinuities. Few of the previously reported multiresolution techniques possess all these properties. Quadtree approaches are domain-dependent representations and require regular sampling [14]. Hierarchical triangulation approaches are more interesting since they are defined in an object-centered reference frame. Among existing methods, ternary triangulations do not meet the equiangularity requirement [3] or are limited to the description of objects without holes [12]. Top-down approaches based on hierarchical Delaunay triangulations lead to equiangular triangulations in 2-D parametric space but do not naturally preserve orientation discontinuties [5, 15]. In the two previous references, it is assumed that the surfaces can be projected on a single parametric grid. It is appropriate to make this assumption for digital terrain modeling. For object modeling, however, this condition limits the applicability of the method to simple objects. The technique proposed in [11] tends to keep orientation discontinuities in the model but does not meet the equiangularity requirement. More recent algorithms reported in [16, 6, 8, 7] are more suited for the representation of 3-D surfaces at different levels of resolution. In the iterative algorithm of Schroeder et al. [16], a vertex Vi is removed from the surface triangulation at a given iteration if the distance between this vertex and the plane approximating its neighboring vertices is smaller than a given threshold. If Vi is effectively removed, the region formed by the triangles having Vi as a vertex is locally retriangulated. The error criterion only considers deviation of the new triangulation from the triangulation created at the previous iteration. Therefore, the global error of a coarse triangulation with

This paper describes a sequential multiresolution surface modeling technique. The initial model consists in a surface triangulation resulting from the integration of a set of range views. Thereafter, a hierarchical triangulation algorithm iteratively removes vertices of the triangulation, always minimizing the retriangulation error. The equiangularity of the surface triangulation is optimized in 3-D space throughout the optimization process. Besides compressing surface information, the proposed technique preserves the topology of the triangulation and surface orientation discontinuities. Experimental models built from complex multipart objects with holes are presented.  1996 Academic Press, Inc.

1. INTRODUCTION

While numerous techniques have been developed to compress 2-D intensity and color images, multiresolution surface triangulation of 3-D objects has not been studied extensively. The compression of high-resolution surface triangulations may find applications in several areas of computer vision. Compression of 3-D surfaces is mandatory in applications that require the real-time display of 3-D objects. In object recognition, the compression can reduce the complexity of the recognition process by representing a complex surface using a small set of significant geometric primitives. Most of the object representation techniques reported in the literature use closed form geometric primitives to describe object parts, such as superquadrics or surface polynomials. Closed form geometric primitives are very compact since they are described by few parameters. They lack, however, the flexibility to represent complex multi-part objects with holes. Multiresolution surface triangulation provides a general representation scheme which can adapt to the most complex free-form surfaces. 1

A preliminary version of this work has been presented at the 1992 IEEE Conference on Computer Vision and Pattern Recognition. 2 To whom correspondence should be addressed, E-mail: laurend@gel. ulaval.ca. 1

1077-3142/96 $18.00 Copyright  1996 by Academic Press, Inc. All rights of reproduction in any form reserved.

2

SOUCY AND LAURENDEAU

respect to the original triangulation cannot be computed. In the work of Hoppe et al. [6], coarser surface triangulations are obtained from an initial high-resolution triangulation by minimizing an energy function. Several parameters must be tuned by the user in order to control the optimization process. An important difference between Hoppe’s work and the one reported here is that he uses a quadratic error norm instead of a maximum error norm. We claim that the maximum error norm is more appropriate to control the degradation of a surface model. The top-down approach reported in [8] refines a rough approximation of a 3-D object by adding triangles to a triangular mesh. It is assumed that the topology of the modeled object can be recovered by fitting a deformable model on scattered 3-D data points. Hoppe et al. [7] have completed their previous work [6] by introducing a model compression technique that fits subdivision surfaces on the 3-D data instead of planes. Using subdivision surfaces improves the quality of surface fitting and reduces the amount of geometric information necessary to represent the model. As they mention however, a preliminary polygon reduction step based on a planar approximation is still useful in their approach since it is much faster than using subdivision surfaces. There are also several applications where a planar facet model is preferable, such as real-time simulation and virtual reality. This paper presents an iterative multi-resolution surface triangulation technique that computes coarse surface triangulations from an initial high-resolution surface triangulation. A sequential optimization process is implemented in order to remove at each iteration the vertex minimizing the retriangulation error. The aim of the sequential algorithm is thus to remove the maximum number of vertices while keeping the triangulation error as low as possible. The equiangularity of surface triangulations is optimized in 3-D space throughout the sequential process. The algorithm also preserves triangulation topology and orientation discontinuities on the surface. In this work, the initial surface triangulation results from the integration of multiple range views [19]. The technique presented here is, however, general and can be applied to any type of surface triangulation. The processed surface triangulations may form closed 3-D volumes but this condition is not mandatory. A surface triangulation embedded in 3-D space will be referred throughout this paper as a 2.5-D triangulation. Section 2 of the paper discusses the importance of optimizing the equiangularity of 2.5-D triangulations. It is demonstrated that a surface triangulation whose equiangularity is optimized in 3-D space maximizes the upper bound on the principal curvatures of the sampled surfaces. A technique to perform the equiangularity optimization is also presented. Section 3 describes the hierarchical triangulation algorithm used to represent a 3-D surface at different levels of resolution. It is shown that the proposed technique

preserves orientation discontinuities on the surface. A time complexity analysis of the multiresolution algorithm is also given. The paper concludes by presenting experimental results obtained on 2.5-D triangulations of complex multipart objects with holes. 2. OPTIMIZATION OF THE EQUIANGULARITY OF A 2.5-D TRIANGULATION

Few researchers have attempted to optimize surface triangulations embedded in 3-D space. At first sight, optimizing the equiangularity of a 2.5-D triangulation seems appropriate since it should intuitively eliminate elongated triangles leading to numerical instability [5]. This section introduces an even more important property of equiangular 2.5-D triangulations. It is demonstrated that equiangular 2.5-D triangulations maximize the upper bound on the principal curvatures of the sampled surfaces. This demonstration shows that optimizing the equiangularity of a 2.5-D triangulation is a good way to achieve the optimal sampling of a 3-D surface. The section also describes a technique to optimize the equiangularity of a 2.5-D triangulation by maximizing the smallest interior angle of each triangle. A. Goals of the Equiangularity Optimization in 3-D Space This section brings a mathematical justification to the optimization of the equiangularity of 2.5-D triangulations. A demonstration shows that minimizing the radii of the circumscribing circles of the triangles of a 2.5-D triangulation maximizes the upper bound on the principal curvatures of the sampled surfaces. Another advantage of equiangular 2.5-D triangulations is that they make the triangulation process invariant to viewpoint. Let us assume a triangle in space formed by points p1 , p2 , and p3 , as illustrated in Fig. 1. These three points define

FIG. 1. The shortest radius of sphere Rs for which the convex surface segment above circle C is entirely visible from the parametric grid of orientation Vi .

HIERARCHICAL TRIANGULATION MODELING

W (xN , yN , zN). The veca plane P with unit normal vector N tor Vi , which is orthogonal to the planar parametric grid onto which the triangle is to be projected (the x–y plane in Fig. 1), is parallel to the z-axis. The circumscribing circle C of the triangle has center pC and radius RC . This circle lies on plane P. The circle fitted on any triplet of points on a spherical surface is also part of the surface. Hence, a spherical surface going through p1 , p2 , and p3 must comprise circle C. In order to find an upper bound on the principal curvatures of the sampled surfaces, we compute the shortest radius of sphere RS for which the convex surface segment above circle C is entirely visible from the viewpoint of the parametric grid of orientation Vi . For any sphere viewed by a parametric grid orthogonal to the zaxis, as in Fig. 1, the occluding contour of the sphere is an equatorial circle parallel to the x–y plane. As demonstrated in Fig. 1, the shortest radius of sphere for which the convex surface segment above circle C is entirely visible is the one for which circle C and the occluding contour of the sphere intersect. This implies that the normal to the sphere at the point on circle C having the smallest z-coordinate should be orthogonal to the z-axis. This point on circle C is called pl(xl , yl , zl) as illustrated in Fig. 1. A second geometrical constraint is that the center of the sphere ps(xs , ys , zs) is located on a straight line orthogonal to plane P and going through the center pC of circle C, as shown again in Fig. 1. The parametric equation of this line is x 5 xc 1 kxN y 5 yc 1 kyN

(1)

z 5 zc 1 kzN . Now zs , the z-coordinate of the center of the sphere, is equal to zl , since the radius joining zl to zs is orthogonal to the z-axis. The coordinates of the center of the sphere are thus

SD SD

xs 5 xc 1 (zl 2 zc)

xN zN

ys 5 yc 1 (zl 2 zc)

yN zN

(2)

zs 5 zl . From Eq. (2), the radius of the sphere can be computed with the procedure given in the Appendix. We simply state the final result which is Rs 5

Rc . zN

(3)

3

The radius Rs of the sphere is equal to the radius of the circumscribing circle divided by the z-component of the unit normal to plane P. Radius Rs is the shortest radius of sphere for which the convex surface segment above circle C is entirely visible from the planar parametric grid of normal Vi . Interestingly enough, Rs is also the shortest radius of the sphere for which the concave surface segment below circle C is entirely visible from the parametric grid. Therefore, any surface segment for which the absolute value of the principal curvatures is smaller than l/Rs at all points ‘‘inside’’ circle C will be visible from the planar parametric grid of normal Vi . By minimizing the radius Rs , one maximizes the upper bound on the principal curvatures of the observed surfaces. Radius Rs may be minimized by either minimizing Rc or maximizing zN . The maximization of zN consists in sampling the surfaces with a grid for which the normal is as close as possible from the surface normal. As for the minimization of Rc , it is obtained by optimizing the equiangularity of the 2.5-D triangulation, as shown in Section 2.B. Dependence to viewpoint is a major problem that restrain the use of 2-D triangulation techniques to triangulate a set of 3-D surface points. Indeed, an optimal equiangular 2-D triangulation, like the Delaunay triangulation, computed in 2-D on a set of 3-D points projected on a plane may lead to different results depending on the orientation of the plane. By optimizing the equiangularity of the surface triangulation in 3-D, the invariance to viewpoint could be theoretically obtained if the optimization function were convex. Unfortunately, as shown below, the technique designed to optimize the equiangularity of a surface triangulation in 3-D is nonconvex. Invariance to viewpoint can thus only be partially achieved. This still gives the present technique an advantage compared to 2-D techniques which do not consider invariance at all. B. A Technique to Optimize the Equiangularity of a 2.5-D Triangulation This section describes a technique to optimize the equiangularity of a surface triangulation embedded in a 3-D space. An iterative algorithm maximizes the smallest interior angle of each triangle of the 2.5-D triangulation. It has been demonstrated in [20] that maximizing the smallest interior angle of a set of 3-D triangles also minimizes the radii of the circumscribing circles of these triangles. In order to optimize a 2.5-D triangulation T, one must be able to modify its triangles. Given a 2.5-D manifold surface triangulation T, a 2.5-D retriangulation is defined as an operation by which a subset of triangles of T is deleted and replaced by another set of triangles. Clearly, a 2.5-D retriangulation must preserve the topology of the 2.5-D triangulation T. To do so, we introduce a constraint for a 2.5-D retriangulation to be attempted. When this

4

SOUCY AND LAURENDEAU

constraint is met, the topology of the surface triangulation is guaranteed to be preserved locally. Let us assume a set of points P in three-dimensional space. Let T be a 2.5-D manifold triangulation of P. A subset of triangles ST of T can be retriangulated if a one-to-one mapping can be defined between ST and a set of triangles SM defined in a 2-D parametric space.

The retriangulation constraint states that the subset of 3-D triangles ST has to be mapped on a 2-D parametric space and that the resulting 2-D triangles of the set SM must not overlap in this 2-D space. Under these conditions, the area covered by SM can be retriangulated using a standard 2-D constrained triangulation technique. The resulting triangles can then be backprojected in 3-D space to replace the triangles in ST . The one-to-one mapping between ST and SM insures that the triangulation topology is preserved since a constrained retriangulation in 2-D space does not modify the topology of the set SM of mapped triangles. In order to optimize the equiangularity of a 2.5-D triangulation, a technique similar to Lawson’s [9] has been devised in a 3-D space. In 2-D space, Lawson proposed to maximize the smallest interior angle of the triangles in order to obtain the most equiangular triangulation. At each iteration, all pairs of adjacent triangles forming a convex quadrilateral are found and the smallest of the six interior angles of the triangles is computed. The diagonal common to both triangles is then permuted. This operation consists in replacing the edge joining the two common vertices by an edge joining the two noncommon vertices. As demonstrated in Fig. 2, the permutation of the common diagonal is only possible when the quadrilateral formed by both triangles is convex. Indeed, when the quadrilateral is not convex, the new diagonal is located outside of the space to be triangulated. If the smallest interior angle of the new triangles obtained by the permutation of the common diagonal is larger than the smallest interior angle of the

previous pair of triangles, the triangulation is modified to take into account the permutation of the diagonal. Sibson [17] has demonstrated that the maximization of the smallest interior angle of the triangles leads to the 2-D Delaunay triangulation. In a 3-D space, the equiangularity of the surface triangulation is optimized by a series of 2.5-D retriangulations which aim to maximize the smallest interior angle of the triangles in 3-D space. When a pair of adjacent triangles is processed, a plane whose orientation averages the orientation of the normals to the triangles is defined. The two 3-D triangles are then mapped orthogonally on this plane. If the two mapped triangles do not overlap and the quadrilateral formed by the mapped triangles is convex, the common diagonal may be permuted. When the smallest interior angle of the 3-D triangles is increased by this permutation of the diagonal, the global 2.5-D triangulation is modified locally by making the permutation effective. This retriangulation process is repeated until the smallest interior angle can no longer be increased. The convergence of this iterative process is guaranteed since the smallest interior angle cannot exceed 608. In 2-D, the process of maximizing the smallest interior angle of the triangles is a concave function. Concave and convex functions share the property that any local maximum or minimum is also a global maximum or minimum [10]. Assuming a set of points on a plane and any initial 2-D triangulation, the optimization process will always lead to a Delaunay triangulation [17]. Unfortunately, in 3-D, the maximization process is not a concave function. With different initial triangulations, different results are obtained at convergence. It is thus possible that the optimization process leads to a local maximum which is not the global maximum. Consequently, it is desirable that the initial 2.5-D triangulation be a good approximation of an equiangular triangulation in 3-D before applying the iterative optimization technique. A good initial approximation should result in a local maximum close to the sought global maximum. In the particular case of this work, 2-D Delau-

FIG. 2. Permutation of the diagonal common to two adjacent triangles. (a) Convex quadrilateral. (b) Nonconvex quadrilateral.

HIERARCHICAL TRIANGULATION MODELING

nay triangulations are already used to optimize subsets of the integrated 2.5-D triangulation [19]. In the more general case, a good initial approximation can be obtained by segmenting the 2.5-D triangulation in small connected subsets of triangles. For each subset of triangles, a retriangulating plane whose orientation approaches the orientations of the normals to the triangles is defined. A good approximation of an equiangular triangulation can then obtained by mapping the subsets on their retriangulating plane and optimizing their 2-D equiangularity. 3. MULTIRESOLUTION SURFACE MODELING ALGORITHM

This section presents the algorithm that computes coarser levels of resolution of a high-resolution 2.5-D triangulation. The criterion that is used to determine the level of compression is the maximum acceptable error of the coarse model with respect to the initial triangulation. Given a maximum error level E that bounds the triangulation error, the algorithm aims at computing the most equiangular 2.5-D triangulation having the smallest number of vertices. A sequential optimization strategy has been devised to reach that goal. At each iteration, the vertex that minimizes the retriangulation error with respect to the initial triangulation is effectively removed from the triangulation. The main advantage of using sequential optimization is that it allows the implementation of a deterministic geometric compression technique that runs in predictable computation time. This section describes in details the sequential optimization process. It is also demonstrated that the proposed method preserves surface orientation discontinuities. Finally, a time complexity analysis of the algorithm is given. A. Hierarchical 2.5-D Triangulation Algorithm A sequential optimization process is an interative process which optimizes a given criterion at each step in order to reach a goal. The criterion optimized by the hierarchical triangulation algorithm is the minimization of the retriangulation error. At each iteration, the vertex that minimizes the retriangulation error is thus effectively removed from the triangulation. This process is interesting since it creates a hierarchy of triangulations ordered by increasing approximation error. The idea of using sequential optimization to build a hierarchy of triangulations has been proposed by De Floriani [5]. However, there are major differences between the work of De Floriani and the one reported here. De Floriani implemented a hierarchical 2-D Delaunay triangulation algorithm that is used to model surfaces that can be projected on a plane. Equiangularity is optimized in 2-D and error vectors are computed along a direction normal to the parametric plane. In the present work, the sequential optimization approach is used to model com-

5

plex surfaces with holes embedded in 3-D space. The equiangularity of the surface triangulation is optimized in 3-D and all error vectors correspond to actual 3-D distances between geometric primitives. As shown in the following, the preservation of surface orientation discontinuities is embedded in the algorithm and is thus performed automatically. In De Floriani’s work, edges to be preserved must be selected by the user. Here is the multiresolution algorithm in a pseudocode form: /* Initialization of the algorithm */ for (i 5 0; i , number vertices; i11) h compute retriangulation vertex i; compute retriangulation error vertex i; j built list ordered vertices; /* Sequential optimization process */ while (maximum error , THRESHOLD) h find vertex minimizing error; find neighbors vertex; remove vertex; delete ancient triangles; inscribe new triangles; for (i 5 0; i , number neighbors; i11) h compute retriangulation neighbor i; compute retriangulation error neighbor i; modify list of vertices; j j At the initialization step of the algorithm, the retriangulation error of each vertex of the 2.5-D triangulation is computed and stored. This computation is achieved by removing the vertex, retriangulating the area affected by this removal, and then computing the minimum 3-D distance between the removed vertex and the local 2.5-D retriangulation. Once all errors are evaluated, a list of all vertices ordered by increasing retriangulation error is built. The sequential optimization process then begins. At each iteration, the vertex minimizing the retriangulation error is effectively removed and the 2.5-D retriangulation resulting from this removal is inscribed in the global 2.5-D triangulation. The retriangulation errors of all neighbors of the removed vertex are then recomputed. Indeed, the local retriangulation caused by the removal of a vertex modifies the retriangulation error of its neighbors in the triangulation. Parameter THRESHOLD is a 3-D distance. The algorithm stops when the smallest retriangulation error of a vertex becomes larger than parameter THRESHOLD which is defined as a 3-D distance. At this phase of the algorithm, the global 2.5-D triangulation has an approxi-

6

SOUCY AND LAURENDEAU

mation error bounded by THRESHOLD. All vertices of the initial triangulation are located at a distance smaller than THRESHOLD from the planar facet model fitted on the 2.5-D triangulation. In the actual implementation, several increasing values for THRESHOLD can be selected by the user. Each time the retriangulation error exceeds one of the values, the 2.5-D triangulation is saved in a file. B. Computation of the Retriangulation Error for an Interior Vertex There are two categories of vertices in the 2.5-D triangulation: interior vertices and contour vertices. Contour vertices are the vertices of edges which are part of only one triangle and which represent the frontiers of the 2.5-D triangulation. No surface information is available beyond the contour edges. Vertices which are not contour vertices are interior vertices of the 2.5-D triangulation. The computation of the retriangulation error for an interior vertex is performed in two steps. The interior vertex is first tentatively removed and the region affected by this removal is retriangulated in 2.5-D. Then, the retriangulation error is evaluated. Figure 3 depicts the retriangulation procedure. Vertex pe , which is a vertex of triangles T1 to T6 , is removed. As stated in Section 2.B, there must exist a 2-D parametric space onto which triangles T1 to T6 can be projected without overlap. Therefore, a triangulating plane is defined whose normal vector is the estimated surface normal at vertex pe . Triangles T1 to T6 are then orthogonally projected on this plane. If the projected triangles do not overlap in 2-D, it is allowed to retriangulate this region in 2-D. The 2-D retriangulation is obtained by computing a constrained Delaunay triangulation of the projected vertices p1 to p6 . Edges p1 –p2 , p2 –p3 , p3 –p4 , p4 –p5 , p5 –p6 , and p6 –p1 are imposed as edges of the 2-D retriangulation. Thereafter, triangle p1 –p5 –p6 , which is located outside the

region formed by triangles T1 to T6 on the triangulating grid, is deleted and the equiangularity of the 2-D retriangulation is optimized in 3-D using the algorithm described in Section 2.B. The final retriangulation is formed by triangles R1 to R4 . Once the retriangulation is obtained, it is possible to evaluate the 3-D error caused by the removal of pe . The error of vertex pe is by definition the smallest distance between pe and the planar facet model fitted on the global 2.5-D triangulation. In order to estimate this error, the smallest 3-D distance between triangles of the 2.5-D retriangulation and vertex pe is computed. This assumes that the triangle of the retriangulation nearest to pe is also the triangle of the global triangulation nearest to pe in 3-D. This estimate of the error of pe cannot be considered as the retriangulation error. There exists a set of vertices PRi , removed at previous iterations, for which the nearest triangles in the global 2.5-D triangulation are part of triangles T1 to T6 . Since these triangles are deleted by the removal of pe , it is mandatory to reevaluate the error of the vertices PRi by finding their nearest triangle in the retriangulation. Let errmax be the largest error among the error of pe and the errors of vertices PRi . The retriangulation error is defined as errmax. In Fig. 3, 10 removed vertices, represented by small black circles, are approximated by a triangle part of triangles T1 to T6 . The error of these 10 vertices and the error of pe are estimated by finding the nearest triangle in the retriangulation. The largest of these 11 computed errors is considered as the retriangulation error resulting from the removal of vertex pe . C. Computation of the Retriangulation Error for a Contour Vertex Three constraints must be met to retriangulate a contour vertex pec . As for the case of an interior vertex, a triangulating plane has to be defined and the set of 3-D triangles Ti having pec as a vertex must not overlap once projected

FIG. 3. Computation of the retriangulation error of interior vertex pe . (a) Initial triangulation. (b) 2-D triangulation. (c) 2.5-D triangulation.

HIERARCHICAL TRIANGULATION MODELING

FIG. 4. Contour vertices which are part of (a) a junction or (b) a triangular hole, cannot be removed from the retriangulation.

on this plane. Two additional topological constraints are imposed for a contour vertex. Vertex pec must only be part of two contour edges and the two other contour vertices pec2 and pec1 linked to pec must not form a contour edge, as depicted in Fig. 4. These constraints avoid the removal of vertex junctions and holes in the triangulation. When these constraints are met, a 2.5-D retriangulation is attempted. However, a retriangulation is not always possible. A retriangulation following the removal of a contour vertex is possible if and only if all interior vertices remain interior vertices and all contour vertices remain contour vertices. This constraint is called ‘‘the border constraint.’’ By replacing the two contour edges pec2 –pec and pec1 – pec by a contour edge jointing pec2 to pec1 , we insure that contour vertices remain contour vertices. To insure that interior vertices remain interior vertices, it is sufficient to verify that no vertex of Ti falls inside the triangle formed by pec2 , pec , pec1 in the 2-D triangulating plane. As shown in Fig. 5a, this condition is always met for concave contours.

7

The convex contour in Fig. 5b presents a case where retriangulation is not possible since an interior vertex falls inside triangle pec2 , pec , pec1 . If pec were removed, this interior vertex would be located outside the triangulation. On the other hand, the convex contour in Fig. 5c) satisfies the border constraint. If the border constraint is respected, the 2.5-D retriangulation may be computed. Two cases of retriangulation are possible, as illustrated in Fig. 6. If the set Ti of triangles having pec as a vertex is a singleton, no retriangulation has to be computed. The single triangle having pec as a vertex is simply deleted. On the other hand, if the set Ti contains several triangles, the set pi of vertices of Ti other than pec , is found. The retriangulation is then estimated by a 2-D constrained Delaunay triangulation of vertices pi computed on the triangulating plane. Frontier edges of the region formed by Ti which do not have pec as a vertex and the new edge joining pec2 and pec1 are imposed as edges of the retriangulation. Thereafter, the equiangularity of this retriangulation is optimized in 3-D. The retriangulation error of contour vertex pec can then be evaluated. In order to constrain the frontiers of the triangulation, the error of a contour vertex pec is defined as the smallest distance between a contour edge of the global 2.5-D triangulation and pec . This error is estimated by computing the 3-D distance between vertex pec and the new contour edge joining pec2 and pec1 . In order to compute the actual retriangulation error associated to pec , the errors of all previously removed vertices, forming the set PRi , which were approximated by triangles in Ti must be reevaluated. For all vertices of PRi which are interior vertices, the error is estimated as the smallest distance between the vertex and the retriangulation. For all contour vertices of PRi , the distance between the vertex and the contour edge pec2 –pec1 is computed instead. The largest of the errors of vertices in PRi and the error of pec is found. This maximal error is the retriangulation error of pec .

FIG. 5. Interior vertices must stay inside of the triangulation when a contour vertex is removed. This constraint is always verifiied when the contour is concave as in (a). The first example of a convex contour in (b) shows a case where retriangulation is not possible while the second example in (c) verifies the constraint.

8

SOUCY AND LAURENDEAU

D. Preservation of Surface Orientation Discontinuities An interesting property of the hierarchical triangulation algorithm is the preservation of orientation discontinuities, e.g., edges on the surface. This property is the result of the sequential nature of the algorithm. For the sake of simplicity, the demonstration will be illustrated in 2-D. The results will then be extended to the 3-D case. The situation in 2-D is illustrated in Fig. 7. Point e is located on an orientation discontinuity. Each side of the discontinuity can be locally represented by a circular arc. All points are measured with an error bounded by Max(err( p)). The 2.5-D triangulation in space becomes a representation by linear segments in the 2-D case. Let us assume two points pa and pb located on a circular arc of radius R as shown in Fig. 8. The error resulting from the approximation of arc AC by a line segment joining pa and pb is bounded by distance «AC between the middle point of the line segment and arc AC. This distance «AC is «AC 5 R 2

!

d 2( pa , pb) R 2 . 4 2

(4)

equal to p21 while p r is p1 . For i 5 0, where i is the number of removed points, we have «la (e)0 5 d(e, p21 p) «la ( p r )0 5 d( p1 , ep2) «la ( p )0 5 d( p21 , ep22). The initial approximation error of a point is the distance between this point and the segment joining the predecessor and the successor of the point. Point e is kept in the linear approximation if «la (e)0 . «la ( p l )0 or

(5)

The bounded interval in (5) is obtained by considering the case where the error vectors on the measured points pa and pb are oriented in a direction orthogonal to the segment pa –pb , and all error vectors of the other measured points on the arc are oriented in the opposite direction. This situation is illustrated in Fig. 8. In Fig. 7, the linear approximation consists of six line segments joining p23 , p22 , p21 , e, p1 , p2 and p3 . The function «la ( p) represents the approximation error caused by the removal of point p. Let p l be the point of the linear approximation closest to the left-hand side of e and p r be the point of the linear approximation closest to the right-hand side of e. At the beginning of the hierarchical algorithm, p l is

FIG. 6. Two cases of retriangulations. In (a), Ti contains one triangle, while in (b), Ti contains several triangles.

(7)

«la (e)0 . «la ( p r )0 and that the following condition is verified for each iteration i for a given level of model error «m : («la (e))i . «la ( p il ) or («la (e))i . «la ( p ri ) or («la (e))i . «m .

The error « resulting from the approximation of arc AC by a line segment is bounded by «AC 1 2 Max(err( p)): « # «AC 1 2 Max(err( p))

(6)

l

(8)

That is, at each iteration, the error caused by the removal of e must be larger than the three following errors: the error caused by the removal of the neighbor on its left, the one caused by the removal of the neighbor on its right, and the threshold of model error. The error «la (e)0 depends on the amplitude of the orientation discontinuity. Errors «la ( p l )0 and «la ( p r )0 are bounded: «la ( p l )0 # Rl 2

!

«la ( p )0 # Rr 2

!

r

R 2l 2 R 2r

d 2( p22 , e) 1 2 Max(err( p)) 4

d 2( p2 , e) 2 1 2 Max(err( p)). 4

(9)

FIG. 7. A curve comprising an orientation discontinuity is modeled by two circular arcs.

9

HIERARCHICAL TRIANGULATION MODELING

FIG. 8. Computation of the maximal distance between the circular arc and the line segment approximating it.

Rl and Rr are the radii of curvature of the circular arcs modeling each side of the discontinuity. Therefore, for Eq. (7) to be verified, the amplitude of the orientation discontinuity must be more significant than the noise level. Furthermore, the radii of curvature Rl and Rr on each side of the discontinuity must be sufficiently large compared to distances d( p22 , e) and d( p2 , e). In practice, this means that surface segments of large curvature need a finer sampling. As for the conditions in Eq. (8), they benefit from a phenomenon of sequential strengthening of the orientation discontinuity. This phenomenon is depicted in Fig. 9. When p l and p r get away from e, the error caused by the removal of e increases. If Eq. (7) is verified, p21 or pl will disappear before e. The discontinuity is then strengthened. If p l or p r can again be removed before e, the discontinuity is still strengthened, and so on. The critical phase for the preservation of orientation discontinuities is at the beginning of the sequential optimization process. If the error caused by the removal of e is sufficiently large, the removal of neighboring points gradually strengthens the discontinuity such that point e will remain in the linear approximation at the following iterations. In 3-D space, an edge on a surface corresponds to the locus of orientation discontinuities between two connected surface segments. Results obtained in 2-D can be extended to the 3-D space. If the curvature of connected surface segments A and B located on each side of the surface edge is not too large compared to the sampling step and if the orientation discontinuities are more significant than the

FIG. 9. Sequential strengthening of an orientation discontinuity.

acquisition noise, points on both side of the surface edge will disappear, strengthening the points nearest to the orientation discontinuities. A new phenomenon then appears. A point e near an orientation discontinuity may be removed if the nearest triangle in the retriangulation caused by the removal of e comprises an edge joining the predecessor and the successor of e on the surface edge. This phenomenon consists in a compression of the edge information. A good approximation of an edge point can only come from two triangulation points neighboring the surface edge. Therefore, the model of a surface edge will be more and more rough as the triangulation error increases. The edge will be represented by a set of vertices approximating the curve which is the locus of orientation discontinuities. The edge preservation property justifies the use of sequential optimization to build a hierarchy of 2.5-D triangulations. Indeed, this property relies heavily on the strengthening phenomenon which is purely sequential in nature. An algorithm which would attempt to remove sets of points simultaneously at each iteration could not guarantee the preservation of surface edges since the strengthening mechanism embedded in the sequential process would not be present.

E. Time Complexity of the Hierarchical Triangulation Algorithm This section presents a time complexity analysis for the sequential optimization process which is the core of the hierarchical triangulation algorithm. An experimental result confirms the validity of the computed function. The average computation time of a 2.5-D local retriangulation is constant at every step of the hierarchical triangulation algorithm. The computation time of the retriangulation error however increases gradually. Let tdpt be the time needed to compute the distance between a vertex and a triangle. Assuming that the retriangulation following the removal of vertex p comprises Nr triangles, the time required to compute the retriangulation error of each removed vertex is equal to tdpt times Nr . To obtain the global retriangulation error, the maximal error among the error of vertex p and the errors of previously removed vertices, which were approximated by triangles having p as a vertex, must be found. Let N be the number of vertices of the initial 2.5-D triangulation and let n be the number of removed vertices. It has been demonstrated in [20] that the number of removed vertices Nrv located in the area affected by a 2.5-D retriangulation can be estimated as follows:

Nrv 5

N . N2n

(10)

10

SOUCY AND LAURENDEAU

FIG. 10. Comparison between theoretical and experimental computation times. The left curve is experimental. The number of removed vertices is on the x-axis, while the y-axis represents the number of seconds. The right graphic shows the theoretical function fitted on the experimental data.

Therefore, the computation time of a retriangulation error for a given value of n is

tre 5

Nrtdpt N . N2n

(11)

If all other computation times are neglected, the time needed to remove a vertex for a given value of n is

tn 5

Nt Nrtdpt N . N2n

(12)

From Eq. (12), it is possible to find a function t(n) which models the computation time needed to go from level 0, where no vertex has been removed, to level n, where n vertices have been removed:

t(n) 5

E t (n) dn 5 N N t n

0

n

t

r dpt N

log

S D

N . N2n

4. EXPERIMENTAL RESULTS

Experimental results of the implemented algorithm are presented in Figs. 11 to 14. The initial triangulated models of the teapot and toy soldier objects were built by integrating multiple range views [19, 20]. The range views were acquired using a synchronized scanners rangefinder developed by the National Research Council of Canada [13]. In the first experiment, three coarse triangulations having maximum errors equal to 0.5, 2.0, and 5.0 mm are computed from an initial surface triangulation of a teapot, shown in Fig. 11. The three resulting compressed triangulations are depicted in Fig. 12 and comprises respectively 3278, 846, and 378 triangles. For a maximum error as low as 0.5 mm, the compression factor on the number of triangles reaches 94%. It can be clearly observed that the topology of the teapot is preserved throughout the compression of

(13)

The function t(n) depends on log(N/(N 2 n)), where N is the number of vertices of the initial triangulation. In Fig. 10, the plot on the left shows an experimental curve while the one on the right illustrates the function of Eq. (13) fitted on the experimental data. The similarity of both graphics supports the validity of the theoretical modeling of time complexity. The part of the algorithm that is most time consuming is thus the computation of the retriangulation errors. The rapid growth of this computation time results from the fact that the errors of all previously removed vertices which were locally approximated by the last deleted triangles have to be reevaluated.

FIG. 11. Flat shaded teapot model built from eight range views (53,821 triangles).

HIERARCHICAL TRIANGULATION MODELING

11

FIG. 12. Coarse teapot triangulations computed by the hierachical triangulation algorithm. (a) Wireframe model of error 0.5 mm (3278 triangles). (b) Wireframe model of error 2.0 mm (846 triangles). (c) Wireframe model of error 5.0 mm (378 triangles).

the initial triangulation. Similarly, orientation discontinuities on the surface of the teapot are preserved and compressed as the allowable maximum error increases. The second experiment is performed on a triangulation of the toy soldier shown in Fig. 13, which was built from eight range views. Two compressed triangulations having maximum errors equal to 0.5 and 2.0 mm are depicted in Fig. 14 using a wireframe and flat shading representation. These two models are made respectively of 5262 and 1449 triangles. These results clearly demonstrate that the hierarchical triangulation algorithm presented in this paper can be used to significantly compress triangulations of complex multipart objects with holes while preserving its essential surface properties. It is worth noting that no hidden tuning parameter has been set to obtain these experimental results. The multiresolution technique is fully autonomous.

5. CONCLUSION

The paper has first described a method to optimize the equiangularity of a surface triangulation embedded in 3-D space. A demonstration has shown that optimizing the equiangularity of a 2.5-D triangulation maximizes the upper bound on the principal curvatures of the sampled surfaces. The main contribution of the paper is a new multiresolution surface modeling technique based on a hierarchical triangulation algorithm. Coarse levels of resolution of an initial high-resolution surface triangulation are obtained by iteratively removing the vertex minimizing the retriangulation error. The proposed algorithm is fully autonomous, optimizes the equiangularity of the triangulation in 3-D, preserves the topology of the initial triangulation and preserves surface orientation discontinuities. Experimental results have shown that the method can

12

SOUCY AND LAURENDEAU

FIG. 13. Flat shaded toy soldier model built from eight range views as viewed from two different viewpoints (70,842 triangles).

FIG. 14. Coarse toy soldier triangulations computed by the hierarchical triangulation algorithm. (a) Wireframe model of error 0.5 mm (5262 triangles). (b) Flat shaded model of error 0.5 mm (5262 triangles). (c) Wireframe model of error 2.0 mm (1449 triangles). (d) Flat shaded model of error 2.0 mm (1449 triangles).

13

HIERARCHICAL TRIANGULATION MODELING

successfully process complex multipart objects with holes. The presented technique will be very useful in numerous 3-D visualization applications. For instance, in the near future, it will be possible to digitize museum artifacts or precious objects and transmit surface models over communication networks. In a visualization context, it is often more important to render the overall shape than to display the finest details. Furthermore, an efficient transmission of the models on communication networks necessitates an important compression of the raw data. The hierarchical triangulation algorithm presented in this paper has demonstrated great compression capabilities combined with excellent surface approximation. In future work, we plan to fit surface patches on the triangulations and interface these models to existing CAD software. APPENDIX: COMPUTATION OF THE RADIUS RS OF THE SPHERE OF SECTION 2.A

From Eq. (2), radius Rs can be computed: Rs 5 Ï(xs 2 xl )2 1 ( ys 2 yl )2 1 (zs 2 zl )2 5

!Az1 B 2 N

Since point pl is the point of circle C for which the z coordinate is the smallest, it is located in the vertical plane formed by the z axis and the vector normal to plane P. The absolute value of (zl 2 zc ) is thus equal to: uzl 2 zc u 5 Rc Ïx 2N 1 y 2N

(19)

Therefore, the value of Rs is: Rs 5

!

R 2c (z 2N 1 x 2N 1 y 2N ) Rc 5 z 2N zN

(20)

ACKNOWLEDGMENTS The authors thank M. Rioux and J. Domey of the National Research Council of Canada (NRCC) for providing the range images used in this work. The research work has been supported by NSERC Grant OGPIN 011. M. Soucy was supported by a NSERC scholarship. The authors are members of the Institute for Robotics and Intelligent Systems (IRIS) and acknowledge the support of the Networks of Centres of Excellence program of the Government of Canada, the Natural Sciences and Engineering Research Council, and the participation of PRECARN Associates Inc. Thanks also to R. Bergevin and the reviewers for many comments that helped to improve the quality of the manuscript.

A 5 z 2N (xc 2 xl )2 1 2zN xN (xc 2 xl )(zl 2 zc ) 1 x 2N (zl 2 zc )2

REFERENCES

B 5 z 2N ( yc 2 yl )2 1 2zN yN ( yc 2 yl )(zl 2 zc ) 1 y 2N (zl 2 zc )2

1. G. Baszenski and L. L. Schumaker, Use of simulated annealing to construct triangular facet surface, in Curves and Surfaces (P.-J. Laurent, A. Le Me´haute´, and L. L. Schumaker, Eds.), pp. 27–32, Academic Press, San Diego, 1991. 2. A. Blake and A. Zisserman, Visual Reconstruction, MIT Press, Cambridge, MA, 1987. 3. L. De Floriani, B. Falcidieno, G. Nagy, and C. Pienovi, A hierarchical structure for surface approximation, Comput. Graphics 8, 1984, 183–193. 4. L. De Floriani and E. Puppo, Constrained Delaunay triangulation for multiresolution surface description, in Proc. of the 9th Intl. Conf. on Pattern Recognition, Roma, Italy, Nov. 1988, pp. 566–569. 5. L. De Floriani, A pyramidal data structure for triangle-based surface description, IEEE Comput. Graphics Appl., Mar. 1989, 67–78. 6. H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, and W. Stuetzle, Mesh optimization, in ACM Computer Graphics Proceedings, 1993, pp. 19–26.

(14) Since we have: z 2N ((xc 2 xl )2 1 ( yc 2 yl )2) 1 (x 2N 1 y 2N )(zl 2 zc )2 5 R 2c z 2N 1 (zl 2 zc )2(1 2 2z 2N ),

(15)

the expression of Rs can be simplified.

Rs 5

!

R 2c z 2N 1 (zl 2 zc )2 1 (zl 2 zc )[22(zl 2 zc )z 2N 1 2xN zN (xc 2 xl ) 1 2yN zN ( yc 2 yl )] z 2N (16)

Rs 5

!

R 2c z 2N 1 (zl 2 zc )2 1 2zN (zl 2 zc ) [(zc 2 zl )zN 1 xN (xc 2 xl ) 1 yN ( yc 2 yl )] z 2N

(17)

The expression between brackets is the dot product between the vector joining pl to pc , located in plane P, and the normal vector to plane P. Since this dot product is zero, Rs becomes: Rs 5

!

R 2c z 2N

1 (zl 2 zc ) z 2N

2

(18)

7. H. Hoppe, T. DeRose, T. Duchamp, M. Halstead, H. Jin, J. McDonald, J. Schweitzer, and W. Stuetzle, Piecewise smooth surface reconstruction, in Proc. of SIGGRAPH 94, pp. 295–302. 8. E. Koh, D. Metaxas, and N. Badler, Hierarchical shape representation using locally adaptive finite elements, in Proc. of ECCV 94, pp. 441–446. 9. C. L. Lawson, Generation of a Triangular Grid with Application to Contour Plotting, Tech. Memo. 299, California Institute of Technology Jet Propulsion Lab., 1972. 10. D. G. Luenberger, Linear and Nonlinear Programming, 2nd ed., Addison–Wesley, Reading, MA, 1984. 11. M. Polis and D. McKeown Jr., Iterative TIN generation from digital elevation models, in Proc. of IEEE Conf. on CVPR, Urbana– Champaign, IL, 1992, pp. 787–790. 12. J. Ponce and O. Faugeras, An object centered hierarchical representa-

14

13. 14. 15.

16.

SOUCY AND LAURENDEAU tion for 3D objects: The prism tree, Comput. Vision Graphics Image Process. 38, Apr. 1987, 1–28. M. Rioux, Laser range finder based on synchronized scanners, Appl. Opt. 23, Nov. 1984, 3837–3844. H. Samet, The Design and Analysis of Spatial Data Structures, Addison–Wesley, Reading, MA, 1989. F. Schmitt, X. Chen, W.-H. Du, and F. Sair, Adaptive G1 approximation of range data using triangular patches, in Curves and Surfaces (P.-J. Laurent, A. Le Me´haute´, and L. L. Schumaker, Eds.), pp. 433–436, Academic Press, San Diego, 1991. W. J. Schroeder, J. A. Zarge, and W. E. Lorensen, Decimation of triangle meshes, ACM Comput. Graphics 26(2), July 1992, 65–70.

17. R. Sibson, Locally equiangular triangulation, Comput. J. 21, 1978, 243–245. 18. M. Soucy, A. Croteau, and D. Laurendeau, A multi-resolution surface model for compact representation of range images, in Proc. of IEEE Intl. Conf. on Robotics and Automation, Nice, France, 1992, pp. 1701– 1706. 19. M. Soucy and D. Laurendeau, Multi-resolution surface modeling from multiple range views, in Proc. of IEEE Conf. on CVPR, Urbana– Champaign, IL, 1992, pp. 348–353. 20. M. Soucy, Modeling the Surfaces of 3-D Objects Using Multiple Range Views, Ph.D. Thesis, Laval University, Dec. 1992. [In French]