目录
- 1.拟合椭圆
- 2.示例代码
爬虫网站自重。
1.拟合椭圆
二次曲线的一般方程为:
A x 2 + B x y + C y 2 + D x + E y + F = 0 Ax^2+Bxy+Cy^2+Dx+Ey+F=0 Ax2+Bxy+Cy2+Dx+Ey+F=0
令:
Δ = B 2 − 4 A C Δ =B^2-4AC Δ=B2−4AC
那么,当 Δ > 0 Δ >0 Δ>0时方程表示双曲线;当 Δ = 0 Δ =0 Δ=0时方程表示抛物线;当 Δ < 0 Δ <0 Δ<0时方程表示椭圆。
即对于椭圆的一般方程:
A x 2 + B x y + C y 2 + D x + E y + F = 0 Ax^2+Bxy+Cy^2+Dx+Ey+F=0 Ax2+Bxy+Cy2+Dx+Ey+F=0
需要有约束条件:
Δ = B 2 − 4 A C < 0 Δ =B^2-4AC<0 Δ=B2−4AC<0
称二元函数 f ( x , y ) f(x,y) f(x,y) 为点 ( x , y ) (x,y) (x,y)到椭圆的代数距离, f ( x , y ) f(x,y) f(x,y)的表达式为:
f ( x , y ) = A x 2 + B x y + C y 2 + D x + E y + F f(x,y)=Ax^2+Bxy+Cy^2+Dx+Ey+F f(x,y)=Ax2+Bxy+Cy2+Dx+Ey+F
如果令向量 a ⃗ \vec{a} a为:
a ⃗ = [ A , B , C , D , E , F ] T \vec{a}=[A,B,C,D,E,F]^T a=[A,B,C,D,E,F]T
令向量 x ⃗ \vec{x} x为:
x ⃗ = [ x 2 , x y , y 2 , x , y , 1 ] \vec{x}=[x^2,xy,y^2,x,y,1] x=[x2,xy,y2,x,y,1]
那么椭圆的一般方程就可以改写为:
a ⃗ ⋅ x ⃗ = 0 \vec{a}\cdot\vec{x}=0 a⋅x=0
即代数距离为: a ⃗ ⋅ x ⃗ \vec{a}\cdot\vec{x} a⋅x。设有 N N N个数据点: ( x i , y i ) , i = 1 , 2 , 3 … N (x_i,y_i),i=1,2,3…N (xi,yi),i=1,2,3…N,那么如果要求最佳拟合的椭圆,便需要求 ∑ i = 1 N ( a ⃗ ⋅ X i ⃗ ) 2 \sum_{i=1}^N(\vec{a}\cdot\vec{X_i})^2 ∑i=1N(a⋅Xi)2的最小值。这可以运用最小二乘法求解,但所得到的结果是圆锥曲线而不是一定椭圆。
所以,如果要获得椭圆解,就需要加入约束条件:
Δ = B 2 − 4 A C < 0 Δ =B^2-4AC<0 Δ=B2−4AC<0
一般地,在一个合适的尺度下,上述形式为不等式的约束条件可以转化为形式为等式的约束条件,即有:
B 2 − 4 A C = − 1 B^2-4AC=-1 B2−4AC=−1
那么,如果定义一个 N × 6 N×6 N×6的矩阵 D D D:
D = [ x 1 2 x 1 y 1 y 1 2 x 1 y 1 1 x 2 2 x 2 y 2 y 2 2 x 2 y 2 1 x 3 2 x 3 y 3 y 3 2 x 3 y 3 1 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ x n 2 x n y n y n 2 x n y n 1 ] D=\left[ \begin{matrix} x_1^2& x_1y_1 & y_1^2 &x_1&y_1&1 \\ x_2^2& x_2y_2 & y_2^2 &x_2&y_2&1 \\ x_3^2& x_3y_3 & y_3^2 &x_3&y_3&1 \\ \vdots & \vdots & \vdots & \vdots& \vdots& \vdots \\ x_n^2& x_ny_n & y_n^2 &x_n&y_n&1 \\ \end{matrix} \right] D= x12x22x32⋮xn2x1y1x2y2x3y3⋮xnyny12y22y32⋮yn2x1x2x3⋮xny1y2y3⋮yn111⋮1
定义一个 6 × 6 6×6 6×6的矩阵 C C C:
C = [ 0 0 2 0 0 0 0 − 1 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ] C=\left[ \begin{matrix} 0& 0 & 2 & 0 & 0 & 0\\ 0 & -1 & 0 & 0 & 0 &0 \\ 2& 0 & 0 & 0 & 0 & 0 \\ 0& 0 & 0 & 0 & 0 & 0 \\ 0& 0 & 0 & 0 & 0 & 0 \\ 0& 0 & 0 & 0 & 0 & 0 \\ \end{matrix} \right] C= 0020000−10000200000000000000000000000
那么,约束条件就可以写为:
a ˉ T C a ˉ = 1 \bar{a}^TC\bar{a}=1 aˉTCaˉ=1
在此条件下,求下式的最小值:
∣ ∣ D a ˉ ∣ ∣ 2 = a ˉ T D T D a ˉ ||D\bar{a}||^2=\bar{a}^TD^TD\bar{a} ∣∣Daˉ∣∣2=aˉTDTDaˉ
构造拉格朗日函数为:
L ( D , λ ) = a ˉ T D T a ˉ − λ ( a ˉ T C a ˉ − 1 ) L(D,\lambda)=\bar{a}^TD^T\bar{a}-\lambda(\bar{a}^TC\bar{a}-1) L(D,λ)=aˉTDTaˉ−λ(aˉTCaˉ−1)
为求极值,可令偏导数为0:
δ L ( D , λ ) δ a ˉ = D T D a ˉ − λ C a ˉ = 0 \frac{\delta L(D,\lambda)}{\delta \bar{a}}=D^TD\bar{a}-\lambda C\bar{a}=0 δaˉδL(D,λ)=DTDaˉ−λCaˉ=0
于是:
D T D a ˉ = λ C a ˉ D^TD\bar{a}=\lambda C\bar{a} DTDaˉ=λCaˉ
如果令:
S = D T D S=D^TD S=DTD
那么,等式就可以化为:
S a ⃗ = λ C a ⃗ S\vec{a}=\lambda C\vec{a} Sa=λCa
如果 S S S可逆,并且令 X = S − 1 C X=S^{-1}C X=S−1C,那么可以将其转化为求取矩阵的特征值的形式:
X a ⃗ = 1 λ a ⃗ X\vec{a}=\frac{1}{\lambda}\vec{a} Xa=λ1a
可以知道,上面的矩阵 X X X为6 阶方阵,其特征多项式为6 阶多项式。
事实上,可以证明,问题的结果对应着 X X X的特征值的绝对值最大的特征向量。
在实际操作中,是可以采用求极限的算法进行求解的。
如果矩阵 X X X有 n n n个不同的特征值: λ 1 , λ 2 , λ 3 , … , λ n λ_1,λ_2,λ_3,…,λ_n λ1,λ2,λ3,…,λn其对应的特征向量分别为:
a ⃗ 1 , a ⃗ 2 , a ⃗ 3 , … , a ⃗ n \vec{a}_1,\vec{a}_2,\vec{a}_3,…,\vec{a}_n a1,a2,a3,…,an
可以知道,这些特征向量线性无关,于是, n n n维向量 x ⃗ \vec{x} x可以表示为:
x ⃗ = k 1 a ⃗ 1 + k 2 a ⃗ 2 + k 3 a ⃗ 3 + . . . + k n a ⃗ n \vec{x}=k_1\vec{a}_1+k_2\vec{a}_2+k_3\vec{a}_3+...+k_n\vec{a}_n x=k1a1+k2a2+k3a3+...+knan
那么,可以得到:
X m x ⃗ = k 1 m a ⃗ 1 + k 2 m a ⃗ 2 + k 3 m a ⃗ 3 + . . . + k n m a ⃗ n X^m\vec{x}=k_1^m\vec{a}_1+k_2^m\vec{a}_2+k_3^m\vec{a}_3+...+k_n^m\vec{a}_n Xmx=k1ma1+k2ma2+k3ma3+...+knman
令 k k k为 k 1 , k 2 , k 3 … k n k_1,k_2,k_3…k_n k1,k2,k3…kn 中的绝对值最大的值,其对应的特征向量为 a ⃗ \vec{a} a,即:
k = m a x ( ∣ k 1 ∣ , ∣ k 2 ∣ , ∣ k 3 ∣ , … , ∣ k n ∣ ) k=max({|k_1|,|k_2|,|k_3|,…,|k_n|}) k=max(∣k1∣,∣k2∣,∣k3∣,…,∣kn∣)
那么,当 m → + ∞ m → + ∞ m→+∞时,有 X m m ⃗ → k m a ⃗ X^m\vec{m}→k^m\vec{a} Xmm→kma。所得的 a ⃗ \vec{a} a即为所求的椭圆的参数所组成的向量。
2.示例代码
在这里插入代码片