quadtrees

I ran into performance issues within some of my simulations. The agents and particles in these simulations tend to do many neighbor comparisons and in the worst case requiring $N^2$ checks. Intuitively, each particle does not need to know information about every other particle to influence its own behavior, with this in mind we can reduce the number of checks a particle has to make. Many simulations apply state changes to colliding particles, others apply constant state changes based on neighboring particles. With the use of a Quadtree we can bring the number of checks down from $N^2$ to $N * log_2(N)$, and in practical terms increase the performance of a system with 1000 particles by more than a 100 times. So what is a Quadtree?

A QuadTree is a tree like data structure that creates quadrilateral regions in the simulating space, where each region can recursively embed another set of quadrilateral regions. How does this help increase performance? It helps by drastically reducing the neighbor search space. Take for example a simulation with 10 particles where the 3 closest neighbors to a particle dictates its behavior. Instead of calculating the 3 closest neighbors by making checks to all particles, we can find the 3 closest neighbors by only making checks to the particles in the same quadrilateral region.

Quadtrees have two primary functions. You can insert into a Quadtree and you can query a Quadtree. Insertion puts the particle into the correct quad region. Each region has a set number of particles it can contain (I tend to use 8 or 16). Once a quad region has reached its max capacity the region subdivides and creates 4 new child regions which increases the max capacity of the parent region. Querying a quad tree is quite simple. It's a traversal from the root quad tree where subtrees that don't intersect with the queried range are skipped providing performance of $log_2(N) + K$ per query where K is the number of points in the region.

Example implementation:

class quadtree {
public:
    const static int QT_CAPACITY = 4;
    quadtree(point2d center, int halfDimension)
        : halfDimension(halfDimension), boundary(center.x, center.y, halfDimension) {}
    ~quadtree() {
        if (tl) delete tl;
        if (tr) delete tr;
        if (bl) delete bl;
        if (br) delete br;
    }

    bool insert(point2d point);
    void subdivide();
    std::vector<point2d> query(rectangle range);


private:
    int halfDimension;
    rectangle boundary;
    std::vector<point2d> points;

    quadtree* tl = nullptr;
    quadtree* tr = nullptr;
    quadtree* bl = nullptr;
    quadtree* br = nullptr;
};

void quadtree::subdivide() {
    this->tl =
        new quadtree({boundary.center.x - halfDimension / 2, boundary.center.y - halfDimension / 2},
                     halfDimension / 2);
    this->tr =
        new quadtree({boundary.center.x + halfDimension / 2, boundary.center.y - halfDimension / 2},
                     halfDimension / 2);
    this->bl =
        new quadtree({boundary.center.x - halfDimension / 2, boundary.center.y + halfDimension / 2},
                     halfDimension / 2);
    this->br =
        new quadtree({boundary.center.x + halfDimension / 2, boundary.center.y + halfDimension / 2},
                     halfDimension / 2);
}

bool quadtree::insert(point2d point) {
    // if the point is not within the boundary don't even try to insert it
    if (!boundary.containsPoint(point)) {
        return false;
    }

    // if there is space left in the quad tree and if it doesn't have subdivision then add it here
    if (points.size() < QT_CAPACITY && tl == nullptr) {
        points.push_back(point);
        return true;
    }

    // other wise subdivide and check children
    if (tl == nullptr) {
        subdivide();
    }

    if (tl->insert(point))
        return true;
    if (tr->insert(point))
        return true;
    if (bl->insert(point))
        return true;
    if (br->insert(point))
        return true;

    // never will get here
    return false;
}

std::vector<point2d> quadtree::query(rectangle range) {
    std::vector<point2d> result;
    if (!boundary.intersects(range)) {
        return result;
    }

    for (int i = 0; i < points.size(); i++) {
        point2d point = points[i];
        if (range.containsPoint(point)) {
            result.push_back(point);
        }
    }

    // if there are no children, terminate
    if (tl == nullptr) {
        return result;
    }

    // otherwise query the children as well and append these points to result
    std::vector<point2d> query = tl->query(range);
    result.insert(result.end(), query.begin(), query.end());

    query = tr->query(range);
    result.insert(result.end(), query.begin(), query.end());

    query = bl->query(range);
    result.insert(result.end(), query.begin(), query.end());

    query = br->query(range);
    result.insert(result.end(), query.begin(), query.end());

    return result;
}

class rectangle {
public:
    point2d center;
    int halfDimension;

    rectangle(int _x, int _y, int _halfDimension)
        : center(_x, _y), halfDimension(_halfDimension) {};
    ~rectangle() = default;

    bool containsPoint(point2d point) {
        const int dist = sqrt(pow(point.x - center.x, 2) + pow(point.y - center.y, 2));

        return dist <= sqrt(pow(halfDimension, 2) + pow(halfDimension, 2));
    }

    bool intersects(rectangle other) {
        const int dist =
            sqrt(pow(center.x - other.center.x, 2) + pow(center.y - other.center.y, 2));

        return dist <= halfDimension + other.halfDimension;
        return false;
    }
};


class point2d {
public:
    int x, y;
    point2d(int x, int y) : x(x), y(y) {}
    ~point2d() {}

};