diff --git a/.travis.yml b/.travis.yml
new file mode 100644
index 0000000..6498482
--- /dev/null
+++ b/.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
diff --git a/CMakeLists.txt b/CMakeLists.txt
index abe2aa0..b17e6f7 100644
--- a/CMakeLists.txt
+++ b/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}")
diff --git a/README.md b/README.md
index 1241e6e..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,8 +7,7 @@ Precise map-matching for public transit schedules (GTFS data).
 ## Requirements
 
  * `cmake`
- * `gcc` >= 4.8
- * `libboost-system` >= 1.56, `libboost-filesystem` >= 1.56, `libboost-geometry` >= 1.56
+ * `gcc` >= 4.8 (may work on lower versions, untested)
 
 ## 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).
 
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/cppgtfs b/src/cppgtfs
index 7c39553..8d25597 160000
--- a/src/cppgtfs
+++ b/src/cppgtfs
@@ -1 +1 @@
-Subproject commit 7c395530d3c5b2d5c147fdabaaeeb0d1babceb04
+Subproject commit 8d255970a97c31891277f63ba8b39a3ea2b7e133
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/pfaedle/eval/Collector.cpp b/src/pfaedle/eval/Collector.cpp
index 7b162ee..334641f 100644
--- a/src/pfaedle/eval/Collector.cpp
+++ b/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>() +
-                                 oldSegs.back().back().get<1>()) /
+  double fac = cos(2 * atan(exp((oldSegs.front().front().getY() +
+                                 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);
 
diff --git a/src/pfaedle/osm/BBoxIdx.cpp b/src/pfaedle/osm/BBoxIdx.cpp
index 8344a55..e5d7c08 100644
--- a/src/pfaedle/osm/BBoxIdx.cpp
+++ b/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>(),
-                                        _root.box.min_corner().get<0>()),
-      util::geo::latLngToWebMerc<float>(_root.box.max_corner().get<1>(),
-                                        _root.box.max_corner().get<0>()));
+      util::geo::latLngToWebMerc<float>(_root.box.getLowerLeft().getY(),
+                                        _root.box.getLowerLeft().getX()),
+      util::geo::latLngToWebMerc<float>(_root.box.getUpperRight().getY(),
+                                        _root.box.getUpperRight().getX()));
 }
 
 // _____________________________________________________________________________
diff --git a/src/pfaedle/osm/OsmBuilder.cpp b/src/pfaedle/osm/OsmBuilder.cpp
index bb3ffd4..455bb1d 100644
--- a/src/pfaedle/osm/OsmBuilder.cpp
+++ b/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);
 }
 
diff --git a/src/pfaedle/router/Router.cpp b/src/pfaedle/router/Router.cpp
index b96aaf6..5cac1cb 100644
--- a/src/pfaedle/router/Router.cpp
+++ b/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>();
-    y += to->pl().getGeom()->get<1>();
+    x += to->pl().getGeom()->getX();
+    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>();
-    y += to->getFrom()->pl().getGeom()->get<1>();
+    x += to->getFrom()->pl().getGeom()->getX();
+    y += to->getFrom()->pl().getGeom()->getY();
     c++;
   }
 
diff --git a/src/pfaedle/router/ShapeBuilder.cpp b/src/pfaedle/router/ShapeBuilder.cpp
index a85286a..ede3c64 100644
--- a/src/pfaedle/router/ShapeBuilder.cpp
+++ b/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>(),
-                                         hop.start->pl().getGeom()->get<1>());
+      FPoint ll = webMercToLatLng<float>(hop.start->pl().getGeom()->getX(),
+                                         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>(),
-                                    hop.end->pl().getGeom()->get<1>());
-        ret->addPoint(ShapePoint(ll.get<1>(), ll.get<0>(), dist, seq));
+        ll = webMercToLatLng<float>(hop.end->pl().getGeom()->getX(),
+                                    hop.end->pl().getGeom()->getY());
+        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>());
-            ret->addPoint(ShapePoint(ll.get<1>(), ll.get<0>(), dist, seq));
+            FPoint ll = webMercToLatLng<float>(cur.getX(), cur.getY());
+            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>());
-            ret->addPoint(ShapePoint(ll.get<1>(), ll.get<0>(), dist, seq));
+            FPoint ll = webMercToLatLng<float>(cur.getX(), cur.getY());
+            ret->addPoint(ShapePoint(ll.getY(), ll.getX(), dist, seq));
             seq++;
             lastDist = dist;
           }
diff --git a/src/util/Misc.h b/src/util/Misc.h
index 1857db9..eecde29 100644
--- a/src/util/Misc.h
+++ b/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)
diff --git a/src/util/geo/BezierCurve.h b/src/util/geo/BezierCurve.h
index 29f669d..f3d25c5 100644
--- a/src/util/geo/BezierCurve.h
+++ b/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"
 }
 }
 
diff --git a/src/util/geo/BezierCurve.tpp b/src/util/geo/BezierCurve.tpp
index 6e7ed14..fb7e6ca 100644
--- a/src/util/geo/BezierCurve.tpp
+++ b/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,
-                         const Point<T>& d)
+BezierCurve<T>::BezierCurve(const Point<T>& a, const Point<T>& b,
+                            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,
-                                 const Point<T>& d) {
-  _xp.a = a.template get<0>();
-  _xp.b = 3.0 * (b.template get<0>() - a.template get<0>());
-  _xp.c = 3.0 * (c.template get<0>() - b.template get<0>()) - _xp.b;
-  _xp.d = d.template get<0>() - a.template get<0>() - _xp.c - _xp.b;
+void BezierCurve<T>::recalcPolynoms(const Point<T>& a, const Point<T>& b,
+                                    const Point<T>& c, const Point<T>& d) {
+  _xp.a = a.getX();
+  _xp.b = 3.0 * (b.getX() - a.getX());
+  _xp.c = 3.0 * (c.getX() - b.getX()) - _xp.b;
+  _xp.d = d.getX() - a.getX() - _xp.c - _xp.b;
 
-  _yp.a = a.template get<1>();
-  _yp.b = 3.0 * (b.template get<1>() - a.template get<1>());
-  _yp.c = 3.0 * (c.template get<1>() - b.template get<1>()) - _yp.b;
-  _yp.d = d.template get<1>() - a.template get<1>() - _yp.c - _yp.b;
+  _yp.a = a.getY();
+  _yp.b = 3.0 * (b.getY() - a.getY());
+  _yp.c = 3.0 * (c.getY() - b.getY()) - _yp.b;
+  _yp.d = d.getY() - a.getY() - _yp.c - _yp.b;
 
   _didRender = false;
 }
diff --git a/src/util/geo/Box.h b/src/util/geo/Box.h
new file mode 100644
index 0000000..f2e77de
--- /dev/null
+++ b/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_
diff --git a/src/util/geo/Geo.h b/src/util/geo/Geo.h
index db34ed4..a020c55 100644
--- a/src/util/geo/Geo.h
+++ b/src/util/geo/Geo.h
@@ -1,49 +1,40 @@
 // Copyright 2016, University of Freiburg,
 // Chair of Algorithms and Data Structures.
 // Authors: Patrick Brosi <brosi@informatik.uni-freiburg.de>
+
 #ifndef UTIL_GEO_GEO_H_
 #define UTIL_GEO_GEO_H_
 
 #define _USE_MATH_DEFINES
 
 #include <math.h>
-#include <boost/geometry.hpp>
-#include <boost/geometry/strategies/transform/matrix_transformers.hpp>
+#include <algorithm>
+#include <cassert>
+#include <iostream>
+#include <sstream>
 #include "util/Misc.h"
+#include "util/geo/Box.h"
+#include "util/geo/Line.h"
+#include "util/geo/Point.h"
+#include "util/geo/Polygon.h"
 
 // -------------------
 // Geometry stuff
 // ------------------
 
-namespace bgeo = boost::geometry;
-
 namespace util {
 namespace geo {
 
-template <typename T>
-using Point = bgeo::model::point<T, 2, bgeo::cs::cartesian>;
-
-template <typename T>
-using Line = bgeo::model::linestring<Point<T>>;
-
-template <typename T>
-using MultiLine = bgeo::model::multi_linestring<Line<T>>;
-
-template <typename T>
-using Polygon = bgeo::model::polygon<Point<T>>;
-
-template <typename T>
-using MultiPolygon = bgeo::model::multi_polygon<Polygon<T>>;
-
-template <typename T>
-using Box = bgeo::model::box<Point<T>>;
-
 // convenience aliases
 
 typedef Point<double> DPoint;
 typedef Point<float> FPoint;
 typedef Point<int> IPoint;
 
+typedef LineSegment<double> DLineSegment;
+typedef LineSegment<float> FLineSegment;
+typedef LineSegment<int> ILineSegment;
+
 typedef Line<double> DLine;
 typedef Line<float> FLine;
 typedef Line<int> ILine;
@@ -52,193 +43,486 @@ typedef Box<double> DBox;
 typedef Box<float> FBox;
 typedef Box<int> IBox;
 
-// _____________________________________________________________________________
-template <typename T>
-inline Line<T> rotate(const Line<T>& geo, double deg, const Point<T>& center) {
-  Line<T> ret;
+typedef Polygon<double> DPolygon;
+typedef Polygon<float> FPolygon;
+typedef Polygon<int> IPolygon;
 
-  bgeo::strategy::transform::translate_transformer<T, 2, 2> translate(
-      -center.template get<0>(), -center.template get<1>());
-  bgeo::strategy::transform::rotate_transformer<bgeo::degree, T, 2, 2> rotate(
-      deg);
-  bgeo::strategy::transform::translate_transformer<T, 2, 2> translateBack(
-      center.template get<0>(), center.template get<1>());
-
-  bgeo::strategy::transform::ublas_transformer<T, 2, 2> translateRotate(
-      prod(rotate.matrix(), translate.matrix()));
-  bgeo::strategy::transform::ublas_transformer<T, 2, 2> all(
-      prod(translateBack.matrix(), translateRotate.matrix()));
-
-  bgeo::transform(geo, ret, all);
-
-  return ret;
-}
-
-// _____________________________________________________________________________
-template <typename T>
-inline MultiLine<T> rotate(const MultiLine<T>& geo, double deg,
-                           const Point<T>& center) {
-  MultiLine<T> ret;
-
-  bgeo::strategy::transform::translate_transformer<T, 2, 2> translate(
-      -center.template get<0>(), -center.template get<1>());
-  bgeo::strategy::transform::rotate_transformer<bgeo::degree, T, 2, 2> rotate(
-      deg);
-  bgeo::strategy::transform::translate_transformer<T, 2, 2> translateBack(
-      center.template get<0>(), center.template get<1>());
-
-  bgeo::strategy::transform::ublas_transformer<T, 2, 2> translateRotate(
-      prod(rotate.matrix(), translate.matrix()));
-  bgeo::strategy::transform::ublas_transformer<T, 2, 2> all(
-      prod(translateBack.matrix(), translateRotate.matrix()));
-
-  bgeo::transform(geo, ret, all);
-
-  return ret;
-}
+const static double EPSILON = 0.00000000001;
+const static double RAD = 0.017453292519943295;  // PI/180
 
 // _____________________________________________________________________________
 template <typename T>
 inline Box<T> pad(const Box<T>& box, double padding) {
-  return Box<T>(Point<T>(box.min_corner().template get<0>() - padding,
-                         box.min_corner().template get<1>() - padding),
-                Point<T>(box.max_corner().template get<0>() + padding,
-                         box.max_corner().template get<1>() + padding));
+  return Box<T>(Point<T>(box.getLowerLeft().getX() - padding,
+                         box.getLowerLeft().getY() - padding),
+                Point<T>(box.getUpperRight().getX() + padding,
+                         box.getUpperRight().getY() + padding));
 }
 
 // _____________________________________________________________________________
 template <typename T>
-inline Line<T> rotate(const Line<T>& geo, double deg) {
-  Point<T> center;
-  bgeo::centroid(geo, center);
-  return rotate(geo, deg, center);
+inline Point<T> centroid(const Point<T> p) {
+  return p;
 }
 
 // _____________________________________________________________________________
 template <typename T>
-inline MultiLine<T> rotate(const MultiLine<T>& geo, double deg) {
-  Point<T> center;
-  bgeo::centroid(geo, center);
-  return rotate(geo, deg, center);
+inline Point<T> centroid(const LineSegment<T> ls) {
+  return Point<T>((ls.first.getX() + ls.second.getX()) / T(2),
+                  (ls.first.getY() + ls.second.getY()) / T(2));
 }
 
 // _____________________________________________________________________________
-template <template <typename> typename Geometry, typename T>
-inline Geometry<T> move(const Geometry<T>& geo, T x, T y) {
-  Geometry<T> ret;
-  bgeo::strategy::transform::translate_transformer<T, 2, 2> translate(x, y);
-  bgeo::transform(geo, ret, translate);
-  return ret;
-}
-
-// TODO: outfactor
-
 template <typename T>
-struct RotatedBox {
-  RotatedBox(const Box<T>& b, double rot, const Point<T>& center)
-      : b(b), rotateDeg(rot), center(center) {}
-  RotatedBox(const Box<T>& b, double rot) : b(b), rotateDeg(rot) {
-    bgeo::centroid(b, center);
+inline Point<T> centroid(const Line<T> ls) {
+  double x = 0, y = 0;
+  for (const auto& p : ls) {
+    x += p.getX();
+    y += p.getY();
   }
+  return Point<T>(x / T(ls.size()), y / T(ls.size()));
+}
 
-  Box<T> b;
-  double rotateDeg;
-  Point<T> center;
+// _____________________________________________________________________________
+template <typename T>
+inline Point<T> centroid(const Polygon<T> ls) {
+  return centroid(ls.getOuter());
+}
 
-  Polygon<T> getPolygon() {
-    Polygon<T> hull;
-    bgeo::convex_hull(b, hull);
-    return rotate(hull, rotateDeg, center);
-  }
-};
+// _____________________________________________________________________________
+template <typename T>
+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) {
+  return p;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Point<T> rotate(Point<T> p, double deg, const Point<T>& c) {
+  deg *= -RAD;
+  double si = sin(deg);
+  double co = cos(deg);
+  p = p - c;
+
+  return Point<T>(p.getX() * co - p.getY() * si,
+                  p.getX() * si + p.getY() * co) +
+         c;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline LineSegment<T> rotate(LineSegment<T> geo, double deg,
+                             const Point<T>& c) {
+  geo.first = rotate(geo.first, deg, c);
+  geo.second = rotate(geo.second, deg, c);
+  return geo;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline LineSegment<T> rotate(LineSegment<T> geo, double deg) {
+  return (geo, deg, centroid(geo));
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Line<T> rotate(Line<T> geo, double deg, const Point<T>& c) {
+  for (size_t i = 0; i < geo.size(); i++) geo[i] = rotate(geo[i], deg, c);
+  return geo;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> rotate(Polygon<T> geo, double deg, const Point<T>& c) {
+  for (size_t i = 0; i < geo.getOuter().size(); i++)
+    geo.getOuter()[i] = rotate(geo.getOuter()[i], deg, 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() {
-  return bgeo::make_inverse<Box<T>>();
+  return Box<T>();
 }
 
 // _____________________________________________________________________________
 template <typename T>
 inline RotatedBox<T> shrink(const RotatedBox<T>& b, double d) {
-  double xd =
-      b.b.max_corner().template get<0>() - b.b.min_corner().template get<0>();
-  double yd =
-      b.b.max_corner().template get<1>() - b.b.min_corner().template get<1>();
+  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;
 
-  Box<T> r(Point<T>(b.b.min_corner().template get<0>() + d,
-                    b.b.min_corner().template get<1>() + d),
-           Point<T>(b.b.max_corner().template get<0>() - d,
-                    b.b.max_corner().template get<1>() - 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; }
 
 // _____________________________________________________________________________
-template <typename Geometry, typename Box>
-inline bool contains(const Geometry& geom, const Box& box) {
-  return bgeo::within(geom, box);
+template <typename T>
+inline bool contains(const Point<T>& p, const Box<T>& box) {
+  return p.getX() >= box.getLowerLeft().getX() &&
+         p.getX() <= box.getUpperRight().getX() &&
+         p.getY() >= box.getLowerLeft().getY() &&
+         p.getY() <= box.getUpperRight().getY();
 }
 
 // _____________________________________________________________________________
 template <typename T>
-inline bool contains(const Point<T>& p1, const Point<T>& q1, const Point<T>& p2,
-                     const Point<T>& q2) {
-  Line<T> a, b;
-  a.push_back(p1);
-  a.push_back(q1);
-  b.push_back(p2);
-  b.push_back(q2);
-
-  return bgeo::covered_by(a, b);
+inline bool contains(const Line<T>& l, const Box<T>& box) {
+  for (const auto& p : l)
+    if (!contains(p, box)) return false;
+  return true;
 }
 
 // _____________________________________________________________________________
 template <typename T>
-inline bool contains(T p1x, T p1y, T q1x, T q1y, T p2x, T p2y, T q2x, T q2y) {
-  Point<T> p1(p1x, p1y);
-  Point<T> q1(q1x, q1y);
-  Point<T> p2(p2x, p2y);
-  Point<T> q2(q2x, q2y);
-
-  return contains(p1, q1, p2, q2);
+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 intersects(const Point<T>& p1, const Point<T>& q1,
-                       const Point<T>& p2, const Point<T>& q2) {
-  /*
-   * checks whether two line segments intersect
-   */
-  Line<T> a, b;
-  a.push_back(p1);
-  a.push_back(q1);
-  b.push_back(p2);
-  b.push_back(q2);
-
-  return !(contains(p1, q1, p2, q2) || contains(p2, q2, p1, q1)) &&
-         bgeo::intersects(a, b);
+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 intersects(T p1x, T p1y, T q1x, T q1y, T p2x, T p2y, T q2x, T q2y) {
-  /*
-   * checks whether two line segments intersect
-   */
-  Point<T> p1(p1x, p1y);
-  Point<T> q1(q1x, q1y);
-  Point<T> p2(p2x, p2y);
-  Point<T> q2(q2x, q2y);
+inline bool contains(const Point<T>& p, const LineSegment<T>& ls) {
+  return fabs(crossProd(p, ls)) < EPSILON && contains(p, getBoundingBox(ls));
+}
 
-  return intersects(p1, q1, p2, q2);
+// _____________________________________________________________________________
+template <typename T>
+inline bool contains(const Point<T>& p, const Line<T>& l) {
+  for (size_t i = 1; i < l.size(); i++) {
+    if (contains(p, LineSegment<T>(l[i - 1], l[i]))) return true;
+  }
+  return false;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool contains(const Point<T>& p, const Polygon<T>& poly) {
+  // see https://de.wikipedia.org/wiki/Punkt-in-Polygon-Test_nach_Jordan
+  int8_t c = -1;
+
+  for (size_t i = 1; i < poly.getOuter().size(); i++) {
+    c *= polyContCheck(p, poly.getOuter()[i - 1], poly.getOuter()[i]);
+    if (c == 0) return true;
+  }
+
+  c *= polyContCheck(p, poly.getOuter().back(), poly.getOuter()[0]);
+
+  return c >= 0;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+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 (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;
+  }
+
+  double d = (b.getX() - a.getX()) * (c.getY() - a.getY()) -
+             (b.getY() - a.getY()) * (c.getX() - a.getX());
+  if (d > 0) return -1;
+  if (d < 0) return 1;
+  return 0;
+}
+
+// _____________________________________________________________________________
+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;
+    }
+  }
+  return true;
+}
+
+// _____________________________________________________________________________
+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;
+    }
+  }
+  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;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool contains(const Box<T>& b, const Polygon<T>& poly) {
+  return contains(b.getLowerLeft(), poly) &&
+         contains(b.getUpperRight(), poly) &&
+         contains(Point<T>(b.getUpperRight().getX(), b.getLowerLeft().getY()),
+                  poly) &&
+         contains(Point<T>(b.getLowerLeft().getX(), b.getUpperRight().getY()),
+                  poly);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool contains(const Polygon<T>& poly, const Box<T>& b) {
+  for (const auto& p : poly.getOuter()) {
+    if (!contains(p, b)) return false;
+  }
+  return true;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool contains(const Polygon<T>& poly, const Line<T>& l) {
+  for (const auto& p : poly.getOuter()) {
+    if (!contains(p, l)) return false;
+  }
+  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) {
+  return intersects(getBoundingBox(ls1), getBoundingBox(ls2)) &&
+         (contains(ls1.first, ls2) || contains(ls1.second, ls2) ||
+          contains(ls2.first, ls1) || contains(ls2.second, ls1) ||
+          ((crossProd(ls1.first, ls2) < 0) ^
+           (crossProd(ls1.second, ls2) < 0)) ||
+          ((crossProd(ls2.first, ls1) < 0) ^ (crossProd(ls2.second, ls1) < 0)));
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Line<T>& ls1, const Line<T>& ls2) {
+  for (size_t i = 1; i < ls1.size(); i++) {
+    for (size_t j = 1; j < ls2.size(); j++) {
+      if (intersects(LineSegment<T>(ls1[i - 1], ls1[i]),
+                     LineSegment<T>(ls2[j - 1], ls2[j])))
+        return true;
+    }
+  }
+
+  return false;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Line<T>& l, const Point<T>& p) {
+  return contains(l, p);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Point<T>& p, const Line<T>& l) {
+  return intersects(l, p);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Polygon<T>& l, const Point<T>& p) {
+  return contains(l, p);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Point<T>& p, const Polygon<T>& l) {
+  return intersects(l, p);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Box<T>& b1, const Box<T>& b2) {
+  return b1.getLowerLeft().getX() <= b2.getUpperRight().getX() &&
+         b1.getUpperRight().getX() >= b2.getLowerLeft().getX() &&
+         b1.getLowerLeft().getY() <= b2.getUpperRight().getY() &&
+         b1.getUpperRight().getY() >= b2.getLowerLeft().getY();
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const LineSegment<T>& ls, const Box<T>& b) {
+  if (intersects(ls, LineSegment<T>(b.getLowerLeft(),
+                                    Point<T>(b.getUpperRight().getX(),
+                                             b.getLowerLeft().getY()))))
+    return true;
+  if (intersects(ls, LineSegment<T>(b.getLowerLeft(),
+                                    Point<T>(b.getLowerLeft().getX(),
+                                             b.getUpperRight().getY()))))
+    return true;
+  if (intersects(ls, LineSegment<T>(b.getUpperRight(),
+                                    Point<T>(b.getLowerLeft().getX(),
+                                             b.getUpperRight().getY()))))
+    return true;
+  if (intersects(ls, LineSegment<T>(b.getUpperRight(),
+                                    Point<T>(b.getUpperRight().getX(),
+                                             b.getLowerLeft().getY()))))
+    return true;
+
+  return contains(ls, b);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const LineSegment<T>& ls, const Polygon<T>& p) {
+  for (size_t i = 1; i < p.getOuter().size(); i++) {
+    if (intersects(LineSegment<T>(p.getOuter()[i - 1], p.getOuter()[i]), ls))
+      return true;
+  }
+
+  if (intersects(LineSegment<T>(p.getOuter().back(), p.getOuter().front()), ls))
+    return true;
+
+  return contains(ls, p);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Polygon<T>& p, const LineSegment<T>& ls) {
+  return intersects(ls, p);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Box<T>& b, const LineSegment<T>& ls) {
+  return intersects(ls, b);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Line<T>& l, const Box<T>& b) {
+  for (size_t i = 1; i < l.size(); i++) {
+    if (intersects(LineSegment<T>(l[i - 1], l[i]), b)) return true;
+  }
+  return false;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Box<T>& b, const Line<T>& l) {
+  return intersects(l, b);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Point<T>& p, const Box<T>& b) {
+  return contains(p, b);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline bool intersects(const Box<T>& b, const Point<T>& p) {
+  return intersects(p, b);
 }
 
 // _____________________________________________________________________________
@@ -268,10 +552,15 @@ inline Point<T> intersection(const Point<T>& p1, const Point<T>& q1,
   /*
    * calculates the intersection between two line segments
    */
-  return intersection(p1.template get<0>(), p1.template get<1>(),
-                      q1.template get<0>(), q1.template get<1>(),
-                      p2.template get<0>(), p2.template get<1>(),
-                      q2.template get<0>(), q2.template get<1>());
+  return intersection(p1.getX(), p1.getY(), q1.getX(), q1.getY(), p2.getX(),
+                      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);
 }
 
 // _____________________________________________________________________________
@@ -281,7 +570,6 @@ inline bool lineIntersects(T p1x, T p1y, T q1x, T q1y, T p2x, T p2y, T q2x,
   /*
    * checks whether two lines intersect
    */
-  double EPSILON = 0.0000001;
   double a = ((q2y - p2y) * (q1x - p1x)) - ((q2x - p2x) * (q1y - p1y));
 
   return a > EPSILON || a < -EPSILON;
@@ -294,10 +582,8 @@ inline bool lineIntersects(const Point<T>& p1, const Point<T>& q1,
   /*
    * checks whether two lines intersect
    */
-  return lineIntersects(p1.template get<0>(), p1.template get<1>(),
-                        q1.template get<0>(), q1.template get<1>(),
-                        p2.template get<0>(), p2.template get<1>(),
-                        q2.template get<0>(), q2.template get<1>());
+  return lineIntersects(p1.getX(), p1.getY(), q1.getX(), q1.getY(), p2.getX(),
+                        p2.getY(), q2.getX(), q2.getY());
 }
 
 // _____________________________________________________________________________
@@ -310,8 +596,7 @@ inline double angBetween(double p1x, double p1y, double q1x, double q1y) {
 // _____________________________________________________________________________
 template <typename T>
 inline double angBetween(const Point<T>& p1, const Point<T>& q1) {
-  return angBetween(p1.template get<0>(), p1.template get<1>(),
-                    q1.template get<0>(), q1.template get<1>());
+  return angBetween(p1.getX(), p1.getY(), q1.getX(), q1.getY());
 }
 
 // _____________________________________________________________________________
@@ -334,41 +619,214 @@ inline double innerProd(double x1, double y1, double x2, double y2, double x3,
 }
 
 // _____________________________________________________________________________
-template <typename G>
-inline double innerProd(const G& a, const G& b, const G& c) {
-  return innerProd(a.template get<0>(), a.template get<1>(),
-                   b.template get<0>(), b.template get<1>(),
-                   c.template get<0>(), c.template get<1>());
+template <typename T>
+inline double innerProd(const Point<T>& a, const Point<T>& b,
+                        const Point<T>& c) {
+  return innerProd(a.getX(), a.getY(), b.getX(), b.getY(), c.getX(), c.getY());
 }
 
 // _____________________________________________________________________________
-template <typename GeometryA, typename GeometryB>
-inline double dist(const GeometryA& p1, const GeometryB& p2) {
-  return bgeo::distance(p1, p2);
+template <typename T>
+inline double crossProd(const Point<T>& a, const Point<T>& b) {
+  return a.getX() * b.getY() - b.getX() * a.getY();
 }
 
 // _____________________________________________________________________________
-template <typename Geometry>
-inline std::string getWKT(Geometry g) {
+template <typename T>
+inline double crossProd(const Point<T>& p, const LineSegment<T>& ls) {
+  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()));
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline double dist(const Point<T>& p1, const Point<T>& p2) {
+  return dist(p1.getX(), p1.getY(), p2.getX(), p2.getY());
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline std::string getWKT(const Point<T>& p) {
   std::stringstream ss;
-  ss << bgeo::wkt(g);
+  ss << "POINT (" << p.getX() << " " << p.getY() << ")";
   return ss.str();
 }
 
 // _____________________________________________________________________________
-template <typename Geometry>
-inline double len(Geometry g) {
-  return bgeo::length(g);
+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 Geometry>
-inline Geometry simplify(Geometry g, double d) {
-  Geometry ret;
-  bgeo::simplify(g, ret, d);
+template <typename T>
+inline std::string getWKT(const Line<T>& l) {
+  std::stringstream ss;
+  ss << "LINESTRING (";
+  for (size_t i = 0; i < l.size(); i++) {
+    if (i) ss << ", ";
+    ss << l[i].getX() << " " << l[i].getY();
+  }
+  ss << ")";
+  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) {
+  return getWKT(Line<T>{l.first, l.second});
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline std::string getWKT(const Box<T>& l) {
+  std::stringstream ss;
+  ss << "POLYGON ((";
+  ss << l.getLowerLeft().getX() << " " << l.getLowerLeft().getY();
+  ss << ", " << l.getUpperRight().getX() << " " << l.getLowerLeft().getY();
+  ss << ", " << l.getUpperRight().getX() << " " << l.getUpperRight().getY();
+  ss << ", " << l.getLowerLeft().getX() << " " << l.getUpperRight().getY();
+  ss << ", " << l.getLowerLeft().getX() << " " << l.getLowerLeft().getY();
+  ss << "))";
+  return ss.str();
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline std::string getWKT(const Polygon<T>& p) {
+  std::stringstream ss;
+  ss << "POLYGON ((";
+  for (size_t i = 0; i < p.getOuter().size(); i++) {
+    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) {
+  return 0;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline double len(const Line<T>& g) {
+  double ret = 0;
+  for (size_t i = 1; i < g.size(); i++) ret += dist(g[i - 1], g[i]);
   return ret;
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline Point<T> simplify(const Point<T>& g, double d) {
+  return g;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline LineSegment<T> simplify(const LineSegment<T>& g, double d) {
+  return g;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+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) {
+  // douglas peucker
+  double maxd = 0;
+  size_t maxi = 0;
+  for (size_t i = 1; i < g.size() - 1; i++) {
+    double dt = distToSegment(g.front(), g.back(), g[i]);
+    if (dt > maxd) {
+      maxi = i;
+      maxd = dt;
+    }
+  }
+
+  if (maxd > d) {
+    auto a = simplify(Line<T>(g.begin(), g.begin() + maxi + 1), d);
+    const auto& b = simplify(Line<T>(g.begin() + maxi, g.end()), d);
+    a.insert(a.end(), b.begin() + 1, b.end());
+
+    return a;
+  }
+
+  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) {
@@ -390,39 +848,37 @@ inline double distToSegment(double lax, double lay, double lbx, double lby,
 template <typename T>
 inline double distToSegment(const Point<T>& la, const Point<T>& lb,
                             const Point<T>& p) {
-  return distToSegment(la.template get<0>(), la.template get<1>(),
-                       lb.template get<0>(), lb.template get<1>(),
-                       p.template get<0>(), p.template get<1>());
+  return distToSegment(la.getX(), la.getY(), lb.getX(), lb.getY(), p.getX(),
+                       p.getY());
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline double distToSegment(const LineSegment<T>& ls, const Point<T>& p) {
+  return distToSegment(ls.first.getX(), ls.first.getY(), ls.second.getX(),
+                       ls.second.getY(), p.getX(), p.getY());
 }
 
 // _____________________________________________________________________________
 template <typename T>
 inline Point<T> projectOn(const Point<T>& a, const Point<T>& b,
                           const Point<T>& c) {
-  if (doubleEq(a.template get<0>(), b.template get<0>()) &&
-      doubleEq(a.template get<1>(), b.template get<1>()))
-    return a;
-  if (doubleEq(a.template get<0>(), c.template get<0>()) &&
-      doubleEq(a.template get<1>(), c.template get<1>()))
-    return a;
-  if (doubleEq(b.template get<0>(), c.template get<0>()) &&
-      doubleEq(b.template get<1>(), c.template get<1>()))
-    return b;
+  if (doubleEq(a.getX(), b.getX()) && doubleEq(a.getY(), b.getY())) return a;
+  if (doubleEq(a.getX(), c.getX()) && doubleEq(a.getY(), c.getY())) return a;
+  if (doubleEq(b.getX(), c.getX()) && doubleEq(b.getY(), c.getY())) return b;
 
   double x, y;
 
-  if (c.template get<0>() == a.template get<0>()) {
+  if (c.getX() == a.getX()) {
     // infinite slope
-    x = a.template get<0>();
-    y = b.template get<1>();
+    x = a.getX();
+    y = b.getY();
   } else {
-    double m = (double)(c.template get<1>() - a.template get<1>()) /
-               (c.template get<0>() - a.template get<0>());
-    double bb = (double)a.template get<1>() - (m * a.template get<0>());
+    double m = (double)(c.getY() - a.getY()) / (c.getX() - a.getX());
+    double bb = (double)a.getY() - (m * a.getX());
 
-    x = (m * b.template get<1>() + b.template get<0>() - m * bb) / (m * m + 1);
-    y = (m * m * b.template get<1>() + m * b.template get<0>() + bb) /
-        (m * m + 1);
+    x = (m * b.getY() + b.getX() - m * bb) / (m * m + 1);
+    y = (m * m * b.getY() + m * b.getX() + bb) / (m * m + 1);
   }
 
   Point<T> ret = Point<T>(x, y);
@@ -440,18 +896,18 @@ template <typename T>
 inline double parallelity(const Box<T>& box, const Line<T>& line) {
   double ret = M_PI;
 
-  double a = angBetween(box.min_corner(),
-                        Point<T>(box.min_corner().template get<0>(),
-                                 box.max_corner().template get<1>()));
-  double b = angBetween(box.min_corner(),
-                        Point<T>(box.max_corner().template get<0>(),
-                                 box.min_corner().template get<1>()));
-  double c = angBetween(box.max_corner(),
-                        Point<T>(box.min_corner().template get<0>(),
-                                 box.max_corner().template get<1>()));
-  double d = angBetween(box.max_corner(),
-                        Point<T>(box.max_corner().template get<0>(),
-                                 box.min_corner().template get<1>()));
+  double a = angBetween(
+      box.getLowerLeft(),
+      Point<T>(box.getLowerLeft().getX(), box.getUpperRight().getY()));
+  double b = angBetween(
+      box.getLowerLeft(),
+      Point<T>(box.getUpperRight().getX(), box.getLowerLeft().getY()));
+  double c = angBetween(
+      box.getUpperRight(),
+      Point<T>(box.getLowerLeft().getX(), box.getUpperRight().getY()));
+  double d = angBetween(
+      box.getUpperRight(),
+      Point<T>(box.getUpperRight().getX(), box.getLowerLeft().getY()));
 
   double e = angBetween(line.front(), line.back());
 
@@ -478,12 +934,6 @@ inline double parallelity(const Box<T>& box, const MultiLine<T>& multiline) {
   return ret / multiline.size();
 }
 
-// _____________________________________________________________________________
-template <typename GeomA, typename GeomB>
-inline bool intersects(const GeomA& a, const GeomB& b) {
-  return bgeo::intersects(a, b);
-}
-
 // _____________________________________________________________________________
 template <typename T, template <typename> typename Geometry>
 inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
@@ -492,18 +942,15 @@ inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
   // 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;
 
-  // rotate in 5 deg steps
+  // rotate in 1 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)) {
+    Box<T> e = getBoundingBox(pol);
+    if (area(tmpBox) > area(e)) {
       tmpBox = e;
       rotateDeg = i;
     }
@@ -512,90 +959,244 @@ inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
   return RotatedBox<T>(tmpBox, -rotateDeg, center);
 }
 
-// _____________________________________________________________________________
-template <typename T>
-inline Box<T> getBoundingBox(Point<T> pol) {
-  Box<T> tmpBox;
-  bgeo::envelope(pol, tmpBox);
-  return tmpBox;
-}
-
-// _____________________________________________________________________________
-template <typename T>
-inline Box<T> getBoundingBox(Line<T> pol) {
-  Box<T> tmpBox;
-  bgeo::envelope(pol, tmpBox);
-  return tmpBox;
-}
-
-// _____________________________________________________________________________
-template <typename T>
-inline Box<T> getBoundingBox(Polygon<T> pol) {
-  Box<T> tmpBox;
-  bgeo::envelope(pol, tmpBox);
-  return tmpBox;
-}
-
 // _____________________________________________________________________________
 template <typename T>
 inline Box<T> extendBox(const Box<T>& a, Box<T> b) {
-  bgeo::expand(b, a);
+  b = extendBox(a.getLowerLeft(), b);
+  b = extendBox(a.getUpperRight(), b);
   return b;
 }
 
 // _____________________________________________________________________________
-template <typename G, typename T>
-inline Box<T> extendBox(G pol, Box<T> b) {
-  Box<T> tmp;
-  bgeo::envelope(pol, tmp);
-  bgeo::expand(b, tmp);
+template <typename T>
+inline Box<T> extendBox(const Point<T>& p, Box<T> b) {
+  if (p.getX() < b.getLowerLeft().getX()) b.getLowerLeft().setX(p.getX());
+  if (p.getY() < b.getLowerLeft().getY()) b.getLowerLeft().setY(p.getY());
+
+  if (p.getX() > b.getUpperRight().getX()) b.getUpperRight().setX(p.getX());
+  if (p.getY() > b.getUpperRight().getY()) b.getUpperRight().setY(p.getY());
   return b;
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline Box<T> getBoundingBox(const Point<T>& p) {
+  return Box<T>(p, p);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Box<T> getBoundingBox(const Line<T>& l) {
+  Box<T> ret;
+  for (const auto& p : l) ret = extendBox(p, ret);
+  return ret;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Box<T> getBoundingBox(const Polygon<T>& pol) {
+  Box<T> ret;
+  for (const auto& p : pol.getOuter()) ret = extendBox(p, ret);
+  return ret;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Box<T> getBoundingBox(const LineSegment<T>& ls) {
+  Box<T> b;
+  b = extendBox(ls.first, b);
+  b = extendBox(ls.second, b);
+  return b;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+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) {
+  for (const auto& p : l) b = extendBox(p, b);
+  return b;
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Box<T> extendBox(const LineSegment<T>& ls, Box<T> b) {
+  b = extendBox(ls.first, b);
+  b = extendBox(ls.second, 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) {
-  T l = std::max(ba.min_corner().template get<0>(),
-                 bb.min_corner().template get<0>());
-  T r = std::min(ba.max_corner().template get<0>(),
-                 bb.max_corner().template get<0>());
-  T b = std::max(ba.min_corner().template get<1>(),
-                 bb.min_corner().template get<1>());
-  T t = std::min(ba.max_corner().template get<1>(),
-                 bb.max_corner().template get<1>());
+  T l = std::max(ba.getLowerLeft().getX(), bb.getLowerLeft().getX());
+  T r = std::min(ba.getUpperRight().getX(), bb.getUpperRight().getX());
+  T b = std::max(ba.getLowerLeft().getY(), bb.getLowerLeft().getY());
+  T t = std::min(ba.getUpperRight().getY(), bb.getUpperRight().getY());
 
   if (l > r || b > t) return 0;
-
   return (r - l) * (t - b);
 }
 
 // _____________________________________________________________________________
 template <typename T, template <typename> typename Geometry>
 inline RotatedBox<T> getFullEnvelope(Geometry<T> pol) {
-  Point<T> center;
-  bgeo::centroid(pol, center);
-
-  Box<T> tmpBox;
-  bgeo::envelope(pol, tmpBox);
+  Point<T> center = centroid(pol);
+  Box<T> tmpBox = getBoundingBox(pol);
   double rotateDeg = 0;
 
-  MultiPolygon<T> ml;
+  std::vector<Polygon<T>> ml;
 
   // 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);
+    Polygon<T> hull = convexHull(pol);
     ml.push_back(hull);
-    Box<T> e;
-    bgeo::envelope(pol, e);
-    if (bgeo::area(tmpBox) > bgeo::area(e)) {
+    Box<T> e = getBoundingBox(pol);
+    if (area(tmpBox) > area(e)) {
       tmpBox = e;
       rotateDeg = i;
     }
   }
 
-  bgeo::envelope(ml, tmpBox);
+  tmpBox = getBoundingBox(ml);
 
   return RotatedBox<T>(tmpBox, rotateDeg, center);
 }
@@ -606,10 +1207,9 @@ inline RotatedBox<T> getOrientedEnvelopeAvg(MultiLine<T> ml) {
   MultiLine<T> orig = ml;
   // get oriented envelope for hull
   RotatedBox<T> rbox = getFullEnvelope(ml);
-  Point<T> center;
-  bgeo::centroid(rbox.b, center);
+  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);
@@ -627,16 +1227,14 @@ inline RotatedBox<T> getOrientedEnvelopeAvg(MultiLine<T> ml) {
 
   // 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);
+  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;
 }
@@ -651,15 +1249,13 @@ inline Line<T> densify(const Line<T>& l, double d) {
   ret.push_back(l.front());
 
   for (size_t i = 1; i < l.size(); i++) {
-    double segd = dist(l[i-1], l[i]);
-    double dx =
-        (l[i].template get<0>() - l[i-1].template get<0>()) / segd;
-    double dy =
-        (l[i].template get<1>() - l[i-1].template get<1>()) / segd;
+    double segd = dist(l[i - 1], l[i]);
+    double dx = (l[i].getX() - l[i - 1].getX()) / segd;
+    double dy = (l[i].getY() - l[i - 1].getY()) / segd;
     double curd = d;
     while (curd < segd) {
-      ret.push_back(Point<T>(l[i-1].template get<0>() + dx * curd,
-                             l[i-1].template get<1>() + dy * curd));
+      ret.push_back(
+          Point<T>(l[i - 1].getX() + dx * curd, l[i - 1].getY() + dy * curd));
       curd += d;
     }
 
@@ -715,29 +1311,29 @@ inline double frechetDist(const Line<T>& a, const Line<T>& b, double d) {
 // _____________________________________________________________________________
 template <typename T>
 inline double accFrechetDistC(const Line<T>& a, const Line<T>& b, double d) {
-
   auto p = densify(a, d);
   auto q = densify(b, d);
 
   std::vector<std::vector<double>> ca(p.size(),
-                                        std::vector<double>(q.size(), 0));
+                                      std::vector<double>(q.size(), 0));
 
-  for (size_t i = 0; i < p.size(); i++) ca[i][0] = std::numeric_limits<double>::infinity();
-  for (size_t j = 0; j < q.size(); j++) ca[0][j] = std::numeric_limits<double>::infinity();
+  for (size_t i = 0; i < p.size(); i++)
+    ca[i][0] = std::numeric_limits<double>::infinity();
+  for (size_t j = 0; j < q.size(); j++)
+    ca[0][j] = std::numeric_limits<double>::infinity();
   ca[0][0] = 0;
 
   for (size_t i = 1; i < p.size(); i++) {
     for (size_t j = 1; j < q.size(); j++) {
-      double d = util::geo::dist(p[i], q[j]) * util::geo::dist(p[i], p[i-1]);
-      ca[i][j] = d + std::min(ca[i-1][j], std::min(ca[i][j-1], ca[i-1][j-1]));
+      double d = util::geo::dist(p[i], q[j]) * util::geo::dist(p[i], p[i - 1]);
+      ca[i][j] =
+          d + std::min(ca[i - 1][j], std::min(ca[i][j - 1], ca[i - 1][j - 1]));
     }
   }
 
   return ca[p.size() - 1][q.size() - 1];
 }
 
-
-
 // _____________________________________________________________________________
 template <typename T>
 inline Point<T> latLngToWebMerc(double lat, double lng) {
@@ -762,8 +1358,8 @@ inline double webMercMeterDist(const G1& a, const G2& b) {
   // euclidean distance on web mercator is in meters on equator,
   // and proportional to cos(lat) in both y directions
 
-  double latA = 2 * atan(exp(a.template get<1>() / 6378137.0)) - 1.5707965;
-  double latB = 2 * atan(exp(b.template get<1>() / 6378137.0)) - 1.5707965;
+  double latA = 2 * atan(exp(a.getY() / 6378137.0)) - 1.5707965;
+  double latB = 2 * atan(exp(b.getY() / 6378137.0)) - 1.5707965;
 
   return util::geo::dist(a, b) * cos((latA + latB) / 2.0);
 }
diff --git a/src/util/geo/Grid.h b/src/util/geo/Grid.h
index adc72da..9156a11 100644
--- a/src/util/geo/Grid.h
+++ b/src/util/geo/Grid.h
@@ -7,6 +7,7 @@
 
 #include <set>
 #include <vector>
+#include <map>
 #include "util/geo/Geo.h"
 
 namespace util {
diff --git a/src/util/geo/Grid.tpp b/src/util/geo/Grid.tpp
index ee59632..3524c87 100644
--- a/src/util/geo/Grid.tpp
+++ b/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 swY = getCellYFromY(box.min_corner().template get<1>());
+  size_t swX = getCellXFromX(box.getLowerLeft().getX());
+  size_t swY = getCellYFromY(box.getLowerLeft().getY());
 
-  size_t neX = getCellXFromX(box.max_corner().template get<0>());
-  size_t neY = getCellYFromY(box.max_corner().template get<1>());
+  size_t neX = getCellXFromX(box.getUpperRight().getX());
+  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 swY = getCellYFromY(box.min_corner().template get<1>());
+  size_t swX = getCellXFromX(box.getLowerLeft().getX());
+  size_t swY = getCellYFromY(box.getLowerLeft().getY());
 
-  size_t neX = getCellXFromX(box.max_corner().template get<0>());
-  size_t neY = getCellYFromY(box.max_corner().template get<1>());
+  size_t neX = getCellXFromX(box.getUpperRight().getX());
+  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,
-                    a.min_corner().template get<1>() - d),
-           Point<T>(a.max_corner().template get<0>() + d,
-                    a.max_corner().template get<1>() + d));
+  Box<T> b(Point<T>(a.getLowerLeft().getX() - d,
+                    a.getLowerLeft().getY() - d),
+           Point<T>(a.getUpperRight().getX() + 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,
-              _bb.min_corner().template get<1>() + y * _cellHeight);
-  Point<T> ne(_bb.min_corner().template get<0>() + (x + 1) * _cellWidth,
-              _bb.min_corner().template get<1>() + (y + 1) * _cellHeight);
+  Point<T> sw(_bb.getLowerLeft().getX() + x * _cellWidth,
+              _bb.getLowerLeft().getY() + y * _cellHeight);
+  Point<T> ne(_bb.getLowerLeft().getX() + (x + 1) * _cellWidth,
+              _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);
 }
diff --git a/src/util/geo/Line.h b/src/util/geo/Line.h
new file mode 100644
index 0000000..dfa29e2
--- /dev/null
+++ b/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_
diff --git a/src/util/geo/Point.h b/src/util/geo/Point.h
new file mode 100644
index 0000000..923e9d4
--- /dev/null
+++ b/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_
diff --git a/src/util/geo/PolyLine.h b/src/util/geo/PolyLine.h
index 5360351..b95f038 100644
--- a/src/util/geo/PolyLine.h
+++ b/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:
diff --git a/src/util/geo/PolyLine.tpp b/src/util/geo/PolyLine.tpp
index 2da1525..8d34c14 100644
--- a/src/util/geo/PolyLine.tpp
+++ b/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 n2 = curP.template get<0>() - lastP.template get<0>();
+    double n1 = lastP.getY() - curP.getY();
+    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.set<1>(lastP.template get<1>() + (n2 * units));
+    lastP.setX(lastP.getX() + (n1 * units));
+    lastP.setY(lastP.getY() + (n2 * units));
 
-    curP.set<0>(curP.template get<0>() + (n1 * units));
-    curP.set<1>(curP.template get<1>() + (n2 * units));
+    curP.setX(curP.getX() + (n1 * 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 n2 = b.template get<1>() - a.template get<1>();
+  double n1 = b.getX() - a.getX();
+  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),
-                  a.template get<1>() + (n2 * p));
+  return Point<T>(a.getX() + (n1 * 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];
-        y += p.template get<1>() * weights[i];
+        x += p.getX() * weights[i];
+        y += p.getY() * weights[i];
       } else {
-        x += p.template get<0>();
-        y += p.template get<1>();
+        x += p.getX();
+        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].set<1>(_line[i].template get<1>() + vy);
+    _line[i].setX(_line[i].getX() + vx);
+    _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 angleY1 = avgP.template get<1>() + sin(angle + M_PI / 2) * length / 2;
+  double angleX1 = avgP.getX() + cos(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 angleY2 = avgP.template get<1>() + sin(angle + M_PI / 2) * -length / 2;
+  double angleX2 = avgP.getX() + cos(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 dy = (b.p.template get<1>() - a.p.template get<1>()) / d;
+  double dx = (b.p.getX() - a.p.getX()) / 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>(),
-                   0.75 * pA.template get<1>() + 0.25 * pB.template get<1>()));
+          Point<T>(0.75 * pA.getX() + 0.25 * pB.getX(),
+                   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>(),
-                   0.25 * pA.template get<1>() + 0.75 * pB.template get<1>()));
+          Point<T>(0.25 * pA.getX() + 0.75 * pB.getX(),
+                   0.25 * pA.getY() + 0.75 * pB.getY()));
     }
 
     smooth.push_back(_line.back());
diff --git a/src/util/geo/Polygon.h b/src/util/geo/Polygon.h
new file mode 100644
index 0000000..d39c9a5
--- /dev/null
+++ b/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_
diff --git a/src/util/geo/output/GeoJsonOutput.tpp b/src/util/geo/output/GeoJsonOutput.tpp
index 9ff6ad7..faec0c9 100644
--- a/src/util/geo/output/GeoJsonOutput.tpp
+++ b/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.template get<1>());
+  _wr.val(p.getX());
+  _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.template get<1>());
+    _wr.val(p.getX());
+    _wr.val(p.template getY());
     _wr.close();
   }
   _wr.close();
diff --git a/src/util/graph/Graph.h b/src/util/graph/Graph.h
index 40146c4..09d0465 100644
--- a/src/util/graph/Graph.h
+++ b/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);
diff --git a/src/util/tests/TestMain.cpp b/src/util/tests/TestMain.cpp
index 3e6fbc2..d5a9220 100644
--- a/src/util/tests/TestMain.cpp
+++ b/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));
-  l.push_back(Point<double>(2.5, 2));
+  l2.push_back(Point<double>(2.5, 1));
+  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() == 1);
+  EXPECT(ret.size() == (size_t)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)));
 }
 
 };