"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Kernel/rtree.h" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

rtree.h  (getdp-3.4.0-source.tgz)	:	rtree.h  (getdp-3.5.0-source.tgz)
/*		/*
TITLE		TITLE
R-TREES: A DYNAMIC INDEX STRUCTURE FOR SPATIAL SEARCHING		R-TREES: A DYNAMIC INDEX STRUCTURE FOR SPATIAL SEARCHING
DESCRIPTION		DESCRIPTION
A C++ templated version of the RTree algorithm.		A C++ templated version of the RTree algorithm.
For more information please read the comments in RTree.h		For more information please read the comments in RTree.h
AUTHORS		AUTHORS
* 1983 Original algorithm and test code by Antonin Guttman and Michael St		* 1983 Original algorithm and test code by Antonin Guttman and Michael
onebraker, UC Berkely		Stonebraker, UC Berkely
* 1994 ANCI C ported from original test code by Melinda Green - melinda@s		* 1994 ANCI C ported from original test code by Melinda Green -
uperliminal.com		melinda@superliminal.com
* 1995 Sphere volume fix for degeneracy problem submitted by Paul Brook		* 1995 Sphere volume fix for degeneracy problem submitted by Paul Brook
* 2004 Templated C++ port by Greg Douglas		* 2004 Templated C++ port by Greg Douglas
LICENSE:		LICENSE:
Entirely free for all uses. Enjoy!		Entirely free for all uses. Enjoy!
FILES		FILES
* RTree.h The C++ templated RTree implementation. Well comment		* RTree.h The C++ templated RTree implementation. Well commented.
ed.		* Test.cpp A simple test program, ported from the original C
* Test.cpp A simple test program, ported from the original C ve		version.
rsion.		* MemoryTest.cpp A more rigourous test to validate memory use.
* MemoryTest.cpp A more rigourous test to validate memory use.		* README.TXT This file.
* README.TXT This file.
TO BUILD		TO BUILD
To build a test, compile only one of the test files with RTree.h.		To build a test, compile only one of the test files with RTree.h.
Both test files contain a main() function.		Both test files contain a main() function.
RECENT CHANGE LOG		RECENT CHANGE LOG
05 Jan 2010		05 Jan 2010
o Fixed Iterator GetFirst() - Previous fix was not incomplete		o Fixed Iterator GetFirst() - Previous fix was not incomplete
03 Dec 2009		03 Dec 2009
o Added Iteartor GetBounds()		o Added Iteartor GetBounds()
o Added Iterator usage to simple test		o Added Iterator usage to simple test
o Fixed Iterator GetFirst() - Thanks Mathew Riek		o Fixed Iterator GetFirst() - Thanks Mathew Riek
o Minor updates for MSVC 2005/08 compilers		o Minor updates for MSVC 2005/08 compilers
*/		*/
#ifndef RTREE_H		#ifndef RTREE_H
#define RTREE_H		#define RTREE_H
#include <algorithm>		#include <algorithm>
// NOTE This file compiles under MSVC 6 SP5 and MSVC .Net 2003 it may not work o		// NOTE This file compiles under MSVC 6 SP5 and MSVC .Net 2003 it may not work
n other compilers without modification.		// on other compilers without modification.
// NOTE These next few lines may be win32 specific, you may need to modify them		// NOTE These next few lines may be win32 specific, you may need to modify them
to compile on other platform		// to compile on other platform
#include <stdio.h>		#include <stdio.h>
#include <math.h>		#include <math.h>
#include <assert.h>		#include <assert.h>
#include <stdlib.h>		#include <stdlib.h>
#define ASSERT assert // RTree uses ASSERT( condition )		#define ASSERT assert // RTree uses ASSERT( condition )
//#ifndef Min		//#ifndef Min
// #define Min std::min		// #define Min std::min
//#endif //Min		//#endif //Min
//#ifndef Max		//#ifndef Max
// #define Max std::max		// #define Max std::max
//#endif //Max		//#endif //Max
//		//
// RTree.h		// RTree.h
//		//
#define RTREE_TEMPLATE template<class DATATYPE, class ELEMTYPE, int NUMDIMS, cla		#define RTREE_TEMPLATE \
ss ELEMTYPEREAL, int TMAXNODES, int TMINNODES>		template <class DATATYPE, class ELEMTYPE, int NUMDIMS, class ELEMTYPEREAL, \
#define RTREE_QUAL RTree<DATATYPE, ELEMTYPE, NUMDIMS, ELEMTYPEREAL, TMAXNODES, T		int TMAXNODES, int TMINNODES>
MINNODES>		#define RTREE_QUAL \
		RTree<DATATYPE, ELEMTYPE, NUMDIMS, ELEMTYPEREAL, TMAXNODES, TMINNODES>
#define RTREE_DONT_USE_MEMPOOLS // This version does not contain a fixed memory
allocator, fill in lines with EXAMPLE to implement one.		#define RTREE_DONT_USE_MEMPOOLS // This version does not contain a fixed memory
#define RTREE_USE_SPHERICAL_VOLUME // Better split classification, may be slower		// allocator, fill in lines with EXAMPLE to
on some systems		// implement one.
		#define RTREE_USE_SPHERICAL_VOLUME // Better split classification, may be slower
		// on some systems
// Fwd decl		// Fwd decl
class RTFileStream; // File I/O helper class, look below for implementation and		class RTFileStream; // File I/O helper class, look below for implementation and
notes.		// notes.
/// \class RTree		/// \class RTree
/// Implementation of RTree, a multidimensional bounding rectangle tree.		/// Implementation of RTree, a multidimensional bounding rectangle tree.
/// Example usage: For a 3-dimensional tree use RTree<Object*, float, 3> myTree;		/// Example usage: For a 3-dimensional tree use RTree<Object*, float, 3> myTree;
///		///
/// This modified, templated C++ version by Greg Douglas at Auran (http://www.au		/// This modified, templated C++ version by Greg Douglas at Auran
ran.com)		/// (http://www.auran.com)
///		///
/// DATATYPE Referenced data, should be int, void*, obj* etc. no larger than siz		/// DATATYPE Referenced data, should be int, void*, obj* etc. no larger than
eof<void*> and simple type		/// sizeof<void*> and simple type ELEMTYPE Type of element such as int or float
/// ELEMTYPE Type of element such as int or float
/// NUMDIMS Number of dimensions such as 2 or 3		/// NUMDIMS Number of dimensions such as 2 or 3
/// ELEMTYPEREAL Type of element that allows fractional and large values such as		/// ELEMTYPEREAL Type of element that allows fractional and large values such as
float or double, for use in volume calcs		/// float or double, for use in volume calcs
///		///
/// NOTES: Inserting and removing data requires the knowledge of its constant Mi		/// NOTES: Inserting and removing data requires the knowledge of its constant
nimal Bounding Rectangle.		/// Minimal Bounding Rectangle.
/// This version uses new/delete for nodes, I recommend using a fixed siz		/// This version uses new/delete for nodes, I recommend using a fixed
e allocator for efficiency.		/// size allocator for efficiency. Instead of using a callback function
/// Instead of using a callback function for returned results, I recommen		/// for returned results, I recommend and efficient pre-sized, grow-only
d and efficient pre-sized, grow-only memory		/// memory array similar to MFC CArray or STL Vector for returning search
/// array similar to MFC CArray or STL Vector for returning search query		/// query result.
result.
///		///
template<class DATATYPE, class ELEMTYPE, int NUMDIMS,		template <class DATATYPE, class ELEMTYPE, int NUMDIMS,
class ELEMTYPEREAL = ELEMTYPE, int TMAXNODES = 8, int TMINNODES = TMAXN		class ELEMTYPEREAL = ELEMTYPE, int TMAXNODES = 8,
ODES / 2>		int TMINNODES = TMAXNODES / 2>
class RTree		class RTree {
{
protected:		protected:
		struct Node; // Fwd decl. Used by other internal structs and iterator
struct Node; // Fwd decl. Used by other internal structs and iterator
public:		public:
// These constant must be declared after Branch and before Node struct		// These constant must be declared after Branch and before Node struct
// Stuck up here for MSVC 6 compiler. NSVC .NET 2003 is much happier.		// Stuck up here for MSVC 6 compiler. NSVC .NET 2003 is much happier.
enum		enum {
{		MAXNODES = TMAXNODES, ///< Max elements in node
MAXNODES = TMAXNODES, ///< Max elements in node		MINNODES = TMINNODES, ///< Min elements in node
MINNODES = TMINNODES, ///< Min elements in node
};		};
public:		public:
RTree();		RTree();
virtual ~RTree();		virtual ~RTree();
/// Insert entry		/// Insert entry
/// \param a_min Min of bounding rect		/// \param a_min Min of bounding rect
/// \param a_max Max of bounding rect		/// \param a_max Max of bounding rect
/// \param a_dataId Positive Id of data. Maybe zero, but negative numbers not		/// \param a_dataId Positive Id of data. Maybe zero, but negative numbers not
allowed.		/// allowed.
void Insert(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], cons		void Insert(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS],
t DATATYPE& a_dataId);		const DATATYPE &a_dataId);
/// Remove entry		/// Remove entry
/// \param a_min Min of bounding rect		/// \param a_min Min of bounding rect
/// \param a_max Max of bounding rect		/// \param a_max Max of bounding rect
/// \param a_dataId Positive Id of data. Maybe zero, but negative numbers not		/// \param a_dataId Positive Id of data. Maybe zero, but negative numbers not
allowed.		/// allowed.
void Remove(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], cons		void Remove(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS],
t DATATYPE& a_dataId);		const DATATYPE &a_dataId);
/// Find all within search rectangle		/// Find all within search rectangle
/// \param a_min Min of search bounding rect		/// \param a_min Min of search bounding rect
/// \param a_max Max of search bounding rect		/// \param a_max Max of search bounding rect
/// \param a_resultCallback Callback function to return result. Callback shou		/// \param a_resultCallback Callback function to return result. Callback
ld return 'true' to continue searching		/// should return 'true' to continue searching \param a_context User context
/// \param a_context User context to pass as parameter to a_resultCallback		/// to pass as parameter to a_resultCallback \return Returns the number of
/// \return Returns the number of entries found		/// entries found
int Search(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], bool		int Search(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS],
a_resultCallback(DATATYPE a_data, void* a_context), void* a_context);		bool a_resultCallback(DATATYPE a_data, void *a_context),
		void *a_context);
/// Remove all entries from tree		/// Remove all entries from tree
void RemoveAll();		void RemoveAll();
/// Count the data elements in this container. This is slow as no internal co		/// Count the data elements in this container. This is slow as no internal
unter is maintained.		/// counter is maintained.
int Count();		int Count();
/// Load tree contents from file		/// Load tree contents from file
bool Load(const char* a_fileName);		bool Load(const char *a_fileName);
/// Load tree contents from stream		/// Load tree contents from stream
bool Load(RTFileStream& a_stream);		bool Load(RTFileStream &a_stream);
/// Save tree contents to file		/// Save tree contents to file
bool Save(const char* a_fileName);		bool Save(const char *a_fileName);
/// Save tree contents to stream		/// Save tree contents to stream
bool Save(RTFileStream& a_stream);		bool Save(RTFileStream &a_stream);
/// Iterator is not remove safe.		/// Iterator is not remove safe.
class Iterator		class Iterator {
{
private:		private:
		enum {
		MAX_STACK = 32
		}; // Max stack size. Allows almost n^32 where n is number of branches in
		// node
enum { MAX_STACK = 32 }; // Max stack size. Allows almost n^32 where n is n		struct StackElement {
umber of branches in node		Node *m_node;
struct StackElement
{
Node* m_node;
int m_branchIndex;		int m_branchIndex;
};		};
public:		public:
		Iterator() { Init(); }
Iterator() { Init(); }		~Iterator() {}
~Iterator() { }
/// Is iterator invalid		/// Is iterator invalid
bool IsNull() { return (m_tos <= 0); }		bool IsNull() { return (m_tos <= 0); }
/// Is iterator pointing to valid data		/// Is iterator pointing to valid data
bool IsNotNull() { return (m_tos > 0); }		bool IsNotNull() { return (m_tos > 0); }
/// Access the current data element. Caller must be sure iterator is not NUL		/// Access the current data element. Caller must be sure iterator is not
L first.		/// NULL first.
DATATYPE& operator*()		DATATYPE &operator*()
{		{
ASSERT(IsNotNull());		ASSERT(IsNotNull());
StackElement& curTos = m_stack[m_tos - 1];		StackElement &curTos = m_stack[m_tos - 1];
return curTos.m_node->m_branch[curTos.m_branchIndex].m_data;		return curTos.m_node->m_branch[curTos.m_branchIndex].m_data;
}		}
/// Access the current data element. Caller must be sure iterator is not NUL		/// Access the current data element. Caller must be sure iterator is not
L first.		/// NULL first.
const DATATYPE& operator*() const		const DATATYPE &operator*() const
{		{
ASSERT(IsNotNull());		ASSERT(IsNotNull());
StackElement& curTos = m_stack[m_tos - 1];		StackElement &curTos = m_stack[m_tos - 1];
return curTos.m_node->m_branch[curTos.m_branchIndex].m_data;		return curTos.m_node->m_branch[curTos.m_branchIndex].m_data;
}		}
/// Find the next data element		/// Find the next data element
bool operator++() { return FindNextData(); }		bool operator++() { return FindNextData(); }
/// Get the bounds for this node		/// Get the bounds for this node
void GetBounds(ELEMTYPE a_min[NUMDIMS], ELEMTYPE a_max[NUMDIMS])		void GetBounds(ELEMTYPE a_min[NUMDIMS], ELEMTYPE a_max[NUMDIMS])
{		{
ASSERT(IsNotNull());		ASSERT(IsNotNull());
StackElement& curTos = m_stack[m_tos - 1];		StackElement &curTos = m_stack[m_tos - 1];
Branch& curBranch = curTos.m_node->m_branch[curTos.m_branchIndex];		Branch &curBranch = curTos.m_node->m_branch[curTos.m_branchIndex];
for(int index = 0; index < NUMDIMS; ++index)		for(int index = 0; index < NUMDIMS; ++index) {
{
a_min[index] = curBranch.m_rect.m_min[index];		a_min[index] = curBranch.m_rect.m_min[index];
a_max[index] = curBranch.m_rect.m_max[index];		a_max[index] = curBranch.m_rect.m_max[index];
}		}
}		}
private:		private:
/// Reset iterator		/// Reset iterator
void Init() { m_tos = 0; }		void Init() { m_tos = 0; }
/// Find the next data element in the tree (For internal use only)		/// Find the next data element in the tree (For internal use only)
bool FindNextData()		bool FindNextData()
{		{
for(;;)		for(;;) {
{		if(m_tos <= 0) { return false; }
if(m_tos <= 0)		StackElement curTos =
{		Pop(); // Copy stack top cause it may change as we use it
return false;
}
StackElement curTos = Pop(); // Copy stack top cause it may change as we
use it
if(curTos.m_node->IsLeaf())		if(curTos.m_node->IsLeaf()) {
{
// Keep walking through data while we can		// Keep walking through data while we can
if(curTos.m_branchIndex+1 < curTos.m_node->m_count)		if(curTos.m_branchIndex + 1 < curTos.m_node->m_count) {
{
// There is more data, just point to the next one		// There is more data, just point to the next one
Push(curTos.m_node, curTos.m_branchIndex + 1);		Push(curTos.m_node, curTos.m_branchIndex + 1);
return true;		return true;
}		}
// No more data, so it will fall back to previous level		// No more data, so it will fall back to previous level
}		}
else		else {
{		if(curTos.m_branchIndex + 1 < curTos.m_node->m_count) {
if(curTos.m_branchIndex+1 < curTos.m_node->m_count)
{
// Push sibling on for future tree walk		// Push sibling on for future tree walk
// This is the 'fall back' node when we finish with the current leve		// This is the 'fall back' node when we finish with the current
l		// level
Push(curTos.m_node, curTos.m_branchIndex + 1);		Push(curTos.m_node, curTos.m_branchIndex + 1);
}		}
// Since cur node is not a leaf, push first of next level to get deepe		// Since cur node is not a leaf, push first of next level to get
r into the tree		// deeper into the tree
Node* nextLevelnode = curTos.m_node->m_branch[curTos.m_branchIndex].m_		Node *nextLevelnode =
child;		curTos.m_node->m_branch[curTos.m_branchIndex].m_child;
Push(nextLevelnode, 0);		Push(nextLevelnode, 0);
// If we pushed on a new leaf, exit as the data is ready at TOS		// If we pushed on a new leaf, exit as the data is ready at TOS
if(nextLevelnode->IsLeaf())		if(nextLevelnode->IsLeaf()) { return true; }
{
return true;
}
}		}
}		}
}		}
/// Push node and branch onto iteration stack (For internal use only)		/// Push node and branch onto iteration stack (For internal use only)
void Push(Node* a_node, int a_branchIndex)		void Push(Node *a_node, int a_branchIndex)
{		{
m_stack[m_tos].m_node = a_node;		m_stack[m_tos].m_node = a_node;
m_stack[m_tos].m_branchIndex = a_branchIndex;		m_stack[m_tos].m_branchIndex = a_branchIndex;
++m_tos;		++m_tos;
ASSERT(m_tos <= MAX_STACK);		ASSERT(m_tos <= MAX_STACK);
}		}
/// Pop element off iteration stack (For internal use only)		/// Pop element off iteration stack (For internal use only)
StackElement& Pop()		StackElement &Pop()
{		{
ASSERT(m_tos > 0);		ASSERT(m_tos > 0);
--m_tos;		--m_tos;
return m_stack[m_tos];		return m_stack[m_tos];
}		}
StackElement m_stack[MAX_STACK]; ///< Stack as we are doing ite		StackElement m_stack[MAX_STACK]; ///< Stack as we are doing iteration
ration instead of recursion		///< instead of recursion
int m_tos; ///< Top Of Stack index		int m_tos; ///< Top Of Stack index
friend class RTree; // Allow hiding of non-public functions while allowing m		friend class RTree; // Allow hiding of non-public functions while allowing
anipulation by logical owner		// manipulation by logical owner
};		};
/// Get 'first' for iteration		/// Get 'first' for iteration
void GetFirst(Iterator& a_it)		void GetFirst(Iterator &a_it)
{		{
a_it.Init();		a_it.Init();
Node* first = m_root;		Node *first = m_root;
while(first)		while(first) {
{		if(first->IsInternalNode() && first->m_count > 1) {
if(first->IsInternalNode() && first->m_count > 1)
{
a_it.Push(first, 1); // Descend sibling branch later		a_it.Push(first, 1); // Descend sibling branch later
}		}
else if(first->IsLeaf())		else if(first->IsLeaf()) {
{		if(first->m_count) { a_it.Push(first, 0); }
if(first->m_count)
{
a_it.Push(first, 0);
}
break;		break;
}		}
first = first->m_branch[0].m_child;		first = first->m_branch[0].m_child;
}		}
}		}
/// Get Next for iteration		/// Get Next for iteration
void GetNext(Iterator& a_it) { ++a_it; }		void GetNext(Iterator &a_it) { ++a_it; }
/// Is iterator NULL, or at end?		/// Is iterator NULL, or at end?
bool IsNull(Iterator& a_it) { return a_it.IsNull(); }		bool IsNull(Iterator &a_it) { return a_it.IsNull(); }
/// Get object at iterator position		/// Get object at iterator position
DATATYPE& GetAt(Iterator& a_it) { return *a_it; }		DATATYPE &GetAt(Iterator &a_it) { return *a_it; }
protected:		protected:
/// Minimal bounding rectangle (n-dimensional)		/// Minimal bounding rectangle (n-dimensional)
struct Rect		struct Rect {
{		ELEMTYPE m_min[NUMDIMS]; ///< Min dimensions of bounding box
ELEMTYPE m_min[NUMDIMS]; ///< Min dimensions of boundin		ELEMTYPE m_max[NUMDIMS]; ///< Max dimensions of bounding box
g box
ELEMTYPE m_max[NUMDIMS]; ///< Max dimensions of boundin
g box
};		};
/// May be data or may be another subtree		/// May be data or may be another subtree
/// The parents level determines this.		/// The parents level determines this.
/// If the parents level is 0, then this is data		/// If the parents level is 0, then this is data
struct Branch		struct Branch {
{		Rect m_rect; ///< Bounds
Rect m_rect; ///< Bounds		union {
union		Node *m_child; ///< Child node
{		DATATYPE m_data; ///< Data Id or Ptr
Node* m_child; ///< Child node
DATATYPE m_data; ///< Data Id or Ptr
};		};
};		};
/// Node for each branch level		/// Node for each branch level
struct Node		struct Node {
{		bool IsInternalNode()
bool IsInternalNode() { return (m_level > 0); } // N		{
ot a leaf, but a internal node		return (m_level > 0);
bool IsLeaf() { return (m_level == 0); } //		} // Not a leaf, but a internal node
A leaf, contains data		bool IsLeaf() { return (m_level == 0); } // A leaf, contains data
int m_count; ///< Count		int m_count; ///< Count
int m_level; ///< Leaf is zero, others posi		int m_level; ///< Leaf is zero, others positive
tive		Branch m_branch[MAXNODES]; ///< Branch
Branch m_branch[MAXNODES]; ///< Branch
};		};
/// A link list of nodes for reinsertion after a delete operation		/// A link list of nodes for reinsertion after a delete operation
struct ListNode		struct ListNode {
{		ListNode *m_next; ///< Next in list
ListNode* m_next; ///< Next in list		Node *m_node; ///< Node
Node* m_node; ///< Node
};		};
/// Variables for finding a split partition		/// Variables for finding a split partition
struct PartitionVars		struct PartitionVars {
{		int m_partition[MAXNODES + 1];
int m_partition[MAXNODES+1];
int m_total;		int m_total;
int m_minFill;		int m_minFill;
int m_taken[MAXNODES+1];		int m_taken[MAXNODES + 1];
int m_count[2];		int m_count[2];
Rect m_cover[2];		Rect m_cover[2];
ELEMTYPEREAL m_area[2];		ELEMTYPEREAL m_area[2];
Branch m_branchBuf[MAXNODES+1];		Branch m_branchBuf[MAXNODES + 1];
int m_branchCount;		int m_branchCount;
Rect m_coverSplit;		Rect m_coverSplit;
ELEMTYPEREAL m_coverSplitArea;		ELEMTYPEREAL m_coverSplitArea;
};		};
Node* AllocNode();		Node *AllocNode();
void FreeNode(Node* a_node);		void FreeNode(Node *a_node);
void InitNode(Node* a_node);		void InitNode(Node *a_node);
void InitRect(Rect* a_rect);		void InitRect(Rect *a_rect);
bool InsertRectRec(Rect* a_rect, const DATATYPE& a_id, Node* a_node, Node** a_		bool InsertRectRec(Rect *a_rect, const DATATYPE &a_id, Node *a_node,
newNode, int a_level);		Node **a_newNode, int a_level);
bool InsertRect(Rect* a_rect, const DATATYPE& a_id, Node** a_root, int a_level		bool InsertRect(Rect *a_rect, const DATATYPE &a_id, Node **a_root,
);		int a_level);
Rect NodeCover(Node* a_node);		Rect NodeCover(Node *a_node);
bool AddBranch(Branch* a_branch, Node* a_node, Node** a_newNode);		bool AddBranch(Branch *a_branch, Node *a_node, Node **a_newNode);
void DisconnectBranch(Node* a_node, int a_index);		void DisconnectBranch(Node *a_node, int a_index);
int PickBranch(Rect* a_rect, Node* a_node);		int PickBranch(Rect *a_rect, Node *a_node);
Rect CombineRect(Rect* a_rectA, Rect* a_rectB);		Rect CombineRect(Rect *a_rectA, Rect *a_rectB);
void SplitNode(Node* a_node, Branch* a_branch, Node** a_newNode);		void SplitNode(Node *a_node, Branch *a_branch, Node **a_newNode);
ELEMTYPEREAL RectSphericalVolume(Rect* a_rect);		ELEMTYPEREAL RectSphericalVolume(Rect *a_rect);
ELEMTYPEREAL RectVolume(Rect* a_rect);		ELEMTYPEREAL RectVolume(Rect *a_rect);
ELEMTYPEREAL CalcRectVolume(Rect* a_rect);		ELEMTYPEREAL CalcRectVolume(Rect *a_rect);
void GetBranches(Node* a_node, Branch* a_branch, PartitionVars* a_parVars);		void GetBranches(Node *a_node, Branch *a_branch, PartitionVars *a_parVars);
void ChoosePartition(PartitionVars* a_parVars, int a_minFill);		void ChoosePartition(PartitionVars *a_parVars, int a_minFill);
void LoadNodes(Node* a_nodeA, Node* a_nodeB, PartitionVars* a_parVars);		void LoadNodes(Node *a_nodeA, Node *a_nodeB, PartitionVars *a_parVars);
void InitParVars(PartitionVars* a_parVars, int a_maxRects, int a_minFill);		void InitParVars(PartitionVars *a_parVars, int a_maxRects, int a_minFill);
void PickSeeds(PartitionVars* a_parVars);		void PickSeeds(PartitionVars *a_parVars);
void Classify(int a_index, int a_group, PartitionVars* a_parVars);		void Classify(int a_index, int a_group, PartitionVars *a_parVars);
bool RemoveRect(Rect* a_rect, const DATATYPE& a_id, Node** a_root);		bool RemoveRect(Rect *a_rect, const DATATYPE &a_id, Node **a_root);
bool RemoveRectRec(Rect* a_rect, const DATATYPE& a_id, Node* a_node, ListNode*		bool RemoveRectRec(Rect *a_rect, const DATATYPE &a_id, Node *a_node,
* a_listNode);		ListNode **a_listNode);
ListNode* AllocListNode();		ListNode *AllocListNode();
void FreeListNode(ListNode* a_listNode);		void FreeListNode(ListNode *a_listNode);
bool Overlap(Rect* a_rectA, Rect* a_rectB);		bool Overlap(Rect *a_rectA, Rect *a_rectB);
void ReInsert(Node* a_node, ListNode** a_listNode);		void ReInsert(Node *a_node, ListNode **a_listNode);
bool Search(Node* a_node, Rect* a_rect, int& a_foundCount, bool a_resultCallba		bool Search(Node *a_node, Rect *a_rect, int &a_foundCount,
ck(DATATYPE a_data, void* a_context), void* a_context);		bool a_resultCallback(DATATYPE a_data, void *a_context),
void RemoveAllRec(Node* a_node);		void *a_context);
		void RemoveAllRec(Node *a_node);
void Reset();		void Reset();
void CountRec(Node* a_node, int& a_count);		void CountRec(Node *a_node, int &a_count);
bool SaveRec(Node* a_node, RTFileStream& a_stream);		bool SaveRec(Node *a_node, RTFileStream &a_stream);
bool LoadRec(Node* a_node, RTFileStream& a_stream);		bool LoadRec(Node *a_node, RTFileStream &a_stream);
Node* m_root; ///< Root of tree		Node *m_root; ///< Root of tree
ELEMTYPEREAL m_unitSphereVolume; ///< Unit sphere constant for		ELEMTYPEREAL m_unitSphereVolume; ///< Unit sphere constant for required number
required number of dimensions		///< of dimensions
};		};
// Because there is not stream support, this is a quick and dirty file I/O helpe		// Because there is not stream support, this is a quick and dirty file I/O
r.		// helper. Users will likely replace its usage with a Stream implementation from
// Users will likely replace its usage with a Stream implementation from their f		// their favorite API.
avorite API.		class RTFileStream {
class RTFileStream		FILE *m_file;
{
FILE* m_file;
public:		public:
		RTFileStream() { m_file = NULL; }
RTFileStream()		~RTFileStream() { Close(); }
{
m_file = NULL;
}
~RTFileStream()
{
Close();
}
bool OpenRead(const char* a_fileName)		bool OpenRead(const char *a_fileName)
{		{
m_file = fopen(a_fileName, "rb");		m_file = fopen(a_fileName, "rb");
if(!m_file)		if(!m_file) { return false; }
{
return false;
}
return true;		return true;
}		}
bool OpenWrite(const char* a_fileName)		bool OpenWrite(const char *a_fileName)
{		{
m_file = fopen(a_fileName, "wb");		m_file = fopen(a_fileName, "wb");
if(!m_file)		if(!m_file) { return false; }
{
return false;
}
return true;		return true;
}		}
void Close()		void Close()
{		{
if(m_file)		if(m_file) {
{
fclose(m_file);		fclose(m_file);
m_file = NULL;		m_file = NULL;
}		}
}		}
template< typename TYPE >		template <typename TYPE> size_t Write(const TYPE &a_value)
size_t Write(const TYPE& a_value)
{		{
ASSERT(m_file);		ASSERT(m_file);
return fwrite((void*)&a_value, sizeof(a_value), 1, m_file);		return fwrite((void *)&a_value, sizeof(a_value), 1, m_file);
}		}
template< typename TYPE >		template <typename TYPE> size_t WriteArray(const TYPE *a_array, int a_count)
size_t WriteArray(const TYPE* a_array, int a_count)
{		{
ASSERT(m_file);		ASSERT(m_file);
return fwrite((void*)a_array, sizeof(TYPE) * a_count, 1, m_file);		return fwrite((void *)a_array, sizeof(TYPE) * a_count, 1, m_file);
}		}
template< typename TYPE >		template <typename TYPE> size_t Read(TYPE &a_value)
size_t Read(TYPE& a_value)
{		{
ASSERT(m_file);		ASSERT(m_file);
return fread((void*)&a_value, sizeof(a_value), 1, m_file);		return fread((void *)&a_value, sizeof(a_value), 1, m_file);
}		}
template< typename TYPE >		template <typename TYPE> size_t ReadArray(TYPE *a_array, int a_count)
size_t ReadArray(TYPE* a_array, int a_count)
{		{
ASSERT(m_file);		ASSERT(m_file);
return fread((void*)a_array, sizeof(TYPE) * a_count, 1, m_file);		return fread((void *)a_array, sizeof(TYPE) * a_count, 1, m_file);
}		}
};		};
RTREE_TEMPLATE		RTREE_TEMPLATE
RTREE_QUAL::RTree()		RTREE_QUAL::RTree()
{		{
ASSERT(MAXNODES > MINNODES);		ASSERT(MAXNODES > MINNODES);
ASSERT(MINNODES > 0);		ASSERT(MINNODES > 0);
// We only support machine word size simple data type eg. integer index or obj		// We only support machine word size simple data type eg. integer index or
ect pointer.		// object pointer. Since we are storing as union with non data branch
// Since we are storing as union with non data branch		ASSERT(sizeof(DATATYPE) == sizeof(void *) || sizeof(DATATYPE) == sizeof(int));
ASSERT(sizeof(DATATYPE) == sizeof(void*) || sizeof(DATATYPE) == sizeof(int));
// Precomputed volumes of the unit spheres for the first few dimensions		// Precomputed volumes of the unit spheres for the first few dimensions
const float UNIT_SPHERE_VOLUMES[] = {		const float UNIT_SPHERE_VOLUMES[] = {
0.000000f, 2.000000f, 3.141593f, // Dimension 0,1,2		0.000000f, 2.000000f, 3.141593f, // Dimension 0,1,2
4.188790f, 4.934802f, 5.263789f, // Dimension 3,4,5		4.188790f, 4.934802f, 5.263789f, // Dimension 3,4,5
5.167713f, 4.724766f, 4.058712f, // Dimension 6,7,8		5.167713f, 4.724766f, 4.058712f, // Dimension 6,7,8
3.298509f, 2.550164f, 1.884104f, // Dimension 9,10,11		3.298509f, 2.550164f, 1.884104f, // Dimension 9,10,11
1.335263f, 0.910629f, 0.599265f, // Dimension 12,13,14		1.335263f, 0.910629f, 0.599265f, // Dimension 12,13,14
0.381443f, 0.235331f, 0.140981f, // Dimension 15,16,17		0.381443f, 0.235331f, 0.140981f, // Dimension 15,16,17
0.082146f, 0.046622f, 0.025807f, // Dimension 18,19,20		0.082146f, 0.046622f, 0.025807f, // Dimension 18,19,20
skipping to change at line 521		skipping to change at line 514
m_unitSphereVolume = (ELEMTYPEREAL)UNIT_SPHERE_VOLUMES[NUMDIMS];		m_unitSphereVolume = (ELEMTYPEREAL)UNIT_SPHERE_VOLUMES[NUMDIMS];
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
RTREE_QUAL::~RTree()		RTREE_QUAL::~RTree()
{		{
Reset(); // Free, or reset node memory		Reset(); // Free, or reset node memory
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::Insert(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMD		void RTREE_QUAL::Insert(const ELEMTYPE a_min[NUMDIMS],
IMS], const DATATYPE& a_dataId)		const ELEMTYPE a_max[NUMDIMS], const DATATYPE &a_dataId)
{		{
#ifdef _DEBUG		#ifdef _DEBUG
for(int index=0; index<NUMDIMS; ++index)		for(int index = 0; index < NUMDIMS; ++index) {
{
ASSERT(a_min[index] <= a_max[index]);		ASSERT(a_min[index] <= a_max[index]);
}		}
#endif //_DEBUG		#endif //_DEBUG
Rect rect;		Rect rect;
for(int axis=0; axis<NUMDIMS; ++axis)		for(int axis = 0; axis < NUMDIMS; ++axis) {
{
rect.m_min[axis] = a_min[axis];		rect.m_min[axis] = a_min[axis];
rect.m_max[axis] = a_max[axis];		rect.m_max[axis] = a_max[axis];
}		}
InsertRect(&rect, a_dataId, &m_root, 0);		InsertRect(&rect, a_dataId, &m_root, 0);
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::Remove(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMD		void RTREE_QUAL::Remove(const ELEMTYPE a_min[NUMDIMS],
IMS], const DATATYPE& a_dataId)		const ELEMTYPE a_max[NUMDIMS], const DATATYPE &a_dataId)
{		{
#ifdef _DEBUG		#ifdef _DEBUG
for(int index=0; index<NUMDIMS; ++index)		for(int index = 0; index < NUMDIMS; ++index) {
{
ASSERT(a_min[index] <= a_max[index]);		ASSERT(a_min[index] <= a_max[index]);
}		}
#endif //_DEBUG		#endif //_DEBUG
Rect rect;		Rect rect;
for(int axis=0; axis<NUMDIMS; ++axis)		for(int axis = 0; axis < NUMDIMS; ++axis) {
{
rect.m_min[axis] = a_min[axis];		rect.m_min[axis] = a_min[axis];
rect.m_max[axis] = a_max[axis];		rect.m_max[axis] = a_max[axis];
}		}
RemoveRect(&rect, a_dataId, &m_root);		RemoveRect(&rect, a_dataId, &m_root);
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
int RTREE_QUAL::Search(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDI		int RTREE_QUAL::Search(const ELEMTYPE a_min[NUMDIMS],
MS], bool a_resultCallback(DATATYPE a_data, void* a_context), void* a_context)		const ELEMTYPE a_max[NUMDIMS],
		bool a_resultCallback(DATATYPE a_data, void *a_context),
		void *a_context)
{		{
#ifdef _DEBUG		#ifdef _DEBUG
for(int index=0; index<NUMDIMS; ++index)		for(int index = 0; index < NUMDIMS; ++index) {
{
ASSERT(a_min[index] <= a_max[index]);		ASSERT(a_min[index] <= a_max[index]);
}		}
#endif //_DEBUG		#endif //_DEBUG
Rect rect;		Rect rect;
for(int axis=0; axis<NUMDIMS; ++axis)		for(int axis = 0; axis < NUMDIMS; ++axis) {
{
rect.m_min[axis] = a_min[axis];		rect.m_min[axis] = a_min[axis];
rect.m_max[axis] = a_max[axis];		rect.m_max[axis] = a_max[axis];
}		}
// NOTE: May want to return search result another way, perhaps returning the n		// NOTE: May want to return search result another way, perhaps returning the
umber of found elements here.		// number of found elements here.
int foundCount = 0;		int foundCount = 0;
Search(m_root, &rect, foundCount, a_resultCallback, a_context);		Search(m_root, &rect, foundCount, a_resultCallback, a_context);
return foundCount;		return foundCount;
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
int RTREE_QUAL::Count()		int RTREE_QUAL::Count()
{		{
int count = 0;		int count = 0;
CountRec(m_root, count);		CountRec(m_root, count);
return count;		return count;
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::CountRec(Node* a_node, int& a_count)		void RTREE_QUAL::CountRec(Node *a_node, int &a_count)
{		{
if(a_node->IsInternalNode()) // not a leaf node		if(a_node->IsInternalNode()) // not a leaf node
{		{
for(int index = 0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{
CountRec(a_node->m_branch[index].m_child, a_count);		CountRec(a_node->m_branch[index].m_child, a_count);
}		}
}		}
else // A leaf node		else // A leaf node
{		{
a_count += a_node->m_count;		a_count += a_node->m_count;
}		}
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::Load(const char* a_fileName)		bool RTREE_QUAL::Load(const char *a_fileName)
{		{
RemoveAll(); // Clear existing tree		RemoveAll(); // Clear existing tree
RTFileStream stream;		RTFileStream stream;
if(!stream.OpenRead(a_fileName))		if(!stream.OpenRead(a_fileName)) { return false; }
{
return false;
}
bool result = Load(stream);		bool result = Load(stream);
stream.Close();		stream.Close();
return result;		return result;
};		};
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::Load(RTFileStream& a_stream)		bool RTREE_QUAL::Load(RTFileStream &a_stream)
{		{
// Write some kind of header		// Write some kind of header
int _dataFileId = ('R'<<0)|('T'<<8)|('R'<<16)|('E'<<24);		int _dataFileId = ('R' << 0) | ('T' << 8) | ('R' << 16) | ('E' << 24);
int _dataSize = sizeof(DATATYPE);		int _dataSize = sizeof(DATATYPE);
int _dataNumDims = NUMDIMS;		int _dataNumDims = NUMDIMS;
int _dataElemSize = sizeof(ELEMTYPE);		int _dataElemSize = sizeof(ELEMTYPE);
int _dataElemRealSize = sizeof(ELEMTYPEREAL);		int _dataElemRealSize = sizeof(ELEMTYPEREAL);
int _dataMaxNodes = TMAXNODES;		int _dataMaxNodes = TMAXNODES;
int _dataMinNodes = TMINNODES;		int _dataMinNodes = TMINNODES;
int dataFileId = 0;		int dataFileId = 0;
int dataSize = 0;		int dataSize = 0;
int dataNumDims = 0;		int dataNumDims = 0;
skipping to change at line 662		skipping to change at line 651
a_stream.Read(dataSize);		a_stream.Read(dataSize);
a_stream.Read(dataNumDims);		a_stream.Read(dataNumDims);
a_stream.Read(dataElemSize);		a_stream.Read(dataElemSize);
a_stream.Read(dataElemRealSize);		a_stream.Read(dataElemRealSize);
a_stream.Read(dataMaxNodes);		a_stream.Read(dataMaxNodes);
a_stream.Read(dataMinNodes);		a_stream.Read(dataMinNodes);
bool result = false;		bool result = false;
// Test if header was valid and compatible		// Test if header was valid and compatible
if( (dataFileId == _dataFileId)		if((dataFileId == _dataFileId) && (dataSize == _dataSize) &&
&& (dataSize == _dataSize)		(dataNumDims == _dataNumDims) && (dataElemSize == _dataElemSize) &&
&& (dataNumDims == _dataNumDims)		(dataElemRealSize == _dataElemRealSize) &&
&& (dataElemSize == _dataElemSize)		(dataMaxNodes == _dataMaxNodes) && (dataMinNodes == _dataMinNodes)) {
&& (dataElemRealSize == _dataElemRealSize)
&& (dataMaxNodes == _dataMaxNodes)
&& (dataMinNodes == _dataMinNodes)
)
{
// Recursively load tree		// Recursively load tree
result = LoadRec(m_root, a_stream);		result = LoadRec(m_root, a_stream);
}		}
return result;		return result;
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::LoadRec(Node* a_node, RTFileStream& a_stream)		bool RTREE_QUAL::LoadRec(Node *a_node, RTFileStream &a_stream)
{		{
a_stream.Read(a_node->m_level);		a_stream.Read(a_node->m_level);
a_stream.Read(a_node->m_count);		a_stream.Read(a_node->m_count);
if(a_node->IsInternalNode()) // not a leaf node		if(a_node->IsInternalNode()) // not a leaf node
{		{
for(int index = 0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{		Branch *curBranch = &a_node->m_branch[index];
Branch* curBranch = &a_node->m_branch[index];
a_stream.ReadArray(curBranch->m_rect.m_min, NUMDIMS);		a_stream.ReadArray(curBranch->m_rect.m_min, NUMDIMS);
a_stream.ReadArray(curBranch->m_rect.m_max, NUMDIMS);		a_stream.ReadArray(curBranch->m_rect.m_max, NUMDIMS);
curBranch->m_child = AllocNode();		curBranch->m_child = AllocNode();
LoadRec(curBranch->m_child, a_stream);		LoadRec(curBranch->m_child, a_stream);
}		}
}		}
else // A leaf node		else // A leaf node
{		{
for(int index = 0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{		Branch *curBranch = &a_node->m_branch[index];
Branch* curBranch = &a_node->m_branch[index];
a_stream.ReadArray(curBranch->m_rect.m_min, NUMDIMS);		a_stream.ReadArray(curBranch->m_rect.m_min, NUMDIMS);
a_stream.ReadArray(curBranch->m_rect.m_max, NUMDIMS);		a_stream.ReadArray(curBranch->m_rect.m_max, NUMDIMS);
a_stream.Read(curBranch->m_data);		a_stream.Read(curBranch->m_data);
}		}
}		}
return true; // Should do more error checking on I/O operations		return true; // Should do more error checking on I/O operations
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::Save(const char* a_fileName)		bool RTREE_QUAL::Save(const char *a_fileName)
{		{
RTFileStream stream;		RTFileStream stream;
if(!stream.OpenWrite(a_fileName))		if(!stream.OpenWrite(a_fileName)) { return false; }
{
return false;
}
bool result = Save(stream);		bool result = Save(stream);
stream.Close();		stream.Close();
return result;		return result;
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::Save(RTFileStream& a_stream)		bool RTREE_QUAL::Save(RTFileStream &a_stream)
{		{
// Write some kind of header		// Write some kind of header
int dataFileId = ('R'<<0)|('T'<<8)|('R'<<16)|('E'<<24);		int dataFileId = ('R' << 0) | ('T' << 8) | ('R' << 16) | ('E' << 24);
int dataSize = sizeof(DATATYPE);		int dataSize = sizeof(DATATYPE);
int dataNumDims = NUMDIMS;		int dataNumDims = NUMDIMS;
int dataElemSize = sizeof(ELEMTYPE);		int dataElemSize = sizeof(ELEMTYPE);
int dataElemRealSize = sizeof(ELEMTYPEREAL);		int dataElemRealSize = sizeof(ELEMTYPEREAL);
int dataMaxNodes = TMAXNODES;		int dataMaxNodes = TMAXNODES;
int dataMinNodes = TMINNODES;		int dataMinNodes = TMINNODES;
a_stream.Write(dataFileId);		a_stream.Write(dataFileId);
a_stream.Write(dataSize);		a_stream.Write(dataSize);
a_stream.Write(dataNumDims);		a_stream.Write(dataNumDims);
skipping to change at line 756		skipping to change at line 735
a_stream.Write(dataMaxNodes);		a_stream.Write(dataMaxNodes);
a_stream.Write(dataMinNodes);		a_stream.Write(dataMinNodes);
// Recursively save tree		// Recursively save tree
bool result = SaveRec(m_root, a_stream);		bool result = SaveRec(m_root, a_stream);
return result;		return result;
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::SaveRec(Node* a_node, RTFileStream& a_stream)		bool RTREE_QUAL::SaveRec(Node *a_node, RTFileStream &a_stream)
{		{
a_stream.Write(a_node->m_level);		a_stream.Write(a_node->m_level);
a_stream.Write(a_node->m_count);		a_stream.Write(a_node->m_count);
if(a_node->IsInternalNode()) // not a leaf node		if(a_node->IsInternalNode()) // not a leaf node
{		{
for(int index = 0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{		Branch *curBranch = &a_node->m_branch[index];
Branch* curBranch = &a_node->m_branch[index];
a_stream.WriteArray(curBranch->m_rect.m_min, NUMDIMS);		a_stream.WriteArray(curBranch->m_rect.m_min, NUMDIMS);
a_stream.WriteArray(curBranch->m_rect.m_max, NUMDIMS);		a_stream.WriteArray(curBranch->m_rect.m_max, NUMDIMS);
SaveRec(curBranch->m_child, a_stream);		SaveRec(curBranch->m_child, a_stream);
}		}
}		}
else // A leaf node		else // A leaf node
{		{
for(int index = 0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{		Branch *curBranch = &a_node->m_branch[index];
Branch* curBranch = &a_node->m_branch[index];
a_stream.WriteArray(curBranch->m_rect.m_min, NUMDIMS);		a_stream.WriteArray(curBranch->m_rect.m_min, NUMDIMS);
a_stream.WriteArray(curBranch->m_rect.m_max, NUMDIMS);		a_stream.WriteArray(curBranch->m_rect.m_max, NUMDIMS);
a_stream.Write(curBranch->m_data);		a_stream.Write(curBranch->m_data);
}		}
}		}
return true; // Should do more error checking on I/O operations		return true; // Should do more error checking on I/O operations
}		}
skipping to change at line 812		skipping to change at line 789
#ifdef RTREE_DONT_USE_MEMPOOLS		#ifdef RTREE_DONT_USE_MEMPOOLS
// Delete all existing nodes		// Delete all existing nodes
RemoveAllRec(m_root);		RemoveAllRec(m_root);
#else // RTREE_DONT_USE_MEMPOOLS		#else // RTREE_DONT_USE_MEMPOOLS
// Just reset memory pools. We are not using complex types		// Just reset memory pools. We are not using complex types
// EXAMPLE		// EXAMPLE
#endif // RTREE_DONT_USE_MEMPOOLS		#endif // RTREE_DONT_USE_MEMPOOLS
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::RemoveAllRec(Node* a_node)		void RTREE_QUAL::RemoveAllRec(Node *a_node)
{		{
ASSERT(a_node);		ASSERT(a_node);
ASSERT(a_node->m_level >= 0);		ASSERT(a_node->m_level >= 0);
if(a_node->IsInternalNode()) // This is an internal node in the tree		if(a_node->IsInternalNode()) // This is an internal node in the tree
{		{
for(int index=0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{
RemoveAllRec(a_node->m_branch[index].m_child);		RemoveAllRec(a_node->m_branch[index].m_child);
}		}
}		}
FreeNode(a_node);		FreeNode(a_node);
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
typename RTREE_QUAL::Node* RTREE_QUAL::AllocNode()		typename RTREE_QUAL::Node *RTREE_QUAL::AllocNode()
{		{
Node* newNode;		Node *newNode;
#ifdef RTREE_DONT_USE_MEMPOOLS		#ifdef RTREE_DONT_USE_MEMPOOLS
newNode = new Node;		newNode = new Node;
#else // RTREE_DONT_USE_MEMPOOLS		#else // RTREE_DONT_USE_MEMPOOLS
// EXAMPLE		// EXAMPLE
#endif // RTREE_DONT_USE_MEMPOOLS		#endif // RTREE_DONT_USE_MEMPOOLS
InitNode(newNode);		InitNode(newNode);
return newNode;		return newNode;
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::FreeNode(Node* a_node)		void RTREE_QUAL::FreeNode(Node *a_node)
{		{
ASSERT(a_node);		ASSERT(a_node);
#ifdef RTREE_DONT_USE_MEMPOOLS		#ifdef RTREE_DONT_USE_MEMPOOLS
delete a_node;		delete a_node;
#else // RTREE_DONT_USE_MEMPOOLS		#else // RTREE_DONT_USE_MEMPOOLS
// EXAMPLE		// EXAMPLE
#endif // RTREE_DONT_USE_MEMPOOLS		#endif // RTREE_DONT_USE_MEMPOOLS
}		}
// Allocate space for a node in the list used in DeletRect to		// Allocate space for a node in the list used in DeletRect to
// store Nodes that are too empty.		// store Nodes that are too empty.
RTREE_TEMPLATE		RTREE_TEMPLATE
typename RTREE_QUAL::ListNode* RTREE_QUAL::AllocListNode()		typename RTREE_QUAL::ListNode *RTREE_QUAL::AllocListNode()
{		{
#ifdef RTREE_DONT_USE_MEMPOOLS		#ifdef RTREE_DONT_USE_MEMPOOLS
return new ListNode;		return new ListNode;
#else // RTREE_DONT_USE_MEMPOOLS		#else // RTREE_DONT_USE_MEMPOOLS
// EXAMPLE		// EXAMPLE
#endif // RTREE_DONT_USE_MEMPOOLS		#endif // RTREE_DONT_USE_MEMPOOLS
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::FreeListNode(ListNode* a_listNode)		void RTREE_QUAL::FreeListNode(ListNode *a_listNode)
{		{
#ifdef RTREE_DONT_USE_MEMPOOLS		#ifdef RTREE_DONT_USE_MEMPOOLS
delete a_listNode;		delete a_listNode;
#else // RTREE_DONT_USE_MEMPOOLS		#else // RTREE_DONT_USE_MEMPOOLS
// EXAMPLE		// EXAMPLE
#endif // RTREE_DONT_USE_MEMPOOLS		#endif // RTREE_DONT_USE_MEMPOOLS
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::InitNode(Node* a_node)		void RTREE_QUAL::InitNode(Node *a_node)
{		{
a_node->m_count = 0;		a_node->m_count = 0;
a_node->m_level = -1;		a_node->m_level = -1;
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::InitRect(Rect* a_rect)		void RTREE_QUAL::InitRect(Rect *a_rect)
{		{
for(int index = 0; index < NUMDIMS; ++index)		for(int index = 0; index < NUMDIMS; ++index) {
{
a_rect->m_min[index] = (ELEMTYPE)0;		a_rect->m_min[index] = (ELEMTYPE)0;
a_rect->m_max[index] = (ELEMTYPE)0;		a_rect->m_max[index] = (ELEMTYPE)0;
}		}
}		}
// Inserts a new data rectangle into the index structure.		// Inserts a new data rectangle into the index structure.
// Recursively descends tree, propagates splits back up.		// Recursively descends tree, propagates splits back up.
// Returns 0 if node was not split. Old node updated.		// Returns 0 if node was not split. Old node updated.
// If node was split, returns 1 and sets the pointer pointed to by		// If node was split, returns 1 and sets the pointer pointed to by
// new_node to point to the new node. Old node updated to become one of two.		// new_node to point to the new node. Old node updated to become one of two.
// The level argument specifies the number of steps up from the leaf		// The level argument specifies the number of steps up from the leaf
// level to insert; e.g. a data rectangle goes in at level = 0.		// level to insert; e.g. a data rectangle goes in at level = 0.
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::InsertRectRec(Rect* a_rect, const DATATYPE& a_id, Node* a_node,		bool RTREE_QUAL::InsertRectRec(Rect *a_rect, const DATATYPE &a_id, Node *a_node,
Node** a_newNode, int a_level)		Node **a_newNode, int a_level)
{		{
ASSERT(a_rect && a_node && a_newNode);		ASSERT(a_rect && a_node && a_newNode);
ASSERT(a_level >= 0 && a_level <= a_node->m_level);		ASSERT(a_level >= 0 && a_level <= a_node->m_level);
int index;		int index;
Branch branch;		Branch branch;
Node* otherNode;		Node *otherNode;
// Still above level for insertion, go down tree recursively		// Still above level for insertion, go down tree recursively
if(a_node->m_level > a_level)		if(a_node->m_level > a_level) {
{
index = PickBranch(a_rect, a_node);		index = PickBranch(a_rect, a_node);
if(index < 0) return false; // Added for Gmsh		if(index < 0) return false; // Added for Gmsh
if (!InsertRectRec(a_rect, a_id, a_node->m_branch[index].m_child, &otherNode		if(!InsertRectRec(a_rect, a_id, a_node->m_branch[index].m_child, &otherNode,
, a_level))		a_level)) {
{
// Child was not split		// Child was not split
a_node->m_branch[index].m_rect = CombineRect(a_rect, &(a_node->m_branch[in		a_node->m_branch[index].m_rect =
dex].m_rect));		CombineRect(a_rect, &(a_node->m_branch[index].m_rect));
return false;		return false;
}		}
else // Child was split		else // Child was split
{		{
a_node->m_branch[index].m_rect = NodeCover(a_node->m_branch[index].m_child		a_node->m_branch[index].m_rect =
);		NodeCover(a_node->m_branch[index].m_child);
branch.m_child = otherNode;		branch.m_child = otherNode;
branch.m_rect = NodeCover(otherNode);		branch.m_rect = NodeCover(otherNode);
return AddBranch(&branch, a_node, a_newNode);		return AddBranch(&branch, a_node, a_newNode);
}		}
}		}
else if(a_node->m_level == a_level) // Have reached level for insertion. Add r		else if(a_node->m_level == a_level) // Have reached level for insertion. Add
ect, split if necessary		// rect, split if necessary
{		{
branch.m_rect = *a_rect;		branch.m_rect = *a_rect;
branch.m_child = (Node*) a_id;		branch.m_child = (Node *)a_id;
// Child field of leaves contains id of data record		// Child field of leaves contains id of data record
return AddBranch(&branch, a_node, a_newNode);		return AddBranch(&branch, a_node, a_newNode);
}		}
else		else {
{
// Should never occur		// Should never occur
ASSERT(0);		ASSERT(0);
return false;		return false;
}		}
}		}
// Insert a data rectangle into an index structure.		// Insert a data rectangle into an index structure.
// InsertRect provides for splitting the root;		// InsertRect provides for splitting the root;
// returns 1 if root was split, 0 if it was not.		// returns 1 if root was split, 0 if it was not.
// The level argument specifies the number of steps up from the leaf		// The level argument specifies the number of steps up from the leaf
// level to insert; e.g. a data rectangle goes in at level = 0.		// level to insert; e.g. a data rectangle goes in at level = 0.
// InsertRect2 does the recursion.		// InsertRect2 does the recursion.
//		//
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::InsertRect(Rect* a_rect, const DATATYPE& a_id, Node** a_root, i		bool RTREE_QUAL::InsertRect(Rect *a_rect, const DATATYPE &a_id, Node **a_root,
nt a_level)		int a_level)
{		{
ASSERT(a_rect && a_root);		ASSERT(a_rect && a_root);
ASSERT(a_level >= 0 && a_level <= (*a_root)->m_level);		ASSERT(a_level >= 0 && a_level <= (*a_root)->m_level);
#ifdef _DEBUG		#ifdef _DEBUG
for(int index=0; index < NUMDIMS; ++index)		for(int index = 0; index < NUMDIMS; ++index) {
{
ASSERT(a_rect->m_min[index] <= a_rect->m_max[index]);		ASSERT(a_rect->m_min[index] <= a_rect->m_max[index]);
}		}
#endif //_DEBUG		#endif //_DEBUG
Node* newRoot;		Node *newRoot;
Node* newNode;		Node *newNode;
Branch branch;		Branch branch;
if(InsertRectRec(a_rect, a_id, *a_root, &newNode, a_level)) // Root split		if(InsertRectRec(a_rect, a_id, *a_root, &newNode, a_level)) // Root split
{		{
newRoot = AllocNode(); // Grow tree taller and new root		newRoot = AllocNode(); // Grow tree taller and new root
newRoot->m_level = (*a_root)->m_level + 1;		newRoot->m_level = (*a_root)->m_level + 1;
branch.m_rect = NodeCover(*a_root);		branch.m_rect = NodeCover(*a_root);
branch.m_child = *a_root;		branch.m_child = *a_root;
AddBranch(&branch, newRoot, NULL);		AddBranch(&branch, newRoot, NULL);
branch.m_rect = NodeCover(newNode);		branch.m_rect = NodeCover(newNode);
branch.m_child = newNode;		branch.m_child = newNode;
AddBranch(&branch, newRoot, NULL);		AddBranch(&branch, newRoot, NULL);
*a_root = newRoot;		*a_root = newRoot;
return true;		return true;
}		}
return false;		return false;
}		}
// Find the smallest rectangle that includes all rectangles in branches of a nod		// Find the smallest rectangle that includes all rectangles in branches of a
e.		// node.
RTREE_TEMPLATE		RTREE_TEMPLATE
typename RTREE_QUAL::Rect RTREE_QUAL::NodeCover(Node* a_node)		typename RTREE_QUAL::Rect RTREE_QUAL::NodeCover(Node *a_node)
{		{
ASSERT(a_node);		ASSERT(a_node);
int firstTime = true;		int firstTime = true;
Rect rect;		Rect rect;
InitRect(&rect);		InitRect(&rect);
for(int index = 0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{		if(firstTime) {
if(firstTime)
{
rect = a_node->m_branch[index].m_rect;		rect = a_node->m_branch[index].m_rect;
firstTime = false;		firstTime = false;
}		}
else		else {
{
rect = CombineRect(&rect, &(a_node->m_branch[index].m_rect));		rect = CombineRect(&rect, &(a_node->m_branch[index].m_rect));
}		}
}		}
return rect;		return rect;
}		}
// Add a branch to a node. Split the node if necessary.		// Add a branch to a node. Split the node if necessary.
// Returns 0 if node not split. Old node updated.		// Returns 0 if node not split. Old node updated.
// Returns 1 if node split, sets *new_node to address of new node.		// Returns 1 if node split, sets *new_node to address of new node.
// Old node updated, becomes one of two.		// Old node updated, becomes one of two.
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::AddBranch(Branch* a_branch, Node* a_node, Node** a_newNode)		bool RTREE_QUAL::AddBranch(Branch *a_branch, Node *a_node, Node **a_newNode)
{		{
ASSERT(a_branch);		ASSERT(a_branch);
ASSERT(a_node);		ASSERT(a_node);
if(a_node->m_count < MAXNODES) // Split won't be necessary		if(a_node->m_count < MAXNODES) // Split won't be necessary
{		{
a_node->m_branch[a_node->m_count] = *a_branch;		a_node->m_branch[a_node->m_count] = *a_branch;
++a_node->m_count;		++a_node->m_count;
return false;		return false;
}		}
else		else {
{
ASSERT(a_newNode);		ASSERT(a_newNode);
SplitNode(a_node, a_branch, a_newNode);		SplitNode(a_node, a_branch, a_newNode);
return true;		return true;
}		}
}		}
// Disconnect a dependent node.		// Disconnect a dependent node.
// Caller must return (or stop using iteration index) after this as count has ch		// Caller must return (or stop using iteration index) after this as count has
anged		// changed
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::DisconnectBranch(Node* a_node, int a_index)		void RTREE_QUAL::DisconnectBranch(Node *a_node, int a_index)
{		{
ASSERT(a_node && (a_index >= 0) && (a_index < MAXNODES));		ASSERT(a_node && (a_index >= 0) && (a_index < MAXNODES));
ASSERT(a_node->m_count > 0);		ASSERT(a_node->m_count > 0);
// Remove element by swapping with the last element to prevent gaps in array		// Remove element by swapping with the last element to prevent gaps in array
a_node->m_branch[a_index] = a_node->m_branch[a_node->m_count - 1];		a_node->m_branch[a_index] = a_node->m_branch[a_node->m_count - 1];
--a_node->m_count;		--a_node->m_count;
}		}
// Pick a branch. Pick the one that will need the smallest increase		// Pick a branch. Pick the one that will need the smallest increase
// in area to accomodate the new rectangle. This will result in the		// in area to accomodate the new rectangle. This will result in the
// least total area for the covering rectangles in the current node.		// least total area for the covering rectangles in the current node.
// In case of a tie, pick the one which was smaller before, to get		// In case of a tie, pick the one which was smaller before, to get
// the best resolution when searching.		// the best resolution when searching.
RTREE_TEMPLATE		RTREE_TEMPLATE
int RTREE_QUAL::PickBranch(Rect* a_rect, Node* a_node)		int RTREE_QUAL::PickBranch(Rect *a_rect, Node *a_node)
{		{
ASSERT(a_rect && a_node);		ASSERT(a_rect && a_node);
bool firstTime = true;		bool firstTime = true;
ELEMTYPEREAL increase;		ELEMTYPEREAL increase;
ELEMTYPEREAL bestIncr = (ELEMTYPEREAL)-1;		ELEMTYPEREAL bestIncr = (ELEMTYPEREAL)-1;
ELEMTYPEREAL area;		ELEMTYPEREAL area;
ELEMTYPEREAL bestArea = (ELEMTYPEREAL)-1;		ELEMTYPEREAL bestArea = (ELEMTYPEREAL)-1;
int best = -1;		int best = -1;
Rect tempRect;		Rect tempRect;
for(int index=0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{		Rect *curRect = &a_node->m_branch[index].m_rect;
Rect* curRect = &a_node->m_branch[index].m_rect;
area = CalcRectVolume(curRect);		area = CalcRectVolume(curRect);
tempRect = CombineRect(a_rect, curRect);		tempRect = CombineRect(a_rect, curRect);
increase = CalcRectVolume(&tempRect) - area;		increase = CalcRectVolume(&tempRect) - area;
if((increase < bestIncr) || firstTime)		if((increase < bestIncr) || firstTime) {
{
best = index;		best = index;
bestArea = area;		bestArea = area;
bestIncr = increase;		bestIncr = increase;
firstTime = false;		firstTime = false;
}		}
else if((increase == bestIncr) && (area < bestArea))		else if((increase == bestIncr) && (area < bestArea)) {
{
best = index;		best = index;
bestArea = area;		bestArea = area;
bestIncr = increase;		bestIncr = increase;
}		}
}		}
return best;		return best;
}		}
// Combine two rectangles into larger one containing both		// Combine two rectangles into larger one containing both
RTREE_TEMPLATE		RTREE_TEMPLATE
typename RTREE_QUAL::Rect RTREE_QUAL::CombineRect(Rect* a_rectA, Rect* a_rectB)		typename RTREE_QUAL::Rect RTREE_QUAL::CombineRect(Rect *a_rectA, Rect *a_rectB)
{		{
ASSERT(a_rectA && a_rectB);		ASSERT(a_rectA && a_rectB);
Rect newRect;		Rect newRect;
for(int index = 0; index < NUMDIMS; ++index)		for(int index = 0; index < NUMDIMS; ++index) {
{		newRect.m_min[index] =
newRect.m_min[index] = std::min(a_rectA->m_min[index], a_rectB->m_min[index]		std::min(a_rectA->m_min[index], a_rectB->m_min[index]);
);		newRect.m_max[index] =
newRect.m_max[index] = std::max(a_rectA->m_max[index], a_rectB->m_max[index]		std::max(a_rectA->m_max[index], a_rectB->m_max[index]);
);
}		}
return newRect;		return newRect;
}		}
// Split a node.		// Split a node.
// Divides the nodes branches and the extra one between two nodes.		// Divides the nodes branches and the extra one between two nodes.
// Old node is one of the new ones, and one really new one is created.		// Old node is one of the new ones, and one really new one is created.
// Tries more than one method for choosing a partition, uses best result.		// Tries more than one method for choosing a partition, uses best result.
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::SplitNode(Node* a_node, Branch* a_branch, Node** a_newNode)		void RTREE_QUAL::SplitNode(Node *a_node, Branch *a_branch, Node **a_newNode)
{		{
ASSERT(a_node);		ASSERT(a_node);
ASSERT(a_branch);		ASSERT(a_branch);
// Could just use local here, but member or external is faster since it is reu		// Could just use local here, but member or external is faster since it is
sed		// reused
PartitionVars localVars;		PartitionVars localVars;
PartitionVars* parVars = &localVars;		PartitionVars *parVars = &localVars;
int level;		int level;
// Load all the branches into a buffer, initialize old node		// Load all the branches into a buffer, initialize old node
level = a_node->m_level;		level = a_node->m_level;
GetBranches(a_node, a_branch, parVars);		GetBranches(a_node, a_branch, parVars);
// Find partition		// Find partition
ChoosePartition(parVars, MINNODES);		ChoosePartition(parVars, MINNODES);
// Put branches from buffer into 2 nodes according to chosen partition		// Put branches from buffer into 2 nodes according to chosen partition
*a_newNode = AllocNode();		*a_newNode = AllocNode();
(*a_newNode)->m_level = a_node->m_level = level;		(*a_newNode)->m_level = a_node->m_level = level;
LoadNodes(a_node, *a_newNode, parVars);		LoadNodes(a_node, *a_newNode, parVars);
ASSERT((a_node->m_count + (*a_newNode)->m_count) == parVars->m_total);		ASSERT((a_node->m_count + (*a_newNode)->m_count) == parVars->m_total);
}		}
// Calculate the n-dimensional volume of a rectangle		// Calculate the n-dimensional volume of a rectangle
RTREE_TEMPLATE		RTREE_TEMPLATE
ELEMTYPEREAL RTREE_QUAL::RectVolume(Rect* a_rect)		ELEMTYPEREAL RTREE_QUAL::RectVolume(Rect *a_rect)
{		{
ASSERT(a_rect);		ASSERT(a_rect);
ELEMTYPEREAL volume = (ELEMTYPEREAL)1;		ELEMTYPEREAL volume = (ELEMTYPEREAL)1;
for(int index=0; index<NUMDIMS; ++index)		for(int index = 0; index < NUMDIMS; ++index) {
{
volume *= a_rect->m_max[index] - a_rect->m_min[index];		volume *= a_rect->m_max[index] - a_rect->m_min[index];
}		}
ASSERT(volume >= (ELEMTYPEREAL)0);		ASSERT(volume >= (ELEMTYPEREAL)0);
return volume;		return volume;
}		}
// The exact volume of the bounding sphere for the given Rect		// The exact volume of the bounding sphere for the given Rect
RTREE_TEMPLATE		RTREE_TEMPLATE
ELEMTYPEREAL RTREE_QUAL::RectSphericalVolume(Rect* a_rect)		ELEMTYPEREAL RTREE_QUAL::RectSphericalVolume(Rect *a_rect)
{		{
ASSERT(a_rect);		ASSERT(a_rect);
ELEMTYPEREAL sumOfSquares = (ELEMTYPEREAL)0;		ELEMTYPEREAL sumOfSquares = (ELEMTYPEREAL)0;
ELEMTYPEREAL radius;		ELEMTYPEREAL radius;
for(int index=0; index < NUMDIMS; ++index)		for(int index = 0; index < NUMDIMS; ++index) {
{		ELEMTYPEREAL halfExtent = ((ELEMTYPEREAL)a_rect->m_max[index] -
ELEMTYPEREAL halfExtent = ((ELEMTYPEREAL)a_rect->m_max[index] - (ELEMTYPEREA		(ELEMTYPEREAL)a_rect->m_min[index]) *
L)a_rect->m_min[index]) * 0.5f;		0.5f;
sumOfSquares += halfExtent * halfExtent;		sumOfSquares += halfExtent * halfExtent;
}		}
radius = (ELEMTYPEREAL)sqrt(sumOfSquares);		radius = (ELEMTYPEREAL)sqrt(sumOfSquares);
// Pow maybe slow, so test for common dims like 2,3 and just use x*x, x*x*x.		// Pow maybe slow, so test for common dims like 2,3 and just use x*x, x*x*x.
if(NUMDIMS == 3)		if(NUMDIMS == 3) { return (radius * radius * radius * m_unitSphereVolume); }
{		else if(NUMDIMS == 2) {
return (radius * radius * radius * m_unitSphereVolume);
}
else if(NUMDIMS == 2)
{
return (radius * radius * m_unitSphereVolume);		return (radius * radius * m_unitSphereVolume);
}		}
else		else {
{
return (ELEMTYPEREAL)(pow(radius, NUMDIMS) * m_unitSphereVolume);		return (ELEMTYPEREAL)(pow(radius, NUMDIMS) * m_unitSphereVolume);
}		}
}		}
// Use one of the methods to calculate retangle volume		// Use one of the methods to calculate retangle volume
RTREE_TEMPLATE		RTREE_TEMPLATE
ELEMTYPEREAL RTREE_QUAL::CalcRectVolume(Rect* a_rect)		ELEMTYPEREAL RTREE_QUAL::CalcRectVolume(Rect *a_rect)
{		{
#ifdef RTREE_USE_SPHERICAL_VOLUME		#ifdef RTREE_USE_SPHERICAL_VOLUME
return RectSphericalVolume(a_rect); // Slower but helps certain merge cases		return RectSphericalVolume(a_rect); // Slower but helps certain merge cases
#else // RTREE_USE_SPHERICAL_VOLUME		#else // RTREE_USE_SPHERICAL_VOLUME
return RectVolume(a_rect); // Faster but can cause poor merges		return RectVolume(a_rect); // Faster but can cause poor merges
#endif // RTREE_USE_SPHERICAL_VOLUME		#endif // RTREE_USE_SPHERICAL_VOLUME
}		}
// Load branch buffer with branches from full node plus the extra branch.		// Load branch buffer with branches from full node plus the extra branch.
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::GetBranches(Node* a_node, Branch* a_branch, PartitionVars* a_pa		void RTREE_QUAL::GetBranches(Node *a_node, Branch *a_branch,
rVars)		PartitionVars *a_parVars)
{		{
ASSERT(a_node);		ASSERT(a_node);
ASSERT(a_branch);		ASSERT(a_branch);
ASSERT(a_node->m_count == MAXNODES);		ASSERT(a_node->m_count == MAXNODES);
// Load the branch buffer		// Load the branch buffer
for(int index=0; index < MAXNODES; ++index)		for(int index = 0; index < MAXNODES; ++index) {
{
a_parVars->m_branchBuf[index] = a_node->m_branch[index];		a_parVars->m_branchBuf[index] = a_node->m_branch[index];
}		}
a_parVars->m_branchBuf[MAXNODES] = *a_branch;		a_parVars->m_branchBuf[MAXNODES] = *a_branch;
a_parVars->m_branchCount = MAXNODES + 1;		a_parVars->m_branchCount = MAXNODES + 1;
// Calculate rect containing all in the set		// Calculate rect containing all in the set
a_parVars->m_coverSplit = a_parVars->m_branchBuf[0].m_rect;		a_parVars->m_coverSplit = a_parVars->m_branchBuf[0].m_rect;
for(int index=1; index < MAXNODES+1; ++index)		for(int index = 1; index < MAXNODES + 1; ++index) {
{		a_parVars->m_coverSplit = CombineRect(
a_parVars->m_coverSplit = CombineRect(&a_parVars->m_coverSplit, &a_parVars->		&a_parVars->m_coverSplit, &a_parVars->m_branchBuf[index].m_rect);
m_branchBuf[index].m_rect);
}		}
a_parVars->m_coverSplitArea = CalcRectVolume(&a_parVars->m_coverSplit);		a_parVars->m_coverSplitArea = CalcRectVolume(&a_parVars->m_coverSplit);
InitNode(a_node);		InitNode(a_node);
}		}
// Method #0 for choosing a partition:		// Method #0 for choosing a partition:
// As the seeds for the two groups, pick the two rects that would waste the		// As the seeds for the two groups, pick the two rects that would waste the
// most area if covered by a single rectangle, i.e. evidently the worst pair		// most area if covered by a single rectangle, i.e. evidently the worst pair
// to have in the same group.		// to have in the same group.
// Of the remaining, one at a time is chosen to be put in one of the two groups.		// Of the remaining, one at a time is chosen to be put in one of the two groups.
// The one chosen is the one with the greatest difference in area expansion		// The one chosen is the one with the greatest difference in area expansion
// depending on which group - the rect most strongly attracted to one group		// depending on which group - the rect most strongly attracted to one group
// and repelled from the other.		// and repelled from the other.
// If one group gets too full (more would force other group to violate min		// If one group gets too full (more would force other group to violate min
// fill requirement) then other group gets the rest.		// fill requirement) then other group gets the rest.
// These last are the ones that can go in either group most easily.		// These last are the ones that can go in either group most easily.
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::ChoosePartition(PartitionVars* a_parVars, int a_minFill)		void RTREE_QUAL::ChoosePartition(PartitionVars *a_parVars, int a_minFill)
{		{
ASSERT(a_parVars);		ASSERT(a_parVars);
ELEMTYPEREAL biggestDiff;		ELEMTYPEREAL biggestDiff;
int group, chosen = -1, betterGroup = -1;		int group, chosen = -1, betterGroup = -1;
InitParVars(a_parVars, a_parVars->m_branchCount, a_minFill);		InitParVars(a_parVars, a_parVars->m_branchCount, a_minFill);
PickSeeds(a_parVars);		PickSeeds(a_parVars);
while (((a_parVars->m_count[0] + a_parVars->m_count[1]) < a_parVars->m_total)		while(
&& (a_parVars->m_count[0] < (a_parVars->m_total - a_parVars->m_minFill))		((a_parVars->m_count[0] + a_parVars->m_count[1]) < a_parVars->m_total) &&
&& (a_parVars->m_count[1] < (a_parVars->m_total - a_parVars->m_minFill)))		(a_parVars->m_count[0] < (a_parVars->m_total - a_parVars->m_minFill)) &&
{		(a_parVars->m_count[1] < (a_parVars->m_total - a_parVars->m_minFill))) {
biggestDiff = (ELEMTYPEREAL) -1;		biggestDiff = (ELEMTYPEREAL)-1;
for(int index=0; index<a_parVars->m_total; ++index)		for(int index = 0; index < a_parVars->m_total; ++index) {
{		if(!a_parVars->m_taken[index]) {
if(!a_parVars->m_taken[index])		Rect *curRect = &a_parVars->m_branchBuf[index].m_rect;
{
Rect* curRect = &a_parVars->m_branchBuf[index].m_rect;
Rect rect0 = CombineRect(curRect, &a_parVars->m_cover[0]);		Rect rect0 = CombineRect(curRect, &a_parVars->m_cover[0]);
Rect rect1 = CombineRect(curRect, &a_parVars->m_cover[1]);		Rect rect1 = CombineRect(curRect, &a_parVars->m_cover[1]);
ELEMTYPEREAL growth0 = CalcRectVolume(&rect0) - a_parVars->m_area[0];		ELEMTYPEREAL growth0 = CalcRectVolume(&rect0) - a_parVars->m_area[0];
ELEMTYPEREAL growth1 = CalcRectVolume(&rect1) - a_parVars->m_area[1];		ELEMTYPEREAL growth1 = CalcRectVolume(&rect1) - a_parVars->m_area[1];
ELEMTYPEREAL diff = growth1 - growth0;		ELEMTYPEREAL diff = growth1 - growth0;
if(diff >= 0)		if(diff >= 0) { group = 0; }
{		else {
group = 0;
}
else
{
group = 1;		group = 1;
diff = -diff;		diff = -diff;
}		}
if(diff > biggestDiff)		if(diff > biggestDiff) {
{
biggestDiff = diff;		biggestDiff = diff;
chosen = index;		chosen = index;
betterGroup = group;		betterGroup = group;
}		}
else if((diff == biggestDiff) && (a_parVars->m_count[group] < a_parVars-		else if((diff == biggestDiff) &&
>m_count[betterGroup]))		(a_parVars->m_count[group] < a_parVars->m_count[betterGroup])) {
{
chosen = index;		chosen = index;
betterGroup = group;		betterGroup = group;
}		}
}		}
}		}
Classify(chosen, betterGroup, a_parVars);		Classify(chosen, betterGroup, a_parVars);
}		}
// If one group too full, put remaining rects in the other		// If one group too full, put remaining rects in the other
if((a_parVars->m_count[0] + a_parVars->m_count[1]) < a_parVars->m_total)		if((a_parVars->m_count[0] + a_parVars->m_count[1]) < a_parVars->m_total) {
{		if(a_parVars->m_count[0] >= a_parVars->m_total - a_parVars->m_minFill) {
if(a_parVars->m_count[0] >= a_parVars->m_total - a_parVars->m_minFill)
{
group = 1;		group = 1;
}		}
else		else {
{
group = 0;		group = 0;
}		}
for(int index=0; index<a_parVars->m_total; ++index)		for(int index = 0; index < a_parVars->m_total; ++index) {
{		if(!a_parVars->m_taken[index]) { Classify(index, group, a_parVars); }
if(!a_parVars->m_taken[index])
{
Classify(index, group, a_parVars);
}
}		}
}		}
ASSERT((a_parVars->m_count[0] + a_parVars->m_count[1]) == a_parVars->m_total);		ASSERT((a_parVars->m_count[0] + a_parVars->m_count[1]) == a_parVars->m_total);
ASSERT((a_parVars->m_count[0] >= a_parVars->m_minFill) &&		ASSERT((a_parVars->m_count[0] >= a_parVars->m_minFill) &&
(a_parVars->m_count[1] >= a_parVars->m_minFill));		(a_parVars->m_count[1] >= a_parVars->m_minFill));
}		}
// Copy branches from the buffer into two nodes according to the partition.		// Copy branches from the buffer into two nodes according to the partition.
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::LoadNodes(Node* a_nodeA, Node* a_nodeB, PartitionVars* a_parVar		void RTREE_QUAL::LoadNodes(Node *a_nodeA, Node *a_nodeB,
s)		PartitionVars *a_parVars)
{		{
ASSERT(a_nodeA);		ASSERT(a_nodeA);
ASSERT(a_nodeB);		ASSERT(a_nodeB);
ASSERT(a_parVars);		ASSERT(a_parVars);
for(int index=0; index < a_parVars->m_total; ++index)		for(int index = 0; index < a_parVars->m_total; ++index) {
{		ASSERT(a_parVars->m_partition[index] == 0 ||
ASSERT(a_parVars->m_partition[index] == 0 || a_parVars->m_partition[index] =		a_parVars->m_partition[index] == 1);
= 1);
if(a_parVars->m_partition[index] == 0)		if(a_parVars->m_partition[index] == 0) {
{
AddBranch(&a_parVars->m_branchBuf[index], a_nodeA, NULL);		AddBranch(&a_parVars->m_branchBuf[index], a_nodeA, NULL);
}		}
else if(a_parVars->m_partition[index] == 1)		else if(a_parVars->m_partition[index] == 1) {
{
AddBranch(&a_parVars->m_branchBuf[index], a_nodeB, NULL);		AddBranch(&a_parVars->m_branchBuf[index], a_nodeB, NULL);
}		}
}		}
}		}
// Initialize a PartitionVars structure.		// Initialize a PartitionVars structure.
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::InitParVars(PartitionVars* a_parVars, int a_maxRects, int a_min		void RTREE_QUAL::InitParVars(PartitionVars *a_parVars, int a_maxRects,
Fill)		int a_minFill)
{		{
ASSERT(a_parVars);		ASSERT(a_parVars);
a_parVars->m_count[0] = a_parVars->m_count[1] = 0;		a_parVars->m_count[0] = a_parVars->m_count[1] = 0;
a_parVars->m_area[0] = a_parVars->m_area[1] = (ELEMTYPEREAL)0;		a_parVars->m_area[0] = a_parVars->m_area[1] = (ELEMTYPEREAL)0;
a_parVars->m_total = a_maxRects;		a_parVars->m_total = a_maxRects;
a_parVars->m_minFill = a_minFill;		a_parVars->m_minFill = a_minFill;
for(int index=0; index < a_maxRects; ++index)		for(int index = 0; index < a_maxRects; ++index) {
{
a_parVars->m_taken[index] = false;		a_parVars->m_taken[index] = false;
a_parVars->m_partition[index] = -1;		a_parVars->m_partition[index] = -1;
}		}
}		}
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::PickSeeds(PartitionVars* a_parVars)		void RTREE_QUAL::PickSeeds(PartitionVars *a_parVars)
{		{
int seed0 = -1, seed1 = -1;		int seed0 = -1, seed1 = -1;
ELEMTYPEREAL worst, waste;		ELEMTYPEREAL worst, waste;
ELEMTYPEREAL area[MAXNODES+1];		ELEMTYPEREAL area[MAXNODES + 1];
for(int index=0; index<a_parVars->m_total; ++index)		for(int index = 0; index < a_parVars->m_total; ++index) {
{
area[index] = CalcRectVolume(&a_parVars->m_branchBuf[index].m_rect);		area[index] = CalcRectVolume(&a_parVars->m_branchBuf[index].m_rect);
}		}
worst = -a_parVars->m_coverSplitArea - 1;		worst = -a_parVars->m_coverSplitArea - 1;
for(int indexA=0; indexA < a_parVars->m_total-1; ++indexA)		for(int indexA = 0; indexA < a_parVars->m_total - 1; ++indexA) {
{		for(int indexB = indexA + 1; indexB < a_parVars->m_total; ++indexB) {
for(int indexB = indexA+1; indexB < a_parVars->m_total; ++indexB)		Rect oneRect = CombineRect(&a_parVars->m_branchBuf[indexA].m_rect,
{		&a_parVars->m_branchBuf[indexB].m_rect);
Rect oneRect = CombineRect(&a_parVars->m_branchBuf[indexA].m_rect, &a_parV
ars->m_branchBuf[indexB].m_rect);
waste = CalcRectVolume(&oneRect) - area[indexA] - area[indexB];		waste = CalcRectVolume(&oneRect) - area[indexA] - area[indexB];
if(waste > worst)		if(waste > worst) {
{
worst = waste;		worst = waste;
seed0 = indexA;		seed0 = indexA;
seed1 = indexB;		seed1 = indexB;
}		}
}		}
}		}
Classify(seed0, 0, a_parVars);		Classify(seed0, 0, a_parVars);
Classify(seed1, 1, a_parVars);		Classify(seed1, 1, a_parVars);
}		}
// Put a branch in one of the groups.		// Put a branch in one of the groups.
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::Classify(int a_index, int a_group, PartitionVars* a_parVars)		void RTREE_QUAL::Classify(int a_index, int a_group, PartitionVars *a_parVars)
{		{
ASSERT(a_parVars);		ASSERT(a_parVars);
ASSERT(!a_parVars->m_taken[a_index]);		ASSERT(!a_parVars->m_taken[a_index]);
ASSERT(a_index >= 0); // Added for Gmsh		ASSERT(a_index >= 0); // Added for Gmsh
a_parVars->m_partition[a_index] = a_group;		a_parVars->m_partition[a_index] = a_group;
a_parVars->m_taken[a_index] = true;		a_parVars->m_taken[a_index] = true;
if (a_parVars->m_count[a_group] == 0)		if(a_parVars->m_count[a_group] == 0) {
{
a_parVars->m_cover[a_group] = a_parVars->m_branchBuf[a_index].m_rect;		a_parVars->m_cover[a_group] = a_parVars->m_branchBuf[a_index].m_rect;
}		}
else		else {
{		a_parVars->m_cover[a_group] = CombineRect(
a_parVars->m_cover[a_group] = CombineRect(&a_parVars->m_branchBuf[a_index].m		&a_parVars->m_branchBuf[a_index].m_rect, &a_parVars->m_cover[a_group]);
_rect, &a_parVars->m_cover[a_group]);
}		}
a_parVars->m_area[a_group] = CalcRectVolume(&a_parVars->m_cover[a_group]);		a_parVars->m_area[a_group] = CalcRectVolume(&a_parVars->m_cover[a_group]);
++a_parVars->m_count[a_group];		++a_parVars->m_count[a_group];
}		}
// Delete a data rectangle from an index structure.		// Delete a data rectangle from an index structure.
// Pass in a pointer to a Rect, the tid of the record, ptr to ptr to root node.		// Pass in a pointer to a Rect, the tid of the record, ptr to ptr to root node.
// Returns 1 if record not found, 0 if success.		// Returns 1 if record not found, 0 if success.
// RemoveRect provides for eliminating the root.		// RemoveRect provides for eliminating the root.
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::RemoveRect(Rect* a_rect, const DATATYPE& a_id, Node** a_root)		bool RTREE_QUAL::RemoveRect(Rect *a_rect, const DATATYPE &a_id, Node **a_root)
{		{
ASSERT(a_rect && a_root);		ASSERT(a_rect && a_root);
ASSERT(*a_root);		ASSERT(*a_root);
Node* tempNode;		Node *tempNode;
ListNode* reInsertList = NULL;		ListNode *reInsertList = NULL;
if(!RemoveRectRec(a_rect, a_id, *a_root, &reInsertList))		if(!RemoveRectRec(a_rect, a_id, *a_root, &reInsertList)) {
{
// Found and deleted a data item		// Found and deleted a data item
// Reinsert any branches from eliminated nodes		// Reinsert any branches from eliminated nodes
while(reInsertList)		while(reInsertList) {
{
tempNode = reInsertList->m_node;		tempNode = reInsertList->m_node;
for(int index = 0; index < tempNode->m_count; ++index)		for(int index = 0; index < tempNode->m_count; ++index) {
{
InsertRect(&(tempNode->m_branch[index].m_rect),		InsertRect(&(tempNode->m_branch[index].m_rect),
tempNode->m_branch[index].m_data,		tempNode->m_branch[index].m_data, a_root, tempNode->m_level);
a_root,
tempNode->m_level);
}		}
ListNode* remLNode = reInsertList;		ListNode *remLNode = reInsertList;
reInsertList = reInsertList->m_next;		reInsertList = reInsertList->m_next;
FreeNode(remLNode->m_node);		FreeNode(remLNode->m_node);
FreeListNode(remLNode);		FreeListNode(remLNode);
}		}
// Check for redundant root (not leaf, 1 child) and eliminate		// Check for redundant root (not leaf, 1 child) and eliminate
if((*a_root)->m_count == 1 && (*a_root)->IsInternalNode())		if((*a_root)->m_count == 1 && (*a_root)->IsInternalNode()) {
{
tempNode = (*a_root)->m_branch[0].m_child;		tempNode = (*a_root)->m_branch[0].m_child;
ASSERT(tempNode);		ASSERT(tempNode);
FreeNode(*a_root);		FreeNode(*a_root);
*a_root = tempNode;		*a_root = tempNode;
}		}
return false;		return false;
}		}
else		else {
{
return true;		return true;
}		}
}		}
// Delete a rectangle from non-root part of an index structure.		// Delete a rectangle from non-root part of an index structure.
// Called by RemoveRect. Descends tree recursively,		// Called by RemoveRect. Descends tree recursively,
// merges branches on the way back up.		// merges branches on the way back up.
// Returns 1 if record not found, 0 if success.		// Returns 1 if record not found, 0 if success.
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::RemoveRectRec(Rect* a_rect, const DATATYPE& a_id, Node* a_node,		bool RTREE_QUAL::RemoveRectRec(Rect *a_rect, const DATATYPE &a_id, Node *a_node,
ListNode** a_listNode)		ListNode **a_listNode)
{		{
ASSERT(a_rect && a_node && a_listNode);		ASSERT(a_rect && a_node && a_listNode);
ASSERT(a_node->m_level >= 0);		ASSERT(a_node->m_level >= 0);
if(a_node->IsInternalNode()) // not a leaf node		if(a_node->IsInternalNode()) // not a leaf node
{		{
for(int index = 0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{		if(Overlap(a_rect, &(a_node->m_branch[index].m_rect))) {
if(Overlap(a_rect, &(a_node->m_branch[index].m_rect)))		if(!RemoveRectRec(a_rect, a_id, a_node->m_branch[index].m_child,
{		a_listNode)) {
if(!RemoveRectRec(a_rect, a_id, a_node->m_branch[index].m_child, a_listN		if(a_node->m_branch[index].m_child->m_count >= MINNODES) {
ode))
{
if(a_node->m_branch[index].m_child->m_count >= MINNODES)
{
// child removed, just resize parent rect		// child removed, just resize parent rect
a_node->m_branch[index].m_rect = NodeCover(a_node->m_branch[index].m		a_node->m_branch[index].m_rect =
_child);		NodeCover(a_node->m_branch[index].m_child);
}		}
else		else {
{
// child removed, not enough entries in node, eliminate node		// child removed, not enough entries in node, eliminate node
ReInsert(a_node->m_branch[index].m_child, a_listNode);		ReInsert(a_node->m_branch[index].m_child, a_listNode);
DisconnectBranch(a_node, index); // Must return after this call as c		DisconnectBranch(
ount has changed		a_node,
		index); // Must return after this call as count has changed
}		}
return false;		return false;
}		}
}		}
}		}
return true;		return true;
}		}
else // A leaf node		else // A leaf node
{		{
for(int index = 0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{		if(a_node->m_branch[index].m_child == (Node *)a_id) {
if(a_node->m_branch[index].m_child == (Node*)a_id)		DisconnectBranch(
{		a_node, index); // Must return after this call as count has changed
DisconnectBranch(a_node, index); // Must return after this call as count
has changed
return false;		return false;
}		}
}		}
return true;		return true;
}		}
}		}
// Decide whether two rectangles overlap.		// Decide whether two rectangles overlap.
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::Overlap(Rect* a_rectA, Rect* a_rectB)		bool RTREE_QUAL::Overlap(Rect *a_rectA, Rect *a_rectB)
{		{
ASSERT(a_rectA && a_rectB);		ASSERT(a_rectA && a_rectB);
for(int index=0; index < NUMDIMS; ++index)		for(int index = 0; index < NUMDIMS; ++index) {
{		if(a_rectA->m_min[index] > a_rectB->m_max[index] ||
if (a_rectA->m_min[index] > a_rectB->m_max[index] ||		a_rectB->m_min[index] > a_rectA->m_max[index]) {
a_rectB->m_min[index] > a_rectA->m_max[index])
{
return false;		return false;
}		}
}		}
return true;		return true;
}		}
// Add a node to the reinsertion list. All its branches will later		// Add a node to the reinsertion list. All its branches will later
// be reinserted into the index structure.		// be reinserted into the index structure.
RTREE_TEMPLATE		RTREE_TEMPLATE
void RTREE_QUAL::ReInsert(Node* a_node, ListNode** a_listNode)		void RTREE_QUAL::ReInsert(Node *a_node, ListNode **a_listNode)
{		{
ListNode* newListNode;		ListNode *newListNode;
newListNode = AllocListNode();		newListNode = AllocListNode();
newListNode->m_node = a_node;		newListNode->m_node = a_node;
newListNode->m_next = *a_listNode;		newListNode->m_next = *a_listNode;
*a_listNode = newListNode;		*a_listNode = newListNode;
}		}
// Search in an index tree or subtree for all data retangles that overlap the ar		// Search in an index tree or subtree for all data retangles that overlap the
gument rectangle.		// argument rectangle.
RTREE_TEMPLATE		RTREE_TEMPLATE
bool RTREE_QUAL::Search(Node* a_node, Rect* a_rect, int& a_foundCount, bool a_re		bool RTREE_QUAL::Search(Node *a_node, Rect *a_rect, int &a_foundCount,
sultCallback(DATATYPE a_data, void* a_context), void* a_context)		bool a_resultCallback(DATATYPE a_data, void *a_context),
		void *a_context)
{		{
ASSERT(a_node);		ASSERT(a_node);
ASSERT(a_node->m_level >= 0);		ASSERT(a_node->m_level >= 0);
ASSERT(a_rect);		ASSERT(a_rect);
if(a_node->IsInternalNode()) // This is an internal node in the tree		if(a_node->IsInternalNode()) // This is an internal node in the tree
{		{
for(int index=0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{		if(Overlap(a_rect, &a_node->m_branch[index].m_rect)) {
if(Overlap(a_rect, &a_node->m_branch[index].m_rect))		if(!Search(a_node->m_branch[index].m_child, a_rect, a_foundCount,
{		a_resultCallback, a_context)) {
if(!Search(a_node->m_branch[index].m_child, a_rect, a_foundCount, a_resu
ltCallback, a_context))
{
return false; // Don't continue searching		return false; // Don't continue searching
}		}
}		}
}		}
}		}
else // This is a leaf node		else // This is a leaf node
{		{
for(int index=0; index < a_node->m_count; ++index)		for(int index = 0; index < a_node->m_count; ++index) {
{		if(Overlap(a_rect, &a_node->m_branch[index].m_rect)) {
if(Overlap(a_rect, &a_node->m_branch[index].m_rect))		DATATYPE &id = a_node->m_branch[index].m_data;
{
DATATYPE& id = a_node->m_branch[index].m_data;		// NOTE: There are different ways to return results. Here's where to
		// modify
// NOTE: There are different ways to return results. Here's where to mo		// if(&a_resultCallback)
dify
//if(&a_resultCallback)
{		{
++a_foundCount;		++a_foundCount;
if(!a_resultCallback(id, a_context))		if(!a_resultCallback(id, a_context)) {
{
return false; // Don't continue searching		return false; // Don't continue searching
}		}
}		}
}		}
}		}
}		}
return true; // Continue searching		return true; // Continue searching
}		}
#undef RTREE_TEMPLATE		#undef RTREE_TEMPLATE
#undef RTREE_QUAL		#undef RTREE_QUAL
#endif //RTREE_H		#endif // RTREE_H
End of changes. 220 change blocks.
573 lines changed or deleted		431 lines changed or added

---