[PYTHON] Approximation by the least squares method of a circle with two fixed points

As with a straight line, you can approximate a circle from points P1, P2, ..., Pn using the least squares method. As a general method, I will refer to other articles, and here I will find a circle that always passes through the start point P1 and end point Pn and that the error is small. Also, try implementing such a function in python (numpy).

policy

As with a straight line, find the difference (residual) between the theoretical value and the measured value for each point, and find a value that reduces the sum of squares.

Generally, three parameters of center coordinate $ (x, y) $ and radius $ r $ are required to represent a circle, but the center of the circle passing through the two points $ P_1 and P_n $ is always the vertical bisector of these two points. It is on the bisector. Once the center is determined, the radius is also determined, so the parameters to be obtained this time can be unified.

Method

It is not impossible to display the center and radius of the circle as parameters using a vector, but it will complicate the calculation, so we will take another method. Let's consider the case of $ P_1 (-a, 0), P_n (a, 0) $ for the time being. At this time, the center is the point $ (0, y_0) $ on the y-axis, the radius $ r ^ 2 = y_0 ^ 2 + a ^ 2 $, and the equation of the circle to be obtained is

x^2+y^2-2yy_0-a^2=0

It will be. Find this $ y_0 $.

The residual for each point is the difference between the distance from each point to the center and the radius, but for the sake of simplicity of calculation, consider the difference after squared in advance. In other words, the residual of each point is calculated by substituting $ (x_k, y_k) $ into the above formula.

x_k^2+y_k^2-2y_ky_0-a^2

is. The sum of squares of this, that is

\sum \{ x_k^2+y_k^2-2y_ky_0-a^2 \} ^2

Find $ y_0 $ that minimizes. Since this formula is a quadratic function of $ y_0 $, $ y_0 $ is calculated by setting the partial derivative of $ y_0 $ to 0.

\begin{align}
\sum -4y_k( x_k^2+y_k^2-2y_ky_0-a^2 ) &= 0 \\
\sum x_k^2y_k+\sum y_k^3-2y_0\sum y_k^2-a^2\sum y_k  &= 0 \\
\end{align}

So $ y_0 $ is:

y_0 = \frac{\sum x_k^2y_k+\sum y_k^3-a^2\sum y_k}{2\sum y_k^2}

Coordinate transformation

Now, to apply this to any point cloud, the point $ P_1 (x_1, y_1) $ goes to $ (-a, 0) $ and the point $ P_n (x_n, y_n) $ for the point cloud. Must be converted so that is moved to $ (a, 0) $. for that purpose,

① Perform translation so that the midpoint $ P_c $ of $ P_1 $ and $ P_n $ moves to the origin. ② Rotate so that $ P_1 $ and $ P_n $ move on the x-axis.

Two steps are required. I will omit the detailed calculation, but you can see that the following conversion should be performed.

\left( \begin{array}{c} x' \\ y' \end{array} \right) = 
 -\frac{1}{a}\left(
    \begin{array}{cc}
      x_1-x_c & -(y_1-y_c)\\
      y_1-y_c & x_1-x_c
    \end{array}
  \right)
\left( \begin{array}{c} x-x_c \\ y-y_c \end{array} \right)

However,x_c,y_cIsP_cCoordinates,aIs|P_1P_c|is. After performing this transformation on all points, approximate it, and then perform the inverse transformation to obtain the center and radius.

Implementation

I implemented it with python.

def approxByArc(points):
	start = points[0]
	end = points[-1]
	sample = points[1:-1]
	
	midpoint = (start + end) /2

	start2 = start - midpoint
	a = np.linalg.norm(start2)

	x = start2[0]
	y = start2[1]

	rotateMatrix = np.array([[x,y],[-y,x]]) * -1 / a

	def conversion(point):
		point = point - midpoint
		p = np.array([[point[0]], [point[1]]])
		p = rotateMatrix @ p
		p = np.array([p[0][0],p[1][0]])
		return p

	sample = np.apply_along_axis(conversion, 1, sample)

	xs = np.array([i[0] for i in sample])
	ys = np.array([i[1] for i in sample])

	p1 = np.sum(xs*xs*ys)
	p2 = np.sum(ys*ys*ys)
	p3 = np.sum(ys) * a * a
	p4 = np.sum(ys*ys) * 2

	y0 = (p1+p2-p3)/p4

	center = np.array([[0],[y0]])
	center = np.linalg.inv(rotateMatrix) @ center

	center = np.array([center[0][0], center[1][0]])
	center = center + midpoint

	radius = np.linalg.norm(center - start)
	return center,radius

Let's see if it works.

points = []
points.append([60,50])
n = 40
for i in range(1,n):
	angle = math.pi / n * i
	radius = 10 + random.uniform(-0.4,0.4)

	x = 50 + radius * math.cos(angle)
	y = 50 + radius * math.sin(angle)

	points.append([x,y])
points.append([40,50])

r,c = approxByArc(np.array(points))
print(r)
print(c)

[50.         49.99319584]
10.00000231483068

It went well.

in conclusion

The method of converting to simple coordinates and then calculating is easier to understand the program, but it is not suitable for manual calculation. But no one manually calculates these complicated steps, so it's okay. maybe.