Merge branch 'release-0.4.6'

CHANGELOG.md

@@ -1,5 +1,11 @@
 # Changelog
 
+## [0.4.6] - 2024-03-11
+### Changed
+- Improved robustness in handling bad polygons
+### Fixed
+- Code crashing when flattening specific cases
+
 ## [0.4.5] - 2024-02-12
 ### Fixed
 - Dependency hotfix

CMakeLists.txt

@@ -4,7 +4,7 @@ project(city4cfd)
 
 #set( CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true )
 
-set(CMAKE_CXX_FLAGS "-O2")
+set(CMAKE_CXX_FLAGS "-O3")
 set(CMAKE_BUILD_TYPE "Release")
 
 if (COMMAND cmake_policy)

src/BoundingRegion.cpp

@@ -116,25 +116,27 @@ void BoundingRegion::calc_bnd_bpg(const Polygon_2& influRegionPoly,
     Config::get().topHeight = elevationMax + hMax * (Config::get().bpgDomainSize.back() - 1);
 
     //-- Blockage ratio handling
-    std::cout << "\nCalculating blockage ratio for flow direction (" << Config::get().flowDirection
-              << ")" << std::endl;
+    if (Config::get().bpgBlockageRatioFlag) {
+        std::cout << "\nCalculating blockage ratio for flow direction (" << Config::get().flowDirection
+                  << ")" << std::endl;
 
 //    double blockRatio = this->calc_blockage_ratio_from_chull(buildings, angle, localPoly);
 //    double blockRatio = this->calc_blockage_ratio_from_ashape(buildings, angle, localPoly);
-    double blockRatio = this->calc_blockage_ratio_comb(buildings, angle, localPoly);
+        double blockRatio = this->calc_blockage_ratio_comb(buildings, angle, localPoly);
 //    double blockRatio = this->calc_blockage_ratio_from_ashape_alt(buildings, angle, localPoly);
 //    double blockRatio = this->calc_blockage_ratio_from_edges(buildings, angle, localPoly);
 
-    std::cout << "    Blockage ratio is: " << blockRatio << std::endl;
-    if (Config::get().bpgBlockageRatioFlag && blockRatio > Config::get().bpgBlockageRatio) {
-        std::cout << "INFO: Blockage ratio is more than " << Config::get().bpgBlockageRatio * 100
-                  << "%. Expanding domain cross section to meet the guideline"
-                  << std::endl;
-        double expRatio = std::sqrt(blockRatio / 0.03);
-        //-- Recalculate the bnd poly and height with new values
-        localPoly.clear();
-        localPoly = this->calc_bnd_poly(candidatePts, hMax, angle, expRatio);
-        Config::get().topHeight = hMax * Config::get().bpgDomainSize.back() * expRatio;
+        std::cout << "    Blockage ratio is: " << blockRatio << std::endl;
+        if (blockRatio > Config::get().bpgBlockageRatio) {
+            std::cout << "INFO: Blockage ratio is more than " << Config::get().bpgBlockageRatio * 100
+                      << "%. Expanding domain cross section to meet the guideline"
+                      << std::endl;
+            double expRatio = std::sqrt(blockRatio / 0.03);
+            //-- Recalculate the bnd poly and height with new values
+            localPoly.clear();
+            localPoly = this->calc_bnd_poly(candidatePts, hMax, angle, expRatio);
+            Config::get().topHeight = hMax * Config::get().bpgDomainSize.back() * expRatio;
+        }
     }
     //-- Return the points back to global coordinates
     for (auto& pt : localPoly) _boundingRegion.push_back(geomutils::rotate_pt(pt, angle));

src/Map3d.cpp

@@ -112,7 +112,7 @@ void Map3d::set_features() {
         _buildingsPtr.push_back(building);
         _allFeaturesPtr.push_back(building);
     }
-    if (Config::get().avoidBadPolys) this->clear_inactives(); // Remove buildings that potentially couldn't be imported
+    this->clear_inactives(); // Remove buildings that potentially couldn't be imported
     //- Imported buildings
     if (!_importedBuildingsJSON.empty()) {
         std::cout << "Importing CityJSON geometries" << std::endl;
@@ -170,6 +170,7 @@ void Map3d::set_features() {
             _allFeaturesPtr.push_back(surfacePoly);
         }
     }
+    this->clear_inactives();
     std::cout << "Polygons read: " << _allFeaturesPtr.size() << std::endl;
 
     //-- Set flat terrain or random thin terrain points
@@ -349,15 +350,13 @@ void Map3d::reconstruct_buildings() {
         auto& b = _buildingsPtr[i];
         if (b->is_active()) this->reconstruct_one_building(b);
     }
-    this->clear_inactives(); // in case of imported-reconstructed fallback
+    this->clear_inactives(); // renumber failed and in case of imported-reconstructed fallback
     // Gather failed reconstructions
-    int failed = 0;
-    for (auto&  b : _buildingsPtr) if (b->has_failed_to_reconstruct()) ++failed;
-    std::cout << "    Number of successfully reconstructed buildings: " << _buildingsPtr.size() - failed << std::endl;
+    std::cout << "    Number of successfully reconstructed buildings: " << _buildingsPtr.size() << std::endl;
     Config::get().logSummary << "Building reconstruction summary: successfully reconstructed buildings: "
-                             << _buildingsPtr.size() - failed << std::endl;
+                             << _buildingsPtr.size() << std::endl;
     Config::get().logSummary << "                                 num of failed reconstructions: "
-                             << failed << std::endl;
+                             << _failedBuildingsPtr.size() << std::endl;
 }
 
 void Map3d::reconstruct_one_building(std::shared_ptr<Building>& building) {
@@ -580,7 +579,7 @@ void Map3d::prep_cityjson_output() { // Temp impl, might change
 };
 
 void Map3d::clear_inactives() {
-    for (unsigned long i = 0; i < _reconstructedBuildingsPtr.size();) {
+    for (int i = 0; i < _reconstructedBuildingsPtr.size();) {
         if (_reconstructedBuildingsPtr[i]->is_active()) ++i;
         else {
             _reconstructedBuildingsPtr.erase(_reconstructedBuildingsPtr.begin() + i);
@@ -596,7 +595,7 @@ void Map3d::clear_inactives() {
                     inactiveBuildingIdxs.push_back(importedBuilding->get_id());
             }
         }
-        for (unsigned long i = 0; i < _importedBuildingsPtr.size();) {
+        for (int i = 0; i < _importedBuildingsPtr.size();) {
             auto it = std::find(inactiveBuildingIdxs.begin(), inactiveBuildingIdxs.end(),
                                 _importedBuildingsPtr[i]->get_id());
             if (it == inactiveBuildingIdxs.end()) ++i;
@@ -606,26 +605,27 @@ void Map3d::clear_inactives() {
             }
         }
     } else {
-        for (unsigned long i = 0; i < _importedBuildingsPtr.size();) {
+        for (int i = 0; i < _importedBuildingsPtr.size();) {
             if (_importedBuildingsPtr[i]->is_active()) ++i;
             else {
                 _importedBuildingsPtr.erase(_importedBuildingsPtr.begin() + i);
             }
         }
     }
-    for (unsigned long i = 0; i < _buildingsPtr.size();) {
-        if (_buildingsPtr[i]->is_active() || _buildingsPtr[i]->has_failed_to_reconstruct()) ++i;
+    for (int i = 0; i < _buildingsPtr.size();) {
+        if (_buildingsPtr[i]->is_active()) ++i;
         else {
+            if (_buildingsPtr[i]->has_failed_to_reconstruct()) _failedBuildingsPtr.push_back(_buildingsPtr[i]);
             _buildingsPtr.erase(_buildingsPtr.begin() + i);
         }
     }
-    for (unsigned long i = 0; i < _surfaceLayersPtr.size();) {
+    for (int i = 0; i < _surfaceLayersPtr.size();) {
         if (_surfaceLayersPtr[i]->is_active()) ++i;
         else {
             _surfaceLayersPtr.erase(_surfaceLayersPtr.begin() + i);
         }
     }
-    for (unsigned long i = 0; i < _allFeaturesPtr.size();) {
+    for (int i = 0; i < _allFeaturesPtr.size();) {
         if (_allFeaturesPtr[i]->is_active()) ++i;
         else {
             _allFeaturesPtr.erase(_allFeaturesPtr.begin() + i);
@@ -634,11 +634,7 @@ void Map3d::clear_inactives() {
 }
 
 BuildingsPtr Map3d::get_failed_buildings() const {
-    BuildingsPtr failedBuildings;
-    for (auto& b : _buildingsPtr) {
-        if (b->has_failed_to_reconstruct()) failedBuildings.push_back(b);
-    }
-    return failedBuildings;
+    return _failedBuildingsPtr;
 }
 
 //-- Templated functions

src/Map3d.h

@@ -57,6 +57,7 @@ class Map3d {
 
     TerrainPtr                  _terrainPtr;
     BuildingsPtr                _buildingsPtr;
+    BuildingsPtr                _failedBuildingsPtr;
     ReconstructedBuildingsPtr   _reconstructedBuildingsPtr;
     ImportedBuildingsPtr        _importedBuildingsPtr;
     SurfaceLayersPtr            _surfaceLayersPtr;

src/PolyFeature.cpp

@@ -89,12 +89,6 @@ PolyFeature::PolyFeature(const Polygon_with_attr& poly, const bool checkSimplici
             this->deactivate();
             return;
         }
-        if (isOuterRing) {
-            if (tempPoly.is_clockwise_oriented()) tempPoly.reverse_orientation();
-            isOuterRing = false;
-        } else {
-            if (tempPoly.is_counterclockwise_oriented()) tempPoly.reverse_orientation();
-        }
         if (checkSimplicity) {
             if (!tempPoly.is_simple()) {
                 if (Config::get().avoidBadPolys) {
@@ -109,6 +103,12 @@ PolyFeature::PolyFeature(const Polygon_with_attr& poly, const bool checkSimplici
                 }
             }
         }
+        if (isOuterRing) {
+            if (tempPoly.is_clockwise_oriented()) tempPoly.reverse_orientation();
+            isOuterRing = false;
+        } else {
+            if (tempPoly.is_counterclockwise_oriented()) tempPoly.reverse_orientation();
+        }
         _poly._rings.push_back(tempPoly);
     }
 }
@@ -222,6 +222,15 @@ bool PolyFeature::flatten_polygon_inner_points(const Point_set_3& pointCloud,
     typedef boost::shared_ptr<CGAL::Polygon_with_holes_2<EPICK>> PolygonPtrWH;
     typedef std::vector<PolygonPtrWH> PolygonPtrVectorWH;
 
+#ifdef CITY4CFD_POLYFEATURE_VERBOSE
+    std::cout << "\nINFO: Flattening a polygon..." << std::endl;
+    int nonSimpleRings = 0;
+    for (auto& ring : _poly.rings()) {
+        if (!ring.is_simple()) ++nonSimpleRings;
+    }
+    std::cout << "Non-simple rings: " << nonSimpleRings << std::endl;
+#endif
+
     std::vector<int>    indices;
     std::vector<double> originalHeights;
     auto building_pt = pointCloud.property_map<std::shared_ptr<Building>>("building_point").first;
@@ -237,9 +246,7 @@ bool PolyFeature::flatten_polygon_inner_points(const Point_set_3& pointCloud,
     std::map<int, std::shared_ptr<Building>> overlappingBuildings; //id-building map
     for (auto& pt3 : subsetPts) {
         Point_2 pt(pt3.x(), pt3.y());
-        if (CGAL::bounded_side_2(_poly._rings.front().begin(),
-                                 _poly._rings.front().end(),
-                                 pt) != CGAL::ON_UNBOUNDED_SIDE) {
+        if (geomutils::point_in_poly_and_boundary(pt, _poly)) {
             auto itIdx = pointCloudConnectivity.find(pt3);
             auto pointSetIt = pointCloud.begin();
             std::advance(pointSetIt, itIdx->second);
@@ -258,37 +265,46 @@ bool PolyFeature::flatten_polygon_inner_points(const Point_set_3& pointCloud,
     //-- If next intersecting a building, clip poly with building
     std::vector<Polygon_with_holes_2> flattenCandidatePolys;
     if (!overlappingBuildings.empty()) {
+#ifdef CITY4CFD_POLYFEATURE_VERBOSE
+        std::cout << "Overlapping buildings: " << overlappingBuildings.size() << std::endl;
+#endif
         // Add clipping polys to set
         CGAL::Polygon_set_2<EPECK> polySet;
         Converter<EPICK, EPECK> to_exact;
         for (auto& intersectBuilding : overlappingBuildings) {
-            polySet.insert(intersectBuilding.second->get_poly().get_exact_outer_boundary());
+            polySet.join(intersectBuilding.second->get_poly().get_exact_outer_boundary());
         }
+#ifdef CITY4CFD_POLYFEATURE_VERBOSE
+        std::cout << "Polyset size? " << polySet.number_of_polygons_with_holes() << std::endl;
+#endif
+        CGAL_assertion(polySet.is_valid());
         // Clip with this polygon
         polySet.complement();
-        polySet.intersection(this->get_poly().get_exact());
+        polySet.intersection(_poly.get_exact());
         // Store this as the new polygon
         std::vector<CGAL::Polygon_with_holes_2<EPECK>> resPolys;
         polySet.polygons_with_holes(std::back_inserter(resPolys));
-#ifdef CITY4CFD_POLYFEATURE_VERBOSE
-        std::cout << "After flatten polygon clipping, total new polygons: " << resPolys.size() << std::endl;
-#endif
         for (auto& flattenBndPoly : resPolys) {
             flattenCandidatePolys.emplace_back(geomutils::exact_poly_to_poly(flattenBndPoly));
         }
     }
-    std::vector<Polygon_2> flattenBndPolys;
+#ifdef CITY4CFD_POLYFEATURE_VERBOSE
+    std::cout << "After flatten polygon clipping, total new polygons: " << flattenCandidatePolys.size() << std::endl;
+#endif
     if (!flattenCandidatePolys.empty()) {
         bool isFirst = true;
         for (auto& poly : flattenCandidatePolys) {
             const double offsetVal = 0.1; // offset value hardcoded
             PolygonPtrVectorWH offset_poly =
                     CGAL::create_interior_skeleton_and_offset_polygons_with_holes_2(offsetVal,
                                                                                     poly.get_cgal_type());
-            if (offset_poly.size() == 1) { // make a check whether the offset is successfully created or not
-                flattenBndPolys.push_back(offset_poly.front()->outer_boundary());
+#ifdef CITY4CFD_POLYFEATURE_VERBOSE
+            std::cout << "Offset polygons created: " << offset_poly.size() << std::endl;
+#endif
+            if (!offset_poly.empty()) { // make a check whether the offset is successfully created
+                poly = *(offset_poly.front());
             } else {
-                std::cout << "Skeleton construction failed!" << std::endl;
+                std::cout << "WARNING: Skeleton construction failed! Some surfaces might not be flattened." << std::endl;
                 return false;
             }
             if (isFirst) {
@@ -300,15 +316,18 @@ bool PolyFeature::flatten_polygon_inner_points(const Point_set_3& pointCloud,
             }
         }
     } else {
-        flattenBndPolys.push_back(_poly.outer_boundary());
+        if (!overlappingBuildings.empty()) {
+            std::cout << "WARNING: Polygon ID" << _id << " flattening failed due to unresolved"
+                                                         " intersections with buildings!" << std::endl;
+            return false;
+        }
+        flattenCandidatePolys.push_back(_poly);
     }
     //-- Collect points that have not been already flattened
     for (auto& pt3 : subsetPts) {
-        for (auto& flattenBndPoly : flattenBndPolys) {
+        for (auto& flattenBndPoly : flattenCandidatePolys) {
             Point_2 pt(pt3.x(), pt3.y());
-            if (CGAL::bounded_side_2(flattenBndPoly.begin(),
-                                     flattenBndPoly.end(),
-                                     pt) != CGAL::ON_UNBOUNDED_SIDE) {
+            if (geomutils::point_in_poly_and_boundary(pt, flattenBndPoly)) {
                 auto itIdx = pointCloudConnectivity.find(pt3);
 
                 auto it = flattenedPts.find(itIdx->second);
@@ -331,16 +350,17 @@ bool PolyFeature::flatten_polygon_inner_points(const Point_set_3& pointCloud,
         flattenedPts[i] = Point_3(pointCloud.point(i).x(), pointCloud.point(i).y(), avgHeight);
     }
     //-- Add additional new segments to constrain
-    if (!flattenBndPolys.empty()) {
-        for (auto& flattenBndPoly : flattenBndPolys) {
-            for (const auto& newEdge: flattenBndPoly.edges()) {
-                ePoint_3 pt1(newEdge.point(0).x(), newEdge.point(0).y(), avgHeight);
-                ePoint_3 pt2(newEdge.point(1).x(), newEdge.point(1).y(), avgHeight);
-                constrainedEdges.emplace_back(pt1, pt2);
+    if (!flattenCandidatePolys.empty()) {
+        for (auto& flattenBndPoly : flattenCandidatePolys) {
+            for (auto& ring: flattenBndPoly.rings()) {
+                for (const auto& newEdge: ring.edges()) {
+                    ePoint_3 pt1(newEdge.point(0).x(), newEdge.point(0).y(), avgHeight);
+                    ePoint_3 pt2(newEdge.point(1).x(), newEdge.point(1).y(), avgHeight);
+                    constrainedEdges.emplace_back(pt1, pt2);
+                }
             }
         }
     }
-
     return true;
 }
 
@@ -425,11 +445,10 @@ void PolyFeature::parse_json_poly(const nlohmann::json& poly, const bool checkSi
             }
         }
         geomutils::pop_back_if_equal_to_front(tempPoly);
-
-        if (_poly._rings.empty()) {
-            if (tempPoly.is_clockwise_oriented()) tempPoly.reverse_orientation();
-        } else {
-            if (tempPoly.is_counterclockwise_oriented()) tempPoly.reverse_orientation();
+        if (tempPoly.size() < 3) { // Sanity check if it is even a polygon
+            std::cout << "WARNING: Skipping import of a zero-area polygon" << std::endl;
+            this->deactivate();
+            return;
         }
         if (checkSimplicity) {
             if (!tempPoly.is_simple()) {
@@ -445,6 +464,11 @@ void PolyFeature::parse_json_poly(const nlohmann::json& poly, const bool checkSi
                 }
             }
         }
+        if (_poly._rings.empty()) {
+            if (tempPoly.is_clockwise_oriented()) tempPoly.reverse_orientation();
+        } else {
+            if (tempPoly.is_counterclockwise_oriented()) tempPoly.reverse_orientation();
+        }
         _poly._rings.push_back(tempPoly);
     }
 }

src/Terrain.cpp

@@ -102,7 +102,7 @@ void Terrain::constrain_features() {
     }
     IO::print_progress_bar(100); std::clog << std::endl;
     // extra edges to constrain when whole polygons couldn't be added
-    if (!_extraConstrainedEdges.empty()) std::cout << "\n    Adding extra constrained edges" << std::endl;
+    if (!_extraConstrainedEdges.empty()) std::cout << "\n    Inserting additional constrained edges" << std::endl;
     for (auto& extraEdge : _extraConstrainedEdges) {
         _cdt.insert_constraint(extraEdge.source(), extraEdge.target());
         ++count;