From c66c0841492a6a8134a8e9873f28154d8d6a5937 Mon Sep 17 00:00:00 2001 From: tatsy Date: Mon, 21 Mar 2016 17:31:32 +0900 Subject: [PATCH] Update hierarchical integrator to consider both indirect and direct lighting at sample points. --- example/pathtracing_example.cc | 4 +- sources/bxdf/bssrdf.cc | 2 +- sources/core/sampling.cc | 81 ++++++++ sources/core/sampling.h | 4 + sources/integrator/hierarchical.cc | 320 +++++++++++++++-------------- sources/integrator/hierarchical.h | 97 +-------- 6 files changed, 261 insertions(+), 247 deletions(-) diff --git a/example/pathtracing_example.cc b/example/pathtracing_example.cc index 0f650b9..54d09d8 100644 --- a/example/pathtracing_example.cc +++ b/example/pathtracing_example.cc @@ -36,8 +36,8 @@ int main(int argc, char **argv) { //VolPathIntegrator integr(camera, sampler); //PathIntegrator integr(camera, sampler); - IrradCacheIntegrator integr(camera, sampler); - //HierarchicalIntegrator integr(camera, sampler); + //IrradCacheIntegrator integr(camera, sampler); + HierarchicalIntegrator integr(camera, sampler); //PPMProbIntegrator integr(camera, sampler); integr.render(scene, params); diff --git a/sources/bxdf/bssrdf.cc b/sources/bxdf/bssrdf.cc index 422223b..8b4e08f 100644 --- a/sources/bxdf/bssrdf.cc +++ b/sources/bxdf/bssrdf.cc @@ -384,7 +384,7 @@ double DiffusionReflectance::Ft(const Vector3d& w) const { double DiffusionReflectance::Fdr() const { if (eta_ >= 1.0) { return -1.4399 / (eta_ * eta_) + 0.7099 / eta_ + 0.6681 + - 0.0636f * eta_; + 0.0636 * eta_; } else { return -0.4399f + 0.7099 / eta_ - 0.3319 / (eta_ * eta_) + 0.0636 / (eta_ * eta_ * eta_); diff --git a/sources/core/sampling.cc b/sources/core/sampling.cc index 7649d1d..46bcd95 100644 --- a/sources/core/sampling.cc +++ b/sources/core/sampling.cc @@ -184,6 +184,87 @@ void sampleUniformHemisphere(const Normal3d& normal, Vector3d* direction, const *direction = (u * cos(t) * z2s + v * sin(t) * z2s + w * sqrt(1.0 - z2)).normalized(); } +void samplePoissonDisk(const Scene& scene, const Point3d& pCamera, + double minDist, std::vector* points) { + std::vector tris; + for (const auto& p : scene.primitives()) { + if (p->material()->isSubsurface()) { + auto ts = p->triangulate(); + tris.insert(tris.end(), ts.begin(), ts.end()); + } + } + + // Initialize utility variables in this method + const auto seed = static_cast(time(0)); + Random rng(seed); + + // Sample random points on trimesh + Bounds3d bounds; + std::vector candPoints; + std::vector candNormals; + for (int i = 0; i < tris.size(); i++) { + const Triangle& tri = tris[i]; + const double A = tri.area(); + const int nSample = static_cast(std::ceil(4.0 * A / (minDist * minDist))); + for (int k = 0; k < nSample; k++) { + double u = rng.nextReal(); + double v = rng.nextReal(); + if (u + v >= 1.0) { + u = 1.0 - u; + v = 1.0 - v; + } + Point3d p = (1.0 - u - v) * tri[0] + u * tri[1] + v * tri[2]; + Normal3d n = (1.0 - u -v) * tri.normal(0) + u * tri.normal(1) + + v * tri.normal(2); + candPoints.push_back(p); + candNormals.push_back(n); + bounds.merge(p); + } + } + + // Create hash grid + const int numCands = static_cast(candPoints.size()); + Vector3d bsize = bounds.posMax() - bounds.posMin(); + const double scale = 1.0 / (2.0 * minDist); + const int numPoints = candPoints.size(); + HashGrid hashgrid; + hashgrid.init(numPoints, scale, bounds); + + RandomQueue que; + for (int i = 0; i < numCands; i++) { + que.push(i); + } + + std::vector sampledIDs; + Vector3d margin(2.0 * minDist, 2.0 * minDist, 2.0 * minDist); + while (!que.empty()) { + int id = que.pop(); + Point3d v = candPoints[id]; + const std::vector& cellvs = hashgrid[v]; + + bool accept = true; + for (int k = 0; k < cellvs.size(); k++) { + if ((candPoints[cellvs[k]] - v).squaredNorm() < minDist * minDist) { + accept = false; + break; + } + } + + if (accept) { + Point3d boxMin = v - margin; + Point3d boxMax = v + margin; + hashgrid.add(id, boxMin, boxMax); + sampledIDs.push_back(id); + } + } + + // Store sampled points + std::vector::iterator it; + for (auto i : sampledIDs) { + points->emplace_back(candPoints[i], candNormals[i]); + } +} + void samplePoissonDisk(const Scene& scene, const Point3d& pCamera, double minDist, std::vector* points) { // Initialize utility variables in this method diff --git a/sources/core/sampling.h b/sources/core/sampling.h index 8731fc1..8e2cb94 100644 --- a/sources/core/sampling.h +++ b/sources/core/sampling.h @@ -63,6 +63,10 @@ SPICA_EXPORTS Vector3d sampleCosineHemisphere(const Point2d& rands); SPICA_EXPORTS inline double cosineHemispherePdf(double cosTheta) { return cosTheta * INV_PI; } SPICA_EXPORTS void sampleUniformHemisphere(const Normal3d& normal, Vector3d* direction, const Point2d& rands); +SPICA_EXPORTS +void samplePoissonDisk(const Scene& scene, const Point3d& pCamera, + double minDist, std::vector* points); + SPICA_EXPORTS void samplePoissonDisk(const Scene& scene, const Point3d& pCamera, double minDist, std::vector* points); diff --git a/sources/integrator/hierarchical.cc b/sources/integrator/hierarchical.cc index 73e00ff..d12c1fa 100644 --- a/sources/integrator/hierarchical.cc +++ b/sources/integrator/hierarchical.cc @@ -24,183 +24,185 @@ namespace spica { -HierarchicalIntegrator::Octree::Octree() - : _root(NULL) - , _numCopies(NULL) - , _parent(NULL) { -} - -HierarchicalIntegrator::Octree::~Octree() { - release(); -} - -HierarchicalIntegrator::Octree::Octree(const Octree& octree) - : _root(NULL) - , _numCopies(NULL) - , _parent(NULL) { - this->operator=(octree); -} - -HierarchicalIntegrator::Octree& -HierarchicalIntegrator::Octree::operator=(const Octree& octree) { - release(); - _numCopies = octree._numCopies; - (*_numCopies) += 1; - _root = octree._root; - _parent = octree._parent; - return *this; -} +struct IrradiancePoint { + IrradiancePoint() + : pos{ 0.0, 0.0, 0.0 } + , area{ 0.0 } + , E{ 0.0 } { + } -void HierarchicalIntegrator::Octree::release() { - if (_numCopies != NULL) { - if ((*_numCopies) == 0) { - deleteNode(_root); - _root = nullptr; - delete _numCopies; - _numCopies = nullptr; - } - else { - (*_numCopies) -= 1; - } + IrradiancePoint(const Point3d& p, double A, const Spectrum& ir) + : pos{ p } + , area{ A } + , E{ ir } { } -} -void HierarchicalIntegrator::Octree::deleteNode( - HierarchicalIntegrator::OctreeNode* node) { - if (node != nullptr) { + Point3d pos; + double area; + Spectrum E; +}; + +struct OctreeNode { + OctreeNode() + : pt{} + , bbox{} + , isLeaf{ false } { for (int i = 0; i < 8; i++) { - deleteNode(node->children[i]); + children[i] = nullptr; } - delete node; - } -} + } -void HierarchicalIntegrator::Octree::construct( - HierarchicalIntegrator* parent, - std::vector& ipoints) { - release(); - - this->_parent = parent; - const int numHitpoints = static_cast(ipoints.size()); - + IrradiancePoint pt; Bounds3d bbox; - for (int i = 0; i < numHitpoints; i++) { - bbox.merge(ipoints[i].pos); + OctreeNode* children[8]; + bool isLeaf; +}; + +class HierarchicalIntegrator::Octree : public Uncopyable { +public: + // Public methods + Octree(double maxError) + : root_{ nullptr } + , maxError_{ maxError } { } - _root = constructRec(ipoints, bbox); -} - -HierarchicalIntegrator::OctreeNode* -HierarchicalIntegrator::Octree::constructRec( - std::vector& ipoints, - const Bounds3d& bbox) { - - // Zero or one child case - if (ipoints.empty()) { - return nullptr; - } else if (ipoints.size() == 1) { - OctreeNode* node = new OctreeNode(); - node->pt = ipoints[0]; - node->bbox = bbox; - node->isLeaf = true; - return node; + ~Octree() { + release(); } - // Divide children into eight groups - Point3d posMid = (bbox.posMin() + bbox.posMax()) * 0.5; - const int numPoints = static_cast(ipoints.size()); - std::vector > childPoints(8); - for (int i = 0; i < numPoints; i++) { - const Point3d& v = ipoints[i].pos; - int id = (v.x() < posMid.x() ? 0 : 4) + - (v.y() < posMid.y() ? 0 : 2) + - (v.z() < posMid.z() ? 0 : 1); - childPoints[id].push_back(ipoints[i]); + void release() { + deleteNode(root_); + root_ = nullptr; } - // Compute child nodes - OctreeNode* node = new OctreeNode(); - for (int i = 0; i < 8; i++) { - Bounds3d childBox; - for (int j = 0; j < childPoints[i].size(); j++) { - childBox.merge(childPoints[i][j].pos); + void deleteNode(OctreeNode* node) { + if (node != nullptr) { + for (int i = 0; i < 8; i++) { + deleteNode(node->children[i]); + node->children[i] = nullptr; + } + delete node; } - node->children[i] = constructRec(childPoints[i], childBox); } - node->bbox = bbox; - node->isLeaf = false; - - // Accumulate child nodes - node->pt.pos = Point3d(0.0, 0.0, 0.0); - node->pt.normal = Normal3d(0.0, 0.0, 0.0); - node->pt.area = 0.0; - - double sumWgt = 0.0; - int nChildren = 0; - for (int i = 0; i < 8; i++) { - if (node->children[i] != nullptr) { - const double weight = node->children[i]->pt.irad.luminance(); - node->pt.pos += weight * node->children[i]->pt.pos; - node->pt.normal += weight * node->children[i]->pt.normal; - node->pt.area += node->children[i]->pt.area; - node->pt.irad += node->children[i]->pt.irad; - sumWgt += weight; - nChildren += 1; + + void construct(HierarchicalIntegrator* parent, + const std::vector& ipoints) { + release(); + + Bounds3d bounds; + for (const auto& p : ipoints) { + bounds.merge(p.pos); } - } - if (sumWgt > 0.0) { - node->pt.pos /= sumWgt; - node->pt.normal /= sumWgt; + root_ = constructRec(ipoints, bounds); } - if (nChildren != 0) { - node->pt.irad /= nChildren; + Spectrum Mo(const SurfaceInteraction& po, + const std::unique_ptr& Rd) const { + return MoRec(root_, po, Rd); } - return node; -} +private: + // Private methods + OctreeNode* constructRec(const std::vector& ipoints, + const Bounds3d& bbox) { + + // Zero or one child case + if (ipoints.empty()) { + return nullptr; + } else if (ipoints.size() == 1) { + OctreeNode* node = new OctreeNode(); + node->pt = ipoints[0]; + node->bbox = bbox; + node->isLeaf = true; + return node; + } -Spectrum HierarchicalIntegrator::Octree::iradSubsurface( - const SurfaceInteraction& po, - const std::unique_ptr& Rd) const { - return iradSubsurfaceRec(_root, po, Rd); -} + // Divide children into eight groups + Point3d posMid = (bbox.posMin() + bbox.posMax()) * 0.5; + const int numPoints = static_cast(ipoints.size()); + std::vector> childPoints(8); + for (int i = 0; i < numPoints; i++) { + const Point3d& v = ipoints[i].pos; + int id = (v.x() < posMid.x() ? 0 : 4) + + (v.y() < posMid.y() ? 0 : 2) + + (v.z() < posMid.z() ? 0 : 1); + childPoints[id].push_back(ipoints[i]); + } -Spectrum HierarchicalIntegrator::Octree::iradSubsurfaceRec( - OctreeNode* node, - const SurfaceInteraction& po, - const std::unique_ptr& Rd) const { - if (node == nullptr) { - return Spectrum(0.0); - } + // Compute child nodes + OctreeNode* node = new OctreeNode(); + for (int i = 0; i < 8; i++) { + Bounds3d childBox; + for (int j = 0; j < childPoints[i].size(); j++) { + childBox.merge(childPoints[i][j].pos); + } + node->children[i] = constructRec(childPoints[i], childBox); + } + node->bbox = bbox; + node->isLeaf = false; + + // Accumulate child nodes + node->pt.pos = Point3d(0.0, 0.0, 0.0); + node->pt.area = 0.0; - double distSquared = (node->pt.pos - po.pos()).squaredNorm(); - double dw = node->pt.area / distSquared; - if (node->isLeaf || (dw < _parent->maxError_ && - !node->bbox.inside(po.pos()))) { - return (*Rd)(node->pt.pos, po.pos()) * node->pt.irad * node->pt.area; - } else { - Spectrum ret(0.0); + double sumWgt = 0.0; + int nChildren = 0; for (int i = 0; i < 8; i++) { if (node->children[i] != nullptr) { - ret += iradSubsurfaceRec(node->children[i], po, Rd); + const double weight = node->children[i]->pt.E.luminance(); + node->pt.pos += weight * node->children[i]->pt.pos; + node->pt.area += node->children[i]->pt.area; + node->pt.E += node->children[i]->pt.E; + sumWgt += weight; + nChildren += 1; } } - return ret; + + if (sumWgt > 0.0) { + node->pt.pos /= sumWgt; + } + node->pt.E /= nChildren; + + return node; } -} + + Spectrum MoRec(OctreeNode* node, + const SurfaceInteraction& po, + const std::unique_ptr& Rd) const { + if (node == nullptr) { + return Spectrum(0.0); + } + + double distSquared = (node->pt.pos - po.pos()).squaredNorm(); + double dw = node->pt.area / distSquared; + if (node->isLeaf || (dw < maxError_ && !node->bbox.inside(po.pos()))) { + return (*Rd)(node->pt.pos, po.pos()) * node->pt.E * node->pt.area; + } else { + Spectrum ret(0.0); + for (int i = 0; i < 8; i++) { + if (node->children[i] != nullptr) { + ret += MoRec(node->children[i], po, Rd); + } + } + return ret; + } + } + +private: + OctreeNode* root_; + double maxError_; + +}; // class Octree HierarchicalIntegrator::HierarchicalIntegrator( const std::shared_ptr& camera, const std::shared_ptr& sampler, double maxError) : SamplerIntegrator{ camera, sampler } - , octree_{} + , octree_{ std::make_unique(maxError) } , dA_{ 0.0 } - , radius_{} - , maxError_{ maxError } { + , radius_{} { } HierarchicalIntegrator::~HierarchicalIntegrator() { @@ -222,7 +224,7 @@ Spectrum HierarchicalIntegrator::Li(const Scene& scene, bool isIntersect = scene.intersect(ray, &isect); // Sample Le which contributes without any loss - if (bounces == 0 || specularBounce) { + if ((bounces == 0 && depth >= 0) || specularBounce) { if (isIntersect) { L += beta * isect.Le(-ray.dir()); } else { @@ -318,17 +320,34 @@ void HierarchicalIntegrator::buildOctree(const Scene& scene, const int threadID = getThreadID(); Spectrum E(0.0); for (int s = 0; s < nSamples; s++) { + // Indirect lighting Vector3d n(points_[i].normal()); Vector3d u, v; vect::coordinateSystem(n, &u, &v); - Vector3d dir = sampleCosineHemisphere(samplers[threadID]->get2D()); + Vector3d dir = sampleCosineHemisphere(samplers[threadID]->get2D()); double pdfDir = cosineHemispherePdf(vect::cosTheta(dir)); if (pdfDir == 0.0) continue; dir = u * dir.x() + v * dir.y() + n * dir.z(); - Ray ray(points_[i].pos(), dir); - E += Li(scene, params, ray, *samplers[threadID], arenas[threadID]) * + Ray ray = points_[i].spawnRay(dir); + E += Li(scene, params, ray, *samplers[threadID], arenas[threadID], -1) * vect::dot(n, dir) / pdfDir; + + // Direct lighting + for (const auto& l : scene.lights()) { + Vector3d wi; + double lightPdf; + VisibilityTester vis; + Spectrum Li = l->sampleLi(points_[i], samplers[threadID]->get2D(), &wi, &lightPdf, &vis); + + if (vect::dot(wi, points_[i].normal()) < 0.0) continue; + if (Li.isBlack() || lightPdf == 0.0) continue; + + Li *= vis.transmittance(scene, *samplers[threadID]); + if (vis.unoccluded(scene)) { + E += Li * vect::absDot(wi, points_[i].normal()) / lightPdf; + } + } } irads[i] = E / nSamples; }); @@ -336,19 +355,18 @@ void HierarchicalIntegrator::buildOctree(const Scene& scene, // Octree construction std::vector iradPoints(numPoints); for (int i = 0; i < numPoints; i++) { - iradPoints[i].pos = points_[i].pos(); - iradPoints[i].normal = points_[i].normal(); - iradPoints[i].area = dA_; - iradPoints[i].irad = irads[i]; + iradPoints[i].pos = points_[i].pos(); + iradPoints[i].area = dA_; + iradPoints[i].E = irads[i]; } - octree_.construct(this, iradPoints); + octree_->construct(this, iradPoints); std::cout << "Octree constructed !!" << std::endl; } Spectrum HierarchicalIntegrator::irradiance(const SurfaceInteraction& po) const { Assertion(po.bssrdf(), "BSSRDF not found!!"); auto Rd = static_cast(po.bssrdf())->Rd(); - const Spectrum Mo = octree_.iradSubsurface(po, Rd); + const Spectrum Mo = octree_->Mo(po, Rd); return (INV_PI * (1.0 - Rd->Fdr())) * Mo; } diff --git a/sources/integrator/hierarchical.h b/sources/integrator/hierarchical.h index 5d2ca56..5502798 100644 --- a/sources/integrator/hierarchical.h +++ b/sources/integrator/hierarchical.h @@ -15,47 +15,6 @@ namespace spica { -struct IrradiancePoint { - Point3d pos; - Normal3d normal; - double area; - Spectrum irad; - - IrradiancePoint() - : pos() - , normal() - , area(0.0) - , irad() - { - } - - IrradiancePoint(const Point3d& pos_, const Normal3d& normal_, - const double area_, const Spectrum& irad_) - : pos(pos_) - , normal(normal_) - , area(area_) - , irad(irad_) - { - } - - IrradiancePoint(const IrradiancePoint& p) - : pos() - , normal() - , area(0.0) - , irad() - { - this->operator=(p); - } - - IrradiancePoint& operator=(const IrradiancePoint& p) { - this->pos = p.pos; - this->normal = p.normal; - this->area = p.area; - this->irad = p.irad; - return *this; - } -}; - /** Irradiance integrator for subsurface scattering objects * @ingroup renderer_module */ @@ -83,55 +42,6 @@ class SPICA_EXPORTS HierarchicalIntegrator : public SamplerIntegrator { int depth = 0) const override; private: - // Private internal classes - struct OctreeNode { - IrradiancePoint pt; - Bounds3d bbox; - OctreeNode* children[8]; - bool isLeaf; - - OctreeNode() - : pt{} - , bbox{} - , isLeaf{ false } { - for (int i = 0; i < 8; i++) { - children[i] = nullptr; - } - } - }; - - class Octree { - private: - OctreeNode* _root; - int* _numCopies; - HierarchicalIntegrator* _parent; - - public: - Octree(); - ~Octree(); - - Octree(const Octree& octree); - Octree& operator=(const Octree& octree); - - void construct(HierarchicalIntegrator* parent, - std::vector& ipoints); - - Spectrum iradSubsurface( - const SurfaceInteraction& intr, - const std::unique_ptr& Rd) const; - - private: - void release(); - void deleteNode(OctreeNode* node); - OctreeNode* constructRec(std::vector& pointers, - const Bounds3d& bbox); - - Spectrum iradSubsurfaceRec( - OctreeNode* node, - const SurfaceInteraction& po, - const std::unique_ptr& Rd) const; - }; - // Private methods Spectrum irradiance(const SurfaceInteraction& po) const; @@ -140,11 +50,12 @@ class SPICA_EXPORTS HierarchicalIntegrator : public SamplerIntegrator { Sampler& sampler); // Private fields - std::vector points_; - Octree octree_; + std::vector points_; + + class Octree; + std::unique_ptr octree_; double dA_; double radius_; - double maxError_; }; // class HierarchicalIntegrator