Add mesh-derived attributes

apps/roofer-app/config.hpp

@@ -1,4 +1,4 @@
 // Copyright (c) 2018-2024 TU Delft 3D geoinformation group, Ravi Peters (3DGI),
 // and Balazs Dukai (3DGI)
 
 // This file is part of roofer (https://github.com/3DBAG/roofer)
@@ -139,7 +139,12 @@
       {"rmse_lod12", "rf_rmse_lod12"},
       {"rmse_lod13", "rf_rmse_lod13"},
       {"rmse_lod22", "rf_rmse_lod22"},
+      {"volume_lod12", "rf_volume_lod12"},
+      {"volume_lod13", "rf_volume_lod13"},
+      {"volume_lod22", "rf_volume_lod22"},
       {"h_ground", "rf_h_ground"},
+      {"slope", "rf_slope"},
+      {"azimuth", "rf_azimuth"},
   };
 };
 

apps/roofer-app/crop_tile.hpp

@@ -49,6 +49,11 @@ bool crop_tile(const roofer::TBox<double>& tile,
     return false;
   }
 
+  if (!pj->data_offset.has_value()) {
+    logger.error("No data offset set after reading inputs");
+    return false;
+  }
+
   // get yoc attribute vector (nullptr if it does not exist)
   auto yoc_vec = attributes.get_if<int>(cfg.yoc_attribute);
   if (!cfg.yoc_attribute.empty() && !yoc_vec) {
@@ -197,7 +202,7 @@ bool crop_tile(const roofer::TBox<double>& tile,
     pt_density.reserve(N_fp);
     is_glass_roof.reserve(N_fp);
     for (unsigned i = 0; i < N_fp; ++i) {
-      h_ground.push_back(ipc.ground_elevations[i]);
+      h_ground.push_back(ipc.ground_elevations[i] + (*pj->data_offset)[2]);
       nodata_r.push_back(ipc.nodata_radii[i]);
       nodata_frac.push_back(ipc.nodata_fractions[i]);
       pt_density.push_back(ipc.pt_densities[i]);
@@ -325,13 +330,13 @@ bool crop_tile(const roofer::TBox<double>& tile,
     {
       BuildingObject& building = output_building_tile.buildings.emplace_back();
       building.attribute_index = i;
+      building.z_offset = (*pj->data_offset)[2];
 
       building.pointcloud =
           input_pointclouds[selected->index].building_clouds[i];
       building.footprint = footprints[i];
       building.h_ground =
-          input_pointclouds[selected->index].ground_elevations[i] +
-          (*pj->data_offset)[2];
+          input_pointclouds[selected->index].ground_elevations[i];
       building.force_lod11 = input_pointclouds[selected->index].lod11_forced[i];
 
       output_building_tile.attributes = attributes;

apps/roofer-app/reconstruct_building.hpp

@@ -23,13 +23,15 @@ enum LOD { LOD12 = 12, LOD13 = 13, LOD22 = 22 };
 
 std::unordered_map<int, roofer::Mesh> extrude(
     roofer::Arrangement_2 arrangement, const float& floor_elevation,
+    const float& z_offset,
 #ifdef RF_USE_VAL3DITY
     std::string& attr_val3dity,
 #endif
-    float& rmse,
+    float& rmse, float& volume,
     roofer::reconstruction::SegmentRasteriserInterface* SegmentRasteriser,
     roofer::reconstruction::PlaneDetectorInterface* PlaneDetector, LOD lod,
-    roofer::ReconstructionConfig* cfg) {
+    RooferConfig* rfcfg) {
+  auto* cfg = &(rfcfg->rec);
   bool dissolve_step_edges = false;
   bool dissolve_all_interior = false;
   bool extrude_LoD2 = true;
@@ -78,9 +80,24 @@ std::unordered_map<int, roofer::Mesh> extrude(
               .with_class_ids(ArrangementExtruder->labels));
 #endif
 
+  auto MeshPropertyCalculator = roofer::misc::createMeshPropertyCalculator();
+  for (auto& [label, mesh] : ArrangementExtruder->multisolid) {
+    mesh.get_attributes().resize(mesh.get_polygons().size());
+    MeshPropertyCalculator->compute_roof_height(
+        mesh, {.z_offset = z_offset,
+               .cellsize = rfcfg->cellsize,
+               .h_50p = rfcfg->n["h_roof_50p"],
+               .h_70p = rfcfg->n["h_roof_70p"],
+               .h_min = rfcfg->n["h_roof_min"],
+               .h_max = rfcfg->n["h_roof_max"]});
+    MeshPropertyCalculator->compute_roof_orientation(
+        mesh, {.slope = rfcfg->n["slope"], .azimuth = rfcfg->n["azimuth"]});
+  }
+
   auto MeshTriangulator =
       roofer::reconstruction::createMeshTriangulatorLegacy();
   MeshTriangulator->compute(ArrangementExtruder->multisolid);
+  volume = MeshTriangulator->volumes.at(0);
   // logger.debug("Completed MeshTriangulator");
 #ifdef RF_USE_RERUN
   rec.log(worldname + "MeshTriangulator",
@@ -113,8 +130,8 @@ std::unordered_map<int, roofer::Mesh> extrude(
   return ArrangementExtruder->multisolid;
 }
 
-void reconstruct_building(BuildingObject& building,
-                          roofer::ReconstructionConfig* cfg) {
+void reconstruct_building(BuildingObject& building, RooferConfig* rfcfg) {
+  auto* cfg = &(rfcfg->rec);
   auto& logger = roofer::logger::Logger::get_logger();
   // split into ground and roof points
   roofer::PointCollection points, points_ground, points_roof;
@@ -314,32 +331,32 @@ void reconstruct_building(BuildingObject& building,
     // logger.debug("LoD={}", cfg->lod);
     if (cfg->lod == 0 || cfg->lod == 12) {
       building.multisolids_lod12 =
-          extrude(arrangement, building.h_ground,
+          extrude(arrangement, building.h_ground, building.z_offset,
 #ifdef RF_USE_VAL3DITY
                   building.val3dity_lod12,
 #endif
-                  building.rmse_lod12, SegmentRasteriser.get(),
-                  PlaneDetector.get(), LOD12, cfg);
+                  building.rmse_lod12, building.volume_lod12,
+                  SegmentRasteriser.get(), PlaneDetector.get(), LOD12, rfcfg);
     }
 
     if (cfg->lod == 0 || cfg->lod == 13) {
       building.multisolids_lod13 =
-          extrude(arrangement, building.h_ground,
+          extrude(arrangement, building.h_ground, building.z_offset,
 #ifdef RF_USE_VAL3DITY
                   building.val3dity_lod13,
 #endif
-                  building.rmse_lod13, SegmentRasteriser.get(),
-                  PlaneDetector.get(), LOD13, cfg);
+                  building.rmse_lod13, building.volume_lod13,
+                  SegmentRasteriser.get(), PlaneDetector.get(), LOD13, rfcfg);
     }
 
     if (cfg->lod == 0 || cfg->lod == 22) {
       building.multisolids_lod22 =
-          extrude(arrangement, building.h_ground,
+          extrude(arrangement, building.h_ground, building.z_offset,
 #ifdef RF_USE_VAL3DITY
                   building.val3dity_lod22,
 #endif
-                  building.rmse_lod22, SegmentRasteriser.get(),
-                  PlaneDetector.get(), LOD22, cfg);
+                  building.rmse_lod22, building.volume_lod22,
+                  SegmentRasteriser.get(), PlaneDetector.get(), LOD22, rfcfg);
     }
   }
 }

apps/roofer-app/roofer-app.cpp

@@ -56,6 +56,7 @@ namespace fs = std::filesystem;
 
 // reconstruct
 #include <roofer/misc/PC2MeshDistCalculator.hpp>
+#include <roofer/misc/MeshPropertyCalculator.hpp>
 #include <roofer/reconstruction/AlphaShaper.hpp>
 #include <roofer/reconstruction/ArrangementBuilder.hpp>
 #include <roofer/reconstruction/ArrangementDissolver.hpp>
@@ -93,6 +94,7 @@ namespace fs = std::filesystem;
 struct BuildingObject {
   roofer::PointCollection pointcloud;
   roofer::LinearRing footprint;
+  float z_offset = 0;
 
   std::unordered_map<int, roofer::Mesh> multisolids_lod12;
   std::unordered_map<int, roofer::Mesh> multisolids_lod13;
@@ -117,6 +119,9 @@ struct BuildingObject {
   float rmse_lod12;
   float rmse_lod13;
   float rmse_lod22;
+  float volume_lod12;
+  float volume_lod13;
+  float volume_lod22;
   std::string val3dity_lod12;
   std::string val3dity_lod13;
   std::string val3dity_lod22;
@@ -842,8 +847,7 @@ int main(int argc, const char* argv[]) {
               auto start = std::chrono::high_resolution_clock::now();
               logger.debug("[reconstructor] starting reconstruction for: {}",
                            building_object_ref.building.jsonl_path.string());
-              reconstruct_building(building_object_ref.building,
-                                   &roofer_cfg.rec);
+              reconstruct_building(building_object_ref.building, &roofer_cfg);
               logger.debug("[reconstructor] finished reconstruction for: {}",
                            building_object_ref.building.jsonl_path.string());
               // TODO: These two seem to be redundant
@@ -1060,6 +1064,8 @@ int main(int argc, const char* argv[]) {
               if (roofer_cfg.rec.lod == 0 || roofer_cfg.rec.lod == 12) {
                 ms12 = &building.multisolids_lod12;
                 attrow.insert(roofer_cfg.n["rmse_lod12"], building.rmse_lod12);
+                attrow.insert(roofer_cfg.n["volume_lod12"],
+                              building.volume_lod12);
 #if RF_USE_VAL3DITY
                 attrow.insert(roofer_cfg.n["val3dity_lod12"],
                               building.val3dity_lod12);
@@ -1068,6 +1074,8 @@ int main(int argc, const char* argv[]) {
               if (roofer_cfg.rec.lod == 0 || roofer_cfg.rec.lod == 13) {
                 ms13 = &building.multisolids_lod13;
                 attrow.insert(roofer_cfg.n["rmse_lod13"], building.rmse_lod13);
+                attrow.insert(roofer_cfg.n["volume_lod13"],
+                              building.volume_lod13);
 #if RF_USE_VAL3DITY
                 attrow.insert(roofer_cfg.n["val3dity_lod13"],
                               building.val3dity_lod13);
@@ -1076,6 +1084,8 @@ int main(int argc, const char* argv[]) {
               if (roofer_cfg.rec.lod == 0 || roofer_cfg.rec.lod == 22) {
                 ms22 = &building.multisolids_lod22;
                 attrow.insert(roofer_cfg.n["rmse_lod22"], building.rmse_lod22);
+                attrow.insert(roofer_cfg.n["volume_lod22"],
+                              building.volume_lod22);
 #if RF_USE_VAL3DITY
                 attrow.insert(roofer_cfg.n["val3dity_lod22"],
                               building.val3dity_lod22);

include/roofer/common/common.hpp

@@ -281,18 +281,12 @@ namespace roofer {
     float* get_data_ptr();
   };
 
-  // struct AttributeVec {
-  //   AttributeVec(std::type_index ttype) : value_type(ttype) {};
-  //   std::vector<std::any> values;
-  //   std::type_index value_type;
-  // };
-
   // class Mesh : public Geometry {
   // use indexed vertices?
   class Mesh {
     std::vector<LinearRing> polygons_;
     std::vector<int> labels_;
-    // std::unordered_map<std::string, AttributeVec>  attributes_;
+    std::vector<AttributeMapRow> attributes_;
 
    public:
     // Mesh() {};
@@ -307,9 +301,8 @@ namespace roofer {
     const std::vector<LinearRing>& get_polygons() const;
     std::vector<int>& get_labels();
     const std::vector<int>& get_labels() const;
-    // std::unordered_map<std::string, AttributeVec>&  get_attributes();
-    // const std::unordered_map<std::string, AttributeVec>&  get_attributes()
-    // const;
+    std::vector<AttributeMapRow>& get_attributes();
+    const std::vector<AttributeMapRow>& get_attributes() const;
   };
 
   struct Image {

include/roofer/io/CityJsonWriter.hpp

@@ -49,10 +49,8 @@ namespace roofer::io {
     size_t written_features_count = 0;
 
     bool prettyPrint_ = false;
-    bool only_output_renamed_ = false;
 
     vec1s key_options;
-    StrMap output_attribute_names;
 
     float translate_x_ = 0.;
     float translate_y_ = 0.;

include/roofer/misc/MeshPropertyCalculator.hpp

@@ -0,0 +1,54 @@
+// Copyright (c) 2018-2024 TU Delft 3D geoinformation group, Ravi Peters (3DGI),
+// and Balazs Dukai (3DGI)
+
+// This file is part of roofer (https://github.com/3DBAG/roofer)
+
+// geoflow-roofer was created as part of the 3DBAG project by the TU Delft 3D
+// geoinformation group (3d.bk.tudelf.nl) and 3DGI (3dgi.nl)
+
+// geoflow-roofer is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any
+// later version. geoflow-roofer is distributed in the hope that it will be
+// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+// Public License for more details. You should have received a copy of the GNU
+// General Public License along with geoflow-roofer. If not, see
+// <https://www.gnu.org/licenses/>.
+
+// Author(s):
+// Ravi Peters
+
+#pragma once
+#include <memory>
+#include <roofer/common/datastructures.hpp>
+
+namespace roofer::misc {
+
+  struct ComputeRoofHeightConfig {
+    float z_offset = 0;
+    float cellsize = 0.5;
+    std::string h_50p = "h_50p";
+    std::string h_70p = "h_70p";
+    std::string h_min = "h_min";
+    std::string h_max = "h_max";
+  };
+
+  struct ComputeRoofOrientationsConfig {
+    std::string slope = "slope";
+    std::string azimuth = "azimuth";
+    std::string roof_type = "roof_type";
+    float is_horizontal_threshold = 4;
+  };
+
+  struct MeshPropertyCalculatorInterface {
+    virtual ~MeshPropertyCalculatorInterface() = default;
+    virtual void compute_roof_height(Mesh& mesh,
+                                     ComputeRoofHeightConfig cfg) = 0;
+    virtual void compute_roof_orientation(
+        Mesh& mesh, ComputeRoofOrientationsConfig cfg) = 0;
+  };
+
+  std::unique_ptr<MeshPropertyCalculatorInterface>
+  createMeshPropertyCalculator();
+}  // namespace roofer::misc

include/roofer/reconstruction/MeshTriangulator.hpp

@@ -39,10 +39,6 @@ namespace roofer::reconstruction {
     vec3f normals;
     vec1i ring_ids;
     vec1f volumes;
-    // add_output("triangles_og", typeid(TriangleCollection));
-    // add_output("segment_ids_og", typeid(vec1i));
-    // add_output("triangles_snapped", typeid(TriangleCollection));
-    // add_output("segment_ids_snapped", typeid(vec1i));
 
     virtual ~MeshTriangulatorInterface() = default;
     virtual void compute(

src/core/common/common.cpp

@@ -422,13 +422,10 @@ namespace roofer {
 
   std::vector<int>& Mesh::get_labels() { return labels_; };
   const std::vector<int>& Mesh::get_labels() const { return labels_; };
-  // std::unordered_map<std::string, AttributeVec>&  Mesh::get_attributes(){
-  //   return attributes_;
-  // };
-  // const std::unordered_map<std::string, AttributeVec>&
-  // Mesh::get_attributes()const {
-  //   return attributes_;
-  // };
+  std::vector<AttributeMapRow>& Mesh::get_attributes() { return attributes_; };
+  const std::vector<AttributeMapRow>& Mesh::get_attributes() const {
+    return attributes_;
+  };
 
   void MultiTriangleCollection::push_back(
       TriangleCollection& trianglecollection) {