Skip to content

Implement Polygon Boolean Operations (Union, Intersection, Difference)#693

Merged
yungyuc merged 1 commit intosolvcon:masterfrom
tigercosmos:polygon_3011
Mar 15, 2026
Merged

Implement Polygon Boolean Operations (Union, Intersection, Difference)#693
yungyuc merged 1 commit intosolvcon:masterfrom
tigercosmos:polygon_3011

Conversation

@tigercosmos
Copy link
Collaborator

@tigercosmos tigercosmos commented Mar 14, 2026

Implement Polygon Boolean Operations (Union, Intersection, Difference)

Summary

Implement polygon boolean operations using trapezoidal decomposition + Y-band sweep with X-interval merge. This approach is based on the algorithm described in:

"A New Method of Applying Polygon Boolean Operations Based on Trapezoidal Decomposition", GIScience & Remote Sensing, Vol 47, No 4 (2010)

The three operations — boolean_union, boolean_intersection, and boolean_difference — are implemented as a single unified algorithm parameterized by a boolean predicate.

Visualization

The visualization code is not included in this PR and is provided only for reference.

Order (left to right, top to down): original outlined polygonss, original colored polygons, the union reseult, the intersection result, the difference result.

image Screenshot 2026-03-14 at 2 56 57 PM Screenshot 2026-03-14 at 2 57 35 PM

Algorithm web explaination

Please use browser to read the web page:
https://gist.github.com/tigercosmos/4a358accfae620ada14a5ffa0776c668

Algorithm Overview

The algorithm has two phases:

Phase Description
1. Decomposition Each polygon is independently decomposed into horizontal-band trapezoids by the existing TrapezoidalDecomposer
2. Boolean Merge The two trapezoid sets are combined using a Y-band sweep with an X-interval sweep inside each band

Phase 2 in Detail

Step 1: Collect critical Y-values
         ├── From trapezoid boundaries (y_bottom, y_top of each trapezoid)
         └── From edge-edge crossings between polygon1 and polygon2 trapezoids

Step 2: For each Y-band [y_low, y_high]:
         ├── Gather X-interval events from spanning trapezoids
         ├── Sort events by X-position
         └── Sweep left-to-right, applying boolean predicate

Step 3: Emit result trapezoids as CCW polygons

Boolean Predicates

Operation Predicate Meaning
Union count_p1 > 0 OR count_p2 > 0 Inside either polygon
Intersection count_p1 > 0 AND count_p2 > 0 Inside both polygons
Difference count_p1 > 0 AND count_p2 == 0 Inside P1 but not P2

Worked Example

Two overlapping squares: P1 = (0,0)-(2,2), P2 = (1,1)-(3,3).

   y=3 - - - - - ┬────────────────┐- - - - - - - - - -
                 │          P2    │        band [2,3]
   y=2 ┬────────────────┐- - - - -│- - - - - - - - - -
       │   P1    │//////│         │        band [1,2]
   y=1 │- - - - -└──────┴─────────┘- - - - - - - - - -
       │                │                  band [0,1]
   y=0 └────────────────┘- - - - - - - - - - - - - - -
       x=0      x=1    x=2       x=3

Y-values: {0, 1, 2, 3} → 3 bands.

Band [1, 2] event sweep:

Events sorted by X:  P1-start@0  P2-start@1  P1-end@2  P2-end@3
                         │            │           │          │
count_p1:                1            1           0          0
count_p2:                0            1           1          0
X Event count_p1 count_p2 Intersection Union Difference
0 P1 start 1 0 - record left=0 record left=0
1 P2 start 1 1 record left=1 (still in) emit [0,1]
2 P1 end 0 1 emit [1,2] (still in) -
3 P2 end 0 0 - emit [0,3] -

Handling Slanted Edges

Each event stores two X-values — one at band bottom, one at band top — because polygon edges can be slanted. The emitted result is a general trapezoid with four corners:

        left_x_top ─────────── right_x_top      ← y_high
       /                              \
      /    result trapezoid            \
     /                                  \
    left_x_bottom ─────────── right_x_bottom    ← y_low

The bottom and top edges may have different widths, naturally handling triangles and trapezoids.

Edge-Crossing Y-Split

When trapezoid edges from different polygons cross within a Y-band, the crossing Y-value is detected and inserted into the critical Y-values set. This subdivides the band so that no two cross-polygon edges swap order within a single sub-band, ensuring correctness.

Before split:              After split:
┌──────────────┐           ┌──────────────┐
│  ╲   ╱       │           │   ╲  ╱       │ band [y_cross, y_top]
│    ╳         │    →      ├────╳─────────┤ ← y_cross inserted
│   ╱ ╲        │           │   ╱ ╲        │ band [y_bottom, y_cross]
│  ╱   ╲       │           │  ╱   ╲       │
└──────────────┘           └──────────────┘

Key Design Decisions

  • Reuses existing TrapezoidalDecomposer — no new decomposition code needed
  • Unified algorithm — all three operations share one function (compute_boolean_with_decomposition) parameterized by a BoolOp enum
  • No custom structs for intermediate data — works directly with TrapezoidPad indices and a source_of() lambda
  • Output as trapezoids — result polygons are emitted as CCW quadrilaterals, matching the existing trapezoidal representation
  • Relative epsilon — all floating-point comparisons use relative thresholds scaled by coordinate magnitude
  • Degenerate filtering — zero-area trapezoids from touching edges are skipped
  • Near-duplicate Y merge — Y-values within relative epsilon are merged to avoid degenerate bands

Files Changed

File Changes
cpp/modmesh/universe/polygon.hpp Added detail::BoolOp enum, detail::compute_boolean_with_decomposition() function, decomposed_trapezoids() accessor, implemented AreaBooleanUnion/Intersection/Difference::compute()
tests/test_polygon3d.py 13 boolean operation tests covering: overlapping squares, non-overlapping, triangles, containment, identical polygons, shared edges, commutativity, partial vertical overlap, concave polygons, CCW winding verification

Test Plan

  • All 52 tests pass for both float32 and float64
  • Overlapping squares: union=7, intersection=1, difference=3
  • Non-overlapping squares: union=2, intersection=0
  • Triangles touching at a point: intersection=0
  • Containment: union=16, intersection=4, difference=12
  • Identical polygons: union=P, intersection=P, difference=empty
  • Shared edge: union=2, intersection=0, difference=1
  • Commutativity: union(A,B)==union(B,A), intersection(A,B)==intersection(B,A)
  • Partial vertical overlap: union=7, intersection=1, difference=5
  • Concave (L-shaped) polygon: intersection=0.75, union=3.25, difference=2.25
  • All result polygons have CCW winding (non-negative signed area)

AI Usage

Used Claude Code with Opus 4.6 to implement the draft. Used Codex with 5.3-codex high effort for review. Fine-tuned the implementation manually with both Claude Code and Codex support. Used AI to generate all documents and the PR description.

@@ -0,0 +1,951 @@
<!DOCTYPE html>
Copy link
Collaborator Author

@tigercosmos tigercosmos Mar 14, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only for PR review purpose. Please use a broswer to open it. Will delete this file once the PR is approved.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please move the extended document in a gist. Make sure this file is not in your commit history.

* 5. Emitting result trapezoids as polygons
*/
template <typename T>
std::shared_ptr<PolygonPad<T>> compute_boolean_with_decomposition(
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Key function. Please check the PR description or ‎docs/polygon_boolean_algorithm.html‎ for the details.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function can be made to a helper class in the future, but the length is OK for now.


# TODO: Verify and implement boolean operations
# once the C++ side is complete.
def _compute_total_area(self, result_pad):
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed previous TODOs.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wonderful.

@tigercosmos
Copy link
Collaborator Author

@yungyuc Could you please take a look at this PR? The C++ implementation is about 400 LOC, but it mainly contains a single critical function and quite a few comment lines. The pytest code is around 300 LOC.

My preference would be to keep everything as is. However, if we need to strictly enforce the 500 LOC limit, I can help trim some comments and test cases.

@yungyuc yungyuc added the geometry Geometry entities and manipulation label Mar 14, 2026
Copy link
Member

@yungyuc yungyuc left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The tests look exciting! Points to address:

  • Move the extended document (‎docs/polygon_boolean_algorithm.html) in a gist.
  • Add all class end marker.
  • Do not use auto when the return type is not written in the RHS.
  • Use an end marker for long functions.

@@ -0,0 +1,951 @@
<!DOCTYPE html>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please move the extended document in a gist. Make sure this file is not in your commit history.

* Critical Y-values: {0, 1, 2, 3}
* Three y-bands: [0,1], [1,2], [2,3]
*
* y=3 - - - - - +----------------+- - - - - - - - - -
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you use a tool to generate the AA? If you did, please leave a comment stating what and how as a reference that may help other contributors.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I asked AI to leave explanation comments for the code. AI generated the AA directly. However, most of the time I need to manually adjust the chart for annotation or alignment.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hand-creation of AA usually took most of the efforts. I guess the manual adjustment did not take too much time.

Union,
Intersection,
Difference
};
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add class end marker.

namespace detail
{

enum class BoolOp
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I prefer BooleanOperation instead of BoolOp, but this is a must-do. Up to the author.

}

// Step 1: Decompose both polygons using the pad's TrapezoidalDecomposer
auto [begin1, end1] = pad->decompose_to_trapezoid(polygon_id1);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The return type is not written in the RHS. Do not use auto here.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use std::tie instead here.

// The source polygon is determined by the index range:
// [begin1, end1) -> polygon1 (source 0)
// [begin2, end2) -> polygon2 (source 1)
auto source_of = [begin1, end1, begin2, end2](size_t idx) -> int
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Many if not all of the lambdas should be made as a member function when this helper function is refactored as a helper class for better testability and performance.

};

// Collect all unique Y values from trapezoid boundaries
std::set<value_type> y_set;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Watch for potential performance hotspot with STL containers in the future.

}

return result;
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use an end marker for long functions.


# TODO: Verify and implement boolean operations
# once the C++ side is complete.
def _compute_total_area(self, result_pad):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wonderful.

result_diff = pad.boolean_difference(polygon1, polygon2)
self.assertEqual(result_diff.num_polygons, 0)

def test_boolean_shared_edge(self):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's good to test for edge touch. Point touch should be tested too, but it can be supplemented in the future.

@yungyuc
Copy link
Member

yungyuc commented Mar 14, 2026

@c1ydehhx please help review if you have time.


// Step 1: Decompose both polygons using the pad's TrapezoidalDecomposer
size_t begin1, end1, begin2, end2;
std::tie(begin1, end1) = pad->decompose_to_trapezoid(polygon_id1);
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed auto.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With the type the lines look clearer. It's a pity the types cannot be written in the same line of the function call.

size_t begin1, end1, begin2, end2;
std::tie(begin1, end1) = pad->decompose_to_trapezoid(polygon_id1);
std::tie(begin2, end2) = pad->decompose_to_trapezoid(polygon_id2);
std::shared_ptr<TrapezoidPad<T>> trap_pad = pad->decomposed_trapezoids();
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed auto.

}

return result;
} /* end function compute_boolean_with_decomposition */
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Mark.

value_type x_lo, x_hi; // X at y_low and y_high
int source;
bool is_start; // true = entering the polygon interval, false = leaving
}; /* end struct Event */
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Mark.

@tigercosmos
Copy link
Collaborator Author

@yungyuc The issues are addressed. Please take a look. Thank you.

  • Move the extended document (‎docs/polygon_boolean_algorithm.html) in a gist. -> Put the gist in the PR description.
  • Add all class end marker. -> Done
  • Do not use auto when the return type is not written in the RHS. -> Done
  • Use an end marker for long functions. -> Done

Copy link
Member

@yungyuc yungyuc left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Additional/remaining points to address:

  • Leave a code comment to reminder to profile for the closure find_crossing in the nested loop. It looks like a performance hotspot.
  • Should be end enum class BooleanOperation, not end enum BooleanOperation.

* Critical Y-values: {0, 1, 2, 3}
* Three y-bands: [0,1], [1,2], [2,3]
*
* y=3 - - - - - +----------------+- - - - - - - - - -
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hand-creation of AA usually took most of the efforts. I guess the manual adjustment did not take too much time.


// Step 1: Decompose both polygons using the pad's TrapezoidalDecomposer
size_t begin1, end1, begin2, end2;
std::tie(begin1, end1) = pad->decompose_to_trapezoid(polygon_id1);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With the type the lines look clearer. It's a pity the types cannot be written in the same line of the function call.

// Find the Y-value where two edges cross within the overlap range.
// Each edge is defined by its X-values at y_lo and y_hi.
// clang-format off
auto find_crossing = [&](value_type edge_a_x_at_bottom, value_type edge_a_x_at_top,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Leave a comment to remind us to profile later and use the timing results to refactor. An issue should be created to track it after the PR is merged.

I am not totally sure about the performance of a closure in the nested loop. My gut feeling is that this is not gonna be fast. It may be hard to profile too.

For the sake of clarity, the nested loop may use a distinct helper class. But it is for future work.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not seeing the comment (code remark) to remind to profile in the future. Please point it out using inline review annotation.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The beginning of the lambda is the primary location to leave a TODO comment for potential performance hotspot.

Union,
Intersection,
Difference
}; /* end enum BooleanOperation */
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should be end enum class BooleanOperation, not end enum BooleanOperation.

Use exact declaration will help grepping.

@tigercosmos
Copy link
Collaborator Author

@yungyuc should be good now. :D

Copy link
Member

@yungyuc yungyuc left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A previous point was not addressed. Please clarify.

// Find the Y-value where two edges cross within the overlap range.
// Each edge is defined by its X-values at y_lo and y_hi.
// clang-format off
auto find_crossing = [&](value_type edge_a_x_at_bottom, value_type edge_a_x_at_top,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not seeing the comment (code remark) to remind to profile in the future. Please point it out using inline review annotation.

return bottom_point.x() + t * (top_point.x() - bottom_point.x());
};

// TODO: potential performance issue: std container
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

TODO

value_type trap_j_right_at_bottom = lerp_x(trap_j_right_bottom, trap_j_right_top, y_lo);
value_type trap_j_right_at_top = lerp_x(trap_j_right_bottom, trap_j_right_top, y_hi);

// TODO: potential performance issue: closure in the nested loop
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

TODO

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not the best place for the TODO comment, but it makes sense to have one note here too. Don't delete it.

@tigercosmos
Copy link
Collaborator Author

@yungyuc I have specified where the TODOs are. Please check the comments above. Thanks!

Copy link
Member

@yungyuc yungyuc left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@yungyuc I have specified where the TODOs are. Please check the comments above. Thanks!

Thank you, now I see them, but one is not at a proper location. If I have time I will modify it myself, but you can do it too. Let's see who's gonna be faster.

value_type trap_j_right_at_bottom = lerp_x(trap_j_right_bottom, trap_j_right_top, y_lo);
value_type trap_j_right_at_top = lerp_x(trap_j_right_bottom, trap_j_right_top, y_hi);

// TODO: potential performance issue: closure in the nested loop
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not the best place for the TODO comment, but it makes sense to have one note here too. Don't delete it.

// Find the Y-value where two edges cross within the overlap range.
// Each edge is defined by its X-values at y_lo and y_hi.
// clang-format off
auto find_crossing = [&](value_type edge_a_x_at_bottom, value_type edge_a_x_at_top,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The beginning of the lambda is the primary location to leave a TODO comment for potential performance hotspot.

Replace the three TODO stubs (union, intersection, difference) with a
working implementation that decomposes both input polygons into
trapezoids, collects critical Y-values including cross-polygon edge
crossings, then sweeps X-intervals per band to emit result trapezoids.

Add comprehensive tests covering overlapping squares, containment,
identical polygons, shared edges, commutativity, concave polygons,
and CCW winding order verification.
@tigercosmos
Copy link
Collaborator Author

@yungyuc I am faster. 🥇

Copy link
Member

@yungyuc yungyuc left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@yungyuc I am faster. 🥇

Glad to lose. Will merge once CI finishes.

@yungyuc yungyuc merged commit 3ea2820 into solvcon:master Mar 15, 2026
14 checks passed
@yungyuc
Copy link
Member

yungyuc commented Mar 15, 2026

@tigercosmos Issue #694 is created for you.

@tigercosmos tigercosmos deleted the polygon_3011 branch March 15, 2026 12:21
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

geometry Geometry entities and manipulation

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants