楕円パラメータの相互変換 (幾何学的 ⇔ 代数的) - 科学と家事とプログラミング (python を中心に)

楕円の表現
- 楕円の幾何学的な表示形式
- 楕円の代数的な表示形式
幾何学的な表示(x0,y0,θ,a,b) ⇒ 代数的な表示(A,B,C,D,E)
代数的な表示(A,B,C,D,E) ⇒ 幾何学的な表示(x0,y0,θ,a,b)
備考
実装例
参考文献

楕円の表現

この記事では、楕円を表す２種類の形式(方程式)と、それらの間の相互変換について述べる。２種類の形式は、一般に canonical form, general form と呼ばれることが多いが、ここでは便宜上、「幾何学的な表示形式」、「代数的な表示形式」と呼ぶことにする。 wikipedia(英語版) の ellipse の項をみると相互変換に関する解析解が掲載されている。それを使っても良いが、式の形や導出がけっこうややこしいので、別解として、比較的簡単に数値的に求める方法を紹介する。

f:id:sken20k:20181105181127p:plain

楕円の幾何学的な表示形式

幾何学的に楕円を表現する時は、パラメータとして $x_0, y_0, \theta, a, b$ を選ぶことが多い。ここで、 $x_0, y_0$ は楕円の中心、 $\theta$ は楕円の傾き、 $a, b$ は楕円の半径(横方向と縦方向)である。 $\cos \theta$ , $\sin \theta$ を $c, s$ で表すと、楕円の方程式は、以下の通り。

$\displaystyle { \frac{ ( c(x - x _ 0) + s(y - y _ 0) )^ 2 }{a^ 2} + \frac{ ( -s(x - x _ 0) + c(y - y _ 0) )^ 2 }{b^ 2} = 1 }$

楕円の代数的な表示形式

代数的に楕円を表現する時は、係数を x, y の冪乗ごとにまとめる。 $x^{2}$ の係数が 1 になるように規格化し、残りの項の係数を $A,B,C,D,E$ とすれば、楕円の方程式(一般的な二次曲線)は、以下の通り。

$\displaystyle { x^{2} + A x y + B y^{2} + C x + D y + E = 0 }$

幾何学的な表示(x0,y0,θ,a,b) ⇒ 代数的な表示(A,B,C,D,E)

幾何学的な表現を展開したあと、 $x, y$ の冪でまとめなおして係数を比較する。

$\begin{aligned} \alpha A &= \gamma \\ \alpha B &= \beta \\ \alpha C &= -2\alpha x_{0} - \gamma y_{0} \\ \alpha D &= -2\beta y_{0} - \gamma x_{0} \\ \alpha E &= \alpha x_{0}^{2} + \beta y_{0}^{2} + \gamma x_{0} y_{0} - 1 \\ \end{aligned}$

ここで、

$\displaystyle{ \alpha = \frac{c^{2}}{a^{2}} + \frac{s^{2}}{b^{2}}, \qquad \beta = \frac{s^{2}}{a^{2}} + \frac{c^{2}}{b^{2}}, \qquad \gamma = 2 c s \left( \frac{1}{a^{2}} - \frac{1}{b^{2}} \right) }$

である。代数的なパラメータ $A,B,C,D,E$ は、上の関係から直ちに求まる。

代数的な表示(A,B,C,D,E) ⇒ 幾何学的な表示(x0,y0,θ,a,b)

$\gamma, \beta$ を消去して、全体を $\alpha$ で割ると以下を得る。

$\begin{aligned} C &= -2 x_{0} - A y_{0} \\ D &= -2 B y_{0} - A x_{0} \\ E &= x_{0}^{2} + B y_{0}^{2} + A x_{0} y_{0} - 1/\alpha \end{aligned}$

$x_{0}, y_{0}$ は、最初の２式を連立して求める。行列形式で表せば、

$\displaystyle{ \begin{pmatrix} C \\ D \end{pmatrix} = \begin{pmatrix} -2 & -A \\ -A & -2B \end{pmatrix} \begin{pmatrix} x_{0} \\ y_{0} \end{pmatrix} }$

$\theta$ は、以下の関係を使って、A,B から直接求める。

$\displaystyle{ \tan(2\theta)=\frac{2cs}{c^{2}-s^{2}}=\frac{\gamma}{\alpha - \beta}=\frac{A}{1 - B} }$

$\alpha$ は、 $x_{0}, y_{0}$ を最後の式に代入して求める。

$\displaystyle{ 1/\alpha = x_{0}^{2} + B y_{0}^{2} + A x_{0} y_{0} - E }$

$a, b$ は、 $\theta$ を $\alpha, \beta$ の定義式に代入して、それらを連立して求める。

$\displaystyle{ \begin{pmatrix} \alpha \\ \beta \end{pmatrix} = \begin{pmatrix} \alpha \\ \alpha B \end{pmatrix} = \begin{pmatrix} c^{2} & s^{2} \\ s^{2} & c^{2} \end{pmatrix} \begin{pmatrix} 1/a^{2} \\ 1/b^{2} \end{pmatrix} }$

備考

「超精度くりこみ法」というのを知ったのをきっかけに、楕円のパラメーター推定をあれこれ試しています。その過程で、楕円パラメータの相互変換が必要になったので、忘備録として記事を書きました。部分的なアークから楕円を推定する方法は、天文の分野でも応用が広そうですね。参考文献に挙げた「３次元～」は、各種アルゴリズムのソースコードが出版社のページで公開されていてます。コードを眺めるだけで動かしていませんが、こういうサービスはうれしいですね。

実装例

def geo_fm_ana5(A, B, C, D, E):
    """
    Parameters
    ----------
    A～E : 楕円方程式の係数 
                 (x^2 + A*x*y + B*y^2 + C*x + D*y + E = 0)
    Returns
    -------
    xc, yc : [-]   中心座標
    th     : [rad] 回転角
    ra, rb : [-]   横半径, 縦半径
    """
    xc, yc = np.linalg.solve([[-2, -A], [-A, -2*B]], [C, D])
    th = np.arctan2(A, 1 - B)/2
    c = np.cos(th)
    s = np.sin(th)
    k = 1/(xc*xc + B*yc*yc + A*xc*yc - E)
    a_inv2, b_inv2 = np.linalg.solve([[c*c, s*s], [s*s, c*c]], [k, k*B])
    ra = np.sqrt(1/a_inv2)
    rb = np.sqrt(1/b_inv2)
    return xc, yc, th, ra, rb

参考文献

最小二乗法による楕円近似 / url1
楕円当てはめの精度比較：最小二乗法から超精度くりこみ法まで url2
金谷健一(他), "３次元コンピュータビジョン計算ハンドブック", 2016, 森北出版
はてなブログで数式を書くurl3