Merge branch 'release-0.4.0' into main

.github/workflows/build.yml

@@ -15,7 +15,7 @@ jobs:
       - name: Install dependencies
         run: |
             sudo apt-get update
-            sudo apt-get install libcgal-dev libboost-all-dev libeigen3-dev libomp-dev
+            sudo apt-get install libcgal-dev libboost-all-dev libeigen3-dev libomp-dev libgdal-dev
       - name: Build
         run: |
             mkdir build && cd build
@@ -28,7 +28,7 @@ jobs:
       - name: Install dependencies
         run: |
             brew update
-            brew install cmake boost cgal eigen libomp
+            brew install cmake boost cgal eigen libomp gdal
       - name: Build
         run: |
             mkdir build && cd build

CHANGELOG.md

@@ -1,5 +1,13 @@
 # Changelog
 
+## [0.4.0] - 2023-05-03
+### Added
+- (breaking) Support for many polygon formats (through GDAL)
+- Optional flag that enforces buildings and terrain intersection
+### Changed
+- Updated LAStools to v2.0.2
+- Minor bugfixes
+
 ## [0.3.0] - 2023-01-18
 ### Added
 - Building surface refinement

CMakeLists.txt

@@ -26,6 +26,13 @@ else()
     return()
 endif ()
 
+# GDAL
+find_package(GDAL REQUIRED 3.0)
+if (NOT GDAL_FOUND)
+    message(SEND_ERROR "City4CFD requires the GDAL library, check the README for configuration")
+endif()
+include_directories(${GDAL_INCLUDE_DIR})
+
 # Eigen
 find_package(Eigen3 3.2.0) #(requires 3.2.0 or greater)
 include(${CMAKE_SOURCE_DIR}/thirdparty/CGAL/cmake/CGAL_Eigen3_support.cmake)
@@ -59,6 +66,7 @@ set_target_properties(
 target_link_libraries(city4cfd
         ${CGAL_LIBRARIES}
         ${CGAL_3RD_PARTY_LIBRARIES}
+        ${GDAL_LIBRARY}
         ${Boost_SYSTEM_LIBRARY}
         ${Boost_FILESYSTEM_LIBRARY}
         ${Boost_LOCALE_LIBRARY}
@@ -67,4 +75,4 @@ target_link_libraries(city4cfd
         ${OpenMP_support}
         )
 
-install(TARGETS city4cfd DESTINATION bin)
+install(TARGETS city4cfd DESTINATION bin)

README.md

@@ -27,7 +27,7 @@ City4CFD is developed by the [3D Geoinformation Research Group](https://3d.bk.tu
 ## Data formats
 **Point clouds** can be imported in LAS/LAZ, TXT/XYZ, or PLY format. We ask separately for ground and building points. While some datasets contain building-ground classification, some do not. Our [point cloud preparation tool](https://github.com/ipadjen/City4CFD_doc/wiki/Point-clouds#automatic-preparation) can extract ground and building points from user-defined classes, or use the [Cloth Simulation Filter](http://ramm.bnu.edu.cn/projects/CSF/) to separate the ground and non-ground points. If you would like to check your points, see if they are classified, or even conduct the filtering and classification yourself, we suggest you use [CloudCompare](https://www.danielgm.net/cc/).
 
-**2D data** (polygons) are imported in GeoJSON format. For all pre-processing related to polygons, including conversion to GeoJSON, you can use [QGIS](https://qgis.org/en/site/).
+**2D data** (polygons) are imported in [GDAL-supported formats](https://gdal.org/drivers/vector/index.html). For all pre-processing related to polygons you can use [QGIS](https://qgis.org/en/site/).
 
 **Geometry import** supports the following formats: OBJ, STL, PLY, OFF, VTP, and CityJSON.
 
@@ -41,18 +41,19 @@ The following libraries are required to build the project:
 - [CGAL](https://www.cgal.org/) version 5
 - Boost >= 1.66
 - Eigen >= 3.2
+- GDAL >= 3.0
 
 *OpenMP* is an optional dependency.
 
 All dependencies are generally available in Linux distributions, e.g. in Debian/Ubuntu/Mint:
 ```
-sudo apt-get install libcgal-dev libboost-all-dev libeigen3-dev libomp-dev
+sudo apt-get install libcgal-dev libboost-all-dev libeigen3-dev libomp-dev libgdal-dev
 ```
 
 In macOS you can install dependencies with Homebrew:
 
 ```
-brew install cmake boost cgal eigen libomp
+brew install cmake boost cgal eigen libomp gdal
 ```
 
 The project uses CMake to generate makefiles, so make sure it is also installed.
@@ -105,4 +106,4 @@ Pađen, Ivan, García-Sánchez, Clara and Ledoux, Hugo (2022). Towards Automatic
 
 ## Acknowledgements
 We would like to acknowledge the authors of the supporting libraries we use in this project:
-[CGAL](https://github.com/CGAL/cgal), [CSF](https://github.com/jianboqi/CSF), [nlohmann/json](https://github.com/nlohmann/json), [LAStools](https://github.com/LAStools), [valijson](https://github.com/tristanpenman/valijson)
+[CGAL](https://github.com/CGAL/cgal), [CSF](https://github.com/jianboqi/CSF), [GDAL](https://github.com/OSGeo/gdal), [LAStools](https://github.com/LAStools), [nlohmann/json](https://github.com/nlohmann/json), [valijson](https://github.com/tristanpenman/valijson)

docker/city4cfd-build-base.dockerfile

@@ -13,4 +13,5 @@ RUN apt-get update && apt-get install -y --no-install-recommends \
     libboost-all-dev   \
     libcgal-dev        \
     libeigen3-dev      \
-    libomp-dev
+    libomp-dev         \
+    libgdal-dev

src/Building.cpp

@@ -40,20 +40,30 @@
 
 Building::Building()
         : PolyFeature(1), _elevation(-global::largnum), _height(-global::largnum),
-          _ptsPtr(std::make_shared<Point_set_3>()) {}
+          _ptsPtr(std::make_shared<Point_set_3>()), _hasFailed(false) {}
 
 Building::Building(const int internalID)
         : PolyFeature(1, internalID), _elevation(-global::largnum), _height(-global::largnum),
-          _ptsPtr(std::make_shared<Point_set_3>()) {}
+          _ptsPtr(std::make_shared<Point_set_3>()), _hasFailed(false) {}
 
 Building::Building(const nlohmann::json& poly)
         : PolyFeature(poly, true, 1), _elevation(-global::largnum), _height(-global::largnum),
-          _ptsPtr(std::make_shared<Point_set_3>()) {}
+          _ptsPtr(std::make_shared<Point_set_3>()), _hasFailed(false) {}
         // 'true' here to check for polygon simplicity
 
 Building::Building(const nlohmann::json& poly, const int internalID)
         : PolyFeature(poly, true, 1, internalID), _elevation(-global::largnum), _height(-global::largnum),
-          _ptsPtr(std::make_shared<Point_set_3>()) {}
+          _ptsPtr(std::make_shared<Point_set_3>()), _hasFailed(false) {}
+        // 'true' here to check for polygon simplicity
+
+Building::Building(const Polygon_with_attr& poly)
+        : PolyFeature(poly, true, 1), _elevation(-global::largnum), _height(-global::largnum),
+          _ptsPtr(std::make_shared<Point_set_3>()), _hasFailed(false) {}
+        // 'true' here to check for polygon simplicity
+
+Building::Building(const Polygon_with_attr& poly, const int internalID)
+        : PolyFeature(poly, true, 1, internalID), _elevation(-global::largnum), _height(-global::largnum),
+          _ptsPtr(std::make_shared<Point_set_3>()), _hasFailed(false) {}
         // 'true' here to check for polygon simplicity
 
 Building::~Building() = default;
@@ -77,6 +87,7 @@ void Building::alpha_wrap(const BuildingsPtr& buildings, Mesh& newMesh) {
     std::vector<std::array<FT, 3>> points;
     std::vector<CGAL_Polygon> polygons;
     for (auto& b : buildings) {
+        if (!b->is_active()) continue; // skip failed reconstructions
         auto& mesh = b->get_mesh();
         for (auto& face: mesh.faces()) {
             CGAL_Polygon p;
@@ -121,18 +132,18 @@ void Building::alpha_wrap(const BuildingsPtr& buildings, Mesh& newMesh) {
     */
 
     //-- Perform CGAL's alpha wrapping
-//    const double relative_alpha = 900.;
-//    const double relative_offset = 15000.;
-//    CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(newMesh);
-//    const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
-//                                         CGAL::square(bbox.ymax() - bbox.ymin()) +
-//                                         CGAL::square(bbox.zmax() - bbox.zmin()));
-//    const double alpha = diag_length / relative_alpha;
-//    const double offset = diag_length / relative_offset;
-//    CGAL::alpha_wrap_3(newMesh, alpha, offset, wrap);
+    const double relative_alpha = 2000.; //1000.
+    const double relative_offset = 7000.; // 12000.
+    CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(newMesh);
+    const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
+                                         CGAL::square(bbox.ymax() - bbox.ymin()) +
+                                         CGAL::square(bbox.zmax() - bbox.zmin()));
+    const double alpha = diag_length / relative_alpha;
+    const double offset = diag_length / relative_offset;
+    CGAL::alpha_wrap_3(newMesh, alpha, offset, newMesh);
 
 //    CGAL::alpha_wrap_3(newMesh, 2, 0.01, newMesh);   // 'coarse'
-    CGAL::alpha_wrap_3(newMesh, 1.5, 0.03, newMesh); // 'medium'
+//    CGAL::alpha_wrap_3(newMesh, 1.5, 0.03, newMesh); // 'medium'
 //    CGAL::alpha_wrap_3(newMesh, 0.7, 0.03, newMesh); // 'fine'
 
 //    CGAL::alpha_wrap_3(newMesh, 0.3, 0.03, newMesh); // that one takes long time
@@ -230,6 +241,14 @@ bool Building::has_self_intersections() const {
     return PMP::does_self_intersect(_mesh);
 }
 
+void Building::mark_as_failed() {
+    _hasFailed = true;
+}
+
+bool Building::has_failed_to_reconstruct() const {
+    return _hasFailed;
+}
+
 void Building::set_to_zero_terrain() {
     //-- Get average footprint height
     std::vector<double> avgRings;
@@ -296,4 +315,4 @@ TopoClass Building::get_class() const {
 
 std::string Building::get_class_name() const {
     return "Building";
-}
+}

src/Building.h

@@ -37,6 +37,8 @@ class Building : public PolyFeature {
     Building(const int internalID);
     Building(const nlohmann::json& poly);
     Building(const nlohmann::json& poly, const int internalID);
+    Building(const Polygon_with_attr& poly);
+    Building(const Polygon_with_attr& poly, const int internalID);
     ~Building();
 
     static void alpha_wrap(const BuildingsPtr& buildings, Mesh& newMesh);
@@ -52,6 +54,8 @@ class Building : public PolyFeature {
     void   translate_footprint(const double h);
     void   check_feature_scope(const Polygon_2& influRegion);
     void   set_clip_flag (const bool flag);
+    void   mark_as_failed();
+    bool   has_failed_to_reconstruct() const;
     bool   has_self_intersections() const;
     void   set_to_zero_terrain();
     double sq_max_dim();
@@ -66,6 +70,7 @@ class Building : public PolyFeature {
     PointSet3Ptr         _ptsPtr;
     double               _elevation;
     double               _height;
+    bool                 _hasFailed;
     bool                 _clip_bottom = Config::get().clip;
 };
 

src/CGALTypes.h

@@ -60,7 +60,7 @@ typedef   EPICK::Vector_2 Vector_2;
 typedef CGAL::Polygon_2<EPICK>              Polygon_2;
 typedef CGAL::Polygon_2<eProjection_traits> Polygon_3;
 
-//-- CGAL's Polygon_with_holes container expanded
+//- CGAL's Polygon_with_holes container expanded
 struct Polygon_with_holes_2 {
     std::vector<Polygon_2> _rings;
 
@@ -118,6 +118,14 @@ struct Polygon_with_holes_2 {
     CGAL::Bbox_2 bbox() const {return _rings.front().bbox();}
 };
 
+//- Polygon_with_holes expanded with attributes from GDAL
+struct Polygon_with_attr {
+    Polygon_with_holes_2 polygon;
+    std::unordered_map<std::string, std::string> attributes;
+//    std::string semanticClass;
+//    std::string id;
+};
+
 //-- Minimum bbox
 struct MinBbox {
     CGAL::Vector_2<EPICK> vec1;
@@ -216,8 +224,9 @@ typedef CGAL::Orthogonal_k_neighbor_search<Traits>                     Neighbor_
 typedef Neighbor_search::Tree                                          SearchTree;
 typedef CGAL::Fuzzy_iso_box<Traits>                                    Fuzzy_iso_box;
 
-//-- Smart pointers with CGAl types
-typedef std::shared_ptr<Point_set_3> PointSet3Ptr;
+//-- Smart pointers with CGAL types
+typedef std::shared_ptr<Point_set_3>                                  PointSet3Ptr;
+typedef std::vector<std::unique_ptr<Polygon_with_attr>>               PolyVecPtr;
 
 //-- Global constants
 namespace global {

src/Config.cpp

@@ -227,6 +227,9 @@ void Config::set_config(nlohmann::json& j) {
     if (j.contains("reconstruct_failed")) {
         reconstructFailed = j["reconstruct_failed"];
     }
+    if (j.contains("intersect_buildings_terrain")) {
+       intersectBuildingsTerrain = j["intersect_buildings_terrain"];
+    }
 
     // Imported buildings
     if (j.contains("import_geometries")) {

src/Config.h

@@ -84,13 +84,15 @@ class Config {
     int       nSmoothIterations                 = 0;
     int       maxSmoothPts                      = -9999;
     bool      flatTerrain                       = false;
+    bool      intersectBuildingsTerrain          = false;
     //- Buildings
     std::string buildingUniqueId;
     std::string lod;
     bool        refineReconstructedBuildings    = false;
     double      buildingPercentile              = -9999.; // Handled by schema
     double      minHeight                       = 2.;
     bool        reconstructFailed               = false;
+    std::string crsInfo;                                  // CRS information of building footprints
     // Height from attributes
     std::string buildingHeightAttribute;
     std::string floorAttribute;
@@ -124,7 +126,6 @@ class Config {
     std::string        logName                  = std::string("log");
     std::ostringstream log;
     std::ostringstream logSummary;
-    std::vector<int>   failedBuildings;
 
     //-- Experimental
     bool       clip                             = false;

src/ImportedBuilding.cpp

@@ -32,8 +32,10 @@
 #include "geomutils.h"
 #include "io.h"
 
-#include <CGAL/alpha_wrap_3.h>
+//#include <CGAL/alpha_wrap_3.h>
 #include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
 #include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
 
 int ImportedBuilding::noBottom = 0;
@@ -280,7 +282,7 @@ void ImportedBuilding::reconstruct() {
     nlohmann::json& geometry = (*_buildingJson)["geometry"][_lodIdx];
 
     _mesh.clear();
-    if (_clip_bottom) {
+    if (_clip_bottom || Config::get().intersectBuildingsTerrain) {
         this->translate_footprint(-5);
     }
     //-- Adjust building height points
@@ -353,21 +355,21 @@ void ImportedBuilding::reconstruct() {
     _mesh = wrap;
      */
 
-    if (_clip_bottom) {
+    if (_clip_bottom || Config::get().intersectBuildingsTerrain) {
         this->translate_footprint(5);
     }
         /*
         //-- Add other surfaces
         std::vector<Mesh::Vertex_index> faceVertices;
-        int facid = -1; //todo temp
+        int facid = -1; // temp
         for (auto& faceLst : faces) {
             if (++facid > 0) std::cout << "YO THERE's A FACE NO: " << facid << std::endl;
 //            if (!(facid > 0)) continue;
             for (auto& face: faceLst) {
                 faceVertices.emplace_back(_mesh.add_vertex(_ptsPtr[face]));
             }
             bool isReconstruct = _mesh.add_face(faceVertices);
-            //todo temp
+            // temp
             if (!isReconstruct) {
                 std::cout << "I have a failed surface!!!" << std::endl;
                 CGAL::Polygon_with_holes_2<EPICK> tempPoly;
@@ -463,4 +465,4 @@ void ImportedBuilding::set_footprint_mesh_connectivity(const std::unordered_map<
         }
         _footprintPtsIdxList.push_back(ringFootprintIdxs);
     }
-}
+}