mirror of https://github.com/CGAL/cgal
transform to same API as existing snap_rounding, remove duplicate, efficient subdivision of segments by vertices
This commit is contained in:
Léo Valque
parent
173ae7b091
commit
012ec64c61
|
|
@ -1,71 +1,46 @@
|
|||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
|
||||
#include <CGAL/Float_snap_rounding_2.h>
|
||||
#include <CGAL/Arr_segment_traits_2.h>
|
||||
#include <CGAL/Surface_sweep_2_algorithms.h>
|
||||
|
||||
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
|
||||
typedef Kernel::Segment_2 Segment_2;
|
||||
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
|
||||
typedef CGAL::Arr_segment_traits_2<Kernel> Traits_2;
|
||||
typedef Traits_2::Curve_2 Segment_2;
|
||||
typedef Kernel::Point_2 Point_2;
|
||||
typedef Kernel::FT FT;
|
||||
|
||||
int main()
|
||||
{
|
||||
std::vector<Point_2> pts;
|
||||
std::vector< std::vector<size_t> > segs;
|
||||
std::vector< Segment_2 > segs;
|
||||
|
||||
Kernel::FT eps(std::pow(2, -60));
|
||||
Kernel::FT one(1.);
|
||||
Kernel::FT two(2.);
|
||||
Kernel::FT e(std::pow(2, -60));
|
||||
|
||||
pts.emplace_back(one-eps, one);
|
||||
pts.emplace_back(-one+eps, -one+2*eps);
|
||||
pts.emplace_back(eps/2, eps/2);
|
||||
pts.emplace_back(one, -one);
|
||||
pts.emplace_back(0, two-eps/2);
|
||||
pts.emplace_back(two, 0);
|
||||
pts.emplace_back(-two+eps, -FT(4.));
|
||||
pts.emplace_back(two, two);
|
||||
pts.emplace_back(-two, two);
|
||||
|
||||
segs.push_back(std::vector<size_t>());
|
||||
segs[0].push_back(0);
|
||||
segs[0].push_back(1);
|
||||
segs.push_back(std::vector<size_t>());
|
||||
segs[1].push_back(2);
|
||||
segs[1].push_back(3);
|
||||
segs.push_back(std::vector<size_t>());
|
||||
segs[2].push_back(4);
|
||||
segs[2].push_back(5);
|
||||
segs.push_back(std::vector<size_t>());
|
||||
segs[3].push_back(4);
|
||||
segs[3].push_back(6);
|
||||
segs.push_back(std::vector<size_t>());
|
||||
segs[4].push_back(7);
|
||||
segs[4].push_back(8);
|
||||
segs.emplace_back(Point_2(1-e, 1), Point_2(-1-e, -1+2*e));
|
||||
segs.emplace_back(Point_2(e/2, e/2), Point_2(1, -1));
|
||||
segs.emplace_back(Point_2(0, 2-e/2), Point_2(2, 0));
|
||||
segs.emplace_back(Point_2(0, 2-e/2), Point_2(-2+e, -4));
|
||||
segs.emplace_back(Point_2(-2, 2), Point_2(2, 2));
|
||||
segs.emplace_back(Point_2(7, 7), Point_2(7+e, 7+e));
|
||||
segs.emplace_back(Point_2(5, 7-e), Point_2(9, 7-e));
|
||||
|
||||
std::cout << "Input" << std::endl;
|
||||
std::cout << "Points" << std::endl;
|
||||
for(auto &p: pts)
|
||||
std::cout << p << std::endl;
|
||||
std::cout << "Segments:" << std::endl;
|
||||
for(auto &seg: segs){
|
||||
for(auto pi : seg)
|
||||
std::cout << pi << " ";
|
||||
std::cout << std::endl;
|
||||
std::cout << seg.source() << " -> " << seg.target() << std::endl;
|
||||
}
|
||||
std::cout << "\n\n" << std::endl;
|
||||
|
||||
CGAL::float_snap_rounding_2(pts, segs);
|
||||
std::vector< std::vector<Point_2 > > out;
|
||||
CGAL::double_snap_rounding_2(segs.begin(), segs.end(), out);
|
||||
|
||||
std::cout << "Output" << std::endl;
|
||||
std::cout << "Points" << std::endl;
|
||||
for(auto &p: pts)
|
||||
std::cout << p << std::endl;
|
||||
std::cout << "Segments:" << std::endl;
|
||||
for(auto &seg: segs){
|
||||
for(auto pi : seg)
|
||||
std::cout << pi << " ";
|
||||
for(auto &poly: out){
|
||||
std::cout << poly[0];
|
||||
for(size_t i=1; i!=poly.size(); ++i)
|
||||
std::cout << " -> " << poly[i];
|
||||
std::cout << std::endl;
|
||||
}
|
||||
|
||||
return(0);
|
||||
return 0;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -13,13 +13,11 @@
|
|||
#ifndef CGAL_FLOAT_SNAP_ROUNDING_2_H
|
||||
#define CGAL_FLOAT_SNAP_ROUNDING_2_H
|
||||
|
||||
#include <CGAL/license/Snap_rounding_2.h>
|
||||
|
||||
|
||||
#include <iostream>
|
||||
#include <CGAL/basic.h>
|
||||
#include <CGAL/enum.h>
|
||||
#include <CGAL/predicates_on_points_2.h>
|
||||
// #include <CGAL/predicates_on_points_2.h>
|
||||
#include <CGAL/intersection_2.h>
|
||||
#include <CGAL/Surface_sweep_2_algorithms.h>
|
||||
#include <CGAL/Box_intersection_d/Box_d.h>
|
||||
|
|
@ -96,85 +94,246 @@ public:
|
|||
#endif
|
||||
}
|
||||
|
||||
template <class PointRange, class PolylineRange>
|
||||
void float_snap_rounding_2(PointRange& points,
|
||||
PolylineRange& segments)
|
||||
template <class PointsRange , class PolylineRange>
|
||||
void double_snap_rounding_2_disjoint(PointsRange &pts,
|
||||
PolylineRange &polylines)
|
||||
{
|
||||
using Point_2 = std::remove_cv_t<typename std::iterator_traits<typename PointRange::const_iterator>::value_type>;
|
||||
using Point_2 = std::remove_cv_t<typename std::iterator_traits<typename PointsRange::iterator>::value_type>;
|
||||
using Kernel = typename Kernel_traits<Point_2>::Kernel;
|
||||
using Segment_2 = typename Kernel::Segment_2;
|
||||
using Vector_2 = typename Kernel::Vector_2;
|
||||
using Comparison_result = typename Kernel::Comparison_result;
|
||||
using Segment_2 = typename Kernel::Segment_2;
|
||||
using Vector_2 = typename Kernel::Vector_2;
|
||||
using PBox = CGAL::Box_intersection_d::Box_with_index_d<double,2>;
|
||||
using SBox = CGAL::Box_intersection_d::Box_with_index_d<double,2>;
|
||||
|
||||
auto round_bound=[](Point_2& p){
|
||||
return std::pow(p.x().approx().sup()-p.x().approx().inf(),2)+std::pow(p.y().approx().sup()-p.y().approx().inf(), 2);
|
||||
auto comp_by_x=[&](size_t ai, size_t bi){
|
||||
return compare(pts[ai].x(),pts[bi].x())==SMALLER;
|
||||
};
|
||||
auto comp_by_y=[&](size_t ai, size_t bi){
|
||||
return compare(pts[ai].y(),pts[bi].y())==SMALLER;
|
||||
};
|
||||
auto comp_by_x_first=[&](size_t ai, size_t bi){
|
||||
Comparison_result res=compare(pts[ai].x(),pts[bi].x());
|
||||
if(res==EQUAL)
|
||||
return compare(pts[ai].y(),pts[bi].y())==SMALLER;
|
||||
return res==SMALLER;
|
||||
};
|
||||
auto comp_by_y_first=[&](size_t ai, size_t bi){
|
||||
Comparison_result res=compare(pts[ai].y(),pts[bi].y());
|
||||
if(res==EQUAL)
|
||||
return compare(pts[ai].x(),pts[bi].x())==SMALLER;
|
||||
return res==SMALLER;
|
||||
};
|
||||
|
||||
|
||||
// Compute the order of the points along the 2 axis
|
||||
// Sorted the points may perform exact computations and thus refine the intervals of the coordinates values
|
||||
// This refine ensures that the order of the points will be preserved by the rounding
|
||||
using Iterator_set_x = typename std::set<size_t, decltype(comp_by_x_first)>::iterator;
|
||||
using Iterator_set_y = typename std::set<size_t, decltype(comp_by_y_first)>::iterator;
|
||||
std::set<size_t, decltype(comp_by_x_first)> p_sort_by_x(comp_by_x_first);
|
||||
std::set<size_t, decltype(comp_by_y_first)> p_sort_by_y(comp_by_y_first);
|
||||
for(size_t i=0; i!=pts.size(); ++i)
|
||||
{
|
||||
p_sort_by_x.insert(i);
|
||||
p_sort_by_y.insert(i);
|
||||
}
|
||||
// std::vector<size_t> p_sort_by_y(pts.size(),0);
|
||||
// std::iota(p_sort_by_x.begin(),p_sort_by_x.end(),0);
|
||||
// std::iota(p_sort_by_y.begin(),p_sort_by_y.end(),0);
|
||||
// std::sort(p_sort_by_x.begin(),p_sort_by_x.end(),comp_by_x);
|
||||
// std::sort(p_sort_by_y.begin(),p_sort_by_y.end(),comp_by_y);
|
||||
|
||||
//Kd-Tree to exhibits pairs
|
||||
std::vector<PBox> points_boxes;
|
||||
std::vector<SBox> segs_boxes;
|
||||
for(size_t i=0; i<points.size(); ++i)
|
||||
points_boxes.emplace_back(points[i].bbox(),i);
|
||||
for(size_t i=0; i<segments.size(); ++i)
|
||||
segs_boxes.emplace_back(points[segments[i][0]].bbox()+points[segments[i][1]].bbox(),i);
|
||||
for(size_t i=0; i<pts.size(); ++i)
|
||||
points_boxes.emplace_back(pts[i].bbox(),i);
|
||||
for(size_t i=0; i<polylines.size(); ++i)
|
||||
segs_boxes.emplace_back(pts[polylines[i][0]].bbox()+pts[polylines[i][1]].bbox(),i);
|
||||
|
||||
auto round_bound=[](Point_2& p){
|
||||
return std::pow(p.x().approx().sup()-p.x().approx().inf(),2)+std::pow(p.y().approx().sup()-p.y().approx().inf(), 2);
|
||||
};
|
||||
auto callback=[&](PBox &bp, SBox &bseg){
|
||||
size_t pi=bp.index();
|
||||
size_t si=bseg.index();
|
||||
size_t si1=segments[bseg.index()][0];
|
||||
size_t si2=segments[bseg.index()][1];
|
||||
size_t si1=polylines[bseg.index()][0];
|
||||
size_t si2=polylines[bseg.index()][1];
|
||||
|
||||
// the point is a vertex of the segment
|
||||
if((pi==si1) || (pi==si2))
|
||||
return;
|
||||
|
||||
Point_2& p= points[bp.index()];
|
||||
Segment_2 seg(points[segments[bseg.index()][0]], points[segments[bseg.index()][1]]);
|
||||
Point_2& p= pts[bp.index()];
|
||||
Segment_2 seg(pts[polylines[bseg.index()][0]], pts[polylines[bseg.index()][1]]);
|
||||
|
||||
// Compute the maximum distance between a point of s and its rounded version
|
||||
// TODO compute one a the beginning for performance
|
||||
// round_bound_si1= (std::max)(points_boxes[si1].bbox().xmax()-points_boxes[si1].bbox().xmin(),
|
||||
// points_boxes[si1].bbox().ymax()-points_boxes[si1].bbox().ymin());
|
||||
// round_bound_si2= (std::max)(points_boxes[si2].bbox().xmax()-points_boxes[si2].bbox().xmin(),
|
||||
// points_boxes[si2].bbox().ymax()-points_boxes[si2].bbox().ymin());
|
||||
// round_bound_s=(std::max)(round_bound_si1,round_bound_si2);
|
||||
double round_bound_s=(std::max)(round_bound(points[si1]), round_bound(points[si2]));
|
||||
double round_bound_s=(std::max)(round_bound(pts[si1]), round_bound(pts[si2]));
|
||||
|
||||
if(possibly(Kernel().compare_squared_distance_2_object()(p, seg, round_bound_s)!=CGAL::LARGER)){
|
||||
points.emplace_back(p.x(), seg.supporting_line().y_at_x(p.x()));
|
||||
segments[si].emplace_back(points.size()-1);
|
||||
pts.emplace_back(p.x(), seg.supporting_line().y_at_x(p.x()));
|
||||
auto pair=p_sort_by_x.insert(pts.size()-1);
|
||||
if(pair.second)
|
||||
p_sort_by_y.insert(pts.size()-1);
|
||||
else
|
||||
pts.pop_back();
|
||||
polylines[si].emplace_back(*pair.first);
|
||||
}
|
||||
};
|
||||
|
||||
CGAL::box_intersection_d(points_boxes.begin(), points_boxes.end(), segs_boxes.begin(), segs_boxes.end(), callback);
|
||||
|
||||
//sort new vertices
|
||||
for(auto &seg: segments)
|
||||
for(auto &polyline: polylines)
|
||||
{
|
||||
if(seg.size()==2)
|
||||
if(polyline.size()==2)
|
||||
continue;
|
||||
Vector_2 ref(points[seg[0]], points[seg[1]]);
|
||||
|
||||
//Sort the points on the polyline along the original vector
|
||||
Vector_2 ref(pts[polyline[0]], pts[polyline[1]]);
|
||||
auto sort_along_ref=[&](size_t pi, size_t qi){
|
||||
Vector_2 v(points[pi], points[qi]);
|
||||
Vector_2 v(pts[pi], pts[qi]);
|
||||
if(is_zero(ref.x()))
|
||||
return is_positive(v.y()*ref.y());
|
||||
return is_positive(v.x()*ref.x());
|
||||
};
|
||||
//Il faut conserver l'ordre des deux sommets d'entrées
|
||||
std::sort(seg.begin(), seg.end(), sort_along_ref);
|
||||
std::sort(polyline.begin(), polyline.end(), sort_along_ref);
|
||||
}
|
||||
|
||||
//round
|
||||
// TODO bug potentiel, l'algo suppose que deux coordonnées en x ne peuvent pas s'inverser mais
|
||||
// l'intervalle d'un nombre pourrait être très petit et l'intervalle d'un autre très grand et provoquer ce genre
|
||||
// d'inversion
|
||||
for(auto &p: points)
|
||||
for(auto &p: pts)
|
||||
p=Point_2(to_double(p.x()), to_double(p.y()));
|
||||
|
||||
//repair some vertices can be on segment
|
||||
//The order of the points on an axis must be preserved for the correctness of the algorithm
|
||||
CGAL_assertion(std::is_sorted(p_sort_by_x.begin(),p_sort_by_x.end(),comp_by_x));
|
||||
CGAL_assertion(std::is_sorted(p_sort_by_y.begin(),p_sort_by_y.end(),comp_by_y));
|
||||
|
||||
//remove duplicate_points
|
||||
std::vector< size_t > unique_points(p_sort_by_x.begin(),p_sort_by_x.end());
|
||||
std::sort(unique_points.begin(),unique_points.end(),comp_by_x_first);
|
||||
std::vector<Point_2> new_pts;
|
||||
std::vector<size_t> old_to_new_index(pts.size());
|
||||
for(size_t i=0; i!=pts.size(); ++i){
|
||||
if(i==0 || (pts[unique_points[i]]!=pts[unique_points[i-1]]))
|
||||
new_pts.push_back(pts[unique_points[i]]);
|
||||
old_to_new_index[unique_points[i]]=new_pts.size()-1;
|
||||
}
|
||||
|
||||
std::swap(pts, new_pts);
|
||||
// Update the polylines by remapping the old indices to new indices
|
||||
for (auto& polyline : polylines) {
|
||||
std::vector<size_t> updated_polyline;
|
||||
for (size_t i=0; i<polyline.size(); ++i) {
|
||||
size_t new_pi=old_to_new_index[polyline[i]];
|
||||
if(i==0 || (new_pi!=updated_polyline[updated_polyline.size()-1]))
|
||||
updated_polyline.push_back(new_pi);
|
||||
assert(new_pi<pts.size());
|
||||
}
|
||||
std::swap(polyline, updated_polyline);
|
||||
}
|
||||
|
||||
//Vertices can be on vertical or horizontal segment, we repair this
|
||||
p_sort_by_x.clear();
|
||||
p_sort_by_y.clear();
|
||||
for(size_t i=0; i!=pts.size(); ++i)
|
||||
{
|
||||
p_sort_by_x.insert(i);
|
||||
p_sort_by_y.insert(i);
|
||||
}
|
||||
// std::vector< size_t > vec_sort_by_x_first(p_sort_by_x.begin(),p_sort_by_x.end());
|
||||
// std::vector< size_t > vec_sort_by_y_first(p_sort_by_y.begin(),p_sort_by_y.end());
|
||||
// std::sort(vec_sort_by_x_first.begin(),vec_sort_by_x_first.end(),comp_by_x_first);
|
||||
// std::sort(vec_sort_by_y_first.begin(),vec_sort_by_y_first.end(),comp_by_y_first);
|
||||
for(auto &poly: polylines){
|
||||
std::vector<size_t> updated_polyline;
|
||||
updated_polyline.push_back(poly[0]);
|
||||
for(size_t i=1; i!=poly.size(); ++i){
|
||||
if(pts[poly[i-1]].x()==pts[poly[i]].x()){
|
||||
Iterator_set_x start, end;
|
||||
if(comp_by_x_first(poly[i-1],poly[i])){
|
||||
start=p_sort_by_x.upper_bound(poly[i-1]);
|
||||
end=p_sort_by_x.lower_bound(poly[i]);
|
||||
} else {
|
||||
start=p_sort_by_x.upper_bound(poly[i]);
|
||||
end=p_sort_by_x.lower_bound(poly[i-1]);
|
||||
}
|
||||
for(auto it=start; it!=end; ++it){
|
||||
updated_polyline.push_back(*it);
|
||||
}
|
||||
}
|
||||
if(pts[poly[i-1]].y()==pts[poly[i]].y()){
|
||||
Iterator_set_y start, end;
|
||||
if(comp_by_y_first(poly[i-1],poly[i])){
|
||||
start=p_sort_by_y.upper_bound(poly[i-1]);
|
||||
end=p_sort_by_y.lower_bound(poly[i]);
|
||||
} else {
|
||||
start=p_sort_by_y.upper_bound(poly[i]);
|
||||
end=p_sort_by_y.lower_bound(poly[i-1]);
|
||||
}
|
||||
for(auto it=start; it!=end; ++it){
|
||||
updated_polyline.push_back(*it);
|
||||
}
|
||||
}
|
||||
updated_polyline.push_back(poly[i]);
|
||||
}
|
||||
std::swap(poly, updated_polyline);
|
||||
}
|
||||
// compute_subcurves(input_begin, input_end, polylines);
|
||||
}
|
||||
|
||||
template <class InputIterator , class OutputContainer>
|
||||
typename OutputContainer::iterator double_snap_rounding_2(InputIterator input_begin,
|
||||
InputIterator input_end,
|
||||
OutputContainer& output)
|
||||
{
|
||||
using Segment_2 = std::remove_cv_t<typename std::iterator_traits<InputIterator>::value_type>;
|
||||
using Polyline = std::remove_cv_t<typename std::iterator_traits<typename OutputContainer::iterator>::value_type>;
|
||||
using Point_2 = typename Default_arr_traits<Segment_2>::Traits::Point_2;
|
||||
|
||||
std::vector<Segment_2> segs;
|
||||
compute_subcurves(input_begin, input_end, std::back_inserter(segs));
|
||||
|
||||
std::set<Point_2> unique_point_set;
|
||||
std::map<Point_2, int> point_to_index;
|
||||
std::vector<Point_2> pts;
|
||||
std::vector< std::vector< size_t> > polylines;
|
||||
|
||||
// Transform range of the segments in the range of points and polyline of indexes
|
||||
for(typename std::vector<Segment_2>::iterator it=segs.begin(); it!=segs.end(); ++it)
|
||||
{
|
||||
const Point_2& p1 = it->source();
|
||||
const Point_2& p2 = it->target();
|
||||
|
||||
// Check and insert the first endpoint if it's not already added
|
||||
if (unique_point_set.find(p1) == unique_point_set.end()) {
|
||||
unique_point_set.insert(p1);
|
||||
pts.push_back(p1);
|
||||
point_to_index[p1] = pts.size() - 1;
|
||||
}
|
||||
if (unique_point_set.find(p2) == unique_point_set.end()) {
|
||||
unique_point_set.insert(p2);
|
||||
pts.push_back(p2);
|
||||
point_to_index[p2] = pts.size() - 1;
|
||||
}
|
||||
}
|
||||
|
||||
for(InputIterator it=input_begin; it!=input_end; ++it)
|
||||
{
|
||||
size_t index1 = point_to_index[it->source()];
|
||||
size_t index2 = point_to_index[it->target()];
|
||||
polylines.push_back({index1, index2});
|
||||
}
|
||||
|
||||
double_snap_rounding_2_disjoint(pts, polylines);
|
||||
for(auto &poly: polylines){
|
||||
Polyline new_line;
|
||||
for(size_t pi: poly){
|
||||
new_line.push_back(pts[pi]);
|
||||
}
|
||||
output.push_back(new_line);
|
||||
}
|
||||
|
||||
return output.begin();
|
||||
}
|
||||
|
||||
} //namespace CGAL
|
||||
|
|
|
|||
|
|
@ -198,7 +198,7 @@ protected:
|
|||
|
||||
/*! Compute intersections between the two given curves.
|
||||
* If the two curves intersect, create a new event (or use the event that
|
||||
* already exits in the intersection point) and insert the curves to the
|
||||
* already exists in the intersection point) and insert the curves to the
|
||||
* event.
|
||||
* \param curve1 The first curve.
|
||||
* \param curve2 The second curve.
|
||||
|
|
|
|||
|
|
@ -76,6 +76,7 @@ public:
|
|||
m_overlapping(overlapping)
|
||||
{}
|
||||
|
||||
|
||||
template <typename CurveIterator>
|
||||
void sweep(CurveIterator begin, CurveIterator end)
|
||||
{
|
||||
|
|
|
|||
Loading…
Reference in New Issue