COMPUTATIONAL GEOMETRY ON THE SPHERE WITH APPLICATION TO AUTOMATED MACHINING Lin-Lin Chen Tony C. Woo Department of Industrial and Operations Engineering The University of Michigan Ann Arbor, Michigan 48105-2117 August 1989 Abstact Based on observations made on the geometry of the cutting tools and the degrees of freedom in 3-, 4-, 5-axis numerical control machines, a new class of geometric algorithms is induced - on the unit Sphere. Spherical algorithms are useful for determining the type of tool, its path, workpiece fixturing, and the type of machine. Basic to these algorithms are four that are presented here: detection of convexity on the sphere, computation for spherical convex hull, determination of the spherical convexity of a union and the intersection of hemispheres. These four algorithms are related by duality and the sharing of partial results. 1

1. Introduction The term "computational geometry" has two interpretations. In his address to the Royal Society, Robin Forrest [4] used the term in the context of computations on analytically complex surfaces for geometric design. As continuity is important to the design of complex surfaces, tangents, normals, and their (dot and cross) products are basic to the computation. The other interpretation stems from an almost orthogonal field in geometry - one that is combinatorially complex. Shamos [16], in his celebrated dissertation, gave a procedural treatise of classical geometric constructs, such as the convex hull. Notably, in the discrete domain of vertices, edges, and polygons, he injected the rigor of algorithm analysis. This paper is an attempt to unify these two fields of research by addressing the algorithmic issues of analytically complex geometric entities. Its roots are humble - numerical control (NC) machining and workpiece setup. But one of the offsprings is rather pleasant - geometric algorithms on the unit sphere, or spherical algorithms, the basics of which are presented here. The applications for which the algorithms are developed involves k-axis NC machines, where k = 3, 4, or 5. Given an object represented by N free-form surfaces, NC applications seek the number of setups (i.e. orientations of the workpieces) and for each setup the tool path. While, it is customary to compute cutter orientations from surface normals, there are two hidden implications. First, the type of cutter is already fixed. If the orientation of the cutter is set to the normal (with the possible difference in the sign of the vectors), the most "efficient" cutter is the flat-end millt. Second, the type of NC machine is also pre-determined. For most free-form surfaces, the wide range of their surface normals can be accommodated only by a 5-axis NC machine. Now, it is not unreasonable to ask the question: "Can a 3-axis (or a 4-axis) NC machine with a different cutter accomplish the same task?" To t If the cutter orientation is orthogonal to the normal, a side mill is implied. 2

be sure, a 5-axis NC machine is more flexible, but is also much more expensive (hence less available) than a 3- or 4-axis machine. Furthermore, the tool path generated for a 3-axis machine is valid for a 4- or 5-axis machine but not vice versa. These observations lead to the following problem definition. Problem MS: Machine Selection. Given a workpiece consisting of N free-form surfaces, determine the minimum number of axes required of the NC machine and the type of cutter. Problem MS as formulated is not without a caveat. It is intuitive that more setups would be required of a 3-axis machine than that of a 4- or 5-axis machine as illustrated in Figure 1.1. Since mounting and dismounting a workpiece is not considered to be productive, it is therefore useful to able to cut as many surfaces as possible in any single setup. This induce a companion problem definition: Problem WS: Workpiece Setup. Given a workpiece consisting of N free-form surfaces, determine an setup orientation such that the maximum number of surfaces are accessible. <Insert Figure 1.1> One of the useful notions from computational geometry is visibility. Two points are visible to each other if the line segment joining them has no intersection [17]. For each point on the surface, it is possible to compute visible rays (to points at infinity). The intersection of the rays for all points on the surface forms a visibility cone for the surface. Thus, if the workpiece is oriented such that the axis of the cutting tool lies in the same direction as one of the rays in the visibility cone, the corresponding surfaces are accessible for machining. Now, consider the geometry of a cutting tool as related to the geometry of a surface. The directly available information from a point on a surface is the 3

1$ too e**te tool adis tup itatong Y (b) 3-axis machining (a) 5-axis machining Figure 1.1 Machine Selection and Workpiece Setup.

outward-pointing normal. Indeed, for a surface, there exists a Gaussian Map which is the intersection of the surface normals with the unit sphere [7]. Without loss of generality, assume that there are three kinds of cutting tools for a milling machine: flat-end, semi-ball-end, and ball-end. With respect to a normal (i.e. a point in the Gaussian Map or GMap), the three kinds of cutting tools can deviate by an angle v = 0, 0 < v < ic/2, and v = 7/2, respectively. This is illustrated in Figure 1.2. Thus, the visibility of a normal is enhanced by an angle of 0 < v < 7/2, depending on the type of cutting tool used. In the subsequent discussion, the ball-end mill, for which v = i/2, is chosen. <Insert Figure 1.2> The visibility map, or VMap, for a given surface can now be constructed. A VMap is an image on the unit sphere whereby every point in the VMap deviates from a corresponding point in a GMap by v < n/2. The construction of a VMap from a GMap is shown in Figure 1.3. For each point pi in the GMap, there corresponds a hemisphere ci (because v = i7/2). If the GMap is a single point (e.g. the north pole), the VMap is a hemisphere (e.g. the Northern Hemisphere). If the GMap consists of two distinct points, the VMap is the intersection of the two corresponding hemispheres. While a GMap is not necessarily (spherically) convex, a VMap is guaranteed to be convex, as it is formed by the intersection of hemispheres that are convex. (More rigorous definitions of spherical convexity are given in Section 2.) <Insert Figure 1.3> Problem WS can now be addressed. A collection of N surfaces is represented by N VMaps on the unit sphere, some of which may overlap. If two VMaps intersect, there exists a common direction along which both surfaces can be machined in a single setup. The degrees of freedom of the 3-, 4-, 5-axis milling machines (Problem MS) are next related to the VMaps. 4

Flat-end Semi-ball-end v 1 0 0 <v<X /2 Figure 12 The v Angle of a Cutting Tool. Ball-end v =~/2 Figure 1.3 Construction of a Visibility Map.

A 4-axis machine is conceptually a 3-axis machine with a rotating table. This extra degree of freedom in rotation can be represented by a great circle, if the rotation covers 21r; otherwise, it is an arc. This situation is illustrated in Figure 1.4(b). For a single setup, if the given arc (great circle) passes through a VMap, the corresponding surface is machinable. The accessibility of surfaces for a 4-axis machine can therefore be formulated as follows: Given N VMaps on the unit sphere, determine an orientation such that a given arc (or a great circle) passes through them (maximally). The utility of a 5-axis machine is now easily representable in spherical geometry. The two additional degrees of freedom in orthogonal rotation correspond to a spherical rectangle as illustrated by Figure 1.4(c). Thus, given a 5-axis machine, the setup of the workpiece can be formulated as a covering problem: Given N VMaps, determine an orientation such that a given spherical rectangle covers them (maximally). <Insert Figure 1.4> By invoking spherical geometry, the preceding description of NC machining and setup induces a structure as shown in Figure 1.5. In this paper, four basic spherical algorithms are presented. They are (1) the detection of convexity (via hemisphericity), (2) the computation for a (spherical) convex hull, (3) the determination of the convexity of a union, and (4) the determination of the intersection of hemispheres. While these four problems may seem disparate, they are shown to be tightly woven via duality and the sharing of partial results. In this order, the algorithms are given, after a brief introduction of terms in spherical geometry and central projection. <Insert Figure 1.5> 5

y (a) 3-axi (b) 4-axis (c) 5-ais Figure 1.4 Spherical representation of NC Machines.

Central Projection Figure 1.5 Problem Structure

2. Definitions and Representations In spherical geometry, the objects considered are sets of points on the d dimensional unit sphere Sd. Let Ed+l denote the d+l dimensional Euclidean space and p = (xl, x2...,xd+1) denote a point in Ed+l, then Sd =p: Ipl =1) In this paper, the discussion will be on S2, the sphere with radius r = 1 in 3 E. The definitions of the principal objects and important concepts are now briefly reviewed [15]: Point. A point p in S2 is a unit vector in E3, which is denoted by a 3-tuple (x1,x2,x3) in this paper. This representation is chosen (over spherical coordinates) such that the correspondence between S2 and E2 through projection can be easily maintained, as will be further explained in Section 3. Antipodal points. Two points p and q are called antipodal if p = -q, where p is called the antipode of q and vice versa. A pair of antipodal points divides the infinite number of great circles containing them into semicircles. Circle and Hemisphere. A circle in S2 is determined by the intersection of the unit sphere with a plane which meets the sphere at more than one point. If the plane contains the origin, the intersection is a great circle; otherwise, it is a small circle. (In this paper, the term circle is used interchangeably with small circle.) The surface of the sphere is partitioned into two hemispheres by a great circle. A closed hemisphere includes the bounding great circle, while an open hemisphere excludes it. The surface of the sphere is partitioned into two portions by a small circle; the smaller of the two is the interior of a small circle. The center of a circle is a point in the interior and equidistant from all the points on the circle. A circle is represented by its center and radius in this paper, as its counterpart in the Euclidean plane. 6

Line. A line L in S2 is a great circle, the intersection of the unit sphere with a plane containing the origin. (The term line and the term great circle are used interchangeably.) Since the radius of a great circle is always equal to 7/2 (as explained under Distance), it can be fully represented by the unit normal of the intersection plane. Line segment. The line segment joining two distinct points p and q in S2, denoted by pq, is the shorter of the two arcs in the great circle containing p and q. This segment is unique, except for the case when p and q are antipodal points where any great semicircle joining p and q is a segment. Distance. The distance between two distinct points p and q in S2 is the length of the line segment (the geodesic) joining p and q. It is defined by the metric: d(p,q) = cos'1 (p.q) <Insert Figure 2.1> Polygon. A (spherical) polygon P is a closed path on the sphere of connected ordered line segments ab, be,..., pq, qa which do not cut across themselves. The points a, b, c,..., p, q are called the vertices of P. The line segments ab, bc,..., pq, qa are called the edges of P. Hemisphericity. A set of points in S2 is said to be hemispherical if it is contained in some open hemisphere [1]. Convexity. A set of points P in S2 is said to be (spherically) convex, if for any two points p and q in P, the segment joining p and q lies entirely in P. Notice that, a closed hemisphere contains antipodal points but it does not contain all the semicircles joining them. Therefore, it is not convex. Quasi-convexity. A set of points P in S2 is said to be quasi-convex, if, for any pair of non-antipodal points p and q in P, the segment pq is also in P. A closed hemisphere is quasi-convex [18]. 7

point p circle (small circle) antipode of p line segment pq p Figure 2.1 Entities on the Spher

Strict convexity. A set of points P in S is said to be strictly convex, if it is nonempty, includes no antipodal pair, and for any two points p and q in P the segment pq is also in P [18]. Convex Hull. A (spherical) convex hull of a set of points P in S2, denoted by SCH(P), is the boundary of the smallest convex set in S2 containing P [14]. <Insert Figure 2.2> 3 Projecting onto the Euclidean Plane Since extensive study in Computational Geometry has been done in 2- and 3-dimensional Euclidean space, it is natural to convert a spherical problem to one in the plane(if possible), apply some planar algorithms, and then convert the planar solution back for the spherical one. One of the methods that converts a spherical problem to a planar one is the central projection (gnomonic projection). However, it does not apply in certain cases. 3.1 Central Projection Central projection maps a pair of antipodal points in S to a point in the two dimensional projective plane P2, which consists of E2 augmented by a line at infinity. By emitting a ray from the origin, a point (x1,x2,x3) in S2 with X3 > 0 is mapped to (x1/3,2/3,1) in the plane X3 = 1, corresponding to the 2 point (xl/x3,x2/x3) in P. By extending the ray in the opposite direction, (-xl,-x2,-x3), the antipode of the point (x1,x2,x3) is also mapped to (x1/x3,x2/x3,1) in X3 = 1 (and thus (xl/x3,x2/x3) in P 2). A point (xl,x2,x3) with X3 = 0 is mapped onto the line at infinity in P2. Central projection therefore is a two-to-one mapping of points in S2 (except for points on the equator) to 2 points in E. As illustrated in Figure 3.1, a pair of antipodal points, p and 8

spherical convex hull Figure 22 Spherical Convexity

-p, in S2 are mapped to a point p' in E2, while a great circle is mapped to a line through central projection. <Insert Figure 3.1> Since central projection is a two-to-one mapping, by keeping track of the sign of X3 for a point, it provides an easy way of going back and forth 2 2 between a geometrical problem in S2 and a corresponding one in E2 Nevertheless, central projection is not a panacea: Not all spherical problems can be solved by applying planar algorithms to the projected counterparts. An example is the problem of finding the smallest enclosing 2 2 circle for a set of points in S2 because a circle in E2 may not be mapped to a circle in S2. It is necessary, then, to review the properties of central projection that allow a spherical problem to be solved in E2 and the limitations that induce new spherical algorithms. 3.2 Properties of Central Projection The merit of central projection is the one-to-one and onto mapping between a line in E2 and a great circle in S2. Lemma 3.2.1 Central projection establishes an one-to-one and onto mapping between great circles in S2 and lines in E2, except for the equator, the great circle that lies on the plane x3 = 0. Proof. Except the equator, any great circle in S2 determines a plane through the origin in E3, which, when intersecting with the plane X3 = 1, determines a line in the plane. Conversely, a line in the plane X3 = 1 together with the origin determine a plane in E3, which, when intersecting with the unit sphere, determines a great circle in S2. Notice that the plane X3 = 1 is itself a space E2. 0 9

Is2 Figure 3.1 Central Projection

By Lemma 3.2.1, central projection preserves incidence since a point 2 p = (xlx2,x3) on a great circle L in S is mapped through central projection to a point p' = (x1/x3,x2/x3) on L', the unique image of L in E2. Obviously, central projection also preserves collinearity. Because of the above properties, the intersection of spherical entities (lines, line segments, polygons, etc.) can be obtained by applying central projection and planar algorithms. Lemma 3.2.2 Central projection maps a (spherically) convex polygon on S2+, the hemisphere of S2 corresponding to x3 > 0, to a convex polygon in E2. Conversely, the inverse of central projection maps a convex polygon in E2 to a (spherically) convex polygon on S2+. Proof. Let P be a convex polygon on S2+ and P' be the image of P through central projection. A convex cone in E3 with the apex at the origin can be constructed by intersecting all the halfspaces each of which is defined by a plane through the origin and an edge of P and contains P. The intersection of this convex cone with the plane X3 = 1 is P', which is convex since it is the intersection of two convex sets in E3. The converse can be shown similarly byconstructing a convex cone in Estruct using convex the edges of P' and the origin. As illustrated in Figure 3.2, a spherically convex polygon not crossing the equator is mapped to a convex polygons in E2 through central projection. On the other hand, a spherical convex polygon crossing the equator is mapped to two open polygons. By Lemma 3.2.2, if a proper rotation is known for a set of points P on 52 such that, after the rotation, the set will be on S2+, then the convex hull for the set can be constructed by applying planar convex hull algorithm to P', the image of P in E2 through central projection. <Insert Figure 3.2> 10

Figure 32 Central Projection of a Polygon

3.3 Limitations of Central Projection For certain problems such as finding the smallest enclosing circle for a set of points in S 2, central projection is not applicable because of the intrinsic nature of these problems. In order to apply central projection (or any other projective mapping), some knowledge of the normal of a good projection plane, e.g. (0,0,1) for x3 = 1, is required. Since a circle in E2 is mapped to a circle in S2 if and only if the normal of the projection plane is the normal of the plane contains the circle in S2, finding such a projection and finding the smallest enclosing circle are essentially the same problem that must be solved on the sphere. Central projection suffers from two limitations. First, central projection does not preserve angle, distance, or area. Therefore, in general, proximity problems, such as the construction of the Voronoi diagram or minimum spanning tree, cannot be solved by applying central projection and planar algorithms. Second, central projection produces largely distorted images for points close to the equator in S2. (This is a problem that needs to be taken care of, because planar algorithms usually do not handle points at infinity.) Hence, it is important to find some initial rotations such that, after applying the rotations to the set of points in S2, no point lies on the equator and therefore none is projected to infinity. It is also essential to find a projective plane such that all projected points are as far away from infinity as possible, so that large distortions and round-off errors can be avoided. Problem P (Projective): Finding Initial Rotations. Given a set of N points P in S2, find some rotations such that, after P is rotated, none of the points in P lies on the equatortt. tt It is possible to formulate a problem of finding the optimal initial rotation: Given a set of N points P in S2, find an optimal rotation such that, after P is rotated, all of the points in P are the farthest away from the equator. This optimal initial rotation can be determined, once the largest empty double wedge on the sphere is found. However,it is shown in [3] that 11

One such algorithm (a slightly modified version due to different dimensionalities) is by Preparata and Muller [13] which finds an initial rotation Ro such that, after applying R0, no point in P is projected to infinity. This algorithm, however, does not try to move the points as far away from the equator as possible. The algorithm proposed here finds two rotations, R1 and R2 that move the points in P as far away from the equator as possible. Let Pi = (il,Xi2,xi3), i = 1,...,N, be points in P. The points Pi =(xil,xi2,xi3) are mapped to Pi' = (0,Xi2,xi3) when parallelly projected onto the plane xl = 0, itself a space E2. The unit circle S1, the image of the unit sphere in the plane x1 = 0, is partitioned into (circular) segments by the unitized pi' and their antipodes in S1. Since the partition is symmetrical with respect to (0,0), only the partition on a semicircle needs to be considered. Arbitrarily, let the semicircle be one of the two determined by the unitized pi' and its antipode. The problem is now transformed to finding the maximum (circular) gap between the points on the semicircle, for which Gonzalez' O(N) algorithm [5, 16] can be used. Once the maximum gap is found, R1 is a rotation that rotates the circular midpoint of the maximum gap to be on the line x3 = 0 (in the plane x1 = 0). After R1 is applied to the unit sphere, the points in P are moved away from the equator, except for the fixed points (1,0,0) and (-1,0,0). Let P' be the result of applying R1 to P. After points in P' are parallelly projected onto the plane x2' = 0, the same procedure is performed again to obtain R2. An illustration of the algorithm is given in Figure 3.3. Given eight points (unit vectors) on the sphere, where point 3, 4, 5, and 6 are on the equator, find a rotation R1 to move points 5 and 6 away from the equator. After the points are parallelly projected to the xi = 0 plane, the largest circular gap is found as shown in Figure 3.3(b). R1 is the transformation that rotates the there are O(N2) empty double wedges. To find the largest empty double wedge, it requires an enumeration, which takes O(N2) time. 12

bisector of the largest gap to the equator. Figure 3.3(c) shows that points 5 and 6 are moved away from the equator after the application of R1. <Insert Figure 3.3> It is noted that, after applying R1 to the unit sphere, the (1,0,0) and (-1,0,0) are the only points that can possibly lie on the plane x3 = 0. The application of R2 to the unit sphere serves to move these two points away from the equator. In addition, each rotation locally moves the points as far away from the equator as possible by moving the circular midpoint of the maximum gap to the equator. This algorithm runs in O(N) time, where N is the number of points in P. Since the procedure for finding R2 is the same as that for R1, only the latter is described here: Algorithm P: Initial Rotations Step 1. Parallelly project the points in P onto the plane x1 = 0. (This can be done in O(N) time.) Step 2. Find the maximum circular gap. Let pi = (XilXi2,xi3), i = 1,...,N, be the unit vectors in P and Pi' = (xi2,xi3) be the image of Pi in the plane x1 = 0. Let Pi' = (x12,X13) be the base vector and the semicircle considered be consisted of all the vectors having a positive cross product with p'. Let qi, i = 1,...,N, represents the degree of the angle between the two vectors Pl' and Pi' (or the its antipode (-xi2,-xi3) if it has a negative cross product with pi'). Then ql = 0. Let qN.1 be equal to c. Now Gonzalez's algorithm can be applied to qi, i = 1,...,N+1, to find the largest gap and its two end points, qj and qk. (This can be done in O(N) time.) Step 3. Compute the locally optimal rotation R1. Let qo be the angle between Pl' and (1,0), a unit vector on x3 = 0. If (1,0) has a positive cross product with pi', R1 consists of a rotation of-qo + (qj+qk)/2; otherwise, R1 consists of a rotation of qO + (qj+qk)/2. (This can be done in 0(1) time.) 13

f co 0 t_ I~ ca o c" -* 0 -.0 Wb

Step 4. Apply R1 to P. (This can be done in O(N) time.) 4. Algorithms Several basic algorithms for solving geometrical problems on the unit sphere are given. They are useful when applied to determining the orientations of the workpiece in NC machining. Problem S1 (Spherical): Detecting Hemisphericity. Given a set P of N points in S2, determine whether P is hemispherical (contained in an open hemisphere.) If P is hemispherical, find a great circle which determines a hemisphere that completely contains P. The detection of hemisphericity for a set P of N points in S2 is basic to spherical algorithms. Since central projection establishes a one-to-one relations between points in S2 and points in E2, if P is hemispherical, it is relatively straight forward to convert a spherical problem for P to a planar counterpart. On the other hand, if P is not hemispherical (or, more precisely, if P contains antipodal points), the convex hull for P is the entire sphere by the definition of convexity. Problem S1 can be solved in E2 through proper transformations. For simplicity, assume that none of the points in P lies on the plane X3 = 0. (If not, Algorithm P described in Section 3.3 can be applied to ensure this.) The points in P can be partitioned into two possibly empty sets P+, those points with X3 > 0, and P-, those with X3 < 0. If either of the two sets is empty, Problem S1 is trivial; P is hemispherical and is contained in a hemisphere determined by the great circle on X3 = 0. If both sets are not empty, P+ and Pare mapped to P'+ and P' in E2, respectively, through central projection. Planar algorithms can now be applied to find the convex hulls of the two planar sets, CH(P +) and CH(P'). Given CH(P +) and CH(P ), it is shown in 14

Theorem 4.1 that the set P is not hemispherical if and only if CH(P ) and CH(P ) intersects. Theorem 4.1 Given a set of N points P, P is not hemispherical if and only if CH(P +) and CH(P ) intersects, where CH(N) denotes the convex hull of a set of points in the euclidean plane and P', P' denote the images of the points on the X3 > 0, X3 < 0 hemispheres, respectively, through central projection. Proof. First, it is shown that, if CH(P +) and CH(P ) intersect, P cannot be hemispherical. There exists a point p' E CH(P +) n CH(P) which is the image of a pair of antipodal points p+ e SCH(P+) and p- E SCH(P-) as illustrated in Figure 4.(a). Since SCH(P) and SCH(P-) are the smallest convex sets containing P+ and P", respectively, any convex set containing P must contain both. Suppose H is a open hemisphere that contains P. Then H contains SCH(P+) and SCH(P-) which implies that it includes at lease a pair of antipodal points. This violates the definition of hemisphericity. Second, it is shown that, if CH(P+) and CH(P-) do not intersect, P is hemispherical. For, if they do not, there exists a separating line SL such that CH(P+) and CH(P') lie on opposite sides of it. Let the great circle corresponding to SL be L. Then S2 is partitioned into 4 regions by L and the equator as shown in Figure 4.1(b). Suppose P exists in the quadrant (L-, PD). If P+ exists in the quadrant (L+, P+), then CH(P+) and CH(P) are on the same side of SL by central projection. Since they are not, P+ can only exist in the quadrant (L-, P+). Hence, P must be hemispherical to the great circle L. Suppose P exists in the quadrant (L+, P-). Same arguments applies and P+ can only exist in the quadrant (L+, P+). Again, P must be hemispherical. O <Insert Figure 4.1> The algorithm for detecting hemisphericity is immediate. 15

tp (a) CHP) and CH(P-) intersect Figure 4.1 Hemisphericity

L',P') ( L-,P) / (L',P-) (L+,P-) (b) CH(P+) and CH(P-) do not intersect

Algorithm Sl: Hemisphericity Step 1. Construct P+ and P by applying the central projection. If either set is empty, the set P is hemispherical. (This can be done in O(N) time.) Step 2. Find CH(P +) and CH(P'). (This can be done in O(NlogN) time using Graham's algorithm [6] or other convex hull algorithms as described in [14].) Step 3. Test if CH(P +) and CH(PD) intersects. (This can be done in O(N) time using the algorithm proposed by O'Rourke et al in [11].) If so, the set P is not hemispherical. otherwise, go to Step 4. Step 4. P is hemispherical. Find a separating line for CH(P +) and CH(P'). (This can be done in O(N) time using Megiddo's algorithm as described in [9,16].) The plane contains the separating line and the origin determines the hemisphere. Problem S2: Constructing Spherical Convex Hull Given a set P of N points in S2, construct its (spherical) convex hull, SCH(P). By definition of convexity, if the set P is not hemispherical, SCH(P) is the entire sphere. Therefore, the first step is to detect whether P is hemispherical by applying Algorithm S1. Then, if P is hemispherical, SCH(P) can be constructed from the by-products of the hemisphericity detection algorithm - the two convex hulls, CH(P') and CH(P ). Algorithm S2: Spherical Convex Hull Step 1. Determine if P is hemispherical using Algorithm S1. (This can be done in O(NlogN) time.) If P is not hemispherical, SCH(P) is the entire sphere; otherwise, go to Step 2. 16

Step 2. P is hemispherical. Construct SCH(P+) from CH(P+) by connecting those points in S2 that give rise to the vertices of CH(P'+), in the order of tracing the boundary of CH(P +). (This can be done in O(N) time.) Similarly, construct SCH(P-) from CH(P-). Step 3. Let L = { x I (a.x) = 0, a = (ai,a2,a3) } be the great circle found at Step 1, which determines a hemisphere that completely contains P. Define a three dimensional rotation R which rotates a = (al,a2,a3) to (0,0,1). (This can be done in constant time.) Then, after SCH(P+) and SCH(P-) are transformed by R, both are completely contained in S2+, the hemisphere corresponds to X3 > 0. (This can be done in O(N) time.) Step 4. Construct CH(P"+) and CH(P"-) by applying central projection. (This requires O(N) time.) Step 5. Find CH(P"), the convex hull of the union of CH(P"+) and CH(P"'). (This can be done in O(N) time using Shamos' algorithm described in [14] or Preparata-Hong's in [12].) Construct SCH(P) from CH(P"). Figure 4.2 illustrates how SCH(P) can be constructed from the by-products of Algorithm S1. After SCH(P+) and SCH(P-) are transformed by the rotation defined in Step 3 of Algorithm S2, both are on the northern hemisphere (above the equator L). Through central projection, the two spherical convex hulls are mapped to two convex hulls, CH(P"+) and CH(P") in E2. Then, the convex hull of the union of CH(P"+) and CH(P"), the shaded region, is found by constructing the support lines. The complexity of this algorithm is O(NlogN), which is dominated by that of Step 1. < Insert Figure 4.2> 17

Figure 42 Constution o Spherical Convex Hul

Problem S3: Detecting Hemisphericity of a Union. Given a spherical convex polygon P with N vertices and a point p in S2, determine if the union of P and p is hemispherical. Problem S3 is fundamental, if a dynamic convex hull algorithm is required. If the union of P and p is found to be hemispherical, its convex hull can be constructed by using projection and planar algorithms; otherwise, its convex hull is the sphere. Without falling back to the O(NlogN) algorithm for Problem S1, S3 is converted to a convex polygon inclusion problem and a simple O(N) algorithm is presented. Algorithm S3: Union Hemisphericity Let the antipode of p be q. Determine whether q is internal to P. Since the (spherical) convex polygon P can be viewed as the intersection of the unit sphere S2 and all the half-spaces, each of which is determined by an edge of P and the origin, q in S2 is internal to P if it is internal to the intersection of all the half-spaces thus defined. Whether q lies in a half-space can be found by taking dot product of q and the normal of the plane defining the half-space. If q is internal to P, the union of P and p includes an antipodal pair, p and q, and, therefore, it is not hemispherical; otherwise, it is hemispherical. This algorithm requires O(N) time because the detection of whether a point lies in a half-space takes constant time and the polygon has N edges. Problem S4: Intersection of Hemispheres. Given N hemispheres H1, H2,...,HN, find their intersection. By observing the duality between a point and a hemisphere on the sphere, It is shown that Problem S4 is the dual of Problem S2. 18

Theorem 4.2 Suppose a hemisphere is defined as Hi = (x I (pi x) >0, Pi= (Pil,Pi2,Pi3) x E S2}, i = 1,...,N. Then, the intersection of hemispheres HI nH2 r"rnHN exists if and only if the set P of points pi, i = 1,...,N, is hemispherical. Proof. Suppose P is hemispherical. The intersection H1 n H2 n rn HN exists since, for any j,k, l~j<N, l~k<N, Hj n Hk is not empty. Suppose P is not hemispherical. Let P* = ( p I (pqj) > 0, p E P } where qj is the antipode of some pj E P. Then P* is the largest subset of P such that the union of P* and q is hemispherical. Since P is not hemispherical, P* is not empty. Let I* = (i I pi E P* ). The intersection of Hi, i E I*, exists since P* is hemispherical. Furthermore, this intersection lies within the open hemisphere Hj' = (x I (p-x) < O, x E S2 which has no intersection with Hj. Therefore, the intersection H1 n H2 n"'n HN is empty. O Theorem 4.3 Let Hi = (x I (pi-x) O0, pi = (pi1,Pi2pi3), x S2, i = 1,...,N, and P be the set of points pi, i = 1,...,N. If P is hemispherical, each vertex of the spherical convex hull SCH(P) corresponds to an edge of the intersection HI n H2 rn"' HN and vice versa. Proof. First, it is shown that the hemisphere determined by any point in SCH(P) that is not a vertex is redundant, i.e. it does not contribute to the boundary of the intersection H1r) H2 n n HN. Let Pk be any point in the interior of the convex hull and Hk = ( x I (pk.x) > 0, x E S2 } be the hemisphere defined by Pk. Suppose Hk contributes to an edge of the intersection. Let the edge be decided by the intersection of Hk and two other hemispheres, Hi and Hj. Since Pk is in the interior of SCH(P), there exists another point Pm E P in the neighborhood of Pk such that the distance between Pm and the edge contributed by H is greater than i/2. Then, the intersection of the hemisphere determined by Pm, Hi, and Hj is smaller than that of Hk, Hi, and Hj. Therefore, Hk does not contribute to the boundary. Now, let pk be any point on an edge of SCH(P) and let pi and pj be the two end points of the edge. The intersection of the hemispheres 19

determined by Pk, Pi, and pj has only two vertices (since they are collinear) and two edges (from pi and pj). Therefore, the hemisphere determined by pk is redundant. Second, it is shown that the hemisphere determined by a vertex of SCH(P) is non-redundant. Let Pk be a vertex of SCH(P) and Hk be the hemisphere determined by Pk. Suppose Hk is redundant. Then, there does not exist x = (xl,x2,x3) such that Pilx1 + Pi2x2 + Pi3X3 > 0, 1li<N, i#k Pklxl + Pk2X2 + Pk3X3 < 0 By Farkas' Lemma [10], there exists some y = (yl,...,Yk-ltYk+l~- ayn) such that P11Y1 +'"+Pk-l,lYk-l + Pk+l,lYk+l + + PN1YN = Pkl P12Y1 +"'+Pk-1,2Yk-1 + Pk+l,2Yk+l + "+ PN2YN = Pk2 P13Y1 +"'+Pk-1,3Yk-1 + Pk+l,3Yk+l + + PN3YN = Pk3 Yl,...,Yk-1,Yk+1,*-.,YN 0 Therefore, Pk can be written as a non-negative linear combination of pi, 1i<N, ink. This contradicts the fact that Pk is a vertex. O Figure 4.6 illustrates the duality between the convex hull of Pi and the intersection of Hi. Each edge of the intersection corresponds to a vertex of the convex hull. Based on Theorem 4.2 and Theorem 4.3, an O(NlogN) algorithm is proposed here for computing the intersection of hemispheres. < Insert Figure 4.3> Algorithm S4: Intersection of Hemispheres Step 1. Determine if P is hemispherical. (This requires O(NlogN) time using Algorithm S1.) If P is not hemispherical, the intersection H1 n H2 rn"" HN is empty; otherwise, go to Step 2. Step 2. Find SCH(P). (This can be done in O(N) time using Algorithm S2.) 20

PI Figure 4.3 Duality between Convex Hull and Intersection

Step 3. Find a point p in the interior of SCH(P). By Theorem 4.3, this point determines a hemisphere that completely contains the intersection H1 n H2 n"'nHN. (This can be done in O(N) time.) Step 4. Traverse SCH(P) and construct the intersection H1 n H2 n'rv HN. For each edge of SCH(P), the intersection of the two corresponding hemispheres determines a pair of antipodal points. One of the pair, which has a positive dot product with p (found at Step 3), is a vertex of H1 n H2 nn"' HN. (This step can be done in O(N) time since the intersection of two hemispheres can be computed in constant time by taking the cross-product of the two unit normals.) 5. Conclusions Unifying two fields of research - one that is analytically complex (but combinatorially simple) and the other combinatorially complex (but analytically simple) - may result in one of two outcomes. If an attempt had been made via the "lowest common denominator" of the two fields, e.g. the machining of planar surfaces, the outcome would have been predictably trivial. This paper reflects an attempt to combine the more appealing aspects of both fields by way of spherical geometry. To avoid duplication of and to give recognition to past work, central projection is adopted in this paper. In computational efficiency, simplicity in structure offers many advantages. Convexity has been employed as the underlying theme in this paper in the form of hemisphericity and its related computations. Hemisphericity is justified by the practical application of NC machining, and, in particular, the geometry of the cutting tools and the degree of freedom in 3-, 4-, 5-axis milling machines. 21

Acknowledgement The authors are grateful to the support of a grant from the Science Research Laboratory, Ford Motor Company. The authors also sincerely appreciates helpful insights and comments from their colleagues: J. Gan, D. S. Kim, S. Y. Chou, K. Tang, J. Y. Park, W. J. Lee, S. Y. Shin, J. Wolter, and D. Dutta. References [1] D. P. Anderson, Efficient algorithms for automatic viewer orientation, Comput. & Graphics, vol. 9, no. 4,407-413, 1985. [2] K. Q. Brown, Geometric transformations for fast geometric algorithms, Ph. D Thesis, Dept. of Computer Science, Carnegie Mellon Univ., Dec. 1979. [3] B. M. Chazelle, L. J. Guibas, and D. T. Lee, The power of geometric duality, BIT 25, 76-90, 1985. [4] A. R. Forrest, Computational Geometry, Proc. Royal Society London, 321 Series 4, 187-195, 1971. [5] T. Gonzalez, Algorithms on sets and related problems, Technical Report, Dept. of Computer Science, Univ. of Oklahoma, 1975. [6] R. L. Graham, An efficient algorithm for determining the convex hull of a finite planar set, Info. Proc. Lett. 1, 132-133, 1972. [7] D. Hilbert and S. Cohn-Vossen, Geometry and the Imagination, New York: Chelsea, 1952. [8] B. K Horn, Extended gaussian image, Proc. of the IEEE, vol. 72, no. 12, 1671-1686, Dec. 1984. [9] N. Megiddo, Linear time algorithm for linear programming in R3 and related problems, SIAM J. of Comput. 12(4), 759-776, Nov. 1983. [10] K. G. Murty, Linear Programming, New York: Wiley, 1983. [11] J. O'Rourke, C. Chien, T. Olson, and D. Naddor, A new linear algorithm for intersecting convex polygons, Computer Graphics and Image Processing 19, 384-391, 1982. 22

[12] F. P. Preparata and S. J. Hong, Convex hulls of finite sets of points in two and three dimensions, Comm. ACM 2(20), 87-93, Feb. 1977. [13] F. P. Preparata and D. E. Muller, Finding the intersection of n half-spaces in time O(nlogn), Theoretical Computer Science 8(1), 45-55, Jan. 1979. [14] F. P. Preparata and M. I. Shamos, Computational Geometry - An Introduction, Springer Verlag, 1985. [15] P. J. Ryan, Euclidean and Non-Euclidean Geometry - An Analytic Approach, New York: Cambridge University Press, 1986. [16] M. I. Shamos, Computational Geometry, Ph.D. Thesis, Dept. of Computer Science, Yale Univ., 1978. [17] S. Y. Shin, Visibility in the Plane and Its Related Problems, Ph.D. Thesis, Dept. of Industrial and Operational Engineering, Univ. of Michigan, Ann Arbor, 1986. [18] J. Stolfi, Primitives for Computational Geometry, Digital Systems Research Center Report 36, 1989. 23

UNIVERSITY OF MICHIGAN 3 9015 04732 7575