diff --git a/GraphicsView/include/CGAL/Qt/Basic_viewer_qt.h b/GraphicsView/include/CGAL/Qt/Basic_viewer_qt.h index 37873145733b..59736fc103d1 100644 --- a/GraphicsView/include/CGAL/Qt/Basic_viewer_qt.h +++ b/GraphicsView/include/CGAL/Qt/Basic_viewer_qt.h @@ -54,8 +54,8 @@ #include #include -#define CGAL_BASIC_VIEWER_INIT_SIZE_X 500 -#define CGAL_BASIC_VIEWER_INIT_SIZE_Y 450 +#define CGAL_BASIC_VIEWER_INIT_SIZE_X 1500 +#define CGAL_BASIC_VIEWER_INIT_SIZE_Y 1500 namespace CGAL { diff --git a/Straight_skeleton_2/doc/Straight_skeleton_2/CGAL/extrude_skeleton.h b/Straight_skeleton_2/doc/Straight_skeleton_2/CGAL/extrude_skeleton.h index 6fcb059ae2d5..7f91ca805d5d 100644 --- a/Straight_skeleton_2/doc/Straight_skeleton_2/CGAL/extrude_skeleton.h +++ b/Straight_skeleton_2/doc/Straight_skeleton_2/CGAL/extrude_skeleton.h @@ -5,18 +5,21 @@ namespace CGAL { * * \brief constructs the straight skeleton-based extrusion of a polygon with holes. * -* Given a polygon with holes and a set of weights, the skeleton extrusion is a volume constructed -* from the weighted straight skeleton by associating a height to the vertices of the skeleton, -* which corresponds to the time at the vertex. The input polygon is placed at `z = 0`. -* -* This function allows cropping the extruded skeleton at a maximum height, using the optional -* `maximum_height()` named parameter. +* Given a polygon with holes and a set of weights (or angles) associated to its edges, +* the skeleton extrusion is a volume constructed from the weighted straight skeleton +* by associating a height to the vertices of the skeleton, which corresponds to the time +* at the vertex. The input polygon is placed at `z = 0`. * * The result is a closed, 2-manifold surface triangle mesh. Note that this mesh can have non-local * self-intersections if a maximal height is provided due to possible (geometric) non-manifold occurences. * -* @tparam PolygonWithHoles must be a model of `SequenceContainer` with value type `InK::Point_2` (e.g. `Polygon_2`) - or a model of `GeneralPolygonWithHoles_2` (e.g. `Polygon_with_holes_2`). +* It is possible to crop the extruded skeleton at a maximum height using the optional +* `maximum_height()` named parameter. A maximum height must be specified if the weights (or angles) +* associated to the edges of the input polygon correspond an outward extrusion, i.e. if no weight +* is greater than zero (or no angle is smaller than `90` degrees). +* +* @tparam Polygon must be a model of `SequenceContainer` with value type `InK::Point_2` (e.g. `Polygon_2`) + or a model of `GeneralPolygonWithHoles_2` (e.g. `Polygon_with_holes_2`). * @tparam PolygonMesh a model of `MutableFaceGraph` * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" * @@ -42,7 +45,7 @@ namespace CGAL { * \cgalParamType{a model of `Range` whose value type is a model of `Range` whose value type is `FT`} * \cgalParamDefault{an empty range (uniform weights are used)} * \cgalParamExtra{Angles are measured in degrees and should be strictly within `0` and `180` degrees -* and should be eitger all acute (inward extrusion) or all obtuse (outward extrusion).} +* and should be either all acute (inward extrusion) or all obtuse (outward extrusion).} * \cgalParamExtra{This parameter is ignored if the `weights` parameter is provided.} * \cgalParamExtra{The conversion to weights involves trigonometry and will be inexact, * even when using a number type with exact square roots.} diff --git a/Straight_skeleton_2/doc/Straight_skeleton_2/Straight_skeleton_2.txt b/Straight_skeleton_2/doc/Straight_skeleton_2/Straight_skeleton_2.txt index b9a2ddb1de77..3a4da692f799 100644 --- a/Straight_skeleton_2/doc/Straight_skeleton_2/Straight_skeleton_2.txt +++ b/Straight_skeleton_2/doc/Straight_skeleton_2/Straight_skeleton_2.txt @@ -142,14 +142,14 @@ The main entry points for straight skeletons are the following functions: \subsection Straight_skeleton_2Weights Weighted Straight Skeletons Weighted straight skeletons are a generalization of straight skeletons: contour edges are assigned -a positive weight, which can be understood as assigning a speed to the wavefront spawned from +a weight, which can be understood as assigning a speed to the wavefront spawned from the contour edge. Vertices no longer move along the angular bisector between two contour edges, but along a weighted bisector. -This \cgal package supports positive multiplicative weights: if the supporting line of a contour edge +This \cgal package supports multiplicative weights: if the supporting line of a contour edge is described through the equation `ax+by+c=0` then the supporting line of the offset edge -at distance `t` is `ax+by+c-t=0`. With a multiplicative weight `w`, the equation becomes `w(ax+by+c)-t=0`. -Therefore, a larger weight implies a faster moving front. +at distance `t` is `ax+by+c-t=0`. With a multiplicative weight `w`, the equation becomes `ax+by+c-w*t=0`. +Therefore, a larger weight implies a slower moving front. \cgalFigureAnchor{SLSWeight}
@@ -282,7 +282,7 @@ to its offset distance, the resulting roof is "correct" in that water will alway to the contour edges (the roof's border), regardless of where it falls on the roof. Laycock and Day \cgalCite{cgal:ld-agrm-03} give an algorithm for roof design based on the straight skeleton. -This \cgal package implements skeleton extrusion for polygons with holes, with support for positive +This \cgal package implements skeleton extrusion for polygons with holes, with support for multiplicative weights. The output is a closed, combinatorially 2-manifold surface triangle mesh. \cgalFigureAnchor{SLSExtrusion} diff --git a/Straight_skeleton_2/examples/Straight_skeleton_2/Create_straight_skeleton_2.cpp b/Straight_skeleton_2/examples/Straight_skeleton_2/Create_straight_skeleton_2.cpp index 5aae2d904f1e..fba376d8c7f1 100644 --- a/Straight_skeleton_2/examples/Straight_skeleton_2/Create_straight_skeleton_2.cpp +++ b/Straight_skeleton_2/examples/Straight_skeleton_2/Create_straight_skeleton_2.cpp @@ -1,9 +1,34 @@ +#define CGAL_SLS_DEBUG_DRAW + +#include +#include +#include + +#define CGAL_SLS_PRINT_QUEUE_BEFORE_EACH_POP +#define CGAL_STRAIGHT_SKELETON_ENABLE_TRACE 10000000 +#define CGAL_STRAIGHT_SKELETON_TRAITS_ENABLE_TRACE 10000000 +#define CGAL_STRAIGHT_SKELETON_VALIDITY_ENABLE_TRACE +#define CGAL_POLYGON_OFFSET_ENABLE_TRACE 10000000 + +void Straight_skeleton_external_trace(std::string m) +{ + std::cout << std::setprecision(17) << m << std::endl << std::endl ; +} + +void Straight_skeleton_traits_external_trace(std::string m) +{ + std::cout << std::setprecision(17) << m << std::endl << std::endl ; +} + + #include +#include +#include #include #include +#include #include -#include #include @@ -11,14 +36,16 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel K ; +typedef K::FT FT ; typedef K::Point_2 Point ; typedef CGAL::Polygon_2 Polygon_2 ; -typedef CGAL::Straight_skeleton_2 Ss ; -typedef boost::shared_ptr SsPtr ; +typedef CGAL::Straight_skeleton_2 Ss ; +typedef boost::shared_ptr SsPtr ; int main() { +#if 0 Polygon_2 poly ; poly.push_back( Point(-1,-1) ) ; poly.push_back( Point(0,-12) ) ; @@ -30,6 +57,7 @@ int main() poly.push_back( Point(-12,0) ) ; assert(poly.is_counterclockwise_oriented()); + CGAL::draw(poly); // You can pass the polygon via an iterator pair SsPtr iss = CGAL::create_interior_straight_skeleton_2(poly.vertices_begin(), poly.vertices_end()); @@ -45,6 +73,21 @@ int main() CGAL::Straight_skeletons_2::IO::print_straight_skeleton(*oss); CGAL::draw(*oss); +#else + Polygon_2 poly ; + poly.push_back( Point(0, 0) ); + // poly.push_back( Point(5, 0) ); + poly.push_back( Point(10, 0) ); + poly.push_back( Point(10, 10) ); + poly.push_back( Point(0, 10) ); + + // You can also use weights to modify the speed of some of the fronts + std::vector weights = {{ 3, 5, 6, 7 }}; + SsPtr iss = CGAL::create_interior_weighted_straight_skeleton_2(poly.vertices_begin(), poly.vertices_end(), + weights.begin(), weights.end()); + CGAL::Straight_skeletons_2::IO::print_straight_skeleton(*iss); + CGAL::draw(*iss); +#endif return EXIT_SUCCESS; } diff --git a/Straight_skeleton_2/examples/Straight_skeleton_2/extrude_skeleton.cpp b/Straight_skeleton_2/examples/Straight_skeleton_2/extrude_skeleton.cpp index 073b50964c09..8555802410e2 100644 --- a/Straight_skeleton_2/examples/Straight_skeleton_2/extrude_skeleton.cpp +++ b/Straight_skeleton_2/examples/Straight_skeleton_2/extrude_skeleton.cpp @@ -1,3 +1,25 @@ +#define CGAL_SLS_DEBUG_DRAW + +#include +#include +#include + +#define CGAL_SLS_PRINT_QUEUE_BEFORE_EACH_POP +#define CGAL_STRAIGHT_SKELETON_ENABLE_TRACE 10000000 +#define CGAL_STRAIGHT_SKELETON_TRAITS_ENABLE_TRACE 10000000 +#define CGAL_STRAIGHT_SKELETON_VALIDITY_ENABLE_TRACE +#define CGAL_POLYGON_OFFSET_ENABLE_TRACE 10000000 + +void Straight_skeleton_external_trace(std::string m) +{ + std::cout << std::setprecision(17) << m << std::endl << std::endl ; +} + +void Straight_skeleton_traits_external_trace(std::string m) +{ + std::cout << std::setprecision(17) << m << std::endl << std::endl ; +} + #include #include #include @@ -170,16 +192,17 @@ int main(int argc, char** argv) Mesh sm; if(use_angles) - CGAL::extrude_skeleton(pwh, sm, CGAL::parameters::angles(speeds).maximum_height(height)); + CGAL::extrude_skeleton(pwh, sm, CGAL::parameters::angles(speeds).maximum_height(height).verbose(true)); else - CGAL::extrude_skeleton(pwh, sm, CGAL::parameters::weights(speeds).maximum_height(height)); + CGAL::extrude_skeleton(pwh, sm, CGAL::parameters::weights(speeds).maximum_height(height).verbose(true)); timer.stop(); std::cout << "Extrusion computation took " << timer.time() << " s." << std::endl; - - CGAL::draw(sm); + std::cout << num_vertices(sm) << " vertices, " << num_faces(sm) << " faces" << std::endl; CGAL::IO::write_polygon_mesh("extruded_skeleton.off", sm, CGAL::parameters::stream_precision(17)); + CGAL::draw(sm); + return EXIT_SUCCESS; } diff --git a/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/IO/print.h b/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/IO/print.h index 33953f24ce39..b425a99fa068 100644 --- a/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/IO/print.h +++ b/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/IO/print.h @@ -31,7 +31,9 @@ template void print_vertex ( HDS_V const& v ) { if(v->has_infinite_time()) - std::cout << "F" << v->id() << " " ; + std::cout << "FN" << v->id() << " " ; + else if (v->is_contour()) + std::cout << "CN" << v->id() << " " ; else std::cout << "N" << v->id() << " " ; } @@ -92,6 +94,7 @@ void print_straight_skeleton( CGAL::Straight_skeleton_2 const& ss ) typedef CGAL::Straight_skeleton_2 Ss ; typedef typename Ss::Vertex_const_handle Vertex_const_handle ; + typedef typename Ss::Vertex_const_iterator Vertex_const_iterator ; typedef typename Ss::Halfedge_const_handle Halfedge_const_handle ; typedef typename Ss::Halfedge_const_iterator Halfedge_const_iterator ; @@ -103,6 +106,32 @@ void print_straight_skeleton( CGAL::Straight_skeleton_2 const& ss ) << " halfedges and " << ss.size_of_faces() << " faces" << std::endl ; + std::cout << "All vertices: " << std::endl; + for ( Vertex_const_iterator v = ss.vertices_begin(); v != ss.vertices_end(); ++v ) + { + print_vertex(v); + std::cout << "@ " << v->time() << std::endl; + } + + std::cout << "All halfedges: " << std::endl; + + for ( Halfedge_const_iterator h = ss.halfedges_begin(); h != ss.halfedges_end(); ++h ) + { + if(h->is_inner_bisector()) + std::cout << "IBH" << h->id() << " " << std::flush ; + else if(h->is_bisector()) + std::cout << "BH" << h->id() << " " << std::flush ; + else + std::cout << "CH" << h->id() << " " << std::flush ; + + print_vertex(h->prev()->vertex()) ; + print_point(h->prev()->vertex()->point()) ; + std::cout << " ==> " << std::flush ; + print_vertex(h->vertex()); + print_point(h->vertex()->point()) ; + std::cout << std::endl; + } + std::cout << "Faces " << std::endl; for ( Halfedge_const_iterator h = ss.halfedges_begin(); h != ss.halfedges_end(); ++h ) { @@ -131,24 +160,6 @@ void print_straight_skeleton( CGAL::Straight_skeleton_2 const& ss ) std::cout << std::endl; } - std::cout << "All halfedges: " << std::endl; - - for ( Halfedge_const_iterator h = ss.halfedges_begin(); h != ss.halfedges_end(); ++h ) - { - if(h->is_inner_bisector()) - std::cout << "IBH" << h->id() << " " << std::flush ; - else if(h->is_bisector()) - std::cout << "BH" << h->id() << " " << std::flush ; - else - std::cout << "CH" << h->id() << " " << std::flush ; - - print_vertex(h->prev()->vertex()) ; - print_point(h->prev()->vertex()->point()) ; - std::cout << " ==> " << std::flush ; - print_vertex(h->vertex()); - print_point(h->vertex()->point()) ; - std::cout << std::endl; - } } } // namespace IO diff --git a/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_aux.h b/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_aux.h index 606f8f27d0d2..0eb287c75eb4 100644 --- a/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_aux.h +++ b/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_aux.h @@ -128,9 +128,9 @@ class Triedge { ss << "{E" ; insert_handle_id(ss,t.e0()) ; - ss << ",E" ; + ss << " E" ; insert_handle_id(ss,t.e1()) ; - ss << ",E" ; + ss << " E" ; insert_handle_id(ss,t.e2()) ; ss << "}" ; return ss ; diff --git a/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_builder_2_impl.h b/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_builder_2_impl.h index f9fc7ca78c41..d92c4db12952 100644 --- a/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_builder_2_impl.h +++ b/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_builder_2_impl.h @@ -106,8 +106,7 @@ Straight_skeleton_builder_2::FindEdgeEvent( Vertex_handle aLNode, { Trisegment_2_ptr lTrisegment = CreateTrisegment(lTriedge,aLNode,aRNode); - CGAL_STSKEL_BUILDER_TRACE(4, "\n[] Considering E" << lTrisegment->e0().mID << " E" << lTrisegment->e1().mID << " E" << lTrisegment->e2().mID - << " Collinearity: " << trisegment_collinearity_to_string(lTrisegment->collinearity()) ) ; + CGAL_STSKEL_BUILDER_TRACE(4, "\nConsidering " << lTrisegment ) ; // The 02 collinearity configuration is problematic: 01 or 12 collinearity has a seed position // giving the point through which the bisector passes. However, for 02, it is not a given. @@ -362,9 +361,9 @@ void Straight_skeleton_builder_2::CollectNewEvents( Vertex_handle aNode CollectSplitEvents(aNode, aPrevEventTriedge) ; EventPtr lLEdgeEvent = FindEdgeEvent( lPrev , aNode, aPrevEventTriedge ) ; - CGAL_STSKEL_BUILDER_TRACE(2, "Done Left " << (lLEdgeEvent ? "Found" : "Not Found")); + CGAL_STSKEL_BUILDER_TRACE(2, "Done; Edge Event on the Left: " << (lLEdgeEvent ? "Found" : "Not Found")); EventPtr lREdgeEvent = FindEdgeEvent( aNode , lNext, aPrevEventTriedge ) ; - CGAL_STSKEL_BUILDER_TRACE(2, "Done Right " << (lREdgeEvent ? "Found" : "Not Found")); + CGAL_STSKEL_BUILDER_TRACE(2, "Done; Edge Event on the Right: " << (lREdgeEvent ? "Found" : "Not Found")); bool lAcceptL = !!lLEdgeEvent ; bool lAcceptR = !!lREdgeEvent ; @@ -672,8 +671,9 @@ void Straight_skeleton_builder_2::HarmonizeSpeeds(boost::mpl::bool_vertex()->point(), lLH->opposite()->vertex()->point(), - lRH->vertex()->point()) == EQUAL ) + lRH->vertex()->point()) == EQUAL ) { return false; // collinear + } // parallel but not collinear, order arbitrarily (but consistently) return K().less_xy_2_object()(lLH->vertex()->point(), lRH->vertex()->point()) ; @@ -699,6 +699,9 @@ void Straight_skeleton_builder_2::HarmonizeSpeeds(boost::mpl::bool_id() << " with " << (*lRes.first)->id() ) ; mTraits.InitializeLineCoeffs(lBorder->id(), (*lRes.first)->id()); + if(lBorder->weight() != (*lRes.first)->weight()) { + throw std::runtime_error("Collinear input edges must have the same weight"); + } } else { @@ -801,7 +804,7 @@ Straight_skeleton_builder_2::LookupOnSLAV ( Halfedge_handle aBorder, Ev CGAL_STSKEL_BUILDER_TRACE ( 3, "Split point found at the " << ( rSite == AT_SOURCE ? "SOURCE vertex" : ( rSite == AT_TARGET ? "TARGET vertex" : "strict inside" ) ) - << " of the offset edge " << vh2str(lPrevN) << " " << vh2str(v) + << " of the offset edge " << vh2str(lPrevN) << " -- " << vh2str(v) ) ; break ; } @@ -1578,7 +1581,7 @@ void Straight_skeleton_builder_2::Propagate() { EventPtr event = mpq.top(); mpq.pop(); - CGAL_STSKEL_BUILDER_TRACE(4, *event); + CGAL_STSKEL_BUILDER_TRACE(4, "MPQ Event: " << *event); } CGAL_STSKEL_BUILDER_TRACE(4, "END MAIN QUEUE --------------------------------------------- "); #endif @@ -1606,13 +1609,18 @@ void Straight_skeleton_builder_2::Propagate() } ++ mStepID ; + + CGAL::draw(*mSSkel); } else { CGAL_STSKEL_BUILDER_TRACE (3,"\nAlready processed") ; } } - else break ; + else + { + break ; + } } mVisitor.on_propagation_finished(); @@ -2410,6 +2418,8 @@ bool Straight_skeleton_builder_2::FinishUp() mVisitor.on_cleanup_started(); + CGAL::draw(*mSSkel); + std::for_each( mSplitNodes.begin() ,mSplitNodes.end () ,[this](Vertex_handle_pair p){ this->MergeSplitNodes(p); } @@ -2433,6 +2443,8 @@ bool Straight_skeleton_builder_2::FinishUp() mVisitor.on_cleanup_finished(); + CGAL::draw(*mSSkel); + // If 'mMaxTime' is sufficiently large, it will be a full skeleton and could be validated as such... if(mMaxTime) // might be a partial skeleton return mSSkel->is_valid(true); diff --git a/Straight_skeleton_2/include/CGAL/Straight_skeleton_builder_2.h b/Straight_skeleton_2/include/CGAL/Straight_skeleton_builder_2.h index bcc7d2b65f4e..a4845520d706 100644 --- a/Straight_skeleton_2/include/CGAL/Straight_skeleton_builder_2.h +++ b/Straight_skeleton_2/include/CGAL/Straight_skeleton_builder_2.h @@ -452,6 +452,8 @@ private : // @todo Should split events always have lower priority than edge events? Comparison_result CompareEventsSupportAngles ( EventPtr const& aA, EventPtr const& aB ) { + CGAL_STSKEL_BUILDER_TRACE(4, "Compare event support angles " << aA->triedge() << " and " << aB->triedge()); + CGAL_precondition ( aA->type() != Event::cEdgeEvent && aB->type() != Event::cEdgeEvent ) ; if(aA->triedge() == aB->triedge()) @@ -503,9 +505,62 @@ private : public: Event_compare ( Self* aBuilder ) : mBuilder(aBuilder) {} + // Priority queue comparison: `A` has higher priority than `B` + // if `operator()(A, B)` is `false`. bool operator() ( EventPtr const& aA, EventPtr const& aB ) const { - return mBuilder->CompareEvents(aA,aB) == LARGER ; + CGAL_STSKEL_BUILDER_TRACE(16, "operator()(" << aA->triedge() << "; " << aB->triedge() << ")"); + + Comparison_result res = mBuilder->CompareEvents(aA,aB); + if(res == EQUAL) + { + CGAL_STSKEL_BUILDER_TRACE(16, "Events have equal time, find earliest seed") ; + CGAL_STSKEL_BUILDER_TRACE(16, " N" << aA->seed0()->id() << " at " << aA->seed0()->time()) ; + CGAL_STSKEL_BUILDER_TRACE(16, " N" << aA->seed1()->id() << " at " << aA->seed1()->time()) ; + CGAL_STSKEL_BUILDER_TRACE(16, " N" << aB->seed0()->id() << " at " << aB->seed0()->time()) ; + CGAL_STSKEL_BUILDER_TRACE(16, " N" << aB->seed1()->id() << " at " << aB->seed1()->time()) ; + + int aContourCounterA = int(aA->seed0()->is_contour()) + int(aA->seed1()->is_contour()) ; + int aContourCounterB = int(aB->seed0()->is_contour()) + int(aB->seed1()->is_contour()) ; + + // 'A' has more contour seeds => 'A' has higher prio => return false + if (aContourCounterA > aContourCounterB) return false ; + if (aContourCounterA < aContourCounterB) return true ; + + Trisegment_2_ptr lTriA, lTriB ; + if (aContourCounterA == 2) // all 4 seeds are contour vertices + { + return aA->seed0()->id() < aB->seed0()->id(); // for determinism + } + else if (aContourCounterA == 1) // compare the non-contour seeds + { + lTriA = aA->seed0()->is_contour() ? aA->seed1()->trisegment() : aA->seed0()->trisegment() ; + lTriB = aB->seed0()->is_contour() ? aB->seed1()->trisegment() : aB->seed0()->trisegment() ; + } + else // no contour vertices, keep for each event the earliest seed + { + lTriA = aA->seed0()->trisegment(); + if (aA->seed0() != aA->seed1()) + { + if (mBuilder->CompareEvents(aA->seed1()->trisegment(), lTriA) == SMALLER) + lTriA = aA->seed1()->trisegment() ; + } + + lTriB = aB->seed0()->trisegment(); + if (aB->seed0() != aB->seed1()) + { + if (mBuilder->CompareEvents(aB->seed1()->trisegment(), lTriB) == SMALLER) + lTriB = aB->seed1()->trisegment() ; + } + } + + // smmallest has priority + return mBuilder->CompareEvents(lTriA, lTriB) == LARGER ; + } + else + { + return res == LARGER ; + } } private: @@ -527,12 +582,7 @@ private : CGAL_precondition( aA->type() != Event::cEdgeEvent || aB->type() != Event::cEdgeEvent ) ; if ( ! mBuilder->AreEventsSimultaneous(aA,aB) ) - { - Comparison_result res = mBuilder->CompareEvents(aA,aB); - if ( res == EQUAL ) - return aA < aB; - return ( res == LARGER ) ; - } + return Event_compare(mBuilder).operator()(aA, aB) ; // @todo cache // There are simultaneous events, we will need to refresh the queue before calling top() // see PopNextSplitEvent() @@ -543,7 +593,7 @@ private : // i.e. `true` if the angle is larger Comparison_result res = mBuilder->CompareEventsSupportAngles(aA, aB); if ( res == EQUAL ) - return aA < aB; + return Event_compare(mBuilder).operator()(aA, aB) ; // @todo cache return ( res == LARGER ) ; } @@ -709,8 +759,6 @@ private : CreateSegment(aTriedge.e2()), aTriedge.e2()->weight() ); - CGAL_STSKEL_BUILDER_TRACE(5,"Trisegment for " << aTriedge << ":\n" << r ) ; - // Consecutive collinear segments must not have the same weight CGAL_assertion_code(if(r->collinearity() == TRISEGMENT_COLLINEARITY_01)) CGAL_assertion_code(if(aTriedge.e0()->weight() != aTriedge.e1()->weight()) {) @@ -854,7 +902,7 @@ private : void AddSplitEvent ( Vertex_handle aV, EventPtr const& aEvent ) { - CGAL_STSKEL_BUILDER_TRACE(2, "V" << aV->id() << " PQ: " << *aEvent); + CGAL_STSKEL_BUILDER_TRACE(2, "V" << aV->id() << " Split PQ: push " << *aEvent); GetVertexData(aV).mSplitEvents.push(aEvent); } @@ -891,6 +939,19 @@ private : Split_event_compare(this, aV)); } +#ifdef CGAL_SLS_PRINT_QUEUE_BEFORE_EACH_POP + CGAL_STSKEL_BUILDER_TRACE(4, "SPLIT QUEUE OF " << aV->id() << " -------------------------------------------------- "); + CGAL_STSKEL_BUILDER_TRACE(4, "Queue size: " << lPQ.size()); + auto spq = lPQ; + while(!spq.empty()) + { + EventPtr event = spq.top(); + CGAL_STSKEL_BUILDER_TRACE(4, "SPQ Event: " << *event); + spq.pop(); + } + CGAL_STSKEL_BUILDER_TRACE(4, "END SPLIT QUEUE OF " << aV->id() << " --------------------------------------------- "); +#endif + rEvent = lPQ.top(); lPQ.pop(); lData.mNextSplitEventInMainPQ = true ; @@ -949,7 +1010,7 @@ private : Comparison_result rResult = aA->triedge() != aB->triedge() ? CompareEvents( aA->trisegment(), aB->trisegment() ) : EQUAL; - CGAL_STSKEL_BUILDER_TRACE(3, "Compare events " << aA->triedge() << " and " << aB->triedge() << " -> " << rResult); + CGAL_STSKEL_BUILDER_TRACE(4, "Compare events " << aA->triedge() << " and " << aB->triedge() << " -> " << rResult); return rResult; } diff --git a/Straight_skeleton_2/include/CGAL/Straight_skeleton_builder_traits_2.h b/Straight_skeleton_2/include/CGAL/Straight_skeleton_builder_traits_2.h index 50266d6a27d5..97b038ebcff7 100644 --- a/Straight_skeleton_2/include/CGAL/Straight_skeleton_builder_traits_2.h +++ b/Straight_skeleton_2/include/CGAL/Straight_skeleton_builder_traits_2.h @@ -468,12 +468,14 @@ class Straight_skeleton_builder_traits_2_implvertex()->point(), lHR->id()); - boost::optional< Line_2 > lL = CGAL_SS_i::compute_weighted_line_coeffC2(lSL, lHL->weight(), mCaches); - boost::optional< Line_2 > lR = CGAL_SS_i::compute_weighted_line_coeffC2(lSR, lHR->weight(), mCaches); + boost::optional< Line_2 > lL = CGAL_SS_i::compute_normalized_line_coeffC2(lSL, mCaches); + boost::optional< Line_2 > lR = CGAL_SS_i::compute_normalized_line_coeffC2(lSR, mCaches); + if ( !lL || !lR ) + return; Vector_2 lVL( lL->b(), - lL->a()) ; Vector_2 lVR(- lR->b(), lR->a()) ; - Vector_2 lVLR = lVL + lVR ; + Vector_2 lVLR = lHR->weight() * lVL + lHL->weight() * lVR ; const Point_2& laP = aNode->point(); Ray_2 bisect_ray(laP, lVLR) ; @@ -499,11 +501,13 @@ class Straight_skeleton_builder_traits_2_impl lSh (s_h, (*h)->id()); boost::optional< Line_2 > lh = CGAL_SS_i::compute_normalized_line_coeffC2(lSh, mCaches); - FT lLambda = - ( lh->a()*laP.x() + lh->b()*laP.y() + lh->c() ) / - ( lh->a()*lVLR.x() + lh->b()*lVLR.y() ) ; + const FT& aL = lL->a(); const FT& bL = lL->b(); const FT& cL = lL->c(); const FT& wL = lHL->weight(); + const FT& aR = lR->a(); const FT& bR = lR->b(); const FT& cR = lR->c(); const FT& wR = lHR->weight(); + const FT& ah = lh->a(); const FT& bh = lh->b(); const FT& ch = lh->c(); - Point_2 lP = laP + lVLR; - FT lBound = lLambda * ( lL->a()*lP.x() + lL->b()*lP.y() + lL->c() ) ; + FT num = -((bh*cR - bR*ch)*aL - (ah*cR - aR*ch)*bL + (ah*bR - aR*bh)*cL) ; + FT den = ah*bL*wR - aL*bh*wR - (ah*bR - aR*bh)*wL ; + FT lBound = num / den ; if(!is_finite(lBound) || !is_positive(lBound)) continue; @@ -781,16 +785,18 @@ class Straight_skeleton_builder_traits_2_implvertex()->point()), lHR->id()); - boost::optional lL = CGAL_SS_i::compute_weighted_line_coeffC2(lSL, lToFiltered(lHL->weight()), mApproximate_traits.mCaches); - boost::optional lR = CGAL_SS_i::compute_weighted_line_coeffC2(lSR, lToFiltered(lHR->weight()), mApproximate_traits.mCaches); - - Target_Point_2 laP = lToFiltered(aNode->point()); + boost::optional lL = CGAL_SS_i::compute_normalized_line_coeffC2(lSL, mApproximate_traits.mCaches); + boost::optional lR = CGAL_SS_i::compute_normalized_line_coeffC2(lSR, mApproximate_traits.mCaches); + if ( !lL || !lR ) + return; // These are weighted direction of the lines supporting the contour segments. // Coefficients for SR are negated because aNode is at the "source" of SR. Target_Vector_2 lVL( lL->b(), - lL->a()) ; Target_Vector_2 lVR(- lR->b(), lR->a()) ; - Target_Vector_2 lVLR = lVL + lVR ; + Target_Vector_2 lVLR = lHR->weight() * lVL + lHL->weight() * lVR ; + + Target_Point_2 laP = lToFiltered(aNode->point()); Target_Ray_2 bisect_ray(laP, lVLR) ; // @todo this should use some kind of spatial searching @@ -815,35 +821,22 @@ class Straight_skeleton_builder_traits_2_implid()); boost::optional lh = CGAL_SS_i::compute_normalized_line_coeffC2(lSh, mApproximate_traits.mCaches); - Target_FT lLambda = - ( lh->a()*laP.x() + lh->b()*laP.y() + lh->c() ) / - ( lh->a()*lVLR.x() + lh->b()*lVLR.y() ) ; + const Target_FT& aL = lL->a(); const Target_FT& bL = lL->b(); const Target_FT& cL = lL->c(); const Target_FT& wL = lToFiltered(lHL->weight()); + const Target_FT& aR = lR->a(); const Target_FT& bR = lR->b(); const Target_FT& cR = lR->c(); const Target_FT& wR = lToFiltered(lHR->weight()); + const Target_FT& ah = lh->a(); const Target_FT& bh = lh->b(); const Target_FT& ch = lh->c(); + + // [[ x == (bh*c0*w1 - b0*ch*w1 - (bh*c1 - b1*ch)*w0)/(ah*b0*w1 - a0*bh*w1 - (ah*b1 - a1*bh)*w0), + // y == -(ah*c0*w1 - a0*ch*w1 - (ah*c1 - a1*ch)*w0)/(ah*b0*w1 - a0*bh*w1 - (ah*b1 - a1*bh)*w0), + // t == -((bh*c1 - b1*ch)*a0 - (ah*c1 - a1*ch)*b0 + (ah*b1 - a1*bh)*c0)/(ah*b0*w1 - a0*bh*w1 - (ah*b1 - a1*bh)*w0) ]] - // Scale it because |d0 + d1| doesn't send aNode to the t=1 line, but to t=lL(aNode+d0+d1)=lR(aNode+d0+d1) - Target_Point_2 lP = laP + lVLR; - Target_FT lBound = lLambda * ( lL->a()*lP.x() + lL->b()*lP.y() + lL->c() ) ; + Target_FT num = -((bh*cR - bR*ch)*aL - (ah*cR - aR*ch)*bL + (ah*bR - aR*bh)*cL) ; + Target_FT den = ah*bL*wR - aL*bh*wR - (ah*bR - aR*bh)*wL ; + Target_FT lBound = num / den ; -#if 0 +#if 1 std::cout << "E" << (*h)->id() << " s_h = " << s_h << std::endl; std::cout << "left/right E" << lHL->id() << " E" << lHR->id() << std::endl; std::cout << "V" << aNode->id() << std::endl; @@ -858,26 +851,17 @@ class Straight_skeleton_builder_traits_2_impl lh1 = CGAL_SS_i::compute_weighted_line_coeffC2(lSR, lToFiltered(lHR->weight()), mApproximate_traits.mCaches); + std::cout << "lL check" << square(lL->a()) + square(lL->b()) << std::endl; + std::cout << "lR check" << square(lR->a()) + square(lR->b()) << std::endl; - std::cout << "lh0 check" << square(lh0->a()) + square(lh0->b()) << std::endl; - std::cout << "lh1 check" << square(lh1->a()) + square(lh1->b()) << std::endl; - std::cout << "l0 time at aNode: " << lh0->a()*laP.x() + lh0->b()*laP.y() + lh0->c() << std::endl; - std::cout << "l1 time at aNode: " << lh1->a()*laP.x() + lh1->b()*laP.y() + lh1->c() << std::endl; - std::cout << "l0 time at aNode + lVLR: " << lh0->a()*(laP + lVLR).x() + lh0->b()*(laP + lVLR).y() + lh0->c() << std::endl; - std::cout << "l1 time at aNode + lVLR: " << lh1->a()*(laP + lVLR).x() + lh1->b()*(laP + lVLR).y() + lh1->c() << std::endl; + std::cout << "lBound: " << lBound << std::endl; auto ipp = FK().intersect_2_object()(s_h, bisect_ray); Target_Point_2* ip = boost::get(&*ipp); - std::cout << "l0 time at inter pt: " << lh0->a()*ip->x() + lh0->b()*ip->y() + lh0->c() << std::endl; - std::cout << "l1 time at inter pt: " << lh1->a()*ip->x() + lh1->b()*ip->y() + lh1->c() << std::endl; - std::cout << "lh-> " << lh->a() << " " << lh->b() << " " << square(lh->a()) + square(lh->b()) << std::endl; - std::cout << "lh time at inter pt: " << lh->a()*ip->x() + lh->b()*ip->y() + lh->c() << std::endl; - std::cout << "lLambda: " << lLambda << std::endl; - std::cout << "lBound: " << lBound << std::endl; - CGAL_assertion(lBound == lh0->a()*ip->x() + lh0->b()*ip->y() + lh0->c()); + std::cout << "l0 time at inter pt: " << (lL->a()*ip->x() + lL->b()*ip->y() + lL->c()) / wL << std::endl; + std::cout << "l1 time at inter pt: " << (lR->a()*ip->x() + lR->b()*ip->y() + lR->c()) / wR << std::endl; + std::cout << "lh at inter pt: " << lh->a()*ip->x() + lh->b()*ip->y() + lh->c() << std::endl; #endif if(!is_finite(lBound) || !is_positive(lBound)) diff --git a/Straight_skeleton_2/include/CGAL/Trisegment_2.h b/Straight_skeleton_2/include/CGAL/Trisegment_2.h index ce7fa34a79c0..aae503b18796 100644 --- a/Straight_skeleton_2/include/CGAL/Trisegment_2.h +++ b/Straight_skeleton_2/include/CGAL/Trisegment_2.h @@ -194,12 +194,14 @@ class Trisegment_2 { const std::string lPadding = std::string(2 * aDepth, ' '); - os << lPadding << "[&: " << &aTri << " ID: " << aTri.id() << "\n" - << lPadding << "\tE" << aTri.e0().mID << " E" << aTri.e1().mID << " E" << aTri.e2().mID << "\n" - << lPadding << "\t" << s2str(aTri.e0()) << " w = " << n2str(aTri.w0()) << ";" << "\n" - << lPadding << "\t" << s2str(aTri.e1()) << " w = " << n2str(aTri.w1()) << ";" << "\n" - << lPadding << "\t" << s2str(aTri.e2()) << " w = " << n2str(aTri.w2()) << ";" << "\n" - << lPadding << "\tCollinearity: " << trisegment_collinearity_to_string(aTri.collinearity()) << "\n" + os << lPadding << "[\n" + << lPadding << " ID: " << aTri.id() << "\n" + << lPadding << " E" << aTri.e0().mID << " E" << aTri.e1().mID << " E" << aTri.e2().mID << "\n" + << lPadding << " " << s2str(aTri.e0()) << " w = " << n2str(aTri.w0()) << ";" << "\n" + << lPadding << " " << s2str(aTri.e1()) << " w = " << n2str(aTri.w1()) << ";" << "\n" + << lPadding << " " << s2str(aTri.e2()) << " w = " << n2str(aTri.w2()) << ";" << "\n" + << lPadding << " Collinearity: " << trisegment_collinearity_to_string(aTri.collinearity()) << "\n" + << lPadding << " Seeds: " << bool(aTri.mChildL) << " " << bool(aTri.mChildR) << " " << bool(aTri.mChildT) << "\n" << lPadding << "]\n" << std::flush; } diff --git a/Straight_skeleton_2/include/CGAL/constructions/Polygon_offset_cons_ftC2.h b/Straight_skeleton_2/include/CGAL/constructions/Polygon_offset_cons_ftC2.h index 78ff0b0368e6..17f9d480a358 100644 --- a/Straight_skeleton_2/include/CGAL/constructions/Polygon_offset_cons_ftC2.h +++ b/Straight_skeleton_2/include/CGAL/constructions/Polygon_offset_cons_ftC2.h @@ -68,8 +68,8 @@ construct_offset_pointC2 ( typename K::FT const& t, { if ( ! CGAL_NTS is_zero(den) ) { - FT numX = t * l1->b() / w0 - t * l0->b() / w1 + l0->b() * l1->c() - l1->b() * l0->c() ; - FT numY = t * l1->a() / w0 - t * l0->a() / w1 + l0->a() * l1->c() - l1->a() * l0->c() ; + FT numX = w0*t*l1->b() - w1*t*l0->b() + l0->b()*l1->c() - l1->b()*l0->c() ; + FT numY = w0*t*l1->a() - w1*t*l0->a() + l0->a()*l1->c() - l1->a()*l0->c() ; x = -numX / den ; y = numY / den ; @@ -92,11 +92,8 @@ construct_offset_pointC2 ( typename K::FT const& t, CGAL_STSKEL_TRAITS_TRACE(" Projected seed point: (" << px << "," << py << ")" ) ; - // When reformulating the progression of the front using the line orthogonal - // to the input segment, the speed (weight) becomes inverted: a large weight - // in the coefficients of the line means the front moves slower - x = px + t * l0->a() / w0 ; - y = py + t * l0->b() / w0 ; + x = px + w0 * t * l0->a() ; + y = py + w0 * t * l0->b() ; ok = CGAL_NTS is_finite(x) && CGAL_NTS is_finite(y) ; } diff --git a/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h b/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h index b968c287195f..cd430a70fcec 100644 --- a/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h +++ b/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h @@ -54,7 +54,7 @@ template Trisegment_collinearity trisegment_collinearity_no_exact_constructions(const Segment_2_with_ID& e0, const Segment_2_with_ID& e1, const Segment_2_with_ID& e2, - Caches& caches ) + Caches& caches) { // 'are_edges_orderly_collinear()' is used to harmonize coefficients, but if the kernel is inexact // we could also have that are_edges_orderly_collinear() returns false, but the computed coefficients @@ -68,7 +68,8 @@ Trisegment_collinearity trisegment_collinearity_no_exact_constructions(const Seg bool is_02 = (l0->a() == l2->a()) && (l0->b() == l2->b()) && (l0->c() == l2->c()); bool is_12 = (l1->a() == l2->a()) && (l1->b() == l2->b()) && (l1->c() == l2->c()); - CGAL_STSKEL_TRAITS_TRACE("coeff equalities: " << is_01 << " " << is_02 << " " << is_12) + CGAL_STSKEL_TRAITS_TRACE("collinearity: [" << e0.id() << ", " << e1.id() << ", " << e2.id() + << "]: " << is_01 << " " << is_02 << " " << is_12); if ( is_01 & !is_02 & !is_12 ) return TRISEGMENT_COLLINEARITY_01; @@ -248,37 +249,6 @@ compute_normalized_line_coeffC2( Segment_2_with_ID const& e, return rRes ; } -// @todo weightless coefficients are stored because we use them sometimes weighted, and sometimes -// inversely weighted (filtering bound). Should we store them weighted also for speed reasons? -template -boost::optional< typename K::Line_2 > compute_weighted_line_coeffC2( Segment_2_with_ID const& e, - typename K::FT const& aWeight, - Caches& aCaches ) -{ - typedef typename K::FT FT ; - typedef typename K::Line_2 Line_2 ; - - CGAL_precondition( CGAL_NTS is_finite(aWeight) && CGAL_NTS is_positive(aWeight) ) ; - - boost::optional< Line_2 > l = compute_normalized_line_coeffC2(e, aCaches); - if( ! l ) - return boost::none ; - - FT a = l->a() * aWeight ; - FT b = l->b() * aWeight ; - FT c = l->c() * aWeight ; - - CGAL_STSKEL_TRAITS_TRACE("\n~~ Weighted line coefficients for E" << e.id() << " " << s2str(e) - << "\nweight=" << n2str(aWeight) - << "\na="<< n2str(a) << "\nb=" << n2str(b) << "\nc=" << n2str(c) - ) ; - - if ( !CGAL_NTS is_finite(a) || !CGAL_NTS is_finite(b) || !CGAL_NTS is_finite(c) ) - return boost::none; - - return cgal_make_optional( K().construct_line_2_object()(a,b,c) ) ; -} - template Rational squared_distance_from_point_to_lineC2( FT const& px, FT const& py, FT const& sx, FT const& sy, FT const& tx, FT const& ty ) { @@ -312,9 +282,6 @@ construct_trisegment ( Segment_2_with_ID const& e0, typedef Trisegment_2 >Trisegment_2 ; typedef typename Trisegment_2::Self_ptr Trisegment_2_ptr ; - CGAL_STSKEL_TRAITS_TRACE("\n~~ Construct trisegment [" << typeid(typename K::FT).name() << "]"); - CGAL_STSKEL_TRAITS_TRACE("Segments E" << e0.id() << " E" << e1.id() << " E" << e2.id()); - Trisegment_collinearity lCollinearity = trisegment_collinearity_no_exact_constructions(e0,e1,e2,aCaches); return Trisegment_2_ptr( new Trisegment_2(e0, w0, e1, w1, e2, w2, lCollinearity, id) ) ; @@ -347,8 +314,7 @@ compute_normal_offset_lines_isec_timeC2 ( Trisegment_2_ptr< Trisegment_2 Optional_line_2 ; - CGAL_STSKEL_TRAITS_TRACE("\n~~ Computing normal offset lines isec time [" << typeid(FT).name() << "]") ; - CGAL_STSKEL_TRAITS_TRACE("Event:\n" << tri); + CGAL_STSKEL_TRAITS_TRACE("\ncompute_normal_offset_lines_isec_timeC2(" << tri->id() << ") [" << typeid(FT).name() << "]") ; FT num(0), den(0) ; @@ -356,55 +322,108 @@ compute_normal_offset_lines_isec_timeC2 ( Trisegment_2_ptr< Trisegment_2 0' being to the left of the line. + // with 't > 0' being to the left of the line. // If 3 such offset lines intersect at the same offset distance, the intersection 't', // or 'time', can be computed solving for 't' in the linear system formed by 3 such equations. - // The result is : // - // sage: var('a0 b0 c0 a1 b1 c1 a2 b2 c2 x y t w0 w1 w2') - // (a0, b0, c0, a1, b1, c1, a2, b2, c2, x, y, t, w0, w1, w2) - // sage: - // sage: eqw0 = w0*a0*x + w0*b0*y + w0*c0 - t == 0 - // sage: eqw1 = w1*a1*x + w1*b1*y + w1*c1 - t == 0 - // sage: eqw2 = w2*a2*x + w2*b2*y + w2*c2 - t == 0 - // sage: - // sage: solve([eqw0,eqw1,eqw2], x,y,t) - // x == (((c1*w1 - c2*w2)*b0 - (b1*w1 - b2*w2)*c0)*w0 - (b2*c1*w2 - b1*c2*w2)*w1) / (((b1*w1 - b2*w2)*a0 - (a1*w1 - a2*w2)*b0)*w0 - (a2*b1*w2 - a1*b2*w2)*w1), - // y == -(((c1*w1 - c2*w2)*a0 - (a1*w1 - a2*w2)*c0)*w0 - (a2*c1*w2 - a1*c2*w2)*w1) / (((b1*w1 - b2*w2)*a0 - (a1*w1 - a2*w2)*b0)*w0 - (a2*b1*w2 - a1*b2*w2)*w1), - // t == -((b2*c1*w2 - b1*c2*w2)*a0*w1 - (a2*c1*w2 - a1*c2*w2)*b0*w1 + (a2*b1*w2 - a1*b2*w2)*c0*w1)*w0/(((b1*w1 - b2*w2)*a0 - (a1*w1 - a2*w2)*b0)*w0 - (a2*b1*w2 - a1*b2*w2)*w1) + // [[ t == ((b2*c1 - b1*c2)*a0 - (a2*c1 - a1*c2)*b0 + (a2*b1 - a1*b2)*c0)/((b2*w1 - b1*w2)*a0 - (a2*w1 - a1*w2)*b0 + (a2*b1 - a1*b2)*w0) ]] bool ok = false ; - Optional_line_2 l0 = compute_weighted_line_coeffC2(tri->e0(), tri->w0(), aCaches) ; - Optional_line_2 l1 = compute_weighted_line_coeffC2(tri->e1(), tri->w1(), aCaches) ; - Optional_line_2 l2 = compute_weighted_line_coeffC2(tri->e2(), tri->w2(), aCaches) ; + Optional_line_2 l0 = compute_normalized_line_coeffC2(tri->e0(), aCaches) ; + Optional_line_2 l1 = compute_normalized_line_coeffC2(tri->e1(), aCaches) ; + Optional_line_2 l2 = compute_normalized_line_coeffC2(tri->e2(), aCaches) ; if ( l0 && l1 && l2 ) { - CGAL_STSKEL_TRAITS_TRACE("coeffs E" << tri->e0().id() << " [" << n2str(l0->a()) << "; " << n2str(l0->b()) << "; " << n2str(l0->c()) << "]" - << "\ncoeffs E" << tri->e1().id() << " [" << n2str(l1->a()) << "; " << n2str(l1->b()) << "; " << n2str(l1->c()) << "]" - << "\ncoeffs E" << tri->e2().id() << " [" << n2str(l2->a()) << "; " << n2str(l2->b()) << "; " << n2str(l2->c()) << "]") ; - - num = (l2->a()*l0->b()*l1->c()) - -(l2->a()*l1->b()*l0->c()) - -(l2->b()*l0->a()*l1->c()) - +(l2->b()*l1->a()*l0->c()) - +(l1->b()*l0->a()*l2->c()) - -(l0->b()*l1->a()*l2->c()); - - den = (-l2->a()*l1->b()) - +( l2->a()*l0->b()) - +( l2->b()*l1->a()) - -( l2->b()*l0->a()) - +( l1->b()*l0->a()) - -( l0->b()*l1->a()); + const FT& a0 = l0->a(); const FT& b0 = l0->b(); const FT& c0 = l0->c(); const FT& w0 = tri->w0(); + const FT& a1 = l1->a(); const FT& b1 = l1->b(); const FT& c1 = l1->c(); const FT& w1 = tri->w1(); + const FT& a2 = l2->a(); const FT& b2 = l2->b(); const FT& c2 = l2->c(); const FT& w2 = tri->w2(); + + CGAL_STSKEL_TRAITS_TRACE("coeffs E" << tri->e0().id() << " [" << n2str(a0) << "; " << n2str(b0) << "; " << n2str(c0) << "]" + << "\ncoeffs E" << tri->e1().id() << " [" << n2str(a1) << "; " << n2str(b1) << "; " << n2str(c1) << "]" + << "\ncoeffs E" << tri->e2().id() << " [" << n2str(a2) << "; " << n2str(b2) << "; " << n2str(c2) << "]") ; + CGAL_STSKEL_TRAITS_TRACE("weight E" << tri->e0().id() << " [" << n2str(w0) << "] " + << "weight E" << tri->e1().id() << " [" << n2str(w1) << "] " + << "weight E" << tri->e2().id() << " [" << n2str(w2) << "]") ; + + num = (b2*c1 - b1*c2)*a0 - (a2*c1 - a1*c2)*b0 + (a2*b1 - a1*b2)*c0 ; + den = (b2*w1 - b1*w2)*a0 - (a2*w1 - a1*w2)*b0 + (a2*b1 - a1*b2)*w0 ; + + // Here we are in a pickle: lines are not collinear but with weights, we could + // get an infinity of solutions still (particularly easy to create with null weights) + if (CGAL_NTS certified_is_zero(den)) { + // If the numerator (which is a 2nd determinant within the augmented matrix) is also 0, + // then we have a consistent system, so an infinity of solution for a fixed value of 't' + if (CGAL_NTS certified_is_zero(num)) { + CGAL_STSKEL_TRAITS_TRACE("Infinite solutions"); + + const auto& e0 = tri->e0(); + const auto& e1 = tri->e1(); + const auto& e2 = tri->e2(); + + CGAL_assertion(are_edges_parallelC2(e0, e1)); + CGAL_assertion(are_edges_parallelC2(e1, e2)); + + // System is consistent, extract t + bool w0_is_zero = CGAL_NTS certified_is_zero(w0); + bool w1_is_zero = CGAL_NTS certified_is_zero(w1); + bool w2_is_zero = CGAL_NTS certified_is_zero(w2); + + // If the intersection is a line, that would mean that all lines are collinear + // and there would be two lines in the same direction, so edges would be collinear + // and we would not be here + CGAL_assertion(!(w0_is_zero && w1_is_zero && w2_is_zero)); + + // If at least one line is not moving, evaluate a moving line at a point of the static line. + if (w0_is_zero && w1_is_zero) { + CGAL_assertion(a0 * a1 < 0); + num = a2*e0.source().x() + b2*e0.source().y() + c2; + den = w0; + } else if (w0_is_zero && w2_is_zero) { + CGAL_assertion(a0 * a2 < 0); + num = a1*e0.source().x() + b1*e0.source().y() + c1; + den = w1; + } else if (w1_is_zero && w2_is_zero) { + CGAL_assertion(a1 * a2 < 0); + num = a0*e1.source().x() + b0*e1.source().y() + c0; + den = w0; + } else if (w0_is_zero) { + num = a2*e0.source().x() + b2*e0.source().y() + c2; + den = w2; + } else if (w1_is_zero) { + num = a1*e0.source().x() + b1*e0.source().y() + c1; + den = w1; + } else if (w2_is_zero) { + num = a0*e1.source().x() + b0*e1.source().y() + c0; + den = w0; + } else { + // We can't have two lines moving in the same direction with the same speeds, + // otherwise they are either identical and we would be collinear (contradiction), + // or there is no solution to the base linear system (also contradiction). + // So if there is the same speed, they move in opposite directions and the only possibility + // is at t=0 + if (CGAL_NTS certified_is_zero(w0 - w1)) { + num = 0; + den = 1; + } else { + num = c1 - c0; + den = w0 - w1; + } + } + } else { + CGAL_assertion_msg(false, "Unexpected case: lines are not collinear, but no valid 't' could be computed."); + } + } else { + CGAL_STSKEL_TRAITS_TRACE("No solution exists (inconsistent system)."); + } ok = CGAL_NTS is_finite(num) && CGAL_NTS is_finite(den); } - CGAL_STSKEL_TRAITS_TRACE("Event time (normal): n=" << num << " d=" << den << " n/d=" << Rational(num,den) ) + CGAL_STSKEL_TRAITS_TRACE("Event time (normal): n=" << num << " d=" << den << " n/d=" << Rational(num,den)); return cgal_make_optional(ok,Rational(num,den)) ; } @@ -423,8 +442,8 @@ compute_oriented_midpoint ( Segment_2_with_ID const& e0, typedef typename K::FT FT ; CGAL_STSKEL_TRAITS_TRACE("Computing oriented midpoint between:" - << "\ne0: " << s2str(e0) - << "\ne1: " << s2str(e1) + << "\n\te0: " << s2str(e0) + << "\n\te1: " << s2str(e1) ); FT delta01 = CGAL::squared_distance(e0.target(),e1.source()); @@ -549,18 +568,17 @@ compute_artifical_isec_timeC2 ( Trisegment_2_ptr< Trisegment_2e0() == tri->e1()); CGAL_precondition(bool(tri->child_l())); - Optional_line_2 l0 = compute_weighted_line_coeffC2(tri->e0(), tri->w0(), aCaches) ; + Optional_line_2 l0 = compute_normalized_line_coeffC2(tri->e0(), aCaches) ; if( !l0 ) return boost::none ; - const Segment_2& contour_seg = tri->e0(); - Direction_2 perp_dir ( contour_seg.source().y() - contour_seg.target().y() , - contour_seg.target().x() - contour_seg.source().x() ) ; boost::optional< typename K::Point_2 > seed = construct_offset_lines_isecC2(tri->child_l(), aCaches ) ; - if(!seed) return boost::none; + const Segment_2& contour_seg = tri->e0(); + const Direction_2 perp_dir(contour_seg.source().y() - contour_seg.target().y(), + contour_seg.target().x() - contour_seg.source().x()); const Ray_2 ray(*seed, perp_dir); const Segment_2& opp_seg = tri->e2(); Boolean inter_exist = K().do_intersect_2_object()(ray, opp_seg); @@ -571,7 +589,6 @@ compute_artifical_isec_timeC2 ( Trisegment_2_ptr< Trisegment_2(&*inter_res)) { // get the segment extremity closest to the seed @@ -588,7 +605,7 @@ compute_artifical_isec_timeC2 ( Trisegment_2_ptr< Trisegment_2(t, FT(1))) ; + return cgal_make_optional(ok, Rational(t, tri->w0())) ; } // Given 3 oriented straight line segments: e0, e1, e2 @@ -601,6 +618,11 @@ compute_artifical_isec_timeC2 ( Trisegment_2_ptr< Trisegment_2 boost::optional< Rational< typename K::FT> > compute_degenerate_offset_lines_isec_timeC2 ( Trisegment_2_ptr< Trisegment_2 > > const& tri, @@ -629,17 +651,17 @@ compute_degenerate_offset_lines_isec_timeC2 ( Trisegment_2_ptr< Trisegment_2collinear_edge(), tri->collinear_edge_weight(), aCaches) ; - Optional_line_2 l1 = compute_weighted_line_coeffC2(tri->other_collinear_edge(), tri->other_collinear_edge_weight(), aCaches) ; - Optional_line_2 l2 = compute_weighted_line_coeffC2(tri->non_collinear_edge(), tri->non_collinear_edge_weight(), aCaches) ; + // else, because l0 is vertical (==> not horizontal), we add the horizontal line: + // + // y(t) = py + w0*t*l0.b + // + // sage: sage: var('a0 b0 c0 a2 b2 c2 x y t w0 w2 px py') + // sage: eq0 = a0*x + b0*y + c0 - w0*t == 0 + // sage: eq2 = a2*x + b2*y + c2 - w2*t == 0 + // sage: eqv = px + w0 * t * a0 == x + // sage: solve([eq0,eq2,eqv], x, y, t) + // [[ x == -(b0*px*w2 - (a0*b2*c0 - a0*b0*c2 + b2*px)*w0)/((a0*a2*b0 - a0^2*b2 + b2)*w0 - b0*w2), + // y == (a0*px*w2 - (a0*a2*c0 - a0^2*c2 + a2*px + c2)*w0 + c0*w2)/((a0*a2*b0 - a0^2*b2 + b2)*w0 - b0*w2), + // t == (a0*b2*px - (a2*px + c2)*b0 + b2*c0)/((a0*a2*b0 - a0^2*b2 + b2)*w0 - b0*w2) ]] + // sage: eqv = py + w0 * t * b0 == y + // sage: solve([eq0,eq2,eqv], x, y, t) + // [[ x == -(b0*py*w2 - (b0*b2*c0 - b0^2*c2 + b2*py + c2)*w0 + c0*w2)/((a2*b0^2 - a0*b0*b2 - a2)*w0 + a0*w2), + // y == (a0*py*w2 - (a2*b0*c0 - a0*b0*c2 + a2*py)*w0)/((a2*b0^2 - a0*b0*b2 - a2)*w0 + a0*w2), + // t == -(a2*b0*py - (b2*py + c2)*a0 + a2*c0)/((a2*b0^2 - a0*b0*b2 - a2)*w0 + a0*w2) ]] + + + Optional_line_2 l0 = compute_normalized_line_coeffC2(tri->collinear_edge(), aCaches) ; + Optional_line_2 l1 = compute_normalized_line_coeffC2(tri->other_collinear_edge(), aCaches) ; + Optional_line_2 l2 = compute_normalized_line_coeffC2(tri->non_collinear_edge(), aCaches) ; Optional_point_2 q = construct_degenerate_seed_pointC2(tri, aCaches); @@ -698,6 +729,8 @@ compute_degenerate_offset_lines_isec_timeC2 ( Trisegment_2_ptr< Trisegment_2collinear_edge_weight() == tri->other_collinear_edge_weight() ) { + const FT& w0 = tri->collinear_edge_weight(); + const FT& w2 = tri->non_collinear_edge_weight(); const FT& l0a = l0->a() ; const FT& l0b = l0->b() ; const FT& l0c = l0->c() ; @@ -705,50 +738,29 @@ compute_degenerate_offset_lines_isec_timeC2 ( Trisegment_2_ptr< Trisegment_2b() ; const FT& l2c = l2->c() ; - // Since l0 and l1 are parallel, we cannot solve the system using: - // l0a*x + l0b*y + l0c = 0 (1) - // l1a*x + l1b*y + l1c = 0 - // l2a*x + l2b*y + l2c = 0 - // Instead, we use the equation of the line orthogonal to l0 (and l1). - // However, rephrasing - // l0a*x + l0b*y + l0c = 0 - // to - // [x, y] = projected_seed + t * N - // requires the norm (l0a² + l0b²) to be exactly '1', which likely isn't the case - // if we are using inexact square roots. In that case, the norm behaves similarly - // to a weight (i.e. speed), and the speed is inverted in the alternate front formulation. - // Equation (1) rewritten with weights is: - // w*l0a*x + w*l0b*y + w*l0c = 0 - // with l0a² + l0b² ~= 1. Extracting the numerical error, we have: - // w'*l0a'*x + w'*l0b'*y + w'*l0c' = 0, - // with l0a'² + l0b'² = 1. - // The orthogonal displacement is rephrased to: - // [x, y] = projected_seed + t / w' * N, - // with w' = weight * (l0a² + l0b²). - const FT sq_w0 = square(l0a) + square(l0b); // l0a and l0b are already *weighted* coefficients - FT num(0), den(0) ; - if ( ! CGAL_NTS is_zero(l0b) ) // Non-vertical + if ( ! CGAL_NTS is_zero(l0b) ) // l0 is not vertical, add the vertical component { - num = ((l2a*l0b - l0a*l2b) * px - l2b*l0c + l0b*l2c) * sq_w0 ; - den = l0a*l0a*l2b - l2b*sq_w0 + l0b*sq_w0 - l0a*l2a*l0b ; + num = l0a*l2b*px - (l2a*px + l2c)*l0b + l2b*l0c ; + den = (l0a*l2a*l0b - l0a*l0a*l2b + l2b)*w0 - l0b*w2 ; CGAL_STSKEL_TRAITS_TRACE("Event time (degenerate, non-vertical) n=" << n2str(num) << " d=" << n2str(den) << " n/d=" << Rational(num,den) ) } - else + else // l0 is vertical, add the horizontal component { - // l0b = 0, and all sq_w0 disappear - num = - l0a*l2b*py - l0a*l2c + l2a*l0c ; - den = l2a - l0a; + num = -(l2a*l0b*py - (l2b*py + l2c)*l0a + l2a*l0c) ; + den = (l2a*l0b*l0b - l0a*l0b*l2b - l2a)*w0 + l0a*w2 ; CGAL_STSKEL_TRAITS_TRACE("Event time (degenerate, vertical) n=" << n2str(num) << " d=" << n2str(den) << " n/d=" << Rational(num,den) ) } - ok = CGAL_NTS is_finite(num) && CGAL_NTS is_finite(den); - return cgal_make_optional(ok, Rational(num,den)) ; + ok = CGAL_NTS is_finite(num) && CGAL_NTS is_finite(den) ; + return cgal_make_optional(ok, Rational(num,den)) ; } else { + // @todo we can't be there, but understand how (where (if)) is parallel + non-collinear handled + // l0 and l1 are collinear but with different speeds, so there cannot be an event. CGAL_STSKEL_TRAITS_TRACE("Event times (degenerate, inequal norms)") CGAL_STSKEL_TRAITS_TRACE("--> Returning 0/0 (no event)"); @@ -812,22 +824,33 @@ construct_normal_offset_lines_isecC2 ( Trisegment_2_ptr< Trisegment_2e0(), tri->w0(), aCaches) ; - Optional_line_2 l1 = compute_weighted_line_coeffC2(tri->e1(), tri->w1(), aCaches) ; - Optional_line_2 l2 = compute_weighted_line_coeffC2(tri->e2(), tri->w2(), aCaches) ; + Optional_line_2 l0 = compute_normalized_line_coeffC2(tri->e0(), aCaches) ; + Optional_line_2 l1 = compute_normalized_line_coeffC2(tri->e1(), aCaches) ; + Optional_line_2 l2 = compute_normalized_line_coeffC2(tri->e2(), aCaches) ; bool ok = false ; if ( l0 && l1 && l2 ) { - FT den = l0->a()*l2->b() - l0->a()*l1->b() - l1->a()*l2->b() + l2->a()*l1->b() + l0->b()*l1->a() - l0->b()*l2->a(); + // [[ x == ((c2*w1 - c1*w2)*b0 - (b2*w1 - b1*w2)*c0 + (b2*c1 - b1*c2)*w0)/((b2*w1 - b1*w2)*a0 - (a2*w1 - a1*w2)*b0 + (a2*b1 - a1*b2)*w0), + // y == -((c2*w1 - c1*w2)*a0 - (a2*w1 - a1*w2)*c0 + (a2*c1 - a1*c2)*w0)/((b2*w1 - b1*w2)*a0 - (a2*w1 - a1*w2)*b0 + (a2*b1 - a1*b2)*w0) ]] + + const FT& a0 = l0->a(); const FT& b0 = l0->b(); const FT& c0 = l0->c(); const FT& w0 = tri->w0(); + const FT& a1 = l1->a(); const FT& b1 = l1->b(); const FT& c1 = l1->c(); const FT& w1 = tri->w1(); + const FT& a2 = l2->a(); const FT& b2 = l2->b(); const FT& c2 = l2->c(); const FT& w2 = tri->w2(); + + CGAL_STSKEL_TRAITS_TRACE("coeffs E" << tri->e0().id() << " [" << n2str(a0) << "; " << n2str(b0) << "; " << n2str(c0) << "] weight [" << n2str(w0) << "]" + << "\ncoeffs E" << tri->e1().id() << " [" << n2str(a1) << "; " << n2str(b1) << "; " << n2str(c1) << "] weight [" << n2str(w1) << "]" + << "\ncoeffs E" << tri->e2().id() << " [" << n2str(a2) << "; " << n2str(b2) << "; " << n2str(c2) << "] weight [" << n2str(w2) << "]") ; + + FT den = (b2*w1 - b1*w2)*a0 - (a2*w1 - a1*w2)*b0 + (a2*b1 - a1*b2)*w0 ; CGAL_STSKEL_TRAITS_TRACE("\tden=" << n2str(den) ) if ( ! CGAL_NTS certified_is_zero(den) ) { - FT numX = l0->b()*l2->c() - l0->b()*l1->c() - l1->b()*l2->c() + l2->b()*l1->c() + l1->b()*l0->c() - l2->b()*l0->c(); - FT numY = l0->a()*l2->c() - l0->a()*l1->c() - l1->a()*l2->c() + l2->a()*l1->c() + l1->a()*l0->c() - l2->a()*l0->c(); + FT numX = (c2*w1 - c1*w2)*b0 - (b2*w1 - b1*w2)*c0 + (b2*c1 - b1*c2)*w0 ; + FT numY = -((c2*w1 - c1*w2)*a0 - (a2*w1 - a1*w2)*c0 + (a2*c1 - a1*c2)*w0) ; CGAL_STSKEL_TRAITS_TRACE("\tnumX=" << n2str(numX) << "\n\tnumY=" << n2str(numY) ) ; @@ -835,8 +858,8 @@ construct_normal_offset_lines_isecC2 ( Trisegment_2_ptr< Trisegment_2e0() == tri->e1()); + CGAL_precondition(!is_zero(tri->w0())); CGAL_precondition(bool(tri->child_l())); const Segment_2& contour_seg = tri->e0(); @@ -931,8 +955,8 @@ construct_degenerate_offset_lines_isecC2 ( Trisegment_2_ptr< Trisegment_2collinear_edge(), tri->collinear_edge_weight(), aCaches) ; - Optional_line_2 l2 = compute_weighted_line_coeffC2(tri->non_collinear_edge(), tri->non_collinear_edge_weight(), aCaches) ; + Optional_line_2 l0 = compute_normalized_line_coeffC2(tri->collinear_edge(), aCaches) ; + Optional_line_2 l2 = compute_normalized_line_coeffC2(tri->non_collinear_edge(), aCaches) ; Optional_point_2 q = construct_degenerate_seed_pointC2(tri, aCaches); @@ -940,68 +964,62 @@ construct_degenerate_offset_lines_isecC2 ( Trisegment_2_ptr< Trisegment_2collinear_edge_weight(), tri->other_collinear_edge_weight()); - if ( res == EQUAL ) + const FT& w0 = tri->collinear_edge_weight(); + const FT& w2 = tri->non_collinear_edge_weight(); + const FT& l0a = l0->a() ; + const FT& l0b = l0->b() ; + const FT& l0c = l0->c() ; + const FT& l2a = l2->a() ; + const FT& l2b = l2->b() ; + const FT& l2c = l2->c() ; + + if ( tri->collinear_edge_weight() == tri->other_collinear_edge_weight() ) { FT px, py ; line_project_pointC2(l0->a(),l0->b(),l0->c(),q->x(),q->y(), px,py); CGAL_STSKEL_TRAITS_TRACE("Degenerate, equal weights " << tri->collinear_edge_weight() ) ; CGAL_STSKEL_TRAITS_TRACE("Seed point: " << p2str(*q) << ". Projected seed point: (" << n2str(px) << "," << n2str(py) << ")" ) ; - const FT& l0a = l0->a() ; - const FT& l0b = l0->b() ; - const FT& l0c = l0->c() ; - const FT& l2a = l2->a() ; - const FT& l2b = l2->b() ; - const FT& l2c = l2->c() ; - // See details in compute_degenerate_offset_lines_isec_timeC2() - const FT sq_w0 = square(l0a) + square(l0b); - - // Note that "* sq_w0" is removed from the numerator expression. - // - // This is because the speed is inverted while representing the front - // progression using the orthogonal normalized vector [l0a, l0b]: P = Q + t/w * V with V normalized. - // However, here l0a & l0b are not normalized but *weighted* coeff, so we need to divide by w0². - // Hence we can just avoid multiplying by w0² in the numerator in the first place. - FT num, den ; + FT num_x, num_y, den ; if ( ! CGAL_NTS is_zero(l0->b()) ) // Non-vertical { - num = ((l2a*l0b - l0a*l2b) * px - l2b*l0c + l0b*l2c) /* * sq_w0 */ ; - den = l0a*l0a*l2b - l2b*sq_w0 + l0b*sq_w0 - l0a*l2a*l0b ; + num_x = -(l0b*px*w2 - (l0a*l2b*l0c - l0a*l0b*l2c + l2b*px)*w0) ; + num_y = (l0a*px*w2 - (l0a*l2a*l0c - l0a*l0a*l2c + l2a*px + l2c)*w0 + l0c*w2) ; + den = ((l0a*l2a*l0b - l0a*l0a*l2b + l2b)*w0 - l0b*w2) ; } else { - num = ((l2a*l0b - l0a*l2b) * py - l0a*l2c + l2a*l0c) /* * sq_w0 */ ; - den = l0a*l0b*l2b - l0b*l0b*l2a + l2a*sq_w0 - l0a*sq_w0; + num_x = -(l0b*py*w2 - (l0b*l2b*l0c - l0b*l0b*l2c + l2b*py + l2c)*w0 + l0c*w2) ; + num_y = (l0a*py*w2 - (l2a*l0b*l0c - l0a*l0b*l2c + l2a*py)*w0) ; + den = ((l2a*l0b*l0b - l0a*l0b*l2b - l2a)*w0 + l0a*w2) ; } - if ( ! CGAL_NTS certified_is_zero(den) && CGAL_NTS is_finite(den) && CGAL_NTS is_finite(num) ) + if ( ! CGAL_NTS certified_is_zero(den) && CGAL_NTS is_finite(den) && + CGAL_NTS is_finite(num_x) && CGAL_NTS is_finite(num_y) ) { - x = px + l0a * num / den ; - y = py + l0b * num / den ; - + x = num_x / den ; + y = num_y / den ; ok = CGAL_NTS is_finite(x) && CGAL_NTS is_finite(y) ; } } else { + CGAL_assertion(false); // this is COLLINEAR, not parallel so we can never be here + CGAL_STSKEL_TRAITS_TRACE("Degenerate, different weights " << n2str(tri->collinear_edge_weight()) << " and " << n2str(tri->other_collinear_edge_weight())); - const FT& l0a = l0->a() ; const FT& l0b = l0->b() ; const FT& l0c = l0->c() ; - const FT& l2a = l2->a() ; const FT& l2b = l2->b() ; const FT& l2c = l2->c() ; - // The line parallel to l0 (and l1) passing through q is: l0a*x + l0b*y + lambda = 0, with const FT lambda = -l0a*q->x() - l0b*q->y(); // The bisector between l0 (l1) and l2 is: - // l0a*x + l0b*y + l0c - t = 0 - // l2a*x + l2b*y + l2c - t = 0 + // l0a*x + l0b*y + l0c - w0*t = 0 + // l2a*x + l2b*y + l2c - w2*t = 0 // The intersection point is thus: - // l0a*x + l0b*y + l0c - t = 0 - // l2a*x + l2b*y + l2c - t = 0 + // l0a*x + l0b*y + l0c - w0*t = 0 + // l2a*x + l2b*y + l2c - w2*t = 0 // l0a*x + l0b*y + lambda = 0 // const FT t = l0c - lambda ; // (3) - (1) @@ -1038,6 +1056,9 @@ construct_offset_lines_isecC2 ( Trisegment_2_ptr< Trisegment_2collinearity() == TRISEGMENT_COLLINEARITY_NONE ? construct_normal_offset_lines_isecC2 (tri, aCaches) : construct_degenerate_offset_lines_isecC2(tri, aCaches); + + CGAL_STSKEL_TRAITS_TRACE("isec = " << (rRes ? p2str(*rRes) : "none")); + aCaches.mPoint_cache.Set(tri->id(), rRes) ; return rRes ; diff --git a/Straight_skeleton_2/include/CGAL/create_offset_polygons_2.h b/Straight_skeleton_2/include/CGAL/create_offset_polygons_2.h index 2a2a481dab07..8f1e868f16d3 100644 --- a/Straight_skeleton_2/include/CGAL/create_offset_polygons_2.h +++ b/Straight_skeleton_2/include/CGAL/create_offset_polygons_2.h @@ -141,7 +141,7 @@ template std::vector< boost::shared_ptr > create_offset_polygons_2 ( FT const& aOffset, Skeleton const& aSs, K const& , Tag_false ) { - static_assert(!(std::is_same::value)); + static_assert(!(std::is_same::value), ""); typedef boost::shared_ptr OutPolygonPtr ; typedef std::vector OutPolygonPtrVector ; @@ -167,7 +167,7 @@ template std::vector< boost::shared_ptr > create_offset_polygons_2 ( FT const& aOffset, Skeleton const& aSs, K const& /*k*/, Tag_true ) { - static_assert(!(std::is_same::value)); + static_assert(!(std::is_same::value), ""); typedef boost::shared_ptr OutPolygonPtr ; typedef std::vector OutPolygonPtrVector ; diff --git a/Straight_skeleton_2/include/CGAL/create_offset_polygons_from_polygon_with_holes_2.h b/Straight_skeleton_2/include/CGAL/create_offset_polygons_from_polygon_with_holes_2.h index 2037a42fa80a..d7671a35fbac 100644 --- a/Straight_skeleton_2/include/CGAL/create_offset_polygons_from_polygon_with_holes_2.h +++ b/Straight_skeleton_2/include/CGAL/create_offset_polygons_from_polygon_with_holes_2.h @@ -101,6 +101,8 @@ create_exterior_skeleton_and_offset_polygons_with_holes_2(const FT& aOffset, std::vector > raw_output = create_exterior_skeleton_and_offset_polygons_2(aOffset, aPoly, ofk, ssk); + CGAL_postcondition(raw_output.size() >= 2); + // filter offset of the outer frame std::swap(raw_output[0], raw_output.back()); raw_output.pop_back(); diff --git a/Straight_skeleton_2/include/CGAL/create_weighted_offset_polygons_2.h b/Straight_skeleton_2/include/CGAL/create_weighted_offset_polygons_2.h index 96783aaa2443..e7cecd95a377 100644 --- a/Straight_skeleton_2/include/CGAL/create_weighted_offset_polygons_2.h +++ b/Straight_skeleton_2/include/CGAL/create_weighted_offset_polygons_2.h @@ -127,7 +127,7 @@ create_partial_exterior_weighted_straight_skeleton_2(const FT& aMaxOffset, typedef typename Kernel_traits::Kernel IK; typedef typename IK::FT IFT; - static_assert((std::is_same::value_type, IFT>::value)); + static_assert((std::is_same::value_type, IFT>::value), ""); boost::shared_ptr > rSkeleton; @@ -170,8 +170,8 @@ create_partial_exterior_weighted_straight_skeleton_2(const FT& aMaxOffset, std::vector holes ; holes.push_back(lPoly) ; - // put a weight large enough such that frame edges are not relevant - const IFT frame_weight = FT(10) * *(std::max_element(aWeightsBegin, aWeightsEnd)); + // weight 0 such that the frame does not play a role + const IFT frame_weight = 0; CGAL_STSKEL_BUILDER_TRACE(4, "Frame weight = " << frame_weight); std::vector lFrameWeights(4, frame_weight); diff --git a/Straight_skeleton_2/include/CGAL/create_weighted_offset_polygons_from_polygon_with_holes_2.h b/Straight_skeleton_2/include/CGAL/create_weighted_offset_polygons_from_polygon_with_holes_2.h index ae1c57f2ab53..6ce84f74f9be 100644 --- a/Straight_skeleton_2/include/CGAL/create_weighted_offset_polygons_from_polygon_with_holes_2.h +++ b/Straight_skeleton_2/include/CGAL/create_weighted_offset_polygons_from_polygon_with_holes_2.h @@ -106,6 +106,8 @@ create_exterior_weighted_skeleton_and_offset_polygons_with_holes_2(const FT& aOf std::vector > raw_output = create_exterior_weighted_skeleton_and_offset_polygons_2(aOffset, aPoly, aWeights, ofk, ssk); + CGAL_postcondition(raw_output.size() >= 2); + // filter offset of the outer frame std::swap(raw_output[0], raw_output.back()); raw_output.pop_back(); diff --git a/Straight_skeleton_2/include/CGAL/create_weighted_straight_skeleton_2.h b/Straight_skeleton_2/include/CGAL/create_weighted_straight_skeleton_2.h index 5dfe7354467e..0328717b139d 100644 --- a/Straight_skeleton_2/include/CGAL/create_weighted_straight_skeleton_2.h +++ b/Straight_skeleton_2/include/CGAL/create_weighted_straight_skeleton_2.h @@ -148,7 +148,7 @@ create_exterior_weighted_straight_skeleton_2(const FT& max_offset, using IK = typename Kernel_traits::Kernel; using IFT = typename IK::FT; - static_assert((std::is_same::value_type, IFT>::value)); + static_assert((std::is_same::value_type, IFT>::value), ""); boost::shared_ptr > skeleton; @@ -190,8 +190,8 @@ create_exterior_weighted_straight_skeleton_2(const FT& max_offset, std::vector holes; holes.push_back(poly); - // put a weight large enough such that frame edges are not relevant - const IFT frame_weight = IFT(10) * *(std::max_element(weights_begin, weights_end)); + // weight 0 such that the frame does not play a role + const IFT frame_weight = 0; CGAL_STSKEL_BUILDER_TRACE(4, "Frame weight = " << frame_weight); std::vector lFrameWeights(4, frame_weight); diff --git a/Straight_skeleton_2/include/CGAL/draw_straight_skeleton_2.h b/Straight_skeleton_2/include/CGAL/draw_straight_skeleton_2.h index 3238ea1e0b1f..8beba39ef0fb 100644 --- a/Straight_skeleton_2/include/CGAL/draw_straight_skeleton_2.h +++ b/Straight_skeleton_2/include/CGAL/draw_straight_skeleton_2.h @@ -13,13 +13,17 @@ #define CGAL_DRAW_SS2_H #include + #include #ifdef CGAL_USE_BASIC_VIEWER #include + +#include #include +#include #include namespace CGAL { @@ -40,19 +44,22 @@ struct DefaultColorFunctorSS2 template class SimpleStraightSkeleton2ViewerQt : public Basic_viewer_qt { - typedef Basic_viewer_qt Base; + typedef Basic_viewer_qt Base; + typedef typename SS2::Vertex_const_handle Vertex_const_handle; typedef typename SS2::Halfedge_const_handle Halfedge_const_handle; + typedef typename SS2::Face_const_handle Face_const_handle; - // typedef typename SS2::Point Point; + typedef typename SS2::Traits::FT FT; + typedef typename SS2::Traits::Point_2 Point; public: /// Construct the viewer. /// @param ass2 the ss2 to view /// @param title the title of the window SimpleStraightSkeleton2ViewerQt(QWidget* parent, const SS2& ass2, - const char* title="Basic SS2 Viewer", - const ColorFunctor& fcolor=ColorFunctor()) : + const char* title="Basic SS2 Viewer", + const ColorFunctor& fcolor = ColorFunctor()) : // First draw: vertices; edges, faces; multi-color; no inverse normal Base(parent, title, true, true, true, false, false), ss2(ass2), @@ -62,32 +69,77 @@ class SimpleStraightSkeleton2ViewerQt : public Basic_viewer_qt } protected: - /* - void compute_face(Facet_const_handle fh) + void compute_face(Halfedge_const_handle eh) { - CGAL::IO::Color c=m_fcolor.run(ss2, fh); - face_begin(c); + // Skip faces with infinite halfedges + Halfedge_const_handle h = eh; + do { + if(h->vertex()->has_infinite_time()) { + return; + } + h = h->next(); + } while(h != eh); - add_point_in_face(fh->vertex(0)->point()); - add_point_in_face(fh->vertex(1)->point()); - add_point_in_face(fh->vertex(2)->point()); + CGAL::IO::Color c(150,170,170); + face_begin(c); + h = eh; + do { + add_point_in_face(h->vertex()->point()); + h = h->next(); + } while(h != eh); face_end(); } - */ + void compute_edge(Halfedge_const_handle eh) { + const bool src_inf = eh->opposite()->vertex()->has_infinite_time(); + const bool dst_inf = eh->vertex()->has_infinite_time(); + + // if there is an inf vertex, put it at the destination + const Point& src_p = src_inf ? eh->vertex()->point() : eh->opposite()->vertex()->point(); + Point dst_p = src_inf ? eh->opposite()->vertex()->point() : eh->vertex()->point(); + + if(src_inf || dst_inf) + dst_p = src_p + 0.1 * (dst_p - src_p); + if(eh->is_bisector()) - add_segment(eh->opposite()->vertex()->point(), eh->vertex()->point(), CGAL::IO::red()); + { + if (src_inf || dst_inf) + add_segment(src_p, dst_p, CGAL::IO::orange()); + else + add_segment(src_p, dst_p, CGAL::IO::red()); + } else - add_segment(eh->opposite()->vertex()->point(), eh->vertex()->point(), CGAL::IO::black()); + { + CGAL_assertion(!src_inf && !dst_inf); + if (eh->weight() == 0) + add_segment(src_p, dst_p, CGAL::IO::violet()); + else + add_segment(src_p, dst_p, CGAL::IO::black()); + } } void print_halfedge_labels(Halfedge_const_handle h) { + // Calculate length of the halfedge + double edge_length = std::sqrt(CGAL::squared_distance(h->opposite()->vertex()->point(), + h->vertex()->point())); + + double radius = 0.01 * edge_length; + + // Generate random point on circle centered at the midpoint + typedef CGAL::Random_points_on_circle_2 > Random_points; + Point midpoint = CGAL::barycenter(h->opposite()->vertex()->point(), 0.75, + h->vertex()->point(), 0.25); + + CGAL::Random rnd(h->id()); + Random_points random_point(radius, rnd); + Point label_pos = midpoint + (*random_point - CGAL::ORIGIN); + std::stringstream label; label << "H" << h->id() << " (V" << h->vertex()->id() << ") "; label << "H" << h->opposite()->id() << " (V" << h->opposite()->vertex()->id() << ") "; - add_text(CGAL::midpoint(h->opposite()->vertex()->point(), h->vertex()->point()), label.str()); + add_text(label_pos, label.str()); } void compute_vertex(Vertex_const_handle vh) @@ -103,9 +155,31 @@ class SimpleStraightSkeleton2ViewerQt : public Basic_viewer_qt } void print_vertex_label(Vertex_const_handle vh) { + // Calculate minimum incident edge length + double min_length = std::numeric_limits::max(); + auto h = vh->halfedge(); + do { + if(!h->vertex()->has_infinite_time() && !h->opposite()->vertex()->has_infinite_time()) { + double len = CGAL::squared_distance(h->vertex()->point(), + h->opposite()->vertex()->point()); + min_length = std::min(min_length, std::sqrt(len)); + } + h = h->next(); + } while(h != vh->halfedge()); + + double radius = 0.001 * min_length; + + // Generate random point on circle + typedef CGAL::Random_points_on_circle_2 > Random_points; + + CGAL::Random rnd(vh->id()); + Random_points random_point(radius, rnd); + Point label_pos = vh->point() + (*random_point - CGAL::ORIGIN); + std::stringstream label; label << "V" << vh->id() << std::ends; - add_text(vh->point(), label.str()); + + add_text(label_pos, label.str()); } void compute_elements() @@ -119,7 +193,11 @@ class SimpleStraightSkeleton2ViewerQt : public Basic_viewer_qt compute_edge(it); print_halfedge_labels(it); } + + if(!it->is_bisector() && !it->is_border()) + compute_face(it); } + for (typename SS2::Vertex_const_iterator it=ss2.vertices_begin(); it!=ss2.vertices_end(); ++it) { compute_vertex(it); diff --git a/Straight_skeleton_2/include/CGAL/predicates/Straight_skeleton_pred_ftC2.h b/Straight_skeleton_2/include/CGAL/predicates/Straight_skeleton_pred_ftC2.h index bdd371ea3704..96b6055b6973 100644 --- a/Straight_skeleton_2/include/CGAL/predicates/Straight_skeleton_pred_ftC2.h +++ b/Straight_skeleton_2/include/CGAL/predicates/Straight_skeleton_pred_ftC2.h @@ -205,7 +205,7 @@ compare_offset_lines_isec_timesC2 ( Trisegment_2_ptr< Trisegment_2 Quotient ; typedef boost::optional Optional_rational ; - CGAL_STSKEL_TRAITS_TRACE("compare_offset_lines_isec_timesC2(\n" << m << "\n" << n << "\n) [" << typeid(FT).name() << "]" ); + CGAL_STSKEL_TRAITS_TRACE("compare_offset_lines_isec_timesC2(" << m->id() << " " << n->id() << ") [" << typeid(FT).name() << "]" ); Uncertain rResult = Uncertain::indeterminate(); @@ -369,8 +369,8 @@ oriented_side_of_event_point_wrt_bisectorC2 ( Trisegment_2_ptr< Trisegment_2 0) ? "POSITIVE" : "NEGATIVE") << " side of reflex bisector" ) ; - } - else + CGAL_STSKEL_TRAITS_TRACE("Point is exactly at bisector"); + + rResult = ON_ORIENTED_BOUNDARY; + } + else + { + // below is independent of weights because weights have the same sign, which only tilts + // the bisector within a quadrant defined by the intersection of the supporting + // lines of the edges + Uncertain smaller = CGAL_NTS certified_is_smaller( validate(l0.a()*l1.b()), validate(l1.a()*l0.b()) ) ; + if ( is_certain(smaller) ) { - rResult = (lCmpResult == LARGER) ? ON_NEGATIVE_SIDE : ON_POSITIVE_SIDE ; - CGAL_STSKEL_TRAITS_TRACE("Event point is on " << ((rResult > 0) ? "POSITIVE" : "NEGATIVE") << " side of convex bisector" ) ; + // Reflex bisector? + if ( smaller ) + { + rResult = (lCmpResult == SMALLER) ? ON_NEGATIVE_SIDE : ON_POSITIVE_SIDE ; + CGAL_STSKEL_TRAITS_TRACE("Event point is on " << ((rResult > 0) ? "POSITIVE" : "NEGATIVE") << " side of reflex bisector" ) ; + } + else + { + rResult = (lCmpResult == LARGER) ? ON_NEGATIVE_SIDE : ON_POSITIVE_SIDE ; + CGAL_STSKEL_TRAITS_TRACE("Event point is on " << ((rResult > 0) ? "POSITIVE" : "NEGATIVE") << " side of convex bisector" ) ; + } } } } @@ -462,6 +494,8 @@ oriented_side_of_event_point_wrt_bisectorC2 ( Trisegment_2_ptr< Trisegment_2 -class Test_polygon_2 : public CGAL::Polygon_2 { +struct Test_polygon_2 : public CGAL::Polygon_2 { typedef CGAL::Polygon_2 Base; + Test_polygon_2() { } Test_polygon_2(const Base&); public: using Base::Base; }; template -class Test_polygon_with_holes_2 : public CGAL::Polygon_with_holes_2 { +struct Test_polygon_with_holes_2 : public CGAL::Polygon_with_holes_2 { typedef CGAL::Polygon_with_holes_2 Base; Test_polygon_with_holes_2(const Base&); public: diff --git a/Straight_skeleton_2/test/Straight_skeleton_2/test_sls_weighted_offset.cpp b/Straight_skeleton_2/test/Straight_skeleton_2/test_sls_weighted_offset.cpp index 40fd4aede6f9..02fb4f56dbaa 100644 --- a/Straight_skeleton_2/test/Straight_skeleton_2/test_sls_weighted_offset.cpp +++ b/Straight_skeleton_2/test/Straight_skeleton_2/test_sls_weighted_offset.cpp @@ -28,19 +28,21 @@ typedef CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt EPECK_w_sqr namespace CGAL { template -class Test_polygon_2 : public CGAL::Polygon_2 { - typedef CGAL::Polygon_2 Base; - Test_polygon_2(const Base&); +struct Test_polygon_2 : public CGAL::Polygon_2 { + typedef CGAL::Polygon_2 Base; + Test_polygon_2() { } + Test_polygon_2(const Base&); public: - using Base::Base; + using Base::Base; }; template -class Test_polygon_with_holes_2 : public CGAL::Polygon_with_holes_2 { - typedef CGAL::Polygon_with_holes_2 Base; - Test_polygon_with_holes_2(const Base&); +struct Test_polygon_with_holes_2 : public CGAL::Polygon_with_holes_2 { + typedef CGAL::Polygon_with_holes_2 Base; + Test_polygon_with_holes_2(); + Test_polygon_with_holes_2(const Base&); public: - using Base::Base; + using Base::Base; }; } // namespace CGAL diff --git a/Straight_skeleton_2/test/Straight_skeleton_2/test_sls_weighted_polygons_with_holes.cpp b/Straight_skeleton_2/test/Straight_skeleton_2/test_sls_weighted_polygons_with_holes.cpp index ffa31c47d5bc..9b3c66a73ea3 100644 --- a/Straight_skeleton_2/test/Straight_skeleton_2/test_sls_weighted_polygons_with_holes.cpp +++ b/Straight_skeleton_2/test/Straight_skeleton_2/test_sls_weighted_polygons_with_holes.cpp @@ -32,6 +32,8 @@ typedef CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt EPECK_w_sqr template void test_kernel(const int hole_n, const int hole_nv, CGAL::Random& rnd) { + std::cout << "\n ==== Test with Kernel: " << typeid(K).name() << " ====" << std::endl; + using FT = typename K::FT; using Point_2 = typename K::Point_2; using Point_3 = typename K::Point_3; @@ -124,8 +126,17 @@ void test_kernel(const int hole_n, const int hole_nv, CGAL::Random& rnd) // CGAL::draw(*ss_ptr); + // randomly switch the weight signs + if(rnd.get_int(0, 2)) + { + std::cout << "using negative weights" << std::endl; + for(std::vector& ws : weights) + for(FT& w : ws) + w = -w; + } + Mesh sm; - bool success = extrude_skeleton(pwh, sm, CGAL::parameters::weights(weights)); + bool success = extrude_skeleton(pwh, sm, CGAL::parameters::weights(weights).verbose(true).maximum_height(hole_n)); assert(success); if(!success) { @@ -135,7 +146,7 @@ void test_kernel(const int hole_n, const int hole_nv, CGAL::Random& rnd) std::cout << num_vertices(sm) << " vertices and " << num_faces(sm) << " faces" << std::endl; -// CGAL::draw(sm); + CGAL::draw(sm); } int main(int argc, char** argv) diff --git a/Straight_skeleton_extrusion_2/include/CGAL/extrude_skeleton.h b/Straight_skeleton_extrusion_2/include/CGAL/extrude_skeleton.h index 38e178d72226..bd2ca370e84a 100644 --- a/Straight_skeleton_extrusion_2/include/CGAL/extrude_skeleton.h +++ b/Straight_skeleton_extrusion_2/include/CGAL/extrude_skeleton.h @@ -20,6 +20,7 @@ #include #include #include + #include #endif #include @@ -85,10 +86,19 @@ inline constexpr FT default_extrusion_height() return (std::numeric_limits::max)(); } +template +typename GeomTraits::Point_3 lift(const typename GeomTraits::Point_2& p, + const typename GeomTraits::FT& h) +{ + return { p.x(), p.y(), h }; +} + +// #define CGAL_SS2_NUDGE_ZERO_WEIGHTS + // @todo Maybe this postprocessing is not really necessary? Do users really care if the point -// is not perfectly above the input contour edge (it generally cannot be anyway if the kernel is not exact except for some -// specific cases)? -#define CGAL_SLS_SNAP_TO_VERTICAL_SLABS +// is not perfectly above the input contour edge (it generally cannot be anyway if the kernel +// is not exact except for some specific cases)? +// #define CGAL_SLS_SNAP_TO_VERTICAL_SLABS #ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS // The purpose of this visitor is to snap back almost-vertical (see preprocessing_weights()) edges @@ -134,20 +144,109 @@ snap_point_to_contour_halfedge_plane(const typename GeomTraits::Point_2& op, } } +template +typename GeomTraits::Point_2 +hds_intersection(typename HDS::Halfedge_const_handle h_1, + typename HDS::Halfedge_const_handle h_2) +{ + using FT = typename GeomTraits::FT; + using Point_2 = typename GeomTraits::Point_2; + using Line_2 = typename GeomTraits::Line_2; + + CGAL_SS_i::Segment_2_with_ID lS1 (h_1->opposite()->vertex()->point(), + h_1->vertex()->point(), + h_1->id()); + CGAL_SS_i::Segment_2_with_ID lS2 (h_2->opposite()->vertex()->point(), + h_2->vertex()->point(), + h_2->id()); + + // @fixme use caches + boost::optional l1 = compute_normalized_line_coeffC2(lS1); + boost::optional l2 = compute_normalized_line_coeffC2(lS2); + CGAL_assertion(bool(l1) && bool(l2)); + + FT denom = l1->a()*l2->b() - l2->a()*l1->b(); + CGAL_assertion(!is_zero(denom)); + FT num_x = (l1->b()*l2->c() - l2->b()*l1->c()); + FT num_y = (l2->a()*l1->c() - l1->a()*l2->c()); + FT ix = num_x / denom; + FT iy = num_y / denom; + + return Point_2(ix, iy); +} + template void snap_skeleton_vertex(typename HDS::Halfedge_const_handle hds_h, - typename HDS::Halfedge_const_handle contour_h, + const typename GeomTraits::FT vertical_weight, std::map& snapped_positions) { - typename HDS::Vertex_const_handle hds_tv = hds_h->vertex(); + using HDS_Vertex_const_handle = typename HDS::Vertex_const_handle; + using HDS_Halfedge_const_handle = typename HDS::Halfedge_const_handle; + + using FT = typename GeomTraits::FT; + using Point_2 = typename GeomTraits::Point_2; + using Line_2 = typename GeomTraits::Line_2; + + HDS_Vertex_const_handle hds_tv = hds_h->vertex(); + if(hds_tv->is_contour()) + return; + + auto insert_res = snapped_positions.emplace(hds_tv->point(), Point_2()); + if(insert_res.second) // snapped position already computed + return; - // this re-applies snapping towards contour_h even if the point was already snapped towards another contour - auto insert_result = snapped_positions.emplace(hds_tv->point(), hds_tv->point()); - insert_result.first->second = snap_point_to_contour_halfedge_plane(insert_result.first->second, contour_h); + // gather incident contour (vertical) halfedges + std::vector vertical_chs; + std::vector non_vertical_chs; - // std::cout << "snap_skeleton_vertex(V" << hds_tv->id() << " pt: " << hds_h->vertex()->point() << ")" - // << " to " << insert_result.first->second << std::endl; + HDS_Halfedge_const_handle curr = hds_h; + do + { + HDS_Halfedge_const_handle ch = curr->defining_contour_edge(); + CGAL_assertion(ch->opposite()->is_border()); + + if(contour_h->weight() == vertical_weight) + vertical_chs.push_back(ch); + else + non_vertical_chs.push_back(ch); + + curr = curr->next()->opposite(); + } + while(curr != hds_h); + + CGAL_assertion(non_vertical_chs.size() + vertical_chs.size() == 3); + + if(vertical_chs.empty()) + return; + + if(vertical_chs.size() == 2) + { + HDS_Halfedge_const_handle non_vertical_ch = non_vertical_chs.front(); + + // The new position is at the intersection of the two vertical lines + const Point_2 sp = hds_intersection(vertical_chs[0], vertical_chs[1]); + + // The time (z coordinate) is the weighted distance to the non-vertical edge + CGAL_SS_i::Segment_2_with_ID seg(non_vertical_ch->opposite()->vertex()->point(), + non_vertical_ch->vertex()->point(), + non_vertical_ch->id()); + + boost::optional line = compute_weighted_line_coeffC2(seg); + CGAL_assertion(bool(line)); + + // Compute weighted distance + const FT a = line->a(), b = line->b(), c = line->c(); + const FT t = a*sp.x() + b*sp.y() + c; + + insert_res.first->second = lift(sp, t); + } + else + { + CGAL_assertion(vertical_chs.size() == 1); + + + } } template @@ -176,7 +275,9 @@ class Skeleton_offset_correspondence_builder_visitor using FT = typename GeomTraits::FT; using Point_2 = typename GeomTraits::Point_2; + using Line_2 = typename GeomTraits::Line_2; + using SS_Vertex_const_handle = typename StraightSkeleton_2::Vertex_const_handle; using SS_Halfedge_const_handle = typename StraightSkeleton_2::Halfedge_const_handle; using HDS = typename StraightSkeleton_2::Base; @@ -199,23 +300,22 @@ class Skeleton_offset_correspondence_builder_visitor { } public: - void on_offset_contour_started() const - { - // std::cout << "~~ new contour ~~" << std::endl; - } + void on_offset_contour_started() const { } // can't modify the position yet because we need arrange_polygons() to still work properly + // + // @fixme on paper one could create a polygon thin-enough w.r.t. the max weight value + // such thatthere is a skeleton vertex that wants to be snapped to two different sides... void on_offset_point(const Point_2& op, SS_Halfedge_const_handle hook) const { - CGAL_assertion(hook->is_bisector()); - -#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS - // @fixme on paper one could create a polygon thin-enough w.r.t. the max weight value such that - // there is a skeleton vertex that wants to be snapped to two different sides... - CGAL_assertion(m_snapped_positions.count(op) == 0); + CGAL_precondition(hook->is_bisector()); HDS_Halfedge_const_handle canonical_hook = (hook < hook->opposite()) ? hook : hook->opposite(); + m_offset_points[canonical_hook] = op; + +#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS + CGAL_precondition(m_snapped_positions.count(op) == 0); SS_Halfedge_const_handle contour_h1 = hook->defining_contour_edge(); CGAL_assertion(contour_h1->opposite()->is_border()); @@ -225,18 +325,30 @@ class Skeleton_offset_correspondence_builder_visitor const bool is_h1_vertical = (contour_h1->weight() == m_vertical_weight); const bool is_h2_vertical = (contour_h2->weight() == m_vertical_weight); - // this can happen when the offset is passing through vertices - m_offset_points[canonical_hook] = op; - - // if both are vertical, it's the common vertex (which has to exist) if(is_h1_vertical && is_h2_vertical) { - CGAL_assertion(contour_h1->vertex() == contour_h2->opposite()->vertex() || - contour_h2->vertex() == contour_h1->opposite()->vertex()); + std::cout << "BOTH VERTICAL" << std::endl; + std::cout << contour_h1->opposite()->vertex()->point() << " " << contour_h1->vertex()->point() << std::endl; + std::cout << contour_h2->opposite()->vertex()->point() << " " << contour_h2->vertex()->point() << std::endl; + if(contour_h1->vertex() == contour_h2->opposite()->vertex()) + { m_snapped_positions[op] = contour_h1->vertex()->point(); - else + } + else if(contour_h2->vertex() == contour_h1->opposite()->vertex()) + { m_snapped_positions[op] = contour_h2->vertex()->point(); + } + else + { + // this can happen for example with outward extrusion and two vertical edges + // that are non consecutive and with a convex region of non-vertical edges in between. + // Their extension will enclose the region. + // + // In this case, the point position will be at the intersection of the supporting lines + // of the two edges. + m_snapped_positions[op] = hds_intersection(contour_h1, contour_h2); + } } else if(is_h1_vertical) { @@ -321,9 +433,14 @@ class Extrusion_builder private: Geom_traits m_gt; + const bool m_verbose; public: - Extrusion_builder(const Geom_traits& gt) : m_gt(gt) { } + Extrusion_builder(const Geom_traits& gt, + const bool verbose) + : m_gt(gt), + m_verbose(verbose) + { } //////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////// @@ -338,7 +455,7 @@ class Extrusion_builder const bool invert_faces = false) { #ifdef CGAL_SLS_DEBUG_DRAW - CGAL::draw(p); + // CGAL::draw(p); #endif CDT cdt; @@ -351,7 +468,8 @@ class Extrusion_builder } catch(const typename CDT::Intersection_of_constraints_exception&) { - std::cerr << "Warning: Failed to triangulate horizontal face" << std::endl; + if(m_verbose) + std::cerr << "Warning: Failed to triangulate horizontal face" << std::endl; return; } @@ -406,7 +524,8 @@ class Extrusion_builder std::rotate(face_points.rbegin(), face_points.rbegin() + 1, face_points.rend()); CGAL_assertion(face_points[0][2] == 0 && face_points[1][2] == 0); - const Vector_3 n = CGAL::cross_product(face_points[1] - face_points[0], face_points[2] - face_points[0]); + const Vector_3 n = CGAL::cross_product(face_points[1] - face_points[0], + face_points[2] - face_points[0]); PK traits(n); PCDT pcdt(traits); @@ -416,7 +535,14 @@ class Extrusion_builder } catch(const typename PCDT::Intersection_of_constraints_exception&) { - std::cerr << "Warning: Failed to triangulate skeleton face" << std::endl; + if(m_verbose) + { + std::cerr << "Warning: Failed to triangulate skeleton face" << std::endl; + for(const auto& p : face_points) { + std::cerr << " " << p << std::endl; + } + } + return; } @@ -474,20 +600,13 @@ class Extrusion_builder continue; HDS_Halfedge_const_handle hds_h = hds_f->halfedge(), done = hds_h; -#ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS - HDS_Halfedge_const_handle contour_h = hds_h->defining_contour_edge(); - CGAL_assertion(hds_h == contour_h); - const bool is_vertical = (contour_h->weight() == vertical_weight); -#endif - do { HDS_Vertex_const_handle hds_tv = hds_h->vertex(); #ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS // this computes the snapped position but does not change the geometry of the skeleton - if(is_vertical && !hds_tv->is_contour()) - snap_skeleton_vertex(hds_h, contour_h, snapped_positions); + snap_skeleton_vertex(hds_h, vertical_weight, snapped_positions); #endif face_points.emplace_back(hds_tv->point().x(), hds_tv->point().y(), hds_tv->time()); @@ -498,7 +617,8 @@ class Extrusion_builder if(face_points.size() < 3) { - std::cerr << "Warning: sm_vs has size 1 or 2: offset crossing face at a single point?" << std::endl; + if(m_verbose) + std::cerr << "Warning: sm_vs has size 1 or 2: offset crossing a face at a vertex?" << std::endl; continue; } @@ -601,7 +721,8 @@ class Extrusion_builder if(face_points.size() < 3) { - std::cerr << "Warning: sm_vs has size 1 or 2: offset crossing face at a single point?" << std::endl; + if(m_verbose) + std::cerr << "Warning: sm_vs has size 1 or 2: offset crossing a face at a vertex?" << std::endl; continue; } @@ -622,6 +743,7 @@ class Extrusion_builder FaceRange& faces) { CGAL_precondition(!is_zero(height)); + CGAL_USE(vertical_weight); const FT abs_height = abs(height); @@ -669,17 +791,22 @@ class Extrusion_builder if(!ss_ptr) { - std::cerr << "Error: encountered an error during skeleton construction" << std::endl; + if(m_verbose) + std::cerr << "Error: encountered an error during skeleton construction" << std::endl; return false; } #ifdef CGAL_SLS_DEBUG_DRAW - // print_straight_skeleton(*ss_ptr); + Straight_skeletons_2::IO::print_straight_skeleton(*ss_ptr); CGAL::draw(*ss_ptr); #endif + std::cout << "built SS" << std::endl; + if(is_default_extrusion_height) { + std::cout << "lateral faces..." << std::endl; + #ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS construct_lateral_faces(*ss_ptr, points, faces, vertical_weight, snapped_positions); #else @@ -691,15 +818,19 @@ class Extrusion_builder #ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS Visitor visitor(*ss_ptr, offset_points, vertical_weight, snapped_positions); #else - Visitor visitor(*ss_ptr, vertical_weight, offset_points); + Visitor visitor(*ss_ptr, offset_points); #endif Offset_builder ob(*ss_ptr, Offset_builder_traits(), visitor); Offset_polygons raw_output; ob.construct_offset_contours(abs_height, std::back_inserter(raw_output)); + std::cout << "horizontal faces..." << std::endl; + Offset_polygons_with_holes output = CGAL::arrange_offset_polygons_2(raw_output); construct_horizontal_faces(output, height, points, faces); + std::cout << "lateral faces..." << std::endl; + #ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS construct_lateral_faces(*ss_ptr, ob, height, points, faces, offset_points, vertical_weight, snapped_positions); #else @@ -715,9 +846,12 @@ class Extrusion_builder } #ifdef CGAL_SLS_SNAP_TO_VERTICAL_SLABS + std::cout << "Apply snapping..." << std::endl; apply_snapping(points, snapped_positions); #endif + std::cout << "Done" << std::endl; + return true; } @@ -732,6 +866,7 @@ class Extrusion_builder { CGAL_precondition(!is_zero(height)); CGAL_precondition(height != default_extrusion_height()); // was checked before, this is just a reminder + CGAL_USE(vertical_weight); const FT abs_height = abs(height); @@ -764,12 +899,13 @@ class Extrusion_builder if(!ss_ptr) { - std::cerr << "Error: encountered an error during outer skeleton construction" << std::endl; + if(m_verbose) + std::cerr << "Error: encountered an error during outer skeleton construction" << std::endl; return false; } #ifdef CGAL_SLS_DEBUG_DRAW - // print_straight_skeleton(*ss_ptr); + Straight_skeletons_2::IO::print_straight_skeleton(*ss_ptr); CGAL::draw(*ss_ptr); #endif @@ -781,6 +917,8 @@ class Extrusion_builder Offset_builder ob(*ss_ptr, Offset_builder_traits(), visitor); ob.construct_offset_contours(abs_height, std::back_inserter(raw_output)); + CGAL_postcondition(raw_output.size() >= 2); + // Manually filter the offset of the outer frame std::swap(raw_output[0], raw_output.back()); raw_output.pop_back(); @@ -813,18 +951,19 @@ class Extrusion_builder CGAL_SS_i::vertices_begin(hole), CGAL_SS_i::vertices_end(hole), std::begin(no_holes), std::end(no_holes), - std::begin(speeds[hole_id]), std::end(speeds[hole_id]), + std::begin(speeds[1 + hole_id]), std::end(speeds[1 + hole_id]), std::begin(no_speeds), std::end(no_speeds), m_gt); if(!ss_ptr) { - std::cerr << "Error: encountered an error during skeleton construction" << std::endl; + if(m_verbose) + std::cerr << "Error: encountered an error during skeleton construction" << std::endl; return EXIT_FAILURE; } #ifdef CGAL_SLS_DEBUG_DRAW - // print_straight_skeleton(*ss_ptr); + Straight_skeletons_2::IO::print_straight_skeleton(*ss_ptr); CGAL::draw(*ss_ptr); #endif @@ -885,17 +1024,27 @@ void convert_angles(AngleRange& angles) { CGAL_precondition(0 < angle && angle < 180); - // @todo should this be an epsilon around 90°? As theta goes to 90°, tan(theta) goes to infinity - // and thus we could get numerical issues (overlfows) if the kernel is not exact - if(angle == 90) + if(angle == 90) { // distinguish to ensure zero is returned return 0; - else - return std::tan(CGAL::to_double(angle * CGAL_PI / 180)); + } else if (angle == 45) { + return 1; + } else if (angle == 135) { + return -1; + } else { + const double rad_angle = CGAL::to_double(angle * CGAL_PI / 180); + return std::cos(rad_angle) / std::sin(rad_angle); // cotan + } }; for(auto& contour_angles : angles) + { for(FT& angle : contour_angles) + { + std::cout << "convert " << angle << " to "; angle = angle_to_weight(angle); + std::cout << angle << std::endl; + } + } } // handle vertical angles (inf speed) @@ -916,6 +1065,8 @@ preprocess_weights(WeightRange& weights) { for(FT& w : contour_weights) { + std::cout << "w = " << w << std::endl; + // '0' means a vertical slab, aka 90° angle (see preprocess_angles()) if(w == 0) continue; @@ -928,12 +1079,12 @@ preprocess_weights(WeightRange& weights) } else if(slope == Slope::INWARD && w < 0) { - std::cerr << "Error: mixing positive and negative weights is not yet supported" << std::endl; + std::cerr << "Error: mixing positive and negative weights is not supported" << std::endl; return {Slope::UNKNOWN, false, FT(-1)}; } else if(slope == Slope::OUTWARD && w > 0) { - std::cerr << "Error: mixing positive and negative weights is not yet supported" << std::endl; + std::cerr << "Error: mixing positive and negative weights is not supported" << std::endl; return {Slope::UNKNOWN, false, FT(-1)}; } @@ -945,17 +1096,16 @@ preprocess_weights(WeightRange& weights) } if(slope == Slope::UNKNOWN) - { - std::cerr << "Warning: all edges vertical?" << std::endl; slope = Slope::VERTICAL; - } +#ifdef CGAL_SS2_NUDGE_ZERO_WEIGHTS // Take a weight which is a large % of the max value to ensure there's no ambiguity // // Since the max value might not be very close to 90°, take the max between of the large-% weight // and the weight corresponding to an angle of 89.9999999° const FT weight_of_89d9999999 = 572957787.3425436; // tan(89.9999999°) - const FT scaled_max = (std::max)(weight_of_89d9999999, 1e3 * max_value); + FT scaled_max = (std::max)(weight_of_89d9999999, 1e3 * max_value); + scaled_max = 1000; for(auto& contour_weights : weights) { @@ -966,7 +1116,11 @@ preprocess_weights(WeightRange& weights) } } + std::cout << "scaled_max = " << scaled_max << std::endl; return {slope, true, scaled_max}; +#else + return {slope, true, 0 /*vertical weight*/}; +#endif } template & contour_weights : weights) { + for(const FT& w : contour_weights) { + std::cout << " " << w; + } + std::cout << std::endl; + } + // End of preprocessing, start the actual skeleton computation if(slope != Slope::INWARD && height == default_extrusion_height()) @@ -1037,7 +1199,7 @@ bool extrude_skeleton(const PolygonWithHoles& pwh, points.reserve(2 * pwh.outer_boundary().size()); // just a reasonnable guess faces.reserve(2 * pwh.outer_boundary().size() + 2*pwh.number_of_holes()); - Extrusion_builder builder(gt); + Extrusion_builder builder(gt, verbose); bool res; if(slope != Slope::OUTWARD) // INWARD or VERTICAL res = builder.inward_construction(pwh, weights, vertical_weight, height, points, faces); @@ -1063,6 +1225,8 @@ bool extrude_skeleton(const PolygonWithHoles& pwh, CGAL_warning(is_valid_polygon_mesh(out) && is_closed(out)); + std::cout << "Finished" << std::endl; + return true; } diff --git a/Straight_skeleton_extrusion_2/test/Straight_skeleton_extrusion_2/CMakeLists.txt b/Straight_skeleton_extrusion_2/test/Straight_skeleton_extrusion_2/CMakeLists.txt index f29dcb4e9296..fbdf8e26116a 100644 --- a/Straight_skeleton_extrusion_2/test/Straight_skeleton_extrusion_2/CMakeLists.txt +++ b/Straight_skeleton_extrusion_2/test/Straight_skeleton_extrusion_2/CMakeLists.txt @@ -5,18 +5,14 @@ cmake_minimum_required(VERSION 3.1...3.23) project(Straight_skeleton_extrusion_2_Tests) find_package(CGAL REQUIRED COMPONENTS Qt5 Core) +find_package(LEDA QUIET) include_directories(BEFORE "include") -# create a target per cppfile -file( - GLOB cppfiles - RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} - ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp) -foreach(cppfile ${cppfiles}) - create_single_source_cgal_program("${cppfile}") -endforeach() +if (CGAL_Core_FOUND OR LEDA_FOUND) + create_single_source_cgal_program("test_sls_extrude.cpp") -if(CGAL_Qt5_FOUND) - target_link_libraries(test_sls_extrude PUBLIC CGAL::CGAL_Basic_viewer) + if(CGAL_Qt5_FOUND) + target_link_libraries(test_sls_extrude PUBLIC CGAL::CGAL_Basic_viewer) + endif() endif() diff --git a/Straight_skeleton_extrusion_2/test/Straight_skeleton_extrusion_2/test_sls_extrude.cpp b/Straight_skeleton_extrusion_2/test/Straight_skeleton_extrusion_2/test_sls_extrude.cpp index eaa971ce880c..4f1f68d8f1b1 100644 --- a/Straight_skeleton_extrusion_2/test/Straight_skeleton_extrusion_2/test_sls_extrude.cpp +++ b/Straight_skeleton_extrusion_2/test/Straight_skeleton_extrusion_2/test_sls_extrude.cpp @@ -1,9 +1,41 @@ +#define CGAL_SLS_DEBUG_DRAW + +#include +#include +#include + +// #define CGAL_SLS_PRINT_QUEUE_BEFORE_EACH_POP +// #define CGAL_STRAIGHT_SKELETON_ENABLE_TRACE 100 +// #define CGAL_STRAIGHT_SKELETON_TRAITS_ENABLE_TRACE 10000000 +// #define CGAL_STRAIGHT_SKELETON_VALIDITY_ENABLE_TRACE +// #define CGAL_POLYGON_OFFSET_ENABLE_TRACE 10000000 + +void Straight_skeleton_external_trace(std::string m) +{ + std::cout << std::setprecision(17) << m << std::endl << std::endl ; +} + +void Straight_skeleton_traits_external_trace(std::string m) +{ + std::cout << std::setprecision(17) << m << std::endl << std::endl ; +} + +#include + #include #include #include #include -#include +#include +#include +#include +#include + +#include +#include +#include +#include #include #include @@ -188,10 +220,246 @@ std::string root_name(const std::string& full_filename) } template -bool test(const char* poly_filename, - const char* angles_filename, - const typename K::FT height, - const typename K::FT expected_volume) +void test_fixed_with_vertical_combinations() +{ + using FT = typename K::FT; + using Point_2 = typename K::Point_2; + using Point_3 = typename K::Point_3; + + using Polygon_2 = CGAL::Polygon_2; + using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; + + using Mesh = CGAL::Surface_mesh; + + auto test_combination = [&](const auto& pwh, + const std::vector >& angles, + const FT max_height) + { + Mesh sm; + bool success = extrude_skeleton(pwh, sm, CGAL::parameters::angles(angles).maximum_height(max_height).verbose(true)); + assert(success); + if(!success) + { + std::cerr << "Error: failed to extrude skeleton" << std::endl; + std::exit(1); + } + + CGAL::draw(sm); + CGAL::IO::write_OFF("last_extrusion.off", sm, CGAL::parameters::stream_precision(17)); + + assert(is_closed(sm)); + }; + + auto test_all_vertical_combinations = [&](const Polygon_2& poly) + { + const std::size_t n = poly.size(); + const FT base_angle = 45; + const FT vertical_angle = 90; + const FT height = 2.0; // @todo multiple heights + + // Generate all possible combinations of vertical/non-vertical edges + for(std::size_t mask = 0; mask < (1u << n); ++mask) + { + std::vector> angles; + std::vector contour_angles(n); + + // Set angles based on the current mask + for(std::size_t i = 0; i < n; ++i) + contour_angles[i] = (mask & (1u << i)) ? vertical_angle : base_angle; + + angles.push_back(contour_angles); + + std::cout << "Testing combination:"; + for(FT angle : contour_angles) + std::cout << " " << angle; + std::cout << std::endl; + + test_combination(poly, angles, height); + } + }; + + std::cout << "\nTesting 'square' combinations:" << std::endl; + std::vector square_vertices = { { 0, 0 }, { 1, 0 }, { 1, 1 }, { 0, 1 } }; + Polygon_2 square { square_vertices.begin(), square_vertices.end() }; + // test_all_vertical_combinations(square); // @tmp + + std::cout << "\nTesting 'touching squares' combinations:" << std::endl; + std::vector touching_squares_vertices = { + { 0, 1 }, { 2, 1 }, { 2, 0 }, { 4, 0 }, { 4, 2 }, { 2, 2 }, { 2, 3 }, { 0, 3 } + }; + Polygon_2 touching_squares { touching_squares_vertices.begin(), touching_squares_vertices.end() }; + test_all_vertical_combinations(touching_squares); +} + +template +void test_fixed() +{ + using FT = typename K::FT; + using Point_2 = typename K::Point_2; + using Point_3 = typename K::Point_3; + + using Polygon_2 = CGAL::Polygon_2; + using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; + + using Mesh = CGAL::Surface_mesh; + + auto test_fixed_polygon = [&](const auto& pwh, + const std::vector >& angles, + const FT max_height, + const FT expected_vol) + { + Mesh sm; + bool success = extrude_skeleton(pwh, sm, CGAL::parameters::angles(angles).verbose(true).maximum_height(max_height)); + assert(success); + if(!success) + { + std::cerr << "Error: failed to extrude skeleton" << std::endl; + std::exit(1); + } + + CGAL::draw(sm); + CGAL::IO::write_OFF("last_extrusion.off", sm, CGAL::parameters::stream_precision(17)); + + FT actual_vol = CGAL::Polygon_mesh_processing::volume(sm); + std::cout << "Expected volume: " << expected_vol << ", actual is: " << actual_vol << std::endl; + if (CGAL::abs(expected_vol - actual_vol) > 1e-6) + { + assert(false); + std::exit(1); + } + }; + + // unit square + std::vector square_vertices = { { 0, 0 }, { 1, 0 }, { 1, 1 }, { 0, 1 } }; + Polygon_2 square { square_vertices.begin(), square_vertices.end() }; + + // unit square, 45° inward extrusion + const std::vector > square_angles_45 = { { 45, 45, 45, 45 } }; + test_fixed_polygon(square, square_angles_45, 1, FT(1)/6); + test_fixed_polygon(square, square_angles_45, -1, FT(1)/6); + + // unit square, 80° inward extrusion + const std::vector > square_angles_80 = { { 80, 80, 80, 80 } }; + test_fixed_polygon(square, square_angles_80, 1000000, 0.9452136366); + test_fixed_polygon(square, square_angles_80, -1000000, 0.9452136366); + test_fixed_polygon(square, square_angles_80, 1, 0.6888009778); + test_fixed_polygon(square, square_angles_80, -1, 0.6888009778); + + // unit square, 10° outward extrusion + // the volume of a truncated pyramid is: 1/3 × h × (a^2 + b^2 + ab) + // and the extruded square has side length: b = a + 2 * (h * tan(alpha - 90°)) + const std::vector > square_angles_m10 = { { 170, 170, 170, 170 } }; + test_fixed_polygon(square, square_angles_m10, 1, 55.2271469426); + test_fixed_polygon(square, square_angles_m10, -1, 55.2271469426); + test_fixed_polygon(square, square_angles_m10, 10, 44028.8396672); + test_fixed_polygon(square, square_angles_m10, -10, 44028.8396672); + + const std::vector > square_angles_m45 = { { 135, 135, 135, 135 } }; + test_fixed_polygon(square, square_angles_m45, 1, FT(13)/3); + test_fixed_polygon(square, square_angles_m45, -1, FT(13)/3); + test_fixed_polygon(square, square_angles_m45, 10, FT(4630)/3); + test_fixed_polygon(square, square_angles_m45, -10, FT(4630)/3); + + // some vertical + const std::vector > square_semi_vertical_1 = { { 90, 45, 45, 45 } }; + test_fixed_polygon(square, square_semi_vertical_1, 1, FT(1)/6 + FT(1)/24); + test_fixed_polygon(square, square_semi_vertical_1, -1, FT(1)/6 + FT(1)/24); + + const std::vector > square_semi_vertical_2 = { { 90, 90, 45, 45 } }; + test_fixed_polygon(square, square_semi_vertical_2, 1, FT(1)/3); + test_fixed_polygon(square, square_semi_vertical_2, -1, FT(1)/3); + + const std::vector > square_semi_vertical_2b = { { 90, 45, 90, 45 } }; + test_fixed_polygon(square, square_semi_vertical_2b, 1, 0.25); + test_fixed_polygon(square, square_semi_vertical_2b, -1, 0.25); + + const std::vector > square_semi_vertical_3 = { { 90, 90, 90, 45 } }; + test_fixed_polygon(square, square_semi_vertical_3, 1, FT(1)/2); + test_fixed_polygon(square, square_semi_vertical_3, -1, FT(1)/2); + + // unit square, 10° inward extrusion + const std::vector > square_angles_10 = { { 10, 10, 10, 10 } }; + test_fixed_polygon(square, square_angles_10, 1, 0.0293878301); // 0.5 * tan(alpha) / 3 + test_fixed_polygon(square, square_angles_10, -1, 0.0293878301); + + // all vertical + const std::vector > square_angles_90 = { { 90, 90, 90, 90 } }; + test_fixed_polygon(square, square_angles_90, 1, 1); + test_fixed_polygon(square, square_angles_90, -1, 1); + test_fixed_polygon(square, square_angles_90, 10, 10); + test_fixed_polygon(square, square_angles_90, -10, 10); + + // larger square, with holes + square_vertices = { { 0, 0 }, { 10, 0 }, { 10, 10 }, { 0, 10 } }; + Polygon_2 large_square { square_vertices.begin(), square_vertices.end() }; + + Polygon_with_holes_2 pwh(large_square); + std::vector hole_centers = { { 2.5, 2.5 }, { 7.5, 2.5 }, { 7.5, 7.5 }, { 2.5, 7.5 } }; + for(const Point_2& hole_center : hole_centers) + { + std::vector hole_pts = { { hole_center.x() - 0.5, hole_center.y() - 0.5 }, + { hole_center.x() - 0.5, hole_center.y() + 0.5 }, + { hole_center.x() + 0.5, hole_center.y() + 0.5 }, + { hole_center.x() + 0.5, hole_center.y() - 0.5 } }; + Polygon_2 hole(hole_pts.begin(), hole_pts.end()); + pwh.add_hole(hole); + } + + // all vertical + const std::vector > large_square_angles_90 = { { 90, 90, 90, 90 }, { 90, 90, 90, 90 }, { 90, 90, 90, 90 }, { 90, 90, 90, 90 }, { 90, 90, 90, 90 } }; + test_fixed_polygon(pwh, large_square_angles_90, 1, 100 - 4 * 1); + test_fixed_polygon(pwh, large_square_angles_90, -1, 100 - 4 * 1); + test_fixed_polygon(pwh, large_square_angles_90, 10, 1000 - 4 * 10); + test_fixed_polygon(pwh, large_square_angles_90, -10, 1000 - 4 * 10); + + // outer vertical, holes outward + const std::vector > large_square_angles_outer_90 = { { 90, 90, 90, 90 }, { 135, 135, 135, 135 }, { 135, 135, 135, 135 }, { 135, 135, 135, 135 }, { 135, 135, 135, 135 } }; + test_fixed_polygon(pwh, large_square_angles_outer_90, 1, 100 - 4 * FT(1)/6); + test_fixed_polygon(pwh, large_square_angles_outer_90, -1, 100 - 4 * FT(1)/6); + test_fixed_polygon(pwh, large_square_angles_outer_90, 10, 1000 - 4 * FT(1)/6); + test_fixed_polygon(pwh, large_square_angles_outer_90, -10, 1000 - 4 * FT(1)/6); + + // outer vertical, holes inwards + const std::vector > large_square_angles_outer_90b = { { 90, 90, 90, 90 }, { 45, 45, 45, 45 }, { 45, 45, 45, 45 }, { 45, 45, 45, 45 }, { 45, 45, 45, 45 } }; + test_fixed_polygon(pwh, large_square_angles_outer_90b, 1, 100 - 4 * FT(13) / 3); // see 45° above + test_fixed_polygon(pwh, large_square_angles_outer_90b, -1, 100 - 4 * FT(13) / 3); + test_fixed_polygon(pwh, large_square_angles_outer_90b, 10, 200 - (4 * FT(62) / 3)); // h = 2 since 45° + test_fixed_polygon(pwh, large_square_angles_outer_90b, -10, 200 - (4 * FT(62) / 3)); + + // outer inwards, holes vertical + const std::vector > large_square_angles_some_90 = { { 10, 10, 10, 10 }, { 90, 90, 90, 90 }, { 90, 90, 90, 90 }, { 90, 90, 90, 90 }, { 90, 90, 90, 90 } }; + test_fixed_polygon(pwh, large_square_angles_some_90, 0.1, 8.5086282194); + test_fixed_polygon(pwh, large_square_angles_some_90, -0.1, 8.5086282194); + test_fixed_polygon(pwh, large_square_angles_some_90, 10, 27.742111630); + test_fixed_polygon(pwh, large_square_angles_some_90, -10, 27.742111630); + + // outer outwardcs, holes vertical + const std::vector > large_square_angles_some_90b = { { 170, 170, 170, 170 }, { 90, 90, 90, 90 }, { 90, 90, 90, 90 }, { 90, 90, 90, 90 }, { 90, 90, 90, 90 } }; + test_fixed_polygon(pwh, large_square_angles_some_90b, 1, 256.310219695 - 4 * 1); + test_fixed_polygon(pwh, large_square_angles_some_90b, -1, 256.310219695 - 4 * 1); + test_fixed_polygon(pwh, large_square_angles_some_90b, 10, 55227.1469423 - 4 * 10); + test_fixed_polygon(pwh, large_square_angles_some_90b, -10, 55227.1469423 - 4 * 10); + + // everyone outwards + const std::vector > large_square_angles_all_170 = { { 170, 170, 170, 170 }, { 170, 170, 170, 170 }, { 170, 170, 170, 170 }, { 170, 170, 170, 170 }, { 170, 170, 170, 170 } }; + test_fixed_polygon(pwh, large_square_angles_all_170, 1, 256.310219695 - 4 * 0.0293878301); // see 10° above + test_fixed_polygon(pwh, large_square_angles_all_170, -1, 256.310219695 - 4 * 0.0293878301); + test_fixed_polygon(pwh, large_square_angles_all_170, 10, 55227.1469423 - 4 * 0.0293878301); + test_fixed_polygon(pwh, large_square_angles_all_170, -10, 55227.1469423 - 4 * 0.0293878301); + + // everyone inwards + const std::vector > large_square_angles_all_45 = { { 45, 45, 45, 45 }, { 45, 45, 45, 45 }, { 45, 45, 45, 45 }, { 45, 45, 45, 45 }, { 45, 45, 45, 45 } }; + test_fixed_polygon(pwh, large_square_angles_all_45, 1, 64); // big truncated pyramid - 4 small truncated pyramids + test_fixed_polygon(pwh, large_square_angles_all_45, -1, 64); + test_fixed_polygon(pwh, large_square_angles_all_45, 10, FT(232)/3); + test_fixed_polygon(pwh, large_square_angles_all_45, -10, FT(232)/3); +} + +template +bool test_dat(const char* poly_filename, + const char* angles_filename, + const typename K::FT height, + const typename K::FT expected_volume) { using FT = typename K::FT; @@ -225,7 +493,9 @@ bool test(const char* poly_filename, else CGAL::extrude_skeleton(pwh, sm, CGAL::parameters::angles(angles).maximum_height(height)); -// assert(CGAL::IO::write_polygon_mesh(root_name(poly_filename) + "_extruded_up.off", sm, CGAL::parameters::stream_precision(17))); +// CGAL::IO::write_polygon_mesh(root_name(poly_filename) + "_extruded_up.off", sm, CGAL::parameters::stream_precision(17)); + + CGAL::draw(sm); FT volume = PMP::volume(sm); const FT rel_eps = 1e-5; @@ -236,7 +506,7 @@ bool test(const char* poly_filename, clear(sm); CGAL::extrude_skeleton(pwh, sm, CGAL::parameters::angles(angles).maximum_height(- height)); -// assert(CGAL::IO::write_polygon_mesh(root_name(poly_filename) + "_extruded_down.off", sm, CGAL::parameters::stream_precision(17))); +// CGAL::IO::write_polygon_mesh(root_name(poly_filename) + "_extruded_down.off", sm, CGAL::parameters::stream_precision(17)); volume = PMP::volume(sm); assert(are_equal(volume, expected_volume, rel_eps, true /*verbose*/)); @@ -244,44 +514,238 @@ bool test(const char* poly_filename, return true; } + +template +void test_random(CGAL::Random& rnd) +{ + std::cout << "\n ==== Test with Kernel: " << typeid(K).name() << " ====" << std::endl; + + using FT = typename K::FT; + using Point_2 = typename K::Point_2; + using Point_3 = typename K::Point_3; + + using Polygon_2 = CGAL::Polygon_2; + using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; + + using Mesh = CGAL::Surface_mesh; + + // ------------------------ test some random polygons ------------------------ + auto generate_random_polygon = [&](CGAL::Random& rnd) -> Polygon_2 + { + typedef CGAL::Random_points_in_square_2 Point_generator; + + Polygon_2 poly; + CGAL::random_polygon_2(10, std::back_inserter(poly), Point_generator(0.25, rnd)); + return poly; + }; + + auto test_polygon = [&](const auto& pwh) + { + std::cout << "== Test Polygon ==" << std::endl; + for(const auto& p : pwh.outer_boundary()) + std::cout << " " << p << std::endl; + + // test all slopes vertical + { + for(int i=0; i<3; ++i) + { + const FT max_height = rnd.get_double(-1, 1); + if(max_height == 0) + continue; + + std::cout << "max_height = " << max_height << std::endl; + + std::vector > weights; + weights.push_back(std::vector(pwh.outer_boundary().size(), 0)); + for(auto hit=pwh.holes_begin(); hit!=pwh.holes_end(); ++hit) + weights.push_back(std::vector(hit->size(), 0)); + + Mesh sm; + bool success = extrude_skeleton(pwh, sm, CGAL::parameters::weights(weights).verbose(true).maximum_height(max_height)); + assert(success); + if(!success) + { + std::cerr << "Error: failed to extrude skeleton" << std::endl; + assert(false); + std::exit(1); + } + + CGAL::draw(sm); + CGAL::IO::write_OFF("last_extrusion.off", sm, CGAL::parameters::stream_precision(17)); + + // check the volume + FT expected_vol = CGAL::abs(pwh.outer_boundary().area() * max_height); + FT actual_vol = CGAL::Polygon_mesh_processing::volume(sm); + if (CGAL::abs(expected_vol - actual_vol) > 1e-6) + { + std::cerr << "Error: expected volume " << expected_vol << " but got " << actual_vol << std::endl; + assert(false); + std::exit(1); + } + + // for each point at z=0, there should be the same point at z=max_height + std::map, std::list > zs; + for (auto v : vertices(sm)) + { + Point_3 p = sm.point(v); + if (p.z() != 0 && p.z() != max_height) + { + std::cerr << "Error: point at z=" << p.z() << " instead of 0 or " << max_height << std::endl; + assert(false); + std::exit(1); + } + + zs[std::make_pair(p.x(), p.y())].push_back(p); + } + + for (const auto& e : zs) + { + if (e.second.size() != 2) + { + std::cerr << "Error: point at (" << e.first.first << ", " << e.first.second << ") has " << e.second.size() << " points" << std::endl; + assert(false); + std::exit(1); + } + + if (e.second.front().z() == e.second.back().z()) + { + std::cerr << "Error: point at (" << e.first.first << ", " << e.first.second << ") has two points at the same height" << std::endl; + assert(false); + std::exit(1); + } + } + } + } + + // test some slopes vertical + { + for(int i=0; i<3; ++i) + { + const FT max_height = rnd.get_double(-1, 1); + if(max_height == 0) + continue; + + std::vector > weights; + weights.push_back(std::vector(pwh.outer_boundary().size(), 0)); + for(auto hit=pwh.holes_begin(); hit!=pwh.holes_end(); ++hit) + weights.push_back(std::vector(hit->size(), 0)); + + for(auto& contour_weights : weights) + { + for(FT& weight : contour_weights) + { + weight = rnd.get_double(0.05, 0.5); + if (!rnd.get_int(0, 3)) { // a third to zero + weight = 0; + } + } + } + + // randomly switch inwards/outwards + bool outwards = false; + if(rnd.get_int(0, 2)) + { + outwards = true; + for(std::vector& ws : weights) + for(FT& w : ws) + w = -w; + } + + Mesh sm; + bool success = extrude_skeleton(pwh, sm, CGAL::parameters::weights(weights) + .verbose(true) + .maximum_height(max_height)); + if(!success) + { + std::cerr << "Error: failed to extrude skeleton" << std::endl; + assert(false); + std::exit(1); + } + + CGAL::draw(sm); + CGAL::IO::write_OFF("last_extrusion.off", sm, CGAL::parameters::stream_precision(17)); + + if(!outwards) + { + success = extrude_skeleton(pwh, sm, CGAL::parameters::weights(weights).verbose(true)); + if(!success) + { + std::cerr << "Error: failed to extrude skeleton" << std::endl; + assert(false); + std::exit(1); + } + + CGAL::draw(sm); + CGAL::IO::write_OFF("last_extrusion.off", sm, CGAL::parameters::stream_precision(17)); + } + } + } + }; + + // Random simple polygon (no holes) + for (int i=0; i<10; ++i) { + Polygon_2 poly = generate_random_polygon(rnd); + Polygon_with_holes_2 pwh(poly); + test_polygon(pwh); + } + + // Random simple polygons with holes + +} + template -void test() +void test_dats() { - test("data/polygon_000.dat", "data/angles_000.dat", 6, 162.37987499999997); - test("data/polygon_001.dat", "data/angles_001.dat", 6, 761.76899999999989); - test("data/polygon_002.dat", "data/angles_002.dat", 22, 15667.658890389464); - test("data/polygon_003.dat", "data/angles_003.dat", 12, 105.79864291299999); - test("data/polygon_004.dat", "data/angles_004.dat", 12, 3119.9357499857151); - test("data/polygon_005.dat", "data/angles_005.dat", 12, 1342.9474791424002); - test("data/polygon_006.dat", "data/angles_006.dat", 12, 249.41520000000008); - test("data/polygon_007.dat", "data/angles_007.dat", 12, 7344.8312073148918); - test("data/polygon_008.dat", "data/angles_008.dat", 12, 7240.6890039677555); - test("data/polygon_009.dat", "data/angles_009.dat", 10, 3704.0787987580375); - test("data/polygon_010.dat", "data/angles_010.dat", 40, 29306.453333333335); - test("data/polygon_011.dat", "data/angles_011.dat", 40, 375866.54633629462); - test("data/polygon_012.dat", "data/angles_012.dat", 12, 4560.0268861722925); - test("data/polygon_013.dat", "data/angles_013.dat", 12, 2221.3622594712501); - test("data/polygon_014.dat", "data/angles_014.dat", 12, 4534.0568515270861); - test("data/polygon_015.dat", "data/angles_015.dat", 12, 1565.5667825255343); - test("data/polygon_016.dat", "data/angles_016.dat", 311, 2518611984.6277928); - test("data/polygon_017.dat", "data/angles_017.dat", 50, 7729166.666666667); - test("data/polygon_018.dat", "data/angles_018.dat", 50, 354166.66666666663); - test("data/polygon_019.dat", "data/angles_019.dat", 311, 92570921.033775225); - test("data/polygon_020.dat", "data/angles_020.dat", 70, 1550161.8131298050); - test("data/polygon_021.dat", "data/angles_021.dat", 70, 37631800.885042846); - test("data/polygon_022.dat", "data/angles_022.dat", 20, 7702444.2118858183); + test_dat("data/polygon_000.dat", "data/angles_000.dat", 6, 162.37987499999997); + test_dat("data/polygon_001.dat", "data/angles_001.dat", 6, 761.76899999999989); + test_dat("data/polygon_002.dat", "data/angles_002.dat", 22, 15667.658890389464); + test_dat("data/polygon_003.dat", "data/angles_003.dat", 12, 105.79864291299999); + test_dat("data/polygon_004.dat", "data/angles_004.dat", 12, 3119.9357499857151); + test_dat("data/polygon_005.dat", "data/angles_005.dat", 12, 1342.9474791424002); + test_dat("data/polygon_006.dat", "data/angles_006.dat", 12, 249.41520000000008); + test_dat("data/polygon_007.dat", "data/angles_007.dat", 12, 7344.8312073148918); + test_dat("data/polygon_008.dat", "data/angles_008.dat", 12, 7240.6890039677555); + test_dat("data/polygon_009.dat", "data/angles_009.dat", 10, 3704.0787987580375); + test_dat("data/polygon_010.dat", "data/angles_010.dat", 40, 29306.453333333335); + test_dat("data/polygon_011.dat", "data/angles_011.dat", 40, 375866.54633629462); + test_dat("data/polygon_012.dat", "data/angles_012.dat", 12, 4560.0268861722925); + test_dat("data/polygon_013.dat", "data/angles_013.dat", 12, 2221.3622594712501); + test_dat("data/polygon_014.dat", "data/angles_014.dat", 12, 4534.0568515270861); + test_dat("data/polygon_015.dat", "data/angles_015.dat", 12, 1565.5667825255343); + test_dat("data/polygon_016.dat", "data/angles_016.dat", 311, 2518611984.6277928); + test_dat("data/polygon_017.dat", "data/angles_017.dat", 50, 7729166.666666667); + test_dat("data/polygon_018.dat", "data/angles_018.dat", 50, 354166.66666666663); + test_dat("data/polygon_019.dat", "data/angles_019.dat", 311, 92570921.033775225); + test_dat("data/polygon_020.dat", "data/angles_020.dat", 70, 1550161.8131298050); + test_dat("data/polygon_021.dat", "data/angles_021.dat", 70, 37631800.885042846); + test_dat("data/polygon_022.dat", "data/angles_022.dat", 20, 7702444.2118858183); } -int main(int, char**) + +template +void test(CGAL::Random& rnd) +{ + test_fixed_with_vertical_combinations(); + // test_fixed(); + // test_dats(); + // test_random(rnd); +} + +int main(int argc, char** argv) { std::cout.precision(17); std::cerr.precision(17); - test(); - test(); - test(); + int seed = (argc > 0) ? std::atoi(argv[0]) : 0; // @fixme std::time(nullptr) + + CGAL::Random rnd(seed); + std::cout << "Seed is " << rnd.get_seed() << std::endl; + + test(rnd); + test(rnd); + test(rnd); - std::cout << "OK" << std::endl; + std::cout << "Done" << std::endl; return EXIT_SUCCESS; }