# Spatial radius search in Boid simulation

I am working on a Boid simulation with 10,000 boids. I've successfully utilized geometry shaders and an array<GLfloat[3], MAX_BOIDS> for vertices to enhance display performance. Currently, I can display 500,000 boids moving randomly at 60 FPS.

However, the neighbor search is a bottleneck. With a naive search, I'm dealing with O(n²) complexity, which is prohibitively slow. Before experimenting with QuadTrees or RTrees, I wanted to benchmark a simple lattice structure. Unfortunately, this approach is already too slow; I can only manage to display 1,000 boids at 60 FPS.

Here the simulator. I can freely move the mouse to highlight seen nearby boids:

How can I optimize this code to improve its efficiency? Perhaps my data structure is not optimal.

Here is the code:

/**
* For a boid simulation, I am testing a bin lattice structure to speed up
* the neighbor search. The idea is to divide the screen into a grid of bins
* and store the boids in the corresponding bin. When querying for neighbors,
* only the bins in a 3x3 square around the boid need to be checked.
*
* Goal is to have a 60 FPS simulation with 10000 boids.
*/
#include <SFML/Graphics.hpp>
#include <cmath>
#include <array>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <span>

#define SHOW_GRID

#define WINDOW_WIDTH 1000
#define WINDOW_HEIGHT 1000

#define NUMBER_BINS 20 // Number of bins should be proportional to sight radius
#define MAX_BOIDS_PER_BIN 100 // Limit number of boids per bin to avoid large queries

#define BOIDS 1000 // To be increased to 10000
#define RADIUS 100 // Radius of the circle around the mouse to query for neighbors

struct Boid {
sf::Vector2f position;
};

/**
* In this implementation I am using a fixed size array for the boids in each bin.
* std::vector is not used because it is slower. I am seeking to optimize the
* neighbor search, so I am not concerned with the memory overhead of a fixed size
* array.
*/
template<size_t WIDTH, size_t HEIGHT, size_t BINS, size_t BIN_SIZE>
class BinLattice {
public:
static constexpr size_t MaxQuerySize = BIN_SIZE * 9;
using Bin = std::pair<std::array<Boid*, BIN_SIZE>, size_t>;
using Lattice = std::array<std::array<Bin, BINS>, BINS>;
using BoidArray = std::array<Boid*, MaxQuerySize>;
using QueryResult = std::pair<BoidArray, size_t>;  // Pair of array and count

BinLattice() {
for (auto &row : lattice) {
for (auto &bin : row) {
bin.first.fill(nullptr);
bin.second = 0;
}
}
}

void addBoid(Boid* boid) {
size_t xIndex = static_cast<size_t>(boid->position.x / (WIDTH / BINS));
size_t yIndex = static_cast<size_t>(boid->position.y / (HEIGHT / BINS));
auto& bin = lattice[xIndex][yIndex];
bin.first[bin.second++] = boid;
}

void clear() {
for (auto &row : lattice) {
for (auto &bin : row) {
bin.second = 0;
}
}
}

void query(sf::Vector2f &position, QueryResult& result, float radius) {
size_t count = 0;

size_t xIndex = static_cast<size_t>(position.x / (WIDTH / BINS));
size_t yIndex = static_cast<size_t>(position.y / (HEIGHT / BINS));

// Query all bins in a 3x3 square around the boid
for (int dx = -1; dx <= 1; ++dx) {
for (int dy = -1; dy <= 1; ++dy) {
size_t newX = std::clamp(xIndex + dx, 0ul, BINS - 1);
size_t newY = std::clamp(yIndex + dy, 0ul, BINS - 1);
auto& bin = lattice[newX][newY];
for (size_t i = 0; i < bin.second; ++i) {
// Check if boid is within radius
if (std::hypot(bin.first[i]->position.x - position.x,
bin.first[i]->position.y - position.y) < radius) {
result.first[count++] = bin.first[i];
}
}
}
}
result.second = count;
}

Lattice lattice;
};

int main() {
std::vector<Boid> boids;
using Lattice = BinLattice<WINDOW_WIDTH, WINDOW_HEIGHT, NUMBER_BINS, MAX_BOIDS_PER_BIN>;
for (int i = 0; i < BOIDS; ++i) {
boids.push_back({{
static_cast<float>(rand() % WINDOW_WIDTH),
static_cast<float>(rand() % WINDOW_HEIGHT)}});
}
Lattice binLattice;
Lattice::QueryResult result;

sf::RenderWindow window(sf::VideoMode(1000, 1000), "Boids");

// Load a font
sf::Font font;

// Setup text
sf::Text text;
text.setFont(font);
text.setCharacterSize(24);
text.setFillColor(sf::Color::White);
text.setPosition(10.f, 10.f);

// Grid lines
sf::RectangleShape hline(sf::Vector2f(WINDOW_WIDTH / NUMBER_BINS, 1.f));
sf::RectangleShape vline(sf::Vector2f(1.f, WINDOW_HEIGHT / NUMBER_BINS));
hline.setFillColor(sf::Color(255, 255, 255, 35));
vline.setFillColor(sf::Color(255, 255, 255, 35));

// Mouse circle
sf::CircleShape spotlight(1.f);
spotlight.setFillColor(sf::Color(255, 255, 255, 35));

// Boid circle
sf::CircleShape boidShape(1.f);
boidShape.setFillColor(sf::Color::Cyan);

// Highlight circle
sf::CircleShape boidSeen(2.f);
boidSeen.setFillColor(sf::Color::Yellow);

sf::Clock frameClock; // Clock to measure frame time
sf::Clock updateClock; // Clock to measure update intervals
float fps = 0.0f;
while (window.isOpen()) {
sf::Event event;
while (window.pollEvent(event)) { if (event.type == sf::Event::Closed) window.close(); }

binLattice.clear();
for (auto& boid : boids) binLattice.addBoid(&boid);

window.clear();

// Display lattice grid with clear gray lines
#ifdef SHOW_GRID

for (size_t i = 0; i < NUMBER_BINS; ++i) {
for (size_t j = 0; j < NUMBER_BINS; ++j) {
hline.setPosition(i * WINDOW_WIDTH / NUMBER_BINS, j * WINDOW_HEIGHT / NUMBER_BINS);
window.draw(hline);
vline.setPosition(i * WINDOW_WIDTH / NUMBER_BINS, j * WINDOW_HEIGHT / NUMBER_BINS);
window.draw(vline);
}
}
#endif

// Get Mouse Position
sf::Vector2i mousePosition = sf::Mouse::getPosition(window);
sf::Vector2f mousePositionFloat = sf::Vector2f(mousePosition.x, mousePosition.y);

// Draw clear alpha circle around mouse
window.draw(spotlight);

for (auto& boid : boids) {
boidShape.setPosition(boid.position);
window.draw(boidShape);
}

binLattice.clear();
for (auto& boid : boids) binLattice.addBoid(&boid);

#if 1 // For each boid (real usecase), very slow < 10 FPS
// *********** SLOW ***********
for (auto& boid : boids) {
for (size_t i = 0; i < result.second; ++i) {
boidSeen.setPosition(result.first[i]->position);
window.draw(boidSeen);
}
}
// *********** /SLOW ***********
#else // Only for mouse position (for profiling): very fast > 120 FPS
for (size_t i = 0; i < result.second; ++i) {
boidSeen.setPosition(result.first[i]->position);
window.draw(boidSeen);
}
#endif

sf::Time frameTime = frameClock.restart();
if (updateClock.getElapsedTime().asSeconds() >= 0.5) {
// Update FPS every 0.5 seconds
fps = 1.0f / frameTime.asSeconds();
updateClock.restart();
}

// Update text
std::stringstream ss;
ss << std::fixed << std::setprecision(2) << fps << " FPS";
text.setString(ss.str());
window.draw(text);

window.display();
}
}


Built on Ubuntu with:

g++ -Os -pg main.cpp -o myapp -lsfml-graphics -lsfml-window -lsfml-system


Curiously gprof doesn't show the bottleneck query method:

granularity: each sample hit covers 4 byte(s) for 25.00% of 0.04 seconds

index % time    self  children    called     name
<spontaneous>
[1]    100.0    0.04    0.00                 main [1]
0.00    0.00  230000/230000      BinLattice<1000ul, 1000ul, 20ul, 100ul>::addBoid(Boid*) [8]
0.00    0.00     230/230         BinLattice<1000ul, 1000ul, 20ul, 100ul>::clear() [9]
0.00    0.00     116/117         std::__cxx11::basic_string<unsigned int, std::char_traits<unsigned int>, std::allocator<unsigned int> >::_M_dispose() [10]
0.00    0.00      11/11          void std::vector<Boid, std::allocator<Boid> >::_M_realloc_insert<Boid>(__gnu_cxx::__normal_iterator<Boid*, std::vector<Boid, std::allocator<Boid> > >, Boid&&) [11]
0.00    0.00       3/3           sf::CircleShape::~CircleShape() [12]
0.00    0.00       2/2           sf::RectangleShape::~RectangleShape() [14]
0.00    0.00       1/1           sf::Text::~Text() [15]
0.00    0.00       1/1           std::_Vector_base<Boid, std::allocator<Boid> >::~_Vector_base() [16]
-----------------------------------------------
0.00    0.00  230000/230000      main [1]
[8]      0.0    0.00    0.00  230000         BinLattice<1000ul, 1000ul, 20ul, 100ul>::addBoid(Boid*) [8]
-----------------------------------------------
0.00    0.00     230/230         main [1]
[9]      0.0    0.00    0.00     230         BinLattice<1000ul, 1000ul, 20ul, 100ul>::clear() [9]

• I know the update didn't really incorporate feedback from the answers but per site policy I have rolled back Rev 2 → 1. Please see What should I do when someone answers my question?. Commented Nov 27, 2023 at 23:54
• @SᴀᴍOnᴇᴌᴀ would it be acceptable to post the updated version using boost on a self-answer? Commented Nov 27, 2023 at 23:56
• Please see point #2 in the numbered list below heading I improved my code based on the reviews. What next? on the link I posted, noting the first bullet point underneath it. Commented Nov 28, 2023 at 0:02
• is BinLattice a view into a vector of Boid if so name it as such for clarity and provide consructors similar to std::string_view. Otherwise use std::shared_ptr. Commented Nov 28, 2023 at 0:17
• Have you considerer compute shaders (or just multipass shaders) to do all or part of the simulation, 10K boids is nothing, 1000K+ boids at 60fps on modern GPU is the expectation. An example using WebGL & multipass shaders Random example shadertoy search "Boids" Commented Nov 29, 2023 at 11:48

Nice problem!

I gotta tell ya, QuadTrees is definitely the appropriate datastructure. And it's just not that hard. Let coord be the bit interleaving of x and y coords. That is, the top two bits are MSB of x and MSB of y. For the next two bits, keep shifting. Last two bits are LSB of x and y. I'm assuming htonl() network order. Now sort in the usual way, and locate a boid at index i. All the boids that are spatially near will also have nearby indexes.

Thank you for the nice manifest constants, they are very clear.

# comments where appropriate

        bin.first[bin.second++] = boid;


It wouldn't hurt to mention to the Gentle reader that .first is a boid and .second is count of boids within a bin. (Similar to the QueryCount comment.)

I see why it isn't necessary, but I'm still a little sad that clear() leaves all the .first pointers untouched.

# edge effect

                size_t newX = std::clamp(xIndex + dx, 0ul, BINS - 1);
size_t newY = std::clamp(yIndex + dy, 0ul, BINS - 1);


Hmmm, interesting. No "index out of bounds" errors. But bins at edge of screen shall be double counted, which seems undesirable. That is, if Percival the pigeon is flying near the center left edge, then result records a pair of pointers to Percival. And if he is flying near one of the corners, the effect is even worse.

# squared vs square root

We pass in delta_x and delta_y:

        if (std::hypot(bin.first[i]->position.x - position.x,
bin.first[i]->position.y - position.y) < radius) {


Finding length of hypoteneuse involves a slow sqrt(). But it would suffice to compare sum of squared deltas against radius². And then a good optimizing compiler would hoist it out of the loop, or you could help it out by assigning a radius_squared = radius * radius temp var.

# no velocity

struct Boid {
sf::Vector2f position;
};


What, no motion vector, no (dx, dy) speed? Ok, I guess it's early days in this prototype, we'll get there before too long.

Here's the query:

        for (auto& boid : boids) {
for (size_t i = 0; i < result.second; ++i) {
boidSeen.setPosition(result.first[i]->position);
window.draw(boidSeen);
}
}


So we're just identifying and highlighting static birds within radius, preparatory to velocity updates, OK. And at that point we can "head toward neighbor" when too far away, and "avoid collision" when too near.

When you are ready for velocity updates, you might consider having each boid maintain pointers to a few nearby boids at several distances, perhaps roughly doubling distances. For many update cycles the identity of a moving boid would be stable, and every so often you would have to pick a new neighbor at an appropriate distance. You could even manage a budget of only picking three such neighbor updates per displayed frame.

I don't understand the nesting there. For each boid, we visit all the spotlight boids? That seems unnecessary, and slow.

# aggregate weights

I was kind of expecting to see "the N-body problem" in the update loop, where we sum up gravitational attraction of our neighbors, or compute some similar potential function.

Distant boids make a negligible contribution, so you can approximate by ignoring them beyond some cutoff threshold. Aggregating their statistics to a bin, and mostly having those numbers influence velocity is an approach worth pursuing. (Maybe avoid aggregation for just the current bin a boid is in.) It might motivate using more than a 3 × 3 bin window.

# square bins

I note in passing that bins have the same aspect ratio as your display window. So you might want e.g. 20 × 15 bins for certain display sizes, if you're applying an isotropic potential function.

This code achieves a subset of its design goals.

It is written in a very clear style. I would be willing to delegate or accept maintenance tasks on it.

• I agree for std::pair this isn't the appropriate way, I could have used std::span instead, but I wanted a fast, simple but ugly thing. But I agree, it is ugly. Commented Nov 27, 2023 at 23:51
• @G.Sliepen, thank you, duly amended. There soo many languages on this stack exchange!
– J_H
Commented Nov 28, 2023 at 16:55

# The slowness comes from drawing

Your query() function is actually quite efficient. The problem comes from drawing the boids. If you remove the calls to boidSeen.setPosition() and window.draw() in the loop that calls query(), then you will have your high FPS back.

Note that you are drawing the same set of seen boids BOIDS times there.

# A 3x3 bin area might not be enough

If you run your program, you'll notice that as you move the mouse, not all boids inside the circle are yellow. That's because if RADIUS is 100, a 3x3 area around the center boid is not enough. You should either restrict the radius, or calculate the number of bins you need to check based on the radius.

# Using fixed-size bins might not work for a simulation

If you have a rather uniform distribution of boids, using a fixed grid of bins, which each can hold a fixed number of boids, is fine. However, if you would run a real simulation, you'll probably notice that there will be areas with very few boids and areas with much more boids, and this varies over time. Your approach might not work then, and instead you might need to use std::vector, std::list, or other dynamically sized data structures instead of std::array.

# Potential to reuse query results

Two boids that are very close together will mostly have the same set of neighbors. You can exploit this and avoid doing a lot of recalculations in query(). Consider for example the case that two boids are in the same bin, and the radius is at least $$\\sqrt 2\$$ times the spatial bin size.

Also consider that with a large enough radius or small enough bins, you can know which bins are fully covered by the circle of neighbors. In that case, you don't need to do the circle test. Furthermore, you can calculate the average position and average heading of all boids in a bin once, and then reuse that many times.

# Square both sides of the circle equation

std::hypot is a nice function, but it performs a square root operation which is quite slow. Instead of doing that, just calculate the difference in position squared, and compare that to the radius of the circle squared.