Application of computational geometry in coordinate measurement

10th Jubilee IEEE International Symposium on Applied Computational Intelligence and Informatics • May 21-23, 2015 • Timişoara, Romania Application of Computational Geometry in Coordinate Measurement Gyula Hermann, György Hermann dept. of Applied Mathematics Óbuda University Budapest, Hungary the least-square solution and terminating after a number of iteration. Abstract—Algorithms for calculating the minimum zone deviation for various geometric forms like straightness, flatness, circularity, spherecity and cylindricity have been successfully established by a number. This paper presents simple and robust algorithms based on computational geometric techniques and discusses them from another point of view. The complexity issues are also mentioned. Keywords—straigthness; flatness; cyrcularity; cylindricity; convex hull; Voronoi diagram The linear programming model used by Carr and Ferreira [4] applying the small displacement screw matrix to linearize the non-linear constrains to solve the problem. Kaiser and Krishnan [13] employ “brute force” to find the minimum zone by simply rotating the initial envelope lines. The algorithm stops as the conditions for straightness definition are met. spherecity; Neural networks and interval regression analysis are used by Suen and Chang [32] to compute the minimum zone straightness and flatness. I. INTRODUCTION The majority of manufactured components are bounded simple geometric elements like straight lines, planes, circles, spheres, cylinders and cones. Coordinate measuring machines are used to capture surface points of these geometric elements. The conformity is determined by fitting the ideal geometric feature to these points using suitable algorithms and compared to the nominal value in order to determine whether the deviation is within the required tolerances. The accuracy of the evaluation process depends on the distribution of the measurement points on the inspected surface element, measurement uncertainty in capturing the surface points and the evaluation algorithms. A. Straigthness in the plane Straightness deviation is the minimal distance between two parallel lines (envelop lines) so that each point of the line must be enclosed by them (Fig.1). The Lp-norm solution minimizes the following function: n Lp = ∑ di Fig. 1. Straightness with given measurement points p The minimum zone straightness has the following properties: i =1 Where n is the number of points captured, di is the distance between the i-th datapoint and the ideal feature and 1≤p≤∞. The L2 norm the so-called least square-fit, minimizes the sum of the squares of the residual errors di ; Min∑|di|2. The Chebyshev best-fit (L∞-norm solution) minimizes the maximum deviation from the ideal feature and results in a Min{Max|di|} objective function. The Chebyshev best-fit problem can be transformed into a non-linear constrained optimization problem. The solution requires iterative process whose convergence depends on the initial guess. The envelop lines pass through at least three measuring points ¾ These measuring points must be in a upper-line/ lower-line/upper-line or a lower-line/ upperline/lower-line sequence The algorithm computing the substitute line and the straightness error consists of two main steps: II. STRAIGTHNESS The minimum zone straightness was found by Murthy and Abdin [22] using simplex search. Their algorithm starts with 978-1-4799-9911-8/15/$31.00 ©2015 IEEE ¾ 511 1. Compute the points of the convex hull, resulting in a list of points in traversal order. 2. Determine the substitute line for this convex point set and calculate the width of the convex hull in one step Gy. Hermann and Gy. Hermann • Application of Computational Geometry in Coordinate Measurement that the computed deviation has the same value in every direction. For the determination of the convex hull of a planar point set a number of algorithms have been invented [25] [31]. For our purpose we have selected the modified incremental algorithm because it has an extension into 3D. According to Zhang et al. [33] that finding the spatial straightness error based on the minimum zone condition cannot be found using a linearized model. They returned to solving the nonlinear model. The determination of the spatial straightness is equivalent to finding the radius of the minimal circumscribed cylinder defined by the coordinates of the measured points. It can be proved that the axis of minimum 3D straightness is parallel to an edge of the convex hull. This leads to the following algorithm: Fig. 2. Determination of the tangency lines from p to Q Compute the vertices and edges of the convex hull set the initial straightness value equal to a large number for i = 1 to ne do for j = 1 to nv do {Compute the projection v’j of all vertices vj onto the plane normal to ei Compute the diameter di of the minimum circumscribed circle of the projected measured point} if (di < straightness) then straightness = di The algorithm starts with the triangle p0,p1,p2. If the point p to be added is inside the triangle, than it is discharged. If it is outside than it is added to the convex hull and those point which can be „seen” from p and are not on the perifery are removed. Therefore two tangent lines from p to Q are determined using the LeftOn test (on which side of the line the point is). If the points to be processed are sorted according their x coordinates the algorithm determines the convex hull in O(nlogn) time. III. FLATNESS Flatness is defined as the orthogonal distance between two parallel planes separated by the specified tolerance within which all surface points must lie (Fig. 2.) After the determination of the convex hull the computation of the substitute line and the width of the convex hull follows. The envelop lines bounding the minimum zone are determined by two points defining one envelope line colinear with one of the edges of the convex hull, the third one is a point with the maximum perpendicular distance to this line. The algorithm does the job in O(n2) time, where n is the number of extreme points on the convex hull: d=0 for i = 1 to n {c=0, j=i+2 while j<2n dist(ei,pj)>d {c= dist(ei,pj), j=j+1} if d<c {d=c, l=i, m=j-l} } Fig. 4. Definition of the flatness with given datapoints The substituting line is defined by the points (pl+pm)/2 and (pl+1+pm)/2. The envelop lines are given respectively by pl, pl+1 and pl,pm-pl+pl+1.. Computational geometric techniques were used by Samuel and Shunmugam [27] used to solve the minimum zone problem for the straightness and flatness problem. But they have not taken the 2-2 case into account. On the contrary the paper by Lee [18] considers the 2-2 model and gives a convex hull based solution. The minimum zone flatness solution has the following properties: B. Spatial straigthness The spatial straightness deviation is defined as the diameter of the minimum cylinder enclosing all measured points on surface. The task is to find the axis of the cylinder and the minimum cycle around the projected measured point on a plane perpendicular to this axis. Fig. 3. Definition of the spatial straightness Huang, Fan and Wu [10] describe a method for spatial straightness error evaluation from the composition of vertical and horizontal planer straightness errors. There is no guarantee 512 ¾ Determining the substitute plane at least four contacting measuring points are used ¾ If two of the four contact points define a line on the upper plane and two on lower plane we speak about a 2-2 model. The line on the upper plane will intersect the line on lower plane when they are projected onto the same plane ¾ If the four contact points divided as three points contacting one plane and one point contacting the other one, the single point must be inside the projected triangle on the other plane (3-1 model). 10th Jubilee IEEE International Symposium on Applied Computational Intelligence and Informatics • May 21-23, 2015 • Timişoara, Romania It is obvious that the measuring points defining the minimum zones are vertex points on the convex hull. Flatness deviation is the distance between the two lines (2-2 model) or the distance between the point and the planes (3-1 model). The overall structure of the algorithm is the same as in two dimensions. If the point to be inserted is inside the convex hull, then it is removed, if not the cone with apex pi and tangent to the convex hull is added. It is clear that faces visible from pi are removed. Fig. 5. Convex hull before and after adding a point Fig. 6. Minimum circumscribed (a), maximum inscribed (b), minimum zone (c) and least-square fit (d) circles The incremental algorithm for the determination of the convex hull is given below: Wang et. al [35] formulated and solved the minimum zone circularity error as a nonlinear optimization problem. Initialize H3 to tetrahedron (p0,p1,p2,p3) for i = 4 to n-1 do for each face f of Hi-1 do Mark f if visible if no faces are visible then dicard pi (it is inside Hi-1) else for each border edge e of Hi-1 do Construct cone face determined by e and pi for each visible face f do Delete f Update Hi Huang [11] solves the minimum zone circularity error using Voronoi diagrams. The plane was divided into a finite number of so-called max regions. Not all intersection point of the farthest and the nearest regions are contenders for the optimum solution. These points can be divided into two categories called X and Y types and only X types participate in the optimal solution. Marking the visible faces and constructing the cones are imbedded in a loop that iterates n times, the complexity of the algorithm is O(n2). Li Xiuming and Shi Zhaoyao [20] followed a different approach. The roundness error is obtained as the distance between two parallel lines. The lines for the maximum inscribed circle and the minimum circumscribed circle can be found based on convex hull and polar coordinates and determined by four critical measured points. The distance between these lines is the roundness error. Having computed the points of the convex hull the next step is the determination of the substitute planes and the width of the convex hull. The boundary planes of the minimum zone are determined by the minimum distance between the face spanned by three points and one point with the maximum perpendicular distance to this face, or in the other case by the distance between two lines each in in a face on the two sides of substitute plane. The strategy described Jiing-Yih Lai and Ing-Hong Che [16] also employs a non-linear transformation to convert a circle into a line. A straightness evaluation was used to obtain appropriate control points and minimum zone deviation. d=0 for i = 1 to n Compute the equation of the face fi from it’s cornerpoints for j = 1 to n Determine the distance of the point pj from the face if dist(pj,fi) > d then d = dist(pj,fi) for i =1 to m for j = 1 to m if dist (ei, ej) > d then d = dist (ei, ej) > d Lei et al. [19] in their paper present a polar coordinate transformation for circularity evaluation. The algorithm allocates a circular region around the center of the least square circle following certain rules and calculates the polar radius for all measuring points by translating the polar coordinate system to each point in the region and obtains minimum circumscribed, maximum inscribed and minimum zone center from comparing each polar radius relative to each polar coordinate system. The complexity of the algorithm is O(m2) time, where m is the number of edges on the convex hull. Samuel and Shunmugam [28] use the convex hull and a heuristic algorithm to find the inner hull. Equi-Distant Voronoi and Equi-Angular diagrams are employed to determine the feature. IV. CIRCULARITY The definition for roundness error given by ANSI [1] and ISO [13] as the minimal radial separation of two concentric circles, the circumscribed and the inscribed circle so, that the ring contains all measured points. The nearest-point Voronoi diagram (denoted by NV) partition the plane by assigning every point in the plane to its 513 Gy. Hermann and Gy. Hermann • Application of Computational Geometry in Coordinate Measurement This leads to a runtime of O(n2) but an average runtime of O(n) can be expected for well distributed sets of sites. The farthest-point Voronoi diagram (denoted by FV) is the set of all points that are farther from the given datapoint (Pi) than from any other measuring point. nearest measuring point (Pi). The Voronoi diagram is the settheoretic intersection of the halfplanes generated by the bisectors of every pointpair: NV ( p i ) = ∩ H ( Pi , P j ) i≠ j The minimum zone circles should pass through at least four points. Three different cases can be distinguished: The regions (point sets) attached in this way to every measuring point form convex regions on the plane. Green and Sibson [8] invented an incremental algorithm for constructing the nearest-point Voronoi diagram. It starts with setting up a simple Voronoi diagram for three data-points, and modifying the diagram by adding measuring points one by one. ¾ Three points determine the circumscribed circle and the fourth point the radius of the inscribed circle. ¾ Three point lie on the inscribed circle and the fourth point determines the radius of the circumscribed circle. ¾ Two point lie both on the circumscribed circle and two points on the inscribed circle. The first step in all cases is the generation of the convex hull of the point set. Hereby the number of points to be considered in the next steps is reduced. Determine the farthestpoint and nearest-point Voronoi diagrams and superimpose them. Fig. 7. Near-point and farthest-point Voronoi diagram Having nearest-point Voronoi diagram NVp-1, we would like to add a new measuring point pk. First, find the measuring point, say pi, whose region contains pk, and draw the perpendicular bisector between pk and pi. The bisector crosses the edge of NV(pi) at two points, point x1 and point x2. Point pk is to the left of the directed line segment x1x2. The line segment x1x2 divides the nearest-point Voronoi polygon NV(pi) into two pieces, the one on the right belonging to the Voronoi region of pk. Thus, we get a nearest-point Voronoi on one edge of the boundary of the Voronoi region of pi. The new edge crosses the boundary at x2, entering the adjacent Voronoi polygon, say NV(pj). Therefore, next draw the edge between pk and pj, and find the point, x3, at which the bisector crosses the boundary of NV(pj). Similarly, find the sequence of segments of perpendicular bisectors of pk and the neighboring sites until we reach the starting point x1. This sequence forms a clockwise boundary of the new Voronoi region. Finally the edge segments of NVp-1 the substructure inside the new Voronoi region are deleted. Fig. 8. Superimposed near-point and farthest point Voronoi diagrams The intersection points of the two diagrams are the candidates for the center of the minimum zone circles. The complexity of finding the minimum zone circularity based on the set of Voronoi diagrams is O(n2). V. SPHERECITY Minimum zone sphericity error is defined as the minimal radial distance between two concentric spheres the circumscribed and the inscribed spheres to the measuring points, so that all points are included in the shell. The farthest-point Voronoi diagram (denoted by FV) is the set of all points that are farther from the given datapoint (Pi) than from any other measuring point. In Huang’s [11] approach the whole space into several regions using Voronoi diagram. The minimum zone spherecity becomes feasible on the vertices of these regions. Calculate the convex hull of the measured points for all p in S add v(p) to K; if n>2 then repeat find p maximizing (radius(before(p),p,next(p)),angle(before(p),p,next(p)); q= before(p); c= center(q,p,next(p)); add c to K; add (c,v(p)) and (c,v(q)) to E; v(q)= c; next(q)=next(p); before(next(q))=q; n=n-1; until n=2; add(v(q),v(next(q)) to E else if n=2 then {S={p1,p2}} add(v(p1),v(p2)) to E Kanada [15] calculated the minimum zone solution using downhill simplex search, a sequential gradient search. The simplex used is a tetrahedron. The object value at a vertex is the difference between the farthest and nearest distances from the measured point. In each iteration step, the vertex with the maximal value will be rejected and replaced by a new measured point in the opposite side of the original tetrahedron. If the objective value of the new point is still the greatest the search is continued in alternate direction, or stopped. The algorithm does not guarantee a global minimum. 514 10th Jubilee IEEE International Symposium on Applied Computational Intelligence and Informatics • May 21-23, 2015 • Timişoara, Romania The minimum zone spheres should pass through at least five point of the dataset. There are four different cases: ¾Four points determine the circumscribed sphere and the fifth point the radius of the inscribed sphere. ¾Three point lie on the circumscribed sphere and two points determines of the inscribed sphere. ¾Two point lie on the circumscribed sphere and the three points on the inscribed sphere. ¾Four points determine the inscribed sphere and one point the radius of the circumscribed sphere. The candidates for being the center of the minimum zone spheres are: ¾The common centers of the farthest-point and nearestpoint Voronoi diagrams Fig. 9. Measured point and the assessment sphere ¾The common centers of the farthest-point and nearestpoint Voronoi diagrams, lying on the edge of a tetrahedron in the other diagram Two papers by Samuel and Shunmugam [29],[30] deals with the spherecity problem. Their solution starts with the construction of the 3D hulls. The outer hull is determined by computational geometric technique, followed by a heuristic approach to find the inner hull. Assesment limacoids are constructed using the concept of 3D equi-angular lines, 3D farthest or nearest equi-angular diagrams. ¾The intersection points of edges of farthest-point and nearest-point Voronoi diagrams In the first case four points determine on the spheres and the fifth the radius of the other one. In the second case four points lie on either the circumscribed or the inscribed sphere depending on whether the center is on the farthest-point or on the nearest-point Voronoi diagram. In the third version the other sphere is defined by three points and the inner one also by three points. The radial distances between the concentric spheres for all candidate centers are calculated and the smallest one is selected. Rossi et. al [23] attack the problem by finding a closed form for the minimum neighborhood of the centroid of the measured sphere. The minimum zone evaluation outlined for circularity can be extended for determining the spherecity error. Let S be a set of point in the 3-dimensional Euclidian space. The set of points which are closer to p, than to any other point in S, form Voronoi cell of a point p in S. The union of all Voronoi cells is the Voronoi diagram. The cells are convex polyhedrons. VI. CYLINDRICITY Delaunay tetrahedralization is the dual to the Voronoi diagram and partitions the space into tetrahedral, where the vertices of the triangles are points p in S. A Voronoi vertex is in the center of the sphere circumscribing its dual tetrahedron. Minimum zone cylindricity is bounded by two coaxial cylinders with a minimal radial separation within which all surface points must lie. The straightness of the median line (the axis of the cylinder) is the diameter of the cylindrical zone containing the median line. There is an incremental algorithm [13] for generating the Delaunay triangulation in O(n2) complexity. This algorithm inserts points one at a time keeping a valid Delaunay triangulation and the tetrahedralization is modified according to the Delaunay criterion in each step. The modification of the diagram is local. The algorithm starts with establishing a “big„ tetrahedron containing all points of S, and hereby the convex hull of S. This has the disadvantage that more tetrahedra are constructed than needed. They are marked easily as they are the ones containing the one of the point of forming the “big„ tetrahedron. By introducing the “big„ tetrahedron a number special cases are avoided. Fig. 10. Cylindricity and straightness error Carr and Ferreira [5] developed a methodology which can be applied for both to the cylindricity and the straightness of the median line. The determination of the minimum zone cylindricity is a non-linear optimization problem which they solved by a sequence of linear programs that converge to a local optimum. The new points are inserted according the algorithm invented by Joe [14]. Next the Delaunay triangulation converted into the 3D near-point Voronoi diagram connecting the center point of circumscribed spheres of the neighboring tetrahedra. The farthest point Voronoi diagram can be derived on a rather similar way. Lai and Chen [16] employ a non-linear transformation to convert a cylinder into a plane and use flatness evaluation scheme to obtain appropriate control points and minimum zone deviation. 515 Gy. Hermann and Gy. Hermann • Application of Computational Geometry in Coordinate Measurement [11] J. Huang “An exact solution for roundness evaluation problem” Computer Aided Design 31(1999) 845-853 [12] ISO/R 1001-1983 Technical drawings – Geometrical Tolerancing – Guidelines, ISO, Geneve [13] B. Joe „Construction of three-dimensional DeLaunay triangulation using local transformations” Computer Aided Geometric Design, 8(1991) 123142 [14] M.J. Kaiser, K.K. Krishnan „Geometry of the minimum zone flatness functional: Planar and spatial case” Precision Engineering 22 (1998)174183 [15] T. Kanada, „Evaluation of spherical form errors – Computation of sphericity by means of minimum zone method and some examinations with using simulated data” Precision Engineering, 17 (1995) 281 – 289 [16] J.-Y. Lai, I.-H. Chen „Minimum zone evaluation of circles and cylinders” International Journal of Machine Tools & Manufacture 36(1996) 435-451 [17] Y.-Z. Lao, H.-W- Leong, F. P. Preparata, G. Sing „Accurate cylindricity evaluation with axis ” Precision Engineering 27 (2003) 429-437 [18] M.K. Lee „A new convex-hull based approach to evaluating flatness tolerance” Computer Aided design 29 (1997) 861-868 [19] Lei Xianqing, Zhang Chunyang, Xue Yujun, Li Jishun „Roundness error evaluation algorithm based on polar coordinate transform” Measurement 44 (2011) 345–350 [20] Li Xiuming, Shi Zhaoyao „Development and application of convex hull in the assessment of roundness error” Int. Journal of Machine Tools & Manufacture (2007) [21] M. Morel, H. Haijtjema “Task Specific Uncertainty Estimation for Roundness Measurement” Proc. of the 3rd euspen International Conference, Eindhoven, The Netherlands, (2002) 513-516 [22] T.S.R. Murty, S.Z. Abdin „Minimum zone evaluation of surfaces” Int. Journal of Machine Tool Design and Research, 20 (1980) 123-136 [23] A. Rossi, S. Chiodi, M. Lanzetta “Minimum centroid neighborhood for minimum zone spherecity” Precision Engineering 38(2014) 337-347 [24] J. O’Rourke „Computational geometry in C” Cambridge University Press 1998 [25] F.P. Preparata, M. Shamos „Computational Geometry” New York Springer Verlag, 1985 [26] A. Rossi, S. Chiodi, M. Lanzetta „ Minimum centroid neighborhood for minimum zone sphericity” Precision Engineering 38 (2014) 337-347 [27] G.L. Samuel, M.S. Shunmugam „Evaluation of straightness and flatness error using computational geometries techniques” Computer Aided design 31 (1999) 829-843 [28] G.L. Samuel, M.S. Shunmugam „Evaluation of circularity from coordinate and form data using computational geometric techniques” Precision Engineering 24 (2000) 251-263 [29] G.L. Samuel, M.S. Shunmugam „Evaluation of sphericity error from coordinate measurement data using computational geometric techniques” Computer methods in applied mechanics and engineering 190 (2001) 6765-6781 [30] G.L. Samuel, M.S. Shunmugam „Evaluation of sphericity error from form data using computational geometric techniques” International Journal of Machine Tools & Manufacture 42 (2002) 415-416 [31] Sariel Har-Peled „An output sensitive algorithm for discrete convex hulls” Computational Geometry 10 (1998) 125-138 [32] D.S. Suen, C.N. Chang „Application of neural network interval regression method for minimum zone straightness and flatness” Precision Engineering 20 (1997) 196-207 [33] Qing Zhang, K.C. Fan, Zhu Li „Evaluation method for spatial straightness errors based on minimum zone condition” Precision Engineering 23 (1999) 264-272 [34] N. Venkaiah, M.S. Shunmugam „Evaluation of form data using computational geometric techniques- Part II: Cylindricity error” Int. Journal of Machine Tools & Manufacture 47 (2007) 1237-1245 [35] M. Wang, S. Hossein Cheraghi, Abu S. M. Masud „Circularity error evaluation theory and algorithm” Precision Engineering 23 (1999) 164176 The method invented by Lao et. al [17] avoids the direct construction of the zone cylinder, but approximate it with a guaranteed accuracy through an iterative process based on the linearization of the problem. Venkaiah and Shunmugam [34] describe the evaluation of cylindricity error using a form tester. Due to misalignment between the component and the instrument axis, a limacon – cylinder has to be used for cylindricity evaluation. The paper established the limacon cylinder using computational geometric technique. Fig.11. Recommended distribution of the measuring point on a cylinder: bird cage, roundness profiles, generatrix, scattered points The evaluation of cylindricty error is influenced by the measuring strategy, the distribution of the measuring points on the surface. In case of regular distribution: bird cage, roundness profile and generatrix an obvious procedure is to estimate the axis by for example finding the minimal enclosing box. One of the principal axis of the box is a good estimation of the cylinder axis we are searching for. Given the axis the minimal and maximal diameters can be calculated using circularity data. However this approach does not give information about the straightness of the axis. For finding the minimum zone cylinders currently is no computational geometry based procedure available. ACKNOWLEDGMENT The authors gratefully acknowledge the stimulating discussion and encouragement of prof. Imre Rudas. REFERENCES [1] ANSI/ASME Y14.5M Dimensioning and Tolerancing. New York, American Society of Mechanical Engineers, 1994 [2] M. de Berg, O. Cheong, M. van Kreveld, M. Overmars “Computational Geometry, Algorithms and applications” Springer Verlag 2008 [3] R. Calvo1, E. Gómez, and R. Domingo „Vectorial Method of Minimum Zone Tolerance for Flatness, Straightness, and their Uncertainty Estimation” Int. Journal of Precision Engineering and Manufacturing Vol. 15, (2014) No. 1, 31-44 [4] K. Carr, P. Ferreira „Verification of form tolerances Part I: Basic issues, flatness, and straightness” Precision Engineering 17 (1995) 144-156 [5] K. Carr, P. Ferreira „Verification of form tolerances Part II: Cylindricity and straightness of a median line” Precision Engineering 17 (1995) 131 [6] O. Devillers, F. P. Preparata “Culling a set of points for roundness or cylindricity evaluation” INRIA Rapport de recherché no 4159 (2001) [7] O. Devillers, F. Preparata. “Evaluating the cylindricity of a nominally cylindrical point set” ACM-SIAM Sympos. Discrete Algorithms, Jan 2000, San Francisco [8] P. J. Green and R. R. Sibson. Computing Dirichlet tessellations in the plane. Comput. J., 21:168{173, 1978. [9] L.J. Gubia, J. Stolfi „Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams”, ACM Transactions on Graphics 4, 74-123 [10] J. Huang, K.C. Fan, J.H. Wu „A new minimum zone method for evaluating straightness errors” Precision Engineering 23(1999) 2-8 516