A Practical Microcylinder Appearance Model for Cloth Rendering の実装(1)

この記事は レイトレアドベントカレンダー2018 22日目の記事です。

大遅刻しましてすみません!!!

If you like english version, Please check this

ushiostarfish.hatenablog.com

布と繊維

布に加工される繊維は多岐にわたりますが、大きな分類としてな、動物繊維・植物繊維・化学繊維といったものが存在し、 例えば以下のようなものがあります。

動物繊維

名前 原料
羊毛(Wool) 羊の毛
カシミヤ(Cashmere) カシミヤ山羊の毛
絹(Silk) 蚕の繭
モヘア(Mohair) アンゴラ山羊の毛

植物繊維

名前 原料
綿(Cotton) ワタ(アオイ科ワタ属)の木綿(もめん)
リネン(Linen) 亜麻
ラミー(Ramie) 苧麻

化学繊維

名前 原料
ポリエステル(Polyester) 石油等
ナイロン(Nylon) 石炭または石油

これら繊維はどれも1次元的な糸の状態です。布は、これらをなんらかの構造で平面の結合状態を作る必要があります。 例えばフェルトは、羊毛などを絡めてシート状にしたものです。このような布を 不織布(ふしょくふ, Nonwoven fabric) と分類します。

f:id:ushiostarfish:20181223175746j:plain

※ 不織布(拡大図)不織布 - Wikipedia より

これとは対照的に、直交する糸を上下させて作る布を 織布(woven cloth) と呼びます。

f:id:ushiostarfish:20181223180828p:plain

※ Silk Crepe de Chine の拡大写真 IMAN SADEGHI, OLEG BISKER, JOACHIM DE DEKEN, HENRIK WANN JENSEN, UC San Diego, "Practical Microcylinder Appearance Model for Cloth Rendering" より

これら以外では、 編物 または ニット(knit) とも呼ばれるものは、糸を屈曲させてループを作って前後左右に繋いだものを言います。

f:id:ushiostarfish:20181223182744p:plain

信州大学繊維学部, はじめて学ぶ繊維 P.140

この構造から、伸縮性が得られるため、Tシャツや靴下やスポーツウェアにも使われます。

f:id:ushiostarfish:20181223183017p:plain

信州大学繊維学部, はじめて学ぶ繊維 P.140

布のシェーディングモデル

布シェーディングについては古くより研究が多くなされています。古いものとしては、Michael Ashikhmin, Simon Premoze, Peter Shirley, "A microfacet-based brdf generator" ではマイクロファセットベースのサテン・ベルベットのモデル化にも触れられています。 大きなものとしては、Irawan, Marschner らが開発したwoven cloth BRDFのモデル(Piti Irawan, "Appearance of Woven Cloth")があり、こちらはラージスケールとスモールスケールをサポートする解析モデルで、Mitsuba にも実装されているそうです。しかしながら大規模すぎてちょっと自分には攻略できませんでした。

これらはサーフェイスベースのアプローチですが、ボリューメトリックなアプローチもあり、Wenzel Jakob, Adam Arbree, Jonathan T. Moon, Kavita Bala, Steve Marschner, "radiative transfer framework for rendering, materials with anisotropic structure" でのscarfの絵はとても美しく心惹かれるものがあります。

f:id:ushiostarfish:20181223192902p:plain

※ Figure 1 (a) Isotropic scattering. Wenzel Jakob, Adam Arbree, Jonathan T. Moon, Kavita Bala, Steve Marschner, "radiative transfer framework for rendering, materials with anisotropic structure" より

いろいろ目移りしますが、自分はボリュームレンダリングには不慣れなので、今回は織布のレンダリングにフォーカスした論文、IMAN SADEGHI, OLEG BISKER, JOACHIM DE DEKEN, HENRIK WANN JENSEN, UC San Diego, "Practical Microcylinder Appearance Model for Cloth Rendering" に着目します。これは織布にフォーカスしたもので、スモールスケールの糸の構造など無視されているものの、布の異方性反射をよく表現でき、計測したBRDFとよく適合するとあります。絵もなかなか綺麗で、Silk Shot Fabricは印象的です。

f:id:ushiostarfish:20181224143716p:plain

※ Figure 24 IMAN SADEGHI, OLEG BISKER, JOACHIM DE DEKEN, HENRIK WANN JENSEN, UC San Diego, "Practical Microcylinder Appearance Model for Cloth Rendering" より

これはたて(経)糸とよこ(緯)糸に違う色が使われているからで、このモデルの特徴がよく表われていると言えます。

A Practical Microcylinder Appearance Model for Cloth Rendering

f:id:ushiostarfish:20181224143513p:plain

※ Figure 24 IMAN SADEGHI, OLEG BISKER, JOACHIM DE DEKEN, HENRIK WANN JENSEN, UC San Diego, "Practical Microcylinder Appearance Model for Cloth Rendering" より

この論文のアプローチは、以下のような流れになっています。

  1. 布を構成する1本の糸についてのBSDF測定を行う
  2. 個別の糸をうまく表現できるモデルを用意して、パラメーターを調整
  3. 糸を垂直に2本並べたBRDFモデルを作る
  4. 織布の織りパターンに応じて可視性や、糸の角度分布ごとにシェーディングを行う
  5. Shadowing and Masking、および角度分布によるエネルギー再分布(論文ではReweighting)を考慮して重みを調整する

だいたいこんな感じです。

座標系

モデルに先立ってまずは座標系を確認します。論文中 Fig. 9 がわかりやすいです。

f:id:ushiostarfish:20181224144148p:plain

※ Figure 9 IMAN SADEGHI, OLEG BISKER, JOACHIM DE DEKEN, HENRIK WANN JENSEN, UC San Diego, "Practical Microcylinder Appearance Model for Cloth Rendering" より

  •  \theta はスレッド t の 法線面からの角度
  •  \phi は 布文脈での法線からの、法線面における角度でありますが、糸単体では法線は面を形成することには注意が必要です

ただ図にある  \theta_i, \theta_r, \phi_i, \phi_r は、個別の糸BSDFモデルでは直接は使われません。代わりに、

 
\displaystyle {
\phi_{ d }=\phi_{ i }-\phi_{ r }\\ \theta _{ h }=\frac { \theta _{ i }+\theta _{ r } }{ 2 } \\ \theta _{ d }=\frac { \theta _{ i }-\theta _{ r } }{ 2 } 
}

といった量が使われます。 それぞれ意味を確認しておきます。



 
\displaystyle {
\phi_{ d }=\phi_{ i }-\phi_{ r }
}

 \phi_d はそれぞれの入射 \phi_i と出射 \phi_r の差です。



 
\displaystyle {
\theta _{ h }=\frac { \theta _{ i }+\theta _{ r } }{ 2 } 
}

 \theta _h は、入射 \theta _i と出射 \theta _r の中間の角度です。これが0に近いほど、鏡面反射成分が強くなることが想像できます。



 
\displaystyle {
\theta _{ d }=\frac { \theta _{ i }-\theta _{ r } }{ 2 } 
}

 \theta _d はそれぞれの入射 \theta_i と出射 \theta_r の差の半分、言い換えると \theta _h  \theta_i の角度とも言えます。

  \phi_d は半分にされていないのに  \theta _d が半分にされているのは、少し混乱しますね。でも論文の表記がそうなのでそれに揃えます。ちなみに Iman Sadeghi氏のPh.D. Thesis、"Controlling the Appearance of Specular Microstructures" では、5.2 Background Theoryでは両方とも半分になっておらず、統一されています。

糸のモデル

これらを使って、糸のモデルを記述します。糸は、誘導体のラフな円筒であると仮定されており、 糸のモデルは以下のように2つの関数で記述します。

 
\displaystyle {
f_{ s }\left( t,\omega _{ i },\omega _{ r } \right) =\frac { f_{ r,s }\left( t,\omega _{ i },\omega _{ r } \right) +f_{ r,v }\left( t,\omega _{ i },\omega _{ r } \right)  }{ \cos ^{ 2 } \theta _{ d } } 
}

ここで  f_{ r,s } は円筒表面で反射する鏡面反射成分、  f_{ r,s } は円筒内部に入って出てくるボリューム拡散成分になります。では分母にある  \cos ^{ 2 } \theta _{ d } は何か?これを少し考えてみます。

分母の  \cos ^{ 2 } \theta _{ d }

仮に円筒表面が理想的になめらかであり、完全な鏡面であったとします。すると鏡面反射する光は、 \theta_i の角度に対して \theta_r = -\theta_i の角度に鏡面反射します。これは もふもふレンダリング(3)のMpの説明の一部 でも軽く紹介していました。改めてビジュアライズすると、

f:id:ushiostarfish:20181224160335g:plain

このことを踏まえて、拡散関数を以下のようにδ関数を使ってモデリングできます。

 
\displaystyle {
S\left( \phi _{ i },\phi _{ r },\theta _{ i },\theta _{ r } \right) =\frac { N\left( \phi _{ i } \right) \delta \left( \theta _{ i }+\theta _{ r } \right)  }{ \cos ^{ 2 }{ \theta  }  } 
}

ここで、 N \left( \phi_i \right)  \phi の密度関数、 \theta = \theta_i = \theta_r とします。

エネルギー保存について考えてみます。

 
\displaystyle {
\int _{ \Omega  }^{  }{ S\left( \phi _{ i },\phi _{ r },\theta _{ i },\theta _{ r } \right) \cos { \theta _{ i } } d\omega  } \\ =\int _{ -\pi  }^{ \pi  } \int _{ -\frac { \pi \:  }{ 2 }  }^{ \frac { \pi  }{ 2 }  } \: S\left( \phi _{ i },\phi _{ r },\theta _{ i },\theta _{ r } \right) \cos ^{ 2 }{ \theta _{ i } } d\theta d\phi \\ =\int _{ -\pi  }^{ \pi  } \int _{ -\frac { \pi \:  }{ 2 }  }^{ \frac { \pi  }{ 2 }  } \frac { N\left( \phi _{ i } \right) \delta \left( \theta _{ i }+\theta _{ r } \right)  }{ \cos ^{ 2 }{ \theta  }  } \cos ^{ 2 }{ \theta _{ i } } d\theta d\phi \\ =\int _{ -\pi \:  }^{ \pi \: \:  } N\left( \phi  \right) d\phi \int _{ -\frac { \pi \:  }{ 2 }  }^{ \frac { \pi  }{ 2 }  } \frac { \delta \left( \theta _{ i }+\theta _{ r } \right)  }{ \cos ^{ 2 }{ \theta  }  } \cos ^{ 2 }{ \theta _{ i } } d\theta \\ =\int _{ -\frac { \pi \:  }{ 2 }  }^{ \frac { \pi  }{ 2 }  } \frac { \delta \left( \theta _{ i }+\theta _{ r } \right)  }{ \cos ^{ 2 }{ \theta  }  } \cos ^{ 2 }{ \theta _{ i } } d\theta \\ =\frac { \cos ^{ 2 }{ \theta _{ i } }  }{ \cos ^{ 2 }{ \theta  }  } \\ =1
}

ここでの係数  \frac {1}{ \cos ^{ 2 } \theta} はうまく機能しています。この項は必要な係数であるのは明らかで、確かに何かを説明しているようです。 ところで角度θは様々な角度をとりますが、拡散の具合はそれぞれで果たして同じだろうか?というのを考えてみますと、例えば極端な2つのケースを考えて見ます。

まずは  \theta = 0 のケース

f:id:ushiostarfish:20181224164108p:plain

次に  \theta \approx 20° のケース

f:id:ushiostarfish:20181224164410p:plain

明らかに前者は後者に比べて、明らかに光線の密度が高いです。 \theta \frac {\pi}{2} に近づくと、  \frac {1}{ \cos ^{ 2 } \theta} はどんどん大きくなっていきますので、つまりこの  \frac {1}{ \cos ^{ 2 } \theta} 項はこのことを説明していたのではないでしょうか。

δ関数からの拡張

Stephen R. Marschner, Henrik Wann Jensen, Mike Cammarano, "Light Scattering from Human Hair Fibers" では、以上のアイディアをもとにδ関数をラフな表面に拡張しています。それは、

 
\displaystyle {
S\left( \phi _{ i },\phi _{ r },\theta _{ i },\theta _{ r } \right) =\frac { N\left( \phi _{ i } \right) M\left( \theta _{ i } \right)  }{ \cos ^{ 2 }{ \theta _{ d } }  } 
}

というものです。 \theta_i = \theta_r がもはや機能しないため、代わりに  \theta_d が使われているのは、近似的ですが仕方のない側面があります。 そんなわけで、糸のモデルが  \theta_d を含むのもこれを参考にしているからということになります。

鏡面反射成分  f_{ r,s }

鏡面反射成分は、以下のような式で表されます。

 
\displaystyle {
f_{ r,s }\left( t,\omega _{ i },\omega _{ { r } } \right) =F_{ r }\left( \eta ,\vec { w } _{ i } \right) \cos  \left( \frac { \phi _{ d } }{ 2 }  \right) g\left( \gamma _{ s },\theta _{ h } \right) 
}

ここで  F_{ r } はフレネル項、 g は正規化ガウス関数で、平均は0です。フレネルは後ほど議論するとして、正規化ガウス関数は、

 
\displaystyle {
g\left( \sigma ,x \right) =\frac { 1 }{ \sqrt { 2\pi { \sigma  }^{ 2 } }  } exp\left( -\frac { { x }^{ 2 } }{ 2{ \sigma  }^{ 2 } }  \right) 
}

であり、 \theta _h は、入射 \theta _i と出射 \theta _r の中間の角度であり、これが0に近づくと鏡面反射成分が強くなることを表現しています。また、分散を調整することでローブの形を調整できるようになっています。

f:id:ushiostarfish:20181224185718g:plain

※プロットではフレネル項は除いています

 \cos  \left( \frac { \phi _{ d } }{ 2 }  \right) はどこから来るのか?

次に、 \cos  \left( \frac { \phi _{ d } }{ 2 }  \right) の項は特徴的です。プロットしてみると、以下のようになります。

f:id:ushiostarfish:20181225233337p:plain

※ desmos : cos(theta/2)

これは、完全になめらかな円筒に入射した光線の反射が周囲にどのように分布するのかを考えたものになります。これについて実際に光線を飛ばして可視化してみますと、次のようになります。

f:id:ushiostarfish:20181224190158p:plain

入射と出射が近いところの密度は高く、離れるほど密度が小さくなっています。 ではこの密度の関係性は実際にはどうなっているか考えてみます。

f:id:ushiostarfish:20181224194405p:plain

まず図のように \phi と点Pのxとの関係を考えて、これは当然、

 
\displaystyle {
x=\frac { 1 }{ 2 } \sin \left( \phi  \right) 
}

となります。ここで、 \phi の箇所に入射してきた光は、 2\phi のところに出射します。

f:id:ushiostarfish:20181224194500p:plain

したがって、

 
\displaystyle {
2\phi =\phi _{ r }\\ \phi =\frac { 1 }{ 2 } \phi _{ r }
}

つまり、 x \phi_{ r } には以下のような関係性があります。

 
\displaystyle {
x=\frac { 1 }{ 2 } \sin \left( \frac { 1 }{ 2 } \phi _{ r } \right) 
}

これを微分して、

 
\displaystyle {
x=\frac { 1 }{ 2 } \sin { \left( \frac { 1 }{ 2 } \phi _{ r } \right)  } \\ \frac { dx }{ d\phi _{ r } } =\frac { 1 }{ 2 } \frac { 1 }{ 2 } \cos  \left( \frac { 1 }{ 2 } \phi _{ r } \right) \\ \frac { dx }{ d\phi _{ r } } =\frac { 1 }{ 4 } \cos  \left( \frac { 1 }{ 2 } \phi _{ r } \right) 
}

ここで  \frac { dx }{ d\phi_{ r } } とは何かと言えば、微小な反射角度あたりの微小なx、元は単位長さ(=1)の部分が反射角にどのように分配されるかを示しています。 つまり先の節 "δ関数からの拡張" のところのNの一つの形としては、

 
\displaystyle {
N\left( \phi _{ i },\phi _{ r } \right) =\frac { 1 }{ 4 } \cos  \left( \frac { 1 }{ 2 } \phi _{ d } \right) 
}

を使うことができることも意味します。 実際、密度関数なら、1回転分の積分が1になるはずで、

 
\displaystyle {
\int _{ -\pi  }^{ \pi  }{ N\left( \phi _{ i },\phi _{ r } \right) d\phi _{ i } } =\int _{ -\pi  }^{ \pi  }{ \frac { 1 }{ 4 } \cos  \left( \frac { 1 }{ 2 } \phi  \right) d\phi  } 
}

は、

 
\displaystyle {
\frac { d }{ d\phi  } \left( \frac { 1 }{ 2 } \sin { \left( \frac { 1 }{ 2 } \phi  \right)  }  \right) =\frac { 1 }{ 4 } \cos  \left( \frac { 1 }{ 2 } \phi  \right) 
}

から、

 
\displaystyle {
\int _{ -\pi  }^{ \pi  }{ \frac { 1 }{ 4 } \cos  \left( \frac { 1 }{ 2 } \phi  \right) d\phi  } \\ =\frac { 1 }{ 2 } \sin { \left( \frac { 1 }{ 2 } \pi  \right)  } -\frac { 1 }{ 2 } \sin { \left( -\frac { 1 }{ 2 } \pi  \right)  } \\ =\frac { 1 }{ 2 } +\frac { 1 }{ 2 } \\ =1
}

と検算で確かめることもできました。 というわけで、鏡面反射成分の式、

 
\displaystyle {
f_{ r,s }\left( t,\omega _{ i },\omega _{ { r } } \right) =F_{ r }\left( \eta ,\vec { w } _{ i } \right) \cos  \left( \frac { \phi _{ d } }{ 2 }  \right) g\left( \gamma _{ s },\theta _{ h } \right) 
}

 \cos  \left( \frac { \phi _{ d } }{ 2 }  \right) はここから来ている、ということがわかりました。

消えた  \frac { 1 }{ 4 }

ところが、この式には  \frac { 1 }{ 4 } がどこにもありません! 実際、Eugene d’Eon, Guillaume Francois, Martin Hill, Joe Letteri, Jean-Marie Aubry, "An Energy-Conserving Hair Reflectance Model" には、

 
\displaystyle {
N\left( \phi _{ i },\phi _{ r } \right) =\frac { 1 }{ 4 } \cos  \left( \frac { 1 }{ 2 } \phi _{ d } \right) \\ M\left( \theta _{ i },\theta _{ r } \right) =g\left( \beta ,\theta _{ h } \right) \\ S\left( \phi _{ i },\phi _{ r },\theta _{ i },\theta _{ r } \right) =\frac { N\left( \phi _{ d } \right) M\left( \theta _{ h } \right)  }{ \cos ^{ 2 }{ \theta _{ d } }  } 
}

のとき、

 
\displaystyle {
\int _{ -\pi  }^{ \pi  } \int _{ -\frac { \pi \:  }{ 2 }  }^{ \frac { \pi  }{ 2 }  } \: S\left( \phi _{ i },\phi _{ r },\theta _{ i },\theta _{ r } \right) \cos ^{ 2 }{ \theta _{ i } } d\theta d\phi 
}

がいくつかの \beta で、 \theta_i を変えながら観察したとき、どうなるのか?という検証がFigure 3にのっています。

f:id:ushiostarfish:20181229114321p:plain

※ Figure 3 Eugene d’Eon, Guillaume Francois, Martin Hill, Joe Letteri, Jean-Marie Aubry, "An Energy-Conserving Hair Reflectance Model"

今回は、自分でも検証のために以下の実装でプロットしてみました。

これをグラフにすると、

f:id:ushiostarfish:20181224231933p:plain

 \beta の単位は度で、横軸は出射角度(ラジアン)、縦軸は出射放射輝度

Total reflectance from cylinders - Google スプレッドシート

これは確かに Eugene d’Eon, Guillaume Francois, Martin Hill, Joe Letteri, Jean-Marie Aubry, "An Energy-Conserving Hair Reflectance Model" の図と一致しており、かなり部分でエネルギーが増幅しているのがわかります。 そしてこれが球面ガウス関数を加工して作る "4.2. An Energy-Conserving Longitudinal Scattering Function Mp" の導入の大きな動機になっています。

つまりここから考えると、本論文の鏡面反射成分  f_{ r,s } はこれの結果のさらに4倍ものエネルギーを増やすポテンシャルがあることになります。 しかしながらこの部分が間違った記述というわけではどうやらなく、こういうもののようです。IMAN SADEGHI, OLEG BISKER, JOACHIM DE DEKEN, HENRIK WANN JENSEN, UC San Diego, "Practical Microcylinder Appearance Model for Cloth Rendering" にはそもそもエネルギー保存に関する記述は無く、この論文の興味の対象では無いのでしょう。ここでも考察はこのくらいにして、次に進みます。

フレネル項

フレネル項は、以下のように定義されています。

 
\displaystyle {
F_{ r }\left( \eta ,\vec { w } _{ i } \right) =F_{ r }\left( \eta ,\cos ^{ -1 }{ \left( \cos { \left( \theta _{ d } \right)  } \cos { \left( \frac { \phi _{ d } }{ 2 }  \right)  }  \right)  }  \right) 
}

この項がどこから来たのか、ちょっとわかりませんでした。しかし、 \theta_{ d }  \frac {\phi _{ d } }{ 2 } も、ハーフアングルまでの角度であり、そのコサインの積ということは、どちらか一方でも0に近づけば、グレージング角で入射したことにする、というような意図は感じます。論文にはマイクロファセットを考慮したフレネル項の導入を検討したが、測定したデータとのマッチングに改善が見られなかった、とありました。 論文には記載がありませんでしたが、今回は誘導体フレネル項の計算には、Bruce Walter, Stephen R. Marschner, Hongsong Li, Kenneth E. Torrance, "Microfacet Models for Refraction through Rough Surfaces" の中で、"A convenient exact formulation for dielectrics with unpolarized light" として書かれている以下のものを使うことにしました。

 
\displaystyle {
{ F\left( \eta _{ t },\eta _{ i },\theta  \right) =\frac { 1 }{ 2 } \frac { \left( g-c \right) ^{ { 2 } } }{ \left( g+c \right) ^{ { 2 } } } \left( 1+\frac { \left( c\left( g+c \right) -1 \right) ^{ { 2 } } }{ \left( c\left( g+c \right) +1 \right) ^{ { 2 } } }  \right)  }\\ g=\sqrt { \frac { { \eta _{ t } }^{ 2 } }{ { \eta _{ i } }^{ 2 } } -1+{ c }^{ 2 } } \\ c=\cos { \theta  } 
}

ボリューム拡散成分  f_{ r,s }

続いて f_{ r,s } です。こちらは以下のように定義されています。

 
\displaystyle {
f_{ r,v }\left( t,\omega _{ i },\omega _{ r } \right) =F\frac { \left( 1-k_{ d } \right) g\left( \gamma _{ v },\theta _{ h } \right) +k_{ d } }{ \cos { \theta _{ i } } +\cos { \theta _{ r } }  } A
}

ここでAは波長依存の吸収係数、ここで色がつきます。この項は鏡面反射成分にはありませんでした。

F 項

Fはフレネル項で論文中では、

 
\displaystyle {
F=F_{ t }\left( \eta ,\vec { \omega _{ i } }  \right) F_{ t }\left( \eta ',\vec { \omega '_{ r } }  \right) 
}

となっています。しかしこれ以上の言及はなかったのですが、おそらく、

 
\displaystyle {
F_{ t }\left( \eta ,\vec { \omega _{ i } }  \right) =1-F_{ r }\left( \eta ,\vec { w } _{ i } \right) 
}

であろうと思いつつ、

 
\displaystyle {
F_{ t }\left( \eta ',\vec { \omega '_{ r } }  \right) 
}

の項は意図がわかりません。とはいえ、シリンダー内部の多重拡散もこの  f_{ r,s } で説明できると仮定すると、

 
\displaystyle {
F=F_{ t }\left( \eta ,\vec { \omega _{ i } }  \right) 
}

としてしまうのもそれほど間違っているとは思えませんが、今回は、

 
\displaystyle {
F={ { F_{ t } }^{ 2 }\left( \eta ,\vec { \omega _{ i } }  \right)  }
}

で実装してみることにしました。

 g\left( \gamma _{ v },\theta _{ h } \right)

次に、 g\left( \gamma _{ v },\theta _{ h } \right) 、この項は分散のパラメータが違うだけで、鏡面反射成分の中にある  g\left( \gamma _{ s },\theta _{ h } \right) と全く同じ働きをします。つまり反射成分だけでなく、一度屈折した光線も出るときにも  \theta に関しては同じふるまいだとモデル化されています。これは真なのでしょうか? これを観察するために、なめらかな誘導体シリンダーへ入射する光線の屈折を追跡してみます。(Houdini vex: tracking refraction · GitHub

f:id:ushiostarfish:20181229123015g:plain

ここで入射光線は \theta = 45° であり、光線の  \sin \theta をプロットしてみると、

f:id:ushiostarfish:20181229123735p:plain

なんと、すべての出射ベクトルは、入射と全く同じ角度で出射しています。ちなみにこれは入射角度を変えても同様のふるまいをします。 実はこのことはすでに、Stephen R. Marschner, Henrik Wann Jensen, Mike Cammarano, "Light Scattering from Human Hair Fibers" の "A Loci of reflections" で触れられています。 "A Loci of reflections" の説明では、誘導体シリンダーの法線を、球体の赤道部分の法線と捉えることでこれを説明しています。

f:id:ushiostarfish:20181229152111p:plain

※ n, n'がある部分がシリンダー法線です。

ここで、球面上の点  v_i を入射ベクトルと考えます。とするとこの入射ベクトル v_i は、赤道上のいろいろな法線と衝突する可能性があります(裏面、つまり v_i n が反対向きでないケースに限りますが)。ここではまず代表的な法線、  v_i の真下にある法線を  n と表記します。

この v_i n の 組み合わせで屈折の振る舞いを考えて、ある相対屈折率では、屈折後は図の  v_t のベクトルになりました。スネルの法則によれば、

 
\displaystyle {
\eta _{ i }\sin  \theta _{ i }=\eta _{ t }\sin  \theta _{ t }
}

の関係性があり、ここで、 \eta =\frac { \eta _{ t } }{ \eta _{ i } } とすると、

 
\displaystyle {
\sin  \theta _{ i }=\eta \sin  \theta _{ t }
}

ここで  \sin\theta_i, \sin\theta_t は何か?と考えてみると、幾何学的にはそれぞれ図の  h_i, h_t だと捉えることができます。 つまり、

 
\displaystyle {
\left| h_{ i } \right| =\eta \left| h_{ t } \right| 
}

という関係性を考える事ができます。 しかし法線は実際には  n だけでなく、赤道の上のどんなベクトルにもなりえます。そこで真下の法線ではない法線  n' について考えてみます。下の図をご覧ください。ここでは、 n' のケースの  \sin\theta にあたる部分、 h'_{ i }, h'_{ t } を足しています。

f:id:ushiostarfish:20181229173945p:plain

このとき、 \left| h'_{ t } \right| を Y 軸に投射した量、

 
\displaystyle {
H=\left| h'_{ t } \right| \left| \frac { h'_{ t } }{ \left| h'_{ t } \right|  } \cdot Y \right| 
}

はどんな値になるのか考えてみます。 まず、やはりスネルの法則によれば、

 
\displaystyle {
\left| h'_{ i } \right| =\eta \left| h'_{ t } \right| 
}

であるので、

 
\displaystyle {
H=\frac { \left| h'_{ i } \right|  }{ \eta  } \left| \frac { h'_{ t } }{ \left| h'_{ t } \right|  } \cdot Y \right| 
}

次に、

 
\displaystyle {
\left| h'_{ i } \right| \left| \frac { h'_{ i } }{ \left| h'_{ i } \right|  } \cdot Y \right| =\left| h_{ i } \right| \\ \left| h'_{ i } \right| ={ \left| \frac { h'_{ i } }{ \left| h'_{ i } \right|  } \cdot Y \right|  }^{ -1 }\left| h_{ i } \right| 
}

であることから、

 
\displaystyle {
H=\frac { 1 }{ \eta  } { \left| \frac { h'_{ i } }{ \left| h'_{ i } \right|  } \cdot Y \right|  }^{ -1 }\left| h_{ i } \right| \left| \frac { h'_{ t } }{ \left| h'_{ t } \right|  } \cdot Y \right| \\ H=\frac { \left| h_{ i } \right|  }{ \eta  } { \left| \frac { h'_{ i } }{ \left| h'_{ i } \right|  } \cdot Y \right|  }^{ -1 }\left| \frac { h'_{ t } }{ \left| h'_{ t } \right|  } \cdot Y \right| 
}

このとき、ベクトル  h'_i と ベクトル  h'_t は、向きがちょうど正反対になります。したがって、

 
\displaystyle {
\left| \frac { h'_{ i } }{ \left| h'_{ i } \right|  } \cdot Y \right| =\left| \frac { h'_{ t } }{ \left| h'_{ t } \right|  } \cdot Y \right| 
}

でなければばならないでしょう。これを使うと、

 
\displaystyle {
H=\frac { \left| h_{ i } \right|  }{ \eta  } 
}

おっと、 \left| h_{ i } \right| =\eta \left| h_{ t } \right| であったことを思い出すと、

 
\displaystyle {
H=\frac { \left| h_{ i } \right|  }{ \eta  } =\frac { \eta \left| h_{ t } \right|  }{ \eta  } =\left| h_{ t } \right| 
}

これの意味するところは、どのような  n' であっても、屈折後のベクトルの赤道面からの垂直距離は屈折率にのみ支配されるということになります。加えて、誘導体シリンダーから入った後に出ることを考えると、屈折率は打ち消し合い、出るときには入るときの  \theta と同じ  \theta で出射する ということもまた言えることになります(正負については反転しますが)。

したがって、ボリューム拡散成分の g\left( \gamma _{ v },\theta _{ h } \right) 項はこの意味で妥当だと考えることができるでしょう。ただし、ラフな誘導体円筒でモデル化しているとはいえ、内部に入り込んだ光は中で拡散を繰り返すことが考えられ、今回のモデルでは、ボリューム拡散成分のgの分散は、鏡面反射とは違う分散を与えられるようになっており、元論文Table Ⅱに記載されているパラメーターも、すべて \gamma_v \gamma_s よりも大きな分散を与えるようになっています。

分母の  \cos { \theta _{ i } } +\cos { \theta _{ r } }

論文では、半無限媒質の多重散乱による拡散反射を考えたもの(Subrahmanyan Chandrasekhar, "Radiative Transfer")だとありますが、すみませんこの項についてはよくわかりませんでした。ただ、この項により測定値とよりマッチするようになったとのことです。

 k_d

 k_d 項 はかなりシンプルで、綿やリネンなどの等方的な拡散を表すために、 g\left( \gamma _{ v },\theta _{ h } \right)  1 k_d で線形補間しています。

以上をプロットしてみると以下のようになります。以上の理解からすると、おおむね直観的な動きだと言えます。

f:id:ushiostarfish:20181229172951g:plain

これで個別の糸のモデル化が完了しました。次回はこれを使ってシェーディングモデルを構築します。

ushiostarfish.hatenablog.com