I am trying to implement quadtree for my Barnes-Hut algorithm implementation. I am not sure that the code I've wrote so far is good implementation - it is dirty, and tend to be slow.

Of course I am also aware that the results might be correct, but I'd be glad if someone can take a look at that and give me a hints. Because the whole code is very big and also includes OpenGL visualization, I am putting here only important parts. I hope that's ok.

Quadtree.h

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

Quadtree.cpp

#define MAX_LEVELS 5

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

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

    if(leaf)
    {
        glm::vec3 newBoundries[8];

        // 0,0,0
        newBoundries[0] = vBoundriesBox[0];
        newBoundries[1] = glm::vec3(vBoundriesBox[0].x + (std::abs(vBoundriesBox[1].x - vBoundriesBox[0].x)/2.0f),vBoundriesBox[0].y,vBoundriesBox[0].z );
        newBoundries[2] = glm::vec3(vBoundriesBox[0].x + (std::abs(vBoundriesBox[1].x - vBoundriesBox[0].x)/2.0f),vBoundriesBox[0].y,vBoundriesBox[0].z - (std::abs(vBoundriesBox[3].z - vBoundriesBox[0].z)/2.0f));
        newBoundries[3] = glm::vec3(vBoundriesBox[0].x,vBoundriesBox[0].y,vBoundriesBox[0].z - (std::abs(vBoundriesBox[3].z - vBoundriesBox[0].z)/2.0f));

        newBoundries[4] = glm::vec3(vBoundriesBox[0].x,vBoundriesBox[0].y - (std::abs(vBoundriesBox[4].y - vBoundriesBox[0].y)/2.0f), vBoundriesBox[0].z);
        newBoundries[5] = glm::vec3(vBoundriesBox[0].x + (std::abs(vBoundriesBox[1].x - vBoundriesBox[0].x)/2.0f),vBoundriesBox[0].y - (std::abs(vBoundriesBox[4].y - vBoundriesBox[0].y)/2.0f),vBoundriesBox[0].z);
        newBoundries[6] = glm::vec3(vBoundriesBox[0].x + (std::abs(vBoundriesBox[1].x - vBoundriesBox[0].x)/2.0f),vBoundriesBox[0].y - (std::abs(vBoundriesBox[4].y - vBoundriesBox[0].y)/2.0f),vBoundriesBox[0].z - (std::abs(vBoundriesBox[3].z - vBoundriesBox[0].z)/2.0f));
        newBoundries[7] = glm::vec3(vBoundriesBox[0].x,vBoundriesBox[0].y - (std::abs(vBoundriesBox[4].y - vBoundriesBox[0].y)/2.0f),vBoundriesBox[0].z - (std::abs(vBoundriesBox[3].z - vBoundriesBox[0].z)/2.0f));

        trees[0][0][0] = new Quadtree(this, newBoundries, level+1);

        // 0,1,0
        newBoundries[0] = glm::vec3(vBoundriesBox[1].x - (std::abs(vBoundriesBox[1].x - vBoundriesBox[0].x)/2.0f),vBoundriesBox[0].y,vBoundriesBox[1][2]);
        newBoundries[1] = vBoundriesBox[1];
        newBoundries[2] = glm::vec3(vBoundriesBox[1].x,vBoundriesBox[1].y,vBoundriesBox[1].z - (std::abs(vBoundriesBox[2][2] - vBoundriesBox[1][2])/2.0f));
        newBoundries[3] = glm::vec3(vBoundriesBox[1].x - (std::abs(vBoundriesBox[1].x - vBoundriesBox[0].x)/2.0f),vBoundriesBox[1][1],vBoundriesBox[1][2] - (std::abs(vBoundriesBox[2].z - vBoundriesBox[1].z)/2.0f));

        newBoundries[4] = glm::vec3(vBoundriesBox[1].x - (std::abs(vBoundriesBox[1].x - vBoundriesBox[0].x)/2.0f),vBoundriesBox[1].y - (std::abs(vBoundriesBox[5].y - vBoundriesBox[1].y)/2.0f),vBoundriesBox[1].z);
        newBoundries[5] = glm::vec3(vBoundriesBox[1].x,vBoundriesBox[1].y - (std::abs(vBoundriesBox[5].y - vBoundriesBox[1].y)/2.0f), vBoundriesBox[1].z);
        newBoundries[6] = glm::vec3(vBoundriesBox[1].x,vBoundriesBox[1].y - (std::abs(vBoundriesBox[5].y - vBoundriesBox[1].y)/2.0f),vBoundriesBox[1].z - (std::abs(vBoundriesBox[2].z - vBoundriesBox[1].z)/2.0f));
        newBoundries[7] = glm::vec3(vBoundriesBox[1].x - (std::abs(vBoundriesBox[1].x - vBoundriesBox[0].x)/2.0f),vBoundriesBox[1].y - (std::abs(vBoundriesBox[5].y - vBoundriesBox[1].y)/2.0f),vBoundriesBox[1].z - (std::abs(vBoundriesBox[2].z - vBoundriesBox[1].z)/2.0f));

        trees[0][1][0] = new Quadtree(this, newBoundries, level+1);

        // 0,0,1
        newBoundries[0] = glm::vec3(vBoundriesBox[3].x,vBoundriesBox[3].y,vBoundriesBox[3].z + (std::abs(vBoundriesBox[3].z - vBoundriesBox[0].z)/2.0f));
        newBoundries[1] = glm::vec3(vBoundriesBox[3].x + (std::abs(vBoundriesBox[2].x - vBoundriesBox[3].x)/2.0f),vBoundriesBox[3].y,vBoundriesBox[3].z + (std::abs(vBoundriesBox[3].z - vBoundriesBox[0].z)/2.0f));
        newBoundries[2] = glm::vec3(vBoundriesBox[3].x + (std::abs(vBoundriesBox[2].x - vBoundriesBox[3].x)/2.0f),vBoundriesBox[3].y,vBoundriesBox[3].z );
        newBoundries[3] = vBoundriesBox[3];

        newBoundries[4] = glm::vec3(vBoundriesBox[3].x,vBoundriesBox[3].y - (std::abs(vBoundriesBox[7].y - vBoundriesBox[3].y)/2.0f),vBoundriesBox[3].z + (std::abs(vBoundriesBox[3].z - vBoundriesBox[0].z)/2.0f));
        newBoundries[5] = glm::vec3(vBoundriesBox[3].x + (std::abs(vBoundriesBox[2].x - vBoundriesBox[3].x)/2.0f),vBoundriesBox[3].y - (std::abs(vBoundriesBox[7].y - vBoundriesBox[3].y)/2.0f),vBoundriesBox[3].z + (std::abs(vBoundriesBox[3].z - vBoundriesBox[0].z)/2.0f));
        newBoundries[6] = glm::vec3(vBoundriesBox[3].x + (std::abs(vBoundriesBox[2].x - vBoundriesBox[3].x)/2.0f),vBoundriesBox[3].y - (std::abs(vBoundriesBox[7].y - vBoundriesBox[3].y)/2.0f),vBoundriesBox[3].z );
        newBoundries[7] = glm::vec3(vBoundriesBox[3].x,vBoundriesBox[3].y - (std::abs(vBoundriesBox[7].y - vBoundriesBox[3].y)/2.0f),vBoundriesBox[3].z);

        trees[0][0][1] = new Quadtree(this, newBoundries, level+1);

        // 0,1,1
        newBoundries[0] = glm::vec3(vBoundriesBox[2].x - (std::abs(vBoundriesBox[2].x - vBoundriesBox[3].x)/2.0f),vBoundriesBox[2].y,vBoundriesBox[2].z + (std::abs(vBoundriesBox[2][2] - vBoundriesBox[1][2])/2.0f));
        newBoundries[1] = glm::vec3(vBoundriesBox[2].x,vBoundriesBox[2].y,vBoundriesBox[2].z + (std::abs(vBoundriesBox[2].z - vBoundriesBox[1].z)/2.0f));
        newBoundries[2] = vBoundriesBox[2];
        newBoundries[3] = glm::vec3(vBoundriesBox[2].x - (std::abs(vBoundriesBox[2].x - vBoundriesBox[3].x)/2.0f),vBoundriesBox[2].y,vBoundriesBox[2].z );

        newBoundries[4] = glm::vec3(vBoundriesBox[2].x - (std::abs(vBoundriesBox[2].x - vBoundriesBox[3].x)/2.0f),vBoundriesBox[2].y - (std::abs(vBoundriesBox[6].y - vBoundriesBox[2].y)/2.0f),vBoundriesBox[2].z + (std::abs(vBoundriesBox[2].z - vBoundriesBox[1].z)/2.0f));
        newBoundries[5] = glm::vec3(vBoundriesBox[2].x,vBoundriesBox[2].y - (std::abs(vBoundriesBox[6].y - vBoundriesBox[2].y)/2.0f),vBoundriesBox[2].z + (std::abs(vBoundriesBox[2].z - vBoundriesBox[1].z)/2.0f));
        newBoundries[6] = glm::vec3(vBoundriesBox[2].x,vBoundriesBox[2].y - (std::abs(vBoundriesBox[6].y - vBoundriesBox[2].y)/2.0f),vBoundriesBox[2].z);
        newBoundries[7] = glm::vec3(vBoundriesBox[2].x - (std::abs(vBoundriesBox[2].x - vBoundriesBox[3].x)/2.0f),vBoundriesBox[2].y - (std::abs(vBoundriesBox[6].y - vBoundriesBox[2].y)/2.0f),vBoundriesBox[2].z);

        trees[0][1][1] = new Quadtree(this, newBoundries, level+1);

        // 1,0,0
        newBoundries[0] = glm::vec3(vBoundriesBox[4].x,vBoundriesBox[4].y + (std::abs(vBoundriesBox[4].y - vBoundriesBox[0].y)/2.0f), vBoundriesBox[4].z);
        newBoundries[1] = glm::vec3(vBoundriesBox[4].x + (std::abs(vBoundriesBox[5].x - vBoundriesBox[4].x)/2.0f),vBoundriesBox[4].y + (std::abs(vBoundriesBox[4].y - vBoundriesBox[0].y)/2.0f),vBoundriesBox[4].z);
        newBoundries[2] = glm::vec3(vBoundriesBox[4].x + (std::abs(vBoundriesBox[5].x - vBoundriesBox[4].x)/2.0f),vBoundriesBox[4].y + (std::abs(vBoundriesBox[4].y - vBoundriesBox[0].y)/2.0f),vBoundriesBox[4].z - (std::abs(vBoundriesBox[7].z - vBoundriesBox[4].z)/2.0f));
        newBoundries[3] = glm::vec3(vBoundriesBox[4].x,vBoundriesBox[4].y + (std::abs(vBoundriesBox[4].y - vBoundriesBox[0].y)/2.0f),vBoundriesBox[4].z - (std::abs(vBoundriesBox[7].z - vBoundriesBox[4].z)/2.0f));

        newBoundries[4] = vBoundriesBox[4];
        newBoundries[5] = glm::vec3(vBoundriesBox[4].x + (std::abs(vBoundriesBox[5].x - vBoundriesBox[4].x)/2.0f),vBoundriesBox[4].y,vBoundriesBox[4].z );
        newBoundries[6] = glm::vec3(vBoundriesBox[4].x + (std::abs(vBoundriesBox[5].x - vBoundriesBox[4].x)/2.0f),vBoundriesBox[4].y,vBoundriesBox[4].z - (std::abs(vBoundriesBox[7].z - vBoundriesBox[4].z)/2.0f));
        newBoundries[7] = glm::vec3(vBoundriesBox[4].x,vBoundriesBox[4].y,vBoundriesBox[4].z - (std::abs(vBoundriesBox[7].z - vBoundriesBox[4].z)/2.0f));

        trees[1][0][0] = new Quadtree(this, newBoundries, level+1);

        // 1,1,0
        newBoundries[0] = glm::vec3(vBoundriesBox[5].x - (std::abs(vBoundriesBox[4].x - vBoundriesBox[5].x)/2.0f),vBoundriesBox[5].y + (std::abs(vBoundriesBox[5].y - vBoundriesBox[1].y)/2.0f),vBoundriesBox[5].z);
        newBoundries[1] = glm::vec3(vBoundriesBox[5].x,vBoundriesBox[5].y + (std::abs(vBoundriesBox[5].y - vBoundriesBox[1].y)/2.0f),vBoundriesBox[5].z);
        newBoundries[2] = glm::vec3(vBoundriesBox[5].x,vBoundriesBox[5].y + (std::abs(vBoundriesBox[5].y - vBoundriesBox[1].y)/2.0f),vBoundriesBox[5].z - (std::abs(vBoundriesBox[5].z - vBoundriesBox[6].z)/2.0f));
        newBoundries[3] = glm::vec3(vBoundriesBox[5].x - (std::abs(vBoundriesBox[4].x - vBoundriesBox[5].x)/2.0f),vBoundriesBox[5].y + (std::abs(vBoundriesBox[5].y - vBoundriesBox[1].y)/2.0f),vBoundriesBox[5].z - (std::abs(vBoundriesBox[5].z - vBoundriesBox[6].z)/2.0f));

        newBoundries[4] = glm::vec3(vBoundriesBox[5].x - (std::abs(vBoundriesBox[4].x - vBoundriesBox[5].x)/2.0f),vBoundriesBox[5].y,vBoundriesBox[5].z);
        newBoundries[5] = vBoundriesBox[5];
        newBoundries[6] = glm::vec3(vBoundriesBox[5].x,vBoundriesBox[5].y,vBoundriesBox[5].z - (std::abs(vBoundriesBox[6].z - vBoundriesBox[5].z)/2.0f));
        newBoundries[7] = glm::vec3(vBoundriesBox[5].x - (std::abs(vBoundriesBox[4].x - vBoundriesBox[5].x)/2.0f),vBoundriesBox[5].y,vBoundriesBox[5].z - (std::abs(vBoundriesBox[5].z - vBoundriesBox[6].z)/2.0f));

        trees[1][1][0] = new Quadtree(this, newBoundries, level+1);

        // 1,0,1
        newBoundries[0] = glm::vec3(vBoundriesBox[7].x,vBoundriesBox[7].y + (std::abs(vBoundriesBox[7].y - vBoundriesBox[3].y)/2.0f),vBoundriesBox[7].z + (std::abs(vBoundriesBox[7].z - vBoundriesBox[4].z)/2.0f));
        newBoundries[1] = glm::vec3(vBoundriesBox[7].x + (std::abs(vBoundriesBox[7].x - vBoundriesBox[6].x)/2.0f),vBoundriesBox[7].y + (std::abs(vBoundriesBox[7].y - vBoundriesBox[3].y)/2.0f),vBoundriesBox[7].z + (std::abs(vBoundriesBox[4].z - vBoundriesBox[7].z)/2.0f));
        newBoundries[2] = glm::vec3(vBoundriesBox[7].x + (std::abs(vBoundriesBox[7].x - vBoundriesBox[6].x)/2.0f),vBoundriesBox[7].y + (std::abs(vBoundriesBox[7].y - vBoundriesBox[3].y)/2.0f),vBoundriesBox[7].z );
        newBoundries[3] = glm::vec3(vBoundriesBox[7].x,vBoundriesBox[7].y + (std::abs(vBoundriesBox[7].y - vBoundriesBox[3].y)/2.0f),vBoundriesBox[7].z );

        newBoundries[4] = glm::vec3(vBoundriesBox[7].x,vBoundriesBox[7].y,vBoundriesBox[7][2] + (std::abs(vBoundriesBox[7].z - vBoundriesBox[4].z)/2.0f));
        newBoundries[5] = glm::vec3(vBoundriesBox[7].x + (std::abs(vBoundriesBox[7].x - vBoundriesBox[6].x)/2.0f),vBoundriesBox[7].y,vBoundriesBox[7].z + (std::abs(vBoundriesBox[7].z - vBoundriesBox[4].z)/2.0f));
        newBoundries[6] = glm::vec3(vBoundriesBox[7].x + (std::abs(vBoundriesBox[7].x - vBoundriesBox[6].x)/2.0f),vBoundriesBox[7].y,vBoundriesBox[7].z );
        newBoundries[7] = vBoundriesBox[7];

        trees[1][0][1] = new Quadtree(this, newBoundries, level+1);

        // 1,1,1
        newBoundries[0] = glm::vec3(vBoundriesBox[6].x - (std::abs(vBoundriesBox[6].x - vBoundriesBox[7].x)/2.0f),vBoundriesBox[6].y + (std::abs(vBoundriesBox[6].y - vBoundriesBox[2].y)/2.0f),vBoundriesBox[6].z + (std::abs(vBoundriesBox[6].z - vBoundriesBox[5].z)/2.0f));
        newBoundries[1] = glm::vec3(vBoundriesBox[6].x,vBoundriesBox[6].y + (std::abs(vBoundriesBox[6].y - vBoundriesBox[2].y)/2.0f),vBoundriesBox[6].z + (std::abs(vBoundriesBox[6].z - vBoundriesBox[5].z)/2.0f));
        newBoundries[2] = glm::vec3(vBoundriesBox[6].x,vBoundriesBox[6].y + (std::abs(vBoundriesBox[6].y - vBoundriesBox[2].y)/2.0f),vBoundriesBox[6].z);
        newBoundries[3] = glm::vec3(vBoundriesBox[6].x - (std::abs(vBoundriesBox[6].x - vBoundriesBox[7].x)/2.0f),vBoundriesBox[6].y + (std::abs(vBoundriesBox[6].y - vBoundriesBox[2].y)/2.0f),vBoundriesBox[6].z);

        newBoundries[4] = glm::vec3(vBoundriesBox[6].x - (std::abs(vBoundriesBox[6].x - vBoundriesBox[7].x)/2.0f),vBoundriesBox[6].y,vBoundriesBox[6].z + (std::abs(vBoundriesBox[6].z - vBoundriesBox[5].z)/2.0f));
        newBoundries[5] = glm::vec3(vBoundriesBox[6].x,vBoundriesBox[6].y,vBoundriesBox[6].z + (std::abs(vBoundriesBox[6].z - vBoundriesBox[5].z)/2.0f));
        newBoundries[6] = vBoundriesBox[6];
        newBoundries[7] = glm::vec3(vBoundriesBox[6].x - (std::abs(vBoundriesBox[6].x - vBoundriesBox[7].x)/2.0f),vBoundriesBox[6].y,vBoundriesBox[6].z );

        trees[1][1][1] = new Quadtree(this, newBoundries, level+1);
        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]->AddObject(obj);
            trees[0][0][0]->combined_weight += obj->weight;
        } else if(contains(trees[0][0][1],obj))
        {
            trees[0][0][1]->AddObject(obj);
            trees[0][0][1]->combined_weight += obj->weight;
        } else if(contains(trees[0][1][0],obj))
        {
            trees[0][1][0]->AddObject(obj);
            trees[0][1][0]->combined_weight += obj->weight;
        } else if(contains(trees[0][1][1],obj))
        {
            trees[0][1][1]->AddObject(obj);
            trees[0][1][1]->combined_weight += obj->weight;
        } else if(contains(trees[1][0][0],obj))
        {
            trees[1][0][0]->AddObject(obj);
            trees[1][0][0]->combined_weight += obj->weight;
        } else if(contains(trees[1][0][1],obj))
        {
            trees[1][0][1]->AddObject(obj);
            trees[1][0][1]->combined_weight += obj->weight;
        } else if(contains(trees[1][1][0],obj))
        {
            trees[1][1][0]->AddObject(obj);
            trees[1][1][0]->combined_weight += obj->weight;
        } else if(contains(trees[1][1][1],obj))
        {
            trees[1][1][1]->AddObject(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 Quadtree::contains(Quadtree* 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;
}

and in my main file I am adding elements into tree like this:

void AddTreeElements(Quadtree* 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])
        {
            tree->AddObject(&node.second);
            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;
        }
    }
}

Because it is octree I am using box to visualize it and test if it's working correctly. Also, boundries array represents this box and this is how I am numbering corners: Box As for the speed I've run sample test by randomly placing new nodes and this is how time/elements plot looks like: enter image description here Is it correct behavior when I am changing MAX_LEVELS value?

Superficial review

  • Naming: If this is an octree, let's call it Octree rather than Quadtree!
  • Naming: Pick a naming convention and stick with it. We have a mix of camel-case, snake-case and pascal-case for members and variables. A particularly confusing name is in for(Element* Element:this->objects), where a variable has the same name as a type.
  • Spelling: "Boundary" and "boundaries" both need an 'a'.
  • Prefer double to float, unless you have a demonstrated problem of storage (in this case, we might need Element to contain float members, but the octree can use a higher precision).
  • Prefer strongly-typed constants rather than macro definitions (although I think MAX_LEVELS shouldn't be a constant at all - see below).
  • The contains() member is a pure lookup, so should be declared const.
  • Avoid owning objects by raw pointers; instead, use smart pointers with automatic memory management.

Missing code

I had to infer the Quadtree and Element class definitions (so I might have come up with something different to what you have):

#include <glm/glm.hpp>

#include <vector>

class Quadtree;

struct Element
{
    glm::vec3 pos;
    float weight;
    Quadtree* parent_group;
};

class Quadtree
{
    unsigned level = 0;
    Quadtree* trees[2][2][2];
    glm::vec3 vBoundriesBox[8];
    bool leaf;
    float combined_weight = 0;
    std::vector<Element*> objects;

public:
    Quadtree(Quadtree* parent,
             const glm::vec3(&)[8],
             unsigned level);
    Quadtree(const Quadtree&) = delete;
    void operator=(const Quadtree&) = delete;

    void AddObject(Element* object);
    void Update();
    bool contains(Quadtree* child, Element* object);
};

And also the constructor:

#include <algorithm>

Quadtree::Quadtree(Quadtree*,
                   const glm::vec3 (&children)[8],
                   unsigned level)
    : level{level}, leaf{true}, objects{}
{
    using std::begin;
    using std::end;
    std::copy(begin(children), end(children),
              begin(vBoundriesBox));
}

Missing tests

This would have been much easier to review with a working test set, using either a C++ test framework or just a simple main().

Re-write

I'm going to start with a similar Element structure, but I don't see a need for an element to know its container here, so I'll remove the intrusive pointer:

struct Element
{
    glm::vec3 pos;
    float weight;
};

Now for the octree. Start with its bounding box, which can be described using two vectors (I'm assuming the usual mathematical right-handed coordinates):

glm::vec3 near_bottom_left;
glm::vec3 far_top_right;

As these don't change, they can be declared const. It's also useful to know the midpoint, so we can allocate points to their quadrants.

We should also have a level member, but I'll change its meaning so that it's then maximum depth permitted below this level. That means we no longer need MAX_LEVELS to tell us how deep the tree can go - it's simply specified as a property of the top node and propagated down until it reaches zero.

We'll want to record the aggregate properties, but I won't actually implement the updating of these as you seem to be comfortable with that part.

double total_weight;
glm::vec3 mean_position;

Finally, to make it an actual tree, we need some children:

Octree* children[8];
std::vector<Element> objects;

For the user interface, we'll need a constructor and some way to add an element:

public:
    explicit Octree(glm::vec3 far_bottom_left,
                    glm::vec3 near_top_right,
                    unsigned max_depth);

    void add(Element e);

The constructor is easy to implement:

Octree::Octree(glm::vec3 far_bottom_left,
               glm::vec3 near_top_right,
               unsigned max_depth)
    : far_bottom_left{far_bottom_left},
      near_top_right{near_top_right},
      middle{0.5f * far_bottom_left + 0.5f * near_top_right},
      level{max_depth}
{}

There's a bit more work for add(), but we can start with the simplest case:

void Octree::add(Element e)
{
    if (objects.empty() && is_leaf()) {
        objects.emplace_back(std::move(e));
        return;
    }
}

Another simple case is at the final level:

    if (level == 0) {
        objects.emplace_back(std::move(e));
        return;
    }

I've assumed an is_leaf method; we'd better define that:

bool Octree::is_leaf() const
{
    // we're a leaf node if all children are null
    return std::all_of(std::begin(children), std::end(children),
                       [](auto& x){return !x;});
}

Now, the more difficult case - delegate to the appropriate child node:

    const bool left = e.pos.x < middle.x;
    const bool down = e.pos.y < middle.y;
    const bool far = e.pos.z < middle.z;
    auto& child = children[4*left + 2*down + far];
    if (!child) {
        // TODO: create a new node
    }
    child.add(std::move(e));

The TODO comment needs filling in. How do we determine the bounds of the new child node? We can start with the bounds of the parent, and move either the left or the right bound to the middle so that the point is still contained, then repeat in the other two dimensions:

        glm::vec3 fbl = far_bottom_left;
        glm::vec3 ntr = near_top_right;
        (left ? fbl : ntr).x = middle.x;
        (down ? fbl : ntr).y = middle.y;
        (far ? fbl : ntr).z = middle.z;
        child = std::make_unique<Octree>(fbl, ntr, level-1);
        // move any directly-held object
        for (auto& o: objects)
            add(std::move(o));
        objects.clear();

Putting it all together

There were a couple of bugs in my thoughts above; having corrected them, we get:

#include <glm/glm.hpp>

#include <algorithm>
#include <iterator>
#include <iosfwd>
#include <memory>
#include <string>
#include <vector>

struct Element
{
    glm::vec3 pos;
    float weight;
};

class Octree
{
    // geometry
    const glm::vec3 far_bottom_left;
    const glm::vec3 near_top_right;
    const glm::vec3 middle;

    // limit to tree depth
    unsigned level;

    // aggregate properties from constituent objects
    double total_weight = 0;
    glm::vec3 mean_position = {};

    // subtree objects - any or all may be null
    // index using (left=4, right=0) + (bottom=2, top=0) + (far=1, near=0)
    std::unique_ptr<Octree> children[8] = {};
    // directly held elements - empty unless all children are null
    std::vector<Element> objects = {};

public:
    explicit Octree(glm::vec3 far_bottom_left,
                    glm::vec3 near_top_right,
                    unsigned max_depth);

    void add(Element e);

    std::ostream& print(std::ostream&, const std::string& indent) const;

    friend std::ostream& operator<<(std::ostream& out, const Octree& o) {
        return o.print(out, {});
    }

private:
    bool is_leaf() const;
};


Octree::Octree(glm::vec3 far_bottom_left,
               glm::vec3 near_top_right,
               unsigned max_depth)
    : far_bottom_left{far_bottom_left},
      near_top_right{near_top_right},
      middle{0.5f * far_bottom_left + 0.5f * near_top_right},
      level{max_depth}
{}


void Octree::add(Element e)
{
    if (level == 0 || objects.empty() && is_leaf()) {
        objects.emplace_back(std::move(e));
        return;
    }
    const bool left = e.pos.x < middle.x;
    const bool down = e.pos.y < middle.y;
    const bool far = e.pos.z < middle.z;
    auto& child = children[4*left + 2*down + far];
    if (!child) {
        glm::vec3 fbl = far_bottom_left;
        glm::vec3 ntr = near_top_right;
        (left ? ntr : fbl).x = middle.x;
        (down ? ntr : fbl).y = middle.y;
        (far ? ntr : fbl).z = middle.z;
        child = std::make_unique<Octree>(fbl, ntr, level-1);
        // move any directly-held object
        auto to_move = std::move(objects);
        objects.clear();
        for (auto& o: to_move)
            add(std::move(o));
    }
    child->add(std::move(e));
}


bool Octree::is_leaf() const
{
    // we're a leaf node if all children are null
    return std::all_of(std::begin(children), std::end(children),
                       [](auto& x){return !x;});
}

To test it, here's some output operators:

std::ostream& operator<<(std::ostream& os, const glm::vec3& pos)
{
    return os << '(' << pos.x
              << ',' << pos.y
              << ',' << pos.z
              << ')';
}

std::ostream& operator<<(std::ostream& os, const Element& e)
{
    return os << e.pos;
}

std::ostream& Octree::print(std::ostream& os, const std::string& indent) const
{
    os << indent
       << "Octree " << far_bottom_left
       << " to " << near_top_right << '\n';
    if (!os) { return os; }
    for (auto const& e: objects) {
        os << indent << "- " << e << '\n';
    }
    if (!os) { return os; }
    for (auto i = 0u;  i < std::size(children);  ++i) {
        if (children[i]) {
            os << indent
               << (i & 1 ? "far" : "near") << ' '
               << (i & 2 ? "bottom" : "top") << ' '
               << (i & 4 ? "left" : "right") << ": \n";
            children[i]->print(os, indent+"  ");
            if (!os) { return os; }
        }
    }
    return os;
}

And a test program to show we're producing the right output, and to measure execution time:

#include <chrono>
#include <iostream>
#include <random>

int main()
{
    {
        // Visual test - do these few points go into the right boxes?
        Octree tree{{0,0,0}, {1,1,1}, 3};

        std::mt19937 gen;
        std::uniform_real_distribution<float> dis{0,1};

        for (auto i = 0u;  i < 20;  ++i) {
            tree.add({{dis(gen), dis(gen), dis(gen)}, 1});
        }

        std::cout << tree;
    }

    std::cout << std::endl;

    {
        // Timing test
        Octree tree{{0,0,0}, {1,1,1}, 5};

        std::mt19937 gen;
        std::uniform_real_distribution<float> dis{0,1};

        auto const starttime = std::chrono::high_resolution_clock::now();

        auto const chunk = 10000;
        for (auto i = 0u;  i < 20;  ++i) {
            for (auto j = 0;  j < chunk;  ++j) {
                tree.add({{dis(gen), dis(gen), dis(gen)}, 1});
            }
            auto const now = std::chrono::high_resolution_clock::now();
            auto const time_in_ms
                 = duration_cast<milliseconds>(now - starttime);
            std::cout << chunk*(i+1) << " points: "
                      << time_in_ms.count() << " ms\n";
        }
    }
}

My results:

Octree (0,0,0) to (1,1,1)
near top right: 
  Octree (0.5,0.5,0.5) to (1,1,1)
  near top right: 
    Octree (0.75,0.75,0.75) to (1,1,1)
    near top right: 
      Octree (0.875,0.875,0.875) to (1,1,1)
      - (0.992881,0.957507,0.996461)
    near top left: 
      Octree (0.75,0.875,0.875) to (0.875,1,1)
      - (0.792207,0.878431,0.959492)
  near top left: 
    Octree (0.5,0.75,0.75) to (0.75,1,1)
    - (0.725839,0.970593,0.98111)
  near bottom left: 
    Octree (0.5,0.5,0.75) to (0.75,0.75,1)
    - (0.503663,0.655741,0.797929)
far top right: 
  Octree (0.5,0.5,0) to (1,1,0.5)
  near top right: 
    Octree (0.75,0.75,0.25) to (1,1,0.5)
    - (0.798106,0.80028,0.297029)
  far top right: 
    Octree (0.75,0.75,0) to (1,1,0.25)
    - (0.964889,0.967695,0.157613)
  near bottom left: 
    Octree (0.5,0.5,0.25) to (0.75,0.75,0.5)
    - (0.740647,0.743132,0.474759)
near bottom right: 
  Octree (0.5,0,0.5) to (1,0.5,1)
  near bottom right: 
    Octree (0.75,0,0.75) to (1,0.25,1)
    near top left: 
      Octree (0.75,0.125,0.875) to (0.875,0.25,1)
      - (0.814724,0.135477,0.905792)
      - (0.835009,0.126987,0.968868)
  far bottom right: 
    Octree (0.75,0,0.5) to (1,0.25,0.75)
    - (0.913376,0.221034,0.632359)
  near top left: 
    Octree (0.5,0.25,0.75) to (0.75,0.5,1)
    - (0.678735,0.398739,0.75774)
far bottom right: 
  Octree (0.5,0,0) to (1,0.5,0.5)
  - (0.957167,0.109862,0.485376)
near top left: 
  Octree (0,0.5,0.5) to (0.5,1,1)
  far top left: 
    Octree (0,0.75,0.5) to (0.25,1,0.75)
    near top right: 
      Octree (0.125,0.875,0.625) to (0.25,1,0.75)
      - (0.211924,0.933993,0.68136)
    near top left: 
      Octree (0,0.875,0.625) to (0.125,1,0.75)
      - (0.112465,0.915736,0.639763)
near bottom left: 
  Octree (0,0,0.5) to (0.5,0.5,1)
  far top right: 
    Octree (0.25,0.25,0.5) to (0.5,0.5,0.75)
    - (0.392227,0.422088,0.655478)
  far bottom right: 
    Octree (0.25,0,0.5) to (0.5,0.25,0.75)
    far top left: 
      Octree (0.25,0.125,0.5) to (0.375,0.25,0.625)
      - (0.278498,0.188382,0.546881)
    far bottom left: 
      Octree (0.25,0,0.5) to (0.375,0.125,0.625)
      - (0.308167,0.0975404,0.547221)
  near top left: 
    Octree (0,0.25,0.75) to (0.25,0.5,1)
    - (0.0357117,0.361294,0.849129)
far bottom left: 
  Octree (0,0,0) to (0.5,0.5,0.5)
  near bottom left: 
    Octree (0,0,0.25) to (0.25,0.25,0.5)
    far top right: 
      Octree (0.125,0.125,0.25) to (0.25,0.25,0.375)
      - (0.173865,0.171187,0.301913)
    near bottom right: 
      Octree (0.125,0,0.375) to (0.25,0.125,0.5)
      - (0.141886,0.00478348,0.421761)
10000 points: 15 ms
20000 points: 27 ms
30000 points: 38 ms
40000 points: 49 ms
50000 points: 59 ms
60000 points: 69 ms
70000 points: 78 ms
80000 points: 88 ms
90000 points: 98 ms
100000 points: 107 ms
110000 points: 118 ms
120000 points: 129 ms
130000 points: 140 ms
140000 points: 152 ms
150000 points: 163 ms
160000 points: 172 ms
170000 points: 181 ms
180000 points: 190 ms
190000 points: 199 ms
200000 points: 208 ms

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.