diff --git a/AABB_tree/include/CGAL/AABB_intersections/Bbox_3_line_3_do_intersect.h b/AABB_tree/include/CGAL/AABB_intersections/Bbox_3_line_3_do_intersect.h index 22dcd8d772b..29610ad1d6d 100644 --- a/AABB_tree/include/CGAL/AABB_intersections/Bbox_3_line_3_do_intersect.h +++ b/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_b = line.point(1); + const Point_3 point = line.point(); + const Vector_3 v = line.to_vector(); return bbox_line_do_intersect_aux( - point_a.x(), point_a.y(), point_a.z(), - point_b.x(), point_b.y(), point_b.z(), + point.x(), point.y(), point.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()) ); } diff --git a/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_line_3_intersection.h b/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_line_3_intersection.h index ac56c56cb62..5ba07be4857 100644 --- a/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_line_3_intersection.h +++ b/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_line_3_intersection.h @@ -31,6 +31,47 @@ namespace CGAL { namespace internal { +template +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 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 = object_cast(&l_bc_obj); - 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(&l_ca_obj); - if ( NULL == l_ca ) - { - CGAL_kernel_assertion(false); - return Object(); - } + const Point_3 l_bc = t3l3_intersection_coplanar_aux(l,b,c,k); + const Point_3 l_ca = t3l3_intersection_coplanar_aux(l,c,a,k); 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 - // triangle in the common plane + // Test whether the 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 @@ -298,8 +324,24 @@ intersection_coplanar(const typename K::Triangle_3 &t, } +template +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(&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); diff --git a/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_ray_3_intersection.h b/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_ray_3_intersection.h index f97f60c95a5..17d2f3b3e58 100644 --- a/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_ray_3_intersection.h +++ b/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_ray_3_intersection.h @@ -26,39 +26,434 @@ namespace CGAL { namespace internal { - - + +template +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 Object -intersection_triangle_ray_aux(const typename K::Ray_3 &r, - const typename K::Segment_3 &s, - const K& k) +t3r3_intersection_coplanar_aux(const typename K::Point_3& a, + const typename K::Point_3& b, + 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_object_3 make_object = k.construct_object_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(); - // typename K::Compare_xyz_3 compare = k.compare_xyz_3_object(); // PA: never used - typename K::Equal_3 equal = k.equal_3_object(); - - const Point_3& p = r.source(); - const Point_3& q1 = s.source(); - 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)); + typename K::Construct_point_on_3 point_on = + k.construct_point_on_3_object(); + + const Point_3& p = point_on(r,0); + + // A ray is not symetric, 2 cases depending on isolated side of c + Orientation cap; + if ( negative_side ) + cap = coplanar_orientation(c,a,p); 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 +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 +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(&obj) ) + return obj; + else + return Object(); } - template 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; - typedef typename K::Ray_3 Ray_3; - typedef typename K::Line_3 Line_3; + + typename K::Construct_vertex_3 vertex_on = + k.construct_vertex_3_object(); + + 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) ) - return Object_3(); + Point_3 d = point_on(ray(a,r.to_vector()),1); + + 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 - // do_intersect between the triangle and the ray - const Object intersection = CGAL::intersection(t.supporting_plane(), - r); - - // intersection is either a point, either a ray - // if it is a ray, then we need to clip it to a segment - if ( const Ray_3* ray = object_cast(&intersection)) - { - typename K::Intersect_3 intersect = k.intersect_3_object(); - typename K::Construct_line_3 f_line = k.construct_line_3_object(); - typename K::Construct_vertex_3 vertex_on = k.construct_vertex_3_object(); - typename K::Construct_object_3 make_object = k.construct_object_3_object(); - typename K::Construct_segment_3 f_segment = k.construct_segment_3_object(); - typename K::Has_on_3 has_on = k.has_on_3_object(); - - const Point_3& t0 = vertex_on(t,0); - const Point_3& t1 = vertex_on(t,1); - const Point_3& t2 = vertex_on(t,2); - - const Line_3 l01 = f_line(t0,t1); - const Line_3 l02 = f_line(t0,t2); - const Line_3 l12 = f_line(t1,t2); - - const Segment_3 s01 = f_segment(t0,t1); - const Segment_3 s02 = f_segment(t0,t2); - const Segment_3 s12 = f_segment(t1,t2); - - const Line_3 lr = ray->supporting_line(); - - const Object_3 inter_01 = intersect(l01,lr); - const Object_3 inter_02 = intersect(l02,lr); - const Object_3 inter_12 = intersect(l12,lr); - - const Point_3* p01 = object_cast(&inter_01); - const Point_3* p02 = object_cast(&inter_02); - const Point_3* p12 = object_cast(&inter_12); - - if ( p01 && has_on(s01, *p01) ) - { - if ( p02 && has_on(s02, *p02) ) - return intersection_triangle_ray_aux(*ray, f_segment(*p01,*p02), k); - else if ( p12 && has_on(s12, *p12) ) - return intersection_triangle_ray_aux(*ray, f_segment(*p01,*p12), k); - else - return make_object(*p01); - } - else if ( p02 && has_on(s02, *p02) ) - { - if ( p12 && has_on(s12, *p12) ) - return intersection_triangle_ray_aux(*ray, f_segment(*p02,*p12), k); - else - return make_object(*p02); - } - else if ( p12 && has_on(s12, *p12) ) - { - return make_object(*p12); - } - - // should not happen - CGAL_kernel_assertion(false); - return Object_3(); - } - else - return intersection; + switch ( abcp ) { + case POSITIVE: + switch ( ray_direction ) { + case POSITIVE: + // the ray lies in the positive open halfspaces defined by the + // triangle's supporting plane + return Object(); + + case NEGATIVE: + // The ray straddles the triangle's plane + // p sees the triangle in counterclockwise 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 a plane parallel to a,b,c support plane + return Object(); + + default: // should not happen. + CGAL_kernel_assertion(false); + return Object(); + } + + case NEGATIVE: + switch ( ray_direction ) { + case POSITIVE: + // The ray straddles the triangle's plane + // q sees the triangle in counterclockwise order + + if ( orientation(q,p,a,b) != POSITIVE + && orientation(q,p,b,c) != POSITIVE + && orientation(q,p,c,a) != POSITIVE ) + return t3r3_intersection_aux(t,r,k); + else + return Object(); + + case NEGATIVE: + // the ray lies in the negative open halfspaces defined by the + // triangle's supporting plane + return Object(); + + case COPLANAR: + // The ray lie in a plane parallel to a,b,c support plane + return Object(); + + default: // should not happen. + CGAL_kernel_assertion(false); + return Object(); + } + + case COPLANAR: // p belongs to the triangle's supporting plane + switch ( ray_direction ) { + case POSITIVE: + // q sees the triangle in counterclockwise order + if ( orientation(q,p,a,b) != POSITIVE + && orientation(q,p,b,c) != POSITIVE + && orientation(q,p,c,a) != POSITIVE ) + return t3r3_intersection_aux(t,r,k); + else + return Object(); + + 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 diff --git a/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_segment_3_intersection.h b/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_segment_3_intersection.h index 8e35ab9c582..06bddb845d9 100644 --- a/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_segment_3_intersection.h +++ b/AABB_tree/include/CGAL/AABB_intersections/Triangle_3_segment_3_intersection.h @@ -33,6 +33,47 @@ namespace CGAL { namespace internal { +template +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 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)); - const Point_3* p_temp = object_cast(&obj); - if ( NULL != p_temp ) - p_side_end_point = *p_temp; + p_side_end_point = t3s3_intersection_coplanar_aux(p,q,b,c,k); } Point_3 q_side_end_point(q); if ( NEGATIVE == coplanar_orientation(c,a,q) ) { - Object obj = intersection(line(p,q),line(c,a)); - const Point_3* p_temp = object_cast(&obj); - if ( NULL != p_temp ) - q_side_end_point = *p_temp; + q_side_end_point = t3s3_intersection_coplanar_aux(p,q,c,a,k); } if ( negative_side ) @@ -114,7 +143,7 @@ t3s3_intersection_coplanar_aux(const typename K::Point_3& a, template 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(&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(&obj) ) + return obj; + else + return Object(); + } else return Object(); diff --git a/AABB_tree/test/AABB_tree/aabb_intersection_test.cpp b/AABB_tree/test/AABB_tree/aabb_intersection_test.cpp index 41412041bed..53a4af87d34 100644 --- a/AABB_tree/test/AABB_tree/aabb_intersection_test.cpp +++ b/AABB_tree/test/AABB_tree/aabb_intersection_test.cpp @@ -31,6 +31,184 @@ #include + + +// ----------------------------------- +// Kernels +// ----------------------------------- +struct Sc_f : public CGAL::Simple_cartesian {}; +struct Sc_d : public CGAL::Simple_cartesian {}; +struct C_f : public CGAL::Cartesian {}; +struct C_d : public CGAL::Cartesian {}; +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 +struct Checker +{ + template + 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(&result) + || NULL != CGAL::object_cast(&result)); + } + } +}; + +template <> +struct Checker +{ + typedef Epic K; + + template + 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(&result) + || NULL != CGAL::object_cast(&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(&result) + || NULL != CGAL::object_cast(&result)); + } + } +}; + +template <> +struct Checker +{ + typedef Epec K; + + template + 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(&result); + const Segment_3* s = CGAL::object_cast(&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 +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 +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 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(bbox), + random_point_in(bbox), + random_point_in(bbox)); + + Plane p = t.supporting_plane(); + + for ( int j=0 ; j<100 ; ++j ) + { + Point a = random_point_in(bbox); + Point b = random_point_in(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 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..." << std::endl ; - bool b = test >(); + // ----------------------------------- + // Test intersection results + // ----------------------------------- + std::cout << "Test precomputed intersection results" << std::endl; + std::cout << "\tTesting with Simple_cartesian..." << std::endl ; + bool b = test(); - std::cout << "Testing with Simple_cartesian..." << std::endl ; - b &= test >(); + std::cout << "\tTesting with Simple_cartesian..." << std::endl ; + b &= test(); - std::cout << "Testing with Cartesian..." << std::endl ; - b &= test >(); + std::cout << "\tTesting with Cartesian..." << std::endl ; + b &= test(); - std::cout << "Testing with Cartesian..." << std::endl ; - b &= test >(); + std::cout << "\tTesting with Cartesian..." << std::endl ; + b &= test(); - std::cout << "Testing with Exact_predicates_inexact_constructions_kernel..." << std::endl ; - b &= test(); + std::cout << "\tTesting with Exact_predicates_inexact_constructions_kernel..." << std::endl ; + b &= test(); - std::cout << "Testing with Exact_predicates_exact_constructions_kernel..." << std::endl ; - b &= test(); + std::cout << "\tTesting with Exact_predicates_exact_constructions_kernel..." << std::endl ; + b &= test(); + + // ----------------------------------- + // Test random intersection + // ----------------------------------- + std::cout << std::endl << "Test random intersections" << std::endl; + std::cout << "\tTesting with Simple_cartesian..." << std::endl ; + random_test(); + + std::cout << "\tTesting with Simple_cartesian..." << std::endl ; + random_test(); + + std::cout << "\tTesting with Cartesian..." << std::endl ; + random_test(); + + std::cout << "\tTesting with Cartesian..." << std::endl ; + random_test(); + + std::cout << "\tTesting with Exact_predicates_inexact_constructions_kernel..." << std::endl ; + random_test(); + + std::cout << "\tTesting with Exact_predicates_exact_constructions_kernel..." << std::endl ; + random_test(); if ( b ) return EXIT_SUCCESS;