Question about Sparse Density Map Marching Cube

[losanc]

Hi,

I have a question about the marching cube step when using sparse density map. In the code here

            // We want to find edges that cross the iso-surface,
            // therefore we can choose to either skip all points above or below the threshold.
            //
            // In most scenes, the sparse density map should contain more entries above than
            // below the threshold, as it contains the whole fluid interior, whereas areas completely
            // devoid of fluid are not part of the density map.
            //
            // Skip points with densities above the threshold to improve efficiency
            if point_value > iso_surface_threshold {
                return;
            }

So it skips all the points with larger values. For points with lower values, it finds the neighbor, and check if it cross the threshold.
But for edged like this

0----------value larger than threshold

So, one side of the edge has value 0. The check will be skipped, because it's not in the density map. On the other side has large enough value, and it will be skipped as well.

In the end, this triangle vertex will be skipped, though it should be a vertex here

I think this is a rare case, but with very small threshold, this can happen

[w1th0utnam3]

Hi, thanks for the message. I think in my implementation of the surface reconstruction pipeline, this case should not happen, as I evaluate the particle kernels with a radius that is basically cellsize * ceil(support_radius/cellsize) (see here). Thus, there should always be a density = 0 ring around the fluid. I guess the comments and documentation of the marching cubes function should be updated to highlight that this is expected.

[losanc]

yes. it does have a density=0 ring. But the grid with density=0 may not explicitly inside the density_map.

For example, if I have only 1 particle at the position (0.5, 0, 0 ), and i want the grid which is all axes are integers. so, cellsize = 1 in this case. If we take support_radius equals to 1 as well. then the actual evaluated radius cellsize * ceil(support_radius/cellsize) is 1 as well.

I will use F(x,y,z) means the value at grid point(x,y,z). So obviously F(0,0,0) has a non zero value. For the point (-1, 0,0), at this line, r_squared < self.kernel_evaluation_radius_sq will returns false. so

*sparse_densities
 .entry(flat_point_index)
.or_insert(R::zero()) += density_contribution;

this part will be skipped for point(-1,0,0)

In the end, sparse_densities contains value at F(0,0,0), but not F(-1,0,0), which is implicitly 0.

If the threshold value is between F(0,0,0) and 0. It won't create a vertex between (-1,0,0) and (0,0,0), but it should. Because when looping over the sparse_densities, F(0,0,0) is skipped by if point_value > iso_surface_threshold , and F(-1,0,0) is not added into sparse_densities

[w1th0utnam3]

Right, I see your point. Even though the self.kernel_evaluation_radius_sq contains a small epsilon, this small epsilon might not be enough to insert an explicit 0 value into the density map in all cases. This can trigger a bug of missing isosurface vertices in your scenario. I think there are two ways to fix this:

1. Reverse the if point_value > iso_surface_threshold { return; } check that you posted in your original message (and adapt the rest of the function accordingly). This way, we could handle cases where edges have a value greater than the threshold on one side and no value at all on the other side. However, this would mean that we would always enter the inner loop of the isosurface vertex creation for all vertices that are inside of the volume of the fluid until we find out that there is no neighbor that is below the threshold. This could be more expensive for scenes with a large interior volume.
2. Ensure that an explicit 0 ring is always added to the density map. This could be done by adjusting the kernel_evaluation_radius to kernel_evaluation_radius = cube_size * (half_supported_cells_real * (R::one() + R::default_epsilon().sqrt()) + R::one()) i.e. adding a 1 x cube_size ring where the kernel is evaluated. This should be the maximum required radius offset to always write an explicit zero if I'm not mistaken? So this will add quite some amount of entries to the density map, proportional to the surface area of the fluid. In terms of performance this will also increase the number of evaluations of the density field for all particles, even in the interior of the fluid.

I'm not entirely sure which approach is preferable in general... Usually I would assume that there are more points of the marching cubes grid than particles to get a nice mesh. But it will also change for scenes with a lot of separated/isolated particles vs. a cube of fluid.

In any case, we should add some tests for this.

[w1th0utnam3]

For now I decided to use the first approach: 919d3dc

Note that the parallel reconstruction using domain decomposition was never affected by this problem, only the "global" reconstruction.

[losanc]

I think both fix works, but it will depend on the scene setup to say which one is better.