So I have been working on a Smooth Particle Hydrodynamics (SPH) simulation and I am trying to implement the Marching Cubes (MC) algorithm to visualize my SPH output. I understand the MC algorithm well enough from reading about it, however usually it is discussed assuming that you already have some scalar field which you can than use to classify your voxels as being either above, below, or on the boundary surface you are trying to create (in my case this amounts to checking if the voxel is above, below, or on the fluid surface). The key part here is that it is easy to check the value of the scalar field at any cube vertex.

However, with my SPH output data I don't have some nice scalar field to use. Instead I have a lot of particles distributed (randomly or haphazardly) in my domain. I know each particles density at these random points, but I don't know the density at the MC grid vertices. For example consider the simple 2d example:

enter image description here

Each of the blue points represent a particle. After performing my SPH solve, I know the position, velocity, and density of each particle at these points. Now, in order to classify each voxel as being above, below, or on the fluid surface I need to know the density at the voxel vertices (i.e. the red points). How is this typically done given you know the density at the blue particles?

My initial thought was to calculate the density at the red points in the same way as is done in the SPH simulation. In this case I would use:

$\rho_{i}(r) = \Sigma_{j}m_{j}W(r -r_{j}, h)$

at each voxel vertex point where $W=0$ when $r > h$. In the picture below this means the green particles would be used in calculating the voxel vertex density value:

enter image description here

Now my concerns/doubts with this approach are as follows:

  1. The $h$ value we use here is probably smaller than the $h$ value we used in the actual SPH simulation because we want a refined voxel grid to allow for refined surface reconstruction.
  2. With the above point in mind, even with a modest $64\times{64}\times{64}$ voxel grid in 3D, this results in over 2 million vertex points (64*64*64*8).
  3. Calculating the density in the actual SPH simulation this way is somewhat expensive, but at least we might only be doing it for O(100k) particles. Here though as just mentioned above we would need to do it for more than 2 million. This I fear would surpass the actual simulation time.

So my questions are as follows:

  1. Is this approach reasonable/correct? If not how should I classify my voxels?
  2. If it is correct, than would you use just a large array of size 64*64*64*8 to store each density value at every voxel vertex?

Note: As mentioned in the title I am using CUDA so keep that in mind.

  • $\begingroup$ What would be the "surface reconstruction" for the example that you present above? Would it have a single interface or the small cluster that is higher would have it's own surface? $\endgroup$
    – nicoguaro
    Mar 27, 2017 at 15:20
  • $\begingroup$ @nicoguaro In the above 2d case it would be an array of line segments that separate the fluid particles from the empty space region. In the case of 3d, it would be composed of triangles that together when rendered form a 2d surface. Also, yes, the small cluster would have its own "surface", i.e. maybe some lines that form a crude circle around the three particles. In particular my surface here is not a bijection if that's causing confusion. $\endgroup$
    – James
    Mar 27, 2017 at 18:20
  • $\begingroup$ @nicoguaro As long as there were voxel vertices that did not get classified as being inside the fluid, than the three particles would form a separate group with there own surface. This allows for fluid sloshing and droplets flying around. $\endgroup$
    – James
    Mar 27, 2017 at 18:26
  • $\begingroup$ What about using a classifier algorithm and then interpolating that result? Some examples using scikit-learn here. $\endgroup$
    – nicoguaro
    Mar 27, 2017 at 18:38

1 Answer 1


It is a different approach but maybe you could do a point cloud surface reconstruction using total variation denoising as in here. They have a demo code that you could use as a base for your implementation, you can get it here. The paper also shows other alternatives, the Poisson surface reconstruction also shows nice results. The poisson surface reconstruction algorithm is implemented in Meshlab.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.