p5.jsのキャンバスで簡単な3D描画

はじめに

　ちょっとした3D描画をしてみようね。楽しそうです。
　カメラの設定：視点[40,60,80], 注視点[140,160,180], 上方向：[0,0,1],
　視野角(fov):PI/3, アスペクト比：1,
　ニアクリップ：10, ファークリップ：1000.
　つまりこんな感じ：

　描画する立体は中心が(100,150,200)にある一辺の長さが40の立方体です。辺は軸に平行で、回転もしていません。コードとしては、こう：

function setup() {
	createCanvas(400,400,WEBGL);
	background(0);
	noFill();
	stroke(255);
	camera(40,60,80,140,160,180,0,0,1);
	perspective(PI/3, 1, 10, 1000);
	translate(100,150,200);
	box(40);
}

そうするとこんな感じの結果になるわけです。これをやってみましょう。

立方体の詳細

　こんな風です。(100,150,200)の上下左右前後に±20ずつ移動させて8つの頂点を用意し、それらを立方体になるようにつなぎます。つまり、こう：

このやり方で辺を描画しても同じ結果になります。やってみるか。

function setup() {
	createCanvas(400,400,WEBGL);
	background(0);
	noFill();
	stroke(255);
	camera(40,60,80,140,160,180,0,0,1);
	perspective(PI/3, 1, 10, 1000);

	const points = [];
	for(let i=-1;i<2;i+=2){for(let j=-1;j<2;j+=2){for(let k=-1;k<2;k+=2){
		points.push(createVector(100+20*i,150+20*j,200+20*k));
	}}}
	const indices = [
		0,1,1,3,3,2,2,0, 4,5,5,7,7,6,6,4, 0,4,1,5,2,6,3,7
	];

	for(let i=0; i<indices.length; i+=2){
		const p = points[indices[i]];
		const q = points[indices[i+1]];
		line(p.x,p.y,p.z, q.x,q.y,q.z);
	}
}

さっきの立方体の辺部分だけ描画するコードです。結果は、こう：

一緒です。まあ当然か。ただ頂点の位置座標が分かんないとこれ以降の議論ができないので必要な過程です。やりたいことは、この画像の算出です。計算で辺の位置を出そうというわけです。そこまでやったら終わりです。

ビュー変換（カメラから見る）

　初めに、3D空間の位置をカメラから見た場合の座標系に移します。つまり、カメラの位置が原点で、直交座標系をカメラ中心に設定した場合の座標の値を出そうというわけ。で、軸の方向はこんな感じで計算しています。

	const v0 = createVector(1,1,1);
	const up = createVector(0,0,1);
	const viewX = v0.cross(up).normalize();
	const viewZ = v0.copy().mult(-1).normalize();
	const viewY = viewZ.cross(viewX).normalize();

v0はカメラから注視点に向かうベクトルです。upは「暫定上方向」で、これを基準に「右」が決まります。カメラの「右」とは、注視点方向、「暫定上方向」の順でベクトルの外積を取って取得します。これがviewX. 次に「前方向」は注視点方向の逆です。カメラから目に向かう方向ですね。これがviewZで、単純に逆にして正規化します。最後に「上方向」（直交するように取る真の上方向）はviewZ cross viewXで用意します。viewYです。これでできました。
　その座標系に移すのは簡単です。視点ベクトルを引いて原点が視点となるようにし、あとは内積を取るだけ。

const applyView = (v) => {
    v.sub(40, 60, 80);
    const x = v.dot(viewX);
    const y = v.dot(viewY);
    const z = v.dot(viewZ);
    v.set(x,y,z);
}

これでビュー座標系になりました。次に、射影です。目標は2次元に落とすことです。なお、透視投影です。

透視投影変換（１）

　射影変換とはすべての座標を-1～1に落とすことです。ただちょっと工夫します。

透視投影では角錐を用意します。視点からの距離、この場合Lであらわされている-zですが、この位置により見える範囲が決まっています。なのでその位置に応じた上下左右の範囲で-1～1となるようにします。たとえばyの場合、視野角の半分である$\theta/2$を使って$y/(L\tan(\theta/2))$を作れば-1～1に落ちるわけですね。

アスペクト比は横幅を縦幅で割ったものです。これを掛けることで、 xも-1～1に落ちます。ここで、これらはLで割って出しているんですが、この処理は非線形なので、行列演算で書けません。そこで、Lの算出は-1倍するだけ、線形ですから、4番目の座標をLにすることで、線形処理のみとする工夫が要求されます。図のように、 xとyは一旦$\tan(\theta/2)$で割ったりアスペクト比を掛けるだけとし、最後に4番目の座標で割ればよいということになります。4番目の座標の算出も線形処理です。あとは3番目の座標を線形処理で出して、これをLで割った時に-1～1となればいいですね。-1になるのは、 zがニアクリップのときです。1になるのは、 zがファークリップのときです。

透視投影変換（２）

ベクトルをx,y,z,1としています。1を追加しました。4つ目の座標はL=-zで確定していますが、3番目をどうするかが未定です。しかし両端が分かっていることと線形であることから算出できます。まず zが-nの場合（注意してください、z軸は逆方向を向いてるので、ニアクリップのとき zは-nで、ファークリップのとき zは-fです）、Lはnですから、これで割って-1とするには変換の結果は-nである必要があります。逆にファークリップの時は-fを fにしなければなりません。つまり、-n,-fを-n, fにする必要があるのです。4番目はそれぞれn,fです。つまり、

M \cdot \begin{pmatrix} -n & -f \\ 1 & 1 \end{pmatrix} = \begin{pmatrix} -n & f \\ n & f \end{pmatrix}

というわけですね。この$M$は2x2行列です。図で言うところの右下です。ここまで落とせばあとは逆行列の計算で$M$は容易に算出できます。できました。

2x2行列の逆行列が分かんない人は高校数学の...あ、やってないのか...？まあいいや、それで、できました。これでうまくいくはずです。

const applyProj = (aspect, fov, n, f, v) => {
    const factor = 1/tan(fov/2);
    const px = v.x * aspect * factor;
    const py = v.y * factor;
    const pz = ((n+f)/(n-f)) * v.z + (2*n*f)/(n-f);
    const pw = -v.z;
    v.set(px/pw, py/pw, pz/pw);
}

最終的に4番目の成分(=L)で割っています。これですべて-1～1に落ちたはずです。描画して確かめましょう。とはいえ今回必要なのは xとyだけですが。

コード全文

3D

function setup() {
	createCanvas(400, 400);

	const v0 = createVector(1,1,1);
	const up = createVector(0,0,1);
	const viewX = v0.cross(up).normalize();
	const viewZ = v0.copy().mult(-1).normalize();
	const viewY = viewZ.cross(viewX).normalize();
	const applyView = (v) => {
		v.sub(40, 60, 80);
		const x = v.dot(viewX);
		const y = v.dot(viewY);
		const z = v.dot(viewZ);
		v.set(x,y,z);
	}
	// 射影
	// ポイントはまず-zが正の数であり、最終的にこれで割るというわけ。その際に
	// xが-1～1となるように、yが-1～1となるようにする。
	// -zが正の数であるならどの範囲にxがあるときz～-zとなるかは簡単にわかる
	// どの範囲にyがあるときz～-zとなるかも簡単にわかる
	// zについてだが-n～-fの範囲にあるとき-zで割って-1～1となるようにする必要があるため
	// 2x2行列を書けた結果が
	// -n,1のとき-n,nとなるよう、-f,1のときf,f... ここが難しい。
	// まず大前提としてzは負の数である
	// それゆえ限界値のときzサイドは-n～-fという値になる。これはいいですね？
	// それで、第4成分は1として行列を掛ける
	// 行列を掛けるのは結局線形変換の逆変換は行列を使うのが一番算出しやすいから。
	// その結果において第4成分が-zになることは決まっていてそれぞれn,fとなるわけで
	// それで割って-1,1になるのだから-n,fで確定するということ
	// ゆえに[-n,1]->[-n,n]で[-f,1]->[f,f]という変換が必要とされるのですよね
	// x,yも同じように考えた方がいいのかもしれないですね
	// -zにtan(fov/2)を掛けたら1になるのがyであってね。だから
	// 1/((-z)*tan(fov/2))を掛けるんですがそれを-zで割って得るんで
	// ひとまず1/tan(fov/2)を掛けようってわけさ。
	// xはそれにaspectを掛けてるだけさ。おわり。
	const applyProj = (aspect, fov, n, f, v) => {
		const factor = 1/tan(fov/2);
		const px = v.x * aspect * factor;
		const py = v.y * factor;
		const pz = ((n+f)/(n-f)) * v.z + (2*n*f)/(n-f);
		const pw = -v.z;
		v.set(px/pw, py/pw, pz/pw);
	}

	const points = [];
	for(let i=-1;i<2;i+=2){for(let j=-1;j<2;j+=2){for(let k=-1;k<2;k+=2){
		points.push(createVector(100+20*i,150+20*j,200+20*k));
	}}}
	const indices = [
		0,1,1,3,3,2,2,0, 4,5,5,7,7,6,6,4, 0,4,1,5,2,6,3,7
	];

	for(const p of points){
		applyView(p);
		applyProj(1, PI/3, 10, 1000, p);
		p.x += 1;
		p.y += 1;
		p.x *= 0.5;
		p.y *= 0.5;
		p.x *= 400;
		p.y *= 400;
	}

	background(0);
	stroke(255);
	strokeWeight(2);

	for(let i=0; i<indices.length; i+=2){
		const p = points[indices[i]];
		const q = points[indices[i+1]];
		line(p.x,p.y, q.x,q.y);
	}
}

　点と辺は改めて3次元ベクトルで用意しています。そのあとでapplyViewとapplyProjを適用しています。最後に、1を足して2で割って0～1とし、400を掛けてキャンバスに合わせています。ほんとはyだけ正規化デバイス座標系との整合性を取るために逆を取らないといけないんですが、p5はいわゆる「y反転」を採用しているのでそのままでいいです。結果はこちら。

サイトを閲覧してコピペではないことを確かめてください。ほんとですよ。まあ絵を作るだけなら zの算出は要らないんですよね...そのうちそういうコードも書けるといいですね。
　以上です。

おわりに

　楽しかったですね。
　ここまでお読みいただき感謝です。