CGAL users discussion list

[Stephane Tayeb <>]

Hi,


 wrote:
> Hi,
> I had to compute the intersection points of a ray and a 3D object
> represented as a 3D complex (for example a tetrahedron mesh in 3D mesh
> generation). So I combined a 3D mesh generation procedure with an AABB
> tree. My query object is a line with fixed direction and crossing a vertex
> of the mesh (I consider parallel rays crossing each vertex for different
> intersection queries) and my set of primitives is a list of triangles
> belonging the 3D complex.
> I suppose, for example, to use the mesh generated in examples
> /Mesh_3/mesh_implicit_sphere.cp
> p (attached code). As I need to consider intersections with both triangles
> on the volume and on the surface I chose an iterator on triangulation
> facets using C3T3 and then I checked whether triangles belong to the 
> complex.
> Is this an efficient procedure to perform intersections with a tetrahedron
> mesh or Polyedron is recommended (see
> examples/AABB_tree/AABB_polyhedron_facet_intersection_example.cpp)?
> Actually I'm wondering about this as I obtain on average only 30
> intersection points for about 30000 triangles. Am I wrong expecting an
> higher number of intersections?

I think 30 is ok. If you have a look at the attached picture (which 
corresponds to the mesh you compute in your code) and mentally draw a 
line through the sphere, 30 intersections seems to be possible. You may 
try to get a more refined mesh to see how the number of intersection 
increases.

> As I need to follow the ray from the "entrance" to the "exit" point in the
> object I sort the list of intersections obtained by means of the tree. The
> question is:
> How can I  keep  the information on the corresponding crossed triangles
> after sorting?
> (without a new location).
> I mean obtain a list of pairs intersection/crossed triangle (vertexes)
> with the new order which starts with the intersection at the object
> boundary and ends when the ray exits from the object. Or, in other words,
> how control the entrance/exit facets of each crossed tetrahedron for each
> couple of consecutives  intersections?
> Thanks in advance for any suggestions!

I'm not sure to really understand what you want to do.

Why don't you store a container of std::pair<Point,Iterator> (instead of 
a Point container) ? Object_and_primitive_id contains the intersection 
Object (oid.first) and also the identifier of the intersected primitive 
(type Iterator in your code).

You could also use a different kind of AABB_primitive, such as Datum is 
Triangle_3 and Id is Tr::Facet_finite_iterator (or Tr::Facet). Thus, you 
can get back the triangulation facet from the intersected 
Object_and_primitive_id. See attached file (I didn't test it, it's just 
to give you an hint).

Hope this helps,
Stéphane.



> #include <list>
> #include <CGAL/AABB_tree.h>
> #include <CGAL/AABB_traits.h>
> #include <CGAL/AABB_triangle_primitive.h>
> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
> #include <CGAL/Mesh_triangulation_3.h>
> #include <CGAL/Mesh_complex_3_in_triangulation_3.h>
> #include <CGAL/Mesh_criteria_3.h>
> #include <CGAL/Implicit_mesh_domain_3.h>
> #include <CGAL/make_mesh_3.h>
> #include <CGAL/iterator.h>
> #include <iostream>
> #include <fstream>
> #include <map>
> #include <vector>
> #include <string>
> #include <CGAL/utility.h>
>
> // Domain
> struct K: public CGAL::Exact_predicates_inexact_constructions_kernel {};
> typedef K::FT FT;
> typedef K::Point_3 Point;
> typedef FT (Function)(const Point&);
> typedef CGAL::Implicit_mesh_domain_3<Function,K> Mesh_domain;
>
> // Triangulation
> typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
> typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
>
> // Criteria
> typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
> typedef Mesh_criteria::Facet_criteria Facet_criteria;
> typedef Mesh_criteria::Cell_criteria Cell_criteria;
>
> typedef C3t3::Triangulation Tr;
> typedef  C3t3::Cell_handle Cell_handle;
> typedef  Tr::Facet Facet;
> typedef  Tr::Finite_facets_iterator Finite_facets_iterator;
> typedef Tr::Finite_vertices_iterator Finite_vertices_iterator;
> typedef Tr::Finite_cells_iterator Finite_cells_iterator;
> typedef Tr::Vertex_handle Vertex_handle;
>
>
> typedef K::Ray_3 Ray;
> typedef K::Line_3 Line;
> typedef K::Direction_3 Direction;
> typedef K::Triangle_3 Triangle;
> typedef K::Segment_3 Segment;
> typedef std::list<Triangle>::iterator Iterator;
> typedef CGAL::AABB_triangle_primitive<K,Iterator> Primitive;
> typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
> typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
> typedef Tree::Object_and_primitive_id Object_and_primitive_id;
> typedef Tree::Primitive_id Primitive_id;
>
>
> // Function
> FT sphere_function (const Point& p)
> {
>         const FT x2=p.x()*p.x();
>         const FT y2=p.y()*p.y();
>         const FT z2=p.z()*p.z();
>
>         return x2+y2+z2-1;
> }
>
> int main()
> {
>         // Domain
>         Mesh_domain domain(sphere_function, K::Sphere_3(CGAL::ORIGIN, 1.));
>
>         // Criteria
>         Facet_criteria facet_criteria(30, 0.1, 0.01); // angle, size, 
> approximation
>         Cell_criteria cell_criteria(3, 0.1); // radius-edge ratio, size
>         Mesh_criteria criteria(facet_criteria, cell_criteria);
>
>         // Mesh generation
>         C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
>
>         // Output
>         std::ofstream medit_file("outmatrix.mesh");
>         c3t3.output_to_medit(medit_file);
>
>         std::cout<< "Number of verticesDT"<<" "<<
> c3t3.triangulation().number_of_vertices() << std::endl;
>         std::cout<< "Number of facetsDT"<<" "<<
> c3t3.triangulation().number_of_finite_facets() << std::endl;
>         std::cout<< "Number of facetsC3T3"<<" "<<c3t3.number_of_facets() <<
> std::endl;
>     std::cout<< "Number of cellsC3T3"<<" "<<c3t3.number_of_cells() <<
> std::endl;
>         std::cout<< "Number of cellsDT"<<"
> "<<c3t3.triangulation().number_of_finite_cells() << std::endl;
>
>
>         // for each facet build a triangle
>         std::list<Triangle> triangles;
>     for(Finite_facets_iterator fi =
> c3t3.triangulation().finite_facets_begin();
>                 fi != c3t3.triangulation().finite_facets_end(); fi++) {
>
>                 Cell_handle cell    = fi->first;
>                 int opposite_vertex_index = fi->second;
>                 if(c3t3.is_in_complex(cell)) {
>                         Triangle 
> t=c3t3.triangulation().triangle(cell,opposite_vertex_index);
>                         triangles.push_back(t);
>                 }
>                 else if (c3t3.is_in_complex(*fi))
>                 {
>                         Triangle t=c3t3.triangulation().triangle(*fi);
>                         triangles.push_back(t);
>                 }
>
>     }
>
>
>     Tree tree(triangles.begin(),triangles.end());
>         Point a(0.0, 0.0, 0.0);
>     Point b(0.7, 0.34, 0.61);
>
>         // set ray and count intersections
>        Ray ray_query(a,b);
>         Direction dir(ray_query);
>
>         std::map<Vertex_handle, int> V;
>         int inum = 1;
>         for( Finite_vertices_iterator vit =
> c3t3.triangulation().finite_vertices_begin();
>                 vit != c3t3.triangulation().finite_vertices_end();
>                 ++vit)
>         {
>                 V[vit] = inum++;
>                 Point c = vit->point();
>                 Line line_vertex(c,dir);
>
>     std::cout << tree.number_of_intersected_primitives(line_vertex)
>         << " intersections(s) with line query" << std::endl;
>
>         // computes all intersections with line query
>     std::list<Object_and_primitive_id> intersections;
>     tree.all_intersections(line_vertex, std::back_inserter(intersections));
>         std::list<Object_and_primitive_id>::iterator it_intersections=
> intersections.begin();
>     std::list<Point> points;
>         std::list<Point>::iterator it_points;
>
>         for(; it_intersections != intersections.end() ; ++it_intersections)
>
>     {
>         // get intersection object
>                 Object_and_primitive_id op = *it_intersections;
>                 CGAL::Object object = op.first;
>         Point point;
>                 Segment segment;
>         if(CGAL::assign(point,object))
>                 {
>                         points.push_back(point);
>                         //std::cout << "intersection object is a point" <<" "<< 
> point.x()<<"
> "<<point.y()<<" "<<point.z()<<std::endl;
>
>                 }
>                 else if(CGAL::assign(segment,object))
>                         std::cout << "intersection object is a segment" 
> <<std::endl;
>
>     }
>
>         //set order in intersection list
>         points.sort();
>                 int ip=1;
>                 for(it_points=points.begin(); it_points != points.end() ; 
> ++it_points) {
>                         std::cout <<  "myorderedlist contains:"<<" " << 
> ip++<<"
> "<<*it_points<<std:: endl;}
>                 std::cout <<"---" << std::endl;
>         //erase multiple intersections
>         points.unique();
>
>         int it=1;
>         for(it_points=points.begin(); it_points != points.end() ; 
> ++it_points) {
>                 std::cout.precision (27);
>                 std::cout << "uniquelist:"<<" " << it++<<" 
> "<<*it_points<<std:: endl;}
>
>         intersections.clear();
>         points.clear();
>
>         }
>
>         return 0;
> }


--
Stephane Tayeb
Software engineer - INRIA Sophia Antipolis
Geometrica Project-Team

template <class Tr>
class AABB_triangulation_triangle_primitive
{
public:
    // types
  typedef typename Tr::Finite_facet_iterator Id; // Id type
  typedef typename Tr::Geom_traits GeomTraits;
  typedef typename GeomTraits::Point_3 Point; // point type
  typedef typename GeomTraits::Triangle_3 Datum; // datum type

private:
  // member data
  Id m_it; // iterator
  Datum m_datum; // 3D triangle

public:
  AABB_triangulation_triangle_primitive() {}
  
  AABB_triangulation_triangle_primitive(Id it)
    : m_it(it)
    , m_datum(it->first->vertex((it->second+1)&3)->point(),
              it->first->vertex((it->second+2)&3)->point(),
              it->first->vertex((it->second+3)&3)->point()){}
  
public:
  Id& id() { return m_it; }
  const Id& id() const { return m_it; }
  Datum& datum() { return m_datum; }
  const Datum& datum() const { return m_datum; }

  /// Returns a point on the primitive
  Point reference_point() const { return m_datum.vertex(0); }
};