ERPT 実装: エネルギー再分布ステップの枠組みと摂動ルーチン

by syoyo

さて、そろそろ ERPT のコアとなる部分を実装していくことにしよう。

エネルギー再分布

パストレにより、まずはデポジットする e_d が求まったはずであるので、
あとはこれをおのおのの変異した経路の対応するピクセルに加算していけばよいはずだ。

e_d はパストレの結果画像バッファから引っこ抜いてきた値を e_ave として、
変異させる回数で割れば求まる。

e_d の計算のコードはこんな感じになるだろうか。

for (j = 0; j height; j++) {
        for (i = 0; i width; i++) {

                vector_t e_ave;
                vector_t e_d;

                // get average energy of this pixel from
                // pass1's path tracing result.
                e_ave.x = render->e_ave[3 * (i + j * render->width) + 0];
                e_ave.y = render->e_ave[3 * (i + j * render->width) + 1];
                e_ave.z = render->e_ave[3 * (i + j * render->width) + 2];

                // e_d is the deposition energy for each mutated
                // path Y.
                e_d.x = e_ave.x / (double)render->nmutations;
                e_d.y = e_ave.y / (double)render->nmutations;
                e_d.z = e_ave.z / (double)render->nmutations;

                energy_redistribution_sampling(i, j, e_d);
        }
}

render->e_ave はパストレの結果の画像バッファをコピーしたものである。
render->nmutation が変異の回数(式 7 の k)である。

energy_redistribution_sampling() は、図 7 相当の処理を行うルーチンであり、
以下のようになる。

void
energy_redistribution_sampling(int x, int y, vector_t e_d)
{
        bool     hit;
        int      k;
        double   u, v;
        vector_t dir;

        path_t   path;
        double   X_f;
        char     str[1024];

        if (fabs(e_d.x) nmutations; k++) {

                u = (double)x + randomMT();
                v = (double)(render->height - 1.0)
                  - (double)y + randomMT();

                get_eyedir(&dir, &render->eyepos, &render->eyedirpos,
                           &render->du, &render->dv, render->T, u, v);

                // create a path, x, int the current pixel.
                hit = path_trace(&path, &render->eyepos, &dir);

                // RGB -> Luminance
                X_f = CIEcol(path.radiance);

                // X_f(x) > 0.0
                if (hit && X_f > 0.0) {
                        equal_deposition_flow(&path,
                                              path.radiance,
                                              render->nchains,
                                              e_d);
                }
        }
}

まず、e_d の輝度値がゼロのとき、つまりピクセルが完全に遮蔽されているような領域だったり、
背景にヒットするような位置である場合は、処理を打ち切る。
(たまたまサンプル数が不十分で、光源へのシャドウレイが全部遮蔽する場合でも
e_d がゼロになる可能性があるが、、、まあそのときはしょうがないということにしよう)

e_d に輝度があれば、まず適当にパストレを行って初期の経路を生成する。
あとはこの経路が同じくどこかに遮蔽されたり背景にヒットするなどされずに
非ゼロな輝度を持つのであれば、 equal_deposition_flow() へと進む。

equal_deposition_flow() はこんな感じ。

void
equal_deposition_flow(path_t *x, vector_t e, int m, vector_t e_d)
{
        int    i, j;
        int    num_chains;
        double e_luminance, e_d_luminance;

        e_luminance    = CIEcol(e);
        e_d_luminance  = CIEcol(e_d);

        num_chains = (int)floor(randomMT() + e_luminance / (m * e_d_luminance));

        render->total_chains += num_chains * render->nchains;
        render->chain_events++;

        for (i = 0; i nchains; j++) {
                        // mutate_path(y, x);
                        // relative_path_density_(y, x);
                        // ごにょごにょ...
                }
        }
}

というわけで、ERPT(MLT でもそうだが) の核心のひとつである、
経路の変異処理、mutate_path()、を実装する部分までやってきた。

変異処理

ERPT では、変異はコースティクス摂動とレンズ摂動の二つしかない
(少なくとも論文ではこの二つしか使っていない)。
つまりはピクセルの位置や光源の方向をちょっとだけずらしてあげる。

また、各変異において、経路の長さは一定で変わらない。

MLT は bidirectional mutation により経路の長さが変わったりするので、
ERPT のほうが変異処理はかんたんである。

まずは、コースティクス摂動とレンズ摂動のために、レイの方向を摂動
(つまりちょっとだけずらす)させる必要があるので、このための基礎コード
を実装することにします。

ERPT でも、摂動処理の基本は Eric Veach の MLT の論文をもとにしている。

これを行うためのコードが、すでに ERPT の第一著者が大学のテクニカルレポートと
して出版している “A practical introduction to Metropolis Light Transport”
に付録として載っている。こちらはパストレを用いて MLT の入門的な実装を
解説している。ERPT と似ている部分も多く、
摂動部分についての解説は ERPT にそのまま流用できるはずだ。

David Cline and Parris K. Egbert,
“A Practical Introduction to Metropolis Light Transport”
BYU Technical Report, 2005
http://rivit.cs.byu.edu/a3dg/publications.php

というわけでこれをそのまま利用すればよい。簡単だ。

void
pixel_offset(double *x, double *y, double r1, double r2)
{
        double phi = randiomMT() * 2.0 * M_PI;
        double r   = r2 * exp(-log(r2 / r1) * randomMT());

        (*x) += r * cos(phi);
        (*y) += r * sin(phi);
}

void
angular_offset(vector_t *N, double theta1, double theta2)
{
        vector_t U, V;
        vector_t axisZ, axisY;

        axisZ.x = 0.0; axisZ.y = 0.0; axisZ.z = 1.0;
        axisY.x = 0.0; axisY.y = 1.0; axisY.z = 0.0;

        // Make a UVN coordinate system from N
        if (fabs(N->x) x = N->x + r * cos(phi) * U.x + r * sin(phi) * V.x;
        N->y = N->y + r * cos(phi) * U.y + r * sin(phi) * V.y;
        N->z = N->z + r * cos(phi) * U.z + r * sin(phi) * V.z;

        vnormalize(N);
}

pixel_offset() がレンズ摂動(画像面のサンプル位置 x, y をずらす)、
angular_offset() がコースティクス摂動(角度をちょっとずらす)となる。

randomMT() は [0.0, 1.0) の範囲で実数ランダム値を返す。

pixel_offset() の r1, r2 は摂動量を調整するパラメータである。
MLT の元論文では、r1 は 0.1、r2 は画像面の面積(横ピクセル数 x 縦ピクセル数)となっている。
A practical 〜 では r1 は 0.1、r2 は横ピクセル数の 10 % に設定している

ちなみに、MLT のレンズ部分経路変異(lens subpath mutation) に相当する処理を行う場合は、
A practical 〜 では pixel_offset() により大きな値を設定、r1 は 1.0、r2 は横ピクセル数の 25 %、
することで代用しているとのこと。

angular_offset() の theta1, theta2 は、MLT も A practical 〜 でも、
theta1 は 0.0001、theta2 は 0.1 に設定しているとのこと。

では、実際に経路を摂動させてみよう…

to be continued…

Advertisements