We have a 3-dimensional grid of voxels (or cells), with individual voxels being of volume $dx\,dy\,dz$ where $dx=dy=dz=1$.

A cone-like surface is defined by some function, $z = f(x, y)$, which in this case is specifically (where $\epsilon$, $w_0$ and $r_0$ are constants):


Over the entire grid, how does one compute the fraction of each voxel's volume lying within the volume enclosed by that cone-like surface?

Obviously values should range from $0\rightarrow1$ and, in this particular case, voxels positioned at large values of $x$ and $y$ and small values of $z$ have none of their volume within the cone, while those lying at large values of $z$, and small values of $x$ and $y$ are completely enclosed by the cone-like surface.

  • $\begingroup$ It's a task for Marching Cube algorithm. Calculate this field in your voxelized grid: $$\mathcal{f}(x,y,z) = z - \Bigg ( \frac{x^{2}+y^{2}}{w_{0}^{2} r_{0}^{-2 \epsilon}} \Bigg )^{\frac{1}{2 \epsilon}}$$ and finally find the surface that belongs to $\mathcal{f}(x,y,z) = 0$. $\endgroup$ – Alone Programmer Apr 1 '20 at 18:54
  • $\begingroup$ According to the MC algorithm, the field defines which voxels lie completely in ($ f(x,y,z)>0 $ for every vertex) and which lie completely out ($f(x,y,z)<0$ for every vertex) of the surface. For the other voxels, a defined grid of $2^8$ 'approximate' surfaces corresponding to different vertex in/out combinations allows for some degree of approximation of fraction of a voxel within the surface. That much I understand, though computing the voxel fractions for each of those possibilities might take a while. I did not however understand your comment to 'find the surface belonging to $f(x,y,z)=0$'? $\endgroup$ – simonp2207 Apr 3 '20 at 9:22

That function looks pretty well behaved/smooth. So assuming your voxels are pretty small compared to the surface curvature, I think you could get away with a signed-distance-like approach, by testing the 8 corners against the surface. This will quickly eliminate the all-corners-inside (1) and all-corners-outside (0) cases.

Near the surface there will be straddling voxels for which this test is inconclusive (corners are a mix of inside and outside). For these you would recurse, bisecting them (in 3D) into eight subvoxels and then sum the results. Of course, some of these subvoxels will still straddle, so recurse again, perhaps up to some fixed/maximum recursion depth. If a subvoxel still straddles all the way down at the leaf of the recursion, give up and call it 0.5.

This is basically like bisection search to find the root of a function, but in 3D to find the "root" of a surface.

This might be a bit slow/inaccurate (kinda depends upon how many voxels you have and what kind of recursion depth you pick), but it looks easy to implement without pulling in a lot of computational geometry. A possible refinement would be to pull in a little computational geometry, and approximate the volume of leaf subvoxels by cutting them with the local tangent plane of the surface into a faceted polyhedron, then compute its volume (probably by breaking it into tetrahedra, which has an easy volume formula). But I would only do this if I was unhappy with the simpler method after I implemented it first.

  • $\begingroup$ Typically my grid will likely contain $\sim10^8-10^9$ voxels, so your AMR-sounding approach may not be too computationally expensive considering it would likely only be computed the once. Thanks for your comment! $\endgroup$ – simonp2207 Apr 3 '20 at 9:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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