Classifying Marching Cube voxels from SPH output data using CUDA

Question

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:

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:

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.

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.