21. メッシュ処理と符号付き距離場

第1部では、面を表す方法を2つ見ました。
ポリゴンメッシュ（第9章、位相は第7章）と、陰関数$f(𝐩)=0$です。
此の章は、其の両方を扱う為の章です。
メッシュについては、モデリングツール、ゲームエンジン、3Dスキャナが裏で動かしているアルゴリズムを見ます。ノイズの有るスキャンの平滑化、LODの為の簡略化、壊れた形状の修復、UVレイアウトの為の平面への展開、そして法線マップへのディテールのベイクです。
陰関数曲面については符号付き距離場（SDF）を扱います。スフィアトレーシングが何故面を見逃さないのかを確かめ、形状を滑らかに融合し、メッシュからSDFを計算し、小さなテクスチャから鮮明な文字を描き、SDFをメッシュに戻します。

ゲームや映画のメッシュで、最初に作られた儘の形で描かれる物は殆ど有りません。
又、SDFはソフトシャドウ、大域照明、衝突判定、文字の描画の為にリアルタイムエンジンの一部になっています。
第9章は第2部に2つの事を約束しました。簡略化の為の二次誤差尺度の導出と、SDFのしっかりした扱いです。其れが此の章です。
内容は検定のシラバスを超えますが、其の中の3つの項目を深めます。ポリゴン曲面の表現、陰関数曲面とマーチングキューブ法、そしてマッピングです。
数学は線型代数シリーズの内積・外積と小さな行列、そして勾配を使います。勾配は途中で説明します。

21.1 離散的な曲面としてのメッシュ

滑らかな曲面には全ての点に法線、曲率、ラプラス作用素が有りますが、三角形メッシュに有るのは頂点、辺、平らな面だけです。
離散微分幾何は、此れらの量のメッシュ版を定義します。其れらは滑らかな場合と同じ様に振る舞い、メッシュが細かくなる程正確になります。
殆どのメッシュアルゴリズムは、其の中の2つ、ラプラシアンと曲率から作られています。
どちらも頂点の1-リング近傍を見ます。辺で繋がった隣の頂点の集合$N(i)$と、周りの面です。第7章のハーフエッジ構造は、正に其れを辿る為の物です。
隣の頂点の数$|N(i)|$が価数で、規則正しい三角形メッシュでは6です。

離散ラプラシアン：アンブレラ演算子

最も単純な離散ラプラシアンは、頂点と隣の頂点の平均を比べます：

隣の頂点が傘の骨の様に頂点の周りに広がるので、アンブレラ演算子と呼ばれます。${𝐜}_i$は其れらの重心です：

此れは第15章（15.5節）のラプラシアンフィルタのメッシュ版です。あちらは画素を4つの隣の画素と比べました。
平らで規則正しい部分では$𝐋=\mathrm{𝟎}$、出っ張りの上では$𝐋$は出っ張りの内側、ほぼ法線の方向を向きます。

計算例。 頂点${𝐯}_i=(0,0.3,0)$が平らな床から0.3上に有り、隣の頂点は$(1,0,0)$、$(0,0,-1)$、$(-1,0,0)$、$(0,0,1)$の4つです。
重心は$\frac{1}{4}(1+0-1+0,0,0-1+0+1)=(0,0,0)$なので、$𝐋({𝐯}_i)=(0,-0.3,0)$です。真下向きで、長さは出っ張りの高さです。

アンブレラ演算子は三角形の形を無視するので、不規則なメッシュでは頂点を横にも引っ張ります。
幾何学的に正しい版は、隣の頂点$j$にコタンジェント重み（Pinkall・Polthier、1993年）$w_{ij}=\frac{1}{2}(\cot {\alpha }_{ij}+\cot {\beta }_{ij})$を掛けます。${\alpha }_{ij}$と${\beta }_{ij}$は辺$ij$の向かいに有る2つの角です。
すると$𝐋$は正確に法線方向を向き、其の長さ（頂点の面積で割った物）は平均曲率、つまり面の曲がり方の平均の2倍を近似します。

角度から曲率を：角度欠損

頂点の周りの面を机の上に平らに並べてみましょう。
平らな面では角の合計が360°になり、隙間無く平らに並びます。
立方体の角では3つの直角の合計が270°しか無く、90°の隙間が残ります。鞍点では角が多過ぎて、紙に皺が寄ってしまいます：

此の足りない角が角度欠損で、頂点に集中した離散的なガウス曲率です：

${\theta }_f$は面$f$の其の頂点での角です。
正なら球の様な点、0なら平ら又は一方向にだけ曲がった点（円柱は平らに広げられます）、負なら鞍の様な点です。

• 立方体の角：$2\pi -3\cdot \frac{\pi }{2}=\frac{\pi }{2}$（90°）。
• 正八面体の頂点、正三角形4つ：$2\pi -4\cdot \frac{\pi }{3}=\frac{2\pi }{3}$（120°）。
• 平らな格子の頂点、正三角形6つ：$2\pi -6\cdot \frac{\pi }{3}=0$。
• 正方形5つが集まる鞍点：$2\pi -5\cdot \frac{\pi }{2}=-\frac{\pi }{2}$。

では閉じたメッシュ全体で足してみましょう。立方体の8つの角で$8\cdot \frac{\pi }{2}=4\pi$、正八面体の6つの頂点で$6\cdot \frac{2\pi }{3}=4\pi$です。
此れはデカルトの定理、離散版のガウス・ボネの定理です：

そして穴の無い閉じたメッシュでは$V-E+F=2$です（オイラーの公式、第7章）。
全曲率は位相だけで決まる定数です。頂点を動かすと曲率は移動しますが、作ったり消したりは出来ません。
トーラス（$V-E+F=0$）の全曲率は0で、外側が正、内側が負です。
実際には、角度欠損は鋭い角の検出、スキャンのノイズの検出、そしてテクスチャアーティストが縁の擦り減りに使う「曲率マップ」の作成に使われます。

21.2 平滑化とフェアリング

3Dスキャン、マーチングキューブ法の面、フォトグラメトリのメッシュには、必ず少しノイズが有ります。
平滑化（フェアリングとも）は、画像のぼかしフィルタの様に、高周波の凸凹を取り除き、全体の形を残します。

ラプラシアンスムージング

全ての頂点を、隣の頂点の重心に向かって割合$\lambda$だけ同時に動かし（新しい位置を全て古い位置から計算してから置き換えます）、此れを繰り返します：

此れがラプラシアンスムージングです。
上の例で$\lambda =0.5$とすると、出っ張りの頂点は$(0,0.3,0)+0.5\cdot (0,-0.3,0)=(0,0.15,0)$に移ります。1回で出っ張りが半分になります。

何故縮むのか

ラプラシアンスムージングには有名な欠点が有ります。ノイズだけでなく、全てが縮むのです。
ノイズの無い半径1の正八角形を考えます。各点の隣の2点は45°離れています。
其の平均は中心からの同じ半直線上、距離$\cos 45\mathrm{°}\approx 0.7071$に有るので、アンブレラベクトルは長さ$1-0.7071=0.2929$で中心を向き、$\lambda =0.5$の1回で半径は次の様に縮みます：

10回後の半径は${0.8536}^{10}\approx 0.205$で、ノイズが全く無かったのに元の5分の1です。

周波数で考えると理由が分かります。
$n$点の閉じたループでは、形$\cos \frac{2\pi jk}{n}$（点$k$、周波数$j$）はアンブレラ演算子で拡大縮小されるだけで、倍率は$-{\kappa }_j$です：

従って平滑化1回で周波数$j$は$1-\lambda {\kappa }_j$倍になります。
八角形の円は最も低い周波数$j=1$で、先程計算した様に${\kappa }_1=0.2929$です。
最もノイズの多い形、ジグザグ（1つおきに点を外へ押し出した物）は$\kappa =2$で、$1-2\lambda$倍になります。$\lambda =0.5$なら1回で消えます。
つまりラプラシアンスムージングはローパスフィルタですが、漏れの有るフィルタです。$j=0$（形全体の移動）以外の全ての周波数が弱められ、形もノイズと一緒に縮みます。

Taubinのλ|μ平滑化

Taubin（1995年）は2つ目のステップで此れを直しました。$\lambda >0$で縮める各ステップの後に、負ので膨らませるステップを入れます：

1組で周波数$\kappa$は次の倍率になります：

此れは$\kappa =0$で1になり、通過帯域周波数

で再び1になります。
${\kappa }_{PB}$より下では形は其の儘通り（実は$f$は1よりほんの少し大きい）、上では凸凹が弱められます。
Taubinは${\kappa }_{PB}$を0.01から0.1の間にする事を勧め、論文の例では$\lambda =0.6307$、${\kappa }_{PB}=0.1$を使っています。此れから$\mu =-0.6732$になります（$\frac{1}{0.6307}-\frac{1}{0.6732}=1.5855-1.4854=0.1001$）。

例えば$\kappa =1$では$(1-0.6307)(1+0.6732)=0.3693\cdot 1.6732=0.6179$です。
細かく標本化された円は残り、ジグザグは1組毎に約40%減ります（マイナスは毎回反転する事を意味します）。
八角形は粗過ぎて、其の「円」が既に高い周波数に入るので、Taubinでも少し縮みます。
±0.05のジグザグが付いた32点の円では、ラプラシアン20回でジグザグは消え、半径は${0.9904}^{20}\approx 0.824$に縮みます。Taubin 20組でもジグザグは消え、半径は保たれます（${1.00066}^{20}\approx 1.013$）。

ラプラシアン10回でノイズは消えますが、100回で曲線は3分の2に縮みます。Taubin 100組は大きさを保ちます（此処では約3%大きくなります。最も低い周波数で$f$が1より少し大きい為です）。
実際のツールは此れに工夫を加えます。頂点が法線方向にだけ動く様にするコタンジェント重み、大きく安定した1ステップの為に線型方程式を解く陰的フェアリング（Desbrunら、1999年）、そして鋭い辺を残す特徴保存フィルタ（第15章のバイラテラルフィルタのメッシュ版）です。
スカルプトソフトの「スムーズ」ブラシはブラシの範囲内のラプラシアンスムージングで、メッシュ処理ツールは生のスキャン用にTaubinの平滑化を備えています。

21.3 二次誤差尺度によるメッシュの簡略化

第9章でエッジ縮約を紹介しました。辺の両端を1つの頂点にまとめ、頂点1個、辺3本、三角形2個を減らす操作です。
其処では本当の問題が残っていました。次にどの辺を縮約し、新しい頂点を何処に置くかです。
GarlandとHeckbertの二次誤差尺度（QEM、1997年）は其の両方に安く上手く答え、今でも殆どの簡略化ツールの基礎です。

考え方

元のメッシュの各頂点は、其の周りの面の平面上に有ります。新しい頂点の良い置き場所は、其れらの平面にまだ近い場所です。
其処で位置の誤差を平面の集合までの距離の2乗の和で測ります。各頂点には其の面の平面を持たせ、2つの頂点がまとまる時は、新しい頂点が両方の集合を受け継ぎます。
「平面」が直線になる2次元では、頂点の周りの誤差は次の様になります：

鋭い角では誤差が全ての方向に急に増えるので、頂点は其の場に固定されます。
緩い曲がりでは曲線に沿ってゆっくり増えるので、頂点は滑る事が出来ます。
平らな部分では直線に沿って誤差が0なので、頂点を取り除いてもコストは掛かりません。
正に簡略化に求める動きです。平らな所から先に減らし、角と折り目は残します。

平面までの距離を行列の積で

平面を単位法線$(a,b,c)$で$ax+by+cz+d=0$と書き、其の数を$𝐪=(a,b,c,d)$にまとめます。
点を同次座標$𝐯=(x,y,z,1)$で表すと、符号付き距離は内積（線型代数シリーズの第4章）になり、其の2乗は行列の「サンドイッチ」になります：

$𝐊=𝐪{𝐪}^{\mathrm{\top }}$は対称な4 × 4行列で、其の平面の基本誤差クアドリックです：

平面が幾つも有る時、距離の2乗の和は

です。
つまり平面が幾つ有っても、1つの対称な4 × 4行列$𝐐$、頂点の二次誤差行列（クアドリック）にまとまります。数は全部で10個です。そして2つの平面の集合をまとめるのは、行列を2つ足すだけです：${𝐐}_{\text{new}}={𝐐}_1+{𝐐}_2$。

計算例：四角錐の頂点のクアドリック

頂点が$(0,1,0)$、底面の角が$(\pm 1,0,\pm 1)$の四角錐を考えます。
頂点と$(1,0,1)$、$(1,0,-1)$を通る面は平面$x+y=1$で、単位法線にすると$𝐪=\frac{1}{\sqrt{2}}(1,1,0,-1)$、

です。
他の側面は$-x+y=1$、$y+z=1$、$y-z=1$です。
4つを足すと、混ざった項は2つずつ打ち消し合い（$x$は1つの平面で$+1$、向かいの平面で$-1$、$z$も同じ）、次の様になります：

誤差は頂点で0になり、全ての方向に増えます（誤差一定の面は楕円体で、角の円の3次元版です）。
頂点を$(0,0.5,0)$まで下げるコストは$2\cdot {0.5}^2=0.5$です。
1つの平面で確かめましょう。$(0,0.5,0)$から$x+y=1$までの距離は$\frac{|0.5-1|}{\sqrt{2}}$、其の2乗は$0.125$で、平面4つで$0.5$になります。

新しい頂点の最適な位置

${𝐯}_1$から${𝐯}_2$への辺を縮約するには、クアドリックを足して$𝐐={𝐐}_1+{𝐐}_2$とし、誤差が最小になる位置$\overline{𝐯}$を探します。
$𝐐$を3 × 3のブロック$𝐀$、列$𝐛$、数$c$に分けると、3次元の点$𝐱$について

です。
此れはお椀形の2次関数なので、最小は勾配が0になる所です：

3 × 3の連立1次方程式です（線型代数シリーズの第10章）。代入し直すと最小の誤差${𝐛}^{\mathrm{\top }}\overline{𝐯}+c$が得られ、此れが辺のコストです。
$𝐀$が特異（又は其れに近い）なら、最適な点は1つに決まりません。例えば平らな領域では、平面上の全ての点の誤差が0です。
其の時の実装は、両端と中点の様な幾つかの候補から一番良い物を選びます。

計算例。 頂点から角$(1,0,1)$への辺を縮約します。
角は面$x+y=1$と$y+z=1$、そして底面$y=0$（単位法線$(0,-1,0)$、$d=0$）に接しています：

方程式$𝐀\overline{𝐯}=-𝐛$は次の通りです：

対称性から$x=z$なので、真ん中の式は$x+4y=3$になります。最初の式から$y=1-3x$、代入すると$x+4-12x=3$なので$x=\frac{1}{11}$、$y=\frac{8}{11}$です。
最適な位置は$\overline{𝐯}=(\frac{1}{11},\frac{8}{11},\frac{1}{11})\approx (0.091,0.727,0.091)$、誤差は${𝐛}^{\mathrm{\top }}\overline{𝐯}+c=-\frac{0.5+24+0.5}{11}+3=\frac{8}{11}\approx 0.727$です。
頂点に縮約するとコストは1、角に縮約すると4、中点$(0.5,0.5,0.5)$なら1.25です。
最適な位置は其の全てより安く、辺の上にさえ有りません。QEMは面を最も良く近似する所に頂点を置くのです。
（此の小さなメッシュで四角錐の頂点を縮約するのは悪い考えで、コストも其れを示しています。平らな領域の辺のコストは丁度0です。）

アルゴリズム

1. 各頂点のクアドリック${𝐐}_i$を其の面の平面から計算します（面積で重み付けする事が多いです）。
2. 各辺について、${𝐐}_1+{𝐐}_2$から最適な位置と其のコストを計算します。
3. 全ての辺をコスト順の優先度付きキュー（最小ヒープ）に入れます。
4. 一番安い辺を最適な位置に縮約し、新しい頂点に${𝐐}_1+{𝐐}_2$を持たせます。
5. 新しい頂点の周りの辺のコストを計算し直し、目標の三角形数になるまでステップ4から繰り返します。

1回の縮約は小さな行列計算が幾つかとヒープの更新だけなので、$n$面のメッシュ全体で$n\log n$程度の時間で済みます。
実際の実装は、縮約の前に確認をします。内部の辺の2つの頂点の共通の隣は丁度2つでなければなりません（さもないと結果が非多様体になります、第7章）。隣の三角形が裏返ってはいけません。開いた境界には、境界の各辺を通り其の面に垂直な重い平面を加えて、輪郭を残します。
UV、法線、色をクアドリックの次元として加えると、継ぎ目や色の境目も残せます（GarlandとHeckbert、1998年、Hoppe、1999年）。

此の地形では、平らな地面が先に数個の大きな三角形になり、小さな四角錐は鋭い辺を保ち、丘は少なく上手く配置された三角形で丸い形を保ちます。

コード：クアドリック

線型代数シリーズのlinalg.hhを使うと、クアドリックは10個の数を持ち、三角形から作られ、足され、評価され、クラメルの公式で最適な位置を求めます：

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

// A quadric Q = sum of q q^T for planes q = (a, b, c, d), stored as its 10 distinct entries.
struct Quadric {
  float a2 = 0, ab = 0, ac = 0, ad = 0, b2 = 0, bc = 0, bd = 0, c2 = 0, cd = 0, d2 = 0;

  // the plane through three points (unit normal n, d = -n.p)
  static Quadric fromTriangle(Vec3 p0, Vec3 p1, Vec3 p2) {
    Vec3 n = normalize(cross(p1 - p0, p2 - p0));
    float a = n.x, b = n.y, c = n.z, d = -dot(n, p0);
    return {a * a, a * b, a * c, a * d, b * b, b * c, b * d, c * c, c * d, d * d};
  }

  Quadric operator+(const Quadric &o) const {
    return {a2 + o.a2, ab + o.ab, ac + o.ac, ad + o.ad, b2 + o.b2,
            bc + o.bc, bd + o.bd, c2 + o.c2, cd + o.cd, d2 + o.d2};
  }

  // v^T Q v with v = (x, y, z, 1): the sum of squared distances to the planes
  float error(Vec3 p) const {
    float x = p.x, y = p.y, z = p.z;
    return a2 * x * x + 2 * ab * x * y + 2 * ac * x * z + 2 * ad * x
         + b2 * y * y + 2 * bc * y * z + 2 * bd * y
         + c2 * z * z + 2 * cd * z + d2;
  }

  // solve A v = -b by Cramer's rule; false if A is (nearly) singular
  bool optimum(Vec3 &v) const {
    float det = a2 * (b2 * c2 - bc * bc) - ab * (ab * c2 - bc * ac) + ac * (ab * bc - b2 * ac);
    if (std::fabs(det) < 1e-6f) return false;
    float rx = -ad, ry = -bd, rz = -cd;
    v.x = (rx * (b2 * c2 - bc * bc) - ab * (ry * c2 - bc * rz) + ac * (ry * bc - b2 * rz)) / det;
    v.y = (a2 * (ry * c2 - bc * rz) - rx * (ab * c2 - bc * ac) + ac * (ab * rz - ry * ac)) / det;
    v.z = (a2 * (b2 * rz - ry * bc) - ab * (ab * rz - ry * ac) + rx * (ab * bc - b2 * ac)) / det;
    return true;
  }
};

int main() {
  // the square pyramid of section 21.3: apex (0, 1, 0), base corners (+-1, 0, +-1)
  Vec3 apex = {0, 1, 0};
  Vec3 c0 = {1, 0, 1}, c1 = {1, 0, -1}, c2 = {-1, 0, -1}, c3 = {-1, 0, 1};
  Quadric qApex = Quadric::fromTriangle(apex, c0, c1) + Quadric::fromTriangle(apex, c1, c2)
                + Quadric::fromTriangle(apex, c2, c3) + Quadric::fromTriangle(apex, c3, c0);
  // the corner c0 touches two side faces and the base square
  Quadric qCorner = Quadric::fromTriangle(apex, c0, c1) + Quadric::fromTriangle(apex, c3, c0)
                  + Quadric::fromTriangle(c0, c3, c2);
  Quadric q = qApex + qCorner;

  printf("apex moved to (0, 0.5, 0): %.3f\n", qApex.error({0, 0.5f, 0}));
  printf("collapse to apex:     %.3f\n", q.error(apex));
  printf("collapse to corner:   %.3f\n", q.error(c0));
  printf("collapse to midpoint: %.3f\n", q.error(0.5f * (apex + c0)));
  Vec3 v;
  if (q.optimum(v)) printf("optimum (%.3f, %.3f, %.3f): %.3f\n", v.x, v.y, v.z, q.error(v));
  return 0;
}

apex moved to (0, 0.5, 0): 0.500
collapse to apex:     1.000
collapse to corner:   4.000
collapse to midpoint: 1.250
optimum (0.091, 0.727, 0.091): 0.727

全ての数が手計算と一致します。
法線の向きは関係有りません。$𝐪$と$-𝐪$は同じ$𝐪{𝐪}^{\mathrm{\top }}$になります。

21.4 プログレッシブメッシュとLODのストリーミング

Hoppeのプログレッシブメッシュ

簡略化を最後まで実行し、全ての縮約を記録します。
逆から読むと、各縮約は頂点分割です。1つの頂点が2つになり、三角形が2つ戻ります。
此れがHoppe（1996年）のプログレッシブメッシュで、粗いベースメッシュの後に頂点分割のリストが続きます：

各矢印が頂点分割1回で、$M^n$が元のメッシュです。
リストの途中までを使っても必ず正しいメッシュになるので、1つのファイルに全ての詳細度が入っています。
20,000三角形のモデルで、ベースメッシュが200三角形なら、分割は$\frac{20000-200}{2}=9900$回で、其々が頂点1個と三角形2個を加えます。
此れで、連続的なLOD（第9章の決まった数個の版ではなく、分割を幾つでも選べる）、ジオモーフ（新しい頂点が親から数フレーム掛けて滑り出し、ポッピングを隠す）、そして段階的な転送（先に$M^0$を送ってから分割を送るので、見る人はすぐに粗いモデルを見られる）が出来ます。
Hoppeの視点依存の詳細化（1997年）は、カメラに近い所、カメラを向いている所、シルエットになる所だけを分割します。

詳細度の選び方：スクリーン空間誤差

QEMのコストは各レベルの幾何学的な誤差を与えます（距離の2乗の和の平方根は距離です）。大事なのは、其れが画面上でどれ位の大きさに見えるかです。
第9章と同じく、垂直視野角60°、高さ1080画素の画面では、距離$d$に有る1 mは次の画素数になります：

計算例。 最大誤差1 cmのレベルを10 mから見ると、ずれは$0.01\times \frac{935.3}{10}\approx 0.94$画素で、見えません。
5 mからでは1.87画素なので、10 mから5 mの間の何処かで、もっと細かいレベルに切り替えるべきです。
エンジンは、物体又は領域毎に、スクリーン空間誤差が約1画素を超えない一番粗いレベルを選びます。

ストリーミングとクラスタLOD

GPUは頂点1つずつではなく大きなまとまりで処理したいので、現代のエンジンはクラスタからLODを作ります。小さな三角形のグループをQEMで簡略化して階層にまとめ、其々をスクリーン空間誤差で決めたレベルで描きます。
違うレベルの隣り合うクラスタはひび割れ無く繋がらなければならないので、クラスタの境界は簡略化の1段階では固定され、次の段階で解放されます。
執筆時点では、Unreal Engine 5のNaniteが最も有名な例です。其のドキュメントには、三角形の階層的なクラスタを描画中に詳細度の間で入れ替え、ひび割れ無く繋ぎ、必要に応じてストリーミングして見える詳細だけをメモリに置く、と説明されています。
粗い物を先に、必要な所だけ細かく、というプログレッシブメッシュの考えは其の儘で、単位がクラスタになっただけです。

21.5 リメッシュとメッシュ修復

等方的リメッシュ

スキャンやマーチングキューブ法の三角形は形が悪いです。細長い三角形、ばらばらな大きさ、3から12までの価数。
リメッシュは、其れらを（ほぼ）同じ面の上の新しい三角形で置き換えます。
等方的リメッシュは、1つの目標の辺の長さ$L$を持つ、ほぼ正三角形のメッシュを目指します。BotschとKobbelt（2004年）の広く使われる版は、4つのステップを10回程繰り返します：

1. $\frac{4}{3}L$より長い辺を全て中点で分割します。
2. $\frac{4}{5}L$より短い辺を全て縮約します。
3. 辺をフリップし（2つの三角形の対角線を入れ替え）、其の4つの頂点の価数が6（境界では4）に近付くなら其れを採用します。
4. 緩和：各頂点を隣の頂点の重心に向けて面に沿ってだけ（接線方向に）動かし、元の面の上に投影し直します。

何故$\frac{4}{3}$なのでしょうか？$L=1$なら、長さ$\frac{4}{3}$の辺は$\frac{1}{3}$長過ぎ、其の半分$\frac{2}{3}$は$\frac{1}{3}$短過ぎます。なので其れより長い辺だけを分割すれば必ず良くなります。
フリップでは、辺の2つの頂点の価数が7、向かいの2つの頂点が5なら、6からのずれは$1+1+1+1=4$です。フリップすると辺の頂点は隣を1つ失い、向かいの頂点は1つ得るので$6,6,6,6$、ずれは0になり、フリップします。
ステップ4は法線方向の成分を取り除いたアンブレラ演算子なので、面は縮みません。

壊れたメッシュの修復

スキャン、CADからの書き出し、古いゲームのメッシュは、第7章の規則を破っている事がよく有ります。
典型的なメッシュ修復の手順は次の通りです。許容値より近い頂点をまとめる（ポリゴンスープを繋がったメッシュに戻す）、面積0の三角形を取り除く、向きを揃える（共有する辺を渡って面を辿り、隣と巻き方が合わない面を裏返す）、非多様体の辺と頂点を別々の多様体の扇に分ける、穴を埋める（各境界ループを三角形で閉じ、周りの面に続く様に其の部分を平滑化する）、自己交差を取り除く（難しく、SDFを計算して新しい面を取り出す事が多いです、21.12節）。

穴を数える。 境界の辺はtwinの無いハーフエッジで、其れを辿ると境界ループが得られます。
数を数えるだけでも多くが分かります。種数$g$（取っ手の数）と境界ループ$b$本を持つ、繋がった向き付け可能な面では

です。
上の面が無い立方体は$V=8$、$E=12$、$F=5$なので、$8-12+5=1=2-0-b$、$b=1$で、穴は1つです。

裏返っていないか？ 面の向きは揃っているが、外向きか内向きか分からない時は、発散定理で符号付き体積を計算します：

各項は、原点から其の三角形までの四面体の符号付き体積の6倍です。
四面体$O=(0,0,0)$、$X=(1,0,0)$、$Y=(0,1,0)$、$Z=(0,0,1)$では、$O$を通る3つの面は0で、外向きに巻いた面$X,Y,Z$が$X\cdot (Y\times Z)=(1,0,0)\cdot (1,0,0)=1$を与えるので、$V=\frac{1}{6}$です。正なので法線は外を向いています。
体積が負なら、メッシュは裏返っていて、全ての面を裏返す必要が有ります。

21.6 メッシュのパラメータ化

パラメータ化とは

テクスチャマッピング（第11章）には、全ての頂点の$(u,v)$、つまり面から平面への写像が必要です。其れを計算するのがパラメータ化で、UV展開の裏に有るアルゴリズムです。
1枚で平らに出来るのは、位相的に円板（境界ループ1本、取っ手無し）の面だけです。なのでメッシュは先に継ぎ目（シーム）で其の様なチャートに切り分けられ、各チャートが平らにされ、チャートがテクスチャアトラスに詰め込まれます。

歪みを測る

曲がった面は伸ばさずに平らには出来ないので、目標は歪みを小さくする事です。
1つの三角形の上では、三角形（其の三角形自身の2次元の座標系）からUV空間への写像はアフィンで、一定の2 × 2のヤコビ行列$𝐉$を持ちます。
其の特異値${\sigma }_1\ge {\sigma }_2$は最大と最小の伸び率です。面上の小さな円は、半径の${\sigma }_1$倍と${\sigma }_2$倍を半軸とする楕円になります。
${\sigma }_1={\sigma }_2$なら写像は等角（角度が保たれ、大きさは変わっても良い）、${\sigma }_1{\sigma }_2=1$なら等積（面積が保たれる）、両方が1なら等長で、此れは平らに広げられる面（円柱、円錐）でしか出来ません。

計算例。 面上の直角三角形で、2本の脚の長さが2と1の物を、UVの三角形$(0,0)$、$(1,0)$、$(0,1)$に写します。
1本目の脚は$\frac{1}{2}$に縮み、2本目は長さが変わらないので、、${\sigma }_1=1$、${\sigma }_2=0.5$です。
$\frac{{\sigma }_1}{{\sigma }_2}=2$なので、市松模様は面上で正方形ではなく1 × 2の長方形に見えます。
${\sigma }_1{\sigma }_2=0.5$なので、此の三角形のテクセル密度（1平方メートル当たりのテクセル数）は歪みの無い三角形の半分です。

Tutteの埋め込み

Tutte（1963年）は、（3連結な平面）グラフの境界を凸多角形に固定し、全ての内部の頂点を隣の頂点の平均の位置に置くと、辺は交差せず、全ての面が凸になる事を証明しました。
円板の形のメッシュでは、境界ループを円又は正方形に固定し、次を解きます：

此れは$𝐋({𝐮}_i)=\mathrm{𝟎}$、つまりUV空間でアンブレラ演算子を0にする事です。
疎な連立1次方程式1つで済み、三角形が裏返る事は決して有りません。パラメータ化が持てる最も価値の有る保証です。
Floaterは、各内部頂点を隣の頂点の凸結合の位置に置く正の重みなら何でも此の保証を保つ事を示し（1997年）、彼の平均値座標の重み（2003年）は単純な平均よりずっと良く形に沿います。

計算例。 正方形のチャートの角が$C_0=(0,0)$、$C_1=(1,0)$、$C_2=(1,1)$、$C_3=(0,1)$、内部の頂点が$𝐮$、$𝐰$で、三角形は$(C_0,C_1,𝐮)$、$(C_1,𝐰,𝐮)$、$(C_1,C_2,𝐰)$、$(C_2,C_3,𝐰)$、$(C_3,𝐮,𝐰)$、$(C_3,C_0,𝐮)$です。
$𝐮$の隣は$C_0,C_1,C_3,𝐰$、$𝐰$の隣は$C_1,C_2,C_3,𝐮$です：

$𝐰=\frac{(2,2)+𝐮}{4}$を最初の式に代入すると$16𝐮=(4,4)+(2,2)+𝐮$なので、$𝐮=\frac{(6,6)}{15}=(0.4,0.4)$、$𝐰=(0.6,0.6)$です。正方形の内側で、裏返った三角形は有りません。
3次元の位置が一度も出て来なかった事に注目して下さい。単純な平均は繋がり方しか使わないので、三角形の大きさが大きく違うメッシュでは写像がひどく歪みます。

等角写像：LSCMとARAP

最小二乗等角写像（LSCM、Lévy・Petitjean・Ray・Maillot、2002年）は、直接${\sigma }_1={\sigma }_2$を目指します。
各三角形の上で、$v$方向が$u$方向を90°回して同じ長さにした物になっていれば、写像は等角です（コーシー・リーマンの方程式）。
全ての三角形で同時に其れが成り立つ事は無いので、LSCMは面積で重み付けした違反の2乗和を最小にします。エネルギーは2次なので、最小は此れも疎な連立1次方程式1つです。位置、回転、大きさを決める為に2つの頂点を固定し、境界は自由に自然な形を見付けます。
自動のテクスチャアトラスの為に作られた方法で、多くのツールの「UV展開」ボタンの裏には其の変種が有ります。
ARAPパラメータ化（as-rigid-as-possible、Liuら、2008年）は、各三角形の写像が回転に出来るだけ近い事（${\sigma }_1={\sigma }_2=1$）を求め、三角形毎の最適な回転（局所ステップ）と全てのUVの連立方程式1つ（大域ステップ）を交互に解きます。

色は$\frac{{\sigma }_1}{{\sigma }_2}$を、1（暗い青緑）から1.7（黄）、3以上（赤）で表しています。
Tutteは正方形の境界を円に押し込み、出っ張りを無視します。角度は殆ど全体で歪みます（平均1.74、4つの角で最悪）。
LSCMは角度をほぼ正しく保ち（平均1.10）、代わりに境界を曲げます。出っ張りの三角形はUV空間で小さくなり、其れが角度を保つ代償です。
UV展開ツールは其の後、チャートを数画素の余白を空けて詰め込み、ミップマップが滲み出さない様にします。

21.7 法線マップとディスプレイスメントマップのベイク

第11章では、法線マップがローポリゴンのモデルを細かく見せる事、そして法線マップは普通、高解像度のスカルプトからベイクされる事を見ました。其の仕組みです。

ベイクのループ

入力は、良いUVを持つローポリゴンのメッシュ（21.6節）と、同じ形のハイポリゴンのメッシュです。
出力のマップの各テクセルについて：

1. テクセルが覆うローポリゴンの面上の点$𝐩$と、補間した法線$𝐍$、接空間の基底$𝐓,𝐁$（第11章）を求めます。
2. 外側の始点から$-𝐍$の方向にレイを飛ばし、ハイポリゴンのメッシュとの交点を求めます。始点は普通ケージの上に取ります。ケージはローポリゴンのメッシュを外側に押し広げ、ハイポリゴンのメッシュを包む様にした複製です。
3. 交点でのハイポリゴンの法線$𝐧$を接空間に変換します：$(𝐧\cdot 𝐓,𝐧\cdot 𝐁,𝐧\cdot 𝐍)$。此れはTBN行列の転置を掛ける事です。
4. 各成分を$\frac{n+1}{2}$として符号化し、保存します。

計算例。 第11章の床の三角形では$𝐓=(1,0,0)$、$𝐁=(0,0,-1)$、$𝐍=(0,1,0)$です。
或るテクセルからのレイが、ワールド空間の法線が$𝐧=(0.6,0.8,0)$の出っ張りでスカルプトに当たります。
接空間では$(0.6,0,0.8)$、符号化すると$(0.8,0.5,0.9)$、8ビットでは$(204,128,230)$です。
第11章で復号したテクセルと全く同じで、復号すると$0.6𝐓+0.8𝐍=(0.6,0.8,0)$に戻ります。ベイクとシェーダーは互いに逆の操作なのです。

ディスプレイスメントマップは、法線の代わりに法線方向の符号付き距離を保存します。
ベイクが高さ−0.05から+0.05を0〜1に写し、スカルプトがローポリゴンの面より0.03上に有るなら、保存する値は$\frac{0.03+0.05}{0.1}=0.8$、8ビットでは204です。1次元の符号付き距離の標本で、此の章の残りの話題です。

上手く行かない所も有ります。ケージがきつ過ぎるとレイがスカルプトの内側から始まり、緩過ぎると違う部分（頬ではなく耳）に当たります。又、ベイクとレンダラーが$𝐓$と$𝐁$を違う方法で計算すると、ライティングに継ぎ目が出ます。
其の為、多くのツールがMikkTSpace（Mikkelsen）に揃えています。例えばglTF 2.0の仕様は、接線が無いファイルではクライアントがMikkTSpaceで接線を計算すべきだとしています。
同じループで、アンビエントオクルージョン（交点から半球に多くのレイを飛ばす）、曲率、厚み、スカルプトの色もベイク出来ます。

21.8 符号付き距離場

定義

第9章で符号付き距離関数を紹介しました。格子やテクスチャに保存した物は普通符号付き距離場と呼ばれ、どちらもSDFと略します。
面$S$を持つソリッド$\mathrm{\Omega }$について：

此のシリーズの何処でもそうである様に、内側が負です。
符号が無ければ符号無し距離場で、開いた面や曲線にも使えます。

性質

SDFを特別にする性質が3つ有ります：

• 勾配の長さが1：勾配が存在する所では$\Vert \mathrm{\nabla }f\Vert =1$です。面から真っ直ぐ1単位離れると、距離は丁度1単位増えます。此れがアイコナール方程式です。
• 勾配が法線：面の上では$\mathrm{\nabla }f$は外向きの単位法線で、何処でも、最も近い面上の点は$𝐩-f(𝐩)\mathrm{\nabla }f(𝐩)$です。
• 1-リプシッツ：$|f(𝐩)-f(𝐪)|\le \Vert 𝐩-𝐪\Vert$。値は動いた距離より速くは変わりません。此のリプシッツ連続性が、スフィアトレーシングを安全にします。

勾配は中心軸の上では存在しません。最も近い面上の点が2つ以上有る点の集まりで、球の中心や箱の内側の対角線などです。其処では場に尾根が有りますが、面からは離れています。

計算例。 単位球$f(𝐩)=\Vert 𝐩\Vert -1$の$𝐩=(3,4,0)$では、$f=5-1=4$、$\mathrm{\nabla }f=\frac{𝐩}{\Vert 𝐩\Vert }=(0.6,0.8,0)$で長さは1、最も近い面上の点は$(3,4,0)-4\cdot (0.6,0.8,0)=(0.6,0.8,0)$で、確かに球の上に有ります。
第9章の$x^2+y^2+z^2-9$も球を表しますが、距離ではありません（勾配$(2x,2y,2z)$が中心からの距離と共に大きくなる）。なので此れらの技は使えません。

正確なSDFと下界

スフィアトレーシングには正確な距離は要りません。
距離の下界（保守的なSDF）とは、符号が正しく、絶対値が本当の距離を決して超えない関数、$|f(𝐩)|\le \mathrm{dist}(𝐩,S)$の事です。
面の上で丁度0になる1-リプシッツな関数は全て此れです。値が1単位当たり高々1しか変わらないなら、距離$d$の点の値が$d$を超える事は出来ません。
シェーダーの「SDF」の殆どは下界で、其れで問題有りません。困るのは距離を大きく見積もる関数です。

解析的なSDF

単純な形には正確な式が有ります：

箱の式では、$𝐪$が点を1つの象限に折り畳み（箱は対称です）、各座標が箱からどれだけはみ出しているかを測ります。
外側では$𝐪$の正の部分が最も近い面、辺、角へのベクトルになり、内側では全ての成分が負で、其の最大値が最も近い面までの近さを表します。
$𝐛=(2,1,1)$では：

左の図（2次元、外側がオレンジ、内側が青、0.2単位毎に線）では、外側の帯は角で丸くなっています。其処では距離が点までの距離だからです。内側の帯には斜めの折り目、中心軸が有ります。

SDFの変換

形を動かしたり回したりするには、逆変換で点を変換します。剛体変換（第5章）は距離を保つので、$f({𝐌}^{-1}𝐩)$は正確です。
$s$倍の一様な拡大縮小は距離も$s$倍するので、$sf(𝐩/s)$です。
一様でない拡大縮小は距離を保ちません。$f(x/2,y,z)$はx方向に2倍長くなりますが、値はもう距離ではありません。最小の倍率を掛けた$min(s_x,s_y,s_z)f(\dots )$は再び下界になります。

21.9 SDFの組み合わせ

和・積・差

第9章のCSGの操作

は、常に内外を正しく与えます。
では結果は距離でしょうか？
和は両方の形の外側では正確です（最も近い面は本当に2つの内の近い方です）が、重なった内側では下界に過ぎません。積と差は内側で正確ですが、外側では下界に過ぎません。
全て1-リプシッツ（1-リプシッツな関数のminとmaxは1-リプシッツ）なので、全てスフィアトレーシングに安全です。

計算例。 2次元の半平面$x\le 0$と$y\le 0$の積（$f_A=x$、$f_B=y$）は、原点に角を持つ4分の1平面です。
$(1,1)$では$max(1,1)=1$ですが、角までの本当の距離は$\sqrt{2}\approx 1.414$です。
下界は小さ過ぎるので、角の近くで数ステップ無駄になりますが、見逃しは決して起きません。
上の中央の図にも同じ効果が見えます。2つの円の和は、場が片方の円からもう片方へ切り替わる所で、外側の帯に折り目が有ります。

スムーズミン

minは2つの形が出会う所に鋭い折り目を作ります。有機的なモデリングでは、第9章のメタボールの様な丸みが欲しくなります。
スムーズミンは$a=b$の近くで融合させます。よく使われる多項式版は、融合の半径$k$で：

$b-a\ge k$なら$h=1$で、結果は丁度$a$です。$a-b\ge k$なら丁度$b$です。
の所でだけ融合し、最も強いのは$a=b$の時で、$h=\frac{1}{2}$、補正は$\frac{k}{4}$です。

$k=0.4$の計算例：

• $a=b=0.2$：$h=0.5$、$\mathrm{smin}=0.2-0.4\cdot 0.25=0.1$。融合が面を此の点の方へ0.1引き寄せます。
• $a=0.7$、$b=0.1$：$h=\mathrm{clamp}(0.5-0.75)=0$、$\mathrm{smin}=0.1=min(a,b)$。継ぎ目から遠ければ何も変わりません。

結果はまだスフィアトレーシングに安全でしょうか？
$\frac{\mathrm{\partial }h}{\mathrm{\partial }a}=-\frac{1}{2k}$と、$h$の定義から出る$a-b=k(1-2h)$を使って微分します：

同様に$\frac{\mathrm{\partial }\mathrm{smin}}{\mathrm{\partial }b}=1-h$です。
従って$\mathrm{\nabla }\mathrm{smin}=h\mathrm{\nabla }a+(1-h)\mathrm{\nabla }b$は長さ高々1の2つのベクトルの重み付き平均で、長さは高々1です。2つのSDFの滑らかな和は1-リプシッツの儘で、正しい下界です。
指数版$-k\ln \left(e^{-a/k}+e^{-b/k}\right)$も同じ性質を持ちますが、場を何処でも少し変えます（$a=b$で$k\ln 2$）。
上の右の図は$k=0.5$の多項式版です。2つの円は滑らかな首で繋がり、外側の帯の折り目も消えています。

丸め、殻、繰り返し、変位

• 丸め：$f(𝐩)-r$は形を$r$だけ太らせ、凸の辺を丸めます。半分の大きさ$𝐛-r$の箱から$r$を引くと、半分の大きさ$𝐛$の角の丸い箱になります。
• 殻：$|f(𝐩)|-t$はソリッドを厚さ$2t$の中空の皮にします。
• 空間の繰り返し：$x\prime =\mathrm{mod}(x+\frac{c}{2},c)-\frac{c}{2}$は、1回分の計算で形を$c$単位毎に繰り返します。$c=1.5$なら、$x=4.0$は$\mathrm{mod}(4.75,1.5)-0.75=0.25-0.75=-0.5$になり、セルの中心4.5から0.5左です。形がセルの中に十分収まる時だけ正確で、そうでなければ隣のセルも調べる必要が有ります。
• 変位：$f(𝐩)+g(𝐩)$は凸凹を加えますが、$g$が1単位当たり1より速く変わると結果は1-リプシッツではなくなります。$g=0.3\sin (10x)$の傾きは最大$3$なので、和は1単位当たり最大$1+3=4$変わります。

21.10 スフィアトレーシング

アルゴリズムと、其れが安全な理由

スフィアトレーシングは、Hart（1996年）が陰関数曲面の描画の為に発表しました。今では多くの人がレイマーチングと呼びますが、其の名前は一定の歩幅で進む方法も含みます。
単位方向$𝐝$のレイ$𝐫(t)=𝐨+t𝐝$について：

を、（当たり）又は$t$が最大値を超える（外れ）まで繰り返します。
面を飛び越えられないのは、$f$が距離の下界だからです。今の点の周りの半径$f(𝐫(t_i))$の球の中には面が無いので、其の球の表面上に有る次の点は、高々面に接するだけです。
Hartは此れをアンバウンディングスフィア（空の球）と呼びました。
GLSLシリーズの第14章でスフィアトレーサーを一歩ずつ完成させます。此処では其の振る舞いを見ます。

左の図のステップ表：床$y=0$と、中心$(5.2,1.4)$、半径1の円なので$f=min(y,\Vert 𝐩-(5.2,1.4)\Vert -1)$です。レイは$(0,1.3)$から方向$(0.9989,-0.0459)$に出発し、で止まります。

最初の3ステップは床が最も近く、歩幅を決めています。其の後は円が代わり、2ステップで到達します（正確な交点は$t=4.249$）。

収束の速さ

正面から当たるレイは、平面に1ステップで届きます。平面までの距離が、レイに沿った距離と丁度同じだからです。
斜めだと遅くなります。平面$f=y$から高さ$h_0$に有り、角度$\theta$で下りる（$d_y=-\sin \theta$）レイは、1ステップで$f\sin \theta$だけ高さを失うので、高さは毎回$1-\sin \theta$倍になります：

$h_0=1$、$\epsilon =0.001$では：

例えば$\frac{\ln 0.001}{\ln 0.9}=\frac{-6.908}{-0.1054}=65.6$なので66ステップです。右の図は$\sin \theta =0.25$で、各円は前の円の0.75倍です。
床に沿ったレイやシルエットを掠めるレイがSDFの描画で一番高く付く部分で、其の為トレーサーはステップ数に上限を設け、「面の近くでステップが尽きた」場合を当たりとして扱います。
過緩和（Keinertら、2014年）は$\omega$を1.2〜1.6程度として$\omega f$だけ進み、新しい空の球が前の球と重ならなければ（歩幅が長過ぎた）、$\omega =1$に戻します。
逆に、$L>1$で$L$-リプシッツでしかない関数（上の変位した形では$L=4$）は、$\frac{f}{L}$ずつ進まなければなりません。

何時止めるか：画素の大きさ

1画素のレイは実は$t$と共に広がる細い円錐なので、一定の$\epsilon$は遠くでは厳し過ぎ、近くでは緩過ぎます。
垂直視野角60°、高さ1080画素の画面では、1画素は約$\frac{2\tan 30\mathrm{°}}{1080}=0.00107$ラジアンなので、$t=10$では円錐の幅は約$0.0107$です。
で止めると「1画素以内」という意味になり、遠くの面は早く止まり、近くの面は正確な儘です（Hartは此の円錐をアンチエイリアスの被覆率の推定にも使いました）。

勾配から法線を

当たった点の法線は正規化した勾配で、小さな$h$による中心差分で推定します：

計算例。 単位球の$𝐩=(0.6,0.8,0)$、$h=0.001$では、$f(0.601,0.8,0)=0.000600$、$f(0.599,0.8,0)=-0.000600$なので、xの差を$2h$で割ると$\frac{0.001200}{0.002}=0.600$です。yは0.800、zは0で、正確な法線になります。
正確なSDFなら正規化前の推定値の長さが既に1に近く、デバッグに便利な確認になります。
中心差分は6回の評価が要ります。四面体の4つの角${𝐤}_j$、$(1,-1,-1)$、$(-1,-1,1)$、$(-1,1,-1)$、$(1,1,1)$で$f$を評価し、${\sum }_j{𝐤}_{j}f(𝐩+h{𝐤}_j)$を使うと、4回で勾配に比例するベクトルが得られます。

ソフトシャドウとアンビエントオクルージョン

SDFはレイが全ての物にどれだけ近くを通るかを知っているので、安くて有名な効果が2つ得られます（どちらもGLSLシリーズの第14章で実装しています）：

• ソフトシャドウ：光源に向かって進みます。各ステップで$\frac{h}{t}$（最も近い物までの距離÷進んだ距離）は、レイが何かを外れる角度の大まかな値で、${min}_i\frac{k{h}_i}{t_i}$を1で抑えた物が、光源の見える割合の推定になります。$k=8$で、$t=2$進んだ所で物から$h=0.1$の所を通ると$\frac{8\cdot 0.1}{2}=0.4$、つまり40%照らされます。$k$が大きい程光源は小さく、影は硬くなります。
• アンビエントオクルージョン：法線に沿った幾つかの距離$\delta$で$f$を標本化します。開けた面なら$f(𝐩+\delta 𝐧)=\delta$で、足りない分$\delta -f$は、近くの何かが空の一部を遮っている事を意味します。

どちらも安い近似です。リアルタイムの影とAO一般は第24章で扱います。

シェーダーのスフィアトレーサー

SDFレンダラーの核です。球、滑らかな和、シーン、そしてマーチングです。

float sdSphere(float3 p, float r) { return length(p) - r; }

float smin(float a, float b, float k) {
  float h = saturate(0.5 + 0.5 * (b - a) / k);
  return lerp(b, a, h) - k * h * (1.0 - h);
}

// the scene: two spheres blended together, standing on the floor y = 0
float map(float3 p) {
  float blob = smin(sdSphere(p - float3(-0.8, 0.5, 0.0), 0.5),
                    sdSphere(p - float3(-0.2, 0.9, 0.1), 0.4), 0.3);
  return min(p.y, blob);
}

// sphere tracing; returns the distance along the ray, or -1.0 for a miss
float trace(float3 ro, float3 rd) {
  float t = 0.0;
  for (int i = 0; i < 128; i++) {
    float h = map(ro + t * rd);
    if (h < 0.001 * t) return t; // closer than about one pixel footprint
    t += h;
    if (t > 50.0) break;
  }
  return -1.0;
}

float sdSphere(vec3 p, float r) { return length(p) - r; }

float smin(float a, float b, float k) {
  float h = clamp(0.5 + 0.5 * (b - a) / k, 0.0, 1.0);
  return mix(b, a, h) - k * h * (1.0 - h);
}

// the scene: two spheres blended together, standing on the floor y = 0
float map(vec3 p) {
  float blob = smin(sdSphere(p - vec3(-0.8, 0.5, 0.0), 0.5),
                    sdSphere(p - vec3(-0.2, 0.9, 0.1), 0.4), 0.3);
  return min(p.y, blob);
}

// sphere tracing; returns the distance along the ray, or -1.0 for a miss
float trace(vec3 ro, vec3 rd) {
  float t = 0.0;
  for (int i = 0; i < 128; i++) {
    float h = map(ro + t * rd);
    if (h < 0.001 * t) return t; // closer than about one pixel footprint
    t += h;
    if (t > 50.0) break;
  }
  return -1.0;
}

float sdSphere(float3 p, float r) { return length(p) - r; }

float smin(float a, float b, float k) {
  float h = saturate(0.5 + 0.5 * (b - a) / k);
  return lerp(b, a, h) - k * h * (1.0 - h);
}

// the scene: two spheres blended together, standing on the floor y = 0
float map(float3 p) {
  float blob = smin(sdSphere(p - float3(-0.8, 0.5, 0.0), 0.5),
                    sdSphere(p - float3(-0.2, 0.9, 0.1), 0.4), 0.3);
  return min(p.y, blob);
}

// sphere tracing; returns the distance along the ray, or -1.0 for a miss
float trace(float3 ro, float3 rd) {
  float t = 0.0;
  for (int i = 0; i < 128; i++) {
    float h = map(ro + t * rd);
    if (h < 0.001 * t) return t; // closer than about one pixel footprint
    t += h;
    if (t > 50.0) break;
  }
  return -1.0;
}

下の図は、此の章の技法でCPU上で描画しました。球と滑らかに融合した角の丸い箱、トーラス、空間の繰り返しによる果てしない柱の列、市松模様の床、$k=8$のソフトシャドウ、5標本のアンビエントオクルージョン、四面体による法線、そして少しの霧です。

三角形は1つも使っていません。

21.11 メッシュからのSDFとSDFの格子

三角形メッシュまでの距離

普通のモデルにSDFの技法を使うには、メッシュのSDFが必要です。符号無しの距離は、全ての三角形について点から三角形までの距離を取った最小値です（三角形の平面に投影し、最も近い点が内側、辺の上、角の何れかを調べます）。遠い三角形はバウンディングボリューム階層（第12章）で飛ばします。
難しいのは符号です：

• レイの偶奇：点から出したレイが面を横切る回数を数え、奇数なら内側です。辺や頂点を掠めるレイは数え間違えるので、実装は複数のレイで多数決を取ります。
• 角度で重み付けした擬法線（Bærentzen・Aanæs、2005年）：最も近い面上の点から問い合わせの点へのベクトルが、其処の法線と逆を向いていれば内側です。最も近い点が辺や頂点の上に有る時は、面法線を角度で重み付けした平均（第9章）を使うと、閉じたメッシュで判定が正しくなります。
• 一般化回転数（Jacobson・Kavan・Sorkine-Hornung、2013年）：点から見た全ての三角形の符号付き立体角を足し、$4\pi$で割ります。閉じたメッシュでは内側で1、外側で0、穴の近くでは中途半端な値になり、0.5で区切っても妥当な答えが得られます。

計算例。 立方体の中心からは、各面が方向の球全体の$\frac{4\pi }{6}$を覆うので、回転数は$\frac{6\cdot 4\pi /6}{4\pi }=1$で、内側です。
面を1つ消すと、レイの偶奇はレイの向きに依って変わりますが、回転数は$\frac{5}{6}\approx 0.83>0.5$で、人間が言う通り内側の儘です。

SDFを格子に保存する

一般の形のSDFには式が無いので、3次元の格子で標本化し、3線型補間（トライリニア補間）（第11章の双線型補間の3次元版）で読み出します。
密な格子は大きくなります。${256}^3$個の標本を2バイトで持つと$33554432$バイト、1物体で32 MiBです。
しかし殆どの場合、面の近くの値だけが大事なので、エンジンはナローバンド（$|f|$が小さい標本だけ）や疎な格子（面の近くにだけ確保し、木やハッシュで探すボクセルの小さなブロック）を使います。
映画のエフェクトで広く使われるOpenVDB（Museth、2013年）は、そんな階層的な疎な格子です。

計算例。 ${256}^3$の格子の中の半径100ボクセルの球で、帯の幅が±3ボクセルなら、帯は半径97と103の間の殻で、$\frac{4}{3}\pi ({103}^3-{97}^3)\approx 754208$ボクセル、格子の$16777216$の約4.5%です。
面は2次元なので、帯は解像度の3乗ではなく2乗で増えます。
格子のSDFは、上の最近接三角形の探索か、面から外へ距離を広げていく距離変換（ファストマーチング法、ファストスイーピング法、GPUではジャンプフラッディング）で作ります。

実際には

執筆時点では、Unreal Engineがよく知られた例です。ドキュメントによれば、Lumenのソフトウェアレイトレーシングはレイの最初の2 mでは各メッシュ自身の距離場を、残りでは其れらをまとめたグローバル距離場をトレースし、エンジンはメッシュの距離場をソフトシャドウとアンビエントオクルージョンにも使います。
物理エンジンはSDFを衝突判定に使い（めり込みの深さと押し出す方向は$f$と$\mathrm{\nabla }f$其の物です、第26章参照）、映画のツールは液体、煙、破壊に疎なレベルセットを使います。

フォントとデカール：SDFテクスチャ

アルファマスクとして保存した文字やデカールは、拡大するとぼやけたり、でこぼこになったりします。被覆率をフィルタリングしても、縁が何処に有るかは分からないからです。
Green（Valve、2007年）は、代わりに距離場を保存する事を提案しました。高解像度の2値画像から距離場を計算し（彼の典型的な例は4096 × 4096から64 × 64）、小さなテクスチャの1チャンネルに保存し、0.5での単純なアルファテストで描きます。
距離の双線型補間は縁の近くでほぼ線型なので、0.5の等値線は大きく拡大しても鮮明な儘です。ValveはTeam Fortress 2で、例えばデカールに此れを使いました。
8ビットのチャンネルは範囲が限られるので、スプレッド$s$（テクセル単位）で0〜1に写します。此のシリーズの符号の約束で、内側が0.5より上になる様にすると：

計算例。 $s=4$テクセルでは、縁の上で$v=0.5$（8ビットで128）、2テクセル内側で$v=0.5+\frac{2}{8}=0.75$（191）、4テクセル以上外側で$v=0$です。

拡大した2つの文字はどちらも32 × 32のテクスチャから作りました。アルファのビットマップは縁が波打ち、距離場は真っ直ぐな儘です。
エフェクトは別のしきい値に過ぎません。縁取りは2つの値の間の帯、グローは更に外側の緩やかな傾斜です。
シェーダーでは、縁を1画面画素に渡って滑らかにしてアンチエイリアスします。1画素がテクスチャの値で幾らに当たるかは、微分fwidthが与えます：

// v: the distance texture's value at this pixel (0.5 on the edge, larger inside)
float sdfCoverage(float v) {
  float w = 0.5 * fwidth(v); // half a screen pixel, in texture value units
  return smoothstep(0.5 - w, 0.5 + w, v);
}

// v: the distance texture's value at this pixel (0.5 on the edge, larger inside)
float sdfCoverage(float v) {
  float w = 0.5 * fwidth(v); // half a screen pixel, in texture value units
  return smoothstep(0.5 - w, 0.5 + w, v);
}

// v: the distance texture's value at this pixel (0.5 on the edge, larger inside)
float sdfCoverage(float v) {
  float w = 0.5 * fwidth(v); // half a screen pixel, in texture value units
  return smoothstep(0.5 - w, 0.5 + w, v);
}

弱点は「A」の頂点に見えます。鋭い角が少し丸くなるのです。1つの距離では、角度を成して出会う2本の縁を表せないからです。
マルチチャンネルSDF（MSDF、Chlumskýら、msdfgenツール）は3つの距離をRGBに保存し、其の中央値を取る事で角を鋭く保ちます。
文字の描画はOpenGLシリーズの第26章で実装します。其処では、エンジンがどんな大きさでも文字を鮮明に保つ方法としてSDFアトラスに触れています。

21.12 SDFからメッシュへ

スカルプトを3Dプリントしたり、シミュレーションしたり、ラスタライズしたりするには、SDFを再びメッシュにする必要が有ります。
第9章のマーチングキューブ法は、符号が変わる格子の辺毎に頂点を置き、256通りの表で其れらを繋ぎます。全ての頂点が格子の辺の上に有るので、箱の角の様にセルの中に有る鋭い特徴は切り落とされます。
2つの双対的な方法は、代わりに頂点をセルの中に置きます：

• サーフェスネット（Gibson、1998年）：面が通る各セルに頂点を1つ、セルの辺との交点の平均の位置に置きます。符号が変わる格子の辺毎に、其の周りのセルの頂点を繋ぎます（3次元では4つのセルで四角形、2次元では2つのセルで線分）。表は要らず、面の形も良くなりますが、角は丸い儘です。
• デュアルコンタリング（Ju・Losasso・Schaefer・Warren、2002年）：繋ぎ方は同じですが、各頂点は交点での法線も使います（$\mathrm{\nabla }f$から。交点と法線の組をエルミートデータと呼びます）。法線${𝐧}_i$を持つ交点${𝐩}_i$は其々平面を決め、頂点は${\sum }_i{\left({𝐧}_i\cdot (𝐱-{𝐩}_i)\right)}^2$を最小にする点です。

此れは平面までの距離の2乗の和、正に21.3節のクアドリックです。論文は此れをQEF（二次誤差関数）と呼び、$𝐀={\sum }_i{𝐧}_i{𝐧}_i^{\mathrm{\top }}$として$𝐀𝐱=-𝐛$を同じ様に解きます。

計算例（2次元）。 セル$[0,1]\times [0,1]$の中に、$x=0.7$と$y=0.6$で囲まれた形の角が有ります（内側はかつ）。
下の辺は$(0.7,0)$で法線$(1,0)$、左の辺は$(0,0.6)$で法線$(0,1)$で横切られます。
マーチングスクエア法は交点を$(0.7,0)$から$(0,0.6)$への線分で結び、角を切り落とします。サーフェスネットは頂点を平均$(0.35,0.3)$に置きます。デュアルコンタリングは$(x-0.7{)}^2+(y-0.6{)}^2$を最小にし、其の最小点$(0.7,0.6)$は角其の物です。

マーチングスクエア法は正方形を八角形に、サーフェスネットは角の丸い正方形にし、デュアルコンタリングは4つの角を全て取り戻します。
代償として、QEFの頂点がセルの外に出たり、面が自己交差したりする事が有るので、実装は頂点をセル内に制限し、結果を多様体に保つ工夫をします。
ボクセルの地形やSDFのスカルプトツールは、よく此れらの双対的な方法を使います。

検定のポイント

此の章はシラバスを超えますが、次の第1部の点は確実にしておきたい所です：

• ポリゴン曲面（第9章）：エッジ縮約は頂点1個、辺3本、三角形2個を減らします。其の逆が頂点分割で、縮約の記録がプログレッシブメッシュになります。ポッピングはLODのレベル間の目に見える跳びです。
• オイラーの公式（第7章）：穴の無い閉じた面では$V-E+F=2$なので、角度欠損の合計は$4\pi$です。
• 陰関数曲面（第9章）：CSGはminとmaxです。SDFの値は面までの距離です。マーチングキューブ法は線型補間で格子の辺の上に頂点を置きます。
• マッピング（第11章）：UV展開は継ぎ目で面を切り開きます。法線マップはハイポリゴンのモデルからベイクし、接空間で保存します。

21.13 まとめ

• メッシュにはラプラシアン（アンブレラ演算子${𝐜}_i-{𝐯}_i$、又はコタンジェント重み）とガウス曲率（角度欠損$2\pi -\sum \theta$、合計は$2\pi (V-E+F)$）が有ります。
• ラプラシアンスムージングは形を縮める漏れの有るローパスフィルタです。Taubinの$\lambda |\mu$のステップは${\kappa }_{PB}=\frac{1}{\lambda }+\frac{1}{\mu }$より下の周波数を通します。
• QEMは頂点の平面のクアドリック$𝐪{𝐪}^{\mathrm{\top }}$を1つの4 × 4行列に足し合わせます。辺のコストは$𝐀\overline{𝐯}=-𝐛$で求まる${𝐯}^{\mathrm{\top }}({𝐐}_1+{𝐐}_2)𝐯$の最小値です。安い物から縮約するので、平らな所が消え、角は残ります。
• プログレッシブメッシュは連続的なLOD、ジオモーフ、ストリーミングを可能にします。エンジンはスクリーン空間誤差でクラスタを詳細化します。
• リメッシュは分割、縮約、フリップ、緩和を繰り返します。修復は頂点をまとめ、向きを揃え、穴を埋めます。オイラーの数え方と符号付き体積が助けになります。
• パラメータ化は円板状のチャートを平らにし、歪みはヤコビ行列の特異値で測ります。Tutteは三角形を決して裏返さず、LSCMは角度を保ち、ARAPは剛体的な三角形を目指します。
• ベイクはケージからハイポリゴンのメッシュへレイを飛ばします。ベイクとレンダラーは接空間の基底で一致しなければなりません。
• SDFの勾配は単位法線で、場は1-リプシッツです。符号が正しい1-リプシッツな下界ならスフィアトレーシングに安全で、min、max、多項式のスムーズミンは其れを保ちます。
• スフィアトレーシングは面を見逃しませんが、掠める角度では遅くなります。法線、ソフトシャドウ、AOは殆どただで得られます。
• メッシュのSDFには頑丈な符号の判定（擬法線、回転数）とナローバンドや疎な格子が必要です。SDFテクスチャは文字を鮮明に保ちます。
• サーフェスネットとデュアルコンタリングは頂点をセルの中に置きます。デュアルコンタリングのQEFは再びQEMのクアドリックで、角を保ちます。

21.14 用語

21.15 課題

1. ★ 正二十面体の頂点（12個の頂点の其々に正三角形が5つ集まる）での角度欠損と、全ての頂点の合計を計算し、$2\pi (V-E+F)$と比べて下さい。次に正三角形が7つ集まる頂点の角度欠損を計算して下さい。其れはどんな点ですか？
2. ★ 半分の大きさ$𝐛=(1,2,0.5)$の箱について、$(3,0,0)$、$(2,3,0)$、$(2,3,1.5)$、$(0,1,0)$でのSDFを計算し、其々で最も近い物（面、辺、角）を答えて下さい。次に、丸めの半径0.2で外側の大きさが同じ角の丸い箱のSDFを$(3,0,0)$で計算して下さい。
3. ★★ 或る頂点が平面$x=0$、$y=0$、$z=0$（其々面1つ）を持ちます。其のクアドリック${𝐐}_1$を書いて下さい。もう1つの頂点は平面$x=1$だけを持ちます。$𝐐={𝐐}_1+{𝐐}_2$を作り、$𝐀\overline{𝐯}=-𝐛$を解き、縮約のコストを求め、結果を言葉で説明して下さい。
4. ★★ $\lambda =0.5$、$\mu =-0.53$のTaubinの平滑化について、${\kappa }_{PB}$と、$\kappa =0.05$、$0.5$、$2$での1組の倍率$f(\kappa )$を計算して下さい。其々の周波数は残りますか、弱められますか、消えますか？
5. ★★★ 指数のスムーズミン$-k\ln (e^{-a/k}+e^{-b/k})$が、$a$と$b$が1-リプシッツなら1-リプシッツである事、そして$min(a,b)-k\ln 2$と$min(a,b)$の間に有る事を示して下さい。次に、$f(𝐩)+0.2\sin (8p_x)\sin (8p_y)\sin (8p_z)$ではスフィアトレーサーが何故歩幅を小さくしなければならないか、最大で何分の1にすれば良いかを説明して下さい。

選択問題

問1 ★ 次のうち、二次誤差尺度（QEM）によるメッシュの簡略化の説明として正しい物はどれか。

問2 ★ 次のうち、符号付き距離場（SDF）とスフィアトレーシングについて正しい物はどれか。

問3 ★★ 多項式のスムーズミンsmin(a, b) = (1 − h) b + h a − k h (1 − h)、h = clamp(0.5 + (b − a) / (2k), 0, 1)を使う時、k = 0.4でのsmin(0.3, 0.1)はどれか。

21.16 実践課題

方眼紙に小さな2次元のSDFのシーンを作り、コンパスでスフィアトレーシングし、其の輪郭をメッシュとして平滑化と簡略化をします。
方眼紙（少なくとも50 × 30マス）、定規、コンパス、鉛筆、2色のペン、電卓を用意して下さい。
1単位 = 4マスとし、$x$は右へ−5から7、$y$は上へ−2から4とします。

シーンは、中心$(0,0)$、半分の大きさ$(2,1)$の箱と、中心$(4,1)$、半径1の円です：

1. 場。 2つの形を描きます。$(0,0)$、$(0,2)$、$(3,0)$、$(3,3)$、$(4,1)$、$(-3,-2)$での$f$を計算して各点の横に書き、内側の点に印を付けます。$(3,0)$と$(-3,-2)$からは面までの最短の線分を描き、長さを測ります。
2. 3本のレイ。 レイAを$(-4,4)$から方向$(0.8,-0.6)$へ、レイBを$(7,4)$から方向$(-0.6,-0.8)$へスフィアトレーシングします。各点で$f_B$、$f_C$、$f$を計算し、コンパスで空の円を描き、レイに沿って円の縁まで進みます。で止めます。$t$、点、距離の表を作ります（小数3桁）。次にレイCを$(-5,2.5)$から方向$(1,0)$へ、円のすぐ上を通るので$t>12$までトレースし、歩幅が最も短くなる所に印を付けます。
3. 法線。 レイBが当たった点で、$h=0.1$の中心差分で法線を推定して正規化し、円の正確な法線$\frac{𝐩-(4,1)}{\Vert 𝐩-(4,1)\Vert }$と比べます。
4. 滑らかな和。 $k=1$で、$(2.5,1)$と$(3,1.8)$での$\mathrm{smin}(f_B,f_C)$を計算し、$min$と比べます。2色目のペンで、滑らかな和が2つの形の間の輪郭をどう変えるかをスケッチします。
5. 輪郭の平滑化。 箱の輪郭を、ノイズの有る8頂点の閉じたメッシュとして扱います。順番に$(2.2,0)$、$(2,1)$、$(0,1.3)$、$(-2,1)$、$(-1.8,0)$、$(-2,-1)$、$(0,-0.8)$、$(2,-1)$です。靴紐の公式$A=\frac{1}{2}|{\sum }_i(x_i{y}_{i+1}-x_{i+1}{y}_i)|$で面積を計算します。$\lambda =0.5$のアンブレラステップを1回（各頂点を2つの隣の中点の方へ）掛けて描き、面積を計算します。次に其の結果に$\mu =-0.53$のステップを1回掛けて描き、もう一度面積を計算します。
6. 輪郭の簡略化。 箱の角と辺の中点に8頂点を置いたノイズの無い輪郭を使います。直線の2次元クアドリック（$𝐪=(a,b,d)$、$𝐊=𝐪{𝐪}^{\mathrm{\top }}$）で、中点$(0,1)$を角$(2,1)$に縮約するコストと、角を中点に縮約するコストを計算します。QEMはどちらを選びますか？4つの中点を全て縮約すると何が残りますか？

自分で確かめる方法：各ステップは前の円の上に乗り、どの円も輪郭を横切りません（接するのは構いません）。レイAとBは箱の上面と円の右上の近くでになって終わり、レイCは抜けていき、円のすぐ上で歩幅が最も短くなります。法線の長さは1に近く、正確な法線と小数2桁程度で一致します。どちらの点でも$\mathrm{smin}\le min$で、2つの距離の差が$k$以上の所では等しくなります。$\lambda$のステップは面積を減らし、$\mu$のステップは元に近い面積に戻します。中点の縮約のコストは0で、角の縮約は0ではないので、QEMは4つの角を残します。