5. 外積

第4章の内積は、二つのベクトルを受け取って一つの数を返しました。
外積は二つのベクトルを受け取って、其の両方から直角に突き出た新しいベクトルを返します。

変な物が欲しいと思うかも知れませんが、此れは3Dグラフィックスで最もよくある問題の一つ、「此の表面はどちらを向いているか?」を解決します。
三角形の法線、カメラの「右」ベクトル、頂点の巡回順序のチェックは全て外積から来ています。
外積は三角形の面積も教えてくれますし、2次元では点が線の左にあるか右にあるかも教えてくれます。其れはGPUが三角形がどのピクセルを覆うかを決める方法でもあります。

5.1 考え方:表面から立ち上がる矢印

机の上に鉛筆を2本、同じ点から別々の方向に向けて置いて下さい。
2本合わせると、机の平らな面が決まります。
其処で3本目の鉛筆を、2本が交わる点から、両方に垂直に真っ直ぐ立てて下さい。
其の3本目の鉛筆が、外積の向きです。

机の上に平らに置かれたaとbと、其処から真っ直ぐ立つa × b

疑問が二つ残ります。3本目の鉛筆は上向きか下向きか、そしてどれ位の長さか?
答えは「右手の法則」と「aとbで張られる平行四辺形の面積」です。
両方とも後で説明します。

先ず一つ重要な制限があります。外積は3次元にしか存在しません。
2次元には突き出る為の三つ目の方向が無く、4次元では二つのベクトルに垂直な方向が無数にあるので、答えが一つに決まりません。
(数学的な珍品として7次元にも似た積がありますが、グラフィックスで必要になる事はありません。)

5.2 公式

𝐚=(ax,ay,az)と𝐛=(bx,by,bz)に対して、外積𝐚×𝐛(「a クロス b」と読みます。英語では「cross product」と言います)は:

𝐚×𝐛=[aybz−azbyazbx−axbzaxby−aybx]

其の長さは:

‖𝐚×𝐛‖=‖𝐚‖‖𝐛‖sinθ

此処でθは𝐚と𝐛の間の角度(0から180°)です。
内積ではcosθでした。内積は二つのベクトルがどれ位一致しているかを、外積はどれ位食い違っているかを測ります。

5.2.1 公式の覚え方

此の公式を眺めて覚えられる人はいません。
あたしが使っている技は「xyzzy」と呼ばれる事もあります。
x成分はyとzだけを使い、「y z 引く z y」と読みます:

cx=aybz−azby(x: yz−zy)

次の成分は、全ての文字をx → y → z → xの順に一つずつずらします:

cy=azbx−axbz(y: zx−xz)

もう一度ずらすと:

cz=axby−aybx(z: xy−yx)

つまり「x = yz − zy」だけを覚えて、文字を回し続ければ良いのです。
どの行でも、計算している成分の文字は右辺に出てきません。

行列式を既に知っている人は、一番上の行に𝐢,𝐣,𝐤を並べた3×3の行列式として外積を見た事があるかも知れません。此れは覚える為の補助に過ぎず、本当の行列式ではありません(一番上の行が数ではなくベクトルだからです)。行列式は第10章でちゃんと学ぶので、此の注が何の事か分からなければ、飛ばして下さい。

5.3 仕組み

5.3.1 計算例

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

cx=aybz−azby=2×6−3×5=12−15=−3cy=azbx−axbz=3×4−1×6=12−6=6cz=axby−aybx=1×5−2×4=5−8=−3

従って𝐚×𝐛=(−3,6,−3)です。

本当に両方に垂直でしょうか?
第4章で判定方法を学びました。垂直なベクトル同士の内積は0です。

(−3,6,−3)⋅(1,2,3)=−3+12−9=0(−3,6,−3)⋅(4,5,6)=−12+30−18=0

どちらも0です。
手で外積を計算した時は、必ず此の確認をする価値があります。殆どの間違いを見つけられます。

5.3.2 何故垂直になるのか

此れは偶然ではありません。
数字の代わりに文字を使って、𝐚と𝐚×𝐛の内積を取ってみましょう:

𝐚⋅(𝐚×𝐛)=ax(aybz−azby)+ay(azbx−axbz)+az(axby−aybx)=axaybz−axazby+ayazbx−ayaxbz+azaxby−azaybx

六つの項を見て下さい。
axaybz(一つ目)とayaxbz(四つ目)は同じ積で符号が逆なので、打ち消し合います。
二つ目と五つ目、三つ目と六つ目も打ち消し合います。
どんなベクトルでも、全部足すと0になります。
𝐛でも同じ事が起こります。
公式は、正に此の打ち消し合いが起こる様に作られているのです。

5.3.3 基底ベクトル

第2章の基底ベクトル𝐢=(1,0,0)、𝐣=(0,1,0)、𝐤=(0,0,1)の外積を取ってみましょう:

𝐢×𝐣=(0×0−0×1,0×0−1×0,1×1−0×0)=(0,0,1)=𝐤

同じ様にして確認出来ます:

𝐢×𝐣=𝐤,𝐣×𝐤=𝐢,𝐤×𝐢=𝐣

逆回りにすると負になります:

𝐣×𝐢=−𝐤,𝐤×𝐣=−𝐢,𝐢×𝐤=−𝐣

i、j、kの巡回:矢印の向きの積は正、逆向きは負

そして、どんなベクトルも自分自身との外積はゼロです:𝐢×𝐢=𝟎。

5.3.4 どちらを向くか:右手の法則

二つのベクトルに垂直な方向は二つあります。其の平面から見て「上」と「下」です。
公式は常に右手の法則で決まる方を選びます:

  1. 右手の指を𝐚の方向に向けます。
  2. 指を𝐛の方へ(近い方の回り方で)曲げます。
  3. 親指が𝐚×𝐛の方向を指します。

右手座標系:右向きのx = i、上向きのy = j、画面から手前を向くz = k = i × j

𝐢×𝐣で試してみましょう。指をx方向の右に向け、y方向の上へ曲げると、親指は画面から手前の自分の方を指します。此のシリーズの右手座標系では+Zです。
此れは𝐢×𝐣=𝐤と一致します。
此れは偶然ではありません。座標系が右手系と呼ばれるのは、正に𝐢×𝐣=𝐤が右手に従うからです。

後の為の注意:公式自体はどのエンジンでも同じですが、左手座標系(DirectXの慣習、Unity、Unreal)では結果は見た目上左手に従います。エンジン間でコードを移植して法線が裏返しになったら、此れが原因である事がよくあります。

5.3.5 順番が大事

𝐚と𝐛を入れ替えると、手を逆に曲げる事になるので親指が反転します:

𝐚×𝐛=−(𝐛×𝐚)

此れを反交換と言います。
先程の数字で:𝐛×𝐚=(4,5,6)×(1,2,3)=(5×3−6×2,6×1−4×3,4×2−5×1)=(3,−6,3)で、丁度(−3,6,−3)の符号を反転した物です。

順番が関係無い普通の掛け算や内積とは大きく違います。
グラフィックスでは、外積の二つの引数を入れ替えると法線が裏返しになります。此れは最もよくあるバグの一つです。

5.3.6 長さ:平行四辺形の面積

次は長さです。
根気があれば全部展開して確認出来る恒等式(ラグランジュの恒等式と呼ばれます)があります:

‖𝐚×𝐛‖2=‖𝐚‖2‖𝐛‖2−(𝐚⋅𝐛)2

少なくとも先程の例で確認しましょう。
‖(−3,6,−3)‖2=9+36+9=54。
右辺:‖𝐚‖2=1+4+9=14、‖𝐛‖2=16+25+36=77、𝐚⋅𝐛=4+10+18=32なので、14×77−322=1078−1024=54。一致します。

次に、内積を第4章の幾何的な形𝐚⋅𝐛=‖𝐚‖‖𝐛‖cosθに置き換えます:

‖𝐚×𝐛‖2=‖𝐚‖2‖𝐛‖2−‖𝐚‖2‖𝐛‖2cos2θ=‖𝐚‖2‖𝐛‖2(1−cos2θ)

第1章の単位円からsin2θ+cos2θ=1(単位円上のピタゴラスの定理です)なので、1−cos2θ=sin2θです。
平方根を取ると(θが0から180°の間ならsinθは負になりません):

‖𝐚×𝐛‖=‖𝐚‖‖𝐛‖sinθ

此の数は何を意味するのでしょうか?
𝐚と𝐛を辺とする平行四辺形を描きます。
其の面積は底辺×高さです。
底辺は‖𝐚‖で、高さは𝐛の先端が底辺の線からどれ位上にあるかで、‖𝐛‖sinθです:

a = 3, 0, 0とb = 2, 2, 0で張られる、高さ2、面積6の平行四辺形

従って、外積の長さは𝐚と𝐛で張られる平行四辺形の面積です。
𝐚と𝐛を辺とする三角形は其の平行四辺形の半分なので、面積は12‖𝐚×𝐛‖です。

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

𝐚×𝐛=(0×0−0×2,0×2−3×0,3×2−0×2)=(0,0,6)

長さは6です。
図で確認:底辺の長さは3で、𝐛の先端はx軸から2上にあるので、平行四辺形の面積は3×2=6です。
公式で確認:‖𝐚‖=3、‖𝐛‖=22、θ=45°、sin45°=22なので、3×22×22=3×2=6。
三角形O、a、bの面積は3です。

5.3.7 平行なベクトルはゼロになる

𝐚と𝐛が同じ(又は丁度反対の)方向を向いていれば、θは0か180°で、sinθ=0なので外積はゼロベクトルです。
納得出来ますね。平行な2本の鉛筆では机の面が決まらないので、「上」がありません。

(1,2,3)×(2,4,6)=(2×6−3×4,3×2−1×6,1×4−2×2)=(0,0,0)

此れは覚えておいて下さい。カメラのコードで実際のバグの原因になります(5.4.3節)。

5.3.8 2次元では?

2次元には本当の外積はありませんが、とても便利な残り物があります。
二つの2次元ベクトルをz=0として3次元に置き、外積を取ります:

(ax,ay,0)×(bx,by,0)=(ay×0−0×by,0×bx−ax×0,axby−aybx)=(0,0,axby−aybx)

z成分だけが残ります。
其の一つの数は、よく「2次元の外積」(又はperp-dot積)と呼ばれます:

𝐚×𝐛=axby−aybx(2次元、スカラー)

其の符号は、𝐚から𝐛へどちらに回るかを表します。正なら反時計回り(左折)、負なら時計回り(右折)、ゼロなら平行です。
5.4.5節で使います。

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

5.4.1 三角形の面法線

三角形は三つの点𝐩0,𝐩1,𝐩2で与えられます。
どちらも𝐩0から始まる二つの辺は、三角形の面の上に平らに乗っています:

𝐞1=𝐩1−𝐩0,𝐞2=𝐩2−𝐩0

其の外積は面に垂直です。此れが面法線です。
正規化(第3章)すれば、ライティング(第4章)に必要な単位法線になります:

𝐧^=𝐞1×𝐞2‖𝐞1×𝐞2‖

計算例:𝐩0=(1,0,0)、𝐩1=(0,1,0)、𝐩2=(0,0,1)。
此の三角形は、箱の角を切り落とす様に、正の領域の角を斜めに横切っています。

𝐞1=(0,1,0)−(1,0,0)=(−1,1,0)𝐞2=(0,0,1)−(1,0,0)=(−1,0,1)𝐞1×𝐞2=(1×1−0×0,0×(−1)−(−1)×1,(−1)×0−1×(−1))=(1−0,0+1,0+1)=(1,1,1)

確認:(1,1,1)⋅(−1,1,0)=0、(1,1,1)⋅(−1,0,1)=0。両方の辺に垂直です。
長さは3≈1.7321なので:

𝐧^=13(1,1,1)≈(0.5774,0.5774,0.5774)

予想通り、原点から離れる方向、「箱の角」の外側を向いています。
おまけとして、三角形の面積は外積の長さの半分です:32≈0.8660。
(此の三角形の各辺の長さは2で、一辺sの正三角形の面積は34s2=34×2=32です。一致します。)

5.4.2 頂点の巡回順序と背面カリング

三角形の二つの面のうち、どちらが「表」でしょうか?
外積は頂点の順番から其れを決めます。
三角形を見た時に頂点が反時計回りに並んでいれば、(𝐩1−𝐩0)×(𝐩2−𝐩0)は自分の方を向きます(右手の法則:指を一つ目の辺に沿わせ、二つ目の辺の方へ曲げると、親指は自分の方を向きます)。

反時計回りの三角形p0, p1, p2と、画面から手前を向く法線

同じ三つの点を逆の順番𝐩0,𝐩2,𝐩1で並べると、外積が反転し(反交換!)、法線は向こう側を向きます。
先程の三角形では𝐞2×𝐞1=(−1,−1,−1)です。

此の順番を巡回順序(ワインディングオーダー)と言い、3Dモデルは普通、全ての三角形が外側から見て反時計回りになる様に書き出されます。
GPUは此れを使ってカメラから背を向けた三角形を飛ばします(背面カリング)。閉じたモデルでは其れらは常に表側の後ろに隠れているので、三角形の約半分はシェーディングする必要がありません。

  • OpenGLでは、表面はデフォルトで反時計回り(glFrontFace(GL_CCW))で、カリングはglEnable(GL_CULL_FACE)とglCullFace(GL_BACK)で有効にします。
  • Vulkanでは、パイプラインのラスタライゼーションステートでfrontFaceとcullModeを設定します。第12章のy軸を反転するVulkan用の投影行列を使えば、VK_FRONT_FACE_COUNTER_CLOCKWISEでOpenGLと同じ規約になります。

GPUは此の判定を、投影後の画面上、つまり2次元で、5.3.8節の2次元の外積を使って行います。三角形の画面上の面積の符号が巡回順序を教えてくれます。
モデルが急に「裏返し」に見えたら(外壁ではなく内壁が見える)、何処かで巡回順序が逆になっています。エクスポーター、鏡映した変換(鏡映で巡回順序が反転する理由は第10章で説明します)、又は入れ替えた外積です。

5.4.3 基底を作る:前方向から右と上を求める

カメラ(又はキャラクター、ミサイル)には、互いに垂直な三つの単位ベクトルが必要です。前、右、上です。
普通分かっているのは、何処を見るべきか、つまり前方向𝐟^と、ワールドで大体どちらが「上」か、𝐮world=(0,1,0)だけです。
残りは外積で作れます:

𝐫^=𝐟^×𝐮world‖𝐟^×𝐮world‖,𝐮^=𝐫^×𝐟^

一つ目の式は、前方向とワールドの上の両方に垂直な「右」ベクトルを作ります。
二つ目の式は、前と右に垂直なカメラ自身の「上」を作ります(カメラが少し下を向いていれば、其の上はワールドの上と違って前に傾きます)。
𝐫^と𝐟^は互いに垂直な単位ベクトルなのでsin90°=1で、其の外積は既に長さ1です。従って二回目の正規化は必要ありません。
此の様な互いに垂直な三つの単位ベクトルを正規直交基底と言います。

計算例:カメラが𝐟^=(0.6,0,−0.8)(主に画面の奥、少し右)を向いています。

𝐟^×𝐮world=(0×0−(−0.8)×1,(−0.8)×0−0.6×0,0.6×1−0×0)=(0.8,0,0.6)‖(0.8,0,0.6)‖=0.64+0.36=1,従って 𝐫^=(0.8,0,0.6)𝐮^=𝐫^×𝐟^=(0×(−0.8)−0.6×0,0.6×0.6−0.8×(−0.8),0.8×0−0×0.6)=(0,0.36+0.64,0)=(0,1,0)

納得ですね。カメラは上にも下にも傾いていないので、其の上はワールドの上です。
単純な𝐟^=(0,0,−1)(−Zを真っ直ぐ見る)の場合は𝐫^=(1,0,0)、つまり+X軸になり、此のシリーズの座標系の通りです。

上から見たカメラの基底:直角を成す前f = 0.6, 0, −0.8と右r = 0.8, 0, 0.6

此れが第11章でビュー行列を作るlookAt関数の核心です。

前方向が上と平行になると何が起こるか。
カメラが真上を向くと𝐟^=(0,1,0)で、前方向とワールドの上が平行になり、5.3.7節で見た通り外積は(0,0,0)になります。
「右」がもうありません。真上の空を見上げると、どう回ったかによってどの方向も右になり得ます。
此のシリーズのnormalizeはゼロベクトルを返し(シェーダーならNaN)、基底が潰れ、画面にはゴミが表示されるか何も表示されません。
よくある対処法は、カメラが決して真上や真下を向かない様にする(ピッチを±89°程度にクランプする、第11章)か、前方向がワールドの上に近付き過ぎたら、(0,0,1)等の別の補助の「上」ベクトルに切り替える事です。

5.4.4 三角形の面積

外積の長さの半分は三角形の面積なので、面積が重要な時はいつでも便利です:

  • メッシュ上に物を散らす:草、パーティクル、ランダムなサンプル点をモデル全体に均等に置くには、各三角形を其の面積に比例した確率で選びます。面積が3と1の二つの三角形があれば、一つ目に75%、二つ目に25%の点を置くべきです。そうしないと、小さい三角形は混み合い、大きい三角形はスカスカに見えます。
  • モデル全体の表面積は、全ての三角形の面積の和です。
  • 重み付き法線:滑らかな頂点法線は、よく頂点の周りの正規化前の面法線の和として作られます。正規化前の外積は「法線×面積の2倍」なので、大きな三角形ほど自動的に影響が大きくなります。追加の計算は要りません。
  • メッシュの掃除:面積が(殆ど)ゼロの三角形は退化していて、使える法線が無いので、普通は削除するべきです。

5.4.5 線の左か右か、三角形の中か

2次元で点𝐚から点𝐛へ向かう線を考えます。
点𝐩は其の左にあるか、右にあるか?
線の向き𝐛−𝐚と𝐩−𝐚の2次元の外積を計算します:

s=(bx−ax)(py−ay)−(by−ay)(px−ax)

正なら左、負なら右、ゼロなら線上です。

計算例:𝐚=(1,1)から𝐛=(5,3)への線なので、𝐛−𝐚=(4,2)です。

点 𝐩−𝐚 s 結果
(2,4) (1,3) 4×3−2×1=10 左
(4,1) (3,0) 4×0−2×3=−6 右
(3,2) (2,1) 4×1−2×2=0 線上

a = 1, 1からb = 5, 3への線の左、右、線上にある点

反時計回りの三角形の三つの辺全てで此れを行うと、点が三つの辺全ての左にある時、丁度其の時だけ、点は三角形の内側にあります。
三角形(0,0),(4,0),(0,4)の場合、点(1,1)は4,8,4(全て正、内側)、(3,3)は12,−8,12(一つが負、外側)になります。
本質的には、此れがGPUのラスタライズの方法です。三角形の近くの各ピクセルについて、此の三つの「エッジ関数」を評価し、三つとも一致すればピクセルを塗ります。
ゲームでも同じ判定を、「マウスは此の2次元の図形の上にあるか?」、AIの操舵(「目標は左か右か、どちらに曲がるか?」)、多角形が凸かどうかのチェック等に使います。

5.4.6 シェーダーでのフラットシェーディング

全ての三角形が一つの平らな色を持つ、カクカクしたローポリの見た目が欲しい時もあります。
メッシュに法線が全く要らない洒落た方法があります。フラグメントシェーダーが外積で面法線を再構築するのです。

ddx/ddy(SlangとHLSL)とdFdx/dFdy(GLSL)は、或る値が此のピクセルから右隣のピクセルへ、そして次の行のピクセルへどれ位変化するかを教えてくれます。
GPUは常にピクセルを2×2の小さなブロックでシェーディングするので、此れを計算出来ます。
値がワールド座標なら、両方の差は三角形の面に沿った小さな一歩なので、其の外積が面法線になります。

struct Scene {
  float3 lightDir;  // unit vector towards the light
  float3 cameraPos;
  float3 albedo;
};
ConstantBuffer<Scene> scene;

[shader("fragment")]
float4 fragmentMain(float3 worldPos : TEXCOORD0) : SV_Target {
  // two steps along the surface -> their cross product is the face normal
  float3 n = normalize(cross(ddx(worldPos), ddy(worldPos)));
  // the sign depends on the window's y direction, so make n face the camera
  if (dot(n, scene.cameraPos - worldPos) < 0.0) n = -n;
  float diffuse = max(dot(n, scene.lightDir), 0.0);
  return float4(scene.albedo * diffuse, 1.0);
}
#version 450

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

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

void main() {
  // two steps along the surface -> their cross product is the face normal
  vec3 n = normalize(cross(dFdx(worldPos), dFdy(worldPos)));
  // the sign depends on the window's y direction, so make n face the camera
  if (dot(n, cameraPos - worldPos) < 0.0) n = -n;
  float diffuse = max(dot(n, lightDir), 0.0);
  outColor = vec4(albedo * diffuse, 1.0);
}
struct Scene {
  float3 lightDir;  // unit vector towards the light
  float3 cameraPos;
  float3 albedo;
};
ConstantBuffer<Scene> scene;

float4 fragmentMain(float3 worldPos : TEXCOORD0) : SV_Target {
  // two steps along the surface -> their cross product is the face normal
  float3 n = normalize(cross(ddx(worldPos), ddy(worldPos)));
  // the sign depends on the window's y direction, so make n face the camera
  if (dot(n, scene.cameraPos - worldPos) < 0.0) n = -n;
  float diffuse = max(dot(n, scene.lightDir), 0.0);
  return float4(scene.albedo * diffuse, 1.0);
}

何故ifがあるのでしょうか?
OpenGLではウィンドウのy軸は上向きですが、Vulkanではフレームバッファのy軸は下向きなので、「次の行のピクセル」が逆方向になり、二つのAPIで外積の符号が反転します。
どちらがどちらかを覚える代わりに、第4章の内積を使った向きの判定で、法線を単純にカメラの方へ反転させています。
フラットシェーディングのもう一つの一般的な方法は、CPUでtriangleNormal(下記)を使って面法線を計算し、補間無しで渡す事です(SlangとHLSLではnointerpolation修飾子、GLSLではflat修飾子)。

5.5 其の他の使い道

外積は、コンピューターグラフィックスが存在する遥か前から物理で使われてきました。

  • トルク:レンチを回す時。ボルトから𝐫の位置を力𝐅で押すと、トルクは𝝉=𝐫×𝐅です。其の向きはボルトが回る軸(右手の法則)で、長さはどれ位強く回るかです。長さ𝐫=(0.5,0,0)メートルのレンチを𝐅=(0,10,0)ニュートンで押すと、𝝉=(0,0,0.5×10−0×0)=(0,0,5)ニュートンメートルで、z軸周りの反時計回りの回転です。sinθの為、レンチに直角に押すのが一番効き、レンチに沿って押しても何も起こりません。
  • 角運動量𝐋=𝐫×𝐩と、回転する物体上の点の速度𝐯=𝝎×𝐫(𝝎、オメガは、回転軸を回転の速さ倍した物)。物理エンジンは剛体の為に毎フレーム両方を使います。
  • ローレンツ力:磁場𝐁の中を速度𝐯で動く電荷は、力q𝐯×𝐁を受けます。力は常に運動に垂直なので、粒子は円を描いて動きます。電気モーター、粒子加速器、昔のブラウン管テレビは全て此れを利用しています。
  • GISと地図:多角形(国、畑、建物の敷地)の面積は、隣り合う角の点の2次元の外積の和の半分です(「靴紐の公式」)。長方形(0,0),(4,0),(4,3),(0,3)の場合:0+12+12+0=24で、其の半分は12、確かに4×3です。符号は角が反時計回り(正)か時計回り(負)かを表し、多くの地図のファイル形式は此れを気にします。
  • ロボット工学:ロボットアームの各関節は軸の周りに回転し、関節が回った時に手がどれ位速く動くかは(関節の軸)×(関節から手へのベクトル)です。
  • 航法と航空:航空機の「右翼」の方向は前×上で、此のシリーズのカメラの基底と全く同じです。

5.6 コード

外積をlinalg.hhのnormalizeの後に追加します:

// 外積 / cross product (chapter 5)
inline Vec3 cross(Vec3 a, Vec3 b) {
  return {
    a.y * b.z - a.z * b.y,
    a.z * b.x - a.x * b.z,
    a.x * b.y - a.y * b.x,
  };
}

公式を一行ずつ書いただけです。「yz − zy」、「zx − xz」、「xy − yx」。
外積は3次元にしか存在しないので、Vec3版しかありません。

三角形の補助関数と2次元の左右判定は小さいので、テストプログラムの中で定義します:

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

// Unit normal of the triangle abc (counter-clockwise = front side).
Vec3 triangleNormal(Vec3 a, Vec3 b, Vec3 c) {
  return normalize(cross(b - a, c - a));
}

// Area of the triangle abc.
float triangleArea(Vec3 a, Vec3 b, Vec3 c) {
  return 0.5f * length(cross(b - a, c - a));
}

// 2D: > 0 if p is left of the line a->b, < 0 if right, 0 if on it.
float side(Vec2 a, Vec2 b, Vec2 p) {
  Vec2 ab = b - a;
  Vec2 ap = p - a;
  return ab.x * ap.y - ab.y * ap.x;
}

void print(const char *name, Vec3 v) {
  printf("%-16s = (%g, %g, %g)\n", name, v.x, v.y, v.z);
}

int main() {
  Vec3 a{1.0f, 2.0f, 3.0f};
  Vec3 b{4.0f, 5.0f, 6.0f};
  Vec3 c = cross(a, b);
  print("cross(a, b)", c);
  print("cross(b, a)", cross(b, a));
  printf("dot(c, a), dot(c, b) = %g, %g\n", dot(c, a), dot(c, b));
  print("i x j", cross({1.0f, 0.0f, 0.0f}, {0.0f, 1.0f, 0.0f}));
  print("a x 2a", cross(a, 2.0f * a));

  // Triangle with counter-clockwise vertices seen from (1, 1, 1)
  Vec3 p0{1.0f, 0.0f, 0.0f};
  Vec3 p1{0.0f, 1.0f, 0.0f};
  Vec3 p2{0.0f, 0.0f, 1.0f};
  Vec3 n = triangleNormal(p0, p1, p2);
  printf("normal           = (%.4f, %.4f, %.4f)\n", n.x, n.y, n.z);
  Vec3 nf = triangleNormal(p0, p2, p1); // swapped winding
  printf("swapped winding  = (%.4f, %.4f, %.4f)\n", nf.x, nf.y, nf.z);
  printf("area             = %.4f\n", triangleArea(p0, p1, p2));

  // Right and up vectors from a forward vector
  Vec3 forward{0.6f, 0.0f, -0.8f};
  Vec3 worldUp{0.0f, 1.0f, 0.0f};
  Vec3 right = normalize(cross(forward, worldUp));
  Vec3 up = cross(right, forward);
  print("right", right);
  print("up", up);
  print("right if f = up", normalize(cross(worldUp, worldUp)));

  // Left or right of the line from (1, 1) to (5, 3)?
  Vec2 l0{1.0f, 1.0f};
  Vec2 l1{5.0f, 3.0f};
  printf("side (2, 4)      = %g\n", side(l0, l1, {2.0f, 4.0f}));
  printf("side (4, 1)      = %g\n", side(l0, l1, {4.0f, 1.0f}));
  printf("side (3, 2)      = %g\n", side(l0, l1, {3.0f, 2.0f}));
  return 0;
}

期待される出力:

cross(a, b)      = (-3, 6, -3)
cross(b, a)      = (3, -6, 3)
dot(c, a), dot(c, b) = 0, 0
i x j            = (0, 0, 1)
a x 2a           = (0, 0, 0)
normal           = (0.5774, 0.5774, 0.5774)
swapped winding  = (-0.5774, -0.5774, -0.5774)
area             = 0.8660
right            = (0.8, -0, 0.6)
up               = (0, 1, 0)
right if f = up  = (0, 0, 0)
side (2, 4)      = 10
side (4, 1)      = -6
side (3, 2)      = 0

はい、-0と表示されています。
浮動小数点数には別の「負のゼロ」があります。rightのy成分は(−0.8)×0−0.6×0として計算され、浮動小数点数では(−0.8)×0は−0になるのです。
全ての比較で0と等しい(-0.0f == 0.0fは真)ので無害で、printfが符号を表示しているだけです。

5.7 まとめ

  • 二つの3次元ベクトルの外積𝐚×𝐛は、両方に垂直な新しいベクトルです。3次元にしか存在しません。
  • 公式:(aybz−azby,azbx−axbz,axby−aybx)。「x = yz − zy」を覚えて文字を回します。
  • 結果を確認しましょう。両方の入力との内積は0になる筈です。
  • 向き:右手の法則。𝐢×𝐣=𝐤、𝐣×𝐤=𝐢、𝐤×𝐢=𝐣。
  • 順番が大事:𝐚×𝐛=−(𝐛×𝐚)。引数を入れ替えると法線が反転します。
  • 長さ:‖𝐚‖‖𝐛‖sinθは平行四辺形の面積で、其の半分が三角形の面積です。平行なベクトルではゼロベクトルになります。
  • 2次元ではaxby−aybxが左と右(反時計回りと時計回り)を区別します。
  • グラフィックス:面法線、巡回順序と背面カリング、カメラの右・上ベクトル(前方向が上と平行な時に注意)、三角形の面積、点が三角形の内側にあるかの判定、フラットシェーディング。

5.8 課題

  1. ★ (2,0,1)×(1,3,0)を計算し、結果が両方のベクトルに垂直である事を確認して下さい。
  2. ★ 𝐢,𝐣,𝐤の規則と分配法則(𝐚+𝐛)×𝐜=𝐚×𝐜+𝐛×𝐜(外積でも成り立ちます)だけを使って(𝐢+𝐣)×𝐤を計算して下さい。次に成分の公式で答えを確認して下さい。
  3. ★ 角が(0,0,0)、(4,0,0)、(0,3,0)の三角形の面積を外積で求め、「底辺×高さ÷2」で確認して下さい。
  4. ★ 点(1,5)は、(0,0)から(3,3)への線の左と右のどちらにありますか?点(4,2)は?
  5. ★★ 三角形(0,0,0)、(2,0,0)、(0,0,−2)にtriangleNormalを使って下さい。此の床の三角形はどちらを向いていますか?(0, 1, 0)と表示される筈です(ゼロの一部が-0と表示される事があります。5.6節を参照)。次に最後の二つの頂点を入れ替えて、(0, -1, 0)と表示される事を確認して下さい。
  6. ★★ side関数を使ってbool pointInTriangle(Vec2 p, Vec2 a, Vec2 b, Vec2 c)を書いて下さい(a、b、cは反時計回りと仮定します)。三角形{0, 0}、{4, 0}、{0, 4}に対して、点{1, 1}、{3, 3}、{0.5, 3}はtrue、false、trueになる筈です。
  7. ★★★ 課題6のpointInTriangleは、三角形が時計回りで与えられると失敗します。三角形自身の辺の2次元の外積を使って巡回順序を検出し、どちらの順番でも動く様にするにはどうすれば良いでしょうか?三つの点が一直線上にある三角形の場合はどうするべきでしょうか?


コメント(0)

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