Merge branch 'release-0.2.0' into main

.dockerignore

@@ -3,3 +3,4 @@
 !/src
 !/thirdparty
 !/docker
+!/tools/prepare_point_cloud

.gitignore

@@ -8,5 +8,3 @@ Mesh*.json
 *.stl
 *.swp
 *.log
-thirdparty/LAStools/bin64/*
-!/thirdparty/LAStools/bin64/serf

CHANGELOG.md

@@ -1,7 +1,15 @@
 # Changelog
 
+## [0.2.0] - 2022-10-07
+### Added
+- Point cloud preparation tool city4cfd_pcprep
+### Changed
+- (breaking) Terrain smoothing overhaul
+### Removed
+- Previous point cloud preparation script
+
 ## [0.1.2] - 2022-09-13
-### Improved
+### Changed
 - Bad triangles handling
 
 ## [0.1.1] - 2022-09-02

CMakeLists.txt

@@ -2,7 +2,6 @@ cmake_minimum_required(VERSION 3.1)
 
 project(city4cfd)
 
-#set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake/Modules/" ${CMAKE_MODULE_PATH})
 #set( CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true )
 
 set(CMAKE_CXX_FLAGS "-O3")
@@ -13,21 +12,16 @@ set(CMAKE_BUILD_TYPE "Release")
 #    cmake_policy(SET CMP0003 NEW)
 #endif()
 
-#==============================
-# Alpha wrap support
-set(ALPHA_WRAP_SUPPORT OFF)
-#==============================
-
 # BOOST
-find_package(Boost 1.66 REQUIRED COMPONENTS filesystem)
+find_package(Boost 1.66 REQUIRED COMPONENTS filesystem locale)
 
 # CGAL
 find_package(CGAL REQUIRED QUIET COMPONENTS)
 if (CGAL_FOUND)
-    if (CGAL_VERSION VERSION_GREATER_EQUAL "5.3")
+    if (CGAL_VERSION VERSION_GREATER_EQUAL "5.5")
         message(STATUS "Found CGAL ${CGAL_VERSION}")
     elseif (CGAL_VERSION VERSION_GREATER_EQUAL "5.0")
-        message(STATUS "Found CGAL version greater than 5.0, but less than 5.3. "
+        message(STATUS "Found CGAL version greater than 5.0, but less than 5.5. "
                 "Proceeding to compile with included CGAL headers.")
         include_directories(${CMAKE_SOURCE_DIR}/thirdparty/CGAL/include)
         set(CGAL_USE_INCLUDED_HEADERS ON)
@@ -37,17 +31,21 @@ if (CGAL_FOUND)
         return()
     endif ()
 endif ()
-if (NOT CGAL_USE_INCLUDED_HEADERS AND ALPHA_WRAP_SUPPORT)
-    message(STATUS "Requested alpha wrap support. Compiling with included CGAL headers.")
-    include_directories( ${CMAKE_SOURCE_DIR}/thirdparty/CGAL/include )
-endif()
+
+# Eigen
+find_package(Eigen3 3.2.0) #(requires 3.2.0 or greater)
+include(${CMAKE_SOURCE_DIR}/thirdparty/CGAL/cmake/CGAL_Eigen3_support.cmake)
 
 # Third-party
 add_subdirectory(${CMAKE_SOURCE_DIR}/thirdparty/LAStools)
+add_subdirectory(${CMAKE_SOURCE_DIR}/thirdparty/CSF/src)
 include_directories(${CMAKE_SOURCE_DIR}/thirdparty)
 include_directories(${CMAKE_SOURCE_DIR}/thirdparty/valijson)
 include_directories(${Boost_INCLUDE_DIR})
 
+# Tools
+add_subdirectory(${CMAKE_SOURCE_DIR}/tools/prepare_point_cloud)
+
 # Creating entries for target: City4CFD
 FILE(GLOB SRC_FILES "src/*.cpp")
 add_executable(city4cfd ${SRC_FILES})
@@ -61,7 +59,9 @@ target_link_libraries(city4cfd
         ${CGAL_3RD_PARTY_LIBRARIES}
         ${Boost_SYSTEM_LIBRARY}
         ${Boost_FILESYSTEM_LIBRARY}
+        ${Boost_LOCALE_LIBRARY}
         LASlib
+        CGAL::Eigen3_support
         )
 
-install(TARGETS city4cfd DESTINATION bin)
+install(TARGETS city4cfd DESTINATION bin)

README.md

@@ -25,7 +25,7 @@ If you happen to use it, feedback is very much appreciated.
 City4CFD is developed by the [3D Geoinformation Research Group](https://3d.bk.tudelft.nl/) at the Delft University of Technology.
 
 ## Data formats
-**Point clouds** can be imported in LAS/LAZ, XYZ, or PLY format. We ask separately for ground and building points. While some datasets contain building-ground classification, some do not. In case your LAS/LAZ file has buildings/ground classified, you can use [point cloud preparation script](https://github.com/tudelft3d/City4CFD/blob/main/tools/prepare_point_cloud.sh) to create separate files. If buildings and terrain are under the same class, or vegetation is not filtered out, we suggest you use [CloudCompare](https://www.danielgm.net/cc/) to prepare points.
+**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/).
 
@@ -34,14 +34,27 @@ City4CFD is developed by the [3D Geoinformation Research Group](https://3d.bk.tu
 **Output** is in the following formats: OBJ, STL, and CityJSON. The ID of each polygon is preserved, and there is a 1-to-1 mapping between the input and the output.
 
 ## Installation
-You can directly compile City4CFD on your system using cmake, run it through a Docker container, or install using Homebrew in case of macOS.
+You can directly compile City4CFD on your system using cmake, run it through a Docker container, or install using Homebrew in the case of macOS.
 
 ### Build from source
 The following libraries are required to build the project:
 - [CGAL](https://www.cgal.org/) version 5
 - Boost >= 1.66
+- Eigen >= 3.2
+
+*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
+```
+
+In macOS you can install dependencies with Homebrew:
+
+```
+brew install cmake boost cgal eigen libomp
+```
 
-Both dependencies are generally available in Linux distributions (Debian/Ubuntu/Mint) as *libcgal-dev* and *libboost-dev*, and in macOS with Homebrew as *cgal* and *boost*.
 The project uses CMake to generate makefiles, so make sure it is also installed.
 
 To build City4CFD, do the following:
@@ -54,14 +67,13 @@ make
 You can speed up compilation by typing *make -j $numcores* where *$numcores* is the number of threads you can spare for compilation.
 
 ### Docker
-We offer built [Docker](https://www.docker.com/) images for every release, available at the [Docker Hub](https://hub.docker.com/r/tudelft3d/city4cfd). Running the docker script for the first time will pull the docker image from the Docker Hub repository.
+We offer built [Docker](https://www.docker.com/) images for every release, available at the [Docker Hub](https://hub.docker.com/r/tudelft3d/city4cfd). Running [the docker script](https://github.com/tudelft3d/City4CFD/tree/main/docker/run) for the first time will pull the docker image from the Docker Hub repository.
 
 ### macOS
 Mac users can install City4CFD through Homebrew:
 
 ```
-brew tap tudelft3d/software
-brew install city4cfd
+brew install tudelft3d/software/city4cfd
 ```
 
 ## Getting started
@@ -76,7 +88,7 @@ in case of building from a source.
 To run through a Docker container, you can use one of the scripts in ```docker/run/```. The script with the extension ```.sh``` can be used in Linux and macOS, the one with the extension ```.ps1``` in Windows Powershell, and the last one with ```.bat``` in Windows Command Prompt. You have to run  a script (you can copy it beforehand) from the root directory of the project (e.g. ```examples/TUD_Campus```), and the arguments are the same as for the compiled executable, e.g.: 
 
 ```
-../../docker/run/city4cfd_run.sh config_bpg.json --output_dir results
+../../docker/run/city4cfd_run.sh city4cfd config_bpg.json --output_dir results
 ```
 
 The script pulls the ```latest``` release from the Docker Hub. For a specific release, replace ```latest``` in the script with the released version tag, e.g. ```0.1.0```. In Linux systems, you will probably have to run the command as a sudo unless you create a 'docker' group and add users to it.
@@ -92,5 +104,5 @@ If you use City4CFD in a scientific context, please cite the following paper:
 Pađen, Ivan, García-Sánchez, Clara and Ledoux, Hugo (2022). Towards Automatic Reconstruction of 3D City Models Tailored for Urban Flow Simulations. *Frontiers in Built Environment*, 8, 2022 [[DOI](https://www.frontiersin.org/articles/10.3389/fbuil.2022.899332)][[BibTeX](https://github.com/tudelft3d/City4CFD/blob/master/CITATION.bib)]
 
 ## Acknowledgements
-We would like to acknowdledge the authors of supporting libraries we use in this project:
-[CGAL](https://github.com/CGAL/cgal), [nlohmann/json](https://github.com/nlohmann/json), [valijson](https://github.com/tristanpenman/valijson), [LAStools](https://github.com/LAStools)
+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)

docker/build-base.sh

@@ -1 +1 @@
-docker build -f city4cfd-build-base.dockerfile -t city4cfd-build-base --no-cache .
+docker build -f city4cfd-build-base.dockerfile -t tudelft3d/city4cfd:build-base --no-cache .

docker/city4cfd-build-base.dockerfile

@@ -11,4 +11,6 @@ RUN apt-get update && apt-get install -y --no-install-recommends \
     build-essential    \
     cmake              \
     libboost-all-dev   \
-    libcgal-dev  
+    libcgal-dev        \
+    libeigen3-dev      \
+    libomp-dev

docker/city4cfd.dockerfile

@@ -50,6 +50,8 @@ RUN bash /tmp/strip-docker-image-export \
     -f /bin/ls \
     -f /data \
     -f /usr/local/share/proj/proj.db \
+    -f /usr/local/bin/city4cfd_pcprep \
+    -f /usr/local/bin/city4cfd_las2las \
     -f /usr/local/bin/city4cfd
 
 #
@@ -70,5 +72,5 @@ COPY --from=builder /export/ /
 
 WORKDIR /data
 
-ENTRYPOINT ["city4cfd"]
-CMD ["--help"]
+#ENTRYPOINT ["city4cfd"]
+#CMD ["--help"]

examples/TUDCampus/config_bpg.json

@@ -30,16 +30,19 @@
   "bnd_type_bpg":      "Rectangle",
   "bpg_blockage_ratio": false,        
   "flow_direction":     [1, 1],
-  "buffer_region":     -20,
+  "buffer_region":      20,
   "reconstruct_boundaries": true,
 
-  "terrain_thinning": 80,
-  "smooth_terrain":   false,
+  "terrain_thinning": 10,
+  "smooth_terrain":   {
+    "iterations": 1,
+    "max_pts": 100000
+  },
 
   "lod": "1.2",
   "building_percentile": 90,
 
-  "edge_max_len": 2,
+  "edge_max_len": 5,
 
   "output_file_name":  "Mesh",
   "output_format":     "obj",

examples/TUDCampus/config_bpg_comments.json

@@ -45,15 +45,18 @@
 
   //-- Reconstruction related
   //- Terrain
-  "terrain_thinning": 90,             // Optional - percentage of randomly removed terrain points
-  "smooth_terrain":   true,           // Optional - smoothing flag
+  "terrain_thinning": 80,             // Optional - percentage of randomly removed terrain points
+  "smooth_terrain": {                 // Optional - smoothing flag
+    "iterations": 1,                  //            number of smoothing iterations. Typical numbers are 1-10
+    "max_pts": 100000                 //            number of points after optimized thinning for smoothing 
+  },  
   "flat_terrain":     false,          // Optional - make terrain flat
   //- Buildings
   "lod": "1.2",
   "building_percentile": 90,
 
   //-- Polygons related
-  "edge_max_len": 2,
+  "edge_max_len": 5,
 
   //-- Output
   "output_file_name":  "Mesh",

examples/TUDCampus/config_import_buildings.json

@@ -36,13 +36,16 @@
   "buffer_region":     5,
   "reconstruct_boundaries": false,
 
-  "terrain_thinning": 80,
-  "smooth_terrain":   false,
+  "terrain_thinning": 10,
+  "smooth_terrain":   {
+    "iterations": 1,
+    "max_pts": 100000
+  },
 
   "lod": "1.2",
   "building_percentile": 90,
 
-  "edge_max_len": 2,
+  "edge_max_len": 5,
 
   "output_file_name":  "Mesh",
   "output_format":     "obj",

examples/TUDCampus/config_manual.json

@@ -30,13 +30,16 @@
   "buffer_region":     5,
   "reconstruct_boundaries": true,
 
-  "terrain_thinning": 80,
-  "smooth_terrain":   false,
+  "terrain_thinning": 10,
+  "smooth_terrain":   {
+    "iterations": 1,
+    "max_pts": 100000
+  },
 
   "lod": "1.2",
   "building_percentile": 90,
 
-  "edge_max_len": 2,
+  "edge_max_len": 5,
 
   "output_file_name":  "Mesh",
   "output_format":     "obj",

src/Building.cpp

@@ -36,6 +36,7 @@
 #include <CGAL/Polygon_mesh_processing/clip.h>
 #include <CGAL/Polygon_mesh_processing/remesh.h>
 #include <CGAL/Mesh_3/dihedral_angle_3.h>
+#include <CGAL/alpha_wrap_3.h>
 
 Building::Building()
         : PolyFeature(1), _height(-global::largnum) {}
@@ -53,6 +54,78 @@ Building::Building(const nlohmann::json& poly, const int internalID)
 
 Building::~Building() = default;
 
+void Building::alpha_wrap(const Buildings& buildings, Mesh& newMesh) {
+    typedef EPICK::FT                 FT;
+    typedef std::array<FT, 3>         Custom_point;
+    typedef std::vector<std::size_t>  CGAL_Polygon;
+
+    //-- Make a single mesh out of all individual buildings
+    std::vector<std::array<FT, 3>> points;
+    std::vector<CGAL_Polygon> polygons;
+    for (auto& b : buildings) {
+        auto& mesh = b->get_mesh();
+        for (auto& face: mesh.faces()) {
+            CGAL_Polygon p;
+            auto vertices = mesh.vertices_around_face(mesh.halfedge(face));
+            for (auto vertex = vertices.begin(); vertex != vertices.end(); ++vertex) {
+                points.push_back(CGAL::make_array<FT>(mesh.point(*vertex).x(),
+                                                      mesh.point(*vertex).y(),
+                                                      mesh.point(*vertex).z()));
+                p.push_back(points.size() - 1);
+            }
+            polygons.push_back(p);
+        }
+    }
+    PMP::repair_polygon_soup(points, polygons, CGAL::parameters::geom_traits(geomutils::Array_traits()));
+    PMP::orient_polygon_soup(points, polygons);
+    PMP::polygon_soup_to_polygon_mesh(points, polygons, newMesh);
+    PMP::triangulate_faces(newMesh);
+
+    /*
+    typedef Mesh::Halfedge_index           halfedge_descriptor;
+    typedef Mesh::Edge_index               edge_descriptor;
+    //-- Set the property map for constrained edges
+    Mesh::Property_map<edge_descriptor,bool> is_constrained =
+            newMesh.add_property_map<edge_descriptor,bool>("e:is_constrained",false).first;
+
+    //-- Detect sharp features
+    for (auto& e : edges(newMesh)) {
+        halfedge_descriptor hd = halfedge(e,newMesh);
+        if (!is_border(e,newMesh)) {
+            double angle = CGAL::Mesh_3::dihedral_angle(newMesh.point(source(hd,newMesh)),
+                                                        newMesh.point(target(hd,newMesh)),
+                                                        newMesh.point(target(next(hd,newMesh),newMesh)),
+                                                        newMesh.point(target(next(opposite(hd,newMesh),newMesh),newMesh)));
+            if (CGAL::abs(angle)<179.5)
+                is_constrained[e]=true;
+        }
+    }
+
+    Mesh wrap;
+    CGAL::alpha_wrap_3(newMesh, 0.1, 0.01, wrap,
+                             CGAL::parameters::edge_is_constrained_map(is_constrained));
+    */
+
+    //-- 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);
+
+//    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, 0.7, 0.03, newMesh); // 'fine'
+
+//    CGAL::alpha_wrap_3(newMesh, 0.3, 0.03, newMesh); // that one takes long time
+//    CGAL::alpha_wrap_3(points, polygons, 0.1, 0.001, newMesh);
+//    newMesh = wrap;
+}
+
 void Building::clip_bottom(const Terrainptr& terrain) {
     if (!_clip_bottom) return;
     if (this->has_self_intersections() && !Config::get().handleSelfIntersect) throw
@@ -81,16 +154,16 @@ void Building::refine() {
     typedef Mesh::Halfedge_index           halfedge_descriptor;
     typedef Mesh::Edge_index               edge_descriptor;
 
-    double target_edge_length = 5; //5;
-    unsigned int nb_iter =  30;   //30;
+    const double target_edge_length = 5; //5;
+    const unsigned int nb_iter =  30;   //30;
 
     PMP::remove_degenerate_faces(_mesh);
     /*
     if (PMP::does_self_intersect(_mesh)) {
         ++config::selfIntersecting;
         PMP::remove_self_intersections(_mesh);
     }
-     */
+    */
 
     //-- Set the property map for constrained edges
     Mesh::Property_map<edge_descriptor,bool> is_constrained =