Volume of a box is zero and it's intersection with a sphere is considered an empty set

[lukablagoje]

Hello,

I'm creating a box within a sphere with the same center:

s0 = Sphere((1,1,1),r=1, res=50,alpha=0.5) b0 = Box((1,1,1),length=1, width=1, height=1,alpha =0.5)

But when I compute the volume of the b0.volume(), I obtain 0. Moreover, when I try Boolean intersection:

i1 = s0.boolean('intersection',b0)

The object i1 doesn't show anything (empty set).

What could be the issue here? Thanks

[marcomusy]

You need to use triangulate, and the two shapes must intersect somewhere, e.g.:

from vedo import *

settings.useDepthPeeling = True

s0 = Sphere((1,1,1), r=1, res=50, alpha=0.5)
b0 = Box((1,1,1), length=3, width=1, height=1, alpha =0.5).triangulate()

i1 = s0.boolean('intersect',b0)
i1.flat()

show([(s0, b0), i1], N=2, axes=1, zoom=0.8)

[lukablagoje]

Thanks for the reply!

I have an additional question which is related. If I want to work with solid spheres and cubes, I assume that I have to convert these meshes to volumes, right? How would I compute and obtain an intersection then (between a solid sphere/ball and a solid cube)?

[marcomusy]

what is the difference in intersecting the surface or a "solid" object? by solid you mean a tetrahedral mesh?

[lukablagoje]

By intersecting a solid object, I mean that I obtain the inside part of the object during the intersection. I'll give an illustrative example, If I had a surface of the human body (skin), I would like to intersect for example a sphere and the upper torso part, in order to obtain the skin and the organs inside (as some data type that I can access, compute it's volume, etc.) with high accuracy.

Is this achievable? How would you approach this?

[marcomusy]

it very much depends on the exact structure of your dataset: is your data made of polygons (a surface mesh) or is it volumetric (tetrahedral mesh or volume)? depending on this there are various examples: https://github.com/marcomusy/vedo/blob/master/examples/advanced/cutWithMesh1.py https://github.com/marcomusy/vedo/blob/master/examples/advanced/cutWithMesh2.py https://github.com/marcomusy/vedo/blob/master/examples/volumetric/tet_cutMesh1.py https://github.com/marcomusy/vedo/blob/master/examples/volumetric/tet_cutMesh2.py https://github.com/marcomusy/vedo/blob/master/examples/volumetric/tet_threshold.py

[lukablagoje]

Thank you for being helpful and responsive. I'll try to describe my project in a nutshell :

• I am trying to determine what's a good way to approach generating and analyzing 3D objects in python. It seems that this library has many required tools that I need.
• I am working with 3D porous data (regions of empty space and objects). Some of them I need to generate myself, for example, a lattice-like structure composed of connected lines and spheres. Some of them will be real datasets of neurons. But the point is that they need to be filled (Instead of being just surfaces).
• At the moment (but not confined to this), I am trying to calculate the fractal dimension of these lattices that I'm creating from spheres and lines. To accomplish this, I am generating a grid of cubes and trying to find their intersections with the lattice.

From what I've gathered and your answers:

• When I generate a sphere, I can "fill" it by converting it to a tetmesh. Afterward, I can find the intersection between that tetmesh and some other mesh (for example a cube). Is this the correct way to approach this problem?

[marcomusy]

Hi, it very much depends on the actual information you want/need to encode into your data. That defines whether you should use 2d or 3d data. If you only deal with tubes (or even 1d lines) and spheres with simple or none internal structure then i would go for 2d polygonal surfaces, otherwise i would go for for 3d tet meshes.

In 2d you can encode information on the surface of the polygons, in 3d you can encode information for all points inside the geometry.

[lukablagoje]

Yes, I see. So If I use volumes (3D data). If I may ask one last question for my initial problem of the sphere and a cube, with the improved perspective. GIven v_1 and v_2:

---

s = Sphere((1,1,1),r=1, res=50,alpha=0.5)

v_1= mesh2Volume(s, spacing=(0.02, 0.02, 0.02)).alpha([0,0.5]).c('blue')

c = Cube((0+1.5,0+1.5,0+1.5),side=1,c=1 ,alpha =0.5)

v_2 = mesh2Volume(c, spacing=(0.02, 0.02, 0.02)).alpha([0,0.5]).c('red')

---

1. How can I find the intersection of v_1 and v_2 (I tried volume operations but haven't succeeded)
2. Given the intersection, how do I access/compute its volume?

Thanks for answering these basic questions, it's my first time working with these kinds of objects.

[marcomusy]

You are welcome to ask as many questions you like :)

What is called Volume in vedo is a dataset which is made out of voxels - the equivalent of pixels in 2d. Think of them as a set of little cubes subdividing a larger cube. Each of these small voxels can contain one or more numbers. Then you have the TetMesh class which fills a region of space of arbitrary shape (not necessarily a cube), made of tetrhaedrons. Each tet can again contain one or more numbers.

What the function mesh2Volume does is to create a Volume from a Mesh (polygonal surface obj). Probably not what you want..

You can have a tetmesh and then cut it with a surface to select a region of interest, whose volume would be the volume of its closed surface (similar to example tet_cutMesh1).

[lukablagoje]

Yes, I see. In that case, I have further simple questions:

1. Given vedo.Volume, vedo.Mesh, vedo.Tetmesh, how can I determine the bounding box for each? Is there a way to access it from the class data?

2. I know that I can compute mesh.volume(). Is same thing possible for vedo.Volume and vedo.Tetmesh?

3. How do I access 3D grid values of voxels of vedo.Volume, for example as 3 dimensional matrix?

4. Is it possible to do intersections with vedo.Volume, to see if the two objects overlap?

5. When I want to convert from Tetmesh to Mesh, I just use tomesh() function. But are there more functions for this? For example for these cases: a)Tetmesh to Vedo.Volume or b)from Vedo.Volume to Tetmesh? c) from Mesh to Tetmesh (I know about delaunay3D, but is there any other way)

Thanks

[marcomusy]

1. yes - all classes have a method called bounds(), see this small example:

   from vedo import *
   disc = Disc(r1=1, r2=1.2)
   mesh = disc.extrude(3, res=50).lineWidth(1)
   mesh.cutWithCylinder([0,0,2], r=0.4, axis='y', invert=True)
   printc("[xmin,xmax, ymin,ymax, zmin,zmax] =", mesh.bounds(), c='yellow')
   show(mesh, mesh.box(), axes=0)

2. in the case of Volume that would be trivial: the volume is always the product of the sides. For TetMesh you can extract the outer polygonal surface and then compute the volume:

   from vedo import *
   tm = TetMesh(dataurl+'limb_ugrid.vtk')
   surf = tm.tomesh(fill=False)
   v = surf.volume()
   printc("volume is", v)
   show(tm, surf, N=2)

3. do vol.getPointArray() or vol.getDataArray()

4. No! that's not how Volumes work, remember that the shape of a Volume is always a box-shape.

5. a) No! see above. b) yes in principle but I never tried that. c) only delauney if mesh is convex, If it's not you may need to use other software to create tetrahedrons inside your surface boundary.

[astroboylrx]

@marcomusy @lukablagoje I am also trying to find the intersection volume of a sphere and a cube. I followed the example in a previous reply

You need to use triangulate, and the two shapes must intersect somewhere, blah blah

However, boolean sometimes work, sometimes boolean hangs forever.

In [1]: import vedo # I also import numpy by default

In [2]: pos = np.array([3.92578125e-01, 4.0234375e-02, 3.9062500e-04])
   ...: r = 0.0006805476689621596
   ...: dx = 0.00078125
   ...: c1 = vedo.Sphere(pos, r=r, res=50, alpha=0.5)
   ...: c2 = vedo.Cube(pos+np.array([dx, 0, 0]), side=dx, alpha=0.5).triangulate().subdivide(2, method=1)

In [3]: cc = c1.boolean("intersect", c2)
   ...: cc.volume()
Out[3]: 1.2907800195096867e-10

In [4]: c3 = vedo.Cube(pos+np.array([dx, dx, 0]), side=dx, alpha=0.5).triangulate().subdivide(2, method=1)

In [5]: cc = c1.boolean("intersect", c3)
# hangs forever

Note the only difference is I shifted c3 with [dx, dx, 0] instead of [dx, 0, 0]. In case you wonder how they look like, the picture below shows c1 and c3. They do intersect with each other, but boolean stucks at somewhere. I tested this script on my Mac and two Linux machines, all with the latest possible version and dependencies (i.e., via pip install -I vedo), boolean always stucks.

I don't quite understand what happened. Any suggestions/ideas would be greatly appreciated!

[marcomusy]

Hi @astroboylrx thanks for reporting, indeed the boolean operation is quite problematic in many cases... it doesn't work very well.. it seems that it is sensitive to the absolute size of the objects involved, e.g.

import numpy as np
import vedo

# pos = np.array([3.92578125e-01, 4.0234375e-02, 3.9062500e-04])
# r = 0.0006805476689621596
# dx = 0.00078125

pos = np.array([0,0,0])
r = 0.2
dx = 0.2

c1 = vedo.Sphere(pos, r=r, res=50, alpha=0.75).lw(1)
c2 = vedo.Cube(pos + [dx, 0, 0], side=dx, alpha=0.75).lw(1)
c2.triangulate().subdivide(2, method=1)

cc = c1.boolean("intersect", c2)
print(cc.volume())

c3 = vedo.Cube(pos+ [dx, dx, 0], side=dx, alpha=0.75).lw(1)
c3.triangulate().subdivide(5, method=1)

cc = c1.boolean("intersect", c3)
print(cc.volume())

vedo.show([[c1, c3], cc], N=2, axes=1)

generates a meaningful result:

So the trick here could be that you put everything back a t the origin and scale to be about 1 before computing the volumes.

(see also https://discourse.vtk.org/t/compute-the-volume-of-the-intersection-between-two-vtkpolydata-objects/8338/5 )

[astroboylrx]

@marcomusy Many thanks for the fast response! I see. Rescaling everything does work as expected. So the issue is due to numerical precision? I wonder why the function hangs like in an endless loop?

[marcomusy]

So the issue is due to numerical precision?

yes, I think that's what happens..

I wonder why the function hangs like in an endless loop?

i've got no good answer there... the vtk boolean operation is know to be a weak point of the vtk library, which might be fixed in the future as independent packages already exist that can be included in vtk.