由圆上三点确定圆心和半径(附Python&Matlab程序)

更多阅读：sppy.site

如何计算曲线   y ( x )   ~y(x)~  y(x) 上的曲率，而曲线是由若干离散点构成。我的第一反应是根据离散点差分得到一阶导数   y ′   ~y'~  y′ 和二阶导数   y ′ ′   ~y''~  y′′ ，然后由下式计算 k = ∣ y ′ ′ ∣ ( 1 + y ′ 2 ) 3 / 2 (0a) \tag{0a} k=\frac{|y''|}{(1+y'^2)^{3/2}} k=(1+y′2)3/2∣y′′∣​(0a)

曲线的凹凸方向则由   y ′ ′   ~y''~  y′′ 的符号确定。

后面，我想到了可以根据曲率的几何意义来计算，即 k = 1 r (0b) \tag{0b} k=\frac{1}{r} k=r1​(0b)

式中， r   r~ r 是该点的曲率半径，可以通过该点及其两个相邻点得到（不共线的三点确定一个圆）。

三个已知点的坐标分别记为   ( x 1 , y 1 ) ~(x_1,y_1)  (x1​,y1​)、 ( x 2 , y 2 ) (x_2,y_2) (x2​,y2​)、 ( x 3 , y 3 ) (x_3,y_3) (x3​,y3​)，圆的一般方程可写为二次多项式，即 A ( x 2 + y 2 ) + B x + C y + D = 0 (1a) \tag{1a} A(x^2+y^2)+Bx+Cy+D=0 A(x2+y2)+Bx+Cy+D=0(1a)

将式 ( 1 a ) (1\mathrm{a}) (1a)变形可得圆的标准方程，即 ( x + B 2 A ) 2 + ( y + C 2 A ) 2 = B 2 + C 2 − 4 A D 4 A 2 (1b) \tag{1b} \bigg(x+\frac{B}{2A}\bigg)^2+\bigg(y+\frac{C}{2A}\bigg)^2=\frac{B^2+C^2-4AD}{4A^2} (x+2AB​)2+(y+2AC​)2=4A2B2+C2−4AD​(1b)

将三个已知点代入式 ( 1 a ) (1\mathrm{a}) (1a)，可得关于   A ~A  A、 B B B、 C C C、 D   D~ D 的齐次线性方程组，即 [ x 2 + y 2 x y 1 x 1 2 + y 1 2 x 1 y 1 1 x 2 2 + y 2 2 x 2 y 2 1 x 3 2 + y 3 2 x 3 y 3 1 ] ⋅ [ A B C D ] = [ 0 0 0 0 ] (2) \tag{2} \begin{bmatrix} x^2+y^2 & x & y & 1 \\[3pt] x^2_1+y^2_1 & x_1 & y_1 & 1 \\[3pt] x^2_2+y^2_2 & x_2 & y_2 & 1 \\[3pt] x^2_3+y^2_3 & x_3 & y_3 & 1 \end{bmatrix} \cdot \begin{bmatrix} A\\ B\\ C\\ D \end{bmatrix} = \begin{bmatrix} 0\\ 0\\ 0\\ 0 \end{bmatrix} ⎣⎢⎢⎢⎡​x2+y2x12​+y12​x22​+y22​x32​+y32​​xx1​x2​x3​​yy1​y2​y3​​1111​⎦⎥⎥⎥⎤​⋅⎣⎢⎢⎡​ABCD​⎦⎥⎥⎤​=⎣⎢⎢⎡​0000​⎦⎥⎥⎤​(2)

显然在三点不共线的前提下，该齐次线性方程组有非零解，其等价于系数矩阵不满秩，即有 ∣ x 2 + y 2 x y 1 x 1 2 + y 1 2 x 1 y 1 1 x 2 2 + y 2 2 x 2 y 2 1 x 3 2 + y 3 2 x 3 y 3 1 ∣ = 0 (3) \tag{3} \begin{vmatrix} x^2+y^2 & x & y & 1 \\[3pt] x^2_1+y^2_1 & x_1 & y_1 & 1 \\[3pt] x^2_2+y^2_2 & x_2 & y_2 & 1 \\[3pt] x^2_3+y^2_3 & x_3 & y_3 & 1 \end{vmatrix} = 0 ∣∣∣∣∣∣∣∣∣​x2+y2x12​+y12​x22​+y22​x32​+y32​​xx1​x2​x3​​yy1​y2​y3​​1111​∣∣∣∣∣∣∣∣∣​=0(3)

将式 ( 3 ) (3) (3)展开，并与式 ( 1 ) (1) (1)对比可得四个系数，即 A = + ∣ x 1 y 1 1 x 2 y 2 1 x 3 y 3 1 ∣ (4a) \tag{4a} A=+\begin{vmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{vmatrix} A=+∣∣∣∣∣∣​x1​x2​x3​​y1​y2​y3​​111​∣∣∣∣∣∣​(4a)

B = − ∣ x 1 2 + y 1 2 y 1 1 x 2 2 + y 2 2 y 2 1 x 3 2 + y 3 2 y 3 1 ∣ (4b) \tag{4b} B=-\begin{vmatrix} x^2_1+y^2_1 & y_1 & 1 \\[3pt] x^2_2+y^2_2 & y_2 & 1 \\[3pt] x^2_3+y^2_3 & y_3 & 1 \end{vmatrix} B=−∣∣∣∣∣∣∣​x12​+y12​x22​+y22​x32​+y32​​y