Flocking simulation is erratic and choppy

Hi all,

I appreciate your time in reading this question. I’m fairly new to Processing and have finally made it into the world of vectors. It’s been funs so far, but I’m having trouble implementing the vectors in a 3D flocking simulation.

More specifically, I’m trying to make a flocking simulation on the surface of a sphere. Think almost like a flock of birds migrating across the globe. I’ve used conversions from spherical to cartesian coordinates, and I’ve applied the alignment, cohesion, and separation forces to the boids. I’ve made it very far and I think I’m super close, but for some reason I can’t get the flocks to behave naturally. Instead they wiggle erratically and move in a choppy format that looks unnatural. I am hoping instead that they would flow more with a natural motion.

I’ve been troubleshooting this for days and am really out of ideas. Hoping someone could take a look and point me in the right direction. Code pasted below. Apologies is it’s not formatted correctly to this forum (still getting the lay of land as I explore coding). Also, this wikipedia article about spherical coordinates is very helpful. I use the ISO convention. It also includes the conversion formulas between spherical and cartesian: Spherical coordinate system - Wikipedia

import peasy.*;
PeasyCam cam;

int boidCount=1000;
Boid[] boids = new Boid[boidCount];

int r=800;
void setup() {
  fullScreen(P3D, 2);
  cam = new PeasyCam(this, 1500);

  for (int i=0; i<boidCount; i++) {
    boids[i]= new Boid();
    boids[i].initialVelocity();
  }
}

void draw() {
  background(255);
  fill(255);
  stroke(0);
  sphere(r-10);
  for (int i=0; i<boids.length; i++) {
    //boids[i].normalVelocityUpdate();
    boids[i].flock();
    boids[i].gravity();
    boids[i].show();
  }
}

class Boid {

  float thetaOne=random(0, PI);
  float phiOne=random(0, TWO_PI);
  PVector position = new PVector(r*sin(thetaOne)*cos(phiOne), r*sin(thetaOne)*sin(phiOne), r*cos(thetaOne));
  float heading=random(0, TWO_PI);

  float thetaTwo;
  float phiTwo;
  float xTwo;
  float yTwo;
  float zTwo;
  PVector positionTwo;

  PVector velocity;
  PVector alignment;
  PVector cohesion;
  PVector separation;
  PVector flockingForces=new PVector(0, 0, 0);

  PVector avgAlignment = new PVector(0, 0, 0);
  PVector avgCohesion= new PVector(0, 0, 0);
  PVector avgSeparation= new PVector(0, 0, 0);

  float perceptionRadius=75;
  float alignmentFactor=0.8;
  float cohesionFactor=0.10;
  float separationFactor=0.3;

  PVector steerAlignment=new PVector(0, 0, 0);
  PVector steerCohesion=new PVector(0, 0, 0);
  PVector steerSeparation=new PVector(0, 0, 0);
  //PVector steeringTwo;
  float maxForce=4;
  float maxSpeed=3;

  PVector positionR;
  PVector positionRForce;


  void initialVelocity() {
    //Calculate initial velocity. To do that, find point tangent to circle and adjacent to starting position. Make velocity vector from current position to tangent position.

    //Calculate theta
    if (position.z>0) {
      thetaOne=atan((sqrt(pow(position.x, 2)+pow(position.y, 2))/position.z));
    } else if (position.z<0) {
      thetaOne=PI+atan((sqrt(pow(position.x, 2)+pow(position.y, 2))/position.z));
    } else if ((position.z==0)&&(sqrt(pow(position.x, 2)+pow(position.y, 2))!=0)) {
      thetaOne=HALF_PI;
    }

    //Calculate phi
    if (position.x>0) {
      phiOne=atan(position.y/position.x);
    } else if ((position.x<0)&&(position.y>=0)) {
      phiOne=atan(position.y/position.x)+PI;
    } else if ((position.x<0)&&(position.y<0)) {
      phiOne=atan(position.y/position.x)-PI;
    } else if ((position.x==0)&&(position.y>0)) {
      phiOne=HALF_PI;
    } else if ((position.x==0)&&(position.y<0)) {
      phiOne=-1*HALF_PI;
    }

    //Using the random heading random assigned, calculate an adjacent theta and phi
    thetaTwo=thetaOne+atan(cos(heading)/r);
    phiTwo=phiOne+atan(-sin(heading)/r);

    //Convert the new theta and phi back to cartesian
    xTwo=r*sin(thetaTwo)*cos(phiTwo);
    yTwo=r*sin(thetaTwo)*sin(phiTwo);
    zTwo=r*cos(thetaTwo);

    //Make velocity based on a vector from current position to adjacent position
    positionTwo=new PVector(xTwo, yTwo, zTwo);
    velocity=positionTwo.sub(position);
    velocity.setMag(maxSpeed);
    position.add(velocity);
  }

   void flock() {


    //Calculate alignment, cohesion, separation

    //ALIGNMENT
    avgAlignment=new PVector(0, 0, 0);
    steerAlignment=new PVector(0, 0, 0);
    int count=0;
    for (int i=0; i<boids.length; i++) {
      float d=PVector.dist(position, boids[i].position);
      if ((boids[i].position!=position) && (d<=perceptionRadius)) {
        avgAlignment.add(boids[i].velocity);
        count++;
      }
    }
    if (count>0) {
      avgAlignment.div(count);
      avgAlignment.limit(maxForce);
      steerAlignment=PVector.sub(avgAlignment, velocity);
      steerAlignment.limit(maxForce);
      steerAlignment.mult(alignmentFactor);
    }



    //COHESION
    count=0;
    avgCohesion=new PVector(0, 0, 0);
    steerCohesion=new PVector(0, 0, 0);
    for (int i=0; i<boids.length; i++) {
      float d=PVector.dist(position, boids[i].position);
      if ((boids[i].position!=position) && (d<=perceptionRadius)) {
        avgCohesion.add(boids[i].position);
        count++;
      }
    }
    if (count>0) {
      avgCohesion.div(count);
      avgCohesion.sub(position);
      avgCohesion.limit(maxForce);
      steerCohesion=PVector.sub(avgCohesion, velocity);
      steerCohesion.limit(maxForce);
      steerCohesion.mult(cohesionFactor);
    }

    //SEPARATION
    count=0;
    avgSeparation=new PVector(0, 0, 0);
    steerSeparation=new PVector(0, 0, 0);
    for (int i=0; i<boids.length; i++) {
      float d=PVector.dist(position, boids[i].position);
      if ((boids[i].position!=position) && (d<=perceptionRadius)) {
        PVector diff=PVector.sub(position, boids[i].position);
        diff.div(pow(d, 0.75));
        avgSeparation.add(diff);
        count++;
      }
    }
    if (count>0) {
      avgSeparation.div(count);
      avgSeparation.sub(velocity);
      avgSeparation.limit(maxForce);
      steerSeparation=PVector.sub(avgSeparation, velocity);
      steerSeparation.limit(maxForce);
      steerSeparation.mult(separationFactor);
    }

    flockingForces.mult(0);                 //flocking forces is the same as "acceleration"
    flockingForces.add(steerAlignment);
    flockingForces.add(steerCohesion);
    flockingForces.add(steerSeparation);
    velocity.add(flockingForces);
    velocity.setMag(maxSpeed);
    position.add(velocity);
  }
  
  
   void gravity() {
    //To keep point on surface and avoid it from flying off or into the sphere, find theta and phi of current point, and then make a new point with the same theta and phi but at only a distance r on the sphere. Make a new vector from current position to position that is on the sphere's surface.


    //Calculate phi
    if (position.z>0) {
      thetaOne=atan((sqrt(pow(position.x, 2)+pow(position.y, 2))/position.z));
    } else if (position.z<0) {
      thetaOne=PI+atan((sqrt(pow(position.x, 2)+pow(position.y, 2))/position.z));
    } else if ((position.z==0)&&(sqrt(pow(position.x, 2)+pow(position.y, 2))!=0)) {
      thetaOne=HALF_PI;
    }

    //    //CALCULATE PHI
    if (position.x>0) {
      phiOne=atan(position.y/position.x);
    } else if ((position.x<0)&&(position.y>=0)) {
      phiOne=atan(position.y/position.x)+PI;
    } else if ((position.x<0)&&(position.y<0)) {
      phiOne=atan(position.y/position.x)-PI;
    } else if ((position.x==0)&&(position.y>0)) {
      phiOne=HALF_PI;
    } else if ((position.x==0)&&(position.y<0)) {
      phiOne=-1*HALF_PI;
    }

    //Cartesian coordinates of new point on the surface
    xTwo=r*sin(thetaOne)*cos(phiOne);
    yTwo=r*sin(thetaOne)*sin(phiOne);
    zTwo=r*cos(thetaOne);
    positionR=new PVector(xTwo, yTwo, zTwo);
    positionRForce=positionR.sub(position);
    velocity.add(positionRForce);
    velocity.setMag(maxSpeed);
    position.add(velocity);
  }

  void show() {
    float thetaDraw = velocity.heading2D() + radians(90);
    float rDraw=2;

    fill(0);
    stroke(255);
    pushMatrix();
    translate(position.x, position.y, position.z);
    rotate(thetaDraw);
    beginShape(TRIANGLES);
    vertex(0, -rDraw*2);
    vertex(-rDraw, rDraw*2);
    vertex(rDraw, rDraw*2);
    endShape();
    popMatrix();
  }
  }

I haven’t tried your code yet, but two things jump out at me. First, never use atan(), instead use atan2() which keeps x and y separate without a potential divide by zero. Second, to keep the boids on the surface of the sphere, instead of converting to spherical coordinates and back, just normalize their position vector, i.e., position.normalize(). Drop the spherical coordinates entirely and just keep everything as 3-D vectors with a final projection onto the sphere.

You could even project onto other surfaces too: https://genart.social/@scdollins/109779417062229176

As to the choppiness, consider what happens when two boids go from greater than one of your cutoff distances to less than. They move along smoothly until they suddenly experience a new force. To keep things smooth, that new force should start at a strength of zero and then only ramp up smoothly as the boids get closer to each other. If the forces don’t go to zero at the threshold, then the boids experience a sudden impulse like a collision and in a flock, they will “bounce” around off of their neighbors.

https://www.youtube.com/watch?v=xiUpAeos168 does a decent job of describing the sort of force functions you might like to use.

Ideally the motion would end up smooth as in https://genart.social/@scdollins/111786009282547893

1 Like

Thank you for the pointers! Doing one final projection would clean up the math and code so much. Thank you!

Good call, never thought about mapping my forces based on distance. Will implement. And thank you for the links. The spherical flow field/flocking gif is awesome!