Specular Polynomial
SIGGRAPH 2024で発表されたCasticsレンダリングの手法、Constraintを多項式表現で表し、高速な計算を行う手法である。
https://scrapbox.io/files/6663f7a099bb51001da2a8a3.png
基本的にcausticsへのサンプリングは間にあるSpecular面におけるCnstraintをパスの関数Cとして表現して、ちゃんと成り立つパスならC=0、そうで無い場合は非ゼロとなるように設計して、C=0を満たす解を探索する形で行われる。
前手法のSMSは解の探索にニュートン法を使っていたが、ニュートン法の欠点である収束の遅さとか苦手な部分があった。
この論文は解の求め方に着目し、Constraint CをPolygonalとして表現し、連立方程式の解を求めるという形で高速化を図った。
あまり知らないけど、Polygonalの連立方程式は高速に解く方法があるらしい
Specularのトラアングルをなんやかんや(Wang et al. 2020等)求めるとする。$ k個のTriangle$ T_1,...,T_kがあるとして、それぞれのTriangleにおける3頂点とその法線を$ P_i =(p_{i,0},p_{i,1},p_{i,2})、$ N_i = (n_{i,0},n_{i,1},n_{i,2})と表す。 この時、Specular経ちをつなぐパスの各頂点を$ x_1,...,x_kとして、各トライアングルにおける重心座標$ U_i =(u_i,v_i)を用いると、重心座標を使えば各頂点は
$ x_i = (1-u_i - v_i) p_{i,0} + u_i p_{i,1} + v_i p_{i,2} = P_i u_i
と表される。
同様に法線も以下のように表すことができる。
$ n_i = (1-u_i - v_i) n_{i,0} + u_i n_{i,1} + v_i n_{i,2} = N_i U_i
$ \hat{n_i} = n_i / |n_i|
ここでは補完で出てきた法線を$ n_i、それを正規化したものを$ \hat{n_i}という風に明確に分けていることに注意（あとあと正規、非正規を考えるので）
Contraint
Specular面では(理想の鏡面を考え)反射屈折の法則にしたがった方向にしか光は行くことはない。その条件を数式にすると、ハーフベクトル$ h_iを使って
$ h_i \times n_i = 0
と表せる。要は$ h_iと$ n_iが完全に一致するときということをあら合わしている。(わざわざクロス積なのは答えの時0だと方程式の解として答えが出てくるので便利だから)
ただし、ここではハーフベクトル$ h_iは屈折も考慮して
$ h_i = \eta_i \hat{d}_i - \eta_{i-1} \hat{d}_{i-1}
$ \eta_iは$ i番目の屈折率を表し、$ d_iは次に進む方向ベクトルを表している。
$ \hat{d_i} = d_i / ||d_i||, d_i = x_{i+1} - x_i
このConstraintが全ての$ x_iで成り立つような$ x_iを解として求める
前手法であるSMSではニュートン法を用いてこの解を求めていたわけである(なので法線の微分とか求める必要があったし、なんか行列積とか必要だった)
この手法では、重心座標$ U_iについてをなんやかんやしてConstraint をPolynomialの方程式に変換して、うまいこと解を求める形でやる。
polynomial
今のConstraintの式には正規化された方向ベクトルがあるので、少なくともRootが式に含まれていてPolynomialな式ではない。これを排除してPolynomialな式へと変換する。
- Reflectionの場合
反射の場合、光路面に対して法線は平行になる。そのため、光路面の法線$ d_{i-1} \times d_iとは垂直な関係があり、式で表すと
$ (d_{i-1} \times d_i) \cdot n_i = 0 \tag{5}
となる。ここでは別に光路面の法線は正規化してもしなくてもConstraintの式は変わらないので、Rootを省いたConstraintが得られた
Directionを展開して頂点形式で書くと
$ (d_{i-1} \times (x_{i+1}-x_{i-1})) \cdot n_i =0 \tag{6}
という形になる。
- Arbitary Constraint
屈折(Refraction)を含めた任意のConstraintはちょっと複雑になる。スネルの法則、
$ \eta_{i-1} \sin{\theta_{i-1}} = \eta_i \sin{\theta_{i}}
と外積の関係性、
$ |n_i| \sin{\theta_i} = |\hat{d_i} \times n_i|
を用いれば、スカラー量のConstraintとして以下の式を得ることができる。
$ \eta_{i-1} |\hat{d}_{i-1} \times n_i| = \eta_i |\hat{d_i} \times n_i|
また、$ \hat{d}_{i-1}, \hat{d}_{i}は1つの光路面上に存在するため、ここから法線と外積取った結果のベクトルは共に平行である
$ \frac{\hat{d}_{i-1} \times n_i}{|\hat{d}_{i-1} \times n_i|} = \frac{\hat{d}_{i} \times n_i}{|\hat{d}_{i} \times n_i|}
これを用いればベクトルのConstraintとして次を得ることができる
$ \eta_{i-1} \hat{d}_{i-1} \times n_i = \eta_i \hat{d} \times n_i(9)
この式は$ \hat{d}_{i-1}, \hat{d}_iという正規化ベクトルを使ってるので中身としてはRootが使われてしまっている。これはどうにかしたい。
- Polynomialize
(9)式はRootが含まれるのでこっからPolynomialな式へと変換する
Rootを取る方法は2種類ある、一つはシンプルに二乗する方法(Square Form)、もう一つドット取ってあげる方法(Productr Form)がある（こっちはReflectionのみ）。
Square Form
(9)式をとりあえずなんかのベクトル$ b(実用上は(1,0,0)^Tとかシンプルなベクトルで)でスカラー値の方程式にする。
$ \eta_{i-1} ((\hat{d}_{i-1} \times n_i) \cdot b) = \eta_{i} ((\hat{d}_i \times n_i) \cdot b)
これに対して、二乗して$ d_i^2 d_{i-1}^2をかけてあげると
$ \eta_{i-1}^2 d_i^2 ((d_{i-1} \times n_i) \cdot b)^2 = \eta_{i}^2 d_{i-1}^2((d_i \times n_i) \cdot b)^2 \tag{10}
という形のConstraintを得られる。
Product Form
Reflectのみ、これは反射の法則から
$ \hat{d}_{i-1} \cdot n_i = - \hat{d}_i \cdot n_i
が成り立つ。ここで$ nに対して垂直な任意のベクトル$ t_i(例えば Triangle tangent)とは$ \hat{d}_{i-1},\hat{d}_{i}については
$ \hat{d}_i \cdot t_i = \hat{d}_{i-1} \cdot t_i
この2つの式をかけて表すと
$ (d_{i-1} \cdot n_i)(d_i \cdot t_i) + (d_{i-1} \cdot t_i)(d_i \cdot n_i) = 0 \tag{13}
と得ることが出る。
- Multivariate Specular Polynomials
こうしてConstraintを作ることは出来たが、各$ U_i全てにあるので全体として2k変数の連立方程式として存在することになる。
また、方程式の次数はrefractionだと6、reflectionなら4となるなる。もし仮にGeometry normalを使うのであればノーマルを定数とみなして4、2に次数を抑えることは可能(果たして使うのかと思うが)。
多変数の連立多項式の解を求める方法はあるみたいだが計算量の観点から、基本的に変数が多いのはお勧めされない。
少ない方が安定かつ早く、一番堅実なsolverは2変数を対称としている
ということで、今回の2k変数のPolynomial Systemをどうにか変数を取り除くこと(variable reduction)で2変数へと落とし込む形で高速化を行う
- variable reduction
Specularというのは(粗さがなければ)レイの進む方向は固定されるため、前のdirectionから一意的に位置と方向ベクトルは決定される。つまり、それを辿っていくとSpecular chainは最初のSpecular面の時点でパスは決定されていることがわかる
なので、本質的には$ U_iというのは$ U_1の関数として表現される。論文ではTriangleのintersectionの関係式を使うと漸化式として、$ U_{i+1},U_{i},U_{i-1}の関係式を示している。この形式をRational Coordinate Mappingと呼ぶ
$ U_{i+1}(U_i,U_{i-1}) = \frac{(\tilde{u}_{i+1}(U_{i},U_{i-1}),(\tilde{v}_{i+1}(U_{i},U_{i-1}))}{\kappa_{i+1} (U_i,U_{i+1})}
各関数については
$ \tilde{u}_{i+1}(U_{i},U_{i-1}) = (\tilde{d}_i \times e_{i+1,2}) \cdot (x_i - p_{i+1,0})
$ \tilde{v}_{i+1}(U_{i},U_{i-1}) = ((x_i - p_{i+1,0}) \times e_{i+1,1}) \cdot \tilde{d_i}
$ \kappa_{i+1}(U_i,U_{i-1}) = (\tilde{d}_i \times e_{i+1,2} ) \cdot e_{i+1,1}
と表される。
ここで使われている$ \tilde{d_i}は次の反射(屈折)方向ベクトル$ \hat{d}_iの非正規化バージョンである。ここでもRootが入るのを防ぐため、こんな感じの非正規方向ベクトルを使っている。
この式を使って各bicentric coordinate $ U_iを$ U_1で展開することで全てのConstraintは$ U_1の2変数関数として落とすこむことができる（式はすごい大変になると思うけど）。
Refrection $ \tilde{d_i}
反射ベクトル$ \hat{d}_iは
$ \hat{d_i} = - 2 (\hat{d}_{i-1} \cdot \hat{n}_i ) \hat{n_i} + \hat{d}_{i-1}
と表されており、Rootを消すため$ n_i^2 \sqrt{d^2_{i-1}}を両辺にかけると$ \tilde{d}_i = n_i^2 \sqrt{d^2_{i-1}}と定義すれば、
$ \tilde{d}_i = -2(d_{i-1} \cdot n_i) n_i + d_{i-1} n_i^2 \tag{20}
という形でRootを含まない方向ベクトルを得ることができる
Refraction $ \tilde{d_i}
屈折ベクトルについては
$ \hat{d}_i = \eta_i' (\hat{d}_{i-1} - (\hat{d}_{i-1} \cdot \hat{n_i}) \hat{n_i} - \sqrt{1 - \eta_i^2(1 - (\hat{d}_{i-1} \cdot \hat{n}_i)^2)}\hat{n}_i
である($ \eta_i'は相対屈折率$ \eta_{i-1} / \eta_i )。これも同様に$ \tilde{d}_i = n_i^2 \sqrt{d^2_{i-1}}とすると、
$ \tilde{d}_i = \eta_i' (d_{i-1}n^2_i - (d_{i-1} \cdot \hat{n_i}) \hat{n_i} - \sqrt{\beta_i}n_i
という式が得られる。
$ \beta_i = n_i^2 d_{i-1}^2 - \eta_i'^2 (n_i^2 d_{i-1}^2 - (d_{i-1} \cdot n_i)^2)
ここで$ \sqrt{\beta_i}というrootが出てしまった。これは式変形ではもうダメらしいので、論文ではこれを有理数近似することで何とかしている。
$ \beta_iは0~1の範囲であり、$ \sqrt{\beta}に対して次のようなrational functionで近似する。
$ \frac{c_{0,i} + c_{1,i} \beta}{d_{0,i} + d_{1,i} \beta}
論文では$ i=0,1,...,5までの係数を事前計算して使っているらしい(supplementalに載ってる?)。数値的には$ 10^{-3}以下の精度で近似できるらしいので大体のケースなら大丈夫とのこと。
これを使えばPolynomal Systemは全て$ U_1の多項式へと展開でき、2変数へと落とし込むことができる(展開たいへんそー)
3.5 Bivariate spacular polynomial
以上の議論で全ての$ U_iはraitional coordinate mappingによって$ U_1で表現することができる。ということでSpecular Polynomialというのは次のように$ U_1に関するConstraint2つに帰着することができる
$ a(U_1) = 0
$ b(U_1) = 0
どうConstraintを取るかについては反射、屈折で分けて指定する必要がある。
1回反射(R)、屈折(T)、２回反射(RR)の場合について論文で例として出されている。
- R
(6)式と(13)式
$ a(U_1) = ((P_1 U_1 - x_0 ) \times (x_2 - x_0) \cdot N_1 U_1 = 0
$ \begin{aligned} a(U_1) = ((P_1 U_1 - x_0) \cdot N_1 U_1)((x_2 - P_1 U_1)\cdot (N_1 U_1 \times e_{1,1})) + \\ ((x_2 - P_1 U_1) \cdot N_1 U_1)((P_1 U_1 - x_0) \cdot (N_1 U_1 \times e_{1,1})) = 0 \end{aligned}
- T
式が大変なので省略、(6)、(10)式が使われている。
- RR
反射が2回生じるとき、(13)式の$ d_{i-1}を式(19)の$ \tilde{d}_{i-1}を置き換え、$ x_2におけるConstraintによってSpecular Polynomialを考えます。
$ a(U_1) = ((P_2 U_2 - P_1 U_1) \times (x_3 - P_1U_1)) \cdot N_2 U_2 =0
$ b(U_1) = (\tilde{d}_1 \cdot n_2)(d_2 \cdot t_2) + (\tilde{d}_1 \cdot t_2) (d_2 \cdot n_2) = 0
Raitonal Coordinateを使うことで最終的にはこのConstraintは$ U_1の式へと展開される
4 efficient polynomial Solver
今回の問題は2変数のPolynomial Systemを解く問題へと帰着された。ここからはシンプルに代数の話に入っていく
今回得られたSpecular Polynomialは持ってる変数が2と結構少ないが、方程式の次数が結構高い(最小でも4次)のでシンプルに解を求めるのはむずい
だけど、 一変数(univariate)であれば高次の方程式をうまく解く方法がいくつか提案されている。
↓
さらに変数を減らして(Variable elimination)１変数に対する方程式を解く形で求めていく
どうやって変数減らすの？
$ U_1 = (u_1,v_1)の2変数多項式である$ a(u_1,v_1),b(u_1,v_1)について、一旦$ u_1の多項式と見なして各係数を$ a_i(v_1),b_i(v_1)という$ v_1の関数で表すとする。ここで$ u_1の最高次数を$ nと表記する。
$ a(u_1,v_1) = \sum_{i=0}^n a_i(u_1) u_1^i \tag{30}
$ b(u_1,v_1) = \sum_{i=0}^n b_i(u_1) u_1^i
ある2つの多項式(1変数)$ f,gに対して共通解があるとき、終結式 Resultantと呼ばれる式$ \mathrm{Res}(f,g)は0になることが知られている(詳細はそっちのメモで書く)。ここで$ a,bが$ u_1の多項式であると考えると、解を持つにはこの終結式$ \mathrm{Res}(a,b)が0にならなくてはならない。 $ \mathrm{Res}(a,b) = 0
この時、各係数は$ v_1の関数として表現されるため終結式というのは$ v_1の関数として表現することができる。これを$ r(v_1)として、
$ \mathrm{Res}(a,b) = r(v_1)
ということになるので、(30)の多項式の終結式は$ u_1を排除した$ v_1の多項式の方程式として現れるということになる!
$ r(v_1) = 0
これが変数を減らせるテクニックである。それで、じゃあ肝心の終結式ってどう求めるのさという話になるが、終結式は Bezout Matrixと呼ばれる行列を使うことで作ることができる。 Bezout Matrixの定義は詳細に書くとして、次数$ nの多項式2つに対して考えられる行列で、今回の場合$ a,bのBezout Matrix $ Bの成分$ b_{ij}を次のように求めることができるとされている。(これは単にBezout Matrixの定義をそのまま使ってるだけ) $ b_{ij} = \sum_{k=0}^{\min(i,n-1-j)} a_{i-k}(v_1) b_{j+1+k}(v_1) - b_{i-k}(v_1) a_{j+1+k}(v_1) \tag{29}
面白い性質としてBezout Matrixの行列式$ \det{B}は元の多項式の終結式に一致することが知られている(というかそう定義されてるっぽい？）
つまり、上記のBezout Matrix$ Bの行列式$ \det Bを求めることで欲しかった終結式$ r(v_1)を求めることができると
いうわけである。
$ \det B(v_i) = r(v_i) = 0 \tag{31}
このような形で$ u_1を排除し、Specular Polynomialを1変数の多項式$ r(v_1)へと変換します（ただし、$ v_1が求まったら$ u_1もそのまま求める必要はある)
- 4. 2
長いことやってきましたが、最終的にやりたいことというの多項式$ r(v_1)の解を求めることになりました。これは(多分式は複雑ですが)単なる1変数の高次多項式でしかないので、色々高速に解く手札があります。
論文では１つの手法ではBisection Solverによって解くことを提案している。Bisection Solverは多分日本語で二分法 Bisection Methodというやつで、多項式の解を正負の判定で求める手法である。 どっちかっていうと計算量の問題になってるのは$ \det Bっぽいです。$ \det Bのラプラシアン展開は1 bounce程度なら大丈夫だけど2bounce以上では計算量的に無理らしい
- ラプラシアン展開? 普通に余因子展開の別名らしい
Laplacianexpansion
これに対してBisection Solverの実行について次のようにやっている。
$ (0,1)を一様に分割し($ v_1は$ (0,1)の範囲であることに留意)、各両端の行列式をガウスの消去法を使って計算
符号が両端で反転している場所だけでBisection Solverを実行する
論文ではこの手法をPiecewise bisenction solverと呼んでるっぽい
論文実装では100個程度の分割と10イテレーションのbisection solverでテストシーンではいい感じになったとのこと（テストシーンだいぶ単純だけどいいのかって感じはあるけど）
Discussion
Piecewise bisection solverでと食うのはいいけど解がクラスタ化されている(解がいっぱいあって細かいみたいな意味?)場合は細かいピースが必要となる。実用的にはそこまで気にする必要はないと述べている（ほんまか？）
代案としては線形化という手法を採用して、固有値問題へと帰着する方法がある。この問題ではQZ分解(?)によって解くことができるとされているが、どうもQZ分解の計算量が$ O(n^6)とえぐい計算量なので大変とのこと。
任意の次数で高速かつ正確なソルバーが欲しいねって言って終わっている。
Pipeline
手法の理論的な説明は以上
実装の流れとしては以下の3つのフェーズで行われる。
Coffient phase
まず、Specular Chain生成し、各トライアングルの頂点と法線を$ u_1,v_1の2変数多項式へと変換して$ a(u_1,v_1),b(u_1,v_1)を取得する。実装上だと各係数の値を求めるって感じになるっぽい(つまり事前に多項式の形自体は手動で求めておく必要がある？）
Eliminate phase
$ a(u_1,v_1),b(u_1,v_1)から$ v_1に関するBezout Matrixを作り、$ v_1に関する一変数多項式を求めて$ v_1の解をPiecewise Bisection Solverで求める。
論文には書かれていないが、多分手順としては
Bezout Matrix$ B(v_1)の生成
$ \lbrack 0,1\rbrackを100分割した区間を用意
各区間$ \lbrack a,b\rbrackで$ B(a),B(b)を求めて、ガウスの消去法で$ \det(B(a)),\det(B(b))を求める。
$ \det(B(a)),\det(B(b))の符号が反対の関係にあったらbisection solverを実行
Bisection Solverの流れは区間を狭める以外は同様に計算していると思う
いくつかの解を求める
という感じなのかと思われる
Path pase
$ v_1の解は求まったとして次はもともとの$ a(u_1,v_1)から$ u_1を求める。これは単純に$ v_1に解を代入して$ a(u_1)|_{v_1} = 0という方程式を解けばよい。（多分ここでもBisection Solverが使われている？この求め方については特に書かれていない）
こうして$ u_1,v_1の解が得られ、解の候補である$ x_1が求められる。$ x_0から$ x_1にパスを飛ばしてSpecularのパスを新たに生成していく。それで遮蔽される場合やパスのConstraintが満たされない場合は解とならないとして廃棄する。(すべての$ u_1の解でこれをやる？）
多分$ x_1以降のパスはBSDFサンプリングに従って進めていき、最後のSpecular面と光源点をつなげた時Constraintが成立するか判定するというような流れになると思われる。
これに合格したらようやくパスの寄与を評価する。
以上がSpecular Polynomialの実装の流れとなる。
Result
Grints Rendering
ラメとか小さい凹凸の光沢をGrintと呼ぶ。問題としてはCausticsとほぼ同等なのでSpecular Renderingの評価でよく使われる。
Path Cutsと比較したところ、Specular Polynomialsは良い結果を出してくれていると述べている（ぼやっとしてるのはBloomが入っているため)
Path Cutsはニュートン法を使用した手法である。基本的にニュートン法で問題となるのはそもそも解に収束しない場合があるというとこ（計算量もそうですが）。Specular Polynomialはニュートン法とは異なり、解をほぼ確実に求められるため、既存の手法では難しい部分もレンダリングが可能であると述べられている。
https://scrapbox.io/files/6663e923f8948b001dda51bd.png
Caustics Rendering
純粋なパストレ、PPG(良く知らない)、SMS、MPG(最近出た学習ベースの方法?)で比較
SMSではニュートン法の影響をもろに受けており、屈折部分では解がうまく求められずアーティファクトが出てしまうという問題があった。（それでもコースティクスがそこそこ出る）
MPGという手法では(あまり知らないのですが）Specular Chainの学習をしてうまいこと回避はしているものの確率の計算が大変みたいで、分散は大きいという特徴があった。
今回の手法ではそれらに比較して、高速かつ正確に出てくれることが画像から読み取れる。テストシーンにおいてはどの手法よりもGTに近い結果となっている。(点光源しか置いていないシーンだけなので、一般的なシーンで使ったらどうなるのか気になる所）
https://scrapbox.io/files/6663e90cac1264001c9ec4f7.png
Validation
今回の手法では多項式に対するSolverをBisectionでやっていたが、Solverをニュートン法で置き換えたらどうなるのかという検証をしている（SMSと同じではないので注意）
結果としてはBisectionalの方が正確な結果になるとのこと(図ではOurって書かれてるのがBisectionという意味)
これに関しては次元が高すぎるせいでニュートン法では解の収束領域がとても小さいためだと考察されている。
SMSのアーティファクトの原因はやっぱりニュートン法によるものっぽいですね
https://scrapbox.io/files/6663edbf23dac6001dc5b358.png
また、Eigenvalue Solverと比較した結果も乗せている。立場としてはEigenvalue Solverはそれらに比較して最も精度よく、ちゃんと解全てを出してくれる一番精度がいい物である(計算量はバカ高いけど)。対称的にBisectionは分割数によって精度が変わるため、場合によっては求められない解も存在する。
ただ、EigenvalueとBisectionを比較しても実用的にはほぼ問題がないぐらいには見た目が変わらないとのこと。なのでBisectionで十分だよねって結論になっている。
Bloomで強調されてるってのがあるのでちょっとGrintの差が大きいですが、Grint以外はほとんど変わらない感じですね
https://scrapbox.io/files/6663f03a7c9be7001c81b769.png
Bisectionの分割数はどのぐらいいのか、という話では次のような図で議論している。下の図はBisectionの分割数で比較したもの、100 pieceぐらいと1000 pieceではほとんど変わらないので100pieceで十分とのこと
https://scrapbox.io/files/6663ef84cb2944001c080a9a.png
Performance
ニュートン法(SMSのことなのかはちょっとわからない）とBisection Solverでパフォーマンスの比較した図が以下の通り、
基本的にCPUだとNewtonより早くなることが示されている。GPUだと10倍ぐらい遅くなっているように見えるが、解を複数個求められるので多分解１つあたりの時間は減っている(のか？)
なんかここだけCPUとGPUがごっちゃだし、比較対象も全部違うのでなんか分かりづらい。RRでCPU、R,TでのGPUの結果が見たい所。
基本的にはBezout matrixの行列式計算と$ v_1の解の計算がネックとなっている。
https://scrapbox.io/files/6663f0ad7c9be7001c81b9b9.png
Discussion and Limitaions
解の計算について
今回の手法の注目点は解の計算の簡略化であった。
Bisection Solverは実用上問題はないものの、Eigen Solverに比べてどうしても求められない解が存在する。
そもそも高次の多項式を求める問題は難しいのでどうにかしたい
ただ、実験ではあまり下に示すように問題の多項式は低次の項の係数が高く、高次の項はあまり影響してこないと述べている。
タスクに合わせたSolverが今後の課題とのこと
https://scrapbox.io/files/6663f2a7d190ff001d3c357f.png
Long Specular Chain
今回の手法は1~2回反射程度の少ない回数に対して行うことを示していたが、一応いくらでも長いSpecular Chainでも使えないことはない。
だけど、とんでもない次数と組み合わせ爆発により多項式が大変なことになる。今のところ３回反射以上でうまくやれる手法は再構築重点的サンプリング(多分 Path Guidingとかのこと)ぐらいしかないので、それとうまく組み合わせたらいいんじゃないかと述べている。
理論的不偏性
屈折の有理数近似以外の式はUnbiusであるものの、解がすべてちゃんと求まるかはSolverに依存する。今回のBisectionだと理論的には無限回の試行が必要である(もしくはEigenvalue Solver)。実用的にはBiusであるもののBisectionで十分であり、もしUnbiusにしたいなら確率的なアプローチと合わせる必要がある（多分シード値みたいなの用意して初期値を変えるとか）
Superfluous solutions
解の候補となるパスは求められるのはいいのですが、結局そのパスがちゃんと遮蔽とかされないのかというのは実際パスを通してみるまではわからず、場合によっては捨てられるというような状況が十分あり得ます。
Grintのシーンでやってみると大体13%ぐらいしか生き残れなかったとのこと(これはどの手法でも言えることなのでしょうがない気はするけど、結構なくなってしまうと感じる)
Glossyな鏡面に対しての扱い
この方法をGlossyに拡張することは簡単らしい
あまりよくわかっていないが、他の手法もGlossyに対応しているらしいのでそちらを見た方がいいかも？
Surface Presentaiton
今回の手法は完全にTriangleを対象とした手法。
他のサーフェイス表現にも使えるようにしたいと述べている
実装
近日公開とのこと