integration-of-complex-surfaces
/
ImplicitSurfaceNetwork
mirror of http://118.195.195.192:3000/ZCWang/ImplicitSurfaceNetwork.git
2
0
Fork 1
Code Issues Projects Releases Wiki Activity
Browse Source

feat: add 2D winding number to determine whether a point inside or outside a p;olygon

feat-integrator
mckay 2 weeks ago
parent
0aee7c297a
commit
c04fba163b
4 changed files with 536 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 60
      surface_integral/include/polygon_winding_number.h
  2. 202
      surface_integral/src/polygon_winding_number.cpp
  3. 265
      surface_integral/test/test_winding_number.cpp
  4. 9
      surface_integral/xmake.lua

60
surface_integral/include/polygon_winding_number.h
View File

@ -0,0 +1,60 @@
// polygon_winding.h
#pragma once
#include <vector>
#include <Eigen/Dense>
#include <string>
#include <iostream>
// Type aliases for convenience
using Point2d = Eigen::Vector2d;
using PolygonRing = std::vector<Point2d>;
/**
* @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<PolygonRing>& 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<PolygonRing>& 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);

202
surface_integral/src/polygon_winding_number.cpp
View File

@ -0,0 +1,202 @@
// polygon_winding.cpp
#include "polygon_winding_number.h"
#include <cmath>
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<PolygonRing>& 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<PolygonRing>& 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;
}

265
surface_integral/test/test_winding_number.cpp
View File

@ -0,0 +1,265 @@
// main.cpp
#include "polygon_winding_number.h"
#include <iostream>
#include <vector>
#include <string>
#include <cassert>
// -----------------------------
// Test Case Definition
// -----------------------------
struct TestCase {
std::string name;
PolygonRing outer_ring;
std::vector<PolygonRing> 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<TestCase> create_test_cases() {
std::vector<TestCase> 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;
}
}

9
surface_integral/xmake.lua
View File

@ -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()

Write Preview
Loading…
Cancel
Save