Merge branch 'refactoring'

.travis.yml

@@ -0,0 +1,24 @@
+language: generic
+sudo: false
+dist: trusty
+
+addons:
+  apt:
+    sources:
+      - ubuntu-toolchain-r-test
+    packages:
+      - cmake
+
+before_script:
+  - mkdir build
+  - cd build
+  - cmake ..
+
+script:
+  - make -j4
+  - make test
+
+notifications:
+  email:
+    on_success: never
+    on_failure: never

CMakeLists.txt

@@ -13,9 +13,6 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/")
 set(EXECUTABLE_OUTPUT_PATH "${CMAKE_SOURCE_DIR}/build")
 
 
-find_package(Boost COMPONENTS system filesystem REQUIRED)
-include_directories("build" ${Boost_INCLUDE_DIRS})
-
 find_package(OpenMP)
 if(OPENMP_FOUND)
 	set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")

README.md

@@ -1,15 +1,4 @@
-```
+![Map-Matched path of a single train through Switzerland](geo/schweiz_mmatched.png?raw=true)
-                 _-====-__-=-__-===-__-=======-__
-               _(                               _)
-            OO(                                )
-         . o  '===-______-===-____-==-__-====='
-      .o
-     . ______          _______________
-   _()_||__|| __o^o___ | [] [] [] [] |
-  (           |      | |             |o
- /-OO----OO""="OO--OO"="OO---------OO"
-############################################################
-```
 
 # pfaedle
 
@@ -18,8 +7,7 @@ Precise map-matching for public transit schedules (GTFS data).
 ## Requirements
 
  * `cmake`
- * `gcc` >= 4.8
+ * `gcc` >= 4.8 (may work on lower versions, untested)
- * `libboost-system` >= 1.56, `libboost-filesystem` >= 1.56, `libboost-geometry` >= 1.56
 
 ## Building and Installation
 
@@ -50,16 +38,29 @@ pfaedle -c <CFG FILE> -x <OSM FILE> <GTFS INPUT FOLDER>
 
 A shape'd version of the input GTFS feed will be written to `./gtfs-out`.
 
-A default configuration file `pfaedle.cfg` can be found in this repo.
-
 By default, shapes are only calculated for trips that don't have a shape in the
 input feed. To drop all existing shapes, use the `-D` flag.
 
+For example, you may generate (and replace existing, see -D parameter) shapes for the GTFS dataset for Freiburg like this:
+
+```
+$ mkdir freiburg_gtfs && cd freiburg_gtfs
+$ wget https://fritz.freiburg.de/csv_Downloads/VAGFR.zip
+$ unzip VAGFR.zip
+$ wget http://download.geofabrik.de/europe/germany/baden-wuerttemberg/freiburg-regbez-latest.osm.bz2
+$ bunzip2 freiburg-regbez-latest.osm.bz2
+$ mkdir gtfs-out
+$ pfaedle -D -c pfaedle.cfg -x freiburg-regbez-latest.osm .
+```
+
+A default configuration file `pfaedle.cfg` can be found in this repo.
+
+
 ## Generating shapes for a specific MOT
 
 To generate shapes only for a specific mot, use the `-m` option. Possible
 values are either `tram`, `bus`, `rail`, `subway`, `ferry`, `funicular`,
-`gondola`, `all`.
+`gondola`, `all` (default).
 
 Multiple values can be specified (comma separated).
 

src/cppgtfs

@@ -1 +1 @@
-Subproject commit 7c395530d3c5b2d5c147fdabaaeeb0d1babceb04
+Subproject commit 8d255970a97c31891277f63ba8b39a3ea2b7e133

src/pfaedle/PfaedleMain.cpp

@@ -154,7 +154,7 @@ int main(int argc, char** argv) {
   if (cfg.evaluate) ecoll.printStats(&std::cout);
 
   if (cfg.feedPaths.size()) {
-    LOG(INFO) << "Writing output GTFS...";
+    LOG(INFO) << "Writing output GTFS to " << cfg.outputPath << " ...";
     ad::cppgtfs::Writer w;
     w.write(&gtfs, cfg.outputPath);
   }

src/pfaedle/config/PfaedleConfig.h

@@ -21,6 +21,7 @@ struct Config {
       : dbgOutputPath("geo"),
         solveMethod("global"),
         evalPath("."),
+        outputPath("gtfs-out"),
         dropShapes(false),
         useHMM(false),
         writeGraph(false),

src/pfaedle/eval/Collector.cpp

@@ -118,8 +118,8 @@ double Collector::add(const Trip* t, const Shape* oldS, const Shape* newS,
   gjout.flush();
   fstr.close();
 
-  double fac = cos(2 * atan(exp((oldSegs.front().front().get<1>() +
+  double fac = cos(2 * atan(exp((oldSegs.front().front().getY() +
-                                 oldSegs.back().back().get<1>()) /
+                                 oldSegs.back().back().getY()) /
                                 6378137.0)) -
                    1.5707965);
 
@@ -401,7 +401,7 @@ std::pair<size_t, double> Collector::getDa(const std::vector<FLine>& a,
   // euclidean distance on web mercator is in meters on equator,
   // and proportional to cos(lat) in both y directions
   double fac =
-      cos(2 * atan(exp((a.front().front().get<1>() + a.back().back().get<1>()) /
+      cos(2 * atan(exp((a.front().front().getY() + a.back().back().getY()) /
                        6378137.0)) -
           1.5707965);
 

src/pfaedle/osm/BBoxIdx.cpp

@@ -31,10 +31,10 @@ bool BBoxIdx::contains(const Point<float>& p) const {
 // _____________________________________________________________________________
 util::geo::Box<float> BBoxIdx::getFullWebMercBox() const {
   return util::geo::FBox(
-      util::geo::latLngToWebMerc<float>(_root.box.min_corner().get<1>(),
+      util::geo::latLngToWebMerc<float>(_root.box.getLowerLeft().getY(),
-                                        _root.box.min_corner().get<0>()),
+                                        _root.box.getLowerLeft().getX()),
-      util::geo::latLngToWebMerc<float>(_root.box.max_corner().get<1>(),
+      util::geo::latLngToWebMerc<float>(_root.box.getUpperRight().getY(),
-                                        _root.box.max_corner().get<0>()));
+                                        _root.box.getUpperRight().getX()));
 }
 
 // _____________________________________________________________________________

src/pfaedle/osm/OsmBuilder.cpp

@@ -991,7 +991,7 @@ double OsmBuilder::webMercDistFactor(const util::geo::FPoint& a) const {
   // euclidean distance on web mercator is in meters on equator,
   // and proportional to cos(lat) in both y directions
 
-  double lat = 2 * atan(exp(a.get<1>() / 6378137.0)) - 1.5707965;
+  double lat = 2 * atan(exp(a.getY() / 6378137.0)) - 1.5707965;
   return cos(lat);
 }
 

src/pfaedle/router/Router.cpp

@@ -121,8 +121,8 @@ NDistHeur::NDistHeur(const RoutingOpts& rOpts,
   size_t c = 0;
   double x = 0, y = 0;
   for (auto to : tos) {
-    x += to->pl().getGeom()->get<0>();
+    x += to->pl().getGeom()->getX();
-    y += to->pl().getGeom()->get<1>();
+    y += to->pl().getGeom()->getY();
     c++;
   }
 
@@ -144,8 +144,8 @@ DistHeur::DistHeur(uint8_t minLvl, const RoutingOpts& rOpts,
   size_t c = 0;
   double x = 0, y = 0;
   for (auto to : tos) {
-    x += to->getFrom()->pl().getGeom()->get<0>();
+    x += to->getFrom()->pl().getGeom()->getX();
-    y += to->getFrom()->pl().getGeom()->get<1>();
+    y += to->getFrom()->pl().getGeom()->getY();
     c++;
   }
 

src/pfaedle/router/ShapeBuilder.cpp

@@ -313,8 +313,8 @@ ad::cppgtfs::gtfs::Shape* ShapeBuilder::getGtfsShape(
   for (const auto& hop : shp.hops) {
     const trgraph::Node* l = hop.start;
     if (hop.edges.size() == 0) {
-      FPoint ll = webMercToLatLng<float>(hop.start->pl().getGeom()->get<0>(),
+      FPoint ll = webMercToLatLng<float>(hop.start->pl().getGeom()->getX(),
-                                         hop.start->pl().getGeom()->get<1>());
+                                         hop.start->pl().getGeom()->getY());
 
       if (dist > -0.5)
         dist += webMercMeterDist(last, *hop.start->pl().getGeom());
@@ -324,7 +324,7 @@ ad::cppgtfs::gtfs::Shape* ShapeBuilder::getGtfsShape(
       last = *hop.start->pl().getGeom();
 
       if (dist - lastDist > 0.01) {
-        ret->addPoint(ShapePoint(ll.get<1>(), ll.get<0>(), dist, seq));
+        ret->addPoint(ShapePoint(ll.getY(), ll.getX(), dist, seq));
         seq++;
         lastDist = dist;
       }
@@ -333,9 +333,9 @@ ad::cppgtfs::gtfs::Shape* ShapeBuilder::getGtfsShape(
       last = *hop.end->pl().getGeom();
 
       if (dist - lastDist > 0.01) {
-        ll = webMercToLatLng<float>(hop.end->pl().getGeom()->get<0>(),
+        ll = webMercToLatLng<float>(hop.end->pl().getGeom()->getX(),
-                                    hop.end->pl().getGeom()->get<1>());
+                                    hop.end->pl().getGeom()->getY());
-        ret->addPoint(ShapePoint(ll.get<1>(), ll.get<0>(), dist, seq));
+        ret->addPoint(ShapePoint(ll.getY(), ll.getX(), dist, seq));
         seq++;
         lastDist = dist;
       }
@@ -351,8 +351,8 @@ ad::cppgtfs::gtfs::Shape* ShapeBuilder::getGtfsShape(
             dist = 0;
           last = cur;
           if (dist - lastDist > 0.01) {
-            FPoint ll = webMercToLatLng<float>(cur.get<0>(), cur.get<1>());
+            FPoint ll = webMercToLatLng<float>(cur.getX(), cur.getY());
-            ret->addPoint(ShapePoint(ll.get<1>(), ll.get<0>(), dist, seq));
+            ret->addPoint(ShapePoint(ll.getY(), ll.getX(), dist, seq));
             seq++;
             lastDist = dist;
           }
@@ -366,8 +366,8 @@ ad::cppgtfs::gtfs::Shape* ShapeBuilder::getGtfsShape(
             dist = 0;
           last = cur;
           if (dist - lastDist > 0.01) {
-            FPoint ll = webMercToLatLng<float>(cur.get<0>(), cur.get<1>());
+            FPoint ll = webMercToLatLng<float>(cur.getX(), cur.getY());
-            ret->addPoint(ShapePoint(ll.get<1>(), ll.get<0>(), dist, seq));
+            ret->addPoint(ShapePoint(ll.getY(), ll.getX(), dist, seq));
             seq++;
             lastDist = dist;
           }

src/util/Misc.h

@@ -5,6 +5,7 @@
 #ifndef UTIL_MISC_H_
 #define UTIL_MISC_H_
 
+#include <cmath>
 #include <chrono>
 
 #define UNUSED(expr) do { (void)(expr); } while (0)

src/util/geo/BezierCurve.h

@@ -47,7 +47,7 @@ class BezierCurve {
   Point<T> valueAt(double t) const;
 };
 
-#include "util/geo/PolyLine.tpp"
+#include "util/geo/BezierCurve.tpp"
 }
 }
 

src/util/geo/BezierCurve.tpp

@@ -4,8 +4,8 @@
 
 // _____________________________________________________________________________
 template <typename T>
-BezierCurve<T>::BezierCurve(const Point<T>& a, const Point<T>& b, const Point<T>& c,
+BezierCurve<T>::BezierCurve(const Point<T>& a, const Point<T>& b,
-                         const Point<T>& d)
+                            const Point<T>& c, const Point<T>& d)
     : _d(dist(a, d)) {
   assert(_d > 0);
   recalcPolynoms(a, b, c, d);
@@ -13,17 +13,17 @@ BezierCurve<T>::BezierCurve(const Point<T>& a, const Point<T>& b, const Point<T>
 
 // _____________________________________________________________________________
 template <typename T>
-void BezierCurve<T>::recalcPolynoms(const Point<T>& a, const Point<T>& b, const Point<T>& c,
+void BezierCurve<T>::recalcPolynoms(const Point<T>& a, const Point<T>& b,
-                                 const Point<T>& d) {
+                                    const Point<T>& c, const Point<T>& d) {
-  _xp.a = a.template get<0>();
+  _xp.a = a.getX();
-  _xp.b = 3.0 * (b.template get<0>() - a.template get<0>());
+  _xp.b = 3.0 * (b.getX() - a.getX());
-  _xp.c = 3.0 * (c.template get<0>() - b.template get<0>()) - _xp.b;
+  _xp.c = 3.0 * (c.getX() - b.getX()) - _xp.b;
-  _xp.d = d.template get<0>() - a.template get<0>() - _xp.c - _xp.b;
+  _xp.d = d.getX() - a.getX() - _xp.c - _xp.b;
 
-  _yp.a = a.template get<1>();
+  _yp.a = a.getY();
-  _yp.b = 3.0 * (b.template get<1>() - a.template get<1>());
+  _yp.b = 3.0 * (b.getY() - a.getY());
-  _yp.c = 3.0 * (c.template get<1>() - b.template get<1>()) - _yp.b;
+  _yp.c = 3.0 * (c.getY() - b.getY()) - _yp.b;
-  _yp.d = d.template get<1>() - a.template get<1>() - _yp.c - _yp.b;
+  _yp.d = d.getY() - a.getY() - _yp.c - _yp.b;
 
   _didRender = false;
 }

src/util/geo/Box.h

@@ -0,0 +1,89 @@
+// 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 geo {
+
+template <typename T>
+class Box {
+ public:
+  // maximum inverse box as default value of box
+  Box()
+      : _ll(std::numeric_limits<T>::max(), std::numeric_limits<T>::max()),
+        _ur(std::numeric_limits<T>::min(), 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; }
+
+  Point<T>& getLowerLeft() { return _ll; }
+  Point<T>& getUpperRight() { 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;
+};
+
+template <typename T>
+class RotatedBox {
+ public:
+  RotatedBox() : _box(), _deg(0), _center() {}
+  RotatedBox(const Box<T>& box)
+      : _box(box),
+        _deg(0),
+        _center(Point<T>(
+            (box.getUpperRight().getX() - box.getLowerLeft().getX()) / T(2),
+            (box.getUpperRight().getY() - box.getLowerLeft().getY()) / T(2))) {}
+  RotatedBox(const Point<T>& ll, const Point<T>& ur)
+      : _box(ll, ur),
+        _deg(0),
+        _center(Point<T>((ur.getX() - ll.getX()) / T(2),
+                         (ur.getY() - ll.getY()) / T(2))) {}
+  RotatedBox(const Box<T>& box, double deg)
+      : _box(box),
+        _deg(deg),
+        _center(Point<T>(
+            (box.getUpperRight().getX() - box.getLowerLeft().getX()) / T(2),
+            (box.getUpperRight().getY() - box.getLowerLeft().getY()) / T(2))) {}
+  RotatedBox(const Point<T>& ll, const Point<T>& ur, double deg)
+      : _box(ll, ur),
+        _deg(deg),
+        _center(Point<T>((ur.getX() - ll.getX()) / T(2),
+                         (ur.getY() - ll.getY()) / T(2))) {}
+  RotatedBox(const Box<T>& box, double deg, const Point<T>& center)
+      : _box(box), _deg(deg), _center(center) {}
+  RotatedBox(const Point<T>& ll, const Point<T>& ur, double deg,
+             const Point<T>& center)
+      : _box(ll, ur), _deg(deg), _center(center) {}
+
+  const Box<T>& getBox() const { return _box; }
+  Box<T>& getBox() { return _box; }
+
+  double getDegree() const { return _deg; }
+  const Point<T>& getCenter() const { return _center; }
+  Point<T>& getCenter() { return _center; }
+
+ private:
+  Box<T> _box;
+  double _deg;
+  Point<T> _center;
+};
+
+}  // namespace geo
+}  // namespace util
+
+#endif  // UTIL_GEO_BOX_H_

src/util/geo/Grid.h

@@ -7,6 +7,7 @@
 
 #include <set>
 #include <vector>
+#include <map>
 #include "util/geo/Geo.h"
 
 namespace util {

src/util/geo/Grid.tpp

@@ -30,9 +30,9 @@ Grid<V, G, T>::Grid(double w, double h, const Box<T>& bbox, bool bValIdx)
       _bb(bbox),
       _hasValIdx(bValIdx) {
   _width =
-      bbox.max_corner().template get<0>() - bbox.min_corner().template get<0>();
+      bbox.getUpperRight().getX() - bbox.getLowerLeft().getX();
   _height =
-      bbox.max_corner().template get<1>() - bbox.min_corner().template get<1>();
+      bbox.getUpperRight().getY() - bbox.getLowerLeft().getY();
 
   if (_width < 0 || _height < 0) {
     _width = 0;
@@ -58,11 +58,11 @@ Grid<V, G, T>::Grid(double w, double h, const Box<T>& bbox, bool bValIdx)
 template <typename V, template <typename> typename G, typename T>
 void Grid<V, G, T>::add(G<T> geom, V val) {
   Box<T> box = getBoundingBox(geom);
-  size_t swX = getCellXFromX(box.min_corner().template get<0>());
+  size_t swX = getCellXFromX(box.getLowerLeft().getX());
-  size_t swY = getCellYFromY(box.min_corner().template get<1>());
+  size_t swY = getCellYFromY(box.getLowerLeft().getY());
 
-  size_t neX = getCellXFromX(box.max_corner().template get<0>());
+  size_t neX = getCellXFromX(box.getUpperRight().getX());
-  size_t neY = getCellYFromY(box.max_corner().template get<1>());
+  size_t neY = getCellYFromY(box.getUpperRight().getY());
 
   for (size_t x = swX; x <= neX && x < _grid.size(); x++) {
     for (size_t y = swY; y <= neY && y < _grid[x].size(); y++) {
@@ -83,11 +83,11 @@ void Grid<V, G, T>::add(size_t x, size_t y, V val) {
 // _____________________________________________________________________________
 template <typename V, template <typename> typename G, typename T>
 void Grid<V, G, T>::get(const Box<T>& box, std::set<V>* s) const {
-  size_t swX = getCellXFromX(box.min_corner().template get<0>());
+  size_t swX = getCellXFromX(box.getLowerLeft().getX());
-  size_t swY = getCellYFromY(box.min_corner().template get<1>());
+  size_t swY = getCellYFromY(box.getLowerLeft().getY());
 
-  size_t neX = getCellXFromX(box.max_corner().template get<0>());
+  size_t neX = getCellXFromX(box.getUpperRight().getX());
-  size_t neY = getCellYFromY(box.max_corner().template get<1>());
+  size_t neY = getCellYFromY(box.getUpperRight().getY());
 
   for (size_t x = swX; x <= neX && x >= 0 && x < _xWidth; x++)
     for (size_t y = swY; y <= neY && y >= 0 && y < _yHeight; y++) get(x, y, s);
@@ -97,10 +97,10 @@ void Grid<V, G, T>::get(const Box<T>& box, std::set<V>* s) const {
 template <typename V, template <typename> typename G, typename T>
 void Grid<V, G, T>::get(const G<T>& geom, double d, std::set<V>* s) const {
   Box<T> a = getBoundingBox(geom);
-  Box<T> b(Point<T>(a.min_corner().template get<0>() - d,
+  Box<T> b(Point<T>(a.getLowerLeft().getX() - d,
-                    a.min_corner().template get<1>() - d),
+                    a.getLowerLeft().getY() - d),
-           Point<T>(a.max_corner().template get<0>() + d,
+           Point<T>(a.getUpperRight().getX() + d,
-                    a.max_corner().template get<1>() + d));
+                    a.getUpperRight().getY() + d));
   return get(b, s);
 }
 
@@ -192,17 +192,17 @@ std::set<std::pair<size_t, size_t> > Grid<V, G, T>::getCells(
 // _____________________________________________________________________________
 template <typename V, template <typename> typename G, typename T>
 Box<T> Grid<V, G, T>::getBox(size_t x, size_t y) const {
-  Point<T> sw(_bb.min_corner().template get<0>() + x * _cellWidth,
+  Point<T> sw(_bb.getLowerLeft().getX() + x * _cellWidth,
-              _bb.min_corner().template get<1>() + y * _cellHeight);
+              _bb.getLowerLeft().getY() + y * _cellHeight);
-  Point<T> ne(_bb.min_corner().template get<0>() + (x + 1) * _cellWidth,
+  Point<T> ne(_bb.getLowerLeft().getX() + (x + 1) * _cellWidth,
-              _bb.min_corner().template get<1>() + (y + 1) * _cellHeight);
+              _bb.getLowerLeft().getY() + (y + 1) * _cellHeight);
   return Box<T>(sw, ne);
 }
 
 // _____________________________________________________________________________
 template <typename V, template <typename> typename G, typename T>
 size_t Grid<V, G, T>::getCellXFromX(double x) const {
-  float dist = x - _bb.min_corner().template get<0>();
+  float dist = x - _bb.getLowerLeft().getX();
   if (dist < 0) dist = 0;
   return floor(dist / _cellWidth);
 }
@@ -210,7 +210,7 @@ size_t Grid<V, G, T>::getCellXFromX(double x) const {
 // _____________________________________________________________________________
 template <typename V, template <typename> typename G, typename T>
 size_t Grid<V, G, T>::getCellYFromY(double y) const {
-  float dist = y - _bb.min_corner().template get<1>();
+  float dist = y - _bb.getLowerLeft().getY();
   if (dist < 0) dist = 0;
   return floor(dist / _cellHeight);
 }

src/util/geo/Line.h

@@ -0,0 +1,28 @@
+// 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 geo {
+
+template <typename T>
+class Line : public std::vector<Point<T>> {
+  using std::vector<Point<T>>::vector;
+};
+
+template <typename T>
+using LineSegment = std::pair<Point<T>, Point<T>>;
+
+template <typename T>
+using MultiLine = std::vector<Line<T>>;
+
+}  // namespace geo
+}  // namespace util
+
+#endif  // UTIL_GEO_LINE_H_

src/util/geo/Point.h

@@ -0,0 +1,48 @@
+// Copyright 2016, University of Freiburg,
+// Chair of Algorithms and Data Structures.
+// Author: Patrick Brosi <brosi@informatik.uni-freiburg.de>
+
+#ifndef UTIL_GEO_POINT_H_
+#define UTIL_GEO_POINT_H_
+
+namespace util {
+namespace geo {
+
+template <typename T>
+class Point {
+ public:
+  Point() : _x(0), _y(0) {}
+  Point(T x, T y) : _x(x), _y(y) {}
+  T getX() const { return _x; }
+  T getY() const { return _y; }
+
+  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());
+  }
+
+  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;
+};
+
+template <typename T>
+using MultiPoint = std::vector<Point<T>>;
+
+}  // namespace geo
+}  // namespace util
+
+#endif  // UTIL_GEO_POINT_H_

src/util/geo/PolyLine.h

@@ -5,6 +5,7 @@
 #ifndef UTIL_GEO_POLYLINE_H_
 #define UTIL_GEO_POLYLINE_H_
 
+#include <cfloat>
 #include <ostream>
 #include <string>
 #include <vector>
@@ -16,6 +17,8 @@ namespace geo {
 static const double MAX_EQ_DISTANCE = 15;
 static const double AVERAGING_STEP = 20;
 
+// legacy code, will be removed in the future
+
 template <typename T>
 struct LinePoint {
   LinePoint() : lastIndex(0), totalPos(-1), p() {}
@@ -45,7 +48,6 @@ struct SharedSegments {
   std::vector<SharedSegment<T>> segments;
 };
 
-// TODO: maybe let this class inherit from a more generic geometry class
 template <typename T>
 class PolyLine {
  public:

src/util/geo/PolyLine.tpp

@@ -2,9 +2,6 @@
 // Chair of Algorithms and Data Structures.
 // Authors: Patrick Brosi <brosi@informatik.uni-freiburg.de>
 
-#include "util/geo/Geo.h"
-#include "util/geo/PolyLine.h"
-
 // _____________________________________________________________________________
 template <typename T>
 PolyLine<T>::PolyLine() {}
@@ -69,11 +66,6 @@ void PolyLine<T>::offsetPerp(double units) {
    *
    * there doesn't seem to be any library which reliably does that,
    * so we do it ourself here until we find one...
-   * boost::geometry only supports buffering a line, resulting in a
-   * polygon. An offsetted line is part of that polygon, but retrieving
-   * it reliably could result in some geometrical diffing hocus pocus which is
-   * bound to go wrong at /some/ point (self intersections, numerical
-   * instability etc)
    */
 
   if (fabs(units) < 0.001) return;
@@ -90,8 +82,8 @@ void PolyLine<T>::offsetPerp(double units) {
   for (size_t i = 1; i < _line.size(); i++) {
     Point<T> curP = _line[i];
 
-    double n1 = lastP.template get<1>() - curP.template get<1>();
+    double n1 = lastP.getY() - curP.getY();
-    double n2 = curP.template get<0>() - lastP.template get<0>();
+    double n2 = curP.getX() - lastP.getX();
     double n = sqrt(n1 * n1 + n2 * n2);
 
     // if n == 0, the segment is effectively a point
@@ -101,11 +93,11 @@ void PolyLine<T>::offsetPerp(double units) {
     n1 = n1 / n;
     n2 = n2 / n;
 
-    lastP.set<0>(lastP.template get<0>() + (n1 * units));
+    lastP.setX(lastP.getX() + (n1 * units));
-    lastP.set<1>(lastP.template get<1>() + (n2 * units));
+    lastP.setY(lastP.getY() + (n2 * units));
 
-    curP.set<0>(curP.template get<0>() + (n1 * units));
+    curP.setX(curP.getX() + (n1 * units));
-    curP.set<1>(curP.template get<1>() + (n2 * units));
+    curP.setY(curP.getY() + (n2 * units));
 
     if (lastIns && befLastIns &&
         lineIntersects(*lastIns, *befLastIns, lastP, curP)) {
@@ -216,7 +208,7 @@ LinePoint<T> PolyLine<T>::getPointAtDist(double atDist) const {
 
   for (size_t i = 1; i < _line.size(); i++) {
     const Point<T>& cur = _line[i];
-    double d = geo::dist(last, cur);
+    double d = geo::dist(*last, cur);
     dist += d;
 
     if (dist > atDist) {
@@ -242,13 +234,13 @@ LinePoint<T> PolyLine<T>::getPointAt(double at) const {
 template <typename T>
 Point<T> PolyLine<T>::interpolate(const Point<T>& a, const Point<T>& b,
                                   double p) const {
-  double n1 = b.template get<0>() - a.template get<0>();
+  double n1 = b.getX() - a.getX();
-  double n2 = b.template get<1>() - a.template get<1>();
+  double n2 = b.getY() - a.getY();
   double n = sqrt(n1 * n1 + n2 * n2);
   n1 = n1 / n;
   n2 = n2 / n;
-  return Point<T>(a.template get<0>() + (n1 * p),
+  return Point<T>(a.getX() + (n1 * p),
-                  a.template get<1>() + (n2 * p));
+                  a.getY() + (n2 * p));
 }
 
 // _____________________________________________________________________________
@@ -307,11 +299,11 @@ PolyLine<T> PolyLine<T>::average(const std::vector<const PolyLine<T>*>& lines,
       const PolyLine* pl = lines[i];
       Point<T> p = pl->getPointAt(a).p;
       if (weighted) {
-        x += p.template get<0>() * weights[i];
+        x += p.getX() * weights[i];
-        y += p.template get<1>() * weights[i];
+        y += p.getY() * weights[i];
       } else {
-        x += p.template get<0>();
+        x += p.getX();
-        y += p.template get<1>();
+        y += p.getY();
       }
     }
     ret << Point<T>(x / total, y / total);
@@ -468,8 +460,8 @@ bool PolyLine<T>::contains(const PolyLine<T>& rhs, double dmax) const {
 template <typename T>
 void PolyLine<T>::move(double vx, double vy) {
   for (size_t i = 0; i < _line.size(); i++) {
-    _line[i].set<0>(_line[i].template get<0>() + vx);
+    _line[i].setX(_line[i].getX() + vx);
-    _line[i].set<1>(_line[i].template get<1>() + vy);
+    _line[i].setY(_line[i].getY() + vy);
   }
 }
 
@@ -627,11 +619,11 @@ PolyLine<T> PolyLine<T>::getOrthoLineAtDist(double d, double length) const {
 
   double angle = angBetween(getPointAtDist(d - 5).p, getPointAtDist(d + 5).p);
 
-  double angleX1 = avgP.template get<0>() + cos(angle + M_PI / 2) * length / 2;
+  double angleX1 = avgP.getX() + cos(angle + M_PI / 2) * length / 2;
-  double angleY1 = avgP.template get<1>() + sin(angle + M_PI / 2) * length / 2;
+  double angleY1 = avgP.getY() + sin(angle + M_PI / 2) * length / 2;
 
-  double angleX2 = avgP.template get<0>() + cos(angle + M_PI / 2) * -length / 2;
+  double angleX2 = avgP.getX() + cos(angle + M_PI / 2) * -length / 2;
-  double angleY2 = avgP.template get<1>() + sin(angle + M_PI / 2) * -length / 2;
+  double angleY2 = avgP.getY() + sin(angle + M_PI / 2) * -length / 2;
 
   return PolyLine(Point<T>(angleX1, angleY1), Point<T>(angleX2, angleY2));
 }
@@ -651,8 +643,8 @@ std::pair<double, double> PolyLine<T>::getSlopeBetween(double ad,
 
   double d = dist(a.p, b.p);
 
-  double dx = (b.p.template get<0>() - a.p.template get<0>()) / d;
+  double dx = (b.p.getX() - a.p.getX()) / d;
-  double dy = (b.p.template get<1>() - a.p.template get<1>()) / d;
+  double dy = (b.p.getY() - a.p.getY()) / d;
 
   return std::pair<double, double>(dx, dy);
 }
@@ -714,11 +706,11 @@ void PolyLine<T>::applyChaikinSmooth(size_t depth) {
       Point<T> pB = _line[i];
 
       smooth.push_back(
-          Point<T>(0.75 * pA.template get<0>() + 0.25 * pB.template get<0>(),
+          Point<T>(0.75 * pA.getX() + 0.25 * pB.getX(),
-                   0.75 * pA.template get<1>() + 0.25 * pB.template get<1>()));
+                   0.75 * pA.getY() + 0.25 * pB.getY()));
       smooth.push_back(
-          Point<T>(0.25 * pA.template get<0>() + 0.75 * pB.template get<0>(),
+          Point<T>(0.25 * pA.getX() + 0.75 * pB.getX(),
-                   0.25 * pA.template get<1>() + 0.75 * pB.template get<1>()));
+                   0.25 * pA.getY() + 0.75 * pB.getY()));
     }
 
     smooth.push_back(_line.back());

src/util/geo/Polygon.h

@@ -0,0 +1,41 @@
+// Copyright 2016, University of Freiburg,
+// Chair of Algorithms and Data Structures.
+// Author: Patrick Brosi <brosi@informatik.uni-freiburg.de>
+
+#ifndef UTIL_GEO_POLYGON_H_
+#define UTIL_GEO_POLYGON_H_
+
+#include <vector>
+#include "./Box.h"
+#include "./Line.h"
+#include "./Point.h"
+
+namespace util {
+namespace geo {
+
+template <typename T>
+class Polygon {
+ public:
+  Polygon() {}
+
+  Polygon(const Line<T>& l) : _outer(l) {}
+  Polygon(const Box<T>& b)
+      : _outer({b.getLowerLeft(),
+                Point<T>(b.getUpperRight().getX(), b.getLowerLeft().getY()),
+                b.getUpperRight(),
+                Point<T>(b.getLowerLeft().getX(), b.getUpperRight().getY())}) {}
+
+  const Line<T>& getOuter() const { return _outer; }
+  Line<T>& getOuter() { return _outer; }
+
+ private:
+  Line<T> _outer;
+};
+
+template <typename T>
+using MultiPolyon = std::vector<Polygon<T>>;
+
+}  // namespace geo
+}  // namespace util
+
+#endif  // UTIL_GEO_LINE_H_

src/util/geo/output/GeoJsonOutput.tpp

@@ -13,8 +13,8 @@ void GeoJsonOutput::print(const Point<T>& p, Attrs attrs) {
   _wr.keyVal("type", "Point");
   _wr.key("coordinates");
   _wr.arr();
-  _wr.val(p.template get<0>());
+  _wr.val(p.getX());
-  _wr.val(p.template get<1>());
+  _wr.val(p.getY());
   _wr.close();
   _wr.close();
   _wr.key("properties");
@@ -36,8 +36,8 @@ void GeoJsonOutput::print(const Line<T>& line, Attrs attrs) {
   _wr.arr();
   for (auto p : line) {
     _wr.arr();
-    _wr.val(p.template get<0>());
+    _wr.val(p.getX());
-    _wr.val(p.template get<1>());
+    _wr.val(p.template getY());
     _wr.close();
   }
   _wr.close();

src/util/graph/Graph.h

@@ -8,6 +8,7 @@
 #include <set>
 #include <string>
 #include <iostream>
+#include <cassert>
 
 #include "util/graph/Edge.h"
 #include "util/graph/Node.h"
@@ -18,7 +19,7 @@ namespace graph {
 template <typename N, typename E>
 class Graph {
  public:
-  ~Graph();
+  virtual ~Graph();
   virtual Node<N, E>* addNd() = 0;
   virtual Node<N, E>* addNd(const N& pl) = 0;
   Edge<N, E>* addEdg(Node<N, E>* from, Node<N, E>* to);

src/util/tests/TestMain.cpp

@@ -129,8 +129,8 @@ CASE("grid") {
   l.push_back(Point<double>(1.5, 2));
 
   Line<double> l2;
-  l.push_back(Point<double>(2.5, 1));
+  l2.push_back(Point<double>(2.5, 1));
-  l.push_back(Point<double>(2.5, 2));
+  l2.push_back(Point<double>(2.5, 2));
 
   g.add(l, 1);
   g.add(l2, 2);
@@ -143,8 +143,7 @@ CASE("grid") {
 
   ret.clear();
   g.getNeighbors(1, 0, &ret);
-  // TODO!
+  EXPECT(ret.size() == (size_t)1);
-  //EXPECT(ret.size() == 1);
 
   ret.clear();
   g.getNeighbors(1, 0.55, &ret);
@@ -165,7 +164,7 @@ CASE("densify") {
   EXPECT(dense.size() == (size_t)10);
 
   for (int i = 0; i < 10; i++) {
-    EXPECT(dense[i].get<0>() == approx(i + 1.0));
+    EXPECT(dense[i].getX() == approx(i + 1.0));
   }
 
   dense = util::geo::simplify(dense, 0.1);
@@ -841,6 +840,352 @@ CASE("nullable") {
 
     EXPECT_THROWS(nullable == voidnull);
   }
+},
+
+// ___________________________________________________________________________
+CASE("geometry") {
+  geo::Point<double> a(1, 2);
+  geo::Point<double> b(2, 3);
+  geo::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(geo::dist(a, a) == approx(0));
+  EXPECT(geo::dist(a, b) == approx(sqrt(2)));
+
+  d = d + d;
+
+  geo::Box<double> box(a, c);
+  EXPECT(geo::contains(a, box));
+  EXPECT(geo::contains(b, box));
+  EXPECT(geo::contains(c, box));
+  EXPECT(!geo::contains(d, box));
+
+  geo::Line<double> line{a, b, c};
+
+  EXPECT(geo::contains(line, box));
+  line.push_back(d);
+  EXPECT(!geo::contains(line, box));
+
+  geo::LineSegment<double> ls{a, b};
+  EXPECT(geo::contains(a, ls));
+  EXPECT(geo::contains(b, ls));
+  EXPECT(!geo::contains(c, ls));
+  EXPECT(geo::contains(a + geo::Point<double>(.5, .5), ls));
+  EXPECT(!geo::contains(a + geo::Point<double>(1.5, 1.5), ls));
+
+  geo::LineSegment<double> lsa{geo::Point<double>(1, 1), geo::Point<double>(2, 2)};
+  geo::LineSegment<double> lsb{geo::Point<double>(1, 2), geo::Point<double>(2, 1)};
+  geo::LineSegment<double> lsc{geo::Point<double>(2.1, 2), geo::Point<double>(3, 3)};
+
+  EXPECT(geo::crossProd(lsa.first, lsb) == approx(-1));
+  EXPECT(geo::crossProd(lsa.second, lsb) == approx(1));
+
+  EXPECT(geo::intersects(lsa, lsb));
+
+  EXPECT(geo::intersects(lsa, lsa));
+  EXPECT(geo::intersects(lsb, lsb));
+  EXPECT(!geo::intersects(lsa, lsc));
+
+  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));
+  EXPECT(geo::contains(geo::Point<double>(2, 3), l));
+
+  geo::Box<double> bbox(geo::Point<double>(1, 1), geo::Point<double>(3, 3));
+  EXPECT(geo::intersects(l, bbox));
+  geo::Line<double> ll{geo::Point<double>(0, 0), geo::Point<double>(4, 4)};
+  EXPECT(geo::intersects(ll, bbox));
+  geo::Line<double> lll{geo::Point<double>(0, 0), geo::Point<double>(0, 4)};
+  EXPECT(!geo::intersects(lll, bbox));
+  geo::Line<double> llll{geo::Point<double>(1.2, 0), geo::Point<double>(1, 2)};
+  EXPECT(geo::intersects(llll, bbox));
+
+  Line<double> l5;
+  l5.push_back(Point<double>(0, 0));
+  l5.push_back(Point<double>(1.5, 2));
+  Box<double> req(Point<double>(.5, 1), Point<double>(1, 1.5));
+
+  EXPECT(geo::getBoundingBox(l5[0]).getLowerLeft().getX() == approx(0));
+  EXPECT(geo::getBoundingBox(l5[0]).getLowerLeft().getY() == approx(0));
+
+  EXPECT(geo::getBoundingBox(l5).getLowerLeft().getX() == approx(0));
+  EXPECT(geo::getBoundingBox(l5).getLowerLeft().getY() == approx(0));
+  EXPECT(geo::getBoundingBox(l5).getUpperRight().getX() == approx(1.5));
+  EXPECT(geo::getBoundingBox(l5).getUpperRight().getY() == approx(2));
+  EXPECT(geo::intersects(geo::getBoundingBox(l5), geo::getBoundingBox(Line<double>{Point<double>(.5, 1), Point<double>(1, 1)})));
+  EXPECT(geo::intersects(l5, Line<double>{Point<double>(.5, 1), Point<double>(1, 1)}));
+  EXPECT(geo::intersects(l5, req));
+
+  Box<double> boxa(Point<double>(1, 1), Point<double>(2, 2));
+  EXPECT(geo::intersects(boxa, Box<double>(Point<double>(1.5, 1.5), Point<double>(1.7, 1.7))));
+  EXPECT(geo::intersects(boxa, Box<double>(Point<double>(0, 0), Point<double>(3, 3))));
+  EXPECT(geo::intersects(boxa, Box<double>(Point<double>(1.5, 1.5), Point<double>(3, 3))));
+  EXPECT(geo::intersects(boxa, Box<double>(Point<double>(0, 0), Point<double>(1.5, 1.5))));
+
+  EXPECT(geo::intersects(Box<double>(Point<double>(1.5, 1.5), Point<double>(1.7, 1.7)), boxa));
+  EXPECT(geo::intersects(Box<double>(Point<double>(0, 0), Point<double>(3, 3)), boxa));
+  EXPECT(geo::intersects(Box<double>(Point<double>(1.5, 1.5), Point<double>(3, 3)), boxa));
+  EXPECT(geo::intersects(Box<double>(Point<double>(0, 0), Point<double>(1.5, 1.5)), boxa));
+
+  Polygon<double> poly({Point<double>(1, 1), Point<double>(3, 2), Point<double>(4, 3), Point<double>(6, 3), Point<double>(5, 1)});
+  EXPECT(geo::getWKT(poly) == "POLYGON ((1 1, 3 2, 4 3, 6 3, 5 1, 1 1))");
+  EXPECT(geo::contains(Point<double>(4, 2), poly));
+  EXPECT(!geo::contains(Point<double>(3, 3), poly));
+  EXPECT(geo::contains(Point<double>(1, 1), poly));
+  EXPECT(geo::contains(Point<double>(3, 2), poly));
+  EXPECT(geo::contains(Point<double>(4, 3), poly));
+  EXPECT(geo::contains(Point<double>(6, 3), poly));
+  EXPECT(geo::contains(Point<double>(5, 1), poly));
+
+  EXPECT(geo::contains(Line<double>{Point<double>(6, 3), Point<double>(5, 1)}, poly));
+  EXPECT(!geo::contains(Line<double>{Point<double>(6, 3), Point<double>(50, 1)}, poly));
+  EXPECT(geo::contains(Line<double>{Point<double>(4, 2), Point<double>(4.5, 2)}, poly));
+  EXPECT(geo::contains(Line<double>{Point<double>(4, 2), Point<double>(5, 1)}, poly));
+
+  Box<double> polybox(Point<double>(1, 1), Point<double>(6, 4));
+  EXPECT(geo::centroid(polybox).getX() == approx(3.5));
+  EXPECT(geo::centroid(polybox).getY() == approx(2.5));
+  EXPECT(geo::contains(poly, polybox));
+  EXPECT(!geo::contains(polybox, poly));
+  Box<double> polybox2(Point<double>(4, 1), Point<double>(5, 2));
+  EXPECT(geo::contains(polybox2, poly));
+  EXPECT(geo::contains(poly, getBoundingBox(poly)));
+
+  Point<double> rotP(2, 2);
+  EXPECT(geo::dist(geo::rotate(rotP, 180, Point<double>(1, 1)), Point<double>(0, 0)) == approx(0));
+  EXPECT(geo::dist(geo::rotate(rotP, 360, Point<double>(1, 1)), rotP) == approx(0));
+
+  Line<double> rotLine({{1, 1}, {3, 3}});
+  EXPECT(geo::rotate(rotLine, 90, Point<double>(2, 2))[0].getX() == approx(1));
+  EXPECT(geo::rotate(rotLine, 90, Point<double>(2, 2))[0].getY() == approx(3));
+  EXPECT(geo::rotate(rotLine, 90, Point<double>(2, 2))[1].getX() == approx(3));
+  EXPECT(geo::rotate(rotLine, 90, Point<double>(2, 2))[1].getY() == approx(1));
+
+  MultiLine<double> multiRotLine({{{1, 1}, {3, 3}}, {{1, 3}, {3, 1}}});
+  EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[0][0].getX() == approx(1));
+  EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[0][0].getY() == approx(3));
+  EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[0][1].getX() == approx(3));
+  EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[0][1].getY() == approx(1));
+  EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[1][0].getX() == approx(3));
+  EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[1][0].getY() == approx(3));
+  EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[1][1].getX() == approx(1));
+  EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[1][1].getY() == approx(1));
+
+  EXPECT(geo::getWKT(multiRotLine) == "MULTILINESTRING ((1 1, 3 3), (1 3, 3 1))");
+
+  EXPECT(geo::contains(multiRotLine[0], geo::move(geo::move(multiRotLine, 1.0, 2.0), -1.0, -2.0)[0]));
+  EXPECT(geo::contains(multiRotLine, geo::getBoundingBox(Line<double>{{1, 1}, {3, 3}, {1, 3}, {3, 1}})));
+
+  EXPECT(geo::contains(getBoundingBox(multiRotLine), geo::getBoundingBox(Line<double>{{1, 1}, {3, 3}, {1, 3}, {3, 1}})));
+  EXPECT(geo::contains(geo::getBoundingBox(Line<double>{{1, 1}, {3, 3}, {1, 3}, {3, 1}}), getBoundingBox(multiRotLine)));
+
+  EXPECT(geo::dist(geo::centroid(rotP), rotP) == approx(0));
+  EXPECT(geo::dist(geo::centroid(rotLine), rotP) == approx(0));
+  EXPECT(geo::dist(geo::centroid(polybox), Point<double>(3.5, 2.5)) == approx(0));
+  EXPECT(geo::dist(geo::centroid(Polygon<double>({{0, 0}, {3, 4}, {4,3}})), Point<double>(7.0/3.0,7.0/3.0)) == approx(0));
+
+  auto polyy = Polygon<double>({{0, 0}, {3, 4}, {4,3}});
+  MultiPolyon<double> mpoly{polyy, polyy};
+
+  EXPECT(geo::getWKT(polyy) == "POLYGON ((0 0, 3 4, 4 3, 0 0))");
+  EXPECT(geo::getWKT(mpoly) == "MULTIPOLYGON (((0 0, 3 4, 4 3, 0 0)), ((0 0, 3 4, 4 3, 0 0)))");
+
+  auto hull = geo::convexHull(Line<double>{{0.1, 3}, {1, 1}, {2, 2}, {4, 4}, {0, 0}, {1, 2}, {3, 1}, {3, 3}});
+  EXPECT(hull.getOuter().size() == size_t(4));
+  EXPECT(hull.getOuter()[0].getX() == approx(0));
+  EXPECT(hull.getOuter()[0].getY() == approx(0));
+  EXPECT(hull.getOuter()[1].getX() == approx(3));
+  EXPECT(hull.getOuter()[1].getY() == approx(1));
+  EXPECT(hull.getOuter()[2].getX() == approx(4));
+  EXPECT(hull.getOuter()[2].getY() == approx(4));
+  EXPECT(hull.getOuter()[3].getX() == approx(0.1));
+  EXPECT(hull.getOuter()[3].getY() == approx(3));
+  EXPECT(geo::contains(geo::convexHull(geo::getBoundingBox(poly)), geo::getBoundingBox(poly)));
+  EXPECT(geo::contains(geo::getBoundingBox(poly), geo::convexHull(geo::getBoundingBox(poly))));
+
+  auto hull2 = geo::convexHull(Line<double>{{0.1, 3}, {1, 1}, {2, 2}, {4, 4}, {0, 0}, {1, 2}, {3, 1}, {3, 3}, {-0.1, 1}});
+  EXPECT(hull2.getOuter().size() == size_t(5));
+  EXPECT(hull2.getOuter()[0].getX() == approx(-.1));
+  EXPECT(hull2.getOuter()[0].getY() == approx(1));
+  EXPECT(hull2.getOuter()[1].getX() == approx(0));
+  EXPECT(hull2.getOuter()[1].getY() == approx(0));
+  EXPECT(hull2.getOuter()[2].getX() == approx(3));
+  EXPECT(hull2.getOuter()[2].getY() == approx(1));
+  EXPECT(hull2.getOuter()[3].getX() == approx(4));
+  EXPECT(hull2.getOuter()[3].getY() == approx(4));
+  EXPECT(hull2.getOuter()[4].getX() == approx(0.1));
+  EXPECT(hull2.getOuter()[4].getY() == approx(3));
+
+  auto hull3 = geo::convexHull(Line<double>{{0.1, 3}, {4, 4}, {0, 0}, {1, 2}, {3, 1}});
+  EXPECT(hull3.getOuter().size() == size_t(4));
+  EXPECT(hull3.getOuter()[0].getX() == approx(0));
+  EXPECT(hull3.getOuter()[0].getY() == approx(0));
+  EXPECT(hull3.getOuter()[1].getX() == approx(3));
+  EXPECT(hull3.getOuter()[1].getY() == approx(1));
+  EXPECT(hull3.getOuter()[2].getX() == approx(4));
+  EXPECT(hull3.getOuter()[2].getY() == approx(4));
+  EXPECT(hull3.getOuter()[3].getX() == approx(0.1));
+  EXPECT(hull3.getOuter()[3].getY() == approx(3));
+
+  hull3 = geo::convexHull(Line<double>{{0.1, 3}, {4, 4}, {2, 1}, {3, 2}, {0, 0}, {1, 2}, {3, 1}});
+  EXPECT(hull3.getOuter().size() == size_t(4));
+  EXPECT(hull3.getOuter()[0].getX() == approx(0));
+  EXPECT(hull3.getOuter()[0].getY() == approx(0));
+  EXPECT(hull3.getOuter()[1].getX() == approx(3));
+  EXPECT(hull3.getOuter()[1].getY() == approx(1));
+  EXPECT(hull3.getOuter()[2].getX() == approx(4));
+  EXPECT(hull3.getOuter()[2].getY() == approx(4));
+  EXPECT(hull3.getOuter()[3].getX() == approx(0.1));
+  EXPECT(hull3.getOuter()[3].getY() == approx(3));
+
+  hull3 = geo::convexHull(Line<double>{{4, 4}, {1, 2}, {2, 1}, {3, 2}, {0.1, 3}, {0, 0}, {1, 2}, {3, 1}});
+  EXPECT(hull3.getOuter().size() == size_t(4));
+  EXPECT(hull3.getOuter()[0].getX() == approx(0));
+  EXPECT(hull3.getOuter()[0].getY() == approx(0));
+  EXPECT(hull3.getOuter()[1].getX() == approx(3));
+  EXPECT(hull3.getOuter()[1].getY() == approx(1));
+  EXPECT(hull3.getOuter()[2].getX() == approx(4));
+  EXPECT(hull3.getOuter()[2].getY() == approx(4));
+  EXPECT(hull3.getOuter()[3].getX() == approx(0.1));
+  EXPECT(hull3.getOuter()[3].getY() == approx(3));
+
+  hull3 = geo::convexHull(Line<double>{{4, 4}, {1, 2}, {3, 1}});
+  EXPECT(hull3.getOuter().size() == size_t(3));
+  EXPECT(hull3.getOuter()[0].getX() == approx(1));
+  EXPECT(hull3.getOuter()[0].getY() == approx(2));
+  EXPECT(hull3.getOuter()[1].getX() == approx(3));
+  EXPECT(hull3.getOuter()[1].getY() == approx(1));
+  EXPECT(hull3.getOuter()[2].getX() == approx(4));
+  EXPECT(hull3.getOuter()[2].getY() == approx(4));
+
+  hull3 = geo::convexHull(Line<double>{{4, 4}, {1, 2}, {3, 10}});
+  EXPECT(hull3.getOuter().size() == size_t(3));
+  EXPECT(hull3.getOuter()[0].getX() == approx(1));
+  EXPECT(hull3.getOuter()[0].getY() == approx(2));
+  EXPECT(hull3.getOuter()[1].getX() == approx(4));
+  EXPECT(hull3.getOuter()[1].getY() == approx(4));
+  EXPECT(hull3.getOuter()[2].getX() == approx(3));
+  EXPECT(hull3.getOuter()[2].getY() == approx(10));
+
+  Line<double> test{{0.3215348546593775, 0.03629583077160248},
+                    {0.02402358131857918, -0.2356728797179394},
+                    {0.04590851212470659, -0.4156409924995536},
+                    {0.3218384001607433, 0.1379850698988746},
+                    {0.11506479756447, -0.1059521474930943},
+                    {0.2622539999543261, -0.29702873322836},
+                    {-0.161920957418085, -0.4055339716426413},
+                    {0.1905378631228002, 0.3698601009043493},
+                    {0.2387090918968516, -0.01629827079949742},
+                    {0.07495888748668034, -0.1659825110491202},
+                    {0.3319341836794598, -0.1821814101954749},
+                    {0.07703635755650362, -0.2499430638271785},
+                    {0.2069242999022122, -0.2232970760420869},
+                    {0.04604079532068295, -0.1923573186549892},
+                    {0.05054295812784038, 0.4754929463150845},
+                    {-0.3900589168910486, 0.2797829520700341},
+                    {0.3120693385713448, -0.0506329867529059},
+                    {0.01138812723698857, 0.4002504701728471},
+                    {0.009645149586391732, 0.1060251100976254},
+                    {-0.03597933197019559, 0.2953639456959105},
+                    {0.1818290866742182, 0.001454397571696298},
+                    {0.444056063372694, 0.2502497166863175},
+                    {-0.05301752458607545, -0.06553921621808712},
+                    {0.4823896228171788, -0.4776170002088109},
+                    {-0.3089226845734964, -0.06356112199235814},
+                    {-0.271780741188471, 0.1810810595574612},
+                    {0.4293626522918815, 0.2980897964891882},
+                    {-0.004796652127799228, 0.382663812844701},
+                    {0.430695573269106, -0.2995073500084759},
+                    {0.1799668387323309, -0.2973467472915973},
+                    {0.4932166845474547, 0.4928094162538735},
+                    {-0.3521487911717489, 0.4352656197131292},
+                    {-0.4907368011686362, 0.1865826865533206},
+                    {-0.1047924716070224, -0.247073392148198},
+                    {0.4374961861758457, -0.001606279519951237},
+                    {0.003256207800708899, -0.2729194320486108},
+                    {0.04310378203457577, 0.4452604050238248},
+                    {0.4916198379282093, -0.345391701297268},
+                    {0.001675087028811806, 0.1531837672490476},
+                    {-0.4404289572876217, -0.2894855991839297}
+
+  };
+  hull3 = geo::convexHull(test);
+  EXPECT(geo::contains(test, hull3));
+  EXPECT(hull3.getOuter().size() == size_t(8));
+  EXPECT(geo::contains(Polygon<double>({{-0.161920957418085, -0.4055339716426413},
+                                        {0.05054295812784038, 0.4754929463150845},
+                                        {0.4823896228171788, -0.4776170002088109},
+                                        {0.4932166845474547, 0.4928094162538735},
+                                        {-0.3521487911717489, 0.4352656197131292},
+                                        {-0.4907368011686362, 0.1865826865533206},
+                                        {0.4916198379282093, -0.345391701297268},
+                                        {-0.4404289572876217,
+                                         -0.2894855991839297}}), hull3));
+  EXPECT(geo::contains(hull3, Polygon<double>({{-0.161920957418085, -0.4055339716426413},
+                                        {0.05054295812784038, 0.4754929463150845},
+                                        {0.4823896228171788, -0.4776170002088109},
+                                        {0.4932166845474547, 0.4928094162538735},
+                                        {-0.3521487911717489, 0.4352656197131292},
+                                        {-0.4907368011686362, 0.1865826865533206},
+                                        {0.4916198379282093, -0.345391701297268},
+                                        {-0.4404289572876217,
+                                         -0.2894855991839297}})));
+
+  hull3 = geo::convexHull(Line<double>{{3, 6}, {8, 10}, {3, 5}, {20, -10}, {-4, 5}, {10, 2}, {5, 1}, {45, 1}, {30, -9}, {3, 14}, {25, -5.5}});
+  EXPECT(hull3.getOuter().size() == size_t(5));
+  EXPECT(hull3.getOuter()[0].getX() == approx(-4));
+  EXPECT(hull3.getOuter()[0].getY() == approx(5));
+  EXPECT(hull3.getOuter()[1].getX() == approx(20));
+  EXPECT(hull3.getOuter()[1].getY() == approx(-10));
+  EXPECT(hull3.getOuter()[2].getX() == approx(30));
+  EXPECT(hull3.getOuter()[2].getY() == approx(-9));
+  EXPECT(hull3.getOuter()[3].getX() == approx(45));
+  EXPECT(hull3.getOuter()[3].getY() == approx(1));
+  EXPECT(hull3.getOuter()[4].getX() == approx(3));
+  EXPECT(hull3.getOuter()[4].getY() == approx(14));
+
+  hull3 = geo::convexHull(Line<double>{{7, 7}, {7, -7}, {-7, -7}, {-7, 7}, {9, 0}, {-9, 0}, {0, 9}, {0, -9}});
+  EXPECT(hull3.getOuter().size() == size_t(8));
+  EXPECT(geo::contains(geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}}), hull3));
+  EXPECT(geo::contains(hull3, geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}})));
+
+  hull3 = geo::convexHull(Line<double>{{7, 7}, {7, -7}, {-7, -7}, {-7, 7}, {9, 0}, {-9, 0}, {0, 9}, {0, -9}, {0, 0}, {1, 2}, {-2, 1}, {-1, -1}, {3, 4}, {4, 3}, {-5, 4}, {6, 5}});
+  EXPECT(hull3.getOuter().size() == size_t(8));
+  EXPECT(geo::contains(geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}}), hull3));
+  EXPECT(geo::contains(hull3, geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}})));
+
+  hull3 = geo::convexHull(Line<double>{{0, 0}, {1, 2}, {-2, 1}, {-1, -1}, {3, 4}, {4, 3}, {-5, 4}, {6, 5}, {7, 7}, {7, -7}, {-7, -7}, {-7, 7}, {9, 0}, {-9, 0}, {0, 9}, {0, -9}, {-8, 0}, {8, 0}, {-7, 0}, {7, 0}, {-6, 0}, {6, 0}, {-5, 0}, {5, 0}, {-4, 0}, {4, 0}, {-3, 0}, {3, 0}, {-2, 0}, {2, 0}, {-1, 0}, {1, 0}, {0, -8}, {0, 8}, {0, -7}, {0, 7}, {0, -6}, {0, 6}, {0, -5}, {0, 5}, {0, -4}, {0, 4}, {0, -3}, {0, 3}, {0, -2}, {0, 2}, {0, -1}, {0, 1}, {1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {1, -1}, {2, -2}, {3, -3}, {4, -4},  {5, -5}, {6, -6}, {-1, 1}, {-2, 2}, {-3, 3}, {-4, 4}, {-5, 5}, {-6, 6}, {-1, -1}, {-2, -2}, {-3, -3}, {-4, -4}, {-5, -5}, {-6, -6}});
+  EXPECT(hull3.getOuter().size() == size_t(8));
+  EXPECT(geo::contains(geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}}), hull3));
+  EXPECT(geo::contains(hull3, geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}})));
+
+  EXPECT(geo::area(geo::Point<double>(1, 2)) == approx(0));
+  EXPECT(geo::area(geo::Line<double>{{1, 2}, {2, 5}}) == approx(0));
+  EXPECT(geo::area(geo::Box<double>({0, 0}, {1, 1})) == approx(1));
+  EXPECT(geo::area(geo::Box<double>({1, 1}, {1, 1})) == approx(0));
+  EXPECT(geo::area(geo::Box<double>({0, 0}, {2, 2})) == approx(4));
+  EXPECT(geo::area(geo::Polygon<double>({{0, 0}, {1, 0}, {1, 1}, {0, 1}})) == approx(1));
+  EXPECT(geo::area(geo::Polygon<double>({{0, 0}, {1, 0}, {1, 1}})) == approx(0.5));
+
+  auto obox = geo::getOrientedEnvelope(geo::Line<double>{{0, 0}, {1, 1}, {1.5, 0.5}});
+  EXPECT(geo::contains(geo::convexHull(obox), geo::Polygon<double>({{0.0, 0.0}, {1.0, 1.0}, {1.5, 0.5}, {0.5, -0.5}})));
+  EXPECT(geo::contains(geo::Polygon<double>({{0.0, 0.0}, {1.0, 1.0}, {1.5, 0.5}, {0.5, -0.5}}), geo::convexHull(obox)));
 }
 
 };