ベクトル・行列の演算を自分で作る

目標

3次元のベクトル演算と4x4の行列の積を自分で作り、それで図形を移動し、回転し、縮小してみます。このラボが終われば、変換行列を読み書きでき、掛ける順序がなぜ重要なのかを、手で確認した状態になります。

なぜ重要なのか

3次元の計算なのに行列が4x4なのは、移動が乗算ではないからです。座標を1つ増やして(x, y, z, 1)と書くと、最後の列がその1と掛け合わされて結果に足され、そうして平行移動が乗算の中に入ります。これを同次座標と呼び、wが1なら点、0なら方向という区別が、ここから出てきます。

そして、行列で表現しておけば、複数の変換を先に掛けて1つにできます。GPUが、頂点ごとにMVP行列を1つだけ掛ける理由です。ただし、行列の積は交換法則が成り立たないので、順序を変えると、物体がその場で回る代わりに、軌道を回ります。この違いを数字で一度見れば、忘れません。

ステップ

/root/linalg/gfxlib.pyにツールボックスを置きます。
/root/linalg/vec.pyに、ベクトル演算を7つ作ります。
/root/linalg/mat.pyに、4x4の行列の積を作ります。
/root/linalg/xform.pyに、移動・拡大縮小・回転の行列を作ります。
掛ける順序を変えた結果を、/root/linalg/out/05-order.txtに書きます。
変換した四角形3つを、/root/linalg/out/square.pngに描きます。
外積で法線と面積を求めて、/root/linalg/out/07-normal.txtに書きます。

参考

すべてのファイルは、/root/linalgの中に置いてください。採点ツールが、そのフォルダーをimportのパスに入れて、関数を直接呼び出してみます。
絵を見るには、nohup python3 -m http.server 8080 -d /root/linalg/out &で立てたあと、Webプレビューでhttp://localhost:8080/square.pngを開いてください。
よくあるミスの1つ目は、mulのインデックスを取り違えて、転置された行列を作ってしまうことです。単位行列で試しても引っかからないので、移動行列で試してください。
よくあるミスの2つ目は、rotate_yの符号を、ほかの軸と同じにそろえてしまうことです。軸の循環の順序のために、y軸だけが逆です。

描画のツールボックスを置く

/root/linalg/gfxlib.pyを例のとおりに保存し、/root/linalg/check.pyでテストパターンを描いて、/root/linalg/out/00-check.pngを作ってください。パターンは、64x64の黒い背景に、(0,0)から(63,63)まで、白(255,255,255)の対角線を引き、そのあとに、(0,32)から(63,32)まで、赤(255,0,0)の横線を重ねて引いたものです。

このラボからは、PNGエンコーダーを作り直しません。最初のラボで手で作ったものと同じコードを、ツールとして配ります。ここで学ぶのは、ファイル形式ではないからです。

ラボのPodにはボリュームがないので、前のラボで作ったファイルが残っていません。そのため、ラボごとに、ツールボックスを置き直すことから始めます。

Canvas(w, h, bg)を作り、line(x0, y0, x1, y1, rgb)で2本の線を引いたあと、write_png(path)で保存してください。交差点(32,32)が赤になるには、横線をあとで引く必要があります。

3次元のベクトル演算

/root/linalg/vec.pyに、add(a,b)、sub(a,b)、scale(a,s)、dot(a,b)、cross(a,b)、length(a)、normalize(a)の7つの関数を作ってください。ベクトルは長さ3のタプルで、dotとlengthだけが実数を返します。

外積の定義を、符号まで正確に書く必要があります。

cross(a, b) = (a1*b2 - a2*b1, a2*b0 - a0*b2, a0*b1 - a1*b0)

cross((1,0,0), (0,1,0))が(0,0,1)になれば、符号は合っています。逆になった場合は、右手座標系が反転しているので、あとで照明が物体の後ろから当たっているように見えてしまいます。

normalizeは、長さが0のベクトルを受け取ることがあります。そのとき0で割らないように、元のベクトルをそのまま返してください。

4x4の行列の積

/root/linalg/mat.pyに、identity()、mul(a,b)、apply(m,v)の3つの関数を作ってください。行列は行優先(長さ4のリスト4つを入れたリスト)で、applyは、点(x,y,z)を列ベクトル(x,y,z,1)とみなして掛けたあと、(x,y,z,w)の4つの成分を返します。

乗算の定義は、out[r][c] = sum(a[r][k] * b[k][c] for k in range(4))です。インデックスの順序を一度取り違えると、転置された行列が出ますが、単位行列で試しても引っかかりません。単位行列は、転置しても自分自身だからです。

そのため、テストは非対称な行列で行ってください。移動行列のように、最後の列にだけ値があるものが良いです。

applyは、m[r][0]*x + m[r][1]*y + m[r][2]*z + m[r][3]*1を、r=0..3について計算します。

移動・拡大縮小・回転の行列

/root/linalg/xform.pyに、translate(tx,ty,tz)、scale(sx,sy,sz)、rotate_x(deg)、rotate_y(deg)、rotate_z(deg)の5つの関数を作ってください。角度は度(degree)で受け取り、右手座標系が基準です。

rotate_zの左上の2x2は、[[cos, -sin], [sin, cos]]です。こうすると、rotate_z(90)が、x軸(1,0,0)をy軸(0,1,0)へ送ります。

同じ規則で、rotate_xはyをzへ、rotate_yはzをxへ送ります。rotate_yだけ、符号が逆に入ります。m[0][2] = sin、m[2][0] = -sinです。軸の循環の順序がx→y→z→xだからですが、ここで符号を統一してしまうと、y軸の回転だけが逆に回ります。

math.radiansで、度をラジアンに変換してください。

掛ける順序が結果を変える

T = translate(2,0,0)、R = rotate_z(90)とおき、点(1,0,0)に、mul(T,R)とmul(R,T)をそれぞれ適用した結果を、/root/linalg/out/05-order.txtに、TR=x,y,zとRT=x,y,zの2行で書いてください。小数点以下6桁まで書きます。

列ベクトルを右に掛ける慣習では、右の行列が先に適用されます。mul(T, R)は、回転を先に行い、移動をあとに行うという意味です。

手で追ってみてください。Rが(1,0,0)をどこへ送るかを先に求め、そこにTを足します。逆の順序は、Tを先に適用してから、Rを掛けます。

一方は、物体がその場で回ったあとに移されたもので、もう一方は、原点から離れたあと、原点を基準に軌道を回ったものです。結果が大きく異なります。

変換した四角形を描く

/root/linalg/plot.pyで、256x256のPNGを/root/linalg/out/square.pngに書いてください。モデル座標の四角形は、(-1,-1,0)、(1,-1,0)、(1,1,0)、(-1,1,0)で、画面座標へは、sx = 128 + 60*x、sy = 128 - 60*yで移します。黒い背景に、四角形を3つ描きます。変換のないものは白(255,255,255)、rotate_z(30)をかけたものは赤(255,60,60)、mul(translate(1.2,-0.8,0), scale(0.5,0.5,1))をかけたものは緑(60,255,60)です。

各四角形は、4つの頂点を順番につないだあと、最後から最初の点に戻る、4本の線分です。gfxlib.Canvasのline(x0, y0, x1, y1, rgb)を使ってください。

画面座標へ移すとき、yに引き算が入っていることに注意してください。数学のyは上に、画像のyは下に向かって増えます。これを抜かすと、緑の四角形が上側に描かれます。

緑の四角形の変換は、縮小を先に行い、移動をあとに行います。mul(translate, scale)の右が先です。

外積で法線と面積を求める

/root/linalg/normal.pyで、2つの三角形の単位法線と面積を求めて、/root/linalg/out/07-normal.txtに、t1_normal=、t1_area=、t2_normal=、t2_area=の4行で書いてください。三角形1はA(0,0,0) B(2,0,0) C(0,3,0)、三角形2はA(1,1,1) B(2,1,1) C(1,1,3)で、法線は、cross(B-A, C-A)を正規化したものです。

外積の長さが、2つの辺が作る平行四辺形の面積なので、三角形の面積は、その半分です。つまり、area = length(cross(B-A, C-A)) / 2です。

法線の符号は、頂点を回る順序で決まります。cross(B-A, C-A)とcross(C-A, B-A)は、正反対の向きです。この順序が、あとで表面と裏面を分ける基準になるので、指示どおりに合わせておいてください。

値は%.6fで書き、ベクトルはカンマでつなぎます。例: t1_normal=0.000000,0.000000,1.000000