{'fraction_non_alphanumeric': 0.03575825799416048, 'fraction_numerical': 0.01458120800647272, 'mean_word_length': 4.9558977582233394, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New Confocal Hyperbola-based Ellipse Fitting with Applications to Estimating Parameters of Mechanical Pipes from Point Clouds

This manuscript presents a new method for fitting ellipses to two-dimensional data using the confocal hyperbola approximation to the geometric distance of points to ellipses. The proposed method was evaluated and compared to established methods on simulated and real-world datasets. First, it was revealed that the confocal hyperbola distance considerably outperforms other distance approximations such as algebraic and Sampson. Next, the proposed ellipse fitting method was compared with five reliable and established methods proposed by Halir, Taubin, Kanatani, Ahn and Szpak. The performance of each method as a function of rotation, aspect ratio, noise, and arc-length were examined. It was observed that the proposed ellipse fitting method achieved almost identical results (and in some cases better) than the gold standard geometric method of Ahn and outperformed the remaining methods in all simulation experiments. Finally, the proposed method outperformed the considered ellipse fitting methods in estimating the geometric parameters of cylindrical mechanical pipes from point clouds. The results of the experiments show that the confocal hyperbola is an excellent approximation to the true geometric distance and produces reliable and accurate ellipse fitting in practical settings.


Introduction: applications of ellipse fitting
Points following elliptic patterns commonly and naturally occur in many real-world datasets. In digital images, the edge points of spheres and circles form ellipses due to projective transformation. In threedimensional point clouds, the intersection of spheres, cones, and cylinders with a plane also delivers points following an elliptic pattern. Reliable and accurate ellipse fitting to data is, hence, commonly an important step towards solving various real-world applications pertaining to shapes such as circles, cylinders and spheres, including, galaxy classification [1], traffic sign detection [2], terrestrial laser scanner (TLS) registration [3], optical instrument calibration [4], geometric as-built building information modeling (BIM) of pipes [5], cylinder fitting of industrial shapes [6], and cylindrical point cloud classification [7] 1 .
To fit an ellipse to a set of points, the ellipse parameters minimizing the sum of the squared distances between the points and an ellipse are estimated (least squares method). Ideally, it is desirable to minimize the function representing the orthogonal distances of points to an ellipse (commonly referred to as the geometric distance). Finding the orthogonal contact point to an ellipse is, however, not an easy task, and most available analytical methods are either numerically unstable, or iterative without guarantee of convergence [8], and consequentially time consuming [9]. The numerically stable methods, such as that proposed by Chernov [8] are not analytically closed-form, and hence, cannot be used as a minimization cost function. As a result, researchers have sought to minimize alternative and less complicated distance functions, such as algebraic [10][11][12][13], and Sampson [9,[14][15][16][17]. The methods minimizing algebraic and Sampson distances, however, share some undesirable properties. Algebraic methods are biased at locations of high curvature (see Section 2). The Sampson distance is not continuous in ℛ 2 (see Sections 2 and 5.1.1 Figure 4) and the methods are only suited for data with "moderate" noise levels [9]. The methods minimizing the algebraic and Sampson distances also require either additional constraints [12] or barrier terms [9] during the minimization of the objective function to guarantee an ellipse. Furthermore, Sampson and algebraic distances do not effectively predict the geometric distance and, as will be revealed in our experiments (see Section 5.1.1 and Figure 4), diverge from the true geometric distance as the point moves farther from the ellipse. This is most likely the reason they are only suited for "moderate" noise levels.
To overcome the aforementioned limitations associated with current ellipse fitting methods, this manuscript presents a new ellipse fitting method that minimizes the sum of the squared confocal hyperbola distances of points to an ellipse. Even though the confocal hyperbola distance approximation was introduced in Rosin [18] with promising geometric properties, its effectiveness as an ellipse fitting method has not yet been examined in the literature. The goals of this manuscript are to: (i) provide a closed-form solution to the confocal hyperbola distance function; (ii) introduce a new algorithm to minimize the sum of the squared confocal hyperbola distances; and (iii) thoroughly investigate its effectiveness in fitting ellipses to twodimensional (2D) data, acquired from simulated image edge points as well as point clouds of cylindrical mechanical pipes.
The remainder of the manuscript is structured as follows. In Section 2, an introduction to current ellipse fitting methods, along with their advantages and limitations is provided. Section 3 includes the detailed explanation of the proposed ellipse fitting algorithm. Section 4 outlines the metrics for validation of the proposed method along with the design and configurations of the various simulation-based and real-world experiments, carried out to assess the applicability of the proposed method compared to state-of-the-art ellipse fitting methods. Sections 5 and 6 detail the results obtained by the simulated, and real-world data in the designed configurations presented in Section 4, respectively. Section 7 discusses the summary of the major findings of this study.

Background: introduction to ellipse fitting
The equation of an ellipse in general form is represented as: where ( , ) are the point coordinates, = ( , , , , ) is the geometric parameter vector of an ellipse, which includes the coordinates of the center, semi major length, semi minor length and rotation angle of the ellipse, respectively. Equation (1) can be reduced to general algebraic form as proposed by Paton [10], which is commonly referred to as the algebraic distance: where ( , ) is the algebraic distance of point ( , ) to the conic, = ( , , , , , ) is the vector of algebraic parameters, = ( 2 , , 2 , , , 1) is the 6-dimensional coordinate set (commonly referred to as the design vector, or for more than one point, design matrix), and (. ) denotes the transpose function. To find a best fit conic to data points, the parameters providing the minimum sum of the squared algebraic distances are required, i.e.: where = is referred to as the scatter matrix. The problem imposed by equation 3 is a linear least squares problem subject to some parameter constraint to prevent the trivial solution, = 0. For instance, Paton [10] proposed to minimize equation 3, subject to ‖ ‖ = 1, where ‖. ‖ denotes the L2-norm. This problem is now reduced to fitting a hyper plane to a six-dimensional point coordinate set = ( 2 , , 2 , , , 1), which is the eigenvector corresponding to the smallest eigenvalue of the 6x6 scatter matrix, . Since this algebraic minimization ignores the correlations between the algebraic parameters 2 , it cannot guarantee an ellipse. In fact, the algebraic parameters must satisfy both of the following conditions to produce an ellipse: where det (. ) denotes the determinant function. To guarantee an ellipse fit, Fitzgibbon [12] proposed to solve equation 3 subject to ∆= −1. By introducing a constraint matrix ( ) for this condition, the problem was reduced to a generalized eigenvalue problem ( = ). Halir [13] observed that the constraint matrix of Fitzgibbon can provide numerically unstable results, especially when the points are exactly on the ellipse. Hence, they proposed a numerically stable method to solve Fitzgibbon's constraint.
To provide some perspective on the meaning of the algebraic distance, Bookstein [11] showed that the algebraic distance is proportional to: where is the distance from the point to the center, and and are shown in Figure 1a. The formulation of Bookstein shows that for the same across all data points (e.g. same noise level), the algebraic distance of the point to the curve is smaller when is larger (i.e. around the major axis). In other words, at locations of high curvature, the algebraic distance under-estimates the distance of the point to the curve. This is referred to as a high curvature bias problem imposed by equation 3 [19]. To solve the curvature bias problem, Agin [14] proposed the minimization of the following distance function: where ∇ denotes the gradient function, and ( , ) is commonly referred to as the Sampson distance of points ( , ) to the ellipse (Sampson 1982). Minimization of the sum of the squared Sampson distances is, however, a non-linear problem. As a result, Agin [14] proposed to approximate equation 6 subject to ∑ �∇ ( , )� =1 2 = 1, which reduces to a linear generalized eigenvalue problem. Later, Taubin [16] independently derived this exact formulation, which is now commonly referred to as the Taubin's method in the literature. Instead of solving the actual non-linear minimization problem of equation 6, Sampson [15] proposed to adopt a reweighting strategy (referred to as the gradient weighted method) using the distance of equation 6 to iteratively improve Bookstein's algebraic conic fitting method. The approximate solution to equation 6 proposed by Taubin and Agin provides a curvature bias correction of the algebraic distance up to the first order. Kanatani [20] went one step further and corrected the algebraic distance approximation up to the second order of the Taylor series expansion of the algebraic distance. Their method also reduces to a generalized eigenvalue minimization problem similar to that of Taubin and Agin. The approximate solutions to equation 6 using Agin's/Taubin's and Kanatani's methods, however, do not necessarily satisfy the conditions of equation 4, and hence, cannot guarantee an ellipse. To this end, Szpak [9] proposed to incorporate a barrier term in the minimization of the Sampson distance, which helps limit the solution space to guarantee an ellipse. Their formulation, however, cannot be reduced to a linear least squares problem, and hence it was solved using the Levenberg-Marquardt algorithm. The Sampson distance, although providing a good approximation to reduce the impact of curvature biases imposed by the algebraic distance function, is undetermined when the point coordinate is on the center of the ellipse (i.e. equation 6 goes to infinity as the point approaches the center of the ellipse). In addition, it bears no direct relationship to the true geometric distance, especially as the point moves farther from the ellipse 3 . Harker and O'Leary [21] pointed out that the Sampson distance is merely the geometric distance of a point to the first order approximation of the algebraic distance function and should not be mistaken with the first order approximation of the geometric distance from the point to the conic. They provided a closed-form solution to the first order approximation of the geometric distance, which was, in fact, also found to be a weighted algebraic distance. However, as pointed out by Uteshev and Goncharova [22], Harker and O'Leary's solution can display singularities in elongated ellipses for points around the major axis. In addition, as will be revealed later in Section 5.1.1, Figure 4, the distance obtained by Haker and O'Leary is undetermined (approaches infinity) for points at the focal points of an ellipse.
The ideal geometric minimization cost function is to minimize the sum of the squared distances represented by the orthogonal projections of the points onto the ellipse (i.e. the geometric distance; Figure 1a; red solid line). Achieving a numerically stable, closed-form analytical solution to the geometric distance is, however, not an easy task [8]. For instance, Safaee-Rad [23] proposed an exact solution to the geometric distance that reduces to finding the roots of a quartic equation. However, analytical methods for finding the roots of a quartic equation are known to be numerically unstable. To this end, Ahn [24] proposed to divide the solution provided by Safaee-Rad into two functions that were iteratively solved using the generalized Gauss-Newton method. However, as pointed by Chernov [8], the method does not guarantee convergence and numerically failed about 3% of the time. Currently, the most accurate and numerically stable process for finding the orthogonal projection of points onto conics is proposed by Chernov [8]. The method converts the orthogonal projection problem to finding one real generalized eigenvalue between the ellipse and its auxiliary conic (which turns out to be a hyperbola; see Figure 1 of [8]). Using one real eigenvalue root, a degenerative conic, comprised of two lines, is produced. The closest orthogonal projection of a point to the ellipse is the point of intersection between the two lines and the ellipse (up to four points), whose distance to the point of interest is the least. However, this elegant method is still not analytically closed-form and, hence, cannot be used as a minimization cost function to estimate the parameters of the best fit ellipse.
Another approximation to the geometric distance is the confocal hyperbola method used in Rosin [18], see Figure 1b. This distance approximation relies on the fact that confocal ellipses and hyperbolas are mutually orthogonal at the intersecting point. Since the confocal hyperbola passing through a given point is nearly linear at the point of intersection, it is a good approximation to the geometric distance of a point to the ellipse. It is worth mentioning that the auxiliary conic to an ellipse at a specific point, as shown by Chernov [8], is also a hyperbola. Hence, the confocal hyperbola distance is speculated to possess similar geometric properties to the true geometric distance of a point to the ellipse. In fact, Rosin had compared 16 different geometric distance approximation methods, including algebraic, Sampson, Harker and O'Leary, and the confocal hyperbola, in terms of properties such as linearity, curvature bias, asymmetry and general goodness [18,25,26]. Overall, the confocal hyperbola was the best in every category, especially for larger noise levels (see Table 2 of [18,25,26]). In Section 5.1, we will demonstrate that the confocal hyperbola method provides almost identical results to that of Chernov [8], which is the current state-of the-art. Although Rosin provided a process to derive the confocal hyperbola distance, it was not presented in general closed-form, imposed singularities for circles, and required the examination of four solutions. As a result, up until now, the confocal hyperbola distance had not been used to fit ellipses to data. In Section 3, we will expand on the formulation provided by Rosin and present an analytically closed-form solution to this geometric distance approximation, which can then be minimized to obtain the best fit ellipse parameters.

Ellipse fitting using confocal hyperbola distance
Various ellipse fitting methods along with their limitations were introduced in Section 2. Here, we will formulate a new ellipse fitting method based on the minimization of the confocal hyperbola distance. To produce a minimization objective function, we shall first derive a closed-form solution to the confocal hyperbola distance by determining the point of intersection between the ellipse and the confocal hyperbola passing through a given point ( , ). Using a translation and rotation, equation 1 can be rewritten as follows: where and are the semi-major and semi-minor lengths of the ellipse, respectively. In this new coordinate system, the problem is now reduced to finding the point of intersection between the ellipse of equation (7) and the confocal hyperbola passing through point ( , ) with equation: where ℎ and ℎ are the semi-major and semi-minor lengths of the hyperbola, respectively. Since the ellipse and hyperbola are confocal, the following relationships are valid: where is the distance of one of the foci to the center. By substituting equation 9 into 8, the following quadratic function is derived: Solving equation 10 for provides two solutions; however, only the following can ensure both and 2 − are positive (this is of course a necessary condition to produce real values for ℎ and ℎ ): Substituting equation (11) into (8) and solving simultaneously with (7) will provide the following simplified solutions to the orthogonal intersecting point ( , ): Since four solutions exists (± , ± ), Rosin [18] proposed to calculate the distances between the point and the four solutions of equation 12 and select the solution with the minimum distance 4 . However, it can be observed that the minimum distance must belong to the solution that is in the same quadrant as the point. Since , ℎ , , ℎ and are all positive, the correct solution can be formulated by simply restricting the point coordinate ( , ) to the first quadrant using an absolute value function. The confocal hyperbola distance in this context can now be defined in closed form as follows: � is a 2D vector of distances in the translated and rotated − coordinate plane (see equation 7), and ‖ . ‖ denotes the L2-norm. In its current form, equation 12 is numerically singular when = (i.e. = 0, when the data follow a circular pattern). However, if the formulations of ℎ and ℎ in equation 11 are rearranged, we will have: therefore, the closed formula for the orthogonal intersecting point of equation 12 will now become: This new formulation shows that when = 0 (i.e. = ), the orthogonal contacting point exists, and hence, does not impose any numerical singularities for circles. Moreover, if equation 15 at = 0 and radius, = = is employed, equation 13 simplifies to: ℎ ( , ) = √ 2 + 2 − , which is the true geometric distance of a point to a circle. Another interesting observation from equation 15 is that when = 0 (points on the minor axis), or when = 0 and | | ≥ (points on the major axis beyond the focal points), even though a confocal hyperbola does not technically exist, the equations of the intersecting point are still deterministic and the distance functions reduce to the true geometric distance through the following simplified formulas: Considering the aforementioned formulations, equations 13, 15 and 16 can now be employed to find the geometric parameter vector = ( , , , , ) by minimizing: where equation 17 is a non-linear least squares problem and must be solved using an iterative heuristic method such as Levenberg-Marquardt. To this end, Algorithm 1: Ellipse Fitting using Confocal Hyperbola is formulated as follows: 1-Find an initial estimate of the geometric parameters = ( , , , , ) of the ellipse using one of the methods presented in Section 2. Here, we will use Halir's numerically stable method, which is fast and guarantees an ellipse, as the initial parameter approximation. It is worth mentioning that there are many variations of the Levenberg-Marquardt algorithm, some include strategies such as parameter-wise re-weighting of the damping parameter [27] as well as geodetic acceleration [28], which may also be used to solve equation 17. The proposed implementation of the Levenberg-Marquardt algorithm was, however, sufficient to achieve the results presented in Sections 5 and 6 in all cases. Another important note is that equation 17 can also be used to derive the best fit hyperbola to a given set of points using Algorithm 1 subject to the parameter vector ℎ = ( , , ℎ , ℎ , ). Fitting a hyperbola to a set of points will, however, not be considered in this manuscript. The closed formulation for the Jacobian, ℎ , of the distance function is provided in Appendix I.

Experimental Design
In this manuscript, two main categories of experiments are designed, namely, simulation, and real-world. Simulation-based experiments were designed to validate the results of the ellipse fitting compared to other state-of-the-art and established methods in data resembling edge points in images. The real-world experiments were designed to evaluate the ability of different ellipse fitting methods in recovering the geometric parameters of cylinders in 3-dimensional point clouds. The two classes are explained in more detail in the following 5 . The summary description of the designed experiments is provided in Table 1.

Metrics for validation of results
For the simulation-based experiments, two metrics are used: (i) the root mean squared error (RMSE) of the best fit ellipse; and (ii) the L2-norm of the estimated best fit parameter vector from the ground truth divided by the L2-norm of the ground truth parameter vector (abbreviated here as P-Error and reported in percentages). The RMSE of the best fit will be calculated using the method of Chernov [8], which is currently the most reliable method for projecting points onto ellipses, and consequentially determining the actual geometric distance from points to the ellipse. For the real-world experiments, the L2-norm of the estimated cylinder parameters from the ground truth parameters, namely center, radius, and axis angle, are used. The ground truth geometric parameters of the cylinders from point clouds are determined using the robust cylinder extraction method proposed by Maalek [5].

Simulation-based experiments
In this section, the method in which elliptic points are generated along with the various simulated configurations used to evaluate different ellipse fitting methods are explained.

Simulating elliptic points
An important consideration for the simulated experiments is the method in which simulated ellipses are generated. Here, our focus for simulation experiments is producing elliptic points resembling edge points in images (three-dimensional (3D) point clouds will be considered in the real-world experiments). The edge points in images contain two important characteristics. First, the larger the ellipse (or any shape for that matter), the higher the number of ellipse edge points will be. In other words, an un-occluded circle with radius of 3 pixels cannot comprise of say 100-pixel points. Second, the discrete pixelated nature of digital images ensures that an extracted edge point will always fall on an image pixel. In other words, even if no random measurement error exists, the edge points representing an ellipse in an image will fall on the closest pixel and not necessarily exactly on the ellipse. The effects of digitization along with those induced by the number of points with respect to the scale of an ellipse must be considered when simulating ellipses. A process is, hence, formulated to automatically determine the pixel points representing the ellipse in images using only the geometric parameters of the ellipse. To this end, Algorithm 2: Simulating Elliptical Edge Points, given the geometric parameter vector = ( , , , , ) , and the polar arc angles for the beginning, , and the end, , of the ellipse, subjected to random measurement error, (0, 2 ), is developed as follows: 1-Find the orthogonal circumscribed rectangle to the ellipse in digital coordinates using the following: Where ( , ) and ( , ) are the minimum and maximum of the limits of the orthogonal circumscribed rectangle in and , respectively, and ⌊. ⌋ and ⌈. ⌉ denote the floor and ceiling functions, respectively. 2-Generate a grid of = − + 1 by = − + 1 points between the rectangular limits of ( , ) and ( , ) with one pixel spacing in the and direction, respectively. 3-Determine the orthogonal contacting point ( , ) of each pixel point ( , ) within the generated grid from the ellipse using the method of Chernov [8]. 4-Find the pixels satisfying the following conditions: 5-Find the orthogonal contacting points ( , ) of step 4 whose polar angles are between and .
6-Subject the ( , ) coordinates satisfying steps 4 and 5 to additive random measurement error (0, 2 ). is in pixel coordinates, and and coordinates of the points are considered independent identically distributed (i.i.d). 7-Find the closest pixel grid to the coordinate points of step 6. Figure 2a and 2b illustrate the simulated pixel coordinates of an ellipse with the geometric parameter vector of (1,3,15,10, 30°) with no random measurement error, and subjected to random measurement error of (0,1), respectively.

Ellipse simulation configurations
Predicated on the proposed methods, two main categories of simulation-based experiments are carried out, namely, point to ellipse distance measure, and ellipse fitting evaluation (further divided into two subcategories), as follows: 1-Point to ellipse distance measure: this experiment is designed to numerically identify the best method amongst various point to ellipse distance measures that reliably and accurately predicts the ground truth geometric distance, obtained by Chernov [8]. To this end, 10,000 ellipses were randomly generated using Algorithm 2, subject to the following conditions: a) rotation angle, randomly chosen between 0 and ; b) aspect ratio, randomly chosen between 1.5 and 4; c) measurement standard deviation, randomly chosen between 0.5 and 15 pixels; d) semi-minor, was chosen as 50 pixels; and e) ellipse spanning arc ( − ), was chosen as 2 (i.e. quarter ellipse).
The ellipse spanning arc was selected as a quarter of an ellipse (rather than say a full ellipse) since the point to ellipse measurements are not expected to drastically change between different quadrants of the same simulated ellipse. For each measurement, the distance of point to ellipse is calculated using the confocal hyperbola, algebraic, Sampson (see Section 2), and Harker and O'Leary [21] . The absolute deviation of the distances for each method from that obtained by Chernov is reported.
2-Overall ellipse fitting evaluation: this experiment is designed to evaluate the effectiveness of the proposed ellipse fitting method using confocal hyperbola (Section 3) compared to the ellipse fitting methods of Ahn [24], Szpak [9], Kanatani [20], Taubin [16], and Halir [13]. Similar to experiment 1, 10,000 randomly generated ellipses were used to determine which ellipse fitting method produces the fit with the lowest RMSE along with P-Error (see Section 4.1). The ellipse parameter configurations were as follows: a) rotation angle, randomly chosen between 0 and ; b) aspect ratio, randomly chosen between 1.5 and 4;
For each randomly-selected configuration, 500 different simulations were carried out to capture the behavior of random measurement errors. The results of the mean, median and 95 th percentile for the 500 simulations in each configuration, achieved by each method, were reported. In addition, the methods of Ahn, Szpak and ours are iterative, whereas the rest are direct. To provide a fair comparison between the different iterative methods, similar iteration-termination criteria are used. For all methods, Halir's ellipse fitting is used to estimate the initial ellipse parameters. The damping parameter (step size for Ahn's method), damping increment, damping decrement, and maximum number of iterations were set as = 0.5, = 10, = 3, and 50, respectively.
3-Parameter specific ellipse fitting evaluation: The last set of simulated experiments involves the assessment of the impact of variables, namely, rotation angle, aspect ratio, measurement standard deviation, and spanning elliptic arc, on the accuracy of the estimated parameters of the best fit ellipse (P-Error) using the methods of Ahn [24], Szpak [9], Kanatani [20], Taubin [16], and Halir [13], and Prasad [29] 6 . To this end, a base parameter setup is used, and only one parameter is changed at a time with the following configurations: a) rotation angle: • Base: Similar to the previous experiment, to correctly capture the behavior of random measurement errors on the performance of each ellipse fitting method, each parameter configuration was simulated 500 times and the mean of the P-Error was recorded.

Real-world experiments
In the early 1800, Dandelin showed that the intersection between a plane and a cylinder produces an ellipse. The geometric parameters of the original cylinder as a function of the parameters of the intersecting ellipse can be retrieved as follows [30]: where , and ⃗ are the radius, a point on the axis (center), and vector of cylinder's axis, respectively, �⃗ and 0 = ( 0 , 0 , 0 ) are the normal vector to and a point on the intersecting plane, respectively, = ( , , 0 ) , and are the center, semi-major length and semi-minor length of the ellipse in the coordinate system of the intersecting plane, respectively, is the rotation matrix 7 that takes the plane's normal vector to vector (0,0,1) , and is the angle between the plane's normal vector and the cylinder's axis. From this formulation, it is possible to approximate the parameters of a cylinder using the parameters of an ellipse created through an intersecting plane. To this end, this experiment is designed to evaluate the ability of different ellipse fitting methods in recovering the geometric parameters in three separate real-world cylindrical point clouds. The selected point clouds are half cylinders of mechanical pipes, acquired using a Leica HDS6100 terrestrial laser scanner (Figures 3a and 3b), as well as 3D reconstructed point cloud using twenty iPhone 11 4K images ( Figure  3c). The point clouds represent three noise levels and point densities with a normalized RMSE 8 of approximately 0.010, 0.015 and 0.020 (for the point clouds of Figures 3a, 3b and 3c, respectively). The outliers, and the ground truth cylinder parameters were automatically estimated using the robust cylinder extraction and fitting proposed in Maalek [5]. For each point cloud, 10,000 random planes passing through a randomly selected point on the cylinder were generated. The closest 50 points on the cylinder to each plane -within 1mm distance from the plane-were identified. The points are then rotated such that the normal to the plane is parallel to (0,0,1) . The best fit ellipse using the methods of Halir [13], Szpak [9], Ahn [24], and ours is applied on the − coordinates of the rotated points 9 . From the best-fit ellipse parameters estimated by each method, the geometric parameters of the original cylinder were computed using equation 21, and the errors between the estimated parameters and the ground truth were calculated.

Real-world experiments
Evaluating the ability of ellipse fitting using algorithm 1, compared to the methods of Halir [13], Szpak [9] and Ahn [24] in recovering geometric parameters of cylinders from point clouds 7 The rotation matrix can be simply obtained using the Rodrigues' rotation formulation. 8 The normalized RMSE is defined in this manuscript as the RMSE divided by the square root of the area of the circle (or ellipse). 9 Here, only methods that guaranteed an ellipse were considered.

Point to ellipse distance measure
In Section 3, a new ellipse fitting method using the confocal hyperbola distance was introduced. Here, we schematically and empirically show the agreement of the confocal hyperbola distance to the ground truth geometric distance, acquired using Chernov's method [8]. The behavior of the confocal hyperbola distance function is also compared to the algebraic, Sampson, and Harker and O'Leary distance approximations.

Visual behavior of distance approximations
Before the numerical results are presented, it is important to also study the behavior of the point to ellipse distance measurements using the confocal hyperbola, algebraic, Sampson, Harker and O'Leary [21], and Chernov [8] (the ground truth), particularly in ill-positioned points relative to the ellipse (i.e. points on the major or minor axes). To this end, consider an ellipse with geometric parameter vector of (0,0,5,3,0) . Figure  4 shows the behavior of the distance functions of points in three different arrangements, namely, points on the major axis (  Figures 4a, 4b and 4c all show that the Sampson distance tends to infinity as the point coordinates approach the center of the ellipse. Figure 4a shows that the Harker and O'Leary's distance also tends to infinity as the point coordinate approaches the focal points of the ellipse (in our case points = ±4 and = 0). Another interesting observation from Figure 4 is that none of the distances obtained using algebraic, Harker and O'Leary, and Sampson, correctly predicts the pattern of the true geometric distance (i.e. Chernov's method shown in blue). In fact, the distances of Harker and O'Leary, Sampson, and algebraic appear to move even farther away from the actual geometric distance as the point moves farther from the ellipse. The confocal hyperbola method on the other hand, is almost indistinguishable from the true geometric distance in the presented three arrangements (the confocal hyperbola, shown in red lies almost exactly underneath Chernov's distance shown in blue). From Figure 4, it can also be observed that

Numerical comparison
As observed, the visual assessment of the distance functions shows the advantages of using the confocal hyperbola method compared to other methods in predicting the behavior of the geometric distance. The results of the simulation presented in the following are, hence, used to numerically quantify the extent of this advantageous performance. To this end, 10,000 different ellipse configurations were simulated as per the properties presented in Section 4.2.2. For each measurement, the absolute difference of the point to ellipse distances using the algebraic, Sampson, Harker and O'Leary and confocal hyperbola, from that of the ground truth geometric distance was calculated. Figure 5a shows the results of the simulated measurements. The horizontal axis represents the absolute difference of each method from the ground truth in pixels, and the vertical axis represents the cumulative probability. It can be clearly observed that the confocal hyperbola achieved significantly better results compared to the remaining methods. To provide some numerical perspective, the mean, median, and 95 th percentile of the absolute distances obtained by each method from the ground truth are presented in Table 2. It can be observed that the 95 th percentile of all measurements using the confocal hyperbola fall within ±0.05 pixels from the ground truth. The confocal hyperbola is, in fact, approximately 34, 73, and 358 times more accurate than the distance approximations of Harker and O'Leary, Sampson, and algebraic (on average). It can also be observed that 50% of the distances using the confocal hyperbola are less than 10 −4 pixels from the ground truth, which is 99, 1136, and 26792 times better than those achieved by the best 50% of Harker and O'Leary, Sampson, and algebraic distances. In fact, the Sampson distance, which has been known as an excellent approximation to the geometric distance [9], is inferior to even Harker and O'Leary, which is also considerably inferior to the confocal hyperbola method.   Figure 5b shows the computation times for each method w.r.t the number of points in the dataset (from 50 to 10,000 points). It was observed that the ground truth method of Chernov was on average approximately 900, 150, 4700, and 12700 times slower than confocal hyperbola, Harker and O'Leary, Sampson, and algebraic methods, respectively. Even though the confocal hyperbola achieved about 5 times slower computation compared to Sampson, it still computed the distance of 10,000 points in around 0.7 milliseconds, while providing highly accurate geometric distance approximations.

Ellipse fitting evaluation
In this section, first, the performance (RMSE and P-Error) of the proposed ellipse fitting using confocal hyperbola in comparison to the methods of Halir [13], Taubin [16], Kanatani [17], Szpak [9], and Ahn [24] was empirically evaluated. To this end, 10,000 random ellipses were simulated using the configurations presented in Section 4.2.2. The mean, median and 95 th percentile of the RMSE and P-Error are provided in Table 3. The bolded numbers represent the best (lowest number) between all methods in each row. As shown, our ellipse fitting provides almost identical (and in some cases slightly better) RMSE and P-Error compared to the gold standard method of Ahn, and outperformed the remaining methods, particularly in the parameter estimation errors. Even though our method outperformed the rest, the difference in the RMSE between our method and the third best performing method, Szpak is still only in the order of 10 −2 pixels, which may be negligible in practical settings. On the other hand, the P-Error, which represents the agreement between the estimated and original parameters, showed around 0.1%, 5.5%, 23.7%, 26.3%, and 107.3% relative improvement compared to Ahn, Szpak, Kanatani, Taubin and Halir, respectively (the average improvement amongst all measurements).
The results, shown in Table 3, provide an overall indication of the expected performance of our method compared to established state-of-the-art methods. It is, however, also important to assess the performance (particularly the P-Error) of all methods with respect to individual variables, such as, rotation angle, aspect ratio, noise level, and spanning arc. The results of these simulations are presented in the following subsections.

Impact of rotation
Here, the impact of change in ellipse rotation angle from zero to , in 20 increments, on the P-Error for the methods of Halir, Taubin, Kanatani, Ahn, Szpak, Prasad, and confocal hyperbola was evaluated. The mean of the P-Error for all rotation angles using the different methods is provided in Table 4. As illustrated, our method achieved the best results followed very closely by Ahn's method. Figure 6 shows the results of the mean of the P-Error for the 500 separate simulations per ellipse configuration vs. the rotation angle. It can be observed that all methods, except for Prasad, remain almost constant with respect to the rotation angle. In other words, the rotation angle does not appear to impact the results of the ellipse fitting for Halir, Taubin, Kanatani, Ahn, Szpak and our method. The method of Prasad, on the other hand, is highly impacted (approximately 60% difference between the lowest and the highest value), especially at rotation angle equal to 2 (i.e. vertical ellipse), which demonstrates that the parameters estimated using Prasad's method not invariant to affine transformation of the data. The reason for this lies in the heart of the objective function of Prasad (see equation 44 of [29]), where the objective function is multiplied by the axis coordinates to enhance numerical stability of their method. Due to the systematic trend imposed through affine transformation of the data, the method of Prasad will no longer be included in the remaining comparisons of this manuscript.

Impact of aspect ratio
The aspect ratio was changed in 20 evenly spaced increments between 1 and 4. The mean of the P-Error for 500 simulations, obtained by the methods of Halir, Taubin, Kanatani, Ahn, Szpak, and confocal hyperbola are shown in Figure 7a. Three important observations were inferred from the provided chart. Firstly, all methods appear to behave similarly up until the aspect ratio of 1.5, where Halir's method starts to produce comparatively biased estimates of the best fit ellipses. Secondly, at around the aspect ratio of 2, the methods of Kanatani and Taubin start to produce comparatively biased results and appear impacted by the elongation of the ellipse. The methods, of Ahn, Szpak and ours appear to perform similarly with ours and Ahn's obtaining slightly better results on average compared to Szpak's (see Table 4). The last important observation is that all methods at the aspect ratio of 1 (i.e. a circle) appear to produce biased estimates of the geometric parameters. In fact, a relatively steep decline can be observed between the aspect ratio of 1 and 1.15. This demonstrates that the available ellipse fitting method may not necessarily be suitable for purely circle fitting applications. It is, hence, desirable to determine whether aspect ratios exist, where a circle fitting method can predict the model parameters more effectively. To this end, the ellipse fitting methods are compared to the reliable hyperaccurate circle fitting of Chernov [31] for aspect ratios between 1 and 1.012 in 0.001 increments. The results are shown in Figure 7b. It can be observed that for very small aspect ratios, it is possible to utilize the hyper accurate circle fitting method to predict the parameters of an ellipse. The circle fitting, however, starts to produce comparatively biased results at the aspect ratio of approximately above 1.005. This is since the semimajor and semi-minor of the simulated ellipses start to differ and a single radius estimation can no longer account for the differences between the semi-minor and semi-major lengths. Therefore, given the configurations used in this section, circle fitting might produce more reliable results for aspect ratio of below 1.005. In such cases, we can set = 0, = 0 and = in our formulations (equations 13, 15, and A1-7) to treat equation 17 as a circle fitting problem by reducing the degrees of freedom from five parameters to three.

Figure 7:
The mean of the parameter estimation errors for 500 simulations in each configuration vs. the aspect ratio: a) comparison of ellipse fitting methods; b) comparison of ellipse fitting with circle fitting method for small aspect ratios

Impact of standard deviation (noise)
Here, the measurement standard deviation is changed from 0 to 5 pixels at 0.25-pixel intervals. The results of the mean parameter error for each method vs. the measurement noise level is shown in Figure 8. Three main observations were drawn from the chart. Firstly, all methods appear to behave similarly up until around 1.5 pixels standard deviation, where Halir's method starts producing biased estimates of the parameters. The second observation is that at around 2.5 pixels standard deviation, the methods of Szpak and Taubin start producing relatively biased results compared to the methods of Kanatani, Ahn and ours. The last observation is that the methods of Ahn, Kanatani and ours are almost identical, and act very close to linear with respect to the standard deviation, with our method achieving an 2 = 0.99 for the best fit line. The method of Kanatani was expected to behave favorably with respect to the standard deviation since the method was specifically designed to include bias corrections to the design and scatter matrices up to a second order error term. Ahn's method and ours appear to behave as favorably as that of Kanatani with respect to changes in measurement noise levels, while also guaranteeing an ellipse, which is one of the main disadvantages of Kanatani's method (see Section 2).  (Table 4). Our method started producing relatively better parameter estimation at arc angles of less than 120° compared to all other methods. Another interesting observation was that the relationship between the estimated parameter errors and arc angle closely followed a power function (e.g. 2 = 0.95 for our method). In fact, using our method (the best performing method), the parameter estimation error for a quarter ellipse and a half ellipse were around 9.0% and 0.7% of the ground truth geometric parameters, respectively. Therefore, careful attention must be given when fitting ellipses to points on elliptic arcs of smaller than half an ellipse, particularly for applications pertaining to high precision metrology.

Computation time and number of iterations
In this section, the relative performance of our method is evaluated in terms of computation time along with the number of iterations to convergence. To this end, the base configuration of the ellipse simulation parameters of Section 4.2.2 was used. The number of points were increased from 200 to 4,000 in 200-point intervals, and at each interval 500 different simulations were performed. Figure 10a shows the mean of the computation time (for the 500 simulations at a given configuration) in logarithmic scale for each method vs. the number of points 10 . As observed, Ahn's, Szpak's and Kanatani's methods achieved the highest computation times, which were approximately 140-, 21-and 1.2-times slower than our method (average over all measurements). Our method on average was around 30 times slower than that of Halir's, which was also 1.4 times slower than Taubin's on average. To provide some perspective, ours, Szpak's and Ahn's took approximately 0.02, 0.66 and 3.91 seconds, respectively, to fitting an ellipse to 4,000 points. Not only did our method considerably outperform the other two iterative methods in terms of computation time, but it also outperformed Kanatani's direct method, particularly for number of points larger than 1,000.
The last comparison presented in this section is the number of iterations till convergence for the three iterative methods. Figure 10b shows the average number of iterations required for convergence vs. the number of points. As observed, Ahn's, Szpak's, and ours required on average approximately 19, 10 and 9 iterations, respectively to converge. An important observation is that Szpak's and ours converged after roughly the same number of iterations, yet our method achieved better computation times on average. This is since Szpak's method requires the evaluation of two objective functions at each iteration (see algorithm 2 of [9]). Not only does Ahn's method require more iterations (about twice as much as ours), which contributes to a higher computation time, but each iteration also takes longer to process. This is since in each iteration Ahn's method requires a separate generalized Gauss-Newton method to acquire the orthogonal contacting point to a given ellipse, a process which our method estimates very effectively through the confocal hyperbola distance function (see equation 13).

Summary of best fit ellipse results
The impact of four important variables, namely, rotation, aspect ratio, noise, and arc angle, on the parameter estimation using established ellipse fitting methods was evaluated. It was observed that the ellipse fitting methods of Halir, Taubin, Kanatani, Szpak, Ahn, and ours were not significantly impacted by the change on the rotation angle. The methods of Szpak, Ahn and ours performed relatively more stable than the other methods with respect to the change in aspect ratio. It was also observed that the circle fitting of Chernov can provide more accurate parameter estimation in small aspect ratios (e.g. for the configurations used in the experiments less than 1.005). The methods of Ahn, Kanatani and ours performed relatively more stable than the other methods as a function of the measurement noise with a close to linear relationship between the estimated parameter errors and measurement noise. The methods of Ahn, Szpak and ours performed more stable than others, especially in arc angles of less than (half ellipse).
To summarize the results, the mean of the parameter errors obtained by each method is presented in Table  4 with the rows representing the results obtained by the methods for each experiment. As observed, for all experiments, ours and Ahn's gold standard geometric fit achieved nearly identical parameter estimation results, and outperformed the methods of Halir, Taubin, Kanatani and Szpak. It is worth noting that the spanning arc for the simulated ellipses in the first three experiments was 2 (full ellipse). When the arc angle was changed, however, our method even outperformed that of Ahn's (Table 4; row 4). The results presented in Tables 3 and 4 show that our method achieved consistently reliable ellipse parameter estimation results, comparable to established geometric fits (Ahn's method), and outperformed reliable algebraic (Halir) and Sampson-based methods (Taubin, Kanatani and Szpak). It was also observed that given the same initial ellipse estimate and convergence criteria, our method converged faster compared to the other two iterative methods (Ahn and Szpak), and consistently converged within 9-10 iterations on average. The results of the simulation experiments demonstrate that the confocal hyperbola distance is an excellent predictor of the true geometric distance and can also be used to produce reliable ellipse fitting results. The bold numbers represent the best method for a given row

Real-world experiments
In this section, the results of fitting ellipses using the methods of Halir [13], Szpak [9], Ahn [24], and ours on 10,000 randomly oriented plane intersections with the cylinder point clouds of Figure 3 are reported. Table  5 shows the results of the mean of the parameter estimation errors for the center, radius, and axis angle of the cylinders of each case w.r.t. the ground truth. Four main observations are made from the results presented in Table 5. Firstly, the estimated center and radius of the cylinders appear to be directly impacted by the normalized RMSE ratio. In other words, the errors in the estimated center and radius (which represents the minor axis of the ellipses) appear to increase as the noise ratio increases. Second, the estimated center and radius using Halir's method is impacted the most (around 3mm and 2.5mm, respectively) with the increase of the normalized RMSE, compared to say our method (which changed around 1.5mm and 0.7mm, respectively). Based on our observations in Section 5.2.3, Halir's method started producing relatively biased results at noise levels above around 1.5 pixels noise (normalized RMSE of around 0.01); hence, it was expected that Halir's method produces more biased results compared to the remaining methods at larger noise levels in cases of Figures 3b and 3c. Halir's method also produced relatively more biased results compared to the remaining methods in each case, which was also expected since the method was shown to be produce considerably biased results for arc angles of less than around ~200° (see Figure 9). The third observation is that the axis angle error, which represents the ratio between the semi-minor and semi-major lengths, does not follow the same trend as the center and radius w.r.t. the relative noise increase. In fact, the ratio between the semi-minor and semi-major lengths appears to improve from case 3b to case 3c with very little change between cases 3a and 3c. Upon further investigation, it was observed that due to the lower point density of case 3b, the average number of points within 1mm of the intersecting planes was 30 in opposed to 45 and 48 points for cases 3a and 3c, respectively, which impacts the estimated ellipse parameters, including the ratio of the estimated semi-minor to semi-major lengths. The fourth observation is that even though Ahn's method provided more accurate radius estimations (corresponding to the semi-minor length) compared to Szpak's method, Szpak's method achieved lower axis angle errors, particularly in cases 3a and 3b. The difference is small (in the order of 10 −1 degrees) and, hence, not a concern in practical settings. However, this shows that in these configurations, Szpak's method provides a slightly better estimation of the ratio between the semi-minor and semi-major lengths compared to Ahn's. The final observation from the data was that our method outperformed all other methods in the parameter estimation of the center, radius, and axis angle errors on average. On average for all measurements, our method achieved 1.5%, 5.7% and 85.3% better center estimation results, 1.7%, 6.0%, and 99.7% better radius estimation results, and 5.2%, 2.9% and 23.7% better cylinder axis estimation results, compared to the methods of Ahn, Szpak, and Halir, respectively. The average of the absolute difference in the estimation of the center, radius, and axis angle using our method were 0.2 mm, 0.2 mm and 0.1 degrees, respectively, compared to Szpak's, and 0.1 mm, 0.1 mm and 0.2 degrees, respectively, compared to Ahn's method, which are still notable in application pertaining to high precision metrology and reliability.
From the results of the simulation experiments presented in Section 5, the methods of Szpak, Ahn and ours behaved similarly in most configurations, other than for noise levels above 2.5 pixels (normalized RMSE of around 0.02 and above), where Szpak's method started to produce relatively biased results. The configurations in the real-world experiment were also in the normalized RMSE range, where Szpak's method was expected to behave similarly to ours and Ahn's. In larger noise levels, however, as shown in Figure 8, ours and Ahn's are expected to perform better. In Szpak [9], it was mentioned that the Sampson based methods are a good approximation for "moderate" noise levels. From our experiments, the 0.02 normalized RMSE ratio might be a good quantification of the definition of "moderate" in this context. Table 5: Results of the mean of the parameter estimation errors for the cases presented in Figure 3 Mean of Cylinder Parameter Errors (Case of Figure 3a)

Conclusions
This manuscript provided a new method for fitting ellipses to 2D points using the confocal hyperbola point to ellipse geometric distance approximation. Ellipse fitting using confocal hyperbola had not been formulated in the literature and, hence, its performance was never investigated. To this end, two comprehensive sets of experiments, including simulation-based and real-world were designed. The effectiveness of the proposed method compared to established state-of-the-art ellipse fitting methods in various configurations was examined. The simulation-based experiments consisted of two main categories, namely, point to ellipse distance approximations, and best fit ellipse evaluation. For point to ellipse geometric distance approximation methods, the confocal hyperbola method outperformed other established geometric distance approximations, namely, algebraic, Sampson, and Harker and O'Leary, and behaved almost identical (in the order of 10 −2 pixels on average) to the ground truth geometric distance, proposed by Chernov [8] . To this end, a new ellipse fitting method using confocal hyperbola was proposed (Algorithm 1) and its effectiveness was evaluated compared to the ellipse fitting methods proposed by Halir [13], Taubin [16], Kanatani [17], Szpak [9], and Ahn [24]. The effects of change of variables, namely, rotation angle, aspect ratio, noise, and spanning arc on the estimated ellipse parameters using each method were thoroughly investigated. Overall, it was observed that the proposed confocal hyperbola-based ellipse fitting method achieved almost identical to (and in some cases outperformed) the gold standard geometric ellipse fit of Ahn, and outperformed the remaining methods in all categories. It was also shown that in all simulations, our method converged after around 9-10 iterations and achieved an average computation speed of roughly 140, 21 and 1.2 times faster than those achieved using the methods of Ahn, Szpak, and Kanatani.
The real-world experiment was designed to determine the ability of different ellipse fitting methods in predicting the parameters of cylindrical point clouds, acquired from mechanical pipes. To this end three cylindrical point clouds with different noise levels, and point densities were considered. The methods of Halir, Szpak, Ahn, and ours, which are all methods guaranteeing ellipses, were considered to estimate ellipse parameters of intersecting planes to these cylinders. The ellipse parameters were then used to recover the original cylinder parameter. The errors in estimating the original cylinder parameters were then recorded. It was observed that our method outperformed the remaining methods in estimating the cylinder's parameters. The method of Halir was clearly less accurate than the remaining methods, which was expected due to the biased behavior of Hair's ellipse parameter estimation in smaller arc angles (approximately less than a semiellipse) and at larger noise levels.
Overall, the comprehensive evaluation in both simulation-based and real-world experiments showed the agreement between the ellipse parameters estimated using the proposed confocal hyperbola ellipse fitting method and the true geometric parameters of the ellipses. The method achieved the best parameter estimation results in all categories of simulation and real-world experiments. We conclude that the confocal hyperbolabased ellipse fitting method shows great promise for fitting ellipses to 2D data in practical settings.