Flatness,Cylindricity and Sphericity Assessment Based on the Seven Classes of Symmetry of the Surfaces

Abstract

Dimensional inspection of a manufactured surface by means of a coordinate measuring machine (CMM) produces a set of Cartesian coordinates. The coordinates are processed to yield the geometric tolerance of the surface. This paper presents a new approach to the evaluation of flatness, cylindricity and sphericity tolerance based on surface invariance with regard to the rigid motions. The proposed algorithm transforms, through homogeneous transformation matrices, the coordinates measured to best fit the reference element of the surface class from which the actual measurements were sampled. The transformation matrix is simplified taking into account the invariance of the sum of the squared normal distances of the measured points from the nominal surface as regards some rigid motions. This invariance is a consequence of the invariance as regards some displacements of the nominal surface from which the data points were sampled. In this way, the number of parameters to be optimised is reduced in comparison with the six parameters characterizing the general homogeneous transform matrix. The methodology was computer implemented and numerical simulations were performed for planes, cylinders, and spheres in order to validate the effectiveness of the approach. The results indicate that the proposed algorithm provides accurate and quick assessments.

1. Introduction

As industry strives towards the increase of product quality, the effects of machining errors on a workpiece need to be notably reduced. This is critical for the manufacturing processes in order to achieve the needed workpiece accuracy. After the machining, it is necessary to verify the possible occurrence of manufacturing errors, in terms of dimensional or geometrical variations of the designed shape. The traditional inspection method of a manufactured part implies several steps: from the inspection plan design with the selection of number and position of sample points to the measurement process, to the identification of an ideal geometric element with known analytical expression that best fits the measured data. Several fitting techniques can be used to determine the substitute geometry, such as the Least Square Method and the Minimum Zone. The final step is the comparison of the evaluated deviations with the imposed tolerances, thus accepting or rejecting the manufactured part.

Nowadays, GPS standards are evolving towards innovative concepts about surfaces classification [1–3]. The modern vision implies a surface classification based on continuous symmetry of geometric objects. This classification is founded on twelve connected Lie subgroups of rigid motion that led to a compact classification of surfaces in seven classes based on their symmetry and invariance with regard to the rigid motions [4–7]. This paper is the first attempt at integrating these principles into a new algorithm for form tolerance assessment.

According to the standards [8, 9] of geometric tolerance, a tolerance applied to a feature of a mechanical part creates a 3D portion of space (the tolerance zone), where the feature may be found, and whose geometry, size, position, and orientation are specified through the type of geometric tolerance assigned, its value, and the datum system chosen, respectively.

The conformance of the manufactured parts to the design specifications is evaluated during the inspection phase. In dimensional inspection, the CMM extracts from the real integral feature a finite set of sampled points and generate the so-called extracted integral feature. The extracted integral feature is elaborated to calculate the associated integral feature (even called substitute feature [10]), that is, the nominal feature that best fits the extracted integral feature, and the corresponding associated derived feature. The deviations of the measured points from the associated integral feature are finally evaluated to yield the form tolerance. For example, to evaluate the straightness error, an ideal straight line has to be established by minimizing, through the choice of a suitable objective function, the maximum deviation between it and the extracted integral feature concerned. The form error is, then, computed as the maximum peak-to-valley distance of the measured points from the associated integral feature. The position error of the actual feature is derived from the estimated parameters of the substitute feature.

In previous works, different objective functions were proposed: sum of the squares of deviations [11], sum of the squares of normal deviations [11, 12], sum of absolute deviations [13, 14], and average deviation [13]. The minimization was carried out based on either direct or random search techniques, Monte Carlo simulation, simplex search, or spiral search techniques. In recent years, a new class of algorithms for the establishment of form tolerance have received much attention [15–21]. This new class, instead of a substitute feature, directly searches the minimum zone of the feature under inspection. For example, in flatness assessment, two parallel planes with minimum distance including the whole cloud of measured points are searched. This method best conforms the ISO standards that define the form tolerance by means of a tolerance zone within which the feature is to be contained. However, this kind of algorithms is very sensitive to asperities, while the methods minimizing an objective function are much less sensitive to asperities and they are mathematically robust [22].

This paper presents a new theoretical approach for the assessment of form tolerances by sampling a set of measured points. The proposed approach is based on the new concepts under development by the ISO/TC 213: in particular, on the classification of surfaces into seven classes of symmetry. Instead of searching for the best fit substitute feature, the proposed method transforms, through homogeneous transformation matrices, the coordinates measured in order to best fit the reference element of the class of the surface from which the actual measurements were sampled. The best fitting transformation matrix is found by means of the minimization of the sum of the squared normal distances of the points measured from the reference element of the class to which the nominal integral feature under inspection belongs. The transformation matrix is simplified taking into account the invariance of the sum of the squared normal distances of the measured points from the nominal surface. This invariance is a consequence of the invariance as regards some displacements of the nominal surface, and subsequently of its datum. In this way, the number of parameters to be optimized is reduced in comparison with the six parameters characterizing the general homogeneous transformation matrix in 3D. This theoretical approach was implemented into numerical algorithms for the assessment of the flatness, cylindricity, and sphericity tolerances. The new algorithms were tested on both simulated data sets and data reported in literature. In the following, the problem of form tolerance assessment is introduced (Section 2). Then, the proposed algorithm is mathematically described (Section 3) and applied to flatness (Section 3.1), cylindricity (Section 3.2), and sphericity (Section 3.3) assessment. Finally, the performances of the proposed algorithm are discussed by means of different sets of data (Section 4).

2. Traditional Approach to Form Tolerance Assessment

The algorithms for the evaluation of the form tolerance found in the literature search for a substitute feature, that is, a feature of the same geometric nature of the surface from which the measurements were sampled. The substitute feature minimizes the deviation between itself and the extracted integral feature.

In general form, an ideal surface can be express as a set of points:

W = {n = {[n_{x}, n_{y}, n_{z}, 1]}^{T} : n_{x}, n_{y}, n_{z} ∊ R, w (n, s) = 0},

(1)

where w is a C¹ function that represents the equation of the nominal surface specified in implicit form.

The function w(n,s) contains internally the set of parameters s that characterize position, orientation, and size of the feature. For example, w for a sphere is

{(x - x_{0})}^{2} + {(y - y_{0})}^{2} + {(z - z_{0})}^{2} = R^{2},

(2)

where the feature parameters are (x₀, y₀, z₀) that individualize the position of the sphere centre and R that specifies its size. w is a C¹ class function when the nominal feature is almost regular (i.e., the normal to the surface at any surface points exists). These types of surfaces are called CMM measurable.

The extracted integral feature is also represented by a set of points:

W_{m} = {r_{i} = {[r_{x i}, r_{y i}, r_{z i}, 1]}^{T} : r_{x i}, r_{y i}, r_{z i} ∊ R},

(3)

where r_i is the ith actual measurement of the m measured points.

The vectorial relationship between the measured points and the substitute feature is

r_{i} = n_{i} + ε_{i},

(4)

where $n_{i} = {[n_{x i}, n_{y i}, n_{z i}, 1]}^{T}$ is the point on the nominal surface such that r_i - n_i is coincident with the surface normal direction. ε_i represents the feature error evaluated along the surface normal direction at the point n_i. It is the sum of two elementary errors: measurement and manufacturing errors.

At this step the methods proposed in the literature proceed to the search of the set of substitute feature parameters, s, by minimizing the deviation between the substitute feature and the extracted integral feature through the minimization of an objective function, that usually has the following expression [16]:

L_{2} (s) = \sum_{i = 1}^{m} {[ε_{i} (s)]}^{2} .

(5)

This function is easily recognized as the optimization function of the normal least squares. Statistically, the least squares fit is the maximum likelihood estimation assuming normally distributed errors [23]. An alternative method proposed by Goch [24] is based on a general L_p norm. Using an L₂ norm for the optimization objective function means to implement the least squares fitting; while an L_∞ norm (Chebyshev norm) means to implement a min-max fit. Goch showed that as the value of p in L_p increases, the fitting algorithm asymptotically approaches a min-max fit.

The minimization of the chosen objective function with respect to the s parameters is in general a complex nonlinear problem, so the objective function usually is linearized to find the solution directly.

3. New Approach to Form Tolerance Assessment

Datum concept is the guiding principle of the proposed approach. Datums are mostly used to position a feature in relation to another one [3, 4]. In Table 1 for each of the seven classes of symmetry, the minimum set of points, lines, and planes that belong to the same automorphism group as the geometric objects in that class is reported. Those points, lines, and planes may substitute each surface within the same class of symmetry in problems of relative positioning.

Table 1:

Seven classes of symmetry. dim(Aut0(S)): dimension of the automorphism group; I: identity motion, [3].

Class of symmetry	Aut₀(S)	Reference element (datum)	dim (Aut₀(S))

Spherical	Rot(3)	point	3
Cylindrical	Tr(1) × Rot(1)	straight line	2
Planar	Tr(2) × Rot(1)	plane	3
Helical	Tr(1) × Rot(1) with pitch μ	helix	1
Revolute	Rot(1)	(point, straight line)	1
Prismatic	Tr(1)	(straight line, plane)	1
Complex	I	(point, straight line, plane)	0

Positioning a cloud of measured points as regards to its nominal integral feature is the same as positioning the cloud of measured points as regards to a reference element of the class of symmetry to which its nominal integral feature belongs. In this way, the deviation of the measured points from a reference element can be defined as an alternative to the deviation defined by (4). The deviation is the distance of the measured points from the reference element of the nominal surface to which the cloud of points has to be fitted, while in (4) the deviation is the error between the measured points and the points on the nominal surface that satisfy the orthogonality condition. This deviation will be characterized for the planar, cylindrical, and spherical surfaces. Thus, instead of finding the parameters of the substitute feature by minimizing the deviation between the substitute feature and the extracted integral feature, the proposed algorithm transforms through an opportune homogeneous matrix the measurement Cartesian coordinates. The best transformation parameters are identified minimizing the distance between the cloud of measured points and a geometric element having the same geometric nature of the reference element of the class of the surface from which the actual measured points were sampled.

By applying a homogeneous coordinate transformation T(t), the extracted integral feature is roto-translated so that the measured points can minimize their distance from the reference element of the nominal surface. This rigid motion of the measured points is carried out in the coordinate frame of the datum element of the class to which the nominal integral feature under inspection belongs. In this paper, the Roll-Pich-Yaw angles are used to describe the spatial rotations and translations. Therefore, the transformation matrix can be described as [6, 7]

T (t) = [\begin{bmatrix} C ϕ C θ & C ϕ S θ S ψ - S ϕ C ψ & C ϕ S θ C ψ + S ϕ S ψ & t_{x} \\ S ϕ C θ & S ϕ S θ S ψ + C ϕ C ψ & S ϕ S θ C ψ γ - C ϕ S ψ & t_{y} \\ - S θ & C θ S ψ & C θ C ψ & t_{z} \\ 0 & 0 & 0 & 1 \end{bmatrix}],

(6)

where C and S mean Cos and Sin, respectively. The transformation matrix is generally characterized by six degrees of freedom that are the six transformation parameters $t = [ψ, θ, ϕ, t_{x}, t_{y}, t_{z}]$ , that is, the three rotation angles ψ, θ, and ϕ around the x-, y-, and z-axis of the datum reference system and the three translations t_x, t_y, and t_z along the same axes.

Considering the new definition of the deviation, stated above as opposed to (4), the objective function to be minimized becomes

Q (t) = \sum_{i = 1}^{m} {[T (t) \cdot r_{i} - n_{i}]}^{2},

(7)

where n_i represents now the point nearest to $T (t) \cdot r {}_{i}$ on the reference element.

The proposed approach presents two advantages. First, the objective function to be minimized has the same transformation parameters, t, for all the nominal integral features, while the classical least squares approach uses parameters, which are specific and different for any different geometric features, [12]. Furthermore, the number of parameters to be minimized, that is the degrees of freedom of the transformation matrix, which in the general case are six, can be reduced considering the invariance properties of the reference element as regards some rigid motions. This invariance produces the invariance of the objective function (7) according to the same rigid motions that leave invariant the nominal feature under inspection belongs.

The proposed method has an iterative implementation, where ε is the chosen precision of the algorithm. For each iteration j and for each measured point, $r_{i}^{j - 1}$ , the nearest point on the reference feature, $n_{i}^{j}$ , is computed. Then, the objective function is minimized searching for the best parameter set t_j:

Q^{j} = \sum_{i = 1}^{m} {(T (t^{j}) r_{i}^{j - 1} - n_{i}^{j})}^{2},

(8)

where $T (t^{j})$ is the transformation matrix computed for the jth set of parameter $t^{j}$ .

The search of the minimum of the function (8) is a nonlinear optimisation problem. It was solved using the numerical well-established method of Levenberg-Marquardt [25] that is very effective when the function to be minimized is a sum of squares.

Finally, the measured points are transformed by $T (t^{j})$ :

r_{i}^{j} = T (t^{j}) \cdot r_{i}^{j - 1},

(9)

and the convergence criterion is checked, $| Q^{j} - Q^{j - 1} | < ε$ , until the convergence is achieved. The best fit transformation, $\tilde{T}$ , is computed by multiplying the transformation matrices found at each step: $\tilde{T} = \prod_{j = 1}^{# of steps} T^{j}$ . The vector of the optimal parameters, extracted from $\tilde{T}$ , is henceforward called $\tilde{t} = [\tilde{ψ}, \tilde{θ}, \tilde{ϕ}, {\tilde{t}}_{x}, {\tilde{t}}_{y}, {\tilde{t}}_{z}]$ .

The proposed approach is applied to assess flatness, cylindricity, and sphericity tolerance in the following paragraph.

3.1. Flatness Assessment

The planar class of symmetry includes the surfaces which are invariant under planar displacements that are the two translations parallel to the plane itself and the rotation around the perpendicular to the plane. The reference element of this class will necessarily be a plane. There are, however, infinite planes, all parallel to each other, that can be chosen as a reference element of a surface of planar class. Then, in flatness assessment the measured points are transformed step by step by the algorithm to best fit a plane. Any plane in the Euclidean space could be chosen, but for the sake of simplicity the x-y plane of the CMM is chosen. This completely characterizes the nearest point of r_i and the structure of the objective function.

It is obvious that the nearest point to r_i on the nominal plane is the point that has the first two coordinates equal to the first two coordinates of r_i and the third equal to zero: $n_{i} = {[r_{x i}^{}, r_{y i}^{}, 0,1]}^{T}$ . Equation (7) can be further simplified considering that the objective function in this case is invariant for rotation and translations of the cloud of measured points around z-axis and along x- and y-axis of the plane reference frame, respectively. This invariance is a consequence of the invariance respect to the same displacements of the nominal surface, and subsequently of its reference element. In this way, the number of parameters to be optimised is reduced from 6 to 3, because it is possible to drop out ϕ, t_x, and t_y from (7), imposing ϕ, t_x, and t_y equal to zero. Then, the generic objective function is replaced for the planar class of symmetry by the following one:

Q (t) = \sum_{i = 1}^{m} {([\begin{bmatrix} C θ & S θ S ψ & S θ C ψ & 0 \\ 0 & C ψ & - S ψ & 0 \\ - S θ & C θ S ψ & C θ C ψ & t_{z} \\ 0 & 0 & 0 & 1 \end{bmatrix}] \cdot r_{i} - [\begin{bmatrix} r_{x i} \\ r_{y i} \\ 0 \\ 1 \end{bmatrix}])}^{2},

(10)

which is simplified to the following form:

Q (ψ, θ, t_{z}) = \sum_{i = 1}^{m} {(- S θ \cdot r_{x i} + C θ S ψ \cdot r_{y i} + C θ C ψ \cdot r_{z i} + t_{z})}^{2} .

(11)

The set of parameters, $(ψ, θ, t_{z})$ , on which the objective function depends, is a reduced version of t, due to the application of (7) to the planar geometry. Then, the reduced set of transformation parameters for flatness assessment is $t = [ψ, θ, t_{z}]$ . The best fit parameters provided by algorithm are $\tilde{t} = [\tilde{ψ}, \tilde{θ}, {\tilde{t}}_{z}]$ . Thus, $T (\tilde{t}) \cdot r_{i}$ , for $i = 1, \dots, m$ , are the measured points transformed to best fit the reference element of the nominal integral feature from which they were sampled.

After the transformation of best fit, it is possible to estimate the flatness tolerance by the application of the tolerance definition supplied by the ASME Y 14.5 standard [9]. ASME Y 14.5 specifies the maximum peak-to-valley distance h_t with reference to the nominal plane as given below

h_{t} = | e_{\max} | - | e_{\min} | .

(12)

On the basis of the previous equation, the flatness tolerance can be estimated directly from the points transformed according to the computed parameters of best fitting. Therefore, the flatness tolerance is given by

T_{f} = \max_{i = 1, . . ., m} {r_{z i}^{\tilde{t}}} - \min_{i = 1, . . ., m} {r_{z i}^{\tilde{t}}},

(13)

where $r_{z i}^{t}$ is the third coordinate of the ith measured point after the best fit transformation, $T (\tilde{t}) \cdot r_{i}$ .

The equation of the best fitting plane in the form $a x + b y + c z + d = 0$ can be obtained transforming the versor normal to the x-y plane, ${[0,0, 1,0]}^{T}$ , through the generalized displacement: $\tilde{t} = [- \tilde{ψ}, - \tilde{θ}, - {\tilde{t}}_{z}]$ .

The product of the homogeneous matrix, representing the reverse best fit transformation, and the normal versor is

[\begin{bmatrix} C \tilde{θ} & S \tilde{θ} S \tilde{ψ} & - S \tilde{θ} C \tilde{ψ} & 0 \\ 0 & C \tilde{ψ} & S \tilde{ψ} & 0 \\ S \tilde{θ} & - C \tilde{θ} S \tilde{ψ} & C \tilde{θ} C \tilde{ψ} & - {\tilde{t}}_{z} \\ 0 & 0 & 0 & 1 \end{bmatrix}] \cdot [\begin{bmatrix} 0 \\ 0 \\ 1 \\ 0 \end{bmatrix}] = [\begin{bmatrix} - S \tilde{θ} C \tilde{ψ} \\ S \tilde{ψ} \\ C \tilde{θ} C \tilde{ψ} \\ 0 \end{bmatrix}] .

(14)

Reminding that the coefficients a, b, and c for a plane are proportional to the components of the versor normal to the plane, their values are obviously derived from (14): $a = - S \tilde{θ} C \tilde{ψ}, b = S \tilde{ψ}, c = C \tilde{θ} C \tilde{ψ}$ .

The coefficient d can be calculated forcing the plane to pass for the origin of the coordinate frame, ${[0,0, 0,1]}^{T}$ , transformed according to $[- \tilde{ψ}, - \tilde{θ}, - {\tilde{t}}_{z}]$ , that is, for the point ${[0,0, 0, - {\tilde{t}}_{z}, 0]}^{T}$ . Then, the plane coefficient d is given by the solution of the following equation:

a \cdot 0 + b \cdot 0 - c \cdot {\tilde{t}}_{z} + \tilde{d} = 0 \Rightarrow \tilde{d} = C \tilde{θ} \cdot C \tilde{ψ} \cdot {\tilde{t}}_{z} .

(15)

3.2. Cylindricity Assessment

The cylindrical class of symmetry encompasses all surfaces invariant under a translation and a rotation around the reference element. Their reference element is their axis of rotation, which is unique.

A straight line is the geometric element from which the proposed algorithm minimizes the distance of the transformed points. For the sake of simplicity the z-axis of the cylinder reference frame is chosen, that is, the straight line of equation x = 0 and y = 0. The nearest point to r_i on the datum is the point $n_{i} = {[0,0, r_{z i}^{}, 1]}^{T}$ . Equation (7) can be simplified considering that the objective function in this case is invariant for rotation and translation of the cloud of measured points around z-axis and along z-axis, respectively. In this way, the number of parameters to be optimised is reduced from 6 to 4, because it is possible to drop out ϕ and t_z from (7), imposing ϕ and t_z equal to zero. Then, the generic objective function, (7), is replaced for the cylindrical class of symmetry by the following one:

\begin{matrix} Q (ψ, θ, t_{x}, t_{y}) & = \sum_{i = 1}^{m} {(C θ \cdot r_{x i} + S θ S ψ \cdot r_{y i} + S θ S ψ \cdot r_{z i} + t_{x})}^{2} \\ + {(C ψ \cdot r_{y i} - C ψ \cdot r_{z i} + t_{y})}^{2} . \end{matrix}

(16)

The set of parameters, $(ψ, θ, t_{x}, t_{y})$ , on which the objective function depends, is a reduced version of t, due to the application of (7) to the cylindrical geometry. Then, the reduced set of transformation parameters for cylindricity assessment is $t = [ψ, θ, t_{x}, t_{y}]$ . The best fit parameters provided by the algorithm are $\tilde{t} = [\tilde{ψ}, \tilde{θ}, {\tilde{t}}_{x}, {\tilde{t}}_{y}]$ . Thus, $T (\tilde{t}) \cdot r_{i}$ , for $i = 1, \dots, m$ , are the measured points transformed to best fit the datum of the nominal integral feature from which they were sampled.

Afterwards the individuation of the best fit set of transformation parameters, the cylindricity of the inspected feature can be assessed following the tolerance definition supplied by the ASME Y 14.5 standard [9]. Thus, the cylindricity of the real feature can be estimated from the points transformed according to the best fit transformation parameter, $r_{i}^{t} = {[r_{x i}^{\tilde{t}}, r_{y i}^{\tilde{t}}, r_{z i}^{\tilde{t}}, 1]}^{T}$ , as

T_{c} = \max_{i = 1, . . ., m} {\sqrt{{(r_{x i}^{\tilde{t}})}^{2} + {(r_{y i}^{\tilde{t}})}^{2}}} - \min_{i = 1, . . ., m} {\sqrt{{(r_{x i}^{\tilde{t}})}^{2} + {(r_{y i}^{\tilde{t}})}^{2}}} .

(17)

It can sometimes be useful for comparison purpose to find the axis of the measured cylinder inside the cylinder reference frame. The position of this axis can be obtained transforming the versor of the z-axis, ${[0,0, 1,0]}^{T}$ , through the generalized displacement: $\tilde{t} = [- \tilde{ψ}, - \tilde{θ}, - {\tilde{t}}_{x}, - {\tilde{t}}_{y}, 0]$ , representing the transformation inverse to the best fit transformation. The product of the corresponding homogeneous matrix and the z-axis versor is

[\begin{bmatrix} C \tilde{θ} & - S \tilde{θ} S \tilde{ψ} & - S \tilde{θ} C \tilde{ψ} & - {\tilde{t}}_{x} \\ 0 & C \tilde{ψ} & S \tilde{ψ} & - {\tilde{t}}_{y} \\ S \tilde{θ} & - C \tilde{θ} S \tilde{ψ} & C \tilde{θ} C \tilde{ψ} & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}] \cdot [\begin{bmatrix} 0 \\ 0 \\ 1 \\ 0 \end{bmatrix}] = [\begin{bmatrix} - S \tilde{θ} C \tilde{ψ} \\ S \tilde{ψ} \\ C \tilde{θ} C \tilde{ψ} \\ 0 \end{bmatrix}] .

(18)

3.3. Sphericity Assessment

The spherical class of symmetry encompasses all surfaces invariant under all rotations about the feature centre (spherical displacement). The reference element of such a surface is a point corresponding to its centre.

The geometric element from which the proposed algorithm minimizes the distance of the transformed points will be, as a consequence, a point. For the sake of simplicity the point $n_{i} = {[0,0, 0,1]}^{T}$ , that is, the origin of the sphere reference system, is chosen. Equation (7) can be simplified considering that the objective function in this case is invariant for rotations of the cloud of measured points around x-, y-, and z-axis. In this way, the number of parameters to be optimised is reduced from 6 to 3, because it is possible to drop out ψ, θ, and ϕfrom (7), imposingψ, θ, and ϕ equal to zero. Then, the objective function becomes

Q (t_{x}, t_{y}, t_{z}) = \sum_{i = 1}^{m} {[\sqrt{{(t_{x} + r_{x i})}^{2} + {(t_{y} + r_{y i})}^{2} + {(t_{z} + r_{z i})}^{2}}]}^{2} .

(19)

The set of parameters, $(t_{x}, t_{y}, t_{z})$ , on which the objective function depends, is a reduced version of t, due to the application of (7) to the spherical geometry. Then, the reduced set of transformation parameters for sphericity assessment is $t = [t_{x}, t_{y}, t_{z}]$ . The best fit parameters provided by the algorithm are $\tilde{t} = [{\tilde{t}}_{x}, {\tilde{t}}_{y}, {\tilde{t}}_{z}]$ . Thus, $T (\tilde{t}) r_{i}$ , for $i = 1, \dots, m$ , are the measured points transformed to best fit the datum of the nominal integral feature from which they were sampled.

Spherical features are feature of size. They are continuous surfaces with possible size, form, and location errors. International Standards, both ISO and ASME [8, 9], do not deal with the sphericity tolerance explicitly and, therefore, the concepts of form tolerance present in these standards have to be employed to the control of the sphericity errors, [18, 26, 27]. The profile tolerancing, for line or surfaces, is considered one of the most versatile controls available. It can be used without datums to control the form of a surface. In addiction, when it is used with one or more datums, it contains a degree of orientation and form control. So profile tolerance is used to control the sphericity errors of the feature. After the best fitting step, the origin of the CMM reference system has become the datum of the measured data. Then, the sphericity tolerance can be calculated from the transformed data by means of the following equation:

\begin{matrix} T_{c} & = \max_{i = 1, . . ., m} {\sqrt{{(r_{x i}^{\tilde{t}})}^{2} + {(r_{y i}^{\tilde{t}})}^{2} + {(r_{z i}^{\tilde{t}})}^{2}}} \\ - \min_{i = 1, . . ., m} {\sqrt{{(r_{x i}^{\tilde{t}})}^{2} + {(r_{y i}^{\tilde{t}})}^{2} + {(r_{z i}^{\tilde{t}})}^{2}}} . \end{matrix}

(20)

From the set of best fitting parameters, $\tilde{t} = [{\tilde{t}}_{x}, {\tilde{t}}_{y}, {\tilde{t}}_{z}]$ , the position of the centre of the measured spherical feature can be obtained transforming the origin of the reference system, ${[0,0, 0,1]}^{T}$ , through the generalized displacement: $\tilde{t} = [- {\tilde{t}}_{x}, - {\tilde{t}}_{y}, - {\tilde{t}}_{z}]$ . This displacement represents the inverse of transformation that best fits the measured points to the origin. The position of the centre in the CMM reference system is given by

[\begin{bmatrix} 0 & 0 & 0 & - {\tilde{t}}_{x} \\ 0 & 0 & 0 & - {\tilde{t}}_{y} \\ 0 & 0 & 0 & - {\tilde{t}}_{z} \\ 0 & 0 & 0 & 1 \end{bmatrix}] \cdot [\begin{bmatrix} 0 \\ 0 \\ 0 \\ 1 \end{bmatrix}] = [\begin{bmatrix} - {\tilde{t}}_{x} \\ - {\tilde{t}}_{y} \\ - {\tilde{t}}_{z} \\ 1 \end{bmatrix}] .

(21)

In this way the location errors of the feature centre can be evaluated.

4. Applications

The above algorithm was implemented on an AMD Athlon, 523 MB, 1800 MHz computer using Mathematica 5.0. The proposed method was tested to verify its correctness using both datasets published in previous papers, for comparative purposes, and datasets created by numerically simulating actual planar, cylindrical, and spherical surfaces. The formers were published in [16, 28] for planar surfaces, in [17, 28] for cylindrical surfaces, and in [29] for spherical surfaces. The latter were generated using an appropriate algorithm similar to that used in [16, 17]. This algorithm samples points from two nominal features at a given distance, so that the form tolerance is a priori known. For example, to simulate measured points from an actual cylinder, two cylinders with the same axis and with a radius difference equal to the imposed cylindricity error are generated. The measured points are extracted from the surfaces of the two coaxial cylinders. Furthermore, a rigid motion is applied to the cloud of measured points; so that it can be positioned and orientated at will into the CMM reference system. The proposed algorithm was tested on both sparse measurements, that is, measured points randomly generated on the surface, and grid points, that is, measured points generated using a grid superimposed on the surface. The results for these two kinds of datasets were approximately identical. Therefore, this paper shows only the results obtained by a grid of measured points. However, the datasets in [16, 17, 28, 29] are sparse measurements. For the presented applications, the adopted convergence criterion was the reduction of the magnitude $| Q^{j} - {Q^{j}}^{- 1} | < ε$ below 10⁻⁵.

The results of these evaluations using the proposed algorithm for each feature class are reported in Tables 2, 3 and 4. In each Table, the values reported in the column Nominal Tolerance are characterized by an indication, reporting the origin of the values. The letter (G) means that the nominal tolerance refers to the relative distance of the two parallel surfaces from which the dataset was generated, (CMM) indicates that the nominal tolerance was evaluated using a least squares technique on a CMM software (the name of the software used can be found in the cited paper), and (MZ) indicates that the nominal value comes from a minimum zone calculation (the method employed is in the corresponding paper).

Table 2:

Flatness results.

	no. of data points	Nominal tolerance (mm)	Proposed algorithm (mm)	Zeiss algorithm (mm)	CPU time (s)

Dataset 1	20	0.1 (G)	0.1151	0.1485	0.391
Dataset 2	40	0.1 (G)	0.1140	0.1137	0.797
Dataset 3	80	0.1 (G)	0.1210	0.1354	2.000
Dataset 4	120	0.1 (G)	0.1095	0.1281	3.984
Dataset 5	240	0.1 (G)	0.1191	0.1270	9.093
Dataset 6	16	0.2 (G)	0.2290	0.2551	0.315
Dataset 7	30	0.05 (G)	0.0563	0.0684	0.532
Dataset 8	48	0.005 (G)	0.0052	0.0056	0.983
Dataset 9 [16]	15	0.0025000 (MZ)	0.0026	0.0029	0.303
Dataset 10 [16]	25	0.0048636 (MZ)	0.0051	0.0059	0.508
Dataset 11 [28]	12	0.0002 (MZ)	0.0004	0.0005	0.252
Dataset 12 [28]	20	0.0027 (MZ)	0.0027	0.0028	0.366
Dataset 13 [28]	15	0.0002 (MZ)	0.0002	0.0002	0.285
Dataset 14 [28]	20	0.0032 (MZ)	0.0032	0.0033	0.401
Dataset 15 [28]	20	0.0025 (MZ)	0.0025	0.0026	0.387

Table 3:

Cylindricity results.

	no. of data points	Nominal tolerance (mm)	Proposed algorithm		Zeiss algorithm		CPU time (s)
	no. of data points	Nominal tolerance (mm)	tolerance (mm)	r (mm)	tolerance (mm)	r (mm)	CPU time (s)

Dataset 1	20	0.1 (G)	0.1499	5.0450	0.1499	5.0450	2.921
Dataset 2	40	0.1 (G)	0.1494	5.0475	0.1495	5.0475	10.969
Dataset 3	80	0.1 (G)	0.1433	5.0512	0.1433	5.0513	45.141
Dataset 4	120	0.1 (G)	0.1246	5.0375	0.1246	5.0375	68.235
Dataset 5	240	0.1 (G)	0.1267	5.0213	0.1306	5.0463	224.236
Dataset 6	40	0.2 (G)	0.3038	5.0800	0.3048	5.0801	14.570
Dataset 7	30	0.05 (G)	0.0733	5.0225	0.0733	5.0225	17.302
Dataset 8	240	0.005 (G)	0.0056	5.0024	0.0068	5.0025	239.420
Dataset 8b	204	0.02 (G)	0.0210	5.0107	0.0231	5.0124	120.657
Dataste 9 [17]	40	0.001 (MZ)	0.0040	30.0008	0.0048	30.0006	18.203
Dataset 10 [17]	20	0.19667 (MZ)	0.2001	60.0031	0.2140	60.0096	5.109
Dataset 11 [28]	16	0.0067 (CMM)	0.0060	0.1491	0.0067	0.1495	4.189
Dataset 12 [28]	18	0.0006 (CMM)	0.0004	0.4946	0.0005	0.4947	4.699
Dataset 13 [28]	16	0.0051 (CMM)	0.0047	0.5353	0.0051	0.5354	4.218
Dataset 14 [28]	24	0.0004 (CMM)	0.0003	0.1998	0.0003	0.1999	5.362
Dataset 15 [28]	18	0.0007 (CMM)	0.0151	0.5018	0.0152	0.5011	4.362

Table 4:

Sphericity results.

	no. of data points	Nominal tolerance (mm)	Proposed algorithm		Zeiss algorithm		CPU time (s)
	no. of data points	Nominal tolerance (mm)	tolerance (mm)	r (mm)	tolerance (mm)	r (mm)	CPU time (s)

Dataset 1	20	0.1 (G)	0.1000	6.100	0.1282	6.060	0.485
Dataset 2	40	0.1 (G)	0.1000	6.000	0.1471	6.045	0.375
Dataset 3	80	0.1 (G)	0.1000	6.100	0.1599	6.053	0.797
Dataset 4	120	0.1 (G)	0.1000	6.100	0.1314	6.052	1.312
Dataset 5	240	0.1 (G)	0.1000	6.000	0.1166	6.048	4.031
Dataset 6	960	0.2 (G)	0.2346	6.095	0.2346	6.095	10.564
Dataset 6b	800	0.1 (G)	0.1039	6.052	0.1040	6.052	7.971
Dataset 7	30	0.05 (G)	0.0500	6.050	0.0853	6.028	0.308
Dataset 8	48	0.005 (G)	0.0050	6.000	0.0072	6.002	0.385
Dataset 9b	800	0.1 (G)	0.1168	6.075	0.1168	6.075	14.954
Dataset 10 [29]	50	0.01 (G)	0.0084	1.004	0.0085	1.004	1.062
Dataset 11 [30]	40	0.008327 (MZ)	0.0085	0.999	0.0091	0.999	0.329
Dataset 12 [30]	36	0.009669 (MZ)	0.0100	50.034	0.0101	50.034	0.299

Datasets 9 and 10 for flatness assessment, Table 2, correspond to datasets 1 and 2 reported in [16] (generated datasets), while datasets 11–15 are the examples 1–5 of [28] (CMM-flatness datasets). Datasets 9 and 10 for cylindricity assessment, Table 3, correspond to the datasets 1 and 2 reported in [17] (generated datasets), while datasets 11–15 correspond to the examples 1–5 in [28] (CMM-cylindricity datasets). Dataset 10 for sphericity assessment, Table 4, corresponds to the dataset presented in [29] (generated) and datasets 11–12 are the datasets 2-3 in [30] (measurements). The other datasets were generated for each class of symmetry. Five of them are characterized by the same form tolerance but with increasing number of measured points to verify the efficiency and the rate of convergence of the algorithm. The remaining datasets have different orientation inside the CMM reference system and different tolerance form of the feature to test the consistency of the algorithm to different feature characteristics. All the datasets marked with b are measurements collected on a portion of a feature. The obtained results were compared with the nominal value of tolerance.

4.1. Flatness Results

The results of the planar cases are presented in Table 2. Datasets 1–5 prove that a linear relationship exists between the number of measured points and the amount of time, in CPU seconds, the flatness algorithm requires to compute the best fit solution, as shown in Figure 1(a). These results indicate that the computational effort is a linear function of the number of measured points.

Efficiency of proposed method: (a) flatness assessment, (b) cylindricity assessment, and (c) sphericity assessment.

Figure 2 shows two views of the best fitting procedure applied to the dataset 6. As it can be seen, the measured points move iteration after iteration towards their datum, that is, the x-y plane.

Best fit of a planar feature (two views of the dataset 6).

For datasets 9 and 10 in the nominal tolerance column, the flatness tolerance computed with the minimum zone method proposed in [16] is reported, because in the cited paper no nominal tolerance is provided for these datasets.

The proposed algorithm gives flatness estimates that are very close to the values reported in the nominal tolerance column. Besides, the results of the proposed algorithm are nearer to those due to minimum zone than the values obtained by Zeiss algorithm [31]. Then, it is proven that the proposed method allows to overestimate flatness less than a traditional least squares algorithm, such as Zeiss, does.

4.2. Cylindricity Results

Table 3 shows the results of the cylindricity assessment. The radius of the assessment feature is also given in Table 3. The amount of time needed to compute the cylindricity tolerance of datasets 1–5 indicates that the computational effort is also in this case a nearly linear function of the number of measured points, as shown in Figure 1(b). However, the CPU time required for cylindricity assessment is remarkably larger than this required for flatness and sphericity assessment due to the major complexity of the algorithm.

Figure 3 shows two views of the best fitting procedure applied to the dataset 8. The measured points, moving iteration after iteration towards their datum, that is, the z-axis, are clearly visible.

Best fit of a cylindrical feature (two views of the dataset 8).

The proposed algorithm gives cylindricity estimates that are very near to the values reported in the nominal tolerance column. Moreover, the results of the proposed algorithm are closer to minimum zone results than the values obtained by Zeiss algorithm. This means that the proposed method allows to overestimate cylindricity less than a traditional least squares algorithm, such as Zeiss, does.

4.3. Spherical Results

Table 4 shows the form errors computed for spherical features. The computational effort is still an almost linear function of the number of data points; see Figure 1(c). The precision of the algorithm for this class of feature is considerably larger than that feasible for the other classes of symmetry because the best fit requires only a translation of the measured points.

Figure 4 shows the best fitting procedure applied to the dataset 6. The measured points, moving iteration after iteration towards their datum, that is, the origin of the reference frame, are clearly visible.

Figure 4

Best fit of a spherical feature (two views of the dataset 6).

The proposed algorithm gives sphericity estimates that are very near to the values reported in the nominal tolerance column. Moreover, the results of the proposed algorithm are nearer to the minimum zone results than the values obtained by Zeiss algorithm. This surely means that the proposed method allows to overestimate sphericity less than a traditional least squares algorithm, such as Zeiss, does.

The algorithm is also applicable to measured data obtained from portion of surface. Figure 5 shows two views of the best fitting procedure applied to the dataset 8b, which was collected from a portion of a cylinder. Figure 6 shows the algorithm applied to the best fit of a portion of a spherical feature (two views of the dataset 6b).

Figure 5

Best fit of a portion of a cylindrical feature (two views of the dataset 8b).

Best fit of a portion of a spherical feature (two views of the dataset 6b).

5. Conclusions

The present work shows a new algorithm to evaluate flatness, cylindricity, and sphericity tolerance that is based on surface invariance with regard to rigid motions. The new algorithm transforms the measured coordinates, through homogeneous transformation matrices, in order to best fit the reference element (datum) of the class of the surface from which the actual measurements were sampled.

The proposed approach was applied to a set of cases in order to validate the effectiveness and the robustness of the approach. The obtained results indicate that the proposed algorithm provides accurate and quick assessment.

In all cases the new method results in a smaller zone values than the classical least squares methods. The proposed algorithm gives flatness, cylindricity, and sphericity estimates that are very near to minimum zone results. Moreover, the results of the proposed algorithm are nearer to those due to minimum zone than the values obtained by Zeiss algorithm. It is demonstrated that the proposed method allows to overestimate form tolerance less than a traditional least squares algorithm, such as Zeiss, does. Then, it is proven that the proposed method is a favourable alternative to the classical least squares algorithms for form tolerance assessment.

References

Homepage of ISO/TC 213, President: Nielsen

Henrik S.

Dr. , 10219 Coral Reef Way, Indianapolis, IN 46256, USA, November 2008, http://isotc213.ds.dk.

ISO/TS 17450:1999, “Geometrical product specification (GPS)—model for geometric specification and verification,” International Organization for Standardization, Geneva, Switzerland, 1999.

Srinivasan

, “A geometrical product specification language based on a classification of symmetry groups,” CAD Computer Aided Design, vol. 31, no. 11, pp. 659–668, 1999.

Desrochers

and Rivière

, “A matrix approach to the representation of tolerance zones and clearances,” International Journal of Advanced Manufacturing Technology, vol. 13, no. 9, pp. 630–636, 1997.

Clemant

Riviere

, and Temmerman

, Cotation Tridimensionelle des Systemes Mecaniques, Ivry-sur-Seine, Paris, France, PYC edition, 1994.

Salomons

O. W.

, Specifying and solving tolerance constraints based on the TTRS concept, computer support in the design of mechanical products, Ph.D. dissertation, University of Twente, Twente, The Netherlands, 1995.

Prisco

and Giorleo

, “Overview of current CAT systems,” Integrated Computer-Aided Engineering, vol. 9, no. 4, pp. 373–387, 2002.

ISO 1101, Technical drawings-Geometrical tolerancing-Tolerancing of form, orientation, location and run out-Generalities, definitions, symbols, indication on drawings.

ANSI Y14.5M-1982, National Standard on Dimensioning and Tolerancing, ASME, New York, NY, USA, 1982.

10.

Capello

and Semeraro

, “The harmonic fitting method for the assessment of the substitute geometry estimate error. Part I: 2D and 3D theory,” International Journal of Machine Tools and Manufacture, vol. 41, no. 8, pp. 1071–1102, 2001.

11.

Murthy

T. S. R.

and Abdin

S. Z.

, “Minimum zone evaluation of surfaces,” International Journal of Machine Tool Design and Research, vol. 20, no. 2, pp. 123–136, 1980.

12.

Shummugam

M. S.

, “On assessment of geometry errors,” International Journal of Production Research, vol. 245, pp. 413–425, 1986.

13.

Shunmugam

M. S.

, “New approach for evaluating form errors of engineering surfaces,” Computer-Aided Design, vol. 19, no. 7, pp. 368–374, 1987.

14.

Wang

, “Minimum zone evaluation of form tolerances,” Manufacturing Review, vol. 5, no. 3, pp. 213–220, 1992.

15.

Roy

and Zhang

, “Establishment of a pair of concentric circles with the minimum radial separation for assessing roundness error,” Computer-Aided Design, vol. 24, no. 3, pp. 161–168, 1992.

16.

Carr

and Ferreira

, “Verification of form tolerances part I: basic issues, flatness, and straightness,” Precision Engineering, vol. 17, no. 2, pp. 131–143, 1995.

17.

Carr

and Ferreira

, “Verification of form tolerances part II: cylindricity and straightness of a median line,” Precision Engineering, vol. 17, no. 2, pp. 144–156, 1995.

18.

Samuel

G. L.

and Shunmugam

M. S.

, “Evaluation of sphericity error from form data using computational geometric techniques,” International Journal of Machine Tools and Manufacture, vol. 42, no. 3, pp. 405–416, 2002.

19.

Damodarasamy

and Anand

, “Evaluation of minimum zone for flatness by normal plane method and simplex search,” IIE Transactions, vol. 31, no. 7, pp. 617–626, 1999.

20.

Lai

J.-Y.

and Chen

I.-H.

, “Minimum zone evaluation of circles and cylinders,” International Journal of Machine Tools and Manufacture, vol. 36, no. 4, pp. 435–451, 1996.

21.

Liu

C.-H.

Jywe

W.-Y.

, and Chen

C.-K.

, “Quality assessment on a conical taper part based on the minimum zone definition using genetic algorithms,” International Journal of Machine Tools and Manufacture, vol. 44, no. 2–3, pp. 183–190, 2004.

22.

Yau

H.-T.

and Menq

C.-H.

, “A unified least-squares approach to the evaluation of geometric errors using discrete measurement data,” International Journal of Machine Tools and Manufacture, vol. 36, no. 11, pp. 1269–1290, 1996.

23.

Choi

and Kurfess

T. R.

, “Dimensional measurement data analysis, part 1: a zone fitting algorithm,” Journal of Manufacturing Science and Engineering, Transactions of the ASME, vol. 121, no. 2, pp. 238–245, 1999.

24.

Goch

and Renker

H. J.

, “Efficient multi-purpose algorithm for approximation and alignment problems in coordinate measurement techniques,” Annals of CIRP, vol. 39, no. 1, pp. 553–556, 1990.

25.

Gill

P. R.

Murray

, and Wright

M. H.

, “The Levenberg-Marquardt method,” in Practical Optimization, pp. 136–137, Academic Press, London, UK, 1981.

26.

Samuel

G. L.

and Shunmugam

M. S.

, “Evaluation of circularity and sphericity from coordinate measurement data,” Journal of Materials Processing Technology, vol. 139, no. 1–3, pp. 90–95, 2003.

27.

Samuel

G. L.

and Shunmugam

M. S.

, “Evaluation of sphericity error from coordinate measurement data using computational geometric techniques,” Computer Methods in Applied Mechanics and Engineering, vol. 190, no. 51–52, pp. 6765–6781, 2001.

28.

Weber

Motavalli

Fallahi

, and Cheraghi

S. H.

, “A unified approach to form error evaluation,” Precision Engineering, vol. 26, no. 3, pp. 269–278, 2002.

29.

Fana

K.-C.

and Lee

J.-C.

, “Analysis of minimum zone sphericity error using minimum potential energy theory,” Precision Engineering, vol. 23, no. 2, pp. 65–72, 1999.

30.

Wen

and Song

, “An immune evolutionary algorithm for sphericity error evaluation,” International Journal of Machine Tools and Manufacture, vol. 44, no. 10, pp. 1077–1084, 2004.

31.

Calypso Guidebook, 2005, 3D Metrology Services, Carl Zeiss 3D Metrology Services, Aalen, 2-25 2-26.