changes

src/util/geo/Box.h

@@ -0,0 +1,39 @@
+// Copyright 2016, University of Freiburg,
+// Chair of Algorithms and Data Structures.
+// Author: Patrick Brosi <brosi@informatik.uni-freiburg.de>
+
+#ifndef UTIL_GEO_BOX_H_
+#define UTIL_GEO_BOX_H_
+
+#include "./Point.h"
+
+namespace util {
+namespace geon {
+
+template <typename T>
+class Box {
+ public:
+  // maximum inverse box as default value of box
+  Box() : _ll(std::numeric_limits<T>::max), _ur(std::numeric_limits<T>::min) {}
+  Box(const Point<T>& ll, const Point<T>& ur) : _ll(ll), _ur(ur) {}
+  const Point<T>& getLowerLeft() const { return _ll; }
+  const Point<T>& getUpperRight() const { return _ur; }
+
+  void setLowerLeft(const Point<T>& ll) { _ll = ll; }
+  void setUpperRight(const Point<T>& ur) { _ur = ur; }
+
+  bool operator==(const Box<T>& b) const {
+    return getLowerLeft() == b.getLowerLeft() &&
+           getUpperRight == b.getUpperRight();
+  }
+
+  bool operator!=(const Box<T>& p) const { return !(*this == p); }
+
+ private:
+  Point<T> _ll, _ur;
+};
+
+}  // namespace geon
+}  // namespace util
+
+#endif  // UTIL_GEO_BOX_H_

src/util/geo/GeoNew.h

@@ -8,6 +8,8 @@
 
 #include <math.h>
 #include "util/Misc.h"
+#include "util/geo/Box.h"
+#include "util/geo/Line.h"
 #include "util/geo/Point.h"
 
 // -------------------
@@ -23,13 +25,13 @@ typedef Point<double> DPoint;
 typedef Point<float> FPoint;
 typedef Point<int> IPoint;
 
-// typedef Line<double> DLine;
+typedef Line<double> DLine;
-// typedef Line<float> FLine;
+typedef Line<float> FLine;
-// typedef Line<int> ILine;
+typedef Line<int> ILine;
 
-// typedef Box<double> DBox;
+typedef Box<double> DBox;
-// typedef Box<float> FBox;
+typedef Box<float> FBox;
-// typedef Box<int> IBox;
+typedef Box<int> IBox;
 
 // typedef Polygon<double> DPolygon;
 // typedef Polygon<float> FPolygon;
@@ -37,7 +39,8 @@ typedef Point<int> IPoint;
 
 // _____________________________________________________________________________
 // template <typename T>
-// inline Line<T> rotate(const Line<T>& geo, double deg, const Point<T>& center) {
+// inline Line<T> rotate(const Line<T>& geo, double deg, const Point<T>& center)
+// {
 // Line<T> ret;
 
 // bgeo::strategy::transform::translate_transformer<T, 2, 2> translate(
@@ -81,13 +84,13 @@ typedef Point<int> IPoint;
 // }
 
 // _____________________________________________________________________________
-// template <typename T>
+template <typename T>
-// inline Box<T> pad(const Box<T>& box, double padding) {
+inline Box<T> pad(const Box<T>& box, double padding) {
-  // return Box<T>(Point<T>(box.min_corner().getX() - padding,
+  return Box<T>(Point<T>(box.getLowerLeft().getX() - padding,
-                         // box.min_corner().getY() - padding),
+                         box.getLowerLeft().getY() - padding),
-                // Point<T>(box.max_corner().getX() + padding,
+                Point<T>(box.getUpperRight().getX() + padding,
-                         // box.max_corner().getY() + padding));
+                         box.getUpperRight().getY() + padding));
-// }
+}
 
 // _____________________________________________________________________________
 template <typename T>
@@ -141,24 +144,24 @@ inline Point<T> move(const Point<T>& geo, T x, T y) {
 // _____________________________________________________________________________
 template <typename T>
 inline Box<T> minbox() {
-  //return bgeo::make_inverse<Box<T>>();
+  return Box<T>();
 }
 
 // _____________________________________________________________________________
 // template <typename T>
 // inline RotatedBox<T> shrink(const RotatedBox<T>& b, double d) {
 // double xd =
-      // b.b.max_corner().getX() - b.b.min_corner().getX();
+// b.b.getUpperRight().getX() - b.b.getLowerLeft().getX();
 // double yd =
-      // b.b.max_corner().getY() - b.b.min_corner().getY();
+// b.b.getUpperRight().getY() - b.b.getLowerLeft().getY();
 
 // if (xd <= 2 * d) d = xd / 2 - 1;
 // if (yd <= 2 * d) d = yd / 2 - 1;
 
-  // Box<T> r(Point<T>(b.b.min_corner().getX() + d,
+// Box<T> r(Point<T>(b.b.getLowerLeft().getX() + d,
-                    // b.b.min_corner().getY() + d),
+// b.b.getLowerLeft().getY() + d),
-           // Point<T>(b.b.max_corner().getX() - d,
+// Point<T>(b.b.getUpperRight().getX() - d,
-                    // b.b.max_corner().getY() - d));
+// b.b.getUpperRight().getY() - d));
 
 // return RotatedBox<T>(r, b.rotateDeg, b.center);
 // }
@@ -167,14 +170,26 @@ inline Box<T> minbox() {
 inline bool doubleEq(double a, double b) { return fabs(a - b) < 0.000001; }
 
 // _____________________________________________________________________________
-// template <typename Geometry, typename Box>
+template <typename T>
-// inline bool contains(const Geometry& geom, const Box& box) {
+inline bool contains(const Point<T>& p, const Box<T>& box) {
-  // return bgeo::within(geom, box);
+  return p.getX() >= box.getLowerLeft().getX() &&
-// }
+         p.getX() <= box.getUpperRight().getX() &&
+         p.getY() >= box.getLowerLeft().getY() &&
+         p.getY() <= box.getUpperRight().getY();
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool contains(const Line<T>& l, const Box<T>& box) {
+  for (const auto& p : l)
+    if (!contains(p, box)) return false;
+  return true;
+}
 
 // _____________________________________________________________________________
 // template <typename T>
-// inline bool contains(const Point<T>& p1, const Point<T>& q1, const Point<T>& p2,
+// inline bool contains(const Point<T>& p1, const Point<T>& q1, const Point<T>&
+// p2,
 // const Point<T>& q2) {
 // Line<T> a, b;
 // a.push_back(p1);
@@ -187,7 +202,8 @@ inline bool doubleEq(double a, double b) { return fabs(a - b) < 0.000001; }
 
 // _____________________________________________________________________________
 // template <typename T>
-// inline bool contains(T p1x, T p1y, T q1x, T q1y, T p2x, T p2y, T q2x, T q2y) {
+// inline bool contains(T p1x, T p1y, T q1x, T q1y, T p2x, T p2y, T q2x, T q2y)
+// {
 // Point<T> p1(p1x, p1y);
 // Point<T> q1(q1x, q1y);
 // Point<T> p2(p2x, p2y);
@@ -215,7 +231,8 @@ inline bool doubleEq(double a, double b) { return fabs(a - b) < 0.000001; }
 
 // _____________________________________________________________________________
 // template <typename T>
-// inline bool intersects(T p1x, T p1y, T q1x, T q1y, T p2x, T p2y, T q2x, T q2y) {
+// inline bool intersects(T p1x, T p1y, T q1x, T q1y, T p2x, T p2y, T q2x, T
+// q2y) {
 /*
  * checks whether two line segments intersect
  */
@@ -254,10 +271,8 @@ inline Point<T> intersection(const Point<T>& p1, const Point<T>& q1,
   /*
    * calculates the intersection between two line segments
    */
-  return intersection(p1.getX(), p1.getY(),
+  return intersection(p1.getX(), p1.getY(), q1.getX(), q1.getY(), p2.getX(),
-                      q1.getX(), q1.getY(),
+                      p2.getY(), q2.getX(), q2.getY());
-                      p2.getX(), p2.getY(),
-                      q2.getX(), q2.getY());
 }
 
 // _____________________________________________________________________________
@@ -280,10 +295,8 @@ inline bool lineIntersects(const Point<T>& p1, const Point<T>& q1,
   /*
    * checks whether two lines intersect
    */
-  return lineIntersects(p1.getX(), p1.getY(),
+  return lineIntersects(p1.getX(), p1.getY(), q1.getX(), q1.getY(), p2.getX(),
-                        q1.getX(), q1.getY(),
+                        p2.getY(), q2.getX(), q2.getY());
-                        p2.getX(), p2.getY(),
-                        q2.getX(), q2.getY());
 }
 
 // _____________________________________________________________________________
@@ -296,8 +309,7 @@ inline double angBetween(double p1x, double p1y, double q1x, double q1y) {
 // _____________________________________________________________________________
 template <typename T>
 inline double angBetween(const Point<T>& p1, const Point<T>& q1) {
-  return angBetween(p1.getX(), p1.getY(),
+  return angBetween(p1.getX(), p1.getY(), q1.getX(), q1.getY());
-                    q1.getX(), q1.getY());
 }
 
 // _____________________________________________________________________________
@@ -321,23 +333,36 @@ inline double innerProd(double x1, double y1, double x2, double y2, double x3,
 
 // _____________________________________________________________________________
 template <typename T>
-inline double innerProd(const Point<T>& a, const Point<T>& b, const Point<T>& c) {
+inline double innerProd(const Point<T>& a, const Point<T>& b,
-  return innerProd(a.template getX(), a.template getY(),
+                        const Point<T>& c) {
-                   b.template getX(), b.template getY(),
+  return innerProd(a.template getX(), a.template getY(), b.template getX(),
-                   c.template getX(), c.template getY());
+                   b.template getY(), c.template getX(), c.template getY());
 }
 
 // _____________________________________________________________________________
 template <typename T>
 inline double dist(const Point<T>& p1, const Point<T>& p2) {
-  return 0;  // TODo
+  return sqrt((p1.getX() - p2.getX()) * (p1.getX() - p2.getX()) +
+              (p1.getY() - p2.getY()) * (p1.getY() - p2.getY()));
 }
 
 // _____________________________________________________________________________
 template <typename T>
-inline std::string getWKT(const Point<T>* p) {
+inline std::string getWKT(const Point<T>& p) {
   std::stringstream ss;
-  // TODO
+  ss << "POINT (" << p.getX() << " " << p.getY() << ")";
+  return ss.str();
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline std::string getWKT(const Line<T>& l) {
+  std::stringstream ss;
+  ss << "LINESTRING (";
+  for (size_t i = 0; i < l.size(); i++) {
+    if (i) ss << ", ";
+    ss << l[i].getX() << " " << l[i].getY();
+  }
   return ss.str();
 }
 
@@ -347,12 +372,26 @@ inline double len(const Point<T>& g) {
   return 0;
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline double len(const Line<T>& g) {
+  double ret = 0;
+  for (size_t i = 1; i < g.size(); i++) ret += dist(g[i - 1], g[i]);
+  return ret;
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline Point<T> simplify(const Point<T>& g, double d) {
   return g;
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline Line<T> simplify(const Line<T>& g, double d) {
+  return g;
+}
+
 // _____________________________________________________________________________
 inline double distToSegment(double lax, double lay, double lbx, double lby,
                             double px, double py) {
@@ -374,24 +413,17 @@ inline double distToSegment(double lax, double lay, double lbx, double lby,
 template <typename T>
 inline double distToSegment(const Point<T>& la, const Point<T>& lb,
                             const Point<T>& p) {
-  return distToSegment(la.getX(), la.getY(),
+  return distToSegment(la.getX(), la.getY(), lb.getX(), lb.getY(), p.getX(),
-                       lb.getX(), lb.getY(),
+                       p.getY());
-                       p.getX(), p.getY());
 }
 
 // _____________________________________________________________________________
 template <typename T>
 inline Point<T> projectOn(const Point<T>& a, const Point<T>& b,
                           const Point<T>& c) {
-  if (doubleEq(a.getX(), b.getX()) &&
+  if (doubleEq(a.getX(), b.getX()) && doubleEq(a.getY(), b.getY())) return a;
-      doubleEq(a.getY(), b.getY()))
+  if (doubleEq(a.getX(), c.getX()) && doubleEq(a.getY(), c.getY())) return a;
-    return a;
+  if (doubleEq(b.getX(), c.getX()) && doubleEq(b.getY(), c.getY())) return 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;
 
@@ -400,13 +432,11 @@ inline Point<T> projectOn(const Point<T>& a, const Point<T>& b,
     x = a.getX();
     y = b.getY();
   } else {
-    double m = (double)(c.getY() - a.getY()) /
+    double m = (double)(c.getY() - a.getY()) / (c.getX() - a.getX());
-               (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) /
+    y = (m * m * b.getY() + m * b.getX() + bb) / (m * m + 1);
-        (m * m + 1);
   }
 
   Point<T> ret = Point<T>(x, y);
@@ -424,18 +454,18 @@ template <typename T>
 inline double parallelity(const Box<T>& box, const Line<T>& line) {
   double ret = M_PI;
 
-  double a = angBetween(box.min_corner(),
+  double a = angBetween(
-                        Point<T>(box.min_corner().getX(),
+      box.getLowerLeft(),
-                                 box.max_corner().getY()));
+      Point<T>(box.getLowerLeft().getX(), box.getUpperRight().getY()));
-  double b = angBetween(box.min_corner(),
+  double b = angBetween(
-                        Point<T>(box.max_corner().getX(),
+      box.getLowerLeft(),
-                                 box.min_corner().getY()));
+      Point<T>(box.getUpperRight().getX(), box.getLowerLeft().getY()));
-  double c = angBetween(box.max_corner(),
+  double c = angBetween(
-                        Point<T>(box.min_corner().getX(),
+      box.getUpperRight(),
-                                 box.max_corner().getY()));
+      Point<T>(box.getLowerLeft().getX(), box.getUpperRight().getY()));
-  double d = angBetween(box.max_corner(),
+  double d = angBetween(
-                        Point<T>(box.max_corner().getX(),
+      box.getUpperRight(),
-                                 box.min_corner().getY()));
+      Point<T>(box.getUpperRight().getX(), box.getLowerLeft().getY()));
 
   double e = angBetween(line.front(), line.back());
 
@@ -498,44 +528,55 @@ inline double parallelity(const Box<T>& box, const Line<T>& line) {
 
 // _____________________________________________________________________________
 template <typename T>
-// inline Box<T> getBoundingBox(Point<T> pol) {
+inline Box<T> getBoundingBox(const Point<T>& p) {
-  // Box<T> tmpBox;
+  return Box<T>(p, p);
-  // bgeo::envelope(pol, tmpBox);
+}
-  // return tmpBox;
-// }
 
 // _____________________________________________________________________________
-// template <typename T>
+template <typename T>
-// inline Box<T> extendBox(const Box<T>& a, Box<T> b) {
+inline Box<T> getBoundingBox(const Line<T>& l) {
-  // bgeo::expand(b, a);
+  Box<T> ret;
-  // return b;
+  for (const auto& p : l) extendBox(p, ret);
-// }
+  return ret;
+}
 
 // _____________________________________________________________________________
-template <typename G, typename T>
+template <typename T>
-inline Box<T> extendBox(G pol, Box<T> b) {
+inline Box<T> extendBox(const Box<T>& a, Box<T> b) {
-  Box<T> tmp;
+  b = extendBox(a.getLowerLeft(), b);
-  bgeo::envelope(pol, tmp);
+  b = extendBox(a.getUpperRight(), b);
-  bgeo::expand(b, tmp);
   return b;
 }
 
 // _____________________________________________________________________________
-// template <typename T>
+template <typename T>
-// inline double commonArea(const Box<T>& ba, const Box<T>& bb) {
+inline Box<T> extendBox(const Point<T>& p, Box<T> b) {
-  // T l = std::max(ba.min_corner().getX(),
+  if (p.getX() < b.getLowerLeft().getX()) b.getLowerLeft().setX(p.getX());
-                 // bb.min_corner().getX());
+  if (p.getY() < b.getLowerLeft().getY()) b.getLowerLeft().setY(p.getY());
-  // T r = std::min(ba.max_corner().getX(),
-                 // bb.max_corner().getX());
-  // T b = std::max(ba.min_corner().getY(),
-                 // bb.min_corner().getY());
-  // T t = std::min(ba.max_corner().getY(),
-                 // bb.max_corner().getY());
 
-  // if (l > r || b > t) return 0;
+  if (p.getX() > b.getUpperRight().getX()) b.getUpperRight().setX(p.getX());
+  if (p.getY() > b.getUpperRight().getY()) b.getUpperRight().setY(p.getY());
+  return b;
+}
 
-  // return (r - l) * (t - b);
+// _____________________________________________________________________________
-// }
+template <typename T>
+inline Box<T> extendBox(const Line<T>& l, Box<T> b) {
+  for (const auto& p : l) b = extendBox(p, b);
+  return b;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline double commonArea(const Box<T>& ba, const Box<T>& bb) {
+  T l = std::max(ba.getLowerLeft().getX(), bb.getLowerLeft().getX());
+  T r = std::min(ba.getUpperRight().getX(), bb.getUpperRight().getX());
+  T b = std::max(ba.getLowerLeft().getY(), bb.getLowerLeft().getY());
+  T t = std::min(ba.getUpperRight().getY(), bb.getUpperRight().getY());
+
+  if (l > r || b > t) return 0;
+  return (r - l) * (t - b);
+}
 
 // _____________________________________________________________________________
 // template <typename T, template <typename> typename Geometry>
@@ -620,14 +661,12 @@ inline Line<T> densify(const Line<T>& l, double d) {
 
   for (size_t i = 1; i < l.size(); i++) {
     double segd = dist(l[i - 1], l[i]);
-    double dx =
+    double dx = (l[i].getX() - l[i - 1].getX()) / segd;
-        (l[i].getX() - l[i-1].getX()) / segd;
+    double dy = (l[i].getY() - l[i - 1].getY()) / segd;
-    double dy =
-        (l[i].getY() - l[i-1].getY()) / segd;
     double curd = d;
     while (curd < segd) {
-      ret.push_back(Point<T>(l[i-1].getX() + dx * curd,
+      ret.push_back(
-                             l[i-1].getY() + dy * curd));
+          Point<T>(l[i - 1].getX() + dx * curd, l[i - 1].getY() + dy * curd));
       curd += d;
     }
 
@@ -683,29 +722,29 @@ inline double frechetDist(const Line<T>& a, const Line<T>& b, double d) {
 // _____________________________________________________________________________
 template <typename T>
 inline double accFrechetDistC(const Line<T>& a, const Line<T>& b, double d) {
-
   auto p = densify(a, d);
   auto q = densify(b, d);
 
   std::vector<std::vector<double>> ca(p.size(),
                                       std::vector<double>(q.size(), 0));
 
-  for (size_t i = 0; i < p.size(); i++) ca[i][0] = std::numeric_limits<double>::infinity();
+  for (size_t i = 0; i < p.size(); i++)
-  for (size_t j = 0; j < q.size(); j++) ca[0][j] = std::numeric_limits<double>::infinity();
+    ca[i][0] = std::numeric_limits<double>::infinity();
+  for (size_t j = 0; j < q.size(); j++)
+    ca[0][j] = std::numeric_limits<double>::infinity();
   ca[0][0] = 0;
 
   for (size_t i = 1; i < p.size(); i++) {
     for (size_t j = 1; j < q.size(); j++) {
       double d = util::geo::dist(p[i], q[j]) * util::geo::dist(p[i], p[i - 1]);
-      ca[i][j] = d + std::min(ca[i-1][j], std::min(ca[i][j-1], ca[i-1][j-1]));
+      ca[i][j] =
+          d + std::min(ca[i - 1][j], std::min(ca[i][j - 1], ca[i - 1][j - 1]));
     }
   }
 
   return ca[p.size() - 1][q.size() - 1];
 }
 
-
-
 // _____________________________________________________________________________
 template <typename T>
 inline Point<T> latLngToWebMerc(double lat, double lng) {

src/util/geo/Line.h

@@ -0,0 +1,21 @@
+// Copyright 2016, University of Freiburg,
+// Chair of Algorithms and Data Structures.
+// Author: Patrick Brosi <brosi@informatik.uni-freiburg.de>
+
+#ifndef UTIL_GEO_LINE_H_
+#define UTIL_GEO_LINE_H_
+
+#include <vector>
+#include "./Point.h"
+
+namespace util {
+namespace geon {
+
+template<typename T>
+using Line = std::vector<Point<T>>;
+
+
+}  // namespace geon
+}  // namespace util
+
+#endif  // UTIL_GEO_LINE_H_

src/util/geo/Point.h

@@ -5,9 +5,6 @@
 #ifndef UTIL_GEO_POINT_H_
 #define UTIL_GEO_POINT_H_
 
-#include <set>
-#include <vector>
-
 namespace util {
 namespace geon {
 
@@ -21,9 +18,20 @@ class Point {
   void setX(T x) { _x = x; }
   void setY(T y) { _y = y; }
 
+  Point<T> operator+(const Point<T>& p) const {
+    return Point<T>(_x + p.getX(), _y + p.getY());
+  }
+
+  bool operator==(const Point<T>& p) const {
+    return p.getX() == _x && p.getY() == _y;
+  }
+
+  bool operator!=(const Point<T>& p) const {
+    return !(*this == p);
+  }
+
  private:
   T _x, _y;
-
 };
 
 }  // namespace geon

src/util/tests/TestMain.cpp

@@ -846,6 +846,43 @@ CASE("nullable") {
 
 // ___________________________________________________________________________
 CASE("geometry") {
+  geon::Point<double> a(1, 2);
+  geon::Point<double> b(2, 3);
+  geon::Point<double> c(4, 5);
+  EXPECT(a.getX() == approx(1));
+  EXPECT(a.getY() == approx(2));
+
+  a.setX(3);
+  EXPECT(a.getX() == approx(3));
+  EXPECT(a.getY() == approx(2));
+
+  a.setY(4);
+  EXPECT(a.getX() == approx(3));
+  EXPECT(a.getY() == approx(4));
+
+  auto d = a + b;
+  EXPECT(d.getX() == approx(5));
+  EXPECT(d.getY() == approx(7));
+
+  a.setX(1);
+  a.setY(2);
+
+  EXPECT(geon::dist(a, a) == approx(0));
+  EXPECT(geon::dist(a, b) == approx(sqrt(2)));
+
+  d = d + d;
+
+  geon::Box<double> box(a, c);
+  EXPECT(geon::contains(a, box));
+  EXPECT(geon::contains(b, box));
+  EXPECT(geon::contains(c, box));
+  EXPECT(!geon::contains(d, box));
+
+  geon::Line<double> line{a, b, c};
+
+  EXPECT(geon::contains(line, box));
+  line.push_back(d);
+  EXPECT(!geon::contains(line, box));
 }
 
 };