Ray Histogram Fusion

概要

Ray Histogram Fusion(RHF)フィルターは、モンテカルロ積分で得たイメージからノイズを除去するフィルターです。

NL-meanフィルタなどのノイズ除去フィルタは、入力にノイズを除去したい画像をそのまま使います。着目しているピクセルの周囲で、似たようなパターンを探し、類似していると認められたものに強くウェイトを掛け、それらの平均を取ることでノイズを除去します。どれ程似ているかを類似度、パターンの単位をパッチと呼びます。

一方、RHFはモンテカルロ積分中に得られたサンプルをヒストグラム化します。類似度はパッチ間のヒストグラムの$\chi^2$距離を利用します。類似度が高いパッチにあるピクセルを、着目しているピクセルの新しいサンプルと考えて再利用します。

デモ

スライダーで、フィルターされる前のイメージと、RHFフィルターを掛けた後のイメージの比較ができます。

周囲のサンプルを自身のサンプルとしてどの程度取り込んだかを可視化すると次のようになります。

オブジェクトのエッジや、ライティング環境が変化する境界付近での追加で得られるサンプル数は減り、光源の映りこみや同一平面が続く場合などはサンプル数が稼げていることが分かります。

実装メモ

  • 「$\chi^2$距離」とは言っていますが、やっていることは$\chi^2$検定です。
  • 直接光(ライトや環境マップ)や、それの完全反射などは追加でサンプルする必要性がないのでスキップする最適化があるかもしれません。
  • パラメーター調整ゲームであることは否めない気はしています。まずヒストグラムを綺麗に使いきれた方が結果は良いですし、閾値kは結果にダイレクトに影響するので丹念に調整する必要があります。
  • 閾値kは「これ以下の差は同じ結果になると思い込もう」というパラメーターなので、大きすぎると当然全てのピクセルを採用してしまい、ただのブラーになります。0にすれば何も取り込まないので何もしないフィルターになります。
  • 十分な情報がある状態、つまりそこそこレンダリング済みではないと、検定の精度は上がりません。このフィルターがやっていることは結局「あとN倍の時間を掛けたとしての結果の近似を得る」ということなので数秒しかレンダリングしていない結果に適用するよりも、何時間、何十時間と掛けた結果に使用して、何百時間掛けたときの結果と似た結果を得るという使い方が正しい気がしています。
  • 参考に上げた大村平氏の書籍は、私が読んだことのある確率統計の入門書の中では群を抜いて分かりやすく、ただの入門書では語られない原理や事実を取り上げていてとても面白いです。(必要なのに)確率統計を避けてきた全ての人にお勧めします。(大村平氏の確率統計関連の書籍は全て面白いので強く推しておきます)。

参考

[1] M. Delbracio, P. Musé, A. Buades, J. Chauvier, N. Phelps, and JM. Morel, “Boosting monte carlo rendering by ray histogram fusion,” ACM Trans. Graph., vol. 33, no. 1, pp. 1–15, 2014.
[2] mdelbra/rhf. https://github.com/mdelbra/rhf
[3] 大村平, 統計のはなし―基礎・応用・娯楽, 2002
[4] 大村平, 統計解析のはなし―データに語らせるテクニック , 2006
[5] 東京大学教養学部統計学教室, 統計学入門, 1991

コード

{% highlight C %} #include #include #include #include #include #include #include #include #include

// #define DEBUG_READ_HIST // #define DEBUG_WRITE_HIST // #define DEBUG_WRITE_RESULT_HDR

#if defined(DEBUG_WRITE_HIST) || defined(DEBUG_READ_HIST) #pragma warning(disable:4996) #define STB_IMAGE_WRITE_IMPLEMENTATION #include “stb_image_write.h” #define STB_IMAGE_IMPLEMENTATION #include “stb_image.h” #endif

const float PI = 3.141592653589f; // float rnd() { static uint64_t x = 123456789; x ^= x » 12; x ^= x « 25; x ^= x » 27; const uint64_t t = x * UINT64_C(2685821657736338717); const float inv = 1.0f / (float)std::numeric_limits<uint64_t>::max(); return (float)t * inv; } // class Vec3 { public: union { struct { float x; float y; float z; }; struct { float r; float g; float b; }; std::array<float, 3> e; }; public: Vec3() :x(0.0f), y(0.0f), z(0.0f) {} ~Vec3() {} Vec3(const Vec3& other) :x(other.x), y(other.y), z(other.z) {} Vec3(float v) :x(v), y(v), z(v) {} Vec3(float ax, float ay, float az) :x(ax), y(ay), z(az) {} float operator [](int32_t index) const { return e[index]; } float& operator [](int32_t index) { return e[index]; } }; static Vec3 operator + (const Vec3& lhs, const Vec3& rhs) { return Vec3( lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z); } static Vec3 operator - (const Vec3& lhs, const Vec3& rhs) { return Vec3( lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z); } static Vec3 operator * (float s, const Vec3& v) { return Vec3(sv.x, sv.y, sv.z); } static Vec3 operator * (const Vec3& v, float s) { return sv; } static Vec3& operator += (Vec3& lhs, const Vec3& rhs) { lhs.x += rhs.x; lhs.y += rhs.y; lhs.z += rhs.z; return lhs; } static float dot(const Vec3& lhs, const Vec3& rhs) { return lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z; } static Vec3 cross(const Vec3& lhs, const Vec3& rhs) { return Vec3( lhs.yrhs.z - lhs.zrhs.y, lhs.zrhs.x - lhs.xrhs.z, lhs.xrhs.y - lhs.yrhs.x); } static Vec3 mul(const Vec3& lhs, const Vec3& rhs) { return Vec3(lhs.x * rhs.x, lhs.y * rhs.y, lhs.z * rhs.z); } static float lenSq(const Vec3& v) { return dot(v, v); } static float len(const Vec3& v) { return sqrtf(lenSq(v)); } static Vec3& saturate(Vec3& v) { v.x = v.x > 1.0f ? 1.0f : v.x < 0.0f ? 0.0f : v.x; v.y = v.y > 1.0f ? 1.0f : v.y < 0.0f ? 0.0f : v.y; v.z = v.z > 1.0f ? 1.0f : v.z < 0.0f ? 0.0f : v.z; return v; } static Vec3 norm(const Vec3& v) { const float s = 1.0f / len(v); return Vec3(s * v.x, s * v.y, s * v.z); } typedef Vec3 Color;

// class Image { public: Image(int32_t w, int32_t h) :width_(w), height_(h), pixels_(w*h) {} const std::vector& pixel() const { return pixels_; } std::vector& pixel() { return pixels_; } const Color& pixel(int32_t x, int32_t y)const { return pixels_[x + width_*y]; } Color& pixel(int32_t x, int32_t y) { return pixels_[x + width_*y]; } void save(const char* filename) { const auto degamma = [](float x) { x = x > 1.0f ? 1.0f : x < 0.0f ? 0.0f : x; return (int32_t)(powf(x, 1.0f / 2.2f) * 255.0f + 0.5f); }; FILE* file = NULL; fopen_s(&file, filename, “w”); fprintf(file, “P3\n%d %d\n%d\n”, width_, height_, 255); for (const auto& p : pixels_) { fprintf(file, “%d %d %d\n”, degamma(p.x), degamma(p.y), degamma(p.z)); } fclose(file); } Color& operator [](int32_t idx) { return pixels_[idx]; } private: int32_t width_; int32_t height_; std::vector pixels_; }; // struct Ray { public: Vec3 o; Vec3 d; float t; public: Ray(const Vec3& ao, const Vec3& ad) :o(ao), d(ad) { t = std::numeric_limits::max(); } }; // enum MaterialType { DIFF, SPEC, REFR }; // struct Sphere { public: float r; Vec3 pos; Vec3 emission; Vec3 diffuse; MaterialType matType; public: Sphere(float ar, const Vec3& apos, const Vec3& aemission, const Vec3& adiffuse, MaterialType amatType) : r(ar), pos(apos), emission(aemission), diffuse(adiffuse), matType(amatType) {} float intersect(const Ray &ray) const { const Vec3 op = pos - ray.o; float b = dot(op, ray.d); float det = b*b - dot(op, op) + r*r; if (det < 0.0f) { return 0; } det = sqrtf(det); const float eps = 1e-4f; float t = b - det; if (t > eps) { return t; } t = b + det; if (t > eps) { return t; } else { return 0.0f; } } }; // static void intersectScene(const Ray& ray, const Sphere** obj, float* t) { // static const Sphere spheres[] = { Sphere(1e3f, Vec3(1e3f + 1.0f, 40.8f, 81.6f), Vec3(0.0f), Vec3(0.75f, 0.25f, 0.25f), DIFF), Sphere(1e3f, Vec3(-1e3f + 99.0f, 40.8f, 81.6f), Vec3(0.0f), Vec3(0.25f, 0.25f, 0.75f), DIFF), Sphere(1e3f, Vec3(50.0f, 40.8f, 1e3f), Vec3(), Vec3(0.75f, 0.75f, 0.75f), DIFF), Sphere(1e3f, Vec3(50.0f, 40.8f, -1e3f + 170.0f), Vec3(), Vec3(), DIFF), Sphere(1e3f, Vec3(50.0f, 1e3f, 81.6f), Vec3(), Vec3(0.75f, 0.75f, 0.75f), DIFF), Sphere(1e3f, Vec3(50.0f, -1e3f + 81.6f, 81.6f), Vec3(), Vec3(0.75f, 0.75f, 0.75f), DIFF), Sphere(16.5f, Vec3(27.0f, 16.5f, 47.0f), Vec3(), Vec3(1.0f, 1.0f, 1.0f)*0.999f, SPEC), Sphere(16.5f, Vec3(73.0f, 16.5f, 78.0f), Vec3(), Vec3(1.0f, 1.0f, 1.0f)*0.999f, REFR), Sphere(600.0f, Vec3(50.0f, 681.6f - 0.27f, 81.6f), Vec3(12.0f, 12.0f, 12.0f), Vec3(), DIFF) }; // *t = std::numeric_limits ::max(); for (const auto& s : spheres) { const float nt = s.intersect(ray); if (nt == 0.0f) { continue; } if (nt < *t) { *t = nt; *obj = &s; } } } // static Vec3 radiance(const Ray& ray, int depth) { const Sphere* obj = NULL; float t; intersectScene(ray, &obj, &t); if (!obj) { return Vec3(0.0f); } const Vec3 x = ray.o + ray.d*t; const Vec3 n = norm(x - obj->pos); const Vec3 nl = dot(n, ray.d) < 0 ? n : -1.0f * n; Vec3 f = obj->diffuse; ++depth; if (depth > 50) { return Vec3(0.0f); } else if (depth > 5) { const float p = f.x > f.y && f.x > f.z ? f.x : f.y > f.z ? f.y : f.z; if (rnd() < p) { f = f*(1 / p); } else { return obj->emission; } } if (obj->matType == DIFF) { const float r1 = 2 * PI*rnd(); const float r2 = rnd(); const float r2s = sqrt(r2); const Vec3 w = nl; const Vec3 u = norm(cross(fabs(w.x) > 0.1f ? Vec3(0.0f, 1.0f, 0.0f) : Vec3(1.0f, 0.0f, 0.0f), w)); const Vec3 v = cross(w, u); const Vec3 d = norm(u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1 - r2)); return obj->emission + mul(f, (radiance(Ray(x, d), depth))); } else if (obj->matType == SPEC) { Vec3 nd = (ray.d - n * 2.0f * dot(n, ray.d)); const Ray nray(x, nd); const Color nr = radiance(nray, depth); return obj->emission + mul(f, nr); } Ray reflRay(x, ray.d - n * 2 * dot(n, ray.d)); bool into = dot(n, nl) > 0.0f; const float nc = 1.0f; const float nt = 1.5f; const float nnt = into ? nc / nt : nt / nc; const float ddn = dot(ray.d, nl); const float cos2t = 1.0f - nnt*nnt*(1 - ddn*ddn); if (cos2t < 0.0f) { return obj->emission + mul(f, radiance(reflRay, depth)); } const Vec3 tdir = norm(ray.d*nnt - n*((into ? 1.0f : -1.0f)*(ddn*nnt + sqrtf(cos2t)))); const float a = nt - nc; const float b = nt + nc; const float R0 = a*a / (b*b); const float c = 1.0f - (into ? -ddn : dot(tdir, n)); const float Re = R0 + (1.0f - R0)*c*c*c*c*c; const float Tr = 1.0f - Re; const float P = 0.25f + 0.5f*Re; const float RP = Re / P; const float TP = Tr / (1 - P); return obj->emission + mul(f, depth > 2 ? (rnd() < P ? radiance(reflRay, depth)*RP : radiance(Ray(x, tdir), depth)*TP) : radiance(reflRay, depth)*Re + radiance(Ray(x, tdir), depth)*Tr); }

// M. Delbracio, P. Musé, A. Buades, J. Chauvier, N. Phelps, and JM. Morel, “Boosting monte carlo rendering by ray histogram fusion,” ACM Trans. Graph., vol. 33, no. 1, pp. 1–15, 2014. class RHFFilter { public: // RHFFilter( int32_t width, int32_t height, float k, int32_t patchSizeHalf, int32_t windowSizeHalf, int32_t numBins = 20, float M = 7.5f, float s = 2.0f, const float gamma = 2.2f) :width_(width), height_(height), patchHalf_(patchSizeHalf), windowHalf_(windowSizeHalf), numBins_(numBins), invM_(1.0f / M), s_(s), invGamma_(1.0f / gamma), k_(k) { pixelInfos.resize(width*height); for (auto& p : pixelInfos) { p.histgram[0].resize(numBins, 0.0f); p.histgram[1].resize(numBins, 0.0f); p.histgram[2].resize(numBins, 0.0f); p.totalWeight = 0.0f; } } // void addPixel(int32_t x, int32_t y, const Color& color) { const int32_t index = x + y*width_; auto& pi = pixelInfos[index]; const float weight = 1.0f; pi.totalWeight += weight; for (int32_t chNo = 0; chNo < 3; ++chNo) { float v = weight * color[chNo]; v = std::max(0.0f, v); v = powf(v, invGamma_) * invM_; v = std::min(v, s_); const float fbin = v * ((float)numBins_ - 2.0f); const float ibinL = floor(fbin); // float wH; int32_t ibL; int32_t ibH; float wbL; float wbH; if ((int32_t)ibinL < numBins_ - 2) { wH = fbin - ibinL; ibL = (int32_t)ibinL; wbL = 1.0f - wH; ibH = (int32_t)ibinL + 1; wbH = wH; } else { wH = (v - 1.0f) / (s_ - 1.0f); ibL = numBins_ - 2; wbL = 1.0f - wH; ibH = numBins_ - 1; wbH = wH; } auto& hist = pi.histgram[chNo]; hist[ibL] += wbL; hist[ibH] += wbH; } // pi.accumColor = pi.accumColor + color; ++pi.numAddPixel; } // void singleScaleRayHistgramFusion(const std::string& fileName ) { // Image inputImage(width_, height_); Image outImage(width_, height_); auto& dst1 = inputImage.pixel(); //auto& dst2 = outImage.pixel(); for (uint32_t pi = 0; pi < dst1.size(); ++pi) { const auto& sp = pixelInfos[pi]; const float iw = 1.0f / sp.totalWeight; const Color col = sp.accumColor * iw; dst1[pi] = col; //dst2[pi] = col; } // const int32_t windowSize = windowHalf_ * 2 + 1; const float weight = 1.0f; // Color v(0.0f); std::vector n; n.resize(pixelInfos.size(), 0); std::vector count; count.resize(pixelInfos.size(), 0); // for (uint32_t pi = 0, pn = pixelInfos.size(); pi < pn; ++pi) { // if (!(pi % (pn/100))) { printf("%f\n", (float)pi / (float)pixelInfos.size()); } // const int32_t x = pi % width_; const int32_t y = pi / width_; // float v = 0.0f; float c = 0.0f; float totalWeight = 0.0f; // const int32_t patchSize = patchHalf_ * 2 + 1; const int32_t patchSize2 = patchSize * patchSize; const int32_t patchHalfAdj = std::min(patchHalf_, std::min(std::min(width_ - 1 - x, height_ - 1 - y), std::min(x, y))); const int32_t patchSizeAdj = patchHalfAdj * 2 + 1; const int32_t patchSizeAdj2 = patchSizeAdj * patchSizeAdj; // const int32_t wxMin = std::max(x - windowHalf_, patchHalfAdj); const int32_t wyMin = std::max(y - windowHalf_, patchHalfAdj); const int32_t wxMax = std::min(x + windowHalf_, width_ - 1 - patchHalfAdj); const int32_t wyMax = std::min(y + windowHalf_, height_ - 1 - patchHalfAdj); // std::vector denoised; denoised.resize(patchSize2, Color(0.0f)); // for (int32_t wy = wyMin; wy < wyMax; ++wy) { for (int32_t wx = wxMin; wx < wxMax; ++wx) { if (wx == x && wy == y) { continue; } // float dx2PxPy = 0.0f; float dx2PxPyNum = 0.0f; for (int32_t chNo = 0; chNo < 3; ++chNo) { for (int32_t pxy = 0; pxy < patchSizeAdj2; ++pxy) { const int32_t px = pxy % patchSizeAdj - patchHalfAdj; const int32_t py = pxy / patchSizeAdj - patchHalfAdj; int32_t nx = x + px; int32_t ny = y + py; int32_t ntx = wx + px; int32_t nty = wy + py; // assert(0 <= nx && nx < width_); assert(0 <= ny && ny < height_); assert(0 <= ntx && ntx < width_); assert(0 <= nty && nty < height_); // const auto& pixelX = pixelInfos[nx + ny*width_]; const auto& pixelY = pixelInfos[ntx + nty*width_]; // float dx2CxCy = 0.0f; float dx2CxCyNum = 0.0f; for (int32_t bn = 0; bn < numBins_; ++bn) { float hx = pixelX.histgram[chNo][bn]; float hy = pixelY.histgram[chNo][bn]; const float nyxrt = sqrtf((float)pixelY.numAddPixel / (float)pixelX.numAddPixel); const float nxyrt = 1.0f / nyxrt; const float tmp = (nyxrt * hx - nxyrt * hy); const float nume = tmp*tmp; const float denom = hx + hy; if (denom <= 1.0f) { continue; } dx2CxCy += nume / denom; dx2CxCyNum += 1.0f; } dx2CxCy /= dx2CxCyNum; dx2PxPy += dx2CxCy; dx2PxPyNum += 1.0f; } } dx2PxPy /= dx2PxPyNum; // if (dx2PxPy < k_) { totalWeight += weight; for (int is = -patchHalfAdj; is <= patchHalfAdj; ++is) { for (int ir = -patchHalfAdj; ir <= patchHalfAdj; ++ir) { const int32_t idx = (patchHalf_ + is) * patchSize + patchHalf_ + ir; const auto& targetPix = inputImage.pixel(wx + ir, wy + is); denoised[idx] += weight * targetPix; } } // n[pi] += 1.0f; } } } // for (int is = -patchHalfAdj; is <= patchHalfAdj; ++is) { for (int ir = -patchHalfAdj; ir <= patchHalfAdj; ++ir) { const int32_t idx = (patchHalf_ + is) * patchSize + patchHalf_ + ir; const auto& targetPix = inputImage.pixel(x + ir, y + is); denoised[idx] += weight * targetPix; } } totalWeight += weight; // const float litlle = 0.00000001f; if (totalWeight > litlle ) { for (int32_t is = -patchHalfAdj; is <= patchHalfAdj; is++) { for (int32_t ir = -patchHalfAdj; ir <= patchHalfAdj; ir++) { const int32_t idx = (patchHalf_ + is) * patchSize + patchHalf_ + ir; const int32_t il = (y + is)*width_ + x + ir; count[il] += 1.0f; const float invTotalWeight = 1.0f / totalWeight; auto& outCol = outImage.pixel()[il]; outCol += denoised[idx] * invTotalWeight; } } } } // for (uint32_t pi = 0; pi < count.size(); ++pi) { auto& outCol = outImage.pixel()[pi]; const float c = count[pi]; if (c > 0.0f) { const float ic = 1.0f / c; outCol = outCol * ic; } else { outCol = inputImage.pixel()[pi]; } } // const auto& ps = outImage.pixel(); const std::string ppmFile = fileName + “.ppm”; outImage.save(ppmFile.c_str()); #if defined(DEBUG_WRITE_RESULT_HDR) const std::string hdrFile = fileName + “.hdr”; stbi_write_hdr(hdrFile.c_str(), width_, height_, 3, (float*)ps.data()); #endif // Image nImage(width_, height_); for (uint32_t pi = 0; pi < n.size(); ++pi) { const float windowSiz2 = (float)(windowSize*windowSize); nImage.pixel()[pi].r = n[pi] / windowSiz2; nImage.pixel()[pi].g = n[pi] > windowSiz2 ? 1.0f : 0.0f; } const std::string nFile = fileName + “_n.ppm”; nImage.save(nFile.c_str()); } #if defined(DEBUG_WRITE_HIST) // void writeHist(const std::string& fileName) { for (int32_t binNo = 0; binNo < numBins_; ++binNo) { // const std::string binName = fileName + “_” + std::to_string(binNo) + “.hdr”; std::vector pixels; for (const auto& p : pixelInfos) { const float r = p.histgram[0][binNo]; const float g = p.histgram[1][binNo]; const float b = p.histgram[2][binNo]; pixels.push_back(Color(r, g, b)); } stbi_write_hdr(binName.c_str(), width_, height_, 3, (float*)pixels.data()); } // Image inputImage(width_, height_); auto& dst = inputImage.pixel(); for (uint32_t pi = 0; pi < dst.size(); ++pi) { const auto& sp = pixelInfos[pi]; const float iw = 1.0f / sp.totalWeight; dst[pi] = sp.accumColor * iw; } const std::string mImageName = fileName + “_ave.hdr”; stbi_write_hdr(mImageName.c_str(), width_, height_, 3, (float*)inputImage.pixel().data()); } #endif #if defined(DEBUG_READ_HIST) // void loadHist(const std::string& fileName) { for (int32_t binNo = 0; binNo < numBins_; ++binNo) { const std::string binName = fileName + “_” + std::to_string(binNo) + “.hdr”; int32_t comp = 3; float* v = stbi_loadf(binName.c_str(), &width_, &height_, &comp, comp); pixelInfos.resize(width_ * height_); for (uint32_t pi = 0; pi < pixelInfos.size(); ++pi) { auto& pinfo = pixelInfos[pi]; auto& hist = pinfo.histgram; hist[0][binNo] = v[pi * 3 + 0]; hist[1][binNo] = v[pi * 3 + 1]; hist[2][binNo] = v[pi * 3 + 2]; } stbi_image_free(v); }

    //
    const std::string mImageName = fileName + "_ave.hdr";
    int32_t comp = 3;
    float* v = stbi_loadf(mImageName.c_str(), &width_, &height_, &comp, comp);
    for (uint32_t pi = 0; pi < pixelInfos.size(); ++pi)
    {
        auto& pix = pixelInfos[pi];
        auto& col = pix.accumColor;
        col.r = v[pi * 3 + 0];
        col.g = v[pi * 3 + 1];
        col.b = v[pi * 3 + 2];
        pix.totalWeight = 1.0f;
        pix.numAddPixel = 128;
    }
    stbi_image_free(v);
}

#endif private: typedef std::vector HistgramPerChannel; struct PixelInfo { std::array<HistgramPerChannel, 3> histgram; float totalWeight; Color accumColor; int32_t numAddPixel; }; std::vector pixelInfos; int32_t width_; int32_t height_; int32_t patchHalf_; int32_t windowHalf_; int32_t numBins_; float invM_; float s_; float invGamma_; float k_;

}; // void main() { const int32_t width = 800; const int32_t height = 450; const int32_t sampNum = 512; const float fw = (float)width; const float fh = (float)height; Image image(width, height); const int32_t patchSizeHalf = 1; const int32_t windowSizeHalf = 6; // or 9 const float k = 0.5f; const int32_t numBins = 20; const float M = 3.0f; RHFFilter rhf(width, height, k, patchSizeHalf, windowSizeHalf, numBins, M ); const Ray cam(Vec3(50.0f, 35.0f, 300.0f), Vec3(0.0f, 0.0f, -1.0f)); const float f = 0.36f; const float aspect = fw / fh; const Vec3 cx = Vec3(aspect*f, 0.0f, 0.0f); const Vec3 cy = norm(cross(cx, cam.d))*f;

#if defined(DEBUG_READ_HIST) rhf.loadHist(“src”); rhf.singleScaleRayHistgramFusion(“result”); return; #endif for (int32_t y = 0; y < height; y++) { printf("%d/%d\n", y, height); for (int32_t x = 0; x < width; x++) { const float invSampleNum = 1.0f / (float)sampNum; Vec3 r(0.0f); for (int32_t si = 0; si < sampNum; ++si) { const float sx = rnd(); const float sy = rnd(); const float fx = ((float)x + sx - 0.5f) / fw - 0.5f; const float fy = -((float)y + sy - 0.5f) / fh + 0.5f; const Vec3 d = fxcx + fycy + cam.d; const Ray ray(cam.o + d * 140.0f, norm(d)); Color nr = radiance(ray, 0); rhf.addPixel(x, y, nr); r = r + nr * invSampleNum; } image.pixel(x, y) = saturate(r); } } image.save(“orig.ppm”); #if defined(DEBUG_WRITE_HIST) rhf.writeHist(“src”); #endif rhf.singleScaleRayHistgramFusion(“result”); }

{% endhighlight %}

自分用メモ

参照実装が少し読みづらいのでここにメモを残しておく。

パッチと呼んでいる物

  • デフォルトサイズは1
  • 平均される小さな領域
  • param->win(ウィンドウではないことに注意)
  • iDWin(パッチの半分サイズ) -> patchHalf_
  • iDWin0(調整済みパッチの半分) -> adjPatchHalf

ブロックと呼んでいるもの

  • デフォルトサイズは6
  • 探索する比較的大きな領域
  • param->bloc,blocなどが相当。
  • iDBloc(ブロックの半分サイズ) -> windowHalf_

読み替え

  • ihwl(パッチサイズ) -> patchSize
  • iwl(パッチ面積) -> patchSize
  • iDWin(パッチの半分サイズ) -> patchHalf_
  • iDWin0(調整済みパッチの半分) -> adjPatchHalf
  • iDBloc(ウィンドウの半分サイズ) -> windowHalf_