# Speed concerns of octree implementation

For couple of days now I am trying to speed up my Force-Directed graph implementation. So far I've implemented Barnes-Hut algorithm that's using octree to decrease number of computations. I've tested it multiple times and number of force related computations is indeed drastically decreased. Below is the plot of computations to number of nodes without Barns-Hut (blue line) and with (red line): Even though now it should be a lot faster, the truth is that in the matter of speed (time) the upgrade is only few percent.

One part I suppose might be causing this is tree creation and elements in tree placing. Because elements are constantly moving I need to recreate tree each loop until some stop condition is reached. But if I will spend lot of time creating tree I will lost there time I've gained on force computation increase. That's my thinking at least. This is how I am adding elements in my main file loop:

void AddTreeElements(Octree* tree, glm::vec3* boundries, Graph& graph)
{
for(auto& node:graph.NodeVector())
{
node.parent_group = nullptr;
if(node.pos[0] < boundries[1][0] && node.pos[0] > boundries[0][0] &&
node.pos[1] > boundries[4][1] && node.pos[1] < boundries[1][1] &&
node.pos[2] < boundries[0][2] && node.pos[2] > boundries[3][2])
{
continue;
}

if(node.pos[0] < boundries[0][0])
{
boundries[0][0] = node.pos[0]-1.0f;
boundries[3][0] = node.pos[0]-1.0f;
boundries[4][0] = node.pos[0]-1.0f;
boundries[7][0] = node.pos[0]-1.0f;
}
else if(node.pos[0] > boundries[1][0])
{
boundries[1][0] = node.pos[0]+1.0f;
boundries[2][0] = node.pos[0]+1.0f;
boundries[5][0] = node.pos[0]+1.0f;
boundries[6][0] = node.pos[0]+1.0f;
}

if(node.pos[1] < boundries[4][1])
{
boundries[4][1] = node.pos[1]-1.0f;
boundries[5][1] = node.pos[1]-1.0f;
boundries[6][1] = node.pos[1]-1.0f;
boundries[7][1] = node.pos[1]-1.0f;
}
else if(node.pos[1] > boundries[0][1])
{
boundries[0][1] = node.pos[1]+1.0f;
boundries[1][1] = node.pos[1]+1.0f;
boundries[2][1] = node.pos[1]+1.0f;
boundries[3][1] = node.pos[1]+1.0f;
}

if(node.pos[2] < boundries[3][2])
{
boundries[2][2] = node.pos[2]-1.0f;
boundries[3][2] = node.pos[2]-1.0f;
boundries[6][2] = node.pos[2]-1.0f;
boundries[7][2] = node.pos[2]-1.0f;
}
else if(node.pos[2] > boundries[0][2])
{
boundries[0][2] = node.pos[2]+1.0f;
boundries[1][2] = node.pos[2]+1.0f;
boundries[4][2] = node.pos[2]+1.0f;
boundries[5][2] = node.pos[2]+1.0f;
}
}
}


What I am doing here is go through all my elements in graph and add them to tree root. Also, I am extending my box that is representing my octree borders for next loop, so all nodes will fit inside.

Fields important to octree structure update are as follows:

Octree* trees[2][2][2];
glm::vec3 vBoundriesBox[8];
bool leaf;
float combined_weight = 0;
std::vector<Element*> objects;


and part of code responsible for update: (in EDIT #2 I've added new, cleaner code)

#define MAX_LEVELS 5

{
this->objects.push_back(object);
}

void Octree::Update()
{
if(this->objects.size()<=1 || level > MAX_LEVELS)
{
for(Element* Element:this->objects)
{
Element->parent_group = this;
}
return;
}

if(leaf)
{
GenerateChildren();
leaf = false;
}

while (!this->objects.empty())
{
Element* obj = this->objects.back();
this->objects.pop_back();
if(contains(trees[0][0][0],obj))
{
trees[0][0][0]->combined_weight += obj->weight;
} else if(contains(trees[0][0][1],obj))
{
trees[0][0][1]->combined_weight += obj->weight;
} else if(contains(trees[0][1][0],obj))
{
trees[0][1][0]->combined_weight += obj->weight;
} else if(contains(trees[0][1][1],obj))
{
trees[0][1][1]->combined_weight += obj->weight;
} else if(contains(trees[1][0][0],obj))
{
trees[1][0][0]->combined_weight += obj->weight;
} else if(contains(trees[1][0][1],obj))
{
trees[1][0][1]->combined_weight += obj->weight;
} else if(contains(trees[1][1][0],obj))
{
trees[1][1][0]->combined_weight += obj->weight;
} else if(contains(trees[1][1][1],obj))
{
trees[1][1][1]->combined_weight += obj->weight;
}
}

for(int i=0;i<2;i++)
{
for(int j=0;j<2;j++)
{
for(int k=0;k<2;k++)
{
trees[i][j][k]->Update();
}
}
}
}

bool Octree::contains(Octree* child, Element* object)
{
if(object->pos[0] >= child->vBoundriesBox[0][0] && object->pos[0] <= child->vBoundriesBox[1][0] &&
object->pos[1] >= child->vBoundriesBox[4][1] && object->pos[1] <= child->vBoundriesBox[0][1] &&
object->pos[2] >= child->vBoundriesBox[3][2] && object->pos[2] <= child->vBoundriesBox[0][2])
return true;
return false;
}


Because I am using pointers to move around tree elements I don't think object creation/destruction is an issue here. The one place I suppose might have impact on speed is this one:

Element* obj = this->objects.back();
this->objects.pop_back();
if(contains(trees[0][0][0],obj))


Although I am not sure how I can ommit/speed it up. Does someone has any suggestions what can be done here?

EDIT #1:

I've done some napkin math and I suppose there is one more place which might be causing major speed decrease. Boundries checking looks like is doing a lot and what I calculated is that the added complexity due to this is in worst case scenario: $$number\_of\_elements*number\_of\_childern*number\_of\_faces*MAX\_LEVELS$$

which in my case is equal to n*240.

Can someone please confirm if my idea is reasonable?

EDIT #2:

Attracted by my napkin math as well as looking at my code design I've decided to modify tree children generation and testing. In the matter of fields I've added those:

struct Bounds
{
glm::vec3 center;
glm::vec3 half_width;
};
typedef struct Bounds Bounds

Bounds bbox;


and as for the code it looks now like this:

void Quadtree::AddObject(Node* object)
{
this->objects.push_back(object);
}

{
if(this->objects.size()<=1 || level > 5)
{
for(Node* node:this->objects)
{
node->parent_group = this;
}
return this->objects.size();
}

if(leaf)
{
for(int i = 0; i<8; i++)
{
Bounds child_bounds = bbox;
child_bounds.center.x += bbox.half_width.x * ((i&4)? 0.5f:-0.5f);
child_bounds.center.y += bbox.half_width.y * ((i&2)? 0.5f:-0.5f);
child_bounds.center.z += bbox.half_width.z * ((i&1)? 0.5f:-0.5f);
child_bounds.half_width = bbox.half_width * 0.5f;
trees[i] = new Quadtree(this, child_bounds, level+1, maxLevel);
}
leaf = false;
}

Node* obj;
while (!this->objects.empty())
{
obj = this->objects.back();
this->objects.pop_back();
int idx = get_children_index(obj->pos);
trees[idx]->combined_weight += obj->weight;
}

int counter = 0;

for(int i=0;i<8;i++)
{
counter += trees[i]->Update();
}
return counter;
}

{
int idx = 0;
if(point.x >= bbox.center.x)
idx |= 1 << 2;
if(point.y >= bbox.center.y)
idx |= 1 << 1;
if(point.z >= bbox.center.z)
idx |= 1 << 0;
return idx;
}


Cleaner, nicer, but unfortunately - not faster very much. What am I still missing? Any idea?

To see what's the difference in comparison to initial Force-Directed graph I've run simple test with still this same dataset. Blue line is initial force-directed graph algorithm, black line is Barnes-Hut with old children creation, red line is current Barnes-Hut implementation. As you can see initial -> Barnes-Hut speed-up is nice, but still I feel it is way too slow for this algorithm.

• It looks like you've updated your code since posting the question. Please have one single piece of code for reviewers to review, and leave old versions to revision history or GitHub.
– anon
Commented May 6, 2016 at 19:41

My first reaction is your dynamic memory allocation:

            trees[i] = new Quadtree(this, child_bounds, level+1, maxLevel);


You are better off popping those objects from a pre-made stack the size of your maximum expected Quadtree object count. Allocating memory is very slow.

My second reaction is your use of std::vector. It may allocate/deallocate big blocks of memory when you push (and pop?). Your usage pattern resembles that of std::list which may serve you better, at least getting close to constant insertion/deletion properties. As for performance the std lib is the fastest there is at doing what it does, but you may want more control over the memory allocation/deallocation part. Do what you wish.

If your parameters are guaranteed to be non-NULL you might want to send them as references instead of pointers. It may or may not improve performance. The compiler may have it easier optimizing your code. And besides you rarely manage cases where pointers may be NULL.

Your contains function should be able to share calculations and results between nodes as the nodes share bounding planes with each other. As it is right now you test 6 planes in best case and 40 planes worst case. By sharing planes you reduce that to a constant 9 for all cases. Probably no big deal.

Always pre-increment loop counters. This will probably have no noticeable effect with decent compiler optimization but it is good coding practice.

A bit off-topic maybe but you may want to look at other platforms to implement the algorithm on such as OpenCL, even if targeting CPU only environments.