diff --git a/.gitignore b/.gitignore
index e5913f35..57b65da9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -43,14 +43,11 @@ vs_build/.vs
 *.lib
 
 # Folders
-/vendor
-/_site
-x64/
-x86/
-
-#Folders
-/vendor
-/_site
+x64
+x86
+vendor
+_site
+.sass-cache
 /resources/convert-to-simulation/simulatorready
 /resources/convert-to-simulation/out.off
 /resources/convert-to-simulation/bin/
diff --git a/CMakeLists.txt b/CMakeLists.txt
index dd76ecaa..207e38d9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -34,9 +34,10 @@ if ( NOT GDAL_FOUND )
   message(SEND_ERROR "3dfier requires the GDAL library")
 endif()
 
-# LIBLAS
-find_package(libLAS REQUIRED)
-find_package (LASzip QUIET)
+# LASLIB
+# find_package(libLAS REQUIRED)
+# find_package (LASzip QUIET)
+find_package (LASlib REQUIRED)
 
 # YamlCpp
 find_package(YamlCpp REQUIRED)
@@ -53,15 +54,15 @@ include( ${CGAL_USE_FILE} )
 
 include_directories( ${CMAKE_SOURCE_DIR}/thirdparty )
 
-include_directories( ${CGAL_INCLUDE_DIR} ${CGAL_3RD_PARTY_INCLUDE_DIR} ${Boost_INCLUDE_DIRS} ${LIBLAS_INCLUDE_DIR} ${LASZIP_INCLUDE_DIR} ${YAMLCPP_INCLUDE_DIR} ${GDAL_INCLUDE_DIR})
+include_directories( ${CGAL_INCLUDE_DIR} ${CGAL_3RD_PARTY_INCLUDE_DIR} ${Boost_INCLUDE_DIRS} ${YAMLCPP_INCLUDE_DIR} ${GDAL_INCLUDE_DIR} ${LASLIB_INCLUDE_DIR})
 link_directories(${YamlCpp_LIBRARY_DIRS})
 
-# build ptinpoly as static lib
-add_library(ptinpoly STATIC thirdparty/ptinpoly.c)
-set_target_properties(
-  ptinpoly
-  PROPERTIES C_STANDARD 11
-)
+# # build ptinpoly as static lib
+# add_library(ptinpoly STATIC thirdparty/ptinpoly.c)
+# set_target_properties(
+#   ptinpoly
+#   PROPERTIES C_STANDARD 11
+# )
 
 # Creating entries for target: 3dfier
 FILE(GLOB SRC_FILES src/*.cpp)
@@ -71,6 +72,6 @@ set_target_properties(
   PROPERTIES CXX_STANDARD 11
 )
 
-target_link_libraries( 3dfier ${CGAL_LIBRARIES} ${CGAL_3RD_PARTY_LIBRARIES} ${GDAL_LIBRARY} ${LIBLAS_LIBRARY} ${LASZIP_LIBRARY} ${YAMLCPP_LIBRARY} Boost::program_options Boost::filesystem Boost::locale ptinpoly)
+target_link_libraries( 3dfier ${CGAL_LIBRARIES} ${CGAL_3RD_PARTY_LIBRARIES} ${GDAL_LIBRARY} ${YAMLCPP_LIBRARY} Boost::program_options Boost::filesystem Boost::locale ${LASLIB_LIBRARY})
 
 install(TARGETS 3dfier DESTINATION bin)
\ No newline at end of file
diff --git a/cmake/Modules/FindLASlib.cmake b/cmake/Modules/FindLASlib.cmake
new file mode 100644
index 00000000..fd9d22a1
--- /dev/null
+++ b/cmake/Modules/FindLASlib.cmake
@@ -0,0 +1,45 @@
+###############################################################################
+#
+# CMake module to search for LASlib library
+#
+# On success, the macro sets the following variables:
+# LASLIB_FOUND       = if the library found
+# LASLIB_LIBRARY   = full path to the library
+# LASLIB_INCLUDE_DIR = where to find the library headers also defined,
+#                       but not for general use are
+# LASLIB_LIBRARY     = where to find the library.
+# LASLIB_VERSION     = version of library which was found, e.g. "1.2.5"
+#
+###############################################################################
+MESSAGE(STATUS "Searching for LASlib library")
+
+IF(LASLIB_INCLUDE_DIR)
+  # Already in cache, be silent
+  SET(LASLIB_FIND_QUIETLY TRUE)
+ENDIF()
+
+# IF(APPLE)
+#   SET(LASLIB_ROOT_DIR /usr/local/)
+#   MESSAGE(STATUS "Trying LasZip using default location LASLIB_ROOT=c:/laszip")
+# ENDIF()
+
+FIND_PATH(LASLIB_INCLUDE_DIR laszip.hpp
+  PATH_PREFIXES laszip
+  PATHS
+  /usr/include
+  /usr/local/include
+  /usr/local/include/LASlib
+  # ${LASLIB_ROOT_DIR}/include
+  NO_DEFAULT_PATH)
+
+SET(LASLIB_NAMES ${LASLIB_LIBRARY} LASlib)
+
+FIND_LIBRARY(LASLIB_LIBRARY
+  NAMES ${LASLIB_NAMES}
+  PATHS
+  /usr/lib
+  /usr/local/lib
+  /usr/local/lib/LASlib
+  ${LASLIB_ROOT_DIR}/lib)
+
+
diff --git a/cmake/Modules/FindLASzip.cmake b/cmake/Modules/FindLASzip.cmake
deleted file mode 100644
index 70d1a4be..00000000
--- a/cmake/Modules/FindLASzip.cmake
+++ /dev/null
@@ -1,56 +0,0 @@
-###############################################################################
-#
-# CMake module to search for LASzip library
-#
-# On success, the macro sets the following variables:
-# LASZIP_FOUND       = if the library found
-# LASZIP_LIBRARIES   = full path to the library
-# LASZIP_INCLUDE_DIR = where to find the library headers also defined,
-#                       but not for general use are
-# LASZIP_LIBRARY     = where to find the PROJ.4 library.
-# LASZIP_VERSION     = version of library which was found, e.g. "1.2.5"
-#
-# Copyright (c) 2009 Mateusz Loskot <mateusz@loskot.net>
-#
-# Module source: http://github.com/mloskot/workshop/tree/master/cmake/
-#
-# Redistribution and use is allowed according to the terms of the BSD license.
-# For details see the accompanying COPYING-CMAKE-SCRIPTS file.
-#
-###############################################################################
-MESSAGE(STATUS "Searching for LASzip ${LASzip_FIND_VERSION}+ library")
-
-IF(LASZIP_INCLUDE_DIR)
-  # Already in cache, be silent
-  SET(LASZIP_FIND_QUIETLY TRUE)
-ENDIF()
-
-IF(WIN32)
-  IF(DEFINED ENV{LASZIP_ROOT})
-    SET(LASZIP_ROOT_DIR $ENV{LASZIP_ROOT})
-    MESSAGE(STATUS "Trying LasZip using environment variable LASZIP_ROOT=$ENV{LASZIP_ROOT}")
-  ELSE()
-    SET(LASZIP_ROOT_DIR c:/laszip)
-    MESSAGE(STATUS "Trying LasZip using default location LASZIP_ROOT=c:/laszip")
-  ENDIF()
-ENDIF()
-
-
-FIND_PATH(LASZIP_INCLUDE_DIR laszip.hpp
-  PATH_PREFIXES laszip
-  PATHS
-  /usr/include
-  /usr/local/include
-  ${LASZIP_ROOT_DIR}/include
-  NO_DEFAULT_PATH)
-
-SET(LASZIP_NAMES ${LASZIP_LIBRARY} laszip)
-
-FIND_LIBRARY(LASZIP_LIBRARY
-  NAMES ${LASZIP_NAMES}
-  PATHS
-  /usr/lib
-  /usr/local/lib
-  ${LASZIP_ROOT_DIR}/lib)
-
-
diff --git a/cmake/Modules/FindlibLAS.cmake b/cmake/Modules/FindlibLAS.cmake
deleted file mode 100644
index 7d6ce712..00000000
--- a/cmake/Modules/FindlibLAS.cmake
+++ /dev/null
@@ -1,64 +0,0 @@
-###############################################################################
-#
-# CMake module to search for libLAS library
-#
-# On success, the macro sets the following variables:
-# LIBLAS_FOUND       = if the library found
-# LIBLAS_LIBRARIES   = full path to the library
-# LIBLAS_INCLUDE_DIR = where to find the library headers also defined,
-#                       but not for general use are
-# LIBLAS_LIBRARY     = where to find the PROJ.4 library.
-# LIBLAS_VERSION     = version of library which was found, e.g. "1.2.5"
-#
-# Copyright (c) 2009 Mateusz Loskot <mateusz@loskot.net>
-#
-# Module source: http://github.com/mloskot/workshop/tree/master/cmake/
-#
-# Redistribution and use is allowed according to the terms of the BSD license.
-# For details see the accompanying COPYING-CMAKE-SCRIPTS file.
-#
-###############################################################################
-MESSAGE(STATUS "Searching for LibLAS ${LibLAS_FIND_VERSION}+ library")
-
-IF(LIBLAS_INCLUDE_DIR)
-  # Already in cache, be silent
-  SET(LIBLAS_FIND_QUIETLY TRUE)
-ENDIF()
-
-IF(WIN32)
-  IF(DEFINED ENV{LIBLAS_ROOT})
-    SET(LIBLAS_ROOT_DIR $ENV{LIBLAS_ROOT})
-    #MESSAGE(STATUS " FindLibLAS: trying LIBLAS using environment variable LIBLAS_ROOT=$ENV{LIBLAS_ROOT}")
-  ELSE()
-    SET(LIBLAS_ROOT_DIR c:/liblas)
-    #MESSAGE(STATUS " FindLibLAS: trying LIBLAS using default location LIBLAS_ROOT=c:/liblas")
-  ENDIF()
-ENDIF()
-
-
-FIND_PATH(LIBLAS_INCLUDE_DIR
-  liblas.hpp
-  PATH_PREFIXES liblas
-  PATHS
-  /usr/include
-  /usr/local/include
-  /tmp/lasjunk/include
-  ${LIBLAS_ROOT_DIR}/include)
-
-find_library(LIBLAS_LIBRARY
-  NAMES liblas.dylib liblas
-  PATHS
-  /usr/lib
-  /usr/local/lib
-  /tmp/lasjunk/lib
-  ${LIBLAS_ROOT_DIR}/lib
-)
-
-
-
-
-
-# Handle the QUIETLY and REQUIRED arguments and set LIBLAS_FOUND to TRUE
-# if all listed variables are TRUE
-#INCLUDE(FindPackageHandleStandardArgs)
-#FIND_PACKAGE_HANDLE_STANDARD_ARGS(libLAS DEFAULT_MSG LIBLAS_LIBRARY LIBLAS_INCLUDE_DIR)
\ No newline at end of file
diff --git a/resources/config_files/myconfig.yml b/resources/config_files/myconfig.yml
index 6f8de9eb..8eb99a64 100644
--- a/resources/config_files/myconfig.yml
+++ b/resources/config_files/myconfig.yml
@@ -96,3 +96,4 @@ options:
   building_radius_vertex_elevation: 3.0
   radius_vertex_elevation: 1.0
   threshold_jump_edges: 0.5
+  threshold_bridge_jump_edges: 0.5
diff --git a/resources/config_files/myconfig_DEFAULTS.yml b/resources/config_files/myconfig_DEFAULTS.yml
index d1704132..9fd1c091 100644
--- a/resources/config_files/myconfig_DEFAULTS.yml
+++ b/resources/config_files/myconfig_DEFAULTS.yml
@@ -36,4 +36,5 @@ input_elevation:
 options:
   building_radius_vertex_elevation: 3.0
   radius_vertex_elevation: 1.0
-  threshold_jump_edges: 0.5
\ No newline at end of file
+  threshold_jump_edges: 0.5
+  threshold_bridge_jump_edges: 0.5
\ No newline at end of file
diff --git a/resources/config_files/myconfig_README.yml b/resources/config_files/myconfig_README.yml
index d7fa1ff1..f3dcc8df 100644
--- a/resources/config_files/myconfig_README.yml
+++ b/resources/config_files/myconfig_README.yml
@@ -90,4 +90,5 @@ options:                                                # Global options
   building_radius_vertex_elevation: 3.0                 # Radius in meters used for point-vertex distance between 3D points and vertices of building polygons, radius_vertex_elevation used when not specified
   radius_vertex_elevation: 1.0                          # Radius in meters used for point-vertex distance between 3D points and vertices of polygons
   threshold_jump_edges: 0.5                             # Threshold in meters for stitching adjacent objects, when the height difference is larger then the threshold a vertical wall is created 
+  threshold_bridge_jump_edges: 0.5                      # Threshold in meters for stitching bridges to adjacent objects, if not specified it falls back to threshold_jump_edges
   extent: xmin, ymin, xmax, ymax                        # Filter the input polygons to this extent
diff --git a/src/Bridge.cpp b/src/Bridge.cpp
index 5baef696..5d81e76d 100644
--- a/src/Bridge.cpp
+++ b/src/Bridge.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Bridge.h b/src/Bridge.h
index 1d0cc388..cf683b64 100644
--- a/src/Bridge.h
+++ b/src/Bridge.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Building.cpp b/src/Building.cpp
index c38eefed..d5a91dfa 100644
--- a/src/Building.cpp
+++ b/src/Building.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Building.h b/src/Building.h
index dc3215fd..159fb9dd 100644
--- a/src/Building.h
+++ b/src/Building.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Forest.cpp b/src/Forest.cpp
index 4c0e39b1..cff77662 100644
--- a/src/Forest.cpp
+++ b/src/Forest.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Forest.h b/src/Forest.h
index c031fe51..1cc3558f 100644
--- a/src/Forest.h
+++ b/src/Forest.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Map3d.cpp b/src/Map3d.cpp
index f6740dda..aef63fd9 100644
--- a/src/Map3d.cpp
+++ b/src/Map3d.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
@@ -51,8 +51,13 @@ Map3d::Map3d() {
   _radius_vertex_elevation = 1.0;
   _building_radius_vertex_elevation = 3.0;
   _threshold_jump_edges = 50;
+  _threshold_bridge_jump_edges = 50;
   _requestedExtent = Box2(Point2(0, 0), Point2(0, 0));
   _bbox = Box2(Point2(999999, 999999), Point2(-999999, -999999));
+  _minxradius = 999999;
+  _maxxradius = 999999;
+  _minyradius = -999999;
+  _maxyradius = -999999;
 }
 
 Map3d::~Map3d() {
@@ -79,6 +84,10 @@ void Map3d::set_threshold_jump_edges(float threshold) {
   _threshold_jump_edges = int(threshold * 100);
 }
 
+void Map3d::set_threshold_bridge_jump_edges(float threshold) {
+  _threshold_bridge_jump_edges = int(threshold * 100);
+}
+
 void Map3d::set_building_include_floor(bool include) {
   _building_include_floor = include;
 }
@@ -155,10 +164,12 @@ Box2 Map3d::get_bbox() {
   return _bbox;
 }
 
-liblas::Bounds<double> Map3d::get_bounds() {
-  double radius = std::max(_radius_vertex_elevation, _building_radius_vertex_elevation);
-  liblas::Bounds<double> bounds(_bbox.min_corner().x() - radius, _bbox.min_corner().y() - radius, _bbox.max_corner().x() + radius, _bbox.max_corner().y() + radius);
-  return bounds;
+bool Map3d::check_bounds(const double xmin, const double xmax, const double ymin, const double ymax) {
+  if ((xmin < _maxxradius || xmax > _minxradius) &&
+    (ymin < _maxyradius || ymax > _minyradius)) {
+    return true;
+  }
+  return false;
 }
 
 bool Map3d::get_cityjson(std::string filename) {
@@ -190,8 +201,7 @@ bool Map3d::get_cityjson(std::string filename) {
     j["vertices"].push_back({std::stod(c[0], NULL), std::stod(c[1], NULL), std::stod(c[2], NULL) });
   }
   std::ofstream o(filename);
-  // o << j.dump(2) << std::endl;      
-  o << j.dump() << std::endl;      
+  o << j.dump() << std::endl;
   return true;
 }
 
@@ -401,7 +411,7 @@ void Map3d::get_obj_per_class(std::wostream& of) {
   of << fs << std::endl;
 }
 
-bool Map3d::get_pdok_output(std::string filename) {
+bool Map3d::get_postgis_output(std::string connstr, bool pdok, bool citygml) {
 #if GDAL_VERSION_MAJOR < 2
   std::cerr << "ERROR: cannot write MultiPolygonZ files with GDAL < 2.0.\n";
   return false;
@@ -409,7 +419,7 @@ bool Map3d::get_pdok_output(std::string filename) {
   if (GDALGetDriverCount() == 0)
     GDALAllRegister();
   GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("PostgreSQL");
-  GDALDataset* dataSource = driver->Create(filename.c_str(), 0, 0, 0, GDT_Unknown, NULL);
+  GDALDataset* dataSource = driver->Create(connstr.c_str(), 0, 0, 0, GDT_Unknown, NULL);
   if (dataSource == NULL) {
     std::cerr << "Starting database connection failed.\n";
     return false;
@@ -425,11 +435,13 @@ bool Map3d::get_pdok_output(std::string filename) {
     std::string layername = f->get_layername();
     if (layers.find(layername) == layers.end()) {
       AttributeMap &attributes = f->get_attributes();
-      //Add additional attribute to list for layer creation
-      attributes["xml"] = std::make_pair(OFTString, "");
-      OGRLayer *layer = create_gdal_layer(driver, dataSource, filename, layername, attributes, f->get_class() == BUILDING);
+      if (pdok) {
+        //Add additional attribute to list for layer creation
+        attributes["xml"] = std::make_pair(OFTString, "");
+      }
+      OGRLayer *layer = create_gdal_layer(driver, dataSource, connstr, layername, attributes, f->get_class() == BUILDING);
       if (layer == NULL) {
-        std::cerr << "ERROR: Cannot open database '" + filename + "' for writing" << std::endl;
+        std::cerr << "ERROR: Cannot open database '" + connstr + "' for writing" << std::endl;
         dataSource->RollbackTransaction();
         GDALClose(dataSource);
         GDALClose(driver);
@@ -451,101 +463,23 @@ bool Map3d::get_pdok_output(std::string filename) {
   int i = 1;
   for (auto& f : _lsFeatures) {
     std::string layername = f->get_layername();
-    //Add additional attribute describing CityGML of feature
-    std::wstring_convert<codecvt<wchar_t, char, std::mbstate_t>>converter;
-    std::wstringstream ss;
-    ss << std::fixed << std::setprecision(3);
-    f->get_citygml_imgeo(ss);
-    std::string gmlAttribute = converter.to_bytes(ss.str());
-    ss.clear();
     AttributeMap extraAttribute = AttributeMap();
-    extraAttribute["xml"] = std::make_pair(OFTString, gmlAttribute);
-
-    if (!f->get_shape(layers[layername], true, extraAttribute)) {
-      return false;
-    }
-    if (i % 1000 == 0) {
-      if (dataSource->CommitTransaction() != OGRERR_NONE) {
-        std::cerr << "Writing to database failed.\n";
-        return false;
-      }  
-      if (dataSource->StartTransaction() != OGRERR_NONE) {
-        std::cerr << "Starting database transaction failed.\n";
-        return false;
+    
+    if (pdok) {
+      //Add additional attribute describing CityGML of feature
+      std::wstring_convert<codecvt<wchar_t, char, std::mbstate_t>>converter;
+      std::wstringstream ss;
+      ss << std::fixed << std::setprecision(3);
+      if (citygml) {
+        f->get_citygml(ss);
       }
-    }
-    i++;
-  }
-  if (dataSource->CommitTransaction() != OGRERR_NONE) {
-    std::cerr << "Writing to database failed.\n";
-    return false;
-  }
-  GDALClose(dataSource);
-  GDALClose(driver);
-  return true;
-#endif
-}
-
-bool Map3d::get_pdok_citygml_output(std::string filename) {
-#if GDAL_VERSION_MAJOR < 2
-  std::cerr << "ERROR: cannot write MultiPolygonZ files with GDAL < 2.0.\n";
-  return false;
-#else
-  if (GDALGetDriverCount() == 0)
-    GDALAllRegister();
-  GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("PostgreSQL");
-  GDALDataset* dataSource = driver->Create(filename.c_str(), 0, 0, 0, GDT_Unknown, NULL);
-  if (dataSource == NULL) {
-    std::cerr << "Starting database connection failed.\n";
-    return false;
-  }
-  if (dataSource->StartTransaction() != OGRERR_NONE) {
-    std::cerr << "Starting database transaction failed.\n";
-    return false;
-  }
-
-  std::unordered_map<std::string, OGRLayer*> layers;
-  // create and write layers first
-  for (auto& f : _lsFeatures) {
-    std::string layername = f->get_layername();
-    if (layers.find(layername) == layers.end()) {
-      AttributeMap &attributes = f->get_attributes();
-      //Add additional attribute to list for layer creation
-      attributes["xml"] = std::make_pair(OFTString, "");
-      OGRLayer *layer = create_gdal_layer(driver, dataSource, filename, layername, attributes, f->get_class() == BUILDING);
-      if (layer == NULL) {
-        std::cerr << "ERROR: Cannot open database '" + filename + "' for writing" << std::endl;
-        dataSource->RollbackTransaction();
-        GDALClose(dataSource);
-        GDALClose(driver);
-        return false;
+      else {
+        f->get_citygml_imgeo(ss);
       }
-      layers.emplace(layername, layer);
+      std::string gmlAttribute = converter.to_bytes(ss.str());
+      ss.clear();
+      extraAttribute["xml"] = std::make_pair(OFTString, gmlAttribute);
     }
-  }
-  if (dataSource->CommitTransaction() != OGRERR_NONE) {
-    std::cerr << "Writing to database failed.\n";
-    return false;
-  }
-  if (dataSource->StartTransaction() != OGRERR_NONE) {
-    std::cerr << "Starting database transaction failed.\n";
-    return false;
-  }
-
-  // create and write features to layers
-  int i = 1;
-  for (auto& f : _lsFeatures) {
-    std::string layername = f->get_layername();
-    //Add additional attribute describing CityGML of feature
-    std::wstring_convert<codecvt<wchar_t, char, std::mbstate_t>>converter;
-    std::wstringstream ss;
-    ss << std::fixed << std::setprecision(3);
-    f->get_citygml(ss);
-    std::string gmlAttribute = converter.to_bytes(ss.str());
-    ss.clear();
-    AttributeMap extraAttribute = AttributeMap();
-    extraAttribute["xml"] = std::make_pair(OFTString, gmlAttribute);
-
     if (!f->get_shape(layers[layername], true, extraAttribute)) {
       return false;
     }
@@ -553,7 +487,7 @@ bool Map3d::get_pdok_citygml_output(std::string filename) {
       if (dataSource->CommitTransaction() != OGRERR_NONE) {
         std::cerr << "Writing to database failed.\n";
         return false;
-      }  
+      }
       if (dataSource->StartTransaction() != OGRERR_NONE) {
         std::cerr << "Starting database transaction failed.\n";
         return false;
@@ -611,7 +545,9 @@ bool Map3d::get_gdal_output(std::string filename, std::string drivername, bool m
         }
         layers.emplace(layername, layer);
       }
-      f->get_shape(layers[layername], true);
+      if (!f->get_shape(layers[layername], true)) {
+        return false;
+      }
     }
     close_gdal_resources(driver, layers);
   }
@@ -736,27 +672,30 @@ const std::vector<TopoFeature*>& Map3d::get_polygons3d() {
   return _lsFeatures;
 }
 
-void Map3d::add_elevation_point(liblas::Point const& laspt) {
+void Map3d::add_elevation_point(LASpoint const& laspt) {
+  //-- only process last returns; 
+  //-- although perhaps not smart for vegetation/forest in the future
+  //-- TODO: always ignore the non-last-return points?
+  if (laspt.return_number != laspt.number_of_returns)
+    return;
+
   std::vector<PairIndexed> re;
-  Point2 minp(laspt.GetX() - _radius_vertex_elevation, laspt.GetY() - _radius_vertex_elevation);
-  Point2 maxp(laspt.GetX() + _radius_vertex_elevation, laspt.GetY() + _radius_vertex_elevation);
+  float x = laspt.get_x();
+  float y = laspt.get_y();
+  Point2 minp(x - _radius_vertex_elevation, y - _radius_vertex_elevation);
+  Point2 maxp(x + _radius_vertex_elevation, y + _radius_vertex_elevation);
   Box2 querybox(minp, maxp);
   _rtree.query(bgi::intersects(querybox), std::back_inserter(re));
-  minp = Point2(laspt.GetX() - _building_radius_vertex_elevation, laspt.GetY() - _building_radius_vertex_elevation);
-  maxp = Point2(laspt.GetX() + _building_radius_vertex_elevation, laspt.GetY() + _building_radius_vertex_elevation);
+  minp = Point2(x - _building_radius_vertex_elevation, y - _building_radius_vertex_elevation);
+  maxp = Point2(x + _building_radius_vertex_elevation, y + _building_radius_vertex_elevation);
   querybox = Box2(minp, maxp);
   _rtree_buildings.query(bgi::intersects(querybox), std::back_inserter(re));
 
   for (auto& v : re) {
     TopoFeature* f = v.second;
     float radius = _radius_vertex_elevation;
-    //-- only process last returns; 
-    //-- although perhaps not smart for vegetation/forest in the future
-    //-- TODO: always ignore the non-last-return points?
-    if (laspt.GetReturnNumber() != laspt.GetNumberOfReturns()) 
-      continue;
 
-    int c = laspt.GetClassification().GetClass();
+    int c = (int)laspt.classification;
     bool bInsert = false;
     bool bWithin = false;
     if (f->get_class() == BUILDING) {
@@ -817,10 +756,9 @@ void Map3d::add_elevation_point(liblas::Point const& laspt) {
         bWithin = true;
       }
     }
- 
     if (bInsert == true) { //-- only insert if in the allowed LAS classes
-      Point2 p(laspt.GetX(), laspt.GetY());
-      f->add_elevation_point(p, laspt.GetZ(), radius, c, bWithin); 
+      Point2 p(x, y);
+      f->add_elevation_point(p, laspt.get_z(), radius, c, bWithin);
     }
   }
 }
@@ -903,7 +841,7 @@ bool Map3d::construct_CDT() {
       p->buildCDT();
     }
     catch (std::exception e) {
-      std::cerr << std::endl << "CDT failed for " << p->get_id() << " (" << p->get_class() << ") with error: " << e.what() << std::endl;
+      std::cerr << std::endl << "CDT failed for object \'" << p->get_id() << "\' (class " << p->get_class() << ") with error: " << e.what() << std::endl;
       return false;
     }
   }
@@ -935,6 +873,12 @@ bool Map3d::construct_rtree() {
       std::min(bg::get<bg::min_corner, 1>(_rtree.bounds()), bg::get<bg::min_corner, 1>(_rtree_buildings.bounds()))),
     Point2(std::max(bg::get<bg::max_corner, 0>(_rtree.bounds()), bg::get<bg::max_corner, 0>(_rtree_buildings.bounds())),
       std::max(bg::get<bg::max_corner, 1>(_rtree.bounds()), bg::get<bg::max_corner, 1>(_rtree_buildings.bounds()))));
+  
+  double radius = std::max(_radius_vertex_elevation, _building_radius_vertex_elevation);
+  _minxradius = std::min(bg::get<bg::min_corner, 0>(_rtree.bounds()), bg::get<bg::min_corner, 0>(_rtree_buildings.bounds())) - radius;
+  _maxxradius = std::min(bg::get<bg::min_corner, 1>(_rtree.bounds()), bg::get<bg::min_corner, 1>(_rtree_buildings.bounds())) + radius;
+  _minyradius = std::max(bg::get<bg::max_corner, 0>(_rtree.bounds()), bg::get<bg::max_corner, 0>(_rtree_buildings.bounds())) - radius;
+  _maxyradius = std::max(bg::get<bg::max_corner, 1>(_rtree.bounds()), bg::get<bg::max_corner, 1>(_rtree_buildings.bounds())) + radius;
   return true;
 }
 
@@ -1135,57 +1079,57 @@ void Map3d::extract_feature(OGRFeature *f, std::string layername, const char *id
   }
 }
 
-//-- http://www.liblas.org/tutorial/cpp.html#applying-filters-to-a-reader-to-extract-specified-classes
 bool Map3d::add_las_file(PointFile pointFile) {
   std::clog << "Reading LAS/LAZ file: " << pointFile.filename << std::endl;
-  std::ifstream ifs;
-  ifs.open(pointFile.filename.c_str(), std::ios::in | std::ios::binary);
-  if (ifs.is_open() == false) {
-    std::cerr << "\tERROR: could not open file: " << pointFile.filename << std::endl;
-    return false;
-  }
-  //-- LAS classes to omit
-  std::vector<liblas::Classification> liblasomits;
-  for (int i : pointFile.lasomits) {
-    liblasomits.push_back(liblas::Classification(i));
-  }
 
-  //-- read each point 1-by-1
-  liblas::ReaderFactory f;
-  liblas::Reader reader = f.CreateWithStream(ifs);
-  liblas::Header const& header = reader.GetHeader();
-
-  //-- check if the file overlaps the polygons
-  liblas::Bounds<double> bounds = header.GetExtent();
-  liblas::Bounds<double> polygonBounds = get_bounds();
-  uint32_t pointCount = header.GetPointRecordsCount();
-  if (polygonBounds.intersects(bounds)) {
-    std::clog << "\t(" << boost::locale::as::number << pointCount << " points in the file)\n";
-    if ((pointFile.thinning > 1)) {
-      std::clog << "\t(skipping every " << pointFile.thinning << "th points, thus ";
-      std::clog << boost::locale::as::number << (pointCount / pointFile.thinning) << " are used)\n";
-    }
-    else
-      std::clog << "\t(all points used, no skipping)\n";
-
-    if (pointFile.lasomits.empty() == false) {
-      std::clog << "\t(omitting LAS classes: ";
-      for (int i : pointFile.lasomits)
-        std::clog << i << " ";
-      std::clog << ")\n";
+  LASreadOpener lasreadopener;
+  lasreadopener.set_file_name(pointFile.filename.c_str());
+  //-- set to compute bounding box
+  lasreadopener.set_populate_header(true);
+  LASreader* lasreader = lasreadopener.open();
+
+  try {
+    //-- check if file is open
+    if (lasreader == 0) {
+      std::cerr << "\tERROR: could not open file: " << pointFile.filename << std::endl;
+      return false;
     }
-    printProgressBar(0);
-    int i = 0;
-    
-    try {
-      while (reader.ReadNextPoint()) {
-        liblas::Point const& p = reader.GetPoint();
+    LASheader header = lasreader->header;
+
+    if (check_bounds(header.min_x, header.max_x, header.min_y, header.max_y)) {
+      //-- LAS classes to omit
+      std::vector<int> lasomits;
+      for (int i : pointFile.lasomits) {
+        lasomits.push_back(i);
+      }
+
+      //-- read each point 1-by-1
+      uint32_t pointCount = header.number_of_point_records;
+
+      std::clog << "\t(" << boost::locale::as::number << pointCount << " points in the file)\n";
+      if ((pointFile.thinning > 1)) {
+        std::clog << "\t(skipping every " << pointFile.thinning << "th points, thus ";
+        std::clog << boost::locale::as::number << (pointCount / pointFile.thinning) << " are used)\n";
+      }
+      else
+        std::clog << "\t(all points used, no skipping)\n";
+
+      if (pointFile.lasomits.empty() == false) {
+        std::clog << "\t(omitting LAS classes: ";
+        for (int i : pointFile.lasomits)
+          std::clog << i << " ";
+        std::clog << ")\n";
+      }
+      printProgressBar(0);
+      int i = 0;
+      while (lasreader->read_point()) {
+        LASpoint const& p = lasreader->point;
         //-- set the thinning filter
         if (i % pointFile.thinning == 0) {
           //-- set the classification filter
-          if (std::find(liblasomits.begin(), liblasomits.end(), p.GetClassification()) == liblasomits.end()) {
+          if (std::find(lasomits.begin(), lasomits.end(), (int)p.classification) == lasomits.end()) {
             //-- set the bounds filter
-            if (polygonBounds.contains(p)) {
+            if (check_bounds(p.X, p.X, p.Y, p.Y)) {
               this->add_elevation_point(p);
             }
           }
@@ -1197,16 +1141,16 @@ bool Map3d::add_las_file(PointFile pointFile) {
       printProgressBar(100);
       std::clog << std::endl;
     }
-    catch (std::exception e) {
-      std::cerr << std::endl << e.what() << std::endl;
-      ifs.close();
-      return false;
+    else {
+      std::clog << "\tskipping file, bounds do not intersect polygon extent\n";
     }
+    lasreader->close();
   }
-  else {
-    std::clog << "\tskipping file, bounds do not intersect polygon extent\n";
+  catch (std::exception e) {
+    std::cerr << std::endl << e.what() << std::endl;
+    lasreader->close();
+    return false;
   }
-  ifs.close();
   return true;
 }
 
@@ -1594,7 +1538,7 @@ void Map3d::stitch_bridges() {
             pis.clear();
             if (!(fadj->get_class() == BRIDGE && fadj->get_top_level()) && fadj->has_point2(ring[i], ringis, pis)) {
               int z = fadj->get_vertex_elevation(ringis[0], pis[0]);
-              if (abs(f->get_vertex_elevation(ringi, i) - z) < _threshold_jump_edges) {
+              if (abs(f->get_vertex_elevation(ringi, i) - z) < _threshold_bridge_jump_edges) {
                 f->set_vertex_elevation(ringi, i, z);
                 if (!(fadj->get_class() == BRIDGE && fadj->get_top_level() == f->get_top_level())) {
                   // Add height to NC
@@ -1766,7 +1710,7 @@ void Map3d::stitch_bridges() {
               int interz = interpolate_height(f, p, ringi, previ, ringi, endCorner.first);
               int prevz = f->get_vertex_elevation(ringi, previ);
               //Allways stich to lower object or if interpolated between corners within threshold or previous within threshold
-              if (stitchz < interz || abs(stitchz - interz) < _threshold_jump_edges || abs(stitchz - prevz) < _threshold_jump_edges) {
+              if (stitchz < interz || abs(stitchz - interz) < _threshold_bridge_jump_edges || abs(stitchz - prevz) < _threshold_bridge_jump_edges) {
                 f->set_vertex_elevation(ringi, pi, stitchz);
               }
               else {
@@ -1798,4 +1742,3 @@ void Map3d::add_allowed_las_class(AllowedLASTopo c, int i) {
 void Map3d::add_allowed_las_class_within(AllowedLASTopo c, int i) {
   _las_classes_allowed_within[c].insert(i);
 }
-
diff --git a/src/Map3d.h b/src/Map3d.h
index 931cb017..1c720596 100644
--- a/src/Map3d.h
+++ b/src/Map3d.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
@@ -56,13 +56,13 @@ class Map3d {
   bool construct_rtree();
   bool threeDfy(bool stitching = true);
   bool construct_CDT();
-  void add_elevation_point(liblas::Point const& laspt);
+  void add_elevation_point(LASpoint const& laspt);
   void cleanup_elevations();
 
   unsigned long get_num_polygons();
   const std::vector<TopoFeature*>&  get_polygons3d();
   Box2 get_bbox();
-  liblas::Bounds<double> get_bounds();
+  bool check_bounds(const double xmin, const double xmax, const double ymin, const double ymax);
 
   void get_citygml(std::wostream& of);
   void get_citygml_multifile(std::string);
@@ -71,8 +71,7 @@ class Map3d {
   bool get_cityjson(std::string filename);
   void get_citygml_imgeo_multifile(std::string ofname);
   void create_citygml_imgeo_header(std::wostream& of);
-  bool get_pdok_output(std::string filename);
-  bool get_pdok_citygml_output(std::string filename);
+  bool get_postgis_output(std::string filename, bool pdok = false, bool citygml = false);
   bool get_gdal_output(std::string filename, std::string drivername, bool multi);
   void get_csv_buildings(std::wostream& of);
   void get_csv_buildings_multiple_heights(std::wostream& of);
@@ -103,6 +102,7 @@ class Map3d {
   void set_radius_vertex_elevation(float radius);
   void set_building_radius_vertex_elevation(float radius);
   void set_threshold_jump_edges(float threshold);
+  void set_threshold_bridge_jump_edges(float threshold);
   void set_requested_extent(double xmin, double ymin, double xmax, double ymax);
 
   void add_allowed_las_class(AllowedLASTopo c, int i);
@@ -133,7 +133,12 @@ class Map3d {
   float       _radius_vertex_elevation;
   float       _building_radius_vertex_elevation;
   int         _threshold_jump_edges; //-- in cm/integer
+  int         _threshold_bridge_jump_edges; //-- in cm/integer
   Box2        _bbox;
+  double      _minxradius;
+  double      _maxxradius;
+  double      _minyradius;
+  double      _maxyradius;
   Box2        _requestedExtent;
 
   //-- storing the LAS allowed for each TopoFeature
@@ -144,7 +149,6 @@ class Map3d {
   NodeColumn                                          _nc_building_walls;
   std::unordered_map<std::string, int>                _bridge_stitches;
   std::vector<TopoFeature*>                           _lsFeatures;
-  std::vector<std::string>                            _allowed_layers;
   bgi::rtree< PairIndexed, bgi::rstar<16> >           _rtree;
   bgi::rtree< PairIndexed, bgi::rstar<16> >           _rtree_buildings;
 
diff --git a/src/Road.cpp b/src/Road.cpp
index 1dcbfcbc..6a5e6c38 100644
--- a/src/Road.cpp
+++ b/src/Road.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Road.h b/src/Road.h
index 75ee95de..ec23520d 100644
--- a/src/Road.h
+++ b/src/Road.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Separation.cpp b/src/Separation.cpp
index 2f11428e..c2260b14 100644
--- a/src/Separation.cpp
+++ b/src/Separation.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Separation.h b/src/Separation.h
index a7f0b045..04e6e59f 100644
--- a/src/Separation.h
+++ b/src/Separation.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Terrain.cpp b/src/Terrain.cpp
index 264c3859..16dfb65c 100644
--- a/src/Terrain.cpp
+++ b/src/Terrain.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/Terrain.h b/src/Terrain.h
index c04d6da6..f8483d9d 100644
--- a/src/Terrain.h
+++ b/src/Terrain.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/TopoFeature.cpp b/src/TopoFeature.cpp
index 85db134d..71e51a97 100644
--- a/src/TopoFeature.cpp
+++ b/src/TopoFeature.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
@@ -373,8 +373,13 @@ bool TopoFeature::get_multipolygon_features(OGRLayer* layer, std::string classNa
 bool TopoFeature::writeAttribute(OGRFeature* feature, OGRFeatureDefn* featureDefn, std::string name, std::string value) {
   int fi = featureDefn->GetFieldIndex(name.c_str());
   if (fi == -1) {
-    std::cerr << "Failed to write attribute " << name << ".\n";
-    return false;
+    // try replace '-' with '_' for postgresql column names
+    std::replace(name.begin(), name.end(), '-', '_');
+    fi = featureDefn->GetFieldIndex(name.c_str());
+    if (fi == -1) {
+      std::cerr << "Failed to write attribute " << name << ".\n";
+      return false;
+    }
   }
   // perform extra character encoding for gdal.
   char* attrcpl = CPLRecode(value.c_str(), "", CPL_ENC_UTF8);
@@ -1180,6 +1185,8 @@ void Boundary3D::detect_outliers(bool flatten){
       // find spikes in roads (due to misclassified lidar points) and fix by averaging between previous and next vertex.
       std::vector<int> idx;
       std::vector<double> x, y, z, coeffs;
+      double x0 = ring[0].x();
+      double y0 = ring[0].y();
       for (int i = 0; i < ring.size(); i++) {
         idx.push_back(i);
         x.push_back(ring[i].x());
@@ -1195,7 +1202,7 @@ void Boundary3D::detect_outliers(bool flatten){
       for (int i = 0; i < niter; i++) {
         // Fit the model
         std::vector<double> correctedvalues;
-        polyfit3d<double>(xtmp, ytmp, ztmp, coeffs, correctedvalues);
+        polyfit3d<double>(xtmp, ytmp, ztmp, x0, y0, coeffs, correctedvalues);
         std::vector<double> residuals, absResiduals;
 
         double sum = 0;
@@ -1244,7 +1251,7 @@ void Boundary3D::detect_outliers(bool flatten){
 
       // get the new values based on the coeffs of the lase equation
       std::vector<double> correctedvalues;
-      polyval3d<double>(x, y, coeffs, correctedvalues);
+      polyval3d<double>(x, y, x0, y0, coeffs, correctedvalues);
       if (flatten) {
         for (int i = 0; i < ring.size(); i++) {
           _p2z[ringi][i] = correctedvalues[i];
diff --git a/src/TopoFeature.h b/src/TopoFeature.h
index 756eed7a..73a3aed9 100644
--- a/src/TopoFeature.h
+++ b/src/TopoFeature.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
@@ -34,7 +34,6 @@
 #include "io.h"
 #include "polyfit.hpp"
 #include "nlohmann-json/json.hpp"
-#include "ptinpoly.h"
 
 class TopoFeature {
 public:
diff --git a/src/Water.cpp b/src/Water.cpp
index f5941610..c0e86070 100644
--- a/src/Water.cpp
+++ b/src/Water.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
@@ -108,10 +108,10 @@ void Water::get_citygml_imgeo(std::wostream& of) {
   std::string attribute;
   if (ondersteunend) {
     if (get_attribute("bgt-type", attribute)) {
-      of << "<wat:class codeSpace=\"http://www.geostandaarden.nl/imgeo/def/2.1#TypeOndersteunendWaterdeel\">" << attribute << "</wat:class>";
+      of << "<wtr:class codeSpace=\"http://www.geostandaarden.nl/imgeo/def/2.1#TypeOndersteunendWaterdeel\">" << attribute << "</wat:class>";
     }
     else if (get_attribute("bgt_type", attribute)) {
-      of << "<wat:class codeSpace=\"http://www.geostandaarden.nl/imgeo/def/2.1#TypeOndersteunendWaterdeel\">" << attribute << "</wat:class>";
+      of << "<wtr:class codeSpace=\"http://www.geostandaarden.nl/imgeo/def/2.1#TypeOndersteunendWaterdeel\">" << attribute << "</wat:class>";
     }
     if (get_attribute("plus-type", attribute)) {
       of << "<imgeo:plus-type codeSpace=\"http://www.geostandaarden.nl/imgeo/def/2.1#TypeOndersteunendWaterdeelPlus\">" << attribute << "</imgeo:plus-type>";
@@ -123,10 +123,10 @@ void Water::get_citygml_imgeo(std::wostream& of) {
   }
   else {
     if (get_attribute("bgt-type", attribute)) {
-      of << "<wat:class codeSpace=\"http://www.geostandaarden.nl/imgeo/def/2.1#TypeWater\">" << attribute << "</wat:class>";
+      of << "<wtr:class codeSpace=\"http://www.geostandaarden.nl/imgeo/def/2.1#TypeWater\">" << attribute << "</wat:class>";
     }
     else if (get_attribute("bgt_type", attribute)) {
-      of << "<wat:class codeSpace=\"http://www.geostandaarden.nl/imgeo/def/2.1#TypeWater\">" << attribute << "</wat:class>";
+      of << "<wtr:class codeSpace=\"http://www.geostandaarden.nl/imgeo/def/2.1#TypeWater\">" << attribute << "</wat:class>";
     }
     if (get_attribute("plus-type", attribute)) {
       of << "<imgeo:plus-type codeSpace=\"http://www.geostandaarden.nl/imgeo/def/2.1#TypeWaterPlus\">" << attribute << "</imgeo:plus-type>";
diff --git a/src/Water.h b/src/Water.h
index 0edeb611..820c23b7 100644
--- a/src/Water.h
+++ b/src/Water.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/definitions.h b/src/definitions.h
index bd0cf00d..866e1cf7 100644
--- a/src/definitions.h
+++ b/src/definitions.h
@@ -2,11 +2,12 @@
 #define __3DFIER__Definitions__
 
 #include <iostream>
+#include <fstream>
 #include <cmath>
 #include <vector>
 #include <unordered_map>
 
-#include <liblas/liblas.hpp>
+#include <lasreader.hpp>
 #include <ogrsf_frmts.h>
 
 #include <boost/geometry.hpp>
diff --git a/src/geomtools.cpp b/src/geomtools.cpp
index aee5fd99..48184094 100644
--- a/src/geomtools.cpp
+++ b/src/geomtools.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
@@ -39,11 +39,14 @@
 #include <unordered_set>
 #include <boost/heap/fibonacci_heap.hpp>
 
+typedef CGAL::Exact_predicates_inexact_constructions_kernel			K;
+
 // fibonacci heap for greedy insertion code
 struct point_error {
-  point_error(int i, double e) : index(i), error(e){}
+  point_error(int i, double e, CGAL::Point_3<K> p) : index(i), error(e), point(p){}
   int index;
   double error;
+  CGAL::Point_3<K> point;
   
   bool operator<(point_error const & rhs) const
   {
@@ -52,8 +55,8 @@ struct point_error {
 };
 typedef boost::heap::fibonacci_heap<point_error> Heap;
 typedef Heap::handle_type heap_handle;
+typedef std::vector<heap_handle> heap_handle_vec;
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel			K;
 typedef CGAL::Projection_traits_xy_3<K>								Gt;
 typedef CGAL::Triangulation_vertex_base_with_id_2<Gt>				Vb;
 struct FaceInfo2
@@ -63,16 +66,16 @@ struct FaceInfo2
   bool in_domain() {
     return nesting_level % 2 == 1;
   }
-  std::vector<heap_handle>* points_inside = nullptr;
+  heap_handle_vec* points_inside = nullptr;
   CGAL::Plane_3<K>* plane = nullptr;
 };
-typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2, Gt>	Fbb;
-typedef CGAL::Constrained_triangulation_face_base_2<Gt, Fbb>		Fb;
-typedef CGAL::Triangulation_data_structure_2<Vb, Fb>				Tds;
-typedef CGAL::Exact_predicates_tag									Itag;
-typedef CGAL::Constrained_Delaunay_triangulation_2<Gt, Tds, Itag>	CDT;
-typedef CDT::Point													Point;
-typedef CGAL::Polygon_2<Gt>											Polygon_2;
+typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2, Gt>  Fbb;
+typedef CGAL::Constrained_triangulation_face_base_2<Gt, Fbb>      Fb;
+typedef CGAL::Triangulation_data_structure_2<Vb, Fb>              Tds;
+typedef CGAL::Exact_predicates_tag                                Itag;
+typedef CGAL::Constrained_Delaunay_triangulation_2<Gt, Tds, Itag> CDT;
+typedef CDT::Point                                                Point;
+typedef CGAL::Polygon_2<Gt>                                       Polygon_2;
 
 struct PointXYHash {
   std::size_t operator()(Point const& p) const noexcept {
@@ -245,7 +248,7 @@ inline double compute_error(Point &p, CDT::Face_handle &face) {
       face->vertex(1)->point(),
       face->vertex(2)->point());
   if(!face->info().points_inside)
-    face->info().points_inside = new std::vector<heap_handle>();
+    face->info().points_inside = new heap_handle_vec();
 
   auto plane = face->info().plane;
   auto interpolate = - plane->a()/plane->c() * p.x() - plane->b()/plane->c()*p.y() - plane->d()/plane->c();
@@ -257,71 +260,133 @@ void greedy_insert(CDT &T, const std::vector<Point3> &pts, double threshold) {
   // assumes all lidar points are inside a triangle
   Heap heap;
 
-  // Convert all elevation points to CGAL points
-  std::vector<Point> cpts;
-  cpts.reserve(pts.size());
-  for (auto& p : pts) {
-    cpts.push_back(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p)));
-  }
-
   // compute initial point errors, build heap, store point indices in triangles
   {
-    std::unordered_set<Point, PointXYHash, PointXYEqual> set;
-    for(int i=0; i<cpts.size(); i++){
-      auto p = cpts[i];
-      // detect and skip duplicate points
-      auto not_duplicate = set.insert(p).second;
-      if(not_duplicate){
-        auto face = T.locate(p);
-        auto e = compute_error(p, face);
-        auto handle = heap.push(point_error(i,e));
+    for (int i = 0; i < pts.size(); i++) {
+      Point p = Point(bg::get<0>(pts[i]), bg::get<1>(pts[i]), bg::get<2>(pts[i]));
+      CDT::Locate_type lt;
+      int li;
+      CDT::Face_handle face = T.locate(p, lt, li);
+      if (lt == CDT::FACE) {
+        double e = compute_error(p, face);
+        auto handle = heap.push(point_error(i, e, p));
         face->info().points_inside->push_back(handle);
       }
+      else {
+        std::cout << "CDT insert; point location not in face but ";
+        if (lt == CDT::VERTEX) {
+          std::cout << "on vertex.";
+        }
+        else if (lt == CDT::EDGE) {
+          std::cout << "on edge.";
+        }
+        else if (lt == CDT::OUTSIDE_CONVEX_HULL) {
+          std::cout << "outside convex hull.";
+        }
+        else if (lt == CDT::OUTSIDE_AFFINE_HULL) {
+          std::cout << "outside affine hull.";
+        }
+        std::cout << " Point; " << std::fixed << std::setprecision(3) << p << std::endl;
+      }
     }
   }
   
   // insert points, update errors of affected triangles until threshold error is reached
   while (!heap.empty() && heap.top().error > threshold){
     // get top element (with largest error) from heap
-    auto maxelement = heap.top();
-    auto max_p = cpts[maxelement.index];
+    point_error maxelement = heap.top();
+    auto max_p = maxelement.point;
 
     // get triangles that will change after inserting this max_p
     std::vector<CDT::Face_handle> faces;
-    T.get_conflicts ( max_p, std::back_inserter(faces) );
+    T.get_conflicts(max_p, std::back_inserter(faces));
+
+    // handle case where max_p somehow coincides with a polygon vertex
+    if (faces.size()==0) {
+        // this should return the already existing vertex
+        auto v = T.insert(max_p);
+
+        // check the incident faces and erase any references to max_p
+        CDT::Face_circulator fcirculator = T.incident_faces(v), done(fcirculator);
+        auto start = fcirculator;
+        do {
+          auto face = *fcirculator;
+          if (face.info().points_inside) {
+            // collect the heap_handles that need to be removed
+            std::vector<heap_handle_vec::iterator> to_erase;
+            for (auto it = face.info().points_inside->begin(); it != face.info().points_inside->end(); ++it) {
+              if (maxelement.index == (**it).index) {
+                to_erase.push_back(it);
+              }
+            }
+            // remove the collected heap_handles
+            for (auto it : to_erase) {
+              face.info().points_inside->erase(it);
+            }
+          }
+        } while (++fcirculator != done);
+        
+        // remove this points from the heap
+        heap.pop();
+        
+        // skip to next iteration
+        continue;
+    }
+
+    // clear info of triangles that just changed, collect points that were inside these triangles
+    heap_handle_vec points_to_update;
 
-    // insert max_p in triangulation
-    auto face_hint = faces[0];
-    auto v = T.insert(max_p, face_hint);
-    face_hint = v->face();
-    
-    // update clear info of triangles that just changed, collect points that were inside these triangles
-    std::vector<heap_handle> points_to_update;
     for (auto face : faces) {
-      if (face->info().plane){
+      if (face->info().plane) {
         delete face->info().plane;
         face->info().plane = nullptr;
       }
-      if (face->info().points_inside) {
-        for (auto h :*face->info().points_inside){
-          if( maxelement.index != (*h).index)
+      if (face->info().points_inside && face->info().points_inside->size() > 0) {
+        for (auto h : *face->info().points_inside) {
+          if (maxelement.index != (*h).index) {
             points_to_update.push_back(h);
+          }
         }
-        face->info().points_inside->clear();
+        heap_handle_vec().swap((*face->info().points_inside));
       }
     }
-    
-    // remove the point we just inserted in the triangulation from the heap
+
+    // insert max_p in triangulation
+    auto face_hint = faces[0];
+    auto v = T.insert(max_p, face_hint);
+    face_hint = v->face();
+
+    // remove this points from the heap
     heap.pop();
 
     // update the errors of affected elevation points
-    for (auto curelement : points_to_update){
-      auto p = cpts[(*curelement).index];
-      auto containing_face = T.locate(p, face_hint);
-      const double e = compute_error(p, containing_face);
-      const point_error new_pe = point_error((*curelement).index, e);
-      heap.update(curelement, new_pe);
-      containing_face->info().points_inside->push_back(curelement);
+    for (heap_handle curelement : points_to_update){
+      point_error element = *curelement;
+      auto p = element.point;
+      CDT::Locate_type lt;
+      int li;
+      CDT::Face_handle containing_face = T.locate(p, lt, li, face_hint);
+      if (lt == CDT::EDGE || lt == CDT::FACE) {
+        element.error = compute_error(p, containing_face);
+        heap.update(curelement, element);
+        containing_face->info().points_inside->push_back(curelement);
+      }
+      else {
+        std::cout << "CDT update; point location not in face but ";
+        if (lt == CDT::VERTEX) {
+          std::cout << "on vertex.";
+          // update error to 0.0 to disable adding point to CDT
+          element.error = 0.0;
+          heap.update(curelement, element);
+        }
+        else if (lt == CDT::OUTSIDE_CONVEX_HULL) {
+          std::cout << "outside convex hull.";
+        }
+        else if (lt == CDT::OUTSIDE_AFFINE_HULL) {
+          std::cout << "outside affine hull.";
+        }
+        std::cout << " Point (index " << element.index << "); " << std::fixed << std::setprecision(3) << p << std::endl;
+      }
     }
   }
 
diff --git a/src/geomtools.h b/src/geomtools.h
index 7453ef24..a15fdc9f 100644
--- a/src/geomtools.h
+++ b/src/geomtools.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
   
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/io.cpp b/src/io.cpp
index 28bccd41..9ae8580f 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
@@ -67,7 +67,7 @@ void get_citygml_namespaces(std::wostream& of) {
   of << "xmlns:app=\"http://www.opengis.net/citygml/appearance/2.0\"\n";
   of << "xmlns:tun=\"http://www.opengis.net/citygml/tunnel/2.0\"\n";
   of << "xmlns:cif=\"http://www.opengis.net/citygml/cityfurniture/2.0\"\n";
-  of << "xsi:schemaLocation=\"http://www.opengis.net/citygml/2.0 ./CityGML_2.0/CityGML.xsd\">\n";
+  of << "xsi:schemaLocation=\"http://www.opengis.net/citygml/2.0 http://schemas.opengis.net/citygml/profiles/base/2.0/CityGML.xsd\">\n";
 }
 
 void get_citygml_imgeo_namespaces(std::wostream& of) {
@@ -136,7 +136,7 @@ void get_extruded_line_gml(std::wostream& of, Point2* a, Point2* b, double high,
   of << bg::get<0>(a) << " " << bg::get<1>(a) << " " << std::setprecision(2) << low << std::setprecision(3) << " ";
   of << bg::get<0>(a) << " " << bg::get<1>(a) << " " << std::setprecision(2) << high << std::setprecision(3) << " ";
   of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << high << std::setprecision(3) << " ";
-  of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << low << std::setprecision(3);
+  of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << low << std::setprecision(3); 
   of << "</gml:posList>";
   of << "</gml:LinearRing>";
   of << "</gml:exterior>";
diff --git a/src/io.h b/src/io.h
index b06bb92d..a66e54d7 100644
--- a/src/io.h
+++ b/src/io.h
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
diff --git a/src/main.cpp b/src/main.cpp
index eb518d73..7ebffe5e 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,7 +1,7 @@
 /*
   3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models.
 
-  Copyright (C) 2015-2018  3D geoinformation research group, TU Delft
+  Copyright (C) 2015-2019  3D geoinformation research group, TU Delft
 
   This file is part of 3dfier.
 
@@ -36,7 +36,7 @@
 #include <boost/program_options.hpp>
 #include <boost/filesystem.hpp>
 
-std::string VERSION = "1.1";
+std::string VERSION = "1.2.2";
 
 bool validate_yaml(const char* arg, std::set<std::string>& allowedFeatures);
 int main(int argc, const char * argv[]);
@@ -52,7 +52,7 @@ int main(int argc, const char * argv[]) {
   std::cout.imbue(loc);
 
   std::string licensewarning =
-    "3dfier Copyright (C) 2015-2018  3D geoinformation research group, TU Delft\n"
+    "3dfier Copyright (C) 2015-2019  3D geoinformation research group, TU Delft\n"
     "This program comes with ABSOLUTELY NO WARRANTY.\n"
     "This is free software, and you are welcome to redistribute it\n"
     "under certain conditions; for details run 3dfier with the '--license' option.\n";
@@ -71,7 +71,6 @@ int main(int argc, const char * argv[]) {
   outputs["Shapefile"] = "";
   outputs["Shapefile-Multifile"] = "";
   outputs["PostGIS"] = "";
-  outputs["PostGIS-Multi"] = "";
   outputs["PostGIS-PDOK"] = "";
   outputs["PostGIS-PDOK-CityGML"] = "";
   outputs["GDAL"] = "";
@@ -96,7 +95,6 @@ int main(int argc, const char * argv[]) {
       ("Shapefile", po::value<std::string>(&outputs["Shapefile"]), "Output ")
       ("Shapefile-Multifile", po::value<std::string>(&outputs["Shapefile-Multifile"]), "Output ")
       ("PostGIS", po::value<std::string>(&outputs["PostGIS"]), "Output ")
-      ("PostGIS-Multi", po::value<std::string>(&outputs["PostGIS-Multi"]), "Output ")
       ("PostGIS-PDOK", po::value<std::string>(&outputs["PostGIS-PDOK"]), "Output ")
       ("PostGIS-PDOK-CityGML", po::value<std::string>(&outputs["PostGIS-PDOK-CityGML"]), "Output ")
       ("GDAL", po::value<std::string>(&outputs["GDAL"]), "Output ")
@@ -129,7 +127,7 @@ int main(int argc, const char * argv[]) {
     }
     if (vm.count("version")) {
       std::cout << "3dfier " << VERSION << std::endl;
-      std::cout << liblas::GetFullVersion() << std::endl;
+      std::cout << "LASlib " << LAS_TOOLS_VERSION << std::endl;
       std::cout << "GDAL " << GDALVersionInfo("--version") << std::endl;
       //-- TODO : put here the date and/or the git-commit?
       return EXIT_SUCCESS;
@@ -353,6 +351,10 @@ int main(int argc, const char * argv[]) {
       map3d.set_building_radius_vertex_elevation(n["building_radius_vertex_elevation"].as<float>());
     if (n["threshold_jump_edges"])
       map3d.set_threshold_jump_edges(n["threshold_jump_edges"].as<float>());
+    if (n["threshold_bridge_jump_edges"])
+      map3d.set_threshold_bridge_jump_edges(n["threshold_bridge_jump_edges"].as<float>());
+    else if (n["threshold_jump_edges"]) // set threshold_jump_edges same for bridge
+      map3d.set_threshold_bridge_jump_edges(n["threshold_jump_edges"].as<float>());
     if (n["stitching"] && n["stitching"].as<std::string>() == "false")
       bStitching = false;
 
@@ -609,7 +611,7 @@ int main(int argc, const char * argv[]) {
     std::string ofname = output.second;
     if (format != "CityGML-Multifile" && format != "CityGML-IMGeo-Multifile" && format != "CityJSON" &&
       format != "Shapefile" && format != "Shapefile-Multifile" &&
-      format != "PostGIS" && format != "PostGIS-Multi" && format != "PostGIS-PDOK" && format != "PostGIS-PDOK-CityGML" &&
+      format != "PostGIS" && format != "PostGIS-PDOK" && format != "PostGIS-PDOK-CityGML" &&
       format != "GDAL") {
       of.open(ofname);
     }
@@ -663,19 +665,15 @@ int main(int argc, const char * argv[]) {
     }
     else if (format == "PostGIS") {
       std::clog << "PostGIS output\n";
-      fileWritten = map3d.get_gdal_output(ofname, "PostgreSQL", false);
-    }
-    else if (format == "PostGIS-Multi") {
-      std::clog << "PostGIS multiple table output\n";
-      fileWritten = map3d.get_gdal_output(ofname, "PostgreSQL", true);
+      fileWritten = map3d.get_postgis_output(ofname, false, false);
     }
     else if (format == "PostGIS-PDOK") {
       std::clog << "PostGIS with IMGeo GML string output\n";
-      fileWritten = map3d.get_pdok_output(ofname);
+      fileWritten = map3d.get_postgis_output(ofname, true,  false);
     }
     else if (format == "PostGIS-PDOK-CityGML") {
       std::clog << "PostGIS with CityGML string output\n";
-      fileWritten = map3d.get_pdok_citygml_output(ofname);
+      fileWritten = map3d.get_postgis_output(ofname, true, true);
     }
     else if (format == "GDAL") { //-- TODO: what is this? a path? how to use?
       if (nodes["output"] && nodes["output"]["gdal_driver"]) {
@@ -704,7 +702,7 @@ std::string print_license() {
   std::string thelicense =
     "======================================================================\n"
     "\n3dfier: takes 2D GIS datasets and '3dfies' to create 3D city models.\n\n"
-    "Copyright (C) 2015-2018  3D geoinformation research group, TU Delft\n\n"
+    "Copyright (C) 2015-2019  3D geoinformation research group, TU Delft\n\n"
     "3dfier is free software: you can redistribute it and/or modify\n"
     "it under the terms of the GNU General Public License as published by\n"
     "the Free Software Foundation, either version 3 of the License, or\n"
@@ -1149,6 +1147,15 @@ bool validate_yaml(const char* arg, std::set<std::string>& allowedFeatures) {
         std::cerr << "\tOption 'options.threshold_jump_edges' invalid.\n";
       }
     }
+    if (n["threshold_bridge_jump_edges"]) {
+      try {
+        boost::lexical_cast<float>(n["threshold_bridge_jump_edges"].as<std::string>());
+      }
+      catch (boost::bad_lexical_cast& e) {
+        wentgood = false;
+        std::cerr << "\tOption 'options.threshold_bridge_jump_edges' invalid.\n";
+      }
+    }
     if (n["stitching"]) {
       std::string s = n["stitching"].as<std::string>();
       if ((s != "true") && (s != "false")) {
diff --git a/thirdparty/polyfit.hpp b/thirdparty/polyfit.hpp
index b3b08709..e749d519 100644
--- a/thirdparty/polyfit.hpp
+++ b/thirdparty/polyfit.hpp
@@ -94,39 +94,37 @@ std::vector<T> polyvalqr(const std::vector<T>& oCoeff, const std::vector<T>& oX)
 }
 
 template<typename T>
-void combineXY(std::vector<T>& oX, std::vector<T>& oY, mathalgo::matrix<T>& oXY) {
+void combineXY(std::vector<T>& oX, std::vector<T>& oY, T x0, T y0, mathalgo::matrix<T>& oXY) {
   size_t nCount = oX.size();
   size_t nCols = 3 * 2; // 3 unknowns in 2 dimensions
   oXY = mathalgo::matrix<T>(nCount, nCols);
-  
-  // normalize x and y matrix
-  for (size_t i = 1; i < nCount; i++) {
-    oX[i] = oX[i] - oX[0];
-    oY[i] = oY[i] - oY[0];
-  }
-  oX[0] = 0;
-  oY[0] = 0;
 
   // create the XY matrix
   for (size_t nRow = 0; nRow < nCount; nRow++) {
+    T x = oX[nRow] - x0;
+    T y = oY[nRow] - y0;
+    if (nRow == 0) {
+      x = 0;
+      y = 0;
+    }
     oXY(nRow, 0) = 1;
-    oXY(nRow, 1) = oX[nRow];
-    oXY(nRow, 2) = oY[nRow];
-    oXY(nRow, 3) = oX[nRow] * oY[nRow];
-    oXY(nRow, 4) = std::pow(oX[nRow], 2);
-    oXY(nRow, 5) = std::pow(oY[nRow], 2);
+    oXY(nRow, 1) = x;
+    oXY(nRow, 2) = y;
+    oXY(nRow, 3) = x * y;
+    oXY(nRow, 4) = x * x;
+    oXY(nRow, 5) = y * y;
   }
 }
 
 // 3D plane fitting
 template<typename T>
-void polyfit3d(std::vector<T>& oX, std::vector<T>& oY, std::vector<T>& oZ, std::vector<T>& coeffs, std::vector<T>& calculated) {
+void polyfit3d(std::vector<T>& oX, std::vector<T>& oY, std::vector<T>& oZ, T x0, T y0, std::vector<T>& coeffs, std::vector<T>& calculated) {
   if (oX.size() != oY.size() || oX.size() != oZ.size())
     throw std::invalid_argument("X and Y or X and Z vector sizes do not match");
 
   size_t nCount = oX.size();
   mathalgo::matrix<T> A;
-  combineXY(oX, oY, A);
+  combineXY(oX, oY, x0, y0, A);
   mathalgo::matrix<T> X(nCount, 1);
 
   // copy z matrix
@@ -161,16 +159,16 @@ void polyfit3d(std::vector<T>& oX, std::vector<T>& oY, std::vector<T>& oZ, std::
 }
 
 template<typename T>
-void polyval3d(std::vector<T>& x, std::vector<T>& y, std::vector<T>& coeff, std::vector<T>& calculated) {
+void polyval3d(std::vector<T>& x, std::vector<T>& y, T x0, T y0, std::vector<T>& coeff, std::vector<T>& calculated) {
   mathalgo::matrix<T> A;
-  combineXY(x, y, A);
-  mathalgo::matrix<double> Y(coeff.size(), 1);
+  combineXY(x, y, x0, y0, A);
+  mathalgo::matrix<T> YT(coeff.size(), 1);
   // build coeff matrix
   for (size_t i = 0; i < coeff.size(); i++) {
-    Y(i, 0) = coeff[i];
+    YT(i, 0) = coeff[i];
   }
   mathalgo::matrix<T> AY;
-  A.multiply(Y, AY);
+  A.multiply(YT, AY);
   calculated = AY.data();
 }
 #endif
\ No newline at end of file
diff --git a/thirdparty/ptinpoly.c b/thirdparty/ptinpoly.c
deleted file mode 100644
index 7c833924..00000000
--- a/thirdparty/ptinpoly.c
+++ /dev/null
@@ -1,651 +0,0 @@
-/* ptinpoly.c - point in polygon inside/outside code.
-
-   by Eric Haines, 3D/Eye Inc, erich@eye.com
-
-   This code contains the following algorithms:
-	grid testing - grid imposed on polygon
-*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include "ptinpoly.h"
-
-#define X	0
-#define Y	1
-
-#ifndef TRUE
-#define TRUE	1
-#define FALSE	0
-#endif
-
-#ifndef HUGE
-#define HUGE	1.797693134862315e+308
-#endif
-
-#ifndef M_PI
-#define M_PI	3.14159265358979323846
-#endif
-
-/* test if a & b are within epsilon.  Favors cases where a < b */
-#define Near(a,b,eps)	( ((b)-(eps)<(a)) && ((a)-(eps)<(b)) )
-
-#define MALLOC_CHECK( a )	if ( !(a) ) {				   \
-				    fprintf( stderr, "out of memory\n" ) ; \
-				    exit(1) ;				   \
-				}
-
-/* ======= Grid algorithm ================================================= */
-
-/* Impose a grid upon the polygon and test only the local edges against the
- * point.
- *
- * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
- * grid resolution _resolution_ and a pointer to a grid structure _p_gs_.
- * Call testing procedure with a pointer to this array and test point _point_,
- * returns 1 if inside, 0 if outside.
- * Call cleanup with pointer to grid structure to free space.
- */
-
-/* Strategy for setup:
- *   Get bounds of polygon, allocate grid.
- *   "Walk" each edge of the polygon and note which edges have been crossed
- *     and what cells are entered (points on a grid edge are always considered
- *     to be above that edge).	Keep a record of the edges overlapping a cell.
- *     For cells with edges, determine if any cell border has no edges passing
- *     through it and so can be used for shooting a test ray.
- *     Keep track of the parity of the x (horizontal) grid cell borders for
- *     use in determining whether the grid corners are inside or outside.
- */
-void GridSetup(pPipoint pgon[], int	numverts, int	resolution, pGridSet p_gs)
-{
-double	*vtx0, *vtx1, *vtxa, *vtxb, *p_gl ;
-int	i, j, gc_clear_flags ;
-double	vx0, vx1, vy0, vy1, gxdiff, gydiff, eps ;
-pGridCell	p_gc, p_ngc ;
-double	xdiff, ydiff, tmax, inv_x, inv_y, xdir, ydir, t_near, tx, ty ;
-double	tgcx, tgcy ;
-int	gcx, gcy, sign_x ;
-int	y_flag, io_state ;
-
-    p_gs->xres = p_gs->yres = resolution ;
-    p_gs->tot_cells = p_gs->xres * p_gs->yres ;
-    p_gs->glx = (double *)malloc( (p_gs->xres+1) * sizeof(double));
-    MALLOC_CHECK( p_gs->glx ) ;
-    p_gs->gly = (double *)malloc( (p_gs->yres+1) * sizeof(double));
-    MALLOC_CHECK( p_gs->gly ) ;
-    p_gs->gc = (pGridCell)malloc( p_gs->tot_cells * sizeof(GridCell));
-    MALLOC_CHECK( p_gs->gc ) ;
-
-    p_gs->minx =
-    p_gs->maxx = pgon[0]->x ;
-    p_gs->miny =
-    p_gs->maxy = pgon[0]->y ;
-
-    /* find bounds of polygon */
-    for ( i = 1 ; i < numverts ; i++ ) {
-	vx0 = pgon[i]->x ;
-	if ( p_gs->minx > vx0 ) {
-	    p_gs->minx = vx0 ;
-	} else if ( p_gs->maxx < vx0 ) {
-	    p_gs->maxx = vx0 ;
-	}
-
-	vy0 = pgon[i]->y ;
-	if ( p_gs->miny > vy0 ) {
-	    p_gs->miny = vy0 ;
-	} else if ( p_gs->maxy < vy0 ) {
-	    p_gs->maxy = vy0 ;
-	}
-    }
-
-    /* add a little to the bounds to ensure everything falls inside area */
-    gxdiff = p_gs->maxx - p_gs->minx ;
-    gydiff = p_gs->maxy - p_gs->miny ;
-    p_gs->minx -= EPSILON * gxdiff ;
-    p_gs->maxx += EPSILON * gxdiff ;
-    p_gs->miny -= EPSILON * gydiff ;
-    p_gs->maxy += EPSILON * gydiff ;
-
-    /* avoid roundoff problems near corners by not getting too close to them */
-    eps = 1e-9 * ( gxdiff + gydiff ) ;
-
-    /* use the new bounds to compute cell widths */
-    TryAgain:
-    p_gs->xdelta =
-	    (p_gs->maxx-p_gs->minx) / (double)p_gs->xres ;
-    p_gs->inv_xdelta = 1.0 / p_gs->xdelta ;
-
-    p_gs->ydelta =
-	    (p_gs->maxy-p_gs->miny) / (double)p_gs->yres ;
-    p_gs->inv_ydelta = 1.0 / p_gs->ydelta ;
-
-    for ( i = 0, p_gl = p_gs->glx ; i < p_gs->xres ; i++ ) {
-	*p_gl++ = p_gs->minx + i * p_gs->xdelta ;
-    }
-    /* make last grid corner precisely correct */
-    *p_gl = p_gs->maxx ;
-
-    for ( i = 0, p_gl = p_gs->gly ; i < p_gs->yres ; i++ ) {
-	*p_gl++ = p_gs->miny + i * p_gs->ydelta ;
-    }
-    *p_gl = p_gs->maxy ;
-
-    for ( i = 0, p_gc = p_gs->gc ; i < p_gs->tot_cells ; i++, p_gc++ ) {
-	p_gc->tot_edges = 0 ;
-	p_gc->gc_flags = 0x0 ;
-	p_gc->gr = NULL ;
-    }
-
-    /* loop through edges and insert into grid structure */
-    vtx0 = pgon[numverts-1] ;
-    for ( i = 0 ; i < numverts ; i++ ) {
-	vtx1 = pgon[i] ;
-
-	if ( vtx0[Y] < vtx1[Y] ) {
-	    vtxa = vtx0 ;
-	    vtxb = vtx1 ;
-	} else {
-	    vtxa = vtx1 ;
-	    vtxb = vtx0 ;
-	}
-
-	/* Set x variable for the direction of the ray */
-	xdiff = vtxb[X] - vtxa[X] ;
-	ydiff = vtxb[Y] - vtxa[Y] ;
-	tmax = sqrt( xdiff * xdiff + ydiff * ydiff ) ;
-
-	/* if edge is of 0 length, ignore it (useless edge) */
-	if ( tmax == 0.0 ) goto NextEdge ;
-
-	xdir = xdiff / tmax ;
-	ydir = ydiff / tmax ;
-
-	gcx = (int)(( vtxa[X] - p_gs->minx ) * p_gs->inv_xdelta) ;
-	gcy = (int)(( vtxa[Y] - p_gs->miny ) * p_gs->inv_ydelta) ;
-
-	/* get information about slopes of edge, etc */
-	if ( vtxa[X] == vtxb[X] ) {
-	    sign_x = 0 ;
-	    tx = HUGE ;
-	} else {
-	    inv_x = tmax / xdiff ;
-	    tx = p_gs->xdelta * (double)gcx + p_gs->minx - vtxa[X] ;
-	    if ( vtxa[X] < vtxb[X] ) {
-		sign_x = 1 ;
-		tx += p_gs->xdelta ;
-		tgcx = p_gs->xdelta * inv_x ;
-	    } else {
-		sign_x = -1 ;
-		tgcx = -p_gs->xdelta * inv_x ;
-	    }
-	    tx *= inv_x ;
-	}
-
-	if ( vtxa[Y] == vtxb[Y] ) {
-	    ty = HUGE ;
-	} else {
-	    inv_y = tmax / ydiff ;
-	    ty = (p_gs->ydelta * (double)(gcy+1) + p_gs->miny - vtxa[Y])
-		* inv_y ;
-	    tgcy = p_gs->ydelta * inv_y ;
-	}
-
-	p_gc = &p_gs->gc[gcy*p_gs->xres+gcx] ;
-
-	vx0 = vtxa[X] ;
-	vy0 = vtxa[Y] ;
-
-	t_near = 0.0 ;
-
-	do {
-	    /* choose the next boundary, but don't move yet */
-	    if ( tx <= ty ) {
-		gcx += sign_x ;
-
-		ty -= tx ;
-		t_near += tx ;
-		tx = tgcx ;
-
-		/* note which edge is hit when leaving this cell */
-		if ( t_near < tmax ) {
-		    if ( sign_x > 0 ) {
-			p_gc->gc_flags |= GC_R_EDGE_HIT ;
-			vx1 = p_gs->glx[gcx] ;
-		    } else {
-			p_gc->gc_flags |= GC_L_EDGE_HIT ;
-			vx1 = p_gs->glx[gcx+1] ;
-		    }
-
-		    /* get new location */
-		    vy1 = t_near * ydir + vtxa[Y] ;
-		} else {
-		    /* end of edge, so get exact value */
-		    vx1 = vtxb[X] ;
-		    vy1 = vtxb[Y] ;
-		}
-
-		y_flag = FALSE ;
-
-	    } else {
-
-		gcy++ ;
-
-		tx -= ty ;
-		t_near += ty ;
-		ty = tgcy ;
-
-		/* note top edge is hit when leaving this cell */
-		if ( t_near < tmax ) {
-		    p_gc->gc_flags |= GC_T_EDGE_HIT ;
-		    /* this toggles the parity bit */
-		    p_gc->gc_flags ^= GC_T_EDGE_PARITY ;
-
-		    /* get new location */
-		    vx1 = t_near * xdir + vtxa[X] ;
-		    vy1 = p_gs->gly[gcy] ;
-		} else {
-		    /* end of edge, so get exact value */
-		    vx1 = vtxb[X] ;
-		    vy1 = vtxb[Y] ;
-		}
-
-		y_flag = TRUE ;
-	    }
-
-	    /* check for corner crossing, then mark the cell we're in */
-	    if ( !AddGridRecAlloc( p_gc, vx0, vy0, vx1, vy1, eps ) ) {
-		/* warning, danger - we have just crossed a corner.
-		 * There are all kinds of topological messiness we could
-		 * do to get around this case, but they're a headache.
-		 * The simplest recovery is just to change the extents a bit
-		 * and redo the meshing, so that hopefully no edges will
-		 * perfectly cross a corner.  Since it's a preprocess, we
-		 * don't care too much about the time to do it.
-		 */
-
-		/* clean out all grid records */
-		for ( i = 0, p_gc = p_gs->gc
-		    ; i < p_gs->tot_cells
-		    ; i++, p_gc++ ) {
-
-		    if ( p_gc->gr ) {
-			free( p_gc->gr ) ;
-		    }
-		}
-
-		/* make the bounding box ever so slightly larger, hopefully
-		 * changing the alignment of the corners.
-		 */
-		p_gs->minx -= EPSILON * gxdiff * 0.24 ;
-		p_gs->miny -= EPSILON * gydiff * 0.10 ;
-
-		/* yes, it's the dreaded goto - run in fear for your lives! */
-		goto TryAgain ;
-	    }
-
-	    if ( t_near < tmax ) {
-		/* note how we're entering the next cell */
-		/* TBD: could be done faster by incrementing index in the
-		 * incrementing code, above */
-		p_gc = &p_gs->gc[gcy*p_gs->xres+gcx] ;
-
-		if ( y_flag ) {
-		    p_gc->gc_flags |= GC_B_EDGE_HIT ;
-		    /* this toggles the parity bit */
-		    p_gc->gc_flags ^= GC_B_EDGE_PARITY ;
-		} else {
-		    p_gc->gc_flags |=
-			( sign_x > 0 ) ? GC_L_EDGE_HIT : GC_R_EDGE_HIT ;
-		}
-	    }
-
-	    vx0 = vx1 ;
-	    vy0 = vy1 ;
-	}
-	/* have we gone further than the end of the edge? */
-	while ( t_near < tmax ) ;
-
-	NextEdge:
-	vtx0 = vtx1 ;
-    }
-
-    /* the grid is all set up, now set up the inside/outside value of each
-     * corner.
-     */
-    p_gc = p_gs->gc ;
-    p_ngc = &p_gs->gc[p_gs->xres] ;
-
-    /* we know the bottom and top rows are all outside, so no flag is set */
-    for ( i = 1; i < p_gs->yres ; i++ ) {
-	/* start outside */
-	io_state = 0x0 ;
-
-	for ( j = 0; j < p_gs->xres ; j++ ) {
-
-	    if ( io_state ) {
-		/* change cell left corners to inside */
-		p_gc->gc_flags |= GC_TL_IN ;
-		p_ngc->gc_flags |= GC_BL_IN ;
-	    }
-
-	    if ( p_gc->gc_flags & GC_T_EDGE_PARITY ) {
-		io_state = !io_state ;
-	    }
-
-	    if ( io_state ) {
-		/* change cell right corners to inside */
-		p_gc->gc_flags |= GC_TR_IN ;
-		p_ngc->gc_flags |= GC_BR_IN ;
-	    }
-
-	    p_gc++ ;
-	    p_ngc++ ;
-	}
-    }
-
-    p_gc = p_gs->gc ;
-    for ( i = 0; i < p_gs->tot_cells ; i++ ) {
-
-	/* reverse parity of edge clear (1==edge clear) */
-	gc_clear_flags = p_gc->gc_flags ^ GC_ALL_EDGE_CLEAR ;
-	if ( gc_clear_flags & GC_L_EDGE_CLEAR ) {
-	    p_gc->gc_flags |= GC_AIM_L ;
-	} else
-	if ( gc_clear_flags & GC_B_EDGE_CLEAR ) {
-	    p_gc->gc_flags |= GC_AIM_B ;
-	} else
-	if ( gc_clear_flags & GC_R_EDGE_CLEAR ) {
-	    p_gc->gc_flags |= GC_AIM_R ;
-	} else
-	if ( gc_clear_flags & GC_T_EDGE_CLEAR ) {
-	    p_gc->gc_flags |= GC_AIM_T ;
-	} else {
-	    /* all edges have something on them, do full test */
-	    p_gc->gc_flags |= GC_AIM_C ;
-	}
-	p_gc++ ;
-    }
-}
-
-int AddGridRecAlloc(pGridCell p_gc, double xa, double ya, double xb, double yb, double eps )
-{
-pGridRec	p_gr ;
-double		slope, inv_slope ;
-
-    if ( Near(ya, yb, eps) ) {
-	if ( Near(xa, xb, eps) ) {
-	    /* edge is 0 length, so get rid of it */
-	    return( FALSE ) ;
-	} else {
-	    /* horizontal line */
-	    slope = HUGE ;
-	    inv_slope = 0.0 ;
-	}
-    } else {
-	if ( Near(xa, xb, eps) ) {
-	    /* vertical line */
-	    slope = 0.0 ;
-	    inv_slope = HUGE ;
-	} else {
-	    slope = (xb-xa)/(yb-ya) ;
-	    inv_slope = (yb-ya)/(xb-xa) ;
-	}
-    }
-
-    p_gc->tot_edges++ ;
-    if ( p_gc->tot_edges <= 1 ) {
-	p_gc->gr = (pGridRec)malloc( sizeof(GridRec) ) ;
-    } else {
-	p_gc->gr = (pGridRec)realloc( p_gc->gr,
-		p_gc->tot_edges * sizeof(GridRec) ) ;
-    }
-    MALLOC_CHECK( p_gc->gr ) ;
-    p_gr = &p_gc->gr[p_gc->tot_edges-1] ;
-
-    p_gr->slope = slope ;
-    p_gr->inv_slope = inv_slope ;
-
-    p_gr->xa = xa ;
-    p_gr->ya = ya ;
-    if ( xa <= xb ) {
-	p_gr->minx = xa ;
-	p_gr->maxx = xb ;
-    } else {
-	p_gr->minx = xb ;
-	p_gr->maxx = xa ;
-    }
-    if ( ya <= yb ) {
-	p_gr->miny = ya ;
-	p_gr->maxy = yb ;
-    } else {
-	p_gr->miny = yb ;
-	p_gr->maxy = ya ;
-    }
-
-    /* P2 - P1 */
-    p_gr->ax = xb - xa ;
-    p_gr->ay = yb - ya ;
-
-    return( TRUE ) ;
-}
-
-/* Test point against grid and edges in the cell (if any).  Algorithm:
- *    Check bounding box; if outside then return.
- *    Check cell point is inside; if simple inside or outside then return.
- *    Find which edge or corner is considered to be the best for testing and
- *	  send a test ray towards it, counting the crossings.  Add in the
- *	  state of the edge or corner the ray went to and so determine the
- *	  state of the point (inside or outside).
- */
-int GridTest(pGridSet p_gs, pPipoint point)
-{
-int	j, count, init_flag ;
-pGridCell	p_gc ;
-pGridRec	p_gr ;
-double	tx, ty, xcell, ycell, bx,by,cx,cy, cornerx, cornery ;
-double	alpha, beta, denom ;
-unsigned short	gc_flags ;
-int	inside_flag ;
-
-    /* first, is point inside bounding rectangle? */
-    if ( ( ty = point->y ) < p_gs->miny ||
-	 ty >= p_gs->maxy ||
-	 ( tx = point->x ) < p_gs->minx ||
-	 tx >= p_gs->maxx ) {
-
-	/* outside of box */
-	inside_flag = FALSE ;
-    } else {
-
-	/* what cell are we in? */
-	ycell = ( ty - p_gs->miny ) * p_gs->inv_ydelta ;
-	xcell = ( tx - p_gs->minx ) * p_gs->inv_xdelta ;
-	p_gc = &p_gs->gc[((int)ycell)*p_gs->xres + (int)xcell] ;
-
-	/* is cell simple? */
-	count = p_gc->tot_edges ;
-	if ( count ) {
-	    /* no, so find an edge which is free. */
-	    gc_flags = p_gc->gc_flags ;
-	    p_gr = p_gc->gr ;
-	    switch( gc_flags & GC_AIM ) {
-	    case GC_AIM_L:
-		/* left edge is clear, shoot X- ray */
-		/* note - this next statement requires that GC_BL_IN is 1 */
-		inside_flag = gc_flags & GC_BL_IN ;
-		for ( j = count+1 ; --j ; p_gr++ ) {
-		    /* test if y is between edges */
-		    if ( ty >= p_gr->miny && ty < p_gr->maxy ) {
-			if ( tx > p_gr->maxx ) {
-			    inside_flag = !inside_flag ;
-			} else if ( tx > p_gr->minx ) {
-			    /* full computation */
-			    if ( ( p_gr->xa -
-				( p_gr->ya - ty ) * p_gr->slope ) < tx ) {
-				inside_flag = !inside_flag ;
-			    }
-			}
-		    }
-		}
-		break ;
-
-	    case GC_AIM_B:
-		/* bottom edge is clear, shoot Y+ ray */
-		/* note - this next statement requires that GC_BL_IN is 1 */
-		inside_flag = gc_flags & GC_BL_IN ;
-		for ( j = count+1 ; --j ; p_gr++ ) {
-		    /* test if x is between edges */
-		    if ( tx >= p_gr->minx && tx < p_gr->maxx ) {
-			if ( ty > p_gr->maxy ) {
-			    inside_flag = !inside_flag ;
-			} else if ( ty > p_gr->miny ) {
-			    /* full computation */
-			    if ( ( p_gr->ya - ( p_gr->xa - tx ) *
-				    p_gr->inv_slope ) < ty ) {
-				inside_flag = !inside_flag ;
-			    }
-			}
-		    }
-		}
-		break ;
-
-	    case GC_AIM_R:
-		/* right edge is clear, shoot X+ ray */
-		inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ;
-
-		/* TBD: Note, we could have sorted the edges to be tested
-		 * by miny or somesuch, and so be able to cut testing
-		 * short when the list's miny > point.y .
-		 */
-		for ( j = count+1 ; --j ; p_gr++ ) {
-		    /* test if y is between edges */
-		    if ( ty >= p_gr->miny && ty < p_gr->maxy ) {
-			if ( tx <= p_gr->minx ) {
-			    inside_flag = !inside_flag ;
-			} else if ( tx <= p_gr->maxx ) {
-			    /* full computation */
-			    if ( ( p_gr->xa -
-				( p_gr->ya - ty ) * p_gr->slope ) >= tx ) {
-				inside_flag = !inside_flag ;
-			    }
-			}
-		    }
-		}
-		break ;
-
-	    case GC_AIM_T:
-		/* top edge is clear, shoot Y+ ray */
-		inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ;
-		for ( j = count+1 ; --j ; p_gr++ ) {
-		    /* test if x is between edges */
-		    if ( tx >= p_gr->minx && tx < p_gr->maxx ) {
-			if ( ty <= p_gr->miny ) {
-			    inside_flag = !inside_flag ;
-			} else if ( ty <= p_gr->maxy ) {
-			    /* full computation */
-			    if ( ( p_gr->ya - ( p_gr->xa - tx ) *
-				    p_gr->inv_slope ) >= ty ) {
-				inside_flag = !inside_flag ;
-			    }
-			}
-		    }
-		}
-		break ;
-
-	    case GC_AIM_C:
-		/* no edge is clear, bite the bullet and test
-		 * against the bottom left corner.
-		 * We use Franklin Antonio's algorithm (Graphics Gems III).
-		 */
-		/* TBD: Faster yet might be to test against the closest
-		 * corner to the cell location, but our hope is that we
-		 * rarely need to do this testing at all.
-		 */
-		inside_flag = ((gc_flags & GC_BL_IN) == GC_BL_IN) ;
-		init_flag = TRUE ;
-
-		/* get lower left corner coordinate */
-		cornerx = p_gs->glx[(int)xcell] ;
-		cornery = p_gs->gly[(int)ycell] ;
-		for ( j = count+1 ; --j ; p_gr++ ) {
-
-		    /* quick out test: if test point is
-		     * less than minx & miny, edge cannot overlap.
-		     */
-		    if ( tx >= p_gr->minx && ty >= p_gr->miny ) {
-
-			/* quick test failed, now check if test point and
-			 * corner are on different sides of edge.
-			 */
-			if ( init_flag ) {
-			    /* Compute these at most once for test */
-			    /* P3 - P4 */
-			    bx = tx - cornerx ;
-			    by = ty - cornery ;
-			    init_flag = FALSE ;
-			}
-			denom = p_gr->ay * bx - p_gr->ax * by ;
-			if ( denom != 0.0 ) {
-			    /* lines are not collinear, so continue */
-			    /* P1 - P3 */
-			    cx = p_gr->xa - tx ;
-			    cy = p_gr->ya - ty ;
-			    alpha = by * cx - bx * cy ;
-			    if ( denom > 0.0 ) {
-				if ( alpha < 0.0 || alpha >= denom ) {
-				    /* test edge not hit */
-				    goto NextEdge ;
-				}
-				beta = p_gr->ax * cy - p_gr->ay * cx ;
-				if ( beta < 0.0 || beta >= denom ) {
-				    /* polygon edge not hit */
-				    goto NextEdge ;
-				}
-			    } else {
-				if ( alpha > 0.0 || alpha <= denom ) {
-				    /* test edge not hit */
-				    goto NextEdge ;
-				}
-				beta = p_gr->ax * cy - p_gr->ay * cx ;
-				if ( beta > 0.0 || beta <= denom ) {
-				    /* polygon edge not hit */
-				    goto NextEdge ;
-				}
-			    }
-			    inside_flag = !inside_flag ;
-			}
-
-		    }
-		    NextEdge: ;
-		}
-		break ;
-	    }
-	} else {
-	    /* simple cell, so if lower left corner is in,
-	     * then cell is inside.
-	     */
-	    inside_flag = p_gc->gc_flags & GC_BL_IN ;
-	}
-    }
-
-    return( inside_flag ) ;
-}
-
-void GridCleanup(pGridSet	p_gs)
-{
-int	i ;
-pGridCell	p_gc ;
-
-    for ( i = 0, p_gc = p_gs->gc
-	; i < p_gs->tot_cells
-	; i++, p_gc++ ) {
-
-	if ( p_gc->gr ) {
-	    free( p_gc->gr ) ;
-	}
-    }
-    free( p_gs->glx ) ;
-    free( p_gs->gly ) ;
-    free( p_gs->gc ) ;
-}
diff --git a/thirdparty/ptinpoly.h b/thirdparty/ptinpoly.h
deleted file mode 100644
index 7a0d26be..00000000
--- a/thirdparty/ptinpoly.h
+++ /dev/null
@@ -1,94 +0,0 @@
-/* ptinpoly.h - point in polygon inside/outside algorithms header file.
- *
- * by Eric Haines, 3D/Eye Inc, erich@eye.com
- */
-
-#ifdef __cplusplus
-extern "C"
-{
-#endif
-/* =========================== System Related ============================= */
-/* SRAN initializes random number generator, if needed */
-#define SRAN()		srand(1)
-/* RAN01 returns a double from [0..1) */
-#define RAN01()		((double)rand() / 32768.0)
-
-/* =========== Grid stuff ================================================= */
-
-#define GR_FULL_VERT	0x01	/* line crosses vertically */
-#define GR_FULL_HORZ	0x02	/* line crosses horizontally */
-
-typedef struct {
-    double	xa,ya ;
-    double	minx, maxx, miny, maxy ;
-    double	ax, ay ;
-    double	slope, inv_slope ;
-} GridRec, *pGridRec;
-
-#define GC_BL_IN	0x0001	/* bottom left corner is in (else out) */
-#define GC_BR_IN	0x0002	/* bottom right corner is in (else out) */
-#define GC_TL_IN	0x0004	/* top left corner is in (else out) */
-#define GC_TR_IN	0x0008	/* top right corner is in (else out) */
-#define GC_L_EDGE_HIT	0x0010	/* left edge is crossed */
-#define GC_R_EDGE_HIT	0x0020	/* right edge is crossed */
-#define GC_B_EDGE_HIT	0x0040	/* bottom edge is crossed */
-#define GC_T_EDGE_HIT	0x0080	/* top edge is crossed */
-#define GC_B_EDGE_PARITY	0x0100	/* bottom edge parity */
-#define GC_T_EDGE_PARITY	0x0200	/* top edge parity */
-#define GC_AIM_L	(0<<10) /* aim towards left edge */
-#define GC_AIM_B	(1<<10) /* aim towards bottom edge */
-#define GC_AIM_R	(2<<10) /* aim towards right edge */
-#define GC_AIM_T	(3<<10) /* aim towards top edge */
-#define GC_AIM_C	(4<<10) /* aim towards a corner */
-#define GC_AIM		0x1c00
-
-#define GC_L_EDGE_CLEAR GC_L_EDGE_HIT
-#define GC_R_EDGE_CLEAR GC_R_EDGE_HIT
-#define GC_B_EDGE_CLEAR GC_B_EDGE_HIT
-#define GC_T_EDGE_CLEAR GC_T_EDGE_HIT
-
-#define GC_ALL_EDGE_CLEAR	(GC_L_EDGE_HIT | \
-				 GC_R_EDGE_HIT | \
-				 GC_B_EDGE_HIT | \
-				 GC_T_EDGE_HIT )
-
-typedef struct {
-    short		tot_edges ;
-    unsigned short	gc_flags ;
-    GridRec		*gr ;
-} GridCell, *pGridCell;
-
-typedef struct {
-    int		xres, yres ;	/* grid size */
-    int		tot_cells ;	/* xres * yres */
-    double	minx, maxx, miny, maxy ;	/* bounding box */
-    double	xdelta, ydelta ;
-    double	inv_xdelta, inv_ydelta ;
-    double	*glx, *gly ;
-    GridCell	*gc ;
-} GridSet, *pGridSet ;
-
-typedef struct {
-  double x, y;
-} Pipoint, *pPipoint;
-
-/* add a little to the limits of the polygon bounding box to avoid precision
-* problems.
-*/
-#define EPSILON		0.00001
-
-/* The following structure is associated with a polygon */
-typedef struct {
-    int		id ;		/* vertex number of edge */
-    int		full_cross ;	/* 1 if extends from top to bottom */
-    double	minx, maxx ;	/* X bounds for bin */
-} Edge, *pEdge ;
-
-void GridSetup(pPipoint pgon[], int	numverts, int	resolution, pGridSet p_gs);
-int AddGridRecAlloc(pGridCell p_gc, double xa, double ya, double xb, double yb, double eps);
-void GridCleanup(pGridSet	p_gs);
-int GridTest(pGridSet p_gs, pPipoint point);
-
-#ifdef __cplusplus
-}
-#endif
\ No newline at end of file
diff --git a/vs_build/3dfier.vcxproj b/vs_build/3dfier.vcxproj
index e341354a..61fcf733 100644
--- a/vs_build/3dfier.vcxproj
+++ b/vs_build/3dfier.vcxproj
@@ -1,4 +1,4 @@
-<?xml version="1.0" encoding="utf-8"?>
+<?xml version="1.0" encoding="utf-8"?>
 <Project DefaultTargets="Build" ToolsVersion="14.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
   <ItemGroup Label="ProjectConfigurations">
     <ProjectConfiguration Include="Debug|x64">
@@ -58,7 +58,7 @@
   </PropertyGroup>
   <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
     <ClCompile>
-      <AdditionalIncludeDirectories>$(LIBLAS_ROOT)\include;$(LASZIP_ROOT)\include;$(BOOST_ROOT);$(OSGEO4W_ROOT)\include;$(YAML-CPP_DIR)\include;$(CGAL_DIR)\include;$(CGAL_DIR)\auxiliary\gmp\include;..\thirdparty;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>$(LASLIB_ROOT)\inc;$(LASZIP_ROOT)\src;$(BOOST_ROOT);$(OSGEO4W_ROOT)\include;$(YAML-CPP_DIR)\include;$(CGAL_DIR)\include;$(CGAL_DIR)\auxiliary\gmp\include;..\thirdparty;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
       <AssemblerListingLocation>$(IntDir)</AssemblerListingLocation>
       <CompileAs>Default</CompileAs>
       <ExceptionHandling>Sync</ExceptionHandling>
@@ -89,8 +89,8 @@
     </Midl>
     <Link>
       <AdditionalOptions> /machine:x64 %(AdditionalOptions)</AdditionalOptions>
-      <AdditionalDependencies>gdal_i.lib;liblas.lib;laszip.lib;libCGAL_Core-vc140-mt-4.8.lib;libCGAL-vc140-mt-4.8.lib;libgmp-10.lib;libyaml-cppmd.lib</AdditionalDependencies>
-      <AdditionalLibraryDirectories>$(LIBLAS_ROOT)\bin;$(LIBLAS_ROOT)\lib;$(LASZIP_ROOT)\lib;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
+      <AdditionalDependencies>gdal_i.lib;LASlib.lib;libCGAL_Core-vc140-mt-4.12.lib;libCGAL-vc140-mt-4.12.lib;libgmp-10.lib;libyaml-cppmd.lib</AdditionalDependencies>
+      <AdditionalLibraryDirectories>$(LASLIB_ROOT)\lib\Release;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
       <GenerateDebugInformation>true</GenerateDebugInformation>
       <IgnoreSpecificDefaultLibraries>%(IgnoreSpecificDefaultLibraries)</IgnoreSpecificDefaultLibraries>
       <ImportLibrary>
@@ -106,7 +106,7 @@
   </ItemDefinitionGroup>
   <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
     <ClCompile>
-      <AdditionalIncludeDirectories>$(LIBLAS_ROOT)\include;$(LASZIP_ROOT)\include;$(BOOST_ROOT);$(OSGEO4W_ROOT)\include;$(YAML-CPP_DIR)\include;$(CGAL_DIR)\include;$(CGAL_DIR)\auxiliary\gmp\include;..\thirdparty;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>$(LASLIB_ROOT)\inc;$(LASZIP_ROOT)\src;$(BOOST_ROOT);$(OSGEO4W_ROOT)\include;$(YAML-CPP_DIR)\include;$(CGAL_DIR)\include;$(CGAL_DIR)\auxiliary\gmp\include;..\thirdparty;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
       <AssemblerListingLocation>$(IntDir)</AssemblerListingLocation>
       <CompileAs>Default</CompileAs>
       <ExceptionHandling>Sync</ExceptionHandling>
@@ -137,8 +137,8 @@
       <ProxyFileName>%(Filename)_p.c</ProxyFileName>
     </Midl>
     <Link>
-      <AdditionalDependencies>gdal_i.lib;liblas.lib;laszip.lib;libCGAL_Core-vc140-mt-4.8.lib;libCGAL-vc140-mt-4.8.lib;libgmp-10.lib;libyaml-cppmd.lib</AdditionalDependencies>
-      <AdditionalLibraryDirectories>$(LIBLAS_ROOT)\bin;$(LIBLAS_ROOT)\lib;$(LASZIP_ROOT)\lib;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
+      <AdditionalDependencies>gdal_i.lib;LASlib.lib;libCGAL_Core-vc140-mt-4.12.lib;libCGAL-vc140-mt-4.12.lib;libgmp-10.lib;libyaml-cppmd.lib</AdditionalDependencies>
+      <AdditionalLibraryDirectories>$(LASLIB_ROOT)\lib\Release;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
       <GenerateDebugInformation>true</GenerateDebugInformation>
       <IgnoreSpecificDefaultLibraries>%(IgnoreSpecificDefaultLibraries)</IgnoreSpecificDefaultLibraries>
       <ImportLibrary>