Merge pull request #10 from bluenote10/feature/refactor_connect_edges

Refactor connect edges

.gitignore

@@ -5,4 +5,6 @@
 .idea/
 .vscode/
 
+# Files generated from debugging
 test_cases.pdf
+*.geojson.generated

lib/src/boolean/compare_segments.rs

@@ -1,6 +1,6 @@
 use super::signed_area::signed_area;
 use super::sweep_event::SweepEvent;
-use num_traits::Float;
+use super::helper::Float;
 use std::cmp::Ordering;
 use std::rc::Rc;
 

lib/src/boolean/compute_fields.rs

@@ -1,6 +1,6 @@
-use super::sweep_event::{EdgeType, SweepEvent};
+use super::sweep_event::{EdgeType, SweepEvent, ResultTransition};
 use super::Operation;
-use num_traits::Float;
+use super::helper::Float;
 use std::rc::Rc;
 
 pub fn compute_fields<F>(event: &Rc<SweepEvent<F>>, maybe_prev: Option<&Rc<SweepEvent<F>>>, operation: Operation)
@@ -15,11 +15,28 @@ where
         } else {
             event.set_in_out(!prev.is_other_in_out(), prev.is_in_out());
         }
+
+        // Connect to previous in result: Only use the given `prev` if it is
+        // part of the result and not a vertical segment. Otherwise connect
+        // to its previous in result if any.
+        if prev.is_in_result() && !prev.is_vertical() {
+            event.set_prev_in_result(prev);
+        } else if let Some(prev_of_prev) = prev.get_prev_in_result() {
+            event.set_prev_in_result(&prev_of_prev);
+        }
     } else {
         event.set_in_out(false, true);
     }
 
-    event.set_in_result(in_result(event, operation));
+    // Determine whether segment is in result, and if so, whether it is an
+    // in-out or out-in transition.
+    let in_result = in_result(event, operation);
+    let result_transition = if !in_result {
+        ResultTransition::None
+    } else {
+        determine_result_transition(&event, operation)
+    };
+    event.set_result_transition(result_transition);
 }
 
 fn in_result<F>(event: &SweepEvent<F>, operation: Operation) -> bool
@@ -40,3 +57,24 @@ where
         EdgeType::NonContributing => false,
     }
 }
+
+fn determine_result_transition<F>(event: &SweepEvent<F>, operation: Operation) -> ResultTransition
+where
+    F: Float,
+{
+    let this_in = !event.is_in_out();
+    let that_in = !event.is_other_in_out();
+    let is_in = match operation {
+        Operation::Intersection => this_in && that_in,
+        Operation::Union        => this_in || that_in,
+        Operation::Xor          => this_in ^  that_in,
+        Operation::Difference   =>
+            // Difference is assymmetric, so subject vs clipping matters.
+            if event.is_subject {
+                this_in && !that_in
+            } else {
+                that_in && !this_in
+            }
+    };
+    if is_in { ResultTransition::OutIn } else { ResultTransition::InOut }
+}

lib/src/boolean/connect_edges.rs

@@ -1,7 +1,6 @@
-use super::sweep_event::SweepEvent;
-use super::Operation;
-use geo_types::{LineString, Polygon};
-use num_traits::Float;
+use super::sweep_event::{SweepEvent, ResultTransition};
+use geo_types::{Coordinate};
+use super::helper::Float;
 use std::collections::HashSet;
 use std::rc::Rc;
 
@@ -30,24 +29,25 @@ where
         }
     }
 
+    // Populate `other_pos` by initializing with index and swapping with other event.
     for (pos, event) in result_events.iter().enumerate() {
-        event.set_pos(pos as i32)
+        event.set_other_pos(pos as i32)
     }
-
     for event in &result_events {
-        if !event.is_left() {
+        if event.is_left() {
             if let Some(other) = event.get_other_event() {
-                let tmp = event.get_pos();
-                event.set_pos(other.get_pos());
-                other.set_pos(tmp);
+                let (a, b) = (event.get_other_pos(), other.get_other_pos());
+                event.set_other_pos(b);
+                other.set_other_pos(a);
             }
         }
     }
 
     result_events
 }
 
-fn next_pos<F>(pos: i32, result_events: &[Rc<SweepEvent<F>>], processed: &mut HashSet<i32>, orig_index: i32) -> i32
+
+fn next_pos<F>(pos: i32, result_events: &[Rc<SweepEvent<F>>], processed: &HashSet<i32>, orig_pos: i32) -> i32
 where
     F: Float,
 {
@@ -73,58 +73,157 @@ where
 
     new_pos = pos - 1;
 
-    while processed.contains(&new_pos) && new_pos >= orig_index as i32 {
+    while processed.contains(&new_pos) && new_pos > orig_pos {
         new_pos -= 1;
     }
     new_pos
 }
 
-pub fn connect_edges<F>(sorted_events: &[Rc<SweepEvent<F>>], operation: Operation) -> Vec<Polygon<F>>
+
+pub struct Contour<F>
+where
+    F: Float
+{
+    /// Raw coordinates of contour
+    pub points: Vec<Coordinate<F>>,
+    /// Contour IDs of holes if any.
+    pub hole_ids: Vec<i32>,
+    /// Contour ID of parent if this contour is a hole.
+    pub hole_of: Option<i32>,
+    /// Depth of the contour. Since the geo data structures don't store depth information,
+    /// this field is not strictly necessary to compute. But it is very cheap to compute,
+    /// so we can add it and see if it has relevance in the future.
+    pub depth: i32,
+}
+
+impl<F> Contour<F>
+where
+    F: Float
+{
+    pub fn new(hole_of: Option<i32>, depth: i32) -> Contour<F> {
+        Contour {
+            points: Vec::new(),
+            hole_ids: Vec::new(),
+            hole_of: hole_of,
+            depth: depth,
+        }
+    }
+
+    /// This logic implements the 4 cases of parent contours from Fig. 4 in the Martinez paper.
+    pub fn initialize_from_context(event: &Rc<SweepEvent<F>>, contours: &mut [Contour<F>], contour_id: i32) -> Contour<F> {
+        if let Some(prev_in_result) = event.get_prev_in_result() {
+            // Note that it is valid to query the "previous in result" for its output contour id,
+            // because we must have already processed it (i.e., assigned an output contour id)
+            // in an earlier iteration, otherwise it wouldn't be possible that it is "previous in
+            // result".
+            let lower_contour_id = prev_in_result.get_output_contour_id();
+            if prev_in_result.get_result_transition() == ResultTransition::OutIn {
+                // We are inside. Now we have to check if the thing below us is another hole or
+                // an exterior contour.
+                let lower_contour = &contours[lower_contour_id as usize];
+                if let Some(parent_contour_id) = lower_contour.hole_of {
+                    // The lower contour is a hole => Connect the new contour as a hole to its parent,
+                    // and use same depth.
+                    contours[parent_contour_id as usize].hole_ids.push(contour_id);
+                    let hole_of = Some(parent_contour_id);
+                    let depth = contours[lower_contour_id as usize].depth;
+                    Contour::new(hole_of, depth)
+                } else {
+                    // The lower contour is an exterior contour => Connect the new contour as a hole,
+                    // and increment depth.
+                    contours[lower_contour_id as usize].hole_ids.push(contour_id);
+                    let hole_of = Some(lower_contour_id);
+                    let depth = contours[lower_contour_id as usize].depth + 1;
+                    Contour::new(hole_of, depth)
+                }
+            } else {
+                // We are outside => this contour is an exterior contour of same depth.
+                let depth = contours[lower_contour_id as usize].depth;
+                Contour::new(None, depth)
+            }
+        } else {
+            // There is no lower/previous contour => this contour is an exterior contour of depth 0.
+            Contour::new(None, 0)
+        }
+    }
+
+    /// Whether a contour is an exterior contour or a hole.
+    /// Note: The semantics of `is_exterior` are in the sense of an exterior ring of a
+    /// polygon in GeoJSON, not to be confused with "external contour" as used in the
+    /// Martinez paper (which refers to contours that are not included in any of the
+    /// other polygon contours; `is_exterior` is true for all outer contours, not just
+    /// the outermost).
+    pub fn is_exterior(&self) -> bool {
+        self.hole_of.is_none()
+    }
+}
+
+
+fn mark_as_processed<F>(processed: &mut HashSet<i32>, result_events: &[Rc<SweepEvent<F>>], pos: i32, contour_id: i32)
+where
+    F: Float,
+{
+    processed.insert(pos);
+    result_events[pos as usize].set_output_contour_id(contour_id);
+}
+
+
+pub fn connect_edges<F>(sorted_events: &[Rc<SweepEvent<F>>]) -> Vec<Contour<F>>
 where
     F: Float,
 {
     let result_events = order_events(sorted_events);
 
-    let mut result: Vec<Polygon<F>> = Vec::new();
+    let mut contours: Vec<Contour<F>> = Vec::new();
     let mut processed: HashSet<i32> = HashSet::new();
 
     for i in 0..(result_events.len() as i32) {
         if processed.contains(&i) {
             continue;
         }
-        let mut contour = LineString::<F>(Vec::new());
+
+        let contour_id = contours.len() as i32;
+        let mut contour = Contour::initialize_from_context(
+            &result_events[i as usize],
+            &mut contours,
+            contour_id,
+        );
+
+        let orig_pos = i; // Alias just for clarity
         let mut pos = i;
-        let initial = result_events[i as usize].point;
 
-        contour.0.push(initial);
+        let initial = result_events[pos as usize].point;
+        contour.points.push(initial);
 
-        while pos >= i {
-            processed.insert(pos);
+        loop {
+            // Loop clarifications:
+            // - An iteration has two kinds of `pos` advancements:
+            //   (A) following a segment via `other_pos`, and
+            //   (B) searching for the next outgoing edge on same point.
+            // - Therefore, the loop contains two "mark pos as processed" steps, using the
+            //   convention that at beginning of the loop, `pos` isn't marked yet.
+            // - The contour is extended after following a segment.
+            // - Hitting pos == orig_pos after search (B) indicates no continuation and
+            //   terminates the loop.
+            mark_as_processed(&mut processed, &result_events, pos, contour_id);
 
-            pos = result_events[pos as usize].get_pos();
-            processed.insert(pos);
-            contour.0.push(result_events[pos as usize].point);
-            pos = next_pos(pos, &result_events, &mut processed, i);
-        }
+            pos = result_events[pos as usize].get_other_pos();              // pos advancement (A)
 
-        if !result_events[i as usize].is_exterior_ring {
-            if result.is_empty() {
-                result.push(Polygon::new(contour, Vec::new()));
-            } else {
-                result
-                    .last_mut()
-                    .expect("Result must not be empty at this point")
-                    .interiors_push(contour);
+            mark_as_processed(&mut processed, &result_events, pos, contour_id);
+            contour.points.push(result_events[pos as usize].point);
+
+            pos = next_pos(pos, &result_events, &processed, orig_pos);      // pos advancement (B)
+
+            if pos == orig_pos {
+                break;
             }
-        } else if operation == Operation::Difference && !result_events[i as usize].is_subject && result.len() > 1 {
-            result
-                .last_mut()
-                .expect("Result must not be empty at this point")
-                .interiors_push(contour);
-        } else {
-            result.push(Polygon::new(contour, Vec::new()));
         }
+
+        // This assert should be possible once the first stage of the algorithm is robust.
+        // debug_assert_eq!(contour.points.first(), contour.points.last());
+
+        contours.push(contour);
     }
 
-    result
+    contours
 }

lib/src/boolean/divide_segment.rs

@@ -1,6 +1,6 @@
 use super::sweep_event::SweepEvent;
 use geo_types::Coordinate;
-use num_traits::Float;
+use super::helper::Float;
 use std::collections::BinaryHeap;
 use std::rc::Rc;
 

lib/src/boolean/fill_queue.rs

@@ -1,5 +1,5 @@
 use geo_types::{LineString, Polygon, Rect};
-use num_traits::Float;
+use super::helper::Float;
 use std::collections::BinaryHeap;
 use std::rc::{Rc, Weak};
 

lib/src/boolean/helper.rs

@@ -1,4 +1,9 @@
 use std::cmp::Ordering;
+use std::fmt::{Debug, Display};
+use num_traits::Float as NumTraitsFloat;
+
+pub trait Float: NumTraitsFloat + Debug + Display {}
+impl<T: NumTraitsFloat + Debug + Display> Float for T {}
 
 #[inline]
 pub fn less_if(condition: bool) -> Ordering {

lib/src/boolean/mod.rs

@@ -1,6 +1,4 @@
-use num_traits::Float;
-
-use geo_types::{Coordinate, MultiPolygon, Polygon, Rect};
+use geo_types::{Coordinate, MultiPolygon, Polygon, LineString, Rect};
 
 pub mod compare_segments;
 pub mod compute_fields;
@@ -14,6 +12,8 @@ mod signed_area;
 pub mod subdivide_segments;
 pub mod sweep_event;
 
+pub use helper::Float;
+
 use self::connect_edges::connect_edges;
 use self::fill_queue::fill_queue;
 use self::subdivide_segments::subdivide;
@@ -110,7 +110,23 @@ where
 
     let sorted_events = subdivide(&mut event_queue, &sbbox, &cbbox, operation);
 
-    MultiPolygon(connect_edges(&sorted_events, operation))
+    let contours = connect_edges(&sorted_events);
+
+    // Convert contours into polygons
+    let polygons: Vec<Polygon<F>> = contours
+        .iter()
+        .filter(|contour| contour.is_exterior())
+        .map(|contour| {
+            let exterior = LineString(contour.points.clone());
+            let mut interios: Vec<LineString<F>> = Vec::new();
+            for hole_id in &contour.hole_ids {
+                interios.push(LineString(contours[*hole_id as usize].points.clone()));
+            }
+            Polygon::new(exterior, interios)
+        })
+        .collect();
+
+    MultiPolygon(polygons)
 }
 
 fn trivial_result<F>(subject: &[Polygon<F>], clipping: &[Polygon<F>], operation: Operation) -> MultiPolygon<F>

lib/src/boolean/possible_intersection.rs

@@ -1,7 +1,7 @@
 use super::divide_segment::divide_segment;
 use super::segment_intersection::{intersection, LineIntersection};
 use super::sweep_event::{EdgeType, SweepEvent};
-use num_traits::Float;
+use super::helper::Float;
 use std::collections::BinaryHeap;
 use std::rc::Rc;
 

lib/src/boolean/segment_intersection.rs

@@ -1,5 +1,5 @@
 use geo_types::Coordinate;
-use num_traits::Float;
+use super::helper::Float;
 
 #[derive(Debug, Clone, Copy, PartialEq)]
 pub enum LineIntersection<F>