diff --git a/README.md b/README.md
index b920872..f583171 100644
--- a/README.md
+++ b/README.md
@@ -1,15 +1,4 @@
-```
-                 _-====-__-=-__-===-__-=======-__
-               _(                               _)
-            OO(                                )
-         . o  '===-______-===-____-==-__-====='
-      .o
-     . ______          _______________
-   _()_||__|| __o^o___ | [] [] [] [] |
-  (           |      | |             |o
- /-OO----OO""="OO--OO"="OO---------OO"
-############################################################
-```
+![Map-Matched path of a single train through Switzerland](geo/schweiz_mmatched.png?raw=true)
 
 # pfaedle
 
@@ -18,7 +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)
 
 ## Building and Installation
 
@@ -49,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).
 
diff --git a/geo/schweiz_mmatched.png b/geo/schweiz_mmatched.png
new file mode 100644
index 0000000..d2c6b17
Binary files /dev/null and b/geo/schweiz_mmatched.png differ
diff --git a/src/pfaedle/PfaedleMain.cpp b/src/pfaedle/PfaedleMain.cpp
index c7d3960..2c54028 100644
--- a/src/pfaedle/PfaedleMain.cpp
+++ b/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);
   }
diff --git a/src/pfaedle/config/PfaedleConfig.h b/src/pfaedle/config/PfaedleConfig.h
index 74047b8..ba9a2a0 100644
--- a/src/pfaedle/config/PfaedleConfig.h
+++ b/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),
diff --git a/src/util/geo/Box.h b/src/util/geo/Box.h
index dfaae15..f2e77de 100644
--- a/src/util/geo/Box.h
+++ b/src/util/geo/Box.h
@@ -41,10 +41,41 @@ class Box {
 template <typename T>
 class RotatedBox {
  public:
-  RotatedBox() : _box() {}
-  RotatedBox(const Box<T>& box) : _box(box), _deg(0), _center() {}
+  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() {}
+      : _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;
diff --git a/src/util/geo/Geo.h b/src/util/geo/Geo.h
index f823296..a020c55 100644
--- a/src/util/geo/Geo.h
+++ b/src/util/geo/Geo.h
@@ -8,6 +8,8 @@
 #define _USE_MATH_DEFINES
 
 #include <math.h>
+#include <algorithm>
+#include <cassert>
 #include <iostream>
 #include <sstream>
 #include "util/Misc.h"
@@ -45,7 +47,7 @@ typedef Polygon<double> DPolygon;
 typedef Polygon<float> FPolygon;
 typedef Polygon<int> IPolygon;
 
-const static double EPSILON = 0.0000001;
+const static double EPSILON = 0.00000000001;
 const static double RAD = 0.017453292519943295;  // PI/180
 
 // _____________________________________________________________________________
@@ -93,6 +95,14 @@ inline Point<T> centroid(const Box<T> box) {
   return centroid(LineSegment<T>(box.getLowerLeft(), box.getUpperRight()));
 }
 
+// _____________________________________________________________________________
+template <typename T, template <typename> typename Geometry>
+inline Point<T> centroid(std::vector<Geometry<T>> multigeo) {
+  Line<T> a;
+  for (const auto& g : multigeo) a.push_back(centroid(g));
+  return centroid(a);
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline Point<T> rotate(const Point<T>& p, double deg) {
@@ -142,12 +152,63 @@ inline Polygon<T> rotate(Polygon<T> geo, double deg, const Point<T>& c) {
   return geo;
 }
 
+// _____________________________________________________________________________
+template <typename T, template <typename> typename Geometry>
+inline std::vector<Geometry<T>> rotate(std::vector<Geometry<T>> multigeo,
+                                       double deg, const Point<T>& c) {
+  for (size_t i = 0; i < multigeo.size(); i++)
+    multigeo[i] = rotate(multigeo[i], deg, c);
+  return multigeo;
+}
+
+// _____________________________________________________________________________
+template <typename T, template <typename> typename Geometry>
+inline std::vector<Geometry<T>> rotate(std::vector<Geometry<T>> multigeo,
+                                       double deg) {
+  auto c = centroid(multigeo);
+  for (size_t i = 0; i < multigeo.size(); i++)
+    multigeo[i] = rotate(multigeo[i], deg, c);
+  return multigeo;
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline Point<T> move(const Point<T>& geo, T x, T y) {
   return Point<T>(geo.getX() + x, geo.getY() + y);
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline Line<T> move(Line<T> geo, T x, T y) {
+  for (size_t i = 0; i < geo.size(); i++) geo[i] = move(geo[i], x, y);
+  return geo;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline LineSegment<T> move(LineSegment<T> geo, T x, T y) {
+  geo.first = move(geo.first, x, y);
+  geo.second = move(geo.second, x, y);
+  return geo;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> move(Polygon<T> geo, T x, T y) {
+  for (size_t i = 0; i < geo.getOuter().size(); i++)
+    geo.getOuter()[i] = move(geo.getOuter()[i], x, y);
+  return geo;
+}
+
+// _____________________________________________________________________________
+template <typename T, template <typename> typename Geometry>
+inline std::vector<Geometry<T>> move(std::vector<Geometry<T>> multigeo, T x,
+                                     T y) {
+  for (size_t i = 0; i < multigeo.size(); i++)
+    multigeo[i] = move(multigeo[i], x, y);
+  return multigeo;
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline Box<T> minbox() {
@@ -155,23 +216,20 @@ inline Box<T> minbox() {
 }
 
 // _____________________________________________________________________________
-// template <typename T>
-// inline RotatedBox<T> shrink(const RotatedBox<T>& b, double d) {
-// double xd =
-// b.b.getUpperRight().getX() - b.b.getLowerLeft().getX();
-// double yd =
-// b.b.getUpperRight().getY() - b.b.getLowerLeft().getY();
+template <typename T>
+inline RotatedBox<T> shrink(const RotatedBox<T>& b, double d) {
+  double xd = b.b.getUpperRight().getX() - b.b.getLowerLeft().getX();
+  double yd = b.b.getUpperRight().getY() - b.b.getLowerLeft().getY();
 
-// if (xd <= 2 * d) d = xd / 2 - 1;
-// if (yd <= 2 * d) d = yd / 2 - 1;
+  if (xd <= 2 * d) d = xd / 2 - 1;
+  if (yd <= 2 * d) d = yd / 2 - 1;
 
-// Box<T> r(Point<T>(b.b.getLowerLeft().getX() + d,
-// b.b.getLowerLeft().getY() + d),
-// Point<T>(b.b.getUpperRight().getX() - d,
-// b.b.getUpperRight().getY() - d));
+  Box<T> r(
+      Point<T>(b.b.getLowerLeft().getX() + d, b.b.getLowerLeft().getY() + d),
+      Point<T>(b.b.getUpperRight().getX() - d, b.b.getUpperRight().getY() - d));
 
-// return RotatedBox<T>(r, b.rotateDeg, b.center);
-// }
+  return RotatedBox<T>(r, b.getDegree(), b.getCenter());
+}
 
 // _____________________________________________________________________________
 inline bool doubleEq(double a, double b) { return fabs(a - b) < 0.000001; }
@@ -199,6 +257,12 @@ inline bool contains(const LineSegment<T>& l, const Box<T>& box) {
   return contains(l.first, box) && contains(l.second, box);
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline bool contains(const Box<T>& b, const Box<T>& box) {
+  return contains(b.getLowerLeft(), box) && contains(b.getUpperRight(), box);
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline bool contains(const Point<T>& p, const LineSegment<T>& ls) {
@@ -236,13 +300,17 @@ inline int8_t polyContCheck(const Point<T>& a, Point<T> b, Point<T> c) {
   if (a.getY() == b.getY() && a.getY() == c.getY())
     return (!((b.getX() <= a.getX() && a.getX() <= c.getX()) ||
               (c.getX() <= a.getX() && a.getX() <= b.getX())));
-  if (a.getY() == b.getY() && a.getX() == b.getX()) return 0;
+  if (fabs(a.getY() - b.getY()) < EPSILON &&
+      fabs(a.getX() - b.getX()) < EPSILON)
+    return 0;
   if (b.getY() > c.getY()) {
     Point<T> tmp = b;
     b = c;
     c = tmp;
   }
-  if (a.getY() <= b.getY() || a.getY() > c.getY()) return 1;
+  if (a.getY() <= b.getY() || a.getY() > c.getY()) {
+    return 1;
+  }
 
   double d = (b.getX() - a.getX()) * (c.getY() - a.getY()) -
              (b.getY() - a.getY()) * (c.getX() - a.getX());
@@ -255,7 +323,9 @@ inline int8_t polyContCheck(const Point<T>& a, Point<T> b, Point<T> c) {
 template <typename T>
 inline bool contains(const Polygon<T>& polyC, const Polygon<T>& poly) {
   for (const auto& p : polyC.getOuter()) {
-    if (!contains(p, poly)) return false;
+    if (!contains(p, poly)) {
+      return false;
+    }
   }
   return true;
 }
@@ -264,7 +334,20 @@ inline bool contains(const Polygon<T>& polyC, const Polygon<T>& poly) {
 template <typename T>
 inline bool contains(const Line<T>& l, const Polygon<T>& poly) {
   for (const auto& p : l) {
-    if (!contains(p, poly)) return false;
+    if (!contains(p, poly)) {
+      return false;
+    }
+  }
+  return true;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool contains(const Line<T>& l, const Line<T>& other) {
+  for (const auto& p : l) {
+    if (!contains(p, other)) {
+      return false;
+    }
   }
   return true;
 }
@@ -298,6 +381,16 @@ inline bool contains(const Polygon<T>& poly, const Line<T>& l) {
   return true;
 }
 
+// _____________________________________________________________________________
+template <typename T, template <typename> typename GeometryA,
+          template <typename> typename GeometryB>
+inline bool contains(const std::vector<GeometryA<T>>& multigeo,
+                     const GeometryB<T>& geo) {
+  for (const auto& g : multigeo)
+    if (!contains(g, geo)) return false;
+  return true;
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline bool intersects(const LineSegment<T>& ls1, const LineSegment<T>& ls2) {
@@ -463,6 +556,13 @@ inline Point<T> intersection(const Point<T>& p1, const Point<T>& q1,
                       p2.getY(), q2.getX(), q2.getY());
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline Point<T> intersection(const LineSegment<T>& s1,
+                             const LineSegment<T>& s2) {
+  return intersection(s1.first, s1.second, s2.first, s2.second);
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline bool lineIntersects(T p1x, T p1y, T q1x, T q1y, T p2x, T p2y, T q2x,
@@ -534,11 +634,10 @@ inline double crossProd(const Point<T>& a, const Point<T>& b) {
 // _____________________________________________________________________________
 template <typename T>
 inline double crossProd(const Point<T>& p, const LineSegment<T>& ls) {
-  LineSegment<T> lss(Point<T>(0, 0),
-                     Point<T>(ls.second.getX() - ls.first.getX(),
-                              ls.second.getY() - ls.first.getY()));
-  return crossProd(lss.second, Point<T>(p.getX() - ls.first.getX(),
-                                        p.getY() - ls.first.getY()));
+  return crossProd(
+      Point<T>(ls.second.getX() - ls.first.getX(),
+               ls.second.getY() - ls.first.getY()),
+      Point<T>(p.getX() - ls.first.getX(), p.getY() - ls.first.getY()));
 }
 
 // _____________________________________________________________________________
@@ -555,6 +654,19 @@ inline std::string getWKT(const Point<T>& p) {
   return ss.str();
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline std::string getWKT(const std::vector<Point<T>>& p) {
+  std::stringstream ss;
+  ss << "MULTIPOINT (";
+  for (size_t i = 0; i < p.size(); i++) {
+    if (i) ss << ", ";
+    ss << "(" << p[i].getX() << " " << p[i].getY() << ")";
+  }
+  ss << ")";
+  return ss.str();
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline std::string getWKT(const Line<T>& l) {
@@ -568,6 +680,26 @@ inline std::string getWKT(const Line<T>& l) {
   return ss.str();
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline std::string getWKT(const std::vector<Line<T>>& ls) {
+  std::stringstream ss;
+  ss << "MULTILINESTRING (";
+
+  for (size_t j = 0; j < ls.size(); j++) {
+    if (j) ss << ", ";
+    ss << "(";
+    for (size_t i = 0; i < ls[j].size(); i++) {
+      if (i) ss << ", ";
+      ss << ls[j][i].getX() << " " << ls[j][i].getY();
+    }
+    ss << ")";
+  }
+
+  ss << ")";
+  return ss.str();
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline std::string getWKT(const LineSegment<T>& l) {
@@ -594,13 +726,35 @@ inline std::string getWKT(const Polygon<T>& p) {
   std::stringstream ss;
   ss << "POLYGON ((";
   for (size_t i = 0; i < p.getOuter().size(); i++) {
-    if (i) ss << ", ";
-    ss << p.getOuter()[i].getX() << " " << p.getOuter()[i].getY();
+    ss << p.getOuter()[i].getX() << " " << p.getOuter()[i].getY() << ", ";
   }
+  ss << p.getOuter().front().getX() << " " << p.getOuter().front().getY();
   ss << "))";
   return ss.str();
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline std::string getWKT(const std::vector<Polygon<T>>& ls) {
+  std::stringstream ss;
+  ss << "MULTIPOLYGON (";
+
+  for (size_t j = 0; j < ls.size(); j++) {
+    if (j) ss << ", ";
+    ss << "((";
+    for (size_t i = 0; i < ls[j].getOuter().size(); i++) {
+      ss << ls[j].getOuter()[i].getX() << " " << ls[j].getOuter()[i].getY()
+         << ", ";
+    }
+    ss << ls[j].getOuter().front().getX() << " "
+       << ls[j].getOuter().front().getY();
+    ss << "))";
+  }
+
+  ss << ")";
+  return ss.str();
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline double len(const Point<T>& g) {
@@ -633,6 +787,12 @@ inline Box<T> simplify(const Box<T>& g, double d) {
   return g;
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline RotatedBox<T> simplify(const RotatedBox<T>& g, double d) {
+  return g;
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline Line<T> simplify(const Line<T>& g, double d) {
@@ -658,6 +818,15 @@ inline Line<T> simplify(const Line<T>& g, double d) {
   return Line<T>{g.front(), g.back()};
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> simplify(const Polygon<T>& g, double d) {
+  auto simple = simplify(g, d);
+  std::rotate(simple.begin(), simple.begin() + simple.size() / 2, simple.end());
+  simple = simplify(simple, d);
+  return Polygon<T>(simple);
+}
+
 // _____________________________________________________________________________
 inline double distToSegment(double lax, double lay, double lbx, double lby,
                             double px, double py) {
@@ -755,43 +924,40 @@ inline double parallelity(const Box<T>& box, const Line<T>& line) {
 }
 
 // _____________________________________________________________________________
-// template <typename T>
-// inline double parallelity(const Box<T>& box, const MultiLine<T>& multiline) {
-// double ret = 0;
-// for (const Line<T>& l : multiline) {
-// ret += parallelity(box, l);
-// }
+template <typename T>
+inline double parallelity(const Box<T>& box, const MultiLine<T>& multiline) {
+  double ret = 0;
+  for (const Line<T>& l : multiline) {
+    ret += parallelity(box, l);
+  }
 
-// return ret / multiline.size();
-// }
+  return ret / multiline.size();
+}
 
 // _____________________________________________________________________________
-// template <typename T, template <typename> typename Geometry>
-// inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
-// // TODO: implement this nicer, works for now, but inefficient
-// // see
-// // https://geidav.wordpress.com/tag/gift-wrapping/#fn-1057-FreemanShapira1975
-// // for a nicer algorithm
+template <typename T, template <typename> typename Geometry>
+inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
+  // TODO: implement this nicer, works for now, but inefficient
+  // see
+  // https://geidav.wordpress.com/tag/gift-wrapping/#fn-1057-FreemanShapira1975
+  // for a nicer algorithm
 
-// Point<T> center;
-// bgeo::centroid(pol, center);
+  Point<T> center = centroid(pol);
+  Box<T> tmpBox = getBoundingBox(pol);
+  double rotateDeg = 0;
 
-// Box<T> tmpBox = getBoundingBox(pol);
-// double rotateDeg = 0;
+  // rotate in 1 deg steps
+  for (int i = 1; i < 360; i += 1) {
+    pol = rotate(pol, 1, center);
+    Box<T> e = getBoundingBox(pol);
+    if (area(tmpBox) > area(e)) {
+      tmpBox = e;
+      rotateDeg = i;
+    }
+  }
 
-// // rotate in 5 deg steps
-// for (int i = 1; i < 360; i += 1) {
-// pol = rotate(pol, 1, center);
-// Box<T> e;
-// bgeo::envelope(pol, e);
-// if (bgeo::area(tmpBox) > bgeo::area(e)) {
-// tmpBox = e;
-// rotateDeg = i;
-// }
-// }
-
-// return RotatedBox<T>(tmpBox, -rotateDeg, center);
-// }
+  return RotatedBox<T>(tmpBox, -rotateDeg, center);
+}
 
 // _____________________________________________________________________________
 template <typename T>
@@ -849,6 +1015,90 @@ inline Box<T> getBoundingBox(const Box<T>& b) {
   return b;
 }
 
+// _____________________________________________________________________________
+template <typename T, template <typename> typename GeometryA>
+inline Box<T> getBoundingBox(const std::vector<GeometryA<T>>& multigeo) {
+  Box<T> b;
+  b = extendBox(multigeo, b);
+  return b;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> convexHull(const Point<T>& p) {
+  return Polygon<T>({p});
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> convexHull(const Box<T>& b) {
+  return Polygon<T>(b);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> convexHull(const LineSegment<T>& b) {
+  return Polygon<T>(Line<T>{b.first, b.second});
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> convexHull(const RotatedBox<T>& b) {
+  auto p = convexHull(b.getBox());
+  p = rotate(p, b.getDegree(), b.getCenter());
+  return p;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline size_t convexHullImpl(const Line<T>& a, size_t p1, size_t p2,
+                             Line<double>* h, uint8_t d) {
+  // quickhull by Barber, Dobkin & Huhdanpaa
+  Point<T> pa;
+  bool found = false;
+  double maxDist = 0;
+  for (const auto& p : a) {
+    double tmpDist = distToSegment((*h)[p1], (*h)[p2], p);
+    double cp = crossProd(p, LineSegment<T>((*h)[p1], (*h)[p2]));
+    if (((cp > 0 && !d) || (cp < 0 && d)) && tmpDist > maxDist) {
+      pa = p;
+      found = true;
+      maxDist = tmpDist;
+    }
+  }
+
+  if (!found) return 0;
+
+  h->insert(h->begin() + p2 + !d, pa);
+  size_t in = 1 + convexHullImpl(a, p1, p2 + !d, h, d);
+  return in + convexHullImpl(a, p2 + in * d + 1 - 2 * d, p2 + in * d, h, d);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> convexHull(const Line<T>& l) {
+  if (l.size() == 2) return convexHull(LineSegment<T>(l[0], l[1]));
+  if (l.size() == 1) return convexHull(l[0]);
+
+  Point<T> left(std::numeric_limits<T>::max(), 0);
+  Point<T> right(std::numeric_limits<T>::min(), 0);
+  for (const auto& p : l) {
+    if (p.getX() <= left.getX()) left = p;
+    if (p.getX() >= right.getX()) right = p;
+  }
+
+  Line<T> hull{left, right};
+  convexHullImpl(l, 0, 1, &hull, 1);
+  convexHullImpl(l, 0, hull.size() - 1, &hull, 0);
+  return Polygon<T>(hull);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> convexHull(const Polygon<T>& p) {
+  return convexHull(p.getOuter());
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline Box<T> extendBox(const Line<T>& l, Box<T> b) {
@@ -864,6 +1114,55 @@ inline Box<T> extendBox(const LineSegment<T>& ls, Box<T> b) {
   return b;
 }
 
+// _____________________________________________________________________________
+template <typename T, template <typename> typename Geometry>
+inline Box<T> extendBox(const std::vector<Geometry<T>>& multigeom, Box<T> b) {
+  for (const auto& g : multigeom) b = extendBox(g, b);
+  return b;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline double area(const Point<T>& b) {
+  UNUSED(b);
+  return 0;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline double area(const LineSegment<T>& b) {
+  UNUSED(b);
+  return 0;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline double area(const Line<T>& b) {
+  UNUSED(b);
+  return 0;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline double area(const Box<T>& b) {
+  return (b.getUpperRight().getX() - b.getLowerLeft().getX()) *
+         (b.getUpperRight().getY() - b.getLowerLeft().getY());
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline double area(const Polygon<T>& b) {
+  double ret = 0;
+  size_t j = b.getOuter().size() - 1;
+  for (size_t i = 0; i < b.getOuter().size(); i++) {
+    ret += (b.getOuter()[j].getX() + b.getOuter()[i].getX()) *
+           (b.getOuter()[j].getY() - b.getOuter()[i].getY());
+    j = i;
+  }
+
+  return fabs(ret / 2.0);
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline double commonArea(const Box<T>& ba, const Box<T>& bb) {
@@ -877,76 +1176,68 @@ inline double commonArea(const Box<T>& ba, const Box<T>& bb) {
 }
 
 // _____________________________________________________________________________
-// template <typename T, template <typename> typename Geometry>
-// inline RotatedBox<T> getFullEnvelope(Geometry<T> pol) {
-// Point<T> center;
-// bgeo::centroid(pol, center);
+template <typename T, template <typename> typename Geometry>
+inline RotatedBox<T> getFullEnvelope(Geometry<T> pol) {
+  Point<T> center = centroid(pol);
+  Box<T> tmpBox = getBoundingBox(pol);
+  double rotateDeg = 0;
 
-// Box<T> tmpBox;
-// bgeo::envelope(pol, tmpBox);
-// double rotateDeg = 0;
+  std::vector<Polygon<T>> ml;
 
-// MultiPolygon<T> ml;
+  // rotate in 5 deg steps
+  for (int i = 1; i < 360; i += 1) {
+    pol = rotate(pol, 1, center);
+    Polygon<T> hull = convexHull(pol);
+    ml.push_back(hull);
+    Box<T> e = getBoundingBox(pol);
+    if (area(tmpBox) > area(e)) {
+      tmpBox = e;
+      rotateDeg = i;
+    }
+  }
 
-// // rotate in 5 deg steps
-// for (int i = 1; i < 360; i += 1) {
-// pol = rotate(pol, 1, center);
-// Polygon<T> hull;
-// bgeo::convex_hull(pol, hull);
-// ml.push_back(hull);
-// Box<T> e;
-// bgeo::envelope(pol, e);
-// if (bgeo::area(tmpBox) > bgeo::area(e)) {
-// tmpBox = e;
-// rotateDeg = i;
-// }
-// }
+  tmpBox = getBoundingBox(ml);
 
-// bgeo::envelope(ml, tmpBox);
-
-// return RotatedBox<T>(tmpBox, rotateDeg, center);
-// }
+  return RotatedBox<T>(tmpBox, rotateDeg, center);
+}
 
 // _____________________________________________________________________________
-// template <typename T>
-// inline RotatedBox<T> getOrientedEnvelopeAvg(MultiLine<T> ml) {
-// MultiLine<T> orig = ml;
-// // get oriented envelope for hull
-// RotatedBox<T> rbox = getFullEnvelope<T>(ml);
-// Point<T> center;
-// bgeo::centroid(rbox.b, center);
+template <typename T>
+inline RotatedBox<T> getOrientedEnvelopeAvg(MultiLine<T> ml) {
+  MultiLine<T> orig = ml;
+  // get oriented envelope for hull
+  RotatedBox<T> rbox = getFullEnvelope(ml);
+  Point<T> center = centroid(rbox.getBox());
 
-// ml = rotate(ml, -rbox.rotateDeg - 45, center);
+  ml = rotate(ml, -rbox.getDegree() - 45, center);
 
-// double bestDeg = -45;
-// double score = parallelity(rbox.b, ml);
+  double bestDeg = -45;
+  double score = parallelity(rbox.b, ml);
 
-// for (double i = -45; i <= 45; i += .5) {
-// ml = rotate(ml, -.5, center);
-// double p = parallelity(rbox.b, ml);
-// if (parallelity(rbox.b, ml) > score) {
-// bestDeg = i;
-// score = p;
-// }
-// }
+  for (double i = -45; i <= 45; i += .5) {
+    ml = rotate(ml, -.5, center);
+    double p = parallelity(rbox.b, ml);
+    if (parallelity(rbox.b, ml) > score) {
+      bestDeg = i;
+      score = p;
+    }
+  }
 
-// rbox.rotateDeg += bestDeg;
+  rbox.rotateDeg += bestDeg;
 
-// // move the box along 45deg angles from its origin until it fits the ml
-// // = until the intersection of its hull and the box is largest
-// Polygon<T> p = rbox.getPolygon();
-// p = rotate(p, -rbox.rotateDeg, rbox.center);
+  // move the box along 45deg angles from its origin until it fits the ml
+  // = until the intersection of its hull and the box is largest
+  Polygon<T> p = convexHull(rbox);
+  p = rotate(p, -rbox.getDegree(), rbox.getCenter());
 
-// Polygon<T> hull;
-// bgeo::convex_hull(orig, hull);
-// hull = rotate(hull, -rbox.rotateDeg, rbox.center);
+  Polygon<T> hull = convexHull(orig);
+  hull = rotate(hull, -rbox.getDegree(), rbox.getCenter());
 
-// Box<T> box;
-// bgeo::envelope(hull, box);
-// rbox = RotatedBox<T>(box, rbox.rotateDeg, rbox.center);
+  Box<T> box = getBoundingBox(hull);
+  rbox = RotatedBox<T>(box, rbox.getDegree(), rbox.getCenter());
 
-// return rbox;
-// }
+  return rbox;
+}
 
 // _____________________________________________________________________________
 template <typename T>
diff --git a/src/util/geo/Line.h b/src/util/geo/Line.h
index d7153d7..dfa29e2 100644
--- a/src/util/geo/Line.h
+++ b/src/util/geo/Line.h
@@ -11,12 +11,16 @@
 namespace util {
 namespace geo {
 
-template<typename T>
-using Line = std::vector<Point<T>>;
+template <typename T>
+class Line : public std::vector<Point<T>> {
+  using std::vector<Point<T>>::vector;
+};
 
-template<typename T>
+template <typename T>
 using LineSegment = std::pair<Point<T>, Point<T>>;
 
+template <typename T>
+using MultiLine = std::vector<Line<T>>;
 
 }  // namespace geo
 }  // namespace util
diff --git a/src/util/geo/Point.h b/src/util/geo/Point.h
index df2b9aa..923e9d4 100644
--- a/src/util/geo/Point.h
+++ b/src/util/geo/Point.h
@@ -39,6 +39,9 @@ class Point {
   T _x, _y;
 };
 
+template <typename T>
+using MultiPoint = std::vector<Point<T>>;
+
 }  // namespace geo
 }  // namespace util
 
diff --git a/src/util/geo/PolyLine.h b/src/util/geo/PolyLine.h
index 89295cc..b95f038 100644
--- a/src/util/geo/PolyLine.h
+++ b/src/util/geo/PolyLine.h
@@ -17,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() {}
@@ -46,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:
diff --git a/src/util/geo/PolyLine.tpp b/src/util/geo/PolyLine.tpp
index ec15b43..8d34c14 100644
--- a/src/util/geo/PolyLine.tpp
+++ b/src/util/geo/PolyLine.tpp
@@ -93,11 +93,11 @@ void PolyLine<T>::offsetPerp(double units) {
     n1 = n1 / n;
     n2 = n2 / n;
 
-    lastP.template set<0>(lastP.getX() + (n1 * units));
-    lastP.template set<1>(lastP.getY() + (n2 * units));
+    lastP.setX(lastP.getX() + (n1 * units));
+    lastP.setY(lastP.getY() + (n2 * units));
 
-    curP.template set<0>(curP.getX() + (n1 * units));
-    curP.template set<1>(curP.getY() + (n2 * units));
+    curP.setX(curP.getX() + (n1 * units));
+    curP.setY(curP.getY() + (n2 * units));
 
     if (lastIns && befLastIns &&
         lineIntersects(*lastIns, *befLastIns, lastP, curP)) {
@@ -460,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].getX() + vx);
-    _line[i].set<1>(_line[i].getY() + vy);
+    _line[i].setX(_line[i].getX() + vx);
+    _line[i].setY(_line[i].getY() + vy);
   }
 }
 
diff --git a/src/util/geo/Polygon.h b/src/util/geo/Polygon.h
index 3b1adac..d39c9a5 100644
--- a/src/util/geo/Polygon.h
+++ b/src/util/geo/Polygon.h
@@ -25,13 +25,16 @@ class Polygon {
                 b.getUpperRight(),
                 Point<T>(b.getLowerLeft().getX(), b.getUpperRight().getY())}) {}
 
-  const std::vector<Point<T>>& getOuter() const { return _outer; }
-  std::vector<Point<T>>& getOuter() { return _outer; }
+  const Line<T>& getOuter() const { return _outer; }
+  Line<T>& getOuter() { return _outer; }
 
  private:
-  std::vector<Point<T>> _outer;
+  Line<T> _outer;
 };
 
+template <typename T>
+using MultiPolyon = std::vector<Polygon<T>>;
+
 }  // namespace geo
 }  // namespace util
 
diff --git a/src/util/tests/TestMain.cpp b/src/util/tests/TestMain.cpp
index cd43ed3..d5a9220 100644
--- a/src/util/tests/TestMain.cpp
+++ b/src/util/tests/TestMain.cpp
@@ -272,38 +272,38 @@ CASE("frechet distance") {
 
 // ___________________________________________________________________________
 CASE("geo box alignment") {
-  // Line<double> a;
-  // a.push_back(Point<double>(1, 1));
-  // a.push_back(Point<double>(1, 2));
+  Line<double> a;
+  a.push_back(Point<double>(1, 1));
+  a.push_back(Point<double>(1, 2));
 
-  // Line<double> b;
-  // b.push_back(Point<double>(1, 2));
-  // b.push_back(Point<double>(2, 2));
+  Line<double> b;
+  b.push_back(Point<double>(1, 2));
+  b.push_back(Point<double>(2, 2));
 
-  // Line<double> c;
-  // c.push_back(Point<double>(2, 2));
-  // c.push_back(Point<double>(2, 1));
+  Line<double> c;
+  c.push_back(Point<double>(2, 2));
+  c.push_back(Point<double>(2, 1));
 
-  // Line<double> d;
-  // d.push_back(Point<double>(2, 1));
-  // d.push_back(Point<double>(1, 1));
+  Line<double> d;
+  d.push_back(Point<double>(2, 1));
+  d.push_back(Point<double>(1, 1));
 
-  // Box<double> box(Point<double>(2, 3), Point<double>(5, 4));
-  // MultiLine<double> ml;
-  // ml.push_back(a);
-  // ml.push_back(b);
-  // ml.push_back(c);
-  // ml.push_back(d);
+  Box<double> box(Point<double>(2, 3), Point<double>(5, 4));
+  MultiLine<double> ml;
+  ml.push_back(a);
+  ml.push_back(b);
+  ml.push_back(c);
+  ml.push_back(d);
 
-  // EXPECT(parallelity(box, ml) == approx(1));
-  // ml = rotate(ml, 45);
-  // EXPECT(parallelity(box, ml) == approx(0));
-  // ml = rotate(ml, 45);
-  // EXPECT(parallelity(box, ml) == approx(1));
-  // ml = rotate(ml, 45);
-  // EXPECT(parallelity(box, ml) == approx(0));
-  // ml = rotate(ml, 45);
-  // EXPECT(parallelity(box, ml) == approx(1));
+  EXPECT(parallelity(box, ml) == approx(1));
+  ml = rotate(ml, 45);
+  EXPECT(parallelity(box, ml) == approx(0));
+  ml = rotate(ml, 45);
+  EXPECT(parallelity(box, ml) == approx(1));
+  ml = rotate(ml, 45);
+  EXPECT(parallelity(box, ml) == approx(0));
+  ml = rotate(ml, 45);
+  EXPECT(parallelity(box, ml) == approx(1));
 },
 
 // ___________________________________________________________________________
@@ -944,7 +944,7 @@ CASE("geometry") {
   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))");
+  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));
@@ -977,8 +977,215 @@ CASE("geometry") {
   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)));
 }
 
 };