diff --git a/surface_integral/include/polygon_winding_number.h b/surface_integral/include/polygon_winding_number.h new file mode 100644 index 0000000..be73913 --- /dev/null +++ b/surface_integral/include/polygon_winding_number.h @@ -0,0 +1,60 @@ +// polygon_winding.h +#pragma once + +#include +#include +#include +#include + +// Type aliases for convenience +using Point2d = Eigen::Vector2d; +using PolygonRing = std::vector; + +/** + * @brief Determines if a point is inside a polygon with holes using the non-zero winding rule. + * + * @param P The query point. + * @param outer_ring The outer boundary of the polygon (should be counter-clockwise). + * @param holes A list of inner rings (holes), each should be clockwise. + * @return true if the point is inside (including boundary), false otherwise. + */ +bool point_in_polygon_with_holes( + const Point2d& P, + const PolygonRing& outer_ring, + const std::vector& holes +); + +/** + * @brief Computes the winding number of a point with respect to a single closed ring. + * + * @param P The query point. + * @param ring A closed polygonal ring (at least 3 vertices). + * @return The winding number (positive for counter-clockwise turns). + */ +int compute_winding_number_ring(const Point2d& P, const PolygonRing& ring); + + +/** + * @brief Validates the polygon input for correct orientation and sufficient vertices. + * + * @param query_point The query point (for logging purposes). + * @param outer_ring The outer boundary ring. + * @param holes The list of hole rings. + * @param context Optional context string for logging. + * @return true if valid, false otherwise. + */ +bool validate_polygon_input( + const Point2d& query_point, + const PolygonRing& outer_ring, + const std::vector& holes, + const std::string& context = "" +); + + +/** + * @brief Computes the signed area of a polygon ring. + * + * @param ring The polygon ring (closed). + * @return The signed area (positive if CCW, negative if CW). + */ +double compute_ring_area(const PolygonRing& ring); \ No newline at end of file diff --git a/surface_integral/src/polygon_winding_number.cpp b/surface_integral/src/polygon_winding_number.cpp new file mode 100644 index 0000000..7e41073 --- /dev/null +++ b/surface_integral/src/polygon_winding_number.cpp @@ -0,0 +1,202 @@ +// polygon_winding.cpp +#include "polygon_winding_number.h" +#include + +const double eps = 1e-8; + +/** + * @brief Computes the contribution of a single edge (A->B) to the winding number. + * + * Uses a rightward horizontal ray from P. + * Only counts crossings where the edge passes **strictly upward** through the ray. + * This avoids double-counting at vertices. + * + * @param P The query point. + * @param A Start point of the edge. + * @param B End point of the edge. + * @return +1 if upward crossing, -1 if downward, 0 otherwise. + */ +int winding_contribution(const Point2d& P, const Point2d& A, const Point2d& B) { + // 1. If the edge is horizontal, no contribution + if (std::abs(A.y() - B.y()) < eps) { + return 0; + } + + // 2. Determine edge direction: does it cross P.y from below to above or vice versa? + bool a_below = (A.y() < P.y() - eps); + bool b_below = (B.y() < P.y() - eps); + + // 3. If both endpoints are on the same side of the ray, no crossing + if (a_below == b_below) { + return 0; + } + + // 4. Compute x-coordinate of intersection with y = P.y + double t = (P.y() - A.y()) / (B.y() - A.y()); + double x_intersect = A.x() + t * (B.x() - A.x()); + + // 5. Only consider intersections to the right of P (or directly on P) + if (x_intersect < P.x() - eps) { + return 0; + } + + // 6. Upward crossing: from below to above → +1 + if (a_below && !b_below) { + return +1; + } + // Downward crossing: from above to below → -1 + else if (!a_below && b_below) { + return -1; + } + + return 0; +} + +/** + * @brief Computes the total winding number of point P with respect to a closed ring. + */ +int compute_winding_number_ring(const Point2d& P, const PolygonRing& ring) { + if (ring.size() < 3) return 0; + int wn = 0; + int n = ring.size(); + for (int i = 0; i < n; ++i) { + const Point2d& A = ring[i]; + const Point2d& B = ring[(i + 1) % n]; + wn += winding_contribution(P, A, B); + } + return wn; +} + +// Helper: Check if P is on segment A-B +bool is_point_on_segment(const Point2d& P, const Point2d& A, const Point2d& B) { + if (P.x() < std::min(A.x(), B.x()) - eps || + P.x() > std::max(A.x(), B.x()) + eps || + P.y() < std::min(A.y(), B.y()) - eps || + P.y() > std::max(A.y(), B.y()) + eps) { + return false; + } + + double cross = (P.y() - A.y()) * (B.x() - A.x()) - (P.x() - A.x()) * (B.y() - A.y()); + if (std::abs(cross) > eps) return false; + + double dot = (P.x() - A.x()) * (B.x() - A.x()) + (P.y() - A.y()) * (B.y() - A.y()); + double len_sq = (B.x() - A.x()) * (B.x() - A.x()) + (B.y() - A.y()) * (B.y() - A.y()); + return (dot >= -eps && dot <= len_sq + eps); +} + +/** + * @brief Checks if a point is inside a polygon with holes using non-zero winding rule. + * + * - Points on outer ring → inside + * - Points on hole boundary → outside + * - Otherwise: wn(outer) + sum(wn(hole)) != 0 → inside + */ +bool point_in_polygon_with_holes( + const Point2d& P, + const PolygonRing& outer_ring, + const std::vector& holes) +{ + // Check boundary first + auto on_ring = [&](const PolygonRing& ring) { + int n = ring.size(); + for (int i = 0; i < n; ++i) { + if (is_point_on_segment(P, ring[i], ring[(i+1)%n])) { + return true; + } + } + return false; + }; + + if (on_ring(outer_ring)) return true; + for (const auto& hole : holes) { + if (on_ring(hole)) return false; + } + + // Compute total winding number + int total_wn = compute_winding_number_ring(P, outer_ring); + for (const auto& hole : holes) { + total_wn += compute_winding_number_ring(P, hole); // hole 是 CW，返回负值 + } + + return total_wn != 0; // 非零即内部 +} + + + +/** + * @brief Validates the polygon input for correct orientation and sufficient vertices. + * + * @param query_point The query point (for logging purposes). + * @param outer_ring The outer boundary ring. + * @param holes The list of hole rings. + * @param context Optional context string for logging. + * @return true if valid, false otherwise. + */ +bool validate_polygon_input( + const Point2d& query_point, + const PolygonRing& outer_ring, + const std::vector& holes, + const std::string& context) +{ + std::string label = context.empty() ? "Polygon" : context; + bool valid = true; + std::string issues = ""; + + // Check outer ring has enough points + if (outer_ring.size() < 3) { + issues += "Outer ring has fewer than 3 vertices."; + valid = false; + } else { + double outer_area = compute_ring_area(outer_ring); + if (outer_area <= 0) { + issues += "Outer ring is not CCW (area = " + std::to_string(outer_area) + ")"; + valid = false; + } + } + + // Check each hole + for (size_t i = 0; i < holes.size(); ++i) { + const auto& hole = holes[i]; + std::string hole_label = "Hole " + std::to_string(i); + + if (hole.size() < 3) { + if (!issues.empty()) issues += ", "; + issues += hole_label + " has fewer than 3 vertices"; + valid = false; + } else { + double hole_area = compute_ring_area(hole); + if (hole_area >= 0) { + if (!issues.empty()) issues += ", "; + issues += hole_label + " is not CW (area = " + std::to_string(hole_area) + ")"; + valid = false; + } + } + } + + // Output issues if invalid + if (!valid) { + std::cerr << "[\033[33mVALIDATION ERROR\033[0m] " << label + << " | Query: " << query_point + << " | Issues: " << issues << "\n"; + } + + return valid; +} + + +/** + * @brief Computes the signed area of a polygon ring. + * + * @param ring The polygon ring (closed). + * @return The signed area (positive if CCW, negative if CW). + */ +double compute_ring_area(const PolygonRing& ring) { + if (ring.size() < 3) return 0.0; + double area = 0.0; + int n = ring.size(); + for (int i = 0; i < n; ++i) { + int j = (i + 1) % n; + area += ring[i].x() * ring[j].y() - ring[j].x() * ring[i].y(); + } + return area * 0.5; +} \ No newline at end of file diff --git a/surface_integral/test/test_winding_number.cpp b/surface_integral/test/test_winding_number.cpp new file mode 100644 index 0000000..05be6c0 --- /dev/null +++ b/surface_integral/test/test_winding_number.cpp @@ -0,0 +1,265 @@ +// main.cpp +#include "polygon_winding_number.h" +#include +#include +#include +#include + +// ----------------------------- +// Test Case Definition +// ----------------------------- + +struct TestCase { + std::string name; + PolygonRing outer_ring; + std::vector holes; + Point2d query_point; + bool expected; // true = inside, false = outside +}; + +// ----------------------------- +// Helper: Print Point +// ----------------------------- +std::ostream& operator<<(std::ostream& os, const Point2d& p) { + return os << "(" << p.x() << ", " << p.y() << ")"; +} + + + +// ----------------------------- +// Helper: Validate Ring Orientation +// ----------------------------- +bool validate_orientation() { + // Test CCW (positive area) + PolygonRing ccw = {{0,0}, {1,0}, {1,1}, {0,1}}; + if (compute_ring_area(ccw) <= 0) { + std::cerr << "Error: CCW ring has non-positive area!\n"; + return false; + } + + // Test CW (negative area) + PolygonRing cw = {{0,0}, {0,1}, {1,1}, {1,0}}; + if (compute_ring_area(cw) >= 0) { + std::cerr << "Error: CW ring has non-negative area!\n"; + return false; + } + return true; +} + +// ----------------------------- +// Validate Test Case: Orientation of outer ring (CCW) and holes (CW) +// ----------------------------- +bool validate_test_case(const TestCase& tc) { + bool valid = true; + std::string status = ""; + + // Validate outer ring orientation (should be CCW → positive area) + double outer_area = compute_ring_area(tc.outer_ring); + if (outer_area <= 0) { + status += "Outer ring not CCW (area = " + std::to_string(outer_area) + ")"; + valid = false; + } + + // Validate holes orientation (should be CW → negative area) + for (size_t i = 0; i < tc.holes.size(); ++i) { + double hole_area = compute_ring_area(tc.holes[i]); + if (hole_area >= 0) { + if (!status.empty()) status += ", "; + status += "Hole " + std::to_string(i) + " not CW (area = " + std::to_string(hole_area) + ")"; + valid = false; + } + } + + if (valid) return true; + + std::cout << "[" + << "\033[33mINVALID\033[0m" // make "INVALID" yellow + << "] " + << tc.name << "\n" + << " Query: " << tc.query_point + << " | Status: \033[33mInvalid Test Case\033[0m" + << " | Reason: " << (valid ? "OK" : status) << "\n"; + + return valid; +} + +// ----------------------------- +// Run Single Test +// ----------------------------- +bool run_test(const TestCase& tc) { + bool result = point_in_polygon_with_holes(tc.query_point, tc.outer_ring, tc.holes); + bool passed = (result == tc.expected); + + std::cout << "[" + << (passed ? "\033[32mPASS\033[0m" : "\033[31mFAIL\033[0m") + << "] " + << tc.name << "\n" + << " Query: " << tc.query_point + << " | Expected: " << (tc.expected ? "Inside" : "Outside") + << " | Got: " << (result ? "Inside" : "Outside") << "\n"; + + return passed; +} + +// ----------------------------- +// Test Case Factory Functions +// ----------------------------- + +std::vector create_test_cases() { + std::vector tests; + + + + // -------------------------------------------------- + // Test 1: Simple hole - point inside hole → Outside + // -------------------------------------------------- + PolygonRing outer1 = {{0,0}, {2,0}, {2,2}, {0,2}}; + PolygonRing hole1 = {{0.5,0.5}, {0.5,1.5}, {1.5,1.5}, {1.5,0.5}}; + tests.push_back({ + "Inside Hole", + outer1, {hole1}, {1.0, 1.0}, false + }); + + // -------------------------------------------------- + // Test 2: Point inside outer, outside all holes → Inside + // -------------------------------------------------- + tests.push_back({ + "Inside Outer, Outside Hole", + outer1, {hole1}, {0.1, 0.1}, true + }); + + // -------------------------------------------------- + // Test 3: Point outside outer ring → Outside + // -------------------------------------------------- + tests.push_back({ + "Outside Outer", + outer1, {hole1}, {3.0, 3.0}, false + }); + + // -------------------------------------------------- + // Test 4: Point on outer ring edge (bottom edge) + // Winding rule: on boundary → considered inside + // -------------------------------------------------- + tests.push_back({ + "On Outer Edge (Bottom)", + outer1, {}, {1.0, 0.0}, true + }); + + // -------------------------------------------------- + // Test 5: Point on outer ring vertex + // -------------------------------------------------- + tests.push_back({ + "On Outer Vertex Offset Slightly Outside", + outer1, {}, {0.0, 0.0-1e-8}, false + }); + + // -------------------------------------------------- + // Test 6: Point on hole edge + // Should be outside (on hole boundary still "outside" the solid) + // -------------------------------------------------- + PolygonRing hole6 = {{1.0,1.0}, {1.0,2.0}, {2.0,2.0}, {2.0,1.0}}; // CW + tests.push_back({ + "On Hole Edge", + {{0,0},{3,0},{3,3},{0,3}}, {hole6}, {1.5, 1.0}, false + }); + + // -------------------------------------------------- + // Test 7: Multiple holes + // -------------------------------------------------- + PolygonRing hole7a = {{0.5,0.5}, {0.5,1.0}, {1.0,1.0}, {1.0,0.5}}; // CW + PolygonRing hole7b = {{1.5,1.5}, {1.5,2.0}, {2.0,2.0}, {2.0,1.5}}; // CW + tests.push_back({ + "Multiple Holes - Inside", + outer1, {hole7a, hole7b}, {0.1, 0.1}, true + }); + + tests.push_back({ + "Multiple Holes - In First Hole", + outer1, {hole7a, hole7b}, {0.7, 0.7}, false + }); + + tests.push_back({ + "Multiple Holes - In Second Hole", + outer1, {hole7a, hole7b}, {1.7, 1.7}, false + }); + + // -------------------------------------------------- + // Test 8: Nested holes (hole inside hole) - NOT standard, but test behavior + // Note: This is not valid in simple polygons, but test robustness + // -------------------------------------------------- + PolygonRing outer8 = {{0,0}, {3,0}, {3,3}, {0,3}}; + PolygonRing hole8a = {{0.5,0.5}, {0.5,2.5}, {2.5,2.5}, {2.5,0.5}}; // large hole + PolygonRing hole8b = {{1.0,1.0}, {1.0,2.0}, {2.0,2.0}, {2.0,1.0}}; // small hole inside hole8a + // This creates a "island" in the hole + tests.push_back({ + "Nested Holes - On Island", + outer8, {hole8a, hole8b}, {1.5, 1.5}, true // inside outer, inside hole8a, but outside hole8b → solid + }); + + // -------------------------------------------------- + // Test 9: Complex L-shaped polygon with hole + // -------------------------------------------------- + PolygonRing outer9 = { + {0,0}, {3,0}, {3,1}, {2,1}, {2,3}, {0,3}, {0,0} // L-shape, explicitly closed and no self-intersection + }; // CCW + PolygonRing hole9 = {{0.5,0.5}, {0.5,1.5}, {1.5,1.5}, {1.5,0.5}}; // CW + tests.push_back({ + "L-Shape with Hole - Inside Leg", + outer9, {hole9}, {0.1, 0.1}, true + }); + tests.push_back({ + "L-Shape with Hole - In Hole", + outer9, {hole9}, {1.0, 1.0}, false + }); + tests.push_back({ + "L-Shape with Hole - In Arm", + outer9, {hole9}, {2.0, 0.5}, true + }); + + // -------------------------------------------------- + // Test 10: Degenerate: Point exactly at intersection of hole and outer? (not possible, but test near) + // -------------------------------------------------- + tests.push_back({ + "Near Degenerate - Close to Corner", + outer1, {hole1}, {0.5000001, 0.5000001}, false // just outside hole + }); + + return tests; +} + +// ----------------------------- +// Main Function +// ----------------------------- +int main() { + std::cout << "Running Polygon Winding Number Tests...\n\n"; + + // Optional: Validate orientation logic + if (!validate_orientation()) { + std::cerr << "Orientation validation failed!\n"; + return 1; + } + + auto test_cases = create_test_cases(); + + int passed = 0; + int total = test_cases.size(); + + for (const auto& tc : test_cases) { + if (!validate_test_case(tc)) { + continue; + } + if (run_test(tc)) { + passed++; + } + } + + std::cout << "\nSummary: " << passed << " / " << total << " tests passed.\n"; + + if (passed == total) { + std::cout << "\033[32m All tests passed!\033[0m\n"; + return 0; + } else { + std::cout << "\033[31m Some tests failed!\033[0m\n"; + return 1; + } +} \ No newline at end of file diff --git a/surface_integral/xmake.lua b/surface_integral/xmake.lua index 55dbe58..576a476 100644 --- a/surface_integral/xmake.lua +++ b/surface_integral/xmake.lua @@ -20,3 +20,12 @@ target("mesh_test") add_deps("frontend") -- data structures from frontend, like polygon_mesh_t add_files("test/test_mesh.cpp") target_end() + +target("polygon_winding_number_test") + set_kind("binary") + add_packages("eigen") + add_includedirs("include") + add_files("src/polygon_winding_number.cpp") + add_files("test/test_winding_number.cpp") +target_end() +