Skip to Content.
Sympa Menu

cgal-discuss - Re: [cgal-discuss] Sort primitives in AABB tree

Subject: CGAL users discussion list

List archive

Re: [cgal-discuss] Sort primitives in AABB tree


Chronological Thread 
  • From: Stephane Tayeb <>
  • To:
  • Subject: Re: [cgal-discuss] Sort primitives in AABB tree
  • Date: Tue, 03 Aug 2010 18:27:38 +0200

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

PNG image


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); }
};



Archive powered by MHonArc 2.6.16.

Top of Page