ガラス・IBL・薄レンズモデル(第7回)

ガラスの実装

ガラスの物理的特性

ガラスの屈折率は本来は光の波長ごとに異なるが、ここでは全波長で同じ屈折率を持つと仮定して話を進めていく。全波長でレンダリングを行うフルスペクトラムレンダラーを実装すれば完全に物理的に正しくなる。

光がガラスに入射すると一部は反射され、残りは屈折してガラス内部を通過していく。この反射の割合を決めるのが フレネルの式(Fresnel Equation) である。

光が屈折すると、光の集光具合が変わるので放射輝度が変化する。

フレネルの式

非常に複雑な式になっている。詳しくは(https://www.wikiwand.com/ja/%E3%83%95%E3%83%AC%E3%83%8D%E3%83%AB%E3%81%AE%E5%BC%8F)

コンピューターグラフィクスでは シュリックの近似式(Schlick's Approximation) というものがよく使われる。光線が法線となす角度を $\theta$ とすると、反射率 $F_r(\theta)$ は

F_r(\theta) = F_0 + (1 - F_0)(1 - \cos{\theta})^5

ここで $F_0$ は光が垂直に入射した場合の反射率であり、空気の屈折率を $n_1$ 、ガラスの屈折率を $n_2$ とすると

F_0 = (\frac{n_1 - n_2}{n_1 + n_2})^2

で与えられる。

材質の屈折率

(http://ww1.tiki.ne.jp/~uri-works/tmp/)で色々な材質の屈折率が見れる。

屈折時の放射輝度の変化

屈折前の放射輝度を $L$ 、屈折後の放射輝度を $L_t$ とすると

L_t = L(\frac{n_2}{n_1})^2

となる。

全反射

屈折率の高い媒質から屈折率の低い媒質に出るときは全反射が起きることがあるので注意する。具体的にはガラスから空気に出る場合には全反射が起きることがある。

屈折方向 $\vec{r}$ は

\vec{r_h} = \frac{n_1}{n_2}\vec{d_h} = \frac{n_1}{n_2}(\vec{d} + (\vec{d}\cdot\vec{n})\vec{n})

として以下の式で与えられるのであった。(第3回宿題)

\vec{r} = \vec{r_h} + \vec{r_p} = \frac{n_1}{n_2}(\vec{d} + (\vec{d}\cdot\vec{n})\vec{n}) -\sqrt{1 - \|\vec{r_h}\|^2}\vec{n}

ルートの中身がマイナスになる場合は屈折方向が存在しないので全反射になる。このときfalseを返すようにすれば屈折と全反射を区別できる。

実装

ガラスの実装の疑似コードを示すと以下にようになる。反射と屈折の選択にはロシアンルーレットを使用している。こうすることで追跡する必要のあるレイは常に一つになる。

double n1, n2; //屈折率
Vec3 hitNormal; //法線

if(空気からガラスにレイが入射している) {
  n1 = 1.0;
  n2 = 1.4;
  hitNormal = ray.hitNormal;
}
else {
  n1 = 1.4;
  n2 = 1.0;
  hitNormal = -ray.hitNormal; //法線を正しい方向にセットする
}

double f = fresnel(ray.direction, hitNormal, n1, n2); //フレネル反射率

//反射
if(rnd() < f) {
  Vec3 refl_direction = reflect(ray.direction, hitNormal); //反射方向
  Ray refl_ray = Ray(ray.hitPos + 0.001*hitNormal, refl_direction); //反射レイ
  return trace(refl_ray, depth + 1); //再帰的に追跡
}
//屈折
else {
  Vec3 refr_direction; //屈折方向
  //屈折した場合
  if(refract(ray.direction, hitNormal, n1, n2, refr_direction)) {
    Ray refr_ray = Ray(ray.hitPos - 0.001*hitNormal, refr_direction); //屈折レイ
    return std::pow(n2/n1, 2.0) * trace(refr_ray, depth + 1);
  }
  //全反射
  else {
    Vec3 refl_direction = reflect(ray.direction, hitNormal); //反射方向
    Ray refl_ray = Ray(ray.hitPos + 0.001*hitNormal, refl_direction); //反射レイ
    return trace(refl_ray, depth + 1);
  }
}

ガラスの外から入射する場合と、ガラスから外に出ていく場合をちゃんと分けて屈折率を変える必要がある。また、ray.hitNormalで取得できる法線は常に物体の外側に向かう方向なので、物体の内側からray.hitNormalを見ると法線の向きが逆になっている。

この事実を使ってレイが物体の外側から入射しているのかを判定することができる。レイの方向を $\vec{r}$ 、ray.hitNormalで取得できる法線を $\vec{n}$ とすると、 $\vec{r} \cdot \vec{n} < 0$ なら外側から入射、 $\vec{r} \cdot \vec{n} > 0$ なら内側から出射している。

bool isInside(const Vec3& r, const Vec3& n) {
  return dot(r, n) > 0;
}

外側か内側かを判定できたら後はフレネルの式より反射率を計算する。rnd()を用いて $[0, 1]$ の一様乱数を発生し、それが反射率より小さかったら反射、そうでなかったら屈折方向を追いかけるようにする。これはまさにロシアンルーレットである。

屈折の場合、全反射が起こる可能性があることに注意する。屈折方向を計算する関数refract()は屈折が起こったときにtrueを返すようにしている。こうすると屈折方向を計算しながら、屈折の場合と全反射の場合をif文を用いて区別することができる。屈折の場合は放射輝度を変化させていることに注意する。

全反射の場合はミラーの実装と同じである。

Image Based Lighting(IBL)

画像を光源として使用する方法を Image Based Rendering(IBL) という。簡単にリアルな見た目のレンダリング画像を生成できるようになるので楽しい。

光源として使用する画像は.jpgや.pngなどの通常の画像ファイルでは駄目で、放射輝度が格納されたHDR画像というものを使用する必要がある。

HDR画像の入手

ここにめっちゃある http://www.hdrlabs.com/sibl/archive.html

HDR画像の読み込み

HDR画像のローダーを自前で実装するのは大変なので、stb_imageという便利なライブラリを使用することにする。シングルファイルなのでこれを必要な場所で#includeするだけで使用できる。使い方はコメントに書いてあるのでそれを参照。

レイの方向から対応するHDR画像の画素を読み込む

sIBLから持ってきたHDR画像をビューワを使って開くと次のように歪んだ画像が表示されるはずだ。

この画像は本来、全天を覆うように球に投影されるように作られているので、それを普通の長方形の画像として表示すると歪んで表示される。

この画像を空の光源として使用するためには、レイの方向に対応する画素を抜き出す必要がある。このためには次のような処理を行えばよい。

レイの方向を球面座標系 $(\phi, \theta)$ に変換する。 $\phi$ は $[0, 2\pi]$ で $\theta$ は $[0, \pi]$ である。
$(u, v)$ に変換する。 $u = \frac{\phi}{2\pi}$ , $v = \frac{\theta}{\pi}$ とする。
画素の座標 $(i, j)$ に変換する。 $i = u*width$ , $j = v*height$ とすればよい。

あとは画素の配列から $(i, j)$ のRGBを読み込むだけである。

stbi_loadfで読み込んだHDR画像は以下のような配列になっている。

R	G	B	R	G	B	...

したがって $(i, j)$ からRを示す配列のインデックス $addr$ を計算するには

addr = 3*i + 3*width*j

とすればよい。すると $(R, G, B)$ は

(R, G, B) = (data[addr], data[addr + 1], data[addr + 2])

となる。

IBLの実装

次のようなクラスを実装するとよいだろう。

class IBL {
  public:
    int width; //横幅
    int height; //縦幅
    float* ibl_data; //HDR画像データ
    
    //ファイル名から読み込む
    IBL(const std::string& filename) {
      int n;
      ibl_data = stbi_loadf(filename.c_str(), &width, &height, &n, 0);
    };
    ~IBL() {
      stbi_image_free(ibl_data);
    };
    
    //レイの方向からの放射輝度を返す
    RGB getColor(const Ray& ray) const {
      //球面座標系(phi, theta)を計算する
      double phi = std::atan2(ray.direction.z, ray.direction.x);
      if(phi < 0) phi += 2*M_PI;
      float theta = std::acos(ray.direction.y);
      
      //(u, v)座標系に直す
      double u = phi/(2*M_PI);
      double v = theta/M_PI;
      
      //画素のインデックスを計算
      int w = (int)(u * width);
      int h = (int)(v * height);
      int adr = 3*w + 3*width*h;
      return RGB(ibl_data[adr], ibl_data[adr + 1], ibl_data[adr + 2]);
    };
}

必要に応じてオフセットを導入して、画像を自由に動かせるようにするとよいだろう。

薄レンズモデル(Thin Lens Model)

今まで使っていたカメラモデルは ピンホールカメラモデル というものであり、全ての物体にピントが合うようになっていた。

イメージセンサーに入射してくる光が必ずピンホールという一点を通過するようになっているので、物体側の点とイメージセンサー上の点が一対一に対応するためである。

しかし、実際のカメラではレンズを用いて光をイメージセンサー上に集光させている。この場合、ピントの合う位置というものが存在し、それ以外の場所では物体側の点がイメージセンサー上の場所に広がりを持って対応するようになるので、ボケ(Bokeh) が発生する。

薄レンズモデル(Thin Lens Model) はごく薄い一枚のレンズを用いたカメラモデルである。本来、レンズでの屈折はレンズ面への入射と出射の両方で起こるが、薄レンズモデルでは一回しか屈折が起こらないと考える。

このような仮定を置くと以下のような性質が得られる。

光軸に平行に入射した光は焦点を通過するように屈折する
焦点を通過して入射してきた光は光軸に平行に出射する
レンズの中心を通る光は屈折せずに直進する

レンズの中心(主点)から焦点までの距離を焦点距離といい、 $f$ で表す。イメージセンサーとレンズ中心までの距離を $a$ で表す。レンズ中心からピントの合う平面までの距離を $b$ で表す。

レンズの式

以下のような関係が成り立つ。

\frac{1}{a} + \frac{1}{b} = \frac{1}{f}

これを使って、 $a$ と $f$ が与えられた場合、レンズ中心からピントの合う平面までの距離は以下の式で求めることができる。

b = \frac{1}{\frac{1}{f} - \frac{1}{a}}

レイの計算方法

ピンホールカメラモデルの場合、イメージセンサーからピンホールに向けてレイを飛ばせば良かったが、薄レンズモデルの場合は、レンズ上の様々な位置を通過して光が入射するので、まずはレンズ上の位置をサンプリングする必要がある。

次にレイの方向を決定する必要があるが、実は屈折の計算をちゃんとやらなくてもレイの方向を決めることができる。

レンズの中心を通る光は方向を変えないで直進するという性質に注目する。直進したレイはレンズ中心から距離 $b$ の位置にある平面と点 $P$ で交わるはずである。

この平面上の点はボケなしでイメージセンサー上に射影されるはずだから、同じイメージセンサー上の点から出たレイは、レンズ上のどの点を通過しても必ずこの点を通るはずである。

これでレイの方向を決定することができた。

点Pを求める計算は次のように行うことができる。カメラの前方向を $\vec{c_f}$ 、イメージセンサーの点 $\vec{I}$ からレンズ中心に向かう方向を $\vec{l}$ とすると、点 $\vec{P}$ は点 $\vec{I}$ から $\vec{l}$ の方向に距離 $\frac{a+b}{\vec{c_f}\cdot\vec{l}}$ だけ離れたところに存在する。よって

\vec{P} = \vec{I} + (\frac{a+b}{\vec{c_f}\cdot\vec{l}})\vec{l}

となる。

まとめると以下のようになる。

レンズ上の点 $\vec{L}$ をサンプリング
点 $\vec{P}$ を計算
点 $\vec{L}$ から方向 $normalize(\vec{P} - \vec{L})$ にレイを飛ばす