lightcuts.cpp

#include "lightcuts.hpp"

template <class PointLightType, class PointLightNodeType>
class PointLightTree
{
public:
	Spectrum samplePosition(Vec3 * wi, bool * isVisible, const Intersection & isect, const Scene & scene, const Sampler & sampler) const { assert(false && "unimplemented"); return Spectrum(0.0); };

	Spectrum sampleDirection(Vec3 * wi, Vec3 * position, const Sampler & sampler) const { assert(false && "unimplemented"); return Spectrum(0.0); }

	Spectrum eval(const Intersection & isect, const Vec3 & wo, const Scene & scene) const
	{
		stat.cutSize = 0;
		return this->evalBranchExhaust(mRoot, isect, wo, scene);
	}

	Spectrum evalBranchExhaust(const shared_ptr<PointLightNodeType> & currentNode, const Intersection & isect, const Vec3 & wo, const Scene & scene) const
	{
		assert(currentNode != nullptr);
		stat.cutSize++;
		if (currentNode->isLeaf())
		{
			stat.geometryTermEval++;
			return currentNode->eval(isect, wo, scene);
		}
		else
		{
			return this->evalBranchExhaust(currentNode->mChild[0], isect, wo, scene) + this->evalBranchExhaust(currentNode->mChild[1], isect, wo, scene);
		}
	}

	Spectrum evalBranch(const shared_ptr<PointLightNodeType> & currentNode, const Intersection & isect, const Vec3 & wo, const Scene & scene) const
	{
		assert(currentNode != nullptr);

		struct HeapNode
		{
			const PointLightNodeType * node;
			Float error;
			Float gvTerm;
			Spectrum contrib;

			HeapNode() {}
		};

		auto CompareLessHeapNode = [] (const HeapNode & a, const HeapNode & b) { return a.error < b.error; };

		std::priority_queue<HeapNode, std::vector<HeapNode>, decltype(CompareLessHeapNode)> heap(CompareLessHeapNode);

		uint32_t cutSize = 0;

		HeapNode rootHeapNode;
		rootHeapNode.node = currentNode->get();
		rootHeapNode.gvTerm = currentNode->geometryVisibility(isect, scene);
		rootHeapNode.contrib = currentNode->eval(rootHeapNode.gvTerm, isect, wo);
		rootHeapNode.error = currentNode->maxError(isect, wo);

		while (heap.size() != 0 && cutSize < 1000)
		{
			const HeapNode top = heap.top();
		}
	}

	void build(const std::vector<shared_ptr<PointLightType>> & pointLights, const Sampler & sampler)
	{
		const size_t numLights = pointLights.size();

		if (numLights == 0) { return; }

		// compute constant c that control relative scaling
		Float c2 = 0.0f;
		{
			Aabb bound;
			for (const shared_ptr<PointLightType> & light : pointLights)
			{
				bound = Aabb::Union(light->mPosition, bound);
			}
			c2 = Aabb::DiagonalLength2(bound);
		}

		// populate active nodes
		std::vector<shared_ptr<PointLightNodeType>> activeNodes(numLights);
		for (size_t i = 0; i < numLights; i++)
		{
			const shared_ptr<PointLightType> & rp = pointLights[i];
			activeNodes[i] = make_shared<PointLightNodeType>(*rp);
		}

		/// TODO:: wait for @kobayashi to unify oriented light and omni light kdtree and 
		while (activeNodes.size() != 1)
		{
			const size_t numActiveNodes = activeNodes.size();
			std::pair<size_t, size_t> minPair;
			Float minDistance = std::numeric_limits<Float>::infinity();
			for (size_t ia = 0;ia < numActiveNodes - 1;ia++)
			{
				const shared_ptr<PointLightNodeType> & rna = activeNodes[ia];
				for (size_t ib = ia + 1;ib < numActiveNodes;ib++)
				{
					const shared_ptr<PointLightNodeType> & rnb = activeNodes[ib];
					Float distance = rna->cost(*rnb, c2);
					if (distance < minDistance)
					{
						minDistance = distance;
						minPair.first = ia;
						minPair.second = ib;
					}
				}
			}

			shared_ptr<PointLightNodeType> & minRnA = activeNodes[minPair.first];
			shared_ptr<PointLightNodeType> & minRnB = activeNodes[minPair.second];
			shared_ptr<PointLightNodeType> newRn = make_shared<PointLightNodeType>(minRnA, minRnB, sampler);

			std::swap(activeNodes[minPair.second], activeNodes.back());
			activeNodes.pop_back();
			std::swap(activeNodes[minPair.first], activeNodes.back());
			activeNodes.pop_back();

			activeNodes.push_back(newRn);
		}
		this->mRoot = activeNodes[0];
	}

	void buildFast(const std::vector<shared_ptr<PointLightType>> & pointLights)
	{
		Rng rng(81521);
		const size_t numPointLights = pointLights.size();

		if (numPointLights == 0) { return; }

		// compute constant c that control relative scaling
		Float diagonal2 = (Float)0.0;
		{
			Aabb bound;
			for (const shared_ptr<PointLightType> & lightPtr : pointLights)
			{
				bound = Aabb::Union(lightPtr->mPosition, bound);
			}
			diagonal2 = Aabb::DiagonalLength2(bound);
		}


		// populate active nodes
		std::vector<const PointLightNodeType *> activeNodes(numPointLights);
		{
			for (size_t i = 0; i < numPointLights; i++)
			{
				const PointLightType & p = *pointLights[i];
				activeNodes[i] = new PointLightNodeType(p);
			}
		}

		// declare kdtree stuffs
		using KDT = LcKdTree<PointLightNodeType>;
		struct QueueElem
		{
			struct NodePair
			{
				shared_ptr<typename KDT::LcKdTreeNode> kdtnode;
				const PointLightNodeType *lightNode;

				bool obsoleted() const
				{
					if (kdtnode->deleted()) { return true; }
					if (kdtnode->lightNode != lightNode) { return true; }
					return false;
				}
			};

			QueueElem(NodePair first, NodePair second, float distance):
				first(first), second(second), distance(distance)
			{
			}

			QueueElem(const QueueElem&) = default;

			NodePair first, second;
			float distance;
			bool operator >(const QueueElem &rhs) const { return distance > rhs.distance; }
		};
		using Queue = std::priority_queue<QueueElem, std::vector<QueueElem>, std::greater<QueueElem>>;

		Queue pairs;
		KDT kdt(activeNodes, diagonal2);

		auto pushNearestToQueue = [&] (shared_ptr<KDT::LcKdTreeNode> p)
		{
			auto nearest = kdt.findNearest(p);
			if (nearest.first)
			{
				pairs.emplace(
					QueueElem::NodePair{p, p->lightNode},
					QueueElem::NodePair{nearest.first, nearest.first->lightNode},
					nearest.second);
			}
		};

		for (auto p = kdt.begin(); p; p = kdt.next(p))
		{
			pushNearestToQueue(p);
		}

		while (!pairs.empty())
		{
			auto pair = pairs.top(); pairs.pop();

			if (pair.first.obsoleted() && pair.second.obsoleted())
			{
				continue;
			}
			else if (pair.first.obsoleted())
			{
				pushNearestToQueue(pair.second.kdtnode);
			}
			else if (pair.second.obsoleted())
			{
				pushNearestToQueue(pair.first.kdtnode);
			}
			else
			{
				const PointLightNodeType * minPlnA = pair.first.kdtnode->lightNode;
				const PointLightNodeType * minPlnB = pair.second.kdtnode->lightNode;
				const PointLightNodeType * newPln = PointLightNodeType::MergeNode(minPlnA, minPlnB, rng);
				auto remain_node = pair.first.kdtnode;
				auto delete_node = pair.second.kdtnode;
				if (remain_node->depth > delete_node->depth)
				{
					std::swap(remain_node, delete_node);
				}
				remain_node->lightNode = newPln;
				delete_node->lightNode = nullptr;
				kdt.cleanup(delete_node);
				if (!pairs.empty())
				{
					if (!kdt.findNearest(remain_node).first)
					{
						kdt.findNearest(remain_node);
					}
				}
				pushNearestToQueue(remain_node);
			}
		}
		auto begin = kdt.begin();
		if (!begin)
		{
			throw std::runtime_error("unreachable! empty tree");
		}
		assert(begin);
		assert(begin->lightNode);
		this->mRoot = begin->lightNode;
	}

	struct
	{
		mutable std::atomic<int> cutSize = 0;
		mutable std::atomic<int> geometryTermEval = 0;
	} stat;

protected:
	shared_ptr<PointLightNodeType> mRoot;
};

lightcuts.hpp

#pragma once

#include <atomic>
#include <queue>
#include <functional>

#include "common\reflectcuts.h"
#include "common\light.h"
#include "common\rng.h"

#include "lightaccels\util\lckdtree.h"

template <typename DerivedType>
class PointLightNode
{
public:
	// must be symmetric a.cost(b) = b.cost(a)
	virtual Float cost(const DerivedType & b, const Float c2) const = 0;

	// for evaluating lightcuts
	virtual Float maxError(const Intersection & isect, const Vec3 & wo) const = 0;

	virtual Float geometryVisibility(const Intersection & isect, const Scene & scene) const = 0;

	// return mRepresent->eval
	virtual Spectrum eval(const Intersection & isect, const Vec3 & wo, const Scene & scene) const = 0;

	// same as eval but don't need to recompute G term
	virtual Spectrum eval(const Float computedG, const Intersection & isect, const Vec3 & wo) const = 0;

	// is leaf
	virtual bool isLeaf() const = 0;
};