10. 行列式と逆行列

もう殆ど何でも行列に出来る様になりました：拡大縮小、回転、せん断、鏡映、そして（第9章からは）平行移動です。
此の章では、作った行列について二つの質問をします。
一つ目：此の行列は空間をどれだけ引き伸ばし、空間を裏返すのか？ 答えは行列式と言う一つの数です。
二つ目：どうやって元に戻すのか？ 答えは逆行列と言う別の行列です。

此の二つの考え方は、グラフィックスの至る所に出て来ます。
カメラのビュー行列は逆行列です。
引き伸ばしたモデルのライティングが正しく見えるのも、逆行列のお陰です。
そして鏡映したモデルは、行列式の符号を確認しないと裏返しに描画されてしまいます。
其の全てがどう動くのか見て行きましょう。

10.1 行列は空間をどれだけ引き伸ばすのか？

第6章で、行列の列は基底ベクトルの行き先だと学びました。
従って2Dでは、辺が $𝐢$ と $𝐣$ の単位正方形（面積1）は、行列の2つの列を辺とする平行四辺形に変わります。

次の行列を例にしましょう：

1列目 $𝐮=(3,1)$ は $𝐢$ の行き先、2列目 $𝐯=(1,2)$ は $𝐣$ の行き先です。
単位正方形は此の平行四辺形になります：

其の面積は5になります。
そして大事なのはここからです：行列は線型なので、全ての図形を同じ倍率で引き伸ばします。
面積2の円は面積10の楕円になります。
100ピクセルを覆うスプライトは、500ピクセルを覆う図形になります。
此の行列の「面積の倍率」が、行列式です。

10.2 公式

2×2行列の行列式は次の通りです：

括弧の代わりに縦棒で書く事もあります：

注意：此の縦棒は「絶対値」の意味ではありません。行列式は負にもなります。

3×3行列では、1行目に沿った余因子展開を使います：1行目の各成分に、其の成分の行と列を消して残った2×2行列式を掛け、符号を交互に変えて足します：

（ここのaからiは単に9つの成分の名前で、基底ベクトルとは関係ありません。）
2Dでは行列式は符号付きの面積の倍率、3Dでは符号付きの体積の倍率です。

10.3 仕組み：2×2の行列式

ad − bcは何処から来るのか

10.1の平行四辺形の面積が本当に $ad-bc=3\cdot 2-1\cdot 1=5$ なのか確かめましょう。
其の周りに、$(0,0)$ から $(4,3)$ までの長方形を描きます。
長方形の面積は $(a+b)(c+d)=4\cdot 3=12$ です。
平行四辺形は、其の長方形から角の6つのピースを除いたものです：

• 2辺が $a=3$ と $c=1$ の三角形が2つ。其れぞれ面積 $\frac{3\cdot 1}{2}=1.5$、合わせて3。
• 2辺が $b=1$ と $d=2$ の三角形が2つ。其れぞれ面積 $\frac{1\cdot 2}{2}=1$、合わせて2。
• 大きさ $b\times c=1\times 1$ の小さな長方形が2つ、合わせて2。

従って面積は $12-3-2-2=5$ です。
同じ引き算を文字で行うと、公式が出て来ます：

行列式のギャラリー

第8章の変換で試してみましょう。

回転の行は ${\sin }^2\theta +{\cos }^2\theta =1$、つまり単位円上のピタゴラスの定理（第1章）を使っています。
回転は面積を決して変えません。紙を回しても大きくはならないので、常識通りですね。
せん断はトランプの山を横に押す様なものです：山は傾きますが、カードの量は同じです。

対角行列（又は三角行列）の行列式は、単に対角成分の積です。
だから倍率 $(2,3,0.5)$ の3Dの拡大縮小行列の行列式は $2\cdot 3\cdot 0.5=3$ です。

符号：裏返し

負の面積とは何を意味するのでしょうか？
其れは、行列が鏡の様に空間の向きを裏返す事を意味します。
元の正方形では、$𝐢$ から $𝐣$ へは反時計回りに回ります。
行列式が正の行列を掛けた後も、1列目から2列目へは反時計回りの方向です。
行列式が負の行列を掛けた後は時計回りになります：図形が鏡像になっているのです。

左手と右手を思い浮かべて下さい：大きさも形も同じですが、どう回転させても片方をもう片方にする事は出来ません。
其れが出来るのは鏡だけで、鏡こそが行列式が負の行列なのです。
此れがレンダリングでとても重要な理由は10.11で見ます。

10.4 仕組み：3×3の行列式

計算例

次の行列の行列式を計算しましょう：

1行目の2、1、0に沿って進みます。
各成分について、其の行と列を消して、残りの2×2行列式を求めます（此れを其の成分の小行列式と言います）：

• 成分2（1行1列）：残りは 、小行列式は $3\cdot 4-1\cdot 0=12$。
• 成分1（1行2列）：残りは 、小行列式は $0\cdot 4-1\cdot 1=-1$。
• 成分0（1行3列）：残りは 、小行列式は $0\cdot 0-3\cdot 1=-3$。

其れらを符号 $+,-,+$ で組み合わせます：

つまり此の行列は、全ての3Dの物体の体積を25倍にし、鏡映はしません。

符号はチェッカーボード模様になっていて、対応する符号を使えばどの行や列に沿って展開しても構いません：

0が一番多い行や列に沿って展開するのがコツです。0の成分一つにつき、2×2行列式一つ分の計算が丸ごと省けます。
（小行列式に其の符号を掛けたものを余因子と言います。逆行列で此の言葉が必要になります。）

符号付き体積とスカラー三重積

3Dでは、単位立方体は行列の3つの列を辺とする斜めの箱（平行六面体）になります。
其の体積の求め方は、第4章と第5章で既に知っています：二辺の外積は両方に垂直で、長さは其の平行四辺形の面積に等しく、3つ目の辺との内積で高さが掛かります。
従って行列式は、列のスカラー三重積でもあります：

例で確かめましょう。
列は $𝐚=(2,0,1)$、$𝐛=(1,3,0)$、$𝐜=(0,1,4)$ です。

全く別の方向から、同じ答えが出ました。
符号はここでも利き手を表します：行列式が正なら3つの列は $𝐢,𝐣,𝐤$ の様な右手系、負なら左手系（鏡映）です。

10.5 行列式の規則

幾つかの規則を知っていると、行列式はずっと扱い易くなります：

• $\det (𝐈)=1$。単位行列は何も変えません。
• $\det (𝐀𝐁)=\det (𝐀)\det (𝐁)$。二つの変換を続けると、倍率が掛け合わされます。
• $\det ({𝐌}^{\mathrm{\top }})=\det (𝐌)$。行と列は対等です。
• n×n行列では $\det (s𝐌)=s^n\det (𝐌)$。3Dの物体を全方向に2倍すると、体積は $2^3=8$ 倍になります。
• 二つの列（又は二つの行）を入れ替えると符号が反転します。二つの軸を入れ替えるのは鏡映です。

積の規則を数字で確かめましょう。
（行列式5）と （行列式 $2\cdot 1-0\cdot 1=2$）を使います。

積の規則からは、鏡についての面白い事実も分かります：鏡一枚なら行列式は負ですが、鏡を二回続けると正になります。
二回鏡映するのは、単なる回転なのです。

行列式がゼロ：空間を平らに潰す

行列式が0の時、行列は平面全体を一本の線（又は一点）に、3D空間全体を平面、線、又は一点に潰してしまいます。
ギャラリーの「潰す」行列を見て下さい：

結果は全て $(1,2)$ の倍数なので、平面の全ての点が一本の線の上に乗ります。
2列目の $(2,4)$ は1列目の $(1,2)$ の丁度2倍なので、「平行四辺形」には幅が全く無いのです。

そして決定的な結果が此れです：情報が失われます。
点 $(2,-1)$ は $(2-2,4-4)=(0,0)$ に行き、$(0,0)$ も同じく $(0,0)$ に行きます。
結果の $(0,0)$ だけを渡されても、どちらから来たのか分かりません。
従って行列式が0の行列は元に戻せません：逆行列が存在しないのです。
此の様な行列を特異行列と言います。

4×4の行列式

同じ余因子の考え方は4×4行列でも使えます（4つの余因子が其れぞれ3×3行列式になります）。但し手計算では大変なので、コードに任せます。
第9章のアフィン行列には近道があります： の行列式は単に $\det (𝐀)$ です。
物体を動かしても体積は変わらないので、納得ですね。
例えば translate({1, 2, 3}) * rotateY(0.4f) * scale({2, 3, 0.5f}) の行列式は $1\cdot 1\cdot 3=3$ です（平行移動が1、回転が1、拡大縮小が3）。

10.6 逆行列：元に戻すボタン

行列 $𝐌$ の逆行列とは、其れを元に戻す行列 ${𝐌}^{-1}$ の事です：

$𝐌$ が $𝐩$ を $𝐩\prime$ にするなら、${𝐌}^{-1}$ は $𝐩\prime$ を $𝐩$ に戻します。
$𝐌$ が左に30°回転するなら、${𝐌}^{-1}$ は右に30°回転します。
$𝐌$ が大きさを2倍にするなら、${𝐌}^{-1}$ は半分にします。

数の割り算と同じ考え方です。
5の逆数は $\frac{1}{5}$ です。何故なら $5\cdot \frac{1}{5}=1$ だからです。
そして0に逆数が無い（0では割れない）のと同じ様に、行列式が0の行列には逆行列がありません。
行列では、「行列式が0でない」事が、丁度「逆行列が存在する」為の条件です。

10.7 公式：2×2の逆行列

2×2行列には短い公式があります：対角の2つの成分を入れ替え、残りの2つの符号を反転し、行列式で割ります。

仕組み

元の行列に、入れ替えて符号を反転した行列を掛けてみましょう：

単位行列の $ad-bc$ 倍になりました。
従って $ad-bc$ で割れば、丁度単位行列になります。
そして行列式が0では駄目な理由も分かります：其れで割らなければならないからです。

計算例

10.1でお馴染みの行列の逆行列を求めましょう：

3と2を入れ替え、二つの1の符号を反転し、5で割ります：

確認：掛け合わせてみます。

では実際に何かを元に戻してみましょう。
$𝐌$ は $(1,1)$ を $(3+1,1+2)=(4,3)$ に送ります。
${𝐌}^{-1}$ で戻るでしょうか？

戻りました。
又、$\det ({𝐌}^{-1})=0.4\cdot 0.6-(-0.2)\cdot (-0.2)=0.24-0.04=0.2=\frac{1}{5}$ である事にも注目して下さい：Mが面積を5倍にするなら、逆行列は5分の1にするのです。

10.8 安い逆行列

一般の逆行列には、かなりの計算が必要です。
然し、ゲームで作る行列は大抵、何で出来ているか分かっています。其の場合、逆行列は殆どタダで手に入ります。

靴下と靴

朝は靴下を履いてから、靴を履きます。
其れを元に戻すには、先に靴を脱いで、其れから靴下を脱ぎます。
行列も同じです：

証明：$(𝐀𝐁)({𝐁}^{-1}{𝐀}^{-1})=𝐀(𝐁{𝐁}^{-1}){𝐀}^{-1}=𝐀𝐈{𝐀}^{-1}=𝐀{𝐀}^{-1}=𝐈$。
（第7章の $(𝐀𝐁{)}^{\mathrm{\top }}={𝐁}^{\mathrm{\top }}{𝐀}^{\mathrm{\top }}$ と同じ、順番の反転です。）

つまり組み合わせの逆行列は、部品の逆行列さえあれば作れます。

部品

• 拡大縮小： $s$ 倍を元に戻すには $\frac{1}{s}$ 倍します。従って $𝐒(s_x,s_y,s_z{)}^{-1}=𝐒(\frac{1}{s_x},\frac{1}{s_y},\frac{1}{s_z})$。予想通り、倍率0は元に戻せません。
• 平行移動： $𝐭$ だけ動かしたのを元に戻すには $-𝐭$ だけ動かします。従って $𝐓(𝐭{)}^{-1}=𝐓(-𝐭)$。
• 回転： $\theta$ の回転を元に戻すには $-\theta$ 回転します。そして更に素敵な言い方があります：${𝐑}^{-1}={𝐑}^{\mathrm{\top }}$。

何故回転の逆行列は転置なのか

回転行列の列は、回転した基底ベクトルです。
回転は長さも角度も変えないので、列は単位長さ（第3章）のままで、互いに垂直（内積0、第4章）のままです。
では積 ${𝐑}^{\mathrm{\top }}𝐑$ を見てみましょう：$(i,j)$ 成分は ${𝐑}^{\mathrm{\top }}$ の i 行目と $𝐑$ の j 列目の内積で、${𝐑}^{\mathrm{\top }}$ の i 行目は $𝐑$ の i 列目そのものです。
つまり全ての成分は二つの列の内積で、列が自分自身と出会えば1、別の列と出会えば0になります。
此れは正に単位行列です。

2Dの回転で確かめましょう（$c=\cos \theta$、$s=\sin \theta$ と書きます）：

列が単位長さで互いに垂直な行列を直交行列と言い、全ての直交行列で逆行列は単に転置です。
回転は行列式1の直交行列、鏡映は行列式−1の直交行列です。
転置には計算が全く要りません。数字を別の順番で読むだけです。だからゲームエンジンは回転が大好きなのです。

モデル行列全体を元に戻す

靴下と靴で全てを組み合わせると：

右から読むと：位置を戻し、回転を戻し、拡大縮小を戻します。
第9章の板で試しましょう：$𝐓(5,0,0)$、$𝐑$ はZ軸周りの90°、$𝐒(2,1,1)$ です。
此のモデル行列は、ローカルの先端 $(1,0,0)$ をワールドの点 $(5,2,0)$ に送りました。
其れを送り返します：

1. $𝐓(-5,0,0)$：$(5,2,0)\to (0,2,0)$。
2. ${𝐑}^{\mathrm{\top }}$（−90°の回転で、$(x,y)$ を $(y,-x)$ に送る）：$(0,2,0)\to (2,0,0)$。
3. $𝐒(\frac{1}{2},1,1)$：$(2,0,0)\to (1,0,0)$。

先端の $(1,0,0)$ に戻りました。

10.9 一般の逆行列：余因子と余因子行列

行列が何で出来ているか分からない時は、一般的な方法が必要です。
2×2の公式は、実は一般的なレシピの小さな例です：

1. 全ての成分について、其の余因子（10.4のチェッカーボードの符号を付けた小行列式）を求めます。
2. 余因子を行列に並べて、転置します。此れを余因子行列（随伴行列）$\mathrm{adj}(𝐌)$ と言います。
3. 行列式で割ります。

2×2で確かめましょう：a の余因子は d、b の余因子は $-c$、c の余因子は $-b$、d の余因子は a です。
転置すると 、10.7の行列と全く同じです。

4×4行列では、linalg.hh は此のレシピを賢い近道と一緒に使っています。
4×4行列の余因子は全て3×3行列式で、3×3行列式は全て2×2行列式で出来ています。
其の2×2行列式の多くは同じものなので、コードは其れらを一度だけ計算します：

• s0 〜 s5 は、上の2行から作れる6つの2×2行列式です（列の組ごとに一つ：0-1、0-2、0-3、1-2、1-3、2-3）。
• c0 〜 c5 は、下の2行から作れる6つの2×2行列式です。

行列式は6つの積の和で、上の2行の列の組を、残りの2列を使う下の2行の組と組み合わせます：s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0。
其の後、逆行列の16個の成分は其れぞれ、此の12個の数から作った余因子を行列式で割ったものです。
其れらの行を覚える必要は全くありません。「余因子を転置して行列式で割る」を効率良くやっているだけだと分かれば十分です。

実用上の注意が二つあります：

• 浮動小数点数では、行列式が丁度0になる事は滅多にありません。ヘッダーは ${10}^{-8}$ より小さい値を0として扱い、省略可能な ok フラグで其れを知らせ、無限大だらけの行列の代わりに単位行列を返します。
• 一般の4×4逆行列には100回近くの掛け算が必要です（ヘッダーの版では94回）。CPUで物体ごと・フレームごとに一回なら何でもありませんが、出来ればシェーダーで頂点ごとに計算するのは避けましょう。

10.10 法線行列

殆どの初心者が混乱する問題がここにあります。
第9章では、法線をモデル行列と $w=0$ で変換しました。
回転と均一な拡大縮小なら其れで正しいのですが、不均一な拡大縮小だと法線がおかしくなります：面に垂直ではなくなってしまうのです。

2Dの計算例

斜め上に向かう面を考えます。面に沿った方向（此れを接線と言います）は $𝐭=(1,1)$ です。
其の法線は $𝐧=(1,-1)$ で、確かに $𝐧\cdot 𝐭=1-1=0$ です。
此のモデルをx方向に2倍に引き伸ばします：。

接線は面の一部なので、単純に引き伸ばされます：$𝐭\prime =𝐒𝐭=(2,1)$。
法線を同じ様に変換すると $𝐒𝐧=(2,-1)$ になり：

垂直ではありません！
面は平らになったのに、素朴に変換した法線は逆の方へ傾いてしまいました。
此の法線でライティングすると見た目がおかしくなり、引き伸ばす程悪化します。

今度は逆行列の転置を使ってみましょう。
で、対角行列は転置しても同じなので：

又垂直になりました。
x方向に引き伸ばすと面は平らになるので、法線のx方向の傾きは小さくならなければなりません。2で割るのは正に其れです。

何故逆行列の転置で上手く行くのか

先ず小さな道具が一つ必要です：内積は行列の積、つまり行と列の積として書けます：

（1×n行列とn×1行列の積は1×1行列、つまり単なる数です。）
接線は $𝐌$ で変換され、法線は未知の行列 $𝐍$ で変換されるとします。
第7章の $(𝐀𝐁{)}^{\mathrm{\top }}={𝐁}^{\mathrm{\top }}{𝐀}^{\mathrm{\top }}$ を使うと：

${𝐍}^{\mathrm{\top }}𝐌=𝐈$ となる様に $𝐍$ を選べば、真ん中が消えて ${𝐧}^{\mathrm{\top }}𝐭=𝐧\cdot 𝐭=0$ が残ります。
つまり ${𝐍}^{\mathrm{\top }}={𝐌}^{-1}$ が必要で、其れは次の意味です：

此れを法線行列と言います。
法線は方向なので、関係するのは左上の3×3の部分（線型部分 $𝐀$）だけです。
アフィンの4×4行列なら、transpose(inverse(model)) の左上3×3は丁度 $({𝐀}^{-1}{)}^{\mathrm{\top }}$ なので、4×4のまま計算して法線に $w=0$ で掛けても大丈夫です。
法線行列は長さを保たないので、結果は必ず正規化して下さい。

省略出来る時

• 回転だけ： $({𝐑}^{-1}{)}^{\mathrm{\top }}=({𝐑}^{\mathrm{\top }}{)}^{\mathrm{\top }}=𝐑$。法線行列はモデル行列そのものです。
• 回転と均一な拡大縮小 $s$：法線行列は $\frac{1}{s}𝐑$ で、$s𝐑$ と同じ方向を向きます。正規化すれば、普通のモデル行列でも同じ答えになります。
• 不均一な拡大縮小（又はせん断）： 本物の法線行列が必要です。

多くのエンジンは、物体の拡大縮小が均一かどうかを確認し、均一なら逆行列の計算を省きます。

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

頂点シェーダーの法線行列

法線行列（transpose(inverse(model))）はCPUで物体ごとに一度計算し、モデル行列と一緒にuniformとして送ります。
GLSLには組み込みの inverse 関数があるので、シェーダーの中で直接 transpose(inverse(mat3(model))) と書いてあるのをよく見掛けます。然し其れだと、全く同じ逆行列を頂点ごとに計算し直す事になります。
uniformとして渡す方が速く、SlangでもGLSLでもHLSLでも同じ様に使えます。

struct Matrices {
  float4x4 model;
  float4x4 view;
  float4x4 proj;
  float4x4 normalMatrix; // transpose(inverse(model)), computed on the CPU
};
ConstantBuffer<Matrices> matrices;

struct VSInput {
  float3 position : POSITION;
  float3 normal   : NORMAL;
};

struct VSOutput {
  float4 position    : SV_Position;
  float3 worldNormal : NORMAL;
};

[shader("vertex")]
VSOutput vertexMain(VSInput input) {
  VSOutput output;
  float4 worldPos = mul(matrices.model, float4(input.position, 1.0));
  output.position = mul(matrices.proj, mul(matrices.view, worldPos));
  // w = 0: a direction. Normalize, because the normal matrix changes lengths.
  output.worldNormal = normalize(mul(matrices.normalMatrix, float4(input.normal, 0.0)).xyz);
  return output;
}

#version 450

layout(set = 0, binding = 0) uniform Matrices {
  mat4 model;
  mat4 view;
  mat4 proj;
  mat4 normalMatrix; // transpose(inverse(model)), computed on the CPU
};

layout(location = 0) in vec3 inPosition;
layout(location = 1) in vec3 inNormal;

layout(location = 0) out vec3 worldNormal;

void main() {
  vec4 worldPos = model * vec4(inPosition, 1.0);
  gl_Position = proj * (view * worldPos);
  // w = 0: a direction. Normalize, because the normal matrix changes lengths.
  worldNormal = normalize((normalMatrix * vec4(inNormal, 0.0)).xyz);
}

struct Matrices {
  float4x4 model;
  float4x4 view;
  float4x4 proj;
  float4x4 normalMatrix; // transpose(inverse(model)), computed on the CPU
};
ConstantBuffer<Matrices> matrices;

struct VSInput {
  float3 position : POSITION;
  float3 normal   : NORMAL;
};

struct VSOutput {
  float4 position    : SV_Position;
  float3 worldNormal : NORMAL;
};

VSOutput vertexMain(VSInput input) {
  VSOutput output;
  float4 worldPos = mul(matrices.model, float4(input.position, 1.0));
  output.position = mul(matrices.proj, mul(matrices.view, worldPos));
  // w = 0: a direction. Normalize, because the normal matrix changes lengths.
  output.worldNormal = normalize(mul(matrices.normalMatrix, float4(input.normal, 0.0)).xyz);
  return output;
}

フラグメントシェーダーでも、もう一度正規化して下さい：三角形の上で単位ベクトル同士を補間すると、長さが1より少し短いベクトルになるからです。

ビュー行列は逆行列

カメラも世界の中の物体の一つで、位置と回転を持つので、モデルと同じ様にワールド行列があります。
然しレンダリングでは逆方向が必要です：全ての物がカメラから見て何処にあるかを知りたいのです。
其れが逆行列です：

例えば、$(0,0,5)$ にあって回転していないカメラのワールド行列は $𝐓(0,0,5)$ なので、ビュー行列は $𝐓(0,0,-5)$ です。
世界の原点はビュー空間で $(0,0,-5)$ に来ます：−Z方向を見ているカメラの5つ前です。
第11章では此れを lookAt でちゃんと作ります。lookAt は正に「回転を転置して、平行移動の符号を反転する」安い逆行列を使っています。

鏡映したモデルは頂点の巻き順が反転する

第5章で、GPUは三角形のどちらの面がカメラを向いているかを、画面上での頂点の巻き順（反時計回りか時計回りか）で決めると学びました。
行列式が負の行列はモデルを鏡映し、鏡映は全ての反時計回りの三角形を時計回りにしてしまいます。
例えば三角形 $(0,0),(1,0),(0,1)$ は反時計回りです：二辺の外積のzは $1\cdot 1-0\cdot 0=1>0$ です。
xを鏡映すると $(0,0),(-1,0),(0,1)$ になり、外積のzは ：時計回りです。

従って右手から左手を作ろうと scale({-1.0f, 1.0f, 1.0f}) を使うと、背面カリングが突然モデルの表側を消してしまい、裏返しに見えます。
直し方は、CPUで行列式を確認して、其の物体だけ表面の設定を入れ替える事です：

// 鏡映された物体か / is this object mirrored?
inline bool isMirrored(const Mat4 &model) {
  return determinant(model) < 0.0f;
}

OpenGLでは、鏡映した物体を描画する間だけ glFrontFace(GL_CCW) の代わりに glFrontFace(GL_CW) を呼びます（又は、普段使っている表面の設定を逆にします）。
Vulkanでは表面の設定はパイプラインのラスタライゼーションステートの一部なので、VK_FRONT_FACE_CLOCKWISE と VK_FRONT_FACE_COUNTER_CLOCKWISE を入れ替えた二つ目のパイプラインを用意するか、ダイナミックステートにして vkCmdSetFrontFace（Vulkan 1.3からコア）を呼びます。

ワールドから物体自身の空間へ

モデル行列の逆行列は、ワールドの位置を物体のローカル空間へ持って行きます。
其れで多くの問題がずっと簡単になります。
例えば：プレイヤーは此の回転したトリガーボックスの中にいるか？
ワールド空間では箱が回転しているので、判定が面倒です。
自分のローカル空間なら、箱は $-𝐡$ から $𝐡$（半分の大きさ）までの単純な軸に沿った箱なので、座標を比べるだけです：

// 回転した箱の中に点があるか / is a point inside a rotated (and moved) box?
bool insideBox(const Mat4 &boxModel, Vec3 halfSize, Vec3 worldPoint) {
  Vec4 p = inverse(boxModel) * Vec4{worldPoint.x, worldPoint.y, worldPoint.z, 1.0f};
  return std::fabs(p.x) <= halfSize.x &&
         std::fabs(p.y) <= halfSize.y &&
         std::fabs(p.z) <= halfSize.z;
}

半分の大きさが $(2,1,1)$ で、Y軸周りに90°回転させて $(10,0,0)$ に置いた箱を使いましょう。長い辺は今、ワールドのZ方向に伸びています。
$(10.5,0,1.5)$ にいるプレイヤーは、$-(10,0,0)$ だけ戻されて $(0.5,0,1.5)$ になり、次にY軸周りに−90°回転を戻されて $(-1.5,0,0.5)$ になります。
3つの座標が全て半分の大きさの内側なので、プレイヤーは中にいます。
もし回転を忘れて $(0.5,0,1.5)$ を其のまま比べたら、$1.5>1$ なので間違って「外」と判定してしまいます。
点 $(12,0,0)$ はローカル空間で $(0,0,2)$ になり、正しく外と判定されます。箱は其の方向には中心から1しか伸びていないからです。

逆行列に出会う其の他の場所

• マウスピッキング： マウスカーソルの下に何があるかを調べるには、投影行列とビュー行列を逆向きに辿ります。$𝐏𝐕$ の逆行列を使います（第12章）。
• シャドウマップとデカール： ワールドの位置を、光源やデカールの投影機の空間へ、其のワールド行列の逆行列で持って行きます。
• 潰れた拡大縮小： 物体を丁度0まで縮めるアニメーションは、行列式を0にします。逆行列（と法線行列）が存在しなくなり、ライティングがNaNになる原因としてよくあります。倍率を小さな正の値で止めるか、其の物体を描画しない様にしましょう。

10.12 其の他の使い道

• 連立方程式を解く： 方程式 $3x+y=9$ と $x+2y=8$ は、お馴染みの と $𝐛=(9,8)$ を使った行列の方程式 $𝐌𝐱=𝐛$ です。すると $𝐱={𝐌}^{-1}𝐛=(0.4\cdot 9-0.2\cdot 8,-0.2\cdot 9+0.6\cdot 8)=(2,3)$。（実際のソルバーは逆行列全体ではなく消去法を使いますが、考え方は同じです。）
• 解が一つだけか： 正方行列の連立方程式 $𝐌𝐱=𝐛$ は、$\det (𝐌)\ne 0$ の時、そして其の時に限り、解を丁度一つ持ちます。行列式が0なら、解が無いか、無限に有るかのどちらかです。
• クラメルの公式： 各未知数は行列式の比で表せます。其の未知数の列を $𝐛$ に置き換えるのです。上の方程式なら $x=\frac{9\cdot 2-1\cdot 8}{5}=2$、$y=\frac{3\cdot 8-9\cdot 1}{5}=3$。2×2や3×3を手で解く時に便利です。
• 暗号： ヒル暗号は、文字のまとまりに鍵行列を掛けて（26を法として）暗号化します。復号には逆行列が必要なので、鍵の行列式は26と公約数を持ってはいけません。
• 電子工学： 回路解析ではキルヒホッフの法則を電圧と電流の行列の方程式にして、逆行列（又は消去法）で解きます。
• 経済学： レオンチェフの産業連関モデルは、各産業がどれだけ生産すべきかを $𝐱=(𝐈-𝐀{)}^{-1}𝐝$ で計算します。
• 統計と機械学習： 線型回帰の古典的な公式には $({𝐗}^{\mathrm{\top }}𝐗{)}^{-1}$ が入っていますし、共分散行列の行列式はデータがどれだけ「広がっているか」を表します。

10.13 コード

此の章で linalg.hh に追加する物です。ヘッダーにある通りです。
どちらの関数も、10.9で説明した12個の2×2行列式を共有しています。

// 行列式（余因子展開） / determinant (cofactor expansion) (chapter 10)
inline float determinant(const Mat4 &m) {
  auto a = [&](int r, int c) { return m.at(r, c); };
  float s0 = a(0,0) * a(1,1) - a(1,0) * a(0,1);
  float s1 = a(0,0) * a(1,2) - a(1,0) * a(0,2);
  float s2 = a(0,0) * a(1,3) - a(1,0) * a(0,3);
  float s3 = a(0,1) * a(1,2) - a(1,1) * a(0,2);
  float s4 = a(0,1) * a(1,3) - a(1,1) * a(0,3);
  float s5 = a(0,2) * a(1,3) - a(1,2) * a(0,3);
  float c5 = a(2,2) * a(3,3) - a(3,2) * a(2,3);
  float c4 = a(2,1) * a(3,3) - a(3,1) * a(2,3);
  float c3 = a(2,1) * a(3,2) - a(3,1) * a(2,2);
  float c2 = a(2,0) * a(3,3) - a(3,0) * a(2,3);
  float c1 = a(2,0) * a(3,2) - a(3,0) * a(2,2);
  float c0 = a(2,0) * a(3,1) - a(3,0) * a(2,1);
  return s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0;
}

// 逆行列 / inverse matrix (chapter 10)
// 行列式が0なら逆行列は存在しない。okにfalseを入れて単位行列を返す。
// If the determinant is 0 there is no inverse: sets ok to false and returns identity.
inline Mat4 inverse(const Mat4 &m, bool *ok = nullptr) {
  auto a = [&](int r, int c) { return m.at(r, c); };
  float s0 = a(0,0) * a(1,1) - a(1,0) * a(0,1);
  float s1 = a(0,0) * a(1,2) - a(1,0) * a(0,2);
  float s2 = a(0,0) * a(1,3) - a(1,0) * a(0,3);
  float s3 = a(0,1) * a(1,2) - a(1,1) * a(0,2);
  float s4 = a(0,1) * a(1,3) - a(1,1) * a(0,3);
  float s5 = a(0,2) * a(1,3) - a(1,2) * a(0,3);
  float c5 = a(2,2) * a(3,3) - a(3,2) * a(2,3);
  float c4 = a(2,1) * a(3,3) - a(3,1) * a(2,3);
  float c3 = a(2,1) * a(3,2) - a(3,1) * a(2,2);
  float c2 = a(2,0) * a(3,3) - a(3,0) * a(2,3);
  float c1 = a(2,0) * a(3,2) - a(3,0) * a(2,2);
  float c0 = a(2,0) * a(3,1) - a(3,0) * a(2,1);
  float det = s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0;

  if (std::fabs(det) < 1e-8f) {
    if (ok) *ok = false;
    return Mat4::identity();
  }
  if (ok) *ok = true;

  float inv = 1.0f / det;
  Mat4 r;
  r.at(0,0) = ( a(1,1) * c5 - a(1,2) * c4 + a(1,3) * c3) * inv;
  r.at(0,1) = (-a(0,1) * c5 + a(0,2) * c4 - a(0,3) * c3) * inv;
  r.at(0,2) = ( a(3,1) * s5 - a(3,2) * s4 + a(3,3) * s3) * inv;
  r.at(0,3) = (-a(2,1) * s5 + a(2,2) * s4 - a(2,3) * s3) * inv;
  r.at(1,0) = (-a(1,0) * c5 + a(1,2) * c2 - a(1,3) * c1) * inv;
  r.at(1,1) = ( a(0,0) * c5 - a(0,2) * c2 + a(0,3) * c1) * inv;
  r.at(1,2) = (-a(3,0) * s5 + a(3,2) * s2 - a(3,3) * s1) * inv;
  r.at(1,3) = ( a(2,0) * s5 - a(2,2) * s2 + a(2,3) * s1) * inv;
  r.at(2,0) = ( a(1,0) * c4 - a(1,1) * c2 + a(1,3) * c0) * inv;
  r.at(2,1) = (-a(0,0) * c4 + a(0,1) * c2 - a(0,3) * c0) * inv;
  r.at(2,2) = ( a(3,0) * s4 - a(3,1) * s2 + a(3,3) * s0) * inv;
  r.at(2,3) = (-a(2,0) * s4 + a(2,1) * s2 - a(2,3) * s0) * inv;
  r.at(3,0) = (-a(1,0) * c3 + a(1,1) * c1 - a(1,2) * c0) * inv;
  r.at(3,1) = ( a(0,0) * c3 - a(0,1) * c1 + a(0,2) * c0) * inv;
  r.at(3,2) = (-a(3,0) * s3 + a(3,1) * s1 - a(3,2) * s0) * inv;
  r.at(3,3) = ( a(2,0) * s3 - a(2,1) * s1 + a(2,2) * s0) * inv;
  return r;
}

ok ポインターは省略可能です：逆行列が存在するか知りたい時は &ok を渡し、存在すると確信している時は省略します。

テストプログラムです。
此の章の行列式を確認し、model * inverse(model) が単位行列になる事を確かめ（浮動小数点の小さな誤差がある為、== ではなく許容誤差で比べます）、第9章の板を元に戻し、10.10の法線を3Dで計算します。
radians() は第1章の物で、tidy は浮動小数点の計算が出したがる -0.00 を隠します。

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

static float tidy(float v) { return std::fabs(v) < 1e-5f ? 0.0f : v; }

static void printMat(const char *name, const Mat4 &m) {
  printf("%s\n", name);
  for (int row = 0; row < 4; ++row) {
    printf("  %7.3f %7.3f %7.3f %7.3f\n",
           tidy(m.at(row, 0)), tidy(m.at(row, 1)), tidy(m.at(row, 2)), tidy(m.at(row, 3)));
  }
}

// 単位行列との最大誤差 / largest difference from the identity matrix
static float errorFromIdentity(const Mat4 &m) {
  Mat4 id = Mat4::identity();
  float worst = 0.0f;
  for (int row = 0; row < 4; ++row)
    for (int col = 0; col < 4; ++col)
      worst = std::fmax(worst, std::fabs(m.at(row, col) - id.at(row, col)));
  return worst;
}

int main() {
  // 1. Determinants
  Mat4 model = translate({1.0f, 2.0f, 3.0f}) * rotateY(0.4f) * scale({2.0f, 3.0f, 0.5f});
  printf("det(model)  = %.3f\n", determinant(model));
  printf("det(mirror) = %.3f\n", determinant(scale({-1.0f, 1.0f, 1.0f})));
  printf("det(flat)   = %.3f\n", determinant(scale({1.0f, 0.0f, 1.0f})));

  // 2. Inverse and the M * M^-1 = I check
  bool ok = false;
  Mat4 inv = inverse(model, &ok);
  printf("ok = %s, error = %s\n", ok ? "true" : "false",
         errorFromIdentity(model * inv) < 1e-5f ? "tiny" : "LARGE");

  inverse(scale({1.0f, 0.0f, 1.0f}), &ok);
  printf("flat matrix: ok = %s\n", ok ? "true" : "false");

  // 3. Undoing a TRS matrix
  Mat4 m = translate({5.0f, 0.0f, 0.0f}) * rotateZ(radians(90.0f)) * scale({2.0f, 1.0f, 1.0f});
  printMat("inverse(TRS)", inverse(m));
  Vec4 back = inverse(m) * Vec4{5.0f, 2.0f, 0.0f, 1.0f};
  printf("back = (%.2f, %.2f, %.2f, %.2f)\n", tidy(back.x), tidy(back.y), tidy(back.z), tidy(back.w));

  // 4. Normal matrix
  Mat4 stretch = scale({2.0f, 1.0f, 1.0f});
  Mat4 normalMatrix = transpose(inverse(stretch));
  Vec4 n = normalMatrix * Vec4{1.0f, -1.0f, 0.0f, 0.0f};
  printf("normal = (%.2f, %.2f, %.2f)\n", tidy(n.x), tidy(n.y), tidy(n.z));
  return 0;
}

期待される出力：

det(model)  = 3.000
det(mirror) = -1.000
det(flat)   = 0.000
ok = true, error = tiny
flat matrix: ok = false
inverse(TRS)
    0.000   0.500   0.000   0.000
   -1.000   0.000   0.000   5.000
    0.000   0.000   1.000   0.000
    0.000   0.000   0.000   1.000
back = (1.00, 0.00, 0.00, 1.00)
normal = (0.50, -1.00, 0.00)

表示された逆行列を見て下さい：左上の部分は ${𝐒}^{-1}{𝐑}^{\mathrm{\top }}$（転置した回転の1行が半分になったもの）で、最後の列 $(0,5,0)$ は $-{𝐒}^{-1}{𝐑}^{\mathrm{\top }}𝐭$ です。
10.8の安い公式を、高い方法で計算した結果と全く同じです。
法線 $(0.5,-1,0)$ はまだ単位長さではありません。シェーダーでは正規化します。

10.14 まとめ

• 行列式は、行列による符号付きの面積（2D）又は体積（3D）の倍率です。2×2なら $ad-bc$、3×3なら余因子展開か、列のスカラー三重積 $𝐚\cdot (𝐛\times 𝐜)$ を使います。
• 行列式が負 = 鏡映（向きが反転）。行列式がゼロ = 空間が平らに潰れ、情報が失われ、逆行列が無い。
• $\det (𝐀𝐁)=\det (𝐀)\det (𝐁)$。回転の行列式は1で、平行移動は行列式を変えません。
• 逆行列は行列を元に戻します：$𝐌{𝐌}^{-1}=𝐈$。2×2なら、対角を入れ替え、残りの符号を反転し、行列式で割ります。
• 靴下と靴：$(𝐀𝐁{)}^{-1}={𝐁}^{-1}{𝐀}^{-1}$。安い部品：拡大縮小は $\frac{1}{s}$、回転は ${𝐑}^{\mathrm{\top }}$、平行移動は $-𝐭$。従って $(𝐓𝐑𝐒{)}^{-1}={𝐒}^{-1}{𝐑}^{\mathrm{\top }}{𝐓}^{-1}$。
• 一般の逆行列は余因子行列（転置した余因子）を行列式で割ったもので、inverse() は失敗を ok で知らせます。
• 法線は法線行列 $({𝐌}^{-1}{)}^{\mathrm{\top }}$（左上の3×3）で変換し、正規化し直します。回転と均一な拡大縮小だけなら、モデル行列で十分です。
• グラフィックスでは：ビュー行列は逆行列、行列式が負なら表面の設定を入れ替える、モデル行列の逆行列はワールドの位置を物体のローカル空間へ持って行く。

10.15 課題

1. ★ と の行列式を計算して下さい。其れぞれ単位正方形の面積をどうしますか？10と0になるはずです。2つ目の行列の列について、0が何を意味するか説明して下さい。
2. ★ の行列式を、1行目に沿った余因子展開で計算し、次に1列目に沿って計算し直して下さい。どちらも22になるはずです。
3. ★ 2×2の公式で の逆行列を求め、元の行列との積が単位行列になる事を確認して下さい。 になるはずです。
4. ★ の逆行列を求めて下さい（行列式が1なので、逆行列は整数だけになります）。掛け算で答えを確認して下さい。
5. ★ 2Dの面の接線が $(1,2)$、法線が $(2,-1)$ です。其れを で拡大縮小します。素朴な法線 $(6,-1)$ が新しい接線に垂直でない事（内積は16）、そして法線行列を使えば垂直な法線が得られる事を示して下さい。
6. ★★ scale、transpose、rotateY、translate だけを使って ${𝐒}^{-1}{𝐑}^{\mathrm{\top }}𝐓(-𝐭)$ を返す Mat4 inverseTRS(Vec3 t, float angleY, Vec3 s) を書いて下さい。t = {1, 2, 3}、angleY = 0.4f、s = {2, 3, 0.5f} に対して、inverse(translate(t) * rotateY(angleY) * scale(s)) と16個の成分全てを比べて下さい。最大の差は0.00001未満になるはずです。
7. ★★★ アフィン行列 の逆行列が である事を示して下さい（第9章のブロックの公式で掛け算します）。次に、$𝐀$ が回転であるカメラの場合に簡単にして下さい。新しい平行移動の列の各成分は内積になります。どのベクトル同士の内積でしょうか？（第11章の lookAt で又出会います。）