fix(layout): improve calculation of band positions for Seg's

Cargo.toml

@@ -52,6 +52,7 @@ contracts-try = "0.7"
 derive-getters.workspace = true
 enum_dispatch = "0.3"
 geo.workspace = true
+log.workspace = true
 petgraph.workspace = true
 rstar.workspace = true
 serde.workspace = true

src/geometry/primitive.rs

@@ -6,7 +6,7 @@ use std::f64::consts::TAU;
 
 use enum_dispatch::enum_dispatch;
 use geo::algorithm::line_measures::{Distance, Euclidean};
-use geo::{point, polygon, Contains, Intersects, Point, Polygon, Rotate};
+use geo::{point, polygon, Contains, Intersects, Line, Point, Polygon, Rotate};
 use rstar::{RTreeObject, AABB};
 
 use crate::{
@@ -148,6 +148,13 @@ impl SegShape {
 
         polygon![p1.0, p2.0, p3.0, p4.0]
     }
+
+    pub(crate) fn middle_line(&self) -> Line {
+        Line {
+            start: self.from.into(),
+            end: self.to.into(),
+        }
+    }
 }
 
 impl MeasureLength for SegShape {

src/layout/layout.rs

@@ -19,7 +19,7 @@ use crate::{
         gear::GearIndex,
         graph::{GetMaybeNet, IsInLayer, MakePrimitive, PrimitiveIndex, PrimitiveWeight},
         loose::LooseIndex,
-        primitive::MakePrimitiveShape,
+        primitive::{GetWeight, MakePrimitiveShape, Primitive},
         rules::AccessRules,
         seg::{
             FixedSegIndex, FixedSegWeight, LoneLooseSegIndex, LoneLooseSegWeight, SegIndex,
@@ -28,16 +28,18 @@ use crate::{
         Cane, Collect, Drawing, DrawingEdit, DrawingException, Infringement,
     },
     geometry::{
+        compound::ManageCompounds,
         edit::ApplyGeometryEdit,
         primitive::{AccessPrimitiveShape, PrimitiveShape, SegShape},
         shape::{AccessShape, Shape},
-        GenericNode,
+        GenericNode, GetSetPos,
     },
     graph::{GenericIndex, GetPetgraphIndex},
     layout::{
         poly::{MakePolygon, Poly, PolyWeight},
         via::{Via, ViaWeight},
     },
+    math::{HomogeneousLine, LineIntersection},
 };
 
 /// Represents a weight for various compounds
@@ -346,11 +348,6 @@ impl<R: AccessRules> Layout<R> {
     pub fn center_of_high_level_node(&self, node: NodeIndex) -> Option<Point> {
         match node {
             NodeIndex::Primitive(primitive) => {
-                use crate::drawing::{
-                    graph::MakePrimitive,
-                    primitive::{GetWeight, Primitive},
-                };
-                use crate::geometry::{compound::ManageCompounds, GetSetPos};
                 if self
                     .drawing()
                     .geometry()
@@ -366,10 +363,7 @@ impl<R: AccessRules> Layout<R> {
                     _ => None,
                 }
             }
-            NodeIndex::Compound(cwidx) => {
-                use crate::geometry::shape::AccessShape;
-                Some(self.node_shape(node).center())
-            }
+            NodeIndex::Compound(cwidx) => Some(self.node_shape(node).center()),
         }
     }
 
@@ -385,15 +379,21 @@ impl<R: AccessRules> Layout<R> {
         assert_ne!(left, right);
         let left_pos = self.node_shape(left).center();
         let right_pos = self.node_shape(right).center();
-        let delta = right_pos - left_pos;
-        let delta_len = delta.y().hypot(delta.x());
-        let fake_seg = PrimitiveShape::Seg(SegShape {
-            from: left_pos,
-            to: right_pos,
-            // TODO: being able to use a zero-width segshape here would be optimal
-            // check if that works
-            width: f64::EPSILON,
-        });
+        let ltr_line = geo::Line {
+            start: left_pos.into(),
+            end: right_pos.into(),
+        };
+        let fake_seg = SegShape {
+            from: left_pos.into(),
+            to: right_pos.into(),
+            width: f64::EPSILON * 16.0,
+        };
+        let mut orig_hline = HomogeneousLine::from(ltr_line);
+        orig_hline.make_normal_unit();
+        let orig_hline = orig_hline;
+        let location_denom = orig_hline.segment_interval(&ltr_line);
+        let location_start = location_denom.start();
+        let location_denom = location_denom.end() - location_denom.start();
 
         let mut bands: Vec<_> = self
             .drawing
@@ -418,12 +418,36 @@ impl<R: AccessRules> Layout<R> {
                 let shape = prim.primitive(&self.drawing).shape();
                 (prim, loose, shape)
             })
-            .filter(|(_, _, shape)| shape.intersects(&fake_seg))
-            .map(|(_, loose, shape)| {
+            .filter_map(|(prim, loose, shape)| {
                 let band_uid = self.drawing.loose_band_uid(loose);
-                // TODO: check that the resulting order is correct
-                let location = crate::math::dot_product(delta, shape.center()) / delta_len;
-                (location, band_uid)
+                let loose_hline = orig_hline.orthogonal_through(&match shape {
+                    PrimitiveShape::Seg(seg) => {
+                        let seg_hline = HomogeneousLine::from(seg.middle_line());
+                        match orig_hline.intersects(&seg_hline) {
+                            LineIntersection::Empty => return None,
+                            LineIntersection::Overlapping => shape.center(),
+                            LineIntersection::Point(pt) => pt,
+                        }
+                    }
+                    _ => {
+                        if !fake_seg.intersects(&shape) {
+                            return None;
+                        }
+                        shape.center()
+                    }
+                });
+                let location = (loose_hline.offset - location_start) / location_denom;
+                log::trace!(
+                    "intersection ({:?}) with {:?} is at {:?}",
+                    band_uid,
+                    shape,
+                    location
+                );
+                if 0.0 <= location && location <= 1.0 {
+                    Some((location, band_uid))
+                } else {
+                    None
+                }
             })
             .collect();
         bands.sort_by(|a, b| f64::total_cmp(&a.0, &b.0));

src/math/mod.rs

@@ -9,6 +9,94 @@ pub use specctra_core::math::{Circle, PointWithRotation};
 mod tangents;
 pub use tangents::*;
 
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub enum LineIntersection {
+    Empty,
+    Overlapping,
+    Point(Point),
+}
+
+#[derive(Clone, Copy, Debug)]
+pub struct HomogeneousLine<T> {
+    // x0*x + y0*y = offset
+    pub x: T,
+    pub y: T,
+    pub offset: T,
+}
+
+impl From<Line> for HomogeneousLine<f64> {
+    fn from(l: Line) -> Self {
+        // the normal vector is perpendicular to the line
+        let normal = geo::point! {
+            x: l.dy(),
+            y: -l.dx(),
+        };
+        Self {
+            x: normal.0.x,
+            y: normal.0.y,
+            offset: perp_dot_product(l.end.into(), l.start.into()),
+        }
+    }
+}
+
+impl HomogeneousLine<f64> {
+    pub fn evaluate_at(&self, pt: Point) -> f64 {
+        self.x * pt.x() + self.y * pt.y() - self.offset
+    }
+
+    pub fn angle(&self) -> f64 {
+        self.y.atan2(self.x)
+    }
+
+    pub fn make_normal_unit(&mut self) {
+        let normal_len = self.y.hypot(self.x);
+        if normal_len > (f64::EPSILON * 16.0) {
+            self.x /= normal_len;
+            self.y /= normal_len;
+        }
+    }
+
+    pub fn intersects(&self, b: &Self) -> LineIntersection {
+        const ALMOST_ZERO: f64 = f64::EPSILON * 16.0;
+        let (mut a, mut b) = (*self, *b);
+        let _ = (a.make_normal_unit(), b.make_normal_unit());
+        let apt = geo::point! { x: a.x, y: a.y };
+        let bpt = geo::point! { x: b.x, y: b.y };
+        let det = perp_dot_product(apt, bpt);
+        let rpx = b.y * a.offset - a.y * b.offset;
+        let rpy = -b.x * a.offset + a.x * b.offset;
+
+        if det.abs() > ALMOST_ZERO {
+            LineIntersection::Point(geo::point! { x: rpx, y: rpy } / det)
+        } else if rpx.abs() <= ALMOST_ZERO && rpy.abs() <= ALMOST_ZERO {
+            LineIntersection::Overlapping
+        } else {
+            LineIntersection::Empty
+        }
+    }
+
+    /// project the point `pt` onto this line, and generate a new line which is orthogonal
+    /// to `self`, and goes through `pt`.
+    #[inline]
+    pub fn orthogonal_through(&self, pt: &Point) -> Self {
+        Self {
+            // recover the original parallel vector
+            x: -self.y,
+            y: self.x,
+            offset: self.x * pt.0.y - self.y * pt.0.x,
+        }
+    }
+
+    pub fn segment_interval(&self, line: &Line) -> core::ops::RangeInclusive<f64> {
+        // recover the original parallel vector
+        let parv = geo::point! {
+            x: -self.y,
+            y: self.x,
+        };
+        dot_product(parv, line.start.into())..=dot_product(parv, line.end.into())
+    }
+}
+
 pub fn intersect_circles(circle1: &Circle, circle2: &Circle) -> Vec<Point> {
     let delta = circle2.pos - circle1.pos;
     let d = Euclidean::distance(&circle2.pos, &circle1.pos);