2次曲線の式から楕円パラメータを導出する



意外とネット上で探しても出てこなかったので、頑張って導出してみようと思います。

2次曲線の式っていうのは$ax^2+bxy+cy^2+dx+ey+f=0$みたいな2次関数のことです。
楕円パラメータっていうのは楕円中心$(x_0,y_0)$、X軸方向半径$a$、Y軸方向半径$b$、回転$\theta$のことのつもりです。
媒介変数とかのことではありません。
![楕円.jpg](https://qiita-image-store.s3.amazonaws.com/0/175507/daa7823a-4641-06da-36e8-3b14443770fe.jpeg)

#2次曲線の式
2次曲線の式は以下の数式で表現することにします。

```math
\left[\begin{array}{c} x & y & f_0 \end{array}\right]
\left[\begin{array}{c} A&B&D\\B&C&E\\D&E&F \end{array}\right]
\left[\begin{array}{c}x\\y\\f_0\end{array}\right]
= Ax^2+2Bxy+Cy^2+2Dxf_0+2Eyf_0+Ff_0^2=0
```

この$f_0$ってのは見慣れないかと思いますが、各項の次元をそろえるための定数です。数値計算的には適当な値に設定することで計算の精度を保つことができます（$F$だけめっちゃ大きくなるみたいなのを防げる）。
基本的には$f_0=1$と考えちゃっていいです。

参考- http://www.iim.cs.tut.ac.jp/~kanatani/papers/newlsellipse.pdf

#結論
先に結論から述べると、

```math
\begin{eqnarray}
x_0 &=& \frac{BE-CD}{AC-B^2}f0 \\
y_0 &=& \frac{BD-AE}{AC-B^2}f0 \\
a^2 &=&\frac{2\left(\left(CD-BE\right)D + \left(AE-BD\right)E - (AC-B^2)F\right)}{(AC-B^2)\left(A+C+\sqrt{(A-C)^2+4B^2}\right)}f_0^2 &=& -2f_0\frac{Dx_0+Ey_0+Ff_0}{A+C+\sqrt{(A-C)^2+4B^2}} \\
b^2 &=&\frac{2\left(\left(CD-BE\right)D + \left(AE-BD\right)E - (AC-B^2)F\right)}{(AC-B^2)\left(A+C-\sqrt{(A-C)^2+4B^2}\right)}f_0^2 &=& -2f_0\frac{Dx_0+Ey_0+Ff_0}{A+C-\sqrt{(A-C)^2+4B^2}} \\
\theta &=& \frac{1}{2}\tan^{-1}\frac{2B}{A-C}
\end{eqnarray}
```

となります。
以下頑張って導出していきます。

#楕円の方程式
楕円の方程式ですが、基本的な楕円の方程式$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$を回転&平行移動させて

```math
\left(\frac{\left(x-x_0\right)\cos\theta + \left(y-y_0\right)\sin\theta}{a}\right)^2 + \left(\frac{-\left(x-x_0\right)\sin\theta + \left(y-y_0\right)\cos\theta}{b}\right)^2 = 1
```
となり、これを$x,y$に関して整理して

```math
\begin{eqnarray}
&&\left(\frac{1}{a^2}\cos^2\theta+\frac{1}{b^2}\sin^2\theta\right)x^2 \\
&+&2\left(\frac{1}{a^2}-\frac{1}{b^2}\right)\sin\theta \cos\theta xy \\
&+&\left(\frac{1}{a^2}\sin^2\theta+\frac{1}{b^2}\cos^2\theta\right)y^2 \\
&-&2\left(\frac{1}{a^2}\cos\theta\left(x_0\cos\theta + y_0\sin\theta\right) - \frac{1}{b^2}\sin\theta\left(-x_0\sin\theta + y_0\cos\theta\right)\right)x \\
&-&2\left(\frac{1}{a^2}\sin\theta\left(x_0\cos\theta + y_0\sin\theta\right) + \frac{1}{b^2}\cos\theta\left(-x_0\sin\theta + y_0\cos\theta\right)\right)y \\
&+&\frac{1}{a^2}\left(x_0\cos\theta + y_0\sin\theta\right)^2 + \frac{1}{b^2}\left(-x_0\sin\theta + y_0\cos\theta\right)^2 - 1 \\
&=&0
\end{eqnarray}
```

となります。

#両式の比較
2次曲線の式と楕円の方程式を比較して、2次曲線の式の各パラメータにはスケールの自由度があることを考えると

```math
\begin{eqnarray}
At &=& \frac{1}{a^2}\cos^2\theta+\frac{1}{b^2}\sin^2\theta \\
Bt &=& \left(\frac{1}{a^2}-\frac{1}{b^2}\right)\sin\theta \cos\theta \\
Ct &=& \frac{1}{a^2}\sin^2\theta+\frac{1}{b^2}\cos^2\theta \\
Dt &=& -\frac{1}{f_0}\left(\frac{1}{a^2}\cos\theta\left(x_0\cos\theta + y_0\sin\theta\right) - \frac{1}{b^2}\sin\theta\left(-x_0\sin\theta + y_0\cos\theta\right)\right) \\
Et &=& -\frac{1}{f_0}\left(\frac{1}{a^2}\sin\theta\left(x_0\cos\theta + y_0\sin\theta\right) + \frac{1}{b^2}\cos\theta\left(-x_0\sin\theta + y_0\cos\theta\right)\right) \\
Ft &=& \frac{1}{f_0^2}\left(\frac{1}{a^2}\left(x_0\cos\theta + y_0\sin\theta\right)^2 + \frac{1}{b^2}\left(-x_0\sin\theta + y_0\cos\theta\right)^2 - 1\right)
\end{eqnarray}
```

ここで、$t$はスケール方向の自由度です。


#楕円パラメータの導出
まず$A$と$C$の式の差をとって

```math
\begin{eqnarray}
At-Ct &=& \frac{1}{a^2}\left(\cos^2\theta-\sin^2\theta\right) + \frac{1}{b^2}\left(\sin^2\theta-\cos^2\theta\right) \\
(A-C)t &=& \left(\frac{1}{a^2}-\frac{1}{b^2}\right)\left(\cos^2\theta-\sin^2\theta\right) \\
\frac{1}{a^2}-\frac{1}{b^2} &=& \frac{A-C}{\cos^2\theta-\sin^2\theta}t
\end{eqnarray}
```

$B$の式に代入して

```math
\begin{eqnarray}
Bt &=& (A-C)\frac{\sin\theta \cos\theta}{\cos^2\theta - \sin^2\theta}t \\
B &=& (A-C)\frac{1}{2}\frac{\sin2\theta}{\cos2\theta} \because 半角公式 \\
\tan2\theta &=& \frac{2B}{A-C} \\
\theta &=& \frac{1}{2}\tan^{-1}\frac{2B}{A-C}
\end{eqnarray}
```

となり$\theta$が求められました。
次に$D$と$E$の式を使って$\frac{1}{b^2}$を消去します。

```math
\begin{eqnarray}
Dt\cos\theta + Et\sin\theta = -\frac{1}{f_0a^2}\left(x_0\cos\theta + y_0\sin\theta\right)
\end{eqnarray}
```
同様に$\frac{1}{a^2}$を消去します。

```math
\begin{eqnarray}
-Dt\sin\theta + Et\cos\theta = -\frac{1}{f_0b^2}\left(-x_0\sin\theta + y_0\cos\theta\right)
\end{eqnarray}
```

この2式から$E$の項を消去して

```math
\begin{eqnarray}
\left(Dt\cos\theta + Et\sin\theta\right)\cos\theta - \left(-Dt\sin\theta + Et\cos\theta\right)\sin\theta &=&
 -\frac{1}{f_0}\left(\frac{1}{a^2}\left(x_0\cos\theta + y_0\sin\theta\right)\cos\theta
 - \frac{1}{b^2}\left(-x_0\sin\theta + y_0\cos\theta\right)\sin\theta\right)\\

Dt\left(\cos^2\theta + \sin^2\theta\right) &=&
 -\frac{1}{f_0}\left(\left(\frac{1}{a^2}\cos^2\theta+\frac{1}{b^2}\sin^2\theta\right)x_0
 + \left(\frac{1}{a^2}-\frac{1}{b^2}\right)\sin\theta \cos\theta y_0\right) \\

Dt &=& -\frac{1}{f_0}\left(Atx_0 + Bty_0\right) \\

Ax_0 + By_0 &=& -Df_0
\end{eqnarray}
```

同様に

```math
\begin{eqnarray}
\left(Dt\cos\theta + Et\sin\theta\right)\sin\theta + \left(-Dt\sin\theta + Et\cos\theta\right)\cos\theta &=&
 -\frac{1}{f_0}\left(\frac{1}{a^2}\left(x_0\cos\theta + y_0\sin\theta\right)\sin\theta
 + \frac{1}{b^2}\left(-x_0\sin\theta + y_0\cos\theta\right)\cos\theta\right)\\

Et\left(\sin^2\theta + \cos^2\theta\right) &=&
 -\frac{1}{f_0}\left(\left(\frac{1}{a^2}-\frac{1}{b^2}\right)\sin\theta \cos\theta x_0
 + \left(\frac{1}{a^2}\sin^2\theta+\frac{1}{b^2}\cos^2\theta\right)y_0\right) \\

Et&=& -\frac{1}{f_0}\left(Btx_0 + Cty_0\right) \\

Bx_0 + Cy_0 &=& -Ef_0
\end{eqnarray}
```

連立方程式

```math
\begin{eqnarray}
Ax_0 + By_0 &=& -Df_0 \\
Bx_0 + Cy_0 &=& -Ef_0
\end{eqnarray}
```

を解いて

```math
\begin{eqnarray}
x_0 &=& \frac{BE-CD}{AC-B^2}f0 \\
y_0 &=& \frac{BD-AE}{AC-B^2}f0
\end{eqnarray}
```

となり、楕円中心$x_0,y_0$を求められました。

$A$と$C$の式を変形して$\frac{1}{b^2}$の項を消去すると

```math
\begin{eqnarray}
At\cos^2\theta - Ct\sin^2\theta &=& \frac{1}{a^2}\left(\cos^4\theta - \sin^4\theta\right) \\

\left(A\cos^2\theta - C\sin^2\theta\right)t &=&
 \frac{1}{a^2}\left(\cos^2\theta - \sin^2\theta\right)\left(\cos^2\theta + \sin^2\theta\right) \\

\left(A\cos^2\theta - C\sin^2\theta\right)t &=& \frac{1}{a^2}\left(\cos^2\theta - \sin^2\theta\right) \\

\frac{1}{a^2} &=& \frac{A\cos^2\theta - C\sin^2\theta}{\cos^2\theta - \sin^2\theta}t
\end{eqnarray}
```
となり、同様に

```math
\begin{eqnarray}
\frac{1}{b^2} &=& \frac{A\sin^2\theta - C\cos^2\theta}{\sin^2\theta - \cos^2\theta}t
\end{eqnarray}
```

また、$F$の式から$t$を求めると

```math
\begin{eqnarray}
t &=& \frac{1}{Ff_0^2}\left(\frac{1}{a^2}\left(x_0\cos\theta + y_0\sin\theta\right)^2 + \frac{1}{b^2}\left(-x_0\sin\theta + y_0\cos\theta\right)^2 - 1 \right) \\
&=& \frac{1}{Ff_0^2}\left(\frac{X^2}{a^2} + \frac{Y^2}{b^2} - 1\right)
\end{eqnarray}
```

となります。ただし$X=x_0\cos\theta + y_0\sin\theta, Y=-x_0\sin\theta + y_0\cos\theta$としました。これを上の$\frac{1}{a^2},\frac{1}{b^2}$の式に代入して

```math
\begin{eqnarray}
\frac{1}{a^2} &=& \frac{A\cos^2\theta - C\sin^2\theta}{\left(\cos^2\theta - \sin^2\theta\right)Ff_0^2}\left(\frac{X^2}{a^2} + \frac{Y^2}{b^2} - 1\right)
 = k_1\left(\frac{X^2}{a^2} + \frac{Y^2}{b^2} - 1\right) \\
\frac{1}{a^2} &=& \frac{A\sin^2\theta - C\cos^2\theta}{\left(\sin^2\theta - \cos^2\theta\right)Ff_0^2}\left(\frac{X^2}{a^2} + \frac{Y^2}{b^2} - 1\right)
 = k_2\left(\frac{X^2}{a^2} + \frac{Y^2}{b^2} - 1\right)
\end{eqnarray}
```

となり、この連立方程式を$a^2,b^2$について解くと

```math
\begin{eqnarray}
a^2 &=& \frac{k_1X^2 + k_2Y^2 - 1}{k_1} \\
b^2 &=& \frac{k_1X^2 + k_2Y^2 - 1}{k_2}
\end{eqnarray}
```

ただしここで、

```math
\begin{eqnarray}
k_1 &=& \frac{A\cos^2\theta - C\sin^2\theta}{\left(\cos^2\theta - \sin^2\theta\right)Ff_0^2} \\
&=& \frac{\frac{1+\cos2\theta}{2}A - \frac{1-\cos2\theta}{2}C}{\frac{1+\cos2\theta}{2} - \frac{1-\cos2\theta}{2}}\frac{1}{Ff_0^2} \because 半角公式\\
&=& \frac{A-C + (A+C)\cos2\theta}{2Ff_0^2\cos2\theta}
\end{eqnarray}
```

同様に

```math
\begin{eqnarray}
k_2 = \frac{C-A + (A+C)\cos2\theta}{2Ff_0^2\cos2\theta}
\end{eqnarray}
```

です。
ここで、

```math
\begin{eqnarray}
X &=& \left(x_0\cos\theta + y_0\sin\theta\right)^2 \\
&=& x_0^2\cos^2\theta + 2x_0y_0\sin\theta\cos\theta + y_0^2\sin^2\theta \\
&=& \frac{1}{2}\left(x_0^2\left(1+\cos2\theta\right) + 2x_0y_0\sin2\theta + y_0^2\left(1-\cos2\theta\right)\right) \\
Y &=& \left(-x_0\sin\theta + y_0\cos\theta\right)^2 \\
&=& x_0^2\sin^2\theta - 2x_0y_0\sin\theta\cos\theta + y_0^2\cos^2\theta \\
&=& \frac{1}{2}\left(x_0^2\left(1-\cos2\theta\right) - 2x_0y_0\sin2\theta + y_0^2\left(1+\cos2\theta\right)\right)
\end{eqnarray}
```
と$k_1,k_2$を代入してまとめると

```math
\begin{eqnarray}
a^2 &=& \frac{2\left(Ax_0^2\cos2\theta + (A-C)x_0y_0\sin2\theta + Cy_0^2\cos2\theta -Ff_0^2\cos2\theta\right)}{A-C+(A+C)\cos2\theta} \\
&=& \frac
{2\left(A(A-C)x_0^2 + 2B(A-C)x_0y_0 + C(A-C)y_0^2 - F(A-C)f_0^2\right)}
{(A-C)\sqrt{(A-C)^2+4B^2} + (A+C)(A-C)} \\
&& \because \tan2\theta = \frac{2B}{A-C} \Rightarrow \sin2\theta = \frac{2B}{\sqrt{(A-C)^2+4B^2}}, \cos2\theta = \frac{A-C}{\sqrt{(A-C)^2+4B^2}} \\
&=&\frac{2\left(Ax_0^2 + 2Bx_0y_0 + Cy_0^2 - Ff_0^2\right)}{A+C+\sqrt{(A-C)^2+4B^2}}
\end{eqnarray}
```

同様に

```math
\begin{eqnarray}
b^2 =\frac{2\left(Ax_0^2 + 2Bx_0y_0 + Cy_0^2 - Ff_0^2\right)}{A+C-\sqrt{(A-C)^2+4B^2}}
\end{eqnarray}
```
となります。
$x_0 = \frac{BE-CD}{AC-B^2}f0, y_0 = \frac{BD-AE}{AC-B^2}f0$だったので、これを代入してまとめると

```math
\begin{eqnarray}
a^2 &=&\frac
{2\left(\left(CD-BE\right)D + \left(AE-BD\right)E - (AC-B^2)F\right)}
{(AC-B^2)\left(A+C+\sqrt{(A-C)^2+4B^2}\right)}f_0^2 \\
b^2 &=&\frac
{2\left(\left(CD-BE\right)D + \left(AE-BD\right)E - (AC-B^2)F\right)}
{(AC-B^2)\left(A+C-\sqrt{(A-C)^2+4B^2}\right)}f_0^2
\end{eqnarray}
```

となり、楕円の径$a,b$を求められます。


#まとめ
最後にもう一度まとめると、

```math
\begin{eqnarray}
x_0 &=& \frac{BE-CD}{AC-B^2}f0 \\
y_0 &=& \frac{BD-AE}{AC-B^2}f0 \\
a^2 &=&\frac{2\left(\left(CD-BE\right)D + \left(AE-BD\right)E - (AC-B^2)F\right)}{(AC-B^2)\left(A+C+\sqrt{(A-C)^2+4B^2}\right)}f_0^2 &=& -2f_0\frac{Dx_0+Ey_0+Ff_0}{A+C+\sqrt{(A-C)^2+4B^2}} \\
b^2 &=&\frac{2\left(\left(CD-BE\right)D + \left(AE-BD\right)E - (AC-B^2)F\right)}{(AC-B^2)\left(A+C-\sqrt{(A-C)^2+4B^2}\right)}f_0^2 &=& -2f_0\frac{Dx_0+Ey_0+Ff_0}{A+C-\sqrt{(A-C)^2+4B^2}} \\
\theta &=& \frac{1}{2}\tan^{-1}\frac{2B}{A-C}
\end{eqnarray}
```

となります。


#感想

とけてよかったとおもいました。
2次方程式が楕円でないときは$a^2$や$b^2$がマイナスになって$a,b$を実数で求められなくなったり、$AC-B^2$が$0$になって$x_0,y_0$が求められなくなったりするみたいです。$a,b$の式はなんかもうちょっと簡単になりそうですが私にはこれが限界でした。楕円の判別は$AC-B^2>0$なので、$AC-B^2>0$で2次曲線が図形になるなら$a^2,b^2>0$になると思うんですが、そのへんの証明もちょっと出来ませんでした。でもまぁ一応は楕円パラメータを求めることが出来たので良かったです。

それではさようなら。