Rendering a particle system as a smooth blob

What is a good algorithm for rendering a particle system as a coherent mesh, in realtime? I am running a smoothed-particle hydrodynamics fluid simulation, and I would like to render the particles not as balls but as a body of water.

I understand that the basis is to find all places where the combined density of the particle system is equal to a certain boundary value, which will yield a surface which is smooth. However, how do you actually do so quickly and translate it into vertices and triangles to render?

Another option is raytracing, raymarching specifically. Just march until the density reaches the boundary value. However, I'm still stuck on how to get the surface normal in this case, as well as how to make it performant. Surely you don't need to iterate through all particles, the far away particles won't have any effect, but I don't know how to spatially cull particles.

1. What are some algorithms for spatially culling far away particles with no local contribution?
2. What are some algorithms for getting the surface normal of a body of water made of particles?
3. What are some algorithms for getting the surface of a body of water made of particles, in the form of a mesh?
• You can use an acceleration structure to avoid raymarching like crazy. As far as rasterization goes, you could use marching cubes. For the normals you can use the numerical derivative at the point of the intersection. Nov 7 '20 at 0:29
• @lightxbulb numerical derivative of what with respect to what? I don't know how I would get the gradient of the surface if I only get an intersection coordinate (using raymarching). Also, isn't marching cubes only applicable to grid structures? Nov 7 '20 at 0:36
• The derivative of the density with respect to the global coordinates. Check out Inigo Quilez's site, or some shadertoy with numerical derivatives. For marching cubes you can make a grid that contains your particles, and estimate the density at each grid point. Nov 7 '20 at 10:29
• You might try searching google for "Particle system isosurface extraction" this brings up several good papers on the subject, the marching cubes algorithm may also be a practical solution..I've seen particle systems composed of metaballs efficiently animated using marching cubes to extract the mesh but haven't (yet) written such a system myself. Nov 9 '20 at 12:04

For culling faraway particles: as lightxbulb said, the keyword to search for is "acceleration structure". There are a variety of options: BVH, octree, kd-tree, uniform grid, or some combination of these. Which one will be most suitable is going to depend on the details of your situation—how large/dense is the simulation, how much you care about performance vs memory vs implementation complexity, and so on.

They all work by subdividing space into volumes in some fashion and keeping track of which particles are in each volume. Then, when you're marching you only need to keep track of which volume you're in, and that will tell you which particles can potentially contribute to the density at that point. For instance if it's a uniform grid, you only need to look at the current cell and possibly some of its immediate neighbors if they're within the maximum particle radius of your current point.

Getting a smoothed mesh with surface normals etc. is usually done by defining a density field for the particle system. This is a function on 3D space that consists of summing up a falloff function for each particle over a defined radius: $$D(\vec r) = \sum_{i \in \text{particles}} f \left(\frac{|\vec r - \vec r_i|}{R_i} \right)$$ where $$f(u)$$ is some function that falls off smoothly from $$f(0) = 1$$ to $$f(1) = 0$$ (for example, a simple cubic falloff: $$f(u) = 1 - 3u^2 + 2u^3$$), and $$\vec r_i, R_i$$ are the position and radius of the $$i$$th particle. Then you set some threshold density value and let the surface be defined by all the points where $$D(\vec r) = \text{threshold}$$.

For getting the surface normal: the surface normal is just the (normalized) gradient vector of the density field. This can be calculated by finite differences between nearby points in the field: \begin{aligned} \nabla D(\vec r) &= \left( \frac{\partial D}{\partial x}, \frac{\partial D}{\partial y}, \frac{\partial D}{\partial z} \right) \\ &\approx \frac{1}{\epsilon} \left( D(\vec r + \epsilon \vec x) - D(\vec r), D(\vec r + \epsilon \vec y) - D(\vec r), D(\vec r + \epsilon \vec z) - D(\vec r) \right) \end{aligned} or analytically by adding up gradient contributions from each particle. \begin{aligned} \nabla D(\vec r) &= \sum_i \nabla \left( f \left(\frac{|\vec r - \vec r_i|}{R_i} \right) \right) \\ &= \sum_i \text{normalize}(\vec r - \vec r_i) \frac{1}{R_i} \frac{df}{du} \end{aligned}

For getting the surface in the form of a mesh: the keyword to search for is "isosurface extraction". Marching cubes is one example of such an algorithm, as lightxbulb mentioned. It evaluates the density field at evenly spaced grid points, uses those values to estimate where the intersection of the surface with the cubes' edges will be, and constructs from that a tiny mesh within each cube. Stitching all those meshes together yields a representation of the full surface. Naturally, the accuracy of it depends on the size of the grid used.

There is also a variant that uses tetahedra instead of cubes, and another similar algorithm called "dual contouring". Here's a nice article showing how these work interactively (in 2D). There's also a great deal of literature on how to improve accuracy, performance, etc. with these kinds of methods.

• Thanks for the great rundown! I still have some trouble seeing how I could easily turn a bunch of random SPH particles into a grid. However, that does seem like the way to go, because a grid allows me to implement marching cubes for the meshification. (The reason I don't do a grid based fluid sim in the first place is because particles allow for a great freedom of simulation domain) Nov 7 '20 at 20:26
• @AnnoyinC You define a density field based on the particles—meaning put some falloff function with a defined radius around each particle, and add up all the functions for all the particles. That gives you an analytic function you can evaluate at any point in space (using an acceleration structure to skip particles too far away to contribute). Then you plop down a grid and evaluate the function and its gradient at each grid point. Nov 7 '20 at 20:48
• So density(x,y,z) = "for each particle{add particle.density*fallof(x,y,z, particle.position)}". This doesn't sound very realtime? Do you know any codebases or papers which implement what you speak of? Nov 7 '20 at 22:55
• Yep. It's certainly possible to make this kind of thing work in realtime, but it's not going to be a trivial job. I'm sure there are SPH implementations out there that you could look at how they do their rendering, or libraries that can automate some of this, but it's not really my area so I don't know of any good examples. Nov 8 '20 at 0:31
• @AnnoyinC Yes $\vec r$ is position in space. $\epsilon$ is just some reasonably small number. $R_i$ is particle radius. $df/du$ is just the derivative of $f(u)$. If you're having trouble with these concepts, you might need to do some studying about vector calculus / multivariable calculus. Nov 10 '20 at 19:21