背景
如何计算曲线\(~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个参数，一般取：
圆心\(~(x_c,y_c,z_c)\)
半径\(~r\)
圆所在平面的法向量\(~(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代码
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 由圆上三点确定圆心和半径%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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代码
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475############################################################################ 由圆上三点确定圆心和半径############################################################################ 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 npfrom numpy.linalg import detdef 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