この記事はレイトレ合宿4!?アドベントカレンダー5週目の記事になります。

レイトレ合宿!?

レイトレ合宿は、お手製のレンダラーを持ち寄る集会です。
今回初参加させてもらうことになりました。よろしくお願いします。
残り１か月ほどで、そろそろ時間が無くなってきましたが、
なかなか作業が進まずやきもきしている次第です。

そんな中！さらに悪いことに
ポケモンgoがリリースされてしまいました。

この記事を書き終えたら、外に出てポケモンを捕まえに行こうと思います。
執筆時、
プレイヤーレベルは７、自分の持っている一番強いモンスターは、ラッタ 273cpです。
しばらくは走り回って健康になろうと思います。

Non-local Means Filter

と余談はこれくらいにして、今回はノイズ除去フィルタとして、Non-local Means、NLMとも略されますが、
このフィルタを試してみようかと思います。
というのも、モンテカルロパストレーシングにおいて、ノイズは死活問題です。
さらに今回、レンダリング制限時間が5分ということで、なかなかにシビアな設定になっています。
PCの性能も上がってきているとはいえ、なかなかきついです。
また、一般にモンテカルロ積分の収束は後半は緩やかになり、なかなか収束しません。
そこで出来上がった画像に後処理としてノイズ除去フィルタを適用することで、最終の画像の品質を底上げすることが可能です。
そこでノイズフィルタリングのとっかかりとして、

OpenCV.jpの解説を参考にしつつ、紹介します。
http://opencv.jp/opencv2-x-samples/non-local-means-filter

ノイズ除去基本戦略

ノイズフィルタ全般にいえることですが、
この手の処理の基本構造として、周囲のピクセルをなんらかの重みでブレンドする、というのがあります。
たとえば、3x3ガウスぼかしなら、


\begin{pmatrix}
\frac{1}{16} & \frac{2}{16} & \frac{1}{16}\\ 
\frac{2}{16} & \frac{4}{16} & \frac{2}{16}\\ 
\frac{1}{16} & \frac{2}{16} & \frac{1}{16}
\end{pmatrix}

の重みで周囲のピクセルを畳み込むことになります。
で、
この重みをどのような重みにするか？
というのが最大の焦点になります。
また、もう一つ忘れてはいけないこととして、重みの合計は1になる、というのは忘れないようにしたいです。
重みの和が1でなくなると、元の画像より暗くなったり、明るくなったりしてしまいます。
（もちろん、重みの和が1なら厳密に画像全体明るさが保存するかというと、そうでもないですが・・・）

このため、うまく重みを決めないと、ただ画像がぼやけるだけという事態にも陥ってしまうリスクがあります。

NLMの重み

NLMの重みの考え方は、

周囲の絵が、真ん中の絵に似ていれば、それは強く混ぜる。
周囲の絵が、真ん中の絵と離れていれば、それはあまり混ぜない。

※真ん中の絵=注目画素周辺, 周囲の絵=注目画素の周囲ピクセル周辺

というものです。
この似たもの同士でのみブレンドする考え方により、ただ画像がボケるのではなく、
エッジがちゃんと残ります。
また類似度を使う分、しっかり模様があるような画像のほうが得意だとあります。

絵がどの程度似たもの同士であるか？というのは、

１、一定のブロック、例えば5x5の小さな画素集合を注目位置から切り出す
２、それのL2ノルムを比較する

という方法をとります。
L2ノルムというのは、N次元ベクトルの距離です。
5x5のrgbなら、5x5x3 = 75次元ベクトルの距離ということですね。
これについては、こちらがわかりやすかったです。（http://mathtrain.jp/lpnorm）

だいたいの概要をつかんだところで、NLMの重みの数式をみます。
NLMの重みは以下の式で計算されます。


w(p, q) = \frac{1}{Z(p)} exp(-\frac{max(\lVert v(p) - v(q))) \rVert_{L_{2}}^2 - 2\sigma^2, 0))}{h^2})

σとhはパラメータです。
これはノイズの標準偏差をいれると効果的だそうです
まず下の式は忘れて、上の式のexp以降を見ましょう。

pが注目画素
qが周辺画素
v(q)が周辺画素の、類似度のためのベクトル(5x5のrgbなら、5x5x3 = 75次元ベクトル)
v(p)が注目画素の、類似度のためのベクトル(5x5のrgbなら、5x5x3 = 75次元ベクトル)


\lVert v(p) - v(q))) \rVert_{L_{2}}

の部分が先ほどの類似度になります。
※L2が重なってみづらいですが、式の中の右肩の2はべき乗です

おっと、Z(p)を説明していませんでした。
これは以下の式になります。


Z(p) = \sum_{ q\in S }^{ } exp(-\frac{max(\lVert v(p) - v(q))) \rVert_{L_{2}}^2 - 2\sigma^2, 0))}{h^2})

おや？ほとんど最初の式と一緒ですね？
というのも、これは重みの和を1にするための正規化係数だからです。全部足して、それで割る。
ただそれだけです。

実装

あとはコードにするだけです。
これをcinder上で実装しました。
ややcinder apiの癖がありますが、おおむね読めるくらいにはなっていると思います。

ただ・・・ちょっと類似度の計算でネストが深いため、
実行速度にやや難があります。


inline cinder::Surface32fRef non_local_means(cinder::Surface32fRef image, float param_h, float sigma) {
    param_h = std::max(0.0001f, param_h);
    sigma = std::max(0.0001f, sigma);

    const int kKernel = 5;
    const int kSupport = 13;
    const int kHalfKernel = kKernel / 2;
    const int kHalfSupport = kSupport / 2;

    int width = image->getWidth();
    int height = image->getHeight();
    auto surface = cinder::Surface32f::create(width, height, false);
    float *src = image->getData();
    float *dst = surface->getData();

    typedef std::array<float, 3 * kKernel * kKernel> Template;

    auto sample_template = [width, height, src, kHalfKernel](int x, int y) {
        Template t;
        int index = 0;
        for (int sx = x - kHalfKernel; sx <= x + kHalfKernel; ++sx) {
            for (int sy = y - kHalfKernel; sy <= y + kHalfKernel; ++sy) {
                int sample_x = sx;
                int sample_y = sy;
                sample_x = std::max(sample_x, 0);
                sample_x = std::min(sample_x, width - 1);
                sample_y = std::max(sample_y, 0);
                sample_y = std::min(sample_y, height - 1);

                float *p = src + (sample_y * width + sample_x) * 3;
                for (int i = 0; i < 3; ++i) {
                    t[index++] = p[i];
                }
            }
        }
        return t;
    };
    auto distance_sqared_template = [](const Template &a, const Template &b) {
        float accum = 0.0f;
        for (int i = 0; i < a.size(); ++i) {
            accum += glm::pow(a[i] - b[i], 2);
        }
        return accum;
    };
    concurrency::parallel_for<int>(0, height, [=](int y) {
        for (int x = 0; x < width; ++x) {
            auto sample = [=](int x, int y) {
                int sample_x = x;
                int sample_y = y;
                sample_x = std::max(sample_x, 0);
                sample_x = std::min(sample_x, width - 1);
                sample_y = std::max(sample_y, 0);
                sample_y = std::min(sample_y, height - 1);

                float *p = src + (sample_y * width + sample_x) * 3;
                return vec3(p[0], p[1], p[2]);
            };

            float *dst_pixel = dst + (y * width + x) * 3;

            auto focus = sample_template(x, y);

            dvec3 sum;
            double sum_weight = 0.0;
            for (int sx = x - kHalfSupport; sx <= x + kHalfSupport; ++sx) {
                for (int sy = y - kHalfSupport; sy <= y + kHalfSupport; ++sy) {
                    auto target = sample_template(sx, sy);
                    auto dist = distance_sqared_template(focus, target);
                    auto arg = -glm::max(dist - 2.0f * sigma * sigma, 0.0f) / (param_h * param_h);
                    auto weight = glm::exp(arg);

                    sum_weight += weight;
                    sum += sample(sx, sy) * weight;
                }
            }
            auto color = sum / sum_weight;

            for (int i = 0; i < 3; ++i) {
                dst_pixel[i] = color[i];
            }
        }
    });

    return surface;
}