ベジェ境界の面光源を境界積分法で解く

0対象と記号

光源の形が曲線で囲まれていても、点 P の照度は光源を細かく分割せずに、境界に沿った1次元の積分だけで求められます。ここでは次の設定で話を進めます。

• 光源:平面上の領域 $S$。境界 $C=\partial S$ は閉じたベジェ曲線($d$ 次を数本。ツールは3次を4本)。輝度 $L$ が一様な完全拡散面で、片面だけが光る。
• 測る点:平面 F 上の点 P。F の単位法線を $\mathbf n$ とする。
• 求める量:P の照度 $E$。輝度 $L$ は光源の性質で、点ごとに変わるのは照度。$E=L\,G$ と書き、$G$(単位は sr)は形状と位置だけで決まる幾何係数。形態係数は $G/\pi$。

1全体の流れ

照度の定義である面積分(2重積分)から出発し、5つのステップで「曲線ごとの1次元積分の和」にします。各ステップの中身は第2節で導きます。

1. 光源全体にわたる面積分
   $$E=L\int_S\frac{\cos\theta_P\cos\theta_Q}{r^2}\,dA$$

   光源上の点 Q すべてについて足し合わせる。2重積分で、分割数に誤差が依存する。

2. 立体角に直す
   $$E=L\int_\Omega\cos\theta\,d\omega=L\,\mathbf n\cdot\int_\Omega\boldsymbol\omega\,d\omega$$

   P から見た光源の広がり $\Omega$(単位球面上の領域)だけで決まる。

3. ストークスの定理で境界の線積分に
   $$\int_\Omega\boldsymbol\omega\,d\omega=\frac12\oint_{\partial\Omega}\boldsymbol\omega\times d\boldsymbol\omega$$

   領域の積分が、領域の縁に沿った積分になる。

4. 縁を曲線 $\mathbf B(t)$ で表す
   $$\mathbf v(t)=\mathbf B(t)-\mathbf P,\qquad \boldsymbol\omega\times d\boldsymbol\omega=\frac{\mathbf v\times\mathbf v'}{|\mathbf v|^{2}}\,dt$$

   P から曲線上の点へのベクトル $\mathbf v$ だけで書ける。

5. 結果:曲線ごとの1次元積分の和
   $$E=\frac L2\oint_C\frac{\mathbf n\cdot\big(\mathbf v\times\mathbf v'\big)}{|\mathbf v|^{2}}\,dt$$

   境界の各曲線について $t$ で積分して足す。

2導出を順に追う

2.1 面積分から立体角へ(ステップ 1→2)

光源上の微小面積 $dA$ を P から見た立体角は $d\omega=\cos\theta_Q\,dA/r^2$ です。したがって $\cos\theta_P\cos\theta_Q/r^2\,dA=\cos\theta_P\,d\omega$ となります。光源側の $\cos\theta_Q$ と $1/r^2$ は、面積を立体角に換算する係数そのものだったわけです。輝度 $L$ は方向によらず一様なので、$L$ を積分の外に出せます。

さらに $\cos\theta_P=\mathbf n\cdot\boldsymbol\omega$ なので、照度は「光源の方向 $\boldsymbol\omega$ を、光源が張る立体角の上で足し合わせたベクトルに、$\mathbf n$ を掛けたもの」になります。なお $\mathbf n\cdot\boldsymbol\omega<0$ の方向(F の水平線より下)は照度に寄与しないので、$\Omega$ は水平線より上の部分に限ります。これが第7節のクリッピングです。

2.2 領域の積分を縁の積分に(ステップ 2→3)

$\displaystyle\int_\Omega\boldsymbol\omega\,d\omega$ を縁の積分にするのが、ストークスの定理です。次のように示せます。

定ベクトル $\mathbf a$ を1つ取り、$\mathbf A(\mathbf x)=\tfrac12\,\mathbf a\times\mathbf x$ とおきます。$\nabla\times\mathbf A=\tfrac12\big(3\mathbf a-\mathbf a\big)=\mathbf a$ です。単位球面(外向き法線 $\boldsymbol\omega$)でストークスの定理を使うと

左辺はスカラー三重積の並べ替えで $\mathbf a\cdot\tfrac12\oint\boldsymbol\omega\times d\boldsymbol\omega$ に等しく、$\mathbf a$ は任意なので、ステップ 3 の式が得られます。

検算:$\Omega$ が上半球なら、縁は赤道 $\boldsymbol\omega=(\cos\varphi,\sin\varphi,0)$ で、$\boldsymbol\omega\times d\boldsymbol\omega=(0,0,1)\,d\varphi$。右辺は $\tfrac12(0,0,2\pi)=(0,0,\pi)$。左辺を直接計算した $\int\cos\theta\,d\omega=\pi$ と一致します。

2.3 曲線で書く(ステップ 3→5)

縁を、パラメータ $t$ のベジェ曲線 $\mathbf B(t)$ の像として表します。P から曲線上の点へのベクトルを $\mathbf v(t)=\mathbf B(t)-\mathbf P$、$r=|\mathbf v|$、$\boldsymbol\omega=\mathbf v/r$ とすると

($\mathbf v\times\mathbf v=\mathbf 0$ なので第2項が消えます。)これをステップ 3 に代入して $\mathbf n$ との内積を取ると、ステップ 5 の式になります。

幾何の意味:$\boldsymbol\omega\times d\boldsymbol\omega=\boldsymbol\gamma\,d\varphi$ と書けます。$d\varphi=|\mathbf v\times\mathbf v'|/|\mathbf v|^2\,dt$ は、P から見た方向が $dt$ の間に回る角、$\boldsymbol\gamma$ はその回転面(P と曲線の接線を含む平面)の単位法線です。つまり

は「P から見た境界の見込み角を、回転面の傾き $\mathbf n\cdot\boldsymbol\gamma$ で重みづけして足したもの」です。回転面が水平($\boldsymbol\gamma=\pm\mathbf n$)なら見込み角がそのまま効き、回転面が $\mathbf n$ に平行なら寄与は 0 です。

3直線辺の場合:Lambert の公式($d=1$)

境界が線分 $\mathbf a\to\mathbf b$($\mathbf a,\mathbf b$ は P から両端へのベクトル)なら、$\mathbf v(t)=\mathbf a+t(\mathbf b-\mathbf a)$ で、次の2つが成り立ちます。

• $\mathbf v\times\mathbf v'=\mathbf a\times(\mathbf b-\mathbf a)=\mathbf a\times\mathbf b$ は $t$ によらない定数。回転面 $\boldsymbol\gamma=\mathbf a\times\mathbf b/|\mathbf a\times\mathbf b|$ が辺の上で一定。
• $|\mathbf v|^2$ は $t$ の2次式。

したがって $\displaystyle\int_0^1\frac{|\mathbf a\times\mathbf b|}{|\mathbf v|^2}\,dt$ は辺全体を見込む角 $\beta$($\mathbf a$ と $\mathbf b$ のなす角)そのものです。辺の寄与は $\beta\,(\mathbf n\cdot\boldsymbol\gamma)$ となり、多角形の $E$ は

(Lambert の公式。頂点の並びの向きは第2.2節の向きに合わせます。)

この式は平方根・逆三角関数だけで書けて、近似誤差がありません。多角形では、これが最後の答えです。

4ベジェ曲線の場合と次数

4.1 被積分関数の次数と、ベジェ形式での表し方

$d$ 次のベジェ曲線 $\mathbf B(t)=\sum_{i=0}^{d}\mathbf P_i\,B_i^d(t)$($B_i^d(t)=\binom di t^{i}(1-t)^{d-i}$ は Bernstein 基底)では、$\mathbf v(t)=\mathbf B(t)-\mathbf P$ は、制御点 $\mathbf p_i=\mathbf P_i-\mathbf P$ を持つ $d$ 次のベジェ曲線です。ステップ 5 の被積分関数を $N(t)/D(t)$ と書きます。

$D$ は $2d$ 次で、$N$ は素朴に数えると $2d-1$ 次です。ところが、$N$ は実際には $2d-2$ 次です。まず、冪の形で確かめます。

冪基底で数える

$\mathbf v(t)=\sum_{k=0}^{d}\mathbf c_k t^{k}$($\mathbf c_k$ は3次元ベクトル)と書くと、$\mathbf v'(t)=\sum_{k=1}^{d}k\,\mathbf c_k t^{k-1}$ で、

です。最高次($j=k=d$)の項は $\mathbf c_d\times\mathbf c_d=\mathbf 0$ で消えます。次の $t^{2d-2}$ の係数は $(d-1)\,\mathbf c_d\times\mathbf c_{d-1}+d\,\mathbf c_{d-1}\times\mathbf c_d=-\mathbf c_d\times\mathbf c_{d-1}$ です。結局 $N$ は高々 $2d-2$ 次です。

ただし、この方法は、ベジェ形式の制御点を冪基底の係数に直し、多項式どうしを掛けています。次数が上がると桁落ちが大きくなります(下の「数値的に安定」)。

ベジェ形式のまま作る(Bernstein 基底の積の公式)

冪基底に直さなくても、$N$ と $D$ は Bernstein 係数のまま作れます。$\mathbf v'$ は、制御点 $\mathbf q_j=d\,(\mathbf p_{j+1}-\mathbf p_j)$ の $d-1$ 次のベジェ曲線(ホドグラフ)です。Bernstein 基底の積は、次の公式で1つの基底にまとまります。

これを、$\mathbf v=\sum_i\mathbf p_iB_i^d$ と $\mathbf v'=\sum_j\mathbf q_jB_j^{d-1}$ の積に使うと、

となります。$N$ も $D$ も、制御点 $(k/m,\ n_k)$ や $(k/2d,\ D_k)$ を持つベジェ関数として、制御点の外積と内積だけで求まります。ツールが球の影で使う6次多項式 $(\mathbf a\cdot\mathbf v)^2-\cos^2\alpha\,|\mathbf v|^2$ も、同じ公式で作っています。

「最高次が消える」は、Bernstein 形式では次数上げとの整合になります。この $N$ は $2d-1$ 次の形ですが、実際は $2d-2$ 次の多項式を次数上げしたものです。次数を1つ下げる漸化式

を最後まで進めると、$b_{2d-2}=n_{2d-1}$ が(丸め誤差の範囲で)成り立ちます。これが、冪基底で「最高次の項が消える」ことに対応します。実験1で、係数を並べて確かめられます。

ベジェ形式で持つ利点は次のとおりです。

• 冪基底への変換が要りません。制御点の差と積の和だけで作れます。
• 数値的に安定です。制御点をランダムに選んだ $N(t)$ を、厳密な有理数の値と比べました(各次数 120 個の多項式、$|N|$ の最大値で正規化)。冪基底で作った場合の最大誤差は、$d=6$ で 4×10^-12、$d=10$ で 5×10^-8、$d=14$ で 2×10^-3 と、次数とともに悪化しました。Bernstein の積の公式では、どの次数でも 10^-15 台(最大 8×10^-16)でした。3次($d=3$)では、冪基底でも最大 4×10^-15 で、実用上の差は出ません。
• 凸包性が使えます。係数の最小値・最大値が、その区間での $N$、$D$ の下限・上限になります。$D_k$ の最小値は、P から曲線までの距離の2乗の下限です(第5.2節)。$n_k$ がすべて同符号なら、被積分関数の符号が区間全体で一定と分かります。
• $N=0$ の根も、ベジェクリッピングで求められます。これは、P から見た方向の回転面が水平になる点($\mathbf n\cdot\boldsymbol\gamma=0$)です。
• $N/D$ は有理ベジェ関数です。$N$ を $2d$ 次に次数上げして係数を $N_k$ とすると、$f(t)=\sum_kN_kB_k^{2d}\big/\sum_kD_kB_k^{2d}$ は、重み $D_k$、制御値 $N_k/D_k$ の有理ベジェ関数です。重みが正なら、凸包性がそのまま使えます。正でない区間は、de Casteljau 法で細分すれば、やがて正になります。

3次のベジェ曲線では、冪基底でも問題は出ません。ベジェ形式の利点が効くのは、高次のベジェ曲線($d\ge5$ ほど)を扱うときと、凸包性で範囲を囲みたいときです。

4.2 次数ごとの整理

$D(t)=|\mathbf v(t)|^2$ は、P が曲線上にない限り実数の $t$ では正で、実根を持ちません。根は $d$ 組の複素共役な対になります。部分分数に展開すれば各対から arctan と log の項が出て、原理的には閉じた式で書けます。しかし、そのために毎回 $2d$ 次方程式を数値的に解く必要があり、P や曲線が変わるたびにやり直しです。$d\ge2$ では、閉じた式を経由するより数値積分のほうが実用的です。線分($d=1$)が特別なのは、分子が定数になって根探しが不要になるからです。

4.3 被積分関数はどんな形か

被積分関数 $f(t)=N/D$ は、$t$ あたりの見込み角の増え方($d\varphi/dt$)に傾き $\mathbf n\cdot\boldsymbol\gamma$ を掛けたものです。P が曲線から離れていれば緩やかですが、P が曲線に近づくと、近い点のまわりで幅がおよそ「距離 ÷ 速さ」の鋭いピークになります(ローレンツ型)。光源が測定点に近いほど、この尖りが問題になります。

4.4 ガウス=ルジャンドル求積と適応分割

滑らかな被積分関数には、ガウス=ルジャンドル求積が向いています。$n$ 点の求積の誤差は $\rho^{-2n}$ 程度($\rho>1$ は被積分関数の解析性で決まる)で、指数的に減ります。$\rho$ は、区間 $[0,1]$ を焦点とするだ円のうち、最も近い複素特異点($D(t)=0$ の複素根)に接するものの大きさで決まります。

• P が曲線から遠い:特異点は実軸から遠く、$\rho$ が大きい。数点で十分。
• P が曲線に近い:特異点が実軸のすぐ近くに来て、$\rho\to1$。同じ点数では精度が出ない。

そこで、区間を二分して $\rho$ を稼ぎます。二分するかどうかは、その区間の曲線が P にどれだけ近いかで決めます。区間の弧長が P までの距離を大きく超えると、被積分関数が区間の中で大きく変わるからです。

距離の下限を、制御点から保証する(近似ベジェクリッピング)

区間の曲線の制御点(P からの相対座標)から、第4.1節の公式で $|\mathbf v|^2$ の Bernstein 係数 $D_k$($2d$ 次)を作ります。凸包性から、その区間で

が保証されます。最近点そのものを求める必要はありません。区間を de Casteljau 法で二分するたびに、下限は真の最小値に近づきます(区間幅の2乗で縮みます)。弧長の上界は、制御多角形の長さです。そこで、次の条件のとき二分します。

端点と中点の3点で距離を見積もる旧方式は、3点の間を曲線が通り抜けると、距離を過大に見積もることがあります。この方式は、距離の下限が保証されます。ただし、保証されるのは距離の下限で、求積の誤差そのものではありません。

もう1つの案は、最近点 $t^*$ を求めることです。$D'(t)=2\,\mathbf v\cdot\mathbf v'=0$ の実根($2d-1$ 次)なので、これもベジェクリッピングで求まります。$t^*$ で区間を分け、$t^*$ に向かって幅を半分ずつ縮めます。ここでも根を厳密に求める必要はなく、許容誤差 $10^{-3}$ で打ち切っても精度は変わりませんでした。

2,400 点で比べた結果です(円盤の半径 0.7〜1.5、高さ 0.08〜3、傾き ±0.6、制御点を ±0.25 ずらした 40 通りの光源に、平面 F 上の点を60個ずつ。基準は、ガウス16点で、弧長が距離の 0.1 倍を超えたら二分した値)。

• 同じ点数(8点)では、保証つき下限は、簡易判定と精度が同じで、時間が約2倍かかります。制御点から下限を作る分の手間です。
• ただし、保証があるので、二分の基準をゆるめ($\beta=2$)、そのぶん点数を上げて(16点)も、安全に使えます。この組み合わせが、簡易判定と同じ時間で、最大誤差を 4 桁小さくしました。ツールの既定にしています。
• 最近点で分ける方法は、根を求める分だけ遅く、この試験では保証つき下限に及びませんでした。ツールでは、選択式で使えます。

4.5 次数が効くもう1つの場所:交点の計算

境界を「水平面で切る」(第7節)、「影の輪郭と交わる」(第9節)には、曲線と平面の交点が要ります。曲線上の点が平面 $\mathbf m\cdot\mathbf v=0$(P を通る平面)にあるという条件は、$\mathbf m\cdot\mathbf v(t)=0$ なので $t$ の $d$ 次方程式です。線分($d=1$)なら1次方程式で済みます。3次以上の実根は、ベジェクリッピング法で求められます(第5節)。

5ベジェクリッピング法が使える箇所

第4.5節の「$d$ 次方程式の実根を求める」場面は、ベジェクリッピング法(Sederberg と Nishita、1990)でまとめて解けます。もとは曲線どうしの交点を求める方法ですが、多項式の実根探しにも使えます。ツールの根探しは、この方法に入れ替えてあります。

5.1 仕組み

多項式を Bernstein 形式 $p(t)=\sum_{i=0}^{n}b_i\,B_i^n(t)$ で持ちます。点 $(i/n,\ b_i)$ を制御点とするベジェ曲線が、ちょうど $y=p(t)$ のグラフです。ベジェ曲線は制御多角形の凸包の中に入る(凸包性)ので、$p(t)=0$ の根は、凸包が $t$ 軸と交わる範囲の中にしかありません。

1. 係数 $b_i$ がすべて同符号なら、根はありません。ここで打ち切ります。
2. そうでなければ、区間を、凸包が $t$ 軸と交わる範囲 $[t_{lo},t_{hi}]$ に縮めます。縮めた区間の係数は de Casteljau 法で作ります。
3. 区間があまり縮まなかったとき(元の7割以上が残るとき)は、二分してそれぞれを続けます(複数の根があるとき)。
4. 区間が十分に小さくなったら、その中点を根とします。

単根の近くでは収束が速く(2次収束)、ツールの試験では平均 5〜8 回の反復で終わりました。二分法は1回に1ビットずつなので、同じ精度に40回ほどかかります。実験3で、1回ごとの様子を見られます。

厳密な根が要らないとき(近似ベジェクリッピング)

ベジェクリッピングは、根や最小値を厳密に求めなくても役に立ちます。

• 打ち切りを粗くする。根を、区間の分割位置のような近似で足りる用途に使うなら、許容誤差を $10^{-3}$ 程度にして早く打ち切れます。
• 根や最小値を求めず、細分した各区間の制御点の最小値・最大値だけを見る。凸包性から、真の最小値の下限(最大値の上限)が保証され、細分するほどきつくなります。第4.4節の分割判定は、この使い方です(実験4)。

5.2 使えるところ、使えないところ

5.3 入れ替えて分かったこと

• 結果は変わりません。旧実装(冪基底に直して、導関数の根で区切って二分法)と、9,000 点で照度係数の最大差は 6×10^-15 でした。
• 速さもほぼ同じです。3次では旧実装も十分速く、64×64 格子で 77 ms が 79 ms でした。効くのは、次数が上がる場合(球の6次)と、根なしの組を打ち切れる影の計算です。
• そのまま作ると根を取りこぼします。区間の端にある根が、丸め誤差で符号が揃って消えます。縮めた区間に、ごくわずかな余白を足して解決しました(3・5・7次で各 1,500 個の多項式が、標本から数えた根の数と全件一致)。
• 内外判定・走査線・標本の判定は、結果も速さも同じで、折れ線への近似が消えました。光源領域の内外判定は、2万点で、256辺の折れ線による判定との不一致が0点でした(時間は 20 ms と 21 ms)。影つき照度は 5,520 点で最大差 0 でした。
• 自己交差の検出は、そのままだと遅くなりました。隣り合う曲線は、共有するアンカーでほぼ重なるため、クリッピングが収束しにくく、円盤で1回 25 ms かかりました。アンカーの近傍どうし(パラメータで 3%×3%)だけを除き、残りを2つの長方形領域に分けて調べると、0.4 ms になりました。400 個の形を、折れ線にした総当たりと比べ、全件一致しました。ただし、1本の曲線のループは、隣り合う部分を調べないので、パラメータで約 0.25 より小さいループは見逃します。

5.4 曲線どうしの交点:ファットライン

曲線 $P$ と $Q$ の交点は、次のように求めます。$P$ の両端を結ぶ直線を中心に、$P$ の制御点をすべて含む帯(ファットライン)を作ります。$Q$ の各制御点までの距離を Bernstein 係数とみなし、その凸包がこの帯と交わる範囲だけに $Q$ の区間を縮めます。次に役割を入れ替えて $P$ を縮め、これを繰り返します。縮まらなければ、長いほうを二分します。

• ツールは、光源の4本の曲線どうしで交点を調べ、交差があれば警告を出します。ランダムに崩した 300 個の形を、折れ線にして総当たりで調べた結果と比べ、全件一致しました。
• 隣り合う曲線は、共有するアンカーで接します。そこを交点と誤検出しないよう、アンカーの近傍(パラメータで 0.002 以内)は除外しています。
• 帯には、丸め誤差に備えて 10^-9 ほどの余白を足しています。

6多角形近似・モンテカルロ法との違い

6.1 3つの方針

折れ線は、1次のベジェ曲線の列です。第4.2節の表のとおり、閉じた式が得られるのは $d=1$ だけです。そこで、光源の形と積分のどちらを近似するか、あるいは面積分を確率的に解くかで、方針が3つに分かれます。

6.2 モンテカルロ法の中身

光源領域 $S$(面積 $A$)に $N$ 個の点を、$S$ の中で一様になるように置きます。それぞれを、面積 $A/N$、半径 $a=\sqrt{A/(\pi N)}$ の小さな円盤光源とみなします。円盤光源の照度は、円盤の中心軸上の点なら解析的に求まります。

($d$ は円盤中心までの距離。円盤の張る円錐が水平線の上にあるとき。)P が軸から外れる分は、円盤を斜めから見る分の $\cos\theta_Q$ を掛けて補います。全体では

です。分母の $d^2+a^2$ が要点です。小円盤を点光源とみなすと分母は $d^2$ で、小円盤の中心が P に近いと1個の値が大きくなり、ばらつきが増えます。円盤の大きさ $a$ を入れると、この発散が抑えられます。代わりに、$a^2/d^2$ 程度の小さな偏りが残りますが、$N$ が増えると円盤が小さくなって消えます。軸から外れる分の補正も、円盤が小さいほど正確です。

• 一様乱数:$N$ を4倍にして誤差がおよそ半分になります。10倍の精度には100倍の点が必要です。
• 準乱数(Halton など):点が領域を均等に埋めるので、誤差がほぼ $1/N$ で減ります。同じ $N$ で一様乱数より1桁近く正確になる例もあります。
• 影:P から各円盤の中心へ光線を飛ばし、遮蔽物に当たれば0とします。円盤ごとの二値なので、半影は $N$ が少ないとまだらになります。

6.3 誤差と手間の比べ方

同じ円盤を、(青)ベジェのまま $n$ 点のガウス求積(分割なし)、(緑)1曲線を $M$ 本の線分で近似した多角形、(紫)$N$ 個の小円盤のモンテカルロ法で解いた誤差を並べます。モンテカルロは、標本を変えた複数回の RMS 誤差です。

横軸は、「評価回数」と「実測時間」を切り替えられます。評価回数は、境界積分が $4n$(4本の曲線 × $n$ 点)、多角形が $4M$(辺の数)、モンテカルロが $N$ です。1回の評価の重さは方式で違うので、実測時間で見ると並びが変わります。

光源が遠い(高さ 3、$x=0$)と、ベジェのままなら 10 点ほどで機械精度(10^-15 程度)に届きます。多角形は 64 分割でも 10^-4 前後、モンテカルロは同じ手間(評価 256 回)で、一様乱数なら 3×10^-3 前後、Halton なら 1×10^-3 前後です。光源を P に近づける(高さ 0.1、$x=1$)と、どの方式も収束が遅くなります。モンテカルロは、小円盤が P に近い分だけばらつきが大きく、一様乱数では 256 個で 0.4 ほど、4096 個でも 0.14 ほどの相対誤差になります。

6.4 収束が速い・遅い理由

• 多角形($1/M^2$):曲率 $\kappa$ の曲線を長さ $h$ の弦で置き換えると、最大のずれは $\kappa h^2/8$ 程度です。立体角のずれはこれに比例し、$h\propto1/M$ から $1/M^2$ になります。分割を10倍にして、やっと2桁です。
• ガウス求積(指数的):被積分関数が滑らかな(解析的な)ので、数点増やすだけで桁が上がります。
• モンテカルロ($1/\sqrt N$):標準偏差は $\sigma/\sqrt N$ です。被積分関数の起伏($\sigma$)が大きい、つまり光源が近いほど遅くなります。

6.5 使い分け

ツールで測った一例です(円盤、高さ 2.2、遮蔽物なし、64×64 格子、ブラウザでの実測。環境や乱数で大きく変わります)。

球の影を入れると差が広がります。P の1点で、誤差 1% に届くまでの時間は、境界積分が 13 µs、モンテカルロが約 7 ms でした(ツールの「精度と時間のグラフ」の結果)。影の縁では、円盤ごとの二値判定がばらつきの原因になります。

7実際の計算手順

ツールが、P ごとに行っている処理の順番です。各手順で「次数」がどこに効くかも併記します。

1. 光源を組み立てる:光源平面の局所座標 $(a,b)$ で12個の制御点を持ち、世界座標に変換して4本の3次ベジェにします。
2. 発光側の判定:P が光源平面の表側にあるか($\mathbf n_L\cdot(\mathbf P-\mathbf O)>0$)。裏側なら $E=0$。
3. 向きの補正:局所座標での面積(シューレース公式)の符号 $\sigma$ で、境界の向きを判定します。
4. 水平面クリッピング:F より下の部分は $\cos\theta_P<0$ で、そのまま入れると誤った(負の)寄与が混ざります。$h(t)=\mathbf n\cdot(\mathbf B(t)-\mathbf P)$ は $d$ 次多項式なので、その実根(ベジェクリッピングで求めます。第5節)で曲線を切り、上側の区間だけ残します。切断で生じた2点は、F 上の直線(弦)でつなぎ、その寄与は角 $\operatorname{atan2}(\mathbf n\cdot(\mathbf a\times\mathbf b),\ \mathbf a\cdot\mathbf b)$ です。P が F 上にあるとき、P の水平面は F そのものなので、この切断はどの P でも同じで、1回だけ計算すれば済みます。
5. 曲線ごとの積分:残った各区間で、適応ガウス求積により $\int n\cdot(v\times v')/|v|^2\,dt$ を求めます。区間を分けるかは、制御点から求めた距離の保証つき下限で決めます(第4.4節)。
6. 合計:$G=\tfrac12\,\kappa\sigma\sum(\text{曲線の積分}+\text{弦の角})$。$\kappa\sigma=\pm1$ は向きの符号です。$E=L\,G$。

計算の中心にあたる部分を抜き出すと、次の形です。$\mathbf n=(0,1,0)$ なので、外積の内積は y 成分だけになります。

// 曲線区間 [ta, tb] の  ∫ n·(v×v') / |v|²  dt   (v = B(t) − P)
function gauss8(c, ta, tb, P) {
  const h = (tb - ta) / 2, m = (ta + tb) / 2;
  let sum = 0;
  for (let k = 0; k < 8; k++) {
    const t = m + h * GX[k];                    // ガウス点
    const v  = B(c, t) - P;                     // P から曲線上の点へ
    const dv = dB(c, t);                        // 接線ベクトル
    const num = v.z * dv.x - v.x * dv.z;        // n·(v×v')   (n = (0,1,0))
    sum += GW[k] * num / dot(v, v);
  }
  return sum * h;
}

// P ひとつ分の境界積分
let s = 0;
for (const piece of clip.pieces) s += integrate(curve[piece.ci], piece.ta, piece.tb, P); // F より上の曲線
for (const chord of clip.chords) s += atan2(n·(a×b), a·b);                              // F 上の弦
const G = 0.5 * KAPPA * sigma * s;              // E = L * G

8数値で追う

円盤の光源で、第7節の手順を数値で追います。P の真下から光源の縁の真下まで動かして、各曲線の積分値がどう変わるか、多角形との差がどうなるかを見てください。

9影の計算

光源と P のあいだに遮蔽物があるとき、境界積分のまま影(半影を含む)を扱えます。ツールが対応しているのは、直方体と球です。ここでは、考え方と、それ以外の形でどこが難しくなるかを整理します。

9.1 考え方

1. P から見た遮蔽物の輪郭 $\Sigma$ は、単位球面上の領域です。P から $\Sigma$ の内側の方向に進むと、遮蔽物に当たります。
2. 光が届く分は $G(S)-G(S\cap\Sigma)$ です。
3. $S\cap\Sigma$ の境界は2種類です。(a) 光源の境界曲線のうち $\Sigma$ の内側の部分。(b) $\Sigma$ の輪郭のうち光源 $S$ の内側の部分。

9.2 直方体(凸多面体)

• 輪郭 $\Sigma$ は、単位球面上の凸多角形で、辺は大円弧(=P から見た直線辺)です。P に向いた面と向いていない面の境界の辺をたどると、向きの揃った閉じた輪が得られます。
• (a) は、境界曲線と「P と $\Sigma$ の辺を含む平面」の交点を求めて区間に分けます。交点は $d$ 次方程式の根で、ベジェクリッピングで解きます。区間が $\Sigma$ の内側かは、区間の中点で判定し、内側の区間を第4節の積分で足します。
• (b) は $\Sigma$ の辺が直線なので、第3節の閉じた式です。
• 多角形近似では光源の辺も直線なので、(a) の交点は1次方程式になります。同じ枠組みで扱えます。
• この方法は、直方体に限りません。面のリストがあれば、角柱や角錐などの凸多面体にそのまま使えます。

9.3 球

球は、P から見ると円錐です。中心方向を $\mathbf a$、半頂角を $\alpha$($\sin\alpha=r/D$、$D$ は P から球の中心までの距離)とすると、$\Sigma$ は $\mathbf a$ を中心とする半頂角 $\alpha$ の円(小円)の内側です。

• (a) の交点:光源境界上の点が円錐の縁にある条件は $(\mathbf a\cdot\mathbf v)^2-\cos^2\alpha\,|\mathbf v|^2=0$ です。$\mathbf v(t)$ は $d$ 次なので、これは $t$ の $2d$ 次多項式(3次曲線なら6次)です。係数は制御点の積(Bernstein 形式の積の公式)で作れるので、ベジェクリッピングで実根を求めます。裏側の円錐($\mathbf a\cdot\mathbf v<0$)の根は除きます。
• (b) の弧:小円は $\boldsymbol\omega(\varphi)=\cos\alpha\,\mathbf a+\sin\alpha(\cos\varphi\,\mathbf e_1+\sin\varphi\,\mathbf e_2)$ と書け、次のとおり $\varphi$ で積分できます。閉じた式で求まります。
  $$\boldsymbol\omega\times d\boldsymbol\omega=\big[\sin^2\alpha\,\mathbf a-\sin\alpha\cos\alpha(\cos\varphi\,\mathbf e_1+\sin\varphi\,\mathbf e_2)\big]\,d\varphi$$
  円全体が光源の内側にあるときは $G=\pi\sin^2\alpha\,(\mathbf n\cdot\mathbf a)$ です。

光線で遮蔽を調べる面積分と比べ、約 4,000 点(うち影の影響あり約 1,000 点)で、差は $G$ で 1〜3×10^-4 でした(1回だけ 6×10^-3 の外れ値が出ましたが、4.8 万点の再試行では再現せず、参照値側の誤差と考えています)。

9.4 それ以外の形は難しいのか

難しい理由は、原理的に不可能だからではなく、(1) の輪郭を求める処理の実装量にあります。非凸の物体は、輪郭の一部が別の部分に隠れます。ドーナツは、外側の輪のほかに、向きの逆の輪(穴)ができます。ベジェ曲面や回転体の輪郭は、$\mathbf N(u,v)\cdot(\mathbf S(u,v)-\mathbf P)=0$ を満たす点の集まりです。これは、1変数の多項式の根ではなく、2変数の陰関数曲線で、曲線としてたどり、隠れる部分を除く必要があります。

9.5 曲面の遮蔽物とベジェクリッピング

遮蔽物がベジェ曲面(球、ドーナツ、回転体をパッチで表したもの)のとき、ベジェクリッピングが効くのは、光線と曲面の交点です(Nishita・Sederberg・Kakimoto、1990)。光線が曲面パッチに当たるかを、パッチの制御点の凸包で絞り込んで判定します。

• モンテカルロ法は、この判定だけで影を扱えます。境界積分のように輪郭を求める必要がありません。
• 面積分による検算(第6節)の遮蔽判定にも使えます。
• 境界積分に必要なのは輪郭で、光線との交点そのものではありません。輪郭を求める部分(2変数の陰関数曲線)には、ベジェクリッピングはそのままでは使えません。
• 回転体は、光線を (r², y) 平面に写すと放物線になり、輪郭の曲線(プロファイル)との交点を、2つの多項式ベジェ曲線の交点としてクリッピングで判定できます。これも、光線ごとの判定です。

ツールは、ドーナツと回転体には対応していません。境界積分で解くには、上の「難しい」行の処理を作る必要があります。

9.6 現実的な選択肢

1. 凸な物体は、凸多面体で近似する。球や楕円体を面で刻むと、箱と同じ処理に載せられます。多角形近似と同様に、面の数を増やすほど誤差が縮まります。
2. 輪郭を折れ線にして、非凸の物体の和集合を作る。輪や穴を持つ多角形に対する処理が要ります。実装量は大きい代わりに、ノイズのない結果が得られます。
3. 併用する。影のない場合の照度を境界積分で正確に求め、遮蔽される分だけをモンテカルロで見積もります。$G(S)$ から、遮蔽された分のモンテカルロ推定を引く形で、ばらつきの大きさは、遮蔽された割合が小さいほど小さくなります。影が小さい範囲では、モンテカルロだけより精度が上がります。
4. すべてモンテカルロにする。形を選ばず、複数の遮蔽物も扱えます。ばらつきは、標本数で抑えます。

2と3は、ツールには入れていません。

10曲面光源について

境界積分が使うのは、光源の面ではなく、P から見た光源の輪郭だけです(ステップ 2 の $\Omega$)。式の導出は、光源が平面であることを使っていません($\mathbf v(t)=\mathbf B(t)-\mathbf P$ は3次元の曲線なら何でもよい)。したがって、光源面が曲がっていても、次の条件がそろえば、境界曲線だけで照度が決まります。

• 光源の輝度が一様で、完全拡散である(ステップ 2 で $\cos\theta_Q$ が消える条件)。
• P から見えるのが光源の表側だけで、光源のどの点も P に向いて発光している。
• 同じ視線が光源に2回当たらない(自己遮蔽がない)。

このとき、曲面の中の形は照度に影響せず、輪郭 $C$ の3次元形状だけが効きます。ただし、ツールの実装は平面光源を前提にしています(発光側の判定、水平面クリッピング、影の内外判定に「境界が同一平面上にある」ことを使っています)。曲面光源を扱うには、これらを3次元の輪郭に対して作り直す必要があります。

11つまずきやすい点