Boids¶

Shared by Antoine Rideau

On this page you will find how to simulate using Castor the flocking behaviour of birds.

Boids - contraction of “bird-oid object” - refers to an artificial life program, developed by Craig Reynolds in 1986, simulating the flocking behaviour of birds.

The complex flocking behaviours arise from three simple rules describing the interaction between each boid:

Separation: steer to avoid crowding local flockmates.
Alignement: steer towards the average heading of local flockmates.
Cohesion: steer to move toward the average position of local flockmates.

The flock is composed of N boids and will fly during nt steps within a perimeter described by dimension

// Parameters
int N = 100;  // Number of boids
int nt = 300; // Number of time steps
matrix<> dimension = {50, 50}; // Permieter lengths

Each boids in the flock is characterized by its position \((x,y)\) and its velocity vector \(\overrightarrow{v}=(dx,dy)\) . Each parameter is stored in a matrix<> and then gathered inside the std::vector Boids.

Furthermore the matrix Out contains these parameters at each steps in order to output them later on for a satisfaying visualization.

// Initialization
std::vector<matrix<>> Boids = { rand(1, N,true) * dimension(0),             // Boids x position
                                rand(1, N) * dimension(1),                  // Boids y position
                                rand(1, N,true) * maxSpeed - maxSpeed / 2,  // Boids x velocity : dx
                                rand(1, N) * maxSpeed - maxSpeed / 2};      // Boids y velocity : dy
matrix<> Out = cat(1, Boids[0], cat(1, Boids[1], cat(1, Boids[2], Boids[3])));

See rand , cat .

_images/boids.png — Neighbourhood of a boid (blue) : zone of separation (red), zone of alignement (yellow) and zone of cohesion (green).¶

So as to apply the set of rules to each boid, it is needed to determine this boid’s flockmates depending on the distance between this boid and the others boids. This is done by the function boidsWithinDistance which returns the index of the boids who are inside the circle of distance radius around the boid of index boidIndex.

matrix<size_t> boidsWithinDistance(std::vector<matrix<>> &Boids, int boidIndex, double distance)
{
    matrix<size_t> I(boidIndex);
    return setdiff(find(sqrt((Boids[0](boidIndex) - Boids[0]) * (Boids[0](boidIndex) - Boids[0]) + (Boids[1](boidIndex) - Boids[1]) * (Boids[1](boidIndex) - Boids[1])) < distance), I);
}

See setdiff , find .

Behaviour¶

Borders¶

There are two boundaries conditions that can be applied :

Boids must stay within the borders.

// Constrain boids to within the borders. If it gets too close to an edge,nudge it back in and reverse its direction.
void stayWithinBorders(std::vector<matrix<>> &Boids, matrix<> border, double margin, double turnFactor)
{
    auto xTooLow = find(Boids[0] < margin);
    auto xTooHigh = find(Boids[0] > border(0) - margin);
    auto yTooLow = find(Boids[1] < margin);
    auto yTooHigh = find(Boids[1] > border(1) - margin);

    if (numel(xTooLow) > 0){Boids[2](xTooLow) = eval(Boids[2](xTooLow)) + turnFactor;}
    if (numel(xTooHigh) > 0){Boids[2](xTooHigh) = eval(Boids[2](xTooHigh)) - turnFactor;}
    if (numel(yTooLow) > 0){Boids[3](yTooLow) = eval(Boids[3](yTooLow)) + turnFactor;}
    if (numel(yTooHigh > 0)){Boids[3](yTooHigh) = eval(Boids[3](yTooHigh)) - turnFactor;}
}

See find , numel , view .

Like Pacman, boids who go out by one side are wrap around to the other side.

//  Checks if boids go out of the window and if so, wraps them around to the other side.
void wrapBorders(std::vector<matrix<>> &Boids, matrix<> border)
{
    auto xTooLow = find(Boids[0] < 0);
    auto xTooHigh = find(Boids[0] > border(0));
    auto yTooLow = find(Boids[1] < 0);
    auto yTooHigh = find(Boids[1] > border(1));

    if (numel(xTooLow) > 0){Boids[0](xTooLow) = eval(Boids[0](xTooLow)) + border(0);}
    if (numel(xTooHigh) > 0){Boids[0](xTooHigh) = eval(Boids[0](xTooHigh)) - border(0);}
    if (numel(yTooLow) > 0){Boids[1](yTooLow) = eval(Boids[1](yTooLow)) + border(1);}
    if (numel(yTooHigh > 0)){Boids[1](yTooHigh) = eval(Boids[1](yTooHigh)) - border(1);}
}

See find , numel , view .

Separation¶

For each boid, once flockmates within separationDistance are known, velocity vector is ajusted toward the opposite directions of flockmates’ positions.

// Boids try to keep a small distance away from other boids
void separation(std::vector<matrix<>> &Boids, double separationDistance, double separationFactor)
{
    for (int i = 0; i < numel(Boids[0]); i++)
    {
        auto FlockMates = boidsWithinDistance(Boids, i, separationDistance);
        if (numel(FlockMates) > 0)
        {
            Boids[2](i) += sum(Boids[0](i) - eval(Boids[0](FlockMates))) * separationFactor;
            Boids[3](i) += sum(Boids[1](i) - eval(Boids[1](FlockMates))) * separationFactor;
        }
    }
}

See sum , view .

Alignement¶

For each boid, once flockmates within alignementDistance are known, velocity vector is ajusted toward the average direction of flockmates’ velocity vector.

// Boids try to match velocity with near boids
void alignement(std::vector<matrix<>> &Boids, double alignementDistance, double alignementFactor)
{
    for (int i = 0; i < numel(Boids[0]); i++)
    {
        auto FlockMates = boidsWithinDistance(Boids, i, alignementDistance);
        if (numel(FlockMates) > 0)
        {
            Boids[2](i) += (sum(eval(Boids[2](FlockMates))) / numel(FlockMates)) * alignementFactor;
            Boids[3](i) += (sum(eval(Boids[3](FlockMates))) / numel(FlockMates)) * alignementFactor;
        }
    }
}

See sum , numel , view .

Cohesion¶

For each boid, once flockmates within cohesionDistance are known, velocity vector is ajusted toward the average position of flockmates.

// Boids try to fly towards the centre of mass of neighbouring boids
void cohesion(std::vector<matrix<>> &Boids, double cohesionDistance, double cohesionFactor)
{
    for (int i = 0; i < numel(Boids[0]); i++)
    {
        auto FlockMates = boidsWithinDistance(Boids, i, cohesionDistance);
        if (numel(FlockMates) > 0)
        {
            Boids[2](i) += ((sum(eval(Boids[0](FlockMates))) / numel(FlockMates)) - Boids[0](i)) * cohesionFactor;
            Boids[3](i) += ((sum(eval(Boids[1](FlockMates))) / numel(FlockMates)) - Boids[1](i)) * cohesionFactor;
        }
    }
}

See sum , numel , view .

Speed limitation¶

Once the rules are applied on the boids, their speed is limited to maxSpeed .

void limitSpeed(std::vector<matrix<>> &Boids, double maxSpeed)
{
    auto Speed = sqrt(Boids[2] * Boids[2] + Boids[3] * Boids[3]);
    matrix<> I = find(Speed > maxSpeed); //Index where velocity > maxSpeed

    if (numel(I) > 0)
    {
        Boids[2](I) = (eval(Boids[2](I)) / eval(Speed(I))) * maxSpeed;
        Boids[3](I) = (eval(Boids[3](I)) / eval(Speed(I))) * maxSpeed;
    }
}

See sqrt , find , view .

Visualisation¶

In order to have a more satisfaying animation of the flying boids Python will be used.

Beforehand, the data stored in Out are written in text file.

writetxt("../", "data.txt", Out);

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.patches as patches
import numpy as np

data = np.loadtxt("./data.txt", skiprows=1)
# print(data)

nt = int(data.shape[0]/4)
N = int(data.shape[1])

xs = [data[i*4, :] for i in range(nt)]
ys = [data[i*4+1, :] for i in range(nt)]
dxs = [data[i*4+2, :] for i in range(nt)]
dys = [data[i*4+3, :] for i in range(nt)]

# Visu initialization
fig, ax = plt.subplots()
p = patches.Rectangle((0, 0), 50, 50, fill=True)
ax.add_artist(p)
ax.set_xlim(0, 50)
ax.set_ylim(0, 50)
ax.set_aspect('equal')

# Length of velocity vector to then normalize its representation
D = [np.sqrt((dxs[0][i])**2+(dys[0][i])**2) for i in range(N)]
boids = [ax.annotate("", xy=(xs[0][i]+dxs[0][i]/D[i], ys[0][i] + dys[0][i]/D[i]), xytext=(xs[0][i], ys[0][i]), arrowprops={
                     "facecolor": "red", 'arrowstyle': 'wedge'})for i in range(N)]


def animate(frame):
    for i in range(N):
        d = np.sqrt((dxs[frame][i])**2+(dys[frame][i])**2)
        pos = np.array([xs[frame][i], ys[frame][i]])

        boids[i].set_position(pos)
        boids[i].xy = pos + (dxs[frame][i]/d, dys[frame][i]/d)

    return boids


# Creating the Animation object
ani = animation.FuncAnimation(
    fig, animate, nt, interval=40, blit=True, repeat_delay=1000)
plt.show()

Flight of 100 boids in a 50-by-50 square.

Code¶

Here is all the code at once, without the functions written above :

int main(int argc, char const *argv[])
{
    // Parameters
    int N = 100;  // Number of boids
    int nt = 300; // Number of time steps
    matrix<> dimension = {50, 50};

    double margin = 2.;
    double turnFactor = 1.;
    double maxSpeed = 1.5;
    double separationDistance = 2.0;
    double separationFactor = 0.05;
    double alignementDistance = 7.5;
    double alignementFactor = 0.10;
    double cohesionDistance = 8.5;
    double cohesionFactor = 0.03;

    // Initialization
    std::vector<matrix<>> Boids = { rand(1, N,true) * dimension(0),              // Boids x position
                                    rand(1, N) * dimension(1),              // Boids y position
                                    rand(1, N,true) * maxSpeed - maxSpeed / 2,   // Boids x velocity : dx
                                    rand(1, N) * maxSpeed - maxSpeed / 2};  // Boids y velocity : dy
    matrix<> Out = cat(1, Boids[0], cat(1, Boids[1], cat(1, Boids[2], Boids[3])));

    disp(Out);

    tic();
    for (int t = 0; t < nt; t++)
    {
        separation(Boids, separationDistance, separationFactor);
        alignement(Boids, alignementDistance, alignementFactor);
        cohesion(Boids, cohesionDistance, cohesionFactor);
        limitSpeed(Boids, maxSpeed);
        stayWithinBorders(Boids, dimension, margin, turnFactor);

        // Update positions
        Boids[0] += Boids[2];
        Boids[1] += Boids[3];

        // wrapBorders(Boids, dimension);

        for (int i = 0; i < 4; i++)
        {
            Out = cat(1, Out, Boids[i]);
        }
    }
    toc();

    writetxt("../", "data.txt", Out);

    return 0;
}

References¶

https://eater.net/boids

https://en.wikipedia.org/wiki/Boids

http://www.red3d.com/cwr/boids/

http://www.vergenet.net/~conrad/boids/pseudocode.html