Computing lines of curvature for implicit surfaces - Semantic Scholar

Computer Aided Geometric Design 26 (2009) 923–940

Contents lists available at ScienceDirect

Computer Aided Geometric Design www.elsevier.com/locate/cagd

Computing lines of curvature for implicit surfaces XiaoPeng Zhang a,b , WuJun Che a,b,∗ , Jean-Claud Paul a,c,d a

LIAMA, CAS Institute of Automation, Beijing 100190, PR China NLPR, CAS Institute of Automation, Beijing 100190, PR China c School of Software, Tsinghua University, Beijing 100084, PR China d INRIA, Le Chesnay Cedex 78153, France b

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 13 August 2008 Received in revised form 2 July 2009 Accepted 13 July 2009 Available online 16 July 2009 Keywords: Implicit surface Line of curvature Principal foliation Principal configuration Umbilical point Visualization

Lines of curvature are important intrinsic characteristics of a curved surface used in a wide variety of geometric analysis and processing. Although their differential attributes have been examined in detail, their global geometric distribution and topological pattern are very difficult to compute over the whole surface because of umbilical points and unstable numerical computation. No studies have yet been carried out on this problem, especially for an implicit surface. In this paper, we present a scheme for computing and visualizing the lines of curvature defined on an implicit surface. A key structure is introduced, conveying significant structure information about lines of curvature to facilitate their investigation, rather than computing their whole net. Our current framework is confined to a collection of manageable structures, consisting of algorithms to locate some seed umbilical points, to compute the lines of curvature through them, and finally to assemble this structure. The numerical implementations are provided in detail and a novel evaluation function measuring the violation of umbilical points in an implicit surface, i.e. indicating how much a point is to be umbilical, is also presented. This paper is the continuation of [Che, W.J., Paul, J.-C., Zhang, X.P., 2007. Lines of curvature and umbilical points for implicit surfaces. Computer Aided Geometric Design 24 (7), 395–409]. © 2009 Elsevier B.V. All rights reserved.

1. Introduction A line of curvature is one of the most important geometric attributes for a curved surface, indicating a direction field for the extreme curvature across the surface. Together with the principal curvatures, lines of curvature (LOC) provide a complete description of variations in curvature of the surface. They are useful tools in surface analysis for exhibiting variations of the principal directions. LOC can guide the analysis of surfaces, widely used in geometric design, shape recognition, polygonization of surfaces, surface rendering etc. (Alliez et al., 2003; Beck et al., 1986; Hallinan et al., 1999; Maekawa et al., 1996; Patrikalakis and Maekawa, 2002). 1.1. Background Given a point p on a smooth manifold surface S in 3D Euclidean space, the normal curvature in a tangent direction at p is the curvature of the intersection curve of S with the plane generated by the tangent vector and the normal to S at p. The curvature depends on the tangent vector. The principal curvatures are the extreme values of normal curvature –

*

Corresponding author. E-mail addresses: [email protected] (X.P. Zhang), [email protected] (W.J. Che), [email protected] (J.-C. Paul).

0167-8396/$ – see front matter doi:10.1016/j.cagd.2009.07.004

© 2009 Elsevier B.V.

All rights reserved.

924

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

maximum or minimum, and their corresponding directions are principal directions. If both principal curvatures are equal, p is an umbilical point and every tangent vector can be thought of as a principal direction. The integral curve of maximal (resp. minimal) principal direction field is called a line of maximal (resp. minimal) curvature. Two principal directions are uniquely determined and are orthogonal at p away from an umbilical point. Therefore, there are always two LOC through non-umbilical p and they cross at a right angle. On S without umbilical points, LOC constitute an orthogonal net, but the existence of umbilical points destroys its elegance – the net patterns are disturbed significantly and could become very complex, leading to a breakdown in the regularity of orthogonality. Therefore, an umbilical point is a singularity with anomalous behavior in the LOC net. The set of umbilical points is denoted by U ; the family of maximal (resp. minimal) LOC by Fi , i = 1, 2, and called the maximal (resp. minimal) principal foliation of S. The triple P = (U , F1 , F2 ) is the principal configuration of S. Historically, principal configuration has been studied from the angle of differential geometry. Its differential geometry is usually based on the parametric representation of surfaces. The study of principal configuration locally around an umbilical point has received considerable attention. See Gutiérrez and Sotomayor (1990), Porteous (2001), Gutiérrez et al. (2004), Sotomayor and Garcia (2007) and references therein. The global structure of principal configuration, however, has been known only for very rare classical examples (Gutiérrez and Sotomayor, 1990). It is difficult to calculate the LOC net to visualize the principal configuration correctly in contrast to principal curvatures and those work provide no methods to obtain it globally. Although there has been some related studies, they mainly concentrate on local computation and explain little about computing the global behavior of LOC over the surface, especially on an implicit surface. To the best of our knowledge, no algorithm is yet able to describe any global differential pattern of LOC with some guarantee. Even their computation with U = ∅ is still impossible. LOC are difficult to obtain because the reconstructed line, with simple numerical integration techniques, diverges from the true one, and the algorithm never ends at the expected point due to computational errors. Relaxation algorithms have been proposed for a closed line, but they can work only if there is some information on how to place an initial closed curve in the surface (Thirion, 1996). An alternative is to generate the field of principal directions via visualization technologies of vector field over its discretization by sampling the surface, but it usually gives us an ambiguous observation for surfaces with a high variation in curvature unless the sampling is dense enough. In order to represent the vector field by its streamlines, i.e. integral curves, we have to choose a set of starting points on the surface in such a way that the density of the points is proportional to the intensity of the field, and then trace each starting point in both positive and negative directions to generate the streamlines, going back to the way by using numerical integration methods. Generally speaking, two aspects need to be taken into account to give an accurate description of P . The first involves the geometric distribution of Fi , i = 1, 2, and the second is the topological patterns among F1 , F2 and U . Their uneven distribution throughout the surface results in the difficulties of robust computation, in particular near an umbilical point. Few studies have dealt with LOC defined on an implicit surface. The primary goal of our paper is to attempt to develop numerical computation techniques for the automatic generation of LOC on an implicit surface. 1.2. Previous research Umbilical points are attributes closely related to LOC. The earliest classification of umbilical points and the generic features of LOC near an umbilical point are due to Darboux (1896) (Maekawa et al., 1996; Porteous, 2001; Sotomayor and Garcia, 2007). Since then, many authors have carried out much research on them, most of which was targeted at local differential properties (Gutiérrez and Sotomayor, 1990; Gutiérrez et al., 2004; Maekawa et al., 1996; Porteous, 2001). The study of the global features of P has also received considerable attention, and see Sotomayor and Garcia (2007) for references. No studies, however, have considered aspects of numerical computation, which are important to understand a given surface in practice since a priori computation is usually infeasible. For a parametric surface, especially for a general free-form surface, some previous work does exist which has investigated the numerical computation of LOC. These methods usually consider two coupled, non-linear differential equations in terms of two parametric coordinates respectively, resulting in a numerical method for a B-spline surface (Beck et al., 1986; Farouki, 1998; Maekawa et al., 1996). Maekawa and Patrikalakis (1994) described a computational method to locate isolated umbilical points on a parametric polynomial surface. But how to derive the LOC net with umbilical points was not fully discussed in these studies. Little work has been carried out on other surface representations. Alliez et al. (2003) reduced LOC calculation on a triangular mesh to a 2D parametric plane by flattening the surface via a discrete conformal parameterization. The umbilical points were found by going over each triangle and solving a 2 × 2 linear system in the tensor field. Their method is still based on the parametric surface form and as some LOC information could be lost, it fails to fit accurately to the computation. Che et al. (2007) studied the explicit representation of principal directions defined on an implicit surface, based upon which the differential equation of LOC and a new criterion for umbilical points were obtained, but their numerical implementation was not carried out. To our knowledge, how to locate umbilical points on an implicit surface has not yet been reported in literature. It is impossible to cover the surprising richness of differential geometry on LOC and umbilical points. The interested readers can refer to related work and textbooks (Beck et al., 1986; Che et al., 2007; Maekawa et al., 1996; Porteous, 2001), and also references therein for basic facts. The study of LOC and their umbilical singularities is described in Gutiérrez and

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

925

Sotomayor (1983, 1990), Gutiérrez et al. (2004), Porteous (2001), Sotomayor and Garcia (2007). Some relevant results are also listed in Section 2. An active area related to the work here presented is the visualization of Vector Field Topology (Delmarcelle and Hesselink, 1994; Helman and Hesselink, 1989, 1991; Scheuermann and Tricoche, 2005). As an ever increasing amount of scientific data sets of vector and tensor challenge scientists and engineers to look for a more manageable representation that may be visualized interactively, Helman and Hesselink (1989, 1991) suggested that a qualitative structure of the global topology based on critical points and streamlines through them can be employed. These algorithms start with a piecewise linear or bilinear interpolation of the data at certain grid positions. Then, one looks for all the saddle points and starts the numerical computation of the streamlines from there (Delmarcelle and Hesselink, 1994). 1.3. Contributions Che et al. (2007) carried out a ground study on the mathematics of LOC on an implicit surface and in this paper we study their numerical implementation and visualization, which do not appear to be explicitly described elsewhere. We address two general problems: identifying umbilical points and generating LOC on an implicit surface as accurately and robustly as possible. But two main challenges are encountered:

• Evaluating differential quantities associated with the implicit surface. • Improving accuracy through the reduction of numerical errors. In this paper, we propose an approach to address the issues above to some extent. We concentrate on numerical computation of a topological structure of P . Our method is structured as follows: 1. A novel criterion is proposed to measure the violation of umbilical points on an implicit surface. Our study shows that it is more efficient than other methods. 2. Novel algorithms are proposed to calculate umbilical points on an implicit surface. 3. A method is proposed for computing a line of curvature. In this method, the integration is performed on the 2-manifold of the implicit surface instead of the parametric domain. 4. A framework, centering on umbilical points, is developed to compute the global structure of principal configuration under some constrained conditions. An important characteristic of our method is that it does not rely on the parametrization of the surface, which makes our approach general enough for applications, and some techniques in this paper can also be applied to other representations, e.g. a parametric surface. In contrast to the scheme proposed by Delmarcelle and Hesselink (1994) for discrete linear tensor field, our framework characterized by several advantages:

• The requirement is eliminated to exhaust all singularities before integration, for it is impossible for a generic smooth manifold surface. Accordingly, a new procedure is designed to detect the rest during LOC integration.

• Bidirectional integration is taken to improve precision in calculation such that an integral curve unnecessarily terminates at a singular point with parabolic sectors. The remainder of this paper is organized as follows. In Section 2, we review essential preliminaries upon which our framework is based and reasonable conditions are imposed upon implicit surfaces handled in this paper, leading to a definition of directed P graph to facilitate LOC computation. An overview of its computation framework are given in Section 3. Greater details are presented in following sections respectively: a new measure function for umbilical points and how to find these points based on it in Section 4; computing Monge form in Section 5 and some useful techniques for integration of LOC in Section 6. Examples are given to illustrate our method and to demonstrate its application in Section 7. Finally we conclude this paper with a discussion in Section 8. 2. Preliminaries 2.1. LOC on an implicit surface The analysis of LOC defined on an implicit surface has already been developed (Che et al., 2007). In this section, we give only a brief review. Let x = (x, y , z) T and dx = (dx, dy , dz) T . Given an implicit surface H (x) = 0, we assume it is regular, i.e. ∇ H = 0. The differential equation of LOC defined on H (x) = 0 is given by



I (dx) = dx T ∗ M ∗ dx = 0 II(dx) = ∇ H · dx = 0

(1)

926

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

where M =

 C

A

B B D E C E F

and A = H y H zx − H z H yx , B = ( H z H xx − H x H zx + H y H zy − H z H y y )/2, C = ( H y H zz − H z H yz + H x H yx −

H y H xx )/2, D = H z H xy − H x H zy , E = ( H x H y y − H y H xy + H z H xz − H x H zz )/2, F = H x H yz − H y H xz . To solve Eqs. (1), let H z = 0 without losing generality, and substituting dz from II(dx) = 0 to I (dx) = 0, we have

III(dx, dy ) = U dx2 + 2V dx dy + W dy 2 = 0

(2)

V = where U = − 2C H x H z + − C H y H z − E H x H z + F H x H y , and W = From Eqs. (1) and (2), we can derive the principal directions, respectively: A H 2z

⎛

T1 = ⎝

dx X

F H x2 ,

(− V + (V −

√

B H 2z

√

V2

⎛

⎞

V 2 − U W )H z U Hz

T2 = ⎝

⎠,

− U W )H x − U H y

(− V − (V +

√

D H 2z

√

V2

V 2 − U W )H z U Hz

− 2E H y H z +

F H 2y ,

respectively.

⎞ ⎠

(3)

− U W )H x − U H y

We write Ti = ( X , Y , Z ) T , i = 1, 2. Thus the differential equation of a line of curvature on H (x) = 0 can be rewritten as = dy = dz , whose a special case is Y Z

dx ds

=

Ti

(4)

Ti 

where s is the parameter of arc length. 2.2. Umbilical points and Monge form We now list some useful results about umbilical points and the patterns of LOC in proximity to them using Monge form. 2.2.1. Umbilical points It is an essential step to find the umbilical points for accurate computation of P . Theorem 2.1. (See Che et al., 2007.) p is an umbilical point if and only if U = V = W = 0. Theorem 2.1 provides the foundation to trace umbilical points in our work. 2.2.2. Monge form To obtain P accurately, we also need to compute limit principal directions of LOC tending to an umbilical point p, which can be determined by the Monge form of the surface S. Suppose T p S is the tangent plane of S at p, then: Definition 2.2. (See Cazals and Pouget, 2005; Porteous, 2001, p. 200.) A limit principal direction of p is a direction in T p S which is tangent to a line of curvature ending at p. For simplicity, we denote such a limit principal direction of p by Dp . Consider a surface of the form as S = [u , v , w (u , v )] T . If the origin of p is umbilical and T p S coincides with the u − v plane, then we can Taylor expand the w-component of S to the Monge form as:

w (u , v ) =

where

κn (0, 0) 2



u2 + v 2 +



+ o u3 + v 3

1 6

w uuu (0, 0)u 3 + 3w uuv (0, 0)u 2 v + 3w uv v (0, 0)uv 2 + w v v v (0, 0) v 3

(5)

κn (0, 0) is the identical normal curvature at the origin. Eq. (5) contains a cubic form w (u , v ) = w uuu (0, 0)u 3 + 3w uuv (0, 0)u 2 v + 3w uv v (0, 0)uv 2 + w v v v (0, 0) v 3

(6)

An umbilical point of U is a singularity of F1 and F2 . Generic umbilical points can be categorized into three types: a star, a monstar and a lemon as in Fig. 1(a)–(c) respectively, also referred to as ordinary umbilical points in Porteous (2001). Their classifying form is just Eq. (6). Other umbilical points exist for degenerative but non-vanishing Eq. (6) (Gutiérrez et al., 2004; Porteous, 2001). If Eq. (6) vanishes, higher order umbilical points need to be considered, but we find that it is discussed quite rarely in literature and are not computationally robust. A special pattern is one through which an infinite number of LOC pass in all directions (Fig. 1(d)), called a perfect star in this paper. Although it is non-generic, a perfect star is often present in geometric surfaces of modeling, especially symmetrical ones. To compute all Dp s, we can further express the Jacobian cubic form of Eq. (6) as









w uuv (0, 0)u 3 − w uuu (0, 0) − 2w uv v (0, 0) u 2 v + w v v v (0, 0) − 2w uuv (0, 0) uv 2 − w uv v (0, 0) v 3 whose root lines are precisely Dp s (Maekawa et al., 1996; Hallinan et al., 1999).

(7)

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

(a) A lemon

(b) A monstar

(c) A star

927

(d) A perfect star

Fig. 1. Patterns of umbilical points: the solid lines correspond to LOC of F1 and the dashed of F2 .

2.3. LOC and umbilical points A compact line of Fi , i = 1, 2, is called a principal cycle. A line of curvature c is either closed, i.e. a cycle, or open. If it is open, then c can be oriented and parameterized by the arc length on its maximal interval I = (ω− ; ω+ ). Its α (c) (resp. ω(c)) limit set is the collection of limit points of sequences c(sn ), convergent in S, with sn tending to ω− (resp. ω+ ). The limit set of c is the set α (c) ∪ ω(c), which may be umbilical points, cycles or other types. LOC are incalculable, but we observe that special lines approaching to an umbilical point divide Fi , i = 1, 2, in its vicinity into groups of different patterns. Definition 2.3. A separatrix1 is a line of curvature approaching an umbilical point p in direction Dp , where p is in its limit set. If it is not unique, one that contacts best with the line collinear with Dp in the vicinity of p is picked up. By Definition 2.3, a lemon has two separatrices, maximal and minimal respectively, while a monstar or a star has six. A perfect star is special, in each D of which a separatrix exists and all of them are either maximal or minimal; its status in U is like an umbilical point in F1 and F2 . Remark 2.4. Special care must be exercised for a monstar through which an infinite number of LOC pass. For a monstar p, three Dp s are contained within a right angle and all LOC in this angle end with the same Dp at it (see Fig. 1(b)). In this case, not all of the LOC relevant to the middle Dp serve as a separatrix but only that contacting best with the in-between root line of Eq. (6). It is clear now that we take Definition 2.3 so as to associate each D with a separatrix tending to this umbilical point. A separatrix of two distinct umbilical points is called an umbilical connection; or an umbilical loop2 if it is of the same umbilical point. A line of curvature can actually reach at infinity. On this occasion, we regard the point at infinity as a special umbilical point. The terminology pertinent to it, e.g. a connection, can be generalized accordingly. Hence, intuitively speaking, each line of curvature usually either starts from an umbilical point, a connection, a loop or a cycle and ends at another one, or has a closed orbit. The principal configuration of a surface has too much diversity to be covered in this paper. The interested reader can refer to Gutiérrez and Sotomayor (1983, 1990), Gutiérrez et al. (2004), Porteous (2001), Sotomayor and Garcia (2007) and references therein for more details. 2.4. P graph Before proceeding further, we need to restrict the class of implicit surfaces, but still covering practically important cases, to introduce P graph: 1. H (x) is continuously differentiable up to fourth order in the vicinity of H (x) = 0, mainly settling for differential analysis of umbilical points (see Sections 4 and 5). 2. The number of umbilical points is finite and they are isolated. 3. The surface is regular, and almost generic, i.e. most umbilical points are ordinary. As in most of the literature for numerical processing of LOC, we mainly deal with ordinary umbilical points; one special type, say a perfect star, is also manipulated, as listed in Fig. 1. 4. The limit set of every line of curvature is the union of umbilical points, connections, loops or cycles.

1

This definition is a little different from that in theory of principal foliation (Gutiérrez and Sotomayor, 1983). See Remark 2.4 too. In the theory of principal foliation such as in Gutiérrez and Sotomayor (1983, 1990), a loop is taken for a special case of a connection whose umbilical points coincide. Here we partition off loops from connections for the sake of clarity. 2

928

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

Fig. 2. Point symbols: Red stands for a lemon, green for a monstar, blue for a star, black for a perfect star, magenta for an initial guess and cyan for a seed point. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

These conditions can be regarded from viewpoint of numerical computation as a generalized variant of those for structurally stable principal configuration (Gutiérrez and Sotomayor, 1983, 1990), which implies that stable principal foliations are described with the set of umbilical points and cycles, and their connection fashions by the separatrices of F1 and F2 . The entire patterns of LOC are dominated by this special structure, capable of offering a clear perception of the structure of P . A separatrix is always related to an umbilical point and a relevant D , so it is naturally oriented. Therefore, we define a directed geometric graph on the surface: Definition 2.5. A graph of principal configuration, or P graph, is a directed graph embedded on a surface satisfying the following conditions: 1. An umbilical point is a simplex vertex. 2. A directed edge is a separatrix at an ordinary umbilical point p oriented in Dp ; p is the initial vertex and the limit set of the other side is a terminal, i.e. umbilical points, connections, loops or cycles. If the terminal vertex is not simplex, it is multiplex. 3. A connection of two ordinary umbilical points or a loop of an ordinary umbilical point is a bidirectional edge. A P graph separates S into simple areas, in each of which principal configuration can be derived from P graph, and thus permits fast extraction of global LOC structures that are directly related to region of interest. By definition, an initial vertex must be an ordinary umbilical point; for a bidirectional edge its vertices are both initial and terminal. Our method takes advantage of bidirection for LOC integration. The primary objective of our work concentrates on the production of the P graph. However, calculating a multiplex vertex is a limit process; no methods are reported to do so on a smooth surface with certain guarantee. Due to the complexity of the issue, developing a framework to extract P graph with multiplex vertices cannot be treated comprehensively in this work yet. We will give some facts of difficulties about it in a nutshell at the last two sections. But here we point out that only basic integration issues of P graph with simplex vertices are studied in this paper, not treating with a multiplex vertex for reasons of technical difficulties. 3. Overview In this section, we will discuss the computational procedures of our framework. For the sake of clarity, we present only a general overview and elaborate on some central computational algorithms in following sections, respectively. 3.1. P graph in implementation The vertex of P graph associated the umbilical point at infinity, however, can never be reached in computation. The treatment is: Its associated line of curvature enters and exits from the domain bounds and results in an intersection point to them; thus this vertex can be replaced in P graph by the intersection point in implementation to keep the graph valid. A P graph can be empty if the surface contains no umbilical points at all. The concerned points relevant to umbilical points in this article are represented in different colors, shown in Fig. 2, where initial guess and seed point will be explained in the following subsections. 3.2. Method outline To accurately compute the P graph of LOC, several aspects should be considered carefully. Our framework centers on umbilical points, and is concerned with the pattern that LOC make in the vicinity of an umbilical point. Fig. 3 presents the procedural flow using an example illustrative of the scheme. The framework is divided into three parts: 1. Location of umbilical points. We reduce it to find the roots of a non-linear equation. Unfortunately, it is still an impossible task to compute all the roots of a generic non-linear function. In our application, the basic idea is to estimate sufficiently good initial guesses of umbilical points and then to converge them to solutions. We propose two approaches to carry this out. Firstly, we polygonize the implicit surface into a sparse mesh and assign to each vertex a differential quantity measuring the violation of an umbilical point; then we go through each vertex to fix local minimal ones serving as initial guesses;

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

(a) Initial guesses

(b) Seed points

929

(c) P graph

Fig. 3. Computation of P graph. In the first step of (a), the input implicit surface is polygonized and those vertices of the resulting sparse mesh, say 3.4K vertices, with local minima are picked as initial guesses. A procedure is then employed to optimize them, resulting in initial umbilical points as seed points, where initial guesses are also displayed in contrast against seed points in (b). The P graph are finally traced based on seed points in (c), where the blue lines correspond to the maximal LOC and the red ones minimal.

1. 2. 3. 4.

Compute initial guesses of umbilical points based on a polygonizer of implicit surfaces. Compute umbilical points based on the guesses above. Calculate D s of ordinary ones and push them into a stack Q as seed points. While ( Q is not empty) Do p=Pop( Q ); For (All unmarked Dp s of p) Mark Dp and integrate from p in direction Dp . If (A new ordinary umbilical point is found during integration) Calculate its D s and push it into Q . If (it is terminal) Mark its D associated with this line of curvature. EndIf EndIf EndFor EndWhile Fig. 4. The basic framework of computation of P graph.

finally we apply an optimization procedure to them and store those points of convergence to the global minimum of zero. This method, of course, cannot detect all umbilical points consequentially, which is the task to be carried out in the second approach, but the resulting points can work as seed points to trace the P graph. In the second approach, the search incorporates the integration of LOC described in the following step (3). During the integration, we make tracks simultaneously for the current solution point of this line of curvature to verify if it is falling into the neighborhood of an umbilical point; if so, a guess is found. 2. Calculation of D s for each ordinary umbilical point separately. For this purpose, a rigid coordinate transformation is applied to set the surface in Monge form. 3. Integration of LOC. A line of curvature is an integral curve of the field of principal directions. Based on Eq. (4), the integration starts from an ordinary umbilical point in one D . When it is propagated step by step along this line of curvature, it is checked whether this point is stepping into a trap neighborhood of another umbilical point, until it reaches a terminal simplex vertex or until it spans the domain bounds. We then integrate all these steps into a basic framework for computing a P graph of LOC (Fig. 4). 4. Finding umbilical points In this section, we explain our approaches to calculating umbilical points. Finding umbilical points on a surface enables us to solve a non-linear system of equations. When H (x) is polynomial, it is a system of polynomial equations, and many efficient methods exist. If H (x) is non-polynomial, as we will explain, this system can be further transformed into finding a global extremum of a function. In general, finding a global extremum is, however, a comparably difficult problem. Standard heuristics is therefore used, i.e. to find local extrema originating from widely varying starting points of the surface. We then pick the most extreme of points and employ the Newton–Raphson method.

930

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

4.1. Non-linear systems of equations of umbilical points According to Theorem 2.1, we need to solve a non-linear system of equations to fix umbilical points:

⎧ U =0 ⎪ ⎨ V =0 ⎪ W ⎩ =0 H =0

(8)

To solve Eqs. (8), we carry out multidimensional root finding by collapsing them into one. Therefore, let f (x) = U 2 + V + W 2 . Then p is an umbilical point if and only if f (p) = 0 subject to H (p) = 0. Since f  0, an umbilical point occurs when the function f (x) has a global minimum of zero. See Section 4.4 for more details. In CAD community, however, polynomial surfaces are popular. On this occasion, Eqs. (8) are then a non-linear polynomial system, for which many global solution methods have been designed to compute all roots. Patrikalakis and Maekawa (2002) distinguished them as three classes: 2

1. Algebraic and hybrid techniques 2. Homotopy (continuation) methods 3. Subdivision methods Class 1 is usually used in systems of computer algebras while Class 3 is the most frequently employed in CAD practice. Once all roots can be found in advance by those methods, the following steps in this section can be omitted to save computation cost as a consequence. Being restricted to the purpose of this paper, we do not spread out this topic further and the interested readers can refer to Patrikalakis and Maekawa (2002) and references therein. 4.2. Criteria for umbilical points In Newton’s method, the Jacobian matrix of the non-linear equations should be computable. However, an evaluation function measuring the violation of umbilical points is also required to confirm whether the current point is close enough to an umbilical point and to stop further iterations if so; the evaluation function should be not only computationally robust but insensitive to the choice of non-zero component of H (x). Therefore, we design two measure functions for iteration of Newton’s method and the violation of umbilical points respectively. 4.2.1. A new evaluation function Although f (x) can serve as a measure function of umbilical points, it depends on the choice of non-zero components of ∇ H . To keep the criterion the same for all points, we need another evaluation function. It has been shown in Che et al. (2007) that I (dx) is an intersecting plane-pair at an umbilical point if M = 0, one of which coincides with the tangent plane at this point, i.e.

I (dx) = II(dx)(y · dx)

(9)

where y is an undetermined vector satisfying Eq. (9). It is obvious that Eq. (9) holds even if M = 0. Therefore Lemma 4.1. p is an umbilical point if and only if Eq. (9) holds. Based on Lemma 4.1, we develop another novel criterion to quantify an umbilical point as follows. Eq. (9) fails to hold true on a non-umbilical point, but an optimal solution in a least-square sense can be defined. Comparing the coefficients on both sides of Eq. (9), we obtain:

A∗y=d

(10)

where

 A=

Hx 0 0

Hy Hx 0

Hz 0 Hx

0 Hy 0

0 Hz Hy

0 0 Hz

T ,

d=(A

B

C

D

E

F )T

The least-square solution of Eq. (10) is y = (A T ∗ A)−1 ∗ A T ∗ d. Lemma 4.1 shows that |A ∗ y − d| = 0 if and only if it is defined on an umbilical point. Let g (x) = |A ∗ y − d|2 = |(A ∗ (A T ∗ A)−1 ∗ A T − I) ∗ d|2 . In this paper, we regard g (x) as the criterion measuring an umbilical point, which is independent of using a non-vanishing component of ∇ H against f (x). On a regular implicit surface, we have

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

931

Fig. 5. A tracing method: pu1 is an initial vertex and pu2 terminal; any gray point can work as an initial guess, e.g. p.

 2  H + H2 + H2  T   x 2 y2 z A ∗ A =  Hx H y   H2 H2 x

H x2 H 2y H x2 + H 2y + H 2z

z

H 2y H 2z

   2 2 >0 H y Hz  2 2 2 H +H +H H x2 H 2z

x

y

z

Thus g (x) is always well defined. Remark 4.2. Both f (x) and g (x) can be used to measure an umbilical point, but they are used in different manners in this paper. We note that the gradient of f (x) is computable but dependent on the non-zero component of ∇ H , while the opposite is true for g (x), i.e. ∇ g is difficult to compute but with good global behavior of measuring an umbilical point. So we make f (x) serve as an evaluation function of Newton’s method while g (x) as the real measure function for a point x 2 tending to be an umbilical point. Although the third measurement h(x) = κ M − κG can also be applied, where κ M is the mean curvature and κG the Gaussian curvature, we will compare their behavior and point out that g usually works better as a measure function via experiments in Section 7. 4.3. Local neighborhood detection We design two approaches for identifying the neighborhood of an umbilical point: a brute-force method and a tracing method. 4.3.1. A brute-force method As its name implies, a brute-force method is used to go over the domain of the surface and fix some points on the surface for initial guesses of umbilical points. The g value is the global minimum of zero at an umbilical point. Thus it is also locally minimal. We polygonize the surface to be a triangular mesh. Then a vertex p of the resulting mesh can be selected as a candidate if g (p)  g (p ) for all adjacent vertices p to p. The proposed polygonizer in Bloomenthal (1994) is used in our work. p, as an initial estimate of an umbilical point, is gradually refined by iteration until it converges to a local minimum, which is then verified if it is a global one. If the mesh is fine enough, compared to the minimum of the distances between umbilical points, it is possible to search for a guess point in the local neighborhood of any umbilical point. But it strongly depends on the surface geometry and cannot be known a priori, so it is difficult to learn in advance whether a starting mesh is sufficiently accurate. Since we usually cannot predict how fine the mesh should be and a sparse mesh is favored here, some umbilical points could be lost. 4.3.2. A tracing method A tracing method is designed based on the fact that the P graph covers the neighborhoods of all umbilical points. Starting from an umbilical point in its one D , the tracing moves with the integration along this line of curvature, as demonstrated in Fig. 5. We keep a close watch on the value of g during the integration advance. If it begins to decrease, the present sample point is falling in the neighborhood of a local minimum. We think now it steps to the trap of an extremum. A trap interval consists of consecutive points along this line of curvature, with g values decreasing firstly and then increasing. Any point in the trap interval can actually work as an initial estimate. But for better convergence, we can select those just before g begins to increase in implementation. The tracing proceeds in a trap-interval-by-trap-interval fashion until the integration terminates. 4.4. Optimization of an umbilical point We now describe how to find the local minima of f (x) subject to the constraint of H (x) = 0. Here, Lagrange multipliers are used to reduce it to a solvable problem without constraints. The Lagrangian, Λ, is defined as

932

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

Λ(x, λ) = f (x) + λ H (x)

(11)

A critical value of Λ occurs when its gradient is zero, that is



∇Λ =

∇x f + λ∇x H



H

=0

(12)

Eq. (12) can be solved by Newton’s method. Before the iteration, we employ the steepest descent method to optimize an initial guess sufficiently close to the root. If a convergence solution is attained, a further validation is required if the g value is close enough to zero to be an umbilical point. Remark 4.3. f can also be defined as f (x) = U 2 + V 2 + W 2 + H 2 . The advantage is that it contains the constraint H (x) = 0 simultaneously and only standard Newton’s method is enough to solve it. We still prefer the former rather than this definition because the value of H (x) does not reflect an umbilical point directly so that it usually imposes an unexpected effect, e.g. local convergence, upon f if f = U 2 + V 2 + W 2 + H 2 is employed. We also find in our experiment that the former one is usually with better global convergence indeed. 5. Monge form Assuming that p is an umbilical point of H (x) = 0, we take p as the origin of the coordinate system, say p–uv w, to represent the surface in the Monge form at p with local coordinates μ = (u , v , w ) T for any point x. We choose unit vectors u, v, w as the axes, where u is a tangent vector in T p , w is the unit normal vector of the surface and v = w × u. Let B = (u, v, w). Then we have x = p + B ∗ μ. Now we represent the implicit surface in the frame p − uv w as

H (p + B ∗ μ) = 0

(13)

To derive the Monge form at p, we need to compute w uuu , w uuv , w uv v , w v v v at this point, which are given by using chain rule on both sides of Eq. (13). Thinking of w as a dependent variable and taking first derivatives with respect to u and v respectively, we obtain

∇H · u ∇H · w ∇H · v wv = − ∇H · w

wu = −

(14) (15)

Taking the derivative of Eq. (14) with respect to u, we have

w uu = −

d∇ H du

∇H · u + w u ddv ·w ∇H · w

(16)

∇H where ddu can be computed as follows:

d∇ H du

= ∇2 H ∗ u + wu∇2 H ∗ w

(17)

Noting that w u (0, 0) = 0, Eqs. (16) and (17) can rewritten at p as, respectively:

w uu (0, 0) = −

 · u  ∇ H · w p d∇ H du

(18)



  = ∇2 H ∗ u  p du p

d∇ H 

(19)

In much the same way as above, we can compute the other intermediate formula as follow:

w v v (0, 0) = −

 · v  , ∇ H · w p d∇ H dv



  = ∇2 H ∗ v  ,  p dv p

d∇ H 



  d2 ∇ H   = ∇ 3 H , v, u + w uv ∇ 2 H ∗ w  , p du dv  p

where [∇ 3 H , v, u] is defined as







   = ∇ 3 H , u, u + w uu ∇ 2 H ∗ w   p

d2 ∇ H  du 2

p



  d2 ∇ H   = ∇ 3 H , v, v + w v v ∇ 2 H ∗ w  p dv 2  p

∇ H x ∗v·u ∇ 2 H y ∗v·u . ∇ 2 H z ∗v·u 2

Finally, third derivatives of w with respect to u and v at p are calculated as following:

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

w uuu (0, 0) = − w uuv (0, 0) = −

w uv v (0, 0) = −

d2 ∇ H du 2

d2 ∇ H du dv

d2 ∇ H dv 2

d ∇H dv 2

p

·u+

∇H w uv ddu

∇H w uu ddv

 · w  

·u+

∇H w uv ddv

∇H w uv ddv

 · w  

∇H w v v ddv

 · w  

·w+ ∇H · w ·w+

∇H · w 2

w v v v (0, 0) = −

 ∇H ∇H · u + w uu ddu · w + w uu ddu · w   ∇H · w

·v+

∇H w v v ddv

·w+

∇H · w

933

p

p

p

Because p is an umbilical point, we also have w uv (0, 0) = 0 and w uu (0, 0) = w v v (0, 0) = κn (0, 0). Thus, they can be rewritten as

 [∇ 3 H , u, u] · u + 3κn ∇ 2 H ∗ u · w  w uuu (0, 0) = −  ∇H · w p  3 2  [∇ H , u, u] · v + κn ∇ H ∗ v · w  w uuv (0, 0) = −  ∇H · w p  [∇ 3 H , u, v] · v + κn ∇ 2 H ∗ u · w  w uv v (0, 0) = −  ∇H · w p  3 2 [∇ H , v, v] · v + 3κn ∇ H ∗ v · w  w v v v (0, 0) = −  ∇H · w p

(20)

(21)

(22)

(23)

Once w uuu , w uuv , w uv v and w v v v at p are obtained, we can compute an angle θ that a Dp makes with u. Therefore the limit vector at p is cos θ u + sin θ v. Remark 5.1. It should be noticed that Eq. (20) needs to be derived from the full representation of Eq. (16) instead of its simplistic form in Eq. (18). Similarly, for the remaining equations, complete derivatives need to be carried out. 6. Integration of LOC Given a starting umbilical point p and its one Dp as the starting tangent vector, Eq. (4) can now be solved as an initial value problem. We numerically solve it by a Runge–Kutta–Fehlberg integration with adaptive step size control in this article. In this section, we will discuss some key issues to integrate a line of curvature. 6.1. Principal directions Expressions of Eqs. (3) need to be adapted for practical numerical calculation in some degenerative situations. Here we take T1 for an example. √ √ Let K = − V + V 2 − U W and L = − V − V 2 − U W , then we have

 T1 =

K Hz U Hz −K H x − U H y

 (24)

Eq. (24) ordinarily works well. An exception is that T1 trends to vanish, i.e. K → 0 and U → 0, causing T1 to be submerged by numerical error. It corresponds to the case that the plane determined by T1 and T2 is parallel with the x–z plane. and not all of U , V and W are zeros simultaneously at a non-umbilical point, To address this issue, we note that UK = W L so we employ another expression for this situation:

 T1 =

W Hz LHz −W H x − L H y



6.2. Integration of LOC We now set out to address several fundamental issues in LOC integration based on Eq. (4).

(25)

934

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

(a) Fix a terminal

(b) Bidirectional integration by 2-piecewise

Fig. 6. Two technical issues of LOC integration.

6.2.1. Terminal umbilical points As stated in Section 4.3.2, the neighborhoods of umbilical points are detected simultaneously during integration of LOC. Once a new umbilical point is met, it is preserved, e.g. pushed into the stack. It will then be checked if it is near the terminal of the integration when the present solution point is very close to it. The verification is intuitive. As illustrated in Fig. 6(a), a solution point of pn is close to a potential terminal umbilical −−−−−−→ −−−→ point pu , and pn−2 and pn−1 are its immediately preceding solution points. We compare the vectors pu pn−1 , pu pn and Dpu even if they are almost collinear; if it is known to be within some tolerance threshold, pu is taken as a terminal of this line of curvature; or else the integration is continued. 6.2.2. Projection operation Because all points p concern the 2-manifold of H (x) = 0, a projection operation is employed to settle them onto H (x) = 0. The projected point p of p is one that realizes the minimum of |p − x|2 on H (x) = 0. The iterative projection involves a non-linear optimization. We can deal with such an optimization problem in much the same way as that of an umbilical point in Section 4.4. 6.2.3. Piecewise integration The basic framework in Fig. 4 could be unsatisfactory for an integral curve covering a long and curved distance, because the cumulative error could lead to a failure to approach the desired terminal umbilical point by virtue of the high density and singular distribution of LOC near it. As delineated in Fig. 6(b), pu0 pu1 is a real line of curvature; but the integration starting from pu0 is lost near pu1 and reaches another point, say p1 , especially when pu1 is a lemon or a monstar. In order for the framework to work more accurately, an alternative strategy is to use piecewise integration, i.e. the integration is performed piecewise by segmenting the integral curve into several parts rather than a whole one. We still cannot, however, estimate the dividing points inside a line of curvature in implementation, with the exception being extreme points, i.e. the initial and the terminal umbilical points respectively. Therefore, we integrate piecewise per segment only by 2 in the experiments of this article. This can be done by reversing the integration starting from pu1 and merging it with the former from pu0 . Because a certain part of the resulting line of curvature from pu1 , say pu1 p0 overlaps that from pu0 (i.e. pu0 p1 ) in their middle parts, or intersects in a sense of smaller than a preset threshold, it involves finding a joint to connect both within their overlapping portion. Suppose p is such an intersecting point, and we replace pu0 pp1 by pu0 ppu1 as the final between pu0 and pu1 . On the other hand, we should distinguish overlapping from the normal intersection between two LOC. Fortunately, it can be done simply and intuitively since the corner of two normal crossing LOC meet at right angles but is zero for the intersection of interest. 7. Numerical results In this section, we demonstrate our method with several examples. In the aforementioned sections, some threshold parameters are described implicitly to control the behavior of the framework. These parameters, of course, cannot be validated everywhere. In our experiments using double-precision arithmetic, they typically are designed as follows, empirically in part: Criterion for an umbilical point. In this article, it is g (x). A precise umbilical point is crucial to trace a separatrix as correctly as possible. The associated arithmetic operations of g are not complicated and very high precision can be archived. But it is confined to that of f in Newton’s method, which is much lower. The threshold value to determine an umbilical point is set to be 10−20 . Criterion for an ordinary umbilical point. We find in our experiments that the above for an umbilical point, though very small already, can impact greatly on the coefficients of Eq. (6) at a perfect star, which are all zeros. Hence, the threshold value for an ordinary umbilical point should be set larger. An umbilical point is ordinary if one of the cubic coefficients is no less than 10−5 . A terminal umbilical point. A threshold angle to determine a terminal in Section 6.2.1 is 5◦ in degree.

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

(a)

(b)

(c)

(d)

935

(e)

(f)

Fig. 7. P graph of a rounded octahedron. The pink point is under consideration for LOC through it. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

(b) Zoom-in of (a)

(c) Zoom-in of (d) (a) A bottom corner

(d) A side face

Fig. 8. The micro-structures of P graph of the RBF surface. At a bottom corner: Three umbilical points located at the corner, with rounded coordinate numbers, from top to bottom as (9.0847, 9.0847, −8.9970) T , (9.0632, 9.0632, −9.0470) T and (9.0156, 9.0156, −9.1532) T , respectively. On one side face: Some maximal or minimal LOC come very close to each other, respectively.

Overlap in piecewise integration. The threshold value for overlap in piecewise integration entirely relies on the P graph of the surface itself. Given the size l of domain bounds or the bounding box of the surface, two segments of two LOC in piecewise integration are considered to intersect if the minimal distance between them is within 10−4 l. In the following experiments, umbilical points are solved by the powerful solver for algebraic equations provided by MAPLE if H (x) is polynomial in terms of x. 7.1. Illustrative examples The first example is a rounded octahedron x4 + y 4 + z4 + 6x2 y 2 + 6 y 2 z2 + 6z2 x2 − 1 = 0. There are 14 umbilical points



1 1 1 ,± √ ,± √ on this surface, i.e. ± √ 4 4 4 21

21

21

T

, (0, 0, ±1) T , (0, ±1, 0) T and (±1, 0, 0) T , as illustrated in Fig. 7(a). Its P graph is

drawn in Fig. 7(b), where umbilical points lying on the coordinate axes are perfect stars, and the others stars. We read off from Fig. 7(b) that the P graph is relatively simple. We can evaluate the distribution of LOC without too much difficulty. We take a pink point in Fig. 7(c) for an example. In this figure, the P subgraphs of maximal and minimal LOC encircling it are displayed in blue and red lines respectively. In order to estimate the line of minimal curvature through this point, we first think of its surrounding P subgraph of minimal LOC, a space quadrilateral whose vertices are all stars. According to the pattern of a star and the symmetry of the surface, we derive that it is also a space quadrilateral shown by a closed red dashed line in Fig. 7(d). As for the line of maximal curvature, we observe that the surrounding P subgraph of maximal LOC is a space quadrilateral too, but two vertices of which are perfect stars. Therefor it is a curve ending at these two perfect stars shown by a blue dashed line in Fig. 7(d). The net of LOC inside the regions of the two subgraphs are illustrated partly in Figs. 7(e) and 7(f) viewed from different angles. The second example is constructed as a RBF implicit surface, a curved tetrahedron with a concave side (see Appendix A for construction detail). The initial guesses are shown in Fig. 3(a), some of which are only estimates of local minima of g with respective to the sparse mesh. Fig. 3(b) demonstrates the resulting real umbilical points as seed points after optimization. Though they fail to contain all umbilical points, the tracing procedure can compensate for it, as shown by the P graph in Fig. 3(c), including four types of umbilical points concerned in this paper. Two classes of micro-structures, however, are not observable in Fig. 3 because their characteristic lengths are much smaller than the scale of the P graph. One class is associated with each bottom corner crowded by 3 umbilical points, say two stars and a monstar, as illustrated in Figs. 8(a)

936

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

Fig. 9. The flattening P graph the RBF surface.

(a) P graph

Fig. 10. The LOC of the RBF surface.

(b) Appearance of a pair of lemons with the annihilation of a perfect star Fig. 11. A rounded octahedron under small perturbation.

and 8(b). The second is with the monstars on its side faces, respectively; see Figs. 8(c) and 8(d). To help understanding this graph, we have pierced the closed RBF surface through the perfect star on the bottom face (see Fig. 8(a)), and then flattened it onto the plane as in Fig. 9. We observe that the points on side faces, e.g. pu1 (refer to Fig. 8(d) too), seem to be in possession of seven degrees instead of six. We argue that pu1 is actually a monstar and the additional edge pu0 pu1 is not a separatrix of pu1 ; pu0 pu1 is a directed edge only from pu0 to pu1 and happens to end at pu1 , leading to the fact that the P graph cannot be traversed only with seed points on the side faces. The net of LOC in the RBF surface is much more complicated, and yet its P graph is also able to assist us in carrying out the LOC analysis of any surface point. For clarity and simplicity, we take a pink point p as an example in the flattened P graph of Fig. 9. The line of maximal curvature is taken into consideration first. The local P subgraph of maximal LOC encircling p properly is an annularity with a hole, whose vertices are either stars or monstars, so that we can infer the line of maximal curvature through p as a cycle enclosing the hole, displayed as a blue dashed line in Fig. 9. Analogously, we can outline the line of minimal curvature through p topologically as the red dashed cycle in Fig. 9. Its corresponding LOC on H (x) = 0 are illustrated by dashed lines in Fig. 10, in white alternating with blue and red respectively. All elements of a P graph concerned in this paper are present to this RBF surface, i.e. four types of umbilical points of U , connections, loops and cycles in F1 and F2 . The third example is a perturbed counterpart of the rounded octahedron in Example 1,3 i.e.

0.9x4 + 1.1 y 4 + 0.8z4 + 5.9x2 y 2 + 6.1 y 2 z2 + 6.2z2 x2 − 0.9 = 0 3

We thank one reviewer for suggesting it as an example.

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

(a) P graph on the surface

937

(b) Flatten graph on the x– y plane

Fig. 12. P graph of a Gaussian blend surface. A vertex at infinity is actually replaced in P graph by the intersection point of a white ball with the domain bounds in implementation. For clarity and perception, the separatrix of a line of minimal curvature is greened.

(a) The simplistic sketch of P graph

(b) A hidden c4

Fig. 13. A sketch map of P graph of the Gaussian blend and its application to a new discovery of c4 , which fails to be calculated by our current framework.

There are 20 umbilical points on this surface and its graph is shown in Fig. 11(a). The perfect stars in Fig. 7(b) die by producing two lemons respectively, as demonstrated in Fig. 11(b). This example also verifies that a perfect star is unstable, as we have learnt. We can estimate its LOC patterns in the fashion similar to Examples 1 and 2. Let us finally consider a surface of a blend between two translated and scaled Gaussian distributions, defined by the following equation:









z − (1 − x)2 exp −x2 − ( y + 0.5)2 + 0.9(1 − y )2 exp −x2 − y 2 = 0 This is an open surface. We restrict the computation domain within (−5, 5) × (−4, 3) × (−5, 5) to calculate its P graph, as exhibited in Fig. 12(a). We flatten it in the parametric space and change color of a red separatrix to be green for clarity and perception as demonstrated in Fig. 12(b), i.e. the line of minimal curvature relevant to pu2 displayed in Fig. 13(a). But this flattening counterpart still gives us an impression upon the whole since the edges of this P graph crowd each other too much to observe its behavior, especially those specious loops with red and green alternating. To make it capable of providing with clear LOC patterns, we redraw the planar outline of P graph in Fig. 13(a) faithfully according to Fig. 12(b). There are four lemons, pu i , i = 1, . . . , 4, on this surface patch; and their relevant separatrices are denoted by small letters from a to g. A separatrix of pu 1 , b, first winds around and away increasingly from pu 1 pu 2 ; then, after gyrating around it nine times, crosses d for the last time and moves away from it finally. Similarly, the separatrix of pu 2 , c, spirals away and finally escapes. To simplify the sketch map, only one round is exhibited for b and c in Fig. 13(a), respectively. Now it is observable that those red-green-alternating loops are not loops at all. Nevertheless, the P graph still fails to unveil one important line of curvature as a boundary of two groups of lines of maximal curvature. We observe that we can read the behavior of lines coming of f — starting from the left side of f , moving upward and then crossing g to the other side and finally passing downward through the spirals of b, c and d between d and f respectively, as illustrated with the lines of c1 and c2 in Fig. 13(b). On the other hand, the calculation

938

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

(a) Edge 1

(b) Edge 2

(d) Edge 4

(e) Edge 5

(c) Edge 3

Fig. 14. Behavior comparison of f (x), g (x) and h(x), where x component is used for f . They are all non-negative on the surface; and f = 0 (resp. g = 0 or h = 0) means an umbilical point of x. They are evaluated on the RBF surface in Example 2, listed in Fig. 9.

shows that a line of curvature coming off d, say c3 , is a spiral winding around d clockwise and outward. Owing to the different patterns of c2 and c3 , there should be another line to separate them. The thick dashed line of c4 is the very line of curvature in Fig. 13(b). Our present framework in this paper, however, cannot extract freely a multiplex vertex like c4 yet. The remarkable LOC patterns of the Gaussian blend surface are the helix structures of minimal LOC around pu 1 pu 2 , say b and c, and of maximal LOC around pu 2 pu 3 , say c3 , respectively. 7.2. Comparison of f (x), g (x) and h(x) We have pointed out in Remark 4.2 that any of f (x), g (x) or h(x) could serve as a measure function for an umbilical point. What are their differences? In this subsection, we compare their behavior using Example 2. For this purpose, we draw graphs over the sampling points along five edges incident with pu in Fig. 9, indexed by 1, 2, 3, 4, and 5, respectively. These graphs are shown in Fig. 14(a)–(e), respectively. By the symmetry of the surface shape, the shapes of Edges 1 and 3 are the same; the front part of Edge 4 is symmetrical to the back part. But f cannot reflect these facts, as indicated in Fig. 14(a), (c) and (d), respectively. Therefore, the global performance of g and h are better than f since they are dependent of the local shape of the surface instead of variable position. On the other hand, h is usually closer to zero than g at non-umbilical points in these figures but the opposite occurs when approaching an umbilical point. Our experiments show that g is usually more efficient in distinguishing umbilical points from non-umbilical ones compared to h, so that we take g for choice. The readers can refer to the electronic version of this article for details of these figures. 8. Discussion and future work We present a framework in this paper to calculate LOC defined on an implicit surface and elaborate on its numerical computation. In particular, we emphasize the generation of the P graph capturing essential topology of principal configuration. Via such a graph, we can take the advantage of 2-piecewise integration to improve the accurate computation of a connection or a loop. P graph proves to be particularly useful for investigating LOC, such as estimating the shape of a line of curvature and correcting its topology. Our framework also fits naturally to parametric surfaces only with adjustment corresponding to computation of its differential quantities in parametric manner, but the latter, of course, will be more efficient since it allows operations on a surface to be performed as if it is flat. It could also be used to guide the further computation of LOC, e.g. providing information on how to place an initial cycle on the surface in relaxation algorithms. A limitation of our method is that the seed points are sought by a brute-force algorithm combining an optimization procedure if H (x) is not polynomial. As pointed out in Section 7.1, it is not always feasible to traverse a P graph from a

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

939

certain set of vertices, so it is critical to select a proper collection of starting vertices for the tracing procedure to search for the P graph. Our method still belongs to heuristics, whose complete success cannot be assured theoretically in any case. Fortunately, this limitation can be abated to some extent in applications by the redundancy of the seed points. Our experiment in Example 2 indicates, for instance, that any ordinary umbilical point, including the monstars on the side faces, can actually be employed as a seed point to trace the whole P graph, due to the tracing of umbilical points based on their trap neighborhoods instead of P graph itself. A main difficulty of LOC computation rests with the reconstructed line diverging from the true one because of computational errors. Via a P graph, we perform the numerical integration by 2-piecewise in opposite directions from two simplex vertices respectively, say a bidirectional edge. The numerical results show that our method can work much better than simple numerical integration techniques. But a potential danger arises from umbilical points beyond the domain bounds, because the algorithms do not handle those points outside them. Thus, our framework is more applicable to closed and symmetrical surfaces. Our framework can also be generalized to higher order umbilical points, although they are unusual from the viewpoint of generic surfaces. One direct approach is to treat a high order umbilical point independently as ordinary; but their related differential geometry needs to be developed accordingly, which is quite difficult, and the stability of high order computation is a problem too. A similar question exists related to degenerative umbilical points, e.g. a degenerative cubic form, which we do not consider in this paper since they are not computationally robust. For instance, after optimizing an initial guess of a perfect star with threshold value for an umbilical point as smaller as 10−20 , the maximal absolute value of coefficients of Eq. (6) at the solution umbilical point can arrive nearly at 10−5 , much larger than 10−20 . This is a reason for setting the threshold value of coefficients as 10−5 . We believe that the challenge to handle more types of umbilical points rather than only of ordinary is to address the issues of the computational precision and robustness. A generic framework for principal configuration is still an open problem. Our current one can still handle P graph with simplex vertices. We close with a few open problems encountered in our study: Problem 1. How to design a surface with prescribed P patterns rather than a trial and error procedure? It is useful to test if a proposed method is robust in practice. Problem 2. How to detect a multiplex vertex in P graph, which is linetype as a connection, a loop, a cycle or their union? The potentially topological multiformity of P determines their difficulties. There still remains much work to be done. Problem 3. How to distinguish between a cycle and a spiral? Sometimes hyperbolic property is imposed on a cycle such that, oriented properly, it is an attractor of nearby lines for stable structure or stable computation (Gutiérrez and Sotomayor, 1983, 1990; Cazals and Pouget, 2005). However, it is a non-trivial task to distinguish a non-hyperbolic cycle and a spiral line of curvature spiralling densely. We see this paper as an important step towards a comprehensive framework for computing and analyzing LOC. Algorithms exploiting the structure of P with multiplex vertices are planned to be developed in the further work. Acknowledgements The authors would like to thank the anonymous reviewers for their constructive comments. This work was supported in part by National Natural Science Foundation of China under Grant Nos. 60672148, 60872120; in part by National High Technology Development 863 Program of China under Grant No. 2008AA01Z301. WuJun Che was also funded jointly by French–Chinese Foundation for Sciences and Their Applications (FFCSA), and Chinese Government Scholarship of China Scholarship Council (CSC). Appendix A. The construction of a RBF implicit surface The implicit function of Example 2 is taken in following form

H (x) = p (x) +

28 



λ i φ |x − x i |

(A.1)

i =1

where p is a polynomial of second degree and the radial basis function is φ(r ) = r 3 . Given a set of non-zero-valued off-surface points xi , total 28, they are divided into two groups, the inner layer and the outer layer, approximately equally spaced on iso-surfaces at −0.2 and 0.2 to that of interest at 0, respectively. The points of the inner layer are (−9, −9, −9), (9, −9, −9), (−9, 9, −9), (9, 9, −9), (−9, −9, 9), (9, −9, 9), (−9, 9, 9), (9, 9, 9), (10.5, 0, 0), (−10.5, 0, 0), (0, 10.5, 0), (0, −10.5, 0), (0, 0, −10.5), (0, 0, 5.25). 1 1 1 1 1 1 ), −(9 + √ ), −(9 + √ )), (9 + √ , −(9 + √ ), −(9 + √ )), (−(9 + The points of the outer layer are (−(9 + √ 5 3

1 √

5 3

), 9 +

1 √

5 3

, −(9 +

1 √

5 3

)), (9 +

1 √

5 3

,9 +

1 √

5 3

, −(9 +

5 3

5 3

1 √

5 3

)), (−(9 +

1 √

5 3

), −(9 +

5 3

1 √

5 3

), 9 +

5 3

5 3

1 √

5 3

), (9 +

1 √

5 3

, −(9 +

1 √ ), 9 + 5 3

940

X.P. Zhang et al. / Computer Aided Geometric Design 26 (2009) 923–940

1 √ ), 5 3

1 1 (−(9 + √ ), 9 + √ ,9 + 5 3 5 3 (0, 0, −10.7), (0, 0, 5.45).

1 √ ), 5 3

(9 +

1 √ ,9 5 3

+

1 √ ,9 5 3

+

1 √ ), 5 3

(10.7, 0, 0), (−10.7, 0, 0), (0, 10.7, 0), (0, −10.7, 0),

The surface of interest is the locus of points where Eq. (A.1) is zero. References Alliez, P., Cohen-Steiner, D., Devillers, O., Lévy, B., Desbrun, M., 2003. Anisotropic polygonal pemeshing. In: Proceedings of ACM SIGGRAPH, pp. 485–493. Beck, J.M., Farouki, R.T., Hinds, J.K., 1986. Surface analysis methods. IEEE Computer Graphics and Applications 6 (12), 18–36. Bloomenthal, J., 1994. An implicit surface polygonizer. Graphics Gems IV, 324–349. Cazals, F., Pouget, M., 2005. Differential topology and geometry of smooth embedded surfaces: selected topics. International Journal of Computational Geometry & Applications 15 (5), 511–536. Che, W.J., Paul, J.-C., Zhang, X.P., 2007. Lines of curvature and umbilical points for implicit surfaces. Computer Aided Geometric Design 24 (7), 395–409. Delmarcelle T., Hesselink L., 1994. The topology of symmetric, second-order tensor fields. In: Proceedings of IEEE Visualization 1994, pp. 140–145. Farouki, R.T., 1998. On integrating lines of curvature. Computer Aided Geometric Design 15 (2), 187–192. Gutiérrez, C., Sotomayor, J., 1983. An approximation theorem for immersions with stable configurations of lines of principal curvature. In: Lecture Notes in Mathematics, vol. 1007. Springer, pp. 332–368. Gutiérrez, C., Sotomayor, J., 1990. Periodic lines of curvature bifurcating from Darbouxian umbilical connections. In: Proceedings of Bifurcations of Planar Vector Fields, 1989. In: Lecture Notes in Mathematics, vol. 1455. Springer-Verlag, pp. 196–229. Gutiérrez, C., Sotomayor, J., Garcia, R., 2004. Bifurcations of umbilic points and related principal cycles. Journal of Dynamics and Differential Equations 16 (2), 321–346. Hallinan, P.W., Gordon, G., Yuille, A.L., Giblin, P., Mumford, D., 1999. Two- and Three-Dimensional Patterns of the Face. A.K. Peters. Helman, J.L., Hesselink, L., 1989. Representation and display of vector field topology in fluid flow data sets. IEEE Computer 22 (8), 27–36. Helman, J.L., Hesselink, L., 1991. Visualizing vector field topology in fluid flows. IEEE Computer Graphics and Applications 11 (3), 36–46. Maekawa, T., Patrikalakis, N.M., 1994. Interrogation of differential geometry properties for design and manufacture. The Visual Computer 10 (4), 216–237. Maekawa, T., Wolter, F.E., Patrikalakis, N.M., 1996. Umbilics and lines of curvature for shape interrogation. Computer Aided Geometric Design 13 (2), 133– 161. Patrikalakis, N.M., Maekawa, T., 2002. Shape Interrogation for Computer Aided Design and Manufacturing. Springer-Verlag. Porteous, I.R., 2001. Geometric Differentiation: For the Intelligence of Curves and Surfaces, second ed. Cambridge University Press. Scheuermann, G., Tricoche, X., 2005. Topological methods for flow visualization. In: Hansen, C.D., Johnson, C.R. (Eds.), The Visualization Handbook. Elsevier. Sotomayor, J., Garcia, R., 2007. Lines of curvature on surfaces, historical comments and recent developments. arXiv:0712.1585v1 [math.DG], http://arxiv. org/abs/0712.1585. Thirion, J.P., 1996. The extremal mesh and the understanding of 3D surfaces. International Journal of Computer Vision 19 (2), 115–128.