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

python 温度計測 湿度計測 DS18B20 USB9097

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

楕円の表現

この記事では、楕円を表す2種類の形式(方程式)と、 それらの間の相互変換について述べる。 2種類の形式は、一般に 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} は、最初の2式を連立して求める。行列形式で表せば、

 \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}
}

備考

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

実装例

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. 金谷健一(他), "3次元コンピュータビジョン計算ハンドブック", 2016, 森北出版
  4. はてなブログで数式を書くurl3

図形を使った解法は、見ためも安心なのです(改)

コンピューターグラフィックスのアイデアが、ピタッとと当てはまる問題というのがあります。 教科書に出てきそうな Bang-Bang 制御の切り替え線の判定も、そんな問題の一つです。 図形の論理演算(and/or)を使って、制御の ON/OFF を判定する例を紹介します。

Bang-Bang 制御

Bang-Bang 制御は、位置  x と速度  \dot{x} に応じて、 一定量の制御を加えるかどうかを判定します。 いろんなパターンがありますが、たとえばこんな感じ(位相面)。

f:id:sken20k:20180525173650j:plain:w200

水色の領域では下向きに、緑色の領域では上向きに制御を加えます(速度を変化させる)。 折れ線になっているので、判定はちょっと面倒ですね。

よくある解法

よくある解法は、直線の交点を求めて x との大小関係で比較直線を選ぶというものでしょう。 そのままロジックを組むと、たくさんの交点を求めるのも面倒だし、 判定も if else がいくつも連なって、ややこしくなりそうです。

図形を使った解法

図形を使った解法はこうです。 比較的単純な領域を定義して、それらの論理和としてもとの領域を表現します。

f:id:sken20k:20180525173956j:plain:w400

制御を加える領域は、上の3つの領域の和になっています。 仕様(図計)を、プログラムに直せば、こんな感じになります。

def judge(x, y):
    y1, y2, y3, y4, y5 = f1(x), f2(x), f3(x),  f4(x), f5(x)
    if (y > y1 and y > y2) or (y>y3 and y>y4) or (y > y5):
        output = -1
   else:
        output = 0
    return output

仕様(図形を使った表現)と実装の対応が良いのも魅力です。 レビューをしっかりすれば、ほぼ間違えることはないのだけれど、 この手のロジックを検証する場合は、作った関数を使って、元の切り替え線(判定領域)を 作図してみるのが確実です(位相面を2次元でサンプリングする)。

おまけ

3次元CG のモデリング(画像生成)では、基本的な形状の組みあわせによって 複雑な形状を表現することがあります。 例えば、レンズの形は、球と球の and 領域として定義できます。 上の、Bang-Bang 制御の判定ロジックは、そこからアイデアをいただきました。 むかし CGや図形処理で遊んでいましたが、そういう経験が結構役立つ場面は多い気がします。 図形処理は、応用範囲が広いということですかね。

多変数最適化の Nelder-Mead は、なかなか使い勝手が良いのです

多変数最適化のアルゴリズムの一種である Nelder-Mead は、なかなか使い勝手の良いアルゴリズムですね。 汎用的に使える(対象が最小二乗に限定されない)ことに加えて、 微分情報を使わないところも魅力です。 評価関数の性格が多少悪くても、上手に最適をみつけてくれる、便利屋さんといったところでしょうか。

Nelder-Mead のお勧めバージョン

Nelder-Mead には、いろんな改良/拡張が存在しますが、 場合分けの見た目の安心感? から、以下の 5 種類の動作を含んだ設計がお勧めです。

  1. reflection
  2. expansion
  3. inside-contraction
  4. outside-contraction
  5. shrinking

細かく言うと、最悪, 2番悪, 最良 の 3点の情報を使って 7 通りに場合分けしたあと、 動作が同じになるところをまとめると、最終的に 5 ケースとなります。 例えば http://adl.stanford.edu/aa222/lecture_notes_files/chapter6_gradfree.pdf

アルゴリズムの妥当性と試験について

この手のアルゴリズムは、経験則の寄せ集め的なところがあって、 絶対的な正解という考え方には馴染まない気がします。 適用限界を明示するのが難しい。 あと設計やコーディングに間違いを含んでいても、 なんとなく動いてしまって、しかも、本来と異なる経路で最適値に収束したりして、 外部仕様を定めにくいですね。つまり、ブラックボックス的な試験設計は難しい。 開発手法(管理)の観点からは、興味深い対象かもしれません(扱いが難しいということ)。 というわけで、ホワイトボックス的な試験に重点を置くのが良さそうです。とにかく全てのパスを通す。 このとき、見やすい図を作って直感的に動作を確認するのは、 そう悪くない方法かと思ったりもします。例を下図に示します。

f:id:sken20k:20180513141435j:plain

適当な評価関数(2次元)を用意して、乱数で初期値を選ぶようにして、 見やすい図が得られるまで乱数のシードをいろいろ試しました。。 スライダーなどの GUI で初期値を選べるようにすると、おもしろいデモが 作れそうですね。

パラメータの一部を固定して評価したいとき

最適化の問題をあつかっていると、一部のパラメータを固定して推定したい。 という場面に出くわすことがあります。 パラメータの組み合わせ毎に最適化対象の関数を定義しはじめると、 評価関数だらけになってしまって、何をやってるのか分からなくなりがちです。 最適化のルーチンの I/F に手を加えて、要素毎に推定する/しないを指定するフラグを追加 しても良いのですが、 python のようなインタープリターを使うのであれば、 固定パラメータを補う関数を動的に生成するのも悪くないですね。 こうすれば、ベーシックな最適化ルーチンに手を加えずシンプルに保ったまま、 目的を達することができます。

落穂ひろい

Nelder-Mead 以外で、実用の見地から押さえておきたい最適化アルゴリズムは、 やっぱり Levenberg-Marquardt でしょう(LM法)。 適用できる問題は最小二乗形式に限定されますが、 変数の数が多めのデリケートな問題でも、良好な結果が得られることが多いですね。 推定誤差や correlation を評価がしやすいところもポイントが高いです。 注意が必要なのは、ダンピングファクターの初期値の調整くらいでしょうか。 実装にもよりますが、ダンピングファクターと収束判定の関係を理解していないと、 使いこなせない場面が、時々あります。 最適化ルーチンをブラックボックス的に利用するのは良いのですが、こういう 特性を理解した上でないと、使いこなせないこともある。 ブラックボックスの表面に、そういう細かな特性を認識できているかどうかで話が大分違ってきそうです。 個々のアルゴリズムや理論の前提条件や適用限界に注意を促したり、それらを身に付けるための方法論が 求められる。そう思う次第です。

剛体回転に関するオイラーの定理は、絵的に理解したいのです

「剛体回転におけるオイラーの定理」について、 天球図を使った直感的な説明を紹介します。 代数的な説明だけだと、なんとなく分かった気がしない...という人におすすめです。

オイラーの定理

オイラーの名の付く定理や公式は、いくつもあって紛らわしいのですが、 その中の一つに「剛体回転におけるオイラーの定理」というのがあります。

【定理】 「点のまわりの回転変位は, その点を通る1つの軸のまわりの回転によって達せられる」

3次元座標の座標回転の言葉で言い直すと、

「原点まわりの任意の座標回転は、ある方向ベクトルを回転の軸とする1回の回転で達成できる」

1回の回転で変換を達成する回転の軸を「オイラー軸」あるいは「オイラー回転軸」と呼び、 その時の回転の角度を「オイラー回転角」と呼んだりします。

オイラー角」と名前が似てますが別物です。「オイラー角」が座標軸回りの3回の回転の組み合わせであるのに対して、「オイラー回転角」は任意方向回りの1回の回転です。

幾何学的な説明

f:id:sken20k:20180427211221j:plain

上の図は、座標系 X-Y-Z から X'-Y'-Z' への回転を表しています。 ひとまず Z を除いて、X と Y に注目した上で、回転軸の成分が保存される (回転軸とXの成す角は、回転軸とX'の成す角と等しい)ことに注目すれば、 回転軸の方向は、 X, X' を底辺とする二等辺(球面)三角形の頂点に存在するはずです。つまり、 X を X' に移す回転軸は X と X' の垂直二等分線のどこかに存在する(球面なので大円のどこか)。 同様に、Y を Y' に移す回転軸は Y と Y' の垂直二等分線のどこかに存在する(球面なので大円のどこか)。 これらを同時に満たす点は2つの大円の交点として求まるので、これを P とする。 (逆方向の解も存在するが、どちらを選んでも議論の一般性は失われないので見やすい方を選ぶ)。

これで、オイラー回転軸の候補 P が手に入りましたが、 X', Y' を同時に満たす回転(角度)が存在するかどうか分かりません。 以下、それを確認します。

XをX'に移す回転角をφx、YをY'に移す回転角をφyとする。 オイラーの定理が成り立つためには、φx とφy が一致すればよい。

Pの成分が座標回転に対して不変なので、 球面三角 P-X-Y と P-X'-Y' は合同であることが分かる。 Pを頂点とする角度から重なり部分を除いた角度がφxとφyなので、 φxとφyは確かに等しい。

Z, Z' についても同様の議論が成り立つ。

これで、オイラーの定例が成り立つことが確認できました。

代数的な説明

先に幾何学的な説明をしましたが、 オイラーの定理は代数的に説明されることの方が多いですね。 座標回転を表す行列の固有値固有ベクトルを使った議論は、 あいまいさなく簡潔に証明できるところが利点と思います。 こちらは、どの教科書にも載っているので説明は省略します。

おまけ

オイラーの定理の絵的な理解は、実践的(応用範囲が広い)知見と思います。 図解した教科書が少ない(気がするので)ここに紹介しました。

補間と言えば akima-spline は知らなきゃもぐりだろう、と言ってみたいのです

スプライン補間でもっともメジャーなのは cubic-spline でしょうか。 ある種の理解しやすさや微分も一緒に求まるという使い良さもあるのですが、 データにギャップがある場合に、補間結果が振動するという やっかいな特性を持っています。 下の図の赤線は、サンプル点(黒点)を cubic-spline で補間した結果です。 ギャップの前後で補間値があばれているのが分かります。

同じ図に、akima-spline で補間した結果を青線で示します。 ギャップの前後も補間値があばれることなく、スムースな補間結果が得られことが分かります。

f:id:sken20k:20180420234011j:plain

参考までに、別の例も載せておきます。 f:id:sken20k:20180420235709j:plain

補間方法には、それぞれ利点/欠点はあるので、 akima-spline が万能というわけではありませんが、 区間ごとに推定した軌道や姿勢をつなぐ場面では、つなぎ目に微小なギャップが 現れることが多いので、そういうとき akima-spline のなめらかさ(オーバーシュートの無さ)は とてもうれしい特性なのでした。

akima は、秋間浩さんという日本人の名前に由来します。

高級な補間公式が、万能とは限らないのです。

補間公式には、ラグランジュ、エルミート、スプライン...と、各種ありますが、 私がよく使うのは、線形補間(直線補間)とAkima-Spline の 2 種類くらいです(JPL の惑星歴のチェビシェフは除く)。 線形補間(直線補間)とは、随分と原始的なものを使うのですね...と思う人がいるかもしれませんが、 高級な補間公式で安全性を保障するのは、結構むずかしかったりします。

高級な補間公式の弱点

高級な補間公式は、(使い方を間違わなければ) 少ないサンプル点で高い精度を実現できますが、 それらは以下のような欠点ももっています。

  • 単調性が保証されない
  • 高階差分を適用するとき桁落ちが発生しやすい

これらの欠点について注意が払われることは、少ないように思います。 数値計算の教科書も、精度の観点から議論が多くを占めて、 単調性が保証されない点に関しては、たとえ記述されていたとしても さらりと一言述べられるだけの事が多いようです。 記述量が少ないため、読み飛ばされることもしばしば。 実運用や解析ツールの設計においては、とても大切な特性ですが、 痛い目に会わないと、なかなか身に付かない知見の一つかと思い、 ここに紹介します。

単調性の保証が大事な時

例えば、軌道情報が時系列データ(点列)として与えられたときに、 距離 R がある閾値を下回るタイミングを求めたいとする。 このとき、単調性が保証されない補間公式を使うと、 偽のイベントがジッターのように出現する可能性がある。 次の図では、本当のタイミングは time が 1.1 を超えあたりにありそうだけど、 高次の補間を使って求めようとしたために 0.9 ~ 1.0 の間に偽の事象が出現しています。

f:id:sken20k:20180417185817j:plain

補間の時に最後の Time = 1.2 の点を使わないのが原因ではないか? と思う人がいるかもしれませんが、Time = 1.0 までの点が揃ったところで、 直近のN 点から事象を判定する場面は、そう無理な仮定ではなかろうと思います。 高次の補間公式で単調性を保証するのは、いろいろとむずかしいので、 とにかく安全を保証したい場面でであれば、単調性が保証できる線形補間を使うことを お勧めする次第です。

上の図は、単調性が保たれない箇所を拡大しています。元のスケールで補間の様子全体を 表示すると以下のようになります(6点を使ったラグランジュ補間)。 f:id:sken20k:20180420183939j:plain

高階差分の桁落ち

高次の補間は、単調性の話以外にも精度の面で問題が生じる場合があります。 具体的には、補間公式には高階の差分演算が現れるため、似たような数値の減算の繰り返し が多くなり、次数が高くなるほど、つまり階数が深くなるほど、桁落ちが生じやすくなるからです。 これに対抗するには、最初のデータの有効桁を増やすしかありません。 時刻を例にとれば、準ユリウス日ユリウス日を時刻に使っていたりすると、 高階差分のために思わぬ誤差が生じることがあるので注意が必要です。 最近では、2000年からの時刻が使われたりすることもありますが、高い精 度が求められる場合には、解析期間の開始時刻を時刻原点に採用するのが無難です。 普段、あまり次数の高い補間公式を使わないようにしていても、DOP853 という 高精度の数値積分のルーチンに組み込まれた densout (稠密出力)などを利用する 場面で、あまり意識することなく、高階の差分を使ってしまうこともあります。

落穂ひろい

高次の補間公式が振動する話は有名ですが、単調性に焦点を当てた議論というのは 少ないように思います。教科書を書いている数値計算の専門家にとっては あたりまえのことであるので、ページ数が限られている教科書においては、 話が精度に偏ってしまうのも仕方ないことかもしれません。あるいは、 数値計算の専門家(研究者)には、実運用のシステムを設計/運用する経験はほとんどないと 思われるので、適用場面で重用になる特性に興味が薄いことも一因かと疑ったりしています。 というわけで、たとえあたりまえの話であっても、ある程度のボリュームを割いて、議論する意味 というのはあるのかなと思う次第です。もちろん、知識の効率的な吸収という観点とは相反するわけだけれども、 ある種の無駄は必要だろうという立場をとります。

少し一般的な話をすると、およそ公式や手法には、それぞれ利点だけでなく欠点もあるものです。 使える場面や条件を明晰に意識できるかどうかが、応用が利くかどうかの分かれ目ですね。 利点だけでなく欠点や制約を2つ3つ列挙できるかどうかが、その手法の理解度を測るバロメータ だと言う人もいます。

Akima Spline については、別記事にて...

緯度経度で分割された星カタログの選択は、なかなかやっかいなのです

カメラ視野を模擬したい

カメラの視野を模擬するとき、 手元の星カタログを全部スキャンすると星の数が多すぎて性能が出ない時があります。 こういう時は、星カタログをいくつかの部分に分けて対応します。 分けたものをサブカタログと呼んだりします。

経度緯度で分割したサブカタログ

サブカタログの分割方法で一番簡単なのは、緯度と経度で分ける方法でしょうか。 さて、分割したからには必要なサブカタログがどれかを調べなければいけません。 簡単なのは、各サブカタログの4隅(極だと3隅)と中心方向が、カメラの視野に入るかどうかを判定することでしょうか。 こういう方法だと、辺の一部に視野が重なる場合などはうまく判定できません。 特に困るのが極付近の判定です。その難しさは、下の図を見てもらえば分かってもらえるでしょうか。

力技による判定は、細かく面でサンプルする方法ですが、判定に時間がかかって残念な感じになってしまいます。 なんとかスマートに判定したいものです。

ビットマップを使ったサブカタログの選択

ここで、ビットマップを使った視野判定の方式を紹介します。

f:id:sken20k:20180330172056j:plain

手順はこうです

  1. 視野より大きめの描画エリアを用意して、カタログの境界線を描画する
  2. 背景色で視野を抜く
  3. 中心付近を視点として、塗りつぶす。
  4. 各サブカタログの中心(代表点)の色を調べる

細かな説明は不要と思いますが、これで視野の端がサブカタログにちょっとだけかかるときも、 ビットマップの分解能で正しく判定できます。実用上は、視野を少し大きく採ればよいでしょう。 けっきょく、面でサンプルしているといえなくもないのですが、サブカタログの境界線と境界線の間は、 ビットマップの塗りつぶし処理だけで代用できているので、ちょっとお得かなと思います。

他の方式

球座標のような特異点を避けるという観点からは、N 面体で分割するのも一案です。 用途によっては、排他的に分割するのではなく、重複を許すと便利そうです。 他にも分割の工夫はありそうですが、最近は使えるメモリーも桁違いに増えたので、 一括で扱った方が良い場面も多いかもしれません。 ビットマップを使った解法というのは、他にも応用が利きそうに思うので、一応紹介しておきます。