Category Archives: Uncategorized

2D大域照明計算のススメ (2D Global Illumination no Susume)

レイトレ合宿3のアドベントカレンダー5週目の記事として, 2Dの大域照明計算についての記事を書きました (PDFです). 文章は英語ですが, 気が向いたら訳すかもしれません.

This article is written for the advent calender of the ray-tracing camp 3. In this article I introduced the 2D global illumination which is useful for the verification of the rendering techniques. In the first part of the article, I briefly introduced the theoretical concepts of 2D global illumination such as 2D radiometry and 2D rendering equation and pointed out its similarities between 3D counterparts. Next I described the implementation details of the 2D global illumination techniques focusing on some tips and pitfalls.

[PDF Link]

Multiple importance samplingの罠

はじめに

レイトレ合宿のアドベントカレンダー8週目の記事として, Multiple importance sampling (MIS) で躓きやすい個所についていくつか実例を挙げて紹介したいと思います.

Multiple importance sampling (MIS) とは, Veachら [1], [2] によって開発されたモンテカルロ積分のための手法で, ひとつの分布から重点的サンプリングすると効率の悪い分布を, 複数の分布からのサンプルを組み合わせて用いることで効率よくサンプリングできます.

ですが, MISがかかわると非常に実装がバグりやすいことで有名です.
Bidirectional path tracing (BPT) [4], [2] がそのひとつで, 正しく実装するためにはMISに対する正しい理解と慎重な実装が求められます.
本記事では実例をあげながらなぜ躓きやすいのか, どういう部分に気を付けて実装すればいいのかについて解説していこうと思います.

MISの概要

モンテカルロ積分によって推定したい積分を

    \[ I=\int_{\Omega} f(x)d\mu(x) \]

とします. このとき, n個の異なる確率分布を考え, その密度関数をp_1,\cdots,p_nとします.
p_iに対する仮定として, {\rm supp}(f)\subseteq \cup_i^n {\rm supp}(p_i) (すべてのf(x)>0となるxで, 少なくともひとつの分布でp_i(x)>0となる) を考えます.
このとき,

    \[ F = \sum_{i=1}^n \frac{1}{n_i} \sum_{j=1}^{n_i} w_i(X_{i,j})\frac{f(X_{i,j})}{p_i(X_{i,j})} \]

multi-sample estimatorと呼びます.
ここで, w_i(x)x\in \cup_i^n {\rm supp}(p_i)に対して定義される重み関数です.

MISを正しく用いるには, 上のFがどのような場合にunbiasedになるのかを把握する必要があります.
重要な特徴として, 重み関数w_iが以下の条件を満たす場合にFはunbiasedとなります.

  • 条件1 : すべてのx\in \cup_i^n {\rm supp}(p_i)に対して, \sum_i^n w_i(x)=1.
  • 条件2 : p_i(x)=0ならばw_i(x)=0.

簡単のため, n=2のときで証明します. 一般の場合もほぼ同じだと思います.
\Omega_1 = {\rm supp}(p_1), \Omega_2 = {\rm supp}(p_2)とします.

    \begin{align*} \mathbb{E}[F] &= \int_{\Omega_1} w_1(x)\frac{f(x)}{p_1(x)}p_1(x)d\mu(x)+ \int_{\Omega_2} w_2(x)\frac{f(x)}{p_2(x)}p_2(x)d\mu(x) \\ &= \int_{\Omega_1} w_1(x)f(x)d\mu(x) + \int_{\Omega_2} w_2(x)f(x)d\mu(x) \\ &= \int_{\Omega_1\setminus\Omega_2} w_1(x)f(x)d\mu(x) \\ &+ \int_{\Omega_1\cap\Omega_2}\left( w_1(x) + w_2(x) \right) f(x)d\mu(x) \\ &+ \int_{\Omega_2\setminus\Omega_1} w_2(x)f(x)d\mu(x) \\ \end{align*}

ここで, x\in\Omega_1\cap\Omega_2に対して条件1より

    \[ \int_{\Omega_1\cap\Omega_2}\left( w_1(x) + w_2(x) \right) f(x)d\mu(x) = \int_{\Omega_1\cap\Omega_2} f(x)d\mu(x). \]

また, x\in\Omega_1\setminus\Omega_2に対して条件2よりw_2(x)=0.
条件1よりw_1(x)=1-w_2(x)=1.
ゆえに,

    \[ \int_{\Omega_1\setminus\Omega_2} w_1(x)f(x)d\mu(x) = \int_{\Omega_1\setminus\Omega_2}f(x)d\mu(x). \]

x\in\Omega_2\setminus\Omega_1のときも同様.
以上により,

    \begin{align*} \mathbb{E}[F] &= \int_{\Omega_1\setminus\Omega_2} f(x)d\mu(x) + \int_{\Omega_1\cap\Omega_2} f(x)d\mu(x) + \int_{\Omega_2\setminus\Omega_1} f(x)d\mu(x) \\ &= \int_{\Omega} f(x)d\mu(x) = I \end{align*}

で, Fはunbinasedです.

重み関数としては, balance heuristics

    \[ w_i(x)=\frac{n_i p_i(x)}{\sum_k n_k p_k(x)} \]

や, power heuristics

    \[ w_i(x)=\frac{(n_i  p_i(x))^\beta}{\sum_k (n_k p_k(x))^\beta} \]

がよく用いられます.

デルタ関数の扱い

本記事のメインパートです.

MISの扱いに関して注意する話題のひとつとして, 組み合わせる確率分布がデルタ関数のケースがあります. 光輸送シミュレーションにおいて, デルタ関数はさまざまな場所にあらわれますが, MISと同時に用いる場合いくつか注意が必要です.

サンプルできない分布は考慮しない

前節ではmulti-sample estimatorがunbiasedになる条件に付いて議論しました.
ここではその条件を拡張して, デルタ関数にも対応させることを考えます.

条件2′ : すべてのx\in \cup_i^n {\rm supp}(p_i)に対して, xp_iから生成される確率が0のとき, w_i(x)=0.

分布関数がゼロなのは「タテがゼロ」のケースで,
デルタ関数の場合は「ヨコがゼロ」のケースです.
いずれもサンプリングできません.

例として, 3つの確率密度関数p_1, p_2, p_2を考えて, そのうちのp_1x_0で定義されるデルタ関数だとしてp_1=\delta(x-x_0)とします. (実はこの書き方は数学的にあやしい : デルタ関数を測度として扱った際にルベーグ測度に対してRadom-Nikodym微分が存在しないので, 密度関数が存在しない. が, 便宜的に広く使われているのでここでは気にしないことにします^^).

ここでp_2(x_0)>0, p_3(x_0)>0ですが, \Omegaは連続なのでp_2, p_3ではX=x_0はサンプリングできません. このような場合, 条件2′を用いてw_2(x_0)=0, w_3(x_0)=0としてp_1, p_2を重みから外す必要があります. (結果として, w_1(x_0) = 1となります.)

デルタ関数の評価値について

実際のデルタ関数は1点で無限大となりますが, サンプリングの際に確率分布の評価値も同時に計算して, デルタ関数の場合には便宜的に有限の値を評価値としておくテクニックがよく用いられています. (私の実装ではcos thetaで割ったものを入れています. PBRTやmitsubaも同じようなことをやっています.)

このデルタ関数の評価値を有限にするテクニックですが, あくまでも便宜的なものなのでそのままの値を用いるとうまくいかないケースがあります. たとえばbalance heuristicsやpower heuristicsを用いている場合に, デルタ関数の評価値として有限の値を用いていると, 他の(非デルタ関数の)評価値の影響により正しい重みの値が計算できません. この場合, [2] p.314 にあるようにフラグを駆使して条件1,2′に従うような重みを計算します.

BPT用のシンプルな重み関数

最後に実際の例として, BPT用のシンプルな重み関数w_{Simple}の実装を示します.
実用上はあまり役に立たない重みですが, 「サンプルできない分布は考慮しない」ことがよくわかると思います.

BPTで用いる光路生成のための確率密度関数をp_{s,t}として, 頂点数nの光路\bar{x}_nに対して,

    \[ w_{Simple}(\bar{x}_n)=\frac{1}{\#\{ (s,t) : s+t=n \wedge P\{  p_{s,t} \text{ generates } \bar{x}_n \} > 0 \}}. \]

与えられたfullPathを生成できる確率分布の数を数えています.
デルタ関数を含む可能性のある場所として

  • s=0, t=nで, ポイントライトが用いられているとき (16行目)
  • s=n, t=0で, ピンホールカメラが用いられているとき (24行目)
  • その他のs,tで, specular BSDFとかdirect lightとかが用いられているとき (35行目)

をチェックしています.

 References

  1. E. Veach and L. Guibas, Optimally combining sampling techniques for Monte Carlo rendering, Procs. of SIGGRAPH 95, pp.419-428, 1995.
  2. E. Veach, Robust Monte Carlo methods for Light Transport Simulation, PhD Thesis, Stanford University, 1997.
  3. J. T. Kajiya, The Rendering Equation, Procs. of SIGGRAPH 86, pp.143-150, 1986.
  4. E. Veach and L. Guibas, Bidirectional estimators for light transport, Procs. of the Fifth Eurographics Workshop on Rendering, pp.147-162, 1994.

 

 

 

IvyBridge RdRand

Introduction

I heard that HW random number generator is implemented in Intel’s CPU since IvyBridge architecture. Today I tried to use it.

How to use?

We can directly access rdrand instruction through assembler but it is easy to use intrinsics offered by the header immintrin.h. The header is same one as AVX intrinsics. For obtaining uniform double random number in the range of [0, 1], we can use the following simple code:

Experiments

Speed test vs. Mersenne twister

This article says rdrand is capable of 7000M instructions / seconds. How fast is it? I compared the speed to MT(std::random implementation) in a simple experiment. Test code:

Execution time (executed in i7-3517U):

  • MT : 0.421 seconds
  • RdRand : 1.322 secods

Trace paths!!

I added simple modification to smallpt (simply replacing erand48 to the above rdrand function). The source is here. And comparison with MT version (note that erand48 is not used in Windows so I used std::random one instead). Source: MT version and rdrand version.

Execution time:

  • MT : 89.107 seconds
  • RdRand : 100.53 seconds

comparison

Seems no difference.

 

VS2012 + googletest

Introduction

Microsoft Visual Studio 2012 introduces a native C++ unit testing tool. It is definitely good thing because we don’t have to write write c++/cli code merely for testings. Still, there is a case when we cannot use VS tools for compatibility reason; in order to make a software platform independent, we’d better to make its testings codes platform independent too. There are many platform independent unit testing frameworks for C++. Among them, I use googletest for this article.

Using googletest with VS2012

Download and build

  • Download latest googletest (I downloaded gtest-1.6.0.zip)
  • Unzip it and open \msvc\gtest.sln. It is VS2003 solution so it requires upgrade.
  • Change number of variadic template parameters (in VS2012, 5 is the default value!). Add _VARIADIC_MAX=10 to PROJECT > Properties > Configuration Properties > C/C++ > Preprocessor > Preprocessor Definitions.
  • Press F7 and build solution both for Debug and Release mode.

Create your project

  • Add references to the built libraries to your project. For example, if the include files and built libraries are located on (your solution directory)\external\gtest-1.6.0\include and \lib. Add $(SolutionDir)external\gtest-1.6.0\include to VC++Directory > Include Directories and $(SolutionDir)external\gtest-1.6.0\lib to VC++Directory > Include Directories respectively, and add *.lib to Linker > Input.
  • Add your test code, e.g.

 Integrate googletest to VS2012 Test Explorer

  • If you don’t see Test Explorer, click TEST > Windows > Test Explorer.
  • Install googletest adapter plugin, and tests in your solution is automatically detected and listed in the Test Explorer window.
  • Click Run All or run selected tests.

VS2012 code coverage tool with googletest

Without any settings, the code coverage tool detects internal code of googletest. The easiest way to avoid it is to wrap your code into a namespace (e.g. mynamespace::*), then tell the code coverage tool not to include the code outside the namespace. 

  • Create .runsettings file in your solution directory (use the template in this page).
  • Add the following snippet to <Functions></Functions>.

  • Click TEST > Test Settings > Select Test Settings File, and load your .runsettings file.