<p>For culling faraway particles: as lightxbulb said, the keyword to search for is &quot;acceleration structure&quot;. 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.</p>
<p>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.</p>
<p>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:
<span class="math-container">$$
D(\vec r) = \sum_{i \in \text{particles}} f \left(\frac{|\vec r - \vec r_i|}{R_i} \right)
$$</span>
where <span class="math-container">$f(u)$</span> is some function that falls off smoothly from <span class="math-container">$f(0) = 1$</span> to <span class="math-container">$f(1) = 0$</span> (for example, a simple cubic falloff: <span class="math-container">$f(u) = 1 - 3u^2 + 2u^3$</span>), and <span class="math-container">$\vec r_i, R_i$</span> are the position and radius of the <span class="math-container">$i$</span>th particle. Then you set some threshold density value and let the surface be defined by all the points where <span class="math-container">$D(\vec r) = \text{threshold}$</span>.</p>
<p>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:
<span class="math-container">$$
\begin{aligned}
\nabla D(\vec r) &amp;= \left( \frac{\partial D}{\partial x}, \frac{\partial D}{\partial y}, \frac{\partial D}{\partial z} \right) \\
&amp;\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}
$$</span>
or analytically by adding up gradient contributions from each particle.
<span class="math-container">$$
\begin{aligned}
\nabla D(\vec r) &amp;= \sum_i \nabla \left( f \left(\frac{|\vec r - \vec r_i|}{R_i} \right) \right) \\
&amp;= \sum_i \text{normalize}(\vec r - \vec r_i) \frac{1}{R_i} \frac{df}{du}
\end{aligned}
$$</span></p>
<p>For getting the surface in the form of a mesh: the keyword to search for is &quot;isosurface extraction&quot;. 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.</p>
<p>There is also a variant that uses tetahedra instead of cubes, and another similar algorithm called &quot;dual contouring&quot;. <a href="https://wordsandbuttons.online/interactive_explanation_of_marching_cubes_and_dual_contouring.html" rel="nofollow noreferrer">Here's a nice article</a> 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.</p>

CC BY-SA 4.0
2020-11-07T20:26:20.870Z
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)
AnnoyinC
CC BY-SA 4.0
2020-11-07T20:48:49.017Z
@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.
Nathan Reed
CC BY-SA 4.0
2020-11-07T22:55:30.430Z
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?
AnnoyinC
CC BY-SA 4.0
2020-11-08T00:31:41.900Z
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.
Nathan Reed
CC BY-SA 4.0
2020-11-10T19:10:29.783Z
@AnnoyinC I added some more explanation of the density field and the gradient.
Nathan Reed
CC BY-SA 4.0
2020-11-10T19:19:30.677Z
Thanks! I assume "vec{r}" means a position vector, and I assume epsilon is just like a "h" in integrals. What is capital R? And lastly, in the "sum each i: normalize(r-r_i)1/R_i df/du", I feel like I do not understand the math behind this yet. Is there a textbook that handles this? df/du means that you derive the entire sum over space relative to fallof f(u)? I have some trouble imagining what it culminates in. It would help me out greatly if you could point me to a class that handles this.
AnnoyinC
CC BY-SA 4.0
2020-11-10T19:21:59.703Z
@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.
Nathan Reed
CC BY-SA 4.0
2020-11-10T19:23:04.233Z
Thank you, I appreciate you taking the time to help me.
AnnoyinC
CC BY-SA 4.0
2020-11-07T04:51:15.923Z
2020-11-10T19:22:29.533Z
2020-11-10T19:22:29.533Z
Nathan Reed
<p>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.</p>
<p>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?</p>
<p>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.</p>
<ol>
<li>What are some algorithms for spatially culling far away particles with no local contribution?</li>
<li>What are some algorithms for getting the surface normal of a body of water made of particles?</li>
<li>What are some algorithms for getting the surface of a body of water made of particles, in the form of a mesh?</li>
</ol>

CC BY-SA 4.0
2020-11-07T00:29:25.757Z
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.
lightxbulb
CC BY-SA 4.0
2020-11-07T00:36:42.207Z
@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?
AnnoyinC
CC BY-SA 4.0
2020-11-07T10:29:28.153Z
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.
lightxbulb
CC BY-SA 4.0
2020-11-09T12:04:32.837Z
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.
pmw1234
<rasterizer><raymarching><particles><fluid-sim>
Rendering a particle system as a smooth blob
AnnoyinC
# Rendering a particle system as a smooth blob

(Asked by: username_0 on 2020-11-06)

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?

**Tags:** rasterizer, raymarching, particles, fluid-sim

> **Comments:**
> **username_1**: 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.  
> **username_0**: @username_1 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?  
> **username_1**: 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.  
> **username_2**: 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.  

---

## 1 answers

### Answer by username_3 on 2020-11-07. Score: 3 | Accepted answer

For culling faraway particles: as username_1 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 username_1 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](https://wordsandbuttons.online/interactive_explanation_of_marching_cubes_and_dual_contouring.html) 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.

> **Comments:**
> **username_0**: 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)  
> **username_3**: @username_0 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.  
> **username_0**: 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?  
> **username_3**: 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.  
> **username_3**: @username_0 I added some more explanation of the density field and the gradient.  
> **username_0**: Thanks! I assume "vec{r}" means a position vector, and I assume epsilon is just like a "h" in integrals. What is capital R? And lastly, in the "sum each i: normalize(r-r_i)1/R_i df/du", I feel like I do not understand the math behind this yet. Is there a textbook that handles this? df/du means that you derive the entire sum over space relative to fallof f(u)? I have some trouble imagining what it culminates in. It would help me out greatly if you could point me to a class that handles this.  
> **username_3**: @username_0 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.  
> **username_0**: Thank you, I appreciate you taking the time to help me.  

---

