### How to Simulate and Render Blobs

Hi everyone,

It's blob time! This is a long one so buckle up.

If you haven't been following along on my tweeter for the last month, I've been working on simulating and rendering blobs. Here's a dump of my progress (easily accessible from the hashtag #vectorboy64). I'll add some captions

This first one is just me trying out position based dynamics and the constraint functions. PBD is a really cool way to simulate physics, and I'm really happy with how flexible it is and how easy it is to implement. I'm definitely going to use it for some other stuff in the future.

I spent a really long time trying to figure out the gradients of some tricky vector math formulas. I have them mostly correct, but there are still spurious rotations like this in the final version. The rotations come from the order of solving each constraint, and aren't really visible in a radially uniform object.

The collisions in this version are using Unity's PolygonCollider2D for collision detection. There are a lot of issues with this and I eventually abandoned this approach. One of the problems is that you have to solve all the blobs sequentially, another is that Unity doesn't give you a distance to edge function that works from within the object (why.. ??)

In this one, I've got it running on the job system, and I've switched to entirely particle-particle collisions. This has some advantages and some disadvantages. The advantage is that it's robust and simpler to implement. The disadvantages are that you need a lot of particles to make sure you don't have any gaps in the surface, otherwise the blobs can pass though each other. It's possible for blobs to pass though each other even if you have a continuous surface if the blobs are moving fast enough. I mostly solve this problem by increasing the radius of the each particle's collision and also just not having my blobs move too fast. Continuous collision detection would be the actual way to solve this problem but it's slow and a lot of work to implement probably.

Just running a lot of blobs at once now. Fun fact: 128 * 128 blobs is too many. It will crash Unity so hard that it will crash Windows. Ah, the joy of native code.

Here's rendering using an ad hoc model, just a power function based on the edge distance. I'll get to how we do the rendering in the main part of the post below.

I'm using the analytic formula for the profile of a liquid under surface tension to calculate the normals.
And this is the final(ish) result. I've tweaked the lighting and added blob nudging when a vectorboy lands on it.

## Position Based Dynamics

PBD is a particle based physics simulation method. Here's the paper if you want it. I found this video lecture to be much more approachable, however. This repo is also a good reference. Setting up the simulation loop is fairly trivial, just implement the pseudocode in the paper. The the core of the constraint update loop requires the gradient of the constraint function, and calculating these gradients is, in my option, the bulk of the work involved in implementing PBD.

The blobs are a set of particles that define the surface of the blob, and the interior is just rendered to make them look solid. With these particles, I implement two primary constraints: distance constraints between every neighboring particle, and an area preservation constraint among all particles. blobs with in-game settings particle radius reduced to show distance constraints

This has all sorts of fun problems like self-intersection and blob-blob penetration, especially if they're moving too fast. For the most part, these issues don't actually create super noticeable artifacts, so I'm not going to fix them. The exception is when a blob gets stuck inside another blob, or folds in on itself, but I am going to try hard to design around that by limiting speeds and forces.

In order to calculate the update for each particle, you have to find the gradient of the constraint function. For the distance constraint, this is given in the paper and video I linked above.

So uh, just use those. One thing that they mention is that this doesn't necessarily preserve angular momentum! Which is where the spurious rotations come from in my simulation.

The area preserving gradient took me a few days to figure out. Mostly because I suck at vector calculus.

If you can't tell, the important bit is that thing in the bottom right corner: (by, -bx). It's kind of hard to explain, but the direction that minimizes the change in area is just the direction pointing inward between each two adjacent particles.

There are two other constraints that are really simple, one constraint between each point and the center of their grid segment. This it to keep each blob in a relatively regular grid shape (which is important for something representing a pixel). The other constraint is collision constraints, which are generated during a collision detection step.

## Collision Detection

Collision detection is done by creating a few AABBs around a set of consecutive particles. In the final gif, each blob has 50 particles, and 5 AABBs. AABB count is a trade off between collision generation and solver time. More AABBs, more time spent generating collisions, but less time in each solver step because each particle has to check less other particles for overlaps. There's probably some spatial hashing garbage I can do on a grid to make this a lot faster, but for now it's good enough.

Now that we have simulated blobs, lets move on to rendering them.

## Rendering

For rendering each blob, I generate a triangle fan from the grid 'center' that spans each blob edge. I started out doing this on the CPU, but moved on to using Graphics.DrawProcedural() for this. Because the center position is fixed... it may actually lie outside the blob, and in that case there will be weird cone shaped artifacts, but that's why I have that center position constraint I mentioned earlier.

Getting SV_InstanceID into the pixel shader was quite an ordeal. When using surface shaders in Unity, they pass all variables from vertex to fragment shader as floats, which is uhh.. pretty not great when you want to pass an integer to the pixel shader. I spent a whole day debugging this thing:

The solution, if you can call it that, is to modify the generated surface shader to pass the int variables as ints like you would expect. It's not really sustainable for iteration purposes, but I don't know what else I can do.

(edit: actually you can pass it as float if you add some epsilon to it, and then round down to the right number in the pixel shader. Still kind of dumb, but it beats having to modify generated shader files)

The last part is calculating the normals on the blob.

## Surface Tension

According to wikipedia: "Surface tension is the tendency of fluid surfaces to shrink into the minimum surface area possible." The blob physics already does this in 2D from a top down perspective, but I also want the blobs to appear as if they do this in 3D. A simple way to do this is to use normal mapping to light the object as if it were 3D. Normals are unit vectors that point out perpendicular to a surface. The tricky bit is calculating what those normals should be.

The first step to calculating the normal is to find the distance from each pixel to the nearest edge. To do this, I use a technique called signed distance fields.

## Signed Distance

Signed distance fields are useful for computing a bunch of effects. In my case, I'm using the distance to the edge to choose the normal to use for lighting. But first, how do we calculate the distance to the closest edge? Ahahahah, well there's probably a fast way to do this, but I just iterate through all the edges. With 50 particles per blob there's only 50 edges, and GPUs are fast so whatever. If this becomes slow, I'll do something else. Distance to the edge visualization, with black being 0 distance from the edge.

## Inverting Surface Tension

So the equation for the profile of a curve under surface tension is this:

This equation is kind of crappy because it gives the profile in terms of h, which is the height. You give it a height and it gives you an x position back. This is not really what you would expect, and won't get us very far if we're trying to get the normal based on how far we are from the edge. To get the thing we want, we have to invert this formula. For the following examples we'll set H = 1 because that's a tuning constant.

Uh oh!

Looking at the graph, it should be obvious that there's no functional inverse to this because it overlaps itself! There's no one to one relationship between input points and output points. Now, we don't really care about the points below the curve, because we're always looking at blobs from directly above, so maybe there's a way to do some kind of piece-wise inversion but I am far too dumb to figure that out.

## Inverse of the Derivative Is Not the Derivative of the Inverse

The derivative of this thing is actually fairly simple, and you can invert it! Oops, that doesn't do what we want.

The inverse is this gnarly thing:

The inverse of the derivative tells you at what h position certain slope occurs... not the slope at a certain x position. This only took me 4 hours to figure out. Please learn from my mistakes.
There are a couple of theorems about the relationship between the derivative and the derivative of the inverse, but as far as I could tell, none of them really easily let you find the derivative of the inverse. The closest thing I found that might work is called Implicit Differentiation, which is what they use to find the derivatives of the inverse trig functions, but they all rely on cute trig identities, which I am not skilled enough to use.

## So, uhh, what now?

Ok, so I'm not good enough to invert this function, or at least get the derivative of the inverse so I can approximate it. Look up table~!

If we take enough samples of the input function, we can construct a piecewise linear model of it, and finding the inverse of that means searching through all the samples to find the X-position which corresponds to any Y position. Do this with uniformly spaced Xs and you've created a look up table for the inverse! This actually works and doesn't have too much error, but for areas of high curvature, it's off by a little bit, but I find this error acceptable Correct normals in red, approximated normals in green.

## Putting it together

Use the signed distance to look up the normal that we calculated from inverting that awful function! Actually, I don't store the normals,  just a lerp factor from pointing directly up to pointing directly out, because the outward direction depends on where on the blob we are and sending less data to the GPU is faster.

Overall, I'm pretty pleased with how it looks. Using Unity's standard shader, with a few directional lights makes it look really pleasing, like skittles or latex paint. If you crank up the metallic, they start to look like gems, or mylar balloons! Or turn down the smoothness and it becomes clay. Another advantage of using the lookup is that I can easily change the edge length - how quickly the curve is applied as the distance goes to zero.

That's about it, thanks for watching my video, remember to smash that like button. Lemme know if you have any questions on the tweeto @pyromuffin

XOXO,
Kelly

1. This looks amazing, gotta try it out.

2. 