* Update (triangle, seg/line/ray) intersection computation

* Enrich (triangle, seg/line/ray) test suite
* Minor formatting change in Bbox_3_line_3_do_intersect.h

AABB_tree/include/CGAL/AABB_intersections/Bbox_3_line_3_do_intersect.h

@@ -121,13 +121,14 @@ namespace internal {
   {
     typedef typename K::FT FT;
     typedef typename K::Point_3 Point_3;
+    typedef typename K::Vector_3 Vector_3;
     
-    const Point_3 point_a = line.point(0);
+    const Point_3 point = line.point();
-    const Point_3 point_b = line.point(1);
+    const Vector_3 v = line.to_vector();
     
     return bbox_line_do_intersect_aux(
-                         point_a.x(), point_a.y(), point_a.z(), 
+                         point.x(), point.y(), point.z(), 
-                         point_b.x(), point_b.y(), point_b.z(),
+                         point.x()+v.x(), point.y()+v.y(), point.z()+v.z(),
                          FT(bbox.xmin()), FT(bbox.ymin()), FT(bbox.zmin()),
                          FT(bbox.xmax()), FT(bbox.ymax()), FT(bbox.zmax()) );
   }

AABB_tree/include/CGAL/AABB_intersections/Triangle_3_line_3_intersection.h

@@ -31,6 +31,47 @@
 namespace CGAL {
 namespace internal {
 
+template <class K>
+typename K::Point_3
+t3l3_intersection_coplanar_aux(const typename K::Line_3& l,
+                               const typename K::Point_3& a,
+                               const typename K::Point_3& b,
+                               const K& k)
+{  
+  // Returns the intersection point between line l and segment [a,b]
+  //
+  // preconditions:
+  //   + l,a,b are coplanar
+  
+  typedef typename K::Point_3 Point_3;
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::FT FT;
+  
+  typename K::Construct_vector_3 vector = 
+    k.construct_vector_3_object();
+  
+  typename K::Construct_cross_product_vector_3 cross_product = 
+    k.construct_cross_product_vector_3_object();
+  
+  typename K::Compute_scalar_product_3 scalar_product = 
+    k.compute_scalar_product_3_object();
+  
+  typename K::Compute_squared_length_3 sq_length = 
+    k.compute_squared_length_3_object();
+  
+  const Point_3 p = l.point();
+  const Vector_3 v = l.to_vector();
+  const Vector_3 ab = vector(a,b);
+  const Vector_3 pa = vector(p,a);
+  
+  const Vector_3 pa_ab = cross_product(pa,ab);
+  const Vector_3 v_ab = cross_product(v,ab);
+  
+  const FT t = scalar_product(pa_ab,v_ab) / sq_length(v_ab);
+  
+  return ( p + t*v );
+}
+  
   
 template <class K>
 Object
@@ -44,12 +85,12 @@ t3l3_intersection_coplanar_aux(const typename K::Point_3& a,
   // This function is designed to clip pq into the triangle abc.
   // Point configuration should be as follows
   //
-  //     +q
-  //     |    +a
   //     |
-  //  +c |       +b
+  //     |    +b
   //     |
-  //     +p
+  //  +c |       +a
+  //     |
+  //     | l
   //
   // We know that c is isolated on the negative side of pq
   
@@ -65,26 +106,13 @@ t3l3_intersection_coplanar_aux(const typename K::Point_3& a,
     k.construct_segment_3_object();
   
   // Let's get the intersection points
-  Object l_bc_obj = intersection(l,line(b,c));
+  const Point_3 l_bc = t3l3_intersection_coplanar_aux(l,b,c,k);
-  const Point_3* l_bc = object_cast<Point_3>(&l_bc_obj);
+  const Point_3 l_ca = t3l3_intersection_coplanar_aux(l,c,a,k);
-  if ( NULL == l_bc )
-  {
-    CGAL_kernel_assertion(false);
-    return Object();
-  }
-  
-  Object l_ca_obj = intersection(l,line(c,a));
-  const Point_3* l_ca = object_cast<Point_3>(&l_ca_obj);
-  if ( NULL == l_ca )
-  {
-    CGAL_kernel_assertion(false);
-    return Object();
-  }
   
   if ( negative_side )
-    return make_object(segment(*l_bc, *l_ca));
+    return make_object(segment(l_bc, l_ca));
   else
-    return make_object(segment(*l_ca, *l_bc));  
+    return make_object(segment(l_ca, l_bc));  
 }
   
   
@@ -140,13 +168,11 @@ intersection_coplanar(const typename K::Triangle_3 &t,
   const Point_3& b = vertex_on(t,k1);
   const Point_3& c = vertex_on(t,k2);
   
-  // Test whether the segment's supporting line intersects the
+  // Test whether the line intersects the triangle in the common plane
-  // triangle in the common plane
   const Orientation pqa = coplanar_orientation(p,q,a);
   const Orientation pqb = coplanar_orientation(p,q,b);
   const Orientation pqc = coplanar_orientation(p,q,c);
   
-  
   switch ( pqa ) {
     // -----------------------------------
     // pqa POSITIVE
@@ -298,8 +324,24 @@ intersection_coplanar(const typename K::Triangle_3 &t,
 }
   
   
+template <class K>
+inline
+Object 
+t3l3_intersection_aux(const typename K::Triangle_3 &t,
+                      const typename K::Line_3 &l,
+                      const K& k)
+{
+  typename K::Intersect_3 intersection =
+    k.intersect_3_object();
   
+  Object obj = intersection(l,t.supporting_plane());
   
+  // Intersection should be a point (because of orientation test done before)
+  if ( NULL == object_cast<typename K::Line_3>(&obj) )
+    return obj;
+  else
+    return Object();
+}
   
   
   
@@ -335,22 +377,23 @@ intersection(const typename K::Triangle_3 &t,
   const Point_3 & p = point_on(l,0);
   const Point_3 & q = point_on(l,1);
   
+
+  
   if ( ( orientation(a,b,c,p) != COPLANAR )
       || ( orientation(a,b,c,q) != COPLANAR ) )
   {
     const Orientation pqab = orientation(p,q,a,b);
     const Orientation pqbc = orientation(p,q,b,c);
-    
     switch ( pqab ) {
       case POSITIVE: 
         if ( pqbc != NEGATIVE && orientation(p,q,c,a) != NEGATIVE )
-          return intersection(l,t.supporting_plane());
+          return t3l3_intersection_aux(t,l,k);
         else
           return Object();
         
       case NEGATIVE:
         if ( pqbc != POSITIVE && orientation(p,q,c,a) != POSITIVE )
-          return intersection(l,t.supporting_plane());
+          return t3l3_intersection_aux(t,l,k);
         else
           return Object();
         
@@ -358,18 +401,18 @@ intersection(const typename K::Triangle_3 &t,
         switch ( pqbc ) {
           case POSITIVE: 
             if ( orientation(p,q,c,a) != NEGATIVE )
-              return intersection(l,t.supporting_plane());
+              return t3l3_intersection_aux(t,l,k);
             else
               return Object();
             
           case NEGATIVE:
             if ( orientation(p,q,c,a) != POSITIVE )
-              return intersection(l,t.supporting_plane());
+              return t3l3_intersection_aux(t,l,k);
             else
               return Object();
               
           case COPLANAR: // pqa or pqb or pqc are collinear
-            return intersection(l,t.supporting_plane());
+            return t3l3_intersection_aux(t,l,k);
           
           default: // should not happen.
             CGAL_kernel_assertion(false);

AABB_tree/include/CGAL/AABB_intersections/Triangle_3_ray_3_intersection.h

@@ -26,39 +26,434 @@
 namespace CGAL {
 namespace internal {
 
-
+  
-
+template <class K>
+typename K::Point_3
+t3r3_intersection_coplanar_aux(const typename K::Point_3& p,
+                               const typename K::Vector_3& v,
+                               const typename K::Point_3& a,
+                               const typename K::Point_3& b,
+                               const K& k)
+{  
+  // Returns the intersection point between line (p,v) and line (a,b)
+  //
+  // preconditions:
+  //   + p,v,a,b are coplanar
+  
+  typedef typename K::Point_3 Point_3;
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::FT FT;
+  
+  typename K::Construct_vector_3 vector = 
+  k.construct_vector_3_object();
+  
+  typename K::Construct_cross_product_vector_3 cross_product = 
+  k.construct_cross_product_vector_3_object();
+  
+  typename K::Compute_scalar_product_3 scalar_product = 
+  k.compute_scalar_product_3_object();
+  
+  typename K::Compute_squared_length_3 sq_length = 
+  k.compute_squared_length_3_object();
+  
+  const Vector_3 ab = vector(a,b);
+  const Vector_3 pa = vector(p,a);
+  
+  const Vector_3 pa_ab = cross_product(pa,ab);
+  const Vector_3 v_ab = cross_product(v,ab);
+  
+  const FT t = scalar_product(pa_ab,v_ab) / sq_length(v_ab);
+  
+  return ( p + t*v );
+}
+  
+  
 template <class K>
 Object
-intersection_triangle_ray_aux(const typename K::Ray_3 &r,
+t3r3_intersection_coplanar_aux(const typename K::Point_3& a,
-                              const typename K::Segment_3 &s,
+                               const typename K::Point_3& b,
-                              const K& k)
+                               const typename K::Point_3& c,
+                               const typename K::Ray_3& r,
+                               const bool negative_side,
+                               const K& k)
 {
+  // This function is designed to clip r into the triangle abc.
+  // Point configuration should be as follows
+  //
+  //     
+  //    p+    +b
+  //     |
+  //  +c |       +a
+  //     | 
+  //     |r
+  //
+  // We know that c is isolated on the negative side of r
+  // but we don't know p position wrt [bc]&[ca]
+  
   typedef typename K::Point_3 Point_3;
+  typedef typename K::Vector_3 Vector_3;
+  
+  typename K::Coplanar_orientation_3 coplanar_orientation = 
+    k.coplanar_orientation_3_object();
+  
+  typename K::Construct_segment_3 segment =
+    k.construct_segment_3_object();
 
-  typename K::Has_on_3 has_on = k.has_on_3_object();
+  typename K::Construct_point_on_3 point_on = 
-  typename K::Construct_object_3 make_object = k.construct_object_3_object();
+    k.construct_point_on_3_object();
-  typename K::Construct_direction_3 direction = k.construct_direction_3_object();
+  
-  typename K::Construct_segment_3 make_segment = k.construct_segment_3_object();
+  const Point_3& p = point_on(r,0);
-  // typename K::Compare_xyz_3 compare = k.compare_xyz_3_object(); // PA: never used
+  
-  typename K::Equal_3 equal = k.equal_3_object();
+  // A ray is not symetric, 2 cases depending on isolated side of c
-
+  Orientation cap;
-  const Point_3& p = r.source();
+  if ( negative_side )
-  const Point_3& q1 = s.source();
+    cap = coplanar_orientation(c,a,p);
-  const Point_3& q2 = s.target();
-
-  if ( ! has_on(s,p) )
-    return make_object(s);
-
-  if ( equal(p,q1) || equal(p,q2) )
-    return make_object(p);
-
-  if ( direction(r) == direction(s) )
-    return make_object(make_segment(p,q2));
   else
-    return make_object(make_segment(p,q1));
+    cap = coplanar_orientation(b,c,p);
+  
+  switch ( cap ) {
+  
+    case NEGATIVE:
+      // p is bellow [c,a]
+      return Object();
+    
+    case COLLINEAR:
+      return make_object(p);
+  
+    case POSITIVE:
+    {
+      // Compute the intersection points between ray and [b,c],[c,a]
+      Vector_3 v = r.to_vector();
+      
+      // Get intersection point at p side
+      Point_3 p_side_end_point(p);
+      Point_3 q_side_end_point;
+      
+      // A ray is not symetric, 2 cases depending on isolated side of c
+      if ( negative_side )
+      {
+        if ( NEGATIVE == coplanar_orientation(b,c,p) )
+          p_side_end_point = t3r3_intersection_coplanar_aux(p,v,b,c,k);
+        
+        // Get other end point (always intersection computation on the unbounded
+        // side of the ray)
+        q_side_end_point = t3r3_intersection_coplanar_aux(p,v,c,a,k);
+      }
+      else
+      {
+        if ( NEGATIVE == coplanar_orientation(c,a,p) )
+          p_side_end_point = t3r3_intersection_coplanar_aux(p,v,c,a,k);
+        
+        // Get other end point (always intersection computation on the unbounded
+        // side of the ray)
+        q_side_end_point = t3r3_intersection_coplanar_aux(p,v,b,c,k);
+      }
+      
+      // Build result
+      return make_object(segment(p_side_end_point, q_side_end_point));
+    }
+    
+    default: // should not happen.
+      CGAL_kernel_assertion(false);
+      return Object();
+  }
+  
+  CGAL_kernel_assertion(false);
+  return Object();
+}
+  
+  
+  
+  
+  
+  
+  
+template <class K>
+Object
+intersection_coplanar(const typename K::Triangle_3 &t, 
+                      const typename K::Ray_3  &r,
+                      const K & k )
+{
+  CGAL_kernel_precondition( ! k.is_degenerate_3_object()(t) ) ;
+  CGAL_kernel_precondition( ! k.is_degenerate_3_object()(r) ) ;
+  
+  typedef typename K::Point_3 Point_3;
+  
+  
+  typename K::Construct_point_on_3 point_on = 
+    k.construct_point_on_3_object();
+  
+  typename K::Construct_vertex_3 vertex_on =
+    k.construct_vertex_3_object();
+  
+  typename K::Coplanar_orientation_3 coplanar_orientation = 
+    k.coplanar_orientation_3_object();
+  
+  typename K::Construct_line_3 line =
+    k.construct_line_3_object();
+  
+  typename K::Construct_segment_3 segment =
+    k.construct_segment_3_object();
+  
+  typename K::Collinear_are_ordered_along_line_3 collinear_ordered =
+    k.collinear_are_ordered_along_line_3_object();
+  
+  
+  const Point_3 & p = point_on(r,0);
+  const Point_3 & q = point_on(r,1);
+  
+  const Point_3 & A = vertex_on(t,0);
+  const Point_3 & B = vertex_on(t,1);
+  const Point_3 & C = vertex_on(t,2);
+  
+  int k0 = 0;
+  int k1 = 1;
+  int k2 = 2;
+  
+  // Determine the orientation of the triangle in the common plane
+  if (coplanar_orientation(A,B,C) != POSITIVE)
+  {
+    // The triangle is not counterclockwise oriented
+    // swap two vertices.
+    std::swap(k1,k2);
+  }
+  
+  const Point_3& a = vertex_on(t,k0);
+  const Point_3& b = vertex_on(t,k1);
+  const Point_3& c = vertex_on(t,k2);
+  
+  // Test whether the ray's supporting line intersects the
+  // triangle in the common plane
+  const Orientation pqa = coplanar_orientation(p,q,a);
+  const Orientation pqb = coplanar_orientation(p,q,b);
+  const Orientation pqc = coplanar_orientation(p,q,c);
+  
+  switch ( pqa ) {
+    // -----------------------------------
+    // pqa POSITIVE
+    // ----------------------------------- 
+    case POSITIVE:
+      switch ( pqb ) {
+        case POSITIVE:
+          switch ( pqc ) {
+            case POSITIVE:
+              // the triangle lies in the positive halfspace
+              // defined by the segment's supporting line.
+              return Object();
+              
+            case NEGATIVE:
+              // c is isolated on the negative side
+              return t3r3_intersection_coplanar_aux(a,b,c,r,true,k);
+              
+            case COLLINEAR:
+              // p,q,c are collinear
+              if ( collinear_ordered(p,c,q) || collinear_ordered(p,q,c) )
+                return make_object(c);
+              else
+                return Object();
+          }
+          
+        case NEGATIVE:
+          if ( POSITIVE == pqc )
+            // b is isolated on the negative side
+            return t3r3_intersection_coplanar_aux(c,a,b,r,true,k);
+          else
+            // a is isolated on the positive side (here mb c could be use as
+            // an endpoint instead of computing an intersection is some cases)
+            return t3r3_intersection_coplanar_aux(b,c,a,r,false,k);
+          
+        case COLLINEAR:
+          switch ( pqc ) {
+            case POSITIVE:
+              // p,q,b are collinear
+              if ( collinear_ordered(p,b,q) || collinear_ordered(p,q,b) )
+                return make_object(b);
+              else
+                return Object();
+              
+            case NEGATIVE:
+              // a is isolated on the positive side (here mb b could be use as
+              // an endpoint instead of computing an intersection)
+              return t3r3_intersection_coplanar_aux(b,c,a,r,false,k);
+              
+            case COLLINEAR:
+              // b,c,p,q are aligned, [p,q]&[b,c] have the same direction
+              if ( collinear_ordered(p,b,c) )
+                return make_object(segment(b,c));
+              else
+                return make_object(segment(p,c));
+          }
+          
+        default: // should not happen.
+          CGAL_kernel_assertion(false);
+          return Object();
+      }
+      
+    // -----------------------------------
+    // pqa NEGATIVE
+    // ----------------------------------- 
+    case NEGATIVE:
+      switch ( pqb ) {
+        case POSITIVE:
+          if ( POSITIVE == pqc )
+            // a is isolated on the negative side
+            return t3r3_intersection_coplanar_aux(b,c,a,r,true,k);
+          else
+            // b is isolated on the positive side (here mb c could be use as
+            // an endpoint instead of computing an intersection, in some cases)
+            return t3r3_intersection_coplanar_aux(c,a,b,r,false,k);
+          
+        case NEGATIVE:
+          switch ( pqc ) {
+            case POSITIVE:
+              // c is isolated on the positive side
+              return t3r3_intersection_coplanar_aux(a,b,c,r,false,k);
+              
+            case NEGATIVE:
+              // the triangle lies in the negative halfspace
+              // defined by the segment's supporting line.
+              return Object();
+              
+            case COLLINEAR:
+              // p,q,c are collinear
+              if ( collinear_ordered(p,c,q) || collinear_ordered(p,q,c) )
+                return make_object(c);
+              else
+                return Object();
+          }
+          
+        case COLLINEAR:
+          switch ( pqc ) {
+            case POSITIVE:
+              // a is isolated on the negative side (here mb b could be use as
+              // an endpoint instead of computing an intersection)
+              return t3r3_intersection_coplanar_aux(b,c,a,r,true,k);       
+              
+            case NEGATIVE:
+              // p,q,b are collinear
+              if ( collinear_ordered(p,b,q) || collinear_ordered(p,q,b) )
+                return make_object(b);
+              else
+                return Object();
+              
+            case COLLINEAR:
+              // b,c,p,q are aligned, [p,q]&[c,b] have the same direction
+              if ( collinear_ordered(p,c,b) )
+                return make_object(segment(c,b));
+              else
+                return make_object(segment(p,b));
+          }
+          
+        default: // should not happen.
+          CGAL_kernel_assertion(false);
+          return Object();
+      }
+      
+    // -----------------------------------
+    // pqa COLLINEAR
+    // -----------------------------------   
+    case COLLINEAR:
+      switch ( pqb ) {
+        case POSITIVE:
+          switch ( pqc ) {
+            case POSITIVE:
+              // p,q,a are collinear
+              if ( collinear_ordered(p,a,q) || collinear_ordered(p,q,a) )
+                return make_object(a);
+              else
+                return Object();
+              
+            case NEGATIVE:
+              // b is isolated on the positive side (here mb a could be use as
+              // an endpoint instead of computing an intersection)
+              return t3r3_intersection_coplanar_aux(c,a,b,r,false,k);
+              
+            case COLLINEAR:
+              // a,c,p,q are aligned, [p,q]&[c,a] have the same direction
+              if ( collinear_ordered(p,c,a) )
+                return make_object(segment(c,a));
+              else
+                return make_object(segment(p,a));
+          }
+          
+        case NEGATIVE:
+          switch ( pqc ) {
+            case POSITIVE:
+              // b is isolated on the negative side (here mb a could be use as
+              // an endpoint instead of computing an intersection)
+              return t3r3_intersection_coplanar_aux(c,a,b,r,true,k);
+              
+            case NEGATIVE:
+              // p,q,a are collinear
+              if ( collinear_ordered(p,a,q) || collinear_ordered(p,q,a) )
+                return make_object(a);
+              else
+                return Object();
+              
+            case COLLINEAR:
+              // a,c,p,q are aligned, [p,q]&[a,c] have the same direction
+              if ( collinear_ordered(p,a,c) )
+                return make_object(segment(a,c));
+              else
+                return make_object(segment(p,c));
+          }
+          
+        case COLLINEAR:
+          switch ( pqc ) {
+            case POSITIVE:
+              // a,b,p,q are aligned, [p,q]&[a,b] have the same direction
+              if ( collinear_ordered(p,a,b) )
+                return make_object(segment(a,b));
+              else
+                return make_object(segment(p,b));
+              
+            case NEGATIVE:
+              // a,b,p,q are aligned, [p,q]&[b,a] have the same direction
+              if ( collinear_ordered(p,b,a) )
+                return make_object(segment(b,a));
+              else
+                return make_object(segment(p,a));
+              
+            case COLLINEAR:
+              // case pqc == COLLINEAR is impossible since the triangle is
+              // assumed to be non flat
+              CGAL_kernel_assertion(false);
+              return Object();
+          }
+          
+        default: // should not happen.
+          CGAL_kernel_assertion(false);
+          return Object();
+          
+      }
+      
+    default:// should not happen.
+      CGAL_kernel_assertion(false);
+      return Object();
+      
+  }
+}
+  
+  
+  
+template <class K>
+inline
+Object 
+t3r3_intersection_aux(const typename K::Triangle_3 &t,
+                      const typename K::Ray_3 &r,
+                      const K& k)
+{
+  typename K::Intersect_3 intersection =
+    k.intersect_3_object();
+  
+  Object obj = intersection(r.supporting_line(),t.supporting_plane());
+  
+  // Intersection should be a point (because of orientation test done before)
+  if ( NULL == object_cast<typename K::Line_3>(&obj) )
+    return obj;
+  else
+    return Object();
 }
-
 
 template <class K>
 Object
@@ -66,81 +461,126 @@ intersection(const typename K::Triangle_3  &t,
              const typename K::Ray_3 &r,
              const K& k)
 {
+  CGAL_kernel_precondition( ! k.is_degenerate_3_object()(t) ) ;
+  CGAL_kernel_precondition( ! k.is_degenerate_3_object()(r) ) ;
+  
   typedef typename K::Point_3 Point_3;
-  typedef typename K::Segment_3 Segment_3;
+  
-  typedef typename K::Object_3 Object_3;
+  typename K::Construct_vertex_3 vertex_on =
-  typedef typename K::Ray_3 Ray_3;
+    k.construct_vertex_3_object();
-  typedef typename K::Line_3 Line_3;
+  
+  typename K::Orientation_3 orientation = 
+    k.orientation_3_object();
+  
+  typename K::Construct_ray_3 ray = 
+    k.construct_ray_3_object();
+  
+  typename K::Construct_point_on_3 point_on = 
+    k.construct_point_on_3_object();
 
+  
+  const Point_3& a = vertex_on(t,0);
+  const Point_3& b = vertex_on(t,1);
+  const Point_3& c = vertex_on(t,2);
+  const Point_3& p = point_on(r,0);
+  const Point_3& q = point_on(r,1);
 
-  if ( !CGAL::do_intersect(t, r) )
+  Point_3 d = point_on(ray(a,r.to_vector()),1);
-    return Object_3();
+  
+  const Orientation ray_direction = orientation(a,b,c,d);
+  const Orientation abcp = orientation(a,b,c,p);
 
-	// TOFIX: here we assume that we have already tested
+  switch ( abcp ) {
-	// do_intersect between the triangle and the ray
+    case POSITIVE: 
-  const Object intersection =  CGAL::intersection(t.supporting_plane(),
+      switch ( ray_direction ) {
-                                                  r);
+        case POSITIVE: 
-
+          // the ray lies in the positive open halfspaces defined by the
-  // intersection is either a point, either a ray
+          // triangle's supporting plane
-  // if it is a ray, then we need to clip it to a segment
+          return Object();
-  if ( const Ray_3* ray = object_cast<Ray_3>(&intersection))
+          
-  {
+        case NEGATIVE:
-    typename K::Intersect_3 intersect = k.intersect_3_object();
+          // The ray straddles the triangle's plane
-    typename K::Construct_line_3 f_line = k.construct_line_3_object();
+          // p sees the triangle in counterclockwise order
-    typename K::Construct_vertex_3 vertex_on = k.construct_vertex_3_object();
+          if ( orientation(p,q,a,b) != POSITIVE
-    typename K::Construct_object_3 make_object = k.construct_object_3_object();
+              && orientation(p,q,b,c) != POSITIVE
-    typename K::Construct_segment_3 f_segment = k.construct_segment_3_object();
+              && orientation(p,q,c,a) != POSITIVE )
-    typename K::Has_on_3 has_on = k.has_on_3_object();
+            return t3r3_intersection_aux(t,r,k);
-
+          else
-    const Point_3& t0 = vertex_on(t,0);
+            return Object();
-    const Point_3& t1 = vertex_on(t,1);
+          
-    const Point_3& t2 = vertex_on(t,2);
+        case COPLANAR:
-
+          // The ray lie in a plane parallel to a,b,c support plane
-    const Line_3 l01 = f_line(t0,t1);
+          return Object();
-    const Line_3 l02 = f_line(t0,t2);
+          
-    const Line_3 l12 = f_line(t1,t2);
+        default: // should not happen.
-
+          CGAL_kernel_assertion(false);
-    const Segment_3 s01 = f_segment(t0,t1);
+          return Object();
-    const Segment_3 s02 = f_segment(t0,t2);
+      }
-    const Segment_3 s12 = f_segment(t1,t2);
+      
-
+    case NEGATIVE:
-    const Line_3 lr = ray->supporting_line();
+      switch ( ray_direction ) {
-
+        case POSITIVE: 
-    const Object_3 inter_01 = intersect(l01,lr);
+          // The ray straddles the triangle's plane
-    const Object_3 inter_02 = intersect(l02,lr);
+          // q sees the triangle in counterclockwise order
-    const Object_3 inter_12 = intersect(l12,lr);
+          
-
+          if ( orientation(q,p,a,b) != POSITIVE
-    const Point_3* p01 = object_cast<Point_3>(&inter_01);
+              && orientation(q,p,b,c) != POSITIVE
-    const Point_3* p02 = object_cast<Point_3>(&inter_02);
+              && orientation(q,p,c,a) != POSITIVE )
-    const Point_3* p12 = object_cast<Point_3>(&inter_12);
+            return t3r3_intersection_aux(t,r,k);
-
+          else
-    if ( p01 && has_on(s01, *p01) )
+            return Object();
-    {
+          
-      if ( p02 && has_on(s02, *p02) )
+        case NEGATIVE:
-        return intersection_triangle_ray_aux(*ray, f_segment(*p01,*p02), k);
+          // the ray lies in the negative open halfspaces defined by the
-      else if ( p12 && has_on(s12, *p12) )
+          // triangle's supporting plane
-        return intersection_triangle_ray_aux(*ray, f_segment(*p01,*p12), k);
+          return Object();
-      else
+          
-        return make_object(*p01);
+        case COPLANAR:
-    }
+          // The ray lie in a plane parallel to a,b,c support plane
-    else if ( p02 && has_on(s02, *p02) )
+          return Object();
-    {
+          
-      if  ( p12 && has_on(s12, *p12) )
+        default: // should not happen.
-        return intersection_triangle_ray_aux(*ray, f_segment(*p02,*p12), k);
+          CGAL_kernel_assertion(false);
-      else
+          return Object();
-        return make_object(*p02);
+      }
-    }
+      
-    else if ( p12 && has_on(s12, *p12) )
+    case COPLANAR: // p belongs to the triangle's supporting plane
-    {
+      switch ( ray_direction ) {
-      return make_object(*p12);
+        case POSITIVE: 
-    }
+          // q sees the triangle in counterclockwise order
-
+          if ( orientation(q,p,a,b) != POSITIVE
-    // should not happen
+              && orientation(q,p,b,c) != POSITIVE
-    CGAL_kernel_assertion(false);
+              && orientation(q,p,c,a) != POSITIVE )
-    return Object_3();
+            return t3r3_intersection_aux(t,r,k);
-  }
+          else
-  else
+            return Object();
-    return intersection;
+          
+        case NEGATIVE:
+          // q sees the triangle in clockwise order
+          if ( orientation(p,q,a,b) != POSITIVE
+              && orientation(p,q,b,c) != POSITIVE
+              && orientation(p,q,c,a) != POSITIVE )
+            return t3r3_intersection_aux(t,r,k);
+          else
+            return Object();
+          
+        case COPLANAR: 
+          // The ray lie in triangle supporting plane
+          return intersection_coplanar(t,r,k);
+          
+        default: // should not happen.
+          CGAL_kernel_assertion(false);
+          return Object();
+      }
+      
+    default: // should not happen.
+      CGAL_kernel_assertion(false);
+      return Object();
+      
+  } 
+  
+  CGAL_kernel_assertion(false);
+  return Object();
 }
 
 } // end namespace internal

AABB_tree/include/CGAL/AABB_intersections/Triangle_3_segment_3_intersection.h

@@ -33,6 +33,47 @@
 namespace CGAL {
 namespace internal {
 
+template <class K>
+typename K::Point_3
+t3s3_intersection_coplanar_aux(const typename K::Point_3& p,
+                               const typename K::Point_3& q,
+                               const typename K::Point_3& a,
+                               const typename K::Point_3& b,
+                               const K& k)
+{  
+  // Returns the intersection point between segment [p,q] and [a,b]
+  //
+  // preconditions:
+  //   + p,q,a,b are coplanar
+  
+  typedef typename K::Point_3 Point_3;
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::FT FT;
+  
+  typename K::Construct_vector_3 vector = 
+    k.construct_vector_3_object();
+  
+  typename K::Construct_cross_product_vector_3 cross_product = 
+    k.construct_cross_product_vector_3_object();
+  
+  typename K::Compute_scalar_product_3 scalar_product = 
+    k.compute_scalar_product_3_object();
+  
+  typename K::Compute_squared_length_3 sq_length = 
+    k.compute_squared_length_3_object();
+  
+  const Vector_3 pq = vector(p,q);
+  const Vector_3 ab = vector(a,b);
+  const Vector_3 pa = vector(p,a);
+  
+  const Vector_3 pa_ab = cross_product(pa,ab);
+  const Vector_3 pq_ab = cross_product(pq,ab);
+  
+  const FT t = scalar_product(pa_ab,pq_ab) / sq_length(pq_ab);
+  
+  return ( p + t*pq );
+}
+  
   
 template <class K>
 Object
@@ -47,12 +88,12 @@ t3s3_intersection_coplanar_aux(const typename K::Point_3& a,
   // This function is designed to clip pq into the triangle abc.
   // Point configuration should be as follows
   //
-  //     +q
-  //     |    +a
-  //     |
-  //  +c |       +b
-  //     |
   //     +p
+  //     |    +b
+  //     |
+  //  +c |       +a
+  //     |
+  //     +q
   //
   // We know that c is isolated on the negative side of pq, but we don't know
   // p position wrt [bc]&[ca] and q position wrt [bc]&[ca]
@@ -62,12 +103,6 @@ t3s3_intersection_coplanar_aux(const typename K::Point_3& a,
   typename K::Coplanar_orientation_3 coplanar_orientation = 
     k.coplanar_orientation_3_object();
   
-  typename K::Intersect_3 intersection =
-    k.intersect_3_object();
-  
-  typename K::Construct_line_3 line =
-    k.construct_line_3_object();
-  
   typename K::Construct_segment_3 segment =
     k.construct_segment_3_object();
   
@@ -89,19 +124,13 @@ t3s3_intersection_coplanar_aux(const typename K::Point_3& a,
     Point_3 p_side_end_point(p);
     if ( NEGATIVE == coplanar_orientation(b,c,p) )
     {
-      Object obj = intersection(line(p,q),line(b,c));
+      p_side_end_point = t3s3_intersection_coplanar_aux(p,q,b,c,k);
-      const Point_3* p_temp = object_cast<Point_3>(&obj);
-      if ( NULL != p_temp )
-        p_side_end_point = *p_temp;
     }
     
     Point_3 q_side_end_point(q);
     if ( NEGATIVE == coplanar_orientation(c,a,q) )
     {
-      Object obj = intersection(line(p,q),line(c,a));
+      q_side_end_point = t3s3_intersection_coplanar_aux(p,q,c,a,k);
-      const Point_3* p_temp = object_cast<Point_3>(&obj);
-      if ( NULL != p_temp )
-        q_side_end_point = *p_temp;
     }
 
     if ( negative_side )
@@ -114,7 +143,7 @@ t3s3_intersection_coplanar_aux(const typename K::Point_3& a,
   
 template <class K>
 Object
-t3s3_intersection_coplanar_aux(const typename K::Point_3& a,
+t3s3_intersection_collinear_aux(const typename K::Point_3& a,
                                const typename K::Point_3& b,
                                const typename K::Point_3& p,
                                const typename K::Point_3& q,
@@ -210,7 +239,6 @@ intersection_coplanar(const typename K::Triangle_3 &t,
   const Orientation pqb = coplanar_orientation(p,q,b);
   const Orientation pqc = coplanar_orientation(p,q,c);
   
-  
   switch ( pqa ) {
     // -----------------------------------
     // pqa POSITIVE
@@ -253,7 +281,7 @@ intersection_coplanar(const typename K::Triangle_3 &t,
               return t3s3_intersection_coplanar_aux(b,c,a,q,p,false,k);
             case COLLINEAR:
               // b,c,p,q are aligned, [p,q]&[b,c] have the same direction
-              return t3s3_intersection_coplanar_aux(b,c,p,q,k);
+              return t3s3_intersection_collinear_aux(b,c,p,q,k);
           }
           
         default: // should not happen.
@@ -302,7 +330,7 @@ intersection_coplanar(const typename K::Triangle_3 &t,
                 return Object();
             case COLLINEAR:
               // b,c,p,q are aligned, [p,q]&[c,b] have the same direction
-              return t3s3_intersection_coplanar_aux(c,b,p,q,k);
+              return t3s3_intersection_collinear_aux(c,b,p,q,k);
           }
           
         default: // should not happen.
@@ -327,7 +355,7 @@ intersection_coplanar(const typename K::Triangle_3 &t,
               return t3s3_intersection_coplanar_aux(c,a,b,q,p,false,k);
             case COLLINEAR:
               // a,c,p,q are aligned, [p,q]&[c,a] have the same direction
-              return t3s3_intersection_coplanar_aux(c,a,p,q,k);
+              return t3s3_intersection_collinear_aux(c,a,p,q,k);
           }
           
         case NEGATIVE:
@@ -342,17 +370,17 @@ intersection_coplanar(const typename K::Triangle_3 &t,
                 return Object();
             case COLLINEAR:
               // a,c,p,q are aligned, [p,q]&[a,c] have the same direction
-              return t3s3_intersection_coplanar_aux(a,c,p,q,k);
+              return t3s3_intersection_collinear_aux(a,c,p,q,k);
           }
           
         case COLLINEAR:
           switch ( pqc ) {
             case POSITIVE:
               // a,b,p,q are aligned, [p,q]&[a,b] have the same direction
-              return t3s3_intersection_coplanar_aux(a,b,p,q,k);
+              return t3s3_intersection_collinear_aux(a,b,p,q,k);
             case NEGATIVE:
               // a,b,p,q are aligned, [p,q]&[b,a] have the same direction
-              return t3s3_intersection_coplanar_aux(b,a,p,q,k);
+              return t3s3_intersection_collinear_aux(b,a,p,q,k);
             case COLLINEAR:
               // case pqc == COLLINEAR is impossible since the triangle is
               // assumed to be non flat
@@ -409,7 +437,6 @@ intersection(const typename K::Triangle_3 &t,
   const Orientation abcp = orientation(a,b,c,p);
   const Orientation abcq = orientation(a,b,c,q);
   
-  
   switch ( abcp ) {
     case POSITIVE:
       switch ( abcq ) {
@@ -423,7 +450,14 @@ intersection(const typename K::Triangle_3 &t,
           if ( orientation(p,q,a,b) != POSITIVE
               && orientation(p,q,b,c) != POSITIVE
               && orientation(p,q,c,a) != POSITIVE )
-            return intersection(s,t.supporting_plane());  
+          { 
+            // The intersection should be a point
+            Object obj = intersection(s.supporting_line(),t.supporting_plane());  
+            if ( NULL == object_cast<typename K::Line_3>(&obj) )
+              return obj;
+            else
+              return Object();
+          }
           else
             return Object();
           
@@ -448,7 +482,13 @@ intersection(const typename K::Triangle_3 &t,
           if ( orientation(q,p,a,b) != POSITIVE
             && orientation(q,p,b,c) != POSITIVE
             && orientation(q,p,c,a) != POSITIVE )
-            return intersection(s,t.supporting_plane());  
+          { 
+            Object obj = intersection(s.supporting_line(),t.supporting_plane());  
+            if ( NULL == object_cast<typename K::Line_3>(&obj) )
+              return obj;
+            else
+              return Object();
+          }
           else
             return Object();
           

AABB_tree/test/AABB_tree/aabb_intersection_test.cpp

@@ -31,6 +31,184 @@
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
 
+
+
+// -----------------------------------
+// Kernels
+// -----------------------------------
+struct Sc_f : public CGAL::Simple_cartesian<float> {};
+struct Sc_d : public CGAL::Simple_cartesian<double> {};
+struct C_f : public CGAL::Cartesian<float> {};
+struct C_d : public CGAL::Cartesian<double> {};
+struct Epic : public CGAL::Exact_predicates_inexact_constructions_kernel {};
+struct Epec : public CGAL::Exact_predicates_exact_constructions_kernel {};
+
+
+
+// -----------------------------------
+// Random intersection tests
+// -----------------------------------
+
+// Checkers (partial specialization for epic & epec)
+template <class K>
+struct Checker
+{
+  template <typename Query>
+  void operator()(const Query& q, const typename K::Triangle_3& t) const
+  {
+    CGAL::Object result = CGAL::intersection(q, t);
+    
+    if ( ! result.empty() )
+    { 
+      assert( CGAL::do_intersect(q, t) );     
+      assert(   NULL != CGAL::object_cast<typename K::Point_3>(&result) 
+             || NULL != CGAL::object_cast<typename K::Segment_3>(&result));
+    }
+  }
+};
+
+template <>
+struct Checker<Epic>
+{
+  typedef Epic K;
+  
+  template <typename Query>
+  void operator()(const Query& q, const typename K::Triangle_3& t) const
+  {
+    CGAL::Object result = CGAL::intersection(q, t);
+    
+    if ( ! result.empty() )
+    { 
+      assert( CGAL::do_intersect(q, t) );
+      assert(   NULL != CGAL::object_cast<typename K::Point_3>(&result) 
+             || NULL != CGAL::object_cast<typename K::Segment_3>(&result));
+    }
+  }
+  
+  void operator()(const K::Line_3& l, const K::Triangle_3& t) const
+  {
+    CGAL::Object result = CGAL::intersection(l, t);
+    
+    if ( ! result.empty() )
+    { 
+      // Here we can't check do_intersect, because there are constructions when
+      // building points on line
+      assert(   NULL != CGAL::object_cast<K::Point_3>(&result) 
+             || NULL != CGAL::object_cast<K::Segment_3>(&result));
+    }
+  }
+};
+
+template <>
+struct Checker<Epec>
+{
+  typedef Epec K;
+  
+  template <typename Query>
+  void operator()(const Query& q, const typename K::Triangle_3& t) const
+  {
+    typedef typename K::Point_3 Point_3;
+    typedef typename K::Segment_3 Segment_3;
+    
+    CGAL::Object result = CGAL::intersection(q, t);
+    
+    if ( ! result.empty() )
+    { 
+      assert( CGAL::do_intersect(q, t) );
+      
+      // Verify answer is correct
+      const Point_3* p = CGAL::object_cast<Point_3>(&result);
+      const Segment_3* s = CGAL::object_cast<Segment_3>(&result);
+      
+      assert(   (NULL!=p && t.has_on(*p) && q.has_on(*p))
+             || (NULL!=s && t.has_on(s->source()) && t.has_on(s->target())
+                         && q.has_on(s->source()) && q.has_on(s->target())) );
+    }
+    else
+    {
+      assert ( !CGAL::do_intersect(q, t) );
+    }
+  }
+};
+
+
+
+// random number generation
+double random_in(const double a,
+                 const double b)
+{
+  double r = rand() / (double)RAND_MAX;
+  return a + (b - a) * r;
+}
+
+template <class K>
+typename K::Point_3 random_point_in(const CGAL::Bbox_3& bbox)
+{
+  typedef typename K::FT FT;
+  FT x = (FT)random_in(bbox.xmin(),bbox.xmax());
+  FT y = (FT)random_in(bbox.ymin(),bbox.ymax());
+  FT z = (FT)random_in(bbox.zmin(),bbox.zmax());
+  return typename K::Point_3(x,y,z);
+}
+
+// random_test()
+template <class K>
+void random_test()
+{
+  typedef typename K::Point_3 Point;
+  typedef typename K::Segment_3 Segment;
+  typedef typename K::Ray_3 Ray;
+  typedef typename K::Line_3 Line;
+  typedef typename K::Triangle_3 Triangle;
+  typedef typename K::Plane_3 Plane;
+  
+  Checker<K> check;
+  
+  double box_size = 1e12;
+  CGAL::Bbox_3 bbox(-box_size,-box_size,-box_size,box_size,box_size,box_size);
+  
+  // Use 10 triangles, 100 queries for each triangle
+  for ( int i=0 ; i<10 ; ++i )
+  {
+    Triangle t(random_point_in<K>(bbox),
+               random_point_in<K>(bbox),
+               random_point_in<K>(bbox));
+    
+    Plane p = t.supporting_plane();
+    
+    for ( int j=0 ; j<100 ; ++j )
+    {
+      Point a = random_point_in<K>(bbox);
+      Point b = random_point_in<K>(bbox);
+      
+      Segment s (a,b);
+      Ray r(a,b);
+      Line l (a,b);
+      
+      check(s,t);
+      check(r,t);
+      check(l,t);
+      
+      // Project points on triangle plane to have degenerate queries
+      Point c = p.projection(a);
+      Point d = p.projection(b);
+      
+      Segment s2 (c,d);
+      Ray r2 (c,d);
+      Line l2 (c,d);
+      
+      check(s2,t);
+      check(r2,t);
+      check(l2,t);
+    }
+  }
+}
+
+
+
+// -----------------------------------
+// Precomputed results test
+// -----------------------------------
 template <class Triangle, class Query, class Result>
 bool test_aux(const Triangle t,
               const Query& q,
@@ -47,7 +225,7 @@ bool test_aux(const Triangle t,
   else
   {
     std::cout << "ERROR: intersection(" << name
-              << ") did not answer the expected result !";  
+    << ") did not answer the expected result !";  
     
     if ( NULL != pr )
       std::cout << " (answer: ["<< *pr << "])";
@@ -208,6 +386,105 @@ bool test()
   b &= test_aux(t,sa8,"t-sa8",p8);
   b &= test_aux(t,sb2,"t-sb2",p2);
   
+  // -----------------------------------
+  // ray queries
+  // -----------------------------------
+  // Edges of t 
+  Ray r12(p1,p2);
+  Ray r21(p2,p1);
+  Ray r13(p1,p3);
+  Ray r23(p2,p3);
+  
+  b &= test_aux(t,r12,"t-r12",s12);
+  b &= test_aux(t,r21,"t-r21",s21);
+  b &= test_aux(t,r13,"t-r13",s13);
+  b &= test_aux(t,r23,"t-r23",s23);
+  
+  // In triangle
+  Point p9_(FT(0.), FT(0.5), FT(0.5));
+  Point p9(FT(0.25), FT(0.375), FT(0.375));
+  
+  Ray r14(p1,p4);
+  Ray r41(p4,p1);
+  Ray r24(p2,p4);
+  Ray r42(p4,p2);
+  Ray r15(p1,p5);
+  Ray r25(p2,p5);
+  Ray r34(p3,p4);
+  Ray r35(p3,p5);
+  Ray r36(p3,p6);
+  Ray r45(p4,p5);
+  Ray r16(p1,p6);
+  Ray r26(p2,p6);
+  Ray r62(p6,p2);
+  Ray r46(p4,p6);
+  Ray r48(p4,p8);
+  Ray r56(p5,p6);
+  Ray r47(p4,p7);
+  Ray r89(p8,p9);
+  Ray r86(p8,p6);
+  Ray r68(p6,p8);
+  Segment r89_res(p8,p9_);
+  
+  b &= test_aux(t,r14,"t-r14",s12);
+  b &= test_aux(t,r41,"t-r41",s41);
+  b &= test_aux(t,r24,"t-r24",s21);
+  b &= test_aux(t,r42,"t-r42",s42);
+  b &= test_aux(t,r15,"t-r15",s15);
+  b &= test_aux(t,r25,"t-r25",s23);
+  b &= test_aux(t,r34,"t-r34",s34);
+  b &= test_aux(t,r35,"t-r35",s32);
+  b &= test_aux(t,r36,"t-r36",s31);
+  b &= test_aux(t,r45,"t-r45",s45);
+  b &= test_aux(t,r16,"t-r16",s13);
+  b &= test_aux(t,r26,"t-r26",s26);
+  b &= test_aux(t,r62,"t-r62",s62);
+  b &= test_aux(t,r46,"t-r46",s46);
+  b &= test_aux(t,r48,"t-r48",s46);
+  b &= test_aux(t,r56,"t-r56",s56);
+  b &= test_aux(t,r47,"t-r47",s45);
+  b &= test_aux(t,r89,"t-t89",r89_res);
+  b &= test_aux(t,r68,"t-r68",s64);
+  b &= test_aux(t,r86,"t-r86",s86);
+  
+  
+  // Outside points (in triangre prane)
+  Ray rAB(pA,pB);
+  Ray rBC(pB,pC);
+  Ray r2E(p2,pE);
+  Ray rE2(pE,p2);
+  Ray r2A(p2,pA);
+  Ray r6E(p6,pE);
+  Ray rB8(pB,p8);
+  Ray rC8(pC,p8);
+  Ray r8C(p8,pC);
+  Ray r1F(p1,pF);
+  Ray rF6(pF,p6);
+  
+  b &= test_aux(t,rAB,"t-rAB",p2);
+  b &= test_aux(t,rBC,"t-rBC",s46);
+  b &= test_aux(t,r2E,"t-r2E",s26);
+  b &= test_aux(t,rE2,"t-rE2",s62);
+  b &= test_aux(t,r2A,"t-r2A",p2);
+  b &= test_aux(t,r6E,"t-r6E",p6);
+  b &= test_aux(t,rB8,"t-rB8",s46);
+  b &= test_aux(t,rC8,"t-rC8",s64);
+  b &= test_aux(t,r8C,"t-r8C",s86);
+  b &= test_aux(t,r1F,"t-r1F",s13);
+  b &= test_aux(t,rF6,"t-rF6",s31);
+  
+  // Outside triangle plane
+  Ray rab(pa,pb);
+  Ray rac(pa,pc);
+  Ray rae(pa,pe);
+  Ray ra8(pa,p8);
+  Ray rb2(pb,p2);
+  
+  b &= test_aux(t,rab,"t-rab",p1);
+  b &= test_aux(t,rac,"t-rac",p6);
+  b &= test_aux(t,rae,"t-rae",p8);
+  b &= test_aux(t,ra8,"t-ra8",p8);
+  b &= test_aux(t,rb2,"t-rb2",p2);
   
   // -----------------------------------
   // Line queries
@@ -224,9 +501,6 @@ bool test()
   b &= test_aux(t,l23,"t-l23",s23);
   
   // In triangle
-  Point p9_(FT(0.), FT(0.5), FT(0.5));
-  Point p9(FT(0.25), FT(0.375), FT(0.375));
-  
   Line l14(p1,p4);
   Line l41(p4,p1);
   Line l24(p2,p4);
@@ -247,7 +521,8 @@ bool test()
   Line l89(p8,p9);
   Line l86(p8,p6);
   Line l68(p6,p8);
-  Segment s89_res(p1,p9_);
+  Segment l89_res(p1,p9_);
+
   
   b &= test_aux(t,l14,"t-l14",s12);
   b &= test_aux(t,l41,"t-l41",s21);
@@ -266,7 +541,7 @@ bool test()
   b &= test_aux(t,l48,"t-l48",s46);
   b &= test_aux(t,l56,"t-l56",s56);
   b &= test_aux(t,l47,"t-l47",s45);
-  b &= test_aux(t,l89,"t-t89",s89_res);
+  b &= test_aux(t,l89,"t-t89",l89_res);
   b &= test_aux(t,l68,"t-l68",s64);
   b &= test_aux(t,l86,"t-l86",s46);
 
@@ -313,25 +588,56 @@ bool test()
 	return b;
 }
 
+
+
+// -----------------------------------
+// Main
+// -----------------------------------
 int main()
 {
-  std::cout << "Testing with Simple_cartesian<float>..." << std::endl ;
+  // -----------------------------------
-  bool b = test<CGAL::Simple_cartesian<float> >();
+  // Test intersection results
+  // -----------------------------------
+  std::cout << "Test precomputed intersection results" << std::endl;
+  std::cout << "\tTesting with Simple_cartesian<float>..." << std::endl ;
+  bool b = test<Sc_f>();
   
-  std::cout << "Testing with Simple_cartesian<double>..." << std::endl ;
+  std::cout << "\tTesting with Simple_cartesian<double>..." << std::endl ;
-	b &= test<CGAL::Simple_cartesian<double> >();
+	b &= test<Sc_d>();
   
-  std::cout << "Testing with Cartesian<float>..." << std::endl ;
+  std::cout << "\tTesting with Cartesian<float>..." << std::endl ;
-	b &= test<CGAL::Cartesian<float> >();
+	b &= test<C_f>();
   
-  std::cout << "Testing with Cartesian<double>..." << std::endl ;
+  std::cout << "\tTesting with Cartesian<double>..." << std::endl ;
-	b &= test<CGAL::Cartesian<double> >();
+	b &= test<C_d>();
   
-  std::cout << "Testing with Exact_predicates_inexact_constructions_kernel..." << std::endl ;
+  std::cout << "\tTesting with Exact_predicates_inexact_constructions_kernel..." << std::endl ;
-  b &= test<CGAL::Exact_predicates_inexact_constructions_kernel>();
+  b &= test<Epic>();
 	
-  std::cout << "Testing with Exact_predicates_exact_constructions_kernel..." << std::endl ;
+  std::cout << "\tTesting with Exact_predicates_exact_constructions_kernel..." << std::endl ;
-  b &= test<CGAL::Exact_predicates_exact_constructions_kernel>();
+  b &= test<Epec>();
+  
+  // -----------------------------------
+  // Test random intersection
+  // -----------------------------------
+  std::cout << std::endl << "Test random intersections" << std::endl;
+  std::cout << "\tTesting with Simple_cartesian<float>..." << std::endl ;
+  random_test<Sc_f>();
+  
+  std::cout << "\tTesting with Simple_cartesian<double>..." << std::endl ;
+	random_test<Sc_d>();
+  
+  std::cout << "\tTesting with Cartesian<float>..." << std::endl ;
+	random_test<C_f>();
+  
+  std::cout << "\tTesting with Cartesian<double>..." << std::endl ;
+	random_test<C_d>();
+  
+  std::cout << "\tTesting with Exact_predicates_inexact_constructions_kernel..." << std::endl ;
+  random_test<Epic>();
+	
+  std::cout << "\tTesting with Exact_predicates_exact_constructions_kernel..." << std::endl ;
+  random_test<Epec>();
   
   if ( b )
     return EXIT_SUCCESS;