止于至善

0%

由圆上三点确定圆心和半径

发表于 更新于 分类于

背景

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

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

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

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

公式

平面三点

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

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

将三个已知点代入式\((1\mathrm{a})\)，可得关于\(~A\)、\(B\)、\(C\)、\(D~\)的齐次线性方程组，即 \[ \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} \]

显然在三点不共线的前提下，该齐次线性方程组有非零解，其等价于系数矩阵不满秩，即有 \[ \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 \]

将式\((3)\)展开，并与式\((1)\)对比可得四个系数，即 \[ \tag{4a} A=+\begin{vmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{vmatrix} \]

\[ \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} \]

\[ \tag{4c} C=+\begin{vmatrix} x^2_1+y^2_1 & x_1 & 1 \\[3pt] x^2_2+y^2_2 & x_2 & 1 \\[3pt] x^2_3+y^2_3 & x_3 & 1 \end{vmatrix} \]

\[ \tag{4d} D=-\begin{vmatrix} x^2_1+y^2_1 & x_1 & y_1 \\[3pt] x^2_2+y^2_2 & x_2 & y_2 \\[3pt] x^2_3+y^2_3 & x_3 & y_3 \end{vmatrix} \]

由式\((1\mathrm{b})\)可得圆心坐标\(~(x_c,y_c)~\)和半径\(~r\)，即 \[ \tag{5} \begin{cases} x_c &=~ -\displaystyle\frac{B}{2A} \\[8pt] y_c &=~ -\displaystyle\frac{C}{2A} \\[8pt] r &=~ \sqrt{\displaystyle\frac{B^2+C^2-4AD}{4A^2}} \end{cases} \]

\(A~\)可能是负数

空间三点

进一步推广到空间中，三个已知点的坐标分别记为\(~(x_1,y_1,z_1)\)、\((x_2,y_2,z_2)\)、\((x_3,y_3,z_3)\)。

空间中确定一个圆至少需要7个参数，一般取：

  1. 圆心\(~(x_c,y_c,z_c)\)
  2. 半径\(~r\)
  3. 圆所在平面的法向量\(~(m,n,p)\)

将三个点所在的球的一般方程写作 \[ \tag{6a} A(x^2+y^2+z^2)+Bx+Cy+Dz+E=0 \]

将式\((6\mathrm{a})\)变形可得球的标准方程，即 \[ \tag{6b} \bigg(x+\frac{B}{2A}\bigg)^2+\bigg(y+\frac{C}{2A}\bigg)^2+\bigg(z+\frac{D}{2A}\bigg)^2=\frac{B^2+C^2+D^2-4AE}{4A^2} \]

将三个已知点代入式\((6\mathrm{a})\)，可得关于\(~A\)、\(B\)、\(C\)、\(D\)、\(E~\)的齐次线性方程组，即 \[ \tag{7} \begin{bmatrix} x^2+y^2+z^2 & x & y & z & 1 \\[3pt] x^2_1+y^2_1+z^2_1 & x_1 & y_1 & z_1 & 1 \\[3pt] x^2_2+y^2_2+z^2_2 & x_2 & y_2 & z_2 & 1 \\[3pt] x^2_3+y^2_3+z^2_3 & x_3 & y_3 & z_3 & 1 \end{bmatrix}_{4\times5} \cdot \begin{bmatrix} A\\ B\\ C\\ D\\ E \end{bmatrix} = \begin{bmatrix} 0\\ 0\\ 0\\ 0 \end{bmatrix} \]

考虑使用式\((3)\)中的做法，但式\((7)\)中的系数矩阵并非方阵。于是利用三个已知点来确定平面\(~\pi\)，其方程为 \[ \tag{8} mx+ny+pz+q=0 \]

将三个已知点代入式\((8)\)，可得关于\(~m\)、\(n\)、\(p\)、\(q~\)的齐次线性方程组，即 \[ \tag{9} \begin{bmatrix} x & y & z & 1 \\ x_1 & y_1 & z_1 & 1 \\ x_2 & y_2 & z_2 & 1 \\ x_3 & y_3 & z_3 & 1 \end{bmatrix} \cdot \begin{bmatrix} m\\ n\\ p\\ q \end{bmatrix} = \begin{bmatrix} 0\\ 0\\ 0\\ 0 \end{bmatrix} \]

同理，在三点不共线的前提下，该齐次线性方程组有非零解，其等价于系数矩阵不满秩。于是将其系数矩阵的行列式展开，并与式\((8)\)对比可得四个系数，即 \[ \tag{10a} m=+\begin{vmatrix} y_1 & z_1 & 1 \\ y_2 & z_2 & 1 \\ y_3 & z_3 & 1 \end{vmatrix} \]

\[ \tag{10b} n=-\begin{vmatrix} x_1 & z_1 & 1 \\ x_2 & z_2 & 1 \\ x_3 & z_3 & 1 \end{vmatrix} \]

\[ \tag{10c} p=+\begin{vmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{vmatrix} \]

\[ \tag{10d} q=-\begin{vmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_2 & y_3 & z_3 \end{vmatrix} \]

由式\((6\mathrm{b})\)可知 \[ \tag{11} x_c=-\frac{B}{2A}, \quad y_c=-\frac{C}{2A}, \quad z_c=-\frac{D}{2A} \]

考虑到圆心\(~(x_c,y_c,z_c)~\)也在平面\(~\pi~\)上，于是将式\((11)\)代入式\((8)\)可得 \[ \tag{12} 2qA-mB-nC-pD=0 \]

将式\((12)\)与式\((7)\)结合，可得 \[ \tag{13} \begin{bmatrix} x^2+y^2+z^2 & x & y & z & 1 \\[3pt] x^2_1+y^2_1+z^2_1 & x_1 & y_1 & z_1 & 1 \\[3pt] x^2_2+y^2_2+z^2_2 & x_2 & y_2 & z_2 & 1 \\[3pt] x^2_3+y^2_3+z^2_3 & x_3 & y_3 & z_3 & 1 \\[3pt] 2q & -m & -n & -p & 0 \end{bmatrix}_{5\times5} \cdot \begin{bmatrix} A\\ B\\ C\\ D\\ E \end{bmatrix} = \begin{bmatrix} 0\\ 0\\ 0\\ 0\\ 0 \end{bmatrix} \]

这样，在三点不共线的前提下，该齐次线性方程组有非零解，其等价于系数矩阵不满秩。于是将其系数矩阵的行列式展开，并与式\((6)\)对比可得五个系数，即 \[ \tag{14a} A=+\begin{vmatrix} x_1 & y_1 & z_1 & 1 \\ x_2 & y_2 & z_2 & 1 \\ x_3 & y_3 & z_3 & 1 \\ -m & -n & -p & 0 \end{vmatrix} \]

\[ \tag{14b} B=-\begin{vmatrix} x^2_1+y^2_1+z^2_1 & y_1 & z_1 & 1 \\[3pt] x^2_2+y^2_2+z^2_2 & y_2 & z_2 & 1 \\[3pt] x^2_3+y^2_3+z^2_3 & y_3 & z_3 & 1 \\[3pt] 2q & -n & -p & 0 \end{vmatrix} \]

\[ \tag{14c} C=+\begin{vmatrix} x^2_1+y^2_1+z^2_1 & x_1 & z_1 & 1 \\[3pt] x^2_2+y^2_2+z^2_2 & x_2 & z_2 & 1 \\[3pt] x^2_3+y^2_3+z^2_3 & x_3 & z_3 & 1 \\[3pt] 2q & -m & -p & 0 \end{vmatrix} \]

\[ \tag{14d} D=-\begin{vmatrix} x^2_1+y^2_1+z^2_1 & x_1 & y_1 & 1 \\[3pt] x^2_2+y^2_2+z^2_2 & x_2 & y_2 & 1 \\[3pt] x^2_3+y^2_3+z^2_3 & x_3 & y_3 & 1 \\[3pt] 2q & -m & -n & 0 \end{vmatrix} \]

\[ \tag{14e} E=+\begin{vmatrix} x^2_1+y^2_1+z^2_1 & x_1 & y_1 & z_1 \\[3pt] x^2_2+y^2_2+z^2_2 & x_2 & y_2 & z_2 \\[3pt] x^2_3+y^2_3+z^2_3 & x_3 & y_3 & z_3 \\[3pt] 2q & -m & -n & -p \end{vmatrix} \]

综上，由式\((10)\)得到\(~m\)、\(n\)、\(p\)、\(q~\)，并由式\((14)\)得到\(~A\)、\(B\)、\(C\)、\(D\)、\(E~\)后，可得圆心坐标\(~(x_c,y_c,z_c)~\)和半径\(~r\)，即 \[ \tag{15} \begin{cases} x_c &=~ -\displaystyle\frac{B}{2A} \\[8pt] y_c &=~ -\displaystyle\frac{C}{2A} \\[8pt] z_c &=~ -\displaystyle\frac{D}{2A} \\[8pt] r &=~ \sqrt{\displaystyle\frac{B^2+C^2+D^2-4AE}{4A^2}} \end{cases} \]

\(A~\)可能是负数

程序

MATLAB代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 由圆上三点确定圆心和半径
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT
% p1   :  - 第一个点坐标, 行向量 1x3
% p2   :  - 第二个点坐标, 行向量 1x3
% p3   :  - 第三个点坐标, 行向量 1x3
% 若输入1x2的行向量, 末位自动补0, 变为1x3的行向量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUT
% pc   :  - 圆心坐标, 行向量 1x3
% r    :  - 半径, 标量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 调用示例1 - 平面上三个点
% [pc1,r1]=points2circle([1,2],[-2,1],[0,-3])
% 调用示例2 - 空间中三个点
% [pc2,r2]=points2circle([1,2,-1],[-2,1,2],[0,-3,-3])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [pc,r]=points2circle(p1,p2,p3)

    % 输入检查
    validateattributes(p1,{'numeric'},{'row'},1);% 行向量
    validateattributes(p2,{'numeric'},{'row'},2);
    validateattributes(p3,{'numeric'},{'row'},3);
    num1=length(p1);num2=length(p2);num3=length(p3);
    if (num1 == num2) && (num2 == num3)
        if num1 == 2
            p1=[p1,0];p2=[p2,0];p3=[p3,0];
        elseif num1 ~= 3
            error('仅支持二维或三维坐标输入');
        end
    else
        error('输入坐标的维数不一致');
    end

    % 共线检查
    temp01=p1-p2;temp02=p3-p2;
    temp03=cross(temp01,temp02);
    temp=(temp03*temp03')/(temp01*temp01')/(temp02*temp02');
    if temp < 10^-6
        error('三点共线, 无法确定圆');
    end

    mat1=[p1,1;p2,1;p3,1];% size = 3x4

    m=+det(mat1(:,2:4));
    n=-det([mat1(:,1),mat1(:,3:4)]);
    p=+det([mat1(:,1:2),mat1(:,4)]);
    q=-det(mat1(:,1:3));

    mat2=[[p1*p1';p2*p2';p3*p3'],mat1;2*q,[-m,-n,-p,0]];% size = 4x5

    A=+det(mat2(:,2:5));
    B=-det([mat2(:,1),mat2(:,3:5)]);
    C=+det([mat2(:,1:2),mat2(:,4:5)]);
    D=-det([mat2(:,1:3),mat2(:,5)]);
    E=+det(mat2(:,1:4));

    pc=-[B,C,D]/2/A;
    r=sqrt(B^2+C^2+D^2-4*A*E)/2/abs(A);

end

Python代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
###########################################################################
# 由圆上三点确定圆心和半径
###########################################################################
# INPUT
# p1   :  - 第一个点坐标, list或者array 1x3
# p2   :  - 第二个点坐标, list或者array 1x3
# p3   :  - 第三个点坐标, list或者array 1x3
# 若输入1x2的行向量, 末位自动补0, 变为1x3的行向量
###########################################################################
# OUTPUT
# pc   :  - 圆心坐标, array 1x3
# r    :  - 半径, 标量
###########################################################################
# 调用示例1 - 平面上三个点
# pc1, r1 = points2circle([1, 2], [-2, 1], [0, -3])
# 调用示例2 - 空间中三个点
# pc2, r2 = points2circle([1, 2, -1], [-2, 1, 2], [0, -3, -3])
###########################################################################
import numpy as np
from numpy.linalg import det

def points2circle(p1, p2, p3):
    p1 = np.array(p1)
    p2 = np.array(p2)
    p3 = np.array(p3)
    num1 = len(p1)
    num2 = len(p2)
    num3 = len(p3)

    # 输入检查
    if (num1 == num2) and (num2 == num3):
        if num1 == 2:
            p1 = np.append(p1, 0)
            p2 = np.append(p2, 0)
            p3 = np.append(p3, 0)
        elif num1 != 3:
            print('\t仅支持二维或三维坐标输入')
            return None
    else:
        print('\t输入坐标的维数不一致')
        return None

    # 共线检查
    temp01 = p1 - p2
    temp02 = p3 - p2
    temp03 = np.cross(temp01, temp02)
    temp = (temp03 @ temp03) / (temp01 @ temp01) / (temp02 @ temp02)
    if temp < 10**-6:
        print('\t三点共线, 无法确定圆')
        return None

    temp1 = np.vstack((p1, p2, p3))
    temp2 = np.ones((3, 1))
    mat1 = np.hstack((temp1, temp2))  # size = 3x4

    m = +det(mat1[:, 1:])
    n = -det(np.delete(mat1, 1, axis=1))
    p = +det(np.delete(mat1, 2, axis=1))
    q = -det(temp1)

    temp3 = np.array([p1 @ p1, p2 @ p2, p3 @ p3]).reshape(3, 1)
    temp4 = np.hstack((temp3, mat1))
    temp5 = np.array([2 * q, -m, -n, -p, 0])
    mat2 = np.vstack((temp4, temp5))  # size = 4x5

    A = +det(mat2[:, 1:])
    B = -det(np.delete(mat2, 1, axis=1))
    C = +det(np.delete(mat2, 2, axis=1))
    D = -det(np.delete(mat2, 3, axis=1))
    E = +det(mat2[:, :-1])

    pc = -np.array([B, C, D]) / 2 / A
    r = np.sqrt(B * B + C * C + D * D - 4 * A * E) / 2 / abs(A)

    return pc, r