A new fitting method for measurement of the curvature radius of a short arc with high precision

The traditional LSF method is based on minimizing the mean square distance from data points to the fitting arc, the objective function is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L=\sum\limits_{i=1}^{n}{d_{i}^{2}=}\sum\limits_{i=1}^{n}{{{(\sqrt{{{({{x}_{i}}-a)}^{2}}+{{({{y}_{i}}-b)}^{2}}}-R)}^{2}}}\nonumber \end{align} \tag{ 1 }$$

where ${{d}_{i}}$ is the geometric distance from the point with the coordinate of $({{x}_{i}},{{y}_{i}})$ to the fitting arc [11–14].

This algorithm is based on a three-parameter circle equation. However, it has been proved that the closed-form solution based on a three-parameter circle equation cannot be obtained in a general case. So, we define the four-parameter equation of circle as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & A({{x}^{2}}+{{y}^{2}})+Bx+Cy+D=0 \nonumber \\ & {\rm s}{\rm .t}{\rm .}\,\,A\ne 0,{{B}^{2}}+{{C}^{2}}-4AD>0. \end{array}\nonumber \end{align} \tag{ 2 }$$

The above four-parameter equation has the following advantages:

• (1)  
  Matrix operation of the equation is much easier and we can obtain the closed-form solution.
• (2)  
  It unifies the equation of the circle and straight line. It gives a line when A  =  0 and a circle when $A\ne 0$.
• (3)  
  The parameters, A, B, C and D, are no longer arbitrarily large and they are reduced to a specific range.
• (4)  
  The objective function is smooth in new four-parameter space and the iterative algorithm ensures convergence of the equation.

The four-parameter equation (2) can also be transformed to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\left(x+\frac{B}{2A} \right)}^{2}}+{{\left(y+\frac{C}{2A} \right)}^{2}}=\frac{{{B}^{2}}+{{C}^{2}}-4AD}{4{{A}^{2}}}.\nonumber \end{align} \tag{ 3 }$$

Then the center coordinate (a, b) and the radius R of circle can be determined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}lllllllllllllllllllllllllllll@{}} a=-\frac{B}{2A} \nonumber \\ b=-\frac{C}{2A} \nonumber \\ R=\sqrt{\frac{{{B}^{2}}+{{C}^{2}}-4AD}{4{{A}^{2}}}} \end{array} \right..\nonumber \end{align} \tag{ 4 }$$

Giving n sampling points on the contour of the arc (x_i, y_i) (i  =  1,2,...,n), the minimizing objective function based on the four-parameter equation of a circle can be defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L=\sum\limits_{i=1}^{n}{{{[A(x_{i}^{2}+y_{i}^{2})+B{{x}_{i}}+C{{y}_{i}}+D]}^{2}}}\nonumber \end{align} \tag{ 5 }$$

where parameters A, B, C and D must satisfy the inequality as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{B}^{2}}+{{C}^{2}}-4AD>0.\nonumber \end{align} \tag{ 6 }$$

So the arc fitting problem can be regarded as the following optimization problem

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \min \,\,L{\rm =}\sum\limits_{i=1}^{n}{{{[A(x_{i}^{2}+y_{i}^{2})+B{{x}_{i}}+C{{y}_{i}}+D]}^{2}}} \nonumber \\ & {\rm s}{\rm .t}{\rm .}\,\,\,\,{{B}^{2}}+{{C}^{2}}-4AD>0. \end{array}\nonumber \end{align} \tag{ 7 }$$

Assuming ${{z}_{i}}=x_{i}^{2}+y_{i}^{2}$, the objective function L can be expressed in matrix form as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L={{{\bf P}}^{{\rm T}}}{\bf MP}\nonumber \end{align} \tag{ 8 }$$

where ${\bf P}={{[A,B,C,D]}^{\operatorname{T}}}$ and $\newcommand{\e}{{\rm e}} {\bf M}=\left[\begin{array}{@{}cccccccccccccccccccc@{}} Mzz & Myz & Mxz & Mz \\ Mxz & Mxy & Mxx & Mx \\ Myz & Mxy & Myy & My \\ Mz & Mx & My & n \end{array} \right]$

The elements in the matrix M are vector moments, for example

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Mxz=\sum\limits_{i=1}^{n}{{{x}_{i}}{{z}_{i}}{\rm ,}\quad Mzz=}\sum\limits_{i=1}^{n}{{{z}_{i}}{{z}_{i}}}{\rm ,}\quad Mz=\sum\limits_{i=1}^{n}{{{z}_{i}}}.\nonumber \end{align}$$

The inequality [6] can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{{\bf P}}^{\operatorname{T}}}{\bf QP}>0\nonumber \end{align} \tag{ 9 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf Q}=\left[\begin{array}{@{}cccccccccccccccccccc@{}} 0 & 0 & 0 & -2 \nonumber \\ 0 & 1 & 0 & 0 \nonumber \\ 0 & 0 & 1 & 0 \nonumber \\ -2 & 0 & 0 & 0 \end{array} \right].\nonumber \end{align}$$

Now the minimizing objective function can be defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & {\rm min}\quad L={{{\bf P}}^{\operatorname{T}}}{\bf MP} \nonumber \\ & {\rm s}{\rm .t}{\rm .}\quad \,\,{{{\bf P}}^{\operatorname{T}}}{\bf QP}>0. \end{array}\nonumber \end{align} \tag{ 10 }$$

This is an optimization problem with inequality constraint. To solve the above problem, a GLF is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{L}^{{\rm *}}}={{{\bf P}}^{\operatorname{T}}}{\bf MP}-\eta {{{\bf P}}^{\operatorname{T}}}{\bf QP}\nonumber \end{align} \tag{ 11 }$$

where η (η  ⩾  0) is the Lagrange multiplier.

Now the objective function can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \underset{P}{\mathop{\min}}\,L{\rm =}\underset{P}{\mathop{\min}}\,\underset{\eta}{\mathop{\max}}\,{{L}^{*}}\nonumber \end{align} \tag{ 12 }$$

where $\newcommand{\e}{{\rm e}} \underset{P}{\mathop{\min}}\,\underset{\eta}{\mathop{\max}}\,{{L}^{*}}$ is defined as the minimax problem of the GLF and it is also called the primal problem. According to the duality of the Lagrange function, the dual problem of the primal problem is a minimax problem, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \underset{\eta}{\mathop{\max}}\,\underset{P}{\mathop{\min}}\,{{L}^{*}}.\nonumber \end{align}$$

The solution of the primal problem can be obtained by analyzing the dual problem.

Differentiating the equation with respect to P gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf MP}=\eta {\bf QP}\nonumber \end{align} \tag{ 13 }$$

where P is the generalized eigenvector and η is the corresponding generalized eigenvalue.

By introducing equation (13) to the objective function (10), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L={{{\bf P}}^{\operatorname{T}}}{\bf MP}{\rm =}\eta {{{\bf P}}^{\operatorname{T}}}{\bf QP}.\nonumber \end{align} \tag{ 14 }$$

Since ${{{\bf P}}^{\operatorname{T}}}{\bf QP}>0$, the minimum of the objective function L corresponds to the smallest nonnegative generalized eigenvalue. Choosing the eigenvector corresponding to the minimum nonnegative generalized eigenvalue ${{[A,B,C,D]}^{T}}$ as the calculation result, the radius and center coordinates of the arc can be calculated by equation (4).

Based on the parameter constraint arc fitting algorithm described above, the radius and center coordinates of the circular arc can be obtained [15]. To further improve its fitting precision, the GD algorithm was adopted. In this method, the original value of the center coordinate is (a₀, b₀), which is used to find the optimal arc parameter R because R can be calculated by the center coordinate.

According to the least square circle fitting, the objective function can be defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L=\sum\limits_{i=1}^{n}{{{(\sqrt{{{({{x}_{i}}-a)}^{2}}+{{({{y}_{i}}-b)}^{2}}}-R)}^{2}}}.\nonumber \end{align} \tag{ 15 }$$

This is a nonlinear problem and it has no closed-form solution. So a GD optimization method can be adopted to solve this objective function.

The gradient is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle d=\left(\frac{\partial L}{\partial a},\frac{\partial L}{\partial b} \right)\nonumber \end{align} \tag{ 16 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial L}{\partial a}=2\sum\limits_{i=1}^{n}{(a-{{x}_{i}})\frac{\sqrt{{{({{x}_{i}}-a)}^{2}}+{{({{y}_{i}}-b)}^{2}}}-R}{\sqrt{{{({{x}_{i}}-a)}^{2}}+{{({{y}_{i}}-b)}^{2}}}}}\nonumber \end{align}$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial L}{\partial b}=2\sum\limits_{i=1}^{n}{(b-{{y}_{i}})\frac{\sqrt{{{({{x}_{i}}-a)}^{2}}+{{({{y}_{i}}-b)}^{2}}}-R}{\sqrt{{{({{x}_{i}}-a)}^{2}}+{{({{y}_{i}}-b)}^{2}}}}.}\nonumber \end{align}$$

The iterative formula is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & a:=a-\lambda \frac{\partial L}{\partial a} \nonumber \\ & b:=b-\lambda \frac{\partial L}{\partial b} \end{array}\nonumber \end{align} \tag{ 17 }$$

where, λ is the step length. Choosing proper step and threshold, the iteration starts and runs until the gradient is less than the threshold.

When the center of the circle (a,b) is determined, the optimal radius R can be solved based on least squares principle.

The flowchart of GD arc fitting is given in figure 1

Zoom In Zoom Out Reset image size

To verify the influence of the numbers of sampling points on the measurement precision of R, we made a further study. Setting the central angle as 10°, the Gauss noise of edge point as 0.4 mm, and the number of point n to be in the range of [100, 1600]. We generate n sampling points and calculate the measurement errors. The results are presented in figure 2. From the figures, we can see that the measurement error of the radius decreases with the increase of n. Our PCF method has the highest accuracy compared with other arc measurement methods under the same conditions.

Zoom In Zoom Out Reset image size

An ideal short arc sampling point set (n  =  1000) with a 14° central angle was generated and the theoretical radius of the arc was 120 mm. A Gauss noise by different standard deviation from 0.0 mm to 0.4 mm was added to the data randomly. The above four algorithms were adopted to calculate the fitting radius, respectively, under different noise standard deviations. The relative deviation of the fitting radius from the theoretical radius was obtained as the final fitting radius. Every algorithm was repeated 100 times under the same noise standard deviation. The average of relative deviations of fitting radius is presented in figure 3. With the increase of the Gauss noise, the related deviations of the LSF method increase rapidly, and the PSSC method shows a similar trend but its fitting precision is much better than that of the LSF method. Both the RC method (the constraint range is in 120 mm  ±  1 mm) and PCF method can keep almost the same fitting precision even the noise increases gradually, and their precisions are much better than that of the PSSC method. By contrast, the fluctuation of the RC method is larger than that of the PCF method. So, the PCF method has better robustness than others.

Zoom In Zoom Out Reset image size

Generating a short arc edge point set (n  =  1000) whose central angle, theoretical radius and range of Gauss noise are 20°, 120 mm and 0–0.08 mm, respectively, and the four algorithms above were adopted to calculate the fitting radius. Under different central angles from 2° to 20°, the relative deviations of fitting radius from the theoretical radius were obtained. Every algorithm repeats 100 times under the same central angle. The average relative deviations are presented in figure 4. With the decreasing of the central angle, the relative deviation of fitting radius with the LSF method increases rapidly when the center angle is less than 20°. The PSSC method shows the similar trend when the central angle is less than 12°. The relative deviation of the PCF method begins to increase only when the center angle is less than 4° and this relative deviation is much less than the other two methods under the same central angle. Only the relative deviation of the RC method always keeps at the same minimum level even if the central angle decreases gradually to 2°. So, the proposed PCF method has better precision without any prior information.

Zoom In Zoom Out Reset image size

Now, we comprehensively study the influence of noise and central angles. Setting the interval of central angle as [6°, 30°], standard deviation interval of Gauss noise of edge point as [0.0 mm, 0.4 mm], we generated 1000 arc edge points for every central angle and Gauss noise. Using the above algorithms, we got the radius. The absolute value of arc radius errors was then computed. The results are presented with a 3D plot and pseudocolor plot in figure 5. All algorithms suffer increasing measurement error with the decreasing of the central angle and the increasing noise. Our PCF method gives more precise measurement than the least squares method and GGS method, and its fitting precision is also better than the Nievergelt method under the same conditions.

Zoom In Zoom Out Reset image size