いろいろな形状の表面に点を一様分布させる

February 22, 2022      

やりたいこと

輻射による熱のやり取りを評価する場合,レイトレースすることで物体間のRadiative Coupling(または形態係数)を求めることが出来る. このとき光線は,物体の表面上で面積に関して一様になるように発生点を指定して,そこからランベルトの余弦則に沿うように光線の方向を決めてやる必要がある. 今回は基本的な表面形状(長方形・三角形・円板・球面・円柱・円錐・放物面)について,0–1の範囲のランダム値からどうやって一様分布を発生させるかを考える.

長方形上の一様分布

下図のような長方形上に一様分布させる場合は,2つのランダム値(q1,q2)(q_1, q_2)を発生させて,以下のように計算すればよい.

q1(p2p1)+q2(p3p1)\begin{equation} q_1 (\vec{p_2} - \vec{p_1}) + q_2 (\vec{p_3} - \vec{p_1}) \end{equation}

三角形上の一様分布

三角形の場合は,発生させた2つのランダム値(q1,q2)(q_1, q_2)の和が1より大きい場合は値を捨てて,和が1以下の組み合わせが得られるまでランダム値を繰り返し発生させる.

q1(p2p1)+q2(p3p1)\begin{equation} q_1 (\vec{p_2} - \vec{p_1}) + q_2 (\vec{p_3} - \vec{p_1}) \end{equation}

円柱側面上の一様分布

円柱側面上に点を一様分布させる場合は,高方向・周方向それぞれ一様に値をとればよい. 高さhhは,0–1の範囲のランダム値qqを用いて以下のように決められる.

h=hmaxq\begin{gather} h = h_\mathrm{max} q \end{gather}

周方向のパラメタφ\varphiは以下のように決められる.

φ=2πq\begin{gather} \varphi = 2\pi q \end{gather}

円板上の一様分布

部分円板上に点を一様にばらまく場合,ある点が半径方向の位置rrの内側に入る確率は次のように表される.

q=rinnerrθstartθendrdrdθrinnerrouterθstartθendrdrdθ=r2rinner2router2rinner2\begin{equation} q = \frac{\int_{r_\mathrm{inner}}^{r} \int_{\theta_\mathrm{start}}^{\theta_\mathrm{end}} r dr d\theta}{\int_{r_\mathrm{inner}}^{r_\mathrm{outer}} \int_{\theta_\mathrm{start}}^{\theta_\mathrm{end}} r dr d\theta} = \frac{r^2 - r_\mathrm{inner}^2}{r_\mathrm{outer}^2 - r_\mathrm{inner}^2} \end{equation}

この式を逆に解いてやると,0–1の値ppに対して半径方向の位置rrを対応させることが出来る.

r=q(router2rinner2)+rinner2\begin{equation} r = \sqrt{q(r_\mathrm{outer}^2 - r_\mathrm{inner}^2) + r_\mathrm{inner}^2} \end{equation}

周方向については,次のように決められる.

θ=q(θendθstart)+θstart\begin{equation} \theta = q (\theta_\mathrm{end} - \theta_\mathrm{start}) + \theta_\mathrm{start} \end{equation}

球面上の一様分布

部分球面上に一様に点をばらまく場合,ある点が緯度方向についてθ\theta以下である確率は次のように表される.

q=θapexθφstartφendsinθdφdθθbaseθendφstartφendsinθdφdθ=[cosθ]θapexθ[cosθ]θapexθbase=cosθapexcosθcosθapexcosθbasewhere  θbase=arccoshbaser,  θapex=arccoshapexr\begin{gather} q = \frac{\int_{\theta_\mathrm{apex}}^{\theta} \int_{\varphi_\mathrm{start}}^{\varphi_\mathrm{end}} \sin \theta d\varphi d\theta}{\int_{\theta_\mathrm{base}}^{\theta_\mathrm{end}} \int_{\varphi_\mathrm{start}}^{\varphi_\mathrm{end}} \sin \theta d\varphi d\theta } = \frac{[-\cos \theta]_{\theta_\mathrm{apex}}^{\theta}}{[-\cos \theta]_{\theta_\mathrm{apex}}^{\theta_\mathrm{base}}} = \frac{\cos \theta_\mathrm{apex}-\cos \theta}{\cos \theta_\mathrm{apex}-\cos \theta_\mathrm{base}} \\ \mathrm{where}~~ \theta_\mathrm{base} = \arccos \frac{h_\mathrm{base}}{r}, ~~ \theta_\mathrm{apex} = \arccos \frac{h_\mathrm{apex}}{r} \end{gather}

この式を逆に解いてやると,0–1の値ppに対して緯度方向のパラメタθ\thetaを対応させることが出来る.

θ=arccos(cosθapexq(cosθapexcosθbase))=arccos((1q)hapexr+qhbaser)\begin{align} \theta &= \arccos (\cos \theta_\mathrm{apex} - q(\cos \theta_\mathrm{apex}-\cos \theta_\mathrm{base})) \\ &= \arccos\left( (1-q) \frac{h_\mathrm{apex}}{r} + q \frac{h_\mathrm{base}}{r} \right) \end{align}

経度方向については,次のように決められる.

φ=q(φendφstart)+φstart\begin{equation} \varphi = q (\varphi_\mathrm{end} - \varphi_\mathrm{start}) + \varphi_\mathrm{start} \end{equation}

円錐面上の一様分布

部分円錐面上に一様に点をばらまく場合,ある点が高さhh以下である確率は次のように表される.

q=0hφstartφend2πrcosθdφdh0hmaxφstartφend2πrcosθdφdh,  where  r=r1tanθ h,  tanθ=r1r2hmax=0hr1tanθ h dh0hmaxr1tanθh dh=[r1htanθ2h2]0h[r1htanθ2h2]0hmax=r1htanθ2h2r1hmaxtanθ2hmax2\begin{align} q &= \frac{\int_{0}^{h} \int_{\varphi_\mathrm{start}}^{\varphi_\mathrm{end}} \frac{2\pi r}{\cos \theta}d\varphi dh}{\int^{h_\mathrm{max}}_{0} \int_{\varphi_\mathrm{start}}^{\varphi_\mathrm{end}} \frac{2\pi r}{\cos \theta}d\varphi dh}, ~~ \mathrm{where}~~ r = r_1 - \tan \theta ~h, ~~ \tan \theta = \frac{r_1 - r_2}{h_\mathrm{max}} \\ &= \frac{\int_{0}^{h} r_1 - \tan \theta~ h ~ dh}{\int^{h_\mathrm{max}}_{0} r_1 - \tan \theta h ~ dh} = \frac{\left[r_1 h - \frac{\tan \theta}{2}h^2 \right]^h_0}{\left[r_1 h - \frac{\tan \theta}{2}h^2 \right]^{h_\mathrm{max}}_0} = \frac{r_1 h - \frac{\tan \theta}{2}h^2}{r_1 h_\mathrm{max} - \frac{\tan \theta}{2}h^2_\mathrm{max}} \end{align}

この式を逆に解いてやると,0–1の値qqに対して高さ方向のパラメタhhを対応させることが出来る. 解は2つ存在するが,高さhhが円錐の頂点を超えないようにしたいので,マイナスの場合を用いる.

h=r1tanθ±(r1tanθ)22qhmaxr1tanθ+qhmax2\begin{equation} h = \frac{r_1}{\tan \theta} \pm \sqrt{\left( \frac{r_1}{\tan \theta} \right)^2 - 2q h_\mathrm{max} \frac{r_1}{\tan \theta} + q h_\mathrm{max}^2} \end{equation}

周方向については,次のように決められる.

φ=q(φendφstart)+φstart\begin{equation} \varphi = q (\varphi_\mathrm{end} - \varphi_\mathrm{start}) + \varphi_\mathrm{start} \end{equation}

放物面上の一様分布

下図のような部分放物面について考えよう.ただし放物面の頂点は原点にあり,放物面の軸が高さ方向の軸に一致するものとする.半径rrと高さhhの関係は以下のように表される.

h=a2r2,  where  a2=hmaxrmax2\begin{equation} h = a^2 r^2, ~~\mathrm{where} ~~ a^2 = \frac{h_\mathrm{max}}{r_\mathrm{max}^2} \end{equation}
r=ha\begin{equation} r = \frac{\sqrt{h}}{a} \end{equation}

接線の傾きをtanθ\tan \thetaで表すとすると,以下のように表される.

tanθ=dhdr=2a2r\begin{equation} \tan{\theta} = \frac{dh}{dr} = 2a^2 r \end{equation}

1+1tan2θ=1sin2θ1+\frac{1}{\tan^2 \theta} = \frac{1}{\sin^2 \theta}を用いると,以下の関係が得られる.

1sinθ=1+1tan2θ=1+14a4r2rsinθ=ha2+h4a6r2=ha2+14a4\begin{gather} \frac{1}{\sin \theta} = \sqrt{1 + \frac{1}{\tan^2 \theta}} = \sqrt{1 + \frac{1}{4a^4r^2}} \\ \frac{r}{\sin \theta} = \sqrt{\frac{h}{a^2} + \frac{h}{4a^6r^2}} = \sqrt{\frac{h}{a^2} + \frac{1}{4a^4}} \end{gather}

一様に点をばらまく場合,ある点が高さhh以下である確率は次のように表される.

q=hminhφstartφendrsinθdφdhhminhmaxφstartφendrsinθdφdh=hminhh+14a2dhhminhmaxh+14a2dh=[23(h+14a2)32]hminh[23(h+14a2)32]hminhmax=(h+14a2)32(hmin+14a2)32(hmax+14a2)32(hmin+14a2)32\begin{align} q &= \frac{\int_{h_\mathrm{min}}^{h} \int_{\varphi_\mathrm{start}}^{\varphi_\mathrm{end}} \frac{r}{\sin \theta}d\varphi dh}{\int^{h_\mathrm{max}}_{h_\mathrm{min}} \int_{\varphi_\mathrm{start}}^{\varphi_\mathrm{end}} \frac{r}{\sin \theta}d\varphi dh} = \frac{\int_{h_\mathrm{min}}^{h} \sqrt{h + \frac{1}{4a^2}} dh}{\int^{h_\mathrm{max}}_{h_\mathrm{min}} \sqrt{h + \frac{1}{4a^2}} dh} \\ &= \frac{\left[ \frac{2}{3}\left(h + \frac{1}{4a^2} \right)^{\frac{3}{2}} \right]^h_{h_\mathrm{min}}}{\left[ \frac{2}{3} \left(h + \frac{1}{4a^2} \right)^{\frac{3}{2}} \right]^{h_\mathrm{max}}_{h_\mathrm{min}}} = \frac{\left(h + \frac{1}{4a^2} \right)^{\frac{3}{2}} - \left(h_\mathrm{min} + \frac{1}{4a^2} \right)^{\frac{3}{2}}}{\left(h_\mathrm{max} + \frac{1}{4a^2} \right)^{\frac{3}{2}} - \left(h_\mathrm{min} + \frac{1}{4a^2} \right)^{\frac{3}{2}}} \end{align}

この式を逆に解いてやると,0–1の値qqに対して高さ方向のパラメタhhを対応させることが出来る.

(h+14a2)32=q{(hmax+14a2)32(hmin+14a2)32}+(hmin+14a2)32h=[q{(hmax+14a2)32(hmin+14a2)32}+(hmin+14a2)32]2314a2\begin{gather} \left(h + \frac{1}{4a^2} \right)^{\frac{3}{2}} = q \left\{ \left(h_\mathrm{max} + \frac{1}{4a^2} \right)^{\frac{3}{2}} - \left(h_\mathrm{min} + \frac{1}{4a^2} \right)^{\frac{3}{2}} \right\} + \left(h_\mathrm{min} + \frac{1}{4a^2} \right)^\frac{3}{2} \\ h = \left[q \left\{ \left(h_\mathrm{max} + \frac{1}{4a^2} \right)^{\frac{3}{2}} - \left(h_\mathrm{min} + \frac{1}{4a^2} \right)^{\frac{3}{2}} \right\} + \left(h_\mathrm{min} + \frac{1}{4a^2} \right)^\frac{3}{2} \right]^\frac{2}{3} - \frac{1}{4a^2} \end{gather}