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.

Recommended Posts

Approximation by the least squares method of a circle with two fixed points
Approximation of distance between two points on the surface of a spheroid (on the surface of the earth)
Calculation of homography matrix by least squares method (DLT method)
Reuse the behavior of the @property method by using a descriptor [16/100]
Hit a method of a class instance with the Python Bottle Web API
Drawing a circle passing through two points
Single regression analysis by least squares method
Try sending the aggregated results of two records by email with pykintone
Find the coefficients of the least squares polynomial
About the accuracy of Archimedean circle calculation method
I tried the least squares method in Python
Approximate a Bezier curve through a specified point using the least squares method in Python
Destroy the intermediate expression of the sweep method with Python
Take a screenshot of the LCD with Python-LEGO Mindstorms
Python points from the perspective of a C programmer
Visualize the characteristic vocabulary of a document with D3.js
Calculate the product of matrices with a character expression?
The copy method of pandas.DataFrame is deep copy by default
Create a 2D array by adding a row to the end of an empty array with numpy
How to plot a lot of legends by changing the color of the graph continuously with matplotlib
Save the output of GAN one by one ~ With the implementation of GAN by PyTorch ~
A simple Python implementation of the k-nearest neighbor method (k-NN)
A network diagram was created with the data of COVID-19.
Measure the importance of features with a random forest tool
Get the id of a GPU with low memory usage
Get UNIXTIME at the beginning of today with a command
Analyze the topic model of becoming a novelist with GensimPy3
Numerical approximation method when the calculation of the derivative is troublesome
The story of making a question box bot with discord.py
A Python script that compares the contents of two directories
Let's prove the addition theorem of trigonometric functions by replacing the function with a function in SymPy (≠ substitution)
Find the white Christmas rate by prefecture with Python and map it to a map of Japan
How to intercept or tamper with the SSL communication of the actual iOS device by a proxy
[Python] How to use the for statement. A method of extracting by specifying a range or conditions.