椭圆

椭圆在高中数学里就开始提到，都是从标准方程开始如： x 2 a 2 + y 2 b 2 = 1 ( a > b > 0 ) \frac{x^2}{a^2}+\frac{y^2}{b^2}=1(a>b>0) a2x2​+b2y2​=1(a>b>0) 其中 a a a 为椭圆的长轴， b b b 为椭圆短轴，通常默认椭圆圆心 ( x c e n t e r , y c e n t e r ) (x_{center},y_{center}) (xcenter​,ycenter​) 为 0 0 0，以及偏移角(这里指长轴相对于 x x x 轴的方位角) θ o f f s e t \theta_{offset} θoffset​ 为 0 0 0 ，如长轴为 1 1 1 ，短轴为 0.5 0.5 0.5 的标准椭圆可以这样画出

如果有一个椭圆的圆心不在坐标原点，而且椭圆的偏移角不为 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

也可以表示为如下的非齐次方程 (后续推导采用) x 2 + A x y + B y 2 + C x + D y + E = 0 x^2+Axy+By^2+Cx+Dy+E=0 x2+Axy+By2+Cx+Dy+E=0

进行平移和旋转变化到标准形式 x ′ ′ 2 a 2 + y ′ ′ 2 b 2 = 1 ( a > b > 0 ) \frac{x''^2}{a^2}+\frac{y''^2}{b^2}=1(a>b>0) a2x′′2​+b2y′′2​=1(a>b>0)

补充：假设椭圆变换到标准椭圆的过程，是先按向量 P = ( x p , y p ) P = ( x_p , y_p) P=(xp​,yp​) 平移到坐标原点得到蓝色椭圆，再旋转 θ r \theta_r θr​ 角度(逆时针)得到红色标准椭圆。

可以表达成 [ x ′ y ′ ] = [ x y ] + [ x p y p ] \begin{gathered}\begin{bmatrix} x' \\ y' \end{bmatrix}=\begin{bmatrix} x \\ y\end{bmatrix}+\begin{bmatrix} x_p \\ y_p\end{bmatrix}\end{gathered} [x′y′​]=[xy​]+[xp​yp​​]​ [ x ′ ′ y ′ ′ ] = [ c o s θ r − s i n θ r s i n θ r c o s θ r ] [ x ′ y ′ ] \begin{gathered}\begin{bmatrix} x'' \\ y'' \end{bmatrix}=\begin{bmatrix} cos\theta_r & -sin\theta_r \\ sin\theta_r & cos\theta_r \end{bmatrix} \begin{bmatrix} x' \\ y'\end{bmatrix}\end{gathered} [x′′y′′​]=[cosθr​sinθr​​−sinθr​cosθr​​][x′y′​]​

第一步，将平移变换代入非齐次方程 [ x y ] = [ x ′ y ′ ] − [ x p y p ] \begin{gathered}\begin{bmatrix} x \\ y \end{bmatrix}=\begin{bmatrix} x' \\ y'\end{bmatrix}-\begin{bmatrix} x_p \\ y_p\end{bmatrix}\end{gathered} [xy​]=[x′y′​]−[xp​yp​​]​ ( x ′ − x p ) 2 + A ( x ′ − x p ) ( y ′ − y p ) + B ( y ′ − y p ) 2 + C ( x ′ − x p ) + D ( y ′ − y p ) + E = 0 (x'-x_p)^2+A(x'-x_p)(y'-y_p)+B(y'-y_p)^2+C(x'-x_p)+D(y'-y_p)+E=0 (x′−xp​)2+A(x′−xp​)(y′−yp​)+B(y′−yp​)2+C(x′−xp​)+D(y′−yp​)+E=0

表达为矩阵形式为 [ 1 A B − 2 x p − A y p + C − A x p − 2 B y p + D x p 2 + A x p y p + B y p 2 − C x p − D y p + E ] T [ x ′ 2 x ′ y ′ y ′ 2 x ′ y ′ 1 ] = 0 \begin{gathered}\begin{bmatrix} 1 \\ A \\ B \\ -2x_p-Ay_p+C \\ -Ax_p-2By_p+D \\ x_p^2+Ax_py_p+By_p^2-Cx_p-Dy_p+E \end{bmatrix}^T\begin{bmatrix} x'^2 \\ x'y' \\ y'^2 \\ x' \\ y' \\ 1\end{bmatrix}=0\end{gathered} ⎣⎢⎢⎢⎢⎢⎢⎡​1AB−2xp​−Ayp​+C−Axp​−2Byp​+Dxp2​+Axp​yp​+Byp2​−Cxp​−Dyp​+E​⎦⎥⎥⎥⎥⎥⎥⎤​T⎣⎢⎢⎢⎢⎢⎢⎡​x′2x′y′y′2x′y′1​⎦⎥⎥⎥⎥⎥⎥⎤​=0​

经过平移椭圆的圆心落在坐标原点，此时 x ′ x' x′ 和 y ′ y' y′ 的系数均为0，因此可以得到如下关系 C = 2 x p + A y p C = 2x_p+Ay_p C=2xp​+Ayp​ D = A x p + 2 B y p D=Ax_p+2By_p D=Axp​+2Byp​可解得 x p = A D − 2 B C A 2 − 4 B x_p = \frac{AD-2BC}{A^2-4B} xp​=A2−4BAD−2BC​ y p = A C − 2 D A 2 − 4 B y_p=\frac{AC-2D}{A^2-4B} yp​=A2−4BAC−2D​

椭圆的圆心 ( x c , y c ) (x_c,y_c) (xc​,yc​) 可表示为 x c e n t e r = − x p x_{center} = -x_p xcenter​=−xp​ y c e n t e r = − y p y_{center}=-y_p ycenter​=−yp​

此时矩阵化简为 [ 1 A B 0 0 − x p 2 − A x p y p − B y p 2 + E ] T [ x ′ 2 x ′ y ′ y ′ 2 x ′ y ′ 1 ] = 0 \begin{gathered}\begin{bmatrix} 1 \\ A \\ B \\ 0 \\ 0 \\ -x_p^2-Ax_py_p-By_p^2+E \end{bmatrix}^T\begin{bmatrix} x'^2 \\ x'y' \\ y'^2 \\ x' \\ y' \\ 1\end{bmatrix}=0\end{gathered} ⎣⎢⎢⎢⎢⎢⎢⎡​1AB00−xp2​−Axp​yp​−Byp2​+E​⎦⎥⎥⎥⎥⎥⎥⎤​T⎣⎢⎢⎢⎢⎢⎢⎡​x′2x′y′y′2x′y′1​⎦⎥⎥⎥⎥⎥⎥⎤​=0​

第二步，将旋转变换代入上式 [ x ′ y ′ ] = [ c o s θ r s i n θ r − s i n θ r c o s θ r ] [ x ′ ′ y ′ ′ ] \begin{gathered}\begin{bmatrix} x' \\ y' \end{bmatrix}=\begin{bmatrix} cos\theta_r & sin\theta_r \\ -sin\theta_r & cos\theta_r \end{bmatrix} \begin{bmatrix} x'' \\ y''\end{bmatrix}\end{gathered} [x′y′​]=[cosθr​−sinθr​​sinθr​cosθr​​][x′′y′′​]​ ( x ′ ′ c o s θ r + y ′ ′ s i n θ r ) 2 + A ( x ′ ′ c o s θ r + y ′ ′ s i n θ r ) ( − x ′ ′ s i n θ r + y ′ ′ c o s θ r ) + B ( − x ′ ′ s i n θ r + y ′ ′ c o s θ r ) 2 − x p 2 − A x p y p − B y p 2 + E = 0 (x''cos\theta_r+y''sin\theta_r)^2+A(x''cos\theta_r+y''sin\theta_r)(-x''sin\theta_r+y''cos\theta_r)\\+B(-x''sin\theta_r+y''cos\theta_r)^2-x_p^2-Ax_py_p-By_p^2+E=0 (x′′cosθr​+y′′sinθr​)2+A(x′′cosθr​+y′′sinθr​)(−x′′sinθr​+y′′cosθr​)+B(−x′′sinθr​+y′′cosθr​)2−xp2​−Axp​yp​−Byp2​+E=0

表达为矩阵形式为 [ c o s 2 θ r − A s i n θ r c o s θ r + B s i n 2 θ r 2 s i n θ r c o s θ r + A ( c o s 2 θ r − s i n 2 θ r ) − 2 B s i n θ r c o s θ r s i n 2 θ r + A s i n θ r c o s θ r + B c o s 2 θ r 0 0 − x p 2 − A x p y p − B y p 2 + E ] T [ x ′ ′ 2 x ′ ′ y ′ ′ y ′ ′ 2 x ′ ′ y ′ ′ 1 ] = 0 \begin{gathered}\begin{bmatrix} cos^2\theta_r-Asin\theta_r cos\theta_r+Bsin^2\theta_r \\ 2sin\theta_r cos\theta_r+A(cos^2\theta_r-sin^2\theta_r)-2Bsin\theta_r cos\theta_r \\ sin^2\theta_r +Asin\theta_r cos\theta_r +Bcos^2\theta_r \\ 0 \\ 0 \\ -x_p^2-Ax_py_p-By_p^2+E \end{bmatrix}^T\begin{bmatrix} x''^2 \\ x''y'' \\ y''^2 \\ x'' \\ y'' \\ 1\end{bmatrix}=0\end{gathered} ⎣⎢⎢⎢⎢⎢⎢⎡​cos2θr​−Asinθr​cosθr​+Bsin2θr​2sinθr​cosθr​+A(cos2θr​−sin2θr​)−2Bsinθr​cosθr​sin2θr​+Asinθr​cosθr​+Bcos2θr​00−xp2​−Axp​yp​−Byp2​+E​⎦⎥⎥⎥⎥⎥⎥⎤​T⎣⎢⎢⎢⎢⎢⎢⎡​x′′2x′′y′′y′′2x′′y′′1​⎦⎥⎥⎥⎥⎥⎥⎤​=0​

经过旋转椭圆的偏移角为0，此时 x ′ ′ y ′ ′ x''y'' x′′y′′ 系数为0，因此可以得到如下关系 2 s i n θ r c o s θ r + A ( c o s 2 θ r − s i n 2 θ r ) − 2 B s i n θ r c o s θ r = 0 2sin\theta_r cos\theta_r+A(cos^2\theta_r-sin^2\theta_r)-2Bsin\theta_r cos\theta_r = 0 2sinθr​cosθr​+A(cos2θr​−sin2θr​)−2Bsinθr​cosθr​=0

t a n 2 θ r = 2 s i n θ r c o s θ r c o s 2 θ r − s i n 2 θ r = A B − 1 tan2\theta_r = \frac{2sin\theta_r cos\theta_r }{cos^2\theta_r-sin^2\theta_r}=\frac{A}{B-1} tan2θr​=cos2θr​−sin2θr​2sinθr​cosθr​​=B−1A​ θ r = 1 2 a r c t a n A B − 1 \theta_r = \frac{1}{2}arctan\frac{A}{B-1} θr​=21​arctanB−1A​ 椭圆的偏移角 θ o f f s e t \theta_{offset} θoffset​ 可表示为 θ o f f s e t = − θ r \theta_{offset}=-\theta_r θoffset​=−θr​

此时矩阵化简为 [ 1 0 s i n 2 θ r + A s i n θ r c o s θ r + B c o s 2 θ r c o s 2 θ r − A s i n θ r c o s θ r + B s i n 2 θ r 0 0 − x p 2 − A x p y p − B y p 2 + E c o s 2 θ r − A s i n θ r c o s θ r + B s i n 2 θ r ] T [ x ′ ′ 2 x ′ ′ y ′ ′ y ′ ′ 2 x ′ ′ y ′ ′ 1 ] = 0 \begin{gathered}\begin{bmatrix} 1 \\ 0 \\ \frac{sin^2\theta_r +Asin\theta_r cos\theta_r +Bcos^2\theta_r}{cos^2\theta_r-Asin\theta_r cos\theta_r+Bsin^2\theta_r} \\ 0 \\ 0 \\ \frac{-x_p^2-Ax_py_p-By_p^2+E }{cos^2\theta_r-Asin\theta_r cos\theta_r+Bsin^2\theta_r}\end{bmatrix}^T\begin{bmatrix} x''^2 \\ x''y'' \\ y''^2 \\ x'' \\ y'' \\ 1\end{bmatrix}=0\end{gathered} ⎣⎢⎢⎢⎢⎢⎢⎢⎡​10cos2θr​−Asinθr​cosθr​+Bsin2θr​sin2θr​+Asinθr​cosθr​+Bcos2θr​​00cos2θr​−Asinθr​cosθr​+Bsin2θr​−xp2​−Axp​yp​−Byp2​+E​​⎦⎥⎥⎥⎥⎥⎥⎥⎤​T⎣⎢⎢⎢⎢⎢⎢⎡​x′′2x′′y′′y′′2x′′y′′1​⎦⎥⎥⎥⎥⎥⎥⎤​=0​

对比标准方程可以得到椭圆的长短轴 a = x p 2 + A x p y p + B y p 2 − E c o s 2 θ r − A s i n θ r c o s θ r + B s i n 2 θ r a = \sqrt{\frac{x_p^2+Ax_py_p+By_p^2-E }{cos^2\theta_r-Asin\theta_r cos\theta_r+Bsin^2\theta_r}} a=cos2θr​−Asinθr​cosθr​+Bsin2θr​xp2​+Axp​yp​+Byp2​−E​ ​ b = x p 2 + A x p y p + B y p 2 − E s i n 2 θ r + A s i n θ r c o s θ r + B c o s 2 θ r b = \sqrt{\frac{x_p^2+Ax_py_p+By_p^2-E }{sin^2\theta_r +Asin\theta_r cos\theta_r +Bcos^2\theta_r}} b=sin2θr​+Asinθr​cosθr​+Bcos2θr​xp2​+Axp​yp​+Byp2​−E​ ​

而反过来，已知椭圆的 a a a ， b b b ， θ o f f s e t \theta_{offset} θoffset​ ， ( x c e n t e r , y c e n t e r ) (x_{center},y_{center}) (xcenter​,ycenter​)，如何表达成一般式。

将之前的旋转平移变换过程 [ x ′ y ′ ] = [ x y ] + [ x p y p ] \begin{gathered}\begin{bmatrix} x' \\ y' \end{bmatrix}=\begin{bmatrix} x \\ y\end{bmatrix}+\begin{bmatrix} x_p \\ y_p\end{bmatrix}\end{gathered} [x′y′​]=[xy​]+[xp​yp​​]​ [ x ′ ′ y ′ ′ ] = [ c o s θ r − s i n θ r s i n θ r c o s θ r ] [ x ′ y ′ ] \begin{gathered}\begin{bmatrix} x'' \\ y'' \end{bmatrix}=\begin{bmatrix} cos\theta_r & -sin\theta_r \\ sin\theta_r & cos\theta_r \end{bmatrix} \begin{bmatrix} x' \\ y'\end{bmatrix}\end{gathered} [x′′y′′​]=[cosθr​sinθr​​−sinθr​cosθr​​][x′y′​]​

代入标准方程 x ′ ′ 2 a 2 + y ′ ′ 2 b 2 = 1 ( a > b > 0 ) \frac{x''^2}{a^2}+\frac{y''^2}{b^2}=1(a>b>0) a2x′′2​+b2y′′2​=1(a>b>0) 即 ( ( x + x p ) c o s θ r − ( y + y p ) s i n θ r ) 2 a 2 + ( ( x + x p ) s i n θ r + ( y + y p ) c o s θ r ) 2 b 2 = 1 \frac{((x+x_p)cos\theta_r-(y+y_p)sin\theta_r)^2}{a^2}+\frac{((x+x_p)sin\theta_r+(y+y_p)cos\theta_r)^2}{b^2}=1 a2((x+xp​)cosθr​−(y+yp​)sinθr​)2​+b2((x+xp​)sinθr​+(y+yp​)cosθr​)2​=1

[ b 2 c o s 2 θ r + a 2 s i n 2 θ r 2 ( a 2 − b 2 ) s i n θ r c o s θ r b 2 s i n 2 θ r + a 2 c o s 2 θ r 2 b 2 ( x p c o s 2 θ r − y p s i n θ r c o s θ r ) + 2 a 2 ( x p s i n 2 θ r + y p s i n θ r c o s θ r ) 2 b 2 ( y p s i n 2 θ r − x p s i n θ r c o s θ r ) + 2 a 2 ( y p c o s 2 θ r + x p s i n θ r c o s θ r ) b 2 ( x p 2 c o s 2 θ r + y p 2 s i n 2 θ r − 2 x p y p s i n θ r c o s θ r ) + a 2 ( x p 2 s i n 2 θ r + y p 2 c o s 2 θ r + 2 x p y p s i n θ r c o s θ r ) − a 2 b 2 ] T [ x 2 x y y 2 x y 1 ] = 0 \begin{gathered}\begin{bmatrix} b^2cos^2\theta_r+a^2sin^2\theta_r \\ 2(a^2-b^2)sin\theta_rcos\theta_r \\ b^2sin^2\theta_r+a^2cos^2\theta_r \\ 2b^2(x_pcos^2\theta_r-y_psin\theta_rcos\theta_r)+2a^2(x_psin^2\theta_r+y_psin\theta_rcos\theta_r) \\ 2b^2(y_psin^2\theta_r-x_psin\theta_rcos\theta_r)+2a^2(y_pcos^2\theta_r+x_psin\theta_rcos\theta_r) \\ b^2(x_p^2cos^2\theta_r+y_p^2sin^2\theta_r-2x_py_psin\theta_rcos\theta_r)+a^2(x_p^2sin^2\theta_r+y_p^2cos^2\theta_r+2x_py_psin\theta_rcos\theta_r)-a^2b^2 \end{bmatrix}^T\begin{bmatrix} x^2 \\ xy \\ y^2 \\ x \\ y \\ 1\end{bmatrix}=0\end{gathered} ⎣⎢⎢⎢⎢⎢⎢⎡​b2cos2θr​+a2sin2θr​2(a2−b2)sinθr​cosθr​b2sin2θr​+a2cos2θr​2b2(xp​cos2θr​−yp​sinθr​cosθr​)+2a2(xp​sin2θr​+yp​sinθr​cosθr​)2b2(yp​sin2θr​−xp​sinθr​cosθr​)+2a2(yp​cos2θr​+xp​sinθr​cosθr​)b2(xp2​cos2θr​+yp2​sin2θr​−2xp​yp​sinθr​cosθr​)+a2(xp2​sin2θr​+yp2​cos2θr​+2xp​yp​sinθr​cosθr​)−a2b2​⎦⎥⎥⎥⎥⎥⎥⎤​T⎣⎢⎢⎢⎢⎢⎢⎡​x2xyy2xy1​⎦⎥⎥⎥⎥⎥⎥⎤​=0​

[ 1 2 ( a 2 − b 2 ) s i n θ r c o s θ r b 2 c o s 2 θ r + a 2 s i n 2 θ r b 2 s i n 2 θ r + a 2 c o s 2 θ r b 2 c o s 2 θ r + a 2 s i n 2 θ r 2 b 2 ( x p c o s 2 θ r − y p s i n θ r c o s θ r ) + 2 a 2 ( x p s i n 2 θ r + y p s i n θ r c o s θ r ) b 2 c o s 2 θ r + a 2 s i n 2 θ r 2 b 2 ( y p s i n 2 θ r − x p s i n θ r c o s θ r ) + 2 a 2 ( y p c o s 2 θ r + x p s i n θ r c o s θ r ) b 2 c o s 2 θ r + a 2 s i n 2 θ r b 2 ( x p 2 c o s 2 θ r + y p 2 s i n 2 θ r − 2 x p y p s i n θ r c o s θ r ) + a 2 ( x p 2 s i n 2 θ r + y p 2 c o s 2 θ r + 2 x p y p s i n θ r c o s θ r ) − a 2 b 2 b 2 c o s 2 θ r + a 2 s i n 2 θ r ] T [ x 2 x y y 2 x y 1 ] = 0 \begin{gathered}\begin{bmatrix} 1 \\ \frac{2(a^2-b^2)sin\theta_rcos\theta_r}{b^2cos^2\theta_r+a^2sin^2\theta_r} \\ \frac{b^2sin^2\theta_r+a^2cos^2\theta_r}{b^2cos^2\theta_r+a^2sin^2\theta_r} \\ \frac{2b^2(x_pcos^2\theta_r-y_psin\theta_rcos\theta_r)+2a^2(x_psin^2\theta_r+y_psin\theta_rcos\theta_r)}{b^2cos^2\theta_r+a^2sin^2\theta_r} \\ \frac{2b^2(y_psin^2\theta_r-x_psin\theta_rcos\theta_r)+2a^2(y_pcos^2\theta_r+x_psin\theta_rcos\theta_r)}{b^2cos^2\theta_r+a^2sin^2\theta_r} \\ \frac{b^2(x_p^2cos^2\theta_r+y_p^2sin^2\theta_r-2x_py_psin\theta_rcos\theta_r)+a^2(x_p^2sin^2\theta_r+y_p^2cos^2\theta_r+2x_py_psin\theta_rcos\theta_r)-a^2b^2}{b^2cos^2\theta_r+a^2sin^2\theta_r} \end{bmatrix}^T\begin{bmatrix} x^2 \\ xy \\ y^2 \\ x \\ y \\ 1\end{bmatrix}=0\end{gathered} ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​1b2cos2θr​+a2sin2θr​2(a2−b2)sinθr​cosθr​​b2cos2θr​+a2sin2θr​b2sin2θr​+a2cos2θr​​b2cos2θr​+a2sin2θr​2b2(xp​cos2θr​−yp​sinθr​cosθr​)+2a2(xp​sin2θr​+yp​sinθr​cosθr​)​b2cos2θr​+a2sin2θr​2b2(yp​sin2θr​−xp​sinθr​cosθr​)+2a2(yp​cos2θr​+xp​sinθr​cosθr​)​b2cos2θr​+a2sin2θr​b2(xp2​cos2θr​+yp2​sin2θr​−2xp​yp​sinθr​cosθr​)+a2(xp2​sin2θr​+yp2​cos2θr​+2xp​yp​sinθr​cosθr​)−a2b2​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​T⎣⎢⎢⎢⎢⎢⎢⎡​x2xyy2xy1​⎦⎥⎥⎥⎥⎥⎥⎤​=0​

[ 1 2 ( a 2 − b 2 ) s i n θ r c o s θ r b 2 c o s 2 θ r + a 2 s i n 2 θ r b 2 s i n 2 θ r + a 2 c o s 2 θ r b 2 c o s 2 θ r + a 2 s i n 2 θ r 2 x p + y p A 2 y p B + x p A x p 2 + y p 2 B + x p y p A − a 2 b 2 b 2 c o s 2 θ r + a 2 s i n 2 θ r ] T [ x 2 x y y 2 x y 1 ] = 0 \begin{gathered}\begin{bmatrix} 1 \\ \frac{2(a^2-b^2)sin\theta_rcos\theta_r}{b^2cos^2\theta_r+a^2sin^2\theta_r} \\ \frac{b^2sin^2\theta_r+a^2cos^2\theta_r}{b^2cos^2\theta_r+a^2sin^2\theta_r} \\ 2x_p+y_pA \\ 2y_pB+x_pA \\ x_p^2+y_p^2B+x_py_pA-\frac{a^2b^2}{b^2cos^2\theta_r+a^2sin^2\theta_r} \end{bmatrix}^T\begin{bmatrix} x^2 \\ xy \\ y^2 \\ x \\ y \\ 1\end{bmatrix}=0\end{gathered} ⎣⎢⎢⎢⎢⎢⎢⎢⎡​1b2cos2θr​+a2sin2θr​2(a2−b2)sinθr​cosθr​​b2cos2θr​+a2sin2θr​b2sin2θr​+a2cos2θr​​2xp​+yp​A2yp​B+xp​Axp2​+yp2​B+xp​yp​A−b2cos2θr​+a2sin2θr​a2b2​​⎦⎥⎥⎥⎥⎥⎥⎥⎤​T⎣⎢⎢⎢⎢⎢⎢⎡​x2xyy2xy1​⎦⎥⎥⎥⎥⎥⎥⎤​=0​

A = 2 ( a 2 − b 2 ) s i n θ r c o s θ r b 2 c o s 2 θ r + a 2 s i n 2 θ r A = \frac{2(a^2-b^2)sin\theta_rcos\theta_r}{b^2cos^2\theta_r+a^2sin^2\theta_r} A=b2cos2θr​+a2sin2θr​2(a2−b2)sinθr​cosθr​​

B = b 2 s i n 2 θ r + a 2 c o s 2 θ r b 2 c o s 2 θ r + a 2 s i n 2 θ r B = \frac{b^2sin^2\theta_r+a^2cos^2\theta_r}{b^2cos^2\theta_r+a^2sin^2\theta_r} B=b2cos2θr​+a2sin2θr​b2sin2θr​+a2cos2θr​​

C = 2 x p + y p A C = 2x_p+y_pA C=2xp​+yp​A

D = 2 y p B + x p A D = 2y_pB+x_pA D=2yp​B+xp​A

E = x p 2 + y p 2 B + x p y p A − a 2 b 2 b 2 c o s 2 θ r + a 2 s i n 2 θ r E = x_p^2+y_p^2B+x_py_pA-\frac{a^2b^2}{b^2cos^2\theta_r+a^2sin^2\theta_r} E=xp2​+yp2​B+xp​yp​A−b2cos2θr​+a2sin2θr​a2b2​

到这里，椭圆的一般方程和标准方程之间的系数转换关系已经全部得出，这里需要注意平移向量和椭圆圆心符号相反，旋转角度和椭圆偏移角符号相反。

椭圆的一般方程 x 2 + A x y + B y 2 + C x + D y + E = 0 x^2+Axy+By^2+Cx+Dy+E=0 x2+Axy+By2+Cx+Dy+E=0

设现在采集到了多个测量点 P i ( x i , y i ) P_i(x_i,y_i) Pi​(xi​,yi​) ，根据最小二乘的原理，所拟合的目标函数为： F ( A , B , C , D , E ) = ∑ i = 1 n ( x i 2 + x i y i A + y i 2 B + x i C + y i D + E ) 2 F(A,B,C,D,E) = \sum_{i=1}^n(x_i^2+x_iy_iA+y_i^2B+x_iC+y_iD+E)^2 F(A,B,C,D,E)=i=1∑n​(xi2​+xi​yi​A+yi2​B+xi​C+yi​D+E)2

为了使 F F F 最小，需使得 F F F 的各项偏导为0，即 ∂ F ∂ A = ∂ F ∂ B = ∂ F ∂ C = ∂ F ∂ D = ∂ F ∂ E = 0 \frac{\partial F}{\partial A} =\frac{\partial F}{\partial B} =\frac{\partial F}{\partial C}=\frac{\partial F}{\partial D}=\frac{\partial F}{\partial E}=0 ∂A∂F​=∂B∂F​=∂C∂F​=∂D∂F​=∂E∂F​=0

可以得到方程： [ ∑ x i 2 y i 2 ∑ x i y i 3 ∑ x i 2 y i ∑ x i y i 2 ∑ x i y i ∑ x i y i 3 ∑ y i 4 ∑ x i y i 2 ∑ y i 3 ∑ y i 2 ∑ x i 2 y i ∑ x i y i 2 ∑ x i 2 ∑ x i y i ∑ x i ∑ x i y i 2 ∑ y i 3 ∑ x i y i ∑ y i 2 ∑ y i ∑ x i y i ∑ y i 2 ∑ x i ∑ y i N ] [ A B C D E ] = − [ ∑ x i 3 y i ∑ x i 2 y i 2 ∑ x i 3 ∑ x i 2 y i ∑ x i 2 ] \begin{gathered}\begin{bmatrix} \sum{x_i^2y_i^2}&\sum{x_iy_i^3}&\sum{x_i^2y_i}&\sum{x_iy_i^2}&\sum{x_iy_i}\\\sum{x_iy_i^3}&\sum{y_i^4}&\sum{x_iy_i^2}&\sum{y_i^3}&\sum{y_i^2}\\\sum{x_i^2y_i}&\sum{x_iy_i^2}&\sum{x_i^2}&\sum{x_iy_i}&\sum{x_i}\\\sum{x_iy_i^2}&\sum{y_i^3}&\sum{x_iy_i}&\sum{y_i^2}&\sum{y_i}\\\sum{x_iy_i}&\sum{y_i^2}&\sum{x_i}&\sum{y_i}&N\end{bmatrix}\begin{bmatrix} A\\B\\C\\D\\E\end{bmatrix}=-\begin{bmatrix} \sum{x_i^3y_i}\\\sum{x_i^2y_i^2}\\\sum{x_i^3}\\\sum{x_i^2y_i}\\\sum{x_i^2}\end{bmatrix}\end{gathered} ⎣⎢⎢⎢⎢⎡​∑xi2​yi2​∑xi​yi3​∑xi2​yi​∑xi​yi2​∑xi​yi​​∑xi​yi3​∑yi4​∑xi​yi2​∑yi3​∑yi2​​∑xi2​yi​∑xi​yi2​∑xi2​∑xi​yi​∑xi​​∑xi​yi2​∑yi3​∑xi​yi​∑yi2​∑yi​​∑xi​yi​∑yi2​∑xi​∑yi​N​⎦⎥⎥⎥⎥⎤​⎣⎢⎢⎢⎢⎡​ABCDE​⎦⎥⎥⎥⎥⎤​=−⎣⎢⎢⎢⎢⎡​∑xi3​yi​∑xi2​yi2​∑xi3​∑xi2​yi​∑xi2​​⎦⎥⎥⎥⎥⎤​​

可以表示为 M 1 [ A B C D E ] = M 2 \begin{gathered}M_1\begin{bmatrix} A\\B\\C\\D\\E\end{bmatrix}=M_2\end{gathered} M1​⎣⎢⎢⎢⎢⎡​ABCDE​⎦⎥⎥⎥⎥⎤​=M2​​ 可以得到 [ A B C D E ] = M 1 − 1 M 2 = M 1 \ M 2 \begin{bmatrix} A\\B\\C\\D\\E\end{bmatrix}=M_1^{-1}M_2=M_1\backslash M_2 ⎣⎢⎢⎢⎢⎡​ABCDE​⎦⎥⎥⎥⎥⎤​=M1−1​M2​=M1​\M2​

解得 A A A， B B B， C C C， D D D， E E E 后再根据前面的推导可以得出椭圆的相关参数。

MATLAB

输入散点坐标，并调用拟合函数