How to quantify the volume of a Volume

[ttsesm]

Hi @marcomusy,

I have the following question, for meshes it is possible to quantify the volume of a mesh by calling the volume() function. Though as I understood this is dependent to the triangulation of the meshes if I am not wrong.

My ultimate goal was to quantify the result of the "intersection" boolean operation between two meshes. But I've noticed that this can become quite slow, and also I know that in vtk the boolean operations are on the weak side of the library. Thus, I thought to do the same using the volumes of two objects. However, I couldn't figure out how to quantify the output.

So for example let's say that we have two boxes:

from vedo import Box, show
vol1 = Box(size=(35,10, 5)).binarize()
vol2 = Box(size=( 5,10,35)).binarize()

vol = vol1.operation("and", vol2)
print(vol)
show(vol1.c('g'), vol2.c('r'), vol.c('magenta'), interactive=True).close()

How can I quantify the vol which is a class Volume(), there is no .volume() function to call on a volume. Then I thought that I can transform it back to mesh by using .tomesh() but this seems to be slow and then calling .volume() gives 0, which I figured out that I can fix by calling the triangulate(). However, the output from that also looks quite high, and then I realized that the volume() output on a mesh probably is related to its triangulation.

Any idea is welcome, thanks.

[marcomusy]

Well I have not tested this but you can try with some specific configuration:

from vedo import *
vol1 = Box(size=(35,10, 5)).binarize()
vol2 = Box(size=( 5,10,35)).binarize()
vol = vol1.operation("and", vol2)
dx, dy, dz = vol.spacing() # voxel size
counts = np.unique(vol.pointdata[0], return_counts=True)
n0, n1 = counts[1]
vol_value = dx*dy*dz * n1
vol1.cmap('g')
vol2.cmap('r')
vol.cmap('m')
show(vol1, vol2, vol, f"volume = {vol_value}", axes=1).close()

[ttsesm]

Hi @marcomusy,

Lovely, this works fine. The point is now I want to apply this on two point clouds.

However, the issue that I face now is that when I try to extract the volume from the point clouds directly I get nothing. I tried the different parameters to the .tovolume() function but it doesn't seem to give me something that I can use:

pcd1.tovolume(kernel='shepard', n=4)
<vedo.volume.Volume object at 0x7fc247ea2800>
pcd1.pointdata
Point Data is empty.

The same with pcd2. Then I tried first to make them as meshes by using the .reconstruct_surface() method and use .binarize() on them and then extract overlap but the reconstructed meshes are really bad.

Any idea how to handle this.

Thanks. pcds.zip

[marcomusy]

Hi, what about this alternative solution

from vedo import *

pcd1 = Points("pcd1.ply").color("blue5")
pcd2 = Points("pcd2.ply").color("red5")

ug1 = pcd1.generate_delaunay3d(radius=0.01)
surf1 = ug1.tomesh().compute_normals()
surf1.color("blue5").alpha(0.1)
ug2 = pcd2.generate_delaunay3d(radius=0.01)
surf2 = ug2.tomesh().compute_normals()
surf2.color("red5").alpha(0.1)

s12 = surf1.boolean("intersect", surf2)
s12.color("green5").lw(1)
vol = s12.volume()

show(pcd1, pcd2, surf1, surf2, s12, f"volume = {vol}", axes=1)

PS : creating a volume from a pointcloud is not possible because there is no face normal which defines the "inside" from the "outside".

[ttsesm]

Using the meshes is good as alternative the only issue that I have though with this approach is that radius needs to change each time because I am loading different point clouds each time so using 0.01 as a fixed number it doesn't always work, e.g.:

that's why I wanted to go volumes instead.

Do you have any idea how I can dynamically specify the radius value for the triangulation, depending on the point cloud that I have each time... :thinking:

[marcomusy]

You can automatize the radius value by histogramming the point relative distances, something like

from vedo import *
from vedo.pyplot import histogram

pcd1 = Points("pcd1.ply").color("blue5")
pcd2 = Points("pcd2.ply").color("red5")

dists1 = []
for p1 in pcd1.coordinates:
    q1 = pcd1.closest_point(p1, n=2)[1]
    dists1.append(mag(p1 - q1))
histo1 = histogram(dists1, bins=25).clone2d()
radius = histo1.mean * 10

tetm1 = pcd1.generate_delaunay3d(radius=radius)
surf1 = tetm1.tomesh().compute_normals()
surf1.color("blue5").alpha(0.1)

tetm2 = pcd2.generate_delaunay3d(radius=radius)
surf2 = tetm2.tomesh().compute_normals()
surf2.color("red5").alpha(0.1)

s12 = surf1.boolean("intersect", surf2)
s12.color("green5").lw(1)
print(s12.volume())

show(pcd1, pcd2, surf1, surf2, s12, histo1, axes=1)

[ttsesm]

Thanks @marcomusy, as always really helpful 👍 For the moment, it seems to do the job so I will be closing this issue.