CGAL - find points inside polyhedron

Question

I'm working in some code using the CGAL library to do geometric stuff.

I have a Polyhedron object and using it I create a mesh domain, following the code found in [the documentation][2].

Using the method CGAL::make_mesh_3(), I create my mesh, which is an object of this [type][3]

Next I need an normal 3d array which is going to represent the space in which the mesh is. I need to go point by point in the array (i,j,k) and decide if the point with these coordinates (after i,j,k are transformed into the right units) is INSIDE the domain defined by the polyhedron.

So what I do is use the method locate[4] from the object Triangulation, to find where the point is. This method returns a Cell_Handle[5] and a Locate_type (from which the location of the point could be know)

The triangulation class has other methods like is_cell() and is_infinite() ...

The point is that if i set my query to know if the point is inside the polyhedron in something like:

if( tr.is_cell(ch)  && !tr.is_infinite(ch) ) {
    return true;
  }else {
    return false;
  }

ch is the Cell_handle returned by the query locate.

It works as long as the polyhedron isn't a complex object. For simple objects (e.g. a cone, sphere etc), I get a pretty output like:

Cone

but if I try a object like a torus, it doesn't work and I get that the hole in the torus is inside the mesh (the points inside the torus give "true" as a result), in the if statement I pasted. But they should not be. I've tried other statements to filter it like:

if ( lt == Tr::INSIDE_CONVEX_HULL && !tr.is_infinite(ch)) {
    return true;
  }

but the result is the same. If the object has holes or the polyhedron is not equal to the CONVEX_HULL, I don't get what I what. Thinking a little, the CONVEX_HULL is, in the best case, what I get.

So how can I decide if a point (x,y,z) lies inside the polyhedron, or inside the mesh?

Is there a way to get this?

Answer 1

If you have located your point inside the mesh, then you have a cell handle containing the point. You can use the member function of the C3T3 is_in_complex(Cell_handle) to know whether this cell is inside or outside the domain.

Note that if you want to directly work on the polyhedron, you can use the functor Side_of_triangle_mesh