PDF の管理方法の整理

今の ray_color() 関数には二つの PDF がハードコードされている:

  1. ライトの形状が関係する p0()
  2. 法線と物体表面の種類が関係する p1()

p0() についてはライト (および多くサンプルしたい hittable) に関する情報を ray_color() 関数に渡し、p1() については material に PDF を尋ねるようにできる。さらに鏡面ベクトルの存在は hit 関数か material クラスが答えるようにする。

拡散反射と鏡面反射

ぜひ表現できるようにしておきたいマテリアルの一つが、ニスを塗った木材のような部分的には理想的な鏡面反射をして部分的には拡散反射をするマテリアルである。マテリアルで鏡面反射レイと拡散反射レイの二つを同時に生成するレンダラもある。ただ私は場合分けが好きではないので、マテリアルで反射と拡散をランダムに選ぶようにする。このアプローチでは PDF を慎重に計算しなければならず、さらに拡散と鏡面に応じて ray_color() の処理を切り替える必要がある。幸い拡散の場合には pdf_value() を呼べば PDF を計算できるので、特別な処理はいらない。

matrial の設計を変更して、関数の引数を struct にまとめるとしよう。hittable の hit_record と同様である:

struct scatter_record {
  ray specular_ray;
  bool is_specular;
  color attenuation;
  shared_ptr<pdf> pdf_ptr;
};

class material {
public:
  virtual ~material() {}

  virtual color emitted(
    const ray& r_in, const hit_record& rec, double u, double v, const point3& p
  ) const {
    return color(0,0,0);
  }

  virtual bool scatter(
    const ray& r_in, const hit_record& rec, scatter_record& srec
  ) const {
    return false;
  }

  virtual double scattering_pdf(
    const ray& r_in, const hit_record& rec, const ray& scattered
  ) const {
    return 0;
  }
};
リスト 3.31 [material.h] material クラスのリファクタリング

lambertian マテリアルは簡単になる:

class lambertian : public material {
public:
  lambertian(shared_ptr<texture> a) : albedo(a) {}

  bool scatter(
    const ray& r_in, const hit_record& rec, scatter_record& srec
  ) const {
    srec.is_specular = false;
    srec.attenuation = albedo->value(rec.u, rec.v, rec.p);
    srec.pdf_ptr = make_shared<cosine_pdf>(rec.normal);
    return true;
  }

  double scattering_pdf(
    const ray& r_in, const hit_record& rec, const ray& scattered
  ) const {
    auto cosine = dot(rec.normal, unit_vector(scattered.direction()));
    return cosine < 0 ? 0 : cosine/pi;
  }

public:
  shared_ptr<texture> albedo;
};
リスト 3.32 [material.h] 新しい lambertian::scatter() メソッド

ray_color() にも変更が必要になる。

color ray_color(
  const ray& r,
  const color& background,
  const hittable& world,
  shared_ptr<hittable> lights,
  int depth
) {
  hit_record rec;

  // 反射回数が一定よりも多くなったら、その時点で追跡をやめる
  if (depth <= 0)
    return color(0,0,0);

  // レイがどのオブジェクトとも交わらないなら、背景色を返す
  if (!world.hit(r, 0.001, infinity, rec))
    return background;

  scatter_record srec;
  color emitted = rec.mat_ptr->emitted(r, rec, rec.u, rec.v, rec.p);
  if (!rec.mat_ptr->scatter(r, rec, srec))
    return emitted;

  auto light_ptr = make_shared<hittable_pdf>(lights, rec.p);
  mixture_pdf p(light_ptr, srec.pdf_ptr);

  ray scattered = ray(rec.p, p.generate(), r.time());
  auto pdf_val = p.value(scattered.direction());

  return emitted
    + srec.attenuation * rec.mat_ptr->scattering_pdf(r, rec, scattered)
               * ray_color(scattered, background, world, lights, depth-1)
               / pdf_val;
}
リスト 3.33 [main.cc] 拡散反射だけを計算する ray_color

ライトに関する情報は main で直接与える:

int main() {
  ...
  camera camera;
  hittable_list world = cornell_box(camera, aspect_ratio);
  shared_ptr<hittable> light =
    make_shared<xz_rect>(213, 343, 227, 332, 554, nullptr);

  for (int j = image_height-1; j >= 0; --j) {
    std::cerr << "\rScanlines remaining: " << j << ' ' << std::flush;
    for (int i = 0; i < image_width; ++i) {
      color pixel_color(0, 0, 0);
      for (int s = 0; s < samples_per_pixel; ++s) {
        auto u = (i + random_double()) / (image_width-1);
        auto v = (j + random_double()) / (image_height-1);
        ray r = camera.get_ray(u, v);
        pixel_color += ray_color(r, background, world, light, max_depth);
      }
      write_color(std::cout, pixel_color, samples_per_pixel);
    }
  }
}
リスト 3.33a [main.cc] 変更後の main

鏡面反射の計算

鏡面反射をするマテリアルは今まで放っておいてきた。また法線を変化させるインスタンスについても何もしていない。しかし scatter_record を使った設計は全体的に整っているので、いずれも修復すればそのまま使える。ここには鏡面マテリアルを直す方法だけを示す。metal は簡単に直せる:

class metal : public material {
public:
  metal(const color& a, double f) : albedo(a), fuzz(f < 1 ? f : 1) {}

  virtual bool scatter(
    const ray& r_in, const hit_record& rec, scatter_record& srec
  ) const {
    vec3 reflected = reflect(unit_vector(r_in.direction()), rec.normal);
    srec.specular_ray = ray(rec.p, reflected+fuzz*random_in_unit_sphere());
    srec.attenuation = albedo;
    srec.is_specular = true;
    srec.pdf_ptr = 0;
    return true;
  }

public:
  color albedo;
  double fuzz;
};
リスト 3.34 [material.h] metal::scatter 関数

fuzz の値が大きいと、このマテリアルは理想的な鏡面反射とは異なる反射レイを生成する。しかし陰的な (PDF を使わない) サンプリングがあるので、このマテリアル以前と同様に正しく動作する。

ray_color() では陰的にサンプルされたレイを生成するケースを追加するだけでよい:

color ray_color(
  const ray& r,
  const color& background,
  const hittable& world,
  shared_ptr<hittable> lights,
  int depth
) {
  hit_record rec;

  // 反射回数が一定よりも多くなったら、その時点で追跡をやめる
  if (depth <= 0)
    return color(0,0,0);

  // レイがどのオブジェクトとも交わらないなら、背景色を返す
  if (!world.hit(r, 0.001, infinity, rec))
    return background;

  scatter_record srec;
  color emitted = rec.mat_ptr->emitted(r, rec, rec.u, rec.v, rec.p);
  if (!rec.mat_ptr->scatter(r, rec, srec))
    return emitted;

  if (srec.is_specular) {
    return srec.attenuation
       * ray_color(srec.specular_ray, background, world, lights, depth-1);
  }

  shared_ptr<pdf> light_ptr = make_shared<hittable_pdf>(lights, rec.p);
  mixture_pdf p(light_ptr, srec.pdf_ptr);

  ray scattered = ray(rec.p, p.generate(), r.time());
  auto pdf_val = p.value(scattered.direction());
  delete srec.pdf_ptr;

  return emitted + srec.attenuation
           * rec.mat_ptr->scattering_pdf(r, rec, scattered)
           * ray_color(scattered, background, world, lights, depth-1)
           / pdf_val;
}
リスト 3.35 [main.cc] ray_color 関数

コーネルボックスの直方体のマテリアルを metal に変更しよう:

hittable_list cornell_box(camera& cam, double aspect) {
  hittable_list world;

  auto red   = make_shared<lambertian>(make_shared<solid_color>(.65, .05, .05));
  auto white = make_shared<lambertian>(make_shared<solid_color>(.73, .73, .73));
  auto green = make_shared<lambertian>(make_shared<solid_color>(.12, .45, .15));
  auto light = make_shared<diffuse_light>(make_shared<solid_color>(15, 15, 15));

  world.add(make_shared<yz_rect>(0, 555, 0, 555, 555, green));
  world.add(make_shared<yz_rect>(0, 555, 0, 555, 0, red));
  world.add(make_shared<flip_face>(
    make_shared<xz_rect>(213, 343, 227, 332, 554, light))
  );
  world.add(make_shared<xz_rect>(0, 555, 0, 555, 555, white));
  world.add(make_shared<xz_rect>(0, 555, 0, 555, 0, white));
  world.add(make_shared<xy_rect>(0, 555, 0, 555, 555, white));

  shared_ptr<material> aluminum =
    make_shared<metal>(color(0.8, 0.85, 0.88), 0.0);
  shared_ptr<hittable> box1 =
    make_shared<box>(point3(0,0,0), point3(165,330,165), aluminum);
  box1 = make_shared<rotate_y>(box1, 15);
  box1 = make_shared<translate>(box1, vec3(265,0,295));
  world.add(box1);

  shared_ptr<hittable> box2 =
    make_shared<box>(point3(0,0,0), point3(165,165,165), white);
  box2 = make_shared<rotate_y>(box2, -18);
  box2 = make_shared<translate>(box2, vec3(130,0,65));
  world.add(box2);

  point3 lookfrom(278, 278, -800);
  point3 lookat(278, 278, 0);
  vec3 vup(0, 1, 0);
  auto dist_to_focus = 10.0;
  auto aperture = 0.0;
  auto vfov = 40.0;
  auto t0 = 0.0;
  auto t1 = 1.0;

  cam = camera(
    lookfrom, lookat, vup, vfov, aspect, aperture, dist_to_focus, t0, t1
  );

  return world;
}
リスト 3.36 [main.cc] アルミニウムマテリアルを追加したコーネルボックス

レンダリングされる画像では、直方体からの反射で照らされる部分 (画像左上) にノイズが多く残る。これは直方体を向くレイを多くサンプルしていないためである。

図 3.16:
鏡面反射と PDF の両方に対応したコーネルボックス

アルミニウムの直方体を PDF に含めることもできる。ただここでは、実装が簡単なガラス球を考えよう。

球オブジェクトのサンプリング

図 3.17:
球に接する円錐

球の外側の点から球を立体角に関して一様にサンプルするとき、球に接する円錐をサンプルすると考えても同じことになる (上図)。図中の \(\theta_{\text{max}}\) が求まっているとする。ランダムな方向の生成で説明したように、\([0,1]\) の乱数 \(r_{2}\) に対して \[ r_2 = \int_{0}^{\theta} 2\pi \cdot f(t) \cdot \sin t\, dt \] が成り立つとき \(\theta\) の PDF は \(f(t)\) となる1。この問題の \(f(t)\) は未知の定数 \(C\) だから \[ r_2 = \int_{0}^{\theta} 2\pi \cdot C \cdot \sin t\, dt \] となる。よって \[ r_2 = 2\pi \cdot C \cdot (1-\cos \theta ) \] すなわち次を得る: \[ \cos \theta = 1 - \frac{r_2}{2 \pi \cdot C} \]

\(r_{2} = 1\) で \(\theta = \theta_{\text{max}}\) となることを使えば \(C\) が求まり、 \[ \cos \theta = 1 + r_2 \cdot (\cos \theta_{\text{max}} - 1) \] を得る。\(\phi\) に関しては以前と同様であり、次の関係が分かる: \[ \begin{aligned} z & = \cos\theta = 1 + r_2 \cdot (\cos \theta_{\text{max}} - 1) \\ x & = \cos\phi \cdot \sin\theta = \cos(2\pi \cdot r_1) \cdot \sqrt{1-z^2} \\ y & = \sin\phi \cdot \sin\theta = \sin(2\pi \cdot r_1) \cdot \sqrt{1-z^2} \end{aligned} \]

続いて \(\theta_{\text{max}}\) を考える。図から \(\sin \theta_{\text{max}} = R / \| \mathbf{c} - \mathbf{p} \|\) が分かるから、 \[ \cos \theta_{\text{max}} = \sqrt{1 - \frac{R^2}{\|\mathbf{c} - \mathbf{p}\|^{2}}} \] である。

さらに方向から PDF を計算できなければならない。方向が球を指すなら、PDF は \(1/({\footnotesize \text{立体角}})\) と表せる。球全体の立体角はいくつだろうか？ この立体角には上で求めた \(C\) が関係する。定義より、立体角は単位球上に射影した領域の面積だから、その値は積分を使って \[ {\footnotesize \text{球の立体角}} = \int_{0}^{2\pi} \int_{0}^{\theta_{\text{max}}} \sin \theta \, d\theta d\phi = 2 \pi \cdot \frac{1}{2 \pi C} = 2 \pi \cdot (1-\cos \theta_{\text{max}}) \] と求まる2。

こういった計算の正しさを確認しておこう。私はこういうとき極端な値を試すことにしている (高校で物理を教えていたホートン先生に感謝する)。球の半径が \(0\) のとき \(\cos\theta_{\text{max}} = 1\) および \({\footnotesize \text{立体角}} = 0\) であり、矛盾はない。また球が \(p\) に接するとき \(\cos \theta_{\text{max}} = 0\) であり、立体角は \(2\pi\) で半球の表面積に等しくなる。ここにも矛盾はない。

球のコードの更新

sphere に PDF に関する関数を二つ追加する:

double sphere::pdf_value(const point3& o, const vec3& v) const {
  hit_record rec;
  if (!this->hit(ray(o, v), 0.001, infinity, rec))
    return 0;

  auto cos_theta_max = sqrt(1 - radius*radius/(center-o).length_squared());
  auto solid_angle = 2*pi*(1-cos_theta_max);

  return  1 / solid_angle;
}

vec3 sphere::random(const point3& o) const {
  vec3 direction = center - o;
  auto distance_squared = direction.length_squared();
  onb uvw;
  uvw.build_from_w(direction);
  return uvw.local(random_to_sphere(radius, distance_squared));
}
リスト 3.37 [sphere.h] PDF に対応した球

次のユーティリティ関数を使っている:

inline vec3 random_to_sphere(double radius, double distance_squared) {
  auto r1 = random_double();
  auto r2 = random_double();
  auto z = 1 + r2*(sqrt(1-radius*radius/distance_squared) - 1);

  auto phi = 2*pi*r1;
  auto x = cos(phi)*sqrt(1-z*z);
  auto y = sin(phi)*sqrt(1-z*z);

  return vec3(x, y, z);
}
リスト 3.38 [pdf.h] ユーティリティ関数 random_to_sphere

最初の試みとしてライトではなく球だけをサンプリングするレンダリングを試す:

shared_ptr<hittable> light_shape =
  make_shared<xz_rect>(213, 343, 227, 332, 554, nullptr);
shared_ptr<hittable> glass_sphere =
  make_shared<sphere>(point3(190, 90, 190), 90, nullptr);

for (int j = image_height-1; j >= 0; --j) {
  std::cerr << "\rScanlines remaining: " << j << ' ' << std::flush;
  for (int i = 0; i < image_width; ++i) {
    color pixel_color(0, 0, 0);
    for (int s=0; s < ns; ++s) {
      auto u = (i + random_double()) / image_width;
      auto v = (j + random_double()) / image_height;
      ray r = cam->get_ray(u, v);
      pixel_color += ray_color(r, background, world, glass_sphere, max_depth);
    }
リスト 3.39 [main.cc] 球だけをサンプリングする

ノイズの多いコーネルボックスが得られるが、球の下の集光模様は上手くレンダリングされている。私の環境ではライトをサンプルするときの五倍の時間がかかった。ガラスに当たるレイは計算コストが高いということだろう。

図 3.18:
新しい PDF 関数を使ったガラス球があるコーネルボックス

オブジェクトのリストに対する PDF

球とライトの両方をサンプルするべきだろう。つまり二つの密度関数の混合密度を使ってサンプルを行いたい。ray_color() 関数に hittable のリストを渡して混合 PDF を構築させることもできるし、hittable_list に PDF 関連の関数を追加することもできる。どちらの方法でも上手く行くと思うが、ここでは hittable_list に関数を追加する方針を取る。

double hittable_list::pdf_value(const point3& o, const vec3& v) const {
  auto weight = 1.0/objects.size();
  auto sum = 0.0;

  for (const auto& object : objects)
    sum += weight * object->pdf_value(o, v);

  return sum;
}

vec3 hittable_list::random(const vec3& o) const {
  auto int_size = static_cast<int>(objects.size());
  return objects[random_int(0, int_size-1)]->random(o);
}
リスト 3.40 [hittable_list.h] 混合密度の作成

ray_color() に渡すリストを main() で作る:

hittable_list lights;
lights.add(make_shared<xz_rect>(213, 343, 227, 332, 554, nullptr));
lights.add(make_shared<sphere>(point3(190, 90, 190), 90, nullptr));
リスト 3.41 [main.cc] シーンの更新

1000 サンプル/ピクセルとすれば以前と同様の整った画像が得られる:

図 3.19:
ガラスとライトの混合 PDF を使ったコーネルボックス

アクネの除去

鋭い読者が 図 3.19 に黒い斑点があることを指摘した。全てのモンテカルロレイトレーサーは次の式をメインループに持つ:

pixel_color = average(たくさんのサンプル)

白または黒のアクネを持つ画像が得られたなら、一つの "悪い" サンプルがピクセル全体をダメにしてしているということだ。そのサンプルは非常に大きい値もしくは NaN (not a number) を持つ。上の画像に含まれるアクネはおそらく NaN であり、私の環境ではおよそ一千万から一億レイに一度起きている。

さて重要な決断をしよう: 隅々まで調査してこのバグを取り除くか、とりあえず出力から NaN を消して今後何も起こらないことを祈るかだ。私は必ず怠惰な方を選択する。浮動小数点の難しさを知っていればなおさらだ。NaN はどのように確認できるだろうか？ NaN について私が常に念頭に置いているのが、NaN はそれ自身と等しくならないという事実である。これを使って write_color 関数を変更し、NaN の要素を \(0\) で置き換える:

void write_color(std::ostream &out, color pixel_color, int samples_per_pixel) {
  auto r = pixel_color.x();
  auto g = pixel_color.y();
  auto b = pixel_color.z();

  // NaN の要素を 0 に置き換える (詳細は『週末レイトレーシング: 余生』を参照)
  if (r != r) r = 0.0;
  if (g != g) g = 0.0;
  if (b != b) b = 0.0;

  // 色の合計をサンプルの数で割り、gamma = 2.0 のガンマ補正を行う
  auto scale = 1.0 / samples_per_pixel;
  r = sqrt(scale * r);
  g = sqrt(scale * g);
  b = sqrt(scale * b);

  // 各成分を [0,255] に変換して出力する
  out << static_cast<int>(256 * clamp(r, 0.0, 0.999)) << ' '
      << static_cast<int>(256 * clamp(g, 0.0, 0.999)) << ' '
      << static_cast<int>(256 * clamp(b, 0.0, 0.999)) << '\n';
}
リスト 3.42 [color.h] NaN に対応した write_color 関数

これでめでたく黒い斑点が消える:

図 3.20:
アクネに対処した write_color 関数を使うコーネルボックス

  1. 訳注: 積分において、原点と球の中心を結ぶ直線を \(z\) 軸としていることに注意。[return]

  2. 訳注: ここでも軸の方向に注意。[return]

広告
このページの文章はテキストゆれないくんを使って校正されました。