座標変換と三次元空間のポリゴン(Rust製ラスタライザ)

from 点を描くところから始めるRust製ソフトウェアラスタライザ

点を描くところから始めるRust製ソフトウェアラスタライザ::三次元空間にあるポリゴンを描画する

https://gyazo.com/919410ffe20f074cf951c8274dbd78ac

4x4行列演算とアフィン変換の定義(Rust製ラスタライザ)で、アフィン変換を実装できたので、カメラを実装していく。

ただし、ここら辺はバグらせやすい要因になるので段階的に開発したいappbird.icon

Ray Tracing in One Weekendの作法に則っている

最終：視点$ \bm pから視線$ \bm Lの方向を向いて、$ \bm uを上とするカメラを作る。

参考資料

床井先生のこちらの資料を参考にさせていただく...

○○くんのために一所懸命書いたものの結局○○くんの卒業に間に合わなかったGLFW による OpenGL 入門

特に、p.114以降 (pdfビューア上ではp.122以降)

p.153

用語などはこちらの資料を参考にする。

検証

$ \bm pを$ (\cos(2\pi t/5), \sin(2\pi t/5), 0)とする。

$ \bm pを$ (0, 0, \sin(2\pi t/5))とする。

$ \bm p = (0, 0, 6)、$ \bm Lを$ (\sin(2\pi t/5), 0, \cos(2\pi t/5))とする。

くるくるしてると良い

$ \bm uを$ (\cos(2\pi t/5), \sin(2\pi t/5), 0)とする。

画面がぐるぐるしていると良い

座標系は？

右手座標系をとる。

カメラが$ (0, 0, 0)を$ (0, 0, z < 0)から見ているとき、上側正

カメラは基本$ +z方向$ (0, 0, 1)を向くとする

この場合、$ xy平面を右向きを$ +x、下向きを$ +yとみる。

通常多くの場合は$ -z方向を向くとするが、これは$ x,y平面を右向き$ +x、上向き$ +yと見れるためであるappbird.icon

（後から$ -z方向に切り替える）

AC.icon 透視変換: (0, 0, 0)から(0, 0, 1)の方向を向いて、(0, 1, 0)を上とするカメラを作る。

点群を正しく変換するジオメトリパイプライン (Geometry Pipeline)の実装が必要

モデルの頂点座標とモデルの座標・角度を持つActorを定義する

モデル（三角ポリゴン）の中心座標$ \bm p

回転（回転軸と角度$ \bm \theta）

自身のモデル変換行列を返す

回転$ \bm\thetaを行った$ R後、$ \bm pだけ並行移動する$ T

行列$ TRを求めれば良い

2. ビュー変換

今は単位行列としてよい

3. 透視変換

資料のPDFとしてはビューア表示においてp.161ぐらいappbird.icon

直交投影ではなく、透視投影を行う。

farを0.1 m, nearを100 mとして投影する。

視錐体

標準視体積に見える範囲の物体を持ってくるappbird.icon

せん断変形 --> 今回は視錐体はビュー変換された地点で直線$ (x, y) = (0, 0)を通っているものと仮定する

ので、ここは単位行列でよい

同時座標系での透視変換を考える

画像乾板に映る像の大きさを考えることになる。

これは、幾何学的に相似を考えればすぐに求められる。

https://scrapbox.io/files/66faa3b9c371ca001ccd9eab.png

$ zで割る分に関しては、同次座標系の$ wへ$ z値を載せればいい

さらに、画像乾板の端っこにあたる座標値は$ -1, 1に正規化しておきたいので

画像乾板の幅$ Wと高さ$ Hで割っておく。

（ピクセル幅ではなく、mで。）

ということで、求める行列は、$ z座標ぶんを除いてここまでわかる。

code:tex

\begin{pmatrix}

x' \\ y' \\ z' \\ w'

\end{pmatrix}

\begin{pmatrix}

near/W & 0 & 0 & 0 \\

0 & near/H & 0 & 0 \\

* & * & * & * \\

0 & 0 & 1 & 0

\end{pmatrix}

\begin{pmatrix}

x \\ y \\ z \\ w

\end{pmatrix}

透視深度

透視投影後の像の深度（カメラに対して、物体がどれぐらいの奥行き方向の距離に存在するか）のことを透視深度という

透視投影では、投影像の大きさが深度 (z 値) に反比例します。透視投影後の深度にも、元の深度の逆数を用います。

この深度は$ \lbrack -1, 1\rbrackの実数で表される。

https://scrapbox.io/files/66faa3f0bbbd59001dde1f67.png

要件

1. 前方面のz値（$ z = N）が-1、後方面のz値（$ z = F）が1になってほしい

2. $ zに反比例してほしい。

これはなぜ？appbird.icon

透視投影後の位置の精度は前方面に近いほど高く、前方面から離れるほど低くなります。このような深度を透視深度といいます。

なるほど、前方面まわりの深度の精度を確保したいからか

$ f(z) = \frac{C_1}{z} + C_2という関数の形を仮定して、このような関数を探してみる。

$ -1 = \frac{C_1}{N} + C_2, $ 1 = \frac{C_1}{F} + C_2を境界条件として満たす必要あり。

$ -N = C_1 + C_2N, $ F = C_1 + C_2F

$ F + N = (F - N) C_2

$ \frac{F + N}{F - N} = C_2

さらに

$ -1 - \frac{ F + N }{ F - N } = \frac{C_1}{N}

$ \frac{ N - F - F - N }{ F - N } = \frac{C_1}{N}

$ - \frac{ 2F }{ F - N } = \frac{C_1}{N}

$ - \frac{ 2FN }{ F - N } = C_1

従って、

$ f(z) = -\frac{2FN}{F-N}\frac{1}{z} + \frac{F + N}{F-N}

同時座標系の$ w値に$ z値が乗ってるので、それも考慮すると

$ zf(z) = -\frac{2FN}{F-N} + \frac{F + N}{F-N}z

https://gyazo.com/c75e2cdbbabff1267c1300f0b35e64a9

ということで、正規化まで含めて求める行列はこうなる。

code:tex

\begin{pmatrix}

x' \\ y' \\ z' \\ w'

\end{pmatrix}

\begin{pmatrix}

N/W & 0 & 0 & 0 \\

0 & N/H & 0 & 0 \\

0 & 0 & \frac{F + N}{F - N} & -\frac{2FN}{F-N} \\

0 & 0 & 1 & 0

\end{pmatrix}

\begin{pmatrix}

x \\ y \\ z \\ w

\end{pmatrix}

ふつう、$ wは$ 1であることに注意したいappbird.icon

実装

モデルからモデル変換行列を提供する

code:rs

pub struct Actor {

pub vertices:Vec4; 3,

position:Vec4,

axis:Vec4,

theta:f64,

}

impl Actor {

fn model_conversion(&self) -> Mat4x4 {

Mat4x4::translate(&self.position)

* Mat4x4::rotation(&self.axis, self.theta)

}

カメラから透視変換行列を提供する

code:rs

pub fn perspective_conversion(&self) -> Mat4x4 {

let n = self.near;

let f = self.far;

let w = self.width;

let h = self.height;

Mat4x4::from_array([

n / w, 0., 0., 0. ,

0., n / h, 0., 0. ,

0., 0., (f+n)/(f-n), -2.*(f*n)/(f - n),

0., 0., 1., 0. ,

])

}

Actorの頂点をモデル変換行列・透視変換行列に通す。

code:rs

pub fn shot(&mut self, canvas:&mut Canvas, actors:Vec<Actor>) {

let perspective = self.perspective_conversion();

let view = Mat4x4::identity();

let pv = perspective*view;

for actor in actors {

let model = actor.model_conversion();

let pvm = &pv * &model;

let projected:Vec4Project; 3 = [

Vec4Project(&pvm * &actor.vertices0),

Vec4Project(&pvm * &actor.vertices1),

Vec4Project(&pvm * &actor.vertices2),

];

self.draw_triangle(canvas, &projected0, &projected1, &projected2, &actor.color);

}

テストappbird.icon

こんな感じの空間を撮ってみる

データ作成にはを使用

https://scrapbox.io/files/66fa0189b60f26001dc50d77.png

https://scrapbox.io/files/66fa0170bdf41c001c3cdcde.png

code:rs

let triangles = [

Vec4::new3d(0., 0., 0.), Vec4::new3d(3., 1., 2.)/5, Vec4::new3d(0., -3., 4.)/5,

Vec4::new3d(0., 4., -5.)/5, Vec4::new3d(-3., -1., -2.)/5, Vec4::new3d(0., 3., 4.)/5,

];

let blue = Vec4::new3d(0.2, 0.2, 0.9);

let green = Vec4::new3d(0.2, 0.9, 0.2);

let red = Vec4::new3d(0.9, 0.2, 0.2);

いけた！！appbird.icon

https://gyazo.com/d76c5de3cb9d2935f03e579cc8f36736

AC.icon 視点$ \bm pから視線$ \bm Lの方向を向いて、$ \bm uを上とするカメラを作る。

まず、視点$ \bm pから。

視点$ \bm pから位置$ \bm xにある点を見ている

== 視点$ 0から位置$ \bm x - \bm pにある点を見ている

code:camera.rs

let perspective = self.perspective_conversion();

let view = self.view_conversion();

let pv = perspective*view;

code:camera.rs

pub fn view_conversion(&self) -> Mat4x4 {

Mat4x4::translate(&-(&self.location))

}

Vec4に同次座標に関する調節

code:Vec4.rs

pub fn newvec(x: f64, y: f64, z: f64) -> Self {

Vec4 { e: x, y, z, 0. }

}

pub fn newpoint(x: f64, y: f64, z: f64) -> Self {

Vec4 { e: x, y, z, 1. }

}

ちょっと三角形を大きくして...

code:main.rs

let world:Vec<Actor> = vec![

Actor{

vertices: Vec4::newpoint(0., 0., 0.), Vec4::newpoint(3., 1., 2.)/3., Vec4::newpoint(0., -3., 4.)/3.,

position: Vec4::newpoint(0., 0., 0.),

axis: Vec4::newvec(0., 0., 0.),

theta: 0.0,

color: blue.clone(), green.clone(), red.clone()

Actor{

vertices: Vec4::newpoint(0., 4., -5.)/3., Vec4::newpoint(-3., -1., -2.)/3., Vec4::newpoint(0., 3., 4.)/3.,

position: Vec4::newvec(0., 0., 0.),

axis: Vec4::newvec(0., 0., 0.),

theta: 0.0,

color: green.clone(), blue.clone(), blue.clone()

}

];

カメラの位置を変える！

code:main.rs

while canvas.update()? {

let t = stopwatch.elapsed_as_sec();

camera.location = Vec4::newpoint(f64::cos(k*t), f64::sin(k*t), 0.)*2.;

テスト

AC.icon $ \bm pを$ (\cos(2\pi t/5), \sin(2\pi t/5), 0)とする。

https://gyazo.com/d0ec56bf5a4441cf8eadf6b3c4587ddd

z-bufferが必要なことがわかるappbird.icon

$ \bm pを$ (0, 0, \sin(2\pi t/5))とする。

WA.icon こっちは動かない...なんでだ....appbird.icon

透視変換を修正するまでの道のり < 点からはじめるラスタライザ

AC.icon bugfix

深度クリッピング < 点からはじめるラスタライザ

ここまでで

物体は$ z \in (-6, 6)のどこかに存在する。

カメラを視点$ (0, 0, 8)から$ (0, 0, 1): z+方向に向けると

https://gyazo.com/1900f8b29e11bf8bf78196b75c320efa

カメラを視点$ (0, 0, -8)から$ (0, 0, 1): z+方向に向けると

https://gyazo.com/6d6a6bc0e3df9ede467dae403ba0f425

となり、確かに深度クリッピングができていると言える。

ViewPortの方で、空間は上向き正としたい。が、画面は下向き正なのが問題appbird.icon

視線$ \bm L、上方向$ \bm uに取り組む。

わからないので資料を見るappbird.icon

$ \bm Lはカメラ座標系の$ z軸の向き$ \bm e_zに対応し、$ Lと$ uは線型独立

$ \bm uと$ \bm Lは必ずしも直交するとは限らない。

が、$ uはy軸の要素を含んでいる

したがって、$ \bm e_x = \bm u \times \bm Lはちょうど$ x軸の向きに対応するとすればよい。

再度$ \bm e_z = \bm e_x \times \bm Lとすれば、直交した$ z軸が得られる。

さて、標準基底を正規直交基底$ \{(\bm e_{x}, \bm e_{y}, \bm e_{z})\}からなる座標系へ写像する場合、このように記述できる。

code:latex

R = \left(

\begin{array}{c}

\begin{array}{c}\\ \bm e_x\\ \\ 0\end{array}

\begin{array}{c}\\ \bm e_y\\ \\ 0\end{array}

\begin{array}{c}\\ \bm e_z\\ \\ 0\end{array}

\begin{array}{c}0\\ 0\\ 0\\ 1 \end{array}

\end{array}

\right)

今回は基底$ \{\bm e_{x}, \bm e_{y}, \bm e_{z}\}を持つ座標系上で表現された点を、標準基底に持っていく必要がある。

したがって、$ Rが回転行列であることから$ R^{-1} = R^Tを計算すればよい。

というわけで三次元の外積、ノルムを定義する。

code:util/vec4.rs

pub fn cross3d(&self, rhs:&Self) -> Vec4 {

let (x1, y1, z1) = (self.x(), self.y(), self.z());

let (x2, y2, z2) = (rhs.x(), rhs.y(), rhs.z());

Vec4::new(

y1*z2 - z1*y2,

z1*x2 - x1*z2,

x1*y2 - y1*x2,

self.w() * rhs.w()

)

}

pub fn normalized3d(&self) -> Vec4 {

let (x, y, z) = (self.x(), self.y(), self.z());

let n = f64::sqrt(x*x + y*y + z*z);

Vec4::new(x/n, y/n, z/n, self.w())

}

続けて今回必要な行列$ R^Tを召喚する

code:util/Mat4x4.rs

pub fn transposed_basis(e1:&Vec4, e2:&Vec4, e3:&Vec4) -> Mat4x4 {

Mat4x4::from_array([

e1.i(0), e1.i(1), e1.i(2), 0.,

e2.i(0), e2.i(1), e2.i(2), 0.,

e3.i(0), e3.i(1), e3.i(2), 0.,

0., 0., 0., 1.,

])

}

これをview_conversionに読み込む

code:camera.rs

pub fn view_conversion(&self) -> Mat4x4 {

let ez = &self.look;

let ex = ez.cross3d(&self.up).normalized3d();

let ey = ez.cross3d(&ex);

Mat4x4::translate(&-(&self.position))

*Mat4x4::transposed_basis(&ex, &ey, ez)

}

$ \bm p = (0, 0, 8)とし、$ \bm L, \bm uを$ (\sin(2\pi t/5), 0, \cos(2\pi t/5)), (0, 1, 0)とする。

くるくるしてると良い

https://gyazo.com/4178d8cef686564b3c9682155783a4ba

なーんか回ってるようには見えませんねえ...？appbird.icon

遠くに行きすぎている感じがするかな、距離感が変動しているように感じられるappbird.icon

側面側だと近く、真ん中側だと遠くに見えてくる

カメラの側面の設定はこんな感じだし、別に変ではないと思うんだけどなあ？

code:camera.rs

let near = 0.1;

let width = 2. * near * f64::tan(PI /4.);

Camera{

near: 0.1,

far: 100.,

width,

height: width * aspect,

position: Vec4::newpoint(0., 0., 6.),

look: Vec4::newvec(0., 0., -1.),

up: Vec4::newvec(0., 1., 0.),

}

あぁ、よく考えたら$ up, ezの順が違いますね...。

code:camera.rs

pub fn view_conversion(&self) -> Mat4x4 {

let ez = &self.look;

let ex = &self.up.cross3d(&ez).normalized3d();

let ey = ez.cross3d(&ex);

Mat4x4::transposed_basis(&ex, &ey, &ez) * Mat4x4::translate(&-(&self.position))

}

いけた！

https://gyazo.com/32397bf7e66e35eabef48059748a7a02

$ \bm uを$ (\cos(2\pi t/5), \sin(2\pi t/5), 0)とする。

https://gyazo.com/94fa1106a441366d2862dd0e021ce70e

まさかの頂点属性関連にバグappbird.icon

三角形の描画においてy基準ソートしていたのが原因か

まあでもカメラとしてはこれで良さそうだappbird.icon

$ \bm p, \bm Lを$ 5(\cos(2\pi t/5),0, \sin(2\pi t/5)), -(\cos(2\pi t/5), 0, \sin(2\pi t/5))とする。

原点のポリゴンを周りからぐるっと見つめる感じになればOK

code:rs

vertices: Vec4::newpoint(0., 0., 0.), Vec4::newpoint(3., 1., 2.), Vec4::newpoint(0., -3., 4.),

vertices: Vec4::newpoint(0., 0., 0.), Vec4::newpoint(-3., -1., -2.), Vec4::newpoint(0., 3., -4.),

頂点をこのように置き直すappbird.icon

https://gyazo.com/669ae920336cc47ffdd56394a7442afd

AC.icon 尻拭い二つappbird.icon

1. まさかの頂点属性関連にバグappbird.icon

y値ソートが原因。元の配列を残さなきゃ

インデックスをソートするように

スキャンライン法をもう一度採用

Barycentric座標

code:canvas/triangle.rs

pub fn draw_triangle(

&mut self,

points: Vec4Screen; 3,

color: &Color; 3

) {

// y基準でソート

// :

if z_segment.and(&bound_z).is_empty() { return; }

// 辺を求める

let mut points_idx:usize; 3 = 0, 1, 2;

points_idx.sort_by(|a, b| pointsa.clone().y().partial_cmp(&pointsb.clone().y()).unwrap_or(std::cmp::Ordering::Equal));

let bottom, middle, top = [

points[points_idx0].to_point2(),

points[points_idx1].to_point2(),

points[points_idx2].to_point2()

];

let lines = [

Line::new(bottom, middle),

Line::new(middle, top),

Line::new(bottom, top),

];

// Barycentric座標

let bound_w = ClosedInterval::between(0.,1.);

let area_abc = area(&points0, &points1, &points2);

if area_abc.abs() < 1e-6 { return; }

let inv_abc = 1./area_abc;

for y in &y_segment.and(&bound_y) {

let edge = if y < middle.y { &lines0 } else { &lines1 };

let i_edge1 = lines2.across_y(y);

let i_edge2 = edge.across_y(y);

let x_segment = i_edge1.or(i_edge2);

for x in &x_segment.and(&bound_x) {

// :

}

こうなると

code:main.rs

camera.position = Vec4::newpoint(0., 0., 8. + 2.*f64::cos(k*t)) ;

camera.look = Vec4::newpoint(0., 0., -1.) ;

camera.up = Vec4::newvec(f64::cos(k*t), f64::sin(k*t), 0.) ;

これで実行してみると

https://gyazo.com/dc4cbe3f53a1c07fdeba3387c3002d72

2. -z方向にカメラを向けたい。

upは本当に画面のupでしょうかappbird.icon

code:camera.rs

position: Vec4::newpoint(0., 0., 8.),

look: Vec4::newvec(0., 0., -1.),

up: Vec4::newvec(0., 1., 0.),

code:main.rs

vertices: Vec4::newpoint(0., 0., 0.), Vec4::newpoint(1., 1., 1.), Vec4::newpoint(0., -3., 4.)

color: red.clone(), red.clone(), red.clone()

vertices: Vec4::newpoint(0., 0., 0.), Vec4::newpoint(-3., -1., -2.), Vec4::newpoint(0., 3., -4.)

color: blue.clone(), blue.clone(), blue.clone()

というわけでこの設定の元で動かしてみる。

https://gyazo.com/2e4bff8fe20b9cb2a3b86ef7927ed810

今日は、座標変換の罠について、知っておこう

https://gyazo.com/c82d60c510b36cd9c65d3c92a131147d

うーん、逆だね。

それじゃ！

$ z方向について思い当たる問題二つ

1. ビュー座標系においてカメラが$ +z方向を向くことを仮定したこと

2. Vec4Projectionのinto_screenメソッドにおいて、$ yの向きが逆転していることを考慮していない点

1. ビュー座標系においてカメラが$ +z方向を向くことを仮定したこと

https://scrapbox.io/files/66faa3f0bbbd59001dde1f67.png

本当は$ -z方向をむくように、$ f(x) = C_1/z + C_2に対して、$ f(-N) = -1, f(-F) = 1の境界条件を設定しないといけない

これは$ +zの時に求めた$ fに対して、$ f(-z)がそれを満たす解である。

まとめると$ (x', y', z')は

$ x' = \frac{N}{-z\cdot W}x, y' = \frac{N}{-z\cdot H}y

$ z' = f(-z) = \frac{2FN}{F-N}\frac{1}{z} + \frac{F + N}{F-N}

($ -zで割る必要がある点に注意。$ z < 0を見ているため、幾何的変換を成り立たせるためにも$ -z > 0で割らないと反転してしまう。)

したがって、割る操作を$ w値に載せることで実現すると

$ w'x' = \frac{N}{W}x, w'y' = \frac{N}{H}y

$ w'z' = -\frac{2FN}{F-N} - \frac{F + N}{F-N}z

$ w' = -z

code:tex

\begin{pmatrix}

x' \\ y' \\ z' \\ w'

\end{pmatrix}

\begin{pmatrix}

N/W & 0 & 0 & 0 \\

0 & N/H & 0 & 0 \\

0 & 0 & \frac{F + N}{F - N} & -\frac{2FN}{F-N} \\

0 & 0 & -1 & 0

\end{pmatrix}

\begin{pmatrix}

x \\ y \\ z \\ w

\end{pmatrix}

となると、視線$ \bm Lはビュー座標系では$ +zと逆の方向を向いていることになる。

というわけで次のように修正

code:rs

pub fn view_conversion(&self) -> Mat4x4 {

let ez = -&self.look; // -つける

let ex = &self.up.cross3d(&ez).normalized3d();

let ey = ez.cross3d(&ex);

Mat4x4::transposed_basis(&ex, &ey, &ez) * Mat4x4::translate(&-(&self.position))

}

これで動かしてみる

https://gyazo.com/d4be0e56d25562c701ef17b7e0009cc7

https://gyazo.com/3ba59c221b69fedbf710545ebd34f704

近づいてきた。ただ上方向に反転しているのでもう一題。

2. Vec4Projectionのinto_screenメソッドにおいて、$ yの向きが逆転していることを考慮していない点

screenはyは下向き正だが、ワールド座標ではのyは上向き正としたい。

ワールド上でyの上側に家の屋根がある場合、何も挟まないと画面下側に屋根が映ってしまう。(---> y値反転が必要)

というわけで、画像をスクリーン座標系に変換する際にもういっちょ。

code:rs

pub fn into_screen(&self, size:&Point2) -> Vec4Screen {

let scale = size.y as f64 / 2.;

let v = self.0.scaled_xy(&scale, &(-scale)) + size.to_vec4() / 2.;

Vec4Screen(v)

}

https://gyazo.com/f25431194d5f7596d2507b69fe37ad58

揃った！

https://gyazo.com/919410ffe20f074cf951c8274dbd78ac