netgen/libsrc/gprim/adtree.hpp
Matthias Hochsteger 16b88e8e67 OCC Mesher Cleanup
2021-11-28 15:14:41 +00:00

1385 lines
31 KiB
C++

#ifndef FILE_ADTREE
#define FILE_ADTREE
/* *************************************************************************/
/* File: adtree.hh */
/* Author: Joachim Schoeberl */
/* Date: 16. Feb. 98 */
/* Redesigned by Wolfram Muehlhuber, May 1998 */
/* *************************************************************************/
namespace netgen
{
/**
Alternating Digital Tree
*/
// #include "../include/mystdlib.h"
// #include "../include/myadt.hpp"
class ADTreeNode
{
public:
ADTreeNode *left, *right, *father;
int dim;
float sep;
float *data;
float *boxmin;
float *boxmax;
int pi;
int nchilds;
ADTreeNode (int adim);
~ADTreeNode ();
friend class ADTree;
};
class ADTreeCriterion
{
public:
ADTreeCriterion() { }
virtual int Eval (const ADTreeNode * node) const = 0;
};
class ADTree
{
int dim;
ADTreeNode * root;
float *cmin, *cmax;
NgArray<ADTreeNode*> ela;
const ADTreeCriterion * criterion;
NgArray<ADTreeNode*> stack;
NgArray<int> stackdir;
int stackindex;
public:
ADTree (int adim, const float * acmin,
const float * acmax);
~ADTree ();
void Insert (const float * p, int pi);
// void GetIntersecting (const float * bmin, const float * bmax,
// NgArray<int> & pis) const;
void SetCriterion (ADTreeCriterion & acriterion);
void Reset ();
int Next ();
void GetMatch (NgArray<int> & matches);
void DeleteElement (int pi);
void Print (ostream & ost) const
{ PrintRec (ost, root); }
void PrintRec (ostream & ost, const ADTreeNode * node) const;
};
class ADTreeNode3
{
public:
ADTreeNode3 *left, *right, *father;
float sep;
float data[3];
int pi;
int nchilds;
ADTreeNode3 ();
void DeleteChilds ();
friend class ADTree3;
static BlockAllocator ball;
void * operator new(size_t);
void operator delete (void *);
};
class ADTree3
{
ADTreeNode3 * root;
float cmin[3], cmax[3];
NgArray<ADTreeNode3*> ela;
public:
ADTree3 (const float * acmin,
const float * acmax);
~ADTree3 ();
void Insert (const float * p, int pi);
void GetIntersecting (const float * bmin, const float * bmax,
NgArray<int> & pis) const;
void DeleteElement (int pi);
void Print (ostream & ost) const
{ PrintRec (ost, root); }
void PrintRec (ostream & ost, const ADTreeNode3 * node) const;
};
/*
// divide each direction
#define ADTN_DIV 10
class ADTreeNode3Div
{
public:
ADTreeNode3Div *father;
ADTreeNode3Div *childs[ADTN_DIV];
float minx, dist;
float data[3];
int pi;
int nchilds;
ADTreeNode3Div ();
void DeleteChilds ();
friend class ADTree3Div;
static BlockAllocator ball;
void * operator new(size_t);
void operator delete (void *);
};
class ADTree3Div
{
ADTreeNode3Div * root;
float cmin[3], cmax[3];
NgArray<ADTreeNode3Div*> ela;
public:
ADTree3Div (const float * acmin,
const float * acmax);
~ADTree3Div ();
void Insert (const float * p, int pi);
void GetIntersecting (const float * bmin, const float * bmax,
NgArray<int> & pis) const;
void DeleteElement (int pi);
void Print (ostream & ost) const
{ PrintRec (ost, root); }
void PrintRec (ostream & ost, const ADTreeNode3Div * node) const;
};
#define ADTN_SIZE 10
// multiple entries
class ADTreeNode3M
{
public:
ADTreeNode3M *left, *right, *father;
float sep;
float data[ADTN_SIZE][3];
int pi[ADTN_SIZE];
int nchilds;
ADTreeNode3M ();
void DeleteChilds ();
friend class ADTree3M;
static BlockAllocator ball;
void * operator new(size_t);
void operator delete (void *);
};
class ADTree3M
{
ADTreeNode3M * root;
float cmin[3], cmax[3];
NgArray<ADTreeNode3M*> ela;
public:
ADTree3M (const float * acmin,
const float * acmax);
~ADTree3M ();
void Insert (const float * p, int pi);
void GetIntersecting (const float * bmin, const float * bmax,
NgArray<int> & pis) const;
void DeleteElement (int pi);
void Print (ostream & ost) const
{ PrintRec (ost, root); }
void PrintRec (ostream & ost, const ADTreeNode3M * node) const;
};
class ADTreeNode3F
{
public:
ADTreeNode3F *father;
ADTreeNode3F *childs[8];
float sep[3];
float data[3];
int pi;
int nchilds;
ADTreeNode3F ();
void DeleteChilds ();
friend class ADTree3F;
static BlockAllocator ball;
void * operator new(size_t);
void operator delete (void *);
};
// fat tree
class ADTree3F
{
ADTreeNode3F * root;
float cmin[3], cmax[3];
NgArray<ADTreeNode3F*> ela;
public:
ADTree3F (const float * acmin,
const float * acmax);
~ADTree3F ();
void Insert (const float * p, int pi);
void GetIntersecting (const float * bmin, const float * bmax,
NgArray<int> & pis) const;
void DeleteElement (int pi);
void Print (ostream & ost) const
{ PrintRec (ost, root); }
void PrintRec (ostream & ost, const ADTreeNode3F * node) const;
};
class ADTreeNode3FM
{
public:
ADTreeNode3FM *father;
ADTreeNode3FM *childs[8];
float sep[3];
float data[ADTN_SIZE][3];
int pi[ADTN_SIZE];
int nchilds;
ADTreeNode3FM ();
void DeleteChilds ();
friend class ADTree3FM;
static BlockAllocator ball;
void * operator new(size_t);
void operator delete (void *);
};
// fat tree
class ADTree3FM
{
ADTreeNode3FM * root;
float cmin[3], cmax[3];
NgArray<ADTreeNode3FM*> ela;
public:
ADTree3FM (const float * acmin,
const float * acmax);
~ADTree3FM ();
void Insert (const float * p, int pi);
void GetIntersecting (const float * bmin, const float * bmax,
NgArray<int> & pis) const;
void DeleteElement (int pi);
void Print (ostream & ost) const
{ PrintRec (ost, root); }
void PrintRec (ostream & ost, const ADTreeNode3FM * node) const;
};
*/
class ADTreeNode6
{
public:
ADTreeNode6 *left, *right, *father;
float sep;
float data[6];
int pi;
int nchilds;
ADTreeNode6 ();
void DeleteChilds ();
friend class ADTree6;
static BlockAllocator ball;
void * operator new(size_t);
void operator delete (void *);
};
class ADTree6
{
ADTreeNode6 * root;
float cmin[6], cmax[6];
NgArray<ADTreeNode6*> ela;
public:
ADTree6 (const float * acmin,
const float * acmax);
~ADTree6 ();
void Insert (const float * p, int pi);
void GetIntersecting (const float * bmin, const float * bmax,
NgArray<int> & pis) const;
void DeleteElement (int pi);
void Print (ostream & ost) const
{ PrintRec (ost, root); }
int Depth () const
{ return DepthRec (root); }
int Elements () const
{ return ElementsRec (root); }
void PrintRec (ostream & ost, const ADTreeNode6 * node) const;
int DepthRec (const ADTreeNode6 * node) const;
int ElementsRec (const ADTreeNode6 * node) const;
void PrintMemInfo (ostream & ost) const;
};
template <int DIM, typename T>
class T_ADTreeNode
{
public:
T_ADTreeNode *left, *right, *father;
float sep;
// float data[DIM];
Point<DIM,float> data;
T pi;
int nchilds;
T_ADTreeNode ()
{
// pi = -1;
SetInvalid(pi);
left = NULL;
right = NULL;
father = NULL;
nchilds = 0;
}
void DeleteChilds (BlockAllocator & ball)
{
if (left)
{
left->DeleteChilds(ball);
ball.Free(left);
left = NULL;
}
if (right)
{
right->DeleteChilds(ball);
ball.Free(right);
right = NULL;
}
}
};
template <int dim, typename T = INDEX>
class T_ADTree
{
T_ADTreeNode<dim,T> * root;
// float cmin[dim], cmax[dim];
Point<dim> cmin, cmax;
// NgArray<T_ADTreeNode<dim>*> ela;
ClosedHashTable<T, T_ADTreeNode<dim,T>*> ela;
BlockAllocator ball{sizeof(T_ADTreeNode<dim,T>)};
public:
T_ADTree (Point<dim> acmin, Point<dim> acmax)
{
cmin = acmin;
cmax = acmax;
root = new (ball.Alloc()) T_ADTreeNode<dim,T>;
root->sep = (cmin[0] + cmax[0]) / 2;
}
~T_ADTree ()
{
root->DeleteChilds(ball);
ball.Free(root);
}
void Insert (Point<dim> p, T pi)
{
T_ADTreeNode<dim,T> *node(NULL);
T_ADTreeNode<dim,T> *next;
int dir;
int lr(0);
Point<dim> bmin = cmin;
Point<dim> bmax = cmax;
next = root;
dir = 0;
while (next)
{
node = next;
if (IsInvalid(node->pi))
{
// memcpy (node->data, p, dim * sizeof(float));
node->data = p;
node->pi = pi;
// if (ela.Size() < pi+1)
// ela.SetSize (pi+1);
ela[pi] = node;
return;
}
if (node->sep > p[dir])
{
next = node->left;
bmax(dir) = node->sep;
lr = 0;
}
else
{
next = node->right;
bmin(dir) = node->sep;
lr = 1;
}
dir++;
if (dir == dim) dir = 0;
}
next = new (ball.Alloc()) T_ADTreeNode<dim,T>;
next->data = p;
next->pi = pi;
next->sep = (bmin[dir] + bmax[dir]) / 2;
// if (ela.Size() < pi+1)
// ela.SetSize (pi+1);
ela[pi] = next;
if (lr)
node->right = next;
else
node->left = next;
next -> father = node;
while (node)
{
node->nchilds++;
node = node->father;
}
}
class inttn {
public:
int dir;
T_ADTreeNode<dim,T> * node;
};
void GetIntersecting (Point<dim> bmin, Point<dim> bmax,
NgArray<T> & pis) const
{
NgArrayMem<inttn,10000> stack(10000);
pis.SetSize(0);
stack[0].node = root;
stack[0].dir = 0;
int stacks = 0;
while (stacks >= 0)
{
T_ADTreeNode<dim,T> * node = stack[stacks].node;
int dir = stack[stacks].dir;
stacks--;
if (!IsInvalid(node->pi)) // != -1)
{
bool found = true;
for (int i = 0; i < dim/2; i++)
if (node->data[i] > bmax[i])
found = false;
for (int i = dim/2; i < dim; i++)
if (node->data[i] < bmin[i])
found = false;
if (found)
pis.Append (node->pi);
/*
if (node->data[0] > bmax[0] ||
node->data[1] > bmax[1] ||
node->data[2] > bmax[2] ||
node->data[3] < bmin[3] ||
node->data[4] < bmin[4] ||
node->data[5] < bmin[5])
;
else
{
pis.Append (node->pi);
}
*/
}
int ndir = (dir+1) % dim;
if (node->left && bmin[dir] <= node->sep)
{
stacks++;
stack[stacks].node = node->left;
stack[stacks].dir = ndir;
}
if (node->right && bmax[dir] >= node->sep)
{
stacks++;
stack[stacks].node = node->right;
stack[stacks].dir = ndir;
}
}
}
void DeleteElement (T pi)
{
T_ADTreeNode<dim,T> * node = ela[pi];
ela.Delete(pi);
SetInvalid(node->pi); // = -1;
node = node->father;
while (node)
{
node->nchilds--;
node = node->father;
}
}
void Print (ostream & ost) const
{ PrintRec (ost, root); }
int Depth () const
{ return DepthRec (root); }
int Elements () const
{ return ElementsRec (root); }
void PrintRec (ostream & ost, const T_ADTreeNode<dim,T> * node) const
{
// if (node->data) // true anyway
{
ost << node->pi << ": ";
ost << node->nchilds << " childs, ";
for (int i = 0; i < dim; i++)
ost << node->data[i] << " ";
ost << endl;
}
if (node->left)
PrintRec (ost, node->left);
if (node->right)
PrintRec (ost, node->right);
}
int DepthRec (const T_ADTreeNode<dim,T> * node) const
{
int ldepth = 0;
int rdepth = 0;
if (node->left)
ldepth = DepthRec(node->left);
if (node->right)
rdepth = DepthRec(node->right);
return 1 + max2 (ldepth, rdepth);
}
int ElementsRec (const T_ADTreeNode<dim,T> * node) const
{
int els = 1;
if (node->left)
els += ElementsRec(node->left);
if (node->right)
els += ElementsRec(node->right);
return els;
}
void PrintMemInfo (ostream & ost) const
{
ost << Elements() << " elements a " << sizeof(ADTreeNode6)
<< " Bytes = "
<< Elements() * sizeof(T_ADTreeNode<dim,T>) << endl;
ost << "maxind = " << ela.Size() << " = " << sizeof(T_ADTreeNode<dim,T>*) * ela.Size() << " Bytes" << endl;
}
};
/*
class ADTreeNode6F
{
public:
ADTreeNode6F * father;
ADTreeNode6F * childs[64];
float sep[6];
float data[6];
int pi;
int nchilds;
ADTreeNode6F ();
void DeleteChilds ();
friend class ADTree6F;
static BlockAllocator ball;
void * operator new(size_t);
void operator delete (void *);
};
class ADTree6F
{
ADTreeNode6F * root;
float cmin[6], cmax[6];
NgArray<ADTreeNode6F*> ela;
public:
ADTree6F (const float * acmin,
const float * acmax);
~ADTree6F ();
void Insert (const float * p, int pi);
void GetIntersecting (const float * bmin, const float * bmax,
NgArray<int> & pis) const;
void DeleteElement (int pi);
void Print (ostream & ost) const
{ PrintRec (ost, root); }
int Depth () const
{ return DepthRec (root); }
void PrintRec (ostream & ost, const ADTreeNode6F * node) const;
int DepthRec (const ADTreeNode6F * node) const;
};
*/
class Point3dTree
{
ADTree3 * tree;
public:
DLL_HEADER Point3dTree (const Point<3> & pmin, const Point<3> & pmax);
DLL_HEADER ~Point3dTree ();
DLL_HEADER void Insert (const Point<3> & p, int pi);
void DeleteElement (int pi)
{ tree->DeleteElement(pi); }
DLL_HEADER void GetIntersecting (const Point<3> & pmin, const Point<3> & pmax,
NgArray<int> & pis) const;
const ADTree3 & Tree() const { return *tree; };
};
template<int dim, typename T=INDEX>
class BoxTree
{
public:
// Number of entries per leaf
static constexpr int N = 100;
struct Node;
struct Leaf
{
Point<2*dim> p[N];
T index[N];
int n_elements;
Leaf() : n_elements(0)
{ }
void Add( ClosedHashTable<T, Leaf*> &leaf_index, const Point<2*dim> &ap, T aindex )
{
p[n_elements] = ap;
index[n_elements] = aindex;
n_elements++;
leaf_index[aindex] = this;
}
};
struct Node
{
union
{
Node *children[2];
Leaf *leaf;
};
double sep;
int level;
Node()
: children{nullptr,nullptr}
{ }
~Node()
{ }
Leaf *GetLeaf() const
{
return children[1] ? nullptr : leaf;
}
};
private:
Node root;
ClosedHashTable<T, Leaf*> leaf_index;
Point<dim> global_min, global_max;
double tol;
size_t n_leaves;
size_t n_nodes;
BlockAllocator ball_nodes;
BlockAllocator ball_leaves;
public:
BoxTree (const Point<dim> & pmin, const Point<dim> & pmax)
: global_min(pmin), global_max(pmax), n_leaves(1), n_nodes(1), ball_nodes(sizeof(Node)), ball_leaves(sizeof(Leaf))
{
root.leaf = (Leaf*) ball_leaves.Alloc(); new (root.leaf) Leaf();
root.level = 0;
tol = 1e-7 * Dist(pmax, pmin);
}
BoxTree (const Box<dim> & box)
: BoxTree(box.PMin(), box.PMax())
{ }
double GetTolerance() { return tol; }
size_t GetNLeaves()
{
return n_leaves;
}
size_t GetNNodes()
{
return n_nodes;
}
template<typename TFunc>
void GetFirstIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
TFunc func=[](auto pi){return false;}) const
{
// static Timer timer("BoxTree::GetIntersecting"); RegionTimer rt(timer);
// static Timer timer1("BoxTree::GetIntersecting-LinearSearch");
ArrayMem<const Node*, 100> stack;
ArrayMem<int, 100> dir_stack;
Point<2*dim> tpmin, tpmax;
for (size_t i : IntRange(dim))
{
tpmin(i) = global_min(i);
tpmax(i) = pmax(i)+tol;
tpmin(i+dim) = pmin(i)-tol;
tpmax(i+dim) = global_max(i);
}
stack.SetSize(0);
stack.Append(&root);
dir_stack.SetSize(0);
dir_stack.Append(0);
while(stack.Size())
{
const Node *node = stack.Last();
stack.DeleteLast();
int dir = dir_stack.Last();
dir_stack.DeleteLast();
if(Leaf *leaf = node->GetLeaf())
{
// RegionTimer rt1(timer1);
for (auto i : IntRange(leaf->n_elements))
{
bool intersect = true;
const auto p = leaf->p[i];
for (int d = 0; d < dim; d++)
if (p[d] > tpmax[d])
intersect = false;
for (int d = dim; d < 2*dim; d++)
if (p[d] < tpmin[d])
intersect = false;
if(intersect)
if(func(leaf->index[i])) return;
}
}
else
{
int newdir = dir+1;
if(newdir==2*dim) newdir = 0;
if (tpmin[dir] <= node->sep)
{
stack.Append(node->children[0]);
dir_stack.Append(newdir);
}
if (tpmax[dir] >= node->sep)
{
stack.Append(node->children[1]);
dir_stack.Append(newdir);
}
}
}
}
void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
NgArray<T> & pis) const
{
pis.SetSize(0);
GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;});
}
void GetIntersecting(const Point<dim> & pmin,
const Point<dim> & pmax,
Array<T> & pis) const
{
pis.SetSize0();
GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;});
}
void Insert (const Box<dim> & box, T pi)
{
Insert (box.PMin(), box.PMax(), pi);
}
void Insert (const Point<dim> & pmin, const Point<dim> & pmax, T pi)
{
// static Timer timer("BoxTree::Insert"); RegionTimer rt(timer);
int dir = 0;
Point<2*dim> p;
for (auto i : IntRange(dim))
{
p(i) = pmin[i];
p(i+dim) = pmax[i];
}
Node * node = &root;
Leaf * leaf = node->GetLeaf();
// search correct leaf to add point
while(!leaf)
{
node = p[dir] < node->sep ? node->children[0] : node->children[1];
dir++;
if(dir==2*dim) dir = 0;
leaf = node->GetLeaf();
}
// add point to leaf
if(leaf->n_elements < N)
leaf->Add(leaf_index, p,pi);
else // assume leaf->n_elements == N
{
// add two new nodes and one new leaf
int n_elements = leaf->n_elements;
ArrayMem<double, N> coords(n_elements);
ArrayMem<int, N> order(n_elements);
// separate points in two halves, first sort all coordinates in direction dir
for (auto i : IntRange(n_elements))
{
order[i] = i;
coords[i] = leaf->p[i][dir];
}
QuickSortI(coords, order);
int isplit = N/2;
Leaf *leaf1 = (Leaf*) ball_leaves.Alloc(); new (leaf1) Leaf();
Leaf *leaf2 = (Leaf*) ball_leaves.Alloc(); new (leaf2) Leaf();
for (auto i : order.Range(isplit))
leaf1->Add(leaf_index, leaf->p[i], leaf->index[i] );
for (auto i : order.Range(isplit, N))
leaf2->Add(leaf_index, leaf->p[i], leaf->index[i] );
Node *node1 = (Node*) ball_nodes.Alloc(); new (node1) Node();
node1->leaf = leaf1;
node1->level = node->level+1;
Node *node2 = (Node*) ball_nodes.Alloc(); new (node2) Node();
node2->leaf = leaf2;
node2->level = node->level+1;
node->children[0] = node1;
node->children[1] = node2;
node->sep = 0.5 * (leaf->p[order[isplit-1]][dir] + leaf->p[order[isplit]][dir]);
// add new point to one of the new leaves
if (p[dir] < node->sep)
leaf1->Add( leaf_index, p, pi );
else
leaf2->Add( leaf_index, p, pi );
ball_leaves.Free(leaf);
n_leaves++;
n_nodes+=2;
}
}
void DeleteElement (T pi)
{
// static Timer timer("BoxTree::DeleteElement"); RegionTimer rt(timer);
Leaf *leaf = leaf_index[pi];
leaf_index.Delete(pi);
auto & n_elements = leaf->n_elements;
auto & index = leaf->index;
auto & p = leaf->p;
for (auto i : IntRange(n_elements))
{
if(index[i] == pi)
{
n_elements--;
if(i!=n_elements)
{
index[i] = index[n_elements];
p[i] = p[n_elements];
}
return;
}
}
}
};
// template <int dim, typename T = INDEX>
// class BoxTree
// {
// T_ADTree<2*dim,T> * tree;
// Point<dim> boxpmin, boxpmax;
// public:
// BoxTree (const Box<dim> & abox)
// {
// boxpmin = abox.PMin();
// boxpmax = abox.PMax();
// Point<2*dim> tpmin, tpmax;
// for (int i = 0; i < dim; i++)
// {
// tpmin(i) = tpmin(i+dim) = boxpmin(i);
// tpmax(i) = tpmax(i+dim) = boxpmax(i);
// }
// tree = new T_ADTree<2*dim,T> (tpmin, tpmax);
// }
//
// BoxTree (const Point<dim> & apmin, const Point<dim> & apmax)
// {
// boxpmin = apmin;
// boxpmax = apmax;
// Point<2*dim> tpmin, tpmax;
// for (int i = 0; i < dim; i++)
// {
// tpmin(i) = tpmin(i+dim) = boxpmin(i);
// tpmax(i) = tpmax(i+dim) = boxpmax(i);
// }
// tree = new T_ADTree<2*dim,T> (tpmin, tpmax);
// }
//
// ~BoxTree ()
// {
// delete tree;
// }
//
// void Insert (const Point<dim> & bmin, const Point<dim> & bmax, T pi)
// {
// Point<2*dim> tp;
//
// for (size_t i = 0; i < dim; i++)
// {
// tp(i) = bmin(i);
// tp(i+dim) = bmax(i);
// }
//
// tree->Insert (tp, pi);
// }
//
// void Insert (const Box<dim> & box, T pi)
// {
// Insert (box.PMin(), box.PMax(), pi);
// }
//
// void DeleteElement (T pi)
// {
// tree->DeleteElement(pi);
// }
//
// void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
// NgArray<T> & pis) const
// {
// Point<2*dim> tpmin, tpmax;
// double tol = Tolerance();
// for (size_t i = 0; i < dim; i++)
// {
// tpmin(i) = boxpmin(i);
// tpmax(i) = pmax(i)+tol;
//
// tpmin(i+dim) = pmin(i)-tol;
// tpmax(i+dim) = boxpmax(i);
// }
//
// tree->GetIntersecting (tpmin, tpmax, pis);
// }
//
//
// double Tolerance() const { return 1e-7 * Dist(boxpmax, boxpmin); } // single precision
// const auto & Tree() const { return *tree; };
// auto & Tree() { return *tree; };
// };
template<int dim, typename T=INDEX, typename TSCAL=double>
class DelaunayTree
{
public:
// Number of entries per leaf
static constexpr int N = 100;
struct Node;
struct Leaf
{
Point<2*dim, TSCAL> p[N];
T index[N];
int n_elements;
int nr;
Leaf() : n_elements(0)
{ }
void Add( Array<Leaf*> &leaves, Array<T> &leaf_index, const Point<2*dim> &ap, T aindex )
{
p[n_elements] = ap;
index[n_elements] = aindex;
n_elements++;
if(leaf_index.Size()<aindex+1)
leaf_index.SetSize(aindex+1);
leaf_index[aindex] = nr;
}
};
struct Node
{
union
{
Node *children[2];
Leaf *leaf;
};
double sep;
int level;
Node()
: children{nullptr,nullptr}
{ }
~Node()
{ }
Leaf *GetLeaf() const
{
return children[1] ? nullptr : leaf;
}
};
private:
Node root;
Array<Leaf*> leaves;
Array<T> leaf_index;
Point<dim> global_min, global_max;
double tol;
size_t n_leaves;
size_t n_nodes;
BlockAllocator ball_nodes;
BlockAllocator ball_leaves;
public:
DelaunayTree (const Point<dim> & pmin, const Point<dim> & pmax)
: global_min(pmin), global_max(pmax), n_leaves(1), n_nodes(1), ball_nodes(sizeof(Node)), ball_leaves(sizeof(Leaf))
{
root.leaf = (Leaf*) ball_leaves.Alloc(); new (root.leaf) Leaf();
root.leaf->nr = 0;
leaves.Append(root.leaf);
root.level = 0;
tol = 1e-7 * Dist(pmax, pmin);
}
DelaunayTree (const Box<dim> & box)
: DelaunayTree(box.PMin(), box.PMax())
{ }
size_t GetNLeaves()
{
return n_leaves;
}
size_t GetNNodes()
{
return n_nodes;
}
template<typename TFunc>
void GetFirstIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
TFunc func=[](auto pi){return false;}) const
{
// static Timer timer("DelaunayTree::GetIntersecting"); RegionTimer rt(timer);
// static Timer timer1("DelaunayTree::GetIntersecting-LinearSearch");
ArrayMem<const Node*, 100> stack;
ArrayMem<int, 100> dir_stack;
Point<2*dim> tpmin, tpmax;
for (size_t i : IntRange(dim))
{
tpmin(i) = global_min(i);
tpmax(i) = pmax(i)+tol;
tpmin(i+dim) = pmin(i)-tol;
tpmax(i+dim) = global_max(i);
}
stack.SetSize(0);
stack.Append(&root);
dir_stack.SetSize(0);
dir_stack.Append(0);
while(stack.Size())
{
const Node *node = stack.Last();
stack.DeleteLast();
int dir = dir_stack.Last();
dir_stack.DeleteLast();
if(Leaf *leaf = node->GetLeaf())
{
// RegionTimer rt1(timer1);
for (auto i : IntRange(leaf->n_elements))
{
bool intersect = true;
const auto p = leaf->p[i];
for (int d = 0; d < dim; d++)
if (p[d] > tpmax[d])
intersect = false;
for (int d = dim; d < 2*dim; d++)
if (p[d] < tpmin[d])
intersect = false;
if(intersect)
if(func(leaf->index[i])) return;
}
}
else
{
int newdir = dir+1;
if(newdir==2*dim) newdir = 0;
if (tpmin[dir] <= node->sep)
{
stack.Append(node->children[0]);
dir_stack.Append(newdir);
}
if (tpmax[dir] >= node->sep)
{
stack.Append(node->children[1]);
dir_stack.Append(newdir);
}
}
}
}
void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
NgArray<T> & pis) const
{
pis.SetSize(0);
GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;});
}
void Insert (const Box<dim> & box, T pi)
{
Insert (box.PMin(), box.PMax(), pi);
}
void Insert (const Point<dim> & pmin, const Point<dim> & pmax, T pi)
{
// static Timer timer("DelaunayTree::Insert"); RegionTimer rt(timer);
int dir = 0;
Point<2*dim> p;
for (auto i : IntRange(dim))
{
p(i) = pmin[i];
p(i+dim) = pmax[i];
}
Node * node = &root;
Leaf * leaf = node->GetLeaf();
// search correct leaf to add point
while(!leaf)
{
node = p[dir] < node->sep ? node->children[0] : node->children[1];
dir++;
if(dir==2*dim) dir = 0;
leaf = node->GetLeaf();
}
// add point to leaf
if(leaf->n_elements < N)
leaf->Add(leaves, leaf_index, p,pi);
else // assume leaf->n_elements == N
{
// add two new nodes and one new leaf
int n_elements = leaf->n_elements;
ArrayMem<TSCAL, N> coords(n_elements);
ArrayMem<int, N> order(n_elements);
// separate points in two halves, first sort all coordinates in direction dir
for (auto i : IntRange(n_elements))
{
order[i] = i;
coords[i] = leaf->p[i][dir];
}
QuickSortI(coords, order);
int isplit = N/2;
Leaf *leaf1 = (Leaf*) ball_leaves.Alloc(); new (leaf1) Leaf();
Leaf *leaf2 = (Leaf*) ball_leaves.Alloc(); new (leaf2) Leaf();
leaf1->nr = leaf->nr;
leaf2->nr = leaves.Size();
leaves.Append(leaf2);
leaves[leaf1->nr] = leaf1;
for (auto i : order.Range(isplit))
leaf1->Add(leaves, leaf_index, leaf->p[i], leaf->index[i] );
for (auto i : order.Range(isplit, N))
leaf2->Add(leaves, leaf_index, leaf->p[i], leaf->index[i] );
Node *node1 = (Node*) ball_nodes.Alloc(); new (node1) Node();
node1->leaf = leaf1;
node1->level = node->level+1;
Node *node2 = (Node*) ball_nodes.Alloc(); new (node2) Node();
node2->leaf = leaf2;
node2->level = node->level+1;
node->children[0] = node1;
node->children[1] = node2;
node->sep = 0.5 * (leaf->p[order[isplit-1]][dir] + leaf->p[order[isplit]][dir]);
// add new point to one of the new leaves
if (p[dir] < node->sep)
leaf1->Add( leaves, leaf_index, p, pi );
else
leaf2->Add( leaves, leaf_index, p, pi );
ball_leaves.Free(leaf);
n_leaves++;
n_nodes+=2;
}
}
void DeleteElement (T pi)
{
// static Timer timer("DelaunayTree::DeleteElement"); RegionTimer rt(timer);
Leaf *leaf = leaves[leaf_index[pi]];
leaf_index[pi] = -1;
auto & n_elements = leaf->n_elements;
auto & index = leaf->index;
auto & p = leaf->p;
for (auto i : IntRange(n_elements))
{
if(index[i] == pi)
{
n_elements--;
if(i!=n_elements)
{
index[i] = index[n_elements];
p[i] = p[n_elements];
}
return;
}
}
}
};
}
#endif