<p>I'm working on a fluid simulation project where two spheres are simulated using smoothed particle hydrodynamics (SPH). I have generated two sets of particles inside and on the surface of a sphere. You can see the result in this picture:</p>
<p><a href="https://i.sstatic.net/lOlW8.png" rel="nofollow noreferrer">scatter plot of spheres' particles</a></p>
<p>I am going to extract the surface of the particles by the marching cubes algorithm. Currently, I am using &quot;scikit-image&quot; implementation of marching cubes in python, and I got this result:</p>
<p><a href="https://i.sstatic.net/izV5Q.png" rel="nofollow noreferrer">Reconstructed surfaces and particles</a></p>
<p>In the above picture, two small spheres are scattered plots of my particles, and two large surfaces are the reconstructed surface that I got from the marching cube algorithm.</p>
<p>I have two problems:</p>
<p>1- As you can see, the surfaces are not smooth enough. So, I'm looking for an efficient method to make these surfaces smoother with a more sphere shape. The function that I used for surface reconstruction of my SHPs is:</p>
<pre><code>Phi(x) = - c + sum_{j in N(x)} of (V_j * W(x - x_j, h))
</code></pre>
<p>where <code>c</code> is the surface threshold, <code>x</code> is the position where we want to evaluate, <code>sum_{j in N(x)}</code> is the sum over all particles <code>j</code> in the neighborhood of point <code>x</code>, <code>V_j</code> is the volume of particle <code>j</code>, <code>W</code> is the SPH kernel function, <code>x_j</code> is the position of particle <code>j</code>, and <code>h</code> is the kernel smoothing length</p>
<p>2- The surfaces are supposed to exactly lie on particles as they are constructed to cover the outer surfaces of the particles. But as you can see in the picture, two surfaces are located somewhere far from the particles. I think my voxelization is not correct.</p>
<p>The code that I have written for surface reconstruction is:</p>
<pre><code>def surface_reconstruction(rezolution=20, surfacelevel=0, bounds=(-4, 4)):
    # First I creat my voxel grid
    x = np.linspace(bounds[0], bounds[1], rezolution)
    y = np.linspace(bounds[0], bounds[1], rezolution)
    z = np.linspace(bounds[0], bounds[1], rezolution)

    xx, yy, zz = np.meshgrid(x, y, z)
    surface = np.zeros((len(xx), len(yy), len(zz)))
    # Then I loop over each voxel's corner and find a value for that point to check whether 
    # it is inside the mesh or no?
    for i in range(len(xx)):
        for j in range(len(yy)):
            for k in range(len(zz)):
                for p1 in env.particles:
                    # p1.X is the position of particle p1. Here, I compute the distance 
                    # between each particle and voxel coordinate
                    dX = p1.X - [xx[i, j, k], yy[i, j, k], zz[i, j, k]]
                    r = np.linalg.norm(dX)
                    # p1.H is a kernel radius of particle p1
                    if r &lt;= p1.H:
                        # p1.vol is particle p1 volume
                        # Wij is a kernel function of SPH
                        surface[i, j, k] += -surfacelevel + p1.vol * Wij(r, p1.H)
   return surface
</code></pre>
<p>Then I apply marching cubes on &quot;surface&quot;</p>
<pre><code>verts, faces, normals, values = skimage.measure.marching_cubes(surface_reconstruction(), 0)
</code></pre>
<p>So, I am wondering</p>
<p>1- How can I make the surfaces smoother?</p>
<p>2- Why surfaces are larger than particles scale and far from particles?</p>

CC BY-SA 4.0
2022-04-23T07:57:45.187Z
Could be wrong but looks to me you compute distance dX between p1 and integer x,y,z which are not in the same space. ie: the grid may be 0-10, 0-500, 0-1000, who knows but your particle position p1 is definitely not in that space. This would cause the results you have shown. Also your marching cubes resolution is far too coarse.
Andrew Wilson
CC BY-SA 4.0
2022-04-24T08:36:08.337Z
@AndrewWilson, yeah, I think my voxelization is not correct. But I do not know how can I correct it
Hassan
<particles><fluid-sim>
How to extract surface of smoothed particle hydrodynamics smoothly?
Hassan
# How to extract surface of smoothed particle hydrodynamics smoothly?

(Asked by: username_0 on 2022-04-22)

I'm working on a fluid simulation project where two spheres are simulated using smoothed particle hydrodynamics (SPH). I have generated two sets of particles inside and on the surface of a sphere. You can see the result in this picture:

[scatter plot of spheres' particles](https://i.sstatic.net/lOlW8.png)

I am going to extract the surface of the particles by the marching cubes algorithm. Currently, I am using "scikit-image" implementation of marching cubes in python, and I got this result:

[Reconstructed surfaces and particles](https://i.sstatic.net/izV5Q.png)

In the above picture, two small spheres are scattered plots of my particles, and two large surfaces are the reconstructed surface that I got from the marching cube algorithm.

I have two problems:

1- As you can see, the surfaces are not smooth enough. So, I'm looking for an efficient method to make these surfaces smoother with a more sphere shape. The function that I used for surface reconstruction of my SHPs is:

```
Phi(x) = - c + sum_{j in N(x)} of (V_j * W(x - x_j, h))

```

where `c` is the surface threshold, `x` is the position where we want to evaluate, `sum_{j in N(x)}` is the sum over all particles `j` in the neighborhood of point `x`, `V_j` is the volume of particle `j`, `W` is the SPH kernel function, `x_j` is the position of particle `j`, and `h` is the kernel smoothing length

2- The surfaces are supposed to exactly lie on particles as they are constructed to cover the outer surfaces of the particles. But as you can see in the picture, two surfaces are located somewhere far from the particles. I think my voxelization is not correct.

The code that I have written for surface reconstruction is:

```
def surface_reconstruction(rezolution=20, surfacelevel=0, bounds=(-4, 4)):
    # First I creat my voxel grid
    x = np.linspace(bounds[0], bounds[1], rezolution)
    y = np.linspace(bounds[0], bounds[1], rezolution)
    z = np.linspace(bounds[0], bounds[1], rezolution)

    xx, yy, zz = np.meshgrid(x, y, z)
    surface = np.zeros((len(xx), len(yy), len(zz)))
    # Then I loop over each voxel's corner and find a value for that point to check whether 
    # it is inside the mesh or no?
    for i in range(len(xx)):
        for j in range(len(yy)):
            for k in range(len(zz)):
                for p1 in env.particles:
                    # p1.X is the position of particle p1. Here, I compute the distance 
                    # between each particle and voxel coordinate
                    dX = p1.X - [xx[i, j, k], yy[i, j, k], zz[i, j, k]]
                    r = np.linalg.norm(dX)
                    # p1.H is a kernel radius of particle p1
                    if r <= p1.H:
                        # p1.vol is particle p1 volume
                        # Wij is a kernel function of SPH
                        surface[i, j, k] += -surfacelevel + p1.vol * Wij(r, p1.H)
   return surface

```

Then I apply marching cubes on "surface"

```
verts, faces, normals, values = skimage.measure.marching_cubes(surface_reconstruction(), 0)

```

So, I am wondering

1- How can I make the surfaces smoother?

2- Why surfaces are larger than particles scale and far from particles?

**Tags:** particles, fluid-sim

> **Comments:**
> **username_1**: Could be wrong but looks to me you compute distance dX between p1 and integer x,y,z which are not in the same space. ie: the grid may be 0-10, 0-500, 0-1000, who knows but your particle position p1 is definitely not in that space. This would cause the results you have shown. Also your marching cubes resolution is far too coarse.  
> **username_0**: @username_1, yeah, I think my voxelization is not correct. But I do not know how can I correct it  

---

## 0 answers

