songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'cgal/6.0.x-branch'

This commit is contained in:
Sébastien Loriot 2024-11-05 17:09:53 +01:00
parent a02519677a bb8eb338a1
commit 3d2d1782cd
8 changed files with 426 additions and 572 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

718
Arrangement_on_surface_2/include/CGAL/Arr_polycurve_traits_2.h
View File

@ -26,7 +26,6 @@
#include <iterator>
#include <type_traits>
#include <variant>
#include <CGAL/basic.h>
@ -41,77 +40,75 @@ namespace CGAL {
template <typename SubcurveTraits_2 = Arr_segment_traits_2<> >
class Arr_polycurve_traits_2 :
public Arr_polycurve_basic_traits_2<SubcurveTraits_2>
{
public Arr_polycurve_basic_traits_2<SubcurveTraits_2> {
public:
typedef SubcurveTraits_2 Subcurve_traits_2;
using Subcurve_traits_2 = SubcurveTraits_2;
private:
typedef Arr_polycurve_basic_traits_2<Subcurve_traits_2> Base;
using Base = Arr_polycurve_basic_traits_2<Subcurve_traits_2>;
public:
/// \name Types inherited from the polycurve basic traits class.
//@{
typedef typename Base::Has_left_category Has_left_category;
typedef typename Base::Has_do_intersect_category
Has_do_intersect_category;
using Has_left_category = typename Base::Has_left_category;
using Has_do_intersect_category = typename Base::Has_do_intersect_category;
typedef typename Base::Left_side_category Left_side_category;
typedef typename Base::Bottom_side_category Bottom_side_category;
typedef typename Base::Top_side_category Top_side_category;
typedef typename Base::Right_side_category Right_side_category;
using Left_side_category = typename Base::Left_side_category;
using Bottom_side_category = typename Base::Bottom_side_category;
using Top_side_category = typename Base::Top_side_category;
using Right_side_category = typename Base::Right_side_category;
typedef typename Base::All_sides_oblivious_category
All_sides_oblivious_category;
using All_sides_oblivious_category =
typename Base::All_sides_oblivious_category;
typedef typename Base::X_monotone_subcurve_2 X_monotone_subcurve_2;
typedef typename Base::Size Size;
typedef typename Base::size_type size_type;
using X_monotone_subcurve_2 = typename Base::X_monotone_subcurve_2;
using Size = typename Base::Size;
using size_type = typename Base::size_type;
typedef typename Base::Point_2 Point_2;
typedef typename Base::X_monotone_curve_2 X_monotone_curve_2;
using Point_2 = typename Base::Point_2;
using X_monotone_curve_2 = typename Base::X_monotone_curve_2;
typedef typename Base::Compare_x_2 Compare_x_2;
typedef typename Base::Compare_xy_2 Compare_xy_2;
typedef typename Base::Construct_min_vertex_2 Construct_min_vertex_2;
typedef typename Base::Construct_max_vertex_2 Construct_max_vertex_2;
typedef typename Base::Is_vertical_2 Is_vertical_2;
typedef typename Base::Compare_y_at_x_2 Compare_y_at_x_2;
typedef typename Base::Compare_y_at_x_left_2 Compare_y_at_x_left_2;
typedef typename Base::Compare_y_at_x_right_2 Compare_y_at_x_right_2;
typedef typename Base::Equal_2 Equal_2;
typedef typename Base::Compare_endpoints_xy_2 Compare_endpoints_xy_2;
typedef typename Base::Construct_opposite_2 Construct_opposite_2;
typedef typename Base::Approximate_2 Approximate_2;
typedef typename Base::Construct_x_monotone_curve_2
Construct_x_monotone_curve_2;
typedef typename Base::Parameter_space_in_x_2 Parameter_space_in_x_2;
typedef typename Base::Parameter_space_in_y_2 Parameter_space_in_y_2;
typedef typename Base::Compare_x_on_boundary_2 Compare_x_on_boundary_2;
typedef typename Base::Compare_x_near_boundary_2 Compare_x_near_boundary_2;
typedef typename Base::Compare_y_on_boundary_2 Compare_y_on_boundary_2;
typedef typename Base::Compare_y_near_boundary_2 Compare_y_near_boundary_2;
typedef typename Base::Is_on_y_identification_2 Is_on_y_identification_2;
typedef typename Base::Is_on_x_identification_2 Is_on_x_identification_2;
using Compare_x_2 = typename Base::Compare_x_2;
using Compare_xy_2 = typename Base::Compare_xy_2;
using Construct_min_vertex_2 = typename Base::Construct_min_vertex_2;
using Construct_max_vertex_2 = typename Base::Construct_max_vertex_2;
using Is_vertical_2 = typename Base::Is_vertical_2;
using Compare_y_at_x_2 = typename Base::Compare_y_at_x_2;
using Compare_y_at_x_left_2 = typename Base::Compare_y_at_x_left_2;
using Compare_y_at_x_right_2 = typename Base::Compare_y_at_x_right_2;
using Equal_2 = typename Base::Equal_2;
using Compare_endpoints_xy_2 = typename Base::Compare_endpoints_xy_2;
using Construct_opposite_2 = typename Base::Construct_opposite_2;
using Approximate_2 = typename Base::Approximate_2;
using Construct_x_monotone_curve_2 =
typename Base::Construct_x_monotone_curve_2;
using Parameter_space_in_x_2 = typename Base::Parameter_space_in_x_2;
using Parameter_space_in_y_2 = typename Base::Parameter_space_in_y_2;
using Compare_x_on_boundary_2 = typename Base::Compare_x_on_boundary_2;
using Compare_x_near_boundary_2 = typename Base::Compare_x_near_boundary_2;
using Compare_y_on_boundary_2 = typename Base::Compare_y_on_boundary_2;
using Compare_y_near_boundary_2 = typename Base::Compare_y_near_boundary_2;
using Is_on_y_identification_2 = typename Base::Is_on_y_identification_2;
using Is_on_x_identification_2 = typename Base::Is_on_x_identification_2;
typedef typename Base::Trim_2 Trim_2;
using Trim_2 = typename Base::Trim_2;
//@}
/// \name Types and functors inherited from the subcurve geometry traits.
//@{
typedef typename Subcurve_traits_2::Has_merge_category Has_merge_category;
typedef typename Subcurve_traits_2::Multiplicity Multiplicity;
typedef typename Subcurve_traits_2::Curve_2 Subcurve_2;
using Has_merge_category = typename Subcurve_traits_2::Has_merge_category;
using Multiplicity = typename Subcurve_traits_2::Multiplicity;
using Subcurve_2 = typename Subcurve_traits_2::Curve_2;
//@}
// Backward compatibility:
typedef Subcurve_2 Segment_2;
using Segment_2 = Subcurve_2;
private:
typedef Arr_polycurve_traits_2<Subcurve_traits_2> Self;
using Self = Arr_polycurve_traits_2<Subcurve_traits_2>;
public:
/*! Default constructor */
@ -128,7 +125,7 @@ public:
/*! A polycurve represents a general continuous piecewise-linear
* curve, without degenerated subcurves.
*/
typedef internal::Polycurve_2<Subcurve_2, Point_2> Curve_2;
using Curve_2 = internal::Polycurve_2<Subcurve_2, Point_2>;
/// \name Basic predicate functors(based on the subcurve traits).
//@{
@ -138,8 +135,7 @@ public:
*/
class Number_of_points_2 : public Base::Number_of_points_2 {
public:
size_type operator()(const Curve_2& cv) const
{
size_type operator()(const Curve_2& cv) const {
size_type num_seg = cv.number_of_subcurves();
return (num_seg == 0) ? 0 : num_seg + 1;
}
@ -161,8 +157,9 @@ public:
//! A functor for subdividing curves into x-monotone curves.
class Make_x_monotone_2 {
protected:
typedef Arr_polycurve_traits_2<Subcurve_traits_2> Polycurve_traits_2;
/*! The traits (in case it has state) */
using Polycurve_traits_2 = Arr_polycurve_traits_2<Subcurve_traits_2>;
//! The traits (in case it has state)
const Polycurve_traits_2& m_poly_traits;
public:
@ -183,12 +180,10 @@ public:
private:
template <typename OutputIterator>
OutputIterator operator_impl(const Curve_2& cv, OutputIterator oi,
Arr_all_sides_oblivious_tag) const
{
typedef std::variant<Point_2, X_monotone_subcurve_2>
Make_x_monotone_subresult;
typedef std::variant<Point_2, X_monotone_curve_2>
Make_x_monotone_result;
Arr_all_sides_oblivious_tag) const {
using Make_x_monotone_subresult =
std::variant<Point_2, X_monotone_subcurve_2>;
using Make_x_monotone_result = std::variant<Point_2, X_monotone_curve_2>;
// If the polycurve is empty, return.
if (cv.number_of_subcurves() == 0) return oi;
@ -314,14 +309,13 @@ public:
return oi;
}
//!
template <typename OutputIterator>
OutputIterator operator_impl(const Curve_2& cv, OutputIterator oi,
Arr_not_all_sides_oblivious_tag) const
{
typedef std::variant<Point_2, X_monotone_subcurve_2>
Make_x_monotone_subresult;
typedef std::variant<Point_2, X_monotone_curve_2>
Make_x_monotone_result;
Arr_not_all_sides_oblivious_tag) const {
using Make_x_monotone_subresult =
std::variant<Point_2, X_monotone_subcurve_2>;
using Make_x_monotone_result = std::variant<Point_2, X_monotone_curve_2>;
// If the polycurve is empty, return.
if (cv.number_of_subcurves() == 0) return oi;
@ -486,7 +480,7 @@ public:
*/
class Push_back_2 : public Base::Push_back_2 {
protected:
typedef Arr_polycurve_traits_2<Subcurve_traits_2> Polycurve_traits_2;
using Polycurve_traits_2 = Arr_polycurve_traits_2<Subcurve_traits_2>;
public:
/*! Constructor. */
@ -522,7 +516,7 @@ public:
*/
class Push_front_2 : public Base::Push_front_2 {
protected:
typedef Arr_polycurve_traits_2<Subcurve_traits_2> Polycurve_traits_2;
using Polycurve_traits_2 = Arr_polycurve_traits_2<Subcurve_traits_2>;
public:
/*! Constructor. */
@ -554,8 +548,9 @@ public:
class Split_2 {
protected:
typedef Arr_polycurve_traits_2<Subcurve_traits_2> Polycurve_traits_2;
/*! The polycurve traits (in case it has state) */
using Polycurve_traits_2 = Arr_polycurve_traits_2<Subcurve_traits_2>;
//! The polycurve traits (in case it has state)
const Polycurve_traits_2& m_poly_traits;
public:
@ -571,8 +566,7 @@ public:
* \pre p lies on cv but is not one of its end-points.
*/
void operator()(const X_monotone_curve_2& xcv, const Point_2& p,
X_monotone_curve_2& xcv1, X_monotone_curve_2& xcv2) const
{
X_monotone_curve_2& xcv1, X_monotone_curve_2& xcv2) const {
const Subcurve_traits_2* geom_traits = m_poly_traits.subcurve_traits_2();
auto min_vertex = geom_traits->construct_min_vertex_2_object();
auto max_vertex = geom_traits->construct_max_vertex_2_object();
@ -670,8 +664,9 @@ public:
class Intersect_2 {
protected:
typedef Arr_polycurve_traits_2<Subcurve_traits_2> Polycurve_traits_2;
/*! The polycurve traits (in case it has state) */
using Polycurve_traits_2 = Arr_polycurve_traits_2<Subcurve_traits_2>;
//! The polycurve traits (in case it has state)
const Polycurve_traits_2& m_poly_traits;
public:
@ -689,15 +684,13 @@ public:
* \return The past-the-end iterator.
*/
template <typename OutputIterator>
OutputIterator
operator()(const X_monotone_curve_2& cv1,
const X_monotone_curve_2& cv2,
OutputIterator oi) const
{
typedef std::pair<Point_2, Multiplicity> Intersection_point;
typedef std::variant<Intersection_point, X_monotone_subcurve_2>
Intersection_base_result;
OutputIterator operator()(const X_monotone_curve_2& cv1,
const X_monotone_curve_2& cv2,
OutputIterator oi) const {
// std::cout << "intersect(" << cv1 << ", " << cv2 << ")\n";
using Intersection_point = std::pair<Point_2, Multiplicity>;
using Intersection_base_result =
std::variant<Intersection_point, X_monotone_subcurve_2>;
const Subcurve_traits_2* geom_traits = m_poly_traits.subcurve_traits_2();
auto cmp_y_at_x = m_poly_traits.compare_y_at_x_2_object();
@ -705,18 +698,21 @@ public:
auto min_vertex = geom_traits->construct_min_vertex_2_object();
auto max_vertex = geom_traits->construct_max_vertex_2_object();
auto intersect = geom_traits->intersect_2_object();
auto cmp_seg_endpts = geom_traits->compare_endpoints_xy_2_object();
auto construct_opposite = geom_traits->construct_opposite_2_object();
auto cmp_endpts = geom_traits->compare_endpoints_xy_2_object();
auto ctr_opposite = geom_traits->construct_opposite_2_object();
auto cmp_xy = m_poly_traits.compare_xy_2_object();
Comparison_result dir1 = cmp_seg_endpts(cv1[0]);
Comparison_result dir2 = cmp_seg_endpts(cv2[0]);
Comparison_result dir1 = cmp_endpts(cv1[0]);
Comparison_result dir2 = cmp_endpts(cv2[0]);
// std::cout << "dir1: " << dir1 << std::endl;
// std::cout << "dir2: " << dir2 << std::endl;
std::vector<X_monotone_subcurve_2> ocv; // Used to represent overlaps.
const bool invert_ocv = ((dir1 == LARGER) && (dir2 == LARGER));
const bool consistent = (dir1 == dir2);
#ifdef CGAL_ALWAYS_LEFT_TO_RIGHT
CGAL_assertion(consistent);
#endif
std::vector<X_monotone_subcurve_2> ocv; // Used to represent overlaps.
const bool invert_ocv = ((dir1 == LARGER) && (dir2 == LARGER));
const std::size_t n1 = cv1.number_of_subcurves();
const std::size_t n2 = cv2.number_of_subcurves();
@ -724,297 +720,295 @@ public:
std::size_t i1 = (dir1 == SMALLER) ? 0 : n1-1;
std::size_t i2 = (dir2 == SMALLER) ? 0 : n2-1;
auto compare_xy = m_poly_traits.compare_xy_2_object();
Comparison_result left_res =
compare_xy(cv1[i1], ARR_MIN_END, cv2[i2], ARR_MIN_END);
Point_2 saved_point;
X_monotone_curve_2 saved_xcv;
bool point_saved = false;
bool xcv_saved = false;
// Save a point
auto update_saved_point = [&](const Point_2& p) {
point_saved = true;
saved_point = p;
};
// Spit out the saved point
auto spit_saved_point = [&]() {
*oi++ = std::make_pair(saved_point, 0);
point_saved = false;
};
// Add a subcurve to the saved x-monotone curve
auto update_saved_xcv = [&](const X_monotone_subcurve_2& sxcv) {
xcv_saved = true;
// We maintain the invariant that if the input curves have opposite
// directions (! consistent), the overalpping curves are directed
// left=>right. This, however, is not guaranteed by all traits for the
// subcurves. (It is guaranteed by some traits, but this is
// insufficient.) Therefore, we need to enforce it. That is, we make
// sure the subcurves are also directed left=>right in this case.
if (! consistent && (cmp_endpts(sxcv) == LARGER)) {
if (invert_ocv) {
saved_xcv.push_front(ctr_opposite(sxcv));
return;
}
saved_xcv.push_back(ctr_opposite(sxcv));
return;
}
if (invert_ocv) {
saved_xcv.push_front(sxcv);
return;
}
saved_xcv.push_back(sxcv);
};
// Spit out the saved x-monotone curve
auto spit_saved_xcv = [&]() {
*oi++ = saved_xcv;
saved_xcv.clear();
xcv_saved = false;
};
auto left_res = cmp_xy(cv1[i1], ARR_MIN_END, cv2[i2], ARR_MIN_END);
if (left_res == SMALLER) {
// cv1's left endpoint is to the left of cv2's left endpoint:
// cv1's left endpoint is to the left of cv2's left endpoint.
// Locate the index i1 of the subcurve in cv1 which contains cv2's
// left endpoint.
i1 = m_poly_traits.locate_impl(cv1, cv2[i2], ARR_MIN_END,
All_sides_oblivious_category());
if (i1 == Polycurve_traits_2::INVALID_INDEX) return oi;
if (equal(max_vertex(cv1[i1]), min_vertex(cv2[i2]))) {
if (((dir1 == SMALLER) && (i1 == n1-1)) ||
((dir1 == LARGER) && (i1 == 0))){
// cv1's right endpoint equals cv2's left endpoint
// Thus we can return this single(!) intersection point
Intersection_point p(max_vertex(cv1[i1]), 0);
*oi++ = p;
return oi;
if (cmp_y_at_x(cv2[i2], ARR_MIN_END, cv1[i1]) == EQUAL) {
const auto& p = min_vertex(cv2[i2]);
update_saved_point(p);
if (equal(max_vertex(cv1[i1]), p)) {
if (dir1 == SMALLER) {
++i1;
if (i1 == n1) {
spit_saved_point();
return oi;
}
}
else {
if (i1 != 0) --i1;
else {
spit_saved_point();
return oi;
}
}
}
dir1 == SMALLER ?
++i1 :
(i1 != 0) ? --i1 : (std::size_t) Polycurve_traits_2::INVALID_INDEX;
left_res = EQUAL;
}
}
else if (left_res == LARGER) {
// cv1's left endpoint is to the right of cv2's left endpoint:
// cv1's left endpoint is to the right of cv2's left endpoint.
// Locate the index i2 of the subcurve in cv2 which contains cv1's
// left endpoint.
i2 = m_poly_traits.locate_impl(cv2, cv1[i1], ARR_MIN_END,
All_sides_oblivious_category());
if (i2 == Polycurve_traits_2::INVALID_INDEX) return oi;
if (equal(max_vertex(cv2[i2]), min_vertex(cv1[i1]))) {
if (((dir2 == SMALLER) && (i2 == n2-1)) ||
((dir2 == LARGER) && (i2 == 0))){
// cv2's right endpoint equals cv1's left endpoint
// Thus we can return this single(!) intersection point
Intersection_point p(max_vertex(cv2[i2]), 0);
*oi++ = p;
return oi;
}
dir2 == SMALLER ?
++i2 :
(i2 != 0) ? --i2 : (std::size_t) Polycurve_traits_2::INVALID_INDEX;
left_res = EQUAL;
}
}
// Check if the left endpoint lies on the other polycurve.
bool left_coincides = (left_res == EQUAL);
bool left_overlap = false;
if (left_res == SMALLER)
left_coincides = (cmp_y_at_x(cv2[i2], ARR_MIN_END, cv1[i1]) == EQUAL);
else if (left_res == LARGER)
left_coincides = (cmp_y_at_x(cv1[i1], ARR_MIN_END, cv2[i2]) == EQUAL);
// The main loop: Go simultaneously over both polycurves.
Comparison_result right_res = left_res;
bool right_coincides = left_coincides;
bool right_overlap = false;
while (((dir1 == SMALLER) && (dir2 == SMALLER) &&
(i1 < n1) && (i2 < n2)) ||
((dir1 != SMALLER) && (dir2 == SMALLER) &&
(i1 != Polycurve_traits_2::INVALID_INDEX) && (i2 < n2)) ||
((dir1 == SMALLER) && (dir2 != SMALLER) && (i1 < n1) &&
(i2 != Polycurve_traits_2::INVALID_INDEX)) ||
((dir1 != SMALLER) && (dir2 != SMALLER) &&
(i1 != Polycurve_traits_2::INVALID_INDEX) &&
(i2 != Polycurve_traits_2::INVALID_INDEX)))
{
right_res = compare_xy(cv1[i1], ARR_MAX_END, cv2[i2], ARR_MAX_END);
right_coincides = (right_res == EQUAL);
if (right_res == SMALLER)
right_coincides =
(cmp_y_at_x(cv1[i1], ARR_MAX_END, cv2[i2]) == EQUAL);
else if (right_res == LARGER)
right_coincides =
(cmp_y_at_x(cv2[i2], ARR_MAX_END, cv1[i1]) == EQUAL);
right_overlap = false;
//! EF: the following code is a bit suspicious. It may erroneously
// assume that the subcurves cannot overlap more than once.
if (! right_coincides && ! left_coincides) {
// Non of the endpoints of the current subcurve of one polycurve
// coincides with the current subcurve of the other polycurve:
// Output the intersection if exists.
std::vector<Intersection_base_result> xections;
intersect(cv1[i1], cv2[i2], std::back_inserter(xections));
for (const auto& xection : xections) {
const X_monotone_subcurve_2* subcv_p =
std::get_if<X_monotone_subcurve_2>(&xection);
if (subcv_p != nullptr) {
ocv.push_back(*subcv_p);
oi = output_ocv (ocv, invert_ocv, oi);
continue;
}
const Intersection_point* p_p =
std::get_if<Intersection_point>(&xection);
if (p_p != nullptr) *oi++ = *p_p;
}
}
else if (right_coincides && left_coincides) {
// An overlap exists between the current subcurves of the
// polycurves: Output the overlapping subcurve.
right_overlap = true;
std::vector<Intersection_base_result> sub_xections;
intersect(cv1[i1], cv2[i2], std::back_inserter(sub_xections));
for (const auto& item : sub_xections) {
const X_monotone_subcurve_2* x_seg =
std::get_if<X_monotone_subcurve_2>(&item);
if (x_seg != nullptr) {
X_monotone_subcurve_2 seg = *x_seg;
// We maintain the variant that if the input curves have opposite
// directions (! consistent), the overalpping curves are directed
// left=>right. This, however, is not guaranteed for the
// subcurves. Therefore, we need to enforce it. That is, we make
// sure the subcurves are also directed left=>right in this case.
if (! consistent && (cmp_seg_endpts(seg) == LARGER))
seg = construct_opposite(seg);
#ifdef CGAL_ALWAYS_LEFT_TO_RIGHT
CGAL_assertion(cmp_seg_endpts(seg) == SMALLER);
#endif
ocv.push_back(seg);
}
const Intersection_point* p_ptr =
std::get_if<Intersection_point>(&item);
if (p_ptr != nullptr) {
// Any point that is not equal to the max_vertex of the
// subcurve should be inserted into oi.
// The max_vertex of the current subcurve (if intersecting)
// will be taken care of as the min_vertex of in the next
// iteration.
if (! equal(p_ptr->first, max_vertex(cv1[i1])))
*oi++ = *p_ptr;
}
}
}
else if (left_coincides && ! right_coincides) {
// std::cout << "Left is coinciding but right is not." << std::endl;
// The left point of the current subcurve of one polycurve
// coincides with the current subcurve of the other polycurve.
if (left_overlap) {
// An overlap occurred at the previous iteration:
// Output the overlapping polycurve.
CGAL_assertion(ocv.size() > 0);
oi = output_ocv (ocv, invert_ocv, oi);
}
else {
// The left point of the current subcurve of one
// polycurve coincides with the current subcurve of the
// other polycurve, and no overlap occurred at the
// previous iteration: Output the intersection
// point. The derivative of at least one of the
// polycurves is not defined at this point, so we give
// it multiplicity 0.
if (left_res == SMALLER) {
Intersection_point p(min_vertex(cv2[i2]), 0);
*oi++ = p;
if (cmp_y_at_x(cv1[i1], ARR_MIN_END, cv2[i2]) == EQUAL) {
const auto& p = min_vertex(cv1[i1]);
update_saved_point(p);
if (equal(max_vertex(cv2[i2]), p)) {
if (dir2 == SMALLER) {
++i2;
if (i2 == n2) {
spit_saved_point();
return oi;
}
}
else {
Intersection_point p(min_vertex(cv1[i1]), 0);
*oi++ = p;
if (i2 != 0) --i2;
else {
spit_saved_point();
return oi;
}
}
}
}
}
else {
CGAL_assertion(left_res == EQUAL);
update_saved_point(min_vertex(cv1[i1]));
}
do {
// std::cout << "i1, i2 = " << i1 <<", " << i2 << std::endl;
std::vector<Intersection_base_result> xsecs;
intersect(cv1[i1], cv2[i2], std::back_inserter(xsecs));
// Iterate over all intersections.
if (! xsecs.empty()) {
auto it = xsecs.begin();
const auto& xsec = *it++;
if (it == xsecs.end()) {
// There is exactly one intersection
const auto* xp_p = std::get_if<Intersection_point>(&xsec);
if (xp_p) {
// The intersection is a point
const auto& xp = xp_p->first;
auto is_min_end_1 = equal(xp, min_vertex(cv1[i1]));
auto is_min_end_2 = equal(xp, min_vertex(cv2[i2]));
if ((is_min_end_1 || is_min_end_2) && (xcv_saved || point_saved)) {
// It is impossible to have a pending point and a pending
// x-monotone curve concurrently
CGAL_assertion((xcv_saved && ! point_saved) ||
(! xcv_saved && point_saved));
if (xcv_saved) spit_saved_xcv();
else spit_saved_point();
}
else {
auto is_max_end_1 = equal(xp, max_vertex(cv1[i1]));
auto is_max_end_2 = equal(xp, max_vertex(cv2[i2]));
if (is_max_end_1 || is_max_end_2) update_saved_point(xp);
else *oi++ = *xp_p;
}
}
else {
// The intersection is an x-monotone curve
auto* sxcv_p = std::get_if<X_monotone_subcurve_2>(&xsec);
CGAL_assertion(sxcv_p != nullptr);
const auto& xp_min = min_vertex(*sxcv_p);
auto is_min_end_1 = equal(xp_min, min_vertex(cv1[i1]));
auto is_min_end_2 = equal(xp_min, min_vertex(cv2[i2]));
const auto& xp_max = max_vertex(*sxcv_p);
auto is_max_end_1 = equal(xp_max, max_vertex(cv1[i1]));
auto is_max_end_2 = equal(xp_max, max_vertex(cv2[i2]));
if (is_min_end_1 || is_min_end_2) {
update_saved_xcv(*sxcv_p);
if (! is_max_end_1 && ! is_max_end_2) spit_saved_xcv();
}
else {
if (xcv_saved) spit_saved_xcv();
if (is_max_end_1 || is_max_end_2) update_saved_xcv(*sxcv_p);
else *oi++ = *sxcv_p;
}
point_saved = false;
}
}
else {
// There is more than one intersection
// Handle the first intersection
const auto* xp_p = std::get_if<Intersection_point>(&xsec);
if (xp_p) {
const auto& xp = xp_p->first;
auto is_min_end_1 = equal(xp, min_vertex(cv1[i1]));
auto is_min_end_2 = equal(xp, min_vertex(cv2[i2]));
if ((is_min_end_1 || is_min_end_2) && (xcv_saved || point_saved)) {
// It is impossible to have a pending point and a pending
// x-monotone curve concurrently
CGAL_assertion((xcv_saved && ! point_saved) ||
(! xcv_saved && point_saved));
if (xcv_saved) spit_saved_xcv();
else spit_saved_point();
}
else *oi++ = *xp_p;
}
else {
const auto* sxcv_p = std::get_if<X_monotone_subcurve_2>(&xsec);
CGAL_assertion(sxcv_p);
const auto& xp = min_vertex(*sxcv_p);
auto is_min_end_1 = equal(xp, min_vertex(cv1[i1]));
auto is_min_end_2 = equal(xp, min_vertex(cv2[i2]));
if (is_min_end_1 || is_min_end_2) {
update_saved_xcv(*sxcv_p);
spit_saved_xcv();
}
else {
if (xcv_saved) spit_saved_xcv();
else *oi++ = *sxcv_p;
}
point_saved = false;
}
// Handle all but the last intersection
for (; it != std::prev(xsecs.end()); ++it) {
const auto& xsec = *it;
xp_p = std::get_if<Intersection_point>(&xsec);
if (xp_p) *oi++ = *xp_p;
else {
const auto* sxcv_p = std::get_if<X_monotone_subcurve_2>(&xsec);
CGAL_assertion(sxcv_p);
*oi++ = *sxcv_p;
}
}
// Handle the last intersection
const auto& xsec = *it;
xp_p = std::get_if<Intersection_point>(&xsec);
if (xp_p) {
const auto& xp = xp_p->first;
auto is_max_end_1 = equal(xp, max_vertex(cv1[i1]));
auto is_max_end_2 = equal(xp, max_vertex(cv2[i2]));
if (is_max_end_1 || is_max_end_2) update_saved_point(xp);
else *oi++ = *xp_p;
}
else {
const auto* sxcv_p = std::get_if<X_monotone_subcurve_2>(&xsec);
CGAL_assertion(sxcv_p);
const auto& xp = max_vertex(*sxcv_p);
auto is_max_end_1 = equal(xp, max_vertex(cv1[i1]));
auto is_max_end_2 = equal(xp, max_vertex(cv2[i2]));
if (is_max_end_1 || is_max_end_2) update_saved_xcv(*sxcv_p);
else *oi++ = *sxcv_p;
}
}
}
// Proceed forward.
if (right_res != SMALLER) {
if (dir2 == SMALLER) ++i2;
else {
if (i2 == 0) i2 = Polycurve_traits_2::INVALID_INDEX;
else --i2;
}
}
// Advance the indices
auto right_res = cmp_xy(cv1[i1], ARR_MAX_END, cv2[i2], ARR_MAX_END);
if (right_res != LARGER) {
if (dir1 == SMALLER)
if (dir1 == SMALLER) {
++i1;
if (i1 == n1) break;
}
else {
if (i1 == 0) i1 = Polycurve_traits_2::INVALID_INDEX;
else --i1;
if (i1 != 0) --i1;
else break;
}
}
left_res = (right_res == SMALLER) ? LARGER :
(right_res == LARGER) ? SMALLER : EQUAL;
if (right_res != SMALLER) {
if (dir2 == SMALLER) {
++i2;
if (i2 == n2) break;
}
else {
if (i2 != 0) --i2;
else break;
}
}
} while (true);
left_coincides = right_coincides;
left_overlap = right_overlap;
} // END of while loop
// Output the remaining overlapping polycurve, if necessary.
if (ocv.size() > 0) {
oi = output_ocv (ocv, invert_ocv, oi);
}
else if (right_coincides) {
typedef std::pair<Point_2,Multiplicity> return_point;
return_point ip;
if (right_res == SMALLER) {
ip = (dir1 == SMALLER) ?
return_point(max_vertex(cv1[i1-1]), 0) :
(i1 != Polycurve_traits_2::INVALID_INDEX) ?
return_point(max_vertex(cv1[i1+1]), 0) :
return_point(max_vertex(cv1[0]), 0);
*oi++ = ip;
}
else if (right_res == LARGER) {
ip = (dir2 == SMALLER) ?
return_point(max_vertex(cv2[i2-1]), 0) :
(i2 != Polycurve_traits_2::INVALID_INDEX) ?
return_point(max_vertex(cv2[i2+1]), 0) :
return_point(max_vertex(cv2[0]), 0);
*oi++ = ip;
}
else if (((i1 > 0) && (dir1 == SMALLER)) ||
((i1 < n1) && (dir1 != SMALLER)) ||
((i1 == Polycurve_traits_2::INVALID_INDEX) &&
(dir1 != SMALLER)))
{
ip = (dir1 == SMALLER) ?
return_point(max_vertex(cv1[i1-1]), 0) :
(i1 != Polycurve_traits_2::INVALID_INDEX) ?
return_point(max_vertex(cv1[i1+1]), 0) :
return_point(max_vertex(cv1[0]), 0);
*oi++ = ip;
}
else {
CGAL_assertion_msg((dir2 == SMALLER && i2 > 0) ||
(dir2 != SMALLER && i2 < n2) ||
(dir2 != SMALLER &&
((i1 == Polycurve_traits_2::INVALID_INDEX) ||
(i2 == Polycurve_traits_2::INVALID_INDEX))),
"Wrong index for xcv2 in Intersect_2 of "
"polycurves.");
ip = (dir2 == SMALLER) ?
return_point(max_vertex(cv2[i2-1]), 0) :
(i2 != Polycurve_traits_2::INVALID_INDEX) ?
return_point(max_vertex(cv2[i2+1]), 0) :
return_point(max_vertex(cv2[0]), 0);
*oi++ = ip;
}
}
if (xcv_saved) spit_saved_xcv();
else if (point_saved) spit_saved_point();
return oi;
}
private:
template <typename OutputIterator>
inline OutputIterator output_ocv
(std::vector<X_monotone_subcurve_2>& ocv, bool invert_ocv, OutputIterator oi) const
{
inline OutputIterator output_ocv(std::vector<X_monotone_subcurve_2>& ocv,
bool invert_ocv, OutputIterator oi) const {
X_monotone_curve_2 curve;
if (invert_ocv)
std::reverse (ocv.begin(), ocv.end());
for (X_monotone_subcurve_2& sc : ocv)
curve.push_back (sc);
*(oi ++) = curve;
if (invert_ocv) std::reverse(ocv.begin(), ocv.end());
for (X_monotone_subcurve_2& sc : ocv) curve.push_back(sc);
*oi++ = curve;
ocv.clear();
return oi;
}
};
/*! Obtain an Intersect_2 functor object. */
Intersect_2 intersect_2_object() const
{ return Intersect_2(*this); }
Intersect_2 intersect_2_object() const { return Intersect_2(*this); }
class Are_mergeable_2 {
protected:
typedef Arr_polycurve_traits_2<Subcurve_traits_2> Polycurve_traits_2;
/*! The polycurve traits (in case it has state) */
using Polycurve_traits_2 = Arr_polycurve_traits_2<Subcurve_traits_2>;
//! The polycurve traits (in case it has state)
const Polycurve_traits_2& m_poly_traits;
public:
/*! Constructor. */
Are_mergeable_2(const Polycurve_traits_2& traits) :
m_poly_traits(traits)
{}
Are_mergeable_2(const Polycurve_traits_2& traits) : m_poly_traits(traits) {}
/*! Check whether it is possible to merge two given x-monotone curves.
* \param cv1 The first curve.
@ -1023,25 +1017,17 @@ public:
* common endpoint and the same orientation;(false) otherwise.
*/
bool operator()(const X_monotone_curve_2& cv1,
const X_monotone_curve_2& cv2) const
{
const X_monotone_curve_2& cv2) const {
const Subcurve_traits_2* geom_traits = m_poly_traits.subcurve_traits_2();
Construct_min_vertex_2 min_vertex =
m_poly_traits.construct_min_vertex_2_object();
Construct_max_vertex_2 max_vertex =
m_poly_traits.construct_max_vertex_2_object();
typename Subcurve_traits_2::Equal_2 equal =
geom_traits->equal_2_object();
typename Subcurve_traits_2::Is_vertical_2 is_seg_vertical =
geom_traits->is_vertical_2_object();
auto min_vertex = m_poly_traits.construct_min_vertex_2_object();
auto max_vertex = m_poly_traits.construct_max_vertex_2_object();
auto equal = geom_traits->equal_2_object();
auto is_seg_vertical = geom_traits->is_vertical_2_object();
Comparison_result dir1 =
m_poly_traits.compare_endpoints_xy_2_object()(cv1);
Comparison_result dir2 =
m_poly_traits.compare_endpoints_xy_2_object()(cv2);
auto dir1 = m_poly_traits.compare_endpoints_xy_2_object()(cv1);
auto dir2 = m_poly_traits.compare_endpoints_xy_2_object()(cv2);
if (dir1 != dir2)
return false;
if (dir1 != dir2) return false;
bool ver1 = is_seg_vertical(cv1[0]);
bool ver2 = is_seg_vertical(cv2[0]);
@ -1073,8 +1059,9 @@ public:
*/
class Merge_2 {
protected:
typedef Arr_polycurve_traits_2<Subcurve_traits_2> Geometry_traits;
/*! The traits (in case it has state) */
using Geometry_traits = Arr_polycurve_traits_2<Subcurve_traits_2>;
//! The traits (in case it has state)
const Geometry_traits& m_poly_traits;
public:
@ -1091,25 +1078,21 @@ public:
*/
void operator()(const X_monotone_curve_2& cv1,
const X_monotone_curve_2& cv2,
X_monotone_curve_2& c) const
{
X_monotone_curve_2& c) const {
CGAL_precondition(m_poly_traits.are_mergeable_2_object()(cv1, cv2));
Construct_min_vertex_2 get_min_v =
m_poly_traits.construct_min_vertex_2_object();
Construct_max_vertex_2 get_max_v =
m_poly_traits.construct_max_vertex_2_object();
Compare_endpoints_xy_2 cmp_seg_endpts =
m_poly_traits.compare_endpoints_xy_2_object();
auto min_vertex = m_poly_traits.construct_min_vertex_2_object();
auto max_vertex = m_poly_traits.construct_max_vertex_2_object();
auto cmp_endpts = m_poly_traits.compare_endpoints_xy_2_object();
Equal_2 equal = m_poly_traits.equal_2_object();
c.clear();
if (// Either both are left-to-right and cv2 is to the right of cv1
((cmp_seg_endpts(cv1)==SMALLER) &&
(equal(get_max_v(cv1),get_min_v(cv2)))) ||
((cmp_endpts(cv1)==SMALLER) &&
(equal(max_vertex(cv1), min_vertex(cv2)))) ||
// or both are right-to-left and cv2 is to the left of cv1
((cmp_seg_endpts(cv1)==LARGER) &&
(equal(get_min_v(cv1), get_max_v(cv2)))))
((cmp_endpts(cv1)==LARGER) &&
(equal(min_vertex(cv1), max_vertex(cv2)))))
{
const std::size_t n1 = cv1.number_of_subcurves();
const std::size_t n2 = cv2.number_of_subcurves();
@ -1148,8 +1131,9 @@ public:
*/
class Construct_curve_2 {
protected:
typedef Arr_polycurve_traits_2<Subcurve_traits_2> Polycurve_traits_2;
/*! The polycurve traits (in case it has state) */
using Polycurve_traits_2 = Arr_polycurve_traits_2<Subcurve_traits_2>;
//! The polycurve traits (in case it has state)
const Polycurve_traits_2& m_poly_traits;
public:
@ -1165,10 +1149,9 @@ public:
* `SubcurveTraits::Point_2` or `SubcurveTraits::Subcurve_2`.
*/
template <typename ForwardIterator>
Curve_2 operator()(ForwardIterator begin, ForwardIterator end) const
{
typedef typename std::iterator_traits<ForwardIterator>::value_type VT;
typedef typename std::is_same<VT, Point_2>::type Is_point;
Curve_2 operator()(ForwardIterator begin, ForwardIterator end) const {
using VT = typename std::iterator_traits<ForwardIterator>::value_type;
using Is_point = typename std::is_same<VT, Point_2>::type;
// Dispatch the range to the appropriate implementation.
return constructor_impl(begin, end, Is_point());
}
@ -1196,8 +1179,7 @@ public:
*/
template <typename ForwardIterator>
Curve_2 constructor_impl(ForwardIterator begin, ForwardIterator end,
std::false_type) const
{
std::false_type) const {
// Range has to contain at least one subcurve
CGAL_precondition(begin != end);
return Curve_2(begin, end);

155
Kinetic_space_partition/include/CGAL/KSP/utils.h
View File

@ -29,17 +29,11 @@
#include <map>
// CGAL includes.
#include <CGAL/Bbox_3.h>
#include <CGAL/centroid.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/Iterator_range.h>
#include <CGAL/convex_hull_2.h>
#include <CGAL/number_utils.h>
#include <CGAL/assertions.h>
#include <CGAL/Cartesian_converter.h>
#include <CGAL/linear_least_squares_fitting_2.h>
#include <CGAL/linear_least_squares_fitting_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
// Boost includes.
@ -70,22 +64,6 @@ decltype(auto) distance(const Point_d& p, const Point_d& q) {
return static_cast<FT>(CGAL::approximate_sqrt(sq_dist));
}
// Project 3D point onto 2D plane.
template<typename Point_3>
typename Kernel_traits<Point_3>::Kernel::Point_2
point_2_from_point_3(const Point_3& point_3) {
return typename Kernel_traits<Point_3>::Kernel::Point_2(
point_3.x(), point_3.y());
}
// Get 3D point from a 2D point.
template<typename Point_2>
typename Kernel_traits<Point_2>::Kernel::Point_3
point_3_from_point_2(const Point_2& point_2) {
return typename Kernel_traits<Point_2>::Kernel::Point_3(
point_2.x(), point_2.y(), typename Kernel_traits<Point_2>::Kernel::FT(0));
}
// Normalize vector.
template<typename Vector_d>
inline const Vector_d normalize(const Vector_d& v) {
@ -96,7 +74,6 @@ inline const Vector_d normalize(const Vector_d& v) {
return v / static_cast<FT>(CGAL::approximate_sqrt(dot_product));
}
// Intersections. Used only in the 2D version.
// For the 3D version, see conversions.h!
template<typename Type1, typename Type2, typename ResultType>
@ -121,40 +98,6 @@ inline const ResultType intersection(const Type1& t1, const Type2& t2) {
return out;
}
// Get boundary points from a set of points.
template<typename Point_2, typename Line_2>
void boundary_points_on_line_2(
const std::vector<Point_2>& input_range,
const std::vector<std::size_t>& indices,
const Line_2& line, Point_2& p, Point_2& q) {
using Traits = typename Kernel_traits<Point_2>::Kernel;
using FT = typename Traits::FT;
using Vector_2 = typename Traits::Vector_2;
FT min_proj_value = (std::numeric_limits<FT>::max)();
FT max_proj_value = -min_proj_value;
const auto ref_vector = line.to_vector();
const auto& ref_point = input_range[indices.front()];
for (const std::size_t index : indices) {
const auto& query = input_range[index];
const auto point = line.projection(query);
const Vector_2 curr_vector(ref_point, point);
const FT value = CGAL::scalar_product(curr_vector, ref_vector);
if (value < min_proj_value) {
min_proj_value = value;
p = point;
}
if (value > max_proj_value) {
max_proj_value = value;
q = point;
}
}
}
// Angles.
// Converts radians to degrees.
@ -202,104 +145,6 @@ angle_2(const Direction_2& dir1, const Direction_2& dir2) {
return CGAL::abs(convert_angle_2(angle_2));
}
// Classes.
template<typename IVertex>
class Indexer {
public:
std::size_t operator()(const IVertex& ivertex) {
const auto pair = m_indices.insert(
std::make_pair(ivertex, m_indices.size()));
const auto& item = pair.first;
const std::size_t idx = item->second;
return idx;
}
void clear() { m_indices.clear(); }
private:
std::map<IVertex, std::size_t> m_indices;
};
template<
typename GeomTraits,
typename InputRange,
typename NeighborQuery>
class Estimate_normals_2 {
public:
using Traits = GeomTraits;
using Input_range = InputRange;
using Neighbor_query = NeighborQuery;
using Kernel = Traits;
using FT = typename Kernel::FT;
using Vector_2 = typename Kernel::Vector_2;
using Line_2 = typename Kernel::Line_2;
using Indices = std::vector<std::size_t>;
using IK = CGAL::Exact_predicates_inexact_constructions_kernel;
using IPoint_2 = typename IK::Point_2;
using ILine_2 = typename IK::Line_2;
using Converter = CGAL::Cartesian_converter<Kernel, IK>;
Estimate_normals_2(
const Input_range& input_range,
const Neighbor_query& neighbor_query) :
m_input_range(input_range),
m_neighbor_query(neighbor_query) {
CGAL_precondition(input_range.size() > 0);
}
void get_normals(std::vector<Vector_2>& normals) const {
normals.clear();
normals.reserve(m_input_range.size());
Indices neighbors;
for (std::size_t i = 0; i < m_input_range.size(); ++i) {
neighbors.clear();
m_neighbor_query(i, neighbors);
const auto line = fit_line(neighbors);
auto normal = line.to_vector();
normal = normal.perpendicular(CGAL::COUNTERCLOCKWISE);
normal = normalize(normal);
normals.push_back(normal);
}
CGAL_assertion(normals.size() == m_input_range.size());
}
private:
const Input_range& m_input_range;
const Neighbor_query& m_neighbor_query;
const Converter m_converter;
const Line_2 fit_line(const Indices& indices) const {
CGAL_assertion(indices.size() > 0);
std::vector<IPoint_2> points;
points.reserve(indices.size());
for (const std::size_t index : indices) {
const auto& point = get(m_neighbor_query.point_map(), index);
points.push_back(m_converter(point));
}
CGAL_assertion(points.size() == indices.size());
ILine_2 fitted_line;
IPoint_2 fitted_centroid;
CGAL::linear_least_squares_fitting_2(
points.begin(), points.end(),
fitted_line, fitted_centroid,
CGAL::Dimension_tag<0>());
const Line_2 line(
static_cast<FT>(fitted_line.a()),
static_cast<FT>(fitted_line.b()),
static_cast<FT>(fitted_line.c()));
return line;
}
};
#endif
} // namespace internal

5
Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h
View File

@ -71,8 +71,6 @@ public:
using Point_3 = typename Kernel::Point_3;
using Index = std::pair<std::size_t, std::size_t>;
/*!
identifies the support of a face in the exported linear cell complex, which is either an input polygon, identified by a non-negative number, a side of the bounding box in the rotated coordinate system or a face of the octree used to partition the scene.
*/
@ -126,6 +124,8 @@ private:
using Triangle_2 = typename Kernel::Triangle_2;
using Transform_3 = CGAL::Aff_transformation_3<Kernel>;
using Index = std::pair<std::size_t, std::size_t>;
using Data_structure = KSP_3::internal::Data_structure<Kernel, Intersection_kernel>;
using IVertex = typename Data_structure::IVertex;
@ -730,6 +730,7 @@ public:
m_partition_nodes[i].m_data->face_to_volumes().clear();
}
if (m_parameters.verbose)
std::cout << "ksp v: " << m_partition_nodes[0].m_data->vertices().size() << " f: " << m_partition_nodes[0].face2vertices.size() << " vol: " << m_volumes.size() << std::endl;
return;

2
Kinetic_surface_reconstruction/examples/Kinetic_surface_reconstruction/ksr_building.cpp
View File

@ -39,7 +39,7 @@ int main(const int, const char**) {
std::vector<Point_3> vtx;
std::vector<std::vector<std::size_t> > polylist;
std::vector<FT> lambdas{0.3, 0.5, 0.7, 0.8, 0.9, 0.95, 0.99};
std::vector<FT> lambdas{0.5, 0.7, 0.8, 0.9};
bool non_empty = false;

11
Kinetic_surface_reconstruction/examples/Kinetic_surface_reconstruction/ksr_parameters.cpp
View File

@ -167,13 +167,6 @@ int main(const int argc, const char** argv) {
// Algorithm.
KSR ksr(point_set, param);
FT max_d, max_dev;
std::size_t num;
ksr.estimate_detection_parameters(max_d, max_dev, num);
std::cout << "d: " << max_d << std::endl;
std::cout << "dev: " << max_dev << std::endl;
std::cout << "num: " << num << std::endl;
Timer timer;
timer.start();
std::size_t num_shapes = ksr.detect_planar_shapes(param);
@ -206,7 +199,7 @@ int main(const int argc, const char** argv) {
FT after_reconstruction = timer.time();
if (polylist.size() > 0)
CGAL::IO::write_polygon_soup("building_c_" + std::to_string(parameters.graphcut_lambda) + (parameters.use_ground ? "_g" : "_") + ".off", vtx, polylist);
CGAL::IO::write_polygon_soup("polylist_" + std::to_string(parameters.graphcut_lambda) + (parameters.use_ground ? "_g" : "_") + ".off", vtx, polylist);
timer.stop();
const FT time = static_cast<FT>(timer.time());
@ -230,7 +223,7 @@ int main(const int argc, const char** argv) {
if (polylist.size() > 0) {
non_empty = true;
CGAL::IO::write_polygon_soup("building_c_" + std::to_string(l) + (parameters.use_ground ? "_g" : "_") + ".off", vtx, polylist);
CGAL::IO::write_polygon_soup("polylist_" + std::to_string(l) + (parameters.use_ground ? "_g" : "_") + ".off", vtx, polylist);
}
}

2
Kinetic_surface_reconstruction/include/CGAL/KSR_3/Graphcut.h
View File

@ -21,8 +21,6 @@
#include <CGAL/boost/graph/alpha_expansion_graphcut.h>
#include <CGAL/boost/graph/Alpha_expansion_MaxFlow_tag.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Cartesian_converter.h>
// Internal includes.

103
Kinetic_surface_reconstruction/include/CGAL/Kinetic_surface_reconstruction_3.h
View File

@ -24,7 +24,6 @@
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Point_set.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/KSP/debug.h>
@ -562,8 +561,6 @@ private:
using Tds = CGAL::Triangulation_data_structure_2<Vbi, Fbi>;
using Delaunay_2 = CGAL::Delaunay_triangulation_2<Kernel, Tds>;
using Delaunay_3 = CGAL::Delaunay_triangulation_3<Kernel>;
typedef CGAL::Linear_cell_complex_traits<3, CGAL::Exact_predicates_exact_constructions_kernel> Traits;
using LCC = CGAL::Linear_cell_complex_for_combinatorial_map<3, 3, Traits, typename KSP::Linear_cell_complex_min_items>;
using Dart_descriptor = typename LCC::Dart_descriptor;
@ -588,8 +585,6 @@ private:
typedef typename CDT::Finite_edges_iterator Finite_edges_iterator;
typedef typename CDT::Finite_faces_iterator Finite_faces_iterator;
//using Visibility = KSR_3::Visibility<Kernel, Intersection_kernel, Point_map, Normal_map>;
using Index = typename KSP::Index;
using Face_attribute = typename LCC::Base::template Attribute_descriptor<2>::type;
using Volume_attribute = typename LCC::Base::template Attribute_descriptor<3>::type;
@ -1716,6 +1711,74 @@ private:
m_face_area_lcc[i] = m_face_area_lcc[i] * 2.0 * m_total_inliers / total_area;
}
FT area(typename LCC::Dart_descriptor face_dart, Plane_3 &pl, std::vector<typename Kernel::Triangle_3> *tris = nullptr) {
std::vector<Point_3> face;
From_exact from_exact;
Dart_descriptor n = face_dart;
do {
face.push_back(from_exact(m_lcc.point(n)));
n = m_lcc.beta(n, 0);
} while (n != face_dart);
pl = from_exact(m_lcc.template info<2>(face_dart).plane);
Delaunay_2 tri;
for (const Point_3& p : face)
tri.insert(pl.to_2d(p));
FT face_area = 0;
// Get area
for (auto fit = tri.finite_faces_begin(); fit != tri.finite_faces_end(); ++fit) {
const typename Kernel::Triangle_2 triangle(
fit->vertex(0)->point(),
fit->vertex(1)->point(),
fit->vertex(2)->point());
face_area += triangle.area();
if (tris)
tris->push_back(typename Kernel::Triangle_3(pl.to_3d(fit->vertex(0)->point()), pl.to_3d(fit->vertex(1)->point()), pl.to_3d(fit->vertex(2)->point())));
}
return face_area;
}
FT volume(typename LCC::Dart_descriptor volume_dart) {
FT x = 0, y = 0, z = 0;
std::size_t count = 0;
From_exact from_exact;
// Collect vertices to obtain point on the inside.
for (auto& fd : m_lcc.template one_dart_per_incident_cell<2, 3>(volume_dart)) {
typename LCC::Dart_descriptor fdh = m_lcc.dart_descriptor(fd);
for (const auto& vd : m_lcc.template one_dart_per_incident_cell<0, 2>(fdh)) {
const auto &p = from_exact(m_lcc.point(m_lcc.dart_descriptor(vd)));
x += p.x();
y += p.y();
z += p.z();
count++;
}
}
Point_3 center(x / count, y / count, z / count);
FT vol = 0;
// Second iteration for computing the area of each face and the volume spanned with the center point.
for (auto& fd : m_lcc.template one_dart_per_incident_cell<2, 3>(volume_dart)) {
typename LCC::Dart_descriptor fdh = m_lcc.dart_descriptor(fd);
Plane_3 plane;
FT a = area(fdh, plane);
Vector_3 n = plane.orthogonal_vector();
FT distance = CGAL::abs((plane.point() - center) * n);
vol += distance * a / 3.0;
}
return vol;
}
void count_volume_votes_lcc() {
// const int debug_volume = -1;
FT total_volume = 0;
@ -1749,8 +1812,6 @@ private:
const auto& point = get(m_point_map, p);
const auto& normal = get(m_normal_map, p);
// count_points++;
for (std::size_t j = 0; j < idx; j++) {
const Vector_3 vec(point, c[j]);
const FT dot_product = vec * normal;
@ -1769,30 +1830,11 @@ private:
}
}
for (auto &d : m_lcc.template one_dart_per_cell<3>()) {
for (auto& d : m_lcc.template one_dart_per_cell<3>()) {
typename LCC::Dart_descriptor dh = m_lcc.dart_descriptor(d);
std::vector<Point_3> volume_vertices;
std::size_t volume_index = m_lcc.template info<3>(dh).volume_id;
// Collect all vertices of volume to calculate volume
for (auto &fd : m_lcc.template one_dart_per_incident_cell<2, 3>(dh)) {
typename LCC::Dart_descriptor fdh = m_lcc.dart_descriptor(fd);
for (const auto &vd : m_lcc.template one_dart_per_incident_cell<0, 2>(fdh))
volume_vertices.push_back(from_exact(m_lcc.point(m_lcc.dart_descriptor(vd))));
}
Delaunay_3 tri;
for (const Point_3& p : volume_vertices)
tri.insert(p);
m_volumes[volume_index] = FT(0);
for (auto cit = tri.finite_cells_begin(); cit != tri.finite_cells_end(); ++cit) {
const auto& tet = tri.tetrahedron(cit);
m_volumes[volume_index] += tet.volume();
}
m_volumes[volume_index] = volume(dh);
total_volume += m_volumes[volume_index];
}
@ -1800,9 +1842,6 @@ private:
// Normalize volumes
for (FT& v : m_volumes)
v /= total_volume;
// for (std::size_t i = 0; i < m_volumes.size(); i++)
// std::cout << i << ": " << m_cost_matrix[0][i] << " o: " << m_cost_matrix[1][i] << " v: " << m_volumes[i] << std::endl;
}
template<typename NamedParameters>
@ -2007,8 +2046,6 @@ private:
}
void map_points_to_faces(const std::size_t polygon_index, const std::vector<Point_3>& pts, std::vector<std::pair<typename LCC::Dart_descriptor, std::vector<std::size_t> > >& face_to_points) {
std::vector<Index> faces;
if (polygon_index >= m_kinetic_partition.input_planes().size())
CGAL_assertion(false);

2
Kinetic_surface_reconstruction/package_info/Kinetic_surface_reconstruction/dependencies
View File

@ -44,6 +44,4 @@ Stream_support
Surface_mesh
Surface_mesh_segmentation
TDS_2
TDS_3
Triangulation_2
Triangulation_3