PAR: A New Representation Scheme for Rotational Parts * Yung-Chia Lee and Kuen-Fang Jack Jea Department of Electrical Engineering and Computer Science The University of Michigan, Ann Arbor, MI 48109-1109 July 28, 1986 Abstract A new representation scheme PAR (Principal Axis Representation) for rotational parts is proposed as an internal representation scheme for Constructive Solid Geometry (CSG). The key idea of PAR is to represent an object by its principal axis and a set of boundary curves. Based on a mathematical.framework, an algorithm is designed to convert a CSG tree into a PAR, which represents the -same object as the CSG tree does but is in an evaluated form. Geometrical properties of parts can then be computed more directly and efficiently from this evaluated representation than from the original CSG tree. In addition to its computational efficiency, the PAR is a unique representation scheme. This work is supported in part by the NSF Grant No. DMC-8504790 and by the SME Grant No. 585-895. i

Contents 1 Introduction 1 2 The Problem and Proposed Approach 3 2.1 Problem...................................... 3 2.2 Basic Ideas................................... 4 3 Formulation of Principal Axis Rcpresentation 5 3.1 Terminology.................................... 5 3.2 Represent Primitive Solids by PAR....................... 8 3.3 Operations on PAR................................ 9 3.4 Relation between PAR and CSG........................ 12 4 Algorithm to Convert CSG Representation to PAR 17 4.1 Data Structure for PAR............................. 17 4.2 Conversion from CSG Representation to PAR................. 17 4.3 Union and Difference of Two PAR's...................... 19 5 Illustrative Examples 27 6 Conclusion and Discussion 29 6.1 Conclusion.................................... 29 6.2 Discussion and Future Work........................... 30

RSD- TR-14-86 1 Introduction Solid modeling techniques have been considered as one of the keys to the integration of computer-aided design and manufacturing (CAD/CAM)[Bee821. Two of the prevailing geometric schemes in solid modeling are Boundary Representation (B-rep) scheme and Conetructive Solid Geometry (CSG) scheme[ReV82]. B-rep represents an object by segmenting its boundary into a finite number of bounded subsets usually called faces or patches, and representing each face by (for example) its bounding edges and vertices. CSG scheme represents objects as constructions or combinations, via the regularized set operators, of solid components. While B-rep in general is in an evaluated form, which explicitly shows what and where each geometric entity is, CSG is not. The much more concise CSG is easy to construct and to check the validity of the object being constructed. The tradeoff between the evaluatedness and conciseness seems to be clear but the choice between B-rep and CSG has never been easy. As a means to integrate design and manufacturing, solid modeling techniques not only display objects but also compute important properties of objects. A more comprehensive understanding of the geometric representation by the computer becomes so important that shape features can then be extracted, parts classified, processes planned, and etc.. Several interesting researches have been done along the direction of enhancing machine's understanding of part geometric representations. Woo[Woo77] focused on cavity recognition in studying problems of transforming volumetric designs of parts into numerical control (NC) descriptions. Staley et al.[StHA831, without specific application emphasized, also dealt with cavity recognition but using syntactic pattern recognition techniques. Grayer(Gra761] compares the part in B-rep with its initial work part in order to compute NC tool paths. In addition to NC path generation, Armstrong[ArCP84] also considered the determination of fixture orientation. For global understanding of part geometry, Woo[Woo821 suggested a convex hull technique which transforms a boundary representation into an expression of alternating sum of volumes. To generate a part code for group technology, Kyprianou(Kyp80] used syntactic pattern recognition techniques to extract features characterized by protruA New Representation Scheme for Rotational Parts 1

RSD- TR-14-86 sions and depressions. Jakubowski[Jak82] also used syntactic methods to develop a part description language primarily for the determination of part shapes. Henderson(Hen84] was able to use PROLOG and expert system techniques to extract features and organize them in a feature graph so as to facilitate automatic process planning. While most of the work mentioned above is based on B-rep, very little has been done on CSG. In CSG, a machine part is represented as a tree with some primitive solids as terminal nodes and the regularized set operations (union and difference) and movement operations as non-terminal nodes. Properties of the machine part can be computed by evaluating the CSG tree. The evaluation procedure is basically a tree traversal routine that combines CSG subtrees into higher level CSG subtrees. The algorithm is simple but the computation involved in calculating the intersection (the major computation of union and difference operations) of two solids are complex and time-consuming. Without evaluation, Lee and Fu[LeK861 proposed an approach to the problem of feature extraction. The approach is based on the spatial relationships among axes of primitive solids. They also suggested a unification process using tree reconstruction to unify feature representations so that properties can be computed more easily. As more general cases for this approach are yet to be studied, there are algorithms that convert CSG into Brep[BoG&2]. Not surprisingly, such algorithms are often called for when certain properties must be computed from a CSG tree[KaO841. However, the conversion often appears to be computationally expensive or, being worse, a waste to simply compute the specific property. There are part families for which CSG should be an excellent representation scheme and, furthermore, a transformation into B-rep can actually be avioded. A such example, rotational parts, is investigated in this study. Instead of B-rep, a new representation scheme PAR (Principal Axis Representation) is proposed as an internal representation scheme for rotational parts that are initially described by CSG. The key idea of PAR is to represent an object by a principal axis and a set of boundary curves. Based on a mathematical framework, an algorithm is designed and implemented to convert a CSG tree into a PAR, which represents the same object as the CSG tree does but is in an evaluated form.'rom PAR, the profile of the part can be efficiently computed. Other properties such as length 2 A New Representation Scheme for Rotational Parts

RSD. TR-14-86 and maximum diameter can also be obtained easily. The conversion algorithm is much simpler than the one converting CSG into B-rep. The main reason is that one dimensional curves (e.g. lines and arcs) are used to characterize the rotational parts and the evaluation of PAR is performed on two dimensional space, e.g. the intersection of two lines on the same plane rather than that of two solids on a three dimensional space. In addition to its computational ease, the PAR is proved to be a unique representation scheme. Being unique, it facilitates feature definition and thus simplifies the tasks of feature extraction. The possible extension of PAR is discussed in Section 6. 2 The Problem and Proposed Approach 2.1 Problem Rotational parts include all parts that are symmetrical with respect to their principal axes. In this study, however, the primitive solids are limited to cylinders, cones, and tori only. The problem is formulated as follows: Given a CSG tree of a machine part, which is 1. axis-symmetrical, and 2. constructed from cylinders, cones and tori in an arbitrary order of combination such as union and difference. can we: 1. by utilizing the property of axis-symmetry, efficiently derive its profile and some geometric properties such as length and diameter so that part classification and process plannin g can possibly be supported? 2. develop an internal scheme which not only supports the above computation but also represents the part uniquely? A New Representation Scheme for Rotational Parts 3

RSD- TR-14-86 2.2 Basic Ideas If a machine part is axis-symmetrical, it can be viewed as a 2iD object and represented by the union of components, each being generated by rotating its boundary curve segments with respect to the principal axis. For example, a cylinder can be represented by rotating a line segment around its axis where the axis and the line segment are parallel to each other and their distance is the radius of the cylinder. Other primitive solids like cones and tori can also be represented in a similar way. Since the machine part is axis-symmetrical, each component solid depicted by a CSG subtree can also be represented by a collection of principal axis segments, each being associated with several pairs of boundary curves. For each boundary curve pair, rotating it around the associated axis segment would generate a volume layer for the corresponding component solid. This idea of representing axis-symmetrical machine parts by a principal axis and a set of bounded curves can be exploited to develop efficient algorithms to evaluate CSG trees and deduce their geometrical properties because the computational complexity of the algorithms can be reduced significantly. For example, in the course of evaluating CSG trees, instead of testing the intersection of two three dimensional solids, we need only to test the intersection of a few one dimensional curves on the same plane, where the number of the ID curves needed to be tested depends on the complexity of the two original 3D solids involved. To evaluate a CSG tree and convert it into this kind of representation, new operations such as union, difference, and movement operated on this new representation scheme must be defined in order to maintain the actual semantics as their counterpart operations do in the CSG tree representation scheme. To deal with the movement operation is easy; if all the coordinates of the segments and boundary curves are relative to a principal axis, then the result of the movement operation is merely applying the movement transformation to the principal axis, and others remain unchanged. But the union and difference operations are not so simple to deal with, as more elaborate processing would be required. After a CSG tree is converted to the new representation, the profile and other geometric properties of the corresponding part can be easily computed. These properties may be 4 A New Representation Scheme for Rotational Parts

RSD. TR-14-86 computed all at the same time after the CSG tree is converted into the new representation. To compute the desired geometric properties directly from the CSG tree, on the other hand, is complex and time-consuming. In summary, the proposed approach includes the following: 1. To define a new representation scheme (Principal Axis Representation) for axissymmetrical machine parts 2. To convert CSG trees into their new representations. 3. To derive desired geometrical properties from this new representation, and it is believed to be easier and more efficient. 4. To prove that the Principal Axis Representation is a unique representation scheme for axis-symmetrical machine parts. 3 Formulation of Principal Axis Representation 3.1 Terminology Definition: Principal Axis A Principal Axis A of an axis-symmetrical object O is a line segment in the three dimensional space such that the object 0 is represented in a way that starts from and ends at the two end points of A, and 0 is symmetrical with respect to A. Definition: Axis Segment An Azis Segment S is a line segment ( i.e. subset ) of a Principal Axis which has two end points. Definition: Principal Axis Coordinate System A curve C and a principal axis A form a coordinate system if they are co-plane and, on this plane, the horizontal axis is the line containing A and its vertical axis is the line A New Representation Scheme for Rotational Parts 5

RSD- TR-14-86 perpendicular to A and passing through the start point of A. This 2-D coordinate system is called principal azi. coordinate system. Definition: Bounded Curve A curve C is bounded at [a bJ with respect to a principal axis A iff in the principal axis coordinate system formed from C and A, there exist two real numbers a and b ( a:b ) such that the curve C is differentiable within the interval (a b). a and b are called the bound points of the curve C, and [a b] is called a pair of bounds of C. Definition: Range of a Bounded Curve If C is a bounded curve within [a b] with respect to a principal axis A, the range of C at x, x E [a b], on the principal axis coordinate system is C(x), i.e. the mapping of C from x. Definition: Arc An Arc is a bounded curve which is also a subset of a circle. Note that: a line segment is a bounded curve. Definition: Bounded Curve Set A Bounded Curve Set is a set of Arca and/or Line Segments. Definition: Principal Axis Representation ( PAR ) If an object O is an axis-symmetrical machine part which can be characterized by a Principal Axis A, then its Principal Azia Representation PAR(O, A) can be defined as a set of tuples (Si, Ci), i = 1,.., n, n E N, such that 1. Si is an Axis Segment, and Si 5 A, U=lI Si = A, Hi, j ~ n, SiUSi is either a point (i.e. their common bound point) if j = i + 1 or 0, otherwise; that is, { Si } is a partition of A. 2. Ci = ( Cii Ij= 1,.., 2mi ) is a Bounded Curve Set. 6 A New Representation Scheme for Rotational Parts

RSD- TR-14-86 3. 3a, b E R such that all Cii's and Si have the same bound [a bJ with respect to A. 4. all Cii's do not intersect within their bound interval (a b). 5. all Ci's are differentiable (i.e. their first derivatives exist ) within (a b), and 3 a curve Ck E C, such that Ck is not differentiable at the bound point a, and 3 a curve C,, E Ci such that Cil is not differentiable at the bound point b. Theorem 1: For a Principal Axis Representation of an axis-symmetrical machine part O and its Principal Axis A, PAR(O, A) = { (Si, CO) I i=l,..,n }. Vi < n, all Cii' ( Cii E C, j = 1,..,2ri ), within (a b) which is the bound of Si with respect to A, form a total ordering. Proof: let's denote Cii(x) to be the range of the bounded curve C,, at x, x E (a b). For Cik and Cim E Ci, k # m, 3z E(a b), Ci,(x) and Ci,(x) E R, and Ci,(x) # Cim(x); otherwise Ck and Cim intersect at x within (a b) and this contradicts to the definition of PAR(O, A). Therefore either Cik(x) < Cim(x) or Cik(x) > Cim(x) is true. And we prove that all Cii(x), Cii E Ci, form a total ordering at some x, x E (a b). Next, we will prove that the total orderings of all Cii (j = 1, 2,.., 2ri) at all x, x E (a b) are consistent. For x1 and ZX2 E (a b), x1 # 72, if the total ordering of all Cii (j = 1, 2,.., 2ri) at x1 and that of all Cii (j = 1, 2,.., 2ri) at 22 are not consistent, then there must exist two different curves, Cik and Cil, such that Cik(zl) < Ci,(zl) and Cik(z2) > Cig(z2). It means that there must exist some value X3 between zl and Z2 such that Cik and Cil intersect at z.X This contradicts to the definition of PAR. Therefore the total orderings of all Cii (j = 1, 2,.., 2ri) at all x, x E (a b), must be consistent, and we may conclude that all Cii (j = 1, 2,.., 2ri) within (a b) form a unique total ordering. Q.E.D. Note: this total ordering will be used to determine the relation among those curves within some bound. Using the concept of Theorem 1, we may define the concept of a layer. Definition: Layer and Layer Set A New Representation Scheme for Rotational Parts 7

RSD. TR-14-86 For PAR(O, A), Si is an Axis Segment and C, is its associated bounded curve set. Cii's (j =1,..,2mi) are those bounded curves in C, and they are sorted in such a sequence: Ci1 > Ci2 > Ci3 >.. > Ci2mi. Then a layer within this Axis Segment S is defined as a pair of curves C;* and Ck+1,denoted as [CikCik+ll, ki =1, 3,..., 2mi-1. The set of layers { [CikC;,+1] I k =1, 3,..., 2mi-1 } is called a layer set. Note: a layer with its bound actually generates a volume of the object by rotating the curves in the layer with respect to its corresponding Axis Segment. 3.2 Represent Primitive Solids by PAR 1. Use PAR to represent cylinder PAR( cylinder, A ) = { (S, C) C = { C1,C2 }, C1 and C2 are line segments; S = A = C2; C1 II C2; I C1 - C2 I = radius of the cylinder; [a b] is the bound for S and C where a and b are the horizontal coordinates of the starting and ending points, respectively, of A in the Principal Axis Coordinate System formed by A and C1. } Notice that the statement'"C II C2' means that the line segment C1 is parallel to the line segment C2, and the statement "I C1 - C2 I" denotes the distance between C1 and C2. 2. Use PAR to represent cone PAR( cone, A ) = { (S, C) I C = { C1,C2 }, C1 and C2 are line segment; S = A = C2; [a b] is the bound for S and C where a and b are the horizontal coordinates of the starting and ending points, respectively, of A in the Principal Axis Coordinate System formed by A and C1; C1 and C2 intersect at b, and the distance of C1 and C2 at a is the radius of the cone. } 3. Use PAR to represent torus PAR( torus, A) = (S, C) I S = A; C = { C1,C2 }, [a b] is the bound for S and C where a and b are the horizontal coordinates of the starting and ending points, respectively, of A in the Principal Axis Coordinate System formed by A and C; C1,C2 are Arcs of the circle with the radius equal to the local radius of the torus and with its center located at (x 8 A New Representation Scheme for Rotational Parts

RSD- TR-14-86 y) in the Principal Axis Coordinate System where z = (a + b)/2, and y is the global radius of the torus, C1 is the upper half and C2 is the lower half of the circle. } 3.3 Operations on PAR Definition: Overlap of Two Layers The overlap of two layers [clelc2] and [c2lc"J with respect to some bound [a bJ is defined as either a layer [c3lc321 with respect to [a bi if C31 = min(cll, C2l), c32 = maX(c12, C22) and c31 > c32; or 0, otherwise. Definition: Layer-Union of Two Layers The layer-union of two layers [Icl C21 and [Cen2 c22] with respect to some bound [a b] is defined as either the layer [C31C321 where cs1 = max(cll, C21), C32 = min(cl2, C22) if the overlap of [Celll2] and [c2212] is not 0; or the two layers [eCllcl2] and [c21c22], otherwise. The result layer or layers are with respect to the same bound[a b]. Definition: Layer-Union of Two Layer Sets The layer-union of two layer sets L1 and L2 is defined as a layer set which is the set union of all layer-union's of two layers, one from L1 and the other from L2, respectively. L1, L2 and the resulting layer set are all with respect to the same bound. Definition: Maximum Union-able Layer Set The maximum union-able layer set of two layer sets L1 and L2 is defined as a layer set which is the transitive closure of applying the layer-union operation to the resulting layers of the layer-union of L1 and L2; that is, if L3 is the layer-union of L1 and L2, then the maximum union-able layer set of L1 and L2 is the results of repeatedly applying the layer-union operation to the layers in Ls until the resulting set is no more changed. L1, L2 and the maximum union-able layer set are all with respect to the same bound. Definition: Union of Two PAR's A New Representation Scheme for Rotational Parts 9

RSD- TR-14-86 For two PAR's, PAR(O1, Al) and PAR(02, A2), if Al and A2 are subsets of a line, then the union of PAR(O1, Al) and PAR(02, A2) is defined as a mapping from PAR x PAR to PAR such that if PAR(O1, Al) = { (Si, C,') I i = 1,..,nl, [alb!] is the bound of S,' }, PAR(02, A2) = { (S2, C) I i = l,..,Ln2, [ab,2 is the bound of Sf2 }, and there exists a PAR(03, A3s){ (S, Ci3) = l1,..,n, [ab?]J is the bound of S$3 }, then 1. 03 = 01 U 02; As= A UA2; 2. S3is a line segment of A3 with bounds [a0]b3, i=l,..,n; and there exists a set of bounds de, m = 1, 2,..,n3 (ns n) which is a partition of the set of bounds [a4b? and there exists a bounded curve set K,3 for each bound [ds e] where K3 = Cil with adjusted bound [d| e" if a] [a2bl] D [de] = 1,..,n2 K,3' = Ck with adjusted bound [d3.e? if a [abl] 2 [de3, 1= l,, K - = maximum union-able layer set of CI and Ca with respect to a bound [d3 if [aWb] ) [d. es ], and [ab] 2 [d e3s ] for some j, k. 3. C3-= Uli=3. KS,(x;, ys are integers) if [d e3 for m = zx,..., yi is a partition of [albJ NOTE: This mapping is from PAR x PAR to PAR, the non-differentiable properties of curves across the bound l[abS1 for the resulting PAR should be preserved by the definition of PAR. Definition: Layer-Difference of Two Layers The layer-difference of two layers [CllCl21 and [c21c22] with respect to some bound [a b] is defined as 1. the layer [elle12l if cll < c22 or C12 > C21; 2. the layer [C22c12] if c21 > cll > c22 > cl2; 3. the layer [CllC2l] if eC > C21 ~ C12 ~ c22; 4. the layers [cllc21] and [c22c12] if Cll > C21 > C22 > ~12; 10 A New Representation Scheme for Rotational Parts

RSD- TR-14-86 5. 0 if C21 > CH > cll> 2 > c; The result layer or layers is with respect to the same bound[a b]. Definition: Layer-Difference of Two Layer Sets The layer-difference of two layer sets LI and L2 (assume L1 = { [cic;+1l I i = 1, 3,..., 2n-1 } is the first operand and L2 = { [didj+l] I i = 1, 3,..., 2m-1 } is the second operand) is defined as a layer set which is the aet union of all the set OVk (k = 1, 2,..., n) where each OVk is the overlap of all the layers 1ki (j = 1, 2,..., m), lkj = the layer-difference of [c2k-lcc2k to [d2j-ld2j]. L1, L2 and the resulting layer set are all with respect to the same bound. Definition: Minimum Difference-able Layer Set The minimum difference-able layer set of two layer sets L1 and L2 with respect to some bound [a bJ is the layer-difference of LI to L2 with respect to the same bound [a b]. Definition: Difference of Two PAR's For two PAR's, PAR(01, Al) and PAR(02, A2), if Al and A2 are subsets of a line, then the difference of PAR(01, Al) and PAR(02, A2) is defined as a mapping from PAR x PAR to PAR such that if PAR(O1, Al) = { (S,, C1) I i = 1,..,n, [alb] is the bound of S }, PAR(02, A2) = { (S2, C2) i = 1,..,n2, [ab2] is the bound of S2 }, and there exists a PAR(03, As) = ( (S3, C) I i = 1,.,n, [ab] is the bound of S,3 }, then 1. 03 = 01- 02;A3 Al; 2. S,3 is a line segment of A3 with bounds [a0bI, i=l,..,n; and there exists a set of bounds I[dses 1 M = 1, 2,..,n3 (ns > n) which is a partition of the set of bounds [a3b?1 and there exists a bounded curve set 4,. for each bound [d3 e3 where K - = C} with adjusted bound [am,.e3 ] if ]! [a1b?] _ [de], I= 1,..,n2 Km = minimun difference-able layer set of Cin and C, with respect to a bound [d? e] if [a=b'] D d3 e3], andf [ datbb] _,e] for some j, k. A New Representation Scheme for Rotational Parts 11

RSD- TR-14-86 3. Cs = U=i K, (zi, y are integers) if (dses, for m = zXi,..., y, is a partition of Definition: Movement of a PAR The movement of a PAR(O, A) is defined as a mapping from PAR to PAR such that the new PAR(O,B) after movement is the same as PAR(O,A) except that A is transformed to B in 3-D space by applying the movement operation specified. Note: The bounds in PAR are with respect to A,and therefore movement-invariant. 3.4 Relation between PAR and CSG In this subsection, the relation between CSG and PAR representation schemes is investigated. It is the basis of our algorithm to convert a CSG tree into a PAR. Theorem 2. For a CSG tree with union as its root, the object represented by this CSG tree is identical to that represented by a PAR which is the union of two PAR's, each corresponding to either the left- or right- subtree of the original CSG tree. Proof: Assume that 03 = 01 U 02 in CSG domain where 03 is the CSG tree, and 01 and 02 are its left- and right- subtrees, respectively. In the PAR domain, assume 01 and 02 are represented by PAR(01,A1) and PAR(02,A2), respectively. Also assume A3 = A1 UA2. Now we want to prove PAR(01,A1)UPAR(02,A2) represents 03. In this proof, we may view an object as a set of points in 3-D space and use the notation V( layers ) to denote the volume generated by rotating the layers with respect to a principal axis in 3-D space. Step 1: To prove 03 5 V( layers of (PAR(O1,A1) U PAR(02,A2)) ). For every point p of 03, p E 01 or p E 02 is true in CSG domain because 03 = 01U 02. But 01 and 02 can be represented by PAR(O1,A1) and PAR(02,A2), respectively. Therefore, p E V( layers of (PAR(O1,A1) ) or p E V( layers of (PAR(02,A2) ). More specifically, on the principal axis coordinate system, the horizontal coordinate of p must lie 12 A New Representation Scheme for Rotational Parts

RSD- TR-14-86 5. 0 if C21 > Cll > C12 2 c22; The result layer or layers is with respect to the same bound[a b]. Definition: Layer-Difference of Two Layer Sets The layer-difference of two layer aeta L1 and L2 (assume L1 = { [cici+1] I i = 1, 3,..., 2n-1 } is the first operand and L2 = { [did+l] I i = 1, 3,..., 2m-1 } is the second operand) is defined as a layer set which is the set union of all the set OVk (k = 1, 2,..., n) where each OVk is the overlap of all the layers lki (i = 1, 2,..., m), 1ki = the layer-difference of (c2k-1C2k] to [d2ild2j]. L1, L2 and the resulting layer set are all with respect to the same bound. Definition: Minimum Difference-able Layer Set The minimum difference-able layer set of two layer sets L1 and L2 with respect to some bound [a b] is the layer-difference of L1 to L2 with respect to the same bound [a b]. Definition: Difference of Two PAR's For two PAR's, PAR(O1, Al) and PAR(02, A2), if Al and A2 are subsets of a line, then the difference of PAR(O1, Al) and PAR(02, A2) is defined as a mapping from PAR x PAR to PAR such that if PAR(O1, Al) = { (S,1, Cl) I i = 1,..,n, [alb] is the bound of S }, PAR(02, A2) = { (S2, C2)I i = 1,..,n2, [a2b2] is the bound of S2 }, and there exists a PAR(03, As) = { (S3, C3) i 1,..,n, [a0b?] is the bound of S3 }, then 1. 03 = 1- 02; A3 5 Al; 2. SV is a line segment of As with bounds [ab, i=l1,..,n; and there exists a set of bounds de, m =i 1, 2,..,n3 (n3 ~ n) which is a partition of the set of bounds [a0bI and there exists a bounded curve set Km, for each bound [ds e]3 where K - = C} with adjusted bound [d me] Jif [a] D [de], =,..,n Km = 0 if v! [abu] _D [dSme ], 1- l,..,nx Ks = minimun difference-able layer set of C}' and C with respect to a bound [d:e] 3 if [abl] _ [d3 e3 ], and [ab2] D >[de3] for some j, k. A New Representation Scheme for Rotational Parts 11

RSD- TR-14-86 3. - U=r K.m, (Zx, wi are integers) if [dse, ] for m = x,,..., y is a partition of [ab3] 1 Definition: Movement of a PAR The movement of a PAR(O, A) is defined as a mapping from PAR to PAR such that the new PAR(O,B) after movement is the same as PAR(O,A) except that A is transformed to B in 3-D space by applying the movement operation specified. Note: The bounds in PAR are with respect to A,and therefore movement-invariant. 3.4 Relation between PAR and CSG In this subsection, the relation between CSG and PAR representation schemes is investigated. It is the basis of our algorithm to convert a CSG tree into a PAR. Theorem 2. For a CSG tree with union as its root, the object represented by this CSG tree is identical to that represented by a PAR which is the union of two PAR's, each corresponding to either the left- or right- subtree of the original CSG tree. Proof: Assume that 03 = 01 U 02 in CSG domain where 03 is the CSG tree, and 01 and 02 are its left- and right- subtrees, respectively. In the PAR domain, assume 01 and 02 are represented by PAR(01,A1) and PAR(02,A2), respectively. Also assume As = A U A2. Now we want to prove PAR(O,A1)UPAR(02,A2) represents 03. In this proof, we may view an object as a set of points in 3-D space and use the notation V( layers ) to denote the volume generated by rotating the layers with respect to a principal axis in 3-D space. Step 1: To prove 03 C V( layers of (PAR(O1,A1) U PAR(02,A2)) ). For every point p of 03, P E 01 or p E 02 is true in CSG domain because 03 = 01 U 02. But 01 and 02 can be represented by PAR(O0,A1) and PAR(02,A2), respectively. Therefore, p E V( layers of (PAR(O1,A1) ) or p E V( layers of (PAR(02,A2) ). More specifically, on the principal axis coordinate system, the horizontal coordinate of p must lie 12 A New Representation Scheme for Rotational Parts

RSD- TR-14-36 in some bounds [aibil, which is a subset of [aJbA] of PAR(Ol,AI) or a subset of [alb2] of PAR(02,A2) for some j or k. 1. if [aibi is a subset of [atbAJ but not a subset of [abJ2l, then p E V( layers specified by Cl of PAR(O,A) ). 2. if (aibi, is a subset of [amb2I but no t a subset of [ab, then p E V( layers specified by Cat of PAR(02,A2) ). 3. if [aibl is both a subset of [a2b6 and a subset of [abJl, then p E V( [CmCm+il) or p E V( [dndn,+l] ) where [Cmcm+lJ and [dndn+lJ are some layers of PAR(O1,A1) and PAR(02,A2), respectively. By the definition of layer-union of two layers and Maximum Union-able Lager Set, p E V( layer-union of [cmc+ll and [dd,,+l ), and therefore p E V( maximum union-able layer set of LI and L2 ) where L1 and L2 are the layer sets with the bound [aibi] that contain [cmcm+i] and [dnd4,+], respectively. Thus, from the definition of the union of two PAR's, p E V( layers of (PAR(O1,A1) U PAR(02,A2)) ). And we prove that 03s 5 V( layers of (PAR(01,A1) U PAR(02,A2))). Step 2: To prove 0s 7 V( layers of (PAR(Oi,AI) U PAR(02,A2)) ). For every point p, p E V( layers of (PAR(O1,A1) U PAR(02,A2)) ), by the definition of the union of two PAR's, there are three possible cases for p: 1. p E V( [cmCm+iJ ) where [cmcm+l], which is bounded by some bound [a b] C [a'blJ but [al2b2], is a layer of PAR(O1,A1). In this case, since PAR(O1,A1) represents 01, p E 01 and therefore p E 01 U 02 = 03. 2. p E V( [dndn+1]. ) where [dnd,,+l], which is bounded by some bound [a b] C [a2b] but X [albfl, is a layer of PAR(02,A2). In this case, since PAR(02,A2) represents 02, p E 02 and therefore p E 01 U 02 = 0.33. p A V( maximum union-able layer set of L1 and L2 ) where L1 and L2 are the sets of layers of PAR(O,A1) and PAR(02,A2) at some bounds [a b], respectively. By the definition of the maximum union-able layer set, there are three possible cases for p: A New Representation Scheme for Rotational Parts 13

RSD.-TR-14-86 (a) p E V( [emcm+l] _ L1 but L2). Thus, p E 01 C (01U 02)= 03. (b) p E V( [dd,+l] C L2 but X L1). Thus, p E 02 C (Q1 U 02) = Os. (c) p E V( [ele+l]1 ) where [el el+,] C (L1 nl 2). In this case, p E (V(L1) n V(L2)). That is p E (O1 n 02) C (01 U 02) = Os. Therefore we prove that 03 2 V( layers of (PAR(O1,AI) U PAR(02,A2)) ). From Step I and Step 2, we prove the theorem. Q.E.D. Notice that the union operator on PAR performs the actual semantics of the union operator on CSG. Theorem 3. For a CSG tree with difference as its root, the object represented by this CSG tree is identical to that represented by a PAR vhich is the difference of two PAR's, which correspond to the left- and right- subtrees of the original CSG tree. Proof: Instead of using the concept of mazimum union-able layer set and the definition of the union of two PAR's, the concept of minimum difference-able layer set and the definition of the difference of two PAR's may be used here to prove this theorem in a similar way as in Theorem 2. Q.E.D. Notice that the difference operator on PAR performs the actual semantics of the difference operator on CSG. Theorem 4. For a CSG tree with movement as its root, the object represented by this CSG tree is identical to that represented by a PAR which is the movement of a PAR, which corresponds to the subtree of the original CSG tree. Proof: By the definition of movement of a PAR, the object represented by the PAR after applying a movement operation is the same as the object represented by the opiginal PAR (i.e. the structure of the object is not changed) except that the location and orientation of the object are changed. In the CSG representation scheme, an object represented by a CSG 14 A New Representation Scheme for Rotational Parts

RSD- TR-14-86 subtree has the same structure as the object represented by the CSG tree after applying a movement operation to the original CSG subtree except these two objects have different location and orientation in 3-D space. Therefore, after the same transformation matrix (i.e. movement operation) is applied to both the PAR and the original CSG subtree which represent the same object, the resulting objects should be identical; that is, they not only have the same structure but also have the same location and orientation in the 3-D space. Q.E.D. Notice that the movement operator on PAR performs the actual semantics of the movement operator on CSG. Theorem 5. An axis-symmetrical machine part represented by a CSG tree can be evaluated on PAR domain and the final resulting PAR after evaluation represents the same object as the CSG tree represents. Proof: A CSG tree is composed of operators (union, difference, movement) as nonterminal nodes and primitive solids (cylinder, cone, torus) as terminal nodes in this study. It is shown in Section 3.2 that these primitive solids can be represented by their corresponding PAR's. If the CSG tree is converted in bottom-up fashion from leaves to root, by Theorem 2, 3 and 4, the final resulting PAR should represent the identical object as the CSG tree represents. Q.E.D. Notice that by this theorem, a given CSG tree can be evaluated in PAR domain to obtain an identical object as the CSG tree represents. This theorem is the theoretical basis of our algorithms in Section 4. Theorem 6. Uniqueness Any axis-symmetrical machine part has a unique PAR representation. Proof: Our theorem may be rephrased as follows: For two representations, PAR(O1,A1) and PAR(02,A2), if 01 = 02 then PAR(O,A1) = PAR(02,A2). A New Representation Scheme for Rotational Parts 15

RSD-TR-14-86 Assume that the object O1 ( or 02 ) is in some 3D coordinate system. Then there must be a unique Principal Axis for the object in this coordinate system, thus Al = A2. Our proof procedure consists of two steps: Step 1: To prove PAR(0,,A1) and PAR(02,A2) have the same Axis Segment bounds with respect to Al and A2, respectively. That is, if PAR(01,A,) = { (S',Ci') I i =1,..,n; S, is Founded by [a b] }, PAR(02,A2) = { (S2,C2) i =,..,m; S is bounded by [a b } then we want to prove n = m, and a=a2 and b-b, for all i= 1,..,n. By the definition of PAR, all the curves within a bound must be differentiable and there exists at least one curve non-differentiable at one of its lower and upper bounds. This means that the bounds on the Principal Axis are uniquely defined by the discontinuity of the first derivatives of the curves specifying the object. Because the object has fixed shape and the bounded curve to specify some part of the shape of an object can be either an arc or a line segment but not both, the bounded curve to describe the object shape and thus the discontinuity points of the bounded curve can be uniquely specified in PAR. Therefore, there is a unique way to partition the Principal Axis in PAR and this unique partition decides a unique set of bounds. Step 2: Within each bound of (Si,C,') and (S,2,C2) of PAR(01,A1 ) and PAR(02,A2), respectively, we want to prove C, = C,2 if S/l = S,. This can be proved by contradition. Assume Ci' # C,2. Since Cl and Ci2, both are composed of a set of layers, if C1 # C,2, there must exist, for some j, a layer [c!cj41 J in C,' and its corresponding layer c)jci+] in C2 such that [el [cc., and generate 11 S) S J),+,],l and generate different volumes for 01 and 02. In this case 01 $ 02, which contradicts to our initial assumption 01 = 02. Therefore, C' = C, if S'-=S2 From the proofs of Stepl and Step2, in PAR, there is a unique set of bounds on the Principal Axis for an object, and on each pair of bounds there is a unique set of layers ( i.e. pairs of curves) to represent the object. Therefore, we conclude that there is only one unique PAR for an axis-symmetrical object and further PAR is a unique representation scheme. Q.E.D. 16 A New Representation Scheme for Rotational Parts

RSD-TR-14-86 4 Algorithm to Convert CSG Representation to PAR In this section, we will describe how to convert a CSG representation to its corresponding Principal Axis Repesentation (PAR). Section 4.1 first describes the data structure for PAR. The algorithm is described in Section 4.2, while Section 4.3 describes how to combine (either union or difference) two PAR's into one. The combination procedure is definitely a key component of the conversion algorithm. 4.1 Data Structure for PAR Since a PAR is composed of a set of tuples (Si, C,), we may use a record to represent each tuple and a link list to link these records together. Being a pair of bounds on the principal axis, each Si can thus be represented by a pair of real numbers and stored in the tuple record. Each Ci consists of a set of curves Cii, and can therefore be represented by a pointer to a link list where each element stores the parameters of one curve pair Ci, and Cii+l, that is, a record representing a layer of the machine part. We will assume the link list for tuples is sorted in the order of their bounds. Within each record for a tuple, the link list for the pairs of curves Cii is also kept sorted in their ordering (from Theorem 1). This makes the concept of layers easy to deal with because each pair of the curves enumerating from the beginning of the sorted sequence is just a layer. 4.2 Conversion from CSG Representation to PAR The conversion algorithm basically traverses the CSG tree from bottom to top. It first converts the CSG leaf nodes to their corresponding PAR's, and then combines (either union or difference) the PAR subtrees into composite PAR subtrees in higher levels. The combination procedure is repeated until the root is visited. The following recurive procedure describes this algorithm. function EvaluateCSG C ( T: CSCtree ): principalaxis_rep; var P, P1. P2: principalaxisrep; A New Representation Scheme for Rotational Parts 17

RSD- TR-14-86 begin case T.type of movement: P: EvaluateCSG ( T.leftchild^ ); Transform ( T.novement, P ); return ( P); union, difference: P1: EvaluateCSC ( T.leftchild^ ); P2:i EvaluateCSGC T.rightchild^ ); P: Combine ( T.type, P1, P2 ); return ( P ); primitive: P:= BuildAxisRep ( T ); return ( P); end; { case } end; { EvaluateCSG } To convert a CSG tree, its root can be passed to this function, which will recursively walk the whole CSG tree and transform it into a PAR. We assume that each node of the CSG tree is either a primitive solid (cylinder, cube, or torus) or an operator (movement, union or difference). Primitives (i.e. the CSG-tree leaf nodes) are converted to corresponding PAR's by BuildAzieRep procedure, which is based on the framework of Section 3.2. Using the data structure representation in Section4.1, the BuildAziaRep procedure simply creates one record to represent the primitive. In this record, the pair of bounds (in terms of the line parameter values of its principal axis) of the primitive solid are stored. A pair of curves representing a layer is also stored within the record. For a cylinder or cone, the outer curve of the pair is the line specifying the outer shape of the primitive solid and the inner curve is just the principal axis. For a torus, this pair is just the outer and inner half circles of the torus. 18 A New Representation Scheme for Rotational Parts

RSD TR-14-86 For the movement node in the CSG tree, its subtree (i.e. left-subtree) is evaluated first, then the transformation matrix stored within the movement node is applied to the evaluated subtree. The transformation matrix is a 4 by 4 matrix denoting the translation and rotation components of the movement. Since all the curves in PAR are specified relative to the principal axis, procedure Transform needs only transform the principal axis, not the whole set of curves. In fact, the transformation matrix is applied only to the two ending points of the principal axis. To deal with the union or difference node, its two subtrees are evaluated first, then the Combine procedure is called to union or difference them together. To implement the definitions of union and difference of two PAR's developed in Section 3.3, the Combine procedure includes four passes. We will describe these four passes in the next section. 4.3 Union and Difference of Two PAR's This algorithm (the procedure Combine) basically employs the split-and-merge paradigm, splitting the two composite axes into segments, computing the union or difference of two segments and then merging all the resulting segments into a new PAR. It has four steps (passes). Step 1. Find New Pairs of Bounds In this stage, the two operand subtrees being combined are represented in the PAR form, each being a link list of segments (records) in the order of the values of the bounded pairs. The function of this step is to compute the new bounded pairs for the resulting PAR. Since the values of the bounds in each operand PAR are in terms of line parameters of its own principal axis, we must first convert these values so that they are in terms of the line parameters of the resulting principal axis. Depending on the operation to be performed, the resulting principal axis can be chosen either as that of the first operand for difference operation or as the union of those of the two operands for the union operation. This selection process can be easily accomplished by simply checking the two ending points of the two operand principal axes. We assume that the two operands have the principal A New Representation Scheme for Rotational Parts 19

RSD. TR-14-86 axes in the same or opposite direction if they are to be combined. (Note that the opposite direction may occur because cones are not symmetrical in its top and buttom and thus have two directions.) After the resulting principal axis is roughly chosen and the values of the bounds are adjusted in term of this resulting axis, we may determine the pairs of all bounds for the resulting PAR. This process is simply the merge-and-sort procedure by repeatedly inputing two values of the bounds from each operand PAR and choosing the smaller value for the new bound. The state transition diagrams of this procedure are shown in Figure 1 and Figure 2. In Figure 1 and Figure 2, we may imagine that there are two stacks for two operand PAR's. Each stack stores the bounds of one operand PAR which are in ascending order with the smallest one on the top of the stack. A pair of bounds is just an even-odd pair of values on the stacks if the top element of the stack is numbered from zero. This implies that a common bound of two neighboring segments has duplicate values on the stack. Let's also assume that a is the top element of the first operand stack A and b is the top element of the second stack B, respectively. Now the following four states in the state transition diagrams can be defined. * State "-O" means a upper bound of a pair in stack A and a upper bound of a pair in stack B have just been popped out. * State "+0" means a lower bound of a pair in stack A and a upper bound of a pair in stack B have just been popped out. * State "+1" means a lower bound of a pair in stack A and a lower bound of a pair in stack B have just been popped out. * State'-1" means a upper bound of a pair in stack A and a lower bound of a pair in stack B have just been popped out. The starting and final state is State "-0". The notation { condition/ action; ction2;.... } is used in the state transition diagrams to indicate that, if the condition is satisfied, the 20 A New Representation Scheme for Rotational Parts

RSD.-TR-14-86 {a<=b / I=pop(A)} {a<=b / u=pop(A); output[I u]} {a>b / a>b / u=pop(B); {a>b / u-pop(B); a>b / unpop(B); I-pop(B)} output[l u]; I=u} output[l u]; IMu} output[I u]; I-u} (a<=b / u=pop(A); output[l u]; Iu}) {a<=b / u=pop(A); output[I u]; l-u) Figure 1: State Transition to Compute the New Bounds for Union. actions specified in the bracket are executed and the new state is pointed to by the arrow sign. Let's also use I and u to denote lower bound and upper bound of a bound pair. Note that these state transition diagrams may generate null bound pairs like [c c], where c is a real number, and they should be discarded. Step 2. Refine the Pairs of Bounds After Step 1, a set of pairs of bounds is obtained for the resulting PAR. For those pairs of bounds which are not subinrtervals common to both operands, the resulting shape within them can be immediately determined by the curves from the original (i.e. operand PAR) bound pairs. For those new bound pairs which are subintervals of both original operand PARs, the story is, however, not so simple. They should be refined. That is, subdivision of them must be considered because the resulting shape within a pair of bounds might be determined by curves from both original PARs. To refine the pairs of bounds which are the subintervals of both original PARs, we compute the intersection points of the curves within them. After sorting these intersecting points, a set of refined pairs of bounds may be created. The resulting curves within these refined pairs of bounds are then determined by the subcurves of the original curves. A New Representation Scheme for Rotational Parts 21

RSD- TR-14-86 {a<=b / I=pop(A)} {a<=b / u=pop(A); output[I u]} {a>b / a>b / I=pop(B)} {a>b / u=pop(B); a>b / u=pop(B); I=pop(B)} output[l u]; I=u} output[l u]; I=u) {a<=b / I=pop(A)) {a<=b / u=pop(A); output[I u]; I=u} Figure 2: State Transition to Compute the New Bounds for Difference. aO al a3 a4 a2 a5 Figure 3: Refine a Pair of Bounds by Computing Intersection Points. 22 A New Representation Scheme for Rotational Parts

RSD- TR-14-86 Consider the example of Figure 3. A cone and a torus are combined ( either union or difference). The cone is specified by lines A and B, and the torus is specified by the upper half circle C and lower half circle D within the bounds [al a2l. Since [al a2] obtained after Step 1 is both a subinterval of the [ao a2l of the cone and of the [al as] of the torus, we compute the intersection points for the curves A, B, C, D. Two intersection points a2 and as are obtained, then we create three refined bound pairs [al as], [as a4], and [a4 a2] to replace the [a, a2]. The curve shape within these three refined intervals can then be uniquely determined, which depends on the operation to be performed. We will discuss it in the next step of the algorithm. Step 3. Apply the Union or Difference Operation In this stage, the resulting curve shape within each pair of bounds can be determined. For those pairs of bounds created from Step 1 but not refined by Step2, that is, they are subintervals that belong to only one original pair of bounds, the curves within the original pair of bounds are directly copied into the newly created bound pair with the two ending points of the curves adjusted so as to be consistent with the new bounds. This works for the union operation. For the difference operation, the bound pairs contributed solely by the second operand are simply thrown away, but those bound pairs from the first operand should be kept and the curves defined within these bound pairs should also be copied. For those refined bound pairs from Step 2, it is assured that there is no curve intersection within them, thus from Theorem 1, the curves within them form a total ordering. This property of total ordering together with the concepts of maximum union-able and minimum difference-able layer aeta is useful for determining which curves from the original ones contribute to the resulting curves. Here we use the concept of layer to determine the resulting curves. A layer is an area bounded by two curves within some interval. Because the curves of layers of the two operands within the bound pair form a total ordering, they are topologically eqivalent to lines which are parallel to the principal axis within that bound pair. At any point within the bound pair, a line, passing this point and perpendicular to the principal axis, intersects all A New Representation Scheme for Rotational Parts 23

RSD- TR-14-86 {a<=b / I=pop(A)} {a<=b / u=pop(A); output[I u]) {a>b / a>b / u=pop(B); {a>b / pop(B)} a>b / pop(B)) I=pop(B)} output[l u]}) -1 {a<=b / pop(A)} +1 {a<=b / pop(A)) {a<=b / pop(A)} Figure 4: State Transition Diagram for Union. the curves within that bound pair. Layers can then be represented by segments (determined by these intersection points) on this line. If the intersection point of this line and the principal axis is assumed to be the zero point, then the coordinates of the original curves and those of segments for layers can thus be uniquely determined in this line coordinate. By using this line coordinate system, computing the union or difference of the line segments on this line is a simple task. It is similar to the procedure of merge-and-sort in Step 1. The coordinates of the line segments (for layers), which are kept sorted, from two operands are compared one by one, and depending on the operation performed (either union or difference), the resulting curves and layers can be determined. The following state transition diagrams (Figure 4 and Figure 5) demonstrate this algorithm. The meaning of the states and symbols is very similar to that in Figure 1 and Figure 2 except that the concept of bound pairs is replaced by the concept of layers and the values of the bounds are replaced by the coordinates of the curves in the line coordinate system. Consider again the example of Figure 3. We assume the cone is operated with the torus, the cone is the first operand (i.e. left subtree) and the torus is the second operand (i.e. right subtree). Within the bounds [al as], lines A and B form a layer for the cone and curves C and D form another layer for the torus. Since the layer bounded by A and B 24 A New Representation Scheme for Rotational Parts

RSD- TR- 14-86 {a<=b / I=pop(A)} {a<=b / u=pop(A); output(I u]) {a>b /pop(B) {a>b / pop(B)} {a>b / u=pop(B); {a>b / I=pop(B)} output(I u]} \/I {a<=b / pop(A)} {a<=b / pop(A)} Figure 5: State Transition Diagram for Difference. covers that bounded by C and D (this can be checked by their line coordinates), if the union operation is performed, the resulting layer would be that bounded by A and B. However, if the difference operation is executed, the results would be the layer bounded by A and C, and the layer bounded by D and B. Similarly, within the bounds [as a4], the result of the union operation is the layer bounded by the curves C and B, and the result of the difference operation is the layer bounded by D and B. Within the bounds [a4 a2J, the result of the union operation is the layer bounded by C and D and the layer bounded by A and B. The result of the difference operation would be the layer bounded by A and B. Step 4. Merge Neighboring Segments If Possible After the processing of the previous steps, the resulting shape of the combined object is obtained. However, it might not be in the PAR form because the curves at one of the two bound points might be differentiable. This may happen if some curves are seperated due to the bound pair refinement and later only those seperated subcurves are saved to be the resulting curves by the selection process of Step 3. In this case, the bound point at which all curves are differentiable is not a really breakpoint and it does not satisfy the A New Representation Scheme for Rotational Parts 25

RSD. TR-14-86 non-differentiable requirement of the PAR definition at the bound points. Therefore further processing is necessary to make the result a PAR. The actual work in this step is to go through the link list and check if the current segment is possible to be merged with its neighboring segments. The conditions for merging two neighboring segments are that 1. two neighboring segments have one common bound point 2. all the curves in one segment are one-to-one connected to all curves in another segment at the common bound point. 3. the two curves from both segments which are connected at the common point must have the common curve type and must belong to (i.e. subcurves) a curve whose curve type is identical to the common curve type. 4. two curves which satisfy condition (3) must be differentiable at the common bound point. After testing the merge condition, if two neighboring segments are mergeable, a segment is deleted from the list and the bound point of the extant segment is extended to cover the range of the deleted segment. Since all the curves in the original two segments are one-toone correspondance and the corresponding curves have the same curve parameters, nothing the bounds of the curves need to be updated. This step makes sure that the result of operating two PAR's is still a PAR. It is very important because it assures that, from Theorem 6, PAR is a unique representational scheme. The uniqueness property of a representational scheme is very useful because it minimizes the redundancy of storing data and eliminates the concern of data inconsistency in a database. Furthermore, it saves programmers' effort of writing code to analyze each of the possible representations for the same object and makes the access of objects in the database more efficient. 26 A New Representation Scheme for Rotational Parts

RSD- TR-14-86 I n I.. Figure 6.1 Figure 6.2 Figure 6.3 Figure 6.4 Figure 6: Examples. 5 Illustrative Examples In this section, the testing results of several examples by our program are shown in Figure 6 and Figure 7. This program implements the algorithms described in Section 4. Also implemented are the algorithms to compute the length and maximum diameter of a symmetrical part and to compute and plot its profile from PAR. The whole program code, which has more than 2,500 lines Pascal code on Apollo computer, can be found in Appendix 1, while Appendix 2 shows the input data for the example "bottle' of Figure 7. Appendix 3 shows all the pictures of the machine parts in the example of Figure 6 and 7. These pictures are generated from a ray tracing program with CSG trees as input data. Figure 6.1 shows the profile of a machine part which is joined together from three cylinders of different sizes. Some part of the middle cylinder originally overlapping with the other two is removed after processing, therefore it is not seen from the profile. Figure 6.2 shows the profile of a machine part which is unioned from two cones and a cylinder. Notice that the two cones are in opposite directions, that is, one of them is rotated 1800 with respect to y-axis from its originial orientation. Figure 6.3 is a profile of a neck. The neck is made from a cylinder, cut off a pipe and further a torus from its outer surface. In Figure 6.4, a pipe and two holes are drilled from a cone. The computation for this object A New Representation Scheme for Rotational Parts 27

RSD- TR-14-86 Figure 7: Example of Bottle. involves many layer difference operations. Its profile clearly displays the layer information. Figure 7 displays the profile of a "bottle". Its length and diameter are also computed by the program (Length = 5.30 units, Diameter = 4.905 units). This "bottle" is composed of 27 primitive solids involing 8 cones, 12 cylinders and 7 tori (see Appendix 2, prmitive file). Its CSG tree contains 53 nodes (object file). Since there is no movement node used in this example, its movement file is empty. Many principal azis union and difference operations as well as layer union and layer difference operations are involved in converting the CSG tree into its PAR representation. Its length, diameter information and profile are then computed from this PAR. To compute a profile of a machine part, its algorithm may take advantages of the PAR representation. Since the layer information which consists of its outer, inner curves, and starting and ending points exists in the PAR representation, the profile of a machine part can be generated by drawing the four boundary curves (two vertical lines for the starting and endinding line boundaries, and two for the inner and outer curves) for each layer in an ezclssive or plotting mode. In such mode, the common boundary of two layers will dissappear, rather than being plotted twice. This algorithm is efficient and easy to implement, it has some problems, however. It 28 A New Representation Scheme for Rotational Parts

RSD- TR-14-86 relies on the exclusive or operation to eliminate the common boundary of two layers, but the exclusive or operation also eliminates the conner points of actual boundary because the intersection point of two boundary curves is plotted twice. In general, a right-turn algorithm [WeA77] can be used to solve this problem. The vertices of layers and their relations in the principal axis coordinate system can be easily obtained from the PAR. Then the rightturn algorithm can be applied to these vertices to get the profile. At this time, we only implement the profile computation using the exclusive or algorithm. Close examination of Figure 6 and Figure 7 reveals the fact of the missing common corner points. 6 Conclusion and Discussion In this section, some concluding remarks are made first in Section 6.1, then several issues about PAR and its future extension are discussed in Section 6.2. 6.1 Conclusion As a new representational scheme toward representational uniqueness, the Principal Axis Representation (PAR) is developed for axis symmetrical machine parts. Based on the CSG, the machine parts can be constructed by applying the regularized set operations union, difference and movement on the primitive solids cylinder, cone and torus. The PAR can represent the machine parts composed of primitives that have the same Principal Axis. The PAR is first defined in this report. The operations for operating ( combining) PAR's: union, difference and movement are then defined. Also shown is how to use PAR to represent primitive solids: cylinder, cone and torus. Based on this formulation, several theorems are proved, which show that a CSG tree can be converted into a PAR representing the same axis-symmetrical object. An algorithm based on the mathematical formulation of PAR to convert a CSG tree has been designed and implemented in this study. It first converts the CSG leaf nodes into their PAR representation and then operates these PAR's (using union, difference, movement operations of PAR ) on the CSG tree from bottom to top until the root of the CSG tree is A New Representation Scheme for Rotational Parts 29

RSD- TR-14-86 visited. Several examples and testing results of this algorithm are shown in Section 5. To support CAD applications, for example, generating the profile of machine parts for Numerical Control and computing the geometrical properties of a machine part such as its length and maximum diameter for part classification and CAPP (Computer Aided Process Planning), the PAR seems to be more efficient than its counterpart: the CSG scheme. This is due to the fact that, to represent an object, the CSG tree is left unevaluated while the PAR is already the result after evaluation. Algorithms that compute the profile and the length as well as diameter of a machine part from PAR have also been implemented. Its results are shown in Section 5. In addition to its computational efficiency, the PAR is proved to possess yet another important property as a representational scheme, that is, representational uniqueness. A unique representatation scheme allows much simpler feature definition and therefore feature eztraction or object recognition because only one representation for a feature or an object is required to deal with[LeK86]. 6.2 Discussion and Euture Work The key ideas of PAR are to represent an object by its principal axis and boundary curves, and to resolve the overlap (intersection) of two composite solids, i.e. to compute the union and difference of layers. At present, PAR works for axis-symmetrical machine parts that are constructed from primitives such as cylinder, cone and torus. To make it more general as a unique representational scheme, several extensions are currently under investigation: * relax the assumption of common principal axis requirement. Based on the PAR, a hierarchical PAR could be defined as a set of PAR's, each having its own common principal axis. But the problem of how to characterize the overlap parts of two PAR's during union and difference operations to ensure that the generalized PAR is also a uniqueness representation needs to be studied. * include more primitives like sphere and cube. In fact, sphere is a special case of torus. But naively integrating sphere into PAR would invalidate the uniqueness property 30 A New Representation Scheme for Rotational Parts

RSD- TR14-86 of the PAR because a sphere has two different representations in PAR, i.e. using sphere primitive or torus primitive. Including cube into PAR relaxes the constraint of axis-symmetry to a large extent. How to characterize the boundary curves, however, requires further study. To make the PAR more useful, several areas of application are also under study: * support feature extraction and object recognition. Since the PAR is representational uniqueness, object features such as round, fillet, keytay, hole can be defined on the PAR more easily than on the CSG tree. Developing algorithms to extract features or recognize object based on the PAR should be simpler because less cases need to be analyzed. * support CAD (Computer Aided Design) applications. In this report, we have shown the profile generation of machine parts for NC (Numerical Control) and the length and diameter computation for machine parts to support part classification and CAPP (Computer Aided Process Planning). Other geometrical properties and geometrical codes [KaO84] could also be computed easily from the PAR. We are working on exploit it now. A New Representation Scheme for Rotational Parts 31

RSD- TR-14-86 References [ArCP841 G.T. Armstrong, G.C. Carey, and A. de Pennington, "Numerical Code Generation from A Geometric Modeling System," in Solid Modeling by Computers: From Theory to Applications, ed. by M.S. Pickett and J.W. Boyse, Plenum Press, 1984. [Bee82] W. Beeby,'The Future of Integrated CAD/CAM Systems: the Boeing Perspective," Computer Graphics and Applications, Vol.2, No.l, Jan.1982, pp51-56. [BoG82] J.W. Boyse and J.E. Gilchrist,'GMSolid: Interactive Modeling for Design and Analysis of Solids," Computer Graphics and Applications, Vol.2, No.2, March 1982, pp27-40. [Gra761 A.R. Grayer, "The Automatic Production of Machined Components Starting from A Stored Geometric Description," PROLAMAT Proceedings, North Holland Publishing Co., 1976. [Hen84] M.R. Henderson, "Extraction of Feature Information from Three Dimensional CAD Data,' Ph.D. Thesis, School of Mechanical Engineering, Purdue University, May, 1984. [Jak82] R. Jakubowski, "Syntactic Characterization of Machine Parts Shapes," Cybernetica and System.: An International Journal, 13, 1982, ppl-24. [KaO84] Y. Kakazu and N. Okino,'Pattern Recognition Approach to GT Code Generation on CSG," Proc. 16th CIRP International Seminar on Manufacturing Systems, Tokyo, 1984, pp. 10-18. [Kyp80] L.K. Kyprianou, "Shape Features in Computer-Aided Design," Ph.D. Thesis, University of Cambridge, Cambridge, England, July, 1980. [LeK86] Y.C. Lee and K.S. Fu, "Machine Understanding of CSG: Extraction and Unification of Manufacturing Features," Computer Graphics and Applications, (to appear as a regular paper in December, 1986) 32 A New Representation Scheme for Rotational Parts

RSD- TR-14-86 [ReV82] A. A. G. Requicha and H. B. Voelcker,'Solid Modeling: A Historical Summary and Contemporary Assessment," Computer Graphics and Applications, Vol.2, No.2, March 1982, pp. 9-24. [StHA83 S.M. Staley, M.R. Henderson, and D.C. Anderson, "Using Syntactic Pattern Recognition to Extract Feature Information from a Solid Geometric Data Base," Computers in Mechanical Engineering, Vol.2, No.2, Sept. 1983, pp61-66. [WeA77] K. Weiler and P. Atherton, "Hidden Surface Removal Using Polygon Area Sorting," Computer Graphics, Vol.11, No.2, March 1977, pp. 9-24. [Woo77] T.C. Woo, "Computer Aided Recognition of Volumetric Designs," in Advances in Computer-Aided Manufacture, edited by D. Mcpherson, North-Holland Publishing Company, 1977, pp121-136. [Woo82] T.C. Woo,'Feature Extraction by Volume Decomposition," Proc. Conference on CAD/CAM Technology in Mechanical Engineering, MIT, Mar. 1982, pp76-94. A New Representation Scheme for Rotational Parts 33

Append/x i. Program Listing

Jul 8 02:16 1986 Page 1 Jul 8 02:16 1986 Page 2 1 57 2 (*58 fupper or lower part of a circle 3 ( This program implements how to convert a CSG * 59 upor~_down =( up, down ) 4 ( input tree into a PAR (Principal Axis Representa- * 60 5 (~tion) that represents the identical object as the *)61 fcurve type 6 ( original CSG tree does. ~ 62 curve-kind ( line, arc ) 7 (* All the algorithms and procedures implemented A)63 8 ( here are based on the definitions of PAR, its A 64[ possible states when two axes are operated 9 ( union, difference and movement operations. A 65 state=( pluszero, minuszero, plusone, iuoe) 10 (A Geometrical properties such as the length and A)66 11 (Adiameter of objects are computed from the result- *)67 (point coordinate on the principal axis coorint 12 (Aing PAR to support part classification and process*) 68 D2_point =record 13 ( planning. Profile is also generated to support A 69 x coord: t value; 14 (~Numerical Control application. A)70 y~coord:y_value; 15 ( )71 end; 16 ( Author K. F. Jack Jea A 72 17 (A Date:6/12/1986 A)73 (actual coordinate in the 3-D space 18 (A)74 D3_point =record 19 (AAAAAAAAAAAAAA*AAAAAAAAAAAA)75 x-coord, 20 76 y-coord, 21 PROGRAM AXIS( input, output ) 77 z coord: real; 22 78 end; 23 %nolist; 79 24 %include'/sys/ins/base. ins. pas'i 80 (chain to link intersection points of two ae 25 % include'/sys/ins/error. ins. pas'; 81 intersect rec~ptr = intersect-rec; 26 %include'/sys/ins/gpr.ins.pas'; 82 27 %include'/sys/ins/time. ins. pas'; 83 (record for storing the intersection points 28 %list; 84 intersect-rec =record 29 85 intersection: t value 30 CONST 86 next-intersect: intersec epr 31 fpart number of the csg tree root 187 end; 32 CSGROOT = -1; 88 33 Pi = 3.1415926; 89( record to store a bounded curve 34 90 curve =record 35 TYPE 91 curve-start, 36 92 curve end:y-value; 37 (normalized parametric value on the axis 193 case curve -type curve-kind of 38 t value =real; 94 line:() 39 95 arc: 40 (value of vertical axis on principal axis coordinate 196 center:Dpit 41 y~value = real; 97 radius:rel 42 98 which_half:upodwn) 43 (disparity of directions of two axes 199 end; 44 direction = (same dir, opposite dir, non co linear); 100 45 101 fpointer to a bounded curve 46 (basic primitive solids currently implemented 1102 curve~ptr - curve; 47 solid-kind =(cylinder, cone, torus );103 48 104 1pointer to a layer 49 1node types of the CSG tree 1105 layer~ptr - layer; 50 node-kind =(primitive, movement, union -op, diff-op );106 51 107 1a layer which is bounded by two curves 52 1axis operation type of two principal axes 1108 layer =record 53 axi~s-_op kind =(union, difference );109 next-layer layer~ptri 54 110 inner-curve, 55 (lower or upper bounds il11 outer-curve:curve-ptr; 56 lower-or upper = ( lower, upper );112 end;

Jul 8 02:16 1986 Page 3 Jul 8 02:16 1986 Page 4 113 169 1 global variables for main program use ] 114 1 pointer to a segment of the axis ] 170 VAR 115 segmentptr = segment; 171 1 original input CSG tree ] 116 172 CSG: CSG_treeptr; 117 f segment record, part of axis ] 173 [ resulting PAR corresponding to CSG above 118 segment = record 174 P: principal_axis; 119 lower bound, 175 1 length and maximum diameter of the object 120 upper_bound: tvalue; 176 L, D: real; 121 next_segment: segmentptr; 177 1 input data for creating the CSG 122 layer_head: layerptr; 178 obj_file, mov file, prim_file: text; 123 end; 179 1 starting and ending points of a combined axis 124 180 min_point, max_point: D3_point; 125 [ information of the principal axis ] 181 126 principal_axis = record 182 1 switches for debugging use ] 127 start_point, 183 DEBUG_INPUT: boolean; 128 end_point: D3_point; 184 DEBUG_EVALUATE_CSG: boolean; 129 segment_head: segment_ptr; 185 DEBUGBUILD_ AXIS: boolean; 130 end; 186 DEBUGLD: boolean; 131 187 DEBUG_DRAW: boolean; 132 [ transformation matrix for movement ] 188 DEBUG_COMPUTE_T: boolean; 133 xformmatrix = record 189 DEBUG_CURVE_CURVE, DEBUG_INTERSECT: boolean; 134 translatex, 190 DEBUG_AXIS_OP, DEBUG_NORMALIZATION, 135 translate_y, 191 DEBUG_PARTITION, DEBUG_MERGE_INTERVAL: boolean; 136 translatez: real; 192 DEBUG_LAYER_DIFFERENCE, DEBUG_COMPUTE_LAYER, 137 rotate x, 193 DEBUG_ADD_LAYER: boolean; 138 rotate y, 194 TRACE_AXIS: boolean; 139 rotate z: real; 195 140 end; 196 procedure dump axis( PX:principal_axis); 141 197 1 dump the principal axis; for debugging purpose ] 142 [ information of a primitive solid ) 198 var 143 solid = record 199 s: segmentptr; 144 base, top: D3point; 200 1: layer_ptr; 145 case solid_type: solid_kind of 201 begin 146 cylinder: ( radius: real ); 202 writeln('@@@(@ Dump Principal Axis @@@ @'); 147 cone: ( radius _b, 203 writeln('start point =', PX.startpoint.xcoord, 148 radiust: real ); 204 PX.startpoint.y_coord, PX.start point.z_coord) 149 torus: ( inner _radius, 205 writeln('end point =', PX.endpoint.xcoord, 150 outer radius:real ); 206 PX.endpoint.y_coord, PX.endpoint.zcoord); 151 end; { solid ] 207 152 208 [ go thru the axis and print its content ] 153 1 pointer to a CSG tree ] 209 s:= PX.segment_head; 154 CSG_tree ptr = CSG_tree; 210 while s <> nil do 155 211 begin 156 { CSG tree node information 212 157 CSG _tree = record 213 writeln('lower and upper bounds=', 158 node_type nodekind; 214 s.lower_bound, s.upper_bound); 159 partno: integer; 215 160 case nodekind of 216 1:= s^.layer_head; 161 primitive: ( prim_solid: solid ); 217 while 1 <> nil do 162 movement: ( move xform_matrix; 218 begin 163 child: CSG_treeptr ); 219 164 union_op, 220 writeln('outer curve, start and end y-value=', 165 diff op: ( leftchild, 221 1 outercurve curve_start, 166 right _child: CSGtreeptr); 222 1 outer_curve curve_end); 167 end; 223 writeln('inner curve, start and end y-value=', 168 224 1 innercurve curvestart,

Jul 8 02:16 1986 Page 5 Jul 8 02:16 1986 Page 6 225 1A.inner_curve-.curveend); 281 t: ml + ( t - fl ) * ( m2 - ml ) / ( f2 - fl ); 226 282 227 1:= 1l.next_layer; 283 end; { modularize_t 228 end; 284 229 285 procedure adjust t value( var P: principal_axis; 230 s:= s^.next_segment; 286 fl, f2, ml, m2: t_value ); 231 end; 287 { This routine use ml and m2 as the factors to adjust the 232 288 t parameters in the principal axis P; formula: 233 writeln('l@@@@ End of Dumping Principal Axis @~@'); 289 t' = ml + ( t - fl ) * ( m2 - ml ) / ( f2 - fl ) ] 234 290 var 235 end; { dumpaxis ] 291 s: segment_ptr; 236 292 1: layer_ptr; 237 293 begin 238 function computet ( ptl, pt2, point: D3_point ): tvalue; 294 239 f calculate the t parameter for point in line segment 295 f go thru each bound and t-value dependent item. ] 240 bounded by ptl and pt2. Note: we do not check the linearity 296 s:= P.segmenthead; 241 of point, ptl and pt2 in this version; it will be refined in 297 while ( s <> nil ) do 242 later version. ] 298 begin 243 const 299 244 { numerical error range 1 300 modularizet ( s.lower_bound, fl, f2, ml, m2 ); 245 EPS = 0.0001; 301 modularize_t ( s^.upperbound, fl, f2, ml, m2 ); 246 var 302 247 t: t_value; 303 1:= s.layer_head; 248 begin 304 while( 1 <> nil ) do 249 305 begin 250 if DEBUG _COMPUTE _T then writeln( 306 if 1.innercurve.curvetype = arc then 251'%%%%%%% In computet, ptl, pt2, point=', 307 modularize t ( 252 ptl.x_coord, ptl.y_coord, ptl.z coord, 308 1.innercurve.center.x_coord, 253 pt2.x_coord, pt2.y_coord, pt2.z coord, 309 fl, f2, ml, m2 ); 254 point.xcoord, point.y_coord, point.z coord ); 310 if l^.outercurve^.curvetype = arc then 255 311 modularizet ( 256 if abs(ptl.x-coord - pt2.x coord) >= EPS then 312 l1.outercurve-.center.x_coord, 257 t:= (point.x-coord- ptl.xcoord) / 313 fl, f2, ml, m2 ); 258 (pt2.x coord - ptl.xcoord) 314 1:= 1.nextlayer; 259 else if abs(ptl.ycoord - pt2.y_coord) >= EPS then 315 end; [ inner while ] 260 t:= (point.y_coord - ptl.y_coord) / 316 261 (pt2.ycoord - ptl.ycoord) 317 s:= s.next_segment; 262 else if abs(ptl.zcoord - pt2.z_coord) >= EPS then 318 263 t:= (point.zcoord - ptl.zcoord) / 319 end; [ while 264 (pt2.zxcoord - ptl.zcoord) 320 265 else writeln('ERROR: ptl and pt2 are the same point,', 321 end; ( adjust t value } 266'in computet procedure'); 322 267 323 function check_direction( P1, P2: principalaxis ):direction; 268 if DEBUG_COMPUTE_T then writeln('%%%%%%% In computet, t =' 324 [ check to see if P1 and P2 are in the same, opposite or 269, t); 325 non-co-linear direction 270 326 const 271 compute_t:= t; 327 EPS = 0.0001; 272 328 var 273 end; [ computet 1 329 xl, yl, zl, x2, y2, z2, dl, d2: real; 274 330 begin 275 procedure modularize_t ( var t: t_value; 331 276 fl, f2, ml, m2: tvalue ); 332 1 get three vector components of P1 and P2 1 277 1 actually modularize t here; 333 xl:= Pl.start_point.xcoord - Pl.end_point.xcoord; 278 t' = ml + ( t - fl ) * ( m2 - ml ) / ( f2 - fl ) ] 334 yl:= Pl.start_point.y_coord - Pl.endpoint.y_coord; 279 begin 335 zl:= Pl.start_point.z coord - Pl.endpoint.z_coord; 280 336 dl:= sqrt( xl * xl + yl * yl + zl * zl );

Jul 8 02:16 1986 Page 7 Jul 8 02:16 1986 Page 8 337 xl:= xl / dl; 393 1.inner_curve^.center.xcoord:= 1.0 - 338 yl: yl / dl; 394 l^.inner_curve.center.x_coord; 339 zl:= zl / dl; 395 340 396 y:= 1.outercurve curve_start; 341 x2:= P2.startpoint.x_coord - P2.end point.xcoord; 397 l^. outer_curve curve_start:= 342 y2:= P2.startpoint.y_coord - P2.end point.y_coord; 398 1.outer_curve curve_end; 343 z2:= P2.start_point.zcoord - P2.endpoint.zcoord; 399 l^.outercurve^.curveend:= y; 344 d2: sqrt( x2 * x2 + y2 * y2 + z2 * z2 ); 400 345 x2:= x2 / d2; 401 if 1.outer_curve^.curve_type = arc then 346 y2: y2 / d2; 402 l^.outer_curve.center.xcoord:= 1.0 - 347 z2: z2 / d2; 403 l1.outer_curve.center.x_coord; 348 404 349 [ check linearity } 405 1:= 1.next_layer; 350 if ( abs(xl - x2) <= EPS) and ( abs(yl - y2) <= EPS) and 406 end; { inner while 351 ( abs(zl - z2) <= EPS) then 407 352 checkdirection:= same_ dir 408 ns:= s.next segment; 353 else if ( abs(xl + x2) <= EPS) and ( abs(yl + y2) <= EPS) 409 s.next_segment:= ps; 354 and ( abs(zl + z2) (= EPS) then 410 ps:= s; 355 checkdirection:= opposite_dir 411 s: ns; 356 else checkdirection:= noncolinear; 412 357 413 end; t while } 358 end; 1 checkdirection } 414 359 415 P.segment_ head:= ps; 360 procedure reversedirection( var P: principal_axis ); 416 361 [ adjust the t-parameters in P; make all the t-value reverse ] 417 end; i reverse_direction } 362 var 418 363 s, ps, ns: segment_ptr; 419 procedure normalization ( var P1, P2: principal_axis; 364 1: layerptr; 420 var ptl, pt2: D3point ); 365 pt: D3_point; 421 { normalize both P1 and P2 to have the same t-parameter 366 t: t_value; 422 basis; note: after difference operation, the principal 367 y: y_value; 423 axis should re-normalized. 368 begin 424 var 369 425 tl, t2: tvalue; 370 pt:= P.startpoint; 426 ml, m2: real; 371 P.startpoint:= P.end_point; 427 dir: direction; 372 P.endpoint:= pt; 428 begin 373 429 374 s:= P.segment_head; 430 if DEBUG_NORMALIZATION then writeln( 375 ps: nil; 431'%%%%% before normalization,Pl.start and end point=, 376 432 Pl.startpointox_coord, Pl.startpoint.y_coord, 377 while ( s <> nil ) do 433 Pl.startpoint. zcoord, Plo.end_pointoxcoord, 378 begin 434 Pl.endpoint.y_coord,Pl.endpoint.zcoord) 379 t:= s lowerbound 435 380 s lowerbound:= 1.0 - s.upper_bound; 436 if DEBUG_NORMALIZATION then writeln( 381 s^ upper_bound:= 1.0 - t; 437'%%%%% before normalization,P2.start and end point=', 382 438 P2.startpoint.x_coord, P2.start point.y_coord, 383 1:= s layer_head; 439 P2.startpoint.z_coord, P2.endpoint.x_coord, 384 whiile ( 1 <> nil ) do 440 P2.endpoint.y_coord,P2.endpoint.zcoord) 385 begin 441 386 442 if DEBUGNORMALIZATION then writeln( 387 y:= 1^ innercurve^.curvestart; 443'%%%%% before normalization,P2.low, up =', 388 l^.innercurve curvestart:= 444 P2.segment_head lower_bound, 389 1.innercurve curveend; 445 P2.segment-head upperbound); 390 l^.innercurve.curve end:= y; 446 391 447 [ check the direction of two axes, ok if the same; 392 if l^.inner curve-.curve type = arc then 448 inverse one if in reverse direction, reject it

Jul 8 02:16 1986 Page 9 Jul 8 02:16 1986 Page 10 449 if not colinear ] 505 this routine is somewhat like getnext function, but it 450 dir:= check_direction( P1, P2. ); 506 will be used by computeintersection and computelayer. 451 case dir of 507 begin 452 samedir: f ok; go thru it } 508 453 if DEBUG_NORMALIZATION then writeln( 509 case flag of 454'%%%%% In normalization,', 510 upper: if p = nil then 455'two axes same direction' ); 511 getcurve:= nil 456 oppositedir: if (Pl.startpoint.zcoord > 512 else begin 457 Pl.end_point.zcoord) then 513 flag:= lower; 458 reversedirection( P1 ) 514 getcurve:= p.outer_curve; 459 else 515 end; 460 reverse_direction( P2 ); 516 lower: begin 461 non_co_linear: writeln(' --— > error,', 517 flag:= upper; 462' two axes do not co-linear'); 518 getcurve:= pA.inner_curve 463 end; [ case ] 519 p:= p.next_layer; 464 520 end; 465 { compute the new starting and ending points 1 521 end; [ case ] 466 tl:= computet( Pl.startpoint, Pl.endpoint, 522 467 P2.start_point ); 523 end; ( getcurve ] 468 t2:= computet( Pl.startpoint, Pl.end_point, 524 469 P2.end_point ); 525 procedure lineline( Cl, C2: curveptr; 470 526 Sl, S2: segment_ptr; 471 if DEBUG_NORMALIZATION then writeln( 527 a, b: tvalue; 472'%%%%% IN normalization, tl, t2 =', tl, t2 ); 528 var Iptr: intersect_recptr ); 473 529 compute the intersection of two lines Cl and C2; 474 if 0.0 <= tl then ptl:= Pl.start_point 530 results are appended to Iptr. 475 else ptl:= P2.startpoint; 531 label 476 if 1.0 >= t2 then pt2:= Pl.endpoint 532 R; 477 else pt2:= P2.end_point; 533 const 478 534 SLOPE EPSILON = 0.001; 479 1 adjust P1 and P2 using the new starting and 535 var 480 ending points 1 536 xl, x2, x3, x4, dx, tl, x: t_value; 481 ml: computet( ptl, pt2, Pl.startpoint ); 537 yl, y2, y3, y4, dy: real; 482 m2:= compute_t( ptl, pt2, Pl.end_point ); 538 p: intersectrecptr; 483 adjust t value( P1, 0.0, 1.0, ml, m2 ); 539 begin 484 540 485 ml:= computet( ptl, pt2, P2.startpoint ); 541 1 get the first line segment } 486 m2:= compute_t( ptl, pt2, P2.endpoint ); 542 xl:= Sl^.lower_bound; 487 adjust t value( P2, 0.0, 1.0, ml, m2 ); 543 x2:= Sl1.upper_bound; 488 544 yl:= Cl.curve_start; 489 if DEBUG_NORMALIZATION then writeln( 545 y2:= Cl.curveend; 490'%%%%% after normalization,P2.low, up =', 546 491 P2.segmenthead^.lower_bound, 547 492 P2.segmenthead^.upper_bound); 548 [ get the second line segment ] 493 549 x3:= S2^.lower_bound; 494 if DEBUG _NORMALIZATION then writeln( 550 x4:= S2.upperbound; 495'%%%%% result of normalization,Pl.start and end point=', 551 y3:= C2.curve_start; 496 Pl.startpoint.x_coord, Pl.startpoint.y_coord, 552 y4: C2.curve_end; 497 Pl.start point.zcoord, Pl.endpoint.xcoord, 553 498 Pl.endpoint.y_coord,Pl.end_point. zcoord); 554 if DEBUGCURVECURVE then 499 555 writeln('$$$$$$ Enter In lineline,', 500 end; [ normalization ] 556'xl,yl,x2,y2,x3,y3,x4,y4', 501 557 xl,yl,x2,y2,x3,y3,x4,y4 ); 502 function getcurve( var p: layer_ptr; 558 503 var flag: lower_orupper ): curveptr; 559 if ( ( xl = x3 ) and ( yl = y3 ) ) and 504 [ get next curve from the layers; if none, then return nil. 560 ( ( x2 = x4 ) and ( y2 = y4 ) ) then

Jul 8 02:16 1986 Page 11 Jul 8 02:16 1986 Page 12 561 f lines coincidence 617 562 begin 618 length:= sqrt( (maxpoint.x_coord - minpointxcoord) 563 { return 619 *(max_point.x_coord - minpoint.xcoord) 564 end 620 +(maxpoint.y_coord - minpoint.ycoord) 565 else if (xl = x2) and (x3 = x4) then 621 *(maxpoint.ycoord - minpoint.ycoord) 566 begin 622 +(maxpoint.z_coord- minpoint.zcoord) 567 {return] 623 *(maxpoint.z_coord - min_point.zcoord) ); 568 end 624 569 else begin 625 f get the first line segment 3 570 if (xl <> x2) and (x3 <> x4) then 626 xl:= Sl^.lower_bound * length; 571 begin 627 x2:= Sl upperbound * length; 572 if abs( (y2 - yl)/(x2 - xl) - 628 yl:= Cl^.curve_start; 573 (y4 -y3)/(x4 - x3) ) <= SLOPEEPSILON 629:= Cl.curve_end; 574 then begin 630 575 goto R; 631 576 end; 632 get the second arc segment ] 577 end; 633 x3 = S2^ lower_bound * length; 578 dy: y4 - y3; 634 x4:= S2 upperbound * length; 579 dx:= x4 - x3; 635 y3:= C2- curve_start; 580 636 y4:= C2 curve_end; 581 if DEBUG CURVECURVE then writeln( 637 xO: C2 center.x_coord * length; 582'$$$$$$ In line line, dx, dy, yl, y2 =',638 yO: C2 center.y_coord; 583 dx, dy, yl, y2 );639 584 640 if DEBUGCURVECURVE then 585 tl:= dy * (x3 - xl) - dx * (y3 - yl); 641 writeln('$$$$$$ In linearc,xl,yl,x2,y2' 586 tl:= tl / ( dy * (x2 - xl) - dx * (y2 - yl) );642 xl,yl,x2,y2 ) 587 x:= xl + tl * ( x2 - xl ); 643 if DEBUGCURVECURVE then 588 if ( a < x ) and ( x < b ) then 644 writeln('$$$$$$ In linearc,xO,yO,x3,y3,x4,y4', 589 begin 645 xO,yO,x3,y3,x4,y4 ); 590 new( p );646 591 p.intersection:= x; 647 if( (xl = x3 ) and ( yl = y3) )and 592 p.next_intersect:= Iptr; 648 ( x2 = x4 ) and ( y2 = y4 ) ) then 593 Iptr:= p; 649 curves intersect at end point) 594 650 begin 595 if DEBUGCURVECURVE then writeln( 651 1 return 596'$$$$$$ In lineline, intersection t = 652 end 597, x );653 else begin 1 computate intersection of line and arc 598 end; 654 dx:= x2 - xl; 599 end 655 dy:= y2 - yl; 600 R: 656 AA:= ( (yo - yl) * dx - (xO - xl) * dy ) / 601 657 C2.radius; 602 end; [ line line ) 658 d:= dx * dx + dy * dy - AA * AA; 603 659 604 procedure line arc( Cl, C2: curveptr; 660 if DEBUGCURVE_CURVE then writeln( 605 Sl, S2: segment_ptr; a, b t_value; 661'$$$$$$ In linearc, dx, dy, AA, d', 606 var Iptr: intersect_recptr ); 662 dx, dy, AA, d ); 607 1 compute the intersection of a line Cl and an arc C2; 663 608 results are appended to Iptr. 664 if d >= 0.0 then 609 } 665 begin 610 var 666 if abs( AA + dy ) <= 0.000001 then 611 xO, xl, x2, x3, x4, dx, x, length: real; 667 begin 612 yO, yl, y2, y3, y4, dy, tl, t2, AA, d: real; 668 theta[l): PI; 613 p: intersectrecptr; 669 theta[2]:- PI; 614 theta, ans: array[l..2] of real; 670 end 615 num, i: integer; 671 else begin 616 begin 672 tl:= (-dx + sqrt(d) ) / ( AA + dy );

Jul 8 02:16 1986 Page 13 Jul 8 02:16 1986 Page 14 673 t2:= (-dx - sqrt(d) ) / ( AA + dy ); 729 +(maxpoint.y_coord - minpoint.ycoord) 674 theta[l]:= 2.0 * arctan( tl ); 730 *(maxpoint.y_coord min_point.ycoord) 675 theta[2]:= 2.0 * arctan( t2 ); 731 +(max_point.zcoord - min_point.zcoord) 676 end; 732 *(max_point.zcoord -minpoint.z-coord) ); 677 733 678 num:= 0; 734 { get the first arc segment ] 679 case C2^.which_half of 735 xl:= Sl^.lower_bound * length; 680 down: { take negative angles ] 736 x2:= Sl^.upperbound * length; 681 for i:= 1 to 2 do 737 yl:= Cl.curve_start; 682 if theta[i] <= 0.0 then 738 y2 C:= cl.curve_end; 683 begin 739 xcl:= Cl^.center.x_coord * length; 684 num:= num + 1; 740 ycl:= Cl^.center.ycoord; 685 ans[num]:= theta[i]; 741 686 end; 742 ( get the second arc segment 687 up: { take positive angles ] 743 x3:= S2^.lower_bound * length; 688 for i:= 1 to 2 do 744 x4:= S2^.upperbound * length; 689 if theta[i] >= 0.0 then 745 y3:= C2^.curve_start; 690 begin 746 y4:= C2^.curve_end; 691 num:= num + 1; 747 xc2: C2^.center.x_coord * length; 692 ans[num]:= theta[i]; 748 yc2:= C2.center.y_coord; 693 end; 749 694 end; [ case ] 750 if ( ( xl = x3 ) and ( yl = y3 ) ) and 695 751 ( (x2 = x4 ) and ( y2 = y4 ) ) then 696 for i:= 1 to num do 752 [curves intersect at end point) 697 begin 753 begin 698 x: xO + cos( ans[i] ) * C2^.radius; 754 { return ] 699 if ( a < (x/length) ) and 755 end 700 ( x < (b/length) ) then 756 else begin { computate intersection of line and arc 701 begin 757 dx: xc2 - xcl; 702 new( p ); 758 dy:= yc2 - ycl; 703 p^.intersection:= x / length; 759 704 p.next_intersect:= Iptr; 760 if DEBUG_EVALUATE_CSG then writeln( 705 Iptr:= p; 761'$$$$$ In arc-arc, C2 radius =', C2.radius ); 706 end; 762 707 end; [ for ] 763 AA:= ( Cl^.radius * Cl^.radius - C2^.radius * 708 end; if 764 C2^.radius - dx * dx - dy * dy ) / 709 end; ( else ] 765 ( 2.0 * C2.radius ); 710 766 d: dx * dx + dy * dy - AA * AA; 711 end; { line arc ] 767 712 768 if d >= 0.0 then 713 procedure arcarc( Cl, C2: curveptr; 769 begin 714 Sl, S2: segment_ptr; a, b: t_value; 770 d:= sqrt( d ); 715 var Iptr: intersect_rec_ptr ); 771 if abs( AA + dy ) <= 0.000001 then 716 f compute the intersection of an arc C1 and an arc C2; 772 begin 717 results are appended to Iptr. 773 theta[l]: PI; 718 ] 774 theta[2]: -PI; 719 var 775 end 720 xcl, xc2, xl, x2, x3, x4, dx, x, length: real; 776 else begin 721 ycl, yc2, yl, y2, y3, y4, dy, tl, t2, AA, d: real; 777 tl:= ( dy + d ) / ( AA + dx ); 722 p: intersect_rec ptr; 778 t2:= ( dy - d ) / ( AA + dx ); 723 theta, ans: array[l..2] of real; 779 theta[l]:= 2.0 * arctan( tl ); 724 num, i: integer; 780 theta[2]: 2.0 * arctan( t2 ) 725 begin 781 end; 726 782 727 length:= sqrt( (maxpoint.xcoord - minpoint.xcoord) 783 num:= 0; 728 *(max_point.xcoord - minpoint.xcoord) 784 case C2.which_half of

Jul 8 02:16 1986 Page 15 Jul 8 02:16 1986 Page 16 785 down (take negative angles 1841 786 for i:= 1 to 2 do 842 while( curvel (> nil ) do 787 if theta~i] <= 0.0 then 843 begin 788 begin 844 P2:S2.layer-head; 789 num:= num, + 1; 845 flag2: upper; 790 ans~num]: theta[i]; 846 curve2: get curve( p2, flag2 ) 791 end; 847 792 up ( take positive angles 848 while curve2 (> nil ) do 793 for i:= 1 to 2 do 849 begin 794 if theta~i] >= 0.0 then 850 case curvel> curve type of 795 begin 851 line: if curve2.curve ye=ln 796 num:= num + 1; 852 then 797 ans[num]: theta~i]; 853 line-line( curvl uv2 798 end; 854 Sl, S2,abpt 799 end; [ case 1855 else 800 856 line-arc( curvcue2 801 for i:= 1 to nun do 857 Sl, s2,abIpr) 802 begin 858 arc: if curve2. curve ye=ln 803 x: xc2 + cos( ans[i])* 859 then 804 C2.radius; 860 line-arc( curv2cuel 805 if ((a < (x/length) )and 861 Sl, S2, a,It 806 (x ( (b/length)) and 862 else 807 (x1l<x and x <x2))and 863 arc-arc( curvel uv2 808 (x3(<x) and x <x4 )then 864 31, s2, abIpr) 809 begin 865 end; [ case 810 new( p ) 866 811 p.intersection:= x/ length; 867 curve2:= get -curve( p2, fla2) 812 p~ next -intersect:=Iptr; 868 end; fcurve2 813 Iptr:= p; 869 814 end; 870 curvel:= get curve( p1, flagl ) 815 end; Ifor 1871 816 end; if) 872 end; Icurvel 817 end; I else 1873 818 874 end; Ielse 819 end;- arc-arc 1875 820 876 end; Iconpute-intersection 821 procedure compute-intersection( 31, S2: segment~ptr.; 877 822 var Iptr:intersect-rec~ptr; 878 procedure sort_intersection( var Iptr: intersect-rctr) 823 a, b: t value); 879 sort the intersection points is ascending order;8241 compute the intersection points of curves in 31 and S2; 880 linear sort is implemented here. 825 results are stored in the list pointed by Iptr.1 881 var 826 var 882 p, q, Il, 12: intersect rec~ptr; 827 p1, p2: layer~ptr,; 883 done: boolean; 828 flagl, flag2: lower-or upper; 884 begin 829 curvel, curve2: curve~ptr; 885 830 begin 886 if DEBUGINTERSECT then 831 887 begin 832 if ( S1 = nil ) or ( S2 =nil) then 888 writeln('#### In sort-intersection, begin dupItli') 833 Iptr:= nil 889 p:= Iptr; 834 else begin 890 while I p(> nil ) do 835 891 begin 836 Iptr:= nil; 892 writeln('#### In sort itrsciointersection,= 837 893, p-A.ntersection ) 838 p1:= S1 layer -head; 894 p:= p.next-intersect; 839 flagl:=upper;- 895 end; 840 curvel: get curve( pl, flagil) 896 writeln('#### In sort-intersection, end dumpprlit)

Jul 8 02:16 1986 Page 17 Jul 8 02:16 1986 Page 18 897 end; 953 q:= I1; 898 954 899 if Iptr <> nil then 955 while ( p <> nil ) do 900 begin 956 begin 901 p: Iptr; 957 if ( p. intersection <= I1A. intersection ) or 902 p:= p.next_intersect; 958 ( p.intersection = qA.intersection ) then 903 while( p <> nil ) do 959 begin 904 begin 960 p:= p-.next intersect; 905 q:=p; 961 q.next_intersect:= p; 906 if q.intersection <= Iptr-.intersection then 962 end 907 begin 963 else if p-.intersection >= I2^.intersection then 908 p:= p-.next_intersect; 964 begin 909 q.next_intersect:= Iptr; 965 p: nil; 910 Iptr:= q; 966 end 911 end 967 else begin 912 else begin 968 q:=p; 913 I2:= Iptr; 969 p:= p.next_intersect; 914 I1:= I2.next_intersect; 970 end; 915 done:= false; 971 end; [ while ] 916 while ( I1 <> nil ) and ( not done ) do 972 917 begin 973 q.next_intersect:= I2; 918 if I1.intersection >= 974 919 q-.intersection then 975 end; [ screen_intersection ] 920 done:= true 976 921 else begin 977 function comnpute_point( t: t_value; C curve_ptr; 922 I2:= I1; 978 xl, x2: t_value ):y_value; 923 I1:= Ii.next_intersect; 979 { given a parameter t, compute the value associated with t 924 end; 980 on the curve C; i.e. C(t) } 925 end; [ while ] 981 var 926 982 m, length: real; 927 p:= p-.next_intersect; 983 y: y_value; 928 q.next_intersect:= I1; 984 begin 929 I2.next_intersect:= q; 985 930 end; 986 case C-.curvetype of 931 end; [ while 1 987 line: begin 932 end; f if ] 988 m:= ( C-.curve_end - C-.curve_start ) / 933 989 ( x2 - xl ); 934 end; [ sort_intersection 990 y: C-.curve_start + m * ( t - xl ) 935 991 compute_point:= y 936 procedure screen_intersection( var Iptr: intersect_recptr; 992 end; 937 a, b: tvalue ); 993 arc: begin 938 1 insert a and b into the intersection list and make sure that 994 length:= sqrt( 939 this list is bounded by a and b, and there are no duplicate 995 (maxpoint.x_coord - minpoint.xcoord) 940 intersection points there. ] 996 * (maxpoint.x_coord - minpoint.xcoord) 941 var 997 + (maxpoint.ycoord - minpoint.y_coord) 942 p, q, I1, I2: intersect_rec ptr; 998 * (maxpoint.ycoord - minpoint.y_coord) 943 begin 999 + (maxpoint.z_coord - minpoint.z_coord) 944 1000 * (maxpoint. zcoord - minpoint.z_coord)); 945 new( I1 ); 1001 946 I1-.intersection:= a; 1002 m:= (t - C.center.x_coord) * 947 Il-.next_intersect:= Iptr; 1003 (t - C-.center.x_coord) * 948 new( I2 ); 1004 length * length; 949 I2.intersection:= b; 1005 m: C-.radius * C-.radius - m; 950 I2.nextintersect:= nil; 1006 m: sqrt(abs( m )); 951 p:= Iptr; 1007 if C^.which half = up then 952 Iptr:= I1; 1008 y:= C.center.ycoord + m

Jul 8 02:16 1986 Page 19 Jul 8 02:16 1986 Page 20 1009 else 1065 begin 1010 y:= C.center.y_coord - m; 1066 1011 computepoint:= y; 1067 [ must save q (i.e. copy list ) before operation 1012 end; 1068 p:= q; 1013 end; { case ] 1069 r:= nil; 1014 1070 while p <> nil do 1015 end; f computepoint 1 1071 begin 1016 1072 if r = nil then 1017 procedure copy( C, D: curveptr ); 1073 begin 1018 ( copy the curve pointed by C to D; 1 1074 new( q ); 1019 begin 1075 r:= q; 1020 1076 end 1021 D curve_type:= C curvetype; 1077 else begin 1022 case C curvetype of 1078 new( r onext_layer ); 1023 line: begin 1079 r:= r.next_layer; 1024 D.curve_start:= C.curvestart; 1080 end; 1025 D^.curve end:= C^.curve end; 1081 1026 end; 1082 r-.inner_curve:= fix_curve( p^inner curve, 1027 arc: begin 1083 low, upp, S ); 1028 D.curve_start:= C.curve-start; 1084 r.outer curve: fix-curve( p.outer curve, 1029 D curve_end:= C.curve_end; 1085 low, upp, S ); 1030 D.center.x_coord:= C.center.x_coord; 1086 r.next_layer:= nil; 1031 D^.center.y_coord:= C.centeroy_coord; 1087 1032 D_.radius:= CA.radius; 1088 p:= p^.next_layer; 1033 D^.whichhalf:= CA.whichhalf; 1089 end; f while ] 1034 end; 1090 1035 end; f case 1 1091 end; { fix bound ] 1036 1092 1037 end; f copy ] 1093 function distance_to_axis( C: curveptr; low, upp: tvalue/ 1038 1094 S: segmentptr): real; 1039 function fixcurve( C: curveptr; low, upp: t_value; 1095 compute the distance from the mid-point of the curve C 1040 S: segmentptr ): curve_ptr; 1096 to principal axis; if C = nil, (means no data) then a 1041 ( fix the start and end points of the curve C ] 1097 minimal distance is assumed. ] 1042 var 1098 const 1043 xl, x2: tvalue; 1099 min distance = -1o0; 1044 D: curveptr; 1100 begin 1045 begin 1101 1046 1102 if C = nil then 1047 new( D ); 1103 distance_to_axis:= min_distance 1048 copy( C, D ); 1104 else begin 1049 xl:= S. lower bound; 1105 1050 x2:= S.upper_bound; 1106 distancetoaxis:= computepoint( 1051 D.curve_start: compute_point( low, C, xl, x2 ); 1107 (low + upp) / 2.0, C, 1052 D.curve-end:= computepoint( upp, C, xl, x2 ); 1108 Sl1owerbound, S.upper-bound ) 1053 1109 end; 1054 fix curve:= D; 1110 1055 1111 end; f distance to axis ] 1056 end; f fix curve ] 1112 1057 1113 function null_layer( a, b: curve ptr ): boolean; 1058 procedure fixbound( var q: layerptr; low, upp: tvalue; 1114 { to test if curve a and b form a null layer 1059 S: segnent_ptr ); 1115 const 1060 { This is for direct copy of intervals ( segments S ). 1116 EPS = 0o0001; 1061 go thru the q chain to fix the strat and end points 1117 var 1062 of the curves. 1 1118 result: boolean; 1063 var 1119 begin 1064 p, r: layerptr; 1120

Jul 8 02:16 1986 Page 21 Jul 8 02:16 1986 Page 22 1121 result:= false; 1177 Sl, S2: segmentptr ); 1122 if ( abs (a-.curve_start - b'.curvestart) <= EPS ) and 1178 1 actually union of layers specified in S1 and S2; both Sl and 1123 ( abs (a'.curveend - b-.curveend) (= EPS ) and 1179 S2 are not nil; new layers are pointed by q and bounded by 1124 ( a.curve_type = bh.curvetype ) then 1180 low and upp; 1125 begin 1181 var 1126 case a-.curvetype of 1182 flagl, flag2: lower_orupper; 1127 line: result:= true; 1183 pl, p2, r: layer_ptr; 1128 arc: result:= (a-.radius = b-.radius) and 1184 C1, C2: curveptr; 1129 (a.which_half = b-.which_half) and 1185 nextstate: state; 1130 (abs(a^.center.xcoord - 1186 dl, d2: real; 1131 b.center.x_coord) (= EPS) and 1187 a, b: curve_ptr; 1132 ( abs(a.center.y_coord - 1188 begin 1133 b-.center.y_coord) <= EPS); 1189 1134 end; ( case ] 1190 flagl:= upper; 1135 end; 1191 flag2:= upper; 1136 null_layer:= result; 1192 p1:= S1.layerhead; 1137 1193 p2:= S2-.layerhead; 1138 end; [ null_layer ] 1194 Cl:= getcurve( p1, flagl ); 1139 1195 C2: getcurve( p2, flag2 ); 1140 procedure add-layer( upperrange, lower_range:curveptr; 1196 dl:= distance_toaxis( Cl, low, upp, S1 ); 1141 var q, r:layer_ptr ); 1197 d2:= distance_toaxis( C2, low, upp, S2 ); 1142 [ add a new layer bounded by [upperrange lowerrange] 1198 q:= nil; 1143 to the q chain; r always points to the head of the q chain. ] 1199 next_state:= minuszero; 1144 begin 1200 1145 1201 repeat 1146 if not nulllayer( upperrange, lower_range ) then 1202 case next_state of 1147 1 null layer will be removed ] 1203 minuszero: 1148 begin 1204 if dl >= d2 then [ + 1149 1205 begin 1150 if DEBUGADD_LAYER then writeln( 1206 a:= fix_curve( C1, low, upp, S1 ) 1151'$$$$ In addlayer, uppercurve(start,end),', 1207 next_state:= pluszero; 1152'lowercurve(start,end)=', 1208 Cl:= getcurve( pl, flagl ); 1153 upperrange.curve_start, upperrange-.curveend, 1209 dl:= distance_to_axis(Cl,low,upp,Sl); 1154 lower_range-.curve_start, lowerrange.curveend ); 1210 end 1155 1211 else 1156 if q = nil then the first one 1212 begin [ 1 1157 begin 1213 a:= fix_curve( C2, low, upp, S2 ); 1158 new( q ); 1214 next_state:= minusone; 1159 q.next_layer:= nil; 1215 C2 = get_curve( p2, flag2 ); 1160 q.outer curve:= upperrange; 1216 d2:= distancetoaxis(C2,low,upp,S2 ); 1161 q.innercurve:= lower_range; 1217 end; 1162 r:= q; 1218 pluszero: 1163 end 1219 if dl >= d2 then [ - 1164 else 1220 begin 1165 begin [ normal process 1 1221 b:= fixcurve( Cl, low, upp, S1 ) 1166 new( r-.nextlayer ); 1222 addlayer( a, b, q, r ); 1167 r:= r-.next_layer; 1223 next_state:= minuszero; 1168 r.outer_curve:= upperrange; 1224 Cl: get_curve( pl, flagl ); 1169 r-.inner curve:= lower_range; 1225 dl:= distance_to_axis(Cl,low,upp,Sl); 1170 r-.nextlayer:= nil; 1226 end 1171 end; 1227 else 1172 end; 1228 begin [ 1 ] 1173 1229 next_state:= plusone; 1174 end; [ addlayer ] 1230 C2:= getcurve( p2, flag2 ); 1175 1231 d2:= distance_to_axis(C2,low,upp,S2 ); 1176 procedure layerunion ( var q: layerptr; low, upp: t_value; 1232 end;

Jul 8 02:16 1986 Page 23 Jul 8 02:16 1986 Page 24 1233 plusone 1289 q: nil; 1234 if dl >= d2 then -1 1290 next state in inuszero; 1235 begin 1291 1236 next-state:=innusone; 1292 repeat 1237 Cl: get curve( p1, flagl ) 1293 1238 dl:=distance to axis(Cl,low,upp,Sl); 1294 if DEBUG_-LAYER_-DIFFERENCE then writeln( 1239 end 1295 ####in layer difference, dl, d2=', dd) 1240 else 1296 1241 begin f0 11297 case next-state of 1242 next state:= pluszero; 1298 minuszero 1243 C2 get-curve( p2, flag2 ) 1299 if dl >= d2 then + 1244 d2:=distance-to-axis(C2,low,upp,32 );1300 begin 1245 end; 1301 a: fix curve( Cl, low, upp31) 1246 minusone: 1302 next-state:= pluszero; 1247 if dl >= d2 then + 1303 Cl: get curve( p1, flagl ) 1248 begin 1304 dl distance-to-axis(C1,lo1 ppS) 1249 next-state:=plusone; 1305 end 1250 Cl: get- curve( p1, flagl ) 1306 else 1251 dl distance-to-axis(Cl,low,upp,Sl); 1307 begin 11 1252 end 1308 next state:= minusone; 1253 else 1309 C2: get-curve( p2, flag2 ) 1254 begin 10 11310 d2:=distance-to-axis(C2,lo1 p,2) 1255 b: fix-curve( C2, low, upp, 32 ) 1311 end; 1256 add layer( a, b, q, r ); 1312 pluszero 1257 next state:= mnunuzero; 1313 if dl >= d2 thenf1258 C2: get- curve( p2, flag2 ) 1314 begin 1259 d2: distance-to-axis(C2,low,upp,S2 ) 1315 b:= fix-curve( Cl, low, upp31) 1260 end; 1316 add_layer( a, b, q, r ); 1261 end; 1case 1 1317 next-state:= minuszero; 1262 until (Cl= nil ) and(C2 =nil); 1318 Cl1: get curve( p1, flagl ) 1263 1319 dl:distance-to-axis(Cl,lowupS) 1264 end;. ( layer-union 11320 end 1265 1321 else 1266 procedure layer-difference( var q: layer~ptr; 1322 begin 11 1267 low, upp: t~value; 1323 b: fix curve( C2, low, upp32) 1268 SI, 32: segmentjptr ) 1324 add layer( a, b, q, r ~ 1269 actually difference of layers specified in 31 and 32; 1325 next-state:= plusone; 1270 both 31 and 32 are not nil; new layers are pointed by q 1326 C2: get curve( p2, flag2 ) 1271 and bounded by low and upp; 1 1327 d2:distance-to-axis(C2,lo1 p,2) 1272 var 1328 end;1273 flagi, flag2 lower -or- upper; 1329 plusone 1274 p1, p2, r: layer~ptr;- 1330 if dl >= d2 thenf1275 Cl, C2: curve-ptr; 1331 begin 1276 next-state: state; 1332 next-state m=iinusone; 1277 dl, d2: real; 1333 Cl1: get curve( p1, flagl ) 1278 a, b: curve~ptr; 1334 dl: distance-to-axis(Cl,lowupS) 1279 begin 1335 end 1280 1336 else 1281 flagl:=upper; 1337 begin 101 1282 flag2: upper;- 1338 a fix curve( C2, low, upp32) 1283 p1:=31 layer-head; 1339 next-state:= pluszero; 1284 p2: 2> layer bhead;' 1340 C2 get curve( p2, flag2 ) 1285 Cl:get curve( p1, flagil) 1341 d2:distantce-to-axis(C2,lo1 p,2) 1286 C2:get curve( p2, flag2 1;1342 end;. 1287 dl: distance -to axis( Cl, low, upp, 31 ) 1343 mninusone 1288 d2: distance to axis( C2, low, upp, 32 ) 1344 if dl >= d2 then +I

Jul 8 02:16 1986 Page 25 Jul 8 02:16 1986 Page 26 1345 begin 1401,'low, upp=',low, upp); 1346 next_state:= plusone; 1402 1347 Cl:= getcurve( pl, flagl ); 1403 if SI = nil then 1348 dl:= distance_to_axis(Cl,low,upp,Sl); 1404 q:= nil 1349 end 1405 else if S2 <> nil then 1350 else 1406 layerdifference( q, low, upp, S, S2 ) 1351 begin { 0 ] 1407 else begin 1352 next_state:= minuszero; 1408 q:= Sl.layerhead; 1353 c2:= getcurve( p2, flag2 ); 1409 fixbound( q, low, upp, S ); 1354 d2:= distance_to_axis(C2,low,upp,S2 ); 1410 end; 1355 end; 1411 end; 1356 end; [ case ] 1412 end; { case ] 1357 until (C1 = nil ) and ( C2 = nil ); 1413 1358 1414 end; { computelayer ] 1359 if DEBUG_LAYER_DIFFERENCE then writeln( 1415 1360'**** In layer_difference, before fix_bound ****'); 1416 procedure addsegment( operation: axisop_kind; 1361 1417 a, b: t_value; var p: segmentjtr; 1362 end; [ layer_differenc ] 1418 S1, S2: segment_ptr ); 1363 1419 [ add a newly created segment [a b] into the principal axis P; 1364 procedure computelayer( operation: axis_op_kind; 1420 where p is the segment pointer. This routine should include 1365 var q: layer_ptr; low, upp: tvalue; 1421 interval intersection computation and layer merge operation.) 1366 S1, S2: segment_ptr ); 1422 var 1367 1 compute the new layers by Si operation S2 and bounded by 1423 Iptr: intersectrec_ptr; 1368 low and upp; results should be sorted according to depth. 1424 low, upp: tvalue; 1369 S1 or S2 may be nil, in this case direct copy the curve is 1425 begin 1370 possible which depends on what operation is applied. New 1426 1371 layer is pointed by q. ] 1427 if a >= b then 1372 var 1428 begin 1373 S: segment_ptr; 1429 { return ] 1 null interval; do nothing ] 1374 begin 1430 end 1375 1431 else begin 1376 case operation of 1432 1377 union: begin 1433 computeintersection( Si, S2, Iptr, a, b ); 1378 1434 sort_intersection( Iptr ); 1379 if DEBUG_COMPUTE_LAYER then writeln( 1435 screen_intersection( Iptr, a, b ) 1380'$$$$ computelayer, union op, low, upp=' 1436 1381,low, upp); 1437 low:= Iptr.intersection; 1382 1438 Iptr:= Iptr.next_intersect; 1383 if ( S1 <> nil ) and ( S2 <> nil ) then 1439 1384 layerunion( q, low, upp, S1, S2 ) 1440 repeat 1385 else if ( S1 = nil ) and ( S2 = nil ) then 1441 upp:= IptrA.intersection; 1386 q:= nil 1442 p.lower_bound: low; 1387 else begin 1443 p upper_bound: upp; 1388 if Si = nil then 1444 computelayer(operation, p^.layerhead, 1389 S:= S2 1445 low, upp,Sl,S2); 1390 else if S2 = nil then 1446 if p.layerhead <> nil then 1391 S:= S1; 1447 begin 1392 q:= S.layerhead; 1448 new( p^.nextsegment ); 1393 fix_bound( q, low, upp, S ); 1449 p:= p nextsegment; 1394 end; 1450 p.next_segment:= nil; 1395 end; 1451 end; 1396 difference: 1452 1397 begin 1453 low:= upp; 1398 1454 Iptr:= Iptr.next_intersect; 1399 if DEBUGCOMPUTE LAYER then writeln( 1455 until Iptr = nil; 1400'$$$$ computelayer, difference op,' 1456

Jul 8 02:16 1986 Page 27 Jul 8 02:16 1986 Page 28 1457 end; 1513 basically it performs the merge step of merge_sort algorithm.] 1458 1514 var 1459 end; f add_segment ] 1515 PX: principal_axis; 1460 1516 S1, S2, NS1, NS2: segmentptr; 1461 function compute D3point( ptl, pt2: D3_point; 1517 flagl, flag2: lower_or_upper; 1462 t: t_value ): D3point; 1518 donel, done2: boolean; 1463 1 given two end points ptl and pt2, and a parameter t, 1519 a, b, al, a2, actual_t, first_t, last_t: t_value; 1464 compute its coordinate 1 1520 ptl, pt2: D3_point; 1465 begin 1521 p: segment_ptr; 1466 1522 next_state: state; 1467 compute D3_oint.x coord:= ptl.x_coord + 1523 begin 1468 t * (pt2.x_coord-ptl.x_coord); 1524 1469 compute D3point.y_coord:= ptl.y_coord + 1525 { need code to initialize PX.start_point and PX.endpoint } 1470 t * (pt2.y_coord-ptl.y_coord); 1526 PX.start_point:= min_point; 1471 compute D3_point.z coord:= ptl.z_coord + 1527 PX.endpoint:= max_point; 1472 t * (pt2.z_coord-ptl.z_coord); 1528 1473 1529 new( PX.segment_head ); 1474 end; f compute D3_point ] 1530 p:= PX.segment_head; 1475 1531 p.next_segment:= nil; 1476 function get next( var flag: lower_or_upper; 1532 NS1:= Pl.segment_head; 1477 var done: boolean; 1533 NS2:= P2.segment_head; 1478 var S, NS: segmentptr ): t_value; 1534 flagl:= lower; 1479 ( get next interval value for use, done is set while no more 1535 flag2:= lower; 1480 can be obtained. S should always point to the interval 1536 donel = false; 1481 currently dealt with. } 1537 done2:= false; 1482 const 1538 1483 max t value = 2.0; 1539 al:= get_next( flagl, donel, Sl, NS1 ); 1484 begin 1540 a2:= get_next( flag2, done2, S2, NS2 ); 1485 1541 next_state:= minuszero; 1486 if done then 1542 1487 get_next:= max t value 1543 repeat 1488 else begin 1544 case nextstate of 1489 case flag of 1545 minuszero: 1490 lower: begin 1546 if al <= a2 then 1491 S:= NS; 1547 begin 1492 flag:= upper; 1548 a:= al; 1493 if S (> nil then 1549 next_state:= pluszero; 1494 get_next:= S^.lower_bound 1550 al:= get_next( flagl, donel, Sl, NSl ); 1495 else begin 1551 end 1496 done:= true; 1552 else begin 1497 get_next:= max t value; 1553 a:= a2; 1498 end; 1554 nextstate:= minusone; 1499 end; 1555 a2:= get_next( flag2, done2, S2, NS2 ); 1500 upper: begin 1556 end, 1501 flag:= lower; 1557 pluszero: 1502 NS:= S next segment; 1558 if al <= a2 then 1503 get_next:= S.upper_bound 1559 begin 1504 end; 1560 b:= al; 1505 end; f case ] 1561 next state:= minuszero; 1506 end; 1562 add_segment(operation, a, b, p, Sl,nil); 1507 1563 al:= get_next( flagl, donel, S1, NS1 ); 1508 end; i get_next ) 1564 end 1509 1565 else begin 1510 function partition ( operation: axis_op_kind; 1566 b:= a2; 1511 P1, P2: principal_axis ): principal axis; 1567 next_state:= plusone; 1512 1 compute new intervals for union/difference operation; 1568 add_segment(operation, a, b, p, S1, nil);

Jul 8 02:16 1986 Page 29 Jul 8 02:16 1986 Page 30 1569 a:= a2; 1625 actual_t:= p^.upperbound; 1570 a2:= getnext( flag2, done2, S2, NS2 ); 1626 end; 1571 end; 1627 p:= p.next_segment; 1572 plusone: 1628 end; 1573 if al <= a2 then 1629 1574 begin 1630 first_t:= PX.segmenthead.lower_bound; 1575 b:= al; 1631 lastt:= actualt; 1576 next_state:= minusone; 1632 1577 addsegment(operation, a, b, p, S1, S2); 1633 [ if operation = difference, should reshape the t 1578 a:= al; 1634 parameters because the first and last several segments 1579 al:= getnext( flagl, donel, S1, NSl ); 1635 might be removed. recalucate the end point and the t 1580 end 1636 parameters. ) 1581 else begin 1637 if (operation = difference ) and 1582 b:= a2; 1638 (( lastt < 1.0 ) or ( first_t > 0.0 )) then 1583 next_state:= pluszero; 1639 begin 1584 addsegment(operation, a, b, p, Sl, S2); 1640 ptl:= PX.startpoint; 1585 a:= a2; 1641 pt2:= PX.endpoint; 1586 a2: get-next( flag2, done2, S2, NS2 ); 1642 PX.start_point:= computeD3_point( 1587 end; 1643 ptl, pt2, first_t ); 1588 minusone: 1644 PX.endpoint:= computeD3_point( 1589 if al <= a2 then 1645 ptl, pt2, lastt ); 1590 begin 1646 adjust-t-value( PX, firstt, lastt, 0.0, 1.0 ); 1591 b:= al; 1647 end; 1592 nextstate:= plusone; 1648 1593 if operation = union then 1649 partition: PX 1594 addsegment( operation, a, b, 1650 1595 p, nil, S2 ); 1651 if DEBUG_PARTITION then writeln( 1596 a:= al; 1652'%%%%% In partition, PX.start and end point=', 1597 al:= get _next( flagl, donel, Sl, NSl ); 1653 PX.start_point.x_coord, PX.startpoint.y_coord, 1598 end 1654 PX.start_point.z_coord, PX.endpoint.xcoord, 1599 else begin 1655 PX.endpoint.ycoord,PX.endpoint.z_coord); 1600 b:= a2; 1656 1601 nextstate:= minuszero; 1657 end; { partition 3 1602 if operation = union then 1658 1603 addsegment( operation, a, b, p, 1659 function smooth_curve( xl, x2: t_value; C1: curveptr 1604 nil, S2 ), 1660 x3, x4: t_value; C2: curveptr ): boolean 1605 a:= a2; 1661 f check to see if the curve Cl and C2 are smoothly connected. 1606 a2:= get next( flag2, done2, S2, NS2 ); 1662 note that Cl and C2 have the same curve type here. 1607 end; 1663 const 1608 end; [ case ] 1664 EPSILON = 0.00001; 1609 until donel and done2; 1665 INFINITY = 99999.99999; 1610 1666 var 1611 { need code to delete the last element of the segment 1667 ml, m2: real; 1612 ( null segment ) ) 1668 begin 1613 p:= PX.segment_head; 1669 1614 if p-.nextsegment = nil then 1670 case Cl^.curve_type of 1615 begin 1671 line: begin { check if they have the same slope 1616 PX.segment_head:= nil; 1672 if abs(xl -x2) <= EPSILON then 1617 actualt: 0.0; 1673 ml:= INFINITY 1618 end; 1674 else ml: ( Cl^.curve_start - 1619 1675 Cl.curveend )/( xl - x2 ); 1620 while ( p <> nil ) do 1676 if abs(x3 -x4) <= EPSILON then 1621 begin 1677 m2:= INFINITY 1622 if p-.next_segment.nextsegment = nil then 1678 else m2:= ( C2-.curve_start - 1623 begin 1679 C2-.curve_end )/( x3 - x4 ); 1624 p.next_segment:= nil; 1680 if abs(ml - m2) <= EPSILON then

Jul 8 02:16 1986 Page 231 Jul 8 02:16 1986 Page 32 1681 smooth-curve:= true 1737 begin 1682 else smooth-curve:=false 1738 a q l1ayer- head; 1683 end;, 1739 b:r.layer- head; 1684 arc: begin 1740 done:= false; 1685 if Cl. radius () C2-.radius then 1741 while not done do 1686 smooth-curve:= false 1742 begin 1687 else if Cl~.which half <> C2.which half 1743 if (a~ inner curve.curve type 0 1688 then smooth-curve:= false 1744 b- inner curve-.curve type) or 1689 else if ( Cl.center.x-coord 0 1745 (a~ outer curve.curve type (> 1690 C2.center.x coord ) or 1746 b.outer curve curve type) 1691 ( Cl.center.y coord <> 1747 then begin 1692 C2.center.y coord ) then 1748 done:= true; 1693 smooth-curve:=false 1749 check-connect:false 1694 else smooth-curve:=true;1750 end 1695 end; 1751 else if (a-.inner-curv&.curve-end ( 1696 end; I case 11752 b-'inner-curve. curve-start)o 1697 1753 (a.outer curve-. curve end ( 1698 end; ( smooth-curve 11754 b-outer-curve. curve-start)he 1699 1755 begin 1700 function check-smooth( q, r segment~ptr; 1756 done:= true;1701 a, b layerjptr): boolean; 1757 check connect: false 1702( within segments q and r, check if layer a and b 1758 end 1703 are smoothly connected. 1759 else if not check -smooth( q,, a, b )te 1704 var 1760 begin 1705 flag: boolean; 1761 done:= true; 1706 xl, x2, x3, x4: t-value; 1762 check-connect: false 1707 begin 1763 end 1708 1764 elsef ok for this layer, let's try nx 1709 xl:=q lower -bound; 1765 begin 1710 x2:=q -upper-bound; 1766 a:a~ next_layer; 1711 x3:=r lower bound; 1767 b b.next -layer; 1712 x4: r upper-bound; 1768 if( a = nil ) and ( b <>nl)o 1713 1769( b = nil ) and ( a <>nl)te 1714 flag:= smooth-curve( xl, x2, a inner-curve, 1770 begin 1715 x3, x4, b -inner -curve ) 1771 done:= true; 1716 if flag then 1772 check-connect: as 1717 flag:= smooth-curve( xl, x2, a outer-curve, 1773 end 1718 x3, x4, b.outer-curve ) 1774 else if (a = nil) and (b=ni)te 1719 check-smooth:flag;1775 begin 1720 1776 done:= true; 1721 end; f check smooth 11777 check connect =tu 1722 1778 end; 1723 function check connect( q, r: segment~ptr): boolean; 1779 end;1724 To check if q and r are connected; check the following 1780 end;[ while 1725 condition: 1781 end; [ else 1726 the same connecting point, the same curve type, 1782 1727 smooth at the connecting point, 1783 end; fcheck-connect 1728 the same number of layers. 1784 1729 var 1785 procedure merge -curve( Cl, C2: curve~ptr ) 1730 a, b: layerjptr; 1786 1 merge the curve Cl and C2; result is in Cl; 1731 done:boolean; 1787 begin 1732 begin 1788 1733 1789 case C1.curve type of 1734 if q~ upper -bound <> r~ lower-bound then 1790 line: Cl. curve-end: C2. curve-end;1735 check-connect:= false 1791 arc: Cl curve-end C2~ curve-end; 1736 else 1792 end; 0- case I

Jul 8 02:16 1986 Page 33 Jul 8 02:16 1986 Page 34 1793 1849 P.end_point.y_coord,P.endpoint.zcoord); 1794 end; f merge_curve ] 1850 1795 1851 end; 1 mergeinterval } 1796 procedure connect( q, r: segmentptr ); 1852 1797 ( merge q and r into a new interval pointed by q. 1853 function axis operation ( operation: axisop kind; 1798 1854 P1, P2: principal_axis ): principal axis; 1799 var 1855 ( union/difference two axes P1 and P2 into a new axis; 1800 p,s: layerptr; 1856 this operation performs the actual union/difference of 1801 begin 1857 two CSG trees; it first calls normalization procedure 1802 1858 to make two trees have t-parameters on the same basis, 1803 { merge segment ] 1859 then applies the merge routine of merge-sort algorithm 1804 q^.upperbound:= r^.upperbound; 1860 to union/difference them together. ] 1805 q.next_segment:= r.next_segment; 1861 var 1806 1862 P: principal_axis; 1807 [ merge each layer ] 1863 begin 1808 p:= q layerhead; 1864 1809 s:= r layerhead; 1865 if DEBUG_AXIS_OP then writeln('%%%% enter normalization'); 1810 while ( p <> nil ) do 1866 1811 begin 1867 normalization ( P1, P2, minpoint, max_point ); 1812 mergecurve( p.inner_curve, s^.innercurve ); 1868 1813 mergecurve( p.outer_curve, s.outer curve ); 1869 if DEBUG_AXIS_OP then writeln('%%%% leave normalization'); 1814 1870 1815 p:= p.next_layer; 1871 if DEBUG_AXIS_OP then writeln('%%%% enter partition'); 1816 s:= s^.nextlayer; 1872 1817 end; [ while ] 1873 P:= partition ( operation, P1, P2 ) 1818 1874 1819 end; ( connect 1875 if DEBUGAXIS_OP then writeln('%%%% leave partition'); 1820 1876 1821 procedure merge_interval ( var P: principalaxis ); 1877 if TRACE_AXIS then dump_axis( P ); 1822 1 try to merge consecutive intervals into a larger one; 1878 1823 this step is important to enforce the uniqueness of 1879 if DEBUG_AXIS_OP then writeln('%%%% enter merge_interval'), 1824 this principal-axis representation. I 1880 1825 var 1881 merge_interval ( P ) 1826 q, r: segmentptr; 1882 1827 begin 1883 if DEBUG_AXIS_OP then writeln('%%%% leave merge_interval'); 1828 1884 1829 q:= P.segmenthead; 1885 if TRACEAXIS then dumnp axis( P ); 1830 r:= q next_segment; 1886 1831 while ( r <> nil ) do 1887 axisoperation:= P; 1832 begin 1888 1833 if checkconnect( q, r ) then 1889 end; ( axisoperation 1 1834 begin 1890 1835 connect( q, r ); 1891 procedure rotx( angle: real; var p: D3_point; 1836 r:= q next_segment; 1892 p0: D3point ); 1837 end 1893 [ rotate the point p about x-axis by amount of angle; 1838 else 1894 since the rotation is relative to local coordinate system 1839 begin 1895 of its own, translation must be done before and after 1840 q:= r; 1896 rotation. ) 1841 r:= q^.next_segment; 1897 var 1842 end; 1898 y, z: real; 1843 end; [ end ) 1899 begin 1844 1900 1845 if DEBUG MERGE INTERVAL then writeln( 1901 [ translation to the origin of local coordinate ] 1846'%%%%% In merge_interval,P.start and end point=', 1902 p.x_coord:= p.x_coord - pO.x_coord; 1847 P.start_point.xcoord, P.startpoint.y_coord, 1903 p.y_coord:= p.y_coord - pO.y_coord; 1848 P. startpoint.zcoord, P.end_point.x_coord, 1904 p.z_coord: p.z_coord - pO.z_coord;

Jul 8 02:16 1986 Page 35 Jul 8 02:16 1986 Page 36 1905 1961 p~ycoord =p~y~coord - pO~ycoord 1906 frotate 11962 pz-coord p-z coord - pO.z~coord 1907 y:p.y_coord; 1963 1908 z:p. z coord; 1964 frotate 1909 angle:= angle * PI / 180.0 f convert to radian degree 1965 x: p.x_coord; 1910 p.y coord: y * cos( angle)- z * sin( angle ) 1966 y: p.y coord1911 p.z_coord: y * sin( angle ) 1 z * cos( angle ) 1967 angle angle A PI / 180.0 f convert to radiandge 1912 1968 p.x~coord: x * cos ( angle) - y * sin( angle ) 1913 f translation back to the global coordinate 1969 p.y_coord: x * sin( angle ) 1 y * cos( angle ) 1914 p.x_coord: p.x_coord + pO.x_coord; 1970 1915 p.y_coord:=p.y~coord + pO.y_coord; 1971 (translation back to the global coordinate 1916 p.z_coord p.z _coord + pO.z_coord;1972 p.x-coord p.x coord + pO~x coord 1917 1973 p.y_coord:=p.y~coord + pO.y_coord 1918 end; f rot-x 11974 p.z_coord:=p.z coord + pO.z_coord 1919 1975 1920 procedure rot~y( angle: real; var p: D3_point; 1976 end; f rot-z 1921 P0: D3_point ); 1977 1922 f rotate the point p about y-axis by amount of angle, 1978 procedure translate( xform: xform matrix; var p:Dpit) 1923 since the rotation is relative to local coordinate system 1979 ftranslate the point p by the amount of the translainpr 1924 of its own, translation must be done before and after 1980 of Xform 1925 rotation. 11981 begin 1926 var 1982 1927 x, z: real; 1983 p.x~coord: p.x_coord + xform~translate_-x; 1928 begin 1984 p.y coord:p~y coord + xform.translatey; 1929 1985 p.z~coord:=p.z_coord +I xform.translatez; 1.930 (translation to the origin of local coordinate 11986 1931 p.x_coord: p.x coord - pO.x_coord; 1987 end; I translate 1932 p.y_coord =p~y~coord - p0O.y_coord;1988 1933 p.z_coord:=p.z coord - pO.z_coord;1989 procedure transform ( xform:xform matrix; 1934 1990 var P:principal axis ) 1935 1rotate 11991 1transform the tree by a transformation matrix; 1936 x:=p.x_coord; 1992 actually the principal axis is transformed; 1937 z p.z_coord; 1993 more specifically only the two end points are 1938 angle:=angle API / 180.0;Iconvert to radian degree 11994 transformed. 1939 p.x_coord:x Acos( angle )+ z A sin( angle ); 1995 1940 p.z_coord:=-x A sin( angle ) + z A cos( angle );1996 begin 1941 1997 1942 1translation back to the global coordinate 11998 f P.startjpoint P.start~point A transform mati; 1943 p.x_coord:=p.x coord + pO.x_coord;1999 1 P. end~point P~end~point A transform mati; 1944 p.y coord p.y coord + pO.y coord;2000 [ note other variables are transformation invarian 1945 p.z_coord:=p.z coord + pO.z_coord,2001 1946 2002 1let's do the rotation first 1947 end; I rot~y 12003 if xform.rotate -x (> 0.0 then 1948 2004 rot-x( xform.rotate-x, P.end~point, P.startpon) 1949 procedure rot-z( angle: real;. var p: D3_point; 2005 1950 P0:. D3_point );- 2006 if xform.rotatey (> 0.'0 then 1951 1rotate the point p about z-axis by amount of angle, 2007 rot~y( xform.rotate-y, P.end~point, P.startpin) 1952 since the rotation is relative to local coordinate system 2008 1953 of its own, translation must be done before and after 2009 if xform.rotate z 0 0.0 then 1954 rotation. )2010 rot-z( xform.rotate-z, P.end~point, P. startjon) 1955 var 2011 1956 x, y:real; 2012 1and then do the translation 1957 begin 2013 if ( xform-translate -x 0 0.0 ) or 1958 2014 ( xform.translatejy 0 0.0 ) or 1959 1translation to the origin of local coordinate 12015 ( xform.translate-z (> 0.0 ) then 1960 p.x_coord:= p.x_coord - pO.x_coord; 2016 begin

Jul 8 02:16 1986 Page 37 Jul 8 02:16 1986 Page 38 2017 translate( xform, P.start_point ); 2073 1.inner curve.curveend:=.0 2018 translate( xform, P.endpoint ); 2074 end; 2019 end; 2075 2020 2076 cone 2021 end; ( transform ] 2077 begin 2022 2078 l^.outercurve^.curve_type:= line 2023 function build_axis ( T: CSGtree ): principalaxis; 2079 1.outer_curve.curvestart = 2024 { convert the primitive solid T to the principalaxis 2080 T.primsolid.radiusb; 2025 type representation ] 2081 l^.outer_curve^.curve_end:= 2026 var 2082 T.primsolid.radiust; 2027 P: principalaxis; 2083 2028 s: segmentptr; 2084 1 inner_curve.curve_type:= line; 2029 1: layerptr; 2085 1. innercurve.curvestart: 0.0 2030 height: real; 2086 1.innercurve.curveend 0.0; 2031 begin 2087 end; 2032 2088 2033 P.startpoint:= T.prim_solid.base; 2089 torus 2034 P.end_point:= T.prim_solid.top; 2090 begin 2035 2091 l1.outer_curve^.curve_type:= arc; 2036 if DEBUG_BUILDAXIS then writeln( 2092 height:= ( T.primsolid.innerradius + 2037'**** In buildaxis,P.start and end point=', 2093 T.prim_solid.outerradius ) / 2.0; 2038 P.startpoint.x_coord, P.start_point.y_coord, 2094 l1.outercurve^.curve_start: height; 2039 P.start_point.z_coord, P.end_point.x coord, 2095 1.outercurve^.curveend:= height; 2040 P.endpoint.y_coord, P.endpoint.zcoord); 2096 1.outercurve^.center.xcoord:= 0.5; 2041 2097 l1.outercurve^.center.y_coord:= height; 2042 if DEBUG_BUILD_AXIS then writeln( 2098 1.outer_curve.radius:= 2043'**** In build_axis,Prim_solid.base and top point=', 2099 T.prim_solid.outer_radius - height; 2044 T.primsolid.base.x_coord,T.primsolid.base.ycoord, 2100 l1.outer_curve^.which_half:= up; 2045 T.primsolid.base.zcoord, T.primsolid.top.xcoord, 2101 2046 T.primsolid.top.y_coord,T.prim_solid.top.z_coord); 2102 l1.innercurve^.curve_type:= arc; 2047 2103 1^.inner curve^.curve start:= height; 2048 new( P.segment_head ); 2104 1^ innercurve-.curve_end:= height; 2049 2105 1.innercurve^.center.xcoord:= 0.5; 2050 s:= P.segment_head; 2106 1.innercurve^.center.ycoord:= height; 2051 s^. lower_bound:= 0.0; 2107 1.inner_curve ^.radius:= height - 2052 sa upper bound:= 1.0; 2108 T.primsolid.inner_radius; 2053 s^.nextsegment:= nil; 2109 l^.innercurve^.which_half: down; 2054 new( s^.layer_head ); 2110 end; 2055 2111 2056 1:= s^ layer_head; 2112 end; { case ] 2057 1.nextlayer:= nil; 2113 2058 new( l^.inner_curve ); 2114 build_axis: P; 2059 new( ^1 outer_curve ); 2115 2060 2116 end; { buildaxis ] 2061 case T.prim_solid.solid_type of 2117 2062 2118 function evaluateCSG ( T: CSGtreeptr ): principal axis; 2063 cylinder: 2119 recursive call itself to evaluate the CSG tree T; 2064 begin 2120 it is basically a tree traversal algorithm. 2065 1.outer_curve^.curve_type:= line; 2121 Our assumption is that the two operands at this point must 2066 1 outercurve.curvestart: 2122 be coaxial although they might not be coaxial at the lower 2067 T.prim_solid.radius; 2123 part of the tree; ] 2068 1-.outer_curve^.curve_end:= 2124 var 2069 T.prim_solid.radius; 2125 P, P1, P2: principal_axis; 2070 2126 begin 2071 l^.innercurve^.curve_type:= line; 2127 2072 1.innercurve.curve_start:= 0.0; 2128 if DEBUG_EVALUATE_CSG then writeln(

Jul 8 02:16 1986 Page 39 Jul 8 02:16 1986 Page 40 2129 in evaluate csg, partno=' T.partno ) 2185 P:= build axis ( T ) 2130 2186 evaluateCSG: P 2131 case T node -type of 2187 2132 movement 2188 if DEBUGEVALUATECSG then writeln( 2133 begin 2189'*** finish evaluate csg, partn 2134 P: evaluateCSG ( T child ) 2190 T-.partno ) 2135 transform( T-.move, P ) 2191 2136 evaluateCSG:P; 2192 end; 2137 2193 end,; case 2138 if DEBUG_ EVALUATE_ CSG then writeln( 2194 2139'*** finish evaluate csg, partno=, 2195 end;- ( evaluateCSG 2140 T-.partno );2196 2141 2197 procedure get input -data( var CSG: CSG -tree~ptr ) 2142 end; 2198( get the input CSG tree f rom outside world; it may b 2143 union-op: 2199 f rom the output of an interactive graphical geometia 2144 begin 2200 modeller; EECS487 term project, for example. 2145 P1: evaluate CSG (T left child );2201 const 2146 P2 evaluate_CSG( T.right child ) 2202 obj~path ='object' 2147 2203 mov-path ='movement;. 2148 if DEBUG_ EVALUATECSG then writeln( 2204 prim~path ='primitives', 2149'##it axis operation for partno=' 2205 2150 T partno );2206 var 2151 2207 openstatus: status$t 2152 P:= axis -operation ( union, P1, P2 ) 2208 2153 evaluateCSG:P; 2209 procedure read -tree( T: CSG-treejptr ) 2154 2210 1 recurvely read data from object, primitives,an 2155 if DEBUG_ EVALUATE_ CSG then writeln( 2211 movement files to construct a CSG tree T 2156'~** finish evaluate-csg, partno, 2212 var 2157 T-.partno ) 2213 partl, part2, part3: integer; 2158 2214 ch:char; 2159 if DEBUG_ EVALUATECSG then writeln( 2215 op arrayfl..9] of char; 2160'H!! result of union,P.start and end point=', 2216 prm: arrayfl..8] of char; 2161 P.startypoint.x -coord,P.start-point.y-coord, 2217 ax, ay, az, x, y, z, xl, yl, zl: real; 2162 P.start~point.z-coord, P.end~point.x-coord, 2218 2163 P.end~point.y~coord,P.endjpoint. zcoord); 2219 begin 2164 2220 2165 end; 2221 readln( obj file, partlch, ch, ch, ch, ch, ch, 2166 diff op: 2222 op, part2, part3 ) 2167 begin 2223 2168 P1I: evaluate_-CSG ( T left-child ) 2224 if DEBUG_ INPUT then writeln( 2169 P2: evaluate_ CSG ( T right child ) 2225 read object, objno =',partl, op), 2170 2226 2171 if DEBUG_ EVALUATECSG then writeln( 2227 if (T.partno C> CSGROOT) and 2172'liii axis operation for partno " 2228 ( T.partno () partl ) then 2173 T partno ) 2229 writeln(' —— > read -tree error,', 2174 2230' input tree structucture incorrect!) 2175 P:= axis -operation ( difference, P1, P2 ) 2231 2176 evaluateCSG:= P; 2232 if op =' union' then T-.node -type:uino 2177 2233 else if op =' differ' then T.node type: ifo 2178 if DEBUG_-,EVALUATECSG then writeln( 2234 else if op ='primitive' then T.node type: rmtv 2179 O** finish evaluate-csg, partno, 2235 else if op ='moved -obj' then T> node -type:=mven 2180 T-.partno ) 2236 else writeln('illegal CSG node type!!); 2181 2237 2182 end;. 2238 1 T-.partno:= partl; 2183 primitive: 2239 2184 begin 2240 case T.node type of

rul 8 02:16 1986 Page 41 Jul 8 02:16 1986 Page 42 2241 unionop, 2297 T^.primsolid.top.z coord:= z+xl; 2242 diffop: 2298 T^.prim_solid.radius_b:= yl; 2243 begin 2299 T^.prim_solidradius_t:= 0.0; 2244 T-rmsi2300 end; 2245 new( TA.leftchild ); 2301 torus: 2246 T. left child.partno:= part2; 2302 begin 2247 readtree( T A.leftchild ); 2303 T^.prim_solid.base.x_coord:= x 2248 2304 T^.prim_solid.base.y_coord:= y 2249 new( T right_child ); 2305 T prim_solid.base.z_coord:= z-yl; 2250 TA.rightchild^.partno:= part3, 2306 T^.prim_solid.top.x_coord:= x; 2251 readtree( T right_child ); 2307 T primsolid.top.y_coord:= y; 2252 2308 T^.prim_solid.top.zcoord:= z+yl; 2253 end; 2309 T.prim_solid.inner_radius:=xl-yli 2254 primitive: 2310 T^.prim_solid.outer_radius:=xl+yl; 2255 begin 2311 end; 2256 readln( primfile, partl, ch, 2312 end; [ case 1 2257 ch, ch, ch, ch, ch, prm, 2313 2258 xl, yl, zl, x, y, z, 2314 if ax <> 0.0 then 2259 ax, ay, az ); 2315 rotx( ax, T.primsolid.top, 2260 2316 T.primsolid.base ); 2261 if DEBUG_INPUT then writeln( 2317 2262'**** read primitive, no =', 2318 if ay 0 0.0 then 2263 partl,'',prm); 2319 rot_y( ay, T primsolid.top, 2264 2320 T.primsolid.base ); 2265 if ( partl <> T.partno ) then 2321 2266 writeln(' —— > read_tree error,', 2322 if az > 0.0 then 2267'input tree structucture incorrect!!'); 2323 rotz( az, T.primsolid.top, 2268 2324 T.primsolid.base ); 2269 f must treat cone specially 1 2325 2270 if prm ='cylinder' then 2326 end; 2271 T prim_solid.solidtype:= cylinder 2327 movement: 2272 else if prm =' cone' then 2328 begin 2273 T.prim_solid.solid_type: cone 2329 readln( movfile, partl, 2274 else if prm =' torus' then 2330 T move.translatex, 2275 TA.prim_solid.solid_type:= torus 2331 T.move.translate_y, 2276 else writeln('illegal CSG', 2332 T.move.translatez, 2277'primitive node type!!'); 2333 T move.rotate_x, 2278 2334 T^.move.rotate_y, 2279 case T.prim_solid.solid_type of 2335 T.move.rotate_z ); 2280 cylinder: 2336 2281 begin 2337 if DEBUG_INPUT then writeln( 2282 T prim_solid.base.x_coord:= x; 2338'**** read movement, no =', partl); 2283 T^.primsolid.base.y_coord:= y; 2339 2284 T^.prim_solid.base.zcoord:= z; 2340 if ( partl <> part3 ) then 2285 T.prim_solid.top.x_coord:= x; 2341 writeln(' —— > read_tree error,', 2286 T.prim_solid.top.y_coord:= y; 2342'input tree structucture incorrect!!'); 2287 T^.primsolid.top.zcoord:= z+xl; 2343 2288 T.primsolid.radius:= yl; 2344 new( T^.child ); 2289 end; 2345 T.child^.partno:= part2; 2290 cone: 2346 read tree( T child ); 2291 begin 2347 end; 2292 T prim_solid.base.x_coord:= x; 2348 end; f case ] 2293 T^.prim_solid.base.y_coord:= y; 2349 2294 T.prim_solid.base.zcoord:= z; 2350 end; ( read tree ] 2295 T^.prim_solid.top.x_coord:= x; 2351 2296 T^ primsolid.top.y_coord:= y; 2352 begin

Jul 8 02:16 1986 Page 43 Jul 8 02:16 1986 Page 44 2353 2409 pAupperbound ) then 2354 [ read data and construct them as a CSG tree; 2410 height:= C.center.ycoord + 2355 there are three data files; one is the movement, 2411 C.radius 2356 another is primitives, the third one is the 2412 else if C^.center.xcoord >= 2357 construction 1 2413 p^.upperbound then 2358 open( obj_file, obj_path,'OLD', openstatus.all );2414 height:= C.curveend 2359 open( mov file, movpath,'OLD', openstatus.all );2415 else 2360 open( primfile, primpath,'OLD', openstatus.all );2416 height:= C^.curve start; 2361 2417 end 2362 reset( obj_file );2418 else ( down 2363 reset( mov file );2419 begin 2364 reset( primfile );2420 if C.curvestart > C^.curveend then 2365 2421 height:= C^ocurve_start 2366 if DEBUG_INPPUT then writeln( 2422 else height:= C.curve_end; 2367'*** open files in getinputdata'); 2423 end; 2368 2424 end; [ case ] 2369 new( CSG );2425 2370 CSG.partno:= CSGROOT; 2426 if height > localheight then 2371 read_tree( CSG );2427 local_height:= height; 2372 2428 2373 if DEBUG_INPUT then writeln( 2429 len:= len + p.upperbound - p.lower_bound; 2374'*** finish get_inputdata'); 2430 q:= p.next_segment; 2375 2431 if q <> nil then 2376 end; [ get_inputdata 2432 begin 2377 2433 if ( q lower_bound - p^.upperbound ) > 00001 2378 procedure compute_length_diameter( PX: principal_axis; 2434 then 2379 var L, D: real ) 2435 begin 1 non-continuous segments ] 2380 { compute the maximal length L of PX; note: 2436 if len > maxlen then 2381 PX might be broken into several segments. Compute the 2437 begin 2382 extremal boundary (diameter D) from the axis 2438 maxlen:= len; 2383 var 2439 maxht: local_height; 2384 p, q: segmentptr; 2440 end; 2385 len, maxlen: t_value; 2441 len:= 0.0; 2386 xl, yl, zl: real; 2442 localheight: 0.0; 2387 maxht, height, localheight: real; 2443 end; 2388 C: curveptr; 2444 end; 2389 begin 2445 p:= q; 2390 2446 end; 2391 maxlen:= 0.0; 2447 2392 len:= 0.0; 2448 if len > maxlen then 2393 maxht:= 0.0; 2449 begin 2394 localheight:= 0.0; 2450 maxlen:= len; 2395 p:= PX.segmenthead; 2451 maxht:= localheight; 2396 2452 end; 2397 while p <> nil do 2453 2398 begin 2454 xl:= PX.startpoint.x_coord - PX.endpoint.xcoord; 2399 C:= p.layer_head.outer_curve; 2455 yl:= PXstartpoint.y_coord - PX.endpointoy_coord; 2400 case C curve_type of 2456 zl:= PX.startpoint.z_coord - PX.endpoint.zcoord; 2401 line. if C^.curve start > C^.curve end then 2457 2402 height:= C.curve_start 2458 if DEBUG_LD then writeln( 2403 else height = C.curve_end; 2459'* Compute L&D, xl, yl, zl, maxlen =', 2404 arc: if C^.which_half = up then 2460 xl, yl, zl, maxlen ); 2405 begin 2461 2406 if ( C.center.x_coord >= 2462 L:= maxlen * sqrt( xl * xl + yl * yl + zl * zl ); 2407 p lowerbound ) and 2463 2408 ( C.center.xcoord <= 2464 D = 2.0 * maxht;

Jul 8 02:16 1986 Page 47 Jul 8 02:16 1986 Page 48 2577 file_name:='output.data', 2633 xxl:= trunc( xl + 0.5 ) + xorigin; 2578 namesize:= 11; 2634 yyl: trunc( yl + 0.5 ) + yorigin; 2579 2635 xx2:= trunc( x2 + 0.5 ) + xorigin; 2580 { open the disk file for external storage 2636 yy2: trunc( y2 + 0.5 ) + y_origin; 2581 2637 2582 gpr_$openbitmap_file( gpr_$create,filename, 2638 call graphic routine to draw aline from 2583 name size,version,window.window size, 2639 ( xxl, yyl ) to ( xx2, yy2 ) 2584 groups,header,attribs,filebm,created,status); 2640 gpr_$move( xxl, yyl, status ); 2585 if status.all <> status_$ok then writeln( 2641 gpr_$1line( xx2, yy2, status ), 2586' —— > In produce_bitmap, Can not open the 2642 2587'bitmap file for the specified file name'); 2643 2588 2644 end; { drawline 3 2589 { set current bitmap for block transfer 2645 2590 2646 procedure draw_arc( xc, yc, radii: real; 2591 gpr_$setbitmap( filebm, status); 2647 half: up_or_down; 2592 if status.all > status_$ok then writeln( 2648 xl, yl, x2, y2: real ); 2593' —-- > In producebitmap, Can not set', 2649 draw an arc of a circle where its center is 2594'the bitmap file'); 2650 at ( xc, yc ) of radius radii; the arc starts from 2595 2651 ( xl, yl ) to ( x2, y2 ) and half specifies either 2596 [ blok transfer data from current bitmap to disk 2652 the upper part or the lower part in the circle 2597 2653 this arc belongs too This routine should interface 2598 gpr_$pixelblt( init bitmap, window, 2654 the graphic routine to draw arc, circle or lines 2599 window.window_base, status); 2655 if approximation is used. } 2600 if status.all <> status_$ok then writeln( 2656 var 2601' —— > In produce_bitmap, Can not Block',2657 x3, y3, dy: real; 2602'transfer bitmap file'); 2658 first, middle, last: gpr_$position_t; 2603 2659 begin 2604 1 set back the original display bitmap )2660 2605 2661 if DEBUG_DRAW then writeln( 2606 gpr_$setbitmap( init_bitmap, status); 2662'&&&& draw arc,xc,yc,xl,yl,x2,y2', 2607 2663 xc,yc,xl,yl,x2,y2); 2608 2664 2609 end; C producebitmapfile 2665 compute the middle point of the arc ] 2610 2666 x3:= ( xl + x2 ) / 2.0; 2611 2667 dy:= sqrt( radii * radii - 2612 procedure terminate_graph; 2668 ( x3 - xc) * ( x3 - xc) ); 2613 ( terminate the graphic environment 2669 if half = up then y3:= yc + dy 2614 probably make a hard copy. 1 2670 else y3:= yc - dy; 2615 begin 2671 2616 2672 1 convert real coordinates to integer coordinates ] 2617 gpr_$terminate( false, status ) 2673 first.x_coord:= trunc( xl + 0.5 ) + x_origin; 2618 2674 first.y_coord:= trunc( yl + 0.5 ) + y_origin 2619 end; [ terminategraph 2675 middle.x_coord = trunc( x3 + 0.5 ) + xorigin; 2620 2676 middle.y_coord:= trunc( y3 + 0.5 ) + y_origin; 2621 2677 last.x_coord = trunc( x2 + 0.5 ) + xorigin; 2622 procedure draw line( xl, yl, x2, y2: real );2678 last.ycoord:= trunc( y2 + 0.5 ) + y_origin; 2623 1 draw a line from ( xl, yl ) to ( x2, y2 ); 2679 2624 this procedure should interface graphic routine 2680 call graphic routine, check Apollo Domain 2625 1 2681 graphic routine ] 2626 var 2682 gpr_$move( first.x coord, first.y_coord, status ) 2627 xxl, yyl, xx2, yy2: integer; 2683 gpr_$arc_3p( middle, last, status ); 2628 begin 2684 2629 2685 2630 if DEBUGDRAW then writeln( 2686 end; [ draw_arc 3 2631'&&&& draw_line,xl,yl,x2,y2',xl,yl,x2,y2); 2687 2632 2688 procedure draw_ curve( C: curveptr; xl, x2: real );

Jul 8 02:16 1986 Page 45 Jul 8 02:16 1986 Page 46 2465 2521 waittime.highl6:= 0; 2466 writeln(' - > length of the part =, L); 2522 wait_time.low32:= 30 * one_second; 2467 writeln( ~ —— > diameter of the part =', D); 2523 time_$wait( time_$relative, wait_time, status); 2468 2524 2469 end; ( computelengthdiameter ] 2525 end; { holdscreen ] 2470 2526 2471 procedure computeprofile ( PX: principal_axis ); 2527 2472 ( obtain profile of part from PX; intend to support 2528 procedure producebitmap_file; 2473 NC application ] 2529 [ store the bitmap of the whole screen into disk, 2474 const 2530 open the disk file for it 2475 { to adjust the output on the center of the screen ) 2531 if file exists or other errors, write error messages 2476 x_origin = 24; 2532 var 2477 yorigin = 512; 2533 window: gpr $windowt 2478 var 2534 ( window of the display bitmap ] 2479 p: segmentptr; 2535 hi-plane: gpr_$planet; 2480 ratio, xl, yl, zl, x_position: real; 2536 groups: integer; 2481 aa, bb: real; 2537 ( the number of groups in external bitmap 2482 status: status_$t; 2538 status: status_$t; 2483 disp_bm_size: gpr_$offsett; 2539 { returned status 2484 init bitmap: gpr_$bitmap_desc_t; 2540 version: gpr_$versiont; 2485 window: gpr_$window_t; 2541 f version number of bitmap file 2486 2542 header: gpr_$bmf_group_headerarrayt; 2487 procedure init_graph; 2543 1 descriptor of external group bitmap header 2488 [ initialize the Apollo Domain graphic environment. 2544 attribs: gpr_$attributedesct; 2489 Clipping to window is implemtented; be careful!! 2545 f attributes which the bitmap will use ] 2490 begin 2546 filebm: gpr_$bitmapdesc_t; 2491 2547 { descriptor of bitmap ] 2492 disp_bm_size.x_size:= 1024; 2548 created: boolean; 2493 dispbm size.y size:= 1024; 2549 1 specifies whether the bitmap file was created 2494 gpr_$init( gpr_$borrow, 1, disp_bm_size, 0, 2550 filename: name_$pname_t; 2495 initbitmap, status ); 2551 I name of the external file for bitmap 2496 2552 namesize: integer; 2497 set window clipping 2553 length of the file name 2498 window.windowbase.x coord:= 0; 2554 2499 window.window_base.ycoord:= 0; 2555 begin 2500 window.windowsize.x size:= 1024; 2556 2501 window.window_size.y size: 1024; 2557 [ set parameters for the window of operation 2502 2558 2503 [ set'exclusive or' raster operation 1 2559 hiplane:= 0; 2504 gpr_$setrasterop( 0, 6, status ); 2560 groups:= 1; 2505 2561 window.window base.x coord: 0; 2506 gpr_$set_clip_window( window, status ); 2562 window.windowbase.ycoord: 0; 2507 gpr_$set clipping active( true, status ); 2563 window.window_size.x_size: 1000; 2508 gpr_$clear( -2, status ); 2564 window.window size.ysize: 800; 2509 2565 with header(0] do 2510 end; [ initgraph 1 2566 begin 2511 2567 nsects:= hiplane + 1; 2512 procedure hold_screen; 2568 pixelsize:= 1; 2513 1 hold screen for a moment; say 30 seconds ] 2569 allocatedsize:= 1; 2514 const 2570 bytesper_line:= 0; 2515 onesecond = 250000; 2571 bytesper_sect: 0; 2516 var 2572 end; [with 2517 wait time: time_$clock_t; 2573 2518 status: status_$t; 2574 gpr-$allocateattributeblock( attribs, status); 2519 begin 2575 2520 2576 ( getfilename( filename, name_size);

Jul 8 02:16 1986 Page 49 Jul 8 02:16 1986 Page 50 2689 fdraw a curve specified by C, 2745 1 2690 xl ans x2 are horizontal bounds. 12746 begin 2691 var 2747 2692 xc, yc: real; 2748 while q (> nil do 2693 begin 2749 begin 2694 2750 draw -curve( q-.outer curve, xl, x2 ) 2695 case C. curve type of 2751 draw -curve( q.inner-curve, xl, x2 ) 2696 line: draw-line( xl, C.curve-start, 2752 draw~_line( xl, q.outer -curve-.curve -start 2697 x2, C curve-end). 2753 xl, q.inner -curve-.curve -starti 2698 arc:begin 2754 draw-line( x2, q.outercre _rv-n 2699 xc: ratio * C.center~x-coord; 2755 x2, q.inner curve curve-end ) 2700 yc ~C~.center.y coord; 2756 2701 draw-arc( xc, yc, C radius, 2757 draw -curve-buddy( q.inner-curve, xl, x2 ) 2702 C which-half, xl, 2758 draw_curve, buddy( q outer-curve, xl, x2 ) 2703 C curve-start, x2, 2759 draw line( xl, -q.outer-curve curve-start 2704 C.curve-end ) 2760 xl, -q inner -curve-.curve- star ) 2705 end; 2761 draw-line( x2, -q-.outer-curve.curve-end, 2706 end; case1 2762 x2, -q inner curve curve end() 2707 2763 2708 end; draw-curve 2764 q:=,q.next layer; 2709 2765 end; 2710 procedure draw curve -buddy( C: curvejptr; xl, x2: real ) 2766 2711 Idraw the counterpart of a curve specified by C; 2767 end; f draw-layer 2712 xl ans x2 are horizontal bounds. 12768 2713 var 2769 2714 D: curve~ptr; 2770 2715 begin 2771 begin 2716 2772 2717 new( D );2773 1initialize the Apollo graphic environment 2718 D.curve type:= C.curve type; 2774 mnit_graph; 2719 case C curve type of 2775 2720 line:begin 2776 xl1: PX.end~point.x coord -PX.start~point~xcod 2721 D curve-start:=- C -curve-start; 2777 yl:=PX.end~point.y~coord -PX.start-pointycrd 2722 D.curve-end: C curve-end; 2778 zl1: PX.end~point.z coord- PX.start~point~zcod 2723 end; 2779 ratio:= sqrt( xl * xl + yl* yl + zl * zl1) 2724 arc: begin 2780 x-position:= 0.0; 2725 D curve- start: C curve-start;- 2781 2726 D curve- end:=- C curve -end; 2782 p:= PX.seqgnent-head; 2727 D.center.x coord:=C.center.x coord; 2783 2728 D. center~ycoord: -C center~ycoord; 2784 while p 0 nil do 2729 D~ radius,:= C -radius; 2785 begin 2730 if C_.which Thalf = up then 2786 x~position:ratio * p. lower -bound; 2731 D which-half:down 2787 xl:= ratio *p. upper-bound; 2732 else 2788 draw-layer( p.layer head, xposition, x; 2733 D.which-half:up; 2789 2734 end; 2790 p:= p next segment; 2735 end; Icase j2791 end; [ while 2736 2792 2737 draw-curve( 13, xl, x2 );2793 1hold the screen for a movement 2738 2794 hold-screen; 2739 end; I draw-curve buddy 12795 2740 2796 Igenerate the bitmap file for dump 2741 2797 produce bitmapf file; 2742 procedure draw Jlayer( q: layer~ptr; xl, x2: real ) 2798 2743 Ito draw the layer curves by going t~hru q chain; 2799 Iterminate the graphic environment 2744 each counterpart of the curve in q must also be drawn; 2800 terminate_graph;

Jul 8 02:16 1986 Page 51 2801 2802 end; f compute~profile 2803 2804 (MAIN PROCEDURE 2805 (show the applications 2806f P is the principal axis; CSG is the CSG tree to be evaluated; 2807 L is the length of this part and D is the diameter of the part; 2808 other properties for various code assignment can also be derived. 2809 2810 begin 2811 2812 (debugging switches 2813 DEBUGINPUT:= false; 2814 DEBUGEVALUATECSG:= false; 2815 DEBUGBUILDAXIS:= false; 2816 DEBUGCOMPUTET:= false; 2817 DEBUGCURVECURVE:= false; 2818 DEBUGINTERSECT:= false; 2819 DEBUGAXISOP:= false; 2820 DEBUGNORMALIZATION:= false; 2821 DEBUGPARTITION:false; 2822 DEBUGADDLAYER:false; 2823 DEBUGCOMPUTELAYER:= false;2824 DEBUGLAYERDIFFERENCE:= false; 2825 DEBUGMERGEINTERVAL:false; 2826 DEBUGLD:= false; 2827 DEBUGDRAW:=false; 2828 TRACE-AXIS:=false; 2829 2830f get input CSG tree 2831 get input data ( CSG ) 2832 2833 ftransform the CSG tree into axis representation 2834 P:= eval 2835 2836 [ find applications 2837 f support part/shape classification and 2838 process planning 2839 compute length diamneter( P, L, D ) 2840 2841 fsupport NC 2842 compute~pr-ofile( P ) 2843 2844 end.f mainI 2845

4Appendix 2 /nput CSG Data for Bottl/e

Jul 8 15:15 1986 bottle.object Page 1 BOTTLE OBJECT 0# ~~OP L_0* R-O#/MV#/dummy -1 differ 99 98 99 union 44 45 44 union 36 37 36 differ 34 35 34 differ 33 30 33 union 27 29 27 union 17 18 17 union 11 12 11 union 5 6 5 differ 3 4 3 union 1 2 1 primitive 1 1 2 primitive 2 2 4 primitive 4 4 6 differ 9 8 9 differ 7 10 7 primitive 7 7 10 primitive 10 10 8 primitive 8 8 12 differ 13 14 13 union 15 16 15 primitive 15 15 16 primitive 16 16 14 primitive 14 14 18 differ 19 20 19 union 21 22 21 differ 23 24 23 primitive 23 23 24 primitive 24 24 22 differ 25 26 25 primitive 25 25 26 primitive 26 26 20 primitive 20 20 29 differ 31 32 31 primitive 31 31 32 primitive 32 32 30 primitive 30 30 35 primitive 35 35 37 differ 40 42 40 differ 38 39 38 primitive 38 38 39 primitive 39 39 42 differ 41 43 41 primitive 41 41 43 primitive 43 43 45 differ 52 53 52 union 50 51 50 union 48 49

48 differ 46 47 46 prlmltLve 46 46 47 primitive 47 47 49 prlmltLve 49 49 51 primitive 51 51 53 primitive 53 53 98 primitive 98 98 PRIMITIVE.................................... O# GC XL YL ZL X Y Z aX aY aZ 1 torus 0.5000 0.2000 0.0000 0.0000 0.0000 0.2000 0.0000 0.0000 0~0000 2 cylinder 0.4000 0.5000 0.0000 0.0000 0.0000 0.0000 0.0000 0~0000 0.0000 4 cylinder 0.4000 0.3000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 7 cone 1.0000 1.2500 0.0000 0.0000 0.0000 1.0000 180.00 0.0000 0.0000 10 cone 0.8400 1.0500 0.0000 0.0000 0.0000 1.0000 180.00 0.0000 0~0000 8 cylinder 0.4000 0.5000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 15 torus 1.2500 0.3000 0.0000 0.0000 0.0000 1.3000 0.0000 0.0000 0.0000 16 cylinder 0.6000 1.2500 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 14 cylinder 0.6000 1.0500 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 23 cone 1.9672 1.9672 0.0000 0.0000 0.0000 2.3172 180.00 0.0000 0.0000 24 cylinder 1.6000 1.2500 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0~0000 25 torus 1.6844 0.4000 0.0000 0.0000 0.0000 2.6000 0.0000 0.0000 0.0000 26 torus 1.6844 0.2500 0.0000 0.0000 0.0000 2.6000 0.0000 0.0000 0.0000 20 cone 1.9086 1.9086 0.0000 0.0000 0.0000 2.4586 180.00 0.0000 0.0000 31 cone 1.9672 1.9672 0.0000 0.0000 0.0000 2.8828 0.0000 0.0000 0.0000 32 cylinder 1.6000 2.0000 0.0000 0.0000 0.0000 3.3172 0.0000 0.0000 0.0000 30 cone 1.9086 1.9086 0.0000 0.0000 0.0000 2.7414 0.0000 0.0000 0.0000 35 cylinder 0.6000 1.6844 0.0000 0.0000 0.0000 2.3000 0.0000 0.0000 0.0000 38 torus 1.6844 0.4000 0.0000 0.0000 0.0000 3.5400 0.0000 0.0000 0.0000 39 torus 1.6844 0.2600 0.0000 0.0000 0.0000 3.5400 0.0000 0.0000 0.0000 41 cylinder 0.8800 2.2500 0.0000 0.0000 0.0000 3.1000 0.0000 0.0000 0.0000 43 cylinder 0.4412 1.6000 0.0000 0.0000 0.0000 3.3172 0.0000 0.0000 0.0000 46 cone 3.2500 2.5000 0.0000 0.0000 0.0000 5.0000 180.00 0.0000 0.0000 47 cylinder 2.0500 1.6000 0.0000 0.0000 0.0000 1.7000 0.0000 0.0000 0.0000 49 torus 2.1525 0.3000 0.0000 0.0000 0.0000 5.0000 0.0000 0.0000 0.0000 51 cylinder 0.3000 2.1525 0.0000 0.0000 0.0000 5.0000 0.0000 0.0000 0.0000 53 cone 3.2200 2.1525 0.0000 0.0000 0.0000 5.0100 180.00 0.0000 0~0000 98 cube 3.0000 3.0000 5.5000 0.0000 -3.000 ~.0000 0.0000 0.0000 0.0000

Appedix Sha~ded Pictures boy Ray Castkig (a) corresponds to Figure 6.1 (b) corresponds to Figure 6.2 (c) corresponds to Figure 6.3 (d) corresponds to Figure 6.4 (e) corresponds to Figure 7. the'bottle"

(b)

UNIVERSITY OF MICHIGAN