mirror of https://github.com/CGAL/cgal
Merge pull request #6779 from sloriot/Intersections_3-better_tetra_tri
Simpler implementation of Tet vs. Tri intersection
This commit is contained in:
Laurent Rineau
commit
4bc403857d
|
|
@ -15,8 +15,7 @@
|
||||||
#ifndef CGAL_INTERNAL_INTERSECTIONS_3_TETRAHEDRON_3_TRIANGLE_3_INTERSECTIONS_H
|
#ifndef CGAL_INTERNAL_INTERSECTIONS_3_TETRAHEDRON_3_TRIANGLE_3_INTERSECTIONS_H
|
||||||
#define CGAL_INTERNAL_INTERSECTIONS_3_TETRAHEDRON_3_TRIANGLE_3_INTERSECTIONS_H
|
#define CGAL_INTERNAL_INTERSECTIONS_3_TETRAHEDRON_3_TRIANGLE_3_INTERSECTIONS_H
|
||||||
|
|
||||||
#include <CGAL/Intersections_3/internal/Plane_3_Tetrahedron_3_intersection.h>
|
#include <CGAL/Intersections_3/internal/Line_3_Plane_3_intersection.h>
|
||||||
#include <CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h>
|
|
||||||
|
|
||||||
#include <CGAL/kernel_basic.h>
|
#include <CGAL/kernel_basic.h>
|
||||||
|
|
||||||
|
|
@ -25,260 +24,12 @@
|
||||||
#include <list>
|
#include <list>
|
||||||
#include <utility>
|
#include <utility>
|
||||||
#include <vector>
|
#include <vector>
|
||||||
|
#include <bitset>
|
||||||
|
|
||||||
namespace CGAL {
|
namespace CGAL {
|
||||||
namespace Intersections {
|
namespace Intersections {
|
||||||
namespace internal {
|
namespace internal {
|
||||||
|
|
||||||
template<typename Segment>
|
|
||||||
void filter_segments(std::list<Segment>& segments)
|
|
||||||
{
|
|
||||||
auto are_equal = [](const Segment& l, const Segment& r) -> bool
|
|
||||||
{
|
|
||||||
return (l == r || l == r.opposite());
|
|
||||||
};
|
|
||||||
|
|
||||||
auto it = std::unique(segments.begin(), segments.end(), are_equal);
|
|
||||||
segments.erase(it, segments.end());
|
|
||||||
}
|
|
||||||
|
|
||||||
// plane going through the ref segment's source, a point above (given by the normal of the input
|
|
||||||
// triangle) and two points (ref_other / query) that need to be ordered
|
|
||||||
template <class K>
|
|
||||||
bool first_comes_first_pt(const typename K::Point_3& ref_source,
|
|
||||||
const typename K::Point_3& ref_z,
|
|
||||||
const typename K::Point_3& ref_other,
|
|
||||||
const typename K::Point_3& query,
|
|
||||||
const K& k)
|
|
||||||
{
|
|
||||||
typename K::Orientation_3 orientation = k.orientation_3_object();
|
|
||||||
|
|
||||||
// points have filtered to remove segments' extremities
|
|
||||||
CGAL_precondition(ref_other != query);
|
|
||||||
|
|
||||||
const Orientation o = orientation(ref_source, ref_z, ref_other, query);
|
|
||||||
CGAL_assertion(o != COPLANAR);
|
|
||||||
|
|
||||||
// ref_other comes first <==> query is on the positive side of the plane
|
|
||||||
return (o == POSITIVE);
|
|
||||||
}
|
|
||||||
|
|
||||||
template <class K, class SegPtVariant>
|
|
||||||
bool first_comes_first(const typename K::Point_3& ref_source,
|
|
||||||
const typename K::Point_3& ref_z,
|
|
||||||
const typename K::Point_3& ref_other,
|
|
||||||
const SegPtVariant& seg_or_pt,
|
|
||||||
const K& k)
|
|
||||||
{
|
|
||||||
typedef typename K::Point_3 Point_3;
|
|
||||||
typedef typename K::Segment_3 Segment_3;
|
|
||||||
|
|
||||||
typedef typename std::list<Segment_3>::iterator SCI;
|
|
||||||
typedef typename std::vector<Point_3>::iterator PCI;
|
|
||||||
|
|
||||||
if(seg_or_pt.which() == 0)
|
|
||||||
{
|
|
||||||
const Segment_3& s = *(boost::get<SCI>(seg_or_pt));
|
|
||||||
return first_comes_first_pt(ref_source, ref_z, ref_other, s.source(), k);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
CGAL_assertion(seg_or_pt.which() == 1);
|
|
||||||
|
|
||||||
const Point_3& p = *(boost::get<PCI>(seg_or_pt));
|
|
||||||
return first_comes_first_pt(ref_source, ref_z, ref_other, p, k);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
template <class K, class SegmentContainer, class PointContainer>
|
|
||||||
typename Intersection_traits<K, typename K::Tetrahedron_3, typename K::Triangle_3>::result_type
|
|
||||||
build_intersection(const typename K::Tetrahedron_3& /*input_tetrahedron*/,
|
|
||||||
const typename K::Triangle_3& input_triangle,
|
|
||||||
PointContainer& points,
|
|
||||||
SegmentContainer& segments,
|
|
||||||
const K& k)
|
|
||||||
{
|
|
||||||
typedef typename Intersection_traits<K, typename K::Tetrahedron_3, typename K::Triangle_3>::result_type result_type;
|
|
||||||
|
|
||||||
typedef typename K::Point_3 Point_3;
|
|
||||||
typedef typename K::Segment_3 Segment_3;
|
|
||||||
typedef typename K::Vector_3 Vector_3;
|
|
||||||
typedef typename K::Triangle_3 Triangle_3;
|
|
||||||
typedef std::vector<Point_3> Poly;
|
|
||||||
|
|
||||||
typedef typename SegmentContainer::iterator SCI;
|
|
||||||
typedef typename PointContainer::iterator PCI;
|
|
||||||
|
|
||||||
// @todo? Could do the 1 segment case with this code too...
|
|
||||||
CGAL_precondition(segments.size() >= 2 && segments.size() <= 4);
|
|
||||||
CGAL_precondition(points.size() <= 2);
|
|
||||||
|
|
||||||
// Constructions @fixme avoidable?
|
|
||||||
const Vector_3 input_triangle_normal = input_triangle.supporting_plane().orthogonal_vector();
|
|
||||||
|
|
||||||
// remove points that are just segments extremities
|
|
||||||
auto is_extremity = [&segments](const Point_3& p) -> bool
|
|
||||||
{
|
|
||||||
for(const Segment_3& s : segments)
|
|
||||||
if(p == s.source() || p == s.target())
|
|
||||||
return true;
|
|
||||||
return false;
|
|
||||||
};
|
|
||||||
points.erase(std::remove_if(points.begin(), points.end(), is_extremity),
|
|
||||||
points.end());
|
|
||||||
|
|
||||||
// Take the first segment as reference, and order the rest to form a convex polygon
|
|
||||||
//
|
|
||||||
// All segments and points involved in the intersection are on the input triangle
|
|
||||||
// and thus everything is coplanar, at least theoretically (the kernel might not provide
|
|
||||||
// exact constructions...)
|
|
||||||
//
|
|
||||||
// Given an arbitrary segment, the code below sorts the other segments and the points
|
|
||||||
// in a ccw order. Using a vector because the number of segments and points is bounded
|
|
||||||
// (max 4 segments and max 2 points) so even if the linear insertion is a little ugly,
|
|
||||||
// it is not expensive anyway.
|
|
||||||
//
|
|
||||||
// Example:
|
|
||||||
/*
|
|
||||||
x p0
|
|
||||||
|
|
||||||
/
|
|
||||||
/
|
|
||||||
s1 / \ s2
|
|
||||||
/ \
|
|
||||||
--------------
|
|
||||||
s0
|
|
||||||
*/
|
|
||||||
//
|
|
||||||
// s0 is chosen as the reference segment
|
|
||||||
// output will be 's0 s2 p0 s1'
|
|
||||||
|
|
||||||
Segment_3& ref_s = segments.front();
|
|
||||||
Point_3 ref_z = ref_s.source() + input_triangle_normal;
|
|
||||||
|
|
||||||
// The reference segment should be such that all other intersection parts are
|
|
||||||
// on the positive side of the plane described by the normal of the triangle and ref_s
|
|
||||||
bool swapped = false;
|
|
||||||
for(SCI slit = std::next(segments.begin()); slit != segments.end(); ++slit)
|
|
||||||
{
|
|
||||||
const Segment_3& other = *slit;
|
|
||||||
|
|
||||||
if(k.orientation_3_object()(ref_s.source(), ref_z, ref_s.target(), other.source()) == CGAL::NEGATIVE ||
|
|
||||||
k.orientation_3_object()(ref_s.source(), ref_z, ref_s.target(), other.target()) == CGAL::NEGATIVE)
|
|
||||||
{
|
|
||||||
ref_s = ref_s.opposite();
|
|
||||||
ref_z = ref_s.source() + input_triangle_normal;
|
|
||||||
swapped = true;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
if(!swapped)
|
|
||||||
{
|
|
||||||
for(PCI plit = points.begin(); plit != points.end(); ++plit)
|
|
||||||
{
|
|
||||||
const Point_3& other = *plit;
|
|
||||||
if(k.orientation_3_object()(ref_s.source(), ref_z, ref_s.target(), other) == CGAL::NEGATIVE)
|
|
||||||
{
|
|
||||||
swapped = true;
|
|
||||||
ref_s = ref_s.opposite();
|
|
||||||
ref_z = ref_s.source() + input_triangle_normal;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
const Point_3& ref_sp = ref_s.source();
|
|
||||||
const Point_3& ref_tp = ref_s.target();
|
|
||||||
|
|
||||||
// Now, order the other parts of the intersection
|
|
||||||
std::list<boost::variant<SCI, PCI> > res_elements; // iterators to the points/segments
|
|
||||||
res_elements.emplace_back(segments.begin());
|
|
||||||
|
|
||||||
for(SCI slit = std::next(segments.begin()); slit != segments.end(); ++slit)
|
|
||||||
{
|
|
||||||
// first, check if the segment is well oriented, meaning its source comes before its target (ccwly)
|
|
||||||
Segment_3& curr_s = *slit;
|
|
||||||
|
|
||||||
if(curr_s.source() == ref_sp || curr_s.target() == ref_tp) // consecutive segments must have consistent orientation
|
|
||||||
{
|
|
||||||
curr_s = curr_s.opposite();
|
|
||||||
}
|
|
||||||
else if(curr_s.source() == ref_tp || curr_s.target() == ref_sp)
|
|
||||||
{
|
|
||||||
// nothing to do here as we know that sp&tp are on the positive side of (normal, ref_s)
|
|
||||||
}
|
|
||||||
else if(first_comes_first_pt(ref_sp, ref_z, curr_s.target(), curr_s.source(), k))
|
|
||||||
{
|
|
||||||
curr_s = curr_s.opposite();
|
|
||||||
}
|
|
||||||
|
|
||||||
// Find where the current segment fit in the final polygon intersection
|
|
||||||
for(auto rit = std::next(res_elements.begin()); ; ++rit)
|
|
||||||
{
|
|
||||||
// always pick the current segment's source to ensure ref_source != ref_other
|
|
||||||
if(rit == res_elements.end() || first_comes_first(ref_sp, ref_z, curr_s.source(), *rit, k))
|
|
||||||
{
|
|
||||||
res_elements.insert(rit, slit);
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
for(PCI plit = points.begin(); plit != points.end(); ++plit)
|
|
||||||
{
|
|
||||||
const Point_3& curr_p = *plit;
|
|
||||||
|
|
||||||
// Find where the current point fits in the boundary of the polygon intersection
|
|
||||||
for(auto rit = std::next(res_elements.begin()); ; ++rit)
|
|
||||||
{
|
|
||||||
if(rit == res_elements.end() || first_comes_first(ref_sp, ref_z, curr_p, *rit, k))
|
|
||||||
{
|
|
||||||
res_elements.insert(rit, plit);
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
CGAL_postcondition(res_elements.size() == points.size() + segments.size());
|
|
||||||
|
|
||||||
// Concatenate points to create the polygonal output
|
|
||||||
Poly res;
|
|
||||||
for(const boost::variant<SCI, PCI>& e : res_elements)
|
|
||||||
{
|
|
||||||
if(const SCI* sci = boost::get<SCI>(&e))
|
|
||||||
{
|
|
||||||
const Segment_3& s = **sci;
|
|
||||||
|
|
||||||
if(res.empty() || s.source() != res.back()) // common extremity for consecutive segments
|
|
||||||
res.push_back(s.source());
|
|
||||||
if(res.empty() || s.target() != res.front())
|
|
||||||
res.push_back(s.target());
|
|
||||||
}
|
|
||||||
else if(const PCI* pci = boost::get<PCI>(&e))
|
|
||||||
{
|
|
||||||
res.push_back(**pci);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
CGAL_assertion(false);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
CGAL_assertion(std::set<Point_3>(res.begin(), res.end()).size() == res.size());
|
|
||||||
CGAL_assertion(res.size() >= 3);
|
|
||||||
|
|
||||||
if(res.size() == 3)
|
|
||||||
{
|
|
||||||
Triangle_3 tr { res[0], res[1], res[2] };
|
|
||||||
return result_type(std::forward<Triangle_3>(tr));
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
return result_type(std::forward<Poly>(res));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
template <class K>
|
template <class K>
|
||||||
typename Intersection_traits<K, typename K::Tetrahedron_3, typename K::Triangle_3>::result_type
|
typename Intersection_traits<K, typename K::Tetrahedron_3, typename K::Triangle_3>::result_type
|
||||||
intersection(const typename K::Tetrahedron_3& tet,
|
intersection(const typename K::Tetrahedron_3& tet,
|
||||||
|
|
@ -286,319 +37,131 @@ intersection(const typename K::Tetrahedron_3& tet,
|
||||||
const K& k)
|
const K& k)
|
||||||
{
|
{
|
||||||
typedef typename Intersection_traits<K, typename K::Tetrahedron_3, typename K::Triangle_3>::result_type result_type;
|
typedef typename Intersection_traits<K, typename K::Tetrahedron_3, typename K::Triangle_3>::result_type result_type;
|
||||||
typedef typename Intersection_traits<K, typename K::Triangle_3, typename K::Triangle_3>::result_type Inter_type;
|
|
||||||
|
|
||||||
CGAL_precondition(!tet.is_degenerate());
|
CGAL_precondition(!tet.is_degenerate());
|
||||||
CGAL_precondition(!tr.is_degenerate());
|
CGAL_precondition(!tr.is_degenerate());
|
||||||
|
|
||||||
typedef typename K::Point_3 Point_3;
|
typedef typename K::Point_3 Point_3;
|
||||||
typedef typename K::Segment_3 Segment_3;
|
typedef typename K::Plane_3 Plane_3;
|
||||||
typedef typename K::Triangle_3 Triangle_3;
|
|
||||||
typedef std::vector<Point_3> Poly;
|
|
||||||
|
|
||||||
typename K::Bounded_side_3 bounded_side = k.bounded_side_3_object();
|
typename K::Construct_plane_3 plane = k.construct_plane_3_object();
|
||||||
typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object();
|
typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object();
|
||||||
typename K::Construct_triangle_3 triangle = k.construct_triangle_3_object();
|
typename K::Construct_triangle_3 triangle = k.construct_triangle_3_object();
|
||||||
|
typename K::Construct_segment_3 segment = k.construct_segment_3_object();
|
||||||
|
typename K::Construct_line_3 line = k.construct_line_3_object();
|
||||||
|
typename K::Oriented_side_3 oriented_side = k.oriented_side_3_object();
|
||||||
|
typename K::Orientation_3 orientation = k.orientation_3_object();
|
||||||
|
|
||||||
std::vector<Bounded_side> vertex_sides(3);
|
|
||||||
|
|
||||||
std::vector<Point_3> points;
|
std::vector<Point_3> res = { vertex(tr,0), vertex(tr,1), vertex(tr,2) };
|
||||||
int inside_points = 0;
|
std::vector<std::bitset<4>> supporting_planes(3); // bitset used to indicate when a point is on a plane
|
||||||
int strictly_inside_points = 0;
|
|
||||||
|
|
||||||
for(int i=0; i<3; ++i)
|
// iteratively clip `tr` with the halfspaces whose intersection form `tet`
|
||||||
|
static constexpr std::array<int8_t, 12> vids = { 1,2,3, 0,3,2, 0,1,3, 1,0,2 };
|
||||||
|
const bool tet_ori_positive = (orientation(tet)==POSITIVE);
|
||||||
|
for (int pid=0; pid<4; ++pid)
|
||||||
{
|
{
|
||||||
vertex_sides[i] = bounded_side(tet, vertex(tr, i));
|
Plane_3 pl = tet_ori_positive
|
||||||
|
? plane(vertex(tet, vids[pid*3]), vertex(tet, vids[pid*3+2]),vertex(tet, vids[pid*3+1]))
|
||||||
|
: plane(vertex(tet, vids[pid*3]), vertex(tet, vids[pid*3+1]),vertex(tet, vids[pid*3+2]));
|
||||||
|
CGAL_assertion(oriented_side(pl, vertex(tet,pid))==ON_POSITIVE_SIDE);
|
||||||
|
|
||||||
if(vertex_sides[i] != ON_UNBOUNDED_SIDE)
|
std::vector<Point_3> current;
|
||||||
++inside_points;
|
std::vector<std::bitset<4>> current_sp;
|
||||||
|
std::vector<Oriented_side> orientations(res.size());
|
||||||
if(vertex_sides[i] == ON_BOUNDED_SIDE)
|
for (std::size_t i=0; i<res.size(); ++i)
|
||||||
{
|
{
|
||||||
++strictly_inside_points;
|
orientations[i]=oriented_side(pl, res[i]);
|
||||||
points.push_back(vertex(tr, i));
|
if (orientations[i]==ON_ORIENTED_BOUNDARY)
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
switch(inside_points)
|
|
||||||
{
|
|
||||||
case 0:
|
|
||||||
{
|
|
||||||
Inter_type intersections[4];
|
|
||||||
std::list<Segment_3> segments;
|
|
||||||
std::vector<std::size_t> seg_ids;
|
|
||||||
for(std::size_t i = 0; i < 4; ++i)
|
|
||||||
{
|
{
|
||||||
const Triangle_3 face = triangle(vertex(tet, (i+1)%4),
|
supporting_planes[i].set(pid);
|
||||||
vertex(tet, (i+2)%4),
|
//workaround for kernels with inexact constructions
|
||||||
vertex(tet, (i+3)%4));
|
//--
|
||||||
intersections[i] = intersection(tr, face, k);
|
if (supporting_planes[i].count()==3)
|
||||||
if(intersections[i])
|
|
||||||
{
|
{
|
||||||
// a face is inside the input tr
|
for (int b=0; i<4; ++b)
|
||||||
if(const Triangle_3* t = boost::get<Triangle_3>(&*intersections[i]))
|
|
||||||
{
|
{
|
||||||
Triangle_3 res = *t;
|
if (!supporting_planes[i].test(b))
|
||||||
return result_type(std::forward<Triangle_3>(res));
|
{
|
||||||
}
|
res[i] = vertex(tet, b);
|
||||||
else if(const Segment_3* s = boost::get<Segment_3>(&*intersections[i]))
|
break;
|
||||||
{
|
}
|
||||||
// get segs and pts to construct poly
|
|
||||||
segments.push_back(*s);
|
|
||||||
seg_ids.push_back(i);
|
|
||||||
}
|
|
||||||
else if(const Point_3* p = boost::get<Point_3>(&*intersections[i]))
|
|
||||||
{
|
|
||||||
points.push_back(*p);
|
|
||||||
}
|
|
||||||
else if(const Poly* p = boost::get<Poly>(&*intersections[i]))
|
|
||||||
{
|
|
||||||
// the input triangle is in the supporting plane of a tet face, return the poly.
|
|
||||||
Poly res = *p;
|
|
||||||
return result_type(std::forward<Poly>(res));
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
//--
|
||||||
|
|
||||||
if(segments.size() > 1)
|
|
||||||
filter_segments(segments);
|
|
||||||
|
|
||||||
// no segments and no inside points, there can still be an intersecting (tet vertex on
|
|
||||||
// an edge|face of the triangle)
|
|
||||||
if(segments.empty())
|
|
||||||
{
|
|
||||||
if(points.empty())
|
|
||||||
return result_type();
|
|
||||||
|
|
||||||
return result_type(std::forward<Point_3>(points.front()));
|
|
||||||
}
|
|
||||||
else if(segments.size() == 1)
|
|
||||||
{
|
|
||||||
// adjacency to an edge, return resulting segment.
|
|
||||||
return result_type(segments.front());
|
|
||||||
}
|
|
||||||
else if(segments.size() > 1)
|
|
||||||
{
|
|
||||||
return build_intersection(tet, tr, points, segments, k);
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
break;
|
|
||||||
case 1:
|
for (std::size_t i=0; i<res.size(); ++i)
|
||||||
case 2: // 1 or 2 inside points
|
|
||||||
{
|
{
|
||||||
Inter_type intersections[4];
|
const bool test_segment = i!=1 || res.size()!=2;
|
||||||
std::list<Segment_3> segments;
|
std::size_t j = (i+1)%res.size();
|
||||||
for(std::size_t i = 0; i < 4; ++i)
|
switch(orientations[j])
|
||||||
{
|
{
|
||||||
const Triangle_3 face = triangle(vertex(tet, (i+1)%4),
|
case ON_POSITIVE_SIDE:
|
||||||
vertex(tet, (i+2)%4),
|
if (test_segment && orientations[i]==ON_NEGATIVE_SIDE)
|
||||||
vertex(tet, (i+3)%4));
|
|
||||||
intersections[i] = intersection(tr, face, k);
|
|
||||||
if(intersections[i])
|
|
||||||
{
|
|
||||||
if(const Triangle_3* t = boost::get<Triangle_3>(&*intersections[i]))
|
|
||||||
{
|
{
|
||||||
Triangle_3 res = *t;
|
current_sp.push_back(supporting_planes[i] & supporting_planes[j]);
|
||||||
return result_type(std::forward<Triangle_3>(res));
|
current_sp.back().set(pid);
|
||||||
}
|
if (current_sp.back().count()==3)
|
||||||
else if(const Segment_3* s = boost::get<Segment_3>(&*intersections[i]))
|
|
||||||
{
|
|
||||||
segments.push_back(*s);
|
|
||||||
}
|
|
||||||
else if(const Point_3* p = boost::get<Point_3>(&*intersections[i]))
|
|
||||||
{
|
|
||||||
points.push_back(*p);
|
|
||||||
}
|
|
||||||
else if(const Poly* p = boost::get<Poly>(&*intersections[i]))
|
|
||||||
{
|
|
||||||
// the input is in a supporting plane of a face
|
|
||||||
Poly res = *p;
|
|
||||||
return result_type(std::forward<Poly>(res));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
if(segments.size() > 1)
|
|
||||||
filter_segments(segments);
|
|
||||||
|
|
||||||
switch(segments.size())
|
|
||||||
{
|
|
||||||
case 0:
|
|
||||||
{
|
|
||||||
// there can only be one point of contact, otherwise by convexity
|
|
||||||
// there would be a full segment on a face (interior segment isn't possible either
|
|
||||||
// because there are at most 2 inside points and an interior segment would also
|
|
||||||
// yield at least a segment on the boundary)
|
|
||||||
return result_type(std::forward<Point_3>(points.front()));
|
|
||||||
}
|
|
||||||
case 1: // 1 segment
|
|
||||||
{
|
|
||||||
const Segment_3& s = segments.front();
|
|
||||||
|
|
||||||
if(strictly_inside_points == 1)
|
|
||||||
{
|
|
||||||
// Doesn't matter whether there is another (non-strictly) inside point: if there is,
|
|
||||||
// it is an extremity of the segment
|
|
||||||
|
|
||||||
const int str_inside_pt_pos =
|
|
||||||
int(std::find(vertex_sides.begin(), vertex_sides.end(), ON_BOUNDED_SIDE) - vertex_sides.begin());
|
|
||||||
CGAL_assertion(str_inside_pt_pos >= 0 && str_inside_pt_pos < 3);
|
|
||||||
|
|
||||||
Triangle_3 res_tr = triangle(vertex(tr, str_inside_pt_pos), s.source(), s.target());
|
|
||||||
return result_type(std::forward<Triangle_3>(res_tr));
|
|
||||||
}
|
|
||||||
else if(strictly_inside_points == 2)
|
|
||||||
{
|
|
||||||
CGAL_assertion(inside_points == 2); // can't be 3 since we're in the 1&2 switch
|
|
||||||
|
|
||||||
Poly res(4);
|
|
||||||
|
|
||||||
// Grab the 2 strictly inside points
|
|
||||||
int id = 0;
|
|
||||||
for(int i=0; i<3; ++i)
|
|
||||||
if(vertex_sides[i] == ON_BOUNDED_SIDE)
|
|
||||||
res[id++] = vertex(tr, i);
|
|
||||||
|
|
||||||
CGAL_assertion(id == 2);
|
|
||||||
|
|
||||||
if((res[0] - res[1]) * (s.source() - s.target()) > 0)
|
|
||||||
{
|
{
|
||||||
res[2] = s.target();
|
for (int b=0; i<4; ++b)
|
||||||
res[3] = s.source();
|
if (!current_sp.back().test(b))
|
||||||
|
{
|
||||||
|
current.push_back(vertex(tet, b));
|
||||||
|
break;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
current.push_back(*CGAL::Intersections::internal::intersection_point(pl, line(res[i], res[j]), k));
|
||||||
res[3] = s.target();
|
|
||||||
res[2] = s.source();
|
|
||||||
}
|
|
||||||
|
|
||||||
return result_type(std::forward<Poly>(res));
|
|
||||||
}
|
}
|
||||||
else if(inside_points == 1) // 1 point on the boundary
|
current.push_back(res[j]);
|
||||||
|
current_sp.push_back(supporting_planes[j]);
|
||||||
|
break;
|
||||||
|
case ON_NEGATIVE_SIDE:
|
||||||
|
if (test_segment && orientations[i]==ON_POSITIVE_SIDE)
|
||||||
{
|
{
|
||||||
CGAL_assertion(strictly_inside_points == 0);
|
current_sp.push_back(supporting_planes[i] & supporting_planes[j]);
|
||||||
|
current_sp.back().set(pid);
|
||||||
// Grab the inside point
|
if (current_sp.back().count()==3)
|
||||||
const int boundary_pt_pos =
|
|
||||||
int(std::find(vertex_sides.begin(), vertex_sides.end(), ON_BOUNDARY) - vertex_sides.begin());
|
|
||||||
CGAL_assertion(boundary_pt_pos >= 0 && boundary_pt_pos < 3);
|
|
||||||
|
|
||||||
const Point_3& boundary_pt = vertex(tr, boundary_pt_pos);
|
|
||||||
if(boundary_pt == s.source() || boundary_pt == s.target())
|
|
||||||
{
|
{
|
||||||
return result_type(s);
|
for (int b=0; i<4; ++b)
|
||||||
|
if (!current_sp.back().test(b))
|
||||||
|
{
|
||||||
|
current.push_back(vertex(tet, b));
|
||||||
|
break;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
current.push_back(*CGAL::Intersections::internal::intersection_point(pl, line(res[i], res[j]), k));
|
||||||
Triangle_3 res_tr = triangle(boundary_pt, s.source(), s.target());
|
|
||||||
return result_type(std::forward<Triangle_3>(res_tr));
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
else // 2 points on the boundary
|
break;
|
||||||
{
|
|
||||||
CGAL_assertion(inside_points == 2 && strictly_inside_points == 0);
|
|
||||||
|
|
||||||
// 2 boundary points and 1 segment, have to distinguish between cases
|
|
||||||
// depending on if the extremities of the segment are triangle extremities
|
|
||||||
|
|
||||||
std::array<int, 2> boundary_pts;
|
|
||||||
std::array<bool, 2> is_boundary_point_an_extremity;
|
|
||||||
|
|
||||||
// Grab the inside points
|
|
||||||
std::size_t id = 0;
|
|
||||||
for(int i=0; i<3; ++i)
|
|
||||||
{
|
|
||||||
if(vertex_sides[i] == ON_BOUNDARY)
|
|
||||||
{
|
|
||||||
boundary_pts[id] = i;
|
|
||||||
|
|
||||||
if(vertex(tr, i) == s.source())
|
|
||||||
is_boundary_point_an_extremity[id] = true;
|
|
||||||
else if(vertex(tr, i) == s.target())
|
|
||||||
is_boundary_point_an_extremity[id] = true;
|
|
||||||
else
|
|
||||||
is_boundary_point_an_extremity[id] = false;
|
|
||||||
|
|
||||||
++id;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
CGAL_assertion(id == 2);
|
|
||||||
|
|
||||||
if(is_boundary_point_an_extremity[0])
|
|
||||||
{
|
|
||||||
if(is_boundary_point_an_extremity[1])
|
|
||||||
{
|
|
||||||
// the segment is composed of the two boundary points
|
|
||||||
return result_type(s);
|
|
||||||
}
|
|
||||||
else // only boundary_pts[0] is an extremity
|
|
||||||
{
|
|
||||||
Triangle_3 res_tr = triangle(s.source(), s.target(), vertex(tr, boundary_pts[1]));
|
|
||||||
return result_type(std::forward<Triangle_3>(res_tr));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
else // boundary_pts[0] is not an extremity
|
|
||||||
{
|
|
||||||
if(is_boundary_point_an_extremity[1]) // only boundary_pts[1] is an extremity
|
|
||||||
{
|
|
||||||
Triangle_3 res_tr = triangle(s.source(), s.target(), vertex(tr, boundary_pts[0]));
|
|
||||||
return result_type(std::forward<Triangle_3>(res_tr));
|
|
||||||
}
|
|
||||||
else // neither boundary points are extremities
|
|
||||||
{
|
|
||||||
Poly res(4);
|
|
||||||
res[0] = vertex(tr, boundary_pts[0]);
|
|
||||||
res[1] = vertex(tr, boundary_pts[1]);
|
|
||||||
|
|
||||||
if((res[0] - res[1]) * (s.source() - s.target()) > 0)
|
|
||||||
{
|
|
||||||
res[2] = s.target();
|
|
||||||
res[3] = s.source();
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
res[3] = s.target();
|
|
||||||
res[2] = s.source();
|
|
||||||
}
|
|
||||||
|
|
||||||
return result_type(std::forward<Poly>(res));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
CGAL_assertion(false);
|
|
||||||
}
|
|
||||||
break;
|
|
||||||
// 2 or 3 segments (and 1 or 2 inside points)
|
|
||||||
case 2:
|
|
||||||
case 3:
|
|
||||||
case 4:
|
|
||||||
{
|
|
||||||
// @todo do that for a single segment too?
|
|
||||||
return build_intersection(tet, tr, points, segments, k);
|
|
||||||
}
|
|
||||||
break;
|
|
||||||
default:
|
default:
|
||||||
// can't have more than 4 segments (1 per tet face)
|
{
|
||||||
CGAL_assertion(false);
|
CGAL_assertion(supporting_planes[j].test(pid));
|
||||||
break;
|
current.push_back(res[j]);
|
||||||
|
current_sp.push_back(supporting_planes[j]);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
break;
|
res.swap(current);
|
||||||
|
supporting_planes.swap(current_sp);
|
||||||
|
|
||||||
case 3:
|
if (res.empty())
|
||||||
{
|
return boost::none;
|
||||||
// the input triangle is entirely contained within the tetrahedron
|
|
||||||
return result_type(tr);
|
|
||||||
}
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
CGAL_assertion(false); // never happens (only 3 pts in a tr)
|
|
||||||
break;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
return result_type();
|
switch(res.size())
|
||||||
|
{
|
||||||
|
case 1:
|
||||||
|
return result_type(res[0]);
|
||||||
|
case 2:
|
||||||
|
return result_type(segment(res[0], res[1]));
|
||||||
|
case 3:
|
||||||
|
return result_type(triangle(res[0], res[1], res[2]));
|
||||||
|
default:
|
||||||
|
return result_type(res);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
template <class K>
|
template <class K>
|
||||||
|
|
|
||||||
|
|
@ -143,6 +143,8 @@ public:
|
||||||
// edge shared, 3rd point outside
|
// edge shared, 3rd point outside
|
||||||
check_intersection(tet, Tr(p(0,1,0), p(1,0,0), P(0.5,0,-100)),
|
check_intersection(tet, Tr(p(0,1,0), p(1,0,0), P(0.5,0,-100)),
|
||||||
S(p(0,1,0), p(1,0,0)));
|
S(p(0,1,0), p(1,0,0)));
|
||||||
|
check_intersection(tet, Tr(P(0.75,0.25,0), p(10,10,10), P(0.25,0.75,0)),
|
||||||
|
S(P(0.75,0.25,0), P(0.25,0.75,0)));
|
||||||
|
|
||||||
// shared edge, 3rd point inside
|
// shared edge, 3rd point inside
|
||||||
check_intersection(tet, Tr(p(0,1,0), p(1,0,0), P(0.25,0.25,0.25)),
|
check_intersection(tet, Tr(p(0,1,0), p(1,0,0), P(0.25,0.25,0.25)),
|
||||||
|
|
@ -166,7 +168,7 @@ public:
|
||||||
|
|
||||||
// face inside tr
|
// face inside tr
|
||||||
check_intersection(tet, Tr(p(0,2,0), p(0,-2,0), p(0,0,2)),
|
check_intersection(tet, Tr(p(0,2,0), p(0,-2,0), p(0,0,2)),
|
||||||
Tr(p(0,0,1), p(0,0,0), p(0,1,0)));
|
Tr(p(0,1,0), p(0,0,0), p(0,0,1)));
|
||||||
|
|
||||||
// on the same plane, polygonal intersection
|
// on the same plane, polygonal intersection
|
||||||
Base::template check_intersection<Poly>(tet, Tr(p(0,-2,0), p(1,1,0), p(1,2,0)));
|
Base::template check_intersection<Poly>(tet, Tr(p(0,-2,0), p(1,1,0), p(1,2,0)));
|
||||||
|
|
@ -189,6 +191,14 @@ public:
|
||||||
// vertex on edge & triangle inside, double segment non-incident
|
// vertex on edge & triangle inside, double segment non-incident
|
||||||
Base::template check_intersection<Poly>(tet, Tr(P(0.25,0,0.25), P(-1,0.5,0.25), P(1.5,0.5,0.25)));
|
Base::template check_intersection<Poly>(tet, Tr(P(0.25,0,0.25), P(-1,0.5,0.25), P(1.5,0.5,0.25)));
|
||||||
|
|
||||||
|
// vertex on face, triangle outside & point intersection
|
||||||
|
Base::check_intersection(tet, Tr(P(-1,1,0.25), P(-1,0,0.25), P(0,0.25,0.25)),
|
||||||
|
P(0, 0.25, 0.25));
|
||||||
|
Base::check_intersection(tet, Tr(P(-1,0,0.25), P(-1,1,0.25), P(0,0.25,0.25)),
|
||||||
|
P(0, 0.25, 0.25));
|
||||||
|
Base::check_intersection(tet, Tr(P(0,0.25,0.25), P(-1,1,0.25), P(-1,0,0.25)),
|
||||||
|
P(0, 0.25, 0.25));
|
||||||
|
|
||||||
// vertex on face, triangle outside & segment intersection
|
// vertex on face, triangle outside & segment intersection
|
||||||
Base::check_intersection(tet, Tr(P(0.5,0,-0.25), P(0.5,0,0.25), P(0.5,-0.5,0)),
|
Base::check_intersection(tet, Tr(P(0.5,0,-0.25), P(0.5,0,0.25), P(0.5,-0.5,0)),
|
||||||
S(P(0.5, 0, 0.25), P(0.5, 0, 0)));
|
S(P(0.5, 0, 0.25), P(0.5, 0, 0)));
|
||||||
|
|
@ -307,12 +317,23 @@ public:
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
void issue_6777()
|
||||||
|
{
|
||||||
|
Tr tri(P(0.191630, -0.331630, -0.370000), P(-0.124185, -0.385815, -0.185000), P(-0.0700000, -0.0700000, 0.00000));
|
||||||
|
Tet tet(P(0, -1, 0), P(-1, 0, 0), P(0, 0, 0), P(0, 0, -1));
|
||||||
|
auto res = intersection(tri, tet);
|
||||||
|
assert(res != boost::none);
|
||||||
|
const std::vector<P> *vps = boost::get<std::vector<P>>(&*res);
|
||||||
|
assert(vps!=nullptr);
|
||||||
|
}
|
||||||
|
|
||||||
void run()
|
void run()
|
||||||
{
|
{
|
||||||
std::cout << "3D Tetrahedron Intersection tests\n";
|
std::cout << "3D Tetrahedron Intersection tests\n";
|
||||||
|
|
||||||
Tet_Tet();
|
Tet_Tet();
|
||||||
Tet_Tr();
|
Tet_Tr();
|
||||||
|
issue_6777();
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue