科学と家事とプログラミング (python を中心に)

python 温度計測 湿度計測 DS18B20 USB9097

トップ > プログラミング

楕円パラメータの相互変換 (幾何学的 ⇔ 代数的)

天文宇宙プログラミング覚書 プログラミング python

楕円の表現

この記事では、楕円を表す２種類の形式(方程式)と、 それらの間の相互変換について述べる。 ２種類の形式は、一般に 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

参考文献

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