00001.tif

PDF

Download

Open document

Transcript (if available)

Content CONVEX HULL OF SURFACE PATCHES: CONSTRUCTION AND APPLICATIONS by Chao-Kuei Hung A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Computer Science) December 1995 Copyright 1996 Chao-Kuei Hung UMI Number: D P 22898 All rights reserved INFORMATION TO ALL U SER S The quality of this reproduction is dep en dent upon the quality of the copy submitted. In the unlikely event that the author did not sen d a com plete manuscript and there are m issing p a g es, th ese will be noted. Also, if material had to be rem oved, a note will indicate the deletion. UMI Dissertation Publishing UMI D P 22898 Published by ProQ uest LLC (2014). Copyright in the Dissertation held by the Author. Microform Edition © ProQ uest LLC. All rights reserved. This work is protected against unauthorized copying under Title 17, United S ta tes C ode ProQ uest LLC. 789 E ast E isenhow er Parkway P.O. Box 1346 Ann Arbor, Ml 4 8 1 0 6 - 1346 UNIVERSITY OF SOUTHERN CALIFORNIA THE GRADUATE SCHOOL UNIVERSITY PARK LOS ANGELES, CALIFORNIA 90007 This dissertation, written by under the direction of hl& Dissertation Committee, and approved by all its members, has been presented to and accepted by The Graduate School, in partial fulfillment of re quirements for the degree of DO CTO R OF PHILOSOPHY Ph.D CpS ’ <\S C.Piao' Ku t\ Hung. Dean of Graduate Studies Date November 28, 1995 DISSERTATION COMMITTEE Chairperson D edication To my parents: whose encouragement made the completion of this thesis possible. A cknow ledgem ents Many ideas in this thesis were first conjectured by Prof. Doug Ierardi, my advisor. Doug is also my oracle for pointers to other information sources, besides being a wonderful source of information himself in subjects ranging from algebra to Linux, from philosophy to Deep Space Nine. The interesting application in stable poses and capture regions would not exist without the helpful discussions with and references from Prof. David Kriegman, who pioneers the study in this field. I am fortunate to have Prof. Francis Bonahon in my committee for consulting mathem atical problems. Prof. Aristides Requicha’s input lead to the discussions on practical issues and finally an implementation, which in turn enhanced my familiarity with C++. I owe my gratitude to Prof. Ken Goldberg, who first posed the convex hull problem to me. He and Anil Rao subsequently provided me with several im portant references on grasping and width functions. Profs. Ram Nevatia and Ming-Deh Huang, along with my other committee members, helped shaping and focusing the thesis. I was first initiated into the world of duality by Prof. Alfred Inselberg, although in a different setting. His “tangential” ideas in parallel coordinates had a subtle but definite impact on the design of the duality in this thesis. I am grateful to him, Avijit Chatterjee, and Yui-Bin Chen for much valuable advice during the early stage of my Ph.D. study. Hsun-Ko Chan and Shih-Hui Lu, a warm and wonderful couple, pushed me gently back to the track when Professor Inselberg had to leave USC and before I found another supportive advisor. Many friends around me, especially Mourad Zerroug, Chih-Wen Hsueh, and dear M artha Escobar, have cooled me down and encouraged me during the typical manic-depression cycles of a Ph.D. student. Finally, it has always been my wish to express my gratitude towards the gurus on the Internet who substantiate the concept of free software. To me, the collaborative nature of such programs as Linux, gnu*, RLaB, PARI/gp, and many others result in their influential and lasting contributions to the computing community. It is they who have been practicing software reuse and citing some of the most relevant papers without necessarily advertising the fact, and on to whose shoulders I wish to climb. C ontents D ed ication ii A cknow ledgem ents iii List O f Figures vii A b stract ix 1 Introduction 1 1.1 Problem Statements and a Brief S urvey....................................................... 1 1.1.1 Convex Hulls of Curved O b je c ts....................................................... 1 1.1.2 Stable Poses and Capture Regions................................................... 2 1.1.3 W idth Function and Tight Bounding Boxes ............................... 3 1.1.4 Practical I s s u e s ..................................................................................... 4 1.2 An Outline ........................................................................................................ 4 2 A N ew Look at th e Planar Problem s 6 2.1 Simple Polygons Bounded by Line S e g m e n ts............................................. 6 2.1.1 P rep ro cessin g ........................................................................................ 8 2.1.2 The Main L o o p .................................................................................... 8 2.1.3 Correctness and C om plexity............................................................. 9 2.2 Simple Polygons Bounded by Algebraic Curve Segments ..................... 10 2.2.1 P rep ro cessin g........................................................................................ 10 2.2.2 The Main L o o p .................................................................................... 12 2.2.3 Correctness and C om plexity............................................................. 13 2.3 Hints and Heuristics in the Design of Algorithms for Higher Dimensions 14 3 C lassical Forms o f D u ality and M ath P relim inaries 19 3.1 Duality in Projective G e o m e try .................................................................... 19 3.1.1 Notations and Terminology ............................................................. 19 3.1.2 Hypersurfaces of the Second O r d e r ................................................ 21 3.2 Pedal Curves and Pedal Surfaces ................................................................ 22 3.2.1 Definitions and P ro p e rtie s ................................................................. 22 3.2.2 Direct Generalizations and a V a r ia n t............................................. 22 3.3 Polar S e t s ............................................................................................................ 24 3.4 Topologically Derived Sets of the Polar S e t ................................................ 25 3.4.1 Partitioning Hyperplanes w.r.t. a S e t ............................................. 25 3.4.2 Topological Properties of the Partitioning ................................... 26 3.4.3 Topologically Derived S e t s ................................................................. 29 3.5 Curvilinear C o m p lex es..................................................................................... 29 4 D uals o f Curved Polyhedra 31 4.1 Object M o d e l...................................................................................................... 31 4.2 Duals of F aces...................................................................................................... 32 4.2.1 Duals of Vertices and F a c e ts .............................................................. 33 4.2.2 Duals of Faces of General D im ensions............................................. 34 4.3 The Cell T heorem ............................................................................................... 34 4.4 Computations and General Properties of D u a l .......................................... 36 4.4.1 General P r o p e rtie s ............................................................................... 37 4.4.2 Implicit F a c e s ......................................................................................... 37 4.4.3 Param etric F a c e s .................................................................................. 41 5 Q uadratic H ypersurface Patches 42 5.1 Non-Degenerate Hosts and the Extent of the D u a l................................... 42 5.2 Generalized Projective D u ality ........................................................................ 43 5.3 Degenerate H o s ts ............................................................................................... 44 5.4 Lower Dimensional H o sts................................................................................. 46 5.5 Degeneracies and Singularities....................................................................... 49 6 C onvex H ull A lgorithm s 52 6.1 Dualizing the I n p u t........................................................................................... 53 6.1.1 Choosing the Coordinate S y s t e m .................................................... 53 6.1.2 Breaking the Faces at S in g u la ritie s................................................ 53 6.1.3 Bounding the Work Space in the Dual S p a c e ................................ 55 6.2 Constructing the Polar S e t .............................................................................. 55 6.3 Recovering the Convex H u l l ........................................................................... 57 6.4 General Curved Polygons in E 2 .................................................................... 62 6.5 Avoiding Hull C o n s tru c tio n ........................................................................... 63 6.6 A Practical A pproach........................................................................................ 64 6.6.1 R a tio n a le ............................................................................................... 64 6.6.2 Precision and Error M e a su re s........................................................... 65 6.6.3 Other Practical Considerations ....................................................... 68 7 Stable P oses and C apture R egions 71 7.1 The Direct A p p ro ach ........................................................................................ 71 7.1.1 Summary of the Case A nalysis........................................................... 72 v 7.1.2 Distinguishing Local Minima from Other Critical Poses .... 73 7.1.3 Coincidences and Singularities.......................................................... 75 7.2 The Dual A pproach......................................................................................... 78 7.2.1 Case A n a ly sis........................................................................................ 78 7.2.2 Coincidences ........................................................................................ 79 7.2.3 Singularities........................................................................................... 80 7.2.4 A Complete A lg o rith m ....................................................................... 82 7.3 Capture Regions ............................................................................................ 82 7.4 Practical C o n sid eratio n s............................................................................... 84 8 W id th Functions and T ight B ounding B oxes 88 8.1 W idth F u n c tio n .................................................................................................. 88 8.2 Planar Tight Bounding B o x e s ....................................................................... 91 9 C onclusions 95 9.1 C o n trib u tio n s................................................................................................... 95 9.2 Future R esearch................................................................................................ 96 R eference List 97 A p p en d ix A Implementation N o t e s ..................................................................................................102 A .l Summary ............................................................................................................ 102 A.2 Object Specification ........................................................................................ 103 A.3 Functionalities and L im ita tio n s .................................................................... 105 A.4 Submatrices: to Copy or to R eferen ce?.......................................................... 106 V I List O f Figures 2.1 Convex hull of a (planar) simple p o ly g o n .................................................... 7 2.2 The cases of the next edge to be examined. The dashed lines indicate the (possible) edges on the polygon; the dotted lines separate the cases. The origin is marked by an “x.” ...................................................... 16 2.3 Convex hull of a curved p o ly g o n .................................................................... 17 2.4 Left: backward edges going beneath—ways to break the snare situa tion. Right: backward edges going above..................................................... 18 3.1 Shaded area represent loci of the perpendicular foots of the straight lines (hyperplanes in E 2) which necessarily cut K ......................................... 28 3.2 p ed H is on the boundary of the pedals of the cutting hyperplanes of K ............................................................................................................................. 28 4.1 A cycloid and its dual It can be decomposed into three smooth open cells—the two end points and the rest of the curve. The two straight lines are images of the vertices. The dual of the rest is the a curve with two asymptotes, which correspond to the two supporting lines from the origin to the curve............................................................................. 35 4.2 An S-shaped figure and its dual....................................................................... 35 4.3 A cardioid and its dual. It need be decomposed into at least two smooth open cells—the cusp the the rest of the curve. The horizontal line is the dual of the cusp. The dual curve begins and ends at the inflection point, at which the horizontal line is the tangent.........................38 4.4 A dimpled surface—the zero set of (4a:2 + 3y2)2 — 4a:2 — 5y2 + 4z2 — 1. 38 4.5 Dual of the dimpled surface. The two “umbrellas” turned inside out correspond to tangent planes to points inside the two dimples. The ap parent sharp cuspidal edges of the umbrellas are duals of two parabolic curves inside the dimples.................................................................................. 39 6.1 The union of two ellipsoids and a truncated cone....................................... 59 6.2 The polar set of the object, recovered as the cell containing the origin in the arrangement of dual surfaces....................................................................59 vii 6.3 The convex hull of the original object, recovered from the dual of the polar dual (previous figure). Note that / j is planar since it is dual to the triple intersection (vertex) in the previous figure. Other new sur faces ( /2, fzi and fi) are developable, and arise from the corresponding edges of the polar set, .................................................................................... 60 6.4 A “big” feature that should not be missed in terms of convex hull construction.......................................................................................................... 66 6.5 A “small” feature that could be missed if the the precision of the approximation is low.......................................................................................... 66 6.6 An e-error circle and its dual hyperbola........................................................ 70 7.1 A PVC pipe fitting.............................................................................................. 74 7.2 Another PVC pipe fitting.................................................................................. 74 7.3 Collinear “two-point” contacts......................................................................... 77 7.4 Singularity caused by symmetry...................................................................... 77 7.5 Equi-potential lines (in the parameter space) of the two-piece PVC pipe example........................................................................................................ 83 7.6 Equi-potential lines (in the parameter space) of the three-piece PVC pipe example................................................................... 83 7.7 The radius function of the three-piece PVC pipe example............................87 8.1 An “object” K consisting of two ellipses and a line segment; shown along with its convex hull................................................................................. 90 8.2 The dual of K . The solid line indicates the boundary of its polar set, K * ........................................................................................................................... 90 8.3 The only discontinuities in the first derivative of the width function occur at either those of the radius of i f as a function of the polar angle, or those of the same function translated by t v . (All three functions are shown.) ........................................................................................................ 92 8.4 The area of the tightest bounding box as a function of the polar angle. Only one period (0 to 7 t / 2) is shown.............................................................. 92 8.5 The tight bounding boxes of maximum and minimum area as found from the previous picture.................................................................................. 94 vm A bstract We solve the problem of constructing the convex hull of piecewise smooth objects in 3-D. This work was motivated by a need to compute the width function of parts in a manufacturing environment where an oblivious grasping algorithm drives a pair of parallel jaw grippers to orient parts on a conveyor belt. We show th at constructing a planar graph describing the convex hull takes randomized expected time 0 ( N 2+t) for arbitrary e > 0, if the input consists of N algebraic faces of bounded degree each bordered by a bounded number of subfaces. We use duality in our construction. This approach allows us to solve two other related problems—finding the stable poses and capture regions of a piecewise smooth object. Such information can be utilized in a model-based vision system to reduce the complexity of the algorithm for recognizing immobile parts. The duality approach provides other advantages as well, such as simplifying the treatm ent of algebraic singularities and unifying the treatm ent of faces of different dimensions. In practice many industrial parts can be described solely in terms of quadratic surface patches cut out by planes. We show that in this case the algebraic overhead in constructing the convex hull is significantly lower than estim ated in the general algorithm. For example, the algebraic complexity of the systems of equations in the stable pose problem is reduced by almost two orders of magnitude compared to existing literature. In this specialized context, it is feasible to adopt a simple scheme for safely processing coincidences, as demonstrated in an implementation in C++ based on the dual approach. ix C hapter 1 Introduction 1.1 P roblem S tatem ents and a B rief Survey Computational geometers have exhaustively studied the problem of constructing convex hulls of discrete point sets in various dimensions. Indeed, Chazelle give optimal algorithms to this problem for arbitrary dimension in [8]. The convex hull of curve segments, however, have received little attention until late eighties. There is little existing literature which addresses the issues in solving the three dimensional, curved version of this problem. We solve the problem of constructing the convex hull of piecewise smooth objects in 3-D using duality. We then apply these methods to two other related problems arising from computer vision and manufacturing. 1.1.1 C onvex H ulls o f C urved O bjects Convex hull algorithms for the curved version of a simple planar polygon have been proposed independently by several researchers. Shaffer and Van Wyk [43], Dobkin and Souvaine [13], and Bajaj and Kim [5] each give an optimal (linear time) algo rithm , although they use different terms such as piecewise smooth Jordan curves and splinegon. Their emphasis is more on the combinatorial aspect than on the algebraic computations, although such issues are also addressed. In contrast and complemen- tarily, the work of Kriegman, Yeh, and Ponce in [32] focus on the algebraic aspect and implementation issues by studying the convex hull of a single closed algebraic curve in the plane. 1 We present a much simpler algorithm for achieving the same optimal results as in [43, 13, 5]. It is shown further that if the input consists of unordered curve segments along with isolated points, the convex hull can still be computed in optimal time Q (N log N). In extending his results to 3-D and to related problems, Kriegman [30, 31] also discusses the algebraic issues. Schwartz and Sharir [45] show that the complexity of the convex hull of a set N spheres is 0(iV2) in the worst case. Boissonnat, Cerezo, Devillers, Duquesne, and Yvinec [6] later present an algorithm for the convex hull of a set of spheres in higher dimensions. They apply duality in a special way which allows the result to be generalized to homothetic objects. However, we have been unable to locate other literature which give a more general result on this problem. We show that constructing a planar graph describing the convex hull takes ran domized expected time 0 ( N 2+e) for arbitrary e > 0, if the input consists of N algebraic faces of bounded degree each bordered by a bounded number of subfaces. Several operations can be performed on the convex hull without its explicit con struction. Point classification against the hull can be effected in time O (N logN ). The supporting cone from a point to the convex hull is shown to have a structural complexity of 0(iV) and can be constructed in tim e 0 ( N log N). 1.1.2 Stab le P oses and C apture R egions Given a description of a piecewise smooth object and the location of its center of mass, its gravitational potential energy is a function of its “pose.” The stable pose problem seeks to determine those poses of the object at which the potential energy function assumes local minima. This problem occurs in model-based vision systems where such information can be used to reduce the complexity of recognizing a static object resting on a plane. Kriegman pioneered the study of stable poses of piecewise smooth objects in [30]. The combinatorial complexity is not given explicitly but as he pointed out, the dominant factor is the algebraic part. Experimental results show that it takes roughly two days on a Sun SPARCStation 1 to process figures similar to Figures 7.1 and 7.2 using homotopy method for root-finding. 2 We show that, as a corollary to the complexity of the convex hull, the number of stable poses of a piecewise smooth object having N faces each with a bounded number of subfaces, is 0 ( N 2+ C ). As a result of theoretical interest, it is shown that they can be computed in the same amount of randomized expected tim e with an algebraic overhead. A direct extension to the stable pose problem is the computation of the capture regions—the regions in the param eter space within each of which the object will eventually converge to one stable pose starting from any pose. A general planner needs this information in controlling a robot manipulator to place an object in the plane in the presence of mechanical uncertainty. In [31] Kriegman extends the idea to computing the capture regions. The analysis is based on stratified Morse theory, under the assumption that the potential energy function is Morse. Implementations for either closed individual algebraic surfaces or general polyhedra are given. The combinatorial time complexity in the latter case is further shown to be 0 ( N 2). We show that the capture regions can be described by a planar graph of size 0(iV 4+e) and can be constructed in corresponding time. 1.1.3 W id th Function and T ight B ounding B oxes The distance between the pair of tight supporting planes which bound a given object is a function of their orientation. Such information is essential to the algorithm which generates oblivious grasping plans to operate a pair of parallel jaw grippers in a parts feeder. Parts orienting and the width function are extensively studied by Mason, Gold berg and Taylor in [36], Goldberg in [20], Rao and Goldberg in [41], Rao, Kriegman, and Goldberg in [40], and by Chen, Hung, and Ierardi in [11]. An efficient algorithm for finding the overall width of a 3-D point set in time 0 ( N 8^5+€) without construct ing the entire width function is given by Chazelle, Edelsbrunner, Guibas, and Sharir in [10] We study the width function and generalize the idea to the description of the area or volume of the tight bounding box as a function of its orientation relative to the object. 3 1.1.4 P ractical Issues The duality approach eases the treatm ent of several practical issues in the 3-D versions of all of the above problems. The tedious analysis of sixteen cases involving combinations of up to three faces of various dimensions is reduced to th at of three cases. We are able to describe how to give an exhaustive listing of singularities up to a fixed degree in the dual space, and give examples of what they correspond to in the primal space. As a practical substitute for the asymptotically efficient algorithms in treating all three problems, we also propose a simple algorithm for sampling points on the convex hull. An intuitive measure of the precision of sampling density is shown to have a simple geometric counterpart in the dual space. Its output may be further processed to produce a planar graph describing the exact boundary of the convex hull, the potential energy function, or the width function, subject to “small” missing features as determined by the precision measure. In practice many industrial parts can be described solely in terms of quadratic surface patches cut out by planes. We show th at in this case the algebraic overhead in these problems is significantly lower than estimated in the general algorithm. Systems of nine quadratic equations in nine variables are replaced by systems of three quadratic equations in three variables, resulting in a reduction of the total degree from 512 to 8. 1.2 A n O utline We first revisit the already-solved two-dimensional version of the convex hull problem in Chapter 2. We introduce concepts which will be elaborated later on in this simpler context and show how they significantly simplify the logic of the optimal algorithms. The mathem atical foundations of the dual approach begin with a review of several forms of traditional duality in Chapter 3. Their relationships are examined and their strengths combined in Chapter 4 to give a definition of inverted pedal, which retains the connection between the dual and the convex hull and at the same time is constructively defined from a computational point of view. The special case of quadratic hypersurface patches cut by hyperplanes is closely studied in Chapter 5. 4 The results here are essential to the significant reductions in the algebraic aspect of algorithms dealing with typical industrial parts. The main algorithm for convex hull construction and several other algorithms directly related to the convex hull are presented in Chapter 6. The first application is studied in Chapter 7, where Kriegman’s discussions and results on stable poses and capture regions are improved upon. We then investigate the second applica tion of convex hull construction—the width functions and tight bounding boxes in Chapter 8. The contributions and future extensions of this thesis are summarized in Chapter 9. Appendix A are the implementation notes of our C++ program. Also presented here is an elegant design of the m atrix class which provides an efficient and neat interface for submatrix manipulations without the need of keeping reference counts. 5 C hapter 2 A N ew Look at th e Planar Problem s To motivated and introduce our approach, in this chapter we first consider the simple 2-D convex hull problem of ordinary and curvilinear simple polygons. The point set version was already studied in [37, 22, 34] and the curved version in [43, 13, 5]. In both cases the problems have long been known to have time complexity 0 (IV), in the latter case with an additional overhead given by the algebraic complexity of the bounding curve segments. The ideas introduced in our versions of the solutions lead to significant simplifi cations of the control flow over the existing algorithms, especially in the curvilinear case. They also lead to the use of and guide in modifying the classical duality def initions and properties, and eventually enable us to provide a general algorithm in higher dimensions. We also use these algorithms to introduce a few notations and give examples of the algebraic issues raised by non-linearity, which we discuss in details in later chapters. 2.1 Sim ple P olygons B ounded by Line Segm ents We begin by the rectilinear version of the convex hull problem. Let K = {Pi, p 2, . . . , p A ,) be the set of vertices of a simple polygon given in counterclockwise direction around its interior. Figure 2.1 shows a C++ implementation of our algo rithm for finding its convex hull. We employ an auxiliary data structure PoinDStack for storing the vertices in a stack. To save memory usage as well as typing, it is actually the indices to the vertices plus a reference to the array of vertices (class Vector<PointD>) which are stored in the stack. Given an index possibly out of 6 in t right_turn(PointD const & a, PointD const & b, PointD const & c) { retu rn ( b -a ) .cross(c-b) < 0; } c lass PointDStack : public Stack<int> { protected: Vector<PointD> const & host; PointD v (in t i) const { return host [(i+host .n e lts O ) % host .n e lts () ] ; } p u b lic : PointDStack(Vector<PointD> const & h) : host(h) {} in t snared_beneath(int i) const { retu rn v ( to p ( ) ).c ro s s ( v ( i) ) <= 0; } in t prev_concave(int i) const { retu rn head->next & & : : right_turn(v(head->next->data), v (to p ()), v ( i) ) ; } in t rig h t_ tu rn (in t i) const { retu rn : :rig h t_ tu rn (v (i-2 ), v ( i - l ) , v(i)); } in t snared_in_pouch(int i) const { retu rn ! prev_concave(i); } >; PointDStack PolygonHull(Vector<PointD> & vertex) { PointDStack CH(vertex); CH.push(i = x.max); 11 s ta rtin g from an extremal vertex CH.push(stop = ++i % = N); / / and ending at the same vertex while (++i % N != stop) { i f (CH.right_turn(i)) while (CH.snared_in_pouch(i)) ++i; else while (CH. snared_beneath(i)) ++i; while (CH.prev_concave(i)) CH.popO; CH.push(i); } CH.popO; / / remove the s ta rtin g point, which was entered twice retu rn CH; Figure 2.1: Convex hull of a (planar) simple polygon 7 array bound (but not “too negative”), its member function v() retrieves the actual coordinates of the vertex by regarding the array as circular. 2.1.1 P reprocessing To find conv K , the convex hull of K , we first choose a convenient coordinate system in which the origin O lies in in t conv K , the interior of conv K. One obvious choice is to use c k = jj J2i Pi as the new origin. It is also easy to show that b # the center of the tight bounding box B of K whose sides parallel the coordinate axes lies in the interior or (with probability 0) on the boundary of conv K. Though not guaranteed to lie in the interior, the latter alternative does not suffer from skew point distributions. We choose |( c k + bff) as the new origin as it is easy to show that it lies in the box concentric with B and half of its size, and that it lies in in t conv K unless all vertices are collinear. We will assume without loss of generality that px is an extremal vertex, which is guaranteed to appear on the boundary of con vK. For ease of presentation we further assume that the coordinate axes are so rotated that px lies on the positive x-axis, although in reality we will not have to actually carry out this rotation. Finally let the polar angles of the vertices be normalized to fall within [0,27r). 2.1.2 T he M ain Loop The algorithm (Figure 2.1) begins with a stack CH containing pa and p2 with p2 on top, and proceeds in the given order of the vertices. We now state the invariants of main loop. At each iteration on entrance to the body of the main loop let p k be the vertex on the top of CH. For every predecessor p, of p f c (that is, i < k), exactly one of the following two statements is true: • Pj is known to lie in in t conv K even with our limited knowledge up to this point; or • Pj lies on the boundary of a convex “sector” beginning with Opx and ending with Opk, and it is contained in CH in sorted order by its polar angle. This is obviously true just prior to the first iteration of the loop. 8 To show that the above is indeed invariant, let us consider what happens until the end of one iteration of the outer loop. For now we ignore the possibility of having a radial segment by just mentioning th at it can be avoided with probability one if we perturb the origin. Then the following cases cover all possibilities at this point: (See also Figure 2.2) 1. p^_iPf c is counter-clockwise (a forward edge). (a) PkPk+i makes too sharp a right turn that it enters the most recent pouch. (b) p A .pfc + 1 makes too sharp a left turn that it goes clockwise (backward) and gets “snared” beneath JbjZiP*.- (c) otherwise. 2. P;-_iPf c is clockwise (a backward edge). (a) PfcP* ; +1 makes too sharp a right turn that it enters the most recent pouch. (b) otherwise. Note th at there is nothing in case 2 which corresponds to case 1(b). On the other hand both cases 1(c) and 2(b) entitle Pj.+1 to be collected into the stack. It is straightforward to verify that the i f part in the main loop processes cases 1(a) and 2(a) alike, whereas the else part processes case 1(b), both discarding succeeding vertices which cannot possibly lie on the hull. Neither of cases 1(c) and 2(b) is affected by either branch. Thus at the end of the i f statem ent, all interior vertices are discarded in each of the 5 cases, and the current vertex is one which we would like to push on to the stack. Finally we proceed to remove the vertices which are now made concave by the addition of the current vertex. When this stops all vertices known to be interior to the hull up to now are removed and the invariant holds again. 2.1.3 C orrectness and C om p lexity The correctness of this algorithm follows from the established loop invariant and the fact that the iterations end when Pi has for a second time been pushed onto CH. At this moment the invariant can be rephrased as saying that vertices in CH are exactly the hull vertices in sorted order, with p a occurring twice. During 9 the algorithm, discarding vertices and stack operations are sufficient to account for all (inner and outer) loop iterations. This number is obviously O(N) since each vertex contributes at most two. The algorithm is therefore Q(N) since each testing operation takes constant time. We note that this algorithm proceeds around the origin, unlike [22] which pro ceeds along some straight line. Thus there is no need to construct the upper and lower hulls separately. Moreover, the earlier assumption th at p 1 is a hull vertex and th at Opi is the polar axis can in reality be removed. Simple modular arithm etic on the vertex indices as performed by the member function v ( ) solves all of these three problems simultaneously. 2.2 Sim ple P olygons B ounded by A lgebraic C urve Segm ents We term a simply-connected planar region bounded by a finite number of seg ments of algebraic curves as a curved polygon. The previous algorithm can be gener alized to find the convex hulls of curved polygons whose edges are given in a circular order around their interiors. Assume th at each curve segment is given as an implicit equation of total degree at most D in x and y. Again we shall reduce the number of possible cases which arise in the algorithm and simplify the logic without introduc ing global flags whose values are altered by function calls and by which the control flow is determined. Figure 2.3 shows the main loop of a C++ implementation. Here the input is an array of Edges instead of V ertexes, and the auxiliary data structure is now a stack of Edges. We leave out the details of the Edge class and only remind the user that it provides enough information for the geometric computations requested in the main loop, which we shall discuss in details promptly. 2.2,1 P reprocessin g Continuing along the line of thought in the rectilinear version, we shall proceed around the origin. To this end, we first break each edge into pieces at inflection points, so th at each resulting piece is either entirely convex (from the origin towards 10 outside) or entirely concave. We further break each piece into yet shorter pieces such th at along each, the points have either monotonically increasing or monotonically decreasing polar angles. We finally break the pieces again if the polar angles of the points on it have a span greater than or equal to 7 r. This last preprocessing ensures an unambiguous definition of the ordering of two points by their polar angles, where the difference between two angles is normalized to be within (— 7r,7r). By how much did we increase the number of input curve segments in doing so? Suppose one input curve segment 7 is given as the zero set of a degree D polynomial f(x,y) = 0. The first step breaks it at no more than D pieces. The second step further breaks 7 at the turning points/ o n 7 , where the tangent line to 7 passes O and satisfies x-j^f + y-§^f, an equation of degree at most D + 1. By Bezout’s theorem, there are at most D(D + 1) such points. We can tighten the bound, however, if we consider / as a projective curve and perform a projective transform ation which brings O to any point at infinity, say the one on the t/-axis. In the new coordinate system the turning points correspond to the extremal points, or folds along the x- axis, which can be no more than D, the maximum number of points of intersection between 7 and an arbitrary straight line. For the last step of breaking pieces longer than 7 r, we dispose of the new coordinate system used only in the previous proof and go back to break 7 at all of its intersections with the polar axis (in the original coordinate system). This creates at most D more breaking points. In practice we break 7 into monotonic segments of angular span no more than 7 t / 4 . In summary the number of segments does not increase by more than a factor of 4D. We make one last surgery on the input. An extremal point is found and inserted as one more vertex which is guaranteed to appear on the hull. To simplify notations, we will thus assume from now on that the input curved polygon K consists of N segments which have already been pre-processed. In the algorithm recorded in each segment edge [i] are s t a r t and sto p , its two end points whose order is determined by tracing the curve in the order of the edges prescribed by the input. It is then meaningful to speak of an edge as being forw ard or backward, depending on whether the normalized difference between the polar angles of s t a r t and sto p is less than or greater than zero. We rotate the indices of the edges so that the extremal vertex serves as the s t a r t vertex of the first edge and sto p vertex of the last edge, both of which are certainly forw ard. 11 2.2.2 T he M ain Loop As before we will proceed counter-clockwise around the interior of the curved poly gon. In the present case, however, it is necessary to push the edges, not just their end points, to the stack. Moreover, new flat edges, namely segments of the common supporting lines between adjacent hull edges, may emerge. On the other hand, we shall see that concave edges can be treated exactly the same way as the convex or flat ones are as far as the control flow of the main loop is concerned. We first push e d g e [l] onto the stack, noting that even if it gets popped off the stack later on, the replacing edge (which is necessarily a supporting line segment) will always have the same s t a r t vertex. If it happens to be concave, we push instead the lid which covers it, without modifying the concave edge itself.1 Since we are pushing an entire edge on to the stack, we have already processed two vertices and the second edge can be treated as the rest of the edges as described below. At each iteration prior to entering the body of the main loop, let 7 be the edge on the top of CH . Note that 7 may or may not be one of the input edges—it may just be a common supporting line of an input edge and a vertex. In any case let edge[k] be the edge in the original input which generated the sto p vertex of 7 . In other words, edge[k] at least shares the same sto p vertex with 7 if they do not in fact coincide. We claim a loop invariant similar to that in the linear case with an additional property. For every predecessor 7, of 7 (i < k), exactly one of the following two statem ents is true: • 7i is known to lie in in t conv K even with our limited knowledge up to this point; or • at least part of 7, - lies on the boundary of a convex “sector” beginning with the ray from O to edge [ 1] . s t a r t and ending with the ray from O to edge [k] .sto p . Furthermore, all edges on the stack are either portions of some forw ard convex edges in the input, or “lids” of forw ard concave edges, or common supporting line segments of the input edges. *It is very important that we should not replace the concave edge by its lid because later on we need to know the tangent direction at its end points for determining whether the next edge goes “above” or “below” it. 12 Supposing that the invariant holds upon entry of the loop body at some iteration, we describe how the algorithm manages to maintain its validity at the end of the iteration. 1. edge[k+ l] is a forw ard edge. 2. edge[k+ l] is backward and goes beneath edge[k]. (Figure 2.4, left.) 3. edge[k+ l] is backward and goes above edge[k]. (Figure 2.4, right.) Case 1 does not require any processing yet at this point. In case 2, we have a snare situation like the linear case. Note that the beneath test is performed against edge[k] and not CurveSeg since the latter may just be a lid whereas it is the tangent of the former at edge [k] . sto p which defines beneath and above. Also note th at the convexity of edge [k] is irrelevant in performing the beneath test. The snare situation is “harder” to break than the linear case because we consider a backward edge to be snared regardless of its polar angle relative to edge[k] . sto p , a n g le (). As always, a backward edge will eventually be covered by a forw ard one and hence can be safely discarded. In case 3, edge[k+l] may or may not get snared into a pouch as in the linear case. If it does, however, it is not necessarily the most recent pouch into which edge[k+l] goes. Thus we skip subsequent backward edges instead, popping off all forw ard edges which were previously on the stack but are now completely shadowed by these backward edges. At this point the current edge is a forw ard one and its starting angle is greater than the starting angle of the topmost edge on CH no m atter which of the three cases we began with earlier. We compute the common supporting line segments of the current edge with the topmost edge on CH and keep popping off those edges which are completely shadowed by the common supports. When we eventually arrive at an edge which is not completely shadowed. After cutting it appropriately, we push the cover as well as the current edge onto the stack if they are not empty. The invariant is restored at this point. 2.2.3 C orrectness and C om plexity Arguments similar to those in the rectilinear case show that this algorithm correct and that the combinatorial complexity is O(N). The algebraic analysis earlier shows 13 th at the the combinatorial complexity is further increased by a factor of AD. Fi nally there is an extra factor incurred by the algebraic computations on the edges during preprocessing and in the main loop, of which finding the common support is the dominant one. In practice this last factor may arise in the form of numerical computations and therefore depends on the desired precision. Note that all the comparison operators for the angles can be unambiguously defined because each edge has a small span of polar angles. The subroutines not shown in the algorithm are purely functional and algebraic, and do not require more information than the explicitly specified parameters. In other words, they have no effects on and are not affected by the control flow of the program. 2.3 H ints and H euristics in th e D esign o f A lgorithm s for H igher D im ensions In comparison to the existing algorithms, we observe that our versions achieve a considerably simplified flow control mostly because the latter proceed around the origin. To put the comparison in a more common perspective, we can also say that the traditional algorithms proceed “around a point at infinity,” much as saying that an orthographic projection is in fact a perspective projection with a point at infinity as its perspective. Consequently we notice the following distinctive features in our algorithms: 1. The coordinate system has to be chosen in such a way th at the origin falls in the interior of the convex hull in order for the algorithms to work. It is, however, irrelevant whether the origin falls in the interior of the polygon itself. The traditional algorithms are not affected by the choice of the straight line with respect to which the upper hull and lower hull are defined. 2. Interesting (or rather, annoying) events occur at points on the curve segments where the tangent to the curves pass the origin. These are therefore places where we break the edges in order to ensure the monotonicity of the polar angles of points on each edge. W ith the traditional algorithms, the corre sponding events occur at the folds, or where the tangents are perpendicular to 14 the straight line separating the upper and lower hull. The idea of monotonicity with respect to a straight line is either implicit or stated in these approaches. 3. In our algorithms as well as some traditional algorithms the edges are also broken at the inflection points. We notice that these features are also present in one or another dual construct in mathematics. It is thus our goal in the next chapter to set out reviewing these concepts as they are related to our problem. Figure 2.2: The cases of the next edge to be examined. The dashed lines indicate the (possible) edges on the polygon; the dotted lines separate the cases. The origin is marked by an “x.” 16 EdgeStack CurvePolygonHu.il (Vector<Edge> & edge) { EdgeStack CH; CH. push (edge [stop = i = x.max] .coverO) ; while ((++i % = M ) != stop) {. i f (edge[i].backward()) i f (e d g e [i].beneath(edge[(i-l+N)%N])) fo r (; ed g e[i].sto p .a n g le () < CH.t o p ( ) . s to p .angle() & & edge[i] .s t a r t . angle() > CH.topO .stop.angle()+M_PI II edge[i].backward(); ++i % = N); else fo r (; edge[i].backward(); ++i % = N) while (edge[i] .stop.angle() < CH.topO .s t a r t .angleO ) CH.pop(); a sse rt(! edge[i].backward() & & edge[i] . s ta r t .angleO > CH.topO .s ta r t .angle() ) ; while ( !CH. isEmptyO) { cover = edge[i].support(CH.top()); i f (cover. s t a r t . angleO > CH.topO .s ta r t .angleO ) break; CH.popO ; > CH. push (CH.popO .cut . a f t e r (cover, s ta r t .angleO )) ; if (cover. s ta r t .angleO < cover.stop.angleO ) CH.push(cover) ; i f (edge[i] .stop.angle() > cover.stop.angleO ) CH. push (edge [i] . cut .before (cover. s to p . a n g le O )); } CH.popO ; retu rn CH; Figure 2.3: Convex hull of a curved polygon 17 Figure 2.4: Left: backward edges going beneath—ways to break the snare situation. Right: backward edges going above. C hapter 3 C lassical Forms o f D uality and M ath P relim inaries The term “dual” in mathematics has a wide variety of definitions in different con texts. In this chapter we review several of them upon which our algorithm is based, and observe the relationships among them. The foundation of our algorithm is es tablished by giving an alternative description of the topologically derived sets of the polar set. Finally we review terminologies from combinatorial topology and specialize towards the algorithmic needs which will arise in later chapters. 3.1 D u ality in P rojective G eom etry 3.1.1 N ota tio n s and T erm inology Conventions for Overloading Notations The d-dimensional projective space P d consists of equivalence classes of (d + 1)- tuples which have at least non-vanishing entry. Two such tuples are equivalent if one is a non-zero scalar multiple of the other. By bracketing the entries of a tuple and displaying it as a column vector, e.g. [£0, £1, • • •, £d)T-, we denote the equivalence class to which the tuple belongs, or the projective point it represents. To avoid introducing an overwhelming number of symbols in later discussions and proofs, we will agree upon the convention of overloading notations as follows. In general let A be a d + 1 by k m atrix enclosed in brackets where k < d + 1. The same symbol 19 A will also be used to denote the linear subspace spanned by its column vectors in contexts where sets are required. Embedding and Improper Points There are an infinite number of ways in which the Euclidean d-space E d may be embedded in P d so that many of the nice properties in P d carry over to E d. One natural choice is to identify each [£o,£i,- G P d for which £0 -= f 0 with (fo ’ to ’ ' ' ‘ ’ !o)T ^ ^ d' ^ P°int identified with some point in E d is called a proper point; and those which left unidentified with any Euclidean point are called improper points. In every embedding of E d in P d, the set of improper points is a (d — 1)- dimensional projective subspace of P d. This set is called the improper hyperplane with respect to the particular embedding. Hyperplanes The hyperplane H = {[£ o > £l > • • • 5 £d] : uo£o + + • • • + Udfd = 0} is uniquely determined by the dT 1 coefficients and also uniquely determines these co efficients up multiplication by a non-zero scalar. Thus H may be uniquely associated with the projective point [£0, £i, • • -, & * ]- We denote H by [u0, v> i,..., «d], the equiva lence class of the (d + l)-tuple displayed as a row vector. The point [no, u%,..., Ud]T and the hyperplane [no, wi,..., uj]T are projective duals to each other. Duality between Linear Bundles and Linear Subspaces Let B be an m by d + 1 m atrix where m < d + 1. The set of hyperplanes whose (d + l)-tuples are linear combinations of the row vectors of B form a linear bundle of dimension r — 1 where r is the rank of the matrix. Accordingly we will overload the symbol B for also representing this linear bundle in a context where a set is required. This notation will not conflict with that for linear subspaces as long as we work in a projective space of fixed dimension d. The intersection of all the hyperplanes belonging to B is the null space of B and is called the carrier or kernel of B (as a set). In general a linear bundle B m < (i+1 dualizes to the linear subspace and vice versa- 20 3.1.2 H ypersurfaces o f th e Second Order A hypersurface of the second order, or simply quadratic hyper surface, is the set of points d {[Co, 6 , : E = 0} i,j=0 for some fixed qt j, 0 < i.j < d with qtJ = q Jt. This hypersurface is uniquely de termined by the real symmetric m atrix whose entries are these qifs; and it can be shown that, in cases we will ever deal with, this hypersurface also uniquely deter mines such a real symmetric m atrix up to multiplication of all entries by a non-zero scalar. We will thus use Q = [qif\ the bracketed version of the m atrix (qK j) to denote its equivalence class, namely those which are its non-zero scalar multiples. The rank and eigenvectors are invariant within each equivalence class and are therefore used to define the rank and eigenvectors of the quadratic hypersurface itself. If the rank of the hypersurface is less than d -f 1, it is said to be degenerate; otherwise it is non-degenerate. Our abuse of notations is completed by using the square m atrix Q for also refer ring to the quadratic hypersurface itself in a set context. In the rest of this subsection Q will always be a real symmetric m atrix which refers to some hypersurface of the second order and its representing matrix. Tangent Hyperplanes and Polars If x is a point on Q , then (Q'x)T represents the tangent hyperplane of Q at x. If, on the other hand, x ^ Q, then all of the points of tangency on Q from x fall on the same hyperplane, namely (<5x)T. In either case, the hyperplane (Qx)T is called the polar of the point x, and the point x is called the pole of the hyperplane (Qx)T, with respect to Q. Dual If Q is non-degenerate, or equivalently the determinant of any representative of Q is non-zero, then Q~l exists and is also real symmetric. In this case Q~~l again represents a hypersurface of the second order, which we call the projective dual of 21 Q. The duals of the points on Q, which are themselves hyperplanes, are exactly the tangent hyperplanes of the projective dual of Q , and vice versa. 3.2 P edal C urves and P edal Surfaces 3.2.1 D efinition s and P rop erties Let 7(t) be a planar curve. The foot of the perpendicular from O to the tangent line to 7 (t) at 7 (t) is given by (7 (i) • N(t))N(t ) where N(t) is the unit normal vector at 7 (t). As t varies, this point traces out a curve called the pedal curve of 7 with respect to O. The derivative of the pedal curve of 7 with respect to the arc length param eter of 7 is S'(t) = -/c ((7 • T )N + (7 • N ) T ) where K(t) and T(t) are the curvature and unit tangent vector of 7 at t, respectively. The second factor vanishes exactly when 7 passes the origin. It is seen that the pedal curve is regular except at those points corresponding to the inflection points of the original curve. Similarly, the pedal surface of a surface in space is the loci of the foots of per pendiculars from O to the tangent planes to the original surface. A parabolic curve on the original surface corresponds to a cuspidal edge on the pedal surface. 3.2.2 D irect G eneralizations and a Variant Pedals in Higher Dimensions and in Degenerate Cases On the one hand, we may generalize the idea of pedal curve and surface to higher dimensions. Define the pedal hypersurface of a smooth hypersurface to be the loci of the foots of perpendiculars from O to the tangent hyperplanes to the original hypersurface. On the other hand we may also agree upon the convention th at in case the pedal fails to be a hypersurface at all, we shall still refer to this construct as the pedal 22 representation of the original hypersurface. In particular the pedal representation of a hyperplane H = {x : n • x = c] is a single point cn ped H = j-j-. \n \ Pedal Representation of a Set of Hyperplanes Let 7t be a set of hyperplanes. We shall use the same notation ped “ H to denote the set {ped H : H € H], knowing that the the context will make the meaning of this notation unambiguous. Inverted Pedal of Hyperplanes The pedal representation of a hyperplane is a bit awkward to work with. If c is non-zero, we may scale the equation so that c = 1. In this case we may further apply the inversion map: inv : E d - {0} — ► E d - {0} 1 to the perpendicular foot and obtain the vector n. We shall refer to this represen tation of H = {x : x • x = 1} by n as the inverted pedal representation of H. This form of duality is employed in linear programming and also in computing the convex hulls of finite point sets. 23 Consider E d as embedded in P d in the usual way. We observe th at [— 1, ni, n2, ..., rid], the projective dual of H, is related to the projective coordinates of the inverted pedal of H by multiplication by the matrix Jd+l 0 0 Id 3.3 Polar Sets Let K be a non-empty subset of E d. The polar set of K is K* = {x £ E d : x • y < 1 for all y £ K} . For example: • if K consists of a single point p, and p 7- O, then K* is the closed half-space {x : p x < 1}; • if K is a closed ball of radius r centered at O. then K* is the closed ball of radius - centered at O: T 1 • if K is the unit cube bounded by X\ = ±1, x 2 = ±1, and x 3 = ±1, then K* is the octahedron with vertices ( ± 1, 0, 0), (0, ± 1, 0), and (0, 0, ± 1). The following general properties of polar sets are useful in our approach to the convex hull problem. Let K ,K i ,K 2, and K a(a £ *4) be nonempty sets. Then: 1. K* is a closed convex set containing O. 2. If Ah C I <2 then A 2 * C K*. 3. (\Ja£AKot)* = Cag^ (K a)*. In words, the polar set of the union is the intersec tion of the polar sets. 4. If A > 0, then (AK)* = (l/\)K*. 24 5. If K is compact and convex, and contains O as an interior point, then K* is compact and contains O as an interior point. 6. If K is closed and convex, and contains O, then (K *)* = K. O b se rv a tio n 1 Let K C E d be non-empty. Then K* C (conv K)*. proof: Suppose y £ K*. Let x £ conv K be arbitrarily given. Then there exist k £ N, Xj, x2, ..., X fc £ K , and 0 < Ai, A2, ..., A * < 1 such that x = x A;X;. Thus k y • X = A iy • Xi i~ l k < i=1 = l. Hence y £ (conv K)*. □ On the other hand by applying (2) to the fact that K C conv K , we see that K* D conv K. Thus we conclude that O b se rv a tio n 2 Let K C E d be non-empty. Then K* = (conv K)* and K** = conv K . 3.4 Topologically D erived Sets o f th e Polar Set 3.4.1 P artition in g H yperplanes w .r.t. a Set Let K C E d be compact and non-empty. Let us consider the set of all hyperplanes Tt in E d. They can be partitioned into three equivalence classes with respect to K as follows. Let H = {x : /(x ) = 0} be a given hyperplane. 1. if /( x ) < 0 for all x £ K or /( x ) > 0 for all x € K, then H is said to strictly bound K ; 25 2. if / takes on a positive value for some x € K and takes on a negative value for some y 6 K , then H is said to cut K; 3. Otherwise the sign of / does not change but does become zero for some choice of x € K and it is said to support K. We will denote the set of strictly bounding hyperplanes, the set of cutting hyper planes, and the set of supporting hyperplanes of K by 7ib(-7Q, Tic(K ), and 7is(K), respectively. It is straightforward to verify that conv K induces exactly the same partitioning of hi as K does, i.e., H b(convA ) = Hh(K) H c(convK) = H C (I<) W8(convJ0 = HS{K) 3.4.2 T opological P rop erties o f th e P artition in g We now show that the pedal representations of the hyperplanes in each class is topologically derived from those of the other two. L e m m a 1 Let K C E d — {0 } be a set containing at least two points. Then ped (7i cK) is open, in the relative topology of E d — {O}. proof: To show that ped (H C K ) is open, consider an arbitrary H € H C (K). We will construct an open neighborhood of ped H which is entirely contained in ped (7i cK). We leave out the technical details giving the precise formulas in the construction of an open ball and instead give an intuitive and pictorial account (Figure 3.1). Let p ,q G K be such that p and q fall on different sides of H. Every hyperplane passing any point on the line segment p q C K cuts K. It is straightforward to verify that their perpendicular foots form the shaded area in the figure. The point ped H is in the interior of that area and hence in the interior of ped 7i cK . □ 26 L e m m a 2 Let K C E d be a closed set such that O € int(conv/F ). Then ped(HhK) Q e x t ped (7icK), in the relative topology of E d — {O}. proof: If K does not have any strictly bounding hyperplane, then neither can it have any supporting hyperplane. Hence H c is the entirety of E d — {O} and ped (HhK) is vacuously contained in the exterior of ped (H C K ). So suppose that HhK is non-empty and let H = {x : /(x ) = 0} € Hh{K) be arbitrarily given and its sign chosen in such a way that x < 0 for all x € K. By assumption f ( 0 ) < 0. Let p ^ K be a point such that f ( 0 ) < f { p) < 0. (Such p always exists; otherwise K wouldn’t be closed.) Also let q be a point such that / ( q) > 0. Repeating an argument similar to th at in the previous proof, we conclude that there is an open neighborhood of ped H which is entirely in the exterior of ped (HC K ). q L e m m a 3 Let K C E d be a closed set such that O € in t(c o n v if). Then ped(HSK) C 5ped (HC K). proof: If K does not have any supporting hyperplane then the statem ent is vacuously true. So suppose that H SK is non-empty and let H = {x : /( x ) = 0} £ H S(K ) be arbitrarily given and its sign chosen in such a way that x < 0 for all x £ K. By assumption f( 0 ) < 0. Figure 3.2 shows th at every open neighborhood of ped H must contain points in ped (H C K ) (those in the open line segment Oped H ), as well as points in the complement of ped(HcK) (those in the closed ray beginning at ped H and pointing away from O). □ P utting all these together, we conclude that if K C E d is closed with O in the interior of its convex hull, then ped (HC K) is open and ped (HhK) and ped (H SK ) are its exterior and boundary, respectively. This follows from the fact th at the three sets of hyperplanes form a partition of H\ while the three topologically derived sets form a partition of E d — {O}, and hence the containment relationships in the last two lemmas are really equalities. 27 Figure 3.1: Shaded area represent loci of the perpendicular foots of the straight lines (hyperplanes in E 2) which necessarily cut K. Figure 3.2: ped if is on the boundary of the pedals of the cutting hyperplanes of K. 28 3.4.3 T opologically D erived Sets L e m m a 4 LetK C E d be such thatO £ in tc o nvK. Then K*—{0} = invped((H b^)U (HsK)). proof: Let p £ Ed — {O}. Then p £ l< <£=> p • x < 1 for all x £ K definition H — {x : p • x = 1} bounds K since O £ i n t conv K < = $ ■ p = in v p e d ((HhK) U ('HSK )) inverted pedal □ Earlier we have seen that the pedal representations of the three sets of hyper planes induced by K form topologically derived sets of one another. This is still the case after the inversion map is applied to these sets since inversion is a homeomor- phism of Ed — {O} onto itself. Now this lemma shows th at the polar set of K is in fact the inverted image of (HhK) U (HSK). We thus conclude with the following T h e o re m 3.1 Let K C E d be closed with O £ i n t conv K. Then i n t (K — {O}) = in v p ed (HhK) d(K — {O}) = inv ped (HSK) ex t (K — {O}) — in v p ed (HC K) . 3.5 C urvilinear C om plexes Linear Complexes and Polyhedra Let Po,Pi,P2> • • • ,Pfc be a linearly independent set of points in E d, k < d. The collection of their convex combinations sk = { ^ f=o*»Pi : = 1 and ti > 29 0 for 0 = 0,1,..., k] defines a dosed k-simplex. Its relative interior is called an open k-simplex. When we use the term k-simplex by itself we shall mean open A;-simplex. The simplex formed by a proper subset of p 0,P n P 2> • • • >Pfc is called a boundary simplex of sk. An m-complex is a collection of disjoint simplexes whose boundary simplexes are also in the collection and in which there is at least one m-simplex. The k-skeleton of an m-complex K is the set of all ^-simplexes in K, 0 < k < m. The point- set union of the simplexes in a complex K is called a polyhedron and denoted \K\. Correspondingly, K is called a triangulation of \K\. Curvilinear Complexes A topological k-simplex is a space homeomorphic to a ^-simplex in a Euclidean space. Topological complexes, topological polyhedra, and triangulation and be defined in terms of their rectilinear counterparts in a natural way. We are interested in topological complexes which are also embedded in a Eu clidean space and whose simplexes are diffeomorphic, rather than just homeomor phic, to the corresponding rectilinear simplexes. We shall use the terms curvilinear simplexes, curvilinear complexes, and curvilinear polyhedra to refer to such objects and the related concepts. 30 C hapter 4 D uals o f Curved P olyhedra 4.1 O bject M odel Input objects we consider in this thesis will be curvilinear polyhedra with additional properties to be stated below. Being topological polyhedra, they are compact. On the other hand, they are not necessarily connected, nor are they necessarily m an ifolds, although each of their triangulating simplexes is indeed a simply connected manifold without boundary. For practical considerations, we make some additional amendments on our de sired curvilinear polyhedron K regarding how they are expected in the input into the algorithms. 1. All diffeomorphisms between the curvilinear simplexes and their defining recti linear counterparts are restrictions of some algebraic functions to the simplexes. 2. Instead of requiring that the input be a description of the triangulating sim plexes of the polyhedron, we will accept descriptions by faces. A k-face is the union of some adjacent fc-simplexes s^. sJ;, ■ • •, sk n of K plus some of their boundary simplexes such that • the entire union is homeomorphic to an open fe-ball, • the homeomorphisms of all faces participating in the union are the re strictions of the same algebraic function / , • this algebraic function has a full rank Jacobian (k ) in the entire union, and 31 • the union can be described as the intersection of algebraic halfspaces with /• 3. Each k-fa.ce has at least (k + 1) subfaces of dimension k —\. A subface of the k- face is a face of K composed of the boundary simplexes of some of sk, s ^ ,..., sk m which do not vanish as a glue. More specifically, each k-f&ce < p k is expected to be described as < p k = {x : /j(x ) = / 2(x) = • • • = / d- fc (x) = 0, 9i(x) < 0,fif2(x) < 0 ,.. .,gm(x) < 0} where m > k and where the Jacobian of / i , / 2, • • ■ , fd-k has a full rank (d — k) everywhere in ipk. If K satisfies the above requirements, we shall call it a curved polyhedron. For computational convenience, we also assume that one sample point on each face is given along with its algebraic description even when the algebraic description alone is sufficient to single out just one simply connected component. We will refer to the algebraic variety defined by /j , / 2, ..., fd-k as the host variety of tpk. A 0-face is also referred to as a vertex, a 1-face as an edge, and a (d — l)-face as a facet. 4.2 D uals o f Faces The polar set duality is directly related to the convex hull of a set K in that one obtains conv K by taking its polar set twice. Its definition, however, does not provide an operational recipe for actually writing a program to compute the convex hull of, say, an industrial part, which we model by a curved polyhedron. The pedal and projective dual of surfaces, on the other hand, are algebraic varieties whose equations can indeed be algorithmically computed. Yet we will have to answer the following questions before we can employ these ideas in solving the convex hull problem. • W hat is a useful definition for the dual of a surface patch as opposed to th at of an entire algebraic surface? 32 • W hat is a useful definition for the dual of a face of a lower dimension? • W hat is a useful definition for the dual of a curved polyhedron once we have definitions for its faces? Here the usefulness of a definition depends on how we will eventually be able to relate the “dual” of a curved polyhedron to its polar set and construct the latter from the former. 4.2.1 D uals o f V ertices and Facets If < p° = {p} is a O-dimensional face on a curved polyhedron K , we define the inverted pedal, or simply dual of < p°, to be < p° = {x : x - p = 1} This is exactly the dual hyperplane in the linear programming literature, and coin cides with the boundary hyperplane of the polar set of cp°. If p d~x is a (d — l)-dimensional face on a curved polyhedron K C E d described by < /-1 = {x : /(x ) = O,0i(x) < 0 ,g 2(x) < 0,.. .,gm(x) < 0} , we define its dual to be ( ,pd~l = {invpedHx : x £ < p} where E x is the tangent hyperplane of the host variety of < p at x. Note that Hx is well-defined because of the assumption about the rank of V /. Neglecting the bound ary constraints imposed by gds, this is just the inversion of the pedal hypersurface of the host variety of ip (to the extent that the host variety is a hypersurface). Figure 4.1 shows a cycloid and its dual. The dual of the cycloid itself (the open 1-face) is the curve with two asymptotes which begins at Ui, goes through two improper points, comes back from the third quadrant, and finally ends at u2. The two end points of the cycloid, being 0-faces themselves, dualize to the two straight 33 lines Vi and v 2. Note that they are tangent to the dual of the curve at the end points Ui and u 2, respectively. Figure 4.2 shows an s-shaped figure. The tangent to the curve at one of the end points Vi happen to pass the origin. In the dual space this is reflected by an asymptotic branch of the dual curve having approaching an improper point from only one side, and vy being its asymptote. 4.2.2 D uals o f Faces o f G eneral D im ension s In general, let tpk = {x : /i(x ) = / 2(x) = .. • = / rf_*(x) = 0, 9i (x) < 0, < / 2(x) < 0,... ,flr m(x) < 0}. where 0 < k < d — 1. Its dual is defined to be the set < p k = {invped H : H has a contact of order 1 or higher at x with ip for some x £ ^ } . This definition subsumes the facet case where k = d — 1. If we consider every hyperplane passing a vertex to have a contact of order 1 with the vertex, then the vertex case where k = 0 may also be subsumed. In any case we observe that a supporting hyperplane of a face (pk necessarily have a contact of order at least one with (pk because of the smoothness assumption for k = 1, 2, ... ,d — 1 and because of the convention for k = 0. 4.3 T he C ell T heorem Let K be a curved polyhedron. The inverted pedal, or simply dual of K is defined to be K = {(p : ip ( E K}. The point set union of K consists of the inversions of the pedals of those hyperplanes which have a contact of order 1 or higher with at least one of IV s faces. 34 Figure 4.1: A cycloid and its dual It can be decomposed into three smooth open cells—the two end points and the rest of the curve. The two straight lines are images of the vertices. The dual of the rest is the a curve with two asymptotes, which correspond to the two supporting lines from the origin to the curve. Figure 4.2: An S-shaped figure and its dual. 35 On the one hand, each of such hyperplanes has a non-empty intersection with K, namely a set containing at least the point of contact. They then therefore necessarily belong to 7icK or 7isK and can never strictly bound K. On the other hand, let us consider an arbitrarily given hyperplane H € H SK. Let x be a point of support of K by H and < p k be the fc-face which x belongs to. If k = 0 then by definition all hyperplanes containing x, including H , belong to K . If 0 < k < d then our object model ensures that < p has a fc-dimensional tangent space at x. Moreover, x is in the relative interior of p. Thus H must contain or equivalently H must have a contact of order at least 1 with p at x, or else it would be a cutting hyperplane. Therefore in this case H also belongs to K. Finally, k cannot be d since any common point of H and < p would be an interior point and cannot be a point of support. We have shown that: inv ped (H S JT) C K C inv ped (HsK) U inv ped (HC K ). Combined with the alternative description of the topologically derived set of K* (Section 3.4.2), we conclude that T h e o re m 4.1 Let K C E d be a curved polyhedron. Then the polar set of K is the closed d-cell containing O in its interior in the arrangement of all hypersurface patches in K . In particular, d(K*) coincides with the boundary of this special cell. 4.4 C om putations and G eneral P rop erties o f D ual In this section we give formulae, properties, and examples of duals, for later reference. We are mainly interested in computing the duals of the host algebraic varieties since the extent of the original faces may be directly modified to serve as the extent of the duals (see below). Thus in the following we will sometimes refer to the “dual of an algebraic variety < p,n which really means the host of the dual of a face hosted by the algebraic variety < p . Many of the properties follow from the definition. Some others are direct corollaries of their counterparts in projective geometry, linear programming, polar set, or pedal curves/surfaces, preserved by the transformations 36 involved. For example, the cusp-inflection duality is a property derived from the corresponding property of a pedal curve, and preserved by inversion, which is a homeomorphism from E d — {O} onto itself and is continuously differentiable up to any order. We will give proofs only in the remaining cases. 4.4.1 G eneral P rop erties • Common supporting hyperplanes of two faces dualize to the intersection of their duals. • The intersection of some finitely many hyperplanes dualizes to the affine hull of the duals of these hyperplanes. • Let < p be a ^-dimensional linear subspace of Ed whose pedal is n. Then its dual < p is the unique (d — k — l)-dimensional linear subspace whose pedal is inv n and which is orthogonal to the spanning vectors of < p . In particular, vertices dualize to hyperplane and in E 3 straight line edges dualize to straight lines. • In E 2, ordinary cusps and ordinary inflection points dualize to each other. (Figure 4.3) • In E 3, a parabolic curve on a surface dualizes to a cuspidal edge on the dual surface. (Figures 4.4 and 4.5, equation taken from [32].) • The envelope of a ^-parameter family of hyperplanes dualizes to a k-face. In particular, let < p be a general developable surface in E 3 and 7 its edge of regression. Then 7 is a developable surface and its edge of regression is < p . • In J53, a ruled surface dualizes to a ruled surface since tangents to a ruling form a family of planes sharing a straight line (namely the ruling), which dualize to a straight line. 4.4.2 Im p licit Faces 37 Figure 4.3: A cardioid and its dual. It need be decomposed into at least two smooth open cells—the cusp the the rest of the curve. The horizontal line is the dual of the cusp. The dual curve begins and ends at the inflection point, at which the horizontal line is the tangent. Figure 4.4: A dimpled surface—the zero set of (4x2 + 3t/2)2 — 4a:2 — 5y 2 + 4z 2 — 1. 38 Figure 4.5: Dual of the dimpled surface. The two “umbrellas” turned inside out correspond to tangent planes to points inside the two dimples. The apparent sharp cuspidal edges of the umbrellas are duals of two parabolic curves inside the dimples. 39 L e m m a 5 Let ip be a k-face whose host variety is defined by the set ofd—k equations /,-(x) = 0, * = 1, 2,... , d - k and let d,-(y) = * = 1, 2,... d — k. V /i(y) • y r p Let A(y) be the (d — k — 1) x d matrix whose i-th row vector is [dj(y) — dd_fc(y)] and let c8 (y) for i = 1,2,...,& +1 be a basis for its null space. Then a set of implicit equations in the variables x representing the host variety of (p is given by eliminating y from (x - d n_*(y)) • c,-(y) = 0, i = 1, 2,..., k + 1 fi{y) = 0, i = 1, 2, . . . , d - k. For example, let ip be the surface in E 3 (p = {(xt,X 2,x 3) : f{x\, X2, X3) = 0} . Then its dual is W ( * 1,S 2,S 3)______ . f , ) = f) Wf(Xi,X2,Xz) ' (xU X2,Xz) The extent of the face (p can be determined as follows. Let x be an arbitrarily given point known to lie on the host of the dual of the face ip. Then the dual hyperplane of x is tangent to the host of ip at some point y. Find all such tangent points and conclude that x € < p if and only if any such point of tangency lies within the extent defining ip on its host surface. In a later chapter we discuss in more details how this procedure can be carried out with less numerical problems. 40 4.4.3 P aram etric Faces Formulae for parametric faces can be similarly derived from the definitions. The extent of the the dual faces is especially easy in this case since the dual is also a param etric face. All of the k parameters in the primal A:-face are also parameters for the dual and retain their extent. The additional (d — k ) parameters of the dual extend over the entirety of R d~k. 1. Let 7 C E 2 be the edge {r(t) : t G 1} where I is an open interval. Then 7 is the param etric curve: f |r fr - ( r .r - )r - 1 \ |rj2|r'|2 — (r • r')2 ' tfc 2. Let < p be a parameterized surface in E 3 given by {x(u, v) : (u , v) G D} where D is open and simply connected in E 2. V = 1 i— XvX ,,X ^ : (m,u) 6 • 3. Let 7 be a parameterized curve in E 3 given by {x(t) : t G D} where / is an open interval. Then 7 is the developable surface f x(* ) |x(t) X X — - 7 T + sx(t) x x'(t) : t G D and 5 G R > . •Wl J 41 C hapter 5 Q uadratic H ypersurface P atches In the previous chapter we have seen that computing the dual of a hypersurface in general involves variable elimination. Such algebraic computations are expensive and often requires the use of exact arithmetic. Since many industrial parts in reality are modeled by low degree surfaces, we would like to investigate if there are less expensive alternatives in these special cases. Fortunately it is the case with quadratic surface patches bounded by planar conic sections, for which the dual constructions can be largely reduced to linear algebra. Although our main focus in this chapter is on quadratic hypersurfaces, the results we obtain are useful in the more general setting where the hosts of the faces are of higher degrees. In the last section we will be able to show that in a practical sense the severity of degeneracy of the duals is determined completely by investigating the partial derivatives of the original algebraic variety. 5.1 N on -D egen erate H osts and th e E xten t o f th e D ual Let H = {x G E d : n x = 1} be a hyperplane in E d where n = (ni, n2, ..., n < * ) G E d was observed to be the inversion of its pedal. We observed in an earlier chapter that if we write n = [l,n1,n2, ..., as a projective point in P d, then the projective dual rF of H is related to n by n' = Jd+in 42 where Jd+1 is the Id+i with the first diagonal entry replaced by — 1. Now consider a non-degenerate quadratic hypersurface {x € P d : x t Q x = 0} where Q = [ql3\, 0 < i,j, < d has non-vanishing determinant. It is straightforward to follow the proofs of the theorems in projective geometry and conclude that Q the inverted pedal of Q is the proper portion of {x £ P d : x TJd+iQ~ 1Jd+ix = 0}. In our notation, we will use JQ~lJ to also denote the dual surface in a set context. The extent of the dual face on the dual host surface can be determined implicitly as in the case for general dual surfaces. In the quadratic case, it can also be deter mined directly in the dual space. Let H be a hyperplane bounding the face < p whose host is Q. The portion of the boundary of C p as determined by H consists of the duals of the tangent hyperplanes taken along the cut. These tangent hyperplanes all intersect at one point, namely p the pole of H with respect to Q. Thus their duals are just the intersection of p with (p. 5.2 G eneralized P rojective D uality The projective duality x t Qx = 0 <=> n T$ -1n = 0 of quadratic hypersurfaces was derived based on the hyperplane-point duality: {x : n • x = 0} n T. When we replaced the above hyperplane-point duality by the pedal duality: {x : n x — 0} (n Jrf+1)T 43 we obtain the conclusion in the previous subsection. In general if we begin with a hyperplane-point duality defined to be {x : n x = 0} RnT where R is a non-singular real symmetric (d + 1) by (d + 1) m atrix, then we obtain the duality x r Qx = 0 4=^ nT(RQR ) 1n = 0. (5.1) in like manner. In this case a linear bundle Bmid+i dualizes to the linear subspace It can be interpreted as dualizing the points, hyperplanes, and non-degenerate quadratic hypersurfaces with respect to the fixed non-degenerate quadratic hyper surface R. The pedal dual is a special case when R = Jd+i, the unit hypersphere. A second interpretation is that we are still taking the ordinary projective dual, except th at we are carrying out all the computations in a new coordinate system in which the imaginary quadratic hypersurface For we can let x = Ry and see that the duality (5.1) taking place in the old co ordinate system translates to the “standard” dual in the new coordinate system expressed in y. 5.3 D egenerate H osts W hen Q has a vanishing determinant, the projective dual is undefined and the tangential dual need be computed directly in this case. takes on the new form d i,j=0 44 Let Q be a degenerate quadratic hypersurface of rank k < d + 1. Write Q as U AUT where j = 0,1, . . . , d the column vectors of U are eigenvectors of Q and where A is the diagonal m atrix composed of Q’s eigenvalues so arranged th at Xi ^ 0 for i — 0,1,..., k and A ^ = 0 for i — k + 1,..., d. Now consider the non degenerate quadratic hypersurface Q' — UA'UT where A' is diagonal with A '- = A , for i = 0,1,..., k and A '- = 1 for i = k + 1,..., d. Suppose first that an arbitrarily given hyperplane Fvecn is tangent to Q. Then nQxo = 0 for some X o € Q- Perform a coordinate transformation x = Uy and consider the new coordinate y 0 = f/Tx 0 of Xo- In the new coordinate system, the quadratic hypersurface takes on the form y TAy = 0 and the given tangent hyperplane is n il. It is straightforward to verify that the tangent hyperplane of A (that is, of Q in the old coordinate system) at each y whose first k + 1 coordinates are the same as y 0, is also n U. Necessarily the last d — k coordinates of nil vanish. This implies that nU is tangent to A' (or, to Q' in the old coordinate system) as well as passes e, for i — k -{ ■ 1^... ,d + 1 where e, has all but one non-zero entry at the ?-th position. Conversely suppose we are given a hyperplane n ' which is tangent to Q', and the last d — k of whose coordinates n'U in the new coordinate system vanish. Let X q be the point of tangency on Q' , i.e., n'Q'Xg = 0. Then it is not difficult to find a point X o such that n ' is tangent to Q at X o- Specifically let Xo be the vector whose coordinates UTx o in the new coordinate system are the same as UTx ' 0 except with the last d — k entries replaced by 0’s. Thus the tangent hyperplanes of A can be characterized as the those of Q' which also pass u*i^+i , ..., u*i ( i. The inverted pedal of Q is therefore the intersection of u ^ + i , . . . , u *td, and Q'. Since Q' is non-singular, Q' = Jd+ 1(Q,)~1Jd+1. The u * /s are linearly independent and dualize to u F j d+1. Hence the intersection of these hyperplanes is a fc-dimensional linear subspace tt C P d. Furthermore, 7 r is not tangent to Q'. because ..., u»> ( / ^ Q'. Thus Q, or equivalently the intersection of 7 r with Q', is a quadratic hypersurface of P k embedded in P d. In other words, the dual of a rank k quadratic hypersurface is (k — 2)-dimensional. 45 5.4 Lower D im ensional H osts Let p be the host of a (d — k — l)-dimensional face (1 < k < d) of a (degenerate or non-degenerate) quadratic hypersurface patch whose host is Q. Let d Hj = {[6 , 6 , • • •, C d] : Y = 0} j = 0, 1,..., k — 1 i=0 be a set of linearly independent hyperplanes whose intersection with Q defines p. Let N be the (d + 1) x k m atrix whose columns are the vectors n*,, which after transposing represents the linear bundle spanned by these hyperplanes according to our convention. Perform QR factorization on N and write it as n = where • The dimensions of U\, C /2, R, and the zero m atrix O are (d + 1) x k , (d + 1) x (d + 1 — k), k x k, and (d + 1 — k) x (d + 1 — k), respectively; • U = (U1U2) is orthogonal; • R is upper triangular and non-singular. W rite Q as Q = (c.c/2) ( ^ , where A of dimension k x k and C of dimension (d — k -f-1) x (d — k + 1) are real symmetric. Then it is easy to verify that the following is another characterization of p equivalent to the original one: {x : U^x = 0 and x TU2C U jx = 0}. Denote the degenerate quadratic hypersurface U2CUJ of rank at most d — k + 1 by Q'- 46 T heorem 5.1 Let p, U2, and C be defined as above. If C is non-singular, then < p = Jd+i,d+\U2C~1U f Jd+i,d+i is a quadratic hypersurface of rank d — k + 1. proof: We will prove a slightly different result in order to keep the notations m an ageable. In this proof only we modify the definition of tangential dual to be taken with respect to the imaginary quadratic hypersurface represented by the identity m atrix, namely the ordinary projective dual, instead of with respect to the unit sphere Jd+i,d+i- Thus let S = U2C~1U f be the quadratic hyper surface which we desire to show to be (p. Once we identify S with (p. we may apply the interpretations in a previous section by actually taking the dual with respect to Jd+i,d+1 and obtain the exact statem ent in the theorem. We first show that < p C S. Let x 0 € < p be arbitrarily chosen and fixed. The tangent hyperplanes of ip at x 0 form a linear bundle of dimension k whose carrier (kernel) is the intersection of the linear subspace N (or equiv alently U] since as matrices their transposes have the same null space) with the hyperplane yffh^CUf. Thus in the dual space these tangent hyperplanes form the linear subspace spanned by the column vectors of U\ and the vec tor L^CT/Jxt). For each arbitrary & -vector w and each arbitrary scalar a, let y = U\W + aU 2C U fx 0. It is easy to verify that yTU2C - 1U fy = 0 and hence conclude that (p C S. Conversely, suppose y0 € S. We will consider L = {ay0 + C/iV : a £ R and v € R k}- in the dual space all at once. It can be verified that L C S and is a linear bundle of dimension k in the primal space. Using the column vectors of U\ and U2 as a basis, it is easy to find the carrier of L in the primal space to be {U2z : YqUJz = 0}. 47 It can also be verified that Zo = C'-1C/Jy0 is a valid choice of z in the above. Furthermore the vector x 0 = U2Z o lies in Q'. The tangent hyperplanes to Q' at x 0 dualize to the span of the vector U2C U j x 0 and of the column vectors U\. Substituting the formula for Xo and Z o into the linear combinations, we see that this linear bundle of dimension k is L' = {aU2U^y0 + Uxw' : a' < E R and w 7 G Rk}}. We would like to show that L = I!. The elements of L for which a = 0 coincide with the elements of L' for which a' = 0. For the other elements in both sets we scale so that a = a' = 1. It now suffices to show that y 0 - [ / 2j72 Ty 0 G {Uiw : w G Rk}. But the span of U\ is just the null space of U2 and it is easily verified th at the above claim is indeed true. □ Finally we mention that if C is singular, we may diagonalize C as VA.VT where the first m diagonal entries vanish, 0 < m < d — k + 1. Let A7 be the same as A with the vanishing diagonal elements replaced by l ’s. If we write Q as ( t w 2) then it is not difficult to see that a combination of the arguments in the degenerate host and lower dimensional host leads to the conclusion that (p is the intersection of JU2V{k')-lVTU lJ with the hyperplanes represented by the first m rows of For d = 3 the only case in which k > 0 and m > 0 corresponds to straight lines and points. Hence we may simply use the linear duals instead and will not have to invoke this result in practice. 48 5.5 D egeneracies and Singularities In this section we study the singularities brought about by dualizing a smooth hy persurface of higher degrees. It turns out to be more convenient to work in the Euclidean space than in the projective space. We first give a formula for computing the dual after a translation of coordinate systems. The duals are to be taken with respect to the unit sphere in the original coordinate system, however. Let p G E d be the origin of the new coordinate system expressed in the old coordinates. Then the duality d {(x 1,x 2, . . . , x d)T : J 2 nixi = 1} (ni,n2, . . . , n d)T in old coordinates translates to d {(X 1, X 2,. . . , X d)T : Y h ni{X i +Pi) = 1} (n l - P l >»2 - P 2, ■•■,«<* -P df ) i=I in new coordinates. Now suppose that we are given H = { p G ,X 2:... , X d)T : i niXi = s} in the new coordinate system to begin with, assuming th at H does not contain the old origin. Let A = —=£------ • Then the dual is An — p in the new S+ L i = l n <Pi coordinate system. Given an algebraic hypersurface < p and a regular point p in 9 2 , we now ask un der what conditions does the point dualize to a singular point, and how bad the singularity is when it does so. More precisely we would like to know the rank of the Jacobian m atrix of the dual with respect to any regular param eterization of < p nearby the given point p. To simplify the computations, we will show the derivation for d = 3, although the idea is the same in higher dimensions. First we rotate the coordinate axes so that the tangent plane at p becomes parallel to the x ^ -p la n e . Obviously the geometry of neither nor the dual is not affected and thus without loss of generality we will assume th at this is the coordinate system initially given to us. Next we translate the origin to p. Nearby the new origin (the given point), the surface can be written as (xi, x2, /(xi, x2))T 49 for some (not necessarily algebraic) function / . The tangent plane corresponding to the parameters {xi,x2) is {(S/i, 2/2,2/s) : {y\ - x i) /i + (j/2 - x 2) / 2 = 2/3 ~ x3}, where /, = are evaluated at the corresponding parameters. We use the formula previously derived to compute its dual g(xi,x2) and find its partial derivative with respect to x\ to be (.f i l C - f 1C1, f 12 C - f 2C1,C 1)T C 2 where C = xtfi + x 2f 2 - f + Pifi + P2/2 ~ P 3 and C\ = -£^C. We note th at /i|( 0,o ) = 1 (0,0) = 0. Thus the above partial derivative evaluates to (— / i l p 3, —fx2ps,pifiL + P2f i ‘ 2)T/C 2 at the corresponding point. Similarly j k g ( x " x >\ Note th at the third row of the Jacobian m atrix 1 ~ / l l p 3 ~ /l2 p 3 ^ - / l 2ps - / 22p3 ^ P l / l l + P2/l2 P l /l2 + P 2 /2 2 y is always a linear combination of the first two (recall that P3 ^ 0) and hence its rank is exactly that of the upper square matrix, which is completely determined by the partial derivatives of up to the second order. It is well known that the pedal of a surface has a cuspidal edge corresponding to the parabolic curves on the original surface. It would have been easy to reach the same conclusion for the inverted pedal by noting that the inversion map is continuously differentiable up to any order and hence preserves the dimension of such singularity. However, the result for the quadratic hypersurfaces together with a generalization of the above calculations enable us to conclude more generally th at the dimension of the dual at a point having k non-Morse terms is d — k — 1. Note th at we do not claim to have categorized different types of singularities. In terms | ° , ° — ^ ( - / i 2 p 3 , - / 2 2 p 3 , P i / i 2 + p 2 / 2 2 ) T . 50 of the computation of the convex hull, the above observation is sufficient for us to specify the algebraic preprocessing of the surfaces for the algorithm. For hosts of lower dimensions (less than d — 1), we observe that eliminating d variables from a system of (d — k) + (k + 1) equations in general leaves us with a single equation. We thus expect that the inverted pedals of most cells, regardless of their dimensions, to be (d — 1)-dimensional. The only degenerate case in E 3 is a straight line edge, which dualizes to a straight line edge. We will therefore not discuss the details of lower dimensional hosts since our focus is on 3-D algorithms and since such edges can be easily detected and discarded. 51 C hapter 6 C onvex H ull A lgorithm s In this chapter we present an asymptotically efficient algorithm for constructing the convex hull of a set of algebraic surface patches and other algorithms directly re lated to convex hulls. We recall that in the 3-D rectilinear case, the presence of the adjacency information on the simple polyhedron does not admit a more efficient algorithm to be developed than a simple one without taking advantage of it. Simi larly, our algorithm will accept as input an unordered set of surface patches (up to their boundaries) and therefore does not rely on the connectedness or homogeneity of the faces. It turns out that it is in fact nearly optimal. The algorithm works in three major steps. After preprocessing, the surface patches are first either discarded (if they are hyperbolic, parabolic, or elliptic but concave) or dualized. Equally im portant, and quite often probably even more so, is to dualize the boundary and dangling curves and the vertices. In the second step we begin with an arrangement of surface patches in the dual space. A number of cell construction algorithms may be chosen from and applied to create the convex cell th at contains the origin in its interior. For the asymptotic analysis we choose a reduction to the construction of the lower envelope. Finally the boundary of the convex hull need to be recovered efficiently from that of the polar set. 52 6.1 D ualizing th e Input 6.1.1 C hoosing th e C oordinate S ystem Application of Theorem 4.1 requires that the convex hull of the object K contain the origin in its interior. To this end, we may arbitrarily pick four or more points which lie on the input surface patches and verify that they are linearly independent. Their average is then in the interior of the convex hull and may be used as the origin of the new coordinate system. It may happen that we can not find any quadru ples of linearly independent points because all of the input surface patches lie on a linear subspace In this case we perform an appropriate coordinate transformation which brings the affine hull of K to one of the coordinate planes and invoke the corresponding algorithm in a lower dimension. In practice we would like not only the aforementioned property to be true, but also th at the origin is as far from the boundary of the convex hull as possible. It is easy to see that if a sphere of radius r > 0 centered at the origin is entirely contained in the convex hull of K , then K* is entirely contained in a sphere of radius ^ centered at the origin in the dual space. Thus we may perform the above choice of quadruples of points several times and compute the ratio of the volume to the area of the tetrahedron they define. This ratio is proportional to the radius of the inscribing sphere of the tetrahedron. We choose the quadruple which gives the largest ratio and use the center of its inscribing circle as the new origin. Let us denote this chosen tetrahedron by To and from now on we will refer to the coordinates and equations of the input surface patches with respect to this new coordinate system. 6.1.2 B reaking th e Faces at Singularities Most of the subsequent algorithms we will invoke assume that they are given “nice” inputs where each face is smooth. We recall that we make the same assumption about the input to the convex hull algorithm. However, in dualizing the faces we may introduce singularities as shown in the previous chapter. We are thus responsible for breaking the dual surfaces at these singularities in order for subsequent algorithms to work. 53 It is relatively straightforward to break each 2-face ip along the curve of support 7 C ip from the origin. This curve corresponds to an improper curve in the dual space. We also have to process the 1-faces in the primal space similarly by breaking them at points at which the tangent pass the origin (which do not exist if we make the general position assumption). Such singularities correspond to the fold, the “least singular” of all singularities. The study of more interesting singularities is a m ajor subject in itself. [7, 18] Fortunately we don’t have to deal with it in its entirety in order to describe the algorithm. Specifically, we have to break each dual surface at singularities into subfaces of lower dimensions, each of which is guaranteed to be regular for its dimension, but we need not worry about how these pieces fit together topologically except up to the subface relationship. In Section 5.5 we see that a point on the dual surface is regular as long as its corresponding point in the primal space has a non-degenerate second fundamental form. More generally, the dimension of the dual surface nearby a point p is at least (and usually exactly) as high as the number of non-vanishing eigenvalues of the second fundamental form at the corresponding point p in the primal space. Here we are using the term “dimension” in a computational context where we require a param eterization of the surface be explicitly established such th at its Jacobian m atrix has the rank equal to the claimed dimension. Thus we introduce artificial vertices in the dual space which correspond in the primal space to flat points [19] where the second fundamental form vanishes com pletely. We further introduce artificial edges in the dual space which correspond in the primal space to the parabolic curves, where the determ inant of the second fundamental form vanishes. In general the original surface is now divided into two or more regions each of which is either hyperbolic or elliptic. We may discard all hyperbolic regions since their contribution to the convex hull, if any, will be mani fested by their boundary curves. If, on the other hand, the original surface happens to be a developable surface, then the entire surface is parabolic and can be discarded for the same reason as hyperbolic regions in ordinary surfaces. Note that we need not introduce these extra lower dimensional faces in the primal space. Specifically, a flat point in the primal space does not contribute more than one tangent plane, if at all, to the convex hull. It does not create a plane in the dual space but is rather simply used to delimit the artificial edges in the dual space. 54 Similarly, a parabolic curve in the primal space contributes to the convex hull no more than the envelope of the tangent planes to the host surface taken along that curve, which forms a developable surface. Thus its dual is not a developable surface but merely a portion of a space curve which delimits the elliptic or hyperbolic regions. We finally observe that one type of singularity is absent in the dual space, namely the folds with respect to the origin. More precisely, a dual surface patch p cannot contain any point at which the tangent passes the origin, for such tangents would correspond to improper points in the primal space. 6.1.3 B ounding th e W ork Space in th e D u al Space Unlike the polar set dual, the compact object K in the primal space may give rise to unbounded portions in the (inverted pedal) dual space even though we have chosen our coordinate system carefully. We will thus have to further cut the dual surfaces and construct a reasonably small work space. For each coordinate axis, for example the s-axis, find its intersection with the boundary of the tetrahedron T0. Clearly the two intersections x = x mjn and x = x max must be such that x mjn < — r and x max > r where r is the radius of the inscribing circle of To- These six intersection points define an octahedron entirely contained in the convex hull. The dual of this octahedron is a rectangular box B whose faces are parallel to the coordinate axes and which contain the polar set. We thus clip each dual surface by this solid. 6.2 C onstructing th e Polar Set Cells in an arrangement of hypersurfaces can be constructed by any of a variety of well known algorithms. For example, a cylindrical decomposition would suffice for this construction. (This approach is similar to the one taken by Kriegman in [32].) However, more specialized and more efficient algorithms are also available for this particular problem [23]. As we are only interested in one particular cell (which is also convex), we may obtain an asymptotically even more efficient algorithm by way of reduction to the computation of the lower envelope, a problem studied in [46]. This randomized 55 algorithm takes as input E a set of N algebraic surface patches of bounded degree without folds, each with a bounded number of similarly constrained subfaces. It then produces the minimization map. This is basically a planar graph with additional geometric information describing the lower envelope of E, or the set of points right above the plane. The algorithm runs in expected time 0(iV2+£) for all e > 0. Let e be the unit vector pointing in the direction of positive 2-axis Consider the following function defined on the halfspace z > 0: fC : x i — ► a _1x + (a — l)e where a = x e . Its effect is to take vectors along the same direction to vectors ending in the same line parallel to e, with their order on the original direction preserved. The intuitive interpretation above also serves to show that if O is in the kernel of a star-shaped cell C in an arrangement of surfaces (which it is in our case), then half of dC fl C is mapped to the lower envelope of the new arrangement. We note that /C and its inverse )C~l : y i—► a y + ae — a 2e where a = y • e, are both algebraic. The transformation produces a valid input for [46]. Specifically, • Each transformed face obviously has a bounded number of subfaces as does the original face. • The transformation preserves the smoothness of the faces and increases the degrees of the host surfaces by a factor of at most two. • Each transformed face < p is fold-free. For let the host < p be given as f ( x i , x2, x 3) = 0 and let (yi, 1 /2, 2 /3) = Then the normal to the transformed face can be found by chain rule to be (y3f Xl, y3f X 2, yxf Xl + y2f X2 +fx3)- It is parallel to the xj/-plane (namely it is a fold) precisely when the tangent to (p passes the origin, which is impossible according to our earlier observation. Constructing K* the polar set of the curved polyhedron is now relatively straight forward. We first need to bound the work space. The transformation K. introduces 56 undesirable poles “along the horizon”—a point with nearly zero z-coordinate is trans formed to a point with a very large z-coordinate. For each face of the box B bound ing the polar set, a pyramid is formed with the face as the base and O as the top. The portion of the arrangement of the dual surface patches within the pyramid is then transformed with respect to the corresponding coordinate plane. It is easy to see th at the transformation of this portion is bounded and the lower envelope can therefore be computed and transformed back. Finally the six patches of dconv K corresponding to the six pyramids are sewn together. As a by-product of this reduction, we have shown that T h e o re m 6.1 The (structural) complexity of the convex hull of a curved polyhedron with N faces is 0 ( N 2+e) fo r any fixed e > 0 provided that the number of subfaces for each face is bounded by a constant. 6.3 R ecovering th e C onvex H ull Regardless of the algorithm we choose to construct the polar set, we can recover the convex hull from the polar set in tim e linear in the number of features on the polar set. For simplicity, we shall assume that the cell construction algorithm has produced a Doubly Connected Edge List of [39] for representing the boundary of a convex curved polyhedron embedded in E 3 as a connected planar graph G. Most existing general-purpose representations are equivalent to this simple structure as far as planar graphs are concerned, including the one generated by the reduction step described above. We could apply the same procedure to K* which was applied to K in order to find Jf**, namely the convex hull. However, we really need not compute the inverted pedal of the 2-faces of K* since they are contained in the input 2-faces. Better yet, we do not even construct the arrangement of the inverted pedal of the faces of K* since we are already dealing with a convex object. Instead, we just modify G to obtain the adjacency information of K** and compute the geometric information as follows. First a new planar graph G' is constructed whose vertices correspond to the vertices, edges, and regions of G. The vertices of G' represent the 2-faces of c o n v /l, the dual of the polar set. An edge is drawn between two vertices of G' if the two 57 corresponding entities in G are incident. These represent the 1-faces of conv K . Clearly G' is planar and the regions of G' represent the 0-faces of conv K . Denote by v the vertex in G that corresponds to v' 6 G '. To recover faces of conv K , we observe that a vertex v' in G' represents 1. part of an input surface patch in the primal space if v is a, region in G. Its geometric information is already present in the input. 2. a new developable surface patch which bridge two original surfaces if v is an edge in G. Its geometric information is computed by the formulaefor inverted pedals. 3. a triple common supporting plane if v is a vertex in G. Its equation is easily found by computing the dual of the corresponding vertex. To recover edges and vertices, we observe the following cases according to what kind of incidence in G an edge e' in G' corresponds to. 1. Vertex-edge adjacencies. In this case, e' represents the ruling in conv K that marks the transition from a plane w to a developable surface a. Its geometric information is easily computed as the dual of the tangent straight line of a at 7 T . 2. Edge-region adjacencies. Here e' represents a curve in conv K that marks the transition from a developable surface a to an input face p. This curve can be considered the intersection of er and < p , and hence the dual of the envelope of the tangent planes of (p along a (that is, their common tangent planes in the dual space). In particular if a can be parameterized then so can this new edge. 3. Vertex-region adjacencies. Now e1 represents a vertex in conv K th at marks the transition point from a plane to an original surface. Its coordinates can be found by dualizing the tangent plane of the region at the vertex. The size (number of vertices) of G' is at most six times th at of G since both are planar. Computing the geometric information in each of the above case contributes only a constant factor determined by the maximum degree of the polynomials. Thus 58 t> l el Figure 6.1: The union of two ellipsoids and a truncated cone. / 4 (planar curve) v\ / 2 (planar curve) / 3 (space curve) Figure 6.2: The polar set of the object, recovered as the cell containing the origin in the arrangement of dual surfaces. 59 t> l ei h Figure 6.3: The convex hull of the original object, recovered from the dual of the polar dual (previous figure). Note that f\ is planar since it is dual to the triple intersection (vertex) in the previous figure. Other new surfaces ( / 2, / 3, and f^) are developable, and arise from the corresponding edges of the polar set. 60 the combinatorial complexity of the convex hull and also that of computing conv K from K* is the same as the complexity of dK*, up to a constant factor. Figure 6.1 shows two ellipsoids and part of a cone in space. The contribution of the cone to the convex hull is limited to its vertex. The ellipsoids have implicit equations e\ : 29a:2 + 36t/2 + 180*2 — 24 xy — 48a; + 144 y — 36 = 0 e2 : 225a;2 + 612y2 + 388*2 - 768yz - 2448y + 1536* + 1548 = 0 . Their inverted pedals are computed to be: ei : 36x2 + 9 y2 + 5 z2 + 24 xy — 20 y — 5 = 0 e2 : 100a;2 - 3y2 + 153*2 + 192yz + 100y - 25 = 0. Figure 6.2 shows the polar set and Figure 6.3 its convex hull. In summary, we have shown that T h e o re m 6.2 Let K be a curved polyhedron with N algebraic faces of bounded de gree each with a bounded number of subfaces. Then conv K can be computed in randomized expected time 0 ( N 2+t) for any e > 0, where the constant depends on the maximum degree of the polynomials in the input as well as on the maximum number of subfaces each face may have. proof: Each step in the series of computations involve processing faces individually (up to their subfaces). Thus each original face is processed at most only a constant number of times. The transformations and preprocessing create no more new faces than a constant number of times the input faces, as determined by the maximum degree of the input polynomials. Hence the complexity of the lower envelop algorithm dominates and we have the desired conclusion. q It is known th at the complexity of the convex hull of N spheres can have com plexity Ll(N2) [6]. Generally one expects a worst case complexity of fl(N2) for algebraic surfaces of bounded degrees. It is also conjectured th at the complexity 61 is actually Q (N \q(N)) where q is a constant depending on the degree of the sur face and Aq(N) — 0 (N log* N) is the length of the maximum Davenport-Schinzel sequence of order q on N symbols. See for example [1] for a short summary and further references. 6.4 G eneral C urved P olygons in E 2 In Chapter 2 we study the convex hull algorithm which takes advantage of the fact th at the input is a simple curved polygon. If the input is a general curved polygon which consists of dangling edges and disconnected components, we may no longer have a linear time algorithm. The special case in which the curved polygon consists of a set of N discrete points serves to give a lower bound of fl(N) for the convex hull problem. In this section we show that this lower bound can be achieved. Let K be a curved polygon consisting of N isolated vertices and/or algebraic edges of degree at most D. As usual, the coordinate system is so chosen th at the origin lies in the interior of the convex hull. As in the 3-D case, we first break the edges at the turning points, namely the folds as seen from the origin. We also break them at the inflection points, which are the 2-D analogs of the parabolic or flat points in 3-D. Finally, we break long edges into pieces with angular spans no more than 7T. The same arguments used in the simple curved polygon case applies here as well and we conclude th at the number of curve segments does not increase by more than a factor of 4D. For each piece or each isolated vertex, we compute its dual. An isolated vertex dualizes to a halfspace of interest to us, namely the one which contains the origin. A curve segment consists of two vertices and an open edge. Its dual as a whole again defines a halfspace delimited no longer by a straight line but rather by one piece of curve segment and at most two rays (fewer if one or both of the vertices in the primal space are turning points). Now th at we have 0 (D N ) halfspaces, we may perform a standard divide-and- conquer technique to find their intersection. The choice of the origin ensures that during the merging process we always have a pair of convex regions with nonempty intersection. Their intersection can be found in time linear in the number of the curve edges bounding the two convex regions. 62 We bound the complexity of the algorithm as follows. We first observe th at the convex hull of a set of k non-intersecting simple curve segments in the plane has a structural complexity of A 2(fc). In other words, there do not exist a pair of the input curve segments a and b such that their portions alternately appear in the convex hull by more than three times (counting both), for otherwise they would intersect each other. Thus the convex hull of k curve segments has a structural complexity of at most A 2(&) = 2k — 1, whether they be algebraic and if so regardless of their degrees. It follows that our algorithm runs in time 0 ( N log N ) where the constant depends on D. To summarize, we have T h e o re m 6.3 The convex hull o f a set of N smooth algebraic curve segments and/or isolated vertices consists of at most 2N — 1 smooth pieces and can be constructed in time 0(iVlog7V) with an overhead determined by the maximum degree of the curves. 6.5 A voiding H ull C onstruction In some applications the convex hull of a curved polyhedron K may not always have to be explicitly constructed. For example one may wish to make a small number of queries as to whether an arbitrarily given point in space lies in the convex hull of a curved polyhedron but otherwise not interested in the convex hull itself. In this case a much simpler algorithm can be used without actually constructing the convex hull. A point p is • in the interior of conv K iff the dual plane p strictly bounds K*; • in the exterior of conv K iff the dual plane p cuts K*\ and • on the boundary of conv K iff the dual plane p supports K*. The dual problem is therefore only slightly more difficult than the 2-D convex hull finding problem. Suppose the input faces have been pre-processed as in the complete algorithm for constructing the convex hull. In the dual space we intersect p with the closure of (p for each face ip. We discard all surface patches which do not intersect p. If none of the dual surface patches remain, we conclude that p strictly bounds K*. 63 Now suppose that some dual surface patches do intersect p. If p supports any of the dual surface patches, we mark p as being potentially on the boundary. In any case, we actually construct the convex “half plane” determined by each intersection. The intersection itself is generally a curve segment with two ends. We always make it into a closed or infinite path by adding two rays tangent to the segment at its end points. Now we may invoke an algorithm similar to the intersection algorithm used in the construction of the general (non-simple) 2-D convex hull, except that we don’t have a fixed point on p which is guaranteed to lie in the intersection. Nevertheless in the merging process, we can take linear time to decide whether an arbitrarily chosen point O' lies within both convex regions. If so, we just proceed to merge their boundaries as before. Otherwise we break the two circular lists of boundary curve segments into at most four sublists, each of which has monotonically increasing or decreasing polar angles with respect to O'. The intersection of the convex regions can still be found in linear time. If the intersection of any pairs turns out to be em pty during this process, then p strictly bounds K*. Now suppose the intersection exists. The algorithm then returns that p lies on the boundary of conv K if it was previously marked as being potentially so; or that p lies outside conv K otherwise. In the latter case where p lies outside the convex hull, it is straightforward to actually construct the supporting cone from p to K . T h e o re m 6.4 Let K be a curved polyhedron with N faces. A point p in space can be classified against conv K in time 0 ( N log N ) without explicit construction of conv K . I f p lies outside the convex hull, the supporting cone from p to conv K has structural complexity 0 (N ) and can be constructed in the same amount of time as the classification problem. 6.6 A P ractical A pproach 6.6.1 R ation ale In many applications limited precision of the output is tolerable as long as the error is predictable, whereas an exact algorithm which is prohibitively slow is to be avoided. The reduction of the convex hull problem to the lower envelope problem—although 64 serves very well its purpose of giving us a satisfactory asymptotic bound—may incur too much algebraic and coding overhead. In this section we consider a numerical approximation for which • an intuitive measure of error can be defined; • the error can be estimated and controlled; • the output can be further processed to generate the “exact” solution, subject to missing “small” features as determined by the error measure; • the coding effort is minimal. The algorithm works as follows. In the dual space we choose a set of nicely spaced directions, for example, ones corresponding to the grid points on the unit sphere with “equally spaced” (in the usual sense) longitudes and latitudes. For each direction we shoot a ray. This ray is intersected with each dual surface to find the intersection point nearest the origin. Thus we have a mesh of sample supporting planes along with the knowledge of the faces they support. This mesh can then be further processed depending on the application which requires the construction of the convex hull. 6.6.2 P recision and Error M easures We choose the number of rays penetrating a unit area on the unit sphere as the measure of precision. Let us explain the intuition behind this choice by considering the two (disconnected) curved polyhedra in Figures 6.4 and 6.5, respectively. In the first case one reasonably asks the algorithm always to report the cone-like area as a new face even if the middle ellipsoid shrinks to the size of the needle. In the second case one should expect a program running on a machine with limited numerical precision to incorrectly report or leave out the “flat island” formed by the middle ellipsoid inside the new surface, even though the area of the reported or missing island may be rather large. Using our proposed measure of precision, the algorithm on input of the first example is not entitled to leave out the new surface, or else it would mean that tangent planes placed at such obvious different angles as those on the opposite sides of the cone are left out or identified. On the other hand, 65 Figure 6.4: A “big” feature that should not be missed in terms of convex hull construction. ’precs.pri’ Figure 6.5: A “small” feature that could be missed if the the precision of the ap proximation is low. 66 the second example does suffer from having nearly coincident planes and should be allowed to miss the big feature in the primal space. Such choice of precision measure can also be interpreted as weighing the sampling density (number of rays per unit area on faces in the primal space) by the Gaussian curvature of the surface since the latter is the limiting ratio between the area of the Gaussian sphere and that of the surface area in the primal space. In other words it realizes the intuition that we need more samples around surface areas which are badly curved, and don’t wish to waste bandwidth on nearly flat faces. The sampling criterion is appropriate for visualization purposes. Rendering algo rithm s based on the normals of the sampled facets usually prefer a “roughly constant change of normal direction” to a “roughly constant size of rendering area.” This technique may be applied to curved polyhedron consisting of higher degree surface patches such as the dimpled surface in Figure 4.4, whose duals are too complicated to derive symbolically. This measure has a geometric meaning in the dual space as well, up to a factor bounded by two positive constants. Let R c and R{ be the radii of the smallest cir cumscribing sphere and largest inscribing sphere centered at the origin, respectively, of jK*. Then the sampling density (number of rays per unit area on the faces in the dual space) on the boundary of the polar set is easily shown to be proportional to th at on the Gaussian sphere to within a factor between and -jg. c i Figure 6.6 shows the precise correspondence between the error in the dual space and th at in the primal space. Consider a small ball of radius e in the dual space representing the range of error. Since inverted pedal is symmetric with respect to the origin, we may assume without loss of generality that the center of the ball lies on one of the positive coordinate axes, say the t/-axis. We further simplify the discussion by looking at only a cross section of the ball where the cutting plane passes the origin, say the xy-plane. Let the center of the boundary circle be (0, h) and assume th at e <. Then the m atrix representation of the circle is B ( h2 - e 2 0 - k \ 0 1 0 \ - h 0 1 / 67 Its dual B \ 1 0 0 - e 2 - h 0 -k h2 — e2 is a hyperbola centered at (0, Points inside B may be thought of as points on smaller balls concentric with B. Thus they dualize to tangent lines to similar- looking hyperbolas whose two sheets are closer than and sit in-between those of B. These lines all fall in the concave area in-between the two sheets of B (Figure 6.6). Conversely, every line falling in that concave area dualizes to a point inside B . Even though the error ball in the dual space dualizes to an unbounded volume except by the enclosing ball in the primal space, the orientation error does indeed approach zero as the error ball shrinks in size. The latter is determined by the angles between the two asymptotes of the hyperbola, namely 2 sin- 1(c//j), and becomes proportional to the radius of the error ball as e — > • 0. 6.6.3 O ther P ractical C onsiderations The output mesh generated by the described algorithm may be further processed to produce a planar graph describing the boundary of the exact convex hull (subject to missing features as explained earlier). We begin by explicitly constructing the planar graph represented by the mesh. For each adjacent pair of vertices 14, V % in the graph which belong to two different dual surfaces, we add an interm ediate vertex along the grid line V{V 2 by computing the intersection of the two host surfaces with the plane containing Vi, V2, and O. Then for each grid for which any boundary grid line has changed, we use a brute-force procedure to add exact edges which correspond to the intersection of the actual dual surfaces, adding new vertices as necessary for triple intersections. Finally we remove redundant edges connecting adjacent vertices which fall in the same dual surface, removing the vertices as well when they become isolated. We may not want to carry out the last step to remove all redundant edges and vertices we could, since certain algebraic operations are simplified when the spans of the surfaces are kept relatively small. An analogous 68 situation is demonstrated in Chapter 2 where we restrict each boundary curve to have an angular span of it or less. In case the hull has to be constructed exactly without missing features, existing algorithms for converting CSG to boundary representations such as [38, 42] may be used in the dual space in place of the near optimal algorithm for convex cell construction. Although the asymptotic worst case behavior is worse, the similarity of these algorithms to their primal counterparts may argue in favor of their use. In CAD and CSG applications, being able to represent the original and de rived surfaces and curves analytically versus numerically means a significant saving in storage. One would even like to have a piecewise parameterized representation available for rendering and other purposes. Fortunately we can do both quite effi ciently and effectively with the duals and convex hulls of curved polyhedron consist ing of quadratic surface patches cut by planes. Using the fact that all dual surface patches are quadratic, it is straightforward to show that the intersection curves in the dual space can be parameterized by solving degree 4 equations. By the re sults in Section 4.4.3, the common supporting developable surfaces can therefore be parameterized as well. 69 Figure 6.6: An e-error circle and its dual hyperbola. 70 C hapter 7 Stable P oses and C apture R egions Given a solid object K , one may wish to find what are the possible orientations at which it can rest on a flat supporting plane. One application may be in manufactur ing environments where an algorithm for producing oblivious grasping plans need to take as input the possible stable poses of a part moving on a conveyor belt. [40] Another application is in model based vision algorithms where the object to be rec ognized rests on a horizontal flat surface such as the floor. More detailed discussions on the latter application can be found in [30]. One may ask a more general question: for each stable pose, what are the initial orientations of the part K starting from which K will eventually converge to that stable pose, assuming that the dynamics are dissipative? We apply inverted pedal to these problems and propose an alternative practical approach which in the quadratic case is far more efficient than the algorithm used in the first implementation (to the best of our knowledge) for solving these problems. [30, 31]. 7.1 T he D irect A pproach Kriegman pioneered the study of stable poses in [30] for piecewise smooth objects. We briefly summarize his approach here, and supplement his discussions by mention ing the difficulties of the direct approach in the verification of local minimality and in dealing with singularities. We will use a slightly different notation. In particular we place the coordinate system so that the center of mass of K coincides with the origin. 71 7.1.1 Sum m ary o f th e C ase A nalysis The conceptually simplest of all possible cases is a pose with only one contact point p located in the relative interior of an elliptic face. A necessary condition for the pose at p to be stable is that S (P )= ° (7.1) p x Vs(p) = 0 This condition characterizes the critical points of the potential energy function. Correspondingly, we will call poses satisfying (7.1) the critical poses. For a critical pose to represent a local minimum as distinguished from local maximum or a saddle point of the potential energy, p also has to satisfy |p| < m in(-|-, — ) (7.2) k\ a? 2 where k\ and k2 are the principle curvatures of surface s at p. Other kinds of critical poses may be analyzed similarly. In each case a set of equations similar to (7.1) may be derived which provide the necessary and sufficient condition for a combination of contact points to determine a critical pose of K . We simply list all possibilities of the set of contact points as discussed in [30] without listing the corresponding equations: • one regular point in a 2-face (discussed above); • two surface points; • one surface point and one edge point; • two edge points; • three surface points; • two surface points and one edge point; • one surface point and two edge points; • three edge points. 72 To make the above a complete list we would have to include cases in which vertices are also involved. The equations involving the vertices are relatively simple and are explicitly excluded from [30]. We will not pursue this issue as our approach will require us to enumerate a different and much smaller set of possibilities. The implementation actually takes triples and pairs of faces/edges which de termine all potential critical points. The homotopy continuation method is then applied to solving each system of equations and find the contact points numerically. The running time is 0 ( N 3dI 9) for an object consisting of N features, each of which is described by algebraic equations of degree at most D. As noted in the paper, it takes about two days to compute the stable poses for a part similar to Figure 7.1 on a Sun SPARCstation 1. 7.1.2 D istin gu ish in g Local M inim a from O ther C ritical P oses Once a critical pose is found, one has to check 1. whether the critical pose is a local minimum and 2. whether the corresponding supporting plane intersects any other portion of K . The second test as described in [30] involves solving less costly systems of equations for each feature on the object (and hence does not contribute an extra factor to the complexity, assuming N is not too huge). For the first test, one need a criterion similar to the inequalities (7.2) for the other cases. In the three-point contact cases the critical pose is a local minima if and only if the supporting plane the projection of the center of mass lies within the triangle formed by the three points of contact. The cases involving two-point contact is less straightforward. We will take the two-surface-point case as an example. Let p x G Si and p 2 € s2 be the points of contact in a critical pose. As given in [30], p x and p 2 satisfy: | 3t -(p,-) = 0 < n<(pi) x n 2(p2) - 0n i(p 1) • (p x - p 2) = 0 1 Pi x P2 • m (p !) = 0 Once we find the solution p x = p x and p 2 = p 2, we have to verify th at a small perturbation of the orientation of the object will move it towards poses with higher 73 opipe.pri Figure 7.2: Another PVC pipe fitting. 74 as opposed to lower potential energy. Intuitively there are two distinct “princi ple” ways of perturbing K out of the critical pose—around the axis determined by Fvecp1Fvecp2 or around the line on the supporting plane which passes the projec tion of O and meets p 1Fvecp2 orthogonally. Other directions of perturbation are possible but it is not clear how to provide a mathematically satisfactory argument which states that it suffices to consider these two kinds of perturbations using the direct approach. Later on we shall argue, taking the dual approach, more rigorously th at this “generally” true, and point out precisely when it fails to be the case. The critical pose is stable with respect to the second kind of perturbation if and only if on the supporting plane, the projection of the center of mass falls between the two points of contact pq and Fvecp2. The other kind of perturbation corresponds to the one-dimensional solution set to (7.3) with the last equation removed. One way to verify the minimality of the solution (px, p2) is to choose a convenient parame terization of the solution curve (in the configuration space) and compute indirectly (using the chain rule) the first and second directional derivatives of the potential energy along the solution curve. It can then be evaluated at run time for each criti cal pose. The resulting formula can be formidable with a system in six variables as such. An alternative is not to compute the second derivative symbolically a priori but to approximate it numerically at run time using the first derivative. 7.1.3 C oincidences and Singularities More often than not, industrial parts are designed with certain symmetry and co planarity. For example, many part are designed so that when they rest on a hori zontal plane, four points or rather four facets come in touch with the plane. It is hard to give a catalog of all kinds of possible singularities. The “stacked ellipsoids” in Figure 7.3 shows another example. Besides the two obvious stable poses of one- point contact, there are two other poses each making a 3-point contact with the supporting plane. In terms of equation solving, we have a one-point contact that gives a critical pose (but not stable by itself). Two such critical solutions together constitute a stable pose. However, a third such critical solution happens to give the same pose. Many common parts consist of cylinders, cones, or other surfaces of revolution, and it is not uncommon for their center of mass to be designed to 75 lie on the axis of symmetry. The case analysis listed in [30] and reviewed above, when supplemented by the other cases involving vertices, give a catalog of cases to which other cases can be reduced. The cost is, however, expensive in terms of the number of coplanar features. For a critical pose consisting of k contact features, the implemented algorithm as described will compute Q(k3) sets of equations and as many number of tests on minimality, only to arrive at the same conclusion for this one critical pose. Then one faces the task of partitioning all of the computed stable poses into equivalent classes, each of which consists of numerically nearby but symbolically distinct solutions. Agrawal studies various kinds of coincidences in his thesis [3] and proposes a gen eral approach for dealing with them. Of particular interest to us in our application here are coincidence by design as exemplified above, and coincidence by theorem. Considering the number of cases, it is at best tedious to modify the algorithm in [30] to take advantage of the proposed solution by symbolically sharing their parameters. If all of the coefficients in the input equations representing the faces are integral (or rational, for th at m atter), such coincidence can be handled as if it were coincidence by theorem. In our context it means that we may use interval arithm etic to com pute the solutions and detect (or reject) the coincidence based on the identity (or the lack of it) of the involved polynomials as suggested by random evaluations. Either way, such overhead will render the direct approach even less attractive in regards to running time. Coincidence and symmetry may bring up singularities in terms of equation solv ing, which require the algorithm to do more than just recognizing the coincidence and solving only a representative set of equations. Consider the toy made of two balls and a connecting rod in Figure 7.4. Suppose that the center of mass lies exactly at the center of symmetry (half way between the balls, on the axis of the rod). This is not a cylinder by itself, which is treated specially in the implementation of [30] and whose symmetry is utilized to speed up computations. Such a two-point contact corresponds to a one-dimensional set of solutions in the configuration space and is likely to cause numerical problems. We have used simple examples to illustrate the variety of possible special cases which may arise due to co-planarity and singularity. W hen the input surfaces are of higher degrees, more complicated and interesting phenomena may conceivably 76 Figure 7.3: Collinear “two-point” contacts. ’ wlt.pri’ ----- Figure 7.4: Singularity caused by symmetry. 77 occur. Compounded with the issue of identifying stable poses among critical poses, and with the number of different cases (recall that we left out the vertex cases), a thorough and reliable implementation based on the direct approach will quickly get out of the hand. 7.2 T he D ual Approach As noted in the paper, one weakness of the direct approach is the lengthy execution tim e spent in solving systems of equations. We will look at the dual approach which not only reduces the number of cases to consider but also admits much faster nu merical methods when the objects under consideration are quadratic surface patches cut by planes. Moreover, the recognition of minimality, coincidences, and singular ities become more intuitive. It is more straightforward to exhaustively enumerate, with respect to inputs of bounded degrees, all possible degenerate cases which may arise. It is also a less formidable task to construct a complete implementation which handles such special cases correctly. 7.2.1 C ase A nalysis We begin by stating the problem in a slightly different way. Given a curved polyhe dron K in space, what are the possible planes which may serve as the floor for K to rest upon stably? The candidates consist of all the supporting planes of K , which may be parameterized by their orientation. In our terms, the stable poses corre spond to the local minima of the potential energy as a function of the orientation of the supporting plane. For our purpose, the potential energy is in turn equivalent to the distance of the supporting plane from the origin. Thus we choose the center of mass to be the origin and find K*. We shall then find the local maxima of the distance |p| from any point p on d(K*) to the origin since this distance is the reciprocal of the distance from the corresponding supporting plane p to the origin. Neglecting coincidences and singularities for the moment, there are three cases which an implementation taking the dual approach will have to consider. Two cases consist of six or more sub-cases whose distinctions 78 are irrelevant to the stable pose algorithm but are implicitly taken care of by the dual mechanism. 1. p is a vertex of dK*. It cannot be the dual of a facet present in the input since facets do not contribute to the convex hull and are not dualized at all. (Recall that by faces we usually mean open faces and that the contribution of a closed facet is accounted for by its O-subfaces.) For similar reasons, it cannot be the intersection of a curve and a face in the dual space. It is generally the intersection of three faces, each of which may be the dual of either a vertex, an edge, or a 2-face. These cases correspond to the three-point contact cases in the direct approach. 2. p is a regular edge point. The edge is generally the intersection of two faces in the dual space, each of which may be the dual of a vertex, an edge, or a face. These cases correspond to the two-point contact cases in the direct approach. 3. p is a regular surface point. In this case the surface < p can only be the dual of a face. It cannot be the dual of an edge for otherwise < p would be a developable surface and it is not difficult to verify that there cannot be a local maximum of the function x 2 + y2 + z2 on such a surface. 7.2.2 C oincidences The various difficulties in identifying the stable poses, in the co-planarity or co linearity, and in singularities of the systems of equations are also present in the dual approach, but in a manner that admits a more systematic treatm ent. Let us first consider the following change of variables in the open halfspace x > 0 of the dual space: and y x + \Jx 2 + y2 z \ /x 2 4-y 2 + V x2 + y 2 + z 2 \Jx2 + y2 + z 2 It is readily verified that this change of variable is bijective and continuously differ entiable up to any order between the open halfspace x > 0 and the corresponding open set { (s ,t,r) : — 1 < s,t < 1 and r > 0}. The family of concentric spheres centered at the origin are mapped to the iso-surfaces r = c for each fixed c in the (s, t, r)-param eter space. Thus dK* D {(x ,y ,z)T : x > 0} corresponds to the graph of a well-defined function r = f ( s , t ) of s and t. Note that / is only piecewise smooth. The stable poses correspond to the local maxima of the graph of / . The three cases in the dual approach remain topologically the same when considered in the pa ram eter space as the change of variables does not alter the smoothness (or roughness for th at m atter) of dK*. Thus the candidate critical poses correspond, in the pa ram eter space, to vertices, to regular curve points at which the tangent is horizontal, or to regular surface points at which V / = 0, All situations involving co-planarity of more than three contact points are accounted for by vertices of degree higher than 3. Similarly, co-linearity of more than two contact points are accounted for by points belonging to an edge which is the intersection of more than two smooth surfaces. 7.2.3 Singularities Singularities corresponding to surfaces of revolutions manifest in the param eter space as horizontal “ridges” along which all poses are equally stable. More generally, a regular surface point at which the Hessian of / does not have a full rank corresponds to a singular pose. • A point with vanishing Hessian corresponds to contact at an umbilical surface point p where the osculating spheres to any surface curves passing p at p are centered at the center of mass. 80 • A point with rank 1 Hessian corresponds to one of the following situations: — contact at a regular surface point p where one of the circles of curvature along the principle directions at p is centered at the center of mass. - contact at a regular edge point p where the circle of curvature at p is centered at the center of mass. In such cases higher differentials up to the maximum degree possible, and only these, need be examined to decide whether the the corresponding critical pose is a stable pose. Similarly, we may have singular poses corresponding to a regular edge point in the param eter space with vanishing second order directional derivative. In this case the osculating plane of the edge point is horizontal. The corresponding 2-point contact in the primal space is less intuitive as explained below. Note th at each of the two contributing features may be either vertices, edges, or 2-faces, although they cannot be both vertices simultaneously. Also recall th at the edge corresponds to, in the primal space, the envelope of the two features, which is a developable surface. Now in the dual space such an edge corresponds to a curve whose osculating sphere at this point is centered at the origin. There are two cases in the dual space: Case A The osculating plane does not pass the origin. Then in the primal space there is a unique cone such that it contains the line of support as a ruling, its vertex lies on the edge of regression of the envelope surface, and its axis passes the center of the mass. This cone has a contact of order two or higher with the envelope and represents a critical pose which may or may not be stable. Case B The osculating plane passes the origin. Then in the primal space there is a unique cylinder such that it contains the line of support as a ruling, its axis passes the center of mass (and hence its axis is horizontal), and its axis passes the center of the mass. This cylinder has a contact of order two or higher with the envelope and represents a critical pose which may or may not be stable. Figure 7.4 shows such an example. In either case, higher order derivatives of the edge in the param eter space must be investigated to determine whether the pose is stable. 8 1 7.2.4 A C om p lete A lgorith m We may now sketch how an implementation can identify the stable poses among all critical poses while taking coincidences and singularities into account. A suitable coincidence resolution scheme such as [3] is first applied to resolve vertex and edge coincidences. Next the parameters corresponding to the critical poses are computed as described in the three major cases. Each critical point p in the param eter space may have zero or more edges incident upon it depending on the case as well as on the number of coincidences at that point. To make sure that p is a local maximum (in the param eter space), we test whether / is “negative definite” around p by iterating over all incident edges as well as all surrounding surface sectors. The test for each edge or sector is performed and singularities are examined by looking at as many higher order directional derivatives or differentials as necessary, until a conclusion is reached or up to the upper bound of the surface degree. In the surface sector case the range of the sector has to be taken into account. 7.3 C apture R egions The technique in the previous section naturally extends to the computation of cap ture regions—the regions in the parameter space within each of which all points correspond to starting poses that will eventually converge to the same stable pose. Kriegman’s treatm ent on this subject is rather detailed for both curved objects and polyhedra. [31] One unsatisfactory aspect, however, is that the analysis is based on the stratified Morse theory, and consequently the functions involved have to be Morse. For curved objects (smooth all over), this means th at his analysis and im plementation would fail when the potential function contains points at which the Hessian m atrix is singular. Thus the singularities discussed in the previous section, which in practice may arise due to design, are explicitly excluded. Also excluded are the coincidences such as co-planarity of more than three contact points and co-linearity of more than two contact points. Finally, piecewise-smooth objects are discussed in [31] but the direct approach does not suggest an efficient algorithm for dealing with such objects. 82 0 .8 /6 0.857 0.837 0.818 0.788 O ./O 'J n wig 0.545 0.525 0.506 0.488 0 .4 6 / 0.447 0.428 0.408 0.389 Figure 7.5: Equi-potential lines (in the param eter space) of the two-piece PVC pipe example. \S \\ t • // / /’ > ^ \ > ^ A t'"/]: W '/mi ’ opipe.rad’ 1.47 1.43 1,39 1.35 1.32 1,28 1.24 1 7 1.17 1.13 1.09 1,05 1.01 0.978 0.938 0.9 0.862 0.824 0.786 0.749 0.711 0.673 0,635 0,587 0.689 0.521 0.483 0.445 0.408 Figure 7.6: Equi-potential lines (in the parameter space) of the three-piece PVC pipe example. 83 Taking the dual approach, however, we are able to deal with all of the problems stated above, at least conceptually. Figures 7.5 and 7.6 show the equi-potential lines of the PVC pipe examples in Figures 7.1 and 7.2, respectively. (The param eterization used is spherical coordinates, not the algebraic change of variables in (7.3).) Finding the capture regions may be accomplished by working in the (.s, t, r)-param eter space, using the algorithm in [31] with appropriate modifications to also process more singular saddles. Theoretically one may use the following asymptotically efficient algorithm which relies on the fact that all input surfaces are algebraic. 1. Find all saddle points, taking into account coincidences and singularities. 2. Find the contours passing the saddle points. Note that these are algebraic curves which intersect only at the computed saddle points and which need not be approximated or traced numerically as done in [31]. 3. Apply a plane sweep algorithm to construct an arrangement consisting of these contours and possibly additional vertical lines. 4. For each peak p (corresponding to a stable pose) locate it in the arrangement. Do a breadth-first search of the adjacent faces starting from the face containing p in order to mark all faces in the arrangement which converge to p. If the input object consists of N features, then the number of smooth pieces which constitute the potential energy function in the param eter space is 0 ( N 2+t) where the constant depends on the maximum degree. Thus there are at most (constant times) this many saddle points. The number of contour curves passing each saddle point is determined by the severity of coincidence and singularity. Each contour curve may consist of up to 0 ( N 2+e) smooth pieces. It is then readily verified that the combinatorial complexity of this algorithm is 0 ( N 4+t) since the plane sweep algorithm need not worry about intersecting contours. 7.4 P ractical C onsiderations We have largely left out the algebraic complexity in the previous discussions in our treatm ent of critical poses and capture regions. Yet in an actual implementation 84 it is the algebraic computations that dominates the computation tim e as pointed out by Kriegman. The transformation introduced in Section 7.2 is algebraic in the relevant direction and does not increase the degrees of the dual surfaces by more than a factor of 5. (In fact a slightly modified variant can be used to achieve a factor of 4.) However, the dualization process may bring a degree D surface to a surface of degree up to D4. Since solving systems of equations is so costly, we might as well evaluate the potential energy all over the places, which can be realized by the ray-shooting tech nique introduced in Section 6.6. The angle between two rays in the dual space is the same as that between the two poses they determine in the primal space. The precision measure discussed there can be used as a param eter the user controls di rectly. If the precise poses have to be computed, the graph describing the potential energy function can be constructed from the mesh. The precise poses can then be computed using the nearby vertices as initial guesses to simple numerical techniques such as Newton’s method, which is much less costly than the homotopy method. A missing feature in this context corresponds to a small capture region, which is likely what the user wouldn’t want anyways. Next let us focus on the special case in which all faces of K are quadratic surface patches bounded by planes. The edge case (typically a two-point contact) can be resolved by parameterizing two variables in the other one (as the solutions of quartic equations), substituting into x 2 + y2+ z2, and finally minimizing a univariate function with radicals. The vertex case (typically a three-point contact) requires us to solve a system of 3 quadratic or linear equations. It may be transformed to a root-finding problem for a univariate function as in the previous case; or one may apply Newton’s method. In the latter case a good initial guess may be found by locating the octant, with respect to the normal forms of each surface, which contains the solution. Now we consider the face case (one-point contact). We will derive the computa tion involved in the general case where the constant term of the quadratic equation does not vanish. The other case is similar but simpler and hence omitted. Let the host surface in question be given as {x : /(x ) = 1 + 2bTx + x t A x = 0}. 85 and let <?(x) = x\ + x\ + x\. Finally let L(xu a r2,z 3, A) = g (x i,x 2,x 3) + X f{xu x 2,x 3) where A is the Lagrange multiplier to be determined. Equating to zero the partial derivatives of L with respect to x \ ,x 2, and x 3, we find that (7 + AA)x + Ab = 0. Expressing x in terms of A and substituting it into the fourth partial derivative of L we obtain a polynomial of degree 6 in A . Thus in each of the three cases, the problem can be reduced to solving a single equation in one variable, albeit non-algebraic. Equivalently, the system in each case is algebraic of total degree at most 8, which is a significant improvement over the pessimistic estim ate of 512 in [30]. 86 ’opipe.rad’ ----- W /jm m m rrnn^ Figure 7.7: The radius function of the three-piece PVC pipe example. 87 C hapter 8 W id th Functions and Tight B ounding B oxes In this chapter we study two other applications of convex hulls of curved polyhe- dra. The ideas and techniques used in constructing the potential energy function of an object can be applied directly to constructing the width function of an object K . Once the width function of K is constructed, it may then be used to find the orientation and dimensions of the tight bounding boxes of K with minimum area, volume, or some other metric. 8.1 W id th Function Let K be a convex curved polyhedron whose interior contains the origin O. Given a fixed direction represented by a vector u, one may determine the distance between the pair of supporting planes of K whose normal vectors are parallel to u. We now study the width function—this distance as a function of the given direction. The need to compute width functions arise in a few applications. In an autom ated manufacturing environment where parts are fed to robotic manipulators, the width function of the planar part is essential to the construction of an oblivious grasping plan [36, 41, 11]. The special case in which the part is a planar polygon is studied in [20, 21] for use in such constructions for orienting parts using parallel jaw grippers. Another application is in the interpretation of the results obtained from a coordinate measuring machine (CMM), an instrument used to obtain the actual coordinates of sampled points on the manufactured surface designed to be flat. The distance between the closest pair of parallel planes enclosing the measured set of points serves 88 as a good measure for “flatness” of the manufactured surface and can be compared with the tolerance. To construct the width function of a curved polyhedron, we first find its polar set either by an ordinary cell construction algorithm or by ray shooting. Each straight line I passing the origin in the dual space represents a family of parallel planes in a fixed direction. It intersects dK* at two points I a and Ib lying on the opposite sides of the origin and representing an enclosing pair of supporting parallel planes of K . The width of K along the direction of I is therefore 1/|A ;| + 1 Clearly the width function is non-smooth only at places where either A or B lies on an edge or coincides with a vertex of dK*. Thus the decomposition (planar graph) of the domain of the “radius” function may be superimposed on a suitable mirror image of its own. The resulting planar graph gives a decomposition of the param eter space such th at the width function is smooth and can be represented symbolically within each face. One convenient choice of such “radius” function is the potential energy function, although the choice of the origin is clearly irrelevant in constructing the width function. The two-dimensional analogs are similar and simpler. Figure 8.1 shows a general curved polygon in the plane, which is disconnected and non-convex, along with its convex hull. Its polar set is depicted in Figure 8.2. The width function is shown in Figure 8.3. In the following table, we summarize the time complexity for constructing a planar graph to describe the width function of a curved polyhedron. planar, simple planar, arbitrary spatial linear Q(N) ©(TVlogiV) 0 ( N 2) log iV curved Q(N) Q{N log N ) 0 ( N 4+‘) For the four planar cases, once the polar set is constructed, it takes only merging two lists whose elements signal the parameters corresponding to the non-smooth points of the width function. The time complexity is obviously optimal for the two simple sub-cases. It is also optimal for the two arbitrary sub-cases because the convex hull can be reconstructed in linear time given the width function. In the spatial linear case, constructing the polar set takes 0 ( N log N). The complexity of the boundary of the point set of polyhedron, however, is linear. Thus the two copies of 89 Figure 8.1: An “object” K consisting of two ellipses and a line segment; shown along with its convex hull. Figure 8.2: The dual of K. The solid line indicates the boundary of its polar set, K*. 90 the boundary can be merged in quadratic time with a logarithmic overhead. For a curved polyhedron in space, the 0 ( N 2+e) structural complexity results in the super- quartic running time of merging two planar subdivisions. Note that the complexity for a curved polyhedron is not completely general in that each 2-face of the input polyhedron is assumed to be bordered by a bounded number of edges. We finally mention th at for an ordinary (non-curved) polyhedron, a faster algorithm exists if one only needs the overall width of the polyhedron and not a complete description of the width function. 8.2 Planar Tight B ounding B oxes Given a curved polygon K in the plane, a parallel pair of supporting lines which enclose K can be thought of as a two-gon which bounds K tightly in a fixed direction. It is natural to extend the idea and consider the area, for example, of the box whose sides are parallel to a given direction and which tightly bounds K . The area is a piecewise smooth function of the the given direction. A complete description of this function may be constructed by superimposing four copies of the “radius” function similar to the potential energy function, each rotated by | more than its predecessor. Figure 8.4 shows this function for the curved polygon in Figure 8.1. From it the tight bounding boxes of maximum and minimum area are found and depicted in Figure 8.5. As is the case with the potential energy function, coincidences and singularities may happen. For example, we expect a bounding box to have generally four contact points with K , although with certain orientations the number of contact points may be more, or there may even be a line segment on K that coincides with one of the edges of the bounding box. Similarly, in order to decide whether the tight bounding box at a certain orientation is locally maximal or minimal, it may be necessary to look at higher derivatives of the area function if both the first and second derivatives vanish. We will not repeat a similar discussion that gives a catalog of the interesting situations in the primal space which correspond to coincidences or singularities in the dual space. It should be clear that by taking the dual approach, the problem of finding the smallest bounding box of a convex K with N is reduced to solving 0 ( N ) instances of non-linear (smooth) optimization problems. It would be at best tedious 91 1 6 14 12 10 8 6 4 2 0 3.14159 6.28319 0 Figure 8.3: The only discontinuities in the first derivative of the width function occur at either those of the radius of K as a function of the polar angle, or those of the same function translated by 7 r. (All three functions are shown.) 200 r -i 180 160 140 120 100 80 60 40 20 0 1.5708 0 Figure 8.4: The area of the tightest bounding box as a function of the polar angle. Only one period (0 to 7 t / 2) is shown. 9 2 to enumerate all possible combinations which give the optimization constraints in the primal space. Two generalizations and their meaningful combinations follow immediately. First, the optimality criterion may be chosen from a wide variety of possibilities such as minimal circumference, tight with maximal shortest edge, tight with minimal longest edge, etc. Each gives rise to a different optimization problem, but does not change the combinatorial structure of the problem. Secondly, the domain of optimization may be all fc-gons whose inner angles in a circular order are the same as a given k-gon, or it may be a reasonable restriction of the above such as all k-gons similar to the given one. In particular minimal bounding squares can be found in the same order of time. In general for a k-gonal problem, the combinatorial complexity will be 0 ( N log k) since the merge has to be performed log k times. Note however, that the increase in the algebraic complexity caused by considering a bounding A:-gon instead of a bounding rectangle or triangle, may or may not be higher than just a factor of k, depending on the numerical method used to solve each instance of the optimization problem. 93 Figure 8.5: The tight bounding boxes of maximum and minimum area as found from the previous picture. 94 C hapter 9 C onclusions 9.1 C ontributions Combinatorial Complexity In this thesis we have shown that the convex hull of a 3-D curved polyhedron has 0 ( N 2+ C ) faces and can be constructed in randomized expected time 0 ( N 2+€), provided each face of the curved polyhedron has a bounded number of subfaces. The graph of the potential energy function of such an object can be decomposed into 0 ( N 2+e) smooth pieces and can be constructed in randomized expected time 0 ( N 2+t). The capture regions of a 3-D curved polyhedron has a structural complex ity of 0(iV 4+t) and can be constructed in deterministic time 0 ( N 4+e). Similarly, the width function of a 3-D curved polyhedron has a structural complexity of 0 ( N 4+e) and can be constructed in deterministic time 0 ( N 4+e). Two operations with respect to the convex hull of an arbitrary 3-D curved poly hedron can be performed more easily and efficiently without explicit construction of the convex hull. The classification of a point again the convex hull of a 3-D curved polyhedron takes tim e 0 ( N log N). Given a curved polyhedron and a point outside its convex hull, the supporting cone from the point to its convex hull has a structural complexity of 0(iV) and can be constructed in time 0 ( N log N). The convex hull of a 2-D curved polygon can be found in time 0(A r log N ).1 The minimal bounding boxes of a 2-D curved polygon can be found in time 0(fV). R eca ll that a curved polygon may consist of more than one connected component and may have dangling 1-faces and/or 0-faces. 95 Algebraic Complexity The new faces in constructing the convex hull of a curved polyhedron made of quadratic surface patches with planar cut have closed form parameterizations. The critical poses of such a curved polyhedron can be found by solving systems of equations each of total degree at most 8. Miscellaneous It is illustrated by an implementation that numerical approximations of most of the above computations can be efficiently (but not optimally in the combinatorial sense) and safely carried out without resorting to sophisticated coincidence m anipu lation strategies. (The implementation does not suggest how to approximate capture regions efficiently.) The existing optimal algorithms for constructing the convex hull of a simple planar polygon, curved or rectilinear, are significantly simplified. The existing discussions of the stable pose and capture region problems are improved by giving a case analysis of coincidences and singularities complete up to a fixed degree. A simple C++ interface is proposed for the safe and efficient manipulations of submatrices without maintaining reference counts. 9.2 Future R esearch Some of the results concerning curved polyhedra are valid only if the curved poly- hedra satisfy the extra constraint that each face has a bounded number of subfaces. The algebraic complexity of the duals of higher degree surfaces and curves is worth studying considering the dramatic difference between the naive estim ate of the de gree 512 system of equations versus the tight bound of 8 after careful analysis of quadratic surfaces using the dual approach. There is room either for improving the algorithmic complexity or for establishing a worst case scenario in the convex hull of a curved polyhedron and its potential energy function. This is also the case with point classification and supporting cone construction without explicit hull construc tion, and especially with the width function and capture region. The prototype implementation allows us to evaluate the dual approach. Two possible extensions are either enhancing the implementation and make it interface 96 with existing solid modelers, or incorporating a numerical solver for dealing with systems of equations. 97 R eference List [1] P. K. Agarwal, M. Sharir, and P. Shor. Sharp upper and lower bounds on the length general Davenport-Schinzel sequences. Journal of Combinatorial Theory, Series A, 52:228-274, 1989. [2] Pankaj K. Agarwal. Intersection and Decomposition Algorithms for Planar Arrangements. Cambridge University Press, 1991. [3] Amitabh Agrawal. A General Approach to the Design of Robust Algorithms for Geometric Modeling. PhD thesis, Computer Science Department, University of Southern California, 1994. [4] D. S. Arnon and B. Buchberger. Algorithms in Real Algebraic Geometry. Aca demic Press, 1988. [5] Chanderjit Bajaj and Myung-Soo Kim. Convex hull of objects bounded by algebraic curves. Algorithmica, pages 533-553, 1991. [6] J. D. Boissonnat, Cerezo, O. Devillers, J. Duquesne, and M. Yvinec. An algo rithm for constructing the convex hull of a set of spheres in dimension d*. In Proceedings of the Fourth Canadian Conderence on Computational Geometry, August 1992. [7] J. W. Bruce and P. J. Giblin. Curves and Singularities. Cambridge University Press, 2nd edition, 1992. [8] Bernard Chazelle. An optimal convex hull algorithm and new results on cut tings. In Proc. 32nd Annu. IEEE Sympos. Found. Comput. Sci., pages 29-38, 1991. [9] Bernard Chazelle. An optimal algorithm for intersecting three-dimensional con vex polyhedra. SIAM Journal on Computing, 21 (4):671— 696, August 1992. [10] Bernard Chazelle, Herbert Edelsbrunner, Leonidas Guibas, and Micha Sharir. Diameter, width, closest line pair, and parameteric searching. Discrete Com putational Geometry, pages 183-196, 1993. 98 11] Yuibin Chen, Chao-Kuei Hung, and Doug Ierardi. Optimal manipulation strate gies for orienting planar workpieces by pushing, submitted to Algorithmica, 1994. 12 13 14 15 16 17 18 19 20 21 22 23 24 K. L. Clarkson and P. W. Shor. Applications of random sampling in computa tional geometry, (ii). Discrete Computational Geometry, 4:387-421, 1989. David P. Dobkin and Diane L. Souvaine. Computational geometry in a curved world. Algorithmica, pages 421-457, 1990. J. J. Dongarra, R. Pozo, and D. W. Walker. Lapack++: A design overview of object-oriented extensions for high performance linear alge bra. Technical report, University of Tennessee, 1993. software available at ftp: / /netlib.att.com /att / C + + /lapack-f+ /. H. Edelsbrunner, L. Guibas, J. Pach, R. Pollack, R. Seidel, and M. Sharir. Arrangements of curves in the plane: topology, combinatorics, and algorithms. Theoretical Computer Science, 92:319-336, 1992. Herbert Edelsbrunner. Algorithms in Combinatorial Geometry. Springer- Verlag, 1987. Michael C. Gemignani. Elementary Topology. Dover, 1990. Republication of the second edition (1972) of the work originally published in 1967 by Addison- Wesley. Robert Gilmore. Catastrophe Theory for Scientists and Engineers. Dover, 1981. A. Goetz. Introduction to Differential Geometry. Addison-Wesley, 1970. K.Y. Goldberg. A complete algorithm for orienting polygonal parts. Algorith mica, 1993. K.Y. Goldberg and A. Rao. Orienting generalized polygons. Technical Report 275, USC Institute for Robotics and Intelligent Systems, 1991. R. L. Graham and F. F. Yao. Finding the convex hull of a simple polygon. Journal of Algorithms, 4(4):324-331, 1983. Joos Heintz, Marie-Francoise Roy, and Pablo Solerno. Description of the con nected components of a semialgebraic set in single exponential time. Discrete and Computational Geometry, pages 121-140, 1994. Roger A. Horn and Charles R. Johnson. Matrix Analysis. University of Cam bridge, 1985. 99 [25] Doug Ierardi and Dexter Kozen. Synthesis of Parallel Algorithms, chapter 15. Morgan Kaufmann, 1993. [26] Alfred Inselberg. Personal communication. [27] Bill Jones, Hebert Edelsbrunner, and Guibas et al. On-line bibliogra phy of computational geometry, anonymous ftp from c s .u s a s k .c a in file pub/geom etry/geom bib.t a r .Z. [28] Andrew Klapper. A lower bound on the complexity of the convex hull problem for simple polyhedra. Information Processing Letters, 25:159-161, 1987. [29] D. J. Kriegman and J. Ponce. On recognizing and positioning curved 3d objects from image contours. In IEEE Transactions on Pattern Analysis and Machine Intelligence, November 1990. [30] David J. Kriegman. Computing stable poses of piecewise smooth objects. In Computer Vision Graphics Image Processing: Image Understanding, 1991. [31] David J. Kriegman. Let them fall where they may: Capture regions of curved objects and polyhedra. Technical Report 9508, Yale University, June 1995. Subm itted to the International Journal of Robotics Research. I [32] David J. Kriegman, Erliang Yeh, and Jean Ponce. Convex hulls of algebraic curves. In SPIE Conference on Curves and Surfaces in Computer Vision and Graphics, November 1992. [33] Steven R. Lay. Convex Sets and Their Applications. John Wiley & Sons, Inc., 1992. original edition 1982. [34] D. T. Lee. On finding the convex hull of a simple polygon. International Journal of Computing and Information Science, 12(2):87-98, 1983. [35] Dinesh Manocha. Solving polynomial systems for curves, surface, and solid modeling. In Proc. of ACM/SIGGRAPH, 1993. [36] M.T. Mason, K.Y. Goldberg, and R.H. Taylor. Planning sequences of squeeze- grasps to orient and grasp polygonal objects. Technical Report CM U -CS-88- 127, Carnegie Mellon University, Computer Science Department, April 1988. [37] D. McCallum and D. Avis. A linear algorithm for finding the convex hull of a simple polygon. Information Processing Letter, 9:201 - 206, 1979. [38] James R. Miller. Incremental boundary evaluation using inference of edge clas sifications. IEEE Computer Graphics and Applications, pages 71-78, January 1993. 100 [39] Franco P. Preparata and Michael Ian Shamos. Computational Geometry: An Introduction. Springer-Verlag, 1985. [40] Anil Rao, David Kriegman, and Ken Goldberg. Complete algorithms for feeding polyhedral parts using pivot grasps. 1994. Submitted to IEEE Transactions on Robotics & Automaton. [41] Anil S. Rao and Kenneth Y. Goldberg. Manipulating algebraic parts in the plane. In IEEE Transactions on Robotics & Automaton, 1993. [42] Aristides A. G. Requicha and Herbert B. Voelcker. Boolean operations in solid modeling: Boundary evaluation and merging algorithms. In Proceedings of the IEEE, 1985. [43] Alejandro A. Schaffer and Christopher J. Van Wyk. Convex hulls of piecewise- smooth jordan curves. Journal of Algorithms, 8:64-94, 1987. [44] O tto Schreier and Emanuel Sperner. Projective Geometry of n Dimensions. Chelsea Publishing Company, New York, 1985. [45] J.T. Schwartz and M. Sharir. On the two dimensional davenport-schinzel prob lem. Journal of Symbolic Computation, 10:371-393, 1990. [46] Micha Sharir. Almost tight upper bounds for lower envelopes in higher di mensions. In Symposium on Foundations of Computer Science, pages 498-507, 1993. 101 A p p en d ix A Im plem entation N otes We have shown that asymptotically efficient algorithms exist for solving many of the problems treated in this thesis. In reality, however, some of these algorithms are either too complicated to implement or simply incur so much algebraic overhead that it dominates the running time. In most cases, the algebraic complexity also forces the implementation to face coincidences and singularities seriously. An implementation based on these “efficient” algorithms thus may not be able to enumerate a small number of such special cases and explicitly treat them separately. Instead, it must adopt some coherent strategy such as that proposed in [3]. Such lim itation, however, is not intrinsic to the dual approach. We now briefly describe a naive and yet very efficient ray-shooting implementation based on the dual approach. It is most useful as a tool for the visualization of the concepts discussed in this thesis, and for use in combination with existing tools which solve (smooth) optimization problems numerically given good initial guesses. We also discuss a C++ technique we adopted in our implementation which is less relevant to the algorithmic issues but is nevertheless very interesting. A .l Sum m ary Platforms and Compiler The program is written in C++. We have compiled and tested it on an IBM PC with i486DX33 processor running Linux and OS/2. The compiler used is gcc in both cases. 102 Dependency The following libraries/packages are used directly or indirectly by the program. L A P A C K The linear algebra package is linked with the program. Schem e This lisp dialect is used to pre-process the input. We only use the Linux version (snow /stk) to actually pre-process the input. Its output is used on other platforms to test the program. G nuplot This plotting package is used to post-process the mesh file created by our program and to create postscript pictures. The following packages are not needed to compile or run our program but they were indispensable tools during the development of the program. R L aB This linear algebra calculator is first used to create a prototype implemen tation for our algorithm and later used to verify the correctness of our m atrix implementation. P A R I/g p We use this symbolic calculator mostly for polynomial computations. A .2 O bject Specification Input to the program are surface patches specified in LISP syntax. The following is a typical description of a surface patch. ( ( c la s s 2) (s c a le #(1.2 1.2 0)) (ro t_ ty p e ezyx) (r o t _an g le #(0 90 0)) (c e n te r #(-0.1 - 1.1 1 .1)) (u .ra n g e #(-180 180)) (v_range #(0 1 .5)) ) 103 The meanings of these fields are explained as follows. class Quadratic surfaces are classified according to their normal forms as in [44]. The possible values are: 0 Cone or plane. 1 Parabolic cylinder, hyperbolic paraboloid, or elliptic paraboloid. 2 Hyperbolic cylinder, elliptic cylinder, hyperboloid of two sheets, hyperboloid of one sheet, or ellipsoid scale Scaling factors along the three principle axes of the surface. rot_type Name of the rotation scheme. rot_angle Rotation parameters to bring the surface from the standard orientation to the desired orientation. center Translation to bring the surface from the standard position to the desired position. u_range The range of the longitudinal parameter. v_range The range of the latitudinal parameter. The description scheme begins with a representative quadratic surface. In our ex ample, a class 2 surface with two positive axis scales and one vanishing axis scale is interpreted as an infinite cone whose axis coincides with the Z-axis. It happens to be a circular cylinder with a 1.2 radius. This cylinder is first rotated using the Euler-Z-Y-X angle specification and its axis become aligned with the Y-axis. Its center is then moved to (— 0.1 — 1.11.1). Finally a patch is cut out from the infinite cylinder as our desired input. The longitudinal param eter specifies the range, in degrees, of the angular param eter around the axis; and the latitudinal param eter specifies the extent in the Z-direction. 104 A .3 F unctionalities and L im itations The program uses ray shooting in the dual space to find the boundary of the polar set. Given an object description file as input, it may output any combination of the following information. m esh o f th e object containing regularly sampled points on the input surface patches m esh o f th e boundary o f th e polar set containing coordinates of points on dK* which are sampled along fixed intervals of longitudes and latitudes m esh of th e radius as a function of the longitudes and latitudes m esh o f th e boundary of convex hull in their Euclidean coordinates. The resolution of the mesh can be specified, although it only has a first-order effect on the output. Currently not all legal inputs to the program are supported. Specifically, the following features are not implemented: 1. cones, planes, hyperboloids, paraboloids, and all parabolic surfaces; 2. cut along longitudinal lines; 3. cut along non-parameter lines; 4. program-chosen origin; All of these features can be added without much difficulty except planes and planar cut along non-parameter lines. The difficulty associated with these two cases is unrelated to the dual approach. In the primal space it is more tedious to create a satisfactory mesh for a planar patch bounded by conic sections. The same is true with surface patches cut along non-parametric lines. The simplicity of the dual approach allows our implementation to take a naive strategy for dealing with coincidences. A macro n e a rly _ z e ro is defined and used consistently in all situations where there may be numerical difficulties. W ithout compiler or manual optimization, it takes about 20 seconds on a 486DX33 running Linux to produce three (the dual, the radius, and the hull) 58 by 29 meshes for the pipe example consisting of six surface patches and ten edges. 105 A .4 Subm atrices: to Copy or to R eference? One of the efficiency issues most commonly raised against the use of C++ is the autom atic invocation of copy constructors for function arguments as well as for their return values. Consider the following simple piece of code for m atrix computation: MatrixDouble A, B, C; C = A * B; If the m atrix operators are not properly defined, there can be as many as four ephemeral matrices which are destroyed almost as soon as they are created. A good design of the class provider will result in the creation and destruction of just one tem porary matrix. Such cost can be reduced to minimum in a series of m atrix computations if the class provides a set of destructive operators/functions which return a reference to one of its arguments, and if the class user makes a judicious use of these operators. Thus, for example, an expression like A *= B += G / = D; will perform three operations without creating any ephemeral objects. Our discussion is based on the assumption that we have agreed upon such design and use of the m atrix class. The question we would like to address is: how do we deal with submatrices efficiently? Do we make a copy of the submatrix and create a new m atrix object, or do we just create a reference to the original m atrix? The first alternative is obviously very inefficient. It is not very pleasant to implement the second alternative, nor is it clear how one should handle user expressions like A.sub(0,0,n/2,n) += A.sub(n/2+l,0 ,n ,n ) ; in which both arguments involve references to the same m atrix. The Lapack++ library [14] basically takes the first alternative but also allows the user to create matrices which represent frequently used references to the submatrices. This is not very satisfactory as the class provider still has to maintain a reference count and the user still suffers from performance penalty when dealing with submatrices which are used only once. We propose a submatrix interface which 106 • avoids unnecessary copying of tlie m atrix entries; • does not require the class provider to m aintain reference counts; • allows nested expressions to contain references to the same matrix; and yet • does not create garbage in the free store. We begin with the declaration of the reference class which has little data content but provides almost all of the m atrix functions: class MatDblRef { / / reference to submatrix frie n d class MatrixDouble; protected: MatrixDouble & h o s t; in t _srO, _scO, _snr, _snc; / / location and size of the submatrix referenced by *th is MatDblRef(MatrixDouble *h, in t i , in t j , in t nr, in t n c ) ; void operator =(MatDblRef const & A) { a s s e r t(0); } public: MatDblRef & pre_mul_eq(MatDblRef const & A, char t = 'nO; MatDblRef & post_mul_eq(MatDblRef const & A, char t = 'n'); ... / / more arithm etic functions and operators }; An object belonging to this class is but a reference to part of a true m atrix (host) containing the actual array of double values. Moreover, ordinary users are not even allowed to create such an object since the only constructor is protected. However, all of the ordinary m atrix operations one would want to use are defined in terms of objects of this class. How could the user use such an object if s/he cannot even create one? The obvious solution is to create a class publicly derived from MatDblRef. The bulky content of an ordinary m atrix—the array of doubles is also defined in this class. class MatrixDouble : public Matrix<double>, public MatDblRef { 107 frie n d class MatDblRef; public: / / constructors MatrixDouble(int m, in t n) : Matrix<double>(m,n), MatDblRef(this, 0, 0, m, n) {} MatrixDouble(MatDblRef const & src) : Matrix<double>(src. _snr, src._snc), MatDblRef(this, 0, 0, src._ sn r, src._snc) { copy(src) ; } / / submatrix access MatDblRef sub(int i , in t j , in t nr, in t nc) { retu rn MatDblRef(this, i,j,nr,nc); } MatDblRef rows(int i , in t nr=l) { retu rn s u b (i,0 ,n r,_ n c o ls ); } MatDblRef c o ls (in t j , in t nc=l) { retu rn sub(0,j,_nrow s,nc); } >; Now an object of class MatrixDouble, once created, is also a reference to the entirety of itself, forever (that is, until it dies). Thus an ordinary user may use all of the arithm etic functions/operators on a MatrixDouble object defined in the parent class as if there were just one ordinary m atrix class. The interesting mechanics begins to work when the user wishes to use subm atri ces, either as the source or as the destination of an operation. The member functions of MatrixDouble for accessing submatrices creates objects of class MatDblRef whose host is the MatrixDouble object in question. The user may use these subm atri ces/references freely in those arithm etic functions/operators defined for MatDblRef in a complicated expression without incurring more overhead than copying 4 integers and one pointer (by which the reference host is implemented) each time. If the user tries to save the result of the final expression, however, a copy will be created out of necessity, by the second constructor in MatrixDouble. Thus the user retains the option of writing efficient expressions without com promising the integrity of the system. W ith our design, there are other subtleties which must be handled carefully, but they are pretty much standard C++ pitfalls and techniques and hence om itted from our discussion. 108

Linked assets