manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Added missing standard vcg header comment

This commit is contained in:
Paolo Cignoni 2014-07-12 05:49:07 +00:00
parent d520fe2f0e
commit c2f98a8595
1 changed files with 442 additions and 419 deletions
vcg/space/index/kdtree
kdtree.h

861
vcg/space/index/kdtree/kdtree.h
View File

@ -1,3 +1,26 @@
/****************************************************************************
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2004 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program 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 (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
#ifndef KDTREE_VCG_H #ifndef KDTREE_VCG_H
#define KDTREE_VCG_H #define KDTREE_VCG_H
@ -11,476 +34,476 @@
namespace vcg { namespace vcg {
template<typename _DataType> template<typename _DataType>
class ConstDataWrapper class ConstDataWrapper
{ {
public: public:
typedef _DataType DataType; typedef _DataType DataType;
inline ConstDataWrapper() inline ConstDataWrapper()
: mpData(0), mStride(0), mSize(0) : mpData(0), mStride(0), mSize(0)
{} {}
inline ConstDataWrapper(const DataType* pData, int size, int stride = sizeof(DataType)) inline ConstDataWrapper(const DataType* pData, int size, int stride = sizeof(DataType))
: mpData(reinterpret_cast<const unsigned char*>(pData)), mStride(stride), mSize(size) : mpData(reinterpret_cast<const unsigned char*>(pData)), mStride(stride), mSize(size)
{} {}
inline const DataType& operator[] (int i) const inline const DataType& operator[] (int i) const
{ {
return *reinterpret_cast<const DataType*>(mpData + i*mStride); return *reinterpret_cast<const DataType*>(mpData + i*mStride);
} }
inline size_t size() const { return mSize; } inline size_t size() const { return mSize; }
protected: protected:
const unsigned char* mpData; const unsigned char* mpData;
int mStride; int mStride;
size_t mSize; size_t mSize;
}; };
template<class StdVectorType> template<class StdVectorType>
class VectorConstDataWrapper :public ConstDataWrapper<typename StdVectorType::value_type> class VectorConstDataWrapper :public ConstDataWrapper<typename StdVectorType::value_type>
{ {
public: public:
inline VectorConstDataWrapper(StdVectorType &vec): inline VectorConstDataWrapper(StdVectorType &vec):
ConstDataWrapper<typename StdVectorType::value_type> ( &(vec[0]), vec.size(), sizeof(typename StdVectorType::value_type)) ConstDataWrapper<typename StdVectorType::value_type> ( &(vec[0]), vec.size(), sizeof(typename StdVectorType::value_type))
{} {}
}; };
template<class MeshType> template<class MeshType>
class VertexConstDataWrapper :public ConstDataWrapper<typename MeshType::CoordType> class VertexConstDataWrapper :public ConstDataWrapper<typename MeshType::CoordType>
{ {
public: public:
inline VertexConstDataWrapper(MeshType &m): inline VertexConstDataWrapper(MeshType &m):
ConstDataWrapper<typename MeshType::CoordType> ( &(m.vert[0].P()), m.vert.size(), sizeof(typename MeshType::VertexType)) ConstDataWrapper<typename MeshType::CoordType> ( &(m.vert[0].P()), m.vert.size(), sizeof(typename MeshType::VertexType))
{} {}
}; };
/** /**
* This class allows to create a Kd-Tree thought to perform the neighbour query (radius search, knn-nearest serach and closest search). * This class allows to create a Kd-Tree thought to perform the neighbour query (radius search, knn-nearest serach and closest search).
* The class implemetantion is thread-safe. * The class implemetantion is thread-safe.
*/ */
template<typename _Scalar> template<typename _Scalar>
class KdTree class KdTree
{ {
public: public:
typedef _Scalar Scalar; typedef _Scalar Scalar;
typedef vcg::Point3<Scalar> VectorType; typedef vcg::Point3<Scalar> VectorType;
typedef vcg::Box3<Scalar> AxisAlignedBoxType; typedef vcg::Box3<Scalar> AxisAlignedBoxType;
typedef HeapMaxPriorityQueue<int, Scalar> PriorityQueue; typedef HeapMaxPriorityQueue<int, Scalar> PriorityQueue;
struct Node struct Node
{ {
union { union {
//standard node //standard node
struct { struct {
Scalar splitValue; Scalar splitValue;
unsigned int firstChildId:24; unsigned int firstChildId:24;
unsigned int dim:2; unsigned int dim:2;
unsigned int leaf:1; unsigned int leaf:1;
}; };
//leaf //leaf
struct { struct {
unsigned int start; unsigned int start;
unsigned short size; unsigned short size;
}; };
}; };
}; };
typedef std::vector<Node> NodeList; typedef std::vector<Node> NodeList;
// return the protected members which store the nodes and the points list // return the protected members which store the nodes and the points list
inline const NodeList& _getNodes(void) { return mNodes; } inline const NodeList& _getNodes(void) { return mNodes; }
inline const std::vector<VectorType>& _getPoints(void) { return mPoints; } inline const std::vector<VectorType>& _getPoints(void) { return mPoints; }
public: public:
KdTree(const ConstDataWrapper<VectorType>& points, unsigned int nofPointsPerCell = 16, unsigned int maxDepth = 64); KdTree(const ConstDataWrapper<VectorType>& points, unsigned int nofPointsPerCell = 16, unsigned int maxDepth = 64);
~KdTree(); ~KdTree();
void doQueryK(const VectorType& queryPoint, int k, PriorityQueue& mNeighborQueue); void doQueryK(const VectorType& queryPoint, int k, PriorityQueue& mNeighborQueue);
void doQueryDist(const VectorType& queryPoint, float dist, std::vector<unsigned int>& points, std::vector<Scalar>& sqrareDists); void doQueryDist(const VectorType& queryPoint, float dist, std::vector<unsigned int>& points, std::vector<Scalar>& sqrareDists);
void doQueryClosest(const VectorType& queryPoint, unsigned int& index, Scalar& dist); void doQueryClosest(const VectorType& queryPoint, unsigned int& index, Scalar& dist);
protected: protected:
// element of the stack // element of the stack
struct QueryNode struct QueryNode
{ {
QueryNode() {} QueryNode() {}
QueryNode(unsigned int id) : nodeId(id) {} QueryNode(unsigned int id) : nodeId(id) {}
unsigned int nodeId; // id of the next node unsigned int nodeId; // id of the next node
Scalar sq; // squared distance to the next node Scalar sq; // squared distance to the next node
}; };
// used to build the tree: split the subset [start..end[ according to dim and splitValue, // used to build the tree: split the subset [start..end[ according to dim and splitValue,
// and returns the index of the first element of the second subset // and returns the index of the first element of the second subset
unsigned int split(int start, int end, unsigned int dim, float splitValue); unsigned int split(int start, int end, unsigned int dim, float splitValue);
void createTree(unsigned int nodeId, unsigned int start, unsigned int end, unsigned int level, unsigned int targetCellsize, unsigned int targetMaxDepth); void createTree(unsigned int nodeId, unsigned int start, unsigned int end, unsigned int level, unsigned int targetCellsize, unsigned int targetMaxDepth);
protected: protected:
AxisAlignedBoxType mAABB; //BoundingBox AxisAlignedBoxType mAABB; //BoundingBox
NodeList mNodes; //kd-tree nodes NodeList mNodes; //kd-tree nodes
std::vector<VectorType> mPoints; //points read from the input DataWrapper std::vector<VectorType> mPoints; //points read from the input DataWrapper
std::vector<unsigned int> mIndices; //points indices std::vector<unsigned int> mIndices; //points indices
}; };
template<typename Scalar> template<typename Scalar>
KdTree<Scalar>::KdTree(const ConstDataWrapper<VectorType>& points, unsigned int nofPointsPerCell, unsigned int maxDepth) KdTree<Scalar>::KdTree(const ConstDataWrapper<VectorType>& points, unsigned int nofPointsPerCell, unsigned int maxDepth)
: mPoints(points.size()), mIndices(points.size()) : mPoints(points.size()), mIndices(points.size())
{ {
// compute the AABB of the input // compute the AABB of the input
mPoints[0] = points[0]; mPoints[0] = points[0];
mAABB.Set(mPoints[0]); mAABB.Set(mPoints[0]);
for (unsigned int i=1 ; i<mPoints.size() ; ++i) for (unsigned int i=1 ; i<mPoints.size() ; ++i)
{ {
mPoints[i] = points[i]; mPoints[i] = points[i];
mIndices[i] = i; mIndices[i] = i;
mAABB.Add(mPoints[i]); mAABB.Add(mPoints[i]);
} }
mNodes.reserve(4*mPoints.size()/nofPointsPerCell); mNodes.reserve(4*mPoints.size()/nofPointsPerCell);
//first node inserted (no leaf). The others are made by the createTree function (recursively) //first node inserted (no leaf). The others are made by the createTree function (recursively)
mNodes.resize(1); mNodes.resize(1);
mNodes.back().leaf = 0; mNodes.back().leaf = 0;
createTree(0, 0, mPoints.size(), 1, nofPointsPerCell, maxDepth); createTree(0, 0, mPoints.size(), 1, nofPointsPerCell, maxDepth);
} }
template<typename Scalar> template<typename Scalar>
KdTree<Scalar>::~KdTree() KdTree<Scalar>::~KdTree()
{ {
} }
/** Performs the kNN query. /** Performs the kNN query.
* *
* This algorithm uses the simple distance to the split plane to prune nodes. * This algorithm uses the simple distance to the split plane to prune nodes.
* A more elaborated approach consists to track the closest corner of the cell * A more elaborated approach consists to track the closest corner of the cell
* relatively to the current query point. This strategy allows to save about 5% * relatively to the current query point. This strategy allows to save about 5%
* of the leaves. However, in practice the slight overhead due to this tracking * of the leaves. However, in practice the slight overhead due to this tracking
* reduces the overall performance. * reduces the overall performance.
* *
* This algorithm also use a simple stack while a priority queue using the squared * This algorithm also use a simple stack while a priority queue using the squared
* distances to the cells as a priority values allows to save about 10% of the leaves. * distances to the cells as a priority values allows to save about 10% of the leaves.
* But, again, priority queue insertions and deletions are quite involved, and therefore * But, again, priority queue insertions and deletions are quite involved, and therefore
* a simple stack is by far much faster. * a simple stack is by far much faster.
* *
* The result of the query, the k-nearest neighbors, are stored into the stack mNeighborQueue, where the * The result of the query, the k-nearest neighbors, are stored into the stack mNeighborQueue, where the
* topmost element [0] is NOT the nearest but the farthest!! (they are not sorted but arranged into a heap). * topmost element [0] is NOT the nearest but the farthest!! (they are not sorted but arranged into a heap).
*/ */
template<typename Scalar> template<typename Scalar>
void KdTree<Scalar>::doQueryK(const VectorType& queryPoint, int k, PriorityQueue& mNeighborQueue) void KdTree<Scalar>::doQueryK(const VectorType& queryPoint, int k, PriorityQueue& mNeighborQueue)
{ {
mNeighborQueue.setMaxSize(k); mNeighborQueue.setMaxSize(k);
mNeighborQueue.init(); mNeighborQueue.init();
mNeighborQueue.insert(0xffffffff, std::numeric_limits<Scalar>::max()); mNeighborQueue.insert(0xffffffff, std::numeric_limits<Scalar>::max());
QueryNode mNodeStack[64]; QueryNode mNodeStack[64];
mNodeStack[0].nodeId = 0; mNodeStack[0].nodeId = 0;
mNodeStack[0].sq = 0.f; mNodeStack[0].sq = 0.f;
unsigned int count = 1; unsigned int count = 1;
while (count) while (count)
{ {
//we select the last node (AABB) inserted in the stack //we select the last node (AABB) inserted in the stack
QueryNode& qnode = mNodeStack[count-1]; QueryNode& qnode = mNodeStack[count-1];
//while going down the tree qnode.nodeId is the nearest sub-tree, otherwise, //while going down the tree qnode.nodeId is the nearest sub-tree, otherwise,
//in backtracking, qnode.nodeId is the other sub-tree that will be visited iff //in backtracking, qnode.nodeId is the other sub-tree that will be visited iff
//the actual nearest node is further than the split distance. //the actual nearest node is further than the split distance.
Node& node = mNodes[qnode.nodeId]; Node& node = mNodes[qnode.nodeId];
//if the distance is less than the top of the max-heap, it could be one of the k-nearest neighbours //if the distance is less than the top of the max-heap, it could be one of the k-nearest neighbours
if (qnode.sq < mNeighborQueue.getTopWeight()) if (qnode.sq < mNeighborQueue.getTopWeight())
{ {
//when we arrive to a lef //when we arrive to a lef
if (node.leaf) if (node.leaf)
{ {
--count; //pop of the leaf --count; //pop of the leaf
//end is the index of the last element of the leaf in mPoints //end is the index of the last element of the leaf in mPoints
unsigned int end = node.start+node.size; unsigned int end = node.start+node.size;
//adding the element of the leaf to the heap //adding the element of the leaf to the heap
for (unsigned int i=node.start ; i<end ; ++i) for (unsigned int i=node.start ; i<end ; ++i)
mNeighborQueue.insert(mIndices[i], vcg::SquaredNorm(queryPoint - mPoints[i])); mNeighborQueue.insert(mIndices[i], vcg::SquaredNorm(queryPoint - mPoints[i]));
} }
//otherwise, if we're not on a leaf //otherwise, if we're not on a leaf
else else
{ {
// the new offset is the distance between the searched point and the actual split coordinate // the new offset is the distance between the searched point and the actual split coordinate
float new_off = queryPoint[node.dim] - node.splitValue; float new_off = queryPoint[node.dim] - node.splitValue;
//left sub-tree //left sub-tree
if (new_off < 0.) if (new_off < 0.)
{ {
mNodeStack[count].nodeId = node.firstChildId; mNodeStack[count].nodeId = node.firstChildId;
//in the father's nodeId we save the index of the other sub-tree (for backtracking) //in the father's nodeId we save the index of the other sub-tree (for backtracking)
qnode.nodeId = node.firstChildId+1; qnode.nodeId = node.firstChildId+1;
} }
//right sub-tree (same as above) //right sub-tree (same as above)
else else
{ {
mNodeStack[count].nodeId = node.firstChildId+1; mNodeStack[count].nodeId = node.firstChildId+1;
qnode.nodeId = node.firstChildId; qnode.nodeId = node.firstChildId;
} }
//distance is inherited from the father (while descending the tree it's equal to 0) //distance is inherited from the father (while descending the tree it's equal to 0)
mNodeStack[count].sq = qnode.sq; mNodeStack[count].sq = qnode.sq;
//distance of the father is the squared distance from the split plane //distance of the father is the squared distance from the split plane
qnode.sq = new_off*new_off; qnode.sq = new_off*new_off;
++count; ++count;
} }
} }
else else
{ {
// pop // pop
--count; --count;
} }
} }
} }
/** Performs the distance query. /** Performs the distance query.
* *
* The result of the query, all the points within the distance dist form the query point, is the vector of the indeces * The result of the query, all the points within the distance dist form the query point, is the vector of the indeces
* and the vector of the squared distances from the query point. * and the vector of the squared distances from the query point.
*/ */
template<typename Scalar> template<typename Scalar>
void KdTree<Scalar>::doQueryDist(const VectorType& queryPoint, float dist, std::vector<unsigned int>& points, std::vector<Scalar>& sqrareDists) void KdTree<Scalar>::doQueryDist(const VectorType& queryPoint, float dist, std::vector<unsigned int>& points, std::vector<Scalar>& sqrareDists)
{ {
QueryNode mNodeStack[64]; QueryNode mNodeStack[64];
mNodeStack[0].nodeId = 0; mNodeStack[0].nodeId = 0;
mNodeStack[0].sq = 0.f; mNodeStack[0].sq = 0.f;
unsigned int count = 1; unsigned int count = 1;
float sqrareDist = dist*dist; float sqrareDist = dist*dist;
while (count) while (count)
{ {
QueryNode& qnode = mNodeStack[count-1]; QueryNode& qnode = mNodeStack[count-1];
Node & node = mNodes[qnode.nodeId]; Node & node = mNodes[qnode.nodeId];
if (qnode.sq < sqrareDist) if (qnode.sq < sqrareDist)
{ {
if (node.leaf) if (node.leaf)
{ {
--count; // pop --count; // pop
unsigned int end = node.start+node.size; unsigned int end = node.start+node.size;
for (unsigned int i=node.start ; i<end ; ++i) for (unsigned int i=node.start ; i<end ; ++i)
{ {
float pointSquareDist = vcg::SquaredNorm(queryPoint - mPoints[i]); float pointSquareDist = vcg::SquaredNorm(queryPoint - mPoints[i]);
if (pointSquareDist < sqrareDist) if (pointSquareDist < sqrareDist)
{ {
points.push_back(mIndices[i]); points.push_back(mIndices[i]);
sqrareDists.push_back(pointSquareDist); sqrareDists.push_back(pointSquareDist);
} }
} }
} }
else else
{ {
// replace the stack top by the farthest and push the closest // replace the stack top by the farthest and push the closest
float new_off = queryPoint[node.dim] - node.splitValue; float new_off = queryPoint[node.dim] - node.splitValue;
if (new_off < 0.) if (new_off < 0.)
{ {
mNodeStack[count].nodeId = node.firstChildId; mNodeStack[count].nodeId = node.firstChildId;
qnode.nodeId = node.firstChildId+1; qnode.nodeId = node.firstChildId+1;
} }
else else
{ {
mNodeStack[count].nodeId = node.firstChildId+1; mNodeStack[count].nodeId = node.firstChildId+1;
qnode.nodeId = node.firstChildId; qnode.nodeId = node.firstChildId;
} }
mNodeStack[count].sq = qnode.sq; mNodeStack[count].sq = qnode.sq;
qnode.sq = new_off*new_off; qnode.sq = new_off*new_off;
++count; ++count;
} }
} }
else else
{ {
// pop // pop
--count; --count;
} }
} }
} }
/** Searchs the closest point. /** Searchs the closest point.
* *
* The result of the query, the closest point to the query point, is the index of the point and * The result of the query, the closest point to the query point, is the index of the point and
* and the squared distance from the query point. * and the squared distance from the query point.
*/ */
template<typename Scalar> template<typename Scalar>
void KdTree<Scalar>::doQueryClosest(const VectorType& queryPoint, unsigned int& index, Scalar& dist) void KdTree<Scalar>::doQueryClosest(const VectorType& queryPoint, unsigned int& index, Scalar& dist)
{ {
QueryNode mNodeStack[64]; QueryNode mNodeStack[64];
mNodeStack[0].nodeId = 0; mNodeStack[0].nodeId = 0;
mNodeStack[0].sq = 0.f; mNodeStack[0].sq = 0.f;
unsigned int count = 1; unsigned int count = 1;
int minIndex = mIndices.size() / 2; int minIndex = mIndices.size() / 2;
Scalar minDist = vcg::SquaredNorm(queryPoint - mPoints[minIndex]); Scalar minDist = vcg::SquaredNorm(queryPoint - mPoints[minIndex]);
minIndex = mIndices[minIndex]; minIndex = mIndices[minIndex];
while (count) while (count)
{ {
QueryNode& qnode = mNodeStack[count-1]; QueryNode& qnode = mNodeStack[count-1];
Node & node = mNodes[qnode.nodeId]; Node & node = mNodes[qnode.nodeId];
if (qnode.sq < minDist) if (qnode.sq < minDist)
{ {
if (node.leaf) if (node.leaf)
{ {
--count; // pop --count; // pop
unsigned int end = node.start+node.size; unsigned int end = node.start+node.size;
for (unsigned int i=node.start ; i<end ; ++i) for (unsigned int i=node.start ; i<end ; ++i)
{ {
float pointSquareDist = vcg::SquaredNorm(queryPoint - mPoints[i]); float pointSquareDist = vcg::SquaredNorm(queryPoint - mPoints[i]);
if (pointSquareDist < minDist) if (pointSquareDist < minDist)
{ {
minDist = pointSquareDist; minDist = pointSquareDist;
minIndex = mIndices[i]; minIndex = mIndices[i];
} }
} }
} }
else else
{ {
// replace the stack top by the farthest and push the closest // replace the stack top by the farthest and push the closest
float new_off = queryPoint[node.dim] - node.splitValue; float new_off = queryPoint[node.dim] - node.splitValue;
if (new_off < 0.) if (new_off < 0.)
{ {
mNodeStack[count].nodeId = node.firstChildId; mNodeStack[count].nodeId = node.firstChildId;
qnode.nodeId = node.firstChildId+1; qnode.nodeId = node.firstChildId+1;
} }
else else
{ {
mNodeStack[count].nodeId = node.firstChildId+1; mNodeStack[count].nodeId = node.firstChildId+1;
qnode.nodeId = node.firstChildId; qnode.nodeId = node.firstChildId;
} }
mNodeStack[count].sq = qnode.sq; mNodeStack[count].sq = qnode.sq;
qnode.sq = new_off*new_off; qnode.sq = new_off*new_off;
++count; ++count;
} }
} }
else else
{ {
// pop // pop
--count; --count;
} }
} }
index = minIndex; index = minIndex;
dist = minDist; dist = minDist;
} }
/** /**
* Split the subarray between start and end in two part, one with the elements less than splitValue, * Split the subarray between start and end in two part, one with the elements less than splitValue,
* the other with the elements greater or equal than splitValue. The elements are compared * the other with the elements greater or equal than splitValue. The elements are compared
* using the "dim" coordinate [0 = x, 1 = y, 2 = z]. * using the "dim" coordinate [0 = x, 1 = y, 2 = z].
*/ */
template<typename Scalar> template<typename Scalar>
unsigned int KdTree<Scalar>::split(int start, int end, unsigned int dim, float splitValue) unsigned int KdTree<Scalar>::split(int start, int end, unsigned int dim, float splitValue)
{ {
int l(start), r(end-1); int l(start), r(end-1);
for ( ; l<r ; ++l, --r) for ( ; l<r ; ++l, --r)
{ {
while (l < end && mPoints[l][dim] < splitValue) while (l < end && mPoints[l][dim] < splitValue)
l++; l++;
while (r >= start && mPoints[r][dim] >= splitValue) while (r >= start && mPoints[r][dim] >= splitValue)
r--; r--;
if (l > r) if (l > r)
break; break;
std::swap(mPoints[l],mPoints[r]); std::swap(mPoints[l],mPoints[r]);
std::swap(mIndices[l],mIndices[r]); std::swap(mIndices[l],mIndices[r]);
} }
//returns the index of the first element on the second part //returns the index of the first element on the second part
return (mPoints[l][dim] < splitValue ? l+1 : l); return (mPoints[l][dim] < splitValue ? l+1 : l);
} }
/** recursively builds the kdtree /** recursively builds the kdtree
* *
* The heuristic is the following: * The heuristic is the following:
* - if the number of points in the node is lower than targetCellsize then make a leaf * - if the number of points in the node is lower than targetCellsize then make a leaf
* - else compute the AABB of the points of the node and split it at the middle of * - else compute the AABB of the points of the node and split it at the middle of
* the largest AABB dimension. * the largest AABB dimension.
* *
* This strategy might look not optimal because it does not explicitly prune empty space, * This strategy might look not optimal because it does not explicitly prune empty space,
* unlike more advanced SAH-like techniques used for RT. On the other hand it leads to a shorter tree, * unlike more advanced SAH-like techniques used for RT. On the other hand it leads to a shorter tree,
* faster to traverse and our experience shown that in the special case of kNN queries, * faster to traverse and our experience shown that in the special case of kNN queries,
* this strategy is indeed more efficient (and much faster to build). Moreover, for volume data * this strategy is indeed more efficient (and much faster to build). Moreover, for volume data
* (e.g., fluid simulation) pruning the empty space is useless. * (e.g., fluid simulation) pruning the empty space is useless.
* *
* Actually, storing at each node the exact AABB (we therefore have a binary BVH) allows * Actually, storing at each node the exact AABB (we therefore have a binary BVH) allows
* to prune only about 10% of the leaves, but the overhead of this pruning (ball/ABBB intersection) * to prune only about 10% of the leaves, but the overhead of this pruning (ball/ABBB intersection)
* is more expensive than the gain it provides and the memory consumption is x4 higher ! * is more expensive than the gain it provides and the memory consumption is x4 higher !
*/ */
template<typename Scalar> template<typename Scalar>
void KdTree<Scalar>::createTree(unsigned int nodeId, unsigned int start, unsigned int end, unsigned int level, unsigned int targetCellSize, unsigned int targetMaxDepth) void KdTree<Scalar>::createTree(unsigned int nodeId, unsigned int start, unsigned int end, unsigned int level, unsigned int targetCellSize, unsigned int targetMaxDepth)
{ {
//select the first node //select the first node
Node& node = mNodes[nodeId]; Node& node = mNodes[nodeId];
AxisAlignedBoxType aabb; AxisAlignedBoxType aabb;
//putting all the points in the bounding box //putting all the points in the bounding box
aabb.Set(mPoints[start]); aabb.Set(mPoints[start]);
for (unsigned int i=start+1 ; i<end ; ++i) for (unsigned int i=start+1 ; i<end ; ++i)
aabb.Add(mPoints[i]); aabb.Add(mPoints[i]);
//bounding box diagonal //bounding box diagonal
VectorType diag = aabb.max - aabb.min; VectorType diag = aabb.max - aabb.min;
//the split "dim" is the dimension of the box with the biggest value //the split "dim" is the dimension of the box with the biggest value
unsigned int dim; unsigned int dim;
if (diag.X() > diag.Y()) if (diag.X() > diag.Y())
dim = diag.X() > diag.Z() ? 0 : 2; dim = diag.X() > diag.Z() ? 0 : 2;
else else
dim = diag.Y() > diag.Z() ? 1 : 2; dim = diag.Y() > diag.Z() ? 1 : 2;
node.dim = dim;
//we divide the bounding box in 2 partitions, considering the average of the "dim" dimension
node.splitValue = Scalar(0.5*(aabb.max[dim] + aabb.min[dim]));
//midId is the index of the first element in the second partition node.dim = dim;
unsigned int midId = split(start, end, dim, node.splitValue); //we divide the bounding box in 2 partitions, considering the average of the "dim" dimension
node.splitValue = Scalar(0.5*(aabb.max[dim] + aabb.min[dim]));
//midId is the index of the first element in the second partition
unsigned int midId = split(start, end, dim, node.splitValue);
node.firstChildId = mNodes.size(); node.firstChildId = mNodes.size();
mNodes.resize(mNodes.size()+2); mNodes.resize(mNodes.size()+2);
{ {
// left child // left child
unsigned int childId = mNodes[nodeId].firstChildId; unsigned int childId = mNodes[nodeId].firstChildId;
Node& child = mNodes[childId]; Node& child = mNodes[childId];
if (midId - start <= targetCellSize || level>=targetMaxDepth) if (midId - start <= targetCellSize || level>=targetMaxDepth)
{ {
child.leaf = 1; child.leaf = 1;
child.start = start; child.start = start;
child.size = midId - start; child.size = midId - start;
} }
else else
{ {
child.leaf = 0; child.leaf = 0;
createTree(childId, start, midId, level+1, targetCellSize, targetMaxDepth); createTree(childId, start, midId, level+1, targetCellSize, targetMaxDepth);
} }
} }
{ {
// right child // right child
unsigned int childId = mNodes[nodeId].firstChildId+1; unsigned int childId = mNodes[nodeId].firstChildId+1;
Node& child = mNodes[childId]; Node& child = mNodes[childId];
if (end - midId <= targetCellSize || level>=targetMaxDepth) if (end - midId <= targetCellSize || level>=targetMaxDepth)
{ {
child.leaf = 1; child.leaf = 1;
child.start = midId; child.start = midId;
child.size = end - midId; child.size = end - midId;
} }
else else
{ {
child.leaf = 0; child.leaf = 0;
createTree(childId, midId, end, level+1, targetCellSize, targetMaxDepth); createTree(childId, midId, end, level+1, targetCellSize, targetMaxDepth);
} }
} }
} }
} }
#endif #endif