Merge pull request #7614 from sloriot/PMP-BF_autorefine

Add autorefine for triangle soup
This commit is contained in:
Sébastien Loriot 2023-12-26 12:11:04 +01:00
commit 6eff89bbc8
22 changed files with 2314 additions and 99 deletions

View File

@ -14,7 +14,7 @@ Release date: October 2023
- **Breaking change**: The usage of `boost::shared_ptr` has been replaced by `std::shared_ptr`. Packages affected are 2D Straight Line Skeleton and Shape Detection.
- **Breaking change**: The usage of `boost::optional` has been replaced by `std::optional`. Packages affected are 2D Straight Line Skeleton, 3D Fast Intersection and Distance Computation (AABB Tree), and the Kernel intersection.
- **Breaking change**: The usage of `boost::variant` has been replaced by `std::variant`. Packages affected are 2D Arrangements, and the Kernel intersection.
- **Breaking chahge**: The file CMake file `UseCGAL.cmake` has been removed from CGAL. Usages of the CMake variables `${CGAL_USE_FILE}` and `${CGAL_LIBRARIES}` must be replaced by a link to the imported target `CGAL::CGAL`, for example: `target_link_library(the_target PRIVATE CGAL::CGAL)`.
- **Breaking change**: The file CMake file `UseCGAL.cmake` has been removed from CGAL. Usages of the CMake variables `${CGAL_USE_FILE}` and `${CGAL_LIBRARIES}` must be replaced by a link to the imported target `CGAL::CGAL`, for example: `target_link_library(the_target PRIVATE CGAL::CGAL)`.
#### 2D Arrangements
@ -45,6 +45,10 @@ Release date: October 2023
- Added the option to use a variable sizing field for `CGAL::Polygon_mesh_processing::isotropic_remeshing()`,
and a sizing function based on a measure of local curvature for adaptive remeshing.
- Added the function `CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel()` that refines a polygon mesh along an isocurve.
- Added the function
`CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` that refines a soup of triangles so that no pair of triangles intersects
in their interiors. Also added, the function `autorefine()` operating directly on a triangle mesh and updating it
using the aforementioned function on a triangle soup.
### [2D Arrangements](https://doc.cgal.org/6.0/Manual/packages.html#PkgArrangementOnSurface2)
- Fixed a bug in the zone construction code applied to arrangements of geodesic arcs on a sphere,

View File

@ -14,6 +14,8 @@
#ifndef CGAL_INTERNAL_INTERSECTIONS_TRIANGLE_3_TRIANGLE_3_INTERSECTION_H
#define CGAL_INTERNAL_INTERSECTIONS_TRIANGLE_3_TRIANGLE_3_INTERSECTION_H
//#define CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
#include <CGAL/Intersection_traits_3.h>
#include <CGAL/Intersections_3/internal/Line_3_Triangle_3_intersection.h>
#include <CGAL/Intersections_3/internal/Line_3_Line_3_intersection.h>
@ -21,6 +23,9 @@
#include <CGAL/kernel_assertions.h>
#include <boost/next_prior.hpp>
#include <boost/container/flat_set.hpp>
#include <list>
#include <map>
#include <vector>
@ -30,53 +35,338 @@ namespace Intersections {
namespace internal{
template <class Kernel>
void intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3& p,
struct Point_on_triangle
{
// triangle points are not stored in this class but are expected
// to always be passed in the same order. For a triangle pqr,
// edge 0 is pq, edge 1 qr and edge 2 rp. Point 0 is p, 1 is q and 2 is r.
//
// (id, -1) point on t1
// (-1, id) point on t2
// (id1, id2) intersection of edges
std::pair<int, int> t1_t2_ids;
boost::container::flat_set<int> extra_t1; // store other ids of edges containing the point
typename Kernel::FT alpha; //
//////
static
inline
const typename Kernel::Point_3&
point_from_id(const typename Kernel::Point_3& p,
const typename Kernel::Point_3& q,
const typename Kernel::Point_3& r,
const Kernel& k,
std::list<typename Kernel::Point_3>& inter_pts)
int id)
{
switch(id)
{
case 0:
return p;
case 1:
return q;
default:
return r;
}
}
Point_on_triangle(int i1, int i2=-1, typename Kernel::FT alpha = 0.) // TODO add global zero()?
: t1_t2_ids(i1,i2)
, alpha(alpha)
{}
// orientation of the current point wrt to edge id1 (p1q1)
Orientation
orientation (const typename Kernel::Point_3& p1, // source of edge edge_id1
const typename Kernel::Point_3& q1, // target of edge edge_id1
const typename Kernel::Point_3& r1,
int edge_id1,
const typename Kernel::Point_3& p2,
const typename Kernel::Point_3& q2,
const typename Kernel::Point_3& r2,
const Kernel& k) const
{
if (t1_t2_ids.first!=-1)
{
if (t1_t2_ids.second==-1)
return (edge_id1==t1_t2_ids.first || (edge_id1+1)%3==t1_t2_ids.first) ? ZERO:POSITIVE; // it is a point on t1
// this is an intersection point
if (t1_t2_ids.first==edge_id1)
return ZERO;
if (t1_t2_ids.first==(edge_id1+1)%3)
{
if (alpha==0) return ZERO;
return alpha>=0 ? POSITIVE:NEGATIVE;
}
CGAL_assertion((t1_t2_ids.first+1)%3==edge_id1);
if (alpha==1) return ZERO;
return alpha<=1?POSITIVE:NEGATIVE;
}
else
{
//this is an input point of t2
typename Kernel::Coplanar_orientation_3 orient = k.coplanar_orientation_3_object();
const typename Kernel::Point_3& query = point_from_id(p2,q2,r2,t1_t2_ids.second);
return orient(p1,q1,r1,query);
}
}
int id1() const { return t1_t2_ids.first; }
int id2() const { return t1_t2_ids.second; }
// construct the intersection point from the info stored
typename Kernel::Point_3
point(const typename Kernel::Point_3& p1,
const typename Kernel::Point_3& q1,
const typename Kernel::Point_3& r1,
const typename Kernel::Point_3& p2,
const typename Kernel::Point_3& q2,
const typename Kernel::Point_3& r2,
const Kernel& k) const
{
if (t1_t2_ids.first==-1)
return point_from_id(p2,q2,r2,t1_t2_ids.second);
if (t1_t2_ids.second==-1)
return point_from_id(p1,q1,r1,t1_t2_ids.first);
return k.construct_barycenter_3_object()(point_from_id(p1,q1,r1,(t1_t2_ids.first+1)%3), alpha, point_from_id(p1,q1,r1,t1_t2_ids.first)) ;
}
};
// the intersection of two triangles is computed by interatively intersection t2
// with halfspaces defined by edges of t1. The following function is called
// for each each on t1 on edge of the current intersection.
// pq is such an edge and p1q1 from t1 defines the halfspace intersection
// we are currently interseted in. We return the intersection point of
// pq with p1q1
template <class Kernel>
Point_on_triangle<Kernel>
intersection(const Point_on_triangle<Kernel>& p,
const Point_on_triangle<Kernel>& q,
int edge_id_t1,
const typename Kernel::Point_3& p1,
const typename Kernel::Point_3& q1,
// const typename Kernel::Point_3& r1,
const typename Kernel::Point_3& p2,
const typename Kernel::Point_3& q2,
const typename Kernel::Point_3& r2,
const Kernel& k)
{
typedef typename std::list<typename Kernel::Point_3>::iterator Iterator;
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " calling intersection: ";
std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) -";
std::cout << " (" << q.id1() << "," << q.id2() << ",[" << q.alpha << "]) || e" << edge_id_t1;
#endif
typename Kernel::Compute_alpha_for_coplanar_triangle_intersection_3 compute_alpha
= k.compute_alpha_for_coplanar_triangle_intersection_3_object();
typedef Point_on_triangle<Kernel> Pot;
switch(p.id1())
{
case -1:
{
switch(q.id1())
{
case -1: // A: (-1, ip2) - (-1, iq2)
{
CGAL_assertion((p.id2()+1)%3 == q.id2() || (q.id2()+1)%3 == p.id2());
// CGAL_assertion(p.extra_t1.empty() && q.extra_t1.empty()); // TMP to see if it's worth implementing special case
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " -- case 1\n";
#endif
typename Kernel::FT alpha = compute_alpha(p1, q1,
Pot::point_from_id(p2, q2, r2, p.id2()),
Pot::point_from_id(p2, q2, r2, q.id2()));
int id2 = (p.id2()+1)%3 == q.id2() ? p.id2() : q.id2();
return Point_on_triangle<Kernel>(edge_id_t1, id2, alpha); // intersection with an original edge of t2
}
default:
if (q.id2()!=-1) // B: (-1, ip2) - (iq1, iq2)
{
if (p.id2() == q.id2() || p.id2() == (q.id2()+1)%3)
{
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " -- case 2\n";
#endif
// points are on the same edge of t2 --> we shorten an already cut edge
typename Kernel::FT alpha = compute_alpha(p1, q1,
Pot::point_from_id(p2, q2, r2, q.id2()),
Pot::point_from_id(p2, q2, r2, (q.id2()+1)%3));
return Point_on_triangle<Kernel>(edge_id_t1, q.id2(), alpha);
}
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " -- case 3\n";
#endif
// point of t1: look for an edge of t1 containing both points
CGAL_assertion( p.extra_t1.count(q.id1())!=0 || p.extra_t1.count(3-q.id1()-edge_id_t1)!=0 );
int eid1 = p.extra_t1.count(q.id1())!=0 ? q.id1() : 3-q.id1()-edge_id_t1;
return Point_on_triangle<Kernel>((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1
}
// C: (-1, ip2) - (iq1, -1)
//vertex of t1, special case t1 edge passed thru a vertex of t2
CGAL_assertion(edge_id_t1 == 2);
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " -- case 4\n";
#endif
CGAL_assertion(q.id1()==1);
CGAL_assertion(!p.extra_t1.empty());
return Point_on_triangle<Kernel>(p.extra_t1.count(0)==1?0:2,-1);
}
}
default:
{
switch(p.id2())
{
case -1:
{
switch(q.id1())
{
case -1: // G: (ip1, -1) - (-1, iq2)
//vertex of t1, special case t1 edge passed thru a vertex of t2
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " -- case 5\n";
#endif
CGAL_assertion(edge_id_t1 == 2);
CGAL_assertion(p.id1()==1);
CGAL_assertion(!q.extra_t1.empty());
return Point_on_triangle<Kernel>(q.extra_t1.count(0)==1?0:2,-1);
default:
{
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " -- case 6\n";
#endif
CGAL_assertion(q.id2()!=-1); // I: (ip1, -1) - (iq2, -1)
//H: (ip1,-1), (iq1, iq2)
CGAL_assertion(edge_id_t1==2);
// p and q are on the same edge of t1
CGAL_assertion(p.id1()==q.id1() || p.id1()==(q.id1()+1)%3);
return Point_on_triangle<Kernel>((q.id1()+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3 , -1);
}
}
}
default:
{
switch(q.id1())
{
case -1: // D: (ip1, ip2) - (-1, iq2)
{
if (q.id2() == p.id2() || q.id2() == (p.id2()+1)%3)
{
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " -- case 7\n";
#endif
// points are on the same edge of t2 --> we shorten an already cut edge
typename Kernel::FT alpha = compute_alpha(p1, q1,
Pot::point_from_id(p2, q2, r2, p.id2()),
Pot::point_from_id(p2, q2, r2, (p.id2()+1)%3));
return Point_on_triangle<Kernel>(edge_id_t1, p.id2(), alpha);
}
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " -- case 8\n";
#endif
// point of t1
//std::cout << "q.extra_t1: "; for(int qet1 : q.extra_t1) std::cout << " " << qet1; std::cout << "\n";
CGAL_assertion( q.extra_t1.count(p.id1())!=0 || q.extra_t1.count(3-p.id1()-edge_id_t1)!=0 );
int eid1 = q.extra_t1.count(p.id1())!=0 ? p.id1() : 3-p.id1()-edge_id_t1;
return Point_on_triangle<Kernel>((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1
}
default:
{
switch(q.id2())
{
case -1: // F: (ip1, ip2) - (iq1, -1)
{
// p and q are on the same edge of t1
CGAL_assertion(q.id1()==p.id1() || q.id1()==(p.id1()+1)%3);
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " -- case 9\n";
#endif
return Point_on_triangle<Kernel>((p.id1()+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3 , -1);
}
default: // E: (ip1, ip2) - (iq1, iq2)
{
if (p.id2()==q.id2())
{
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " -- case 10\n";
#endif
typename Kernel::FT alpha = compute_alpha(p1, q1,
Pot::point_from_id(p2, q2, r2, q.id2()),
Pot::point_from_id(p2, q2, r2, (q.id2()+1)%3));
return Point_on_triangle<Kernel>(edge_id_t1, q.id2(), alpha);
}
// we are intersecting an edge of t1
CGAL_assertion(p.id1()==q.id1() || edge_id_t1==2);
int eid1 = p.id1()==q.id1() ? p.id1() : 1;
return Point_on_triangle<Kernel>((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1
}
}
}
}
}
}
}
}
}
template <class Kernel>
void intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3& p1,
const typename Kernel::Point_3& q1,
const typename Kernel::Point_3& r1,
int edge_id,
const typename Kernel::Point_3& p2,
const typename Kernel::Point_3& q2,
const typename Kernel::Point_3& r2,
const Kernel& k,
std::list<Point_on_triangle<Kernel>>& inter_pts)
{
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " cutoff using e" << edge_id << ": "
<< to_double(p1.x()) << " " << to_double(p1.y()) << " " << to_double(p1.z()) << " "
<< to_double(q1.x()) << " " << to_double(q1.y()) << " " << to_double(q1.z()) << "\n";
#endif
typedef typename std::list<Point_on_triangle<Kernel>>::iterator Iterator;
if(inter_pts.empty())
return;
typename Kernel::Coplanar_orientation_3 orient = k.coplanar_orientation_3_object();
typename Kernel::Construct_line_3 line = k.construct_line_3_object();
//orient(p1,q1,r1,r1) is POSITIVE
std::map<const Point_on_triangle<Kernel>*,Orientation> orientations; // TODO skip map
for (Point_on_triangle<Kernel>& pot : inter_pts)
{
orientations[ &pot ]=pot.orientation(p1,q1,r1,edge_id,p2,q2,r2,k);
if (pot.id1()==-1 && orientations[ &pot ]==COLLINEAR)
pot.extra_t1.insert(edge_id);
}
//orient(p,q,r,r) is POSITIVE
std::map<const typename Kernel::Point_3*,Orientation> orientations;
for (Iterator it=inter_pts.begin();it!=inter_pts.end();++it)
orientations[ &(*it) ]=orient(p,q,r,*it);
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " Orientations:";
for (const Point_on_triangle<Kernel>& pot : inter_pts)
std::cout << " " << orientations[ &pot ];
std::cout << "\n";
#endif
CGAL_kernel_assertion_code(int pt_added = 0);
CGAL_kernel_assertion_code(int pt_added = 0;)
Iterator prev = std::prev(inter_pts.end());
const typename Kernel::Point_3* prev = &(*std::prev(inter_pts.end()));
Iterator stop = inter_pts.size() > 2 ? inter_pts.end() : std::prev(inter_pts.end());
for(Iterator it=inter_pts.begin(); it!=stop; ++it)
{
const typename Kernel::Point_3& curr = *it;
Orientation or_prev = orientations[prev],
or_curr = orientations[&curr];
Orientation or_prev = orientations[&(*prev)],
or_curr = orientations[&(*it)];
if((or_prev == POSITIVE && or_curr == NEGATIVE) ||
(or_prev == NEGATIVE && or_curr == POSITIVE))
{
typename Intersection_traits<Kernel, typename Kernel::Line_3, typename Kernel::Line_3>::result_type
obj = intersection(line(p,q), line(*prev,curr), k);
Point_on_triangle<Kernel> new_pt = intersection(*prev, *it, edge_id, p1, q1, p2, q2, r2, k);
// assert "not empty"
CGAL_kernel_assertion(bool(obj));
const typename Kernel::Point_3* inter = intersect_get<typename Kernel::Point_3>(obj);
CGAL_kernel_assertion(inter != nullptr);
prev = &(*inter_pts.insert(it,*inter));
orientations[prev] = COLLINEAR;
CGAL_kernel_assertion_code(++pt_added;)
prev = inter_pts.insert(it,new_pt);
orientations[&(*prev)] = COLLINEAR;
CGAL_kernel_assertion_code(++pt_added);
}
prev = &(*it);
prev = it;
}
CGAL_kernel_assertion(pt_added<3);
@ -96,35 +386,77 @@ intersection_coplanar_triangles(const typename K::Triangle_3& t1,
const typename K::Triangle_3& t2,
const K& k)
{
const typename K::Point_3& p = t1.vertex(0),
q = t1.vertex(1),
r = t1.vertex(2);
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
auto to_string = [](const typename K::Triangle_3& t)
{
std::stringstream sstr;
sstr << "4 "
<< to_double(t[0].x()) << " " << to_double(t[0].y()) << " " << to_double(t[0].z()) << " "
<< to_double(t[1].x()) << " " << to_double(t[1].y()) << " " << to_double(t[1].z()) << " "
<< to_double(t[2].x()) << " " << to_double(t[2].y()) << " " << to_double(t[2].z()) << " "
<< to_double(t[0].x()) << " " << to_double(t[0].y()) << " " << to_double(t[0].z()) << "\n";
return sstr.str();
};
std::list<typename K::Point_3> inter_pts;
inter_pts.push_back(t2.vertex(0));
inter_pts.push_back(t2.vertex(1));
inter_pts.push_back(t2.vertex(2));
std::cout << "intersection_coplanar_triangles\n";
std::ofstream("/tmp/t1.polylines.txt") << std::setprecision(17) << to_string(t1) << "\n";
std::ofstream("/tmp/t2.polylines.txt") << std::setprecision(17) << to_string(t2) << "\n";
#endif
const typename K::Point_3& p1 = t1.vertex(0),
q1 = t1.vertex(1),
r1 = t1.vertex(2);
const typename K::Point_3& p2 = t2.vertex(0),
q2 = t2.vertex(1),
r2 = t2.vertex(2);
std::list<Point_on_triangle<K>> inter_pts;
inter_pts.push_back(Point_on_triangle<K>(-1,0));
inter_pts.push_back(Point_on_triangle<K>(-1,1));
inter_pts.push_back(Point_on_triangle<K>(-1,2));
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
auto print_points = [&]()
{
for(auto p : inter_pts) std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) "; std::cout <<"\n";
};
std::cout << " ipts size: " << inter_pts.size() << "\n";
print_points();
#endif
//intersect t2 with the three half planes which intersection defines t1
intersection_coplanar_triangles_cutoff(p,q,r,k,inter_pts); //line pq
intersection_coplanar_triangles_cutoff(q,r,p,k,inter_pts); //line qr
intersection_coplanar_triangles_cutoff(r,p,q,k,inter_pts); //line rp
intersection_coplanar_triangles_cutoff(p1,q1,r1,0,p2,q2,r2,k,inter_pts); //line pq
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " ipts size: " << inter_pts.size() << "\n";
print_points();
#endif
intersection_coplanar_triangles_cutoff(q1,r1,p1,1,p2,q2,r2,k,inter_pts); //line qr
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " ipts size: " << inter_pts.size() << "\n";
print_points();
#endif
intersection_coplanar_triangles_cutoff(r1,p1,q1,2,p2,q2,r2,k,inter_pts); //line rp
#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
std::cout << " ipts size: " << inter_pts.size() << "\n";
print_points();
#endif
auto point = [&](const Point_on_triangle<K>& pot){ return pot.point(p1,q1,r1,p2,q2,r2,k); };
switch(inter_pts.size())
{
case 0:
return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>();
case 1:
return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(*inter_pts.begin());
return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(point(*inter_pts.begin()));
case 2:
return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(
k.construct_segment_3_object()(*inter_pts.begin(), *std::next(inter_pts.begin())) );
k.construct_segment_3_object()(point(*inter_pts.begin()), point(*std::next(inter_pts.begin()))) );
case 3:
return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(
k.construct_triangle_3_object()(*inter_pts.begin(), *std::next(inter_pts.begin()), *std::prev(inter_pts.end())) );
k.construct_triangle_3_object()(point(*inter_pts.begin()), point(*std::next(inter_pts.begin())), point(*std::prev(inter_pts.end()))) );
default:
return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(
std::vector<typename K::Point_3>(inter_pts.begin(),inter_pts.end()));
std::vector<typename K::Point_3>(boost::make_transform_iterator(inter_pts.begin(), point),
boost::make_transform_iterator(inter_pts.end(), point)));
}
}

View File

@ -96,6 +96,13 @@ void test_coplanar_triangles(){
assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
obj=CGAL::intersection(t2,t1);
assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
// TK10 case C'
t1=Triangle(Point(88.7921, 89.0007, 1.25), Point(88.1912, 88.3997, 1.25), Point(89.8224, 90.031, 1.25));
t2=Triangle(Point(88.0497, 88.2583, 1.25), Point(82.9292, 81.8747, 1.25), Point(91.1726, 91.3812, 1.25));
obj=CGAL::intersection(t1,t2);
assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
obj=CGAL::intersection(t2,t1);
assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
//Intersection is a point
//edges are collinear, one vertex in common
t1=Triangle( Point(0,0,0),Point(0,1,0),Point(1,0,0) );
@ -153,6 +160,13 @@ void test_coplanar_triangles(){
assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
obj=CGAL::intersection(t2,t1);
assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
// TK10 case D
t1=Triangle(Point(-34.893700000000003, -16.0351, 3.1334899999999998e-12), Point(-34.893700000000003, -18.5351, 3.1334899999999998e-12), Point(-42.393700000000003, -16.0351, 3.1334899999999998e-12));
t2=Triangle(Point(-34.893700000000003, -32.0351, 3.1334899999999998e-12), Point(-34.893700000000003, -9.7851400000000002, 3.1334899999999998e-12), Point(-31.643699999999999, -17.201799999999999, 3.1334899999999998e-12));
obj=CGAL::intersection(t1,t2);
assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
obj=CGAL::intersection(t2,t1);
assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
//Intersection is a polygon
//David's star
t1=Triangle( Point(0,0,0),Point(1,0,0),Point(0.5,1.5,0) );
@ -181,6 +195,51 @@ void test_coplanar_triangles(){
obj=CGAL::intersection(t2,t1);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
// TK10 case A
t1=Triangle(Point(3.74861, 12.4822, 14.0112), Point(5.40582, 12.4822, 15.6895), Point(5.37748, 12.4822, 15.7206));
t2=Triangle(Point(5.49972, 12.4822, 13.491), Point(5.27627, 12.4822, 15.8106), Point(5.32119, 12.4822, 15.8126));
obj=CGAL::intersection(t1,t2);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
obj=CGAL::intersection(t2,t1);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
// TK10 case C
t1=Triangle(Point(5, -94.6659, 3.85175), Point(5, -94.5682, 3.08638), Point(5, -94.8182, 3.08638));
t2=Triangle(Point(5, -94.4317, 3.76399), Point(5, -97.6182, 3.08638), Point(5, -94.5659, 2.99682));
obj=CGAL::intersection(t1,t2);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
obj=CGAL::intersection(t2,t1);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
// TK10 case E
t1=Triangle(Point(-955.858, -45.032, -0.016), Point(-955.856, -45.032, -0.004), Point(-955.856, -45.032, -0.002));
t2=Triangle(Point(-955.856, -45.032, 0.006), Point(-955.854, -45.032, -0.002), Point(-955.876, -45.032, -0.034));
obj=CGAL::intersection(t1,t2);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
obj=CGAL::intersection(t2,t1);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
// TK10 case F
t1=Triangle(Point(141.172, 20.576, 155.764), Point(141.172, 20.588, 155.766), Point(141.172, 20.59, 155.766));
t2=Triangle(Point(141.172, 20.602, 155.768), Point(141.172, 20.594, 155.766), Point(141.172, 20.574, 155.764));
obj=CGAL::intersection(t1,t2);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
obj=CGAL::intersection(t2,t1);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
// TK10 case D
t1=Triangle(Point(152.864, 126.324, 0.950001), Point(152.77, 126.483, 0.950001), Point(153.072, 125.973, 0.950001));
t2=Triangle(Point(153.322, 125.551, 0.950001), Point(152.218, 127.415, 0.950001), Point(153.66, 124.768, 0.950001));
obj=CGAL::intersection(t1,t2);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
obj=CGAL::intersection(t2,t1);
assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
//Intersection is empty
t1=Triangle( Point(0,0,0),Point(0,1,0),Point(1,0,0) );
t2=Triangle( Point(-0.1,-0.1,0),Point(-0.1,-0.9,0),Point(-1,-0.1,0) );

View File

@ -2224,6 +2224,107 @@ namespace CommonKernelFunctors {
}
};
template <typename K>
class Construct_planes_intersection_point_3
{
typedef typename K::Plane_3 Plane;
typedef typename K::Point_3 Point;
typename K::Construct_plane_3 construct_plane;
public:
typedef Point result_type;
Point
operator()(const Point& p1, const Point& q1, const Point& r1,
const Point& p2, const Point& q2, const Point& r2,
const Point& p3, const Point& q3, const Point& r3) const
{
Plane plane1 = construct_plane(p1, q1, r1);
Plane plane2 = construct_plane(p2, q2, r2);
Plane plane3 = construct_plane(p3, q3, r3);
const auto res = typename K::Intersect_3()(plane1, plane2, plane3);
CGAL_assertion(res!=std::nullopt);
const Point* e_pt = std::get_if<Point>(&(*res));
CGAL_assertion(e_pt!=nullptr);
return *e_pt;
}
Point
operator()(const Plane& plane1, const Plane& plane2, const Plane& plane3) const
{
const auto res = typename K::Intersect_3()(plane1, plane2, plane3);
CGAL_assertion(res!=std::nullopt);
const Point* e_pt = std::get_if<Point>(&(*res));
CGAL_assertion(e_pt!=nullptr);
return *e_pt;
}
};
template <typename K>
class Construct_coplanar_segments_intersection_point_3
{
typedef typename K::Segment_3 Segment;
typedef typename K::Point_3 Point;
typename K::Construct_segment_3 construct_segment;
public:
typedef Point result_type;
Point
operator()(const Point& p1, const Point& q1,
const Point& p2, const Point& q2) const
{
Segment s1 = construct_segment(p1, q1);
Segment s2 = construct_segment(p2, q2);
const auto res = typename K::Intersect_3()(s1, s2);
CGAL_assertion(res!=std::nullopt);
const Point* e_pt = std::get_if<Point>(&(*res));
CGAL_assertion(e_pt!=nullptr);
return *e_pt;
}
Point
operator()(const Segment& s1, const Segment& s2) const
{
const auto res = typename K::Intersect_3()(s1, s2);
CGAL_assertion(res!=std::nullopt);
const Point* e_pt = std::get_if<Point>(&(*res));
CGAL_assertion(e_pt!=nullptr);
return *e_pt;
}
};
template <typename K>
class Compute_alpha_for_coplanar_triangle_intersection_3
{
typedef typename K::Point_3 Point_3;
typedef typename K::Vector_3 Vector_3;
public:
typedef typename K::FT result_type;
result_type
operator()(const Point_3& p1, const Point_3& p2, // segment 1
const Point_3& p3, const Point_3& p4) const // segment 2
{
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();
const Vector_3 v1 = vector(p1, p2);
const Vector_3 v2 = vector(p3, p4);
CGAL_assertion(K().coplanar_3_object()(p1,p2,p3,p4));
const Vector_3 v3 = vector(p1, p3);
const Vector_3 v3v2 = cross_product(v3,v2);
const Vector_3 v1v2 = cross_product(v1,v2);
const typename K::FT sl = K().compute_squared_length_3_object()(v1v2);
CGAL_assertion(!certainly(is_zero(sl)));
const typename K::FT t = ((v3v2.x()*v1v2.x()) + (v3v2.y()*v1v2.y()) + (v3v2.z()*v1v2.z())) / sl;
return t; // p1 + (p2-p1) * t
}
};
template <typename K>
class Construct_point_on_2
{

View File

@ -398,6 +398,12 @@ CGAL_Kernel_cons(Construct_plane_3,
construct_plane_3_object)
CGAL_Kernel_cons(Construct_plane_line_intersection_point_3,
construct_plane_line_intersection_point_3_object)
CGAL_Kernel_cons(Construct_planes_intersection_point_3,
construct_planes_intersection_point_3_object)
CGAL_Kernel_cons(Construct_coplanar_segments_intersection_point_3,
construct_coplanar_segments_intersection_point_3_object)
CGAL_Kernel_cons(Compute_alpha_for_coplanar_triangle_intersection_3,
compute_alpha_for_coplanar_triangle_intersection_3_object)
CGAL_Kernel_cons(Construct_point_on_2,
construct_point_on_2_object)
CGAL_Kernel_cons(Construct_point_on_3,

View File

@ -0,0 +1,27 @@
/// \ingroup PkgPolygonMeshProcessingConcepts
/// \cgalConcept
///
/// The concept `PMPAutorefinementVisitor` defines the requirements for the visitor
/// used in `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` to track
/// the creation of new triangles.
///
/// \cgalRefines{CopyConstructible}
/// \cgalHasModelsBegin
/// \cgalHasModels{CGAL::Polygon_mesh_processing::Autorefinement::Default_visitor}
/// \cgalHasModelsEnd
class PMPAutorefinementVisitor{
public:
/// @name Functions called only if at least one intersection has been found
/// @{
/// called when the final number of output triangles is known, `nbt` being the total number of triangles in the output.
void number_of_output_triangles(std::size_t nbt);
/// called for triangle with no intersection, `tgt_id` is the position in the triangle container after calling
/// `autorefine_triangle_soup()`, while `src_id` was its position before calling the function.
void verbatim_triangle_copy(std::size_t tgt_id, std::size_t src_id);
/// called for each subtriangle created from a triangle with intersection, `tgt_id` is the position in the triangle container after calling
/// `autorefine_triangle_soup()` of the subtriangle, while `src_id` was the position of the original support triangle before calling the function.
void new_subtriangle(std::size_t tgt_id, std::size_t src_id);
/// @}
};

View File

@ -113,6 +113,8 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
- `CGAL::Polygon_mesh_processing::surface_intersection()`
- `CGAL::Polygon_mesh_processing::clip()`
- `CGAL::Polygon_mesh_processing::split()`
- `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()`
- `CGAL::Polygon_mesh_processing::autorefine()`
\cgalCRPSection{Meshing Functions}
- `CGAL::Polygon_mesh_processing::remesh_planar_patches()`

View File

@ -229,7 +229,7 @@ it can try to equalize the angles between incident edges, or (and) move vertices
Border vertices are considered constrained and do not move at any step of the procedure. No vertices are inserted at any time.
Angle and area smoothing algorithms are based on Surazhsky and Gotsman \cgalCite{cgal:sg-hqct-04}.
Since area smoothing considers only areas as a smoothing criterion, it may result in long and skinny
triangles. To paliate this phenomenon, area smoothing is followed by an (optional) step of Delaunay-based edge flips.
triangles. To palliate this phenomenon, area smoothing is followed by an (optional) step of Delaunay-based edge flips.
In any case, area smoothing is guaranteed to improve the spatial distribution of the vertices
over the area that is being smoothed. A simple example can be found in \ref Polygon_mesh_processing/mesh_smoothing_example.cpp.
@ -449,7 +449,7 @@ depicted in red); (iii) `clip_volume=true`: clipping of the volume bounded by th
\cgalFigureEnd
\cgalFigureBegin{coref_clip_compact, clip_compact.png}
Clipping a cube with a halfspace: compacity of the clipper (`clip_volume=false` in both cases).
Clipping a cube with a halfspace: compactivity of the clipper (`clip_volume=false` in both cases).
From left to right: (i) initial cube and the plane defining the clipping halfspace,
note that a whole face of the cube (2 triangles) is exactly contained in the plane;
(ii) `use_compact_clipper=true`: clipping of the surface mesh with a compact halfspace: coplanar faces are part of the output;
@ -879,6 +879,19 @@ vertices at identical positions, can be used to repair this configuration.
\section PMPGeometricRepair Geometric Repair
****************************************
\subsection PMPAutoref Self-intersection Resolution (Autorefinement) in Triangle Soups
Given a soup of triangles, a self-intersection is defined as the intersection of two triangles from the soup
such that the intersection is not defined by the convex hull of one, two or three shared vertices.
In other words, it is an intersection that happens in the interior of one of the two triangles, or in the interior
of one their edges, except if identical points are associated to different vertices of the triangle soup which
would then also includes overlaps of duplicated points.
The function `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` provides a way to refine a triangle soup
using the intersections of the triangles from the soup. In particular, if some points are duplicated they will be
merged. Note that if a kernel with exact predicates but inexact constructions is used, some new self-intersections
might be introduced due to rounding issues of points coordinates.
To guarantee that the triangle soup is free from self-intersections, a kernel with exact constructions must be used.
\subsection PMPRemoveCapsNeedles Removal of Almost Degenerate Triangle Faces
Triangle faces of a mesh made up of almost collinear points are badly shaped elements that
might not be desirable to have in a mesh. The function

View File

@ -47,5 +47,6 @@
\example Polygon_mesh_processing/remesh_planar_patches.cpp
\example Polygon_mesh_processing/remesh_almost_planar_patches.cpp
\example Polygon_mesh_processing/sample_example.cpp
\example Polygon_mesh_processing/soup_autorefinement.cpp
*/
*/

View File

@ -53,6 +53,8 @@ create_single_source_cgal_program("cc_compatible_orientations.cpp")
create_single_source_cgal_program("hausdorff_distance_remeshing_example.cpp")
create_single_source_cgal_program("hausdorff_bounded_error_distance_example.cpp")
create_single_source_cgal_program("isotropic_remeshing_with_custom_sizing_example.cpp")
create_single_source_cgal_program("triangle_mesh_autorefinement.cpp")
create_single_source_cgal_program("soup_autorefinement.cpp")
find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater)
include(CGAL_Eigen3_support)
@ -133,6 +135,7 @@ if(TARGET CGAL::TBB_support)
target_link_libraries(self_intersections_example PUBLIC CGAL::TBB_support)
target_link_libraries(hausdorff_distance_remeshing_example PUBLIC CGAL::TBB_support)
target_link_libraries(hausdorff_bounded_error_distance_example PUBLIC CGAL::TBB_support)
target_link_libraries(soup_autorefinement PUBLIC CGAL::TBB_support)
create_single_source_cgal_program("corefinement_parallel_union_meshes.cpp")
target_link_libraries(corefinement_parallel_union_meshes PUBLIC CGAL::TBB_support)

View File

@ -0,0 +1,38 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/autorefinement.h>
#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <CGAL/Real_timer.h>
#include <boost/container/small_vector.hpp>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char** argv)
{
const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off")
: std::string(argv[1]);
std::vector<Point> input_points;
std::vector<boost::container::small_vector<std::size_t, 3>> input_triangles;
CGAL::IO::read_polygon_soup(filename, input_points, input_triangles);
PMP::repair_polygon_soup(input_points, input_triangles);
PMP::triangulate_polygons(input_points, input_triangles);
CGAL::Real_timer t;
t.start();
PMP::autorefine_triangle_soup(input_points, input_triangles,
CGAL::parameters::concurrency_tag(CGAL::Parallel_if_available_tag()));
t.stop();
std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << " in " << t.time() << " sec." << std::endl;
CGAL::IO::write_polygon_soup("autorefined.off", input_points, input_triangles, CGAL::parameters::stream_precision(17));
return 0;
}

View File

@ -0,0 +1,29 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/autorefinement.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef CGAL::Surface_mesh<Point> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char** argv)
{
Mesh mesh;
const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off")
: std::string(argv[1]);
CGAL::IO::read_polygon_mesh(filename, mesh);
PMP::autorefine(mesh);
CGAL::IO::write_polygon_mesh("autorefined.off", mesh, CGAL::parameters::stream_precision(17));
return 0;
}

View File

@ -28,8 +28,8 @@ struct Non_manifold_feature_map
typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor;
typedef dynamic_edge_property_t<std::size_t> Edge_to_id_tag;
typedef dynamic_vertex_property_t<std::size_t> Vertex_to_id_tag;
typedef typename boost::property_map<PolygonMesh, Edge_to_id_tag>::type Edge_to_nm_id;
typedef typename boost::property_map<PolygonMesh, Vertex_to_id_tag>::type Vertex_to_nm_id;
typedef typename boost::property_map<PolygonMesh, Edge_to_id_tag>::const_type Edge_to_nm_id;
typedef typename boost::property_map<PolygonMesh, Vertex_to_id_tag>::const_type Vertex_to_nm_id;
Edge_to_nm_id e_nm_id;
Vertex_to_nm_id v_nm_id;
std::vector< std::vector<edge_descriptor> > non_manifold_edges;
@ -39,7 +39,7 @@ struct Non_manifold_feature_map
{}
template <class Vpm>
Non_manifold_feature_map(PolygonMesh& pm, Vpm vpm)
Non_manifold_feature_map(const PolygonMesh& pm, Vpm vpm)
: e_nm_id(get(Edge_to_id_tag(), pm))
, v_nm_id(get(Vertex_to_id_tag(), pm))
{
@ -99,6 +99,14 @@ struct Non_manifold_feature_map
}
}
}
void clear()
{
non_manifold_edges.clear();
non_manifold_vertices.clear();
e_nm_id = Edge_to_nm_id();
v_nm_id = Vertex_to_nm_id();
}
};
} } // end of CGAL::Polygon_mesh_processing

File diff suppressed because it is too large Load Diff

View File

@ -106,6 +106,7 @@ if(TARGET CGAL::TBB_support)
target_link_libraries(test_pmp_distance PUBLIC CGAL::TBB_support)
target_link_libraries(orient_polygon_soup_test PUBLIC CGAL::TBB_support)
target_link_libraries(self_intersection_surface_mesh_test PUBLIC CGAL::TBB_support)
target_link_libraries(test_autorefinement PUBLIC CGAL::TBB_support)
else()
message(STATUS "NOTICE: Intel TBB was not found. Tests will use sequential code.")
endif()

View File

@ -1,28 +1,28 @@
data-autoref/test_01.off 0 0 4 1 4 0
data-autoref/test_02.off 1 5 13 1 5 0
data-autoref/test_03.off 4 8 18 0 18 0
data-autoref/test_04.off 1 5 17 0 17 0
data-autoref/test_05.off 1 5 17 0 17 0
data-autoref/test_06.off 3 55 141 1 76 0
data-autoref/test_07.off 1 4 10 1 4 0
data-autoref/test_08.off 1 31 87 1 52 0
data-autoref/test_09.off 1 4 5 1 4 0
data-autoref/test_10.off 1 3 13 0 13 0
data-autoref/test_11.off 1 3 12 0 12 0
data-autoref/test_12.off 2 2 11 1 11 0
data-autoref/test_13.off 1 5 16 1 8 0
data-autoref/test_14.off 1 5 16 1 12 0
data-autoref/test_15.off 3 8 16 1 12 0
data-autoref/test_16.off 1 2 6 1 4 0
data-autoref/test_17.off 1 2 6 1 4 0
data-autoref/triple_inter_exception/triple.off 0 0 0 0 0 1
data-autoref/triple_inter_exception/cubes_cpln_1.off 0 0 0 0 0 1
data-autoref/triple_inter_exception/cubes_cpln_2.off 0 0 0 0 0 1
data-autoref/triple_inter_exception/cubes_cpln_3.off 0 0 0 0 0 1
data-autoref/open_01.off 1 65 1377 1 1313 0
data-autoref/open_02.off 1 33 595 1 562 0
data-autoref/cpln_01.off 18 42 48 1 30 0
data-autoref/cpln_02.off 28 56 40 1 24 0
data-autoref/cpln_03.off 15 35 42 1 27 0
data-autoref/four_cubes.off 12 94 184 1 78 0
data-autoref/spiral.off 7 14 26 0 26 0
data-autoref/test_01.off 0 0 4 1 4 0 4 2
data-autoref/test_02.off 1 5 13 1 5 0 10 17
data-autoref/test_03.off 4 8 18 0 18 0 14 20
data-autoref/test_04.off 1 5 17 0 17 0 13 20
data-autoref/test_05.off 1 5 17 0 17 0 13 20
data-autoref/test_06.off 3 55 141 1 76 0 92 248
data-autoref/test_07.off 1 4 10 1 4 0 8 12
data-autoref/test_08.off 1 31 87 1 52 0 57 148
data-autoref/test_09.off 1 4 5 1 4 0 4 4
data-autoref/test_10.off 1 3 13 0 13 0 10 13
data-autoref/test_11.off 1 3 12 0 12 0 9 12
data-autoref/test_12.off 2 2 11 1 11 0 9 6
data-autoref/test_13.off 1 5 16 1 8 0 16 22
data-autoref/test_14.off 1 5 16 1 12 0 16 22
data-autoref/test_15.off 3 8 16 1 12 0 16 24
data-autoref/test_16.off 1 2 6 1 4 0 6 2
data-autoref/test_17.off 1 2 6 1 4 0 6 2
data-autoref/triple_inter_exception/triple.off 0 0 0 0 0 1 15 18
data-autoref/triple_inter_exception/cubes_cpln_1.off 0 0 0 0 0 1 66 224
data-autoref/triple_inter_exception/cubes_cpln_2.off 0 0 0 0 0 1 54 196
data-autoref/triple_inter_exception/cubes_cpln_3.off 0 0 0 0 0 1 61 204
data-autoref/open_01.off 1 65 1377 1 1313 0 1317 2622
data-autoref/open_02.off 1 33 595 1 562 0 565 1056
data-autoref/cpln_01.off 18 42 48 1 30 0 30 88
data-autoref/cpln_02.off 28 56 40 1 24 0 24 72
data-autoref/cpln_03.off 15 35 42 1 27 0 27 76
data-autoref/four_cubes.off 12 94 184 1 78 0 106 352
data-autoref/spiral.off 7 14 26 0 26 0 19 31

View File

@ -3,6 +3,10 @@
#include <CGAL/Polygon_mesh_processing/intersection.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/autorefinement.h>
#include <CGAL/Polygon_mesh_processing/self_intersections.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <fstream>
#include <sstream>
@ -13,23 +17,64 @@ typedef CGAL::Surface_mesh<K::Point_3> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
template <class TriangleMesh>
struct My_visitor :
struct My_exp_visitor :
public CGAL::Polygon_mesh_processing::Corefinement::Default_visitor<TriangleMesh>
{
void after_subface_creations(TriangleMesh&){++(*i);}
My_visitor()
My_exp_visitor()
: i (new int(0) )
{}
std::shared_ptr<int> i;
};
void test(const char* fname, std::size_t nb_polylines, std::size_t total_nb_points,
struct My_visitor
{
My_visitor(std::size_t nb_input, std::size_t expected_nb_output)
: nb_input(nb_input)
, expected_nb_output(expected_nb_output)
{}
~My_visitor()
{
for(std::size_t i=0; i<tgt_check.size(); ++i)
{
assert( tgt_check[i]==1 );
}
}
void number_of_output_triangles(std::size_t nbt)
{
tgt_check.assign(expected_nb_output, 0);
assert(nbt==expected_nb_output);
}
void verbatim_triangle_copy(std::size_t tgt_id, std::size_t src_id)
{
assert(src_id<nb_input);
assert(tgt_id<expected_nb_output);
assert(tgt_check.size()==expected_nb_output && tgt_check[tgt_id]==0);
tgt_check[tgt_id]=1;
}
void new_subtriangle(std::size_t tgt_id, std::size_t src_id)
{
assert(src_id<nb_input);
assert(tgt_id<expected_nb_output);
assert(tgt_check.size()==expected_nb_output && tgt_check[tgt_id]==0);
tgt_check[tgt_id]=1;
}
std::size_t nb_input;
std::size_t expected_nb_output;
std::vector<int> tgt_check;
};
void test_coref_based(const char* fname, std::size_t nb_polylines, std::size_t total_nb_points,
std::size_t nb_vertices_after_autorefine, bool all_fixed, std::size_t nb_vertices_after_fix,
bool triple_intersection)
{
std::cout << "Running tests on " << fname << "\n";
std::cout << "Running tests (coref based) on " << fname << "\n";
std::ifstream input(fname);
Mesh mesh;
@ -41,7 +86,7 @@ void test(const char* fname, std::size_t nb_polylines, std::size_t total_nb_poin
input.close();
std::size_t nb_vertices_before_autorefine = num_vertices(mesh);
// Testing surface_self_intersection()
// Testing PMP::experimental::surface_self_intersection()
try{
std::vector< std::vector<K::Point_3> >polylines;
PMP::experimental::surface_self_intersection(mesh, std::back_inserter(polylines));
@ -57,9 +102,9 @@ void test(const char* fname, std::size_t nb_polylines, std::size_t total_nb_poin
assert( triple_intersection );
}
// Testing autorefine()
// Testing PMP::experimental::autorefine()
try{
My_visitor<Mesh> visitor;
My_exp_visitor<Mesh> visitor;
PMP::experimental::autorefine(mesh,
CGAL::parameters::visitor(visitor));
mesh.collect_garbage();
@ -72,7 +117,7 @@ void test(const char* fname, std::size_t nb_polylines, std::size_t total_nb_poin
assert( triple_intersection );
}
// Testing autorefine_and_remove_self_intersections()
// Testing PMP::experimental::autorefine_and_remove_self_intersections()
try{
input.open(fname);
mesh.clear();
@ -89,10 +134,42 @@ void test(const char* fname, std::size_t nb_polylines, std::size_t total_nb_poin
}
}
template <class Tag>
void test(const char* fname, std::size_t nb_vertices_after_autorefine, std::size_t expected_nb_output, Tag tag)
{
std::cout << "Running tests on " << fname;
if (std::is_same_v<Tag, CGAL::Sequential_tag>)
std::cout << " (Sequential)\n";
else
std::cout << " (Parallel)\n";
std::vector<K::Point_3> points;
std::vector< std::vector<std::size_t> > triangles;
if (!CGAL::IO::read_polygon_soup(fname, points, triangles))
{
std::cerr << " Input mesh is not a valid file." << std::endl;
exit(EXIT_FAILURE);
}
// Testing autorefine()
My_visitor visitor(triangles.size(), expected_nb_output);
PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::visitor(visitor).concurrency_tag(tag));
assert( nb_vertices_after_autorefine==points.size());
assert( expected_nb_output==triangles.size());
assert( !PMP::does_triangle_soup_self_intersect(points, triangles) );
// CGAL::IO::write_polygon_soup("/tmp/debug.off", points, triangles);
}
int main(int argc, const char** argv)
{
// file nb_polylines total_nb_points nb_vertices_after_autorefine all_fixed nb_vertices_after_fix triple_intersection
for (int i=0;i<(argc-1)/7; ++i)
test(argv[1+7*i], atoi(argv[1+7*i+1]), atoi(argv[1+7*i+2]),
atoi(argv[1+7*i+3]), atoi(argv[1+7*i+4])==0?false:true, atoi(argv[1+7*i+5]), atoi(argv[1+7*i+6])==0?false:true);
for (int i=0;i<(argc-1)/9; ++i)
{
test_coref_based(argv[1+9*i], atoi(argv[1+9*i+1]), atoi(argv[1+9*i+2]),
atoi(argv[1+9*i+3]), atoi(argv[1+9*i+4])==0?false:true, atoi(argv[1+9*i+5]), atoi(argv[1+9*i+6])==0?false:true);
test(argv[1+9*i], atoi(argv[1+9*i+7]), atoi(argv[1+9*i+8]), CGAL::Sequential_tag());
#ifdef CGAL_LINKED_WITH_TBB
test(argv[1+9*i], atoi(argv[1+9*i+7]), atoi(argv[1+9*i+8]), CGAL::Parallel_tag());
#endif
}
}

View File

@ -133,7 +133,11 @@ target_link_libraries(
qt6_wrap_ui( repairUI_FILES RemoveNeedlesDialog.ui SelfSnapDialog.ui)
polyhedron_demo_plugin(repair_polyhedron_plugin Repair_polyhedron_plugin ${repairUI_FILES} KEYWORDS PMP)
target_link_libraries(repair_polyhedron_plugin PUBLIC scene_points_with_normal_item scene_surface_mesh_item)
target_link_libraries(repair_polyhedron_plugin PUBLIC scene_points_with_normal_item scene_surface_mesh_item scene_polygon_soup_item)
if(TARGET CGAL::TBB_support)
target_link_libraries(repair_polyhedron_plugin PUBLIC CGAL::TBB_support)
endif()
if(TARGET CGAL::Eigen3_support)
qt6_wrap_ui(isotropicRemeshingUI_FILES Isotropic_remeshing_dialog.ui)

View File

@ -1,6 +1,7 @@
#include <QtCore/qglobal.h>
#include "Scene_surface_mesh_item.h"
#include "Scene_polygon_soup_item.h"
#include "Scene_points_with_normal_item.h"
#include <CGAL/Three/Scene_interface.h>
#include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
@ -19,6 +20,8 @@
#include <CGAL/Polygon_mesh_processing/repair_degeneracies.h>
#include <CGAL/Polygon_mesh_processing/merge_border_vertices.h>
#include <CGAL/Polygon_mesh_processing/internal/Snapping/snap.h>
#include <CGAL/Polygon_mesh_processing/autorefinement.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include "ui_RemoveNeedlesDialog.h"
#include "ui_SelfSnapDialog.h"
@ -52,8 +55,9 @@ public:
actionDuplicateNMVertices = new QAction(tr("Duplicate Non-Manifold Vertices"), mw);
actionExtractNMVertices = new QAction(tr("Extract Non-Manifold Vertices"), mw);
actionMergeDuplicatedVerticesOnBoundaryCycles = new QAction(tr("Merge Duplicated Vertices on Boundary Cycles"), mw);
actionAutorefine = new QAction(tr("Autorefine Mesh"), mw);
actionAutorefineAndRMSelfIntersections = new QAction(tr("Autorefine and Remove Self-Intersections"), mw);
actionAutorefine = new QAction(tr("Autorefine Mesh (Deprecated)"), mw);
actionNewAutorefine = new QAction(tr("Autorefine"), mw);
actionAutorefineAndRMSelfIntersections = new QAction(tr("Autorefine and Remove Self-Intersections (Deprecated)"), mw);
actionRemoveNeedlesAndCaps = new QAction(tr("Remove Needles And Caps"));
actionSnapBorders = new QAction(tr("Snap Boundaries"));
@ -65,6 +69,7 @@ public:
actionExtractNMVertices->setObjectName("actionExtractNMVertices");
actionMergeDuplicatedVerticesOnBoundaryCycles->setObjectName("actionMergeDuplicatedVerticesOnBoundaryCycles");
actionAutorefine->setObjectName("actionAutorefine");
actionNewAutorefine->setObjectName("actionNewAutorefine");
actionAutorefineAndRMSelfIntersections->setObjectName("actionAutorefineAndRMSelfIntersections");
actionRemoveNeedlesAndCaps->setObjectName("actionRemoveNeedlesAndCaps");
actionSnapBorders->setObjectName("actionSnapBorders");
@ -77,6 +82,7 @@ public:
actionExtractNMVertices->setProperty("subMenuName", "Polygon Mesh Processing/Repair");
actionMergeDuplicatedVerticesOnBoundaryCycles->setProperty("subMenuName", "Polygon Mesh Processing/Repair");
actionAutorefine->setProperty("subMenuName", "Polygon Mesh Processing/Repair/Experimental");
actionNewAutorefine->setProperty("subMenuName", "Polygon Mesh Processing/Repair");
actionAutorefineAndRMSelfIntersections->setProperty("subMenuName", "Polygon Mesh Processing/Repair/Experimental");
actionSnapBorders->setProperty("subMenuName", "Polygon Mesh Processing/Repair/Experimental");
@ -93,16 +99,29 @@ public:
<< actionExtractNMVertices
<< actionMergeDuplicatedVerticesOnBoundaryCycles
<< actionAutorefine
<< actionNewAutorefine
<< actionAutorefineAndRMSelfIntersections
<< actionRemoveNeedlesAndCaps
<< actionSnapBorders;
}
bool applicable(QAction*) const
bool applicable(QAction* action) const
{
if (action!=actionNewAutorefine)
{
int item_id = scene->mainSelectionIndex();
return qobject_cast<Scene_surface_mesh_item*>(scene->item(item_id));
}
for (Scene_interface::Item_id index : scene->selectionIndices())
{
if (qobject_cast<Scene_surface_mesh_item*>(scene->item(index)))
return true;
if (qobject_cast<Scene_polygon_soup_item*>(scene->item(index)))
return true;
}
return false;
}
template <typename Item>
void on_actionRemoveIsolatedVertices_triggered(Scene_interface::Item_id index);
template <typename Item>
@ -120,6 +139,8 @@ public:
template <typename Item>
void on_actionAutorefine_triggered(Scene_interface::Item_id index);
template <typename Item>
void on_actionNewAutorefine_triggered(const std::vector<Scene_interface::Item_id>& indices);
template <typename Item>
void on_actionAutorefineAndRMSelfIntersections_triggered(Scene_interface::Item_id index);
public Q_SLOTS:
@ -131,6 +152,7 @@ public Q_SLOTS:
void on_actionExtractNMVertices_triggered();
void on_actionMergeDuplicatedVerticesOnBoundaryCycles_triggered();
void on_actionAutorefine_triggered();
void on_actionNewAutorefine_triggered();
void on_actionAutorefineAndRMSelfIntersections_triggered();
void on_actionRemoveNeedlesAndCaps_triggered();
void on_actionSnapBorders_triggered();
@ -144,6 +166,7 @@ private:
QAction* actionExtractNMVertices;
QAction* actionMergeDuplicatedVerticesOnBoundaryCycles;
QAction* actionAutorefine;
QAction* actionNewAutorefine;
QAction* actionAutorefineAndRMSelfIntersections;
QAction* actionRemoveNeedlesAndCaps;
QAction* actionSnapBorders;
@ -421,6 +444,78 @@ void Polyhedron_demo_repair_polyhedron_plugin::on_actionAutorefine_triggered()
QApplication::restoreOverrideCursor();
}
template <typename Item>
void Polyhedron_demo_repair_polyhedron_plugin::on_actionNewAutorefine_triggered(const std::vector<Scene_interface::Item_id>& indices)
{
namespace PMP = CGAL::Polygon_mesh_processing;
Polygon_soup::Points points;
Polygon_soup::Polygons polygons;
if (indices.size()==1)
{
if (Scene_surface_mesh_item* smi_ptr = qobject_cast<Scene_surface_mesh_item*>(scene->item(indices[0])))
PMP::polygon_mesh_to_polygon_soup(*smi_ptr->polyhedron(), points, polygons);
else if (Scene_polygon_soup_item* spi_ptr = qobject_cast<Scene_polygon_soup_item*>(scene->item(indices[0])))
{
points = spi_ptr->points();
polygons = spi_ptr->polygons();
}
}
else
{
for (Scene_interface::Item_id id : indices)
{
Polygon_soup::Points l_points;
Polygon_soup::Polygons l_polygons;
if (Scene_surface_mesh_item* smi_ptr = qobject_cast<Scene_surface_mesh_item*>(scene->item(id)))
PMP::polygon_mesh_to_polygon_soup(*smi_ptr->polyhedron(), l_points, l_polygons);
else if (Scene_polygon_soup_item* spi_ptr = qobject_cast<Scene_polygon_soup_item*>(scene->item(id)))
{
l_points = spi_ptr->points();
l_polygons = spi_ptr->polygons();
}
std::size_t offset=points.size();
points.insert(points.end(), l_points.begin(), l_points.end());
std::size_t psize=polygons.size();
polygons.insert(polygons.end(), l_polygons.begin(), l_polygons.end());
for (std::size_t i=psize; i<polygons.size(); ++i)
for(std::size_t& id : polygons[i])
id+=offset;
}
}
PMP::triangulate_polygons(points, polygons);
PMP::autorefine_triangle_soup(points, polygons,
CGAL::parameters::concurrency_tag(CGAL::Parallel_if_available_tag()));
Scene_polygon_soup_item* new_item = new Scene_polygon_soup_item();
new_item->load(points, polygons);
QString name = scene->item(indices[0])->name();
for (std::size_t k=1; k<indices.size(); ++k)
name += " + " + scene->item(indices[k])->name();
new_item->setName(name+" autorefined");
scene->addItem(new_item);
new_item->invalidateOpenGLBuffers();
Q_EMIT new_item->itemChanged();
}
void Polyhedron_demo_repair_polyhedron_plugin::on_actionNewAutorefine_triggered()
{
std::vector<Scene_interface::Item_id> indices;
for (Scene_interface::Item_id index : scene->selectionIndices())
{
if (qobject_cast<Scene_surface_mesh_item*>(scene->item(index)))
indices.push_back(index);
else if (qobject_cast<Scene_polygon_soup_item*>(scene->item(index)))
indices.push_back(index);
}
QApplication::setOverrideCursor(Qt::WaitCursor);
on_actionNewAutorefine_triggered<Scene_surface_mesh_item>(indices);
QApplication::restoreOverrideCursor();
}
template <typename Item>
void Polyhedron_demo_repair_polyhedron_plugin::on_actionAutorefineAndRMSelfIntersections_triggered(Scene_interface::Item_id index)
{

View File

@ -661,7 +661,7 @@ struct Boolean_property_map
return pm.set_ptr->count(k) != 0;
}
friend void put(Boolean_property_map<Set>& pm, const key_type& k, bool v)
friend void put(Boolean_property_map<Set> pm, const key_type& k, bool v)
{
CGAL_assertion(pm.set_ptr!=nullptr);
if (v)

View File

@ -47,6 +47,7 @@ CGAL_add_named_parameter(implementation_tag_t, implementation_tag, implementatio
CGAL_add_named_parameter(prevent_unselection_t, prevent_unselection, prevent_unselection)
CGAL_add_named_parameter(verbose_t, verbose, verbose)
CGAL_add_named_parameter(concurrency_tag_t, concurrency_tag, concurrency_tag)
// List of named parameters used for IO
CGAL_add_named_parameter(vertex_normal_output_iterator_t, vertex_normal_output_iterator, vertex_normal_output_iterator)

View File

@ -330,11 +330,11 @@ public:
#if 1
template <class Segment_2>
static decltype(auto) get_source(const Segment_2& segment){
return segment.source();
return segment[0];
}
template <class Segment_2>
static decltype(auto) get_target(const Segment_2& segment){
return segment.target();
return segment[1];
}
static const Point& get_source(const Constraint& cst){