Optionally ignore direction. #9

Implement by forcing orientation of the entire linestring. Is this
brittle?

rust/cli/src/main.rs

@@ -5,7 +5,9 @@ use clap::{Arg, ArgAction, Command};
 use geo::GeodesicLength;
 use geojson::{Feature, FeatureCollection, GeoJson};
 
-use overline::{aggregate_properties, feature_to_line_string, overline, Aggregation, Output};
+use overline::{
+    aggregate_properties, feature_to_line_string, overline, Aggregation, Options, Output,
+};
 
 fn main() -> Result<()> {
     let mut args = Command::new("overline")
@@ -24,6 +26,9 @@ fn main() -> Result<()> {
     let input_path = args.remove_one::<String>("FILE").unwrap();
     let output_path = args.remove_one::<String>("output").unwrap();
 
+    // TODO Add a flag
+    let options = Options::default();
+
     println!("Reading and deserializing {input_path}");
     let mut now = Instant::now();
     let geojson: GeoJson = std::fs::read_to_string(input_path)?.parse()?;
@@ -36,7 +41,7 @@ fn main() -> Result<()> {
 
     println!("Running overline on {} line-strings", input.len());
     now = Instant::now();
-    let output = overline(&input);
+    let output = overline(&input, options);
     println!("... took {:?}", now.elapsed());
 
     if let Some(sum_property) = args.get_one::<String>("summary") {

rust/overline/src/lib.rs

@@ -1,5 +1,6 @@
 use std::collections::HashMap;
 
+use geo::{HaversineBearing, Winding};
 use geojson::Feature;
 use ordered_float::NotNan;
 use rayon::prelude::*;
@@ -20,33 +21,61 @@ pub struct Output {
     pub indices: Vec<usize>,
 }
 
+#[derive(Default)]
+pub struct Options {
+    /// TODO describe. may reverse input geometry.
+    pub ignore_direction: bool,
+}
+
 /// Ignores anything aside from LineStrings. Returns LineStrings chopped up to remove overlap, with
 /// exactly one property -- `indices`, an array of numbers indexing the input that share that
 /// geometry.
-pub fn overline(input: &Vec<Feature>) -> Vec<Output> {
+pub fn overline(input: &Vec<Feature>, options: Options) -> Vec<Output> {
+    // Extract LineStrings from input
+    let input_linestrings: Vec<Option<geo::LineString<f64>>> = input
+        .par_iter()
+        .map(|f| {
+            feature_to_line_string(f).map(|mut linestring| {
+                if options.ignore_direction {
+                    // TODO Do we need to project to euclidean first?
+                    println!("input is cw? {:?}", linestring.winding_order());
+                    linestring.make_cw_winding();
+
+                    let pt1 = linestring.0[0];
+                    let pt2 = *linestring.0.last().unwrap();
+                    let bearing = geo::Point::from(pt1).haversine_bearing(geo::Point::from(pt2));
+                    if bearing < 0.0 {
+                        linestring.0.reverse();
+                    }
+                }
+                linestring
+            })
+        })
+        .collect();
+
     // For every individual (directed) line segment, record the index of inputs matching there
     let mut line_segments: HashMap<(HashedPoint, HashedPoint), Vec<usize>> = HashMap::new();
-    for (idx, input) in input.iter().enumerate() {
-        if let Some(geom) = feature_to_line_string(input) {
-            for line in geom.lines().map(HashedPoint::new_line) {
+    for (idx, maybe_linestring) in input_linestrings.iter().enumerate() {
+        if let Some(ref linestring) = maybe_linestring {
+            for line in linestring.lines().map(HashedPoint::new_line) {
                 line_segments.entry(line).or_insert_with(Vec::new).push(idx);
             }
         }
     }
     // Then look at each input, accumulating points as long all the indices match
-    input
+    input_linestrings
         .par_iter()
         .enumerate()
-        .flat_map(|(idx, input_feature)| {
+        .flat_map(|(idx, maybe_linestring)| {
             let mut intermediate_output = Vec::new();
 
-            // This state is reset as we look through this input's points
-            let mut pts_so_far = Vec::new();
-            let mut indices_so_far = Vec::new();
-            let mut keep_this_output = false;
+            if let Some(ref linestring) = maybe_linestring {
+                // This state is reset as we look through this input's points
+                let mut pts_so_far = Vec::new();
+                let mut indices_so_far = Vec::new();
+                let mut keep_this_output = false;
 
-            if let Some(geom) = feature_to_line_string(input_feature) {
-                for line in geom.lines() {
+                for line in linestring.lines() {
                     // The segment is guaranteed to exist
                     let indices = &line_segments[&HashedPoint::new_line(line)];
 
@@ -76,13 +105,14 @@ pub fn overline(input: &Vec<Feature>) -> Vec<Output> {
                     // below, since we process input in order.
                     keep_this_output = indices_so_far.iter().all(|i| *i >= idx);
                 }
-            }
-            // This input ended; add to output if needed
-            if !pts_so_far.is_empty() && keep_this_output {
-                intermediate_output.push(Output {
-                    geometry: pts_so_far.into(),
-                    indices: indices_so_far,
-                });
+
+                // This input ended; add to output if needed
+                if !pts_so_far.is_empty() && keep_this_output {
+                    intermediate_output.push(Output {
+                        geometry: pts_so_far.into(),
+                        indices: indices_so_far,
+                    });
+                }
             }
 
             intermediate_output