更多阅读：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​​y1​y2​y3​​111​∣∣∣∣∣∣∣​(4b)

C = + ∣ x 1 2 + y 1 2 x 1 1 x 2 2 + y 2 2 x 2 1 x 3 2 + y 3 2 x 3 1 ∣ (4c) \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} C=+∣∣∣∣∣∣∣​x12​+y12​x22​+y22​x32​+y32​​x1​x2​x3​​111​∣∣∣∣∣∣∣​(4c)

D = − ∣ x 1 2 + y 1 2 x 1 y 1 x 2 2 + y 2 2 x 2 y 2 x 3 2 + y 3 2 x 3 y 3 ∣ (4d) \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} D=−∣∣∣∣∣∣∣​x12​+y12​x22​+y22​x32​+y32​​x1​x2​x3​​y1​y2​y3​​∣∣∣∣∣∣∣​(4d)

由式 ( 1 b ) (1\mathrm{b}) (1b)可得圆心坐标   ( x c , y c )   ~(x_c,y_c)~  (xc​,yc​) 和半径   r ~r  r，即
{ x c =   − B 2 A y c =   − C 2 A r =   B 2 + C 2 − 4 A D 4 A 2 (5) \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} ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧​xc​yc​r​= −2AB​= −2AC​= 4A2B2+C2−4AD​ ​​(5)

A   A~ A 可能是负数

空间三点

进一步推广到空间中，三个已知点的坐标分别记为   ( x 1 , y 1 , z 1 ) ~(x_1,y_1,z_1)  (x1​,y1​,z1​)、 ( x 2 , y 2 , z 2 ) (x_2,y_2,z_2) (x2​,y2​,z2​)、 ( x 3 , y 3 , z 3 ) (x_3,y_3,z_3) (x3​,y3​,z3​)。

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

1. 圆心   ( x c , y c , z c ) ~(x_c,y_c,z_c)  (xc​,yc​,zc​)
2. 半径   r ~r  r
3. 圆所在平面的法向量   ( m , n , p ) ~(m,n,p)  (m,n,p)

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

将式 ( 6 a ) (6\mathrm{a}) (6a)变形可得球的标准方程，即
( x + B 2 A ) 2 + ( y + C 2 A ) 2 + ( z + D 2 A ) 2 = B 2 + C 2 + D 2 − 4 A E 4 A 2 (6b) \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} (x+2AB​)2+(y+2AC​)2+(z+2AD​)2=4A2B2+C2+D2−4AE​(6b)

将三个已知点代入式 ( 6 a ) (6\mathrm{a}) (6a)，可得关于   A ~A  A、 B B B、 C C C、 D D D、 E   E~ E 的齐次线性方程组，即
[ x 2 + y 2 + z 2 x y z 1 x 1 2 + y 1 2 + z 1 2 x 1 y 1 z 1 1 x 2 2 + y 2 2 + z 2 2 x 2 y 2 z 2 1 x 3 2 + y 3 2 + z 3 2 x 3 y 3 z 3 1 ] 4 × 5 ⋅ [ A B C D E ] = [ 0 0 0 0 ] (7) \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} ⎣⎢⎢⎢⎡​x2+y2+z2x12​+y12​+z12​x22​+y22​+z22​x32​+y32​+z32​​xx1​x2​x3​​yy1​y2​y3​​zz1​z2​z3​​1111​⎦⎥⎥⎥⎤​4×5​⋅⎣⎢⎢⎢⎢⎡​ABCDE​⎦⎥⎥⎥⎥⎤​=⎣⎢⎢⎡​0000​⎦⎥⎥⎤​(7)

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

将三个已知点代入式 ( 8 ) (8) (8)，可得关于   m ~m  m、 n n n、 p p p、 q   q~ q 的齐次线性方程组，即
[ 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 ] ⋅ [ m n p q ] = [ 0 0 0 0 ] (9) \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} ⎣⎢⎢⎡​xx1​x2​x3​​yy1​y2​y3​​zz1​z2​z3​​1111​⎦⎥⎥⎤​⋅⎣⎢⎢⎡​mnpq​⎦⎥⎥⎤​=⎣⎢⎢⎡​0000​⎦⎥⎥⎤​(9)

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

n = − ∣ x 1 z 1 1 x 2 z 2 1 x 3 z 3 1 ∣ (10b) \tag{10b} n=-\begin{vmatrix} x_1 & z_1 & 1 \\ x_2 & z_2 & 1 \\ x_3 & z_3 & 1 \end{vmatrix} n=−∣∣∣∣∣∣​x1​x2​x3​​z1​z2​z3​​111​∣∣∣∣∣∣​(10b)

p = + ∣ x 1 y 1 1 x 2 y 2 1 x 3 y 3 1 ∣ (10c) \tag{10c} p=+\begin{vmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{vmatrix} p=+∣∣∣∣∣∣​x1​x2​x3​​y1​y2​y3​​111​∣∣∣∣∣∣​(10c)

q = − ∣ x 1 y 1 z 1 x 2 y 2 z 2 x 2 y 3 z 3 ∣ (10d) \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} q=−∣∣∣∣∣∣​x1​x2​x2​​y1​y2​y3​​z1​z2​z3​​∣∣∣∣∣∣​(10d)

由式 ( 6 b ) (6\mathrm{b}) (6b)可知
x c = − B 2 A , y c = − C 2 A , z c = − D 2 A (11) \tag{11} x_c=-\frac{B}{2A}, \quad y_c=-\frac{C}{2A}, \quad z_c=-\frac{D}{2A} xc​=−2AB​,yc​=−2AC​,zc​=−2AD​(11)

考虑到圆心   ( x c , y c , z c )   ~(x_c,y_c,z_c)~  (xc​,yc​,zc​) 也在平面   π   ~\pi~  π 上，于是将式 ( 11 ) (11) (11)代入式 ( 8 ) (8) (8)可得
2 q A − m B − n C − p D = 0 (12) \tag{12} 2qA-mB-nC-pD=0 2qA−mB−nC−pD=0(12)

将式 ( 12 ) (12) (12)与式 ( 7 ) (7) (7)结合，可得
[ x 2 + y 2 + z 2 x y z 1 x 1 2 + y 1 2 + z 1 2 x 1 y 1 z 1 1 x 2 2 + y 2 2 + z 2 2 x 2 y 2 z 2 1 x 3 2 + y 3 2 + z 3 2 x 3 y 3 z 3 1 2 q − m − n − p 0 ] 5 × 5 ⋅ [ A B C D E ] = [ 0 0 0 0 0 ] (13) \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} ⎣⎢⎢⎢⎢⎢⎢⎡​x2+y2+z2x12​+y12​+z12​x22​+y22​+z22​x32​+y32​+z32​2q​xx1​x2​x3​−m​yy1​y2​y3​−n​zz1​z2​z3​−p​11110​⎦⎥⎥⎥⎥⎥⎥⎤​5×5​⋅⎣⎢⎢⎢⎢⎡​ABCDE​⎦⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎡​00000​⎦⎥⎥⎥⎥⎤​(13)

这样，在三点不共线的前提下，该齐次线性方程组有非零解，其等价于系数矩阵不满秩。于是将其系数矩阵的行列式展开，并与式 ( 6 ) (6) (6)对比可得五个系数，即
A = + ∣ 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 ∣ (14a) \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} A=+∣∣∣∣∣∣∣∣​x1​x2​x3​−m​y1​y2​y3​−n​z1​z2​z3​−p​1110​∣∣∣∣∣∣∣∣​(14a)

B = − ∣ x 1 2 + y 1 2 + z 1 2 y 1 z 1 1 x 2 2 + y 2 2 + z 2 2 y 2 z 2 1 x 3 2 + y 3 2 + z 3 2 y 3 z 3 1 2 q − n − p 0 ∣ (14b) \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} B=−∣∣∣∣∣∣∣∣∣​x12​+y12​+z12​x22​+y22​+z22​x32​+y32​+z32​2q​y1​y2​y3​−n​z1​z2​z3​−p​1110​∣∣∣∣∣∣∣∣∣​(14b)