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..b920872 100644
--- a/README.md
+++ b/README.md
@@ -19,7 +19,6 @@ Precise map-matching for public transit schedules (GTFS data).
 
  * `cmake`
  * `gcc` >= 4.8
- * `libboost-system` >= 1.56, `libboost-filesystem` >= 1.56, `libboost-geometry` >= 1.56
 
 ## Building and Installation
 
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/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.tpp b/src/util/geo/BezierCurve.tpp
index da27f83..fb7e6ca 100644
--- a/src/util/geo/BezierCurve.tpp
+++ b/src/util/geo/BezierCurve.tpp
@@ -15,15 +15,15 @@ BezierCurve<T>::BezierCurve(const Point<T>& a, const Point<T>& b,
 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;
+  _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/Geo.h b/src/util/geo/Geo.h
index 27a1391..5568992 100644
--- a/src/util/geo/Geo.h
+++ b/src/util/geo/Geo.h
@@ -6,11 +6,13 @@
 
 #define _USE_MATH_DEFINES
 
+#include <sstream>
 #include <math.h>
 #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
@@ -37,11 +39,11 @@ typedef Box<double> DBox;
 typedef Box<float> FBox;
 typedef Box<int> IBox;
 
-const static double EPSILON = 0.0000001;
+typedef Polygon<double> DPolygon;
+typedef Polygon<float> FPolygon;
+typedef Polygon<int> IPolygon;
 
-// typedef Polygon<double> DPolygon;
-// typedef Polygon<float> FPolygon;
-// typedef Polygon<int> IPolygon;
+const static double EPSILON = 0.0000001;
 
 // _____________________________________________________________________________
 // template <typename T>
@@ -224,6 +226,32 @@ inline bool intersects(const LineSegment<T>& ls1, const LineSegment<T>& ls2) {
           ((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 Box<T>& b1, const Box<T>& b2) {
@@ -256,15 +284,39 @@ inline bool intersects(const LineSegment<T>& ls, const Box<T>& b) {
   return contains(ls, b);
 }
 
+// _____________________________________________________________________________
+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;
+    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);
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline Point<T> intersection(T p1x, T p1y, T q1x, T q1y, T p2x, T p2y, T q2x,
@@ -378,8 +430,7 @@ inline double crossProd(const Point<T>& p, const LineSegment<T>& ls) {
 // _____________________________________________________________________________
 template <typename T>
 inline double dist(const Point<T>& p1, const Point<T>& p2) {
-  return sqrt((p1.getX() - p2.getX()) * (p1.getX() - p2.getX()) +
-              (p1.getY() - p2.getY()) * (p1.getY() - p2.getY()));
+  return dist(p1.getX(), p1.getY(), p2.getX(), p2.getY());
 }
 
 // _____________________________________________________________________________
@@ -399,6 +450,34 @@ inline std::string getWKT(const Line<T>& l) {
     if (i) ss << ", ";
     ss << l[i].getX() << " " << l[i].getY();
   }
+  ss << ")";
+  return ss.str();
+}
+
+// _____________________________________________________________________________
+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++) {
+    if (i) ss << ", ";
+    ss << p.getOuter()[i].getX() << " " << p.getOuter()[i].getY();
+  }
+  ss << "))";
   return ss.str();
 }
 
@@ -424,10 +503,41 @@ inline Point<T> simplify(const Point<T>& g, double d) {
 
 // _____________________________________________________________________________
 template <typename T>
-inline Line<T> simplify(const Line<T>& g, double d) {
+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 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()};
+}
+
 // _____________________________________________________________________________
 inline double distToSegment(double lax, double lay, double lbx, double lby,
                             double px, double py) {
@@ -453,6 +563,13 @@ inline double distToSegment(const Point<T>& la, const Point<T>& lb,
                        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,
@@ -591,7 +708,7 @@ inline Box<T> getBoundingBox(const Point<T>& p) {
 template <typename T>
 inline Box<T> getBoundingBox(const Line<T>& l) {
   Box<T> ret;
-  for (const auto& p : l) extendBox(p, ret);
+  for (const auto& p : l) ret = extendBox(p, ret);
   return ret;
 }
 
diff --git a/src/util/geo/Grid.tpp b/src/util/geo/Grid.tpp
index ece7481..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.getUpperRight().template getX() - bbox.getLowerLeft().template getX();
+      bbox.getUpperRight().getX() - bbox.getLowerLeft().getX();
   _height =
-      bbox.getUpperRight().template getY() - bbox.getLowerLeft().template getY();
+      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.getLowerLeft().template getX());
-  size_t swY = getCellYFromY(box.getLowerLeft().template getY());
+  size_t swX = getCellXFromX(box.getLowerLeft().getX());
+  size_t swY = getCellYFromY(box.getLowerLeft().getY());
 
-  size_t neX = getCellXFromX(box.getUpperRight().template getX());
-  size_t neY = getCellYFromY(box.getUpperRight().template getY());
+  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.getLowerLeft().template getX());
-  size_t swY = getCellYFromY(box.getLowerLeft().template getY());
+  size_t swX = getCellXFromX(box.getLowerLeft().getX());
+  size_t swY = getCellYFromY(box.getLowerLeft().getY());
 
-  size_t neX = getCellXFromX(box.getUpperRight().template getX());
-  size_t neY = getCellYFromY(box.getUpperRight().template getY());
+  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.getLowerLeft().template getX() - d,
-                    a.getLowerLeft().template getY() - d),
-           Point<T>(a.getUpperRight().template getX() + d,
-                    a.getUpperRight().template getY() + 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.getLowerLeft().template getX() + x * _cellWidth,
-              _bb.getLowerLeft().template getY() + y * _cellHeight);
-  Point<T> ne(_bb.getLowerLeft().template getX() + (x + 1) * _cellWidth,
-              _bb.getLowerLeft().template getY() + (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.getLowerLeft().template getX();
+  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.getLowerLeft().template getY();
+  float dist = y - _bb.getLowerLeft().getY();
   if (dist < 0) dist = 0;
   return floor(dist / _cellHeight);
 }
diff --git a/src/util/geo/Point.h b/src/util/geo/Point.h
index de6bfd6..942a361 100644
--- a/src/util/geo/Point.h
+++ b/src/util/geo/Point.h
@@ -11,6 +11,7 @@ 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; }
diff --git a/src/util/geo/PolyLine.h b/src/util/geo/PolyLine.h
index 5360351..89295cc 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>
diff --git a/src/util/geo/PolyLine.tpp b/src/util/geo/PolyLine.tpp
index 7ace68a..ec15b43 100644
--- a/src/util/geo/PolyLine.tpp
+++ b/src/util/geo/PolyLine.tpp
@@ -66,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;
@@ -87,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
@@ -98,11 +93,11 @@ void PolyLine<T>::offsetPerp(double units) {
     n1 = n1 / n;
     n2 = n2 / n;
 
-    lastP.template set<0>(lastP.template get<0>() + (n1 * units));
-    lastP.template set<1>(lastP.template get<1>() + (n2 * units));
+    lastP.template set<0>(lastP.getX() + (n1 * units));
+    lastP.template set<1>(lastP.getY() + (n2 * units));
 
-    curP.template set<0>(curP.template get<0>() + (n1 * units));
-    curP.template set<1>(curP.template get<1>() + (n2 * units));
+    curP.template set<0>(curP.getX() + (n1 * units));
+    curP.template set<1>(curP.getY() + (n2 * units));
 
     if (lastIns && befLastIns &&
         lineIntersects(*lastIns, *befLastIns, lastP, curP)) {
@@ -213,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) {
@@ -239,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));
 }
 
 // _____________________________________________________________________________
@@ -304,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);
@@ -465,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].set<0>(_line[i].getX() + vx);
+    _line[i].set<1>(_line[i].getY() + vy);
   }
 }
 
@@ -624,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));
 }
@@ -648,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);
 }
@@ -711,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..2841830
--- /dev/null
+++ b/src/util/geo/Polygon.h
@@ -0,0 +1,32 @@
+// 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 "./Point.h"
+#include "./Line.h"
+
+namespace util {
+namespace geo {
+
+template <typename T>
+class Polygon {
+ public:
+  // maximum inverse box as default value of box
+  Polygon(const Line<T>& l) : _outer(l) {}
+
+  const std::vector<Point<T>>& getOuter() const { return _outer; }
+  std::vector<Point<T>>& getOuter() {return _outer; }
+
+ private:
+  std::vector<Point<T>> _outer;
+};
+
+
+}  // 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/tests/TestMain.cpp b/src/util/tests/TestMain.cpp
index e439d8d..b15c2b5 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);
@@ -907,6 +906,42 @@ CASE("geometry") {
   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));
 }
 
 };