4. 内積

此処まででベクトルの足し算、拡大縮小、長さの測定、正規化が出来る様になりました。
然し、とても自然な疑問にはまだ答えられません。
二つのベクトルはどれ位同じ方向を向いているのか?

内積は、正に其の疑問に一つの数で答えます。
恐らくグラフィックスプログラミング全体で最もよく使われる演算です。
3Dゲームで見た事のあるライティングされたピクセルは全て、少なくとも一回の内積で計算されています。其れに加えて、角度、射影、反射、「其れは前にあるか?」の判定、壁ずり移動も内積で出来ます。
此のシリーズで一つの章だけ覚えるなら、此の章にして下さい。

4.1 考え方:二つの矢印はどれ位一致しているか?

ロープでソリを引いている所を想像して下さい。
真っ直ぐ前に引けば、力は全てソリを動かすのに使われます。
ロープが斜め上を向いていると、力の一部は前に動かす代わりにソリを持ち上げる為に無駄になります。
真上に引くと、ソリは全く前に進みません。

ソリを真っ直ぐ前、斜め上、真上に引いた時に前に進む量

内積が測るのは、此の「片方の矢印のどれ位がもう片方の矢印に沿っているか」です。
矢印が同じ方向を向いていれば大きな正の数、直角なら0、反対方向なら負の数になります。

4.2 公式

二つのベクトル𝐚と𝐛の内積は𝐚⋅𝐛と書きます(「a ドット b」と読みます。英語では「dot product」と言います)。
結果はベクトルではなく、一つの数(スカラー)です。
公式は二つあり、此の章の魔法は、其の二つが常に同じ答えになる事です。

代数的な公式:対応する成分同士を掛けて、全部足します。

𝐚⋅𝐛=axbx+ayby+azbz

2次元ではzの項を省くだけです。一般には、第1章のシグマ記号で:

𝐚⋅𝐛=∑i=1naibi=a1b1+a2b2+⋯+anbn

幾何的な公式:二つの長さと、矢印の間の角度θ(シータ)のコサインを掛けます。

𝐚⋅𝐛=‖𝐚‖‖𝐛‖cosθ

一つ目の公式はコンピューターが計算する方法です。
二つ目の公式は其の意味です。

4.3 仕組み

4.3.1 計算してみる

𝐚=(1,2,3)、𝐛=(4,−5,6)とします。

𝐚⋅𝐛=1×4+2×(−5)+3×6=4−10+18=12

掛け算三回と足し算二回です。
とても安い計算で、其れがGPUが毎秒何十億回も出来る理由の一つです。

4.3.2 性質

代数的な公式から直接、幾つかの規則が出てきます。各成分の普通の掛け算と足し算で成り立つ規則だからです:

  • 交換法則:𝐚⋅𝐛=𝐛⋅𝐚。axbx=bxaxなので、順番は関係ありません。
  • 分配法則:𝐚⋅(𝐛+𝐜)=𝐚⋅𝐛+𝐚⋅𝐜。
  • スカラーは外に出せる:(s𝐚)⋅𝐛=s(𝐚⋅𝐛)。
  • 自分自身との内積は長さの二乗:𝐚⋅𝐚=ax2+ay2+az2=‖𝐚‖2。

最後の規則が此の章と第3章を繋ぎます。
𝐚=(1,2,3)の場合:𝐚⋅𝐚=1+4+9=14で、‖𝐚‖=14です。
つまり、第3章の長さの二乗は単にdot(v, v)です。

4.3.3 何故二つの公式は一致するのか

此れは此の章で最も重要な「何故」なので、ゆっくり進みましょう。
三つのステップで説明します。

ステップ1:内積は長さだけで決まる。
‖𝐚−𝐛‖2を見てみます。
最後の性質から、此れは(𝐚−𝐛)⋅(𝐚−𝐛)に等しく、分配法則を使って(x−y)2=x2−2xy+y2と同じ様に展開出来ます:

‖𝐚−𝐛‖2=𝐚⋅𝐚−𝐚⋅𝐛−𝐛⋅𝐚+𝐛⋅𝐛=‖𝐚‖2−2𝐚⋅𝐛+‖𝐛‖2

並べ替えると:

𝐚⋅𝐛=‖𝐚‖2+‖𝐛‖2−‖𝐚−𝐛‖22

先程の数字で確認しましょう。
‖𝐚‖2=14、‖𝐛‖2=16+25+36=77、𝐚−𝐛=(−3,7,−3)なので‖𝐚−𝐛‖2=9+49+9=67です。
すると14+77−672=242=12。先程と同じです。

何故此れが役に立つのでしょうか?
𝐚と𝐛を原点から始まる矢印として描きます。
二つの先端を結ぶ矢印𝐚−𝐛と合わせると三角形が出来て、上の公式は其の三辺の長さしか使っていません。

原点からのベクトルaとb、其の間の角度θ、先端を結ぶa − bが作る三角形

ステップ2:図を回しても内積は変わらない。
三角形全体を回転させても、辺の長さは変わりません。
従って、内積も変わりません。
と言う事は、図を最も都合の良い位置に回しても良いのです。𝐚が正のx軸上に来る位置にします。

ステップ3:都合の良い位置で計算する。
此れで𝐚=(‖𝐚‖,0)です。
そして𝐛はx軸と角度θを成すので、第1章の単位円(第3章の通り𝐛の長さ倍したもの)から𝐛=(‖𝐛‖cosθ,‖𝐛‖sinθ)です。
代数的な公式を使います:

𝐚⋅𝐛=‖𝐚‖×‖𝐛‖cosθ+0×‖𝐛‖sinθ=‖𝐚‖‖𝐛‖cosθ

完成です。
3次元でも同じです。原点から始まる二つの矢印は常に一つの平面上にあり、其の平面をxy平面になるまで回せます。

おまけ:幾何的な公式をステップ1に代入すると、‖𝐚−𝐛‖2=‖𝐚‖2+‖𝐛‖2−2‖𝐚‖‖𝐛‖cosθとなり、此れは三角法の「余弦定理」です。θ=90°ではコサインが0になり、ピタゴラスの定理になります。

4.3.4 符号が向きを教えてくれる

長さは決して負にならないので、𝐚⋅𝐛の符号はcosθの符号と同じです。
単位円から、其の意味が分かります:

内積 角度 意味
>0 90°未満 大体同じ方向
=0 丁度90° 垂直
<0 90°より大きい 大体反対方向

鋭角・直角・鈍角での内積の符号

其の為、内積は垂直なベクトルの標準的な判定方法です。0なら直角です。
例:(1,2,2)⋅(2,1,−2)=2+2−4=0なので、数字を見ただけでは分かり難いですが、此の二つのベクトルは垂直です。

4.3.5 角度を求める

二つの公式が同じ数になると分かっているので、角度について解く事が出来ます:

cosθ=𝐚⋅𝐛‖𝐚‖‖𝐛‖,θ=arccos(𝐚⋅𝐛‖𝐚‖‖𝐛‖)

此処でarccos(アークコサイン、コードではacos)はコサインの逆関数です。−1から1の間のコサインの値を与えると、0からπの間の角度を返します。

計算例:𝐚=(1,1,0)、𝐛=(0,1,1)。

𝐚⋅𝐛=1×0+1×1+0×1=1‖𝐚‖=2,‖𝐛‖=2cosθ=122=12θ=arccos(12)=π3≈1.0472 rad=60°

2次元の別の例:(2,0)と(1,1)ではcosθ=222=12≈0.7071で、45°です。

本当に便利な特別な場合があります。両方が単位ベクトルなら、内積はコサイン其の物です。両方の長さが1だからです:

𝐚^⋅𝐛^=cosθ

其の為、グラフィックスのコードは内積を取る前に全てを正規化します。そうすれば掛け算三回の値段でコサインが手に入ります。

クランプの落とし穴。
浮動小数点数は正確ではありません。
同じ方向を向いた二つの単位ベクトルの内積が、1ではなく1.0000001になる事があります。
acosは−1から1の値にしか定義されていないので、acos(1.0000001)はNaNを返します。第3章で、NaNがどれだけ有害か既に知っていますね。
acosを呼ぶ前に、必ず値を[−1, 1]にクランプして下さい。

実際には、角度其の物が必要無い事もよくあります。
「角度は30°より小さいか?」を知りたいなら、代わりに単位ベクトルの内積をcos30°≈0.866と比較します。
コサインが大きい程角度は小さいので、高価なacosを節約出来ます。

4.3.6 射影:あるベクトルの、別のベクトル上の影

ソリの話に戻りましょう。「引く力のどれ位が前に向かうか」は、引く力の前方向への射影と呼ばれます。
𝐛に向かって真上から光を当てる所を想像して下さい。𝐛の直線上に落ちる𝐚の影が射影です。

a = 4, 2のb = 3, 4への射影2.4, 3.2と、垂直な部分1.6, −1.2

スカラー射影は、其の影の符号付きの長さです。
𝐚の長さ×cosθで、此れは𝐛の単位ベクトルとの内積です:

影の長さ=‖𝐚‖cosθ=𝐚⋅𝐛^

ベクトル射影は、影を矢印にした物です。其の長さ×𝐛の向き:

proj𝐛𝐚=(𝐚⋅𝐛^)𝐛^=𝐚⋅𝐛𝐛⋅𝐛𝐛

二つ目の形は、𝐛⋅𝐛=‖𝐛‖2なので平方根を使わずに済みます。

計算例:𝐚=(4,2)を𝐛=(3,4)に射影します。

𝐚⋅𝐛=4×3+2×4=12+8=20‖𝐛‖=5,𝐛^=(0.6,0.8)𝐚⋅𝐛^=205=4proj𝐛𝐚=4×(0.6,0.8)=(2.4,3.2)

残りの部分𝐚−proj𝐛𝐚=(4−2.4,2−3.2)=(1.6,−1.2)は、𝐛に垂直な筈です。
確認:(1.6,−1.2)⋅(3,4)=4.8−4.8=0。正しいです。

つまり、全てのベクトルは𝐛に沿った部分と、其れに垂直な部分に分けられます。
此の分解が、次の二つの技の基礎です。

4.3.7 成分を取り除く:壁ずり移動

キャラクターが斜めに壁に向かって歩いた時、ぴたっと止まって欲しくはありません。
壁に沿って滑って欲しいのです。
言い換えると、移動のうち壁に沿った部分は残し、壁に向かう部分は捨てます。

壁の「壁から出る」方向は、其の単位法線𝐧^です。
速度𝐯のうち𝐧^に沿った部分は射影(𝐯⋅𝐧^)𝐧^なので、其れを引きます:

𝐯slide=𝐯−(𝐯⋅𝐧^)𝐧^

上から見た壁ずり移動:速度v、壁の法線、残った滑る向き

計算例:プレイヤーが𝐯=(4,0,−3)(右前方)で動いていて、右側に法線が左を向いた壁𝐧^=(−1,0,0)があります。

𝐯⋅𝐧^=4×(−1)+0+(−3)×0=−4𝐯slide=(4,0,−3)−(−4)(−1,0,0)=(4,0,−3)−(4,0,0)=(0,0,−3)

横向きの部分は消えて、プレイヤーは壁に沿って前に歩き続けます。
内積が負である事に注目して下さい。プレイヤーが壁の中へ向かって動いているという意味です。
𝐯⋅𝐧^≥0なら、プレイヤーは壁から離れているので、速度は其のままにしておくべきです。

4.3.8 反射

跳ねるボール、鏡に当たるレーザー、水面に映る空。表面の単位法線𝐧^に対して、向き𝐝を鏡の様に反転させる必要があります。
壁ずりでは法線方向の部分を一回取り除きました。
反射では二回取り除きます。一回目で表面まで、もう一回で反対側へ出てきます。

𝐫=𝐝−2(𝐝⋅𝐧^)𝐧^

向きd = 1, −1が法線0, 1の床で反射してr = 1, 1になる様子

計算例1:𝐝=(1,−1)で動くボールが床𝐧^=(0,1)に当たります。
𝐝⋅𝐧^=−1なので、𝐫=(1,−1)−2(−1)(0,1)=(1,−1)+(0,2)=(1,1)。
横の動きは残り、縦の動きが反転します。正に跳ね返りの見た目です。

計算例2:真下に向かう光線𝐝=(0,−1,0)が、傾いた表面𝐧^=(0.6,0.8,0)に当たります(0.36+0.64=1なので単位ベクトルです)。

𝐝⋅𝐧^=0×0.6+(−1)×0.8+0×0=−0.8𝐫=(0,−1,0)−2(−0.8)(0.6,0.8,0)=(0,−1,0)+(0.96,1.28,0)=(0.96,0.28,0)

確認:‖𝐫‖2=0.9216+0.0784=1。
鏡なので当然ですが、反射は長さを変えません。
此の公式は𝐧^が単位ベクトルの時だけ正しく動きます。長い法線だと引き過ぎてしまいます。

GLSL、Slang、HLSLには、全く此の公式のreflect(d, n)が組み込まれています。
法線は正規化されている必要があり、dは表面に向かう向き(入射方向)です。

4.4 グラフィックスプログラミングでの使い方

4.4.1 ランバート拡散反射

つや消しの表面(紙、粘土、壁)は、光が正面から当たる時に最も明るく、光が浅い角度で入る程暗くなります。
ヨハン・ハインリッヒ・ランベルトが1760年に此れを記述しました。明るさはcosθに比例し、θは表面の法線と光源への向きの間の角度です。

光が正面から当たる場合と浅い角度で当たる場合。浅い角度では同じ光がより広い面積に広がる

何故コサインなのでしょうか?
一定の幅の光線が浅い角度で表面に当たると、より広い面積に広がるので、1平方センチメートル当たりの光が少なくなります。
そして両方のベクトルが単位ベクトルなら、コサインは単なる内積です。
光が表面の裏側にあると内積は負になります。負の光は意味が無いので、0にクランプします:

diffuse=max(𝐧^⋅𝐥^,0)

数字で:床𝐧^=(0,1,0)に、𝐥^=(0.6,0.8,0)(光源への向き)から光が来ます。
𝐧^⋅𝐥^=0+0.8+0=0.8なので、床は光の明るさの80%を受け取ります。
光が床の下にあれば𝐥^=(0.6,−0.8,0)で、内積は−0.8になり、maxが其れを0に変えます。全く照らされません。

小さな、完全なフラグメントシェーダーです:

struct Scene {
  float3 lightDir;   // unit vector from the surface towards the light
  float3 lightColor;
  float3 albedo;     // surface color
};
ConstantBuffer<Scene> scene;

[shader("fragment")]
float4 fragmentMain(float3 normal : NORMAL) : SV_Target {
  float3 n = normalize(normal); // interpolation made it shorter (chapter 3)
  float diffuse = max(dot(n, scene.lightDir), 0.0);
  float3 color = scene.albedo * scene.lightColor * diffuse;
  return float4(color, 1.0);
}
#version 450

layout(location = 0) in vec3 normal;
layout(location = 0) out vec4 outColor;

layout(set = 0, binding = 0) uniform Scene {
  vec3 lightDir;   // unit vector from the surface towards the light
  vec3 lightColor;
  vec3 albedo;     // surface color
};

void main() {
  vec3 n = normalize(normal); // interpolation made it shorter (chapter 3)
  float diffuse = max(dot(n, lightDir), 0.0);
  vec3 color = albedo * lightColor * diffuse;
  outColor = vec4(color, 1.0);
}
struct Scene {
  float3 lightDir;   // unit vector from the surface towards the light
  float3 lightColor;
  float3 albedo;     // surface color
};
ConstantBuffer<Scene> scene;

float4 fragmentMain(float3 normal : NORMAL) : SV_Target {
  float3 n = normalize(normal); // interpolation made it shorter (chapter 3)
  float diffuse = max(dot(n, scene.lightDir), 0.0);
  float3 color = scene.albedo * scene.lightColor * diffuse;
  return float4(color, 1.0);
}

注意点が二つあります。
一つ目に、nとlightDirは同じ座標空間(普通はワールド空間かビュー空間。空間については第11章で説明します)になければなりません。
二つ目に、光の向きは表面から光源へ向かいます。エンジンが「光が進む向き」として保存しているなら先に符号を反転させて下さい。そうしないと、照らされる面と暗い面が全部逆になります。
実際のレンダラーでは、照らされない側が真っ黒にならない様に、小さな一定の「環境光(アンビエント)」の項も加えます。

4.4.2 ブリン・フォンの鏡面ハイライト

つやのある表面には、光が目の方向へ反射される所に明るい点(ハイライト)が見えます。
ブリン・フォンのモデルは、此れをハーフベクトル𝐡^で計算します。光源への向き𝐥^と視点への向き𝐯^の丁度中間の向きです:

𝐡^=𝐥^+𝐯^‖𝐥^+𝐯^‖,specular=max(𝐧^⋅𝐡^,0)s

表面の法線が丁度ハーフベクトルを向いていれば、表面は光と目の間の完全な鏡の様に振る舞い、ハイライトは最大の強さになります。
指数s(「シャイニネス」)はハイライトがどれ位早く弱まるかを決めます。大きなsは小さく鋭いハイライト(磨いた金属)、小さなsは広いハイライト(プラスチック)になります。

𝐧^=(0,1,0)、𝐥^=(0.6,0.8,0)での計算例:

  • 視点が鏡の位置、𝐯^=(−0.6,0.8,0):𝐥^+𝐯^=(0,1.6,0)なので𝐡^=(0,1,0)、𝐧^⋅𝐡^=1、specular = 1。
  • 視点が真上、𝐯^=(0,1,0):𝐥^+𝐯^=(0.6,1.8,0)で、長さは0.36+3.24=3.6≈1.8974。
  • 従って𝐡^≈(0.3162,0.9487,0)で、𝐧^⋅𝐡^=1.83.6≈0.9487。其の二乗は丁度3.243.6=0.9です。
  • s=8なら0.94878=0.94=0.6561。s=32なら0.916≈0.1853。

同じ角度でも、よりつやのある表面ではハイライトの殆どが既に消えています。

float blinnPhong(float3 n, float3 l, float3 v, float shininess) {
  float3 h = normalize(l + v);
  return pow(max(dot(n, h), 0.0), shininess);
}
float blinnPhong(vec3 n, vec3 l, vec3 v, float shininess) {
  vec3 h = normalize(l + v);
  return pow(max(dot(n, h), 0.0), shininess);
}
float blinnPhong(float3 n, float3 l, float3 v, float shininess) {
  float3 h = normalize(l + v);
  return pow(max(dot(n, h), 0.0), shininess);
}

4.4.3 リムライト

リムライトは物体の輪郭を光らせる効果で、キャラクター、幽霊、バリア等によく使われます。
輪郭とは表面がカメラから背を向ける所、つまり法線が視線方向に垂直な所なので、𝐧^⋅𝐯^は0に近くなります。
物体の真ん中、こちらを向いた所では𝐧^⋅𝐯^は1に近いです。
従って1−𝐧^⋅𝐯^は真ん中で0、端で1になり、正にリムのマスクです:

float rim(float3 n, float3 v, float power) {
  return pow(1.0 - saturate(dot(n, v)), power);
}
float rim(vec3 n, vec3 v, float power) {
  return pow(1.0 - clamp(dot(n, v), 0.0, 1.0), power);
}
float rim(float3 n, float3 v, float power) {
  return pow(1.0 - saturate(dot(n, v)), power);
}

例えば𝐯^=(0,1,0)の場合、カメラを向いた法線(0,1,0)は1−1=0、法線(0.6,0.8,0)は1−0.8=0.2、シルエット上の法線(1,0,0)は1−0=1になります。
powerを大きくするとリムが細く、小さくすると太くなります。

4.4.4 向きの判定:前か後ろか?

内積の符号は「同じ側か反対側か」を教えてくれるので、何かが何かの前にあるかを調べる標準的な方法です。

敵は自分の前にいるか?
自分の前方向𝐟^と、自分から敵への矢印𝐞−𝐩を取ります。
内積が正なら敵は前に、負なら後ろにいます。
長さは符号を変えないので、此の判定には正規化すら必要ありません。

例:自分は原点に立って−Zを向いています、𝐟^=(0,0,−1)。
(3,0,−4)にいる敵:(3,0,−4)⋅(0,0,−1)=4>0、前にいます。
(1,0,2)にいる敵:(1,0,2)⋅(0,0,−1)=−2<0、後ろにいます。

敵は視野の円錐の中にいるか?
ステルスゲームの見張りは左右に45°ずつ見えるとします。
今度はコサインが欲しいので正規化します。
自分から敵への矢印を正規化した物を𝐝^とします。
𝐟^⋅𝐝^>cos45°≈0.7071なら、敵は円錐の中にいます。

(3,0,−4)にいる敵の場合:長さは5なので𝐝^=(0.6,0,−0.8)、𝐟^⋅𝐝^=0.8です。
0.8>0.7071なので、見張りは敵を見つけます(実際の角度はarccos(0.8)≈36.87°ですが、計算する必要はありませんでした)。
cos45°は毎フレームではなく、ゲーム開始時に一回だけ計算しておきます。

上から見た左右45°ずつの見張りの視野と、其の中の3, 0, −4にいる敵

此の三角形はカメラを向いているか?
法線𝐧を持つ表面上の点𝐩は、𝐧⋅(𝐜−𝐩)>0なら𝐜にあるカメラの方を向いています。
閉じたモデルの裏側は決して見えないので、描画を省略出来ます(「背面カリング」)。
実際のGPUは代わりに画面上での三角形の頂点の巡回順序で判断します。其れは第5章で見ますが、考え方は同じです。

4.4.5 シェーダーでの反射と壁ずり

組み込みのreflectは鏡の様な効果に使われます。カメラから表面への向きを計算し、法線で反射させ、其の結果で空のテクスチャ(キューブマップ)を参照します。

float3 mirrorDirection(float3 cameraPos, float3 worldPos, float3 normal) {
  float3 d = normalize(worldPos - cameraPos); // incoming: from the eye to the surface
  return reflect(d, normalize(normal));
}
vec3 mirrorDirection(vec3 cameraPos, vec3 worldPos, vec3 normal) {
  vec3 d = normalize(worldPos - cameraPos); // incoming: from the eye to the surface
  return reflect(d, normalize(normal));
}
float3 mirrorDirection(float3 cameraPos, float3 worldPos, float3 normal) {
  float3 d = normalize(worldPos - cameraPos); // incoming: from the eye to the surface
  return reflect(d, normalize(normal));
}

壁ずり移動は普通、CPU上の物理やキャラクターコントローラーのコードにありますが、上と同じ射影の公式です。下のC++のテストプログラムで見られます。

4.5 其の他の使い道

内積は、科学と工学全体で最もよく使われる演算の一つです。

  • 物理の仕事:力𝐅が物体を𝐝だけ動かす時の仕事はW=𝐅⋅𝐝です。ソリを𝐅=(30,40)ニュートンで𝐝=(5,0)メートル引くと、W=150+0=150ジュールです。幾何的な公式でも:‖𝐅‖=50、cosθ=0.6で、50×5×0.6=150。同じ答えです。
  • コサイン類似度:検索エンジン、レコメンドシステム、AIの「埋め込み(エンベディング)」は、文書、商品、文を長いベクトルに変換し、cosθ=𝐚^⋅𝐛^で比較します。三つの文書で「猫」「犬」「車」の単語を数えると𝐚=(2,1,0)、𝐛=(4,2,0)、𝐜=(0,1,3)。𝐚と𝐛の類似度は丁度1(同じ話題で、片方が長いだけ)、𝐚と𝐜は1510≈0.1414(似ていない)です。
  • 統計:二つのデータ系列の相関係数は、各系列から其の平均を引いた後のコサイン類似度です。
  • 重み付き和:買い物のレシートは内積です。値段(120,250,80)円の商品を(3,2,5)個ずつ買うと、360+500+400=1260円です。重み(0.2,0.3,0.5)と点数(80,90,70)の最終成績は16+27+35=78です。
  • ニューラルネットワーク:人工ニューロン一つは、重みと入力の内積(に定数を足した物)を計算します。行列版は第7章で見ます。
  • オーディオ:信号を或る周波数のサイン波と比較するのは内積です。其れがフーリエ変換の核となる考え方で、音にどの周波数(音程)が含まれているかを教えてくれます。

4.6 コード

内積を、三つ全てのベクトルの大きさについてlinalg.hhに追加します。
第2章のスカラー倍の直後に置いて下さい:

// 内積 / dot product (chapter 4)
inline float dot(Vec2 a, Vec2 b) { return a.x * b.x + a.y * b.y; }
inline float dot(Vec3 a, Vec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
inline float dot(Vec4 a, Vec4 b) { return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; }

Vec4版は第6章で重要になります。其処では行列の各行とベクトルの内積を取ります。

角度、射影、反射の補助関数は短いので、ヘッダーではなくテストプログラムの中で定義します。
reflectはシェーダーの組み込み関数と同じく、単位法線を受け取ります。
degrees()は第1章の補助関数です。

#include "linalg.hh"
#include <cstdio>

// Angle between two vectors, in radians.
float angleBetween(Vec3 a, Vec3 b) {
  float c = dot(normalize(a), normalize(b));
  if (c > 1.0f) c = 1.0f;   // rounding can push c slightly
  if (c < -1.0f) c = -1.0f; // outside [-1, 1]; acos would return NaN
  return std::acos(c);
}

// Vector projection of a onto b (b must not be zero).
Vec3 project(Vec3 a, Vec3 b) {
  return (dot(a, b) / dot(b, b)) * b;
}

// Reflect direction d off a surface with unit normal n.
Vec3 reflect(Vec3 d, Vec3 n) {
  return d - 2.0f * dot(d, n) * n;
}

int main() {
  Vec3 a{1.0f, 2.0f, 3.0f};
  Vec3 b{4.0f, -5.0f, 6.0f};
  printf("dot(a, b)       = %g\n", dot(a, b));
  printf("dot(a, a)       = %g\n", dot(a, a));

  float angle = angleBetween({1.0f, 1.0f, 0.0f}, {0.0f, 1.0f, 1.0f});
  printf("angle           = %.4f rad = %.2f deg\n", angle, degrees(angle));

  Vec3 p = project({4.0f, 2.0f, 0.0f}, {3.0f, 4.0f, 0.0f});
  printf("projection      = (%g, %g, %g)\n", p.x, p.y, p.z);

  // Wall sliding: remove the part of v that goes into the wall
  Vec3 v{4.0f, 0.0f, -3.0f};
  Vec3 wall{-1.0f, 0.0f, 0.0f};
  Vec3 slide = v - dot(v, wall) * wall;
  printf("slide           = (%g, %g, %g)\n", slide.x, slide.y, slide.z);

  Vec3 r = reflect({0.0f, -1.0f, 0.0f}, {0.6f, 0.8f, 0.0f});
  printf("reflect         = (%.2f, %.2f, %.2f)\n", r.x, r.y, r.z);

  // Lambert: brightness = max(dot(n, l), 0)
  Vec3 n{0.0f, 1.0f, 0.0f};
  Vec3 l{0.6f, 0.8f, 0.0f};
  float diffuse = dot(n, l);
  printf("diffuse         = %g\n", diffuse > 0.0f ? diffuse : 0.0f);

  // Is the enemy in front of me, and inside a 45 degree view cone?
  Vec3 forward{0.0f, 0.0f, -1.0f};
  Vec3 toEnemy = Vec3{3.0f, 0.0f, -4.0f} - Vec3{0.0f, 0.0f, 0.0f};
  float c = dot(forward, normalize(toEnemy));
  printf("cos to enemy    = %g -> %s\n", c, c > 0.7071f ? "inside cone" : "outside cone");

  // Shopping: quantities . prices
  Vec3 quantity{3.0f, 2.0f, 5.0f};
  Vec3 price{120.0f, 250.0f, 80.0f};
  printf("total           = %g yen\n", dot(quantity, price));
  return 0;
}

期待される出力:

dot(a, b)       = 12
dot(a, a)       = 14
angle           = 1.0472 rad = 60.00 deg
projection      = (2.4, 3.2, 0)
slide           = (0, 0, -3)
reflect         = (0.96, 0.28, 0.00)
diffuse         = 0.8
cos to enemy    = 0.8 -> inside cone
total           = 1260 yen

reflectの細かい点:2.0f * dot(d, n) * nは左から右へ計算されるので、先ず二つのfloatが掛けられ、残りは第2章のfloat * Vec3演算子が行います。

4.7 まとめ

  • 内積𝐚⋅𝐛は、対応する成分を掛けて足し合わせます。結果は数です。
  • 幾何的には𝐚⋅𝐛=‖𝐚‖‖𝐛‖cosθ。内積は長さだけで決まり、長さは回転しても変わらないので、二つの公式は一致します。
  • 𝐚⋅𝐚=‖𝐚‖2。
  • 正なら大体同じ方向、0なら垂直、負なら反対方向。
  • 単位ベクトル同士なら、内積は角度のコサインです。acosの前に[−1, 1]にクランプするか、もっと良いのはコサイン同士を直接比較する事です。
  • 射影:(𝐚⋅𝐛^)𝐛^。壁ずり:𝐯−(𝐯⋅𝐧^)𝐧^。反射:𝐝−2(𝐝⋅𝐧^)𝐧^。
  • グラフィックス:ランバートのmax(dot(n, l), 0)、ハーフベクトルを使うブリン・フォン、リムライト、前後の判定、視野の円錐。

4.8 課題

  1. ★ (2,−1,4)⋅(3,5,1)と(1,2)⋅(−2,1)を計算して下さい。垂直なのはどちらの組ですか?
  2. ★ (1,0)と(1,1)の間の角度、及び(1,0)と(−1,1)の間の角度を度で求めて下さい。
  3. ★ (3,1)を(0,2)に射影し、次に(1,1)に射影して下さい。どちらの場合も、残りの部分が垂直である事を確認して下さい。
  4. ★ (3,−2,1)で動くボールが床𝐧^=(0,1,0)で跳ね返ります。新しい向きは?代わりに壁ずりを使うとどうなりますか?
  5. ★★ 両方の入力を正規化してmax(dot, 0)を返すC++の関数float lambert(Vec3 n, Vec3 l)を書いて下さい。n = {0, 2, 0}、l = {3, 4, 0}なら0.8、l = {3, -4, 0}なら0と表示される筈です。
  6. ★★ bool inViewCone(Vec3 pos, Vec3 forward, Vec3 target, float halfAngle)(角度はラジアン)を書いて下さい。原点で{0, 0, -1}を向き、半角30°の見張りの場合、目標{1, 0, -3}、{2, 0, -2}、{0, 0, 5}はtrue、false、falseになる筈です。
  7. ★★★ 内積だけを使って、𝐯−(𝐯⋅𝐧^)𝐧^が常に𝐧^に垂直である理由を説明して下さい(ヒント:𝐧^との内積を取り、𝐧^⋅𝐧^=1を使います)。次に、反射𝐫の長さが𝐝と同じである事を示して下さい。


コメント(0)

まだコメントはありません。