From 713262f756beafab0846b04ecded980b38c61de4 Mon Sep 17 00:00:00 2001
From: Patrick Brosi <brosi@informatik.uni-freiburg.de>
Date: Wed, 8 Aug 2018 17:36:28 +0200
Subject: [PATCH] fix bug in line segment intersection

---
 src/util/geo/Geo.h          | 14 +++++++-------
 src/util/tests/TestMain.cpp |  2 ++
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/util/geo/Geo.h b/src/util/geo/Geo.h
index dc52d4a..570ba66 100644
--- a/src/util/geo/Geo.h
+++ b/src/util/geo/Geo.h
@@ -398,9 +398,9 @@ inline bool intersects(const LineSegment<T>& ls1, const LineSegment<T>& ls2) {
   return intersects(getBoundingBox(ls1), getBoundingBox(ls2)) &&
          (contains(ls1.first, ls2) || contains(ls1.second, ls2) ||
           contains(ls2.first, ls1) || contains(ls2.second, ls1) ||
-          ((crossProd(ls1.first, ls2) < 0) ^
-           (crossProd(ls1.second, ls2) < 0)) ||
-          ((crossProd(ls2.first, ls1) < 0) ^ (crossProd(ls2.second, ls1) < 0)));
+          (((crossProd(ls1.first, ls2) < 0) ^
+           (crossProd(ls1.second, ls2) < 0)) &&
+          ((crossProd(ls2.first, ls1) < 0) ^ (crossProd(ls2.second, ls1) < 0))));
 }
 
 // _____________________________________________________________________________
@@ -953,7 +953,7 @@ inline Point<T> projectOn(const Point<T>& a, const Point<T>& b,
   if (doubleEq(a.getX(), c.getX()) && doubleEq(a.getY(), c.getY())) return a;
   if (doubleEq(b.getX(), c.getX()) && doubleEq(b.getY(), c.getY())) return b;
 
-  double x, y;
+  T x, y;
 
   if (c.getX() == a.getX()) {
     // infinite slope
@@ -963,8 +963,8 @@ inline Point<T> projectOn(const Point<T>& a, const Point<T>& b,
     double m = (double)(c.getY() - a.getY()) / (c.getX() - a.getX());
     double bb = (double)a.getY() - (m * a.getX());
 
-    x = (m * b.getY() + b.getX() - m * bb) / (m * m + 1);
-    y = (m * m * b.getY() + m * b.getX() + bb) / (m * m + 1);
+    x = (m * b.getY() + b.getX() - m * bb) / (m * m + 1.0);
+    y = (m * m * b.getY() + m * b.getX() + bb) / (m * m + 1.0);
   }
 
   Point<T> ret = Point<T>(x, y);
@@ -1017,7 +1017,7 @@ inline double parallelity(const Box<T>& box, const MultiLine<T>& multiline) {
     ret += parallelity(box, l);
   }
 
-  return ret / multiline.size();
+  return ret / static_cast<float>(multiline.size());
 }
 
 // _____________________________________________________________________________
diff --git a/src/util/tests/TestMain.cpp b/src/util/tests/TestMain.cpp
index dbf1eb5..4e5ec62 100644
--- a/src/util/tests/TestMain.cpp
+++ b/src/util/tests/TestMain.cpp
@@ -976,6 +976,8 @@ CASE("geometry") {
   EXPECT(geo::intersects(lsb, lsb));
   EXPECT(!geo::intersects(lsa, lsc));
 
+  EXPECT(!geo::intersects(geo::Point<double>(871569.2, 6104550.4), geo::Point<double>(871581.2, 6104536), geo::Point<double>(871580.3, 6104541.3), geo::Point<double>(871625.7, 6104510.1)));
+
   geo::Line<double> l{geo::Point<double>(1, 1), geo::Point<double>(2, 2), geo::Point<double>(2, 4)};
   EXPECT(!geo::contains(geo::Point<double>(1, 2), l));
   EXPECT(geo::contains(geo::Point<double>(2, 2), l));