update util lib

CMakeLists.txt

@@ -18,7 +18,7 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/")
 set(EXECUTABLE_OUTPUT_PATH "${CMAKE_SOURCE_DIR}/build")
 
 # set compiler flags, see http://stackoverflow.com/questions/7724569/debug-vs-release-in-cmake
-set(CMAKE_CXX_FLAGS                "-Ofast -fno-signed-zeros -fno-trapping-math -Wall -Wno-format-extra-args -Wextra -Wformat-nonliteral -Wformat-security -Wformat=2 -Wextra -Wno-implicit-fallthrough -pedantic")
+set(CMAKE_CXX_FLAGS                "-Ofast -fno-signed-zeros -fno-trapping-math -Wall -Wno-format-extra-args -Wextra -Wformat-nonliteral -Wformat-security -Wformat=2 -Wextra -Wno-implicit-fallthrough -pedantic -Wno-keyword-macro")
 set(CMAKE_CXX_FLAGS_SANITIZE       "-Og -g -fsanitize=address -fsanitize=leak -fsanitize=undefined -DLOGLEVEL=3 -DPFAEDLE_DBG=1")
 set(CMAKE_CXX_FLAGS_PROFILE        "-g -pg -DLOGLEVEL=3 -DPFAEDLE_DBG=1")
 set(CMAKE_CXX_FLAGS_DEBUG          "-Og -g -DLOGLEVEL=3 -DPFAEDLE_DBG=1")

src/util/3rdparty/MurmurHash3.cpp

@@ -31,13 +31,11 @@
 
 #define FORCE_INLINE inline __attribute__((always_inline))
 
-inline uint32_t rotl32 ( uint32_t x, int8_t r )
+inline uint32_t rotl32(uint32_t x, int8_t r) {
-{
   return (x << r) | (x >> (32 - r));
 }
 
-inline uint64_t rotl64 ( uint64_t x, int8_t r )
+inline uint64_t rotl64(uint64_t x, int8_t r) {
-{
   return (x << r) | (x >> (64 - r));
 }
 
@@ -52,21 +50,14 @@ inline uint64_t rotl64 ( uint64_t x, int8_t r )
 // Block read - if your platform needs to do endian-swapping or can only
 // handle aligned reads, do the conversion here
 
-FORCE_INLINE uint32_t getblock32 ( const uint32_t * p, int i )
+FORCE_INLINE uint32_t getblock32(const uint32_t *p, int i) { return p[i]; }
-{
-  return p[i];
-}
 
-FORCE_INLINE uint64_t getblock64 ( const uint64_t * p, int i )
+FORCE_INLINE uint64_t getblock64(const uint64_t *p, int i) { return p[i]; }
-{
-  return p[i];
-}
 
 //-----------------------------------------------------------------------------
 // Finalization mix - force all bits of a hash block to avalanche
 
-FORCE_INLINE uint32_t fmix32 ( uint32_t h )
+FORCE_INLINE uint32_t fmix32(uint32_t h) {
-{
   h ^= h >> 16;
   h *= 0x85ebca6b;
   h ^= h >> 13;
@@ -78,8 +69,7 @@ FORCE_INLINE uint32_t fmix32 ( uint32_t h )
 
 //----------
 
-FORCE_INLINE uint64_t fmix64 ( uint64_t k )
+FORCE_INLINE uint64_t fmix64(uint64_t k) {
-{
   k ^= k >> 33;
   k *= BIG_CONSTANT(0xff51afd7ed558ccd);
   k ^= k >> 33;
@@ -91,9 +81,7 @@ FORCE_INLINE uint64_t fmix64 ( uint64_t k )
 
 //-----------------------------------------------------------------------------
 
-void MurmurHash3_x86_32 ( const void * key, int len,
+void MurmurHash3_x86_32(const void *key, int len, uint32_t seed, void *out) {
-                          uint32_t seed, void * out )
-{
   const uint8_t *data = (const uint8_t *)key;
   const int nblocks = len / 4;
 
@@ -107,8 +95,7 @@ void MurmurHash3_x86_32 ( const void * key, int len,
 
   const uint32_t *blocks = (const uint32_t *)(data + nblocks * 4);
 
-  for(int i = -nblocks; i; i++)
+  for (int i = -nblocks; i; i++) {
-  {
     uint32_t k1 = getblock32(blocks, i);
 
     k1 *= c1;
@@ -127,12 +114,19 @@ void MurmurHash3_x86_32 ( const void * key, int len,
 
   uint32_t k1 = 0;
 
-  switch(len & 3)
+  switch (len & 3) {
-  {
+    case 3:
-  case 3: k1 ^= tail[2] << 16;
+      k1 ^= tail[2] << 16;
-  case 2: k1 ^= tail[1] << 8;
+      // fall through
-  case 1: k1 ^= tail[0];
+    case 2:
-          k1 *= c1; k1 = ROTL32(k1,15); k1 *= c2; h1 ^= k1;
+      k1 ^= tail[1] << 8;
+      // fall through
+    case 1:
+      k1 ^= tail[0];
+      k1 *= c1;
+      k1 = ROTL32(k1, 15);
+      k1 *= c2;
+      h1 ^= k1;
   };
 
   //----------
@@ -147,9 +141,8 @@ void MurmurHash3_x86_32 ( const void * key, int len,
 
 //-----------------------------------------------------------------------------
 
-void MurmurHash3_x86_128 ( const void * key, const int len,
+void MurmurHash3_x86_128(const void *key, const int len, uint32_t seed,
-                           uint32_t seed, void * out )
+                         void *out) {
-{
   const uint8_t *data = (const uint8_t *)key;
   const int nblocks = len / 16;
 
@@ -168,28 +161,47 @@ void MurmurHash3_x86_128 ( const void * key, const int len,
 
   const uint32_t *blocks = (const uint32_t *)(data + nblocks * 16);
 
-  for(int i = -nblocks; i; i++)
+  for (int i = -nblocks; i; i++) {
-  {
     uint32_t k1 = getblock32(blocks, i * 4 + 0);
     uint32_t k2 = getblock32(blocks, i * 4 + 1);
     uint32_t k3 = getblock32(blocks, i * 4 + 2);
     uint32_t k4 = getblock32(blocks, i * 4 + 3);
 
-    k1 *= c1; k1  = ROTL32(k1,15); k1 *= c2; h1 ^= k1;
+    k1 *= c1;
+    k1 = ROTL32(k1, 15);
+    k1 *= c2;
+    h1 ^= k1;
 
-    h1 = ROTL32(h1,19); h1 += h2; h1 = h1*5+0x561ccd1b;
+    h1 = ROTL32(h1, 19);
+    h1 += h2;
+    h1 = h1 * 5 + 0x561ccd1b;
 
-    k2 *= c2; k2  = ROTL32(k2,16); k2 *= c3; h2 ^= k2;
+    k2 *= c2;
+    k2 = ROTL32(k2, 16);
+    k2 *= c3;
+    h2 ^= k2;
 
-    h2 = ROTL32(h2,17); h2 += h3; h2 = h2*5+0x0bcaa747;
+    h2 = ROTL32(h2, 17);
+    h2 += h3;
+    h2 = h2 * 5 + 0x0bcaa747;
 
-    k3 *= c3; k3  = ROTL32(k3,17); k3 *= c4; h3 ^= k3;
+    k3 *= c3;
+    k3 = ROTL32(k3, 17);
+    k3 *= c4;
+    h3 ^= k3;
 
-    h3 = ROTL32(h3,15); h3 += h4; h3 = h3*5+0x96cd1c35;
+    h3 = ROTL32(h3, 15);
+    h3 += h4;
+    h3 = h3 * 5 + 0x96cd1c35;
 
-    k4 *= c4; k4  = ROTL32(k4,18); k4 *= c1; h4 ^= k4;
+    k4 *= c4;
+    k4 = ROTL32(k4, 18);
+    k4 *= c1;
+    h4 ^= k4;
 
-    h4 = ROTL32(h4,13); h4 += h1; h4 = h4*5+0x32ac3b17;
+    h4 = ROTL32(h4, 13);
+    h4 += h1;
+    h4 = h4 * 5 + 0x32ac3b17;
   }
 
   //----------
@@ -202,47 +214,95 @@ void MurmurHash3_x86_128 ( const void * key, const int len,
   uint32_t k3 = 0;
   uint32_t k4 = 0;
 
-  switch(len & 15)
+  switch (len & 15) {
-  {
+    case 15:
-  case 15: k4 ^= tail[14] << 16;
+      k4 ^= tail[14] << 16;
-  case 14: k4 ^= tail[13] << 8;
+      // fall through
-  case 13: k4 ^= tail[12] << 0;
+    case 14:
-           k4 *= c4; k4  = ROTL32(k4,18); k4 *= c1; h4 ^= k4;
+      k4 ^= tail[13] << 8;
-
+      // fall through
-  case 12: k3 ^= tail[11] << 24;
+    case 13:
-  case 11: k3 ^= tail[10] << 16;
+      k4 ^= tail[12] << 0;
-  case 10: k3 ^= tail[ 9] << 8;
+      k4 *= c4;
-  case  9: k3 ^= tail[ 8] << 0;
+      k4 = ROTL32(k4, 18);
-           k3 *= c3; k3  = ROTL32(k3,17); k3 *= c4; h3 ^= k3;
+      k4 *= c1;
-
+      h4 ^= k4;
-  case  8: k2 ^= tail[ 7] << 24;
+      // fall through
-  case  7: k2 ^= tail[ 6] << 16;
+    case 12:
-  case  6: k2 ^= tail[ 5] << 8;
+      k3 ^= tail[11] << 24;
-  case  5: k2 ^= tail[ 4] << 0;
+      // fall through
-           k2 *= c2; k2  = ROTL32(k2,16); k2 *= c3; h2 ^= k2;
+    case 11:
-
+      k3 ^= tail[10] << 16;
-  case  4: k1 ^= tail[ 3] << 24;
+      // fall through
-  case  3: k1 ^= tail[ 2] << 16;
+    case 10:
-  case  2: k1 ^= tail[ 1] << 8;
+      k3 ^= tail[9] << 8;
-  case  1: k1 ^= tail[ 0] << 0;
+      // fall through
-           k1 *= c1; k1  = ROTL32(k1,15); k1 *= c2; h1 ^= k1;
+    case 9:
+      k3 ^= tail[8] << 0;
+      k3 *= c3;
+      k3 = ROTL32(k3, 17);
+      k3 *= c4;
+      h3 ^= k3;
+      // fall through
+    case 8:
+      k2 ^= tail[7] << 24;
+      // fall through
+    case 7:
+      k2 ^= tail[6] << 16;
+      // fall through
+    case 6:
+      k2 ^= tail[5] << 8;
+      // fall through
+    case 5:
+      k2 ^= tail[4] << 0;
+      k2 *= c2;
+      k2 = ROTL32(k2, 16);
+      k2 *= c3;
+      h2 ^= k2;
+      // fall through
+    case 4:
+      k1 ^= tail[3] << 24;
+      // fall through
+    case 3:
+      k1 ^= tail[2] << 16;
+      // fall through
+    case 2:
+      k1 ^= tail[1] << 8;
+      // fall through
+    case 1:
+      k1 ^= tail[0] << 0;
+      k1 *= c1;
+      k1 = ROTL32(k1, 15);
+      k1 *= c2;
+      h1 ^= k1;
   };
 
   //----------
   // finalization
 
-  h1 ^= len; h2 ^= len; h3 ^= len; h4 ^= len;
+  h1 ^= len;
+  h2 ^= len;
+  h3 ^= len;
+  h4 ^= len;
 
-  h1 += h2; h1 += h3; h1 += h4;
+  h1 += h2;
-  h2 += h1; h3 += h1; h4 += h1;
+  h1 += h3;
+  h1 += h4;
+  h2 += h1;
+  h3 += h1;
+  h4 += h1;
 
   h1 = fmix32(h1);
   h2 = fmix32(h2);
   h3 = fmix32(h3);
   h4 = fmix32(h4);
 
-  h1 += h2; h1 += h3; h1 += h4;
+  h1 += h2;
-  h2 += h1; h3 += h1; h4 += h1;
+  h1 += h3;
+  h1 += h4;
+  h2 += h1;
+  h3 += h1;
+  h4 += h1;
 
   ((uint32_t *)out)[0] = h1;
   ((uint32_t *)out)[1] = h2;
@@ -252,9 +312,8 @@ void MurmurHash3_x86_128 ( const void * key, const int len,
 
 //-----------------------------------------------------------------------------
 
-void MurmurHash3_x64_128 ( const void * key, const int len,
+void MurmurHash3_x64_128(const void *key, const int len, const uint32_t seed,
-                           const uint32_t seed, void * out )
+                         void *out) {
-{
   const uint8_t *data = (const uint8_t *)key;
   const int nblocks = len / 16;
 
@@ -269,18 +328,27 @@ void MurmurHash3_x64_128 ( const void * key, const int len,
 
   const uint64_t *blocks = (const uint64_t *)(data);
 
-  for(int i = 0; i < nblocks; i++)
+  for (int i = 0; i < nblocks; i++) {
-  {
     uint64_t k1 = getblock64(blocks, i * 2 + 0);
     uint64_t k2 = getblock64(blocks, i * 2 + 1);
 
-    k1 *= c1; k1  = ROTL64(k1,31); k1 *= c2; h1 ^= k1;
+    k1 *= c1;
+    k1 = ROTL64(k1, 31);
+    k1 *= c2;
+    h1 ^= k1;
 
-    h1 = ROTL64(h1,27); h1 += h2; h1 = h1*5+0x52dce729;
+    h1 = ROTL64(h1, 27);
+    h1 += h2;
+    h1 = h1 * 5 + 0x52dce729;
 
-    k2 *= c2; k2  = ROTL64(k2,33); k2 *= c1; h2 ^= k2;
+    k2 *= c2;
+    k2 = ROTL64(k2, 33);
+    k2 *= c1;
+    h2 ^= k2;
 
-    h2 = ROTL64(h2,31); h2 += h1; h2 = h2*5+0x38495ab5;
+    h2 = ROTL64(h2, 31);
+    h2 += h1;
+    h2 = h2 * 5 + 0x38495ab5;
   }
 
   //----------
@@ -291,32 +359,66 @@ void MurmurHash3_x64_128 ( const void * key, const int len,
   uint64_t k1 = 0;
   uint64_t k2 = 0;
 
-  switch(len & 15)
+  switch (len & 15) {
-  {
+    case 15:
-  case 15: k2 ^= ((uint64_t)tail[14]) << 48;
+      k2 ^= ((uint64_t)tail[14]) << 48;
-  case 14: k2 ^= ((uint64_t)tail[13]) << 40;
+      // fall through
-  case 13: k2 ^= ((uint64_t)tail[12]) << 32;
+    case 14:
-  case 12: k2 ^= ((uint64_t)tail[11]) << 24;
+      k2 ^= ((uint64_t)tail[13]) << 40;
-  case 11: k2 ^= ((uint64_t)tail[10]) << 16;
+      // fall through
-  case 10: k2 ^= ((uint64_t)tail[ 9]) << 8;
+    case 13:
-  case  9: k2 ^= ((uint64_t)tail[ 8]) << 0;
+      k2 ^= ((uint64_t)tail[12]) << 32;
-           k2 *= c2; k2  = ROTL64(k2,33); k2 *= c1; h2 ^= k2;
+      // fall through
-
+    case 12:
-  case  8: k1 ^= ((uint64_t)tail[ 7]) << 56;
+      k2 ^= ((uint64_t)tail[11]) << 24;
-  case  7: k1 ^= ((uint64_t)tail[ 6]) << 48;
+      // fall through
-  case  6: k1 ^= ((uint64_t)tail[ 5]) << 40;
+    case 11:
-  case  5: k1 ^= ((uint64_t)tail[ 4]) << 32;
+      k2 ^= ((uint64_t)tail[10]) << 16;
-  case  4: k1 ^= ((uint64_t)tail[ 3]) << 24;
+      // fall through
-  case  3: k1 ^= ((uint64_t)tail[ 2]) << 16;
+    case 10:
-  case  2: k1 ^= ((uint64_t)tail[ 1]) << 8;
+      k2 ^= ((uint64_t)tail[9]) << 8;
-  case  1: k1 ^= ((uint64_t)tail[ 0]) << 0;
+      // fall through
-           k1 *= c1; k1  = ROTL64(k1,31); k1 *= c2; h1 ^= k1;
+    case 9:
+      k2 ^= ((uint64_t)tail[8]) << 0;
+      k2 *= c2;
+      k2 = ROTL64(k2, 33);
+      k2 *= c1;
+      h2 ^= k2;
+      // fall through
+    case 8:
+      k1 ^= ((uint64_t)tail[7]) << 56;
+      // fall through
+    case 7:
+      k1 ^= ((uint64_t)tail[6]) << 48;
+      // fall through
+    case 6:
+      k1 ^= ((uint64_t)tail[5]) << 40;
+      // fall through
+    case 5:
+      k1 ^= ((uint64_t)tail[4]) << 32;
+      // fall through
+    case 4:
+      k1 ^= ((uint64_t)tail[3]) << 24;
+      // fall through
+    case 3:
+      k1 ^= ((uint64_t)tail[2]) << 16;
+      // fall through
+    case 2:
+      k1 ^= ((uint64_t)tail[1]) << 8;
+      // fall through
+    case 1:
+      k1 ^= ((uint64_t)tail[0]) << 0;
+      k1 *= c1;
+      k1 = ROTL64(k1, 31);
+      k1 *= c2;
+      h1 ^= k1;
   };
 
   //----------
   // finalization
 
-  h1 ^= len; h2 ^= len;
+  h1 ^= len;
+  h2 ^= len;
 
   h1 += h2;
   h2 += h1;
@@ -332,4 +434,3 @@ void MurmurHash3_x64_128 ( const void * key, const int len,
 }
 
 //-----------------------------------------------------------------------------
-

src/util/Misc.h

@@ -9,8 +9,10 @@
 #include <cstring>
 #include <chrono>
 #include <sstream>
+#include <iomanip>
 #include <immintrin.h>
 #include <iostream>
+#include <sstream>
 #include <vector>
 #include <unistd.h>
 #include <sys/types.h>
@@ -114,6 +116,21 @@ inline bool isFloatingPoint(const std::string& str) {
   return ss.eof() && ! ss.fail();
 }
 
+// _____________________________________________________________________________
+inline std::string formatFloat(double f, size_t digits) {
+  std::stringstream ss;
+  ss << std::fixed << std::setprecision(digits) << f;
+	std::string ret = ss.str();
+
+  if (ret.find('.') != std::string::npos) {
+		auto p = ret.find_last_not_of('0');
+    if (ret[p] == '.') return ret.substr(0, p);
+    return ret.substr(0, p + 1);
+  }
+
+  return ret;
+}
+
 // _____________________________________________________________________________
 inline double atof(const char* p, uint8_t mn) {
   // this atof implementation works only on "normal" float strings like

src/util/geo/Box.h

@@ -24,6 +24,9 @@ class Box {
   Point<T>& getLowerLeft() { return _ll; }
   Point<T>& getUpperRight() { return _ur; }
 
+  Point<T> getUpperLeft() { return {_ll.getX(), _ur.getY()}; }
+  Point<T> getLowerRight() { return {_ur.getX(), _ll.getY()}; }
+
   void setLowerLeft(const Point<T>& ll) { _ll = ll; }
   void setUpperRight(const Point<T>& ur) { _ur = ur; }
 

src/util/geo/Geo.h

@@ -10,6 +10,7 @@
 #include <math.h>
 #include <algorithm>
 #include <cassert>
+#include <iomanip>
 #include <iostream>
 #include <sstream>
 #include <stdexcept>
@@ -742,6 +743,19 @@ inline double angBetween(const Point<T>& p1, const Point<T>& q1) {
   return angBetween(p1.getX(), p1.getY(), q1.getX(), q1.getY());
 }
 
+// _____________________________________________________________________________
+template <typename T>
+inline double angBetween(const Point<T>& p1, const MultiPoint<T>& points) {
+  double sinSum = 0;
+  double cosSum = 0;
+  for (auto q1 : points) {
+    double a = angBetween(p1.getX(), p1.getY(), q1.getX(), q1.getY());
+    sinSum += sin(a);
+    cosSum += cos(a);
+  }
+  return atan2(sinSum / points.size(), cosSum / points.size());
+}
+
 // _____________________________________________________________________________
 inline double dist(double x1, double y1, double x2, double y2) {
   return sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
@@ -932,19 +946,19 @@ inline Line<T> lineFromWKT(std::string wkt) {
 // _____________________________________________________________________________
 template <typename T>
 inline std::string getWKT(const Point<T>& p) {
-  std::stringstream ss;
+  return std::string("POINT (") + formatFloat(p.getX(), 6) + " " +
-  ss << std::fixed << std::setprecision(6) << "POINT (" << p.getX() << " " << p.getY() << ")";
+         formatFloat(p.getY(), 6) + ")";
-  return ss.str();
 }
 
 // _____________________________________________________________________________
 template <typename T>
 inline std::string getWKT(const std::vector<Point<T>>& p) {
   std::stringstream ss;
-  ss << std::fixed << std::setprecision(6) << "MULTIPOINT (";
+  ss << "MULTIPOINT (";
   for (size_t i = 0; i < p.size(); i++) {
     if (i) ss << ", ";
-    ss << "(" << p[i].getX() << " " << p[i].getY() << ")";
+    ss << "(" << formatFloat(p.getX(), 6) << " " << formatFloat(p.getY(), 6)
+       << ")";
   }
   ss << ")";
   return ss.str();
@@ -954,10 +968,10 @@ inline std::string getWKT(const std::vector<Point<T>>& p) {
 template <typename T>
 inline std::string getWKT(const Line<T>& l) {
   std::stringstream ss;
-  ss << std::fixed << std::setprecision(6) << "LINESTRING (";
+  ss << "LINESTRING (";
   for (size_t i = 0; i < l.size(); i++) {
     if (i) ss << ", ";
-    ss << l[i].getX() << " " << l[i].getY();
+    ss << formatFloat(l[i].getX(), 6) << " " << formatFloat(l[i].getY(), 6);
   }
   ss << ")";
   return ss.str();
@@ -967,14 +981,15 @@ inline std::string getWKT(const Line<T>& l) {
 template <typename T>
 inline std::string getWKT(const std::vector<Line<T>>& ls) {
   std::stringstream ss;
-  ss << std::fixed << std::setprecision(6) << "MULTILINESTRING (";
+  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 << formatFloat(ls[j][i].getX(), 6) << " "
+         << formatFloat(ls[j][i].getY(), 6);
     }
     ss << ")";
   }
@@ -994,11 +1009,16 @@ template <typename T>
 inline std::string getWKT(const Box<T>& l) {
   std::stringstream ss;
   ss << "POLYGON ((";
-  ss << l.getLowerLeft().getX() << " " << l.getLowerLeft().getY();
+  ss << formatFloat(l.getLowerLeft().getX(), 6) << " "
-  ss << ", " << l.getUpperRight().getX() << " " << l.getLowerLeft().getY();
+     << formatFloat(l.getLowerLeft().getY(), 6);
-  ss << ", " << l.getUpperRight().getX() << " " << l.getUpperRight().getY();
+  ss << ", " << formatFloat(l.getUpperRight().getX(), 6) << " "
-  ss << ", " << l.getLowerLeft().getX() << " " << l.getUpperRight().getY();
+     << formatFloat(l.getLowerLeft().getY(), 6);
-  ss << ", " << l.getLowerLeft().getX() << " " << l.getLowerLeft().getY();
+  ss << ", " << formatFloat(l.getUpperRight().getX(), 6) << " "
+     << formatFloat(l.getUpperRight().getY(), 6);
+  ss << ", " << formatFloat(l.getLowerLeft().getX(), 6) << " "
+     << formatFloat(l.getUpperRight().getY(), 6);
+  ss << ", " << formatFloat(l.getLowerLeft().getX(), 6) << " "
+     << formatFloat(l.getLowerLeft().getY(), 6);
   ss << "))";
   return ss.str();
 }
@@ -1009,9 +1029,11 @@ 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 << formatFloat(p.getOuter()[i].getX(), 6) << " "
+       << formatFloat(p.getOuter()[i].getY(), 6) << ", ";
   }
-  ss << p.getOuter().front().getX() << " " << p.getOuter().front().getY();
+  ss << formatFloat(p.getOuter().front().getX(), 6) << " "
+     << formatFloat(p.getOuter().front().getY(), 6);
   ss << "))";
   return ss.str();
 }
@@ -1026,11 +1048,11 @@ inline std::string getWKT(const std::vector<Polygon<T>>& ls) {
     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 << formatFloat(ls[j].getOuter()[i].getX(), 6) << " "
-         << ", ";
+         << formatFloat(ls[j].getOuter()[i].getY(), 6) << ", ";
     }
-    ss << ls[j].getOuter().front().getX() << " "
+    ss << formatFloat(ls[j].getOuter().front().getX(), 6) << " "
-       << ls[j].getOuter().front().getY();
+       << formatFloat(ls[j].getOuter().front().getY(), 6);
     ss << "))";
   }
 
@@ -1224,7 +1246,8 @@ inline double parallelity(const Box<T>& box, const MultiLine<T>& multiline) {
 
 // _____________________________________________________________________________
 template <template <typename> class Geometry, typename T>
-inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
+inline RotatedBox<T> getOrientedEnvelope(std::vector<Geometry<T>> pol,
+                                         double step) {
   // TODO: implement this nicer, works for now, but inefficient
   // see
   // https://geidav.wordpress.com/tag/gift-wrapping/#fn-1057-FreemanShapira1975
@@ -1235,8 +1258,8 @@ inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
   double rotateDeg = 0;
 
   // rotate in 1 deg steps
-  for (int i = 1; i < 360; i += 1) {
+  for (double i = step; i < 360; i += step) {
-    pol = rotate(pol, 1, center);
+    pol = rotate(pol, step, center);
     Box<T> e = getBoundingBox(pol);
     if (area(tmpBox) > area(e)) {
       tmpBox = e;
@@ -1247,6 +1270,50 @@ inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
   return RotatedBox<T>(tmpBox, -rotateDeg, center);
 }
 
+// _____________________________________________________________________________
+template <template <typename> class Geometry, typename T>
+inline RotatedBox<T> getOrientedEnvelope(std::vector<Geometry<T>> pol) {
+  return getOrientedEnvelope(pol, 1);
+}
+
+// _____________________________________________________________________________
+template <template <typename> class Geometry, typename T>
+inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol, double step) {
+  std::vector<Geometry<T>> mult{pol};
+  return getOrientedEnvelope(mult, step);
+}
+
+// _____________________________________________________________________________
+template <template <typename> class Geometry, typename T>
+inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
+  std::vector<Geometry<T>> mult{pol};
+  return getOrientedEnvelope(mult, 1);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> buffer(const Line<T>& line, double d, size_t points) {
+  // so far only works correctly if the input polygon is convex
+  MultiPoint<T> pSet;
+  for (const auto& p : line) {
+    Point<T> anchor{p.getX() + d, p.getY()};
+    double deg = 0;
+    pSet.push_back(p);
+    for (size_t i = 0; i < points; i++) {
+      pSet.push_back(rotate(anchor, deg, p));
+      deg += 360 / (1.0 * points);
+    }
+  }
+
+  return convexHull(pSet);
+}
+
+// _____________________________________________________________________________
+template <typename T>
+inline Polygon<T> buffer(const Polygon<T>& pol, double d, size_t points) {
+  return buffer(pol.getOuter(), d, points);
+}
+
 // _____________________________________________________________________________
 template <typename T>
 inline Box<T> extendBox(const Box<T>& a, Box<T> b) {
@@ -1799,7 +1866,6 @@ inline double haversine(const Point<T>& a, const Point<T>& b) {
   return haversine(a.getY(), a.getX(), b.getY(), b.getX());
 }
 
-
 // _____________________________________________________________________________
 template <typename T>
 inline Line<T> densify(const Line<T>& l, double d) {
@@ -1963,8 +2029,8 @@ inline double accFrechetDistCHav(const Line<T>& a, const Line<T>& b, double d) {
 
   for (size_t i = 1; i < p.size(); i++) {
     for (size_t j = 1; j < q.size(); j++) {
-      float d =
+      float d = util::geo::haversine(p[i], q[j]) *
-          util::geo::haversine(p[i], q[j]) * util::geo::haversine(p[i], p[i - 1]);
+                util::geo::haversine(p[i], p[i - 1]);
       ca[i * q.size() + j] =
           d + std::min(ca[(i - 1) * q.size() + j],
                        std::min(ca[i * q.size() + (j - 1)],
@@ -2042,7 +2108,6 @@ inline double latLngLen(const Line<T>& g) {
   return ret;
 }
 
-
 // _____________________________________________________________________________
 template <typename G>
 inline double webMercDistFactor(const G& a) {

src/util/geo/PolyLine.h

@@ -71,8 +71,6 @@ class PolyLine {
   double distTo(const PolyLine<T>& g) const;
   double distTo(const Point<T>& p) const;
 
-  SharedSegments<T> getSharedSegments(const PolyLine<T>& pl, double dmax) const;
-
   double getLength() const;
 
   // return point at dist

src/util/geo/PolyLine.tpp

@@ -188,9 +188,9 @@ PolyLine<T> PolyLine<T>::getSegment(const LinePoint<T>& start,
   ret << end.p;
 
   // find a more performant way to clear the result of above
-  // ret.simplify(0);
+  ret.simplify(0);
 
-  // assert(ret.getLine().size());
+  assert(ret.getLine().size());
 
   return ret;
 }
@@ -478,120 +478,6 @@ void PolyLine<T>::move(double vx, double vy) {
   }
 }
 
-// _____________________________________________________________________________
-template <typename T>
-SharedSegments<T> PolyLine<T>::getSharedSegments(const PolyLine<T>& pl,
-                                                 double dmax) const {
-  /**
-   * Returns the segments this polyline share with pl
-   * atm, this is a very simple distance-based algorithm
-   */
-  double STEP_SIZE = 2;
-  double MAX_SKIPS = 4;
-  double MIN_SEG_LENGTH = 0.1;  // dmax / 2;  // make this configurable!
-
-  SharedSegments<T> ret;
-
-  if (distTo(pl) > dmax) return ret;
-
-  bool in = false, single = true;
-  double curDist = 0;
-  double curTotalSegDist = 0;
-  size_t skips;
-
-  LinePoint<T> curStartCand, curEndCand, curStartCandCmp, curEndCandCmp;
-
-  double comp = 0, curSegDist = 0;
-  double length = getLength(), plLength = pl.getLength();
-
-  for (size_t i = 1; i < _line.size(); ++i) {
-    const Point<T>& s = _line[i - 1];
-    const Point<T>& e = _line[i];
-
-    bool lastRound = false;
-
-    double totalDist = dist(s, e);
-    while (curSegDist <= totalDist) {
-      const auto& curPointer = interpolate(s, e, curSegDist);
-
-      if (pl.distTo(curPointer) <= dmax) {
-        LinePoint<T> curCmpPointer = pl.projectOn(curPointer);
-        LinePoint<T> curBackProjectedPointer = projectOn(curCmpPointer.p);
-
-        skips = 0;
-
-        if (in) {
-          curEndCand = curBackProjectedPointer;
-          curEndCandCmp = curCmpPointer;
-
-          if (curEndCand.totalPos < curStartCand.totalPos) {
-            curEndCand = curStartCand;
-          }
-
-          single = false;
-
-          comp = fabs(curStartCand.totalPos * length -
-                      curEndCand.totalPos * length) /
-                 fabs(curStartCandCmp.totalPos * plLength -
-                      curEndCandCmp.totalPos * plLength);
-        } else {
-          in = true;
-          curStartCand = curBackProjectedPointer;
-          curStartCandCmp = curCmpPointer;
-        }
-      } else {
-        if (in) {
-          skips++;
-          if (skips > MAX_SKIPS) {  // TODO: make configurable
-            if (comp > 0.8 && comp < 1.2 && !single &&
-                (fabs(curStartCand.totalPos * length -
-                      curEndCand.totalPos * length) > MIN_SEG_LENGTH &&
-                 fabs(curStartCandCmp.totalPos * plLength -
-                      curEndCandCmp.totalPos * plLength) > MIN_SEG_LENGTH)) {
-              ret.segments.push_back(
-                  SharedSegment<T>(std::pair<LinePoint<T>, LinePoint<T>>(
-                                       curStartCand, curStartCandCmp),
-                                   std::pair<LinePoint<T>, LinePoint<T>>(
-                                       curEndCand, curEndCandCmp)));
-
-              // TODO: only return the FIRST one, make this configuralbe
-              return ret;
-            }
-
-            in = false;
-            single = true;
-          }
-        }
-      }
-
-      if (curSegDist + STEP_SIZE > totalDist && !lastRound) {
-        lastRound = true;
-        double finalStep = totalDist - curSegDist - 0.0005;
-        curSegDist += finalStep;
-        curDist += finalStep;
-      } else {
-        curSegDist += STEP_SIZE;
-        curDist += STEP_SIZE;
-      }
-    }
-
-    curSegDist = curSegDist - totalDist;
-    curTotalSegDist += totalDist;
-  }
-
-  if (comp > 0.8 && comp < 1.2 && in && !single &&
-      (fabs(curStartCand.totalPos * length - curEndCand.totalPos * length) >
-           MIN_SEG_LENGTH &&
-       fabs(curStartCandCmp.totalPos * plLength -
-            curEndCandCmp.totalPos * plLength) > MIN_SEG_LENGTH)) {
-    ret.segments.push_back(SharedSegment<T>(
-        std::pair<LinePoint<T>, LinePoint<T>>(curStartCand, curStartCandCmp),
-        std::pair<LinePoint<T>, LinePoint<T>>(curEndCand, curEndCandCmp)));
-  }
-
-  return ret;
-}
-
 // _____________________________________________________________________________
 template <typename T>
 std::set<LinePoint<T>, LinePointCmp<T>> PolyLine<T>::getIntersections(

src/util/graph/DirGraph.h

@@ -28,7 +28,7 @@ class DirGraph : public Graph<N, E> {
   Node<N, E>* addNd(const N& pl);
   Edge<N, E>* addEdg(Node<N, E>* from, Node<N, E>* to, const E& p);
 
-  Node<N, E>* mergeNds(Node<N, E>* a, Node<N, E>* b);
+  virtual Node<N, E>* mergeNds(Node<N, E>* a, Node<N, E>* b);
 
 };
 

src/util/graph/EDijkstra.h

@@ -114,16 +114,18 @@ class EDijkstra : public ShortestPath<EDijkstra> {
   using SettledInit = tsl::robin_map<Edge<N, E>*, RouteEdgeInit<N, E, C>>;
 
   template <typename N, typename E, typename C>
-  using SettledInitNoRes = tsl::robin_map<Edge<N, E>*, RouteEdgeInitNoRes<N, E, C>>;
+  using SettledInitNoRes =
+      tsl::robin_map<Edge<N, E>*, RouteEdgeInitNoRes<N, E, C>>;
 
   template <typename N, typename E, typename C>
   using PQInit = radix_heap::pair_radix_heap<C, RouteEdgeInit<N, E, C>>;
 
   template <typename N, typename E, typename C>
-  using PQInitNoRes = radix_heap::pair_radix_heap<C, RouteEdgeInitNoRes<N, E, C>>;
+  using PQInitNoRes =
+      radix_heap::pair_radix_heap<C, RouteEdgeInitNoRes<N, E, C>>;
 
   template <typename N, typename E, typename C>
-  static C shortestPathImpl(const std::set<Edge<N, E>*> from,
+  static C shortestPathImpl(const std::set<Edge<N, E>*>& from,
                             const std::set<Edge<N, E>*>& to,
                             const util::graph::CostFunc<N, E, C>& costFunc,
                             const util::graph::HeurFunc<N, E, C>& heurFunc,
@@ -148,6 +150,13 @@ class EDijkstra : public ShortestPath<EDijkstra> {
                             const util::graph::HeurFunc<N, E, C>& heurFunc,
                             EList<N, E>* resEdges, NList<N, E>* resNodes);
 
+  template <typename N, typename E, typename C>
+  static C shortestPathImpl(const std::set<Node<N, E>*>& from,
+                            const std::set<Node<N, E>*>& to,
+                            const util::graph::CostFunc<N, E, C>& costFunc,
+                            const util::graph::HeurFunc<N, E, C>& heurFunc,
+                            EList<N, E>* resEdges, NList<N, E>* resNodes);
+
   template <typename N, typename E, typename C>
   static std::unordered_map<Edge<N, E>*, C> shortestPathImpl(
       const std::set<Edge<N, E>*>& from,
@@ -161,6 +170,15 @@ class EDijkstra : public ShortestPath<EDijkstra> {
       std::unordered_map<Edge<N, E>*, EList<N, E>*> resEdges,
       std::unordered_map<Edge<N, E>*, NList<N, E>*> resNodes);
 
+  template <typename N, typename E, typename C>
+  static std::unordered_map<Edge<N, E>*, C> shortestPathImpl(
+      const std::set<Edge<N, E>*>& from, const std::set<Edge<N, E>*>& to,
+      const std::unordered_map<Edge<N, E>*, C>& initCosts,
+      const util::graph::CostFunc<N, E, C>& costFunc,
+      const util::graph::HeurFunc<N, E, C>& heurFunc,
+      std::unordered_map<Edge<N, E>*, EList<N, E>*> resEdges,
+      std::unordered_map<Edge<N, E>*, NList<N, E>*> resNodes);
+
   template <typename N, typename E, typename C>
   static std::unordered_map<Edge<N, E>*, C> shortestPathImpl(
       const std::set<Edge<N, E>*>& from, const std::set<Edge<N, E>*>& to,
@@ -202,11 +220,10 @@ class EDijkstra : public ShortestPath<EDijkstra> {
                                PQInit<N, E, C>& pq);
 
   template <typename N, typename E, typename C>
-  static inline void relaxInitNoResEdgs(RouteEdgeInitNoRes<N, E, C>& cur,
+  static inline void relaxInitNoResEdgs(
-                               const std::set<Edge<N, E>*>& to, C stall,
+      RouteEdgeInitNoRes<N, E, C>& cur, const std::set<Edge<N, E>*>& to,
-                               const util::graph::CostFunc<N, E, C>& costFunc,
+      C stall, const util::graph::CostFunc<N, E, C>& costFunc,
-                               const util::graph::HeurFunc<N, E, C>& heurFunc,
+      const util::graph::HeurFunc<N, E, C>& heurFunc, PQInitNoRes<N, E, C>& pq);
-                               PQInitNoRes<N, E, C>& pq);
 
   template <typename N, typename E, typename C>
   static void relaxInv(RouteEdge<N, E, C>& cur,

src/util/graph/EDijkstra.tpp

@@ -47,7 +47,30 @@ C EDijkstra::shortestPathImpl(Edge<N, E>* from, const std::set<Node<N, E>*>& to,
 
 // _____________________________________________________________________________
 template <typename N, typename E, typename C>
-C EDijkstra::shortestPathImpl(const std::set<Edge<N, E>*> from,
+C EDijkstra::shortestPathImpl(const std::set<Node<N, E>*>& from,
+                              const std::set<Node<N, E>*>& to,
+                              const util::graph::CostFunc<N, E, C>& costFunc,
+                              const util::graph::HeurFunc<N, E, C>& heurFunc,
+                              EList<N, E>* resEdges, NList<N, E>* resNodes) {
+  std::set<Edge<N, E>*> frEs;
+  std::set<Edge<N, E>*> toEs;
+
+  for (auto n : from) {
+    frEs.insert(n->getAdjListIn().begin(), n->getAdjListIn().end());
+  }
+
+  for (auto n : to) {
+    toEs.insert(n->getAdjListIn().begin(), n->getAdjListIn().end());
+  }
+
+  C cost = shortestPathImpl(frEs, toEs, costFunc, heurFunc, resEdges, resNodes);
+
+  return cost;
+}
+
+// _____________________________________________________________________________
+template <typename N, typename E, typename C>
+C EDijkstra::shortestPathImpl(const std::set<Edge<N, E>*>& from,
                               const std::set<Edge<N, E>*>& to,
                               const util::graph::CostFunc<N, E, C>& costFunc,
                               const util::graph::HeurFunc<N, E, C>& heurFunc,

src/util/graph/Graph.h

@@ -5,10 +5,10 @@
 #ifndef UTIL_GRAPH_GRAPH_H_
 #define UTIL_GRAPH_GRAPH_H_
 
+#include <cassert>
+#include <iostream>
 #include <set>
 #include <string>
-#include <iostream>
-#include <cassert>
 
 #include "util/graph/Edge.h"
 #include "util/graph/Node.h"
@@ -25,7 +25,7 @@ class Graph {
   Edge<N, E>* addEdg(Node<N, E>* from, Node<N, E>* to);
   virtual Edge<N, E>* addEdg(Node<N, E>* from, Node<N, E>* to, const E& p) = 0;
   Edge<N, E>* getEdg(Node<N, E>* from, Node<N, E>* to);
-  const Edge<N, E>* getEdg(Node<N, E>* from, Node<N, E>* to) const;
+  const Edge<N, E>* getEdg(const Node<N, E>* from, const Node<N, E>* to) const;
 
   virtual Node<N, E>* mergeNds(Node<N, E>* a, Node<N, E>* b) = 0;
 
@@ -43,7 +43,7 @@ class Graph {
 };
 
 #include "util/graph/Graph.tpp"
-}
+}  // namespace graph
-}
+}  // namespace util
 
 #endif  // UTIL_GRAPH_GRAPH_H_

src/util/graph/Graph.tpp

@@ -22,8 +22,7 @@ const std::set<Node<N, E>*>& Graph<N, E>::getNds() const {
 
 // _____________________________________________________________________________
 template <typename N, typename E>
-typename std::set<Node<N, E>*>::iterator Graph<N, E>::delNd(
+typename std::set<Node<N, E>*>::iterator Graph<N, E>::delNd(Node<N, E>* n) {
-    Node<N, E>* n) {
   return delNd(_nodes.find(n));
 }
 
@@ -69,10 +68,10 @@ Node<N, E>* Graph<N, E>::sharedNode(const Edge<N, E>* a, const Edge<N, E>* b) {
   return r;
 }
 
-
 // _____________________________________________________________________________
 template <typename N, typename E>
-const Edge<N, E>* Graph<N, E>::getEdg(Node<N, E>* from, Node<N, E>* to) const {
+const Edge<N, E>* Graph<N, E>::getEdg(const Node<N, E>* from,
+                                      const Node<N, E>* to) const {
   for (auto e : from->getAdjList()) {
     if (e->getOtherNd(from) == to) return e;
   }

src/util/graph/ShortestPath.h

@@ -79,6 +79,16 @@ class ShortestPath {
                                resNodes);
   }
 
+  template <typename N, typename E, typename C>
+  static C shortestPath(const std::set<Node<N, E>*>& from,
+                        const std::set<Node<N, E>*>& to,
+                        const CostFunc<N, E, C>& costFunc) {
+    EList<N, E>* el = 0;
+    NList<N, E>* nl = 0;
+    return D::shortestPathImpl(from, to, costFunc, ZeroHeurFunc<N, E, C>(),
+                               el, nl);
+  }
+
   template <typename N, typename E, typename C>
   static C shortestPath(const std::set<Node<N, E>*> from,
                         const std::set<Node<N, E>*>& to,

src/util/graph/UndirGraph.h

@@ -30,7 +30,7 @@ class UndirGraph : public Graph<N, E> {
   Node<N, E>* addNd(const N& pl);
   Edge<N, E>* addEdg(Node<N, E>* from, Node<N, E>* to, const E& p);
 
-  Node<N, E>* mergeNds(Node<N, E>* a, Node<N, E>* b);
+  virtual Node<N, E>* mergeNds(Node<N, E>* a, Node<N, E>* b);
 
 };
 

src/util/graph/UndirGraph.tpp

@@ -42,12 +42,14 @@ Edge<N, E>* UndirGraph<N, E>::addEdg(Node<N, E>* from, Node<N, E>* to,
 template <typename N, typename E>
 Node<N, E>* UndirGraph<N, E>::mergeNds(Node<N, E>* a, Node<N, E>* b) {
   for (auto e : a->getAdjListOut()) {
+    if (e->getFrom() != a) continue;
     if (e->getTo() != b) {
       addEdg(b, e->getTo(), e->pl());
     }
   }
 
   for (auto e : a->getAdjListIn()) {
+    if (e->getTo() != a) continue;
     if (e->getFrom() != b) {
       addEdg(e->getFrom(), b, e->pl());
     }

src/util/graph/radix_heap.h

@@ -128,7 +128,7 @@ class pair_radix_heap {
   void push(key_type key, const value_type &value) {
     unsigned_key_type x = encoder_type::encode(key);
     if (last_ > x) {
-      std::cerr << "Not monotone: " << last_ << " vs " << x << std::endl;
+      std::cerr << "PQ: not monotone: " << last_ << " vs " << x << std::endl;
       x = last_;
     }
     ++size_;
@@ -140,7 +140,7 @@ class pair_radix_heap {
   void push(key_type key, value_type &&value) {
     unsigned_key_type x = encoder_type::encode(key);
     if (last_ > x) {
-      std::cerr << "Not monotone: " << last_ << " vs " << x << std::endl;
+      std::cerr << "PQ: not monotone: " << last_ << " vs " << x << std::endl;
       x = last_;
     }
     ++size_;

src/util/json/Writer.h

@@ -30,9 +30,9 @@ struct Null {};
 struct Val {
   enum VAL_T { UINT, INT, FLOAT, STRING, ARRAY, DICT, BOOL, JSNULL };
   VAL_T type;
-  int i;
+  int i = 0;
-  uint64_t ui;
+  uint64_t ui = 0;
-  double f;
+  double f = 0;
   std::string str;
   std::vector<Val> arr;
   std::map<std::string, Val> dict;

src/util/tests/TestMain.cpp

@@ -1843,7 +1843,6 @@ int main(int argc, char** argv) {
   TEST(geo::getWKT(geo::segment(DLine{{0, 0}, {0, 1}, {0, 2}}, 0.5, 1)), ==, "LINESTRING (0 1, 0 2)");
 }
 
-
   // inversion count
   std::vector<int> test = {2, 1};
   TEST(inversions(test), ==, 1);
@@ -1880,4 +1879,13 @@ int main(int argc, char** argv) {
 
   test = {9, 8, 7, 6, 5, 4, 3, 2, 1};
   TEST(inversions(test), ==, 8 + 7 + 6 + 5 + 4 + 3 + 2 + 1);
+
+  // nice float formatting
+	TEST(formatFloat(15.564, 3), ==, "15.564");
+	TEST(formatFloat(15.564, 0), ==, "16");
+	TEST(formatFloat(15.0000, 10), ==, "15");
+	TEST(formatFloat(15.0100, 10), ==, "15.01");
+	TEST(formatFloat(0.0000, 10), ==, "0");
+	TEST(formatFloat(-1.0000, 10), ==, "-1");
+	TEST(formatFloat(-15.000001, 10), ==, "-15.000001");
 }