mingzailao
/
mitsuba
mirror of https://github.com/Quartenia/mitsuba.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

implemented a photon map performance testcase

metadata
Wenzel Jakob 2011-08-28 03:09:35 -04:00
parent 2ba07a76fd
commit 749d30ab7f
7 changed files with 792 additions and 165 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

560
include/mitsuba/core/kdtree.h
View File

@ -20,58 +20,80 @@
#define __KDTREE_H
#include <mitsuba/core/aabb.h>
#include <boost/foreach.hpp>
#include <mitsuba/core/timer.h>
MTS_NAMESPACE_BEGIN
/**
* \brief Basic point node for use with \ref TKDTree.
* \brief Simple kd-tree node for use with \ref PointKDTree.
*
* \tparam PointType Underlying point data type (e.g. \ref TPoint3<float>)
* \tparam DataRecord Custom payload to be attached to each node
* This class is an example of how one might write a space-efficient kd-tree
* node that is compatible with the \ref PointKDTree class. The implementation
* supports associating a custom data record with each node and works up to
* 16 dimensions.
*
* \tparam _PointType Underlying point data type (e.g. \ref TPoint3<float>)
* \tparam _DataRecord Custom storage that should be associated with each
* tree node
*
* \ingroup libcore
* \sa PointKDTree
* \sa LeftBalancedKDNode
*/
template <typename PointType, typename DataRecord> struct BasicKDNode {
typedef PointType point_type;
template <typename _PointType, typename _DataRecord> struct SimpleKDNode {
typedef _PointType PointType;
typedef _DataRecord DataRecord;
typedef uint32_t IndexType;
typedef typename PointType::value_type Scalar;
enum {
ELeafFlag = 0x10,
EAxisMask = 0x0F
};
static const bool leftBalancedLayout = false;
PointType position;
uint32_t right;
uint16_t flags;
uint16_t axis;
DataRecord value;
IndexType right;
DataRecord data;
uint8_t flags;
/// Initialize a KD-tree node
inline BasicKDNode() : position((Float) 0),
right(0), flags(0), axis(0), value() { }
inline SimpleKDNode() : position((Scalar) 0),
right(0), data(), flags(0) { }
/// Initialize a KD-tree node with the given data record
inline BasicKDNode(const DataRecord &value) : position((Float) 0),
right(0), flags(0), axis(0), value(value) { }
inline SimpleKDNode(const DataRecord &data) : position((Scalar) 0),
right(0), data(data), flags(0) { }
/// Given the current node's index, return the index of the right child
inline uint32_t getRightIndex(uint32_t curIndex) { return right; }
/// Given the current node's index, return the index of the right child (const version)
inline const uint32_t getRightIndex(uint32_t curIndex) const { return right; }
inline IndexType getRightIndex(IndexType self) const { return right; }
/// Given the current node's index, set the right child index
inline void setRightIndex(uint32_t curIndex, uint32_t value) { right = value; }
inline void setRightIndex(IndexType self, IndexType value) { right = value; }
/// Given the current node's index, return the index of the left child
inline uint32_t getLeftIndex(uint32_t curIndex) { return curIndex + 1; }
/// Given the current node's index, return the index of the left child (const version)
inline const uint32_t getLeftIndex(uint32_t curIndex) const { return curIndex + 1; }
inline IndexType getLeftIndex(IndexType self) const { return self + 1; }
/// Given the current node's index, set the left child index
inline void setLeftIndex(uint32_t curIndex, uint32_t value) {
if (value != curIndex+1) SLog(EError, "Not supported!"); }
inline void setLeftIndex(IndexType self, IndexType value) {
#if defined(MTS_DEBUG)
if (value != self+1)
SLog(EError, "SimpleKDNode::setLeftIndex(): Internal error!");
#endif
}
/// Check whether this is a leaf node
inline bool isLeaf() const { return flags & 1; }
inline bool isLeaf() const { return flags & (uint8_t) ELeafFlag; }
/// Specify whether this is a leaf node
inline void setLeaf(bool value) { if (value) flags |= 1; else flags &= ~1; }
inline void setLeaf(bool value) {
if (value)
flags |= (uint8_t) ELeafFlag;
else
flags &= (uint8_t) ~ELeafFlag;
}
/// Return the split axis associated with this node
inline uint16_t getAxis() const { return axis; }
/// Set the split axis associated with this node
inline void setAxis(uint16_t value) { axis = value; }
inline uint16_t getAxis() const { return flags & (uint8_t) EAxisMask; }
/// Set the split flags associated with this node
inline void setAxis(uint8_t axis) { flags = (flags & (uint8_t) ~EAxisMask) | axis; }
/// Return the position associated with this node
inline const PointType &getPosition() const { return position; }
@ -79,32 +101,126 @@ template <typename PointType, typename DataRecord> struct BasicKDNode {
inline void setPosition(const PointType &value) { position = value; }
/// Return the data record associated with this node
inline DataRecord &getValue() { return value; }
inline DataRecord &getData() { return data; }
/// Return the data record associated with this node (const version)
inline const DataRecord &getValue() const { return value; }
inline const DataRecord &getData() const { return data; }
/// Set the data record associated with this node
inline void setValue(const DataRecord &val) { value = val; }
inline void setData(const DataRecord &val) { data = val; }
};
/**
* \brief Left-balanced kd-tree node for use with \ref PointKDTree.
*
* This class provides a basic kd-tree node layout that can be used with
* the \ref PointKDTree class and the \ref PointKDTree::ELeftBalanced
* tree construction heuristic. The advantage of a left-balanced tree is
* that there is no need to store node pointers within nodes.
*
* \tparam _PointType Underlying point data type (e.g. \ref TPoint3<float>)
* \tparam _DataRecord Custom storage that should be associated with each
* tree node
*
* \ingroup libcore
* \sa PointKDTree
* \sa SimpleKDNode
*/
template <typename _PointType, typename _DataRecord> struct LeftBalancedKDNode {
typedef _PointType PointType;
typedef _DataRecord DataRecord;
typedef uint32_t IndexType;
typedef typename PointType::value_type Scalar;
enum {
ELeafFlag = 0x10,
EAxisMask = 0x0F
};
static const bool leftBalancedLayout = true;
PointType position;
DataRecord data;
uint8_t flags;
/// Initialize a KD-tree node
inline LeftBalancedKDNode() : position((Scalar) 0), data(), flags(0) { }
/// Initialize a KD-tree node with the given data record
inline LeftBalancedKDNode(const DataRecord &data) : position((Scalar) 0),
data(data), flags(0) { }
/// Given the current node's index, return the index of the left child
inline IndexType getLeftIndex(IndexType self) const { return 2 * self + 1; }
/// Given the current node's index, set the left child index
inline void setLeftIndex(IndexType self, IndexType value) {
#if defined(MTS_DEBUG)
if (value != 2*self + 1)
SLog(EError, "LeftBalancedKDNode::setLeftIndex(): Internal error!");
#endif
}
/// Given the current node's index, return the index of the right child
inline IndexType getRightIndex(IndexType self) const { return 2 * self + 2; }
/// Given the current node's index, set the right child index
inline void setRightIndex(IndexType self, IndexType value) {
#if defined(MTS_DEBUG)
if (value != 0 && value != 2*self + 2)
SLog(EError, "LeftBalancedKDNode::setRightIndex(): Internal error!");
#endif
}
/// Check whether this is a leaf node
inline bool isLeaf() const { return flags & (uint8_t) ELeafFlag; }
/// Specify whether this is a leaf node
inline void setLeaf(bool value) {
if (value)
flags |= (uint8_t) ELeafFlag;
else
flags &= (uint8_t) ~ELeafFlag;
}
/// Return the split axis associated with this node
inline uint16_t getAxis() const { return flags & (uint8_t) EAxisMask; }
/// Set the split flags associated with this node
inline void setAxis(uint8_t axis) { flags = (flags & (uint8_t) ~EAxisMask) | axis; }
/// Return the position associated with this node
inline const PointType &getPosition() const { return position; }
/// Set the position associated with this node
inline void setPosition(const PointType &value) { position = value; }
/// Return the data record associated with this node
inline DataRecord &getData() { return data; }
/// Return the data record associated with this node (const version)
inline const DataRecord &getData() const { return data; }
/// Set the data record associated with this node
inline void setData(const DataRecord &val) { data = val; }
};
/**
* \brief Generic multi-dimensional kd-tree data structure for point data
*
* Organizes a list of point data in a hierarchical manner. For data
* with spatial extents, \ref GenericKDTree and \ref ShapeKDTree will be
* more appropriate.
* This class organizes point data in a hierarchical manner so various
* types of queries can be performed efficiently. It supports several
* heuristics for building ``good'' trees, and it is oblivious to the
* data layout of the nodes themselves.
*
* \tparam KDNode Underlying node data structure. See \ref BasicKDNode as
* Note that this class is meant for point data only --- for things that
* have some kind of spatial extent, the classes \ref GenericKDTree and
* \ref ShapeKDTree will be more appropriate.
*
* \tparam _NodeType Underlying node data structure. See \ref SimpleKDNode as
* an example for the required public interface
*
* \ingroup libcore
* \see SimpleKDNode
*/
template <typename KDNode> class TKDTree {
template <typename _NodeType> class PointKDTree {
public:
typedef KDNode node_type;
typedef typename KDNode::point_type point_type;
typedef typename point_type::value_type value_type;
typedef typename point_type::vector_type vector_type;
typedef TAABB<point_type> aabb_type;
typedef _NodeType NodeType;
typedef typename NodeType::PointType PointType;
typedef typename NodeType::IndexType IndexType;
typedef typename PointType::value_type Scalar;
typedef typename PointType::vector_type VectorType;
typedef TAABB<PointType> AABBType;
/// Supported tree construction heuristics
enum EHeuristic {
@ -122,9 +238,11 @@ public:
/**
* \brief Choose the split plane by optimizing a cost heuristic
* based on the ratio of voxel volumes. Note that the implementation
* here is not particularly optimized, and it furthermore runs in
* time O(n (log n)^2) instead of O(n log n)
* based on the ratio of voxel volumes.
*
* Note that Mitsuba's implementation of this heuristic is not
* particularly optimized --- the tree construction construction
* runs time O(n (log n)^2) instead of O(n log n).
*/
EVoxelVolume
};
@ -132,9 +250,11 @@ public:
/// Result data type for k-nn queries
struct SearchResult {
Float distSquared;
uint32_t index;
IndexType index;
inline SearchResult(Float distSquared, uint32_t index)
inline SearchResult() {}
inline SearchResult(Float distSquared, IndexType index)
: distSquared(distSquared), index(index) { }
std::string toString() const {
@ -164,16 +284,28 @@ public:
* \brief Create an empty KD-tree that can hold the specified
* number of points
*/
inline TKDTree(size_t nodes, EHeuristic heuristic = ESlidingMidpoint)
inline PointKDTree(size_t nodes = 0, EHeuristic heuristic = ESlidingMidpoint)
: m_nodes(nodes), m_heuristic(heuristic), m_depth(0) { }
/// Clear the kd-tree node array
inline void clear() { m_nodes.clear(); }
/// Resize the kd-tree node array
inline void resize(size_t size) { m_nodes.resize(size); }
/// Reserve a certain amount of memory for the kd-tree node array
inline void reserve(size_t size) { m_nodes.reserve(size); }
/// Append a kd-tree node to the node array
inline void push_back(const NodeType &node) { m_nodes.push_back(node); }
/// Return one of the KD-tree nodes by index
inline KDNode &operator[](size_t idx) { return m_nodes[idx]; }
inline NodeType &operator[](size_t idx) { return m_nodes[idx]; }
/// Return one of the KD-tree nodes by index (const version)
inline const KDNode &operator[](size_t idx) const { return m_nodes[idx]; }
inline const NodeType &operator[](size_t idx) const { return m_nodes[idx]; }
/// Return the AABB of the underlying point data
inline const aabb_type &getAABB() const { return m_aabb; }
inline const AABBType &getAABB() const { return m_aabb; }
/// Return the depth of the constructed KD-tree
inline size_t getDepth() const { return m_depth; }
/// Return the size of the kd-tree
@ -181,14 +313,45 @@ public:
/// Construct the KD-tree hierarchy
void build() {
ref<Timer> timer = new Timer();
SLog(EInfo, "Building a %i-dimensional kd-tree over " SIZE_T_FMT " data points",
PointType::dim, m_nodes.size());
m_aabb.reset();
BOOST_FOREACH(KDNode &node, m_nodes) {
m_aabb.expandBy(node.getPosition());
}
for (size_t i=0; i<m_nodes.size(); ++i)
m_aabb.expandBy(m_nodes[i].getPosition());
int aabbTime = timer->getMilliseconds();
timer->reset();
/* Instead of shuffling around the node data itself, only modify
an indirection table initially. Once the tree construction
is done, this table will contain a indirection that can then
be applied to the data in one pass */
std::vector<IndexType> indirection(m_nodes.size());
for (size_t i=0; i<m_nodes.size(); ++i)
indirection[i] = i;
m_depth = 0;
build(1, m_nodes.begin(), m_nodes.end());
if (NodeType::leftBalancedLayout) {
std::vector<IndexType> permutation(m_nodes.size());
buildLB(0, 1, indirection.begin(), indirection.begin(),
indirection.end(), permutation);
permute_inplace(&m_nodes[0], permutation);
} else {
build(1, indirection.begin(), indirection.begin(), indirection.end());
permute_inplace(&m_nodes[0], indirection);
}
int constructionTime = timer->getMilliseconds();
timer->reset();
int permutationTime = timer->getMilliseconds();
SLog(EInfo, "Done after %i ms (breakdown: aabb: %i ms, build: %i ms, permute: %i ms). ",
aabbTime + constructionTime + permutationTime, aabbTime, constructionTime, permutationTime);
}
/**
@ -197,23 +360,28 @@ public:
* \param p Search position
* \param k Maximum number of search results
* \param results Index list of search results
* \param searchRadius Maximum search radius (this can be used to
* restrict the knn query to a subset of the data)
* \param sqrSearchRadius
* Specifies the squared maximum search radius. This parameter can be used
* to restrict the k-nn query to a subset of the data -- it that is not
* desired, simply set it to positive infinity. After the query
* finishes, the parameter value will correspond to the (potentially lower)
* maximum query radius that was necessary to ensure that the number of
* results did not exceed \c k.
* \return The number of used traversal steps
*/
size_t nnSearch(const point_type &p, size_t k, std::vector<SearchResult> &results,
Float searchRadius = std::numeric_limits<Float>::infinity()) const {
uint32_t *stack = (uint32_t *) alloca((m_depth+1) * sizeof(uint32_t));
size_t nnSearch(const PointType &p, size_t k, std::vector<SearchResult> &results,
Float &sqrSearchRadius) const {
// IndexType *stack = (IndexType *) alloca((m_depth+1) * sizeof(IndexType));
IndexType stack[30];
size_t index = 0, stackPos = 1, traversalSteps = 0;
bool isHeap = false;
Float distSquared = searchRadius*searchRadius;
stack[0] = 0;
results.clear();
results.reserve(k+1);
while (stackPos > 0) {
const KDNode &node = m_nodes[index];
const NodeType &node = m_nodes[index];
++traversalSteps;
int nextIndex;
@ -222,26 +390,31 @@ public:
Float distToPlane = p[node.getAxis()]
- node.getPosition()[node.getAxis()];
uint32_t first, second;
bool searchBoth = distToPlane*distToPlane <= distSquared;
IndexType first, second;
bool searchBoth = distToPlane*distToPlane <= sqrSearchRadius;
if (distToPlane > 0) {
first = node.getRightIndex(index);
if (NodeType::leftBalancedLayout && first >= m_nodes.size())
first = 0;
second = searchBoth ? node.getLeftIndex(index) : 0;
} else {
first = node.getLeftIndex(index);
second = searchBoth ? node.getRightIndex(index) : 0;
if (NodeType::leftBalancedLayout && second >= m_nodes.size())
second = 0;
}
if (first != 0 && second != 0) {
if (first && second) {
nextIndex = first;
stack[stackPos++] = second;
} else if (first != 0) {
nextIndex = first;
} else if (second != 0) {
nextIndex = second;
} else {
nextIndex = stack[--stackPos];
nextIndex = first | second;
if (!nextIndex)
nextIndex = stack[--stackPos];
}
} else {
nextIndex = stack[--stackPos];
@ -250,7 +423,7 @@ public:
/* Check if the current point is within the query's search radius */
const Float pointDistSquared = (node.getPosition() - p).lengthSquared();
if (pointDistSquared < distSquared) {
if (pointDistSquared < sqrSearchRadius) {
/* Switch to a max-heap when the available search
result space is exhausted */
if (results.size() < k) {
@ -272,7 +445,7 @@ public:
results.pop_back();
/* Reduce the search radius accordingly */
distSquared = results[0].distSquared;
sqrSearchRadius = results[0].distSquared;
}
}
index = nextIndex;
@ -280,27 +453,42 @@ public:
return traversalSteps;
}
/**
* \brief Run a k-nearest-neighbor search query without any
* search radius threshold
*
* \param p Search position
* \param k Maximum number of search results
* \param results Index list of search results
* \return The number of used traversal steps
*/
inline size_t nnSearch(const PointType &p, size_t k, std::vector<SearchResult> &results) {
Float searchRadiusSqr = std::numeric_limits<Float>::infinity();
return nnSearch(p, k, results, searchRadiusSqr);
}
/**
* \brief Execute a search query and run the specified functor on them,
* while potentially modifying nodes within the search radius
*
* The functor must have an operator() implementation, which accepts
* a \a KDNode as its argument.
* a \a NodeType as its argument.
*
* \param p Search position
* \param functor Functor to be called on each search result
* \param searchRadius Search radius
* \return The number of used traversal steps
*/
template <typename Functor> size_t executeModifier(const point_type &p,
template <typename Functor> size_t executeModifier(const PointType &p,
Float searchRadius, Functor &functor) {
uint32_t *stack = (uint32_t *) alloca((m_depth+1) * sizeof(uint32_t));
IndexType *stack = (IndexType *) alloca((m_depth+1) * sizeof(IndexType));
size_t index = 0, stackPos = 1, traversalSteps = 0;
Float distSquared = searchRadius*searchRadius;
stack[0] = 0;
while (stackPos > 0) {
KDNode &node = m_nodes[index];
NodeType &node = m_nodes[index];
++traversalSteps;
int nextIndex;
@ -309,7 +497,7 @@ public:
Float distToPlane = p[node.getAxis()]
- node.getPosition()[node.getAxis()];
uint32_t first, second;
IndexType first, second;
bool searchBoth = distToPlane*distToPlane <= distSquared;
if (distToPlane > 0) {
@ -349,22 +537,22 @@ public:
* \brief Execute a search query and run the specified functor on them
*
* The functor must have an operator() implementation, which accepts
* a constant reference to a \a KDNode as its argument.
* a constant reference to a \a NodeType as its argument.
*
* \param p Search position
* \param functor Functor to be called on each search result
* \param searchRadius Search radius
* \return The number of used traversal steps
*/
template <typename Functor> size_t executeQuery(const point_type &p,
template <typename Functor> size_t executeQuery(const PointType &p,
Float searchRadius, Functor &functor) const {
uint32_t *stack = (uint32_t *) alloca((m_depth+1) * sizeof(uint32_t));
IndexType *stack = (IndexType *) alloca((m_depth+1) * sizeof(IndexType));
size_t index = 0, stackPos = 1, traversalSteps = 0;
Float distSquared = searchRadius*searchRadius;
stack[0] = 0;
while (stackPos > 0) {
const KDNode &node = m_nodes[index];
const NodeType &node = m_nodes[index];
++traversalSteps;
int nextIndex;
@ -373,7 +561,7 @@ public:
Float distToPlane = p[node.getAxis()]
- node.getPosition()[node.getAxis()];
uint32_t first, second;
IndexType first, second;
bool searchBoth = distToPlane*distToPlane <= distSquared;
if (distToPlane > 0) {
@ -397,7 +585,7 @@ public:
} else {
nextIndex = stack[--stackPos];
}
/* Check if the current point is within the query's search radius */
const Float pointDistSquared = (node.getPosition() - p).lengthSquared();
@ -418,8 +606,8 @@ public:
* \param searchRadius Search radius
* \return The number of used traversal steps
*/
size_t search(const point_type &p, Float searchRadius, std::vector<uint32_t> &results) const {
uint32_t *stack = (uint32_t *) alloca((m_depth+1) * sizeof(uint32_t));
size_t search(const PointType &p, Float searchRadius, std::vector<IndexType> &results) const {
IndexType *stack = (IndexType *) alloca((m_depth+1) * sizeof(IndexType));
size_t index = 0, stackPos = 1, traversalSteps = 0;
Float distSquared = searchRadius*searchRadius;
stack[0] = 0;
@ -427,7 +615,7 @@ public:
results.clear();
while (stackPos > 0) {
const KDNode &node = m_nodes[index];
const NodeType &node = m_nodes[index];
++traversalSteps;
int nextIndex;
@ -436,7 +624,7 @@ public:
Float distToPlane = p[node.getAxis()]
- node.getPosition()[node.getAxis()];
uint32_t first, second;
IndexType first, second;
bool searchBoth = distToPlane*distToPlane <= distSquared;
if (distToPlane > 0) {
@ -473,76 +661,142 @@ public:
}
protected:
struct CoordinateOrdering : public std::binary_function<KDNode, KDNode, bool> {
struct CoordinateOrdering : public std::binary_function<IndexType, IndexType, bool> {
public:
inline CoordinateOrdering(int axis) : m_axis(axis) { }
inline bool operator()(const KDNode &n1, const KDNode &n2) const {
return n1.getPosition()[m_axis] < n2.getPosition()[m_axis];
inline CoordinateOrdering(std::vector<NodeType> &nodes, int axis)
: m_nodes(nodes), m_axis(axis) { }
inline bool operator()(const IndexType &i1, const IndexType &i2) const {
return m_nodes[i1].getPosition()[m_axis] < m_nodes[i2].getPosition()[m_axis];
}
private:
std::vector<NodeType> &m_nodes;
int m_axis;
};
struct LessThanOrEqual : public std::unary_function<KDNode, bool> {
struct LessThanOrEqual : public std::unary_function<IndexType, bool> {
public:
inline LessThanOrEqual(int axis, value_type value) : m_axis(axis), m_value(value) { }
inline bool operator()(const KDNode &n1) const {
return n1.getPosition()[m_axis] <= m_value;
inline LessThanOrEqual(std::vector<NodeType> &nodes, int axis, Scalar value)
: m_nodes(nodes), m_axis(axis), m_value(value) { }
inline bool operator()(const IndexType &i) const {
return m_nodes[i].getPosition()[m_axis] <= m_value;
}
private:
std::vector<NodeType> &m_nodes;
int m_axis;
value_type m_value;
Scalar m_value;
};
/**
* Given a number of entries, this method calculates the number of nodes
* nodes on the left subtree of a left-balanced tree. There are two main
* cases here:
*
* 1) It is possible to completely fill the left subtree
* 2) It doesn't work - the last level contains too few nodes, e.g :
* O
* / \
* O O
* /
* O
*
* The function assumes that "count" > 1.
*/
inline IndexType leftSubtreeSize(IndexType count) const {
/* Layer 0 contains one node */
IndexType p = 1;
/* Traverse downwards until the first incompletely
filled tree level is encountered */
while (2*p <= count)
p *= 2;
/* Calculate the number of filled slots in the last level */
IndexType remaining = count - p + 1;
if (2*remaining < p) {
/* Case 2: The last level contains too few nodes. Remove
overestimate from the left subtree node count and add
the remaining nodes */
p = (p >> 1) + remaining;
}
return p - 1;
}
void buildLB(IndexType idx, size_t depth,
typename std::vector<IndexType>::iterator base,
typename std::vector<IndexType>::iterator rangeStart,
typename std::vector<IndexType>::iterator rangeEnd,
typename std::vector<IndexType> &permutation) {
m_depth = std::max(depth, m_depth);
size_t count = rangeEnd-rangeStart;
SAssert(count > 0);
if (count == 1) {
/* Create a leaf node */
m_nodes[*rangeStart].setLeaf(true);
permutation[idx] = *rangeStart;
return;
}
typename std::vector<IndexType>::iterator split
= rangeStart + leftSubtreeSize(count);
int axis = m_aabb.getLargestAxis();
std::nth_element(rangeStart, split, rangeEnd,
CoordinateOrdering(m_nodes, axis));
NodeType &splitNode = m_nodes[*split];
splitNode.setAxis(axis);
splitNode.setLeaf(false);
permutation[idx] = *split;
/* Recursively build the children */
Scalar temp = m_aabb.max[axis],
splitPos = splitNode.getPosition()[axis];
m_aabb.max[axis] = splitPos;
buildLB(2*idx+1, depth+1, base, rangeStart, split, permutation);
m_aabb.max[axis] = temp;
if (split+1 != rangeEnd) {
temp = m_aabb.min[axis];
m_aabb.min[axis] = splitPos;
buildLB(2*idx+2, depth+1, base, split+1, rangeEnd, permutation);
m_aabb.min[axis] = temp;
}
}
void build(size_t depth,
typename std::vector<KDNode>::iterator rangeStart,
typename std::vector<KDNode>::iterator rangeEnd) {
typename std::vector<IndexType>::iterator base,
typename std::vector<IndexType>::iterator rangeStart,
typename std::vector<IndexType>::iterator rangeEnd) {
m_depth = std::max(depth, m_depth);
if (rangeEnd-rangeStart <= 0) {
SLog(EError, "Internal error!");
} else if (rangeEnd-rangeStart == 1) {
size_t count = rangeEnd-rangeStart;
SAssert(count > 0);
if (count == 1) {
/* Create a leaf node */
rangeStart->setLeaf(true);
m_nodes[*rangeStart].setLeaf(true);
return;
}
int axis = 0;
typename std::vector<KDNode>::iterator split;
typename std::vector<IndexType>::iterator split;
switch (m_heuristic) {
case EBalanced: {
/* Split along the median */
split = rangeStart + (rangeEnd-rangeStart)/2;
split = rangeStart + count/2;
axis = m_aabb.getLargestAxis();
std::nth_element(rangeStart, split, rangeEnd, CoordinateOrdering(axis));
std::nth_element(rangeStart, split, rangeEnd,
CoordinateOrdering(m_nodes, axis));
};
break;
case ELeftBalanced: {
size_t treeSize = rangeEnd-rangeStart;
/* Layer 0 contains one node */
size_t p = 1;
/* Traverse downwards until the first incompletely
filled tree level is encountered */
while (2*p <= treeSize)
p *= 2;
/* Calculate the number of filled slots in the last level */
size_t remaining = treeSize - p + 1;
if (2*remaining < p) {
/* Case 2: The last level contains too few nodes. Remove
overestimate from the left subtree node count and add
the remaining nodes */
p = (p >> 1) + remaining;
}
split = rangeStart + leftSubtreeSize(count);
axis = m_aabb.getLargestAxis();
split = rangeStart + (p - 1);
std::nth_element(rangeStart, split, rangeEnd,
CoordinateOrdering(axis));
CoordinateOrdering(m_nodes, axis));
};
break;
@ -550,11 +804,11 @@ protected:
/* Sliding midpoint rule: find a split that is close to the spatial median */
axis = m_aabb.getLargestAxis();
value_type midpoint = (value_type) 0.5f
Scalar midpoint = (Scalar) 0.5f
* (m_aabb.max[axis]+m_aabb.min[axis]);
size_t nLT = std::count_if(rangeStart, rangeEnd,
LessThanOrEqual(axis, midpoint));
LessThanOrEqual(m_nodes, axis, midpoint));
/* Re-adjust the split to pass through a nearby point */
split = rangeStart + nLT;
@ -565,23 +819,25 @@ protected:
--split;
std::nth_element(rangeStart, split, rangeEnd,
CoordinateOrdering(axis));
CoordinateOrdering(m_nodes, axis));
};
break;
case EVoxelVolume: {
Float bestCost = std::numeric_limits<Float>::infinity();
for (int dim=0; dim<point_type::dim; ++dim) {
std::sort(rangeStart, rangeEnd, CoordinateOrdering(dim));
for (int dim=0; dim<PointType::dim; ++dim) {
std::sort(rangeStart, rangeEnd,
CoordinateOrdering(m_nodes, dim));
size_t numLeft = 1, numRight = rangeEnd-rangeStart-2;
aabb_type leftAABB(m_aabb), rightAABB(m_aabb);
size_t numLeft = 1, numRight = count-2;
AABBType leftAABB(m_aabb), rightAABB(m_aabb);
Float invVolume = 1.0f / m_aabb.getVolume();
for (typename std::vector<KDNode>::iterator it = rangeStart+1; it != rangeEnd; ++it) {
for (typename std::vector<IndexType>::iterator it = rangeStart+1;
it != rangeEnd; ++it) {
++numLeft; --numRight;
leftAABB.max[dim] = it->getPosition()[dim];
rightAABB.min[dim] = it->getPosition()[dim];
Float pos = m_nodes[*it].getPosition()[dim];
leftAABB.max[dim] = rightAABB.min[dim] = pos;
Float cost = (numLeft * leftAABB.getVolume()
+ numRight * rightAABB.getVolume()) * invVolume;
@ -593,38 +849,40 @@ protected:
}
}
std::nth_element(rangeStart, split, rangeEnd,
CoordinateOrdering(axis));
CoordinateOrdering(m_nodes, axis));
};
break;
}
value_type splitPos = split->getPosition()[axis];
split->setAxis(axis);
if (split+1 != rangeEnd)
split->setRightIndex(rangeStart - m_nodes.begin(), (uint32_t) (split + 1 - m_nodes.begin()));
else
split->setRightIndex(rangeStart - m_nodes.begin(), 0);
NodeType &splitNode = m_nodes[*split];
splitNode.setAxis(axis);
splitNode.setLeaf(false);
split->setLeftIndex(rangeStart - m_nodes.begin(), rangeStart + 1 - m_nodes.begin());
split->setLeaf(false);
if (split+1 != rangeEnd)
splitNode.setRightIndex(rangeStart - base, split + 1 - base);
else
splitNode.setRightIndex(rangeStart - base, 0);
splitNode.setLeftIndex(rangeStart - base, rangeStart + 1 - base);
std::iter_swap(rangeStart, split);
/* Recursively build the children */
value_type temp = m_aabb.max[axis];
Scalar temp = m_aabb.max[axis],
splitPos = splitNode.getPosition()[axis];
m_aabb.max[axis] = splitPos;
build(depth+1, rangeStart+1, split+1);
build(depth+1, base, rangeStart+1, split+1);
m_aabb.max[axis] = temp;
if (split+1 != rangeEnd) {
temp = m_aabb.min[axis];
m_aabb.min[axis] = splitPos;
build(depth+1, split+1, rangeEnd);
build(depth+1, base, split+1, rangeEnd);
m_aabb.min[axis] = temp;
}
}
protected:
std::vector<KDNode> m_nodes;
aabb_type m_aabb;
std::vector<NodeType> m_nodes;
AABBType m_aabb;
EHeuristic m_heuristic;
size_t m_depth;
};

20
include/mitsuba/core/util.h
View File

@ -160,20 +160,28 @@ template<typename T> inline T endianness_swap(T value) {
* in linear time without requiring additional memory. This is based on
* the fact that each permutation can be decomposed into a disjoint set
* of permutations, which can then be applied individually.
*
* \param data
* Pointer to the data that should be permuted
* \param perm
* Input permutation vector having the same size as \c data. After
* the function terminates, this vector will be set to the
* identity permutation.
*/
template <typename T> void permute_inplace(T *values, std::vector<size_t> &perm) {
template <typename DataType, typename IndexType> void permute_inplace(
DataType *data, std::vector<IndexType> &perm) {
for (size_t i=0; i<perm.size(); i++) {
if (perm[i] != i) {
/* The start of a new cycle has been found. Save
the value at this position, since it will be
overwritten */
size_t j = i;
T curval = values[i];
IndexType j = i;
DataType curval = data[i];
do {
/* Shuffle backwards */
size_t k = perm[j];
values[j] = values[k];
IndexType k = perm[j];
data[j] = data[k];
/* Also fix the permutations on the way */
perm[j] = j;
@ -183,7 +191,7 @@ template <typename T> void permute_inplace(T *values, std::vector<size_t> &perm)
} while (perm[j] != i);
/* Fix the final position with the saved value */
values[j] = curval;
data[j] = curval;
perm[j] = j;
}
}

134
include/mitsuba/render/photon.h
View File

@ -20,7 +20,7 @@
#define __PHOTON_H
#include <mitsuba/core/serialization.h>
#include <mitsuba/core/aabb.h>
#include <mitsuba/core/kdtree.h>
MTS_NAMESPACE_BEGIN
@ -171,9 +171,141 @@ protected:
static bool createPrecompTables();
};
/**
* \sa NewPhoton
*/
struct PhotonData {
#if defined(SINGLE_PRECISION) && SPECTRUM_SAMPLES == 3
uint8_t power[4]; //!< Photon power stored in Greg Ward's RGBE format
#else
Spectrum power; //!< Accurate spectral photon power representation
#endif
uint8_t theta; //!< Discretized photon direction (\a theta component)
uint8_t phi; //!< Discretized photon direction (\a phi component)
uint8_t thetaN; //!< Discretized surface normal (\a theta component)
uint8_t phiN; //!< Discretized surface normal (\a phi component)
uint16_t depth; //!< Photon depth (number of preceding interactions)
};
/** \brief Memory-efficient photon representation for use with
* \ref PointKDTree
*
* Requires 24 bytes of memory when Mitsuba is compiled with single
* precision and RGB-based color spectra.
*
* \ingroup librender
*/
struct MTS_EXPORT_RENDER NewPhoton : public SimpleKDNode<Point, PhotonData> {
friend class PhotonMap;
public:
/// Dummy constructor
inline NewPhoton() { }
/// Construct from a photon interaction
NewPhoton(const Point &pos, const Normal &normal,
const Vector &dir, const Spectrum &power,
uint16_t depth);
/// Unserialize from a binary data stream
NewPhoton(Stream *stream);
/// @}
// ======================================================================
/// Return the depth (in # of interactions)
inline int getDepth() const {
return data.depth;
}
/**
* Convert the photon direction from quantized spherical coordinates
* to a floating point vector value. Precomputation idea based on
* Jensen's implementation.
*/
inline Vector getDirection() const {
return Vector(
m_cosPhi[data.phi] * m_sinTheta[data.theta],
m_sinPhi[data.phi] * m_sinTheta[data.theta],
m_cosTheta[data.theta]
);
}
/**
* Convert the normal direction from quantized spherical coordinates
* to a floating point vector value.
*/
inline Normal getNormal() const {
return Normal(
m_cosPhi[data.phiN] * m_sinTheta[data.thetaN],
m_sinPhi[data.phiN] * m_sinTheta[data.thetaN],
m_cosTheta[data.thetaN]
);
}
/// Convert the photon power from RGBE to floating point
inline Spectrum getPower() const {
#if defined(SINGLE_PRECISION) && SPECTRUM_SAMPLES == 3
Spectrum result;
result.fromRGBE(data.power);
return result;
#else
return data.power;
#endif
}
/// Serialize to a binary data stream
inline void serialize(Stream *stream) const {
position.serialize(stream);
#if defined(SINGLE_PRECISION) && SPECTRUM_SAMPLES == 3
stream->write(data.power, 8);
#else
data.power.serialize(stream);
stream->writeUChar(data.phi);
stream->writeUChar(data.theta);
stream->writeUChar(data.phiN);
stream->writeUChar(data.thetaN);
#endif
stream->writeUShort(data.depth);
stream->writeUChar(flags);
}
/// Return a string representation (for debugging)
std::string toString() const {
std::ostringstream oss;
oss << "Photon[pos = " << getPosition().toString() << ","
<< ", power = " << getPower().toString()
<< ", direction = " << getDirection().toString()
<< ", normal = " << getNormal().toString()
<< ", axis = " << getAxis()
<< ", depth = " << getDepth()
<< "]";
return oss.str();
}
protected:
// ======================================================================
/// @{ \name Precomputed lookup tables
// ======================================================================
static Float m_cosTheta[256];
static Float m_sinTheta[256];
static Float m_cosPhi[256];
static Float m_sinPhi[256];
static Float m_expTable[256];
static bool m_precompTableReady;
/// @}
// ======================================================================
/// Initialize the precomputed lookup tables
static bool createPrecompTables();
};
#if defined(SINGLE_PRECISION) && SPECTRUM_SAMPLES == 3
/* Compiler sanity check */
BOOST_STATIC_ASSERT(sizeof(Photon) == 24);
//BOOST_STATIC_ASSERT(sizeof(NewPhoton) == 24);
#endif
MTS_NAMESPACE_END

84
src/librender/photon.cpp
View File

@ -31,6 +31,15 @@ Float Photon::m_expTable[256];
bool Photon::m_precompTableReady = Photon::createPrecompTables();
Float NewPhoton::m_cosTheta[256];
Float NewPhoton::m_sinTheta[256];
Float NewPhoton::m_cosPhi[256];
Float NewPhoton::m_sinPhi[256];
Float NewPhoton::m_expTable[256];
bool NewPhoton::m_precompTableReady = NewPhoton::createPrecompTables();
bool Photon::createPrecompTables() {
for (int i=0; i<256; i++) {
Float angle = (Float) i * ((Float) M_PI / 256.0f);
@ -109,4 +118,79 @@ Photon::Photon(const Point &p, const Normal &normal,
#endif
}
bool NewPhoton::createPrecompTables() {
for (int i=0; i<256; i++) {
Float angle = (Float) i * ((Float) M_PI / 256.0f);
m_cosPhi[i] = std::cos(2.0f * angle);
m_sinPhi[i] = std::sin(2.0f * angle);
m_cosTheta[i] = std::cos(angle);
m_sinTheta[i] = std::sin(angle);
m_expTable[i] = std::ldexp((Float) 1, i - (128+8));
}
m_expTable[0] = 0;
return true;
}
NewPhoton::NewPhoton(Stream *stream) {
position = Point(stream);
#if defined(SINGLE_PRECISION) && SPECTRUM_SAMPLES == 3
stream->read(data.power, 8);
#else
data.power = Spectrum(stream);
data.phi = stream->readUChar();
data.theta = stream->readUChar();
data.phiN = stream->readUChar();
data.thetaN = stream->readUChar();
#endif
data.depth = stream->readUShort();
flags = stream->readUChar();
}
NewPhoton::NewPhoton(const Point &p, const Normal &normal,
const Vector &dir, const Spectrum &P,
uint16_t _depth) {
if (!P.isValid())
SLog(EWarn, "Creating an invalid photon with power: %s", P.toString().c_str());
/* Possibly convert to single precision floating point
(if Mitsuba is configured to use double precision) */
position = p;
data.depth = _depth;
flags = 0;
/* Convert the direction into an approximate spherical
coordinate format to reduce storage requirements */
data.theta = (uint8_t) std::min(255,
(int) (std::acos(dir.z) * (256.0 / M_PI)));
int tmp = std::min(255,
(int) (std::atan2(dir.y, dir.x) * (256.0 / (2.0 * M_PI))));
if (tmp < 0)
data.phi = (uint8_t) (tmp + 256);
else
data.phi = (uint8_t) tmp;
if (normal.isZero()) {
data.thetaN = data.phiN = 0;
} else {
data.thetaN = (uint8_t) std::min(255,
(int) (std::acos(normal.z) * (256.0 / M_PI)));
tmp = std::min(255,
(int) (std::atan2(normal.y, normal.x) * (256.0 / (2.0 * M_PI))));
if (tmp < 0)
data.phiN = (uint8_t) (tmp + 256);
else
data.phiN = (uint8_t) tmp;
}
#if defined(SINGLE_PRECISION) && SPECTRUM_SAMPLES == 3
/* Pack the photon power into Greg Ward's RGBE format */
P.toRGBE(data.power);
#else
data.power = P;
#endif
}
MTS_NAMESPACE_END

2
src/librender/photonmap.cpp
View File

@ -183,7 +183,7 @@ void PhotonMap::quickPartition(photon_iterator left, photon_iterator right,
}
/**
* Given a number of entries, this method calculates the maximum amount of
* Given a number of entries, this method calculates the number of nodes
* nodes on the left subtree of a left-balanced tree. There are two main
* cases here:
*

54
src/tests/test_kd.cpp
View File

@ -26,8 +26,8 @@ MTS_NAMESPACE_BEGIN
class TestKDTree : public TestCase {
public:
MTS_BEGIN_TESTCASE()
MTS_DECLARE_TEST(test01_sutherlandHodgman)
MTS_DECLARE_TEST(test02_bunnyBenchmark)
// MTS_DECLARE_TEST(test01_sutherlandHodgman)
// MTS_DECLARE_TEST(test02_bunnyBenchmark)
MTS_DECLARE_TEST(test03_pointKDTree)
MTS_END_TESTCASE()
@ -127,9 +127,11 @@ public:
Log(EInfo, "");
}
}
void test03_pointKDTree() {
typedef TKDTree< BasicKDNode<Point2, Float> > KDTree2;
typedef PointKDTree< SimpleKDNode<Point2, Float> > KDTree2;
typedef PointKDTree< LeftBalancedKDNode<Point2, Float> > KDTree2Left;
size_t nPoints = 50000, nTries = 20;
ref<Random> random = new Random();
@ -138,7 +140,7 @@ public:
for (size_t i=0; i<nPoints; ++i) {
kdtree[i].setPosition(Point2(random->nextFloat(), random->nextFloat()));
kdtree[i].setValue(random->nextFloat());
kdtree[i].setData(random->nextFloat());
}
std::vector<KDTree2::SearchResult> results, resultsBF;
@ -179,8 +181,48 @@ public:
Point2 p(random->nextFloat(), random->nextFloat());
nTraversals += kdtree.search(p, 0.05, results2);
}
Log(EInfo, "Average number of traversals for a radius=0.05 search query = " SIZE_T_FMT, nTraversals / nTries);
Log(EInfo, "Average number of traversals for a radius=0.05 search query = " SIZE_T_FMT, nTraversals / nTries);
}
Log(EInfo, "Testing the left-balanced kd-tree construction heuristic with left-balanced nodes");
KDTree2Left kdtree(nPoints, KDTree2Left::ELeftBalanced);
for (size_t i=0; i<nPoints; ++i) {
kdtree[i].setPosition(Point2(random->nextFloat(), random->nextFloat()));
kdtree[i].setData(random->nextFloat());
}
std::vector<KDTree2Left::SearchResult> results, resultsBF;
ref<Timer> timer = new Timer();
kdtree.build();
Log(EInfo, "Construction time = %i ms, depth = %i", timer->getMilliseconds(), kdtree.getDepth());
for (int k=1; k<=10; ++k) {
size_t nTraversals = 0;
for (size_t it = 0; it < nTries; ++it) {
Point2 p(random->nextFloat(), random->nextFloat());
nTraversals += kdtree.nnSearch(p, k, results);
resultsBF.clear();
for (size_t j=0; j<nPoints; ++j)
resultsBF.push_back(KDTree2Left::SearchResult((kdtree[j].getPosition()-p).lengthSquared(), (uint32_t) j));
std::sort(results.begin(), results.end(), KDTree2Left::SearchResultComparator());
std::sort(resultsBF.begin(), resultsBF.end(), KDTree2Left::SearchResultComparator());
for (int j=0; j<k; ++j)
assertTrue(results[j] == resultsBF[j]);
}
Log(EInfo, "Average number of traversals for a %i-nn query = " SIZE_T_FMT, k, nTraversals / nTries);
}
std::vector<uint32_t> results2;
size_t nTraversals = 0;
for (size_t it = 0; it < nTries; ++it) {
Point2 p(random->nextFloat(), random->nextFloat());
nTraversals += kdtree.search(p, 0.05, results2);
}
Log(EInfo, "Average number of traversals for a radius=0.05 search query = " SIZE_T_FMT, nTraversals / nTries);
Log(EInfo, "Normal node size = " SIZE_T_FMT " bytes", sizeof(KDTree2::NodeType));
Log(EInfo, "Left-balanced node size = " SIZE_T_FMT " bytes", sizeof(KDTree2Left::NodeType));
}
};

103
src/tests/test_phot.cpp Normal file
View File

@ -0,0 +1,103 @@
/*
This file is part of Mitsuba, a physically based rendering system.
Copyright (c) 2007-2011 by Wenzel Jakob and others.
Mitsuba is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License Version 3
as published by the Free Software Foundation.
Mitsuba is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <mitsuba/core/plugin.h>
#include <mitsuba/core/kdtree.h>
#include <mitsuba/core/fresolver.h>
#include <mitsuba/render/testcase.h>
#include <mitsuba/render/gatherproc.h>
MTS_NAMESPACE_BEGIN
class TestPhotonMap : public TestCase {
public:
MTS_BEGIN_TESTCASE()
MTS_DECLARE_TEST(test01_photonmap)
MTS_END_TESTCASE()
void test01_photonmap() {
ref<FileResolver> fres = Thread::getThread()->getFileResolver();
fres->addPath("scenes/sponza2");
ref<Scene> scene = loadScene("scenes/sponza2/sponza.xml");
scene->configure();
scene->initialize();
size_t nPhotons = 1000000;
ref<Scheduler> sched = Scheduler::getInstance();
ref<GatherPhotonProcess> proc = new GatherPhotonProcess(
GatherPhotonProcess::ESurfacePhotons, nPhotons,
1000, 10, 10, false, true, NULL);
ref<Sampler> sampler = static_cast<Sampler *>(
PluginManager::getInstance()->createObject(Properties("halton")));
int sceneResID = sched->registerResource(scene);
int cameraResID = sched->registerResource(scene->getCamera());
int qmcSamplerID = sched->registerResource(sampler);
proc->bindResource("scene", sceneResID);
proc->bindResource("camera", cameraResID);
proc->bindResource("sampler", qmcSamplerID);
sched->schedule(proc);
sched->wait(proc);
ref<PhotonMap> pmap = proc->getPhotonMap();
pmap->setScaleFactor(1 / (Float) proc->getShotParticles());
PointKDTree<NewPhoton> phot(pmap->getPhotonCount(), PointKDTree<NewPhoton>::EVoxelVolume);
for (size_t i=0; i<pmap->getPhotonCount(); ++i) {
const Photon &photon = pmap->getPhoton(i+1);
phot[i] = NewPhoton(photon.getPosition(),
photon.getNormal(), photon.getDirection(),
photon.getPower(), photon.getDepth());
}
pmap->balance();
phot.build();
ref<Random> random = new Random();
PhotonMap::search_result results[101];
std::vector<PointKDTree<NewPhoton>::SearchResult> results2(101);
size_t nQueries = 100000;
Point *queryPositions = new Point[nQueries];
for (size_t i=0; i<nQueries; ++i)
queryPositions[i] = pmap->getPhoton(random->nextUInt(pmap->getPhotonCount())+1).getPosition();
ref<Timer> timer = new Timer();
Log(EInfo, "Testing query performance (photon map) ..");
size_t nResults1 = 0, nResults2=0;
for (size_t i=0; i<nQueries; ++i) {
Float searchRadius = std::numeric_limits<Float>::infinity();
nResults1 += pmap->nnSearch(queryPositions[i], searchRadius, 100, results);
}
Log(EInfo, "Done. Took %i ms, got " SIZE_T_FMT " results", timer->getMilliseconds(), nResults1);
timer->reset();
Log(EInfo, "Testing query performance .. ");
for (size_t i=0; i<nQueries; ++i) {
Float searchRadius = std::numeric_limits<Float>::infinity();
phot.nnSearch(queryPositions[i], 100, results2, searchRadius);
nResults2 += results2.size();
}
Log(EInfo, "Done. Took %i ms, got " SIZE_T_FMT " results", timer->getMilliseconds(), nResults2);
timer->reset();
}
};
MTS_EXPORT_TESTCASE(TestPhotonMap, "Testcase for kd-tree related code")
MTS_NAMESPACE_END