From 633fff03ca9a051df0814e5437ba149d35d43917 Mon Sep 17 00:00:00 2001 From: Kenny Weiss Date: Mon, 29 Nov 2021 09:04:04 -0800 Subject: [PATCH] WIP: Fixing curved polygon intersection algorithm for special cases --- .../detail/intersect_bezier_impl.hpp | 21 +- .../detail/intersect_curved_poly_impl.hpp | 304 +++++++++++++----- src/axom/primal/operators/intersect.hpp | 14 +- src/axom/primal/operators/orientation.hpp | 4 +- .../tests/primal_curved_polygon_intersect.cpp | 42 ++- 5 files changed, 282 insertions(+), 103 deletions(-) diff --git a/src/axom/primal/operators/detail/intersect_bezier_impl.hpp b/src/axom/primal/operators/detail/intersect_bezier_impl.hpp index b467cb302c..f9f846d6f4 100644 --- a/src/axom/primal/operators/detail/intersect_bezier_impl.hpp +++ b/src/axom/primal/operators/detail/intersect_bezier_impl.hpp @@ -152,6 +152,7 @@ bool intersect_bezier_curves(const BezierCurve &c1, s_scale *= scaleFac; + // Note: we want to find all intersections, so don't short-circuit if(intersect_bezier_curves(c2, c3, tp, @@ -185,11 +186,11 @@ bool intersect_bezier_curves(const BezierCurve &c1, return foundIntersection; } -template -bool intersect_2d_linear(const Point &a, - const Point &b, - const Point &c, - const Point &d, +template +bool intersect_2d_linear(const Point &a, + const Point &b, + const Point &c, + const Point &d, T &s, T &t) { @@ -199,18 +200,16 @@ bool intersect_2d_linear(const Point &a, // Note: Uses exact floating point comparisons since the subdivision algorithm // provides both sides of the line segments for interior curve points. - AXOM_STATIC_ASSERT(NDIMS == 2); - // compute signed areas of endpoints of segment (c,d) w.r.t. segment (a,b) - auto area1 = twoDcross(a, b, c); - auto area2 = twoDcross(a, b, d); + const auto area1 = twoDcross(a, b, c); + const auto area2 = twoDcross(a, b, d); // early return if both have same orientation, or if d is collinear w/ (a,b) if(area2 == 0. || (area1 * area2) > 0.) return false; // compute signed areas of endpoints of segment (a,b) w.r.t. segment (c,d) - auto area3 = twoDcross(c, d, a); - auto area4 = area3 + area1 - area2; // equivalent to twoDcross(c,d,b) + const auto area3 = twoDcross(c, d, a); + const auto area4 = area3 + area1 - area2; // equivalent to twoDcross(c,d,b) // early return if both have same orientation, or if b is collinear w/ (c,d) if(area4 == 0. || (area3 * area4) > 0.) return false; diff --git a/src/axom/primal/operators/detail/intersect_curved_poly_impl.hpp b/src/axom/primal/operators/detail/intersect_curved_poly_impl.hpp index 4249c369c2..d8ae24bc21 100644 --- a/src/axom/primal/operators/detail/intersect_curved_poly_impl.hpp +++ b/src/axom/primal/operators/detail/intersect_curved_poly_impl.hpp @@ -12,6 +12,8 @@ #ifndef AXOM_PRIMAL_INTERSECTION_CURVED_POLYGON_IMPL_HPP_ #define AXOM_PRIMAL_INTERSECTION_CURVED_POLYGON_IMPL_HPP_ +#include "axom/core.hpp" + #include "axom/primal/geometry/BoundingBox.hpp" #include "axom/primal/geometry/OrientedBoundingBox.hpp" #include "axom/primal/geometry/Point.hpp" @@ -22,8 +24,6 @@ #include "axom/primal/geometry/BezierCurve.hpp" #include "axom/primal/geometry/CurvedPolygon.hpp" -#include "axom/core/utilities/Utilities.hpp" - #include "axom/primal/operators/squared_distance.hpp" #include "axom/primal/operators/detail/intersect_bezier_impl.hpp" @@ -43,21 +43,21 @@ int isContained(const CurvedPolygon& p1, const CurvedPolygon& p2, double sq_tol = 1e-10); -template +template class DirectionalWalk { public: + static constexpr int NDIMS = 2; using CurvedPolygonType = CurvedPolygon; - using EdgeLabels = std::vector; + using IndexArray = std::vector; public: /*! \class EdgeIntersectionInfo * * \brief For storing intersection points between edges of \a CurvedPolygon instances so they can be easily sorted by parameter value using std::sort */ - class EdgeIntersectionInfo + struct EdgeIntersectionInfo { - public: T myTime; // parameter value of intersection on curve on first CurvePolygon int myEdge; // index of curve on first CurvedPolygon T otherTime; // parameter value of intersection on curve on second CurvedPolygon @@ -78,7 +78,7 @@ class DirectionalWalk NON_JUNCTION, /// Not a junction, e.g. a vertex of the original APPLIED, /// Use after junction has been applied CROSS, /// Typical case: edges of polygons intersect - GRAZE /// Atypical case: polygons intersecte at a terminal vertex + GRAZE /// Atypical case: polygons intersect at a terminal vertex }; /*! @@ -90,23 +90,31 @@ class DirectionalWalk using EdgeIndex = int; static const EdgeIndex INVALID_EDGE_INDEX = -1; static const int INVALID_JUNCTION_INDEX = -1; - - // indices of polygon edges leading into junction - EdgeIndex edgeIndex[2] {INVALID_EDGE_INDEX, INVALID_EDGE_INDEX}; - // describes the junction type and status - JunctionState junctionState {JunctionState::UNINITIALIZED}; - JunctionIndex index {INVALID_JUNCTION_INDEX}; + static const int NON_JUNCTION_INDEX = 0; bool isActive() const { return junctionState > JunctionState::APPLIED; } - EdgeIndex currentEdgeIndex(bool active) const { return edgeIndex[active]; } - EdgeIndex nextEdgeIndex(bool active) const { return edgeIndex[active] + 1; } + EdgeIndex currentEdgeIndex(bool active) const + { + return incomingEdgeIndex[active]; + } + EdgeIndex nextEdgeIndex(bool active) const + { + return incomingEdgeIndex[active] + 1; + } bool operator==(const Junction& other) const { return index == other.index; } bool operator!=(const Junction& other) const { return !(*this == other); } + + public: + // indices of polygon edges leading into junction + EdgeIndex incomingEdgeIndex[2] {INVALID_EDGE_INDEX, INVALID_EDGE_INDEX}; + // describes the junction type and status + JunctionState junctionState {JunctionState::UNINITIALIZED}; + JunctionIndex index {INVALID_JUNCTION_INDEX}; }; /*! @@ -119,10 +127,10 @@ class DirectionalWalk using VertexIndex = int; using EdgeIndex = int; - PolygonEdge(int polygonID, const EdgeLabels& labels) + PolygonEdge(int polygonID, const IndexArray& endpointJunctionIndices) : m_polygonID(polygonID) - , m_labels(labels) - , m_numEdges(labels.size()) + , m_endpointJunctionIndices(endpointJunctionIndices) + , m_numEdges(endpointJunctionIndices.size()) { SLIC_ASSERT_MSG(m_numEdges > 0, "Polygon " << polygonID << " was empty"); } @@ -130,8 +138,8 @@ class DirectionalWalk EdgeIndex getIndex() const { return m_index; } void setIndex(EdgeIndex idx) { m_index = (idx % m_numEdges); } - int getStartLabel() const { return m_labels[prevIndex()]; } - int getEndLabel() const { return m_labels[m_index]; } + int getStartLabel() const { return m_endpointJunctionIndices[prevIndex()]; } + int getEndLabel() const { return m_endpointJunctionIndices[m_index]; } void advance() { m_index = nextIndex(); } bool isJunction() const { return getEndLabel() > 0; } @@ -158,17 +166,17 @@ class DirectionalWalk EdgeIndex m_index; const int m_polygonID; - const EdgeLabels& m_labels; + const IndexArray& m_endpointJunctionIndices; const int m_numEdges; }; public: - explicit DirectionalWalk(bool verbose = false) : m_verbose(verbose) { } + explicit DirectionalWalk(bool verbose = true) : m_verbose(verbose) { } /*! * Splits edges of each polygon based on the intersections with the other polygon * When the two polygons intersect, the split polygons are returned in \a psplit - * The edges are labeled by the types of intersections in \a edgelabels + * The edges are labeled by the types of intersections in \a endpointJunctionIndices * and \a orientation contains the relative orientation of the first intersection */ int splitPolygonsAlongIntersections(const CurvedPolygonType& p1, @@ -176,8 +184,10 @@ class DirectionalWalk double sq_tol) { // We store intersection information for each edge in EdgeIntersectionInfo structures - std::vector> E1IntData(p1.numEdges()); - std::vector> E2IntData(p2.numEdges()); + using EdgeIntersectionInfoArray = + std::vector>; + EdgeIntersectionInfoArray p1IntersectionData(p1.numEdges()); + EdgeIntersectionInfoArray p2IntersectionData(p2.numEdges()); EdgeIntersectionInfo firstinter; // Need to do orientation test on first intersection // Find all intersections and store @@ -209,8 +219,10 @@ class DirectionalWalk for(int k = 0; k < edgeIntersections; ++k, ++numinters) { - E1IntData[i].push_back({p1times[k], i, p2times[k], j, numinters + 1}); - E2IntData[j].push_back({p2times[k], j, p1times[k], i, numinters + 1}); + p1IntersectionData[i].push_back( + {p1times[k], i, p2times[k], j, numinters + 1}); + p2IntersectionData[j].push_back( + {p2times[k], j, p1times[k], i, numinters + 1}); if(m_verbose) { @@ -239,18 +251,18 @@ class DirectionalWalk for(int i = 0; i < p1.numEdges(); ++i) { - std::sort(E1IntData[i].begin(), E1IntData[i].end()); + std::sort(p1IntersectionData[i].begin(), p1IntersectionData[i].end()); } for(int i = 0; i < p2.numEdges(); ++i) { - std::sort(E2IntData[i].begin(), E2IntData[i].end()); + std::sort(p2IntersectionData[i].begin(), p2IntersectionData[i].end()); } psplit[0] = p1; psplit[1] = p2; - edgelabels[0].reserve(p1.numEdges() + numinters); - edgelabels[1].reserve(p2.numEdges() + numinters); + endpointJunctionIndices[0].reserve(p1.numEdges() + numinters); + endpointJunctionIndices[1].reserve(p2.numEdges() + numinters); junctions = std::vector(numinters + 1); @@ -260,8 +272,8 @@ class DirectionalWalk SLIC_INFO("Poly2 before split: " << psplit[1]); } - splitPolygon(psplit[0], E1IntData, edgelabels[0], 0); - splitPolygon(psplit[1], E2IntData, edgelabels[1], 1); + this->splitPolygon(0, p1IntersectionData); + this->splitPolygon(1, p2IntersectionData); if(m_verbose) { @@ -279,14 +291,31 @@ class DirectionalWalk * * This function must be called after splitPolygonsAlongIntersections * which inserts the intersection points into each polygon and generates - * the \a edgelabels structure. + * the \a endpointJunctionIndices structure. * * The results will be a vector of polygons inserted into OUT parameter \a pnew */ void findIntersectionRegions(std::vector& pnew) { - PolygonEdge currentEdge[2] = {PolygonEdge(0, edgelabels[0]), - PolygonEdge(1, edgelabels[1])}; + PolygonEdge currentEdge[2] = {PolygonEdge(0, endpointJunctionIndices[0]), + PolygonEdge(1, endpointJunctionIndices[1])}; + + if(m_verbose) + { + SLIC_INFO("At start of 'findIntersectionRegions', there are " + << junctions.size() << " junctions:"); + for(int i = 0; i < junctions.size(); ++i) + { + SLIC_INFO( + axom::fmt::format("\tJunction {}: incomingEdge[0]: {}, " + "incomingEdge[1]: {}, state: {}, index: {} ", + i, + junctions[i].incomingEdgeIndex[0], + junctions[i].incomingEdgeIndex[1], + junctions[i].junctionState, + junctions[i].index)); + } + } // We use junctionIndex to loop through the active junctions starting with index 1 int junctionIndex = 1; @@ -309,6 +338,15 @@ class DirectionalWalk { break; } + // switch(startJunction->junctionState) + // { + // case JunctionState::CROSS: + // break; + // case JunctionState::GRAZE: + // break; + // default: + // break; + // } else { startJunction->junctionState = JunctionState::APPLIED; @@ -324,17 +362,28 @@ class DirectionalWalk if(m_verbose) { - SLIC_INFO("" - << "Starting with junction " << startJunction->index - << "\n\t edges[0] {in: " << startJunction->currentEdgeIndex(0) - << " -- " << psplit[0][startJunction->currentEdgeIndex(0)] - << "; out: " << startJunction->nextEdgeIndex(0) << " -- " - << psplit[0][startJunction->nextEdgeIndex(0)] << "}" - << "\n\t edges[1] {in: " << startJunction->currentEdgeIndex(1) - << " -- " << psplit[1][startJunction->currentEdgeIndex(1)] - << "; out: " << startJunction->nextEdgeIndex(1) << " -- " - << psplit[1][startJunction->nextEdgeIndex(1)] << "}" - << "\n\t active: " << (active ? 1 : 0)); + const int p0Edges = psplit[0].numEdges(); + const int p1Edges = psplit[1].numEdges(); + const int p0CurrIdx = startJunction->currentEdgeIndex(0); + const int p0NextIdx = startJunction->nextEdgeIndex(0) % p0Edges; + const int p1CurrIdx = startJunction->currentEdgeIndex(1); + const int p1NextIdx = startJunction->nextEdgeIndex(1) % p1Edges; + + SLIC_INFO( + axom::fmt::format("Starting with junction {}" + "\n\t edges[0] {} -> {} {{in: {}; out: {} }}" + "\n\t edges[1] {} -> {} {{in: {}; out: {} }}" + "\n\t active: {}", + startJunction->index, + p0CurrIdx, + p0NextIdx, + psplit[0][p0CurrIdx], + psplit[0][p0NextIdx], + p1CurrIdx, + p1NextIdx, + psplit[1][p1CurrIdx], + psplit[1][p1NextIdx], + (active ? 1 : 0))); } CurvedPolygonType aPart; // Tracks the edges of the current intersection polygon @@ -408,65 +457,138 @@ class DirectionalWalk } private: - /// Splits a CurvedPolygon \a p1 at every intersection point stored in \a IntersectionData - /// \a edgeLabels stores intersection id for new vertices and 0 for original vertices - void splitPolygon(CurvedPolygon& polygon, - std::vector>& IntersectionData, - std::vector& edgelabels, - int polygonID) + /// Splits the polygon with id \a polygonID at every intersection point in \a edgeIntersections. + /// Uses internal array \a edgeLabels to store ids for vertices (junction indices or 0 for original vertices) + void splitPolygon( + int polygonID, + std::vector>& allEdgeIntersections) { using axom::utilities::isNearlyEqual; + CurvedPolygonType& polygon = psplit[polygonID]; + IndexArray& junctionIndices = endpointJunctionIndices[polygonID]; + + bool fixupBeginning = false; + int fixupBeginningJunctionIdx = Junction::INVALID_JUNCTION_INDEX; + int addedIntersections = 0; - const int numEdges = polygon.numEdges(); - for(int i = 0; i < numEdges; ++i) // foreach edge + const int numOrigEdges = polygon.numEdges(); + for(int i = 0; i < numOrigEdges; ++i) // foreach edge { - edgelabels.push_back(0); // mark start current vertex as 'original' - const int nIntersect = IntersectionData[i].size(); + // mark this edge's endpoint as 'original' + junctionIndices.push_back(Junction::NON_JUNCTION_INDEX); + double previous_tj = 0.; + auto& curEdgeIntersections = allEdgeIntersections[i]; + const int nIntersect = curEdgeIntersections.size(); for(int j = 0; j < nIntersect; ++j) // foreach intersection on this edge { // split edge at parameter t_j - const double t_j = IntersectionData[i][j].myTime; + const double t_j = curEdgeIntersections[j].myTime; const int edgeIndex = i + addedIntersections; + if(m_verbose) + { + SLIC_INFO( + fmt::format("i {}, j {}, added {}, t_j {}, previous t_j {}, " + "edgeIndex {} (sz: {}), polygon size {}", + i, + j, + addedIntersections, + t_j, + previous_tj, + edgeIndex, + junctionIndices.size(), + polygon.numEdges())); + } + + const bool nearlyZero = isNearlyEqual(t_j, 0.); + // TODO: Handle case where t_j is nearly 0. + if(!nearlyZero) + { + polygon.splitEdge(edgeIndex, t_j); + + // update edge label + const int idx = curEdgeIntersections[j].numinter; + junctionIndices.insert(junctionIndices.begin() + edgeIndex, idx); + junctions[idx].incomingEdgeIndex[polygonID] = edgeIndex; + //if(junctions[idx].junctionState == JunctionState::UNINITIALIZED) + //{ + junctions[idx].junctionState == JunctionState::CROSS; + //} + junctions[idx].index = idx; + + ++addedIntersections; + previous_tj = t_j; + } + else + { + const int idx = curEdgeIntersections[j].numinter; - polygon.splitEdge(edgeIndex, t_j); + if(edgeIndex > 0) + { + junctionIndices[edgeIndex] = idx; + junctions[idx].incomingEdgeIndex[polygonID] = edgeIndex; + } + else + { + fixupBeginning = true; + fixupBeginningJunctionIdx = idx; + } - // update edge label - const int label = IntersectionData[i][j].numinter; - edgelabels.insert(edgelabels.begin() + edgeIndex, label); - junctions[label].edgeIndex[polygonID] = edgeIndex; - junctions[label].junctionState = - JunctionState::CROSS; // TODO: Figure out what junction type needs to be - junctions[label].index = label; + junctions[idx].junctionState = + JunctionState::CROSS; //JunctionState::GRAZE; + junctions[idx].index = idx; + } - ++addedIntersections; + if(m_verbose) + { + SLIC_INFO(axom::fmt::format("junctionIndices: {}\n junctions:", + axom::fmt::join(junctionIndices, " "))); + for(const auto& jj : junctions) + { + SLIC_INFO(fmt::format( + "\t{{index: {}, state: {}, edgeIndex[0]: {}, edgeIndex[1]: {} }}", + jj.index, + jj.junctionState, + jj.incomingEdgeIndex[0], + jj.incomingEdgeIndex[1])); + } + } // update remaining intersections on this edge; special case if already at end of curve if(!isNearlyEqual(1., t_j)) { for(int k = j + 1; k < nIntersect; ++k) { - const double t_k = IntersectionData[i][k].myTime; - IntersectionData[i][k].myTime = (t_k - t_j) / (1.0 - t_j); + const double t_k = curEdgeIntersections[k].myTime; + curEdgeIntersections[k].myTime = (t_k - t_j) / (1.0 - t_j); } } else { for(int k = j + 1; k < nIntersect; ++k) { - IntersectionData[i][k].myTime = 1.; + curEdgeIntersections[k].myTime = 1.; } } } } + + // If the first vertex of the first edge was a junction, we need to fix it up using the current last edge + if(fixupBeginning) + { + const int edgeIndex = junctionIndices.size() - 1; + junctionIndices[edgeIndex] = fixupBeginningJunctionIdx; + junctions[fixupBeginningJunctionIdx].incomingEdgeIndex[polygonID] = + edgeIndex; + } } public: // Objects to store completely split polygons (split at every intersection point) and vector with unique id for each intersection and zeros for corners of original polygons. CurvedPolygonType psplit[2]; // The two completely split polygons will be stored in this array - EdgeLabels edgelabels[2]; // 0 for curves that end in original vertices, unique id for curves that end in intersection points + IndexArray endpointJunctionIndices[2]; // 0 for curves that end in original vertices, unique id for curves that end in intersection points std::vector junctions; int numinters {0}; bool orientation {false}; @@ -482,10 +604,10 @@ class DirectionalWalk * \param [in] sq_tol tolerance parameter for the base case of intersect_bezier_curve * \param [out] pnew vector of type CurvedPolygon holding CurvedPolygon objects representing boundaries of intersection regions. */ -template -bool intersect_polygon(const CurvedPolygon& p1, - const CurvedPolygon& p2, - std::vector>& pnew, +template +bool intersect_polygon(const CurvedPolygon& p1, + const CurvedPolygon& p2, + std::vector>& pnew, double sq_tol) { // Intersection is empty if either of the two polygons are empty @@ -494,7 +616,7 @@ bool intersect_polygon(const CurvedPolygon& p1, return false; } - DirectionalWalk walk; + DirectionalWalk walk; int numinters = walk.splitPolygonsAlongIntersections(p1, p2, sq_tol); @@ -637,6 +759,36 @@ bool orient(const BezierCurve& c1, const BezierCurve& c2, T s, T t) return (orientation > 0); } +// when interesection is at a vertex of the original polygon(s), +// determine if the in/out edges of second polygon is on the same +// side of the first polygon +template +int getJunctionIntersectionType(const BezierCurve& p0In, + const BezierCurve& p0Out, + const BezierCurve& p1In, + const BezierCurve& p1Out) +{ + using SegmentType = primal::Segment; + const auto p0Tangents[2] = {-p0In.dt(1.).unitVector(), + p0Out.dt(0.).unitVector()}; + const SegmentType p0Seg {p0Tangents[0], p0Tangents[1]}; + + const auto p1Tangents[2] = {-p1In.dt(1.).unitVector(), + p1Out.dt(0.).unitVector()}; + + const int orientIn = primal::orientation(p1Tangents[0], p0Seg); + const int orientOut = primal::orientation(p1Tangents[1], p0Seg); + + if(orientIn == orientOut) + { + return (orientIn == ON_POSITIVE_SIDE) ? 1 : -1; + } + else + { + return (orientOut == ON_POSITIVE_SIDE) ? 0 : -1; + } +} + } // namespace detail } // namespace primal } // namespace axom diff --git a/src/axom/primal/operators/intersect.hpp b/src/axom/primal/operators/intersect.hpp index 9ab5f8b624..caad3a63c0 100644 --- a/src/axom/primal/operators/intersect.hpp +++ b/src/axom/primal/operators/intersect.hpp @@ -505,10 +505,10 @@ bool intersect(const OrientedBoundingBox& b1, * Thus, component curves do not intersect at \f$ s==1 \f$ or at * \f$ t==1 \f$. */ -template -bool intersect(const CurvedPolygon& p1, - const CurvedPolygon& p2, - std::vector>& pnew, +template +bool intersect(const CurvedPolygon& p1, + const CurvedPolygon& p2, + std::vector>& pnew, double tol = 1E-8) { // for efficiency, linearity check actually uses a squared tolerance @@ -546,9 +546,9 @@ bool intersect(const CurvedPolygon& p1, * contain their first endpoint, but not their last endpoint. Thus, the * curves do not intersect at \f$ s==1 \f$ or at \f$ t==1 \f$. */ -template -bool intersect(const BezierCurve& c1, - const BezierCurve& c2, +template +bool intersect(const BezierCurve& c1, + const BezierCurve& c2, std::vector& sp, std::vector& tp, double tol = 1E-8) diff --git a/src/axom/primal/operators/orientation.hpp b/src/axom/primal/operators/orientation.hpp index d00d665e5d..0aae6cef49 100644 --- a/src/axom/primal/operators/orientation.hpp +++ b/src/axom/primal/operators/orientation.hpp @@ -60,7 +60,7 @@ inline int orientation(const Point& p, const Triangle& tri) { orient = ON_BOUNDARY; } - else if(det < 0.0f) + else if(det < 0.) { // outside orient = ON_POSITIVE_SIDE; @@ -109,7 +109,7 @@ inline int orientation(const Point& p, const Segment& seg) // collinear orient = ON_BOUNDARY; } - else if(det < 0.0f) + else if(det < 0.) { // outside, clockwise, to the right orient = ON_POSITIVE_SIDE; diff --git a/src/axom/primal/tests/primal_curved_polygon_intersect.cpp b/src/axom/primal/tests/primal_curved_polygon_intersect.cpp index b70924ebdd..22efe85a5e 100644 --- a/src/axom/primal/tests/primal_curved_polygon_intersect.cpp +++ b/src/axom/primal/tests/primal_curved_polygon_intersect.cpp @@ -170,7 +170,7 @@ void checkIntersection( //Compute intersection using algorithm with tolerance of eps intersect(bPolygon1, bPolygon2, intersectionPolys, eps); //Check that expected number of intersection regions are found - EXPECT_EQ(expbPolygon.size(), intersectionPolys.size()); + ASSERT_EQ(expbPolygon.size(), intersectionPolys.size()); //Check that expected intersection curves are found to within test_eps const int nPolygons = expbPolygon.size(); @@ -244,6 +244,32 @@ primal::CurvedPolygon createPolygon( return bPolygon; } +//---------------------------------------------------------------------------------- + +TEST(primal_curvedpolygon, detail_intersection_type) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test intersection of two squares (single region)"); + + std::vector orders = {1, 1; + + std::vector CP = {PointType {-1, 0}, + PointType {0, 0}, + PointType {1, 0}}; + CurvedPolygonType bPolygon1 = createPolygon(CP, orders); + + std::vector CP2 = {PointType {1, -1}, + PointType {0, 0}, + PointType {-1, -1}}; + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders); + + + + //---------------------------------------------------------------------------------- TEST(primal_curvedpolygon, intersection_squares) @@ -382,7 +408,7 @@ TEST(primal_curvedpolygon, intersections_triangle_rectangle) // No intersection case: Triangle is below the rectangle { CurvedPolygonType bTriangle = createPolygon(triPts, triOrders); - primal::detail::DirectionalWalk walk(bVerbose); + primal::detail::DirectionalWalk walk(bVerbose); int nIntersections = walk.splitPolygonsAlongIntersections(bRectangle, bTriangle, SQ_EPS); EXPECT_EQ(0, nIntersections); @@ -401,7 +427,7 @@ TEST(primal_curvedpolygon, intersections_triangle_rectangle) triPts[2][1] = 0; CurvedPolygonType bTriangle = createPolygon(triPts, triOrders); - primal::detail::DirectionalWalk walk(bVerbose); + primal::detail::DirectionalWalk walk(bVerbose); int nIntersections = walk.splitPolygonsAlongIntersections(bRectangle, bTriangle, SQ_EPS); EXPECT_EQ(1, nIntersections); @@ -429,7 +455,7 @@ TEST(primal_curvedpolygon, intersections_triangle_rectangle) triPts[2][1] = 1; CurvedPolygonType bTriangle = createPolygon(triPts, triOrders); - primal::detail::DirectionalWalk walk(bVerbose); + primal::detail::DirectionalWalk walk(bVerbose); int nIntersections = walk.splitPolygonsAlongIntersections(bRectangle, bTriangle, SQ_EPS); EXPECT_EQ(2, nIntersections); @@ -460,7 +486,7 @@ TEST(primal_curvedpolygon, intersections_triangle_rectangle) triPts[2][1] = 2; CurvedPolygonType bTriangle = createPolygon(triPts, triOrders); - primal::detail::DirectionalWalk walk(bVerbose); + primal::detail::DirectionalWalk walk(bVerbose); int nIntersections = walk.splitPolygonsAlongIntersections(bRectangle, bTriangle, SQ_EPS); EXPECT_EQ(3, nIntersections); @@ -486,7 +512,7 @@ TEST(primal_curvedpolygon, intersections_triangle_rectangle) triPts[2][1] = 3; CurvedPolygonType bTriangle = createPolygon(triPts, triOrders); - primal::detail::DirectionalWalk walk(bVerbose); + primal::detail::DirectionalWalk walk(bVerbose); int nIntersections = walk.splitPolygonsAlongIntersections(bRectangle, bTriangle, SQ_EPS); EXPECT_EQ(4, nIntersections); @@ -698,7 +724,7 @@ TEST(primal_curvedpolygon, doubleIntersection) double walkIntersectionArea = 0.; { const bool verbose = true; - primal::detail::DirectionalWalk walk(verbose); + primal::detail::DirectionalWalk walk(verbose); int numIntersections = walk.splitPolygonsAlongIntersections(bPolygon1, bPolygon2, EPS * EPS); @@ -873,6 +899,7 @@ TEST(primal_curvedpolygon, adjacent_squares) } // Sharing a vertex + if(false) { std::vector CP1 = {PointType {-1, 0}, PointType {0, 0}, @@ -911,6 +938,7 @@ TEST(primal_curvedpolygon, adjacent_squares) } } + if(false) { std::vector expIntersections; intersect(bPolygon2, bPolygon1, expIntersections, EPS);