Piecewise B-spline Surfaces Fitting to Arbitrary Triangle Meshes

Piecewise B-spline Surfaces Fitting to Arbitrary Triangle Meshes L.Y. Zhang, R.R. Zhou, J.Y. Zhu(l), X. Wu College of Mechanical and Electrical Engineering, Nanjing University of Aeronautics 8, Astronautics, Nanjing 210016, P. R. of China Abstract Due to technical progress in shape acquiring from range data, triangle mesh representations of models can now be obtained conveniently. However, triangle meshes are not the final representation form in many applications, for reasons of compactness, appearance, modification, and manufacturability etc. In this paper, we present a solution for fitting smooth B-spline surfaces to triangle meshes of arbitrary topology. The procedure consists of two major phases: (1) partitioning triangle meshes of arbitrary topology into quadrilateral patches, (2) fitting tensor-product B-spline surfaces to each quadrilateral patch and keeping patches continuous across shared boundaries and corners. One merit of our solution is that it makes a reasonable tradeoff between automation and human control in complex model segmentation and surface fitting. Experiments in the paper show that our solution is practical and effective. Keywords: Reverse engineering, Surface reconstruction, Spline

1 INTRODUCTION Triangle mesh representation of models has many advantages, especially the flexibility to model arbitrary topology type. Due to technical progress in shape acquiring from range data, triangle mesh representation of models can now be obtained conveniently. However, B-spline(or more generally NURBS, Non Uniform Rational B-Splines) representation is even more attractive for engineering design, due to its characteristics of compact, accurate representation and global parametrization, and becomes the de facto standard representation of surfaces in almost all commercial CAD packages. But a single tensor product B-spline surface can model only rectangle topology (the cylinder and torus can be considered to be special cases of rectangles). Techniques are needed for fitting B-spline surfaces to triangle meshes of arbitrary topology. To solve this problem, quadrilateral patches network of the mesh must be constructed first. This is generally known to be difficult. Second, an algorithm should be employed to fit piecewise B-spline surfaces to the quadrilateral network, keeping patches continuous across shared boundaries and corners. There has been certain work on segmentation of arbitrary triangle meshes [I]-[7]. Fischer[l] put forward an approach for extracting a continuous 3D boundary of an open surface by using image processing techniques. Her method, however, can not be used to the reconstruction of models with arbitrary topologies. Guo [3] and Bajaj [4] made use of 3D a-shape [8] ,which is quite expensive as Guo [3] pointed out, to capture the topology of the object. Krishnamurthy [5], Milroy [6] etc. implemented the model segmentation of arbitrary topology by manually planning and interactively specifying boundary points or boundary curves. Unfortunately, it is a difficult task for users to manually plan and delineate all the corners and the boundary curves of the segmented patches, especially for complex or wraparound models. It is quite difficult to keep the overall segmentation consistent. In contrast, Eck [7] presented a scheme to automatically construct both a network of quadrilateral patches and a parametrization of the data points over these patches. Ecks method is thoroughly based on topology constraints, without human intervention at all. It has both

advantages and disadvantages because users have no freedom to control the result (geometric behavior) of segmentation. It demands an eclectic scheme, which can automatically provide a rough segmentation, and allows for further human interaction to meet the user's desire to the maximum extent. In this paper, we present a semiautomated solution to segment the arbitrary triangle meshes into four-sided patches, which makes a reasonable tradeoff between automation and human control. Then we fit a B-spline surface to each patch, taking continuity condition into consideration. 2

SEGMENTATION OF TRIANGLE MESHES

2.1 Constraining corner points In the beginning of our segmentation approach, we allow users to interactively specify some points on the triangle mesh M as the corners of the B-spline patches. This step is optional. Users can constrain the obvious common points of several natural surfaces (if exist) as corner points, or randomly distribute some corner points, or specify no corner points at all. All the other necessary corner points will be automatically computed in the next step (see section 2.2). The purpose of this step is to augment the freedom of human control and to accelerate the following disk topology tiles partitioning procedure. If users label the corner points, it is recommended that the corner points should be distribute evenly, and there need not be many. 2.2

Disk topology tiles partitioning

A disk topology tile is a set of contiguous facets of M surrounding a corner point as illustrated in Figurel(a), where the solid circles represent the corner points and the bold lines are boundaries of the disk tiles. Denote the corner points as set C, which is initialized with the user specified points (more accurately speaking, the centroid points of the triangle facets that the user specified points are located in) in the last step. If users did not define any corner points, our algorithm initializes C with the centroid of a single randomly chosen facet in M. The key idea of this step is as follows: given a set C(c,,c2, ... ,cn) of corner points, for each facet fk in M, if

where d(c,, fk) denotes the geodesic distance between c, and the centroid point of fk, facet fk belongs to the tile corresponding to c,. However, The calculation of exact geodesic distance in dense mesh surface involves high space and time costs, and it is unnecessary in our research context, so we take an approximate solution in which the topological distance is calculated instead. Take the triangle mesh M as a graph G, which consists of vertices and edges of the mesh, with edge length given by Euclidean distance. We model the approximate geodesic distance calculation as the shortest path problem in the dual graph of G, denoted Gd. The nodes of Gd represent faces in M, and the edges of the graph connect nodes corresponding to adjacent faces. The cost of edges in Gd is set to be the distance between centroids of the corresponding faces. We calculate shortest paths in the graph using the Dijkstra's algorithm (most textbooks on graph theory, e.g. [9], have detailed descriptions on the algorithm). Taking the topological distance to a corner point as the priority of a facet, we establish a priority queue for each corner point. Using the priority queues, the algorithm dynamically grows the tiles from the corner points simultaneously, i.e. adds the nearest facet in the prior queues to the corresponding tile every time. The growth procedure continues until the tiles cover M or at least one of the tiles is no longer in disk topology. The criteria of disk topology is as follows: if the three vertices of the most recently added facet are already in the tile, that tile is no longer a disk. The algorithm will add that facet to the corner point set C and resume the tile growth procedure. If the newly added facet is on the boundary of M (If M is open) and there were no boundary facets in the tile before, the algorithm will also add the boundary facet to C and restart the growth procedure. When the growth procedure is completed, the algorithm checks if the following two situations exist in the tiles: (1) two adjacent tiles share more than one boundary, or (2) more than four tiles meet at one point. For situation ( I ) , add one corner point at the midpoint of one of the shared boundaries to set C, and restart the tile growth procedure. For situation (2), add the corresponding point to C, and resume the tile growth procedure

Figurel. Illustration of the corner points, disk topology tiles and the rough segmentation. We implemented the above method by referring to the Voronoi diagram construction approach presented by Eck et a/[7],[10]. But they imposed three more convergence conditions in the above iteration procedure to meet the requirements of thoroughly automated quadrilateral patches construction. Because the above procedure uses a greedy search, it can not theoretically guarantee convergence. More convergence conditions are more likely to result in algorithm divergence. The convergence conditions we used, fewer than that in [7],[10], obviously improve the convergence rate. Experiments show that the above disk topology tiles partitioning algorithm converges in all our test cases.

2.3 Constructing the rough segmentation The disk topology tiles constructed in the last step can also be represented by a graph, denoted Gt. In this step, we calculate the dual graph of Gt to construct the rough segmentation of M as shown in Figurel(b). The dashed lines, i.e. the dual graph of the disk tiles, constitute the boundaries of the rough segmentation. Eck et a/ [7] used the harmonic mapping method to calculate the dual graph of Voronoi diagram. For dense triangle meshes, harmonic mapping is expensive and can not guarantee the embedding of 3D meshes into 2D convex polygons. Instead, our approach consists of three steps. First, move each corner point to the vertex nearest to it on M. For each boundary tile, the corner point must be moved to the nearest vertex on the boundary of M in order to make sure the rough segmentation completely covers M. Second, calculate the shortest path between the corner points of each pair of adjacent disk topology tiles. Here again, the shortest path is in the criterion of topological distance and can be calclated using Dijkstra's algorithm by taking M as a graph G. All the shortest paths constitute the boundaries set in the rough segmentation. Third, smooth the boundaries. Figure2 illustrates the method we used in the smoothing process. Let A and B be two corner points. The bold successive edges in Figure 2(a) demonstrate the shortest path between A and B. C is a vertex in the middle of the path. In Figure 2(b), the plane P is defined by the three points A, B, and C. The relatively smooth path (the bold lines) is obtained by intersecting P with the corresponding facets in M.

Figure 2 : Smoothing the shortest path The rough segmentation, which is the dual graph of the disk tiles, consists of triangle and quadrilateral patches, for we have constrained that no more than four tiles meet at one point. 2.4 Adjusting the segmentation The last two procedures, the disk topology tiles partitioning and the rough segmentation construction, actually run as a single process. The automatically generated rough segmentation serves as the basis for quadrilateral patches partitioning. To improve the segmentation further and get the final quadrilateral patches that are suitable for piecewise B-spline surface fitting, we have developed a set of semi-automated facilities. Merginghubdividing patches. As mentioned above, there exist triangle and quadrilateral patches in the rough segmentation. In our solution, users could select two adjacent triangle patches to merge into a quadrilateral one as shown in Figure 3(a). There are three adjacent triangle patches at most around a triangle patch Ti. Which one of the three triangle patches will be merged with Ti is up to the users. Since the number of patches in the rough segmentation is typically not large, the interactive merging process is not difficult for users. Compared to the approach of automated matching of adjacent triangle patches, our

interactive method is more intuitive and reliable, and users have much freedom to control the final segmentation result. If some quadrilateral patches are too large or there are triangle patches that cannot be merged, we provide users a splitting tool as shown in Figure 3(b). The key ingredient of this tool is also calculating the shortest paths and then smooth them, as stated in the last section. We use the winged-edge structure to represent the topological structure of the patches. When users merge or subdivide the patches, the algorithm automatically updates the winged-edge structure to maintain the adjacency relations of the patches.

deletes the unwanted segment and constructs a new segment by connecting the specified points. If no points are specified except the segment end points, the algorithm automatically connects the end points. In some situations, especially when a corner point or a patch boudary curve is dragged far from the original lacation, the plane being moved could probably intersect the mesh more than once. When this occurs, the algorithm calculate first intersection segment and highlight it to warn the users to select a new location. If user ignore the warning, the algorithm will automatically set the boundary curve to the original location. 3 PIECEWISE BSPLINE SURFACES FITTING The equation of a tensor product B-spline surface can be written as

( b) Figure 3: Merging and subdividing patches.

where Nl,k(u),i=O,l,...m ; N,Jv), j=O,l, _ _ n _ , are B-spline basis functions of order k+l and 1+1 respectively, and d,,, denotes the control points (Textbooks on geometric modeling, e.g. [ I I ] , have detailed descriptions of B-spline definition). Consider the problem of fitting a B-spline surface to a given set of points P={p,, ...,p,} on a quadrilateral patch. If we can specify the corresponding parameters u, and v, of each point p,, we can substitute p I ( ~ l , ~into I ) equation (1) to generate a system of linear equations r n n

Modifying the corner points. This facility is for dragging a corner point on M to a new position as illustrated in Figure 4. Modifying a corner point will change all the boundaries adjacent to it. Again, we calculate the intersection between a plane and the facets in M to determine the corresponding boundary curves after modification. For a corner point with even degree as shown in Figure 4(a), the new corner position A and each pair of opposite corner points, e.g. B and C, determine the plane P that is used to calculate the new boundary. For a corner point with odd degree as shown in Figure 4(b), project the middle point of the straight line connecting the new corner position A and the adjacent corner point, e.g. B, to M. The projected point, A , and B define the plane to be used for the calculation of boundary A B .

Figure 4: Modifying the corner points Modifying the boundary curves. We modify the entire path of a patch boundary by specifying a point that the new boundary path will go through. We then use the plane defined by this point and the two end points of the boundary to intersect with M to get the new boundary. If part of a boundary is not satisfying, we have also developed a facility to locally modify the boundary curves. First, specify two end points of the unwanted segment and the points on M that the new segment will pass through. The algorithm then

where the d,,, are unknowns. In the context of this research, the number of p, is greater than that of d,,,, so the system is over-constrained. We can get the control points d,,, by solving the linear system using the least squares method. Therefore, the first thing is to associate a parameter pair (u,,v,)to each fitting point p,. 3.1 Parametrization of data points Consider an individual quadrilateral patch Q. We take the entire vertex of the triangles in Q as the fitting point set P. If P is not sufficiently dense, the centroid of each triangle and/or the midpoints of each edge in Q can be added to the set P. Obviously, the fitting points in P are unorganized. Without taking the gridded resampling approaches [5], [7], we directly calculate the parametrization of the unorganized points. First, construct a biconic Bezier surface S using the four boundaries of Q. Then, project each point p, in P onto S and the parameter (u,,v,) of the projected point on S is taken to be the initial parametrization of p,. Ref. [I21 had a more detail discussion on the parametrization of unorganized points. In the following surface fitting procedure, an iterative method is used to correct the initial parametrization. 3.2 Surface fitting We use end-point interpolating, uniform, cubic B-spline, which is widely used in the modeling industry. Every patch has the same number of control points, which is specified by the users. The knot vectors in the u and v directions of each patch are also set to be the same. The surface fitting method consists of the following iterative steps: Stepl: According to the parametrization of point set P, solve linear equations (2) using the least squares method for every patches. In this procedure, we enforce shared end points for each pair of adjacent patches to construct

piecewise position-continuous (i.e. Co continuity) B-spline s u rfaces. Step2: For each patch, project its fitting point set P on the fitted surface to get a more accurate parametrization of point set P and return to stepl. In most cases, the iteration converges very quickly. We initially impose Co continuity in the least squares procedure rather than lower or higher continuity for the following reasons. First, Co continuity is obligatory. Generally speaking, no one would want gaps to exist between ad'acent patches. Second, for some creases or corners, Cd IS adequate. If we enforce G1 or higher continuity in the above least squares process, we will lose these characteristic lines or points. Furthermore, G1 or even C2 continuity is inadequate in some applications, e.g. in constructing a car body surface. So, we prefer post-adjusting approaches (e.g. [6]) to meet higher continuity requirements wherever necessary. Compared to the surface spline approach [7], the fitting accuracy in our surface fitting method can be controlled by interactively changing spline resolution (the number of control points in each patch) rather than by subdividing the quadrilateral patches numerous times. 4 EXAMPLES Figure 5 shows an example. We first obtained the points on the cavum of a diesel engine by a laser range scanner and then reconstructed the triangle mesh shown in Figure 5(a). Figure 5(b) shows the diagram of the disk topoplogy tiles. We do not show this diagram to users when the program runs. It is for illustration purpose only. Figure 5(c) is the automatically generated skeleton segmentation. Figure 5(d) illustrates the quadrilateral patches and the iso-curves of the piecewise B-spline surfaces. Each step runs at almost real time on our PentiumIII733, 256MB memory computers.

(d)

Figure 5:B-spline surfaces fitting to a triangle mesh. 5 SUMMARY AND FUTURE WORK We have presented a solution for fitting piecewise Bspline surfaces to triangle meshes. The characteristics of the method are: (1) the meshes can be arbitrary, including closed or open with boundaries. (2) The segmentation process makes a reasonable tradeoff between automation and human control. In the future, we would like to develop a semi-automated characteristic lines recognition procedure to obtain more optimal segmentation.

6 AC KNOWLEDGM ENTS This work was partly sponsored by National Science Foundation of China grant 59905013, China Hi-Tech Project Foundation grant 863-51 1-942-022, Science Foundation of Jiangsu, China, grant BK2001408 and China Aeronautic Foundation grant 00H52069. 7

REFERENCES Fischer, A,, et a/, 1999, Utilizing image processing techniques for 3D reconstruction of laser-scanned data, Annals of the CIRP, 48/1: 99-102. Fischer, A,, et a/, 1997, Modeling in Reverse Engineering for Injection Molding Analysis of 3D Thin Objects, Annals of the CIRP, 46/1: 97-100 Guo, B.,1997, Surface reconstruction: from points to splines, Computer-Aided Design, 29(4): 269-277. Bajaj, C.L. , et a/, 1995, Automatic reconstruction of surfaces and scalar fields from 3D scans. SIGGRAPH'95 Proceedings, 109-1 17. Krishnamurthy, V., Levoy, M., 1996, Fitting smooth surfaces to dense polygon meshes, SIGGRAPH'96 Proceedings, 313-324. Milroy, M. J., et a/, 1995, G1 continuity of B-spline surface patches in reverse engineering, ComputerAided Design, 27(6): 471-478. Eck, M., Hoppe, H., 1996, Automatic reconstruction of B-spline surfaces of arbitrary topological type. SIGGRAPH'96 Proceedings, 325-333. Edelsbrunner, H. , et a/, 1994, Three-dimensional alpha shapes, ACM trans. on Graph., 13(1): 43-72. Aho, A. , et a/, 1983, Data structure and algorithms, Addison-Wesley, Reading, Mass. Eck, M., et a/, 1995, Multiresolution analysis of arbitrary meshes, SIGGRAPH'95 Proceedings: 173182. Farin, G., 1992, Curves and Surfaces for Computer Aided Design, 3rded. ,Academic Press. Ma, W., Kruth, J.P., 1995, Parameterization of randomly measured points for least squares fitting of B-spline curves and surfaces, Computer-Aided Design, 27(9): 663-675.