4x4行列演算とアフィン変換の定義(Rust製ラスタライザ)
点を描くところから始めるRust製ソフトウェアラスタライザ::4x4行列演算・4次元ベクトルを実装する
AC.icon 行列Mat4x4の実装
目次
Mat4x4, Vec4
Mat4x4, Mat4x4
回転行列の実装
任意の軸基準での回転の実装
スケーリング
並行移動
テスト
行列のデータ構造
e: [[f64; 4]; 4]で保持する。
行列積周りをさっぱりと定義したい
行列の積の定義ベタ書きはあんまりやりたくない
あと、こういう感じの記述がよく出てくるようになった（ベクトルの各要素に番号0, 1, 2, 3と振っているため）
code:rs
Vec4::new(
element(0),
element(1),
element(2),
element(3),
)
この書き方を抽象化したい
ということで、まずVec4に次のメソッドを実装
code:rs
/** i番目の要素がelement(i)であるようなベクトルを構成する。
* ただし、iについて0はx, 1はy, 2はz, 3はw座標をそれぞれ示す。
*/
pub fn construct<F>(element:F) -> Self
where
F:Fn(usize) -> f64
{
Vec4::new(
element(0),
element(1),
element(2),
element(3),
)
}
そうすると、さまざまなものの記述が楽になる。
二つのベクトルをとってそれぞれの要素の和を取るaddとかの記述がすごい楽になる。
code:rs
/** 要素ごとの積 */
pub fn hadamard(&self, rhs:&Self) -> Self { Vec4::bin(self, rhs, |a, b|{a * b}) }
/** 内積 */
pub fn dot(&self, rhs:&Self) -> f64 { self.hadamard(rhs).fold() }
impl ops::Add for &Vec4 {
type Output = Vec4;
fn add(self, rhs: Self) -> Self::Output {
Vec4::construct(|i| self.ei + rhs.ei) }
}
同様に、Mat4x4でも同じような関数を定義する。
code:rs
fn construct<F>(element:F) -> Mat4x4
where
F: Fn(usize, usize) -> f64
{
for r in 0..4 {
for c in 0..4 {
}
}
Mat4x4{e}
}
行・列ベクトルを取り出す関数も定義しておく。
（これもconstruct関数のおかげで定義は簡単に書ける）
ここ、ちょっとオーバーヘッドになりそう
（ベタ書きで実装すると比べると、わざわざVec4オブジェクトでラップしてるのがパフォーマンス上懸念）
コンパイラはVec4でラップするところも果たして展開してくれるだろうか...？appbird.icon
code:rs
/** 行列のr行目の行ベクトルを取り出す。 */
fn row(&self, r:usize) -> Vec4 {
Vec4::construct(|i| self.eri) }
/** 行列のc列目の列ベクトルを取り出す。 */
fn col(&self, c:usize) -> Vec4 {
Vec4::construct(|i| self.eic) }
すると、行列積の実装本体はワンラインで書ける（し、行列積の意味もすごく取りやすくなる。）
code:rs
impl Mul<Vec4> for Mat4x4 {
type Output = Vec4;
fn mul(self, rhs: Vec4) -> Self::Output {
Vec4::construct(|r| self.row(r).dot(&rhs))
}
}
impl Mul for Mat4x4 {
type Output = Mat4x4;
fn mul(self, rhs: Mat4x4) -> Self::Output {
Mat4x4::construct(|r, c| self.row(r).dot(&rhs.col(c)))
}
}
結構満足appbird.icon
ゼロコスト抽象化がうまく働いてくれているなら、こう抽象化してもさほどコストがかさまないと期待する
回転行列の実装
任意の軸基準での回転の実装
いくつかの方針は考えられる。
---> 今からもういっちょクォータニオンを実装するのはちょいしんどい
行列で済むなら済ませたい
ちゃんとやるなら導出すべきだけど、線形代数の理解が目的ではないので今回はこれを認める。
実装として
fn rotation(axis:Vec4, theta:f64) -> Mat4x4
ベタ書きでもいいけど、せっかくなら構造的に描きたい
まず、w要素に関しては計算に介入させないようにする
(r, c)が(3, 3)なら$ 1
それ以外の3行目・3列目の要素は$ 0にする
それ以外の部分
code:tex
\left[\begin{array}{ccc}
c+n_1^2(1-c) & n_1 n_2(1-c)-n_3 s & n_1 n_3(1-c)+n_2 s \\
n_2 n_1(1-c)+n_3 s & c+n_2^2(1-c) & n_2 n_3(1-c)-n_1 s \\
n_3 n_1(1-c)-n_2 s & n_3 n_2(1-c)+n_1 s & c+n_3^2(1-c)
\end{array}\right]
まず、対角成分は
$ a_{ii} = c + n_in_i (1 - c)
非対角成分は
$ a_{ij} = g\cdot n_{3-i-j}s + n_in_j (1 - c)
符号$ gは$ (i, j)が$ (2, 1, 0, 2, 1, 0)の連続部分列であれば$ 1、そうでなければ$ -1
code:rs
/** Rodriguesの回転行列を使ってaxis軸周りに角度thetaだけ点を回転させる行列を作る */
fn rotation(axis:Vec4, theta: f64) -> Mat4x4 {
let (cos, sin) = (f64::cos(theta), f64::sin(theta));
let lcos = 1. - cos;
Self::construct(
|r, c|
if (r, c) == (3, 3) { 1. }
else if r == 3 || c == 3 { 0. }
else if r == c {
} else {
let sign = if
(r, c) == (2, 1) ||
(r, c) == (1, 0) ||
(r, c) == (0, 2)
{ 1. } else { -1. };
}
)
}
スケーリング・並行移動
construct関数のおかげでなんとかappbird.icon
code:rs
/** s倍するアフィン変換行列 */
pub fn scale(s:Vec4) -> Mat4x4 {
Self::construct(
|r, c|
if r != c { 0. }
// # FIXED 4になってた
else if r != 3 { s.i(r) }
else { 1. }
)
}
/** v方向に移動させるアフィン変換行列を作る */
pub fn translate(v:Vec4) -> Mat4x4 {
Self::construct(
|r, c|
if r == c { 1. }
else if c == 3 { v.i(r) }
else { 0. }
)
}
適当に小さな点を打つメソッドをCanvasに実装した。その上で、点群と線分がアフィン変換を受けて変化するかを見る。
code:rs
pub fn draw_point(&mut self, center:&Point2, color: &Color) {
let half_width = 3;
for x in -half_width .. half_width {
for y in -half_width .. half_width {
if x*x + y*y < half_width*half_width {
self.draw_pixel(&Point2::new(center.x + x, center.y + y), color);
}
}
}
}
毎フレーム画面をリセットすることにするappbird.icon
長方形の頂点$ (0.5, 0, 0), (0.5, 0.25, 0), (0, 0, 0), (0, 0.25, 0)にある点を描画する。
https://scrapbox.io/files/66f965186563ec001d6d5142.png
$ (0.5, 0.25)だけ移動する。
https://scrapbox.io/files/66f97aa83aeab5001d411b60.png
constructの実装が間違ってることに気づいた
あと、渡すw座標が1.じゃないといけない
$ z座標中心で$ \pi/6\cdot tだけ回転する。
https://scrapbox.io/files/66f97bc393a251001dd0119a.png
回転方向は時計回り
これは、このシステムが$ y軸下向き正にとっているため
下向き負にとった場合は反時計回りになっている必要がある
$ (0.8, 0.5)倍にスケーリングする。
https://scrapbox.io/files/66f97bfc4431ba001c1e8799.png
$ (0.8,0.5)倍にスケーリングして、$ \pi/6\cdot tだけ回転させた後、$ (0.5, 0.25)だけ移動する。
https://scrapbox.io/files/66f97d8f54a623001cc5c2c6.png
この図で言う右下の頂点中心に回転できていればOK