URCHN Arkipelago Roach Crowd Sim

Roach Crowd Sim

From URCHN Arkipelago

Rough concept

Similar to boids except that the cockroachs' bodies are parented to an empty/root bone which is shrinkwrapped to a low poly mesh (offset some distance from the floor) which marks the area the roach's are constrained to.

The system uses empties to input data. Custom properties define whether an empty represents an aligner, a control or a boid (roach). Controls (which can be of type 'follow' 'attract' 'repel' 'herd' 'panic' 'sleep' etc) are empties which affect boids within a certain radius during the simulation time. Proximity data is refreshed on a frame by frame basis. Aligners are used during setup to quickly orientate and configure a large number of boids using a weighted linear fall off to interpolate starting direction and speed between multiple aligners.

The roachs have an emotional state which is adjusted each frame by various influences (including how crowded the flock is, whether they have been involved in a collision, if they have strayed from the flock, or if they are under the influence of a controller). Emotional states decay/normalize over time.

The roachs have a friend/foe database. Roachs they have collided with slide to the foe end of the scale, wheras roachs they have herded with (moved parallel to) become their friends.


Python is incredibly slow. With 500 boids there are 500x499 inter-boid distances (or separations) to measure. These are stored in a half-matrix (1/2 n (n + 1)) indexed by [boid with lower index of pair][boid with higher index of pair]. The separations are also indexed in another dict by their value ie dict[separation distance]=[index a][index b]. Both sortings of the separations are required so although the script is basically holding duplicate data, this will be faster than reindexing the half-matrix by separation size each time we need to find the pair which is closest together, or a list of all pairs which have collided.

Having to regenerate all the separations every frame would make the script run extremely slowly. Instead, we can only regenerate the critical separations (ie those inter-roach distances which are small enough to mean that either a collision has taken place, or that a roach can 'see' another roach - the user defines a 'sight distance' for the roaches). The sight distance will obviously be larger than the 2 * roach radius (the collision distance). As we know the maximum distance a roach can travel in a single frame (defined by the user), we can calculate how often separation data will need to be refreshed.

Ignoring the roach's radius for simplicity's sake, if a roach can see 1 unit away, and the maximum speed of a roach is 0.1 units, then any pairs of roaches whose separation is less than / equal to 1.1 units will need to be refreshed after 1 frame. After 2 frames a roach which was 1.2 units away could now be within the sight distance (as its max speed 0.1*2 frames = 0.2... + sight distance of 1 = 1.2). After 3 frames a roach which was 1.3 units away... and so on.

Once data has been refreshed, it lasts for a set number of frames before it needs to be refreshed again:

  • roachs within distance 0-1.1 need regenerating every frame
  • roachs within distance 1-1.1.2 need regenerating every other frame
  • roachs within distance 1-2.1.3 need regenerating every third frame
  • roachs within distance 1-3.1.4 need regenerating every fourth frame
  • roachs within distance 1-4.1.5 need regenerating every fifth frame

So for example on the 6th frame, the range from 0-1.3 and 1.5-1.6 needs refreshing, but the data from 1.3-1.5 is still valid as no roaches within this range or radii could possibly be within the sight distance even if they have moved at maximum speed.

My python code currently keeps track of the age of the separation data so only the bare minimum of data gets refreshed each frame saving a vast amount of computation time.


Overlap Collisions: (where the bounding sphere for one roach overlaps with the bounding sphere of another) These are resolved by freezing one roach and moving the other along the vector from the centre of the frozen roach to the centre of the unfrozen roach until the unfrozen roach no longer overlaps. The reason for freezing rather than standard Newtonian physics calculations is that when multiple roachs collide in a single frame, calculations can simply begin at the most locked in roach and work outwards from there, freezing each roach as it is processed. The unfrozen roach's rotation is pointed away from the frozen roach instantaneously. The rotation of the roachs represents the direction of their movement NOT the orientation of their bodies - ie a roach can face forwards but walk backwards out of a collision. The orientation of their bodies can be determined after the script has run by smoothing over the rotation curves.

Early Warning System: Roachs have a 'last minute repel' whereby if they are getting very close to another roach they begin to gradually change direction and rotation away from it. The rate of change is determined by how close they are to an overlap collision. This is intended to prevent the overlap collisions from ever executing, however it will not force a roach out of a collision in the same way the overlap collisions algorithm will.

Escape Route Searching: Roachs look around them and determine the bearing of all their neighbours (projected into the roach's local XY plane to take account for when roachs are not on a perfectly flat surface). They look for the largest gap or escape route possible by finding the largest of the difference between subsequent sorted bearings of neighbours. For example if there is a roach due east of them and a roach due south then the best escape route is to the north west. The roach will gradually change rotation (the rate of change determined by the global simulation speed) towards this bearing.

Surface Constraints

Failed Attempt: Adding a shrinkwrap constraint for every object worked, but excessively slowly. It seems that each constraint builds an acceleration tree but perhaps that these trees do not talk to each other (so trees are generated per constrained object rather than per mesh). Another problem is that animation keys affect whether a constrained object's matrix (which contains the visual location, rotation, scale etc in the matrix_world.translation_part() .rotation_part() etc...). This may be a bug in blender. The only stable way to convert a desired loc/rot/sc into a visual loc/rot/sc is:

  • 1. set the desired loc/rot/sc on boid_X
  • 2. key the visual loc/rot/sc on boid_X
  • 3. read back the visual loc/rot/sc from boid_X
  • 4. check the returned loc/rot/sc from boid_X does not exceed any limits imposed by the sim (max, min speed, rotation limits etc)
  • 5. set the adjusted loc/rot/sc on boid_X
  • 6. key the loc/rot/sc on boid_X

It seems that in keying using bpy.ops.anim.keyframe_insert() the fcurves are regenerated (or some other slow operation takes place) wheras using object.keyframe_insert() is much faster. However object.keyframe_insert() does not consistently update the matrix (and visual loc/rot/sc) and so is useless for this purpose.

Proposed Method: At script start create a single empty constrained to the surface (constrained_empty). This empty must not have any keys on it. The boid empties must be free of constraints. The method for keying will be as follows:

  • 1. set the desired loc/rot/sc for boid_X on constrained_empty without keying anything
  • 2. may need to update the matrix using scene.update() and marking constrained_empty for update
  • 3. read back the visual loc/rot/sc from constrained_empty.matrix_world
  • 4. check the returned loc/rot/sc does not exceed limits (as before)
  • 5. set the adjusted loc/rot/sc on boid_X
  • 6. key boid_X using boid_X.keyframe_insert() and specifying relevant data paths.

Fingers crossed this method will key much faster by using boid_X.keyframe_insert() specifying to not regenerate fcurves, and by having only one constrained empty (resulting in memory savings etc).

Surface Orientation

Failed Attempt: Initially I proposed having two empties parented to each boid (one at local position (1,0,0) and the other at (0,1,0)) and constrained to the surface to gather enough information to set the local x and y axes orientations of each boid. This suffered from the same speed reduction problems as the above method for surface constraints.

Proposed Method 1: Without using any constraints gather the last two direction vectors for each boid, crossing these gives a z-up so long as:

  • 1. they are crossed in a clockwise direction (use global z-up or last z-up to determine which way is clockwise).
  • 2. the angle between them is greater than 0 and as close to 90 as possible.
  • 3. the variation in surface tangent is not too severe.

This method should be significantly faster than the failed attempt and proposed method 2.

Proposed Method 2: As the failed attempt, but only have local axes empties on the one constrained_empty rather than on every boid.