Optimizing the Particle Life: Spatial Partitioning with QuadTree (Part 2)
TypeScript Particle Life
Now that the decision is made, let’s get to work. I’ll start with writing a basic JavaScript (or, to be precise, TypeScript that will later compile to JS) implementation based on the Python code I created in the previous part. To improve the code, I’ll need to add the following:
- Space wrapping around the edges.
- Space partitioning for optimization purposes.
- Allowing the user to adjust the parameters.
The code below is part of a Vue component. VueJS is great when you need a lightweight, simple, and robust tool for rapid development. While the whole application structure is out of the scope of this article, please have a look at this GitHub repo to see how I organized my project.
I will only explain in detail the parts not covered in the previous article, so make sure you read it first. A very important disclaimer: comments regarding optimization (e.g. typed vs regular arrays) only refer to code running in Chrome browser. Performance will vary from environment to environment, and it would have to be taken into consideration in real-life products.
const canvas = document.getElementById("pl") as HTMLCanvasElement; //Grab a reference to Canvas element and cast its type as HTMLCanvasElement to have access to all its properties and methods
const context = (canvas).getContext("2d"); //Create 2D rendering context for our canvas
//Make our canvas full-screen
canvas!.width = document.documentElement.clientWidth
canvas!.height = document.documentElement.clientHeight
//Since our x,y coordinates have values between (0,1), we need screenRatio when drawing particles, so that we don't render stretched/distorted animation.
const screenRatio = canvas!.width / canvas!.height
//Let's store our systems parameters in localStorage, so that they persist between user's session. If no 'parameters' object is in localStorage, fallback to default config.
const parameters = window.localStorage.getItem('parameters') ? JSON.parse(window.localStorage.getItem('parameters')!) : defaultConfig
const m = parameters.colorsNumber;
const n = parameters.particlesCountJS;
const radius = parameters.interactionRadiusJS;
const fr = Math.pow(0.5, 20);
const matrix = createRandomMatrix(m);
//using typed arrays would give us no performance advantage
let colors: number[] = [];
let positionsX: number[] = [];
let positionsY: number[] = [];
let velocitiesX: any[] = [];
let velocitiesY: any[] = [];
for (let i = 0; i < n; i++) {
colors[i] = Math.floor(Math.random() * m);
positionsX[i] = Math.random();
positionsY[i] = Math.random();
velocitiesX[i] = 0;
velocitiesY[i] = 0;
}
renderFrame();
function renderFrame() {
//Clear previous positions on every frame
context!.fillStyle = "black";
context!.fillRect(0, 0, canvas!.clientWidth, canvas!.clientHeight);
for (let i = 0; i < n; i++) {
//beginPath() method is called before drawing every particle, so that they may be filled with different colors.
context!.beginPath();
let x = positionsX[i] * canvas!.clientWidth;
let y = positionsY[i] * canvas!.clientHeight
//Adjusting position for two cases:
//a) mobile screens -> x<y
//b) desktops -> x>y
x = screenRatio < 1 ? x * (1 / screenRatio) : x;
y = screenRatio > 1 ? y * screenRatio : y;
if (!(x < 1) && !(y < 1) && !(x > canvas!.clientWidth - 1) && !(y > canvas!.clientHeight - 1)) { //avoid blinking at edges due to how space wrapping is done
//Drawing the actual particle
context!.arc(x, y, 1, 0, 2 * Math.PI); //js compilers are smart enough for me not to have to move math.pi to var
context!.fillStyle = `hsl(${360 * (colors[i] / m)},100%,50%)`;
context!.fill();
}
}
updateParticles();
// Native window method to request animation frame. Note: The frequency of calls aims to match the display refresh rate.
// This method passes DOMHighResTimeStamp as an argument to callback, so that we can calculate how much our animation should progress in one frame
// Otherwise there will be different speed for animation on screens with different refresh rates
requestAnimationFrame(renderFrame);
}
function updateParticles() {
for (let i = 0; i < n; i++) {
let totalForceX = 0;
let totalForceY = 0;
for (let j = 0; j < n; j++) {
if (i === j) continue;
let x = positionsX[j] - positionsX[i];
let y = positionsY[j] - positionsY[i];
let d = Math.sqrt(x * x + y * y); // No performance advantage in removing math.sqrt. 10x faster than math.hypot.
if (d < radius) {
let f = force(
d / radius,
matrix[colors[i]][colors[j]])
totalForceX += (x / d) * f;
totalForceY += (y / d) * f;
}
velocitiesX[i] += totalForceX;
velocitiesY[i] += totalForceY;
velocitiesX[i] *= fr;
velocitiesY[i] *= fr;
positionsX[i] += velocitiesX[i];
positionsY[i] += velocitiesY[i];
x = positionsX[i];
y = positionsY[i];
// Prevent particles from running away from the viewport.
// If x or y is > 1, substract 1 from it, so that it moves to other part of the screen
// This only handles particle position - it's not interacting with
// particles from the other side of the screen. In other words, particle at position (0,0) with radius of 0.1
// will not interact with particle at position (1.0,1.0). This gives us interesting visual effect (see for GPU.js examples screenshots)
// but also causes some issues, like particles aggregating at the edges
x = x > 1 ? x - 1 : x;
y = y > 1 ? y - 1 : y;
x = x < 0 ? x + 1 : x;
y = y < 0 ? y + 1 : y;
positionsX[i] = x;
positionsY[i] = y;
}
}
}
//Instead of using one force function allow user to choose from several options.
//Graphs of every function (assuming the force value between two particles to be 1.0 and a repulsive radius of 1/4 of the max radius)
//available on live site of the project
function force(d: number, f: number): number {
let repRadius = parameters.repulsiveRadius;
switch (parameters.forceFunction) {
case '1':
return forceVariantI(d, f, repRadius);
case '2':
return forceVariantII(d, f, repRadius);
case '3':
return forceVariantIII(d, f, repRadius);
default:
return defaultForce(d, f, repRadius)
}
}
function forceVariantI(d: number, f: number, repRadius: number): number {
if (d <= repRadius) {
return d / repRadius - 1
} else {
return f * (1 - Math.abs(2 * d - 1 - repRadius) / 1 - repRadius);
}
}
function forceVariantII(d: number, f: number, repRadius: number): number {
if (d <= repRadius) {
return (d / repRadius - 1) + (f * (1 - d))
} else {
return f * (1 - d)
}
}
function forceVariantIII(d: number, f: number, repRadius: number): number {
return defaultForce(d,f,repRadius)
}
function defaultForce(d: number, f: number, repRadius: number): number {
if (d <= repRadius) {
return d / repRadius - 1
} else {
return f * (1 - d)
}
}
QuadTree
Now that my code is running in the browser, it's time for optimization. We have a set of points in 2D space that can only interact with one another within a given radius. To see which particles fall within this range, we must examine every particle in the system. This is the root cause of the inefficiency of our program.
So, how do we retrieve particles that only reside in the surroundings of the particle we’re calculating the velocity of? We would need a data structure that stores regions of space and their corresponding particles. This way, we could just query regions within a given area and thus get all particles living there. Here comes the QuadTree!
Consider the root node, which will translate to our entire 2D plane:
First, let’s examine if our node has more particles than its carrying capacity. Let’s assume the carrying capacity for every sector is 4. Below, we have 5 particles. The number of particles exceeds the carrying capacity of our sector, so we need to subdivide it. We’ll divide one sector into 4 smaller ones:
.png)
Now, let’s repeat the first step for every one of our new nodes. In this case, none of the child nodes exceed the carrying capacity, so we stop here, and this is our final quadtree. Had it not been the case, we’d repeat the division recursively, creating smaller and smaller sections for every node that exceeds its capacity. This makes it a very good data structure for our simulation, where attraction between particles will usually lead to an uneven distribution of particles. Only crowded sectors get subdivided.
Now, every sector needs the following properties:
- boundary - a region of space it covers. We can represent it as a rectangle with a center point (x,y) and its height and width;
- capacity - carrying capacity;
- points - an array of points, given the sector is not subdivided;
- divided - indicates if the sector is subdivided;
Of course, we also need at least methods to insert, delete, and query particles. These methods will work recursively. For example, to implement query() method, we need to:
- Check if the queried area intersects with the root node
- If it doesn’t return an empty array
- If it does and the root node is not subdivided, return all points belonging to this node
- If it does and the node is subdivided, repeat the same for every child node
This way, we can efficiently traverse our tree and return only relevant points. I decided to use ready implementation, as I don’t like reinventing the tree - thanks to The coding train for sharing his implementation
Drawbacks
Spatial partitioning can already take our Particle Life to another level, but it is far from an ideal solution. In the end, what limits us is the number of particles in the range that we need to test. If we want to keep things efficient, the radius needs to be small enough to contain a limited number of particles. This way, we create a zoomed-out, uniform version of the particle life. It’s an interesting optimization exercise, but nothing more. The more particles we want, the smaller the radius we need to set, and the less we see.
What if we wanted to see how, let’s say, 100000 particles in range interact with one another? To achieve this, we need to think in a different way. Instead of sequentially calculating forces for every particle (which would require 10 billion checks for every frame), let’s use the GPU to calculate force for every particle in parallel. So it happens that I was reading about the new webGPU API a few days ago, and I think it would be a nice programming experiment to implement our particle life this way.
Further considerations
If we were to continue developing our javascript implementation, there are more things we could do to optimize the code. Our system is very dynamic, and everything is moving with every frame, forcing us to reconstruct the tree. It happens to be quite an expensive operation.
To further improve this, we could draw from data structure used in computer simulations, verlet list. Instead of querying particles within a maximum interaction radius, we add a padding of X%. Then, we consider how many repeated position updates we can calculate before particles move outside of our adjusted radius. This way, we can rebuild our quadtree not on every frame, but once in every few frames. The tricky part is finding a balance between an increased radius (and hence more particles to loop over) and the frequency of rebuilding the quadtree. This will vary depending on our system parameters.
Particle Life online
I deployed my code to Netlify, where you can play with an unoptimized JS version, quadtree enhanced code, and webGPU 3D particles. Remember to increase the parameters carefully if you’re not on a good machine, otherwise your browser will freeze. The code is available here.
See you in the next part, where I’ll take things to the next level by implementing and optimizing the webGPU version!