20170214

matplotlibのaxes3Dで、aspectの設定がうまくできなかったので、行列とかをいじった。

まだ、フーリエ変換関係の動画を作ろうとしている。

  • DFT計算詳細
  • iDFT計算詳細
  • FFT
  • 2DDFT
  • i2DDFT

については、詳細動画作りたい。作ってて楽しいし。

その過程で、フーリエ変換の結果を、3Dで図示する必要があった。(Re, Im, kの3軸)

matplotlibのaxes3dを使った。

from mpl_toolkits.mplot3d.axes3d import Axes3D
from mpl_toolkits.mplot3d import proj3d

で、3Dプロットをしている。

この機能は、重なった部分を描画すると、最後に書いた形状が、常に手前に表示されるという未解決のバグがあるが、複雑な形状じゃなければさほど問題にならない。(もし問題があれば、mayaviを使う

軸のaspectがうまくいかなかった。

何も指定しないと、自動的に、余ったスペースにフィットする感じにプロットが拡大される。
しかし、縦横比が気に入らない。
Re・Im比が気に入らない.png

アスペクト比を設定する、


ax.set_aspect('equal')

は、x:y:z全てを、1:1:1にする。Nが大きい場合、64:1:1とかになるから、問題がある。
だから、

  • 複素平面を±1〜2の範囲で描く。
  • k軸を、その1.5〜2倍の大きさで描く。

事を目指した。
2Dプロットの場合は、

ax.set_aspect()

の引数に、数値を入れれば指定できた。数値は一個しか取れない。x:yの比なので、一個でも表現できる。
3Dの場合にも、同じメソッドは存在する。なぜかこちらも数値は一個しか取れない。x:y:zの比は、一個では表現でき…る?無理だと思うんだけど…

  • aspectメソッド自体はあっても、2Dプロット用のためか、x,yの二つしか設定できなかった。
  • x,y,zのバランスが整えられない。

という結論で、アスペクトをいじるための、他の。方法がほしくなった。

aspectの解決策

ここを参考にして、matplotlibのaxes3dのaspectを調節した

ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([1.2, 0.5, 1.5, 1]))

このコマンドだけで、下の図のように、Re・Im(縦・横)比は、なんとなく解決できた。

上記方法だと、プロット全体が、画面の中心からずれる。

しかし、余白が目立つ。どうにもシフトしてしまうようだ。

get_projで生まれる余白.png
こんな感じになる。
いや、もともと、この位の隙間が最初からあったんだっけ?
とにかく、余白が気に入らない。

あ、コレ、アフィン行列だ。

アスペクトを直すときに使ったコマンドをもう一度載せる。

ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([1.2, 0.5, 1.5, 1]))

get_projには、Axes3D.get_proj(ax)と、アフィン行列の、行列の積、を計算する、lambda式を入れる。
ちなみに、上記の、

np.diag([1.2, 0.5, 1.5, 1])

の結果は、

>>> np.diag([1.2, 0.5, 1.5, 1])
array([[ 1.2, 0. , 0. , 0. ],
[ 0. , 0.5, 0. , 0. ],
[ 0. , 0. , 1.5, 0. ],
[ 0. , 0. , 0. , 1. ]])

$$\rm aff= \left[ \begin{array}{ccc}1.2&0&0&0 \\\\ 0&0.5&0&0 \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]$$

この行列は、アフィン変換の拡大縮小行列ですな。この行列と、図中の座標の、行列の積(np.dotsで計算している部分。行列の積)を計算する関数(lambda式)が、get_projに入るという流れ。

アフィン変換のwikipediaは、難しくて読めないけど、「回転、拡大・縮小、平行移動、ミラー、せん断行列の事が、アフィン変換」というのが、自分の認識。多分、それ以外にも使うことがあるから、wikiでは、ああなるんだと思う。フーリエ変換も、そのせいで意味不明だし。まぁ、wikipediaは、あくまでも百科事典で、教科書じゃないから。

拡大縮小の行列は、3Dなら、

$$\left[ \begin{array}{ccc} x'\\\\ y' \\\\ z' \\\\ 1 \end{array}\right] = \left[ \begin{array}{ccc}a&0&0&0 \\\\ 0&b&0&0 \\\\ 0&0&c&0 \\\\ 0&0&0&1 \end{array}\right]\left[ \begin{array}{ccc}x\\\\ y \\\\ z \\\\ 1 \end{array}\right]$$

なはず。
3Dなのに、4D用意しているのは、平行移動行列の時とか、

$$\left[ \begin{array}{ccc} x'\\\\ y' \\\\ z' \\\\ 1\end{array}\right] = \left[ \begin{array}{ccc}1&0&0&d \\\\ 0&1&0&h \\\\ 0&0&1&l \\\\ 0&0&0&1 \end{array}\right]\left[ \begin{array}{ccc}x\\\\ y \\\\ z \\\\ 1 \end{array}\right]$$

だから、z軸動かしたい時は、3行目の4列目がないと困る。などがある。拡大・縮小や、回転などは、3*3で足りる変換もあるが、4*4が必要な他の変換等を、複数組み合わせて使うときに、同じサイズの行列だと、あらかじめかけ算して、ひとまとめにしてから計算することが可能だから、4*4にする。

  • そこは、マークしてるんだおれは。その内やろうとは、思っているんだ。
  • 円周上の点の座標を、回転行列に入れて、計算とかしたんだ。
  • その時、アフィン行列という分野を知ったんだ。
  • だから、アフィン行列のLaTeXは、2D 用から3D用まで、全部用意しているんだ。
  • あとは、フーリエへの動画を作り終わるだけなんだ。
  • そしたら、アフィン行列のアニメを作るんだ。
  • その後、レンダリング方程式のアニメも作るんだ。

んだんだ。

行列のかけ算、アフィン行列。この辺は、フーリエ変換とも直接関係する部分なので、存じておりました。
フーリエ変換は、フーリエ変換行列と、データの、行列のかけ算の計算だから。

np.diagは、対角項を指定した、配列(行列)を作る。

np.diagは、対角項というか、対角行列。
斜めのやつですな。
英語なら、diagonal matrixとか、この場合は、diagonal arrayとか。

当然、正方行列。
wikipediaに、そう載ってた

理由はともあれ、解決策

縦横比の指定がアフィンの拡大行列を使って行われていたので、
平行移動も、アフィンで行うことにする。

$$\left[ \begin{array}{ccc} x'\\\\ y' \\\\ z' \\\\ 1\end{array}\right] = \left[ \begin{array}{ccc}1&0&0&d \\\\ 0&1&0&h \\\\ 0&0&1&l \\\\ 0&0&0&1 \end{array}\right]\left[ \begin{array}{ccc}x\\\\ y \\\\ z \\\\ 1 \end{array}\right]$$

アフィン同士は、どう処理すればいいんだっけか。

3D回転の時は、x軸中心に回すやつと、y軸中心に回すやつの、行列の積を出して、使った覚えがある。



去年の年末、せっせと作っていた。

せっかくなので、アップして貼っておく。
しばらくしたら、もう少し、追加の情報を入れて、アップしなおす。

動画では、リングをくるくる回すために…

  1. リングを作るために、各点のx,y座標を、(cosθ, sin θ) 0<θ<2*πで作って、
  2. 回転行列(変数θが、for構文で動く)との行列のかけ算を計算して、回転後のx,y座標を作って、
  3. 計算結果を3Dプロットする、

という作戦で作った。

最後の、グチャーとしたのが、二軸回転のトラジェクトリ。
これは、二つの回転行列の、行列のかけ算で出した。

二軸の回転角θ、φは、別々に設定した。

上記は、たまたま、同種のアフィン行列(y,z,軸の回転行列)二つの組み合わせだった。

別種類の行列同士、(例えば、今回の”回転行列”と、”平行移動”)でも、行列のかけ算をするべきだろうか?

でも、めんどくせえ。
この場合、単純に行列の和でいいのかな?そっちなら、簡単だし、そっちであってほしい。
というか、インデキシングした配列に代入すればいっか…な?
どうでしょう。
「どのみち結果は、一緒」って事も、ありそう。

両方試せばいっか!

結局は、lambda式の中の、行列のかけ算の右側の引数を、うまい具合に設定できれば、なんだっていい。

複数のアフィン変換をする場合、それらの行列同士の積をするべきだろうが、

  • 代入
  • アフィン同士の積

の両方を試した。

代入で試してみた。

$$\left[ \begin{array}{ccc} x'\\\\ y' \\\\ z' \\\\ 1\end{array}\right] = \left[ \begin{array}{ccc}1.2&0&0&\fbox{-0.1} \\\\ 0&0.5&0&\fbox{0.1} \\\\ 0&0&1.5&\fbox{0} \\\\ 0&0&0&1 \end{array}\right]\left[ \begin{array}{ccc}x\\\\ y \\\\ z \\\\ 1 \end{array}\right]$$

四角で囲ったセルに、移動させたい数値を代入した。python的には下の様になった。

aff = np.diag([1.2, 0.5, 1.5, 1])
aff[0][3] = -0.1#x軸担当
aff[1][3] = 0.1#y軸担当
ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), aff)

ん?動かない。

aff_shift_x_-0p1.png


極端に大きなな数字いれてやれ。

aff[0][3] = -30

aff_shift_x_-30y_0.png

あ、動いてますね。
わーい。

x,y,z3軸いじれるので、書き直すけど、

aff = np.diag([1.2, 0.5, 1.5, 1])
aff[0][3] = 0#x軸担当
aff[1][3] = 0#y軸担当
aff[2][3] = 0#z軸担当
ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), aff)

こういう事ですわ。

  • xの場合
コードの変えた部分行列結果
なし$$\left[ \begin{array}{ccc}1.2&0&0&\fbox{0} \\\\ 0&0.5&0&0 \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]$$aff_shift_x_0y_0.png
aff[0][3] = -30#x軸担当$$\left[ \begin{array}{ccc}1.2&0&0&\fbox{-30} \\\\ 0&0.5&0&0 \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]$$aff_shift_x_-30y_0.png
aff[0][3] = -60#x軸担当$$\left[ \begin{array}{ccc}1.2&0&0&\fbox{-60} \\\\ 0&0.5&0&0 \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]$$aff_shift_x_-60y_0.png
  • yの場合
コードの変えた部分行列結果
なし$$\left[ \begin{array}{ccc}1.2&0&0&0 \\\\ 0&0.5&0&\fbox{0} \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]$$aff_shift_x_0y_0.png
aff[1][3] = 2#y軸担当$$\left[ \begin{array}{ccc}1.2&0&0&0 \\\\ 0&0.5&0&\fbox{2} \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]$$aff_shift_x_0y_2.png
aff[1][3] = 5#y軸担当$$\left[ \begin{array}{ccc}1.2&0&0&0 \\\\ 0&0.5&0&\fbox{5} \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]$$aff_shift_x_0y_5.png
  • zの場合
コードの変えた部分行列結果
なし$$\left[ \begin{array}{ccc}1.2&0&0&0 \\\\ 0&0.5&0&0 \\\\ 0&0&1.5&\fbox{0} \\\\ 0&0&0&1 \end{array}\right]$$aff_shift_x_0y_0.png
aff[2][3] = 2#z軸担当$$\left[ \begin{array}{ccc}1.2&0&0&0 \\\\ 0&0.5&0&0 \\\\ 0&0&1.5&\fbox{2} \\\\ 0&0&0&1 \end{array}\right]$$aff_shift_x_-0y_0z_2.png
aff[2][3] = 5#z軸担当$$\left[ \begin{array}{ccc}1.2&0&0&0 \\\\ 0&0.5&0&0 \\\\ 0&0&1.5&\fbox{5} \\\\ 0&0&0&1 \end{array}\right]$$aff_shift_x_-0y_0z_5.png

行列のかけ算は、代入と答えが違う…訳ないだろう。

「もう、動いてくれたし、別に、やらんでもいいかなぁ〜」とか思っちゃったけど、やる。

# aff1拡大・縮小行列
aff1 = np.diag([1.2, 0.5, 1.5, 1])

# aff2平行移動行列
aff2=np.diag([1., 1., 1., 1.])
aff2[0][3] = -20.0
aff2[1][3] = 0.5

$$\left[ \begin{array}{ccc}1.2&0&0&0 \\\\ 0&0.5&0&0 \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]\cdot\left[ \begin{array}{ccc}1&0&0&-20 \\\\ 0&1&0&0.5 \\\\ 0&0&1&0 \\\\ 0&0&0&1 \end{array}\right]$$

>>> np.dot(aff1, aff2)
array([[ 1.2 , 0. , 0. , -24. ],
[ 0. , 0.5 , 0. , 0.25],
[ 0. , 0. , 1.5 , 0. ],
[ 0. , 0. , 0. , 1. ]])

!?、代入と、微妙に答えが違うぜ。

いや、行列のかけ算は、順番が大事なんだぜ。
行列の掛け算では、左右入れ替えると、

  • 答えが変わる
  • そもそも、計算できなくなる。

とか、起こり得る。
今回は、正方行列なので、左右を入れ替えても、計算はできる。
しかし、結果は変わる。

>>> np.dot(aff2,aff)
array([[ 1.2, 0. , 0. , -20. ],
[ 0. , 0.5, 0. , 0.5],
[ 0. , 0. , 1.5, 5. ],
[ 0. , 0. , 0. , 1. ]])
もどった。というか、代入した場合と、同じになった。
逆の場合も、間違いではないんだろう。処理の順番によって、結果が変わるだけ。
下の表に、行列のかけ算をまとめた。
左側右側答え
$$\left[ \begin{array}{ccc}1.2&0&0&0 \\\\ 0&0.5&0&0 \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]$$$$\left[ \begin{array}{ccc}1&0&0&-20 \\\\ 0&1&0&0.5 \\\\ 0&0&1&0 \\\\ 0&0&0&1 \end{array}\right]$$$$\left[ \begin{array}{ccc} 1.2&0&0&-24 \\\\ 0&0.5&0&0.25 \\\\ 0&0&1.5&0 \\\\ 0&0&0&1  \end{array}\right]$$
$$\left[ \begin{array}{ccc}1&0&0&-20 \\\\ 0&1&0&0.5 \\\\ 0&0&1&0 \\\\ 0&0&0&1 \end{array}\right]$$$$\left[ \begin{array}{ccc}1.2&0&0&0 \\\\ 0&0.5&0&0 \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]$$$$\left[ \begin{array}{ccc} 1.2&0&0&-20 \\\\ 0&0.5&0&0.5 \\\\ 0&0&1.5&0 \\\\ 0&0&0&1 \end{array}\right]$$

この順番が、どこから来るのか分からないけど、順番注意なようですな。
拡大縮小→平行移動
平行移動→拡大縮小
では、意味が違うようです。

表の、下側の結果は、代入したのと、同じになっているし、代入した結果で事足りてる。

自分の中の結論としては、

  • 代入の方が、行数少ない。
  • 複数組み合わせて難しいことするなら行列のかけ算使う。しかし順番注意。

アフィンで、軸を動かせれば、何かしら使えそう

アフィンの部分を利用すれば、

  • ゆらゆら軸を揺らしながら、プロットしたり、
  • どっくんどっくん軸を脈動させながら、プロットしたり、

特殊なアニメーションができそう。動かない方が見やすいけど。
どういうシーンで、それを使うのかは、別にして…。
これからも、是非、是非、無駄なことをやっていきたいと思う。

ん…?結局グラフは、完成させないまま終わってしまった。凄い!8000文字近く書いた。

posted by yuchan at 07:00 | Comment(0) | python
この記事へのコメント
コメントを書く
お名前:

メールアドレス:

ホームページアドレス:

コメント: