diff --git a/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_collision_detector_2.h b/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_collision_detector_2.h index 9a9a1a58d5e..36cacd84772 100644 --- a/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_collision_detector_2.h +++ b/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_collision_detector_2.h @@ -9,40 +9,48 @@ namespace CGAL { // Tests whether two polygons P and Q overlap for different translations of Q. template -class AABB_collision_detector_2 { +class AABB_collision_detector_2 +{ public: - typedef typename Kernel_::Point_2 Point_2; - typedef typename Kernel_::Vector_2 Vector_2; - typedef typename CGAL::Polygon_2 Polygon_2; - typedef typename Polygon_2::Edge_const_iterator Edge_iterator; - typedef AABB_segment_2_primitive Tree_segment_2; - typedef AABB_traits_2 Tree_traits; - typedef AABB_tree Tree_2; + typedef typename Kernel_::Point_2 Point_2; + typedef typename Kernel_::Vector_2 Vector_2; + typedef typename CGAL::Polygon_2 Polygon_2; + typedef typename Polygon_2::Edge_const_iterator Edge_iterator; + typedef AABB_segment_2_primitive + Tree_segment_2; + typedef AABB_traits_2 Tree_traits; + typedef AABB_tree Tree_2; public: - AABB_collision_detector_2(const Polygon_2 &p, const Polygon_2 &q) - : m_stationary_tree((p.edges_begin()), (p.edges_end())), m_translating_tree((q.edges_begin()), (q.edges_end())), m_p(q), m_q(p) { + AABB_collision_detector_2(const Polygon_2 &p, const Polygon_2 &q) + : m_stationary_tree((p.edges_begin()), (p.edges_end())), + m_translating_tree((q.edges_begin()), (q.edges_end())), m_p(q), m_q(p) + { + } + + // Returns true iff the polygons' boundaries intersect or one polygon is completely inside of the other one. Q is translated by t. + bool check_collision(const Point_2 &t) + { + if (m_stationary_tree.do_intersect(m_translating_tree, t)) + { + return true; } - // Returns true iff the polygons' boundaries intersect or one polygon is completely inside of the other one. Q is translated by t. - bool check_collision(const Point_2 &t) { - if (m_stationary_tree.do_intersect(m_translating_tree, t)) { - return true; - } - - Polygon_2 translated_p = transform(typename Kernel_::Aff_transformation_2(Translation(), Vector_2(ORIGIN, t)), m_p); - return (translated_p.has_on_bounded_side(*(m_q.vertices_begin())) || m_q.has_on_bounded_side(*(translated_p.vertices_begin()))); - } + Polygon_2 translated_p = transform(typename Kernel_::Aff_transformation_2( + Translation(), Vector_2(ORIGIN, t)), m_p); + return (translated_p.has_on_bounded_side(*(m_q.vertices_begin())) + || m_q.has_on_bounded_side(*(translated_p.vertices_begin()))); + } private: - Tree_2 m_stationary_tree; - Tree_2 m_translating_tree; - const Polygon_2 &m_p; - const Polygon_2 &m_q; + Tree_2 m_stationary_tree; + Tree_2 m_translating_tree; + const Polygon_2 &m_p; + const Polygon_2 &m_q; }; } // namespace CGAL diff --git a/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_segment_2_primitive.h b/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_segment_2_primitive.h index 3d93fbbceee..2593db39b06 100644 --- a/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_segment_2_primitive.h +++ b/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_segment_2_primitive.h @@ -5,40 +5,46 @@ namespace CGAL { // Wraps around a Segment_2 and provides its iterator as Id template -class AABB_segment_2_primitive { +class AABB_segment_2_primitive +{ public: - typedef Iterator_ Id; - typedef typename GeomTraits::Segment_2 Datum; - typedef typename GeomTraits::Point_2 Point; - typedef ContainerType Container; + typedef Iterator_ Id; + typedef typename GeomTraits::Segment_2 Datum; + typedef typename GeomTraits::Point_2 Point; + typedef ContainerType Container; - AABB_segment_2_primitive() {} + AABB_segment_2_primitive() {} - AABB_segment_2_primitive(Id it) : m_it(it) { - } + AABB_segment_2_primitive(Id it) : m_it(it) + { + } - AABB_segment_2_primitive(const AABB_segment_2_primitive &primitive) { - m_it = primitive.id(); - } + AABB_segment_2_primitive(const AABB_segment_2_primitive &primitive) + { + m_it = primitive.id(); + } - const Id &id() const { - return m_it; - } + const Id &id() const + { + return m_it; + } - const Datum datum() const { - return *m_it; - } + const Datum datum() const + { + return *m_it; + } - // Return a point on the primitive - Point reference_point() const { - return datum().source(); - } + // Return a point on the primitive + Point reference_point() const + { + return datum().source(); + } private: - Id m_it; + Id m_it; }; diff --git a/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_traits_2.h b/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_traits_2.h index 9db860c143a..7ed7c4a7a2e 100644 --- a/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_traits_2.h +++ b/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/AABB_traits_2.h @@ -4,164 +4,200 @@ namespace CGAL { template -class AABB_traits_2 { +class AABB_traits_2 +{ public: - typedef AABB_primitive_ Primitive; - typedef typename Primitive::Id Id; - typedef typename Primitive::Datum Datum; - typedef typename Primitive::Container Container; + typedef AABB_primitive_ Primitive; + typedef typename Primitive::Id Id; + typedef typename Primitive::Datum Datum; + typedef typename Primitive::Container Container; - typedef typename GeomTraits::Point_2 Point; - typedef typename GeomTraits::Vector_2 Vector_2; - typedef typename CGAL::Bbox_2 Bounding_box; + typedef typename GeomTraits::Point_2 Point; + typedef typename GeomTraits::Vector_2 Vector_2; + typedef typename CGAL::Bbox_2 Bounding_box; - typedef typename std::pair Object_and_primitive_id; - typedef typename std::pair Point_and_primitive_id; + typedef typename std::pair Object_and_primitive_id; + typedef typename std::pair Point_and_primitive_id; - // Types for AABB_tree - typedef typename GeomTraits::FT FT; - typedef typename GeomTraits::Point_2 Point_3; - typedef typename GeomTraits::Circle_2 Sphere_3; - typedef typename GeomTraits::Iso_rectangle_2 Iso_cuboid_3; - typedef typename GeomTraits::Construct_center_2 Construct_center_3; - typedef typename GeomTraits::Construct_iso_rectangle_2 Construct_iso_cuboid_3; - typedef typename GeomTraits::Construct_min_vertex_2 Construct_min_vertex_3; - typedef typename GeomTraits::Construct_max_vertex_2 Construct_max_vertex_3; - typedef typename GeomTraits::Compute_squared_radius_2 Compute_squared_radius_3; - typedef typename GeomTraits::Cartesian_const_iterator_2 Cartesian_const_iterator_3; - typedef typename GeomTraits::Construct_cartesian_const_iterator_2 Construct_cartesian_const_iterator_3; + // Types for AABB_tree + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Point_2 Point_3; + typedef typename GeomTraits::Circle_2 Sphere_3; + typedef typename GeomTraits::Iso_rectangle_2 Iso_cuboid_3; + typedef typename GeomTraits::Construct_center_2 Construct_center_3; + typedef typename GeomTraits::Construct_iso_rectangle_2 Construct_iso_cuboid_3; + typedef typename GeomTraits::Construct_min_vertex_2 Construct_min_vertex_3; + typedef typename GeomTraits::Construct_max_vertex_2 Construct_max_vertex_3; + typedef typename GeomTraits::Compute_squared_radius_2 Compute_squared_radius_3; + typedef typename GeomTraits::Cartesian_const_iterator_2 + Cartesian_const_iterator_3; + typedef typename GeomTraits::Construct_cartesian_const_iterator_2 + Construct_cartesian_const_iterator_3; - AABB_traits_2(const Point &translation_point): m_translation_point(translation_point) { - m_interval_x = Interval_nt(to_interval(translation_point.x())); - m_interval_y = Interval_nt(to_interval(translation_point.y())); - }; + AABB_traits_2(const Point &translation_point): m_translation_point( + translation_point) + { + m_interval_x = Interval_nt(to_interval(translation_point.x())); + m_interval_y = Interval_nt(to_interval(translation_point.y())); + }; - AABB_traits_2() { - m_translation_point = Point(0, 0); - m_interval_x = Interval_nt(0); - m_interval_y = Interval_nt(0); - }; + AABB_traits_2() + { + m_translation_point = Point(0, 0); + m_interval_x = Interval_nt(0); + m_interval_y = Interval_nt(0); + }; - Interval_nt get_interval_x() const { - return m_interval_x; + Interval_nt get_interval_x() const + { + return m_interval_x; + } + + Interval_nt get_interval_y() const + { + return m_interval_y; + } + + Point get_translation_point() const + { + return m_translation_point; + } + + // Put the n/2 smallest primitives in the front, the n/2 largest primitives + // in the back. They are compared along the bbox' longest axis. + class Sort_primitives + { + public: + template + void operator()(PrimitiveIterator first, + PrimitiveIterator beyond, + const Bounding_box &bbox) const + { + PrimitiveIterator middle = first + (beyond - first) / 2; + + if (bbox.xmax()-bbox.xmin() >= bbox.ymax()-bbox.ymin()) + { + std::nth_element(first, middle, beyond, less_x); // sort along x + } + else + { + std::nth_element(first, middle, beyond, less_y); // sort along y + } + } + }; + + Sort_primitives sort_primitives_object() const + { + return Sort_primitives(); + } + + // Computes the bounding box of a set of primitives + class Compute_bbox + { + public: + template + Bounding_box operator()(ConstPrimitiveIterator first, + ConstPrimitiveIterator beyond) const + { + Bounding_box bbox = first->datum().bbox(); + + for (++first; first != beyond; ++first) + { + bbox = bbox + first->datum().bbox(); + } + + return bbox; + } + }; + + Compute_bbox compute_bbox_object() const + { + return Compute_bbox(); + } + + class Do_intersect + { + + private: + + AABB_traits_2 *m_traits; + + public: + + Do_intersect(AABB_traits_2 *_traits): m_traits(_traits) {} + + bool operator()(const Bounding_box &q, const Bounding_box &bbox) const + { + Interval_nt x1 = Interval_nt(q.xmin(), q.xmax()); + Interval_nt y1 = Interval_nt(q.ymin(), q.ymax()); + Interval_nt x2 = Interval_nt(bbox.xmin(), + bbox.xmax()) + m_traits->get_interval_x(); + Interval_nt y2 = Interval_nt(bbox.ymin(), + bbox.ymax()) + m_traits->get_interval_y(); + + return x1.do_overlap(x2) && y1.do_overlap(y2); } - Interval_nt get_interval_y() const { - return m_interval_y; + bool operator()(const Primitive &q, const Bounding_box &bbox) const + { + Interval_nt x1 = Interval_nt(q.datum().bbox().xmin(), + q.datum().bbox().xmax()); + Interval_nt y1 = Interval_nt(q.datum().bbox().ymin(), + q.datum().bbox().ymax()); + Interval_nt x2 = Interval_nt(bbox.xmin(), + bbox.xmax()) + m_traits->get_interval_x(); + Interval_nt y2 = Interval_nt(bbox.ymin(), + bbox.ymax()) + m_traits->get_interval_y(); + + return x1.do_overlap(x2) && y1.do_overlap(y2); } - Point get_translation_point() const { - return m_translation_point; + bool operator()(const Bounding_box &q, const Primitive &pr) const + { + Datum tr_pr = pr.datum().transform(typename GeomTraits::Aff_transformation_2( + Translation(), + Vector_2(ORIGIN, m_traits->get_translation_point()))); + + return do_overlap(q, tr_pr.bbox()); } - // Put the n/2 smallest primitives in the front, the n/2 largest primitives - // in the back. They are compared along the bbox' longest axis. - class Sort_primitives { - public: - template - void operator()(PrimitiveIterator first, - PrimitiveIterator beyond, - const Bounding_box &bbox) const { - PrimitiveIterator middle = first + (beyond - first) / 2; + bool operator()(const Primitive &q, const Primitive &pr) const + { + Datum tr_pr = pr.datum().transform(typename GeomTraits::Aff_transformation_2( + Translation(), Vector_2(ORIGIN, m_traits->get_translation_point()))); - if (bbox.xmax()-bbox.xmin() >= bbox.ymax()-bbox.ymin()) { - std::nth_element(first, middle, beyond, less_x); // sort along x - } else { - std::nth_element(first, middle, beyond, less_y); // sort along y - } - } - }; + if (!do_overlap(q.datum().bbox(), tr_pr.bbox())) + { + return false; + } - Sort_primitives sort_primitives_object() const { - return Sort_primitives(); + return do_intersect(q.datum(), tr_pr); } + }; - // Computes the bounding box of a set of primitives - class Compute_bbox { - public: - template - Bounding_box operator()(ConstPrimitiveIterator first, - ConstPrimitiveIterator beyond) const { - Bounding_box bbox = first->datum().bbox(); - - for (++first; first != beyond; ++first) { - bbox = bbox + first->datum().bbox(); - } - - return bbox; - } - }; - - Compute_bbox compute_bbox_object() const { - return Compute_bbox(); - } - - class Do_intersect { - - private: - - AABB_traits_2 *m_traits; - - public: - - Do_intersect(AABB_traits_2 *_traits): m_traits(_traits) {} - - bool operator()(const Bounding_box &q, const Bounding_box &bbox) const { - Interval_nt x1 = Interval_nt(q.xmin(), q.xmax()); - Interval_nt y1 = Interval_nt(q.ymin(), q.ymax()); - Interval_nt x2 = Interval_nt(bbox.xmin(), bbox.xmax()) + m_traits->get_interval_x(); - Interval_nt y2 = Interval_nt(bbox.ymin(), bbox.ymax()) + m_traits->get_interval_y(); - - return x1.do_overlap(x2) && y1.do_overlap(y2); - } - - bool operator()(const Primitive &q, const Bounding_box &bbox) const { - Interval_nt x1 = Interval_nt(q.datum().bbox().xmin(), q.datum().bbox().xmax()); - Interval_nt y1 = Interval_nt(q.datum().bbox().ymin(), q.datum().bbox().ymax()); - Interval_nt x2 = Interval_nt(bbox.xmin(), bbox.xmax()) + m_traits->get_interval_x(); - Interval_nt y2 = Interval_nt(bbox.ymin(), bbox.ymax()) + m_traits->get_interval_y(); - - return x1.do_overlap(x2) && y1.do_overlap(y2); - } - - bool operator()(const Bounding_box &q, const Primitive &pr) const { - Datum tr_pr = pr.datum().transform(typename GeomTraits::Aff_transformation_2(Translation(), - Vector_2(ORIGIN, m_traits->get_translation_point()))); - - return do_overlap(q, tr_pr.bbox()); - } - - bool operator()(const Primitive &q, const Primitive &pr) const { - Datum tr_pr = pr.datum().transform(typename GeomTraits::Aff_transformation_2(Translation(), Vector_2(ORIGIN, m_traits->get_translation_point()))); - - if (!do_overlap(q.datum().bbox(), tr_pr.bbox())) { - return false; - } - - return do_intersect(q.datum(), tr_pr); - } - }; - - Do_intersect do_intersect_object() { - return Do_intersect(this); - } + Do_intersect do_intersect_object() + { + return Do_intersect(this); + } private: - Point m_translation_point; - Interval_nt m_interval_x; - Interval_nt m_interval_y; + Point m_translation_point; + Interval_nt m_interval_x; + Interval_nt m_interval_y; - // Comparison functions - static bool less_x(const Primitive &pr1, const Primitive &pr2) { - return pr1.reference_point().x() < pr2.reference_point().x(); - } + // Comparison functions + static bool less_x(const Primitive &pr1, const Primitive &pr2) + { + return pr1.reference_point().x() < pr2.reference_point().x(); + } - static bool less_y(const Primitive &pr1, const Primitive &pr2) { - return pr1.reference_point().y() < pr2.reference_point().y(); - } + static bool less_y(const Primitive &pr1, const Primitive &pr2) + { + return pr1.reference_point().y() < pr2.reference_point().y(); + } }; } // namespace CGAL diff --git a/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/Minkowski_sum_by_reduced_convolution_2.h b/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/Minkowski_sum_by_reduced_convolution_2.h index 81c6da4738f..a91df2f3d32 100644 --- a/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/Minkowski_sum_by_reduced_convolution_2.h +++ b/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/Minkowski_sum_by_reduced_convolution_2.h @@ -16,331 +16,393 @@ namespace CGAL { template -class Minkowski_sum_by_reduced_convolution_2 { +class Minkowski_sum_by_reduced_convolution_2 +{ private: - typedef Kernel_ Kernel; - typedef Container_ Container; + typedef Kernel_ Kernel; + typedef Container_ Container; - // Basic types: - typedef CGAL::Polygon_2 Polygon_2; - typedef typename Kernel::Point_2 Point_2; - typedef typename Kernel::Vector_2 Vector_2; - typedef typename Kernel::Direction_2 Direction_2; - typedef typename Kernel::Triangle_2 Triangle_2; - typedef typename Kernel::FT FT; + // Basic types: + typedef CGAL::Polygon_2 Polygon_2; + typedef typename Kernel::Point_2 Point_2; + typedef typename Kernel::Vector_2 Vector_2; + typedef typename Kernel::Direction_2 Direction_2; + typedef typename Kernel::Triangle_2 Triangle_2; + typedef typename Kernel::FT FT; - // Segment-related types: - typedef Arr_segment_traits_2 Traits_2; - typedef typename Traits_2::X_monotone_curve_2 Segment_2; - typedef std::list Segment_list; - typedef Arr_default_dcel Dcel; - typedef std::pair State; + // Segment-related types: + typedef Arr_segment_traits_2 Traits_2; + typedef typename Traits_2::X_monotone_curve_2 Segment_2; + typedef std::list Segment_list; + typedef Arr_default_dcel Dcel; + typedef std::pair State; - // Arrangement-related types: - typedef Arrangement_with_history_2 Arrangement_history_2; - typedef typename Arrangement_history_2::Halfedge_handle Halfedge_handle; - typedef typename Arrangement_history_2::Face_iterator Face_iterator; - typedef typename Arrangement_history_2::Face_handle Face_handle; - typedef typename Arrangement_history_2::Ccb_halfedge_circulator Ccb_halfedge_circulator; - typedef typename Arrangement_history_2::Originating_curve_iterator Originating_curve_iterator; + // Arrangement-related types: + typedef Arrangement_with_history_2 Arrangement_history_2; + typedef typename Arrangement_history_2::Halfedge_handle Halfedge_handle; + typedef typename Arrangement_history_2::Face_iterator Face_iterator; + typedef typename Arrangement_history_2::Face_handle Face_handle; + typedef typename Arrangement_history_2::Ccb_halfedge_circulator + Ccb_halfedge_circulator; + typedef typename Arrangement_history_2::Originating_curve_iterator + Originating_curve_iterator; - // Function object types: - typename Kernel::Construct_translated_point_2 f_add; - typename Kernel::Construct_vector_2 f_vector; - typename Kernel::Construct_direction_2 f_direction; - typename Kernel::Orientation_2 f_orientation; - typename Kernel::Compare_xy_2 f_compare_xy; - typename Kernel::Counterclockwise_in_between_2 f_ccw_in_between; + // Function object types: + typename Kernel::Construct_translated_point_2 f_add; + typename Kernel::Construct_vector_2 f_vector; + typename Kernel::Construct_direction_2 f_direction; + typename Kernel::Orientation_2 f_orientation; + typename Kernel::Compare_xy_2 f_compare_xy; + typename Kernel::Counterclockwise_in_between_2 f_ccw_in_between; public: - Minkowski_sum_by_reduced_convolution_2() { - // Obtain kernel functors - Kernel ker; - f_add = ker.construct_translated_point_2_object(); - f_vector = ker.construct_vector_2_object(); - f_direction = ker.construct_direction_2_object(); - f_orientation = ker.orientation_2_object(); - f_compare_xy = ker.compare_xy_2_object(); - f_ccw_in_between = ker.counterclockwise_in_between_2_object(); + Minkowski_sum_by_reduced_convolution_2() + { + // Obtain kernel functors + Kernel ker; + f_add = ker.construct_translated_point_2_object(); + f_vector = ker.construct_vector_2_object(); + f_direction = ker.construct_direction_2_object(); + f_orientation = ker.orientation_2_object(); + f_compare_xy = ker.compare_xy_2_object(); + f_ccw_in_between = ker.counterclockwise_in_between_2_object(); + } + + template + void operator()(const Polygon_2 &pgn1, const Polygon_2 &pgn2, + Polygon_2 &outer_boundary, OutputIterator holes) const + { + + Timer timer; // TODO: remove when optimization is done + timer.start(); + + CGAL_precondition(pgn1.is_simple()); + CGAL_precondition(pgn2.is_simple()); + CGAL_precondition(pgn1.orientation() == COUNTERCLOCKWISE); + CGAL_precondition(pgn2.orientation() == COUNTERCLOCKWISE); + + timer.stop(); + std::cout << timer.time() << " s: Preconditions" << std::endl; + timer.reset(); + timer.start(); + + // Initialize collision detector. It operates on pgn2 and on the inversed pgn1: + const Polygon_2 inversed_pgn1 = transform(Aff_transformation_2(SCALING, + -1), pgn1); + AABB_collision_detector_2 collision_detector(pgn2, + inversed_pgn1); + + timer.stop(); + std::cout << timer.time() << " s: AABB init" << std::endl; + timer.reset(); + timer.start(); + + // Compute the reduced convolution + Segment_list reduced_convolution; + build_reduced_convolution(pgn1, pgn2, reduced_convolution); + + timer.stop(); + std::cout << timer.time() << " s: Convolution" << std::endl; + timer.reset(); + timer.start(); + + // Insert the segments into an arrangement + Arrangement_history_2 arr; + insert(arr, reduced_convolution.begin(), reduced_convolution.end()); + + timer.stop(); + std::cout << timer.time() << " s: Arrangement" << std::endl; + timer.reset(); + timer.start(); + + // Trace the outer loop and put it in 'outer_boundary' + get_outer_loop(arr, outer_boundary); + + timer.stop(); + std::cout << timer.time() << " s: Outer Loop" << std::endl; + timer.reset(); + timer.start(); + + // Check for each face whether it is a hole in the M-sum. If it is, add it to 'holes'. + for (Face_iterator face = arr.faces_begin(); face != arr.faces_end(); ++face) + { + handle_face(arr, face, holes, collision_detector); } - template - void operator()(const Polygon_2 &pgn1, const Polygon_2 &pgn2, - Polygon_2 &outer_boundary, OutputIterator holes) const { - - Timer timer; // TODO: remove when optimization is done - timer.start(); - - CGAL_precondition(pgn1.is_simple()); - CGAL_precondition(pgn2.is_simple()); - CGAL_precondition(pgn1.orientation() == COUNTERCLOCKWISE); - CGAL_precondition(pgn2.orientation() == COUNTERCLOCKWISE); - - timer.stop(); - std::cout << timer.time() << " s: Preconditions" << std::endl; - timer.reset(); - timer.start(); - - // Initialize collision detector. It operates on pgn2 and on the inversed pgn1: - const Polygon_2 inversed_pgn1 = transform(Aff_transformation_2(SCALING, -1), pgn1); - AABB_collision_detector_2 collision_detector(pgn2, inversed_pgn1); - - timer.stop(); - std::cout << timer.time() << " s: AABB init" << std::endl; - timer.reset(); - timer.start(); - - // Compute the reduced convolution - Segment_list reduced_convolution; - build_reduced_convolution(pgn1, pgn2, reduced_convolution); - - timer.stop(); - std::cout << timer.time() << " s: Convolution" << std::endl; - timer.reset(); - timer.start(); - - // Insert the segments into an arrangement - Arrangement_history_2 arr; - insert(arr, reduced_convolution.begin(), reduced_convolution.end()); - - timer.stop(); - std::cout << timer.time() << " s: Arrangement" << std::endl; - timer.reset(); - timer.start(); - - // Trace the outer loop and put it in 'outer_boundary' - get_outer_loop(arr, outer_boundary); - - timer.stop(); - std::cout << timer.time() << " s: Outer Loop" << std::endl; - timer.reset(); - timer.start(); - - // Check for each face whether it is a hole in the M-sum. If it is, add it to 'holes'. - for (Face_iterator face = arr.faces_begin(); face != arr.faces_end(); ++face) { - handle_face(arr, face, holes, collision_detector); - } - - timer.stop(); - std::cout << timer.time() << " s: Holes" << std::endl; - } + timer.stop(); + std::cout << timer.time() << " s: Holes" << std::endl; + } private: - // Builds the reduced convolution using a fiber grid approach. For each - // starting vertex, try to add two outgoing next states. If a visited - // vertex is reached, then do not explore further. This is a BFS-like - // iteration beginning from each vertex in the first column of the fiber - // grid. - void build_reduced_convolution(const Polygon_2 &pgn1, const Polygon_2 &pgn2, Segment_list &reduced_convolution) const { - unsigned int n1 = pgn1.size(); - unsigned int n2 = pgn2.size(); + // Builds the reduced convolution using a fiber grid approach. For each + // starting vertex, try to add two outgoing next states. If a visited + // vertex is reached, then do not explore further. This is a BFS-like + // iteration beginning from each vertex in the first column of the fiber + // grid. + void build_reduced_convolution(const Polygon_2 &pgn1, const Polygon_2 &pgn2, + Segment_list &reduced_convolution) const + { + unsigned int n1 = pgn1.size(); + unsigned int n2 = pgn2.size(); - // Init the direcions of both polygons - std::vector p1_dirs = directions_of_polygon(pgn1); - std::vector p2_dirs = directions_of_polygon(pgn2); + // Init the direcions of both polygons + std::vector p1_dirs = directions_of_polygon(pgn1); + std::vector p2_dirs = directions_of_polygon(pgn2); - // Contains states that were already visited - boost::unordered_set visited_states; + // Contains states that were already visited + boost::unordered_set visited_states; - // Init the queue with vertices from the first column - std::queue state_queue; - for (int i = n1-1; i >= 0; --i) { - state_queue.push(State(i, 0)); - } - - while (state_queue.size() > 0) { - State curr_state = state_queue.front(); - state_queue.pop(); - - int i1 = curr_state.first; - int i2 = curr_state.second; - - // If this state was already visited, skip it - if (visited_states.count(curr_state) > 0) { - continue; - } - visited_states.insert(curr_state); - - int next_i1 = (i1+1) % n1; - int next_i2 = (i2+1) % n2; - int prev_i1 = (n1+i1-1) % n1; - int prev_i2 = (n2+i2-1) % n2; - - // Try two transitions: From (i,j) to (i+1,j) and to (i,j+1). Add - // the respective segments, if they are in the reduced convolution. - for(int step_in_pgn1 = 0; step_in_pgn1 <= 1; step_in_pgn1++) { - int new_i1, new_i2; - if (step_in_pgn1) { - new_i1 = next_i1; - new_i2 = i2; - } else { - new_i1 = i1; - new_i2 = next_i2; - } - - // If the segment's direction lies counterclockwise in between - // the other polygon's vertex' ingoing and outgoing directions, - // the segment belongs to the full convolution. - bool belongs_to_convolution; - if (step_in_pgn1) { - belongs_to_convolution = f_ccw_in_between(p1_dirs[i1], p2_dirs[prev_i2], p2_dirs[i2]) || - p1_dirs[i1] == p2_dirs[i2]; - } else { - belongs_to_convolution = f_ccw_in_between(p2_dirs[i2], p1_dirs[prev_i1], p1_dirs[i1]) || - p2_dirs[i2] == p1_dirs[prev_i1]; - } - - if (belongs_to_convolution) { - state_queue.push(State(new_i1, new_i2)); - - // Only edges added to convex vertices can be on the M-sum's boundary. - // This filter only leaves the *reduced* convolution. - bool convex; - if (step_in_pgn1) { - convex = is_convex(pgn2[prev_i2], pgn2[i2], pgn2[next_i2]); - } else { - convex = is_convex(pgn1[prev_i1], pgn1[i1], pgn1[next_i1]); - } - - if (convex) { - Point_2 start_point = get_point(i1, i2, pgn1, pgn2); - Point_2 end_point = get_point(new_i1, new_i2, pgn1, pgn2); - reduced_convolution.push_back(Segment_2(start_point, end_point)); - } - } - } - } + // Init the queue with vertices from the first column + std::queue state_queue; + for (int i = n1-1; i >= 0; --i) + { + state_queue.push(State(i, 0)); } - // Returns a sorted list of the polygon's edges - std::vector directions_of_polygon(const Polygon_2 &p) const { - std::vector directions; - unsigned int n = p.size(); + while (state_queue.size() > 0) + { + State curr_state = state_queue.front(); + state_queue.pop(); - for (int i = 0; i < n-1; ++i) { - directions.push_back(f_direction(f_vector(p[i], p[i+1]))); + int i1 = curr_state.first; + int i2 = curr_state.second; + + // If this state was already visited, skip it + if (visited_states.count(curr_state) > 0) + { + continue; + } + visited_states.insert(curr_state); + + int next_i1 = (i1+1) % n1; + int next_i2 = (i2+1) % n2; + int prev_i1 = (n1+i1-1) % n1; + int prev_i2 = (n2+i2-1) % n2; + + // Try two transitions: From (i,j) to (i+1,j) and to (i,j+1). Add + // the respective segments, if they are in the reduced convolution. + for(int step_in_pgn1 = 0; step_in_pgn1 <= 1; step_in_pgn1++) + { + int new_i1, new_i2; + if (step_in_pgn1) + { + new_i1 = next_i1; + new_i2 = i2; } - directions.push_back(f_direction(f_vector(p[n-1], p[0]))); - - return directions; - } - - bool is_convex(const Point_2 &prev, const Point_2 &curr, const Point_2 &next) const { - return f_orientation(prev, curr, next) == LEFT_TURN; - } - - // Returns the point corresponding to a state (i,j). - Point_2 get_point(int i1, int i2, const Polygon_2 &pgn1, const Polygon_2 &pgn2) const { - - return f_add(pgn1[i1], Vector_2(Point_2(ORIGIN), pgn2[i2])); - } - - // Put the outer loop of the arrangement in 'outer_boundary' - void get_outer_loop(Arrangement_history_2 &arr, Polygon_2 &outer_boundary) const { - Ccb_halfedge_circulator circ_start = *(arr.unbounded_face()->holes_begin()); - Ccb_halfedge_circulator circ = circ_start; - - do { - outer_boundary.push_back(circ->source()->point()); - } while (--circ != circ_start); - } - - // Check whether the face is on the M-sum's border. Add it to 'holes' if it is. - template - void handle_face(const Arrangement_history_2 &arr, const Face_handle face, OutputIterator holes, AABB_collision_detector_2 &collision_detector) const { - - // If the face contains holes, it can't be on the Minkowski sum's border - if (face->holes_begin() != face->holes_end()) { - return; + else + { + new_i1 = i1; + new_i2 = next_i2; } - Ccb_halfedge_circulator start = face->outer_ccb(); - Ccb_halfedge_circulator circ = start; - - // The face needs to be orientable - do { - if (!do_original_edges_have_same_direction(arr, circ)) { - return; - } - } while (++circ != start); - - // When the reversed polygon 1, translated by a point inside of this face, collides with polygon 2, this cannot be a hole - Point_2 inner_point = get_point_in_face(face); - if (collision_detector.check_collision(inner_point)) { - return; + // If the segment's direction lies counterclockwise in between + // the other polygon's vertex' ingoing and outgoing directions, + // the segment belongs to the full convolution. + bool belongs_to_convolution; + if (step_in_pgn1) + { + belongs_to_convolution = f_ccw_in_between(p1_dirs[i1], p2_dirs[prev_i2], + p2_dirs[i2]) || + p1_dirs[i1] == p2_dirs[i2]; + } + else + { + belongs_to_convolution = f_ccw_in_between(p2_dirs[i2], p1_dirs[prev_i1], + p1_dirs[i1]) || + p2_dirs[i2] == p1_dirs[prev_i1]; } - // At this point, the face is a real hole, add it to 'holes' - Polygon_2 pgn_hole; - circ = start; + if (belongs_to_convolution) + { + state_queue.push(State(new_i1, new_i2)); - do { - pgn_hole.push_back(circ->source()->point()); - } while (--circ != start); + // Only edges added to convex vertices can be on the M-sum's boundary. + // This filter only leaves the *reduced* convolution. + bool convex; + if (step_in_pgn1) + { + convex = is_convex(pgn2[prev_i2], pgn2[i2], pgn2[next_i2]); + } + else + { + convex = is_convex(pgn1[prev_i1], pgn1[i1], pgn1[next_i1]); + } - *holes = pgn_hole; - ++holes; + if (convex) + { + Point_2 start_point = get_point(i1, i2, pgn1, pgn2); + Point_2 end_point = get_point(new_i1, new_i2, pgn1, pgn2); + reduced_convolution.push_back(Segment_2(start_point, end_point)); + } + } + } + } + } + + // Returns a sorted list of the polygon's edges + std::vector directions_of_polygon(const Polygon_2 &p) const + { + std::vector directions; + unsigned int n = p.size(); + + for (int i = 0; i < n-1; ++i) + { + directions.push_back(f_direction(f_vector(p[i], p[i+1]))); + } + directions.push_back(f_direction(f_vector(p[n-1], p[0]))); + + return directions; + } + + bool is_convex(const Point_2 &prev, const Point_2 &curr, + const Point_2 &next) const + { + return f_orientation(prev, curr, next) == LEFT_TURN; + } + + // Returns the point corresponding to a state (i,j). + Point_2 get_point(int i1, int i2, const Polygon_2 &pgn1, + const Polygon_2 &pgn2) const + { + + return f_add(pgn1[i1], Vector_2(Point_2(ORIGIN), pgn2[i2])); + } + + // Put the outer loop of the arrangement in 'outer_boundary' + void get_outer_loop(Arrangement_history_2 &arr, + Polygon_2 &outer_boundary) const + { + Ccb_halfedge_circulator circ_start = *(arr.unbounded_face()->holes_begin()); + Ccb_halfedge_circulator circ = circ_start; + + do + { + outer_boundary.push_back(circ->source()->point()); + } + while (--circ != circ_start); + } + + // Check whether the face is on the M-sum's border. Add it to 'holes' if it is. + template + void handle_face(const Arrangement_history_2 &arr, const Face_handle face, + OutputIterator holes, AABB_collision_detector_2 + &collision_detector) const + { + + // If the face contains holes, it can't be on the Minkowski sum's border + if (face->holes_begin() != face->holes_end()) + { + return; } - // Check whether the convolution's original edge(s) had the same direction as the arrangement's half edge - bool do_original_edges_have_same_direction(const Arrangement_history_2 &arr, const Halfedge_handle he) const { - Originating_curve_iterator segment_itr; + Ccb_halfedge_circulator start = face->outer_ccb(); + Ccb_halfedge_circulator circ = start; - for (segment_itr = arr.originating_curves_begin(he); segment_itr != arr.originating_curves_end(he); ++segment_itr) { - if (f_compare_xy(segment_itr->source(), segment_itr->target()) == (Comparison_result)he->direction()) { - return false; - } - } + // The face needs to be orientable + do + { + if (!do_original_edges_have_same_direction(arr, circ)) + { + return; + } + } + while (++circ != start); - return true; + // When the reversed polygon 1, translated by a point inside of this face, collides with polygon 2, this cannot be a hole + Point_2 inner_point = get_point_in_face(face); + if (collision_detector.check_collision(inner_point)) + { + return; } - // Return a point in the face's interior by finding a diagonal - Point_2 get_point_in_face(const Face_handle face) const { - Ccb_halfedge_circulator current_edge = face->outer_ccb(); - Ccb_halfedge_circulator next_edge = current_edge; - next_edge++; + // At this point, the face is a real hole, add it to 'holes' + Polygon_2 pgn_hole; + circ = start; - Point_2 a, v, b; - - // Move over the face's vertices until a convex corner is encountered: - do { - a = current_edge->source()->point(); - v = current_edge->target()->point(); - b = next_edge->target()->point(); - - current_edge++; - next_edge++; - } while (!is_convex(a, v, b)); - - Triangle_2 ear(a, v, b); - FT min_distance = -1; - Point_2 min_q; - - // Of the remaining vertices, find the one inside of the "ear" with minimal distance to v: - while (++next_edge != current_edge) { - Point_2 q = next_edge->target()->point(); - if (ear.has_on_bounded_side(q)) { - FT distance = squared_distance(q, v); - if (min_distance == -1 || distance < min_distance) { - min_distance = distance; - min_q = q; - } - } - } - - // If there was no vertex inside of the ear, return it's centroid. - // Otherwise, return a point between v and min_q. - if (min_distance == -1) { - return centroid(ear); - } else { - return midpoint(v, min_q); - } + do + { + pgn_hole.push_back(circ->source()->point()); } + while (--circ != start); + + *holes = pgn_hole; + ++holes; + } + + // Check whether the convolution's original edge(s) had the same direction as the arrangement's half edge + bool do_original_edges_have_same_direction(const Arrangement_history_2 &arr, + const Halfedge_handle he) const + { + Originating_curve_iterator segment_itr; + + for (segment_itr = arr.originating_curves_begin(he); + segment_itr != arr.originating_curves_end(he); ++segment_itr) + { + if (f_compare_xy(segment_itr->source(), + segment_itr->target()) == (Comparison_result)he->direction()) + { + return false; + } + } + + return true; + } + + // Return a point in the face's interior by finding a diagonal + Point_2 get_point_in_face(const Face_handle face) const + { + Ccb_halfedge_circulator current_edge = face->outer_ccb(); + Ccb_halfedge_circulator next_edge = current_edge; + next_edge++; + + Point_2 a, v, b; + + // Move over the face's vertices until a convex corner is encountered: + do + { + a = current_edge->source()->point(); + v = current_edge->target()->point(); + b = next_edge->target()->point(); + + current_edge++; + next_edge++; + } + while (!is_convex(a, v, b)); + + Triangle_2 ear(a, v, b); + FT min_distance = -1; + Point_2 min_q; + + // Of the remaining vertices, find the one inside of the "ear" with minimal distance to v: + while (++next_edge != current_edge) + { + Point_2 q = next_edge->target()->point(); + if (ear.has_on_bounded_side(q)) + { + FT distance = squared_distance(q, v); + if (min_distance == -1 || distance < min_distance) + { + min_distance = distance; + min_q = q; + } + } + } + + // If there was no vertex inside of the ear, return it's centroid. + // Otherwise, return a point between v and min_q. + if (min_distance == -1) + { + return centroid(ear); + } + else + { + return midpoint(v, min_q); + } + } }; } // namespace CGAL