@@ -1036,6 +1036,18 @@ impl Vector {
10361036 false
10371037 }
10381038
1039+ /// Compute signed area of a face using the shoelace formula.
1040+ /// Negative area indicates clockwise winding, positive indicates counterclockwise.
1041+ fn compute_signed_area ( & self , face : & FoundSubpath ) -> f64 {
1042+ let mut area = 0.0 ;
1043+ for edge in & face. edges {
1044+ let start_pos = self . point_domain . positions ( ) [ edge. start ] ;
1045+ let end_pos = self . point_domain . positions ( ) [ edge. end ] ;
1046+ area += ( end_pos. x - start_pos. x ) * ( end_pos. y + start_pos. y ) ;
1047+ }
1048+ area * 0.5
1049+ }
1050+
10391051 /// Find all minimal closed regions (faces) in a branching mesh vector.
10401052 pub fn find_closed_regions ( & self ) -> Vec < FoundSubpath > {
10411053 let mut regions = Vec :: new ( ) ;
@@ -1122,6 +1134,18 @@ impl Vector {
11221134 }
11231135 }
11241136
1137+ // Filter out outer (counterclockwise) faces, keeping only inner (clockwise) faces
1138+ regions. retain ( |face| {
1139+ // Always include 2-edge faces (floating point issues with area calculation)
1140+ if face. edges . len ( ) == 2 {
1141+ return true ;
1142+ }
1143+
1144+ let area = self . compute_signed_area ( face) ;
1145+ // Keep clockwise faces (negative area), exclude counterclockwise (positive area)
1146+ area < 0.0
1147+ } ) ;
1148+
11251149 regions
11261150 }
11271151
@@ -1140,62 +1164,45 @@ impl Vector {
11401164 let target = from_point;
11411165 let mut prev_segment = start_segment;
11421166
1143- let mut iteration = 0 ;
11441167 let max_iterations = adjacency. len ( ) * 2 ;
11451168
11461169 // Follow the "rightmost" edge at each vertex to find minimal face
1147- loop {
1148- iteration += 1 ;
1170+ for _iteration in 1 ..=max_iterations {
1171+ let neighbors = adjacency . get ( & current ) ? ;
11491172
1150- if iteration > max_iterations {
1151- return None ;
1152- }
1173+ // Since neighbors are pre-sorted by angle, find the index of the edge we came from
1174+ // and take the next edge in the sorted circular list
1175+ let prev_index = neighbors. iter ( )
1176+ . position ( |( seg, _next, _rev) | * seg == prev_segment) ?;
11531177
1154- let neighbors = adjacency . get ( & current ) ? ;
1155- // Find the next edge in counterclockwise order (rightmost turn)
1156- let prev_direction = self . point_domain . positions ( ) [ current ] - self . point_domain . positions ( ) [ path . last ( ) ? . start ] ;
1178+ // The next edge in the sorted list is the minimal face edge (rightmost turn)
1179+ // Wrap around using modulo to handle circular list
1180+ let next = neighbors [ ( prev_index + 1 ) % neighbors . len ( ) ] ;
11571181
1158- let angle_between = |v1 : DVec2 , v2 : DVec2 | -> f64 {
1159- let angle = v2. y . atan2 ( v2. x ) - v1. y . atan2 ( v1. x ) ;
1160- if angle < 0.0 { angle + 2.0 * std:: f64:: consts:: PI } else { angle }
1161- } ;
1182+ let ( next_segment, next_point, next_reversed) = next;
11621183
1163- let next = neighbors. iter ( ) . filter ( |( seg, _next, _rev) | * seg != prev_segment) . min_by ( |a, b| {
1164- let dir_a = self . point_domain . positions ( ) [ a. 1 ] - self . point_domain . positions ( ) [ current] ;
1165- let dir_b = self . point_domain . positions ( ) [ b. 1 ] - self . point_domain . positions ( ) [ current] ;
1166- let angle_a = angle_between ( prev_direction, dir_a) ;
1167- let angle_b = angle_between ( prev_direction, dir_b) ;
1168- angle_a. partial_cmp ( & angle_b) . unwrap_or ( std:: cmp:: Ordering :: Equal )
1169- } ) ?;
1184+ // Check if we've created a cycle (might not be back to start)
1185+ if path. iter ( ) . any ( |e| e. end == next_point && e. id != next_segment) {
1186+ return None ;
1187+ }
11701188
1171- let ( next_segment, next_point, next_reversed) = * next ;
1189+ path . push ( HalfEdge :: new ( next_segment, current , next_point, next_reversed) ) ;
11721190
1191+ // Check if we've completed the face
11731192 if next_point == target {
1174- // Completed the cycle
1175- path. push ( HalfEdge :: new ( next_segment, current, next_point, next_reversed) ) ;
1176-
11771193 // Mark all half-edges as used
11781194 for edge in & path {
11791195 used_half_edges. insert ( ( edge. id , edge. reverse ) ) ;
11801196 }
1181-
11821197 return Some ( FoundSubpath :: new ( path) ) ;
11831198 }
11841199
1185- // Check if we've created a cycle (might not be back to start)
1186- if path. iter ( ) . any ( |e| e. end == next_point && e. id != next_segment) {
1187- return None ;
1188- }
1189-
1190- path. push ( HalfEdge :: new ( next_segment, current, next_point, next_reversed) ) ;
11911200 prev_segment = next_segment;
11921201 current = next_point;
1193-
1194- // Prevent infinite loops
1195- if path. len ( ) > adjacency. len ( ) {
1196- return None ;
1197- }
11981202 }
1203+
1204+ // If we exit the loop without returning, the path didn't close
1205+ None
11991206 }
12001207}
12011208
0 commit comments