ブログ移行しました -> MochiMochi 3D

この記事はレイトレ(Raytracing) - Qiita Advent Calendar 2025 - Qiitaの一日目の記事です。モンテカルロレイトレーシングにおけるNEE(Next Event Estimation)の名称の由来について調べてみた記事です

NEE(Next Event Estimation)

モンテカルロレイトレーシングにおけるNEE(Next Event Estimation)とは光源点と衝突点を接続するような方向をサンプリングするモンテカルロ法のテクニックの一つです。これ自体はモンテカルロライトレーシングではとてもメジャーな手法です。詳しい理論や実装に関してはshocker氏の記事や私の記事を参照してください。

rayspace.xyz

kinakomoti321.hatenablog.com

さて、NEEを勉強したことがある人はきっと一度はこう感じたことがあるのではないでしょうか。

「名前が意味不明すぎる」

何がNext Eventで、どこをEstimationしているんだと疑問に思ったことはあるはずです。私も勉強したての頃からずっと疑念を抱いておりました。

ということで、今回の記事はnext event estimationという名称がどこから来たのか何が由来なのかを調査してみました。

そもそものNEEとは何なのか

とりあえず、PBRのバイブルであるpbrt-4th editionを参照してみます。13,2.3 Incremental Path ConstructionにNEEに対する言及が見つかりました。

pbr-book.org

With path tracing, the last vertex of the path, which is on the surface of a light source, gets special treatment. Rather than being sampled incrementally, it is sampled from a distribution that is just over the surfaces of the lights. (Sampling the last vertex in this way is often referred to as next event estimation (NEE), after a Monte Carlo technique with that name.)

簡単に翻訳すると

パストレーシングでは光源の表面上にあるパスの最終頂点は特別な扱いを受ける。この最終頂点は逐次的にサンプルされるのではなく、光源上の分布からサンプルされる。(この方法による最終頂点のサンプリングはモンテカルロ法の手法にちなんでnext event estimation(NEE)と呼ばれることが多いです。)

この文面からNEEについて以下のことが読み取れます。

3DCGとは無関係に元からNEEというモンテカルロ法のテクニックが存在した
NEEの本質的な所は経路の最終頂点を強制的に決定するところにある

NEEはよく直接光のlight samplingと混合されがちなイメージがありますが(私も思っていました)、どうやらそういうわけではなくあくまでNEEは「最終頂点を決定する手法」のようです。なので、恐らく光源点以外を最終頂点としてサンプリングしたとしても、それもNEEと呼ぶことができると思います(なので NEE $\neq$ Light sampling)。

NEEの起源

pbrt-4th editionではFurther Renderingの項にNEEの元論文としてAdjoint and importance in Monte Carlo application[Coveyou, et al. 1967]が提示されています。

pbr-book.org

この論文を見てみると228ページに確かにNext-Event Estimationの文字が確認できました(課金必須)。一応この論文より前にNEEの名前が使われていないかをGoogle Scholarで確認してみましたが、特に見当たらなかったので、NEEという名称の起源はこの論文であることは正しそうです。

Google Scholarの検索結果 [Coveyou, et al. 1967]以前の論文は発見できず

ここで語られるNext-Event Estimationについて論文を引用しつつちょっと見ていこうと思います。まずバックグラウンドですが、この論文は3DCGではなく、中性子のtransportの計算を目的としたもので(おそらく原子力周りの研究っぽいです)、それにおける効率的なモンテカルロ法を提案するというものです*1。

簡単に問題設定を説明します。この論文の目的は空間内に粒子(ここでは中性子)が散乱、吸収するなどして飛び回ってる状況において、空間中の中性子の量を推定することです(最終的にはセンサーに入る量みたいなのをシミュレーションするみたいです)。

このモデルでは粒子に関する何らかの事象をイベント eventと呼びます。このモデルでは初期(initial)と消失(lethal)、遷移(transition)の3つのイベントを考えます。

初期イベントというのは単に初期条件で、最初に粒子は位置$x$に対して $S(x)$ によって分布していることにします。これはソース関数と呼び、レンダリングで言う所の光源みたいなものです。これは分布であるため、

$\displaystyle{ S(x) \ge 0 \tag{1.1} }$

$\displaystyle{ \int S(x) dx = 1 \tag{1.2} }$

を満たします(あんま議論の本質ではないので気にしなくていいです)。

対照的に消失と遷移は粒子に対する作用であり、レンダリングの文脈で言う吸収と散乱に相当します。このイベントの表現として、ある位置 $x$ に対して消失する確率を $\alpha(x)$ 、遷移が起きる確率を $\sigma(x)$ と書きます。

遷移はある種の散乱みたいなもので、ある位置 $x$ から $y$ へと粒子が移動する確率密度を遷移カーネル(transition karnel)と呼び、 $K(x,y)$ と書きます。遷移カーネル $K(x,y)$ は次のような性質を持ちます。

$\displaystyle{ K(x,y) \ge 0 \tag{1.3} }$

$\displaystyle{ \sigma(x) = \int K(x,y) dy \tag{1.4} }$

行ってしまえば遷移カーネル $K$ はBSDF的な存在です。

こうして様々なイベントを経てある地点における粒子の確率密度を求めることにします。これをイベント密度 $E(x)$ と呼び、次のような方程式で求められます。

$\displaystyle{ E(x) = S(x) + \int K(x,y) E(y) dy \tag{1.5} }$

ここから粒子の割合が分かるわけなので、センサーに検出される粒子数など空間における粒子を把握することが可能になります。論文の大目的はこれの数値計算をすることです。

イベント密度 $E(x)$ の数値計算法としてモンテカルロ法が当時既に使われていたらしく、効率的な計算手法を求めていたというバックグラウンドがありました。この論文はそれに対してNEE(などの様々な手法)を提案したという流れになっています。

ではNext-Event Estimationがどんなアルゴリズムだったか、論文の記述を追ってみます。まず、estimatorを $\eta$ と定義し、経路 $(x_0, x_1, ..., x_n)$ における各点の寄与(contribution)を $P$ として次のように表現することにします。

$\displaystyle{ \eta = \sum_{k=0}^n P(x_k) \tag{1.6} }$

寄与 $P$ はいわゆる評価値とpdfで割ったアレです。具体的な値は考えず、とりあえずそのように表現しているという認識でOKです。期待値 $\langle \eta \rangle$ の値をここでは $\lambda$ と表記します(論文ではイベント密度に関わらず一般的なモンテカルロ法の議論をしています、イベント密度に適用した場合は $\lambda$ は $E(x)$ です)

$\displaystyle{ \lambda = \langle \eta \rangle \tag{1.7} }$

これに対して、Next-Event Estimationは次のように説明されています

The principle title of this section is enclosed in quotation marks because the name, though em bedded in the literature, is a thorough misnomer. The method differs from other estimation meth ods, not because the desired answers are esti mated statistically, since all Monte Carlo results are, but because of the way in which these estimates are made. The secondary title is more descriptive of the technique. In the most usual form of the method, the estimator $\bar{\eta}$ chosen to replace

$\displaystyle{ \eta = \sum_{k=0}^n P(x_k) \tag{1.8} }$

for a particle of biography $(x_0, x_1, ... , x_n)$ is given by

$\displaystyle{ \bar{\eta} = P(x_0) + \sum_{k=0}^{n} W_1(x_k) \tag{1.9} }$

with

$\displaystyle{ W_1(x) = \int K(x,y) P(y) dy \tag{1.10} }$

as defined in Eq. (25). We will see that $\bar{\eta}$ . The choice of this estimator is based upon the quite sound intuition that replacing the contribu tion $P(x)$ of an event at $x$ by $W_i(x)$ (the expected value of this contribution averaged over all pos sible next events) may reduce the variance of the estimate. The contribution of initial events thus neglected is accounted for by the term $P(x_0)$ . The procedure here is that of Fig. 7.

これをまとめますと、NEEとはestimator $\eta$ の $x_0$ 以降の寄与を $W_1$ という積分値に置き換えた新しいestimator $\bar{\eta}$ を使う手法です。この手法のコアは $W_1$ というウェイトで、定義を見るとあらゆる $y$ から $x$ に来る寄与 $K(x,y)P(y)$ を計算したものであることが分かります。

論文では「この推定量の選択は、ある事象 $x$ における寄与 $P(x)$ を $W_i(x)$ （この寄与の期待値を全ての可能な次事象について平均したもの）で置き換えることで推定量の分散を低減できるという、非常に妥当な直観に基づいている」(Deepl翻訳)と述べています。パストレしている方ならわかると思いますが、この手の問題のvarianceの原因は高い寄与を持つ地点がとても限られており、サンプルを取ることが上手くできないことにあります。ここにおけるNEEは平均値(積分値)を使うことでvarianceを抑えよう！という直感的な考え方に基づいて考案されたものっぽいです。

それでvarianceが低くなるのは納得ですが、そもそも平均値に積分が入ってるやろどうやって計算すんねんという話です。ですが、どうも論文に具体的な $W_1$ の計算方法が記載されていませんでした。アルゴリズムの解説を見ても単に「 $W_1$ を足す」としか書いていません。事前計算するなり、新しくestimatorを作るなりで $W_1$ の計算をする感じだったのかと思います。

この辺は論文でも「varianceが減少する保証があるかどうかは分からない」と述べていたりとこの辺は結構曖昧というか議論が十分にされていない感じがあります。記述もシンプルなので恐らく紹介程度のノリで書かれていたのかなと思います。

Although it seems quite plausible that this proce dure is variance-reducing in general, no proof of this was found, and, indeed, it may well be untrue. It is true that in many problems the use of this technique pays handsome dividends in reduced variance.

$W_1$ の式を見ると次のイベントの寄与を全て確認しているわけですから、ある種の推定をしていると解釈してもよさそうです。つまり、次(next)のイベント(event)の推定(estimation)をしているという意味でNext-Event Estimationという名称がつけられたのだと推測できます。確かにこの文脈でつけるなら正しい名称のように感じます。こうしてNEEの名称が生まれたわけですね。

あと、重要な点として、おそらくNEEのアイディア自体はこの論文が初出ではないと思われます。NEEの冒頭に「このセクションの主要タイトルが引用符で囲まれているのは、その名称が文献に定着しているにもかかわらず、完全に誤った呼称であるためである。」とあり、別の名前(Statistical Estimation)で呼ばれていたことが伺えます。従って、[Coveyou, et al. 1967]はNext-Event Estimationという名称の起源ではあるが、アルゴリズムの開祖ではないと思われます。

こうして調べてみると、[Coveyou, et al. 1967]のNEEは現在のNEEと比べると全然内容が大きく違うように見えます。少なくともNEEで次の寄与の積分値を取得しているというわけはありませんし、pbrtの説明のように直接的にパスをつなぐみたいなことをは特に明確な言及はありません。次は現在と過去のNEEがどのように繋がっているのかを調べていこうと思います。

第二種フレッドホルム積分方程式

ちょっと話は変わりますが、イベント密度の方程式はレンダリング方程式とそのまま対応付けることができます。こうしたレンダリング方程式や(式番号)のような自身が積分として含まれているような再帰的な積分方程式はこうした方程式はちゃんと名前があり、第二種フレッドホルム積分方程式(Fredholm Integral Equation of Second Kind)と呼ばれています。

$\displaystyle{ f(x) = g(x) + \int K(x,y) f(y) dy \tag{2.1} }$

この $K$ のことをカーネル(kernel)と言い、先ほどの例では遷移カーネル、レンダリング方程式ではBSDF+cosine項に相当します。

つまるところ、[Coveyou, et al. 1967]の問題とレンダリング方程式は本質的には同じであるわけです。よって、分野は違えど[Coveyou, et al. 1967]で提案されたNEEはレンダリング方程式にも応用可能です(抽象度で言えばイベント密度の方程式の方が高いですし)。ということで、このNEEが3DCGにも持ち込まれ、なんやかんやあることで現在のNEEのアルゴリズムに落ち着いた、という歴史的な経緯が推察できます。

3DCG分野におけるNEE

初期のNEEは原子物理学と分野がかなり離れていますが、これがいつ頃3DCGに応用されたのかが気になるところです。Google Scholarで調べてみる限り、Adjoint Equations and Random Walks for Illumination Computation[Pattanaik, Mudur, 1995]と言う論文が初めて明確にNEEをレンダリング方程式に応用したと思われます。

[Pattanaik, Mudur, 1995] 5.1 Next-Event Estimationより引用

[Pattanaik, Mudur, 1995]の式を見る限り、手法としては[Coveyou, et al. 1967]のNEEをそのままレンダリング方程式に応用したという感じです。私たちが知るレンダリング方程式のestimatorは(2)の方になります( (1)は多分光源->カメラのパストレーシングの式)。 $\epsilon$ は光源からの放射輝度、 $f_r$ はBSDFを表しています。ここでNEEのウェイトとなる $L^1$ はdirect lightの積分値として定義されてます。

しかしながら、[Coveyou, et al. 1967]と同様に具体的に積分値をどうやって計算するのかは特に提示されていないようでした。また、(2)の式を見てもらえるとわかるんですが、direct lightの積分値をただ足していく形で寄与が構成されており、かなり怪しい式をしています。実際はindirect lightとして計算する必要があるため、前の反射の値が入ってくるはずなのですが、ここではそう言ったことはされてません。恐らくはPDFによってその辺が打ち消される、という仮定を置いている...のかもしれません(後に紹介する論文もそうだった)。(当時のモンテカルロ感がいまいち把握できない)

ただここで重要なのが、ウェイトの中に使われている被積分関数が光源 $\epsilon$ とカーネルで表されており、これ以上の展開を必要としません。[Coveyou, et al. 1967]ではこの辺かなりぼかして書かれていたのでよくわかりませんでしたが、これでちょっと明確になりました。あるパス $(x_0, x_1, x_2, ... , x_n)$ が与えられたら、NEEは各点のlocal directionを計算してそれを加算するということを行うアルゴリズムであることが少なくともここの式から読み取れます(暗黙的にindirect lightの係数がかかったうえで)。

こうしてみるとかなりNEEっぽさが出てきました。もう少し明確なpdfを使った議論をしてみれば繋がりが見えてきそうです。Google Scholarを探ってみると、Mathematical Models and Monte Carlo Algorithms for Physically Based Rendering[Lafortune, 1996]という論文が分かりやすくNEEを定式化しているのを見つけました。

ここではまず、一般的な第二種フレッドホルム積分方程式

$\displaystyle{ f(x) = g(x) + \int K(x,y) f(y) dy \tag{3.1} }$

に対するNEEの適用を考えます。今回の解説では論文に従いつつ、分かりやすさのためにいくつか変更を加えて話していきます。この論文では次のようなtex: hを定義して、NEEを構築します。

$\displaystyle{ h(x) = \int K(x,y) f(y) dy \tag{3.2} }$

この $h$ は再帰的に展開すれば次のように書くことができます。

$\displaystyle{ \begin{align} h(x) &= \int K(x,y) f(y) dy \\ &= \int K(x,y) (g(y) + h(y)) dy\\ &= \int K(x,y) g(y) dy + \int K(x,y) h(y) dy \tag{3.3}\\ \end{align} }$

ある関数 $F$ のestimatorをここの節では $\langle F \rangle$ と表記することにします(前の節とは表記がごっちゃになって申し訳ないですが、論文の表記を重視します)。 $f(x)$ のestimator $\langle f(x) \rangle$ は

$\displaystyle{ \langle f(x) \rangle = g(x) + \langle h(x) \rangle \tag{3.4} }$

となり、 $h(x)$ のestimator $\langle h(x) \rangle$ は以下のように書きます

$\displaystyle{ \begin{align} \langle h(x) \rangle &= \langle K(x,y) g(y) \rangle + \langle K(x,y) h(y) \rangle \tag{3.5} \end{align} }$

ここで、第二項目のestimatorについて次のように変形できます(期待値は変わらないので大丈夫...たぶん)

$\displaystyle{ \langle K(x,y) h(y) \rangle = \langle K(x,y) \langle h(y) \rangle \rangle \tag{3.6} }$

それぞれのestimatorは独立なので各々サンプリング戦略 $P_1, P_1'$ を持ち、それぞれのサンプルを $\xi_1,\zeta_1$ とします。こうした時、estimatorは具体的に

と書くことができます。これを再帰的に繰り返すと

$\displaystyle{ \begin{align} \langle K(x,y) g(y) \rangle + \langle K(x,y) h(y) \rangle &= \frac{K(x,\zeta_1) g(\zeta_1)}{P_1'(\zeta_1)} + \frac{K(x,\xi_1)}{P_1(\xi_1)} \langle h(\xi_1) \rangle \\ &= \frac{K(x,\zeta_1) g(\zeta_1)}{P_1'(\zeta_1)} + \frac{K(x,\xi_1)}{P_1(\xi_1)} \frac{K(\xi_1,\zeta_2) g(\zeta_2)}{P_2'(\zeta_2)} +\frac{K(x,\xi_1)}{P_1(\xi_1)} \frac{K(\xi_1,\xi_2)}{P_2(\xi_2)} \langle h(\xi_2) \rangle \\ &= \frac{K(x,\zeta_1) g(\zeta_1)}{P_1'(\zeta_1)} + ... + \lbrace \prod_{k=0}^n \frac{K(\xi_{k-1},\xi_k)}{P_k(\xi_k)} \rbrace \frac{K(x,\zeta_i) g(\zeta_i)}{P_i'(\zeta_i)} + ... \tag{3.8} \end{align} }$

従って、 $f$ のestimatorは次のように書くことができます。

$\displaystyle{ \begin{align} \langle f(x) \rangle =g(x) + \sum_{n=0}^{\infty} \lbrace \prod_{k=0}^n \frac{K(\xi_{k-1},\xi_k)}{P_k(\xi_k)} \rbrace \frac{K(x,\zeta_i) g(\zeta_i)}{P_i'(\zeta_i)} \tag{3.9} \end{align} }$

NEEのウェイト $W_1$ を[Pattanaik, Mudur, 1995]に従って定義すると

$\displaystyle{ \begin{align} W_1(x) = \int K(x,y) g(y) dy \tag{3.10} \end{align} }$

であるわけですから、 $W_1$ を同様にモンテカルロ法をするとしてそのestimator $\langle W_1 \rangle$ を使えば

$\displaystyle{ \begin{align} \langle f(x) \rangle =g(x) + \sum_{n=0}^{\infty} \lbrace \prod_{k=0}^n \frac{K(\xi_{k-1},\xi_k)}{P_k(\xi_k)} \rbrace \langle W_1(\xi_k) \rangle \tag{3.11} \end{align} }$

とも表現することができます。つまり $W_1$ の計算はモンテカルロ法で一緒くたに計算するということを意味しています。

ここで重要なのが $W_1$ のestimatorが取るサンプル戦略は完全に独立しているという所です。はっきり言って $W_1$ をモンテカルロ法でただ解くとしたらそもそもNEEをやる意味がありません(それはナイーブにモンテカルロ法をしているのと何も変わらない)。しかし、ここでは $W_1$ の計算にはサンプリング戦略 $P_i$ とは異なる $P_i'$ を使って、独立したサンプル $\zeta_i$ を使用しています。

パス $(x, \xi_1, \xi_2, ... , \xi_n)$ がサンプルされたとして、各地点における $W_1$ はまた別のサンプリング戦略によって生成されたサンプル $\zeta_i$ を使って計算されます。そして、各ウェイトを適切な係数をかけて総和を取れば(3.11)を求めることができます。

ここで $W_1$ のサンプリング戦略にlight samplingを使用すれば、これは確かに現在のNEEのアルゴリズムであることが分かります。従って、[Coveyou, et al. 1967]は現在3DCGで使用されているNEEの起源であることが確認できます。

また、ここからNEEの本質的な所は「評価するestimatorとパスを選択するestimatorが独立している」ことで、NEEのアルゴリズムで言う所の「移動して、サンプル取って寄与を計算して、また移動する」部分が原義のNEEであると言えます。パスを大量に作っているともとれるので、パスの終端という視点で見ればpbrtの説明も正しいと言えます。

結論

Next Event Estimationの名称は物理学的な文脈で使われる"event"に由来していた。

おわり

ここまで見て頂きありがとうございました。

レイトレ合宿でNEEって名前変じゃない？みたいな話題が出てきたのが今回の記事のきっかけでした。前々よりこの疑問が頭の片隅にずっとあったので、せっかくなので調査をしてみた、というのが今回の記事です。

調べてみて、まさか原子物理学と関わっているとは思ってもいなかったので結構驚きました。色々見てみると意外とRendering Equationと中性子の散乱が同じ問題であり、なかなか面白い繋がり方しているな～と感じました(まあ、対象が光子から中性子に変わっただけでしかないですしね)。最近はこういう意外な歴史の繋がりを見るのが好きなので、なかなか楽しんで調査できました。

今回の記事に載せた内容はかなり要点だけ抜き出した感じではあるので、実際はもっと深い歴史があると思います。実際1980年代ぐらいの原子物理学の論文でちょこちょこNEEの名前がヒットしていて、恐らくそっちの方で色々発展とかがあったんだろうな～と推測しています(大半が課金しないと閲覧できなかったので断念)。3DCGの方でも資料が残ってないだけで、多分色々あった気がします(もう人に聞きに行く必要がありそう)。とりあえず自分が納得できる範囲でやりました。

正直これが正しいのかよくわからないので、誰か偉い人がまとめてくんないかな～

Reference

Adjoint and importance in Monte Carlo application[Coveyou, et al. 1967]
Adjoint Equations and Random Walks for Illumination Computation[Pattanaik, Mudur, 1995]
Mathematical Models and Monte Carlo Algorithms for Physically Based Rendering[Lafortune, 1996]
https://pbr-book.org/4ed/Light_Transport_I_Surface_Reflection/Path_Tracing
https://en.wikipedia.org/wiki/Fredholm_integral_equation

*1:後でわかったんですがモンテカルロ法自体が中性子のtransportを解くために考案されたらしいです。なのでこっちが実は本場