CGAL users discussion list

["Sebastien Loriot (GeometryFactory)" <>]

Add this in your for loop after

"Line line_vertex(vertexpoint,dir);"


std::vector<Cell_handle> incident_cells;
bool is_on_chull=false;
int nb_inter=0;

c3t3.triangulation().incident_cells(vit,std::back_inserter(incident_cells));
for (std::vector<Cell_handle>::iterator cit=incident_cells.begin();
                                        cit!=incident_cells.end();
                                        ++cit)
{
    if (c3t3.triangulation().is_infinite(*cit)){
      is_on_chull=true;
      continue;
    }
    int index=(*cit)->index(vit);
    K::Triangle_3 triangle((*cit)->vertex((index+1)%4)->point(),
                           (*cit)->vertex((index+2)%4)->point(),
                           (*cit)->vertex((index+3)%4)->point());
    if( CGAL::do_intersect(line_vertex,triangle) ) ++nb_inter;
  }
  CGAL_assertion( (is_on_chull && nb_inter==1) ||
                  (!is_on_chull && nb_inter==2) );

You should not have any assertion failure.

In your code, I would say that the intersection does not show up
because the intersection is computed using a kernel with inexact
constructions.

Sebastien.



 wrote:
> Hi,
> I cannot reproduce the behaviour you suggested:
>
> <Suppose that the line does not intersect the opposite facets on an edge.
> <Given a vertex, they are two cases:
>
> <1) if the vertex v is not incident to the infinite vertex:
> <the set of finite cells incident to v is a topological ball with
> <the point of v in the interior the ball. So in that case you must
> <intersect exactly two facets opposite to v.
>
> <2) if the vertex v is incident to the infinite vertex:
> <then the set of finite cells incident to v is a topological ball with
> <the point of v on the boundary of that ball. Due to the convexity of
> <any small neighborhood of v, the number of intersected facets opposite
> <to v ranges from 0 to 1 (it depends on the direction).
>
> <S.
>
> I have some single intersections even though the vertex is not incident to
> the infinite vertex. (The result doesn't change using a Facet_iterator
> instead of a Finite_facet_iterator for Id and testing
> c3t3.is_in_complex(cell) or considering only cells incident to the
> vertex). The case 2 works fine.
> Thank you for your attention.
>
> #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/AABB_tree.h>
> #include <CGAL/AABB_traits.h>
> #include <list>
> #include <CGAL/iterator.h>
> #include <iostream>
> #include <fstream>
> #include <map>
> #include <vector>
> #include <string>
> #include <utility>
> #include <cassert>
>
>         template <class Tr>
>         class AABB_triangulation_triangle_primitive
>         {
>         public:
>                 // types
>                 typedef typename Tr::Finite_facets_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); }
>         };
>
> // 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/domain
> 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;
>
> // Triangulation
> typedef  C3t3::Triangulation Tr;
> typedef  C3t3::Cell_handle Cell_handle;
> typedef  Tr::Finite_vertices_iterator Finite_vertices_iterator;
> typedef  Tr::Vertex_handle Vertex_handle;
>
> //Geometric objects
> typedef K::Ray_3 Ray;
> typedef K::Line_3 Line;
> typedef K::Direction_3 Direction;
> typedef K::Segment_3 Segment;
>
> //AABB Tree
> typedef AABB_triangulation_triangle_primitive<Tr> Primitive;
> typedef CGAL::AABB_traits<K,Primitive> Traits;
> typedef CGAL::AABB_tree<Traits> Tree;
> typedef Tree::Object_and_primitive_id Object_and_primitive_id;
> typedef Tree::Primitive_id Id;
> //typedef Tr::Finite_facets_iterator Id;
> typedef K::Triangle_3 Datum;
>
> // To avoid verbose function and named parameters call
> using namespace CGAL::parameters;
>
> // Function
> FT sphere_function (const Point& p)
> { return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-1; }
>
> int main()
> {
> // Domain
> Mesh_domain domain(sphere_function, K::Sphere_3(CGAL::ORIGIN, 2.));
>
> // Criteria
> Mesh_criteria criteria(facet_angle=30, facet_size=0.3,
> facet_distance=0.15,
>            cell_radius_edge=2, cell_size=0.3);
>
>
> // Mesh generation
> //C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
> // Mesh generation and optimization
>         C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
>                                                                               
>   lloyd(time_limit=30),
>                                                                               
>   no_perturb(),
>                                                                               
>   exude(time_limit=30, sliver_bound=10));
> // Output
> std::ofstream medit_file("out.mesh");
> c3t3.output_to_medit(medit_file);
>         Tree
> tree(c3t3.triangulation().finite_facets_begin(),c3t3.triangulation().finite_facets_end());
>         Point a(0.3, 0.8, 0.1);
>         Point b(0.4, 0.5, 0.7);
>         Ray ray_query(a,b);
>         Direction dir(ray_query);
>         for(Finite_vertices_iterator vit
> =c3t3.triangulation().finite_vertices_begin();
>                 vit != c3t3.triangulation().finite_vertices_end();
>                 ++vit)
>         {
>                 Point vertexpoint = vit->point();
>                 Line line_vertex(vertexpoint,dir);
>                 std::cout << tree.number_of_intersected_primitives(line_vertex) 
> <<
> "INTERSECTED PRIMITIVES"<<std::endl;
>                 // compute and store 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<std::pair<Point,Id> > myvector;
>                 std::list<std::pair<Point,Id> >::iterator it;
>                 for(; it_intersections != intersections.end() ; 
> it_intersections++)  {
>                         // get intersection object
>                         Object_and_primitive_id op = *it_intersections;
>                         CGAL::Object object = op.first;
>                         Id id=op.second;
>                         Point point;
>                         Segment segment;
>                         if(CGAL::assign(point,object))
>                         {
>                                 myvector.push_back(make_pair(point,id));
>                         }
>                         else if(CGAL::assign(segment,object))
>                         {
>                                 std::cout << "intersection object is a segment" 
> <<std::endl;
>                         }
>                 }
>                 //convex hull
>                 std::list<Vertex_handle>  verticeschull;
>                 std::list<Vertex_handle>::iterator it_inf;
>                 
> c3t3.triangulation().incident_vertices(c3t3.triangulation().infinite_vertex(),
> std::back_inserter(verticeschull));
>                 bool infinitevert = false;
>                 for ( it_inf=verticeschull.begin() ; it_inf != 
> verticeschull.end();
> it_inf++ ){
>                         if (vertexpoint == (*it_inf)->point())
>                                 infinitevert= true;
>                 }
>                 int ipr=1;
>           for ( it=myvector.begin() ; it != myvector.end(); it++ ){
>                   Cell_handle cell = (*((*it).second)).first;
>                     int opposite_vertex_index =(*it).second->second;
>                     Point
> vertexoppositepoint=(cell->vertex(opposite_vertex_index)->point());
>                                                 if((vertexpoint == 
> vertexoppositepoint)&&(infinitevert == true)) {
>                                 std::cout<<" "<<"squared distance"<<"
> "<<CGAL::squared_distance(vertexpoint,Point(CGAL::ORIGIN))<<"
> "<<"current vertex"<<" "<< vertexpoint<<" "<<"on BOUNDARY"<<std::endl;
>                                 //if(c3t3.is_in_complex(cell)==true) {
>                                         std::cout << ipr++<<"th 
> intersection/coordinates"<<"
> "<<std::setprecision(15)<<(*it).first << " => " << std::endl;
>                                                                         for(int 
> i = 0; i < 4; i++){
>                                                 const  Vertex_handle& vh 
> =cell->vertex(i);
>                                                 std::cout <<"tetrahedron vertex 
> coordinates:"<<"
> "<<vh->point().x()<<" "<<vh->point().y()<<" "<<vh->point().z()<<"
> "<<"squared distance"<<"
> "<<CGAL::squared_distance((vh->point()),Point(CGAL::ORIGIN))<<std::endl;
>                                         }
>                                         //}
>                         }
>                         else if((vertexpoint == 
> vertexoppositepoint)&&(infinitevert == false))
> {
>                                 std::cout<<"current vertex"<< " "<< vertexpoint<<" 
> "<<" inside the
> VOLUME"<<" "<<std::endl;
>                                 std::cout<<"squared distance"<<"
> "<<CGAL::squared_distance(vertexpoint,Point(CGAL::ORIGIN))<<std::endl;
>                                 //if(c3t3.is_in_complex(cell)==true) {
>                                 std::cout << ipr++<<"th 
> intersection/coordinates"<<"
> "<<std::setprecision(15)<<(*it).first << " => " << std::endl;
>                                 for(int i = 0; i < 4; i++){
>                                         const  Vertex_handle& vh 
> =cell->vertex(i);
>                                         std::cout <<"tetrahedron vertex 
> coordinates:"<<"
> "<<vh->point().x()<<" "<<vh->point().y()<<" "<<vh->point().z()<<"
> "<<"squared distance"<<"
> "<<CGAL::squared_distance((vh->point()),Point(CGAL::ORIGIN))<<std::endl;
>                                 }
>                                 //      }
>                                                 }
>                 }
>         }
> return 0;
> }