\newcommand{\expected}{\mathbb{E}}
\newcommand{\variance}{\mathbb{V}}

この記事はレイトレ合宿8のアドベントカレンダー2週目の記事として書かれました。

はじめに

皆さんこんにちは、@Shocker_0x15です。今年2022年、3年ぶりに開催されることになったレイトレ合宿ですが、なんと開催場所が沖縄！ですね。割りと現地参加するハードルが高くなったのにも関わらず新規参加者も多く今から既に楽しみにしております。2018年頃のGPUよりレイトレーシング支援ハードウェアが搭載され始めており、SIGGRAPHなどの学会でもリアルタイムレンダリングでの使用を視野に入れた手法の提案が数多くなされています。レイトレの波を感じますね。ハードウェアがどういう形に収束していくのかはわかりませんが、ソフトウェア面におけるモンテカルロレイトレーシング(的な手法)の重要度はずっと増していくことでしょう。

さて、今回のレイトレ合宿では本戦に加えてデノイズの質を競うらしいサブ部門があります。また本戦でも今年のレギュレーションでは、動画作品の提出が義務付けられているため、こちらでもデノイザーを使って時間方向の安定性を確保することが技術面では結構重要になるのではないでしょうか。このような背景もあり、今回のアドベントカレンダーではデノイザーについて取り上げてみます。モンテカルロレイトレーシング用のデノイザーには大別してMachine Learningベースのものと、そうでないもの(典型的には画像空間におけるフィルタリング)の2種類があります。リアルタイムレンダリングでは現状後者に属する画像空間フィルタリングベースのものが多く使われている気がします。その中でも様々な手法の基となっていると思われる"Spatiotemporal Variance-Guided Filtering"、通称SVGFのアイディアについて紹介します。元論文は以下のリンク先にあります。

当記事は実装の詳細を細かに説明するよりはSVGFのアイディアの概要、使われているテクニックの紹介としたいと思います。
代わりにSVGFを実装したものを以下のリポジトリに置いております。

フィルタリングというアイディア

モンテカルロレイトレーシングは、エリアライトによる照明や光沢面の反射、間接照明など様々な表現を統一されたアルゴリズムで実現できる強みがあるのですが、一方でレンダリング画像は確率的な分散を持ち視覚的なノイズとして知覚されるため、これを気にならなくなるレベルまで抑えるには大量のサンプリングを行う必要があります。しかしサンプル数¹の増加はそのままレンダリング時間の増加に繋がります。1/60[秒]といった時間内で処理を終える必要のあるリアルタイムレンダリングではもっと問題です。2022年現在だとシーンの規模やもちろんGPUグレードにもよるのですが、1ピクセルあたり1つのパストレーシングサンプル(複数レイからなる)、つまり1[spp]でもフレーム中のかなりの時間を占めます。1[spp]で得られる画像は当ページトップの左側のようになり、そのままでは視覚的なノイズが多すぎて使い物になりません。

典型的なシーンを考えると、物体の境界やマテリアルの境界、影の境界以外における空間的な近傍同士では(時間をかけた場合に得られる)レンダリング結果の輝度・色は似ています。そこで特性が似ている物体表面同士で分散の多いノイジーなライティング結果をフィルタリングして実質のサンプル数を大きくすることでレンダリング品質を向上させられそう、というアイディアが登場することになります。

画像空間フィルタリング

直接3Dな(ワールド)空間でのフィルタリングはメリットもありますが、面倒なので一般的な画像処理のように画像空間中でのフィルタリングを考えます。通常のガウシアンフィルターなどのように物体表面の特性を考えずに何でもかんでもフィルタリングしてしまうと、レンダリング画像が明らかにボケボケになるのでバイラテラルフィルターのような選択的なウェイト、edge-stopping関数を考慮します。具体例として次のような要素に基づいたウェイトを考えます。

ピクセル間の距離
これは通常のボックスフィルターやガウシアンフィルターのウェイトです。
ワールド空間的な距離
近傍ピクセルはワールド空間でも近いことが多いが、必ずしも正しくはないのでデプスバッファーの値などを使います。
法線ベクトルの類似度
法線ベクトルの異なる面間でフィルタリングしてしまうと物体のエッジが見えなくなるので考慮に加えます。
輝度の類似度
輝度差があまりに大きいケースは影の境界だったりするのでウェイトを小さくします。

画像空間以外のフィルタリング

これにはPath Space Filteringといったアイディアなどが含まれますが、当記事では紹介しません。

SVGF

パイプライン概要

純粋なSVGFのパイプラインは次のようになります。

G-Bufferの生成。
アルベドや法線、オブジェクトIDやマテリアルID、モーションベクターを出力します。
パストレーシングなどによるレンダリング
G-Bufferからプライマリーヒットの物体表面を復元、セカンダリーレイ以降をトレースしてレンダリング、ノイジーなレンダリング結果をバッファーに記録します。
また各ピクセルのモーメントやサンプル数も記録。
各ピクセルの分散を推定
モーメントやサンプル数を用いてピクセルの分散をバッファーに記録。
À-trousフィルターを何度か繰り返す
各ピクセルの色と分散を、デプスや法線、輝度の類似度、分散に基づいてウェイトを制御されたÀ-trousフィルターにかける。

アルベド除去

論文ではAlbedo demodulation(とmodulation)と書かれています。SVGFにピクセルごとのライティング結果を直接入力するのではなく、プライマリーヒットのアルベドで割ったものを入力、デノイズ後に再度掛け直します。これによってプライマリーヒットにおけるテクスチャーのディテールを保持することができます。ランバート反射の場合はアルベドとは単に反射率パラメターのことになりますが、複雑なBRDFの場合は何をアルベドとして使うか曖昧です。おそらく視線ベクトルに対するDirectional-Hemispherical Reflectance(の推定値)であったり、スペキュラー成分しか無い場合には垂直反射率などが好ましくなると思います。

以下はディフューズとスペキュラー両方の成分を持つBRDFに関してアルベド除去を行った例です。複合的なBRDFに関して完璧なアルベド除去は難しいので微妙に色や模様が残っている箇所もありますが、床や壁の模様を除去できていることがわかります。レンガ状のブロック境界が見えているのはバンプマップを行っていることが理由です。

ピクセルごとの分散の推定

フィルタリング半径を大きくすることで実質のサンプル数をより稼ぐことができるのですが、一方でedge-stopping関数があるといってもどうしてもフィルタリング結果はある程度ぼけたものになってしまいます。そもそも分散が小さいケースではそこまでフィルタリング半径を大きくする必要がありません。極端な話、分散がゼロならフィルタリングするメリットはありません。そこで、ピクセルごとの分散に従って実質のフィルタリング半径をコントロールしたくなります。

分散は以下の有名な式にしたがって求めることができます。
$$ \variance[X] = \expected[X^2] - \expected[X]^2 $$

式自体は非常に簡単なのですが $ \expected[X^2] $ も $ \expected[X] $ も1フレーム1サンプルの場合求めることができません。(正確さはともかく)最低でも2サンプル必要です。しかし、カメラが静止している(かつピクセル中のサンプル位置も同じ)場合を考えれば、前フレームの同じピクセルが持っていたサンプルは現在のフレームと同じ条件下で評価したサンプルと言えるので、分散の評価に使えそうです。そこで各ピクセルにおいてノイジーなサンプルを得た後、以下のコード片のようにピクセルにおける輝度のモーメント2つとサンプル数をトラッキングします。

RGB color = render(pix); // レンダリングで求めた確率的なピクセルカラー
RGB albedo = albedoBuffer[pix];

// 物体表面色に依存しないライティングに関してデノイズするために
// アルベドでライティング結果をdemodulation。
RGB demColor = color / albedo;

float luminance = calcLuminance(demColor);
float sqLuminance = luminance * luminance;

// 前フレームのモーメントとサンプル数を読み込む。
MomentPair_SampleCount prevMS = prev_momentPair_sampleCount_buffer[pix];

// E[X], E[X^2] を求める。
uint32_t newSampleCount = prevMS.sampleCount + 1;
float curWeight = 1.0f / newSampleCount;
float avgLuminance = (1 - curWeight) * prevMS.firstMoment + curWeight * luminance;
float avgSqLuminance = (1 - curWeight) * prevMS.secondMoment + curWeight * sqLuminance;

// モーメントとサンプル数を書き込む。
MomentPair_SampleCount curMS;
curMS.firstMoment = avgLuminance;
curMS.secondMoment = avgSqLuminance;
curMS.sampleCount = newSampleCount;
cur_momentPair_sampleCount_buffer[pix] = curMS;

分散を求めるコードは書くまでも無いかもしれませんが以下のようになります。

MomentPair_SampleCount curMS = prev_momentPair_sampleCount_buffer[pix];
float variance = curMS.secondMoment - curMS.firstMoment * curMS.firstMoment;

RGBそれぞれについて分散を求めてフィルタリングするほうが正確ですが、SVGFでは代わりに輝度の分散を使用します。

カメラや物体が動いている場合

カメラや物体が動いている場合、前フレームのデータから現在処理中のピクセル座標でデータを読むと、異なる物体表面のサンプルになってしまいます。そこでG-Buffer生成時に計算しておいたモーションベクターを用いて前フレームの対応ピクセルを計算、同ピクセルからサンプルのモーメントとサンプル数のデータを読みます。同ピクセルに対応する物体表面が法線やデプスを用いた類似度テストに通らない場合は読んだデータを棄却して現在フレームのデータだけで結果を計算します(前フレームのサンプル数ゼロと等価)。論文では前フレームの対応ピクセルを1つに決定するのではなく、近傍4ピクセルのデータを読み出しバイリニア補間した結果を前フレームのデータとしています。

//MomentPair_SampleCount prevMS = prev_momentPair_sampleCount_buffer[pix];
MomentPair_SampleCount prevMS = getPrevious_MomentPair_SampleCount(pix - motionVector);

このように毎フレーム複数サンプル評価するのではなく、過去フレームと合わせて時間的に分散を求めている処理がSVGFのSpatiotemporalの所以です。

一方で、類似度テストにより前フレームのデータを棄却した場合や、そもそも1フレーム目には1サンプル分しかデータが無いため、分散を求めることができません。正確に言えば求めることはできるもののゼロになってしまいます。分散が小さいほどフィルタリング半径を絞るので、これは適切なフィルタリングができない可能性があります。
そこでSVGFは有効なサンプル数が少ないケースにおいて、Spatiotemporalたる所以である分散の空間的な推定へのフォールバックを行います。

MomentPair_SampleCount curMS = cur_momentPair_sampleCount_buffer[pix];
float firstMoment = curMS.firstMoment;
float secondMoment = curMS.secondMoment;

// 空間的な分散推定へのフォールバック
if (curMS.sampleCount < 4) {
    constexpr float filterKernel[] = {
        0.00598, 0.060626, 0.241843, 0.383103, 0.241843, 0.060626, 0.00598
    };
    constexpr float centerWeight = filterKernel[3] * filterKernel[3];
    float sumFirstMoments = centerWeight * firstMoment;
    float sumSecondMoments = centerWeight * secondMoment;

    float depth = depthBuffer[pix];
    float3 normal = normalBuffer[pix];
    ... // 実際にはもう少し付随情報を読む。

    float sumWeights = 0.0f;
    for (int i = -3; i <= 3; ++i) {
        float hy = filterKernel[i];
        for (int j = -3; j <= 3; ++j) {
            if (i == 0 && j == 0)
                continue;

            float hx = filterKernel[j];
            
            int2 nbPix = int2(pix.x + j, pix.y + i);
            float nbDepth = depthBuffer[nbPix];
            float3 nbNormal = normalBuffer[nbPix];

            // モーメントのフィルタリングにもノイズ除去自体と同様のedge-stopping関数を使う。
            float wz = calcDepthWeight(depth, nbDepth, ...);
            float wn = calcNormalWeight(normal, nbNormal, ...);
            float weight = hx * hy * wz * wn;

            MomentPair_SampleCount nbMS = cur_momentPair_sampleCount_buffer[nbPix];
            sumFirstMoments = weight * nbMS.firstMoment;
            sumSecondMoments = weight * nbMS.secondMoment;
            sumWeights += weight;
        }
    }
    firstMoment = sumFirstMoments / sumWeights;
    secondMoment = sumSecondMoments / sumWeights;
}

float variance = secondMoment - firstMoment * firstMoment;
varianceBuffer[pix] = variance;

累積移動平均と指数移動平均

上記ではモーメントを求める際に以下のような単純な平均(累積移動平均, Cumulative Moving Average, CMA)を用いました。

float curWeight = 1.0f / sampleCount;
float avg = (1 - curWeight) * prevData + curWeight * curData;

が、これではアニメーションしているケースでも過去のフレームを延々とひきずることになり好ましくありません。そこで代わりに次のような指数移動平均(Exponential Moving Average, EMA)を用います。

float curWeight = 0.2f;
float avg = (1 - curWeight) * prevData + curWeight * curData;

0.2という数字は実際には(0, 1)の範囲のパラメターです。数字が小さいほどアニメーションがない場合には高品質な結果に収束していきますが、一方で応答性が悪くなります。SVGFの論文では0.2を用いているようです。EMAは過去フレームのウェイトが指数関数的に小さくなっていくのであまり古い情報を引きずらない利点がありますが、一方でサンプル数が極めて少ない場合にも新たなサンプルのウェイトが0.2などにしかならないため、これはこれで応答性が悪くなります。そこでCMAとEMAのハイブリッドな手法も使われます。

constexpr uint32_t threshold = 5;
float curWeight = 1.0f / threshold; // Exponential Moving Average
if (sampleCount < threshold) // Cumulative Moving Average
    curWeight = 1.0f / sampleCount;
float avg = (1 - curWeight) * prevData + curWeight * curData;

以下はシーン中をうさぎのモデルが前後に往復アニメーションしている際の最終レンダリングとサンプル数、分散を可視化したものです。手前から奥へとうさぎが進んだ後、折り返して手前に戻り始めた直後をキャプチャーしています。サンプル数は1~255をピクセルの明るさに変換(255以上は真っ白)、分散は厳密には標準偏差を適当なスケールをかけて可視化しています。うさぎの向こう側はDisocclusionが発生した直後でサンプル数が少ないことがわかります。分散の可視化中で対応する箇所を見てみると空間的な分散推定へフォールバックしたことにより、画像の他の部分に比べてぼやけていることがわかります。

À-trousフィルター

分散の大きいところでは、実質のサンプル数を稼ぐためフィルタリング半径を大きくしたいのですが単純に大きいフィルター半径を使うと非常に重くなるため、SVGFではÀ-trousフィルターと呼ばれるものを使用します。À-trousフィルターは多段階のフィルターとなっており上図のように、各段階でサンプルする周辺ピクセルの数は変えず、サンプルする距離を前段階の2倍に増加させます。つまりフィルターの半径中でサンプルしない点が数多くあります(À-trousはフランス語で英語だと"with holes"と意味らしい)。フィルターカーネル自体は論文では5x5のガウス(風？)フィルターが提案されていますが、3x3のボックスフィルターでも十分高品質かつ高速に動作しました。
以下にコード片を掲載します。前述の分散推定と似ていますが、隣接ピクセルへの距離がÀ-trousフィルターの段階に応じて変わるのと、edge-stopping関数にSVGFのコアである輝度に応じたウェイトが含まれています。関数の詳細は論文を参照してください。また分散も直接使うのではなく、安定化のため3x3の範囲でフィルタリングしたものを使うことが提案されています。

RGB noisyLighting = noisyLightingBuffer[pix];
float variance = varianceBuffer[pix];
float luminance = calcLuminance(noisyLighting);
float depth = depthBuffer[pix];
float3 normal = normalBuffer[pix];
... // 実際にはもう少し付随情報を読む。

// 安定化のために3x3の範囲で分散をフィルターにかける。
for (int i = -1; i <= 1; ++i) {
    for (int j = -1; j <= 1; ++j) {
        ...
    }
}
float localMeanStdDev = std::sqrt(localMeanVars);

// ここではKernelは2Dの全要素を含んだ1次元の配列で実装した。
constexpr float weights[] = {
    1, 1, 1,
    1, 1, 1,
    1, 1, 1,
};
constexpr int2 offsets[] = {
    int2(-1, -1), int2( 0, -1), int2( 1, -1),
    int2(-1,  0), int2( 0,  0), int2( 1,  0),
    int2(-1,  1), int2( 0,  1), int2( 1,  1),
};
constexpr uint32_t kernelSize = 9;
constexpr uint32_t centerIndex = kernelSize / 2;

constexpr float centerWeight = weights(centerIndex);
RGB sumLightings = centerWeight * noisyLighting;
float sumVars = centerWeight * centerWeight * variance;
float sumWeights = 0.0f;
for (int i = 0; i < kernelSize; ++i) {
    if (i == centerIndex)
        continue;

    // 隣接ピクセルへのオフセットはステージごとに倍になる。
    float h = weights[i];
    int2 offset = offsets[i] << filterStageIndex;
    
    int2 nbPix = int2(pix.x + offset.x, pix.y + offset.y);
    float nbDepth = depthBuffer[nbPix];
    float3 nbNormal = normalBuffer[nbPix];
    RGB nbNoisyLighting = noisyLightingBuffer[nbPix];
    float nbLuminance = calcLuminance(nbNoisyLighting);
    float nbVar = varianceBuffer[nbPix];

    float wz = calcDepthWeight(depth, nbDepth, ...);
    float wn = calcNormalWeight(normal, nbNormal);
    float wl = calcLuminanceWeight(luminance, nbLuminance, localMeanStdDev);

    float weight = h * wz * wn * wl;

    sumLightings += weight * nbNoisyLighting;
    sumVars += weight * weight * nbVar;
    sumWeights += weight;
}

RGB filteredLighting = sumLightings / sumWeights;
float filteredVar = sumVars / (sumWeights * sumWeights);

filteredLightingBuffer[pix] = filteredLighting;
filteredVarianceBuffer[pix] = filteredVar;

Temporal AccumulationとTemporal Anti-Aliasing

SVGFによる視覚的なノイズの低減の威力はなかなか素晴らしいですが、実際にはピュアにSVGFだけではそこまでの効果を発揮できません。SVGFのコアとはほぼ直交する概念であるTemporal AccumulationとTemporal Anti-Aliasingを組み合わせることで安定したレンダリング結果を得ることができます。SVGFの論文では安定化のため、Temporal Accumulationで読み出す前フレームのノイジーなライティング結果として、À-trousの1段階目の出力を使用することが提案されています。
これらのテクニックと組み合わせたSVGFのパイプラインは次のようになります。

G-Bufferの生成。
パストレーシングなどによるレンダリング
Temporal Accumulation
処理方法は前フレームのモーメント・サンプルカウントをモーションベクターを使って取得したのと同じです。
各ピクセルの分散を推定
À-trousフィルターを何度か繰り返す
1段階目の出力を次フレームのTemporal Accumulationで使えるように書き出しておく。
Temporal Anti-Aliasing

以下に1sppのパストレーシング、(+ Temporal Accumulation)、(+ SVGF)、(+ Temporal AA)の画像を載せます。各段階で徐々に品質が向上していっていることがわかるでしょう。(本当は動画比較もすべきかもしれませんが)

その他

直接照明と間接照明の分離

論文ではSVGFのパイプラインを直接照明成分と間接照明成分で分離することが提案されています。成分ごとに適切な分散推定やフィルタリングができるので、結果の品質は向上するのですが、デノイズのコストが単純に倍近くになるのでコストパフォーマンスは非常に悪い印象があります。これをやるくらいならDiffuseとSpecularに分けて何か別々の処理を加えたほうが良さそうです。

Path Space Regularization

名前が仰々しいですが、Specular-Specularな光輸送経路などはパストレーシングでは高い分散になりがちなので、セカンダリーヒット以降はサーフェスのsmoothnessを下げる(例えば半分)ことで分散を抑えます。もちろんこれによってレンダリングの結果は変わってしまうのですが、二次反射以降なので大きくは目立たないことが期待されています。

制約

確率的な一次の可視性

カメラから見て一次の可視点において分散推定を行っているため、そもそも一次の可視点すら確率的になるようなエフェクト、つまりDepth of Field (被写界深度)やモーションブラー、Stochastic Transparencyと併せて使うことができません。

スペキュラー反射

アルベド除去やモーションベクターはプライマリーヒットに関する情報であるため、シャープな鏡面に映る像がかなりボケてしまったり、カメラ移動時に過去フレームの情報を引きずりゴースト現象が出ます。

光源環境の変化

SVGF(というよりNaiveなTemporal Accumulation)では光源の移動や新たな光源の出現・消失などを考慮に入れていないため、光源環境が変化した場合にレンダリング結果が即座に変化するのではなく、遅れて徐々に変化するようになります。SVGFの発展手法の多くはこの部分に手を加えて応答性を改善しています。

相関のあるサンプルの取り扱い

昨今ReSTIRを始めとして、そもそもデノイズよりも前のレンダリングステップの時点で時間・空間的なサンプルの再利用を行い品質向上を図る手法が数多く提案されています。この場合各フレームで得られるサンプルが前フレームと相関を持っているため、分散を同じ手順で求めると適切ではないフィルタリング半径などを使ってしまいます。

さいごに

SVGFを改善した手法やSVGFに基づいた手法へのリンクを示して、この記事を終わりにします。

Adaptive-SVGF
Gradient Estimation for Real-time Adaptive Temporal Filtering
一部のピクセルで前フレームと同じ乱数列を用いてライトトランスポートパスの時間的な変化を求め、Temporal Accumulationのブレンド率を制御することでゴーストを抑えます。
SLGF
Spatiotemporal Luminance-Guided Filter
A-SVGFほどは複雑な実装は行わずに光源環境への追従を考慮。
ReLAX
ReLAX: A Denoiser Tailored to Work with the ReSTIR Algorithm
オリジナルのSVGFと比較して、デノイズ対象のDiffuseとSpecular信号への分離、Temporal Accumulationの改善、Disocclusionのケア、ファイアフライの明示的な除去、フィルタリングにおけるウェイトの追加等さまざまな改善が加えられています。

レンダリング界隈ではサンプルサイズではなくサンプル数と呼ぶ誤用慣習があるみたいなのでサンプル数と呼びます。 ↩