Three-Point Inverse and Forward Kinematic Algorithms for Circle Measurement from Distributed Displacement Sensor Network

Automatic fitting of an arc center and radius is a quality problem frequently encountered when manufacturing a mechanical component. Due to the complexity of the measurement, validating each manufactured component via inspection is not feasible or economical. This paper introduces a new validation procedure for measuring arcs from distributed sensors. The goal of this proposed measurement process is to improve measurement throughput (i.e., parts measured per unit of time) and reduce measurement errors associated with hardware and algorithms. This proposed model develops a three-point inverse kinematic algorithm (TPIK) accompanied by a calibration master to obtain the relative location of the measurement system by solving a set of six non-linear equations. This technique allows deployment of a high accuracy gauge systems that in general, reduces machine and algorithm errors. The direct fitting is validated by using mathematical, CAD, and experimental models. Furthermore, a modified definition for the roundness index is introduced based on the proposed forward and inverse algorithms. The simulations examine the roundness index in relation to the measurement precision, sampling angle, nominal radius, and part variation. A benefit of this proposed method is accurate and rapid inspection of the radii and elimination of the human error associated with part loading variation during conventional radii measurement. The rapid, accurate inspection and corresponding reduction in human error make this method an excellent process for inspection of large quantities of components.


Introduction
Fitting surface geometry to a part is the process of estimating the interpolated geometry between the real data points collected using metrology instruments on real components. This method is an efficient way to develop a geometric model of the feature. The problem of fitting a circular arc (say, for a fillet or a section of a spherical mirror or cylinder) has motivated a large amount of literature in science and engineering with applications to microwave measurement [1], computer vision [2], design tolerance [3], metrology and inspection of mechanical parts [4,5], archaeology [6], geodesy [7], and many others. For example, the manufacturing precision for a camshaft in an internal combustion engine significantly correlates to the performance of the engine [8].
In recent years, a significant issue has been identified with automated inspection techniques that utilize measurement equipment such as coordinate measuring machines (CMMs) [3]. For example, to determine the diameter of a cylinder, a CMM is programmed to acquire points on the surface and then optimize the parameters for best-fit geometry. This type of measurement can be affected i.
Three-point inverse kinematic algorithm (TPIK): a calibration procedure that references the local coordinates of each sensor into one global coordinate system (G). This produces system parameters whose values are based on a numerical solution of mixed non-linear equations. Further discussion to follow in Section 2.2. ii.
Three-point forward kinematic algorithm (TPFK): a closed form solution that algebraically fits the transformed data points through the circumference of a circle. Further discussion to follow in Section 2.1.
3 coordinates of each sensor into one global coordinate system (G). This produces system parameters whose values are based on a numerical solution of mixed non-linear equations. Further discussion to follow in Section 2.2. ii.
Three-point forward kinematic algorithm (TPFK): a closed form solution that algebraically fits the transformed data points through the circumference of a circle. Further discussion to follow in Section 2.1.
The TPFK algorithm uses the standard representation of circle with a radius and a center. The goal is to extract the radius regardless of where the center is located relative to the global coordinate system. The TPIK algorithm suggests deployment of a calibration master of known geometry to estimate location of the sensors. Additionally, TPIK provides a mean for an operator to calibrate the overall offset caused by the measurement system, which could be due to linear or non-linear shifts in sensor reading, or unknown relative movements of sensors. Therefore, it is recommended to carry out periodic calibration to avoid accumulation drifts in the measurements.

Three-Point Forward Kinematics (TPFK) Algorithm-Recovery of Circle Geometry
The TPFK algorithm is a deterministic method that calculates the radius and the center ( , ) of a circle from point coordinates measured by spatially distributed displacement sensors { , , } shown in Figure 1. The sensors are arranged to measure the vertical displacements { , , } relative to their local coordinates. The position of the two local coordinate sensors { , } are measured with respect to a global coordinate system { }, which we chose to fix at sensor { }. Therefore, the vector positions of the points located at the arc { , , } can be expressed in reference to the global coordinate system as { (0, ), ( , + ), ( , + )}. The radius is defined by the linear distance between any point on the circle and its center ( , ), thus can be expresses as: where the linear distance between a point on a circle ( , ) and its center is simply obtained from three data points.
The circle center ( , ) is referenced to the global coordinate. Therefore, rewriting the Equation (1) yields: Figure 1. Three-point method: Measurement of the circle geometry by using three arbitrary points on an arc {P 1 , P 2 , P 3 } using global reference G located at s 1 . The sensors are arranged in parallel to reduce system complexity and set-up error.
The TPFK algorithm uses the standard representation of circle with a radius and a center. The goal is to extract the radius regardless of where the center is located relative to the global coordinate system. The TPIK algorithm suggests deployment of a calibration master of known geometry to estimate location of the sensors. Additionally, TPIK provides a mean for an operator to calibrate the overall offset caused by the measurement system, which could be due to linear or non-linear shifts in sensor reading, or unknown relative movements of sensors. Therefore, it is recommended to carry out periodic calibration to avoid accumulation drifts in the measurements.

Three-Point Forward Kinematics (TPFK) Algorithm-Recovery of Circle Geometry
The TPFK algorithm is a deterministic method that calculates the radius R and the center C(a, b) of a circle from point coordinates measured by spatially distributed displacement sensors {s 1 , s 2 , s 3 } shown in Figure 1. The sensors are arranged to measure the vertical displacements y 1 , y 2 , y 3 relative to their local coordinates. The position of the two local coordinate sensors {s 2 , s 3 } are measured with respect to a global coordinate system {G}, which we chose to fix at sensor {s 1 }. Therefore, the vector positions of the points located at the arc {P 1 , P 2 , P 3 } can be expressed in reference to the global coordinate system as P 1 (0, y 1 ), P 2 (x 2o , y 2o + y 2 ), P 3 (x 3o , y 3o + y 3 ) . The radius is defined by the linear distance between any point on the circle and its center C(a, b), thus can be expresses as: where the linear distance between a point on a circle P(x, y) and its center is simply obtained from three data points. The circle center C(a, b) is referenced to the global coordinate. Therefore, rewriting the Equation (1) yields: and Upon expanding and simplifying the Equations (3) and (4), the results is a linear equation in Upon solving for → C(a, b) = M −1 B, a point can be used in the circle such as P 1 (0, y 1 ) to calculate circle radius from Equation (2).

Three-Point Inverse Kinematics (TPIK) Algorithm-Calibration Procedures
The TPFK problem discussed in the previous section requires prior knowledge of the position S p x 2o , y 2o , x 3o , y 3o , which is the origin of the sensors {s 2 , s 3 }. At the origin, one assumes that all initial displacement y 1 , y 2 , y 3 values are reset or forced to read zero. Often, for various reasons, it is impractical to obtain the S values by using direct measurement. Barriers to direct, efficient measurement include a lack of proper tooling, desired level of measurement accuracy is not obtainable, excessive set-up time, set up is not conducive for the production environment, measurement system idle time attributed to quality protocol, and change of operator and/or the part. As a result, there is a need to conduct a periodic calibration procedure that resets the errors because of the controlled or the uncontrolled factors that result in changing the value of S. These factors could be attributed to the mechanical movement of sensors, retooling, thermal expansion, changes in the characteristics of the sensors, and many other unpredictable factors. Therefore, in this section a systematic calibration procedure is proposed that inversely calculates S p x 2o , y 2o , x 3o , y 3o , which must remain unchanged during measurement. The S p represents part of the system parameters that are later plugged into Equation (5) to calculate the center of the circle → C(a, b), and its radius by using Equation (2). It should be noted that TPFK as well as TPIK (will be shown later in this section) do not require information about the center of the calibration master.
A calibration master was proposed and depicted in Figure 2a. The calibration master is designed with a minimum of two different concentric radii. For example, The calibration master in Figure 2a is fabricated with three known concentric radii {R c1 , R c2 , R c3 }, whose design values depend on the size of the part being inspected and the range of the sensors used. The fabrication of a good calibration master with precise {R c2 , R c3 } is crucial to the success of the calibration and can be achieved either by additive manufacturing, such as 3D printing, or subtractive fabrication such as machining the calibration master from a thick cylinder shown in Figure 2b. However, the fabrication cannot guarantee accurate design values, therefore it is necessary to measure the radii {R c2 , R c3 } after fabrication at given operating conditions. The value of R c1 is not part of the calculation-as we will see later-however, it is used in the calibration to initialize the location of the sensors, S. The center of the calibration master C c relative to any global coordinate system G (G was chosen to be placed at sensor s 1 as shown in Figure 1) should remain unchanged during the calibration procedure. To better understand the design requirement of the calibration master in relation to the inspected part, it was assumed that the minimum sensing range was ∆ . Thus, the nominal radius value of the inspected part was ~ , and the calibration master should be ( − ) < ∆ and > > . The sensors were arranged in parallel such that each arc in the calibration master could be measured at once by all three sensors as shown in Figure 2c. The sensors were made adjustable in the y-direction to allow inspection of a wider range of parts. The shared center of the calibration master must be tightly secured in a fixed rotating axis; allowing the measurement of displacement for every arc segment while their shared center remains unchanged.
Following the discussion above, the derivation of the TPIK algorithm can be shown by using the following example: Assume that spring-loaded probes in Figure 2c are brought on contact and then pushed against the segment of the calibration master. This is to maximize the reading range of the sensors. Each sensor is fastened onto the fixture and then reset to read zero displacement. This step has established { , , , }, however the values are yet to be found. Equation (5) is applied to evaluate the center of the calibration master by either using segment or . Substituting the known value yields the following matrix equation: where { , , } displacement are measured on the segment. The above equation breaks down into its element form: In addition, substitute the center and radius of the segment in equation (2) to yield: and write as: Similarly, for , we can write: To better understand the design requirement of the calibration master in relation to the inspected part, it was assumed that the minimum sensing range was ∆L min . Thus, the nominal radius value of the inspected part was~R a , and the calibration master should be (R c3 − R c1 ) < ∆L min and R c3 > R a > R c1 . The sensors were arranged in parallel such that each arc in the calibration master could be measured at once by all three sensors as shown in Figure 2c. The sensors were made adjustable in the y-direction to allow inspection of a wider range of parts. The shared center C c of the calibration master must be tightly secured in a fixed rotating axis; allowing the measurement of displacement for every arc segment while their shared center remains unchanged.
Following the discussion above, the derivation of the TPIK algorithm can be shown by using the following example: Assume that spring-loaded probes in Figure 2c are brought on contact and then pushed against the R c1 segment of the calibration master. This is to maximize the reading range of the sensors. Each sensor is fastened onto the fixture and then reset to read zero displacement. This step has established S p x 2o , y 2o , x 3o , y 3o , however the values are yet to be found. Equation (5) is applied to evaluate the center of the calibration master C c by either using segment R c2 or R c3 . Substituting the known R c2 value yields the following matrix equation: x 2 2o − y 2 12 + y 2 22 + 2y 22 y 2o + y 2 2o x 2 3o − y 2 12 + y 2 32 + 2y 32 y 3o + y 2 3o B c2 (6) where y 12 , y 22 , y 32 displacement are measured on the R c2 segment. The above equation breaks down into its element form: In addition, substitute the center and radius of the R c2 segment in Equation (2) to yield: and write as: Similarly, for R c3 , we can write: x 2 2o − y 2 13 + y 2 23 + 2y 23 y 2o + y 2 2o x 2 3o − y 2 13 + y 2 33 + 2y 33 y 3o + y 2 where y 13 , y 23 , y 33 displacement is measured on the R c3 segment. The above equation breaks down into its element form: In addition, substitute the center and radius of the R c3 segment in Equation (2) to yield: which is expressed as: It must be noted that C c2 (1) = C c3 (1) = a c , and C c2 (2) = C c3 (2) = b c . Therefore Equations (7), (9), (11), and (13) provide six non-linear equations with a total of six unknowns S x 2o , y 2o , x 3o , y 3o, , a c , b c . A successful solution will provide an indirect way for calculating system parameter S = x 2o , y 2o , x 3o , y 3o, , which is sufficient for the TPFK algorithm. However, the subset solution {a c , b c } is unnecessary for the TPFK algorithm, and is only needed to complete the computation of TPIK algorithm. Likewise, the center of the inspected part in the TPFK algorithm is implicit in the computation of the radius where it can be introduced to the measurement system from any location as long as the sensors can read any three points on the surface. Consequently, fitting a circle by using these proposed methods reduces handling errors and reduces the measurement time.
Finally, it should be noted that there is no restriction on how the sensors are initialized, i.e., it is possible to replace the circular shape R c1 in the calibration master with a flat surface to level the offsets between the sensors { y 1o − y 3o ≈ y 1o − y 2o ≈ 0}. This could help the numerical solution of the TPIK algorithm to converge faster to a solution as will be observed later in the simulation sections.

Roundness Index Based on TPIK and TPFK Algorithms
Roundness was defined in the index in this paper by the average nominal radius R a that minimizes the square summation of residual errors for every arc-segment in Figure 3a whose radius R 3p,i is computed by using TPIK and TPFK algorithms. This measure is given by: where N is the total number of arc-segments. The optimization of Equation (14) with respect to R a , simply gives the average R a = N i=1 R 3p,i /N for the best fitted circle in Figure 3b. Similarly, the where a 3p,i and b 3p,i are obtained for each arc-segment in the TPFK algorithm. Roundness index is defined by β = 1 − √ W/N/R a , where a perfect unit circle has a value of 1. Alternatively, one can define β * = 1 − √ W/N/R e when the exact nominal radius R e is known. It should be noted that in ideal scenarios, R a = R e when all arc profiles has no perturbation, as will be shown in later sections. The computation of a β or β * requires a onetime calibration via the TPIK algorithm, and then the application of TPFK algorithm to every arc-segment around the particular part being inspected. The characteristics of the roundness as function of perturbation pattern(s) will be presented in later section. One advantage of this definition is that the roundness is well suited for the automatic inspection because it does not require a common center C(a, b) for all arc-segments, and can be achieved by various sampling patterns. Examples of sampling patterns encompasses non-overlapping or overlapping arc-segmentations as shown in the three cases in Figure 4. Numerical studies will be conducted in later sections to show the effect of the sampling angle θ on the roundness index.
Examples of sampling patterns encompasses non-overlapping or overlapping arc-segmentations as shown in the three cases in Figure 4. Numerical studies will be conducted in later sections to show the effect of the sampling angle on the roundness index.

Validation of the TPIK Algorithm by Mathematical Simulation
In this section, a mathematical model was constructed to simulate the y-displacements between sensor coordinates and parametric circles. For simplicity, it was assumed that the sensors were initialized at the same level, i.e., = = = 0. The calibration master was represented with parametric equations ( ( ), ( )) that simulate two concentric circles of a known center and two given radii { , } as shown in Figure 5. To find the y-displacements between the sensor coordinates and the parametric equations, substitute the x-location in = cos (( − )⁄ ) to find t, and then compute displacements from = − ( ). Typically, for a circle, there are two solutions obtained for t, and one must select the solution, which corresponds to the smallest . The y-displacement results are summarized in Table 1. 7 Examples of sampling patterns encompasses non-overlapping or overlapping arc-segmentations as shown in the three cases in Figure 4. Numerical studies will be conducted in later sections to show the effect of the sampling angle on the roundness index.

Validation of the TPIK Algorithm by Mathematical Simulation
In this section, a mathematical model was constructed to simulate the y-displacements between sensor coordinates and parametric circles. For simplicity, it was assumed that the sensors were initialized at the same level, i.e., = = = 0. The calibration master was represented with parametric equations ( ( ), ( )) that simulate two concentric circles of a known center and two given radii { , } as shown in Figure 5. To find the y-displacements between the sensor coordinates and the parametric equations, substitute the x-location in = cos (( − )⁄ ) to find t, and then compute displacements from = − ( ). Typically, for a circle, there are two solutions obtained for t, and one must select the solution, which corresponds to the smallest . The y-displacement results are summarized in Table 1.

Validation of the TPIK Algorithm by Mathematical Simulation
In this section, a mathematical model was constructed to simulate the y-displacements between sensor coordinates and parametric circles. For simplicity, it was assumed that the sensors were initialized at the same level, i.e., y 1o = y 2o = y 3o = 0. The calibration master was represented with parametric equations (x(t), y(t)) that simulate two concentric circles of a known center and two given radii {R c2 , R c3 } as shown in Figure 5. To find the y-displacements between the sensor coordinates and the parametric equations, substitute the x-location in t = cos −1 ((a c − x)/R) to find t, and then compute displacements from y = b c − R sin(t). Typically, for a circle, there are two solutions obtained for t, and one must select the solution, which corresponds to the smallest y. The y-displacement results are summarized in Table 1.   (7), (9), (11), and (13). It is important to reemphasize here that the solution set is sufficient for TPFK algorithm, but the solution set is necessary to compute . The TPIK algorithm is comprised of non-linear equations, which were intentionally solved iteratively using optimization techniques. The 'fsolve' command in MATLAB [19] was utilized along with the 'trust-region-dogleg' algorithm and the 'optimplotfirstorderopt' option to plot first-order optimality at each iteration. This numerical method requires an initial estimation for a solution, S. The validity of the solution depends on the convergence criterion of the selected optimization method and the final values of the functions in Equations (7), (9), (11), and (13). The numerical iterations of the optimization method might prematurely stop producing multiple final values depending on the convergence criterion and initial estimation. Therefore, the convergence of the solution was studied and several initial estimations benchmarked relative to true solution that could be obtained from the parametric model. Table 2 is a summary of the studies conducted for several initial estimations. In the absence of information about the true solution, success was evaluated form the numerical solution using its convergence to a final solution and the residual that is defined from Equations (7), (9), (11), and (13) by: The Ers value could be used as an indicator to evaluate the difference between the numerical solution and the true solution if the solution exists and is unique. While the existence and uniqueness of a real solution for both the TPIK algorithm and the optimization method was not investigated, the researcher purposely limited this paper to the numerical comparison. The numerical values of in Table 2 are rounded to 3 decimal points and displayed for each trial when the optimization algorithm converges or stops. For example, trials 1, 4, 5, 8, and 9 in Table 2 produced acceptable results because Mathematical simulation of calibration procedures: the circles are constructed for (a c , b c ) = (10, 40) mm and two radii R = {38, 36} mm by using two parametric functions In practice, it is desired to recover the solution values S x 2o , x 3o , y 2o , y 3o , a c , b c by using the displacement information D y 12 , y 22 , y 32 , y 13 , y 23 , y 33 in the TPIK Equations (7), (9), (11), and (13). It is important to reemphasize here that the solution set S p is sufficient for TPFK algorithm, but the S solution set is necessary to compute S p . The TPIK algorithm is comprised of non-linear equations, which were intentionally solved iteratively using optimization techniques. The 'fsolve' command in MATLAB [19] was utilized along with the 'trust-region-dogleg' algorithm and the 'optimplotfirstorderopt' option to plot first-order optimality at each iteration. This numerical method requires an initial estimation for a solution, S. The validity of the solution depends on the convergence criterion of the selected optimization method and the final values of the functions in Equations (7), (9), (11), and (13). The numerical iterations of the optimization method might prematurely stop producing multiple final values depending on the convergence criterion and initial estimation. Therefore, the convergence of the solution was studied and several initial estimations benchmarked relative to true solution that could be obtained from the parametric model. Table 2 is a summary of the studies conducted for several initial estimations. In the absence of information about the true solution, success was evaluated form the numerical solution using its convergence to a final solution and the residual that is defined from Equations (7), (9), (11), and (13) by: The E rs value could be used as an indicator to evaluate the difference between the numerical solution and the true solution if the solution exists and is unique. While the existence and uniqueness of a real solution for both the TPIK algorithm and the optimization method was not investigated, the researcher purposely limited this paper to the numerical comparison. The numerical values of S in Table 2 are rounded to 3 decimal points and displayed for each trial when the optimization algorithm converges or stops. For example, trials 1, 4, 5, 8, and 9 in Table 2 produced acceptable results because they converged quickly with very small E rs values. While trials 6, 7, and 10 failed to produce acceptable results due to their non-convergence and very large E rs values. On the other hand, trials 2 and 3 did not converge, however they produced adequate final values because E rs were considerably small when compared to trials 6, 7, and 10. These studies also show that the numerical solution in trial 8, which has the fastest convergence and smallest E rs , had the least difference from the true solution.

Validation of the TPFK Algorithm by Mathematical Simulation
As previously discussed, the TPFK algorithm requires that system parameters S p x 2o , x 3o , y 2o , y 3o are previously known or computed from the TPIK algorithm. It is assumed that the displacement measurements of the inspected part are the same as the values in Table 2. The TPFK algorithm performance was compared by using two sets of system parameters: true solution and the computed solution. Table 3 shows that the computed {a, b, R 3p } values for both radii agree closely with the true radii values.

Validation of the TPIK Algorithm by CAD Simulation
A 3D assembly CAD model was constructed using SolidWorks software [20] and a dimensioning accuracy of ±0.01 mm. The model consists of three identical sensor probes with round tips, a fixture, a disc-part of unknown radius (inspected part), and a calibration master of three known radii as shown in Figure 6. The mating constraints are applied such that the probes can slide, and the calibration master can freely rotate. The tip of each probe is brought to contact at the largest radius (R c1 = 40 mm) in the calibration master. The probes location is fixed, and the calibration master is rotated about its center to measure displacements corresponding to each radius in the calibration master. The tip of sensor 1 is used as a reference in the Cartesian coordinate system to measure the displacements y 12 , y 22 , y 32 and { y 13 , y 23 , y 33 for R c2 and R c3 , respectively. The measurements are recorded in Table 4. The TPIK algorithm was used to compute (a c , b c ), (x 2o , y 2o ), (x 3o , y 3o ) for several initial estimations and the measurement resolution of the displacements results are shown in Table 5. A comparison of these values with the true dimensions extracted from the 2D-CAD drawing of the assembly is shown in Figure 7. In this simulation, the measurement resolution of the displacement, n, was achieved by rounding dimensions in Table 4.

10
of the assembly is shown in Figure 7. In this simulation, the measurement resolution of the displacement, n, was achieved by rounding dimensions in Table 4. Table 4. TPIK algorithm based on CAD simulation: displacements measured between sensor coordinates and the calibration master for two radii by using direct CAD dimensioning tools rounded to n = 6 decimal places.
Sensor  It is observed the convergence to a solution in trial 3 of Table 5 was unsuccessful despite the high measurement resolution used in the computation. This could be attributed to using initial estimations that are far-off from the true values or because the optimization algorithm might not perform adequately when it searches for non-zero values. Therefore, to improve the convergence, it is recommended to initialize the y-locations of the sensors 2 and 3 { , } to zero by replacing the large radius in the calibration master with a flat surface. However, trials 1 and 2 did not converge,    It is observed the convergence to a solution in trial 3 of Table 5 was unsuccessful despite the high measurement resolution used in the computation. This could be attributed to using initial estimations that are far-off from the true values or because the optimization algorithm might not perform adequately when it searches for non-zero values. Therefore, to improve the convergence, it is recommended to initialize the y-locations of the sensors 2 and 3 { , } to zero by replacing the large radius in the calibration master with a flat surface. However, trials 1 and 2 did not converge, It is observed the convergence to a solution in trial 3 of Table 5 was unsuccessful despite the high measurement resolution used in the computation. This could be attributed to using initial estimations that are far-off from the true values or because the optimization algorithm might not perform adequately when it searches for non-zero values. Therefore, to improve the convergence, it is recommended to initialize the y-locations of the sensors 2 and 3 y 2o , y 3o to zero by replacing the large radius in the calibration master with a flat surface. However, trials 1 and 2 did not converge, yet produced an accurate solution when compared to the true values. This might be due to a relatively small E rs value. The effect of the measurement resolution was studied using the initial estimations to set the true values. This was carried out for n = {6,5,4,2} in trials 4-8, which confirms that high measurement resolution was not always necessary to attain an acceptable solution. Similarly, a low-resolution measurement may successfully converge to a solution, but is far off from the true values. For example, trial 5 was rounded to n = 5 decimal places, and the optimization had successfully converged to the true values with a small E rs . The simulation was repeated using the same initial estimation with n = {4,3,2}. Interestingly, trial 7 was simulated for n = 3 (trial 5) and had not only converged quicker than that of n = 5 (trial 3), but also had a smaller E rs value. However, the solution for trial 7 was farther off from the true value when evaluated against trial 5.

Validation of the TPFK Algorithm by CAD Simulation
It was presented in previous sections that the TPFK algorithm computed the radius of the part independent of how it was positioned. In this study, the true location of the sensors was obtained {x 2o = 10.000, x 3o = 20.000, y 2o = −1.175, y 3o = −0.014} mm and the true inspected part geometry by using CAD dimensioning tools in Figure 6 where all values were approximated to n = 3 decimal points. The first group of simulations was conducted for inspected parts whose true radii R e ranged between 38.000 mm to 39.500 mm and with all of their true centers fixed at (10.010, 38.788) mm. The TPFK algorithm was implemented to compute the "fit" of their radii and centers based on the displacement measurement and the sensor location.
The absolute error between the true and computed radius |R − R e | was calculated and the results are summarized in Table 6. The error increased as the radius of the inspected part radius approached 40 mm. Specifically, when the y-displacements measurements approached zero, the TPFK became sensitive to rounding. This could be attributed to sensors that were initialized for a calibration master whose radius was at 40 mm. The second group of study was conducted for an inspected part of a fixed radius of 38 mm, however, the inspected part was introduced to the sensors from several arbitrary locations. It was determined that the computed radius of the inspected part was in close proximity with the true value for the arbitrary locations. Additionally, the |R − R e | error followed the same findings as in the first study group where the error generally increased when the displacement measurement vector decreased.

Study of the Roundness Index Based on the Mathematical Model
The fundamental aim of this section was to determine how well the proposed TPIK and TPFK algorithms work in computing the roundness index β of the inspected part. It was assumed that the inspected part had perturbation around the circler profile. Additionally, it was desired to observe the effect of the sampling pattern on the roundness index. Furthermore, it was assumed that the inspected part is fixed at the center C e (a, b) and was able to rotate in increments of δθ with a total rotational angle θ measured from a given reference. The perturbed circle was modeled by using the nominal radius R e and a sinusoidal function as shown in Figure 8. The actual radius value of a point at an angle of ϕ can be represented by a piecewise smooth curve in reference to its polar coordinate while it is not in rotation (θ = 0): where R p is the amplitude of the fluctuation of a sinusoidal function, and ω(θ) is the frequency function. The first step is to estimate the number of sample points required for scanning a fluctuating profile. Theoretically, there are infinite points along the length of the inspected part profile. The total arc length L of R(ϕ) is computed over the entire interval 2π: where is the rotational angle 360 ≥ ≥ 0 with increment. Figure 9a simulates several fluctuating circles by using a non-zero frequency and fluctuation amplitude. The sensor displacements , , and shown in Figure 9b represent the intersection between the vertical line passing between each sensor and the inspected part. It is desired to observe how the roundness index varies with the frequency, fluctuation amplitude , nominal radius , and measurement resolution n. Table 7 shows an example of the simulation carried out for = 20 , frequency = 5 rad, = 38/40 mm, = 38 mm, a = 10 mm, b = 40 mm, and n = 3. The sensor locations shown in Figure 8b were set to x = 0, x = 10.000 mm, x = 20.000 mm, y = 0, y = 0 mm, and y = 0. The simulation was run as follows: The initial step was to solve equation 18 Figure 8b where { ( ), ( ), ( )} represent the displacement distances between sensors and the part surface, i.e., { , , }. The second step was to apply the TPFK algorithm to calculate the corresponding part radius and its center from these three displacement vectors. The previous steps were repeated for every ; the results are shown in Table  7. The average center of the inspected part was , = (9.9994, 39.9821) mm, and the equivalent average radius was = 37.9992 mm. These values were obtained for the inspected part whose center was (10,40) mm, nominal radius = 38 mm, and fluctuation amplitude = 38/40 mm. Therefore, stochastically, the average radius could be used for describing the radius of the inspected part, and the fluctuation of , values could be used as an indication of roundness. Comprehensive studies of the relationship between roundness and the fluctuation ratio / were simulated for several conditions in Figure 10. For example, Figure 10a shows a simulation of the roundness index for three sampling increments ={0.5°,10°,60°} and = 50 rad. It is evident that the computation of the roundness index was sensitive to the sampling increment when the fluctuation amplitude was large, but when the fluctuation was small, all roundness indices converged to the same final value regardless of the value of . Let the desired linear increment ∆L between any two data points on the surface of the inspected part be used to sample the entire length L. This generates approximately ∼ L/∆L sample points over the entire length. Suppose a circle (is utilized) whose circumference is equal to the length L, or has a radius of R L = L/2π, then the average angular increment that uniformly dissects this circle is δθ * ∆L/R L radian, or δθ * 2π∆L/L radian. The significance of δθ * is to provide the programmer a reference value for selecting δθ. Likewise, the reference value can be used to calculate the spacing between two points when a uniform rotational increment sampling is imposed.
The function is parameterized to simplify the analysis in Equation (16) as follows: where θ is the rotational angle 360 o ≥ θ ≥ 0 with δθ increment. Figure 9a simulates several fluctuating circles by using a non-zero frequency and fluctuation amplitude. The sensor displacements y 1 , y 2 , and y 3 shown in Figure 9b represent the intersection between the vertical line passing between each sensor and the inspected part. It is desired to observe how the roundness index varies with the frequency, fluctuation amplitude R p , nominal radius R e , and measurement resolution n. Table 7 shows an example of the simulation carried out for δθ = 20 o , frequency ω = 5 rad, R p = 38/40 mm, R e = 38 mm, a = 10 mm, b = 40 mm, and n = 3. The sensor locations shown in Figure 9b were set to x 1o = 0, x 2o = 10.000 mm, x 3o = 20.000 mm, y 1o = 0, y 2o = 0 mm, and y 3o = 0. The simulation was run as follows: The initial step was to solve Equation (18) at every angle θ to find the solution set t for which x(t 1 ) = 0, x(t 2 ) = x 2o = 10 mm, and (t 3 ) = x 3o = 20 mm. The solution requires a numerical search algorithm such that y(t) < b ∀t ∈ {t 1 , t 2 , t 3 }. Thus, the solution lies on the lower half of the part shown in Figure 9b where y(t 1 ), y(t 2 ), y(t 3 ) represent the displacement distances between sensors and the part surface, i.e., y 1 , y 2 , y 3 . The second step was to apply the TPFK algorithm to calculate the corresponding part radius and its center from these three displacement vectors. The previous steps were repeated for every θ; the results are shown in Table 7. The average center of the inspected part was a, b = (9.9994, 39.9821) mm, and the equivalent average radius was R a = 37.9992 mm. These values were obtained for the inspected part whose center was (10,40) mm, nominal radius R e = 38 mm, and fluctuation amplitude R p = 38/40 mm. Therefore, stochastically, the average radius R a could be used for describing the radius of the inspected part, and the fluctuation of R 3p,i values could be used as an indication of roundness. Comprehensive studies of the relationship between roundness β and the fluctuation ratio R e /R p were simulated for several conditions in Figure 10. For example, Figure 10a shows a simulation of the roundness index for three sampling increments δθ ={0.5 • ,10 • ,60 • } and ω = 50 rad. It is evident that the computation of the roundness index was sensitive to the sampling increment when the fluctuation amplitude was large, but when the fluctuation was small, all roundness indices converged to the same final value regardless of the value of δθ.  Table 7. Calculation of Roundness for R p = 38/40 mm, R e = 38 mm, and frequency ω = 5 rad. The system parameters were set to x 1o = 0, x 2o = 10.000 mm, x 3o = 20.000 mm, y 1o = 0, y 2o = 0 mm, and y 3o = 0.
True Part Center is Fixed and Rotates about a = 10 mm and b = 20 mm by Sampling Increment of δθ = 20 • , and Frequency ω = 5 rad.
Probes Measurement (mm), n = 3 Calculated (mm) This implies that in such cases, fewer sampling points could be used to calculate the roundness of the inspected part. Under the same conditions as applied to Figure 10a, the effect of the fluctuation ratio R e /R p to the ratio between average radius to the exact radius R a /R e was observed with results shown in Figure 10b. It was determined that the R a converged to R e as the fluctuation amplitude became small, however, this convergence was rapid when the frequency of the fluctuation decreased as noted by comparing Figure 10b,d.
In addition, a comparison between Figure 10a,c show that the computation of the roundness index was less sensitive to δθ when ω was small. Finally, the relation between the roundness index and the measurement resolution of the displacement was studied for n = {6,5,4,3,2} and were all found to be identical. This indicates that it could be possible to evaluate the roundness of the inspected part using a less expensive sensor with lower resolution.
Since the fluctuation had a symmetric pattern, the computation of average center {a = 9.9994, b = 39.9821} mm was close to the exact value {a = 10, b = 40} mm. Nevertheless, this did not imply that the difference between the average and exact values was predestined to worsen when the rotation increment increased. However, such a difference would also depend on the effect of the scanning resolution and the sample point distribution as will be presented in the subsequent section.

Study of Roundness Index Based on CAD Model
The "equation driven curve" function in SolidWorks was used to plot the parametric Equation (18). Figure 11a shows a solid model assembly comprised of three probes and a fluctuating part. To simplify the measurement in CAD software, the researcher simulated the fluctuation profile using Equation (18) in 2D drawing. Figure 11b shows an example of inspected part with R e /R p = 100, and ω = 20 rad/s. The center of the inspected part rotated freely. The displacements { y 1 , y 2 , y 3 were measured at every ∆θ = 2 o by using the dimensioning tools in 2D drawing. Subsequently, the {a, b, R 3p } were measured at every increment using the TPFK algorithm with results in Table 8.
Note that for increment ∆θ = 2 • , the measurements values repeated after every 16 • , therefore the average values of {a, b, R a } within the first J = 9 rotations would produce similar values if calculated from sample size M × J rotations, where M was any real integer. Therefore, one could reduce the cost of data collection when the surface pattern of the inspected part is known to have a uniform repetition. The model in Figure 11 produced an average radius R a = 45.1944 mm and an average center of (a, b) = (7.7273, 46.1588) mm. Although the sampling increment was small ∆θ = 2, the average values differed from the exact value (a, b) = (7.873513, 43.48817) mm, and were also biased. Unlike the case study in Table 7 (ω = 5 rad/s), which produced smaller differences despite the large sampling increment used (∆θ = 20). This could suggest that a smaller sampling increment be used when the fluctuation frequency is high.
16 average values differed from the exact value (a, b) = (7.873513, 43.48817) mm, and were also biased. Unlike the case study in Table 7 ( = 5 rad/s), which produced smaller differences despite the large sampling increment used (∆ = 20). This could suggest that a smaller sampling increment be used when the fluctuation frequency is high.

Validation of the TPIK and TPFK Algorithms Based on the Empirical Study
A radial load bearing was mounted on a customized 3D printed fixture as shown in Figure 12 Three micrometers {S1, S2, S3} were used to measure the displacements from the initialized. coordinates into the surface of the bearing. Researchers applied adhesive tape of known thickness on the bearing surface to create multiple radii for both calibration and inspected part experiments.
The calibration procedures require two known radii sharing the same center, therefore, the base radius of the bearing = 40.00 mm was utilized to measure the first group of displacements, and also 2.13 mm thick tape was applied on another section of the bearing surface to create the second radius = 42.13 mm and obtain the second group of displacements. The average and standard

Validation of the TPIK and TPFK Algorithms Based on the Empirical Study
A radial load bearing was mounted on a customized 3D printed fixture as shown in Figure 12 Three micrometers {S 1 , S 2 , S 3 } were used to measure the displacements from the initialized. coordinates into the surface of the bearing. Researchers applied adhesive tape of known thickness on the bearing surface to create multiple radii for both calibration and inspected part experiments.
deviation sd measurements are summarized in Table 9. Likewise, following the calibration procedures an inspected part with radius of = 41.03 mm was created by adding a 1.03 mm thick tape on a section on the bearing. Following this, the corresponding displacements were measured and results are shown in Table 10. A summary of all measurements is shown in the error plot in Figure 13.  = .

Sensor
Avg./SD (mm) Measurement (mm) S1 ( , sd1) (13.10, 0.03) S2 ( , sd2) (2.85, 0.01) S3 ( , sd3) (13.69,0.11) S1 S2 S3 Sensor Number The calibration procedures require two known radii sharing the same center, therefore, the base radius of the bearing R c2 = 40.00 mm was utilized to measure the first group of displacements, and also 2.13 mm thick tape was applied on another section of the bearing surface to create the second radius R c3 = 42.13 mm and obtain the second group of displacements. The average y and standard deviation s d measurements are summarized in Table 9. Likewise, following the calibration procedures an inspected part with radius of R e = 41.03 mm was created by adding a 1.03 mm thick tape on a section on the bearing. Following this, the corresponding displacements were measured and results are shown in Table 10. A summary of all measurements is shown in the error plot in Figure 13. The researcher used the displacement values of the calibration master in Table 9 to recover the system parameters S x 2o , x 3o , y 2o , y 3o , a c , b c , by using the TPIK algorithm. While the measurement of S was difficult to obtain and inaccurate in practice, it roughly measured the sensors location and the center of the master calibration part S{26,53,1,1,26,45} mm, and then these values were applied as the initial estimation in the algorithm to calculate S based on data in Table 9. The final values obtained from the TPIK algorithm after 106 iterations was S{25.30, 53.32, 1.12, 1.04, 25.71, 45.02} mm with a total numerical error E rs = 0.0071 mm. The final values were applied in the TPFK algorithm to calculate the radius of the inspected part using the data in Table 10. The estimated {a, b, R 3p } values were {25.72, 44.75, 40.78} mm, respectively. This indicates that the error in measuring the inspected part radius was er = R 3p − R e = 0.25 mm or 0.61% percentage error, which represents acceptable agreement between values given the simplicity of the experiment setup used in this study.

Conclusions and Future Work
This paper developed practical methods to measure radius and roundness of a circular arc from three distributed gauges. The new proposed TPIK algorithm provides a calibration mechanism that resets the measurement error resulting from machine and human errors by computing the location of distributed gauges in reference to a global Cartesian coordinate. The TPIK algorithm is a non-stochastic fitting method but requires numerical optimization techniques to determine a solution. The solution was investigated based on iterative optimization techniques that depend on an initial estimation. While this paper was not concerned with the approach used in finding a solution, however, it was determined that the final value solution was sensitive to the initial estimation. The existence and uniqueness of solution have not yet been investigated mathematically, but it was verified geometrically, that there must be at least one solution to satisfy the six non-linear TPIK equations. In addition, a TPFK algorithm was offered, which is a closed form solution to compute the geometry of the circular arc from three distributed gauges. This approach is a modification of the direct algebraic method that uses three points to fit a circle, however, the algorithm itself contains estimated parameters from the TPIK algorithm. Due to simplicity, accuracy, and computational efficiency of the TPFK algorithm, it was not possible to compute the radius and the center of a part being inspected, but the method was to compute the roundness of a fluctuating surface. The proposed roundness index was based on aforementioned algorithms and demonstrated the capability to present the degree of circularity of a part. It was determined that a successful implementation of roundness index depends on the selection of a proper scanning rate for the waviness attributes of the part being inspected, which includes fluctuation amplitude and frequency. The degree of accuracy of measurement did not play a major role when the roundness was computed, however, the sample point distributions of the fluctuation had an additional impact. Future work will investigate numerical methods for solving the TPIK algorithms and the corresponding performance tradeoffs.