add option to output basic terrain

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 - Repair of non-manifold and self-intersecting roof junctions during arrangement snapping. This prevents certain geometric errors like non-watertight meshes, and highly non-planar faces that could prevriously happen in complex roof arrangements. See #168.
+- Optional per-tile triangulated terrain output as CityJSON `TINRelief` features, see the `--terrain` option.
 
 ## Changed
 - The pointcloud-insufficient test is now an absolute per-building density floor (`min-building-density`, points/m²) instead of a tile-relative `mean - 2·std` outlier test. The old test depended the entire tile that was being processed, so the same building could be flagged differently between runs with different tiling — potentially skipping point-rich buildings and leaving them without geometry. The decision is now deterministic and depends only on the building's own data.

apps/roofer-app/config.hpp

@@ -25,7 +25,7 @@
 #include <concepts>
 #include <thread>
 #include <algorithm>
 #include "toml.hpp"
 #include "fmt/format.h"
 #include "fmt/ranges.h"
 
@@ -86,6 +86,27 @@
   std::vector<fileExtent> file_extents;
 };
 
+inline std::optional<size_t> select_terrain_pointcloud(
+    const std::vector<InputPointcloud>& input_pointclouds) {
+  std::optional<size_t> selected;
+  for (size_t i = 0; i < input_pointclouds.size(); ++i) {
+    const auto& candidate = input_pointclouds[i];
+    if (candidate.select_only_for_date) continue;
+    if (!selected.has_value()) {
+      selected = i;
+      continue;
+    }
+
+    const auto& current = input_pointclouds[*selected];
+    if (candidate.quality < current.quality ||
+        (candidate.quality == current.quality &&
+         candidate.date > current.date)) {
+      selected = i;
+    }
+  }
+  return selected;
+}
+
 struct RooferConfigHandler;
 
 struct RooferConfig {
@@ -111,6 +132,7 @@
   float min_building_density = 1.0;
   float terrain_grid_cellsize = 10.0;
   int terrain_grid_search_radius = 3;
+  std::string terrain_nodata_mode = "fill_small_gaps";
   int lod11_fallback_area = 69000;
   float lod11_fallback_density = 5;
   roofer::arr2f tilesize = {1000, 1000};
@@ -133,6 +155,7 @@
   // crop output
   bool split_cjseq = false;
   bool omit_metadata = false;
+  bool output_terrain = false;
   roofer::arr3d cj_scale = {0.001, 0.001, 0.001};
   std::optional<roofer::arr3d> cj_translate;
   std::string building_toml_file_spec =
@@ -440,6 +463,13 @@
              "Number of terrain grid cells to search around a building when "
              "its local fallback cells do not contain terrain points.",
              cfg_.terrain_grid_search_radius, {check::HigherOrEqualTo<int>(0)});
+    crop.add("terrain-nodata-mode",
+             "How missing terrain grid samples are handled during output "
+             "triangulation: `complete_quads`, `local_triangles`, or "
+             "`fill_small_gaps`.",
+             cfg_.terrain_nodata_mode,
+             {check::OneOf<std::string>(
+                 {"complete_quads", "local_triangles", "fill_small_gaps"})});
     crop.add(
         "lod11-fallback-area",
         "LoD 1.1 fallback threshold area in square meters. If the area of the "
@@ -559,6 +589,10 @@
     output.add("omit-metadata",
                "Omit metadata line in CityJSONSequence output.",
                cfg_.omit_metadata);
+    output.add("terrain",
+               "Write one triangulated TINRelief feature per tile from the "
+               "highest-quality pointcloud source.",
+               cfg_.output_terrain);
     output.add("cj-scale", "Scaling applied to CityJSON output vertices",
                cfg_.cj_scale);
     output

apps/roofer-app/crop_tile.hpp

@@ -19,6 +19,7 @@
 // Author(s):
 // Ravi Peters
 #pragma once
+
 bool crop_tile(const roofer::TBox<double>& tile,
                std::vector<InputPointcloud>& input_pointclouds,
                BuildingTile& output_building_tile, const RooferConfig& cfg,
@@ -33,6 +34,14 @@ bool crop_tile(const roofer::TBox<double>& tile,
   auto vector_ops = roofer::misc::createVector2DOpsGEOS(*pj);
   auto LASWriter = roofer::io::createLASWriter(*pj);
 
+  const auto terrain_pointcloud =
+      cfg.output_terrain ? select_terrain_pointcloud(input_pointclouds)
+                         : std::nullopt;
+  if (cfg.output_terrain && !terrain_pointcloud.has_value()) {
+    logger.warning(
+        "Terrain output requested, but no eligible pointcloud source exists");
+  }
+
   // logger.info("region_of_interest.has_value()? {}",
   // region_of_interest.has_value()); if(region_of_interest.has_value())
   vector_reader->layer_name = cfg.layer_name;
@@ -83,7 +92,9 @@ bool crop_tile(const roofer::TBox<double>& tile,
       polygon_extent.pmax[0], polygon_extent.pmax[1], polygon_extent.pmax[2]));
 
   // Crop all pointclouds
-  for (auto& ipc : input_pointclouds) {
+  for (size_t ipc_index = 0; ipc_index < input_pointclouds.size();
+       ++ipc_index) {
+    auto& ipc = input_pointclouds[ipc_index];
     logger.info("Cropping pointcloud {}...", ipc.name);
 
     auto intersecting_files = ipc.rtree->query(polygon_extent_untransformed);
@@ -103,8 +114,45 @@ bool crop_tile(const roofer::TBox<double>& tile,
          .building_class = ipc.bld_class,
          .use_acquisition_year = static_cast<bool>(yoc_vec),
          .terrain_grid_cellsize = cfg.terrain_grid_cellsize,
-         .terrain_grid_search_radius = cfg.terrain_grid_search_radius});
+         .terrain_grid_search_radius = cfg.terrain_grid_search_radius,
+         .retain_terrain_grid = terrain_pointcloud == ipc_index});
     ipc.min_ground_elevation = PointCloudCropper->get_min_terrain_elevation();
+    if (terrain_pointcloud == ipc_index) {
+      const auto* terrain_grid = PointCloudCropper->get_terrain_grid();
+      if (terrain_grid != nullptr) {
+        auto nodata_mode = roofer::io::TerrainNoDataMode::COMPLETE_QUADS;
+        if (cfg.terrain_nodata_mode == "local_triangles") {
+          nodata_mode = roofer::io::TerrainNoDataMode::LOCAL_TRIANGLES;
+        } else if (cfg.terrain_nodata_mode == "fill_small_gaps") {
+          nodata_mode = roofer::io::TerrainNoDataMode::FILL_SMALL_GAPS;
+        }
+        auto triangles =
+            roofer::io::triangulateTerrainGrid(*terrain_grid, nodata_mode);
+        if (triangles.empty()) {
+          logger.warning(
+              "No complete terrain grid quads found for tile {} using "
+              "pointcloud {}",
+              output_building_tile.id, ipc.name);
+        } else {
+          TerrainData terrain;
+          terrain.components =
+              roofer::io::splitTerrainConnectedComponents(triangles);
+          terrain.attributes.insert("rf_pc_source", ipc.name);
+          terrain.attributes.insert("rf_pc_quality", ipc.quality);
+          terrain.attributes.insert("rf_pc_date", ipc.date);
+          terrain.attributes.insert("rf_ground_class", ipc.grnd_class);
+          terrain.attributes.insert("rf_terrain_grid_cellsize",
+                                    cfg.terrain_grid_cellsize);
+          terrain.attributes.insert("rf_terrain_aggregation",
+                                    std::string("cell_minimum"));
+          terrain.attributes.insert("rf_terrain_interpolation",
+                                    std::string("cell_center_linear"));
+          terrain.attributes.insert("rf_terrain_nodata_mode",
+                                    cfg.terrain_nodata_mode);
+          output_building_tile.terrain = std::move(terrain);
+        }
+      }
+    }
     if (ipc.date != 0) {
       logger.info("Overriding acquisition year from config file");
       std::fill(ipc.acquisition_years.begin(), ipc.acquisition_years.end(),

apps/roofer-app/roofer-app.cpp

@@ -143,6 +143,11 @@ struct BuildingObject {
   // bool was_skipped;  // b3_reconstructie_onvolledig;
 };
 
+struct TerrainData {
+  std::vector<std::vector<roofer::LinearRing>> components;
+  roofer::AttributeMapRow attributes;
+};
+
 /**
  * @brief Indicates the progress of the object through the whole roofer process.
  */
@@ -204,6 +209,7 @@ struct BuildingTile {
   std::unique_ptr<roofer::misc::projHelperInterface> proj_helper;
   // extent
   roofer::TBox<double> extent;
+  std::optional<TerrainData> terrain;
 
   std::vector<std::pair<Progress, size_t>> count_progresses() const;
 };
@@ -937,6 +943,29 @@ int main(int argc, const char* argv[]) {
             }
           }
 
+          if (building_tile.terrain.has_value()) {
+            const auto terrain_id = fmt::format("terrain-{}", building_tile.id);
+            if (handler.cfg_.split_cjseq) {
+              const int minx = building_tile.extent.min()[0];
+              const int miny = building_tile.extent.min()[1];
+              const auto terrain_path =
+                  fs::path(handler.cfg_.output_path) /
+                  fmt::format("{:06d}_{:06d}.terrain.city.jsonl", minx, miny);
+              fs::create_directories(terrain_path.parent_path());
+              std::ofstream terrain_ofs(terrain_path);
+              CityJsonWriter->write_metadata(terrain_ofs, project_srs.get(),
+                                             building_tile.extent,
+                                             {.identifier = terrain_id});
+              CityJsonWriter->write_tin_relief_feature(
+                  terrain_ofs, terrain_id, building_tile.terrain->components,
+                  building_tile.terrain->attributes);
+            } else {
+              CityJsonWriter->write_tin_relief_feature(
+                  ofs, terrain_id, building_tile.terrain->components,
+                  building_tile.terrain->attributes);
+            }
+          }
+
           for (auto& building : building_tile.buildings) {
             if (handler.cfg_.split_cjseq) {
               fs::create_directories(building.jsonl_path.parent_path());

docs/cli_application.md

@@ -50,7 +50,9 @@ Output directory. The building models will be written to a CityJSONSequence file
 ```
 
 ## Output format
-The output of the {program}`roofer` CLI application are [CityJSONSequence](https://www.cityjson.org/cityjsonseq/) files. These are a JSON Lines files that contain a sequence of CityJSON features, each feature represents all the information for one building.
+The output of the {program}`roofer` CLI application are [CityJSONSequence](https://www.cityjson.org/cityjsonseq/) files. These are JSON Lines files that contain a sequence of CityJSON features. Each feature represents one building, with an optional per-tile `TINRelief` terrain feature when `--terrain` is enabled.
+
+For terrain output, `--terrain-nodata-mode` controls whether only complete grid quads are triangulated, three-sample local triangles are retained, or isolated missing samples are interpolated before local triangulation.
 
 > [!TIP]
 >  In the near future we expect to introduce a converter program to convert the CityJSONSequence files to other formats like GPKG, OBJ and others.

include/roofer/io/CityJsonWriter.hpp

@@ -78,6 +78,10 @@ namespace roofer::io {
         const std::unordered_map<int, Mesh>* geometry_lod13,
         const std::unordered_map<int, Mesh>* geometry_lod22,
         const AttributeMapRow attributes) = 0;
+    virtual void write_tin_relief_feature(
+        std::ostream& output_stream, const std::string& id,
+        const std::vector<std::vector<LinearRing>>& components,
+        const AttributeMapRow& attributes) = 0;
 
     // virtual void write(const std::string& source, const LinearRing&
     // footprints,

include/roofer/io/StreamCropper.hpp

@@ -21,10 +21,17 @@
 
 #pragma once
 #include <memory>
+#include <roofer/common/Raster.hpp>
 #include <roofer/common/datastructures.hpp>
 #include <roofer/misc/projHelper.hpp>
 
 namespace roofer::io {
+  enum class TerrainNoDataMode {
+    COMPLETE_QUADS,
+    LOCAL_TRIANGLES,
+    FILL_SMALL_GAPS,
+  };
+
   struct PointCloudCropperConfig {
     float cellsize = 50.0;
     float buffer = 1.0;
@@ -43,6 +50,7 @@ namespace roofer::io {
     bool use_acquisition_year = true;
     float terrain_grid_cellsize = 10.0;
     int terrain_grid_search_radius = 3;
+    bool retain_terrain_grid = false;
   };
   struct PointCloudCropperInterface {
     roofer::misc::projHelperInterface& pjHelper;
@@ -61,6 +69,7 @@ namespace roofer::io {
         PointCloudCropperConfig cfg = PointCloudCropperConfig{}) = 0;
 
     virtual std::optional<float> get_min_terrain_elevation() const = 0;
+    virtual const RasterTools::Raster* get_terrain_grid() const = 0;
     // virtual void process(
     //     std::string filepaths, std::vector<LinearRing>& polygons,
     //     std::vector<LinearRing>& buf_polygons,
@@ -71,4 +80,10 @@ namespace roofer::io {
 
   std::unique_ptr<PointCloudCropperInterface> createPointCloudCropper(
       roofer::misc::projHelperInterface& pjh);
+
+  std::vector<LinearRing> triangulateTerrainGrid(
+      const RasterTools::Raster& terrain_grid,
+      TerrainNoDataMode nodata_mode = TerrainNoDataMode::COMPLETE_QUADS);
+  std::vector<std::vector<LinearRing>> splitTerrainConnectedComponents(
+      const std::vector<LinearRing>& triangles);
 }  // namespace roofer::io

include/roofer/reconstruction/PlaneDetector.hpp

@@ -51,7 +51,8 @@ namespace roofer::reconstruction {
 
     // limits
     bool with_limits = false;
-    // Deterministic complexity backstop: abort region growing (and fall back to LoD 1.1) once this many planes are detected.
+    // Deterministic complexity backstop: abort region growing (and fall back to
+    // LoD 1.1) once this many planes are detected.
     int limit_n_regions = 100;
   };
 

src/extra/io/CityJsonWriter.cpp

@@ -452,6 +452,54 @@ namespace roofer::io {
 
       write_to_stream(outputJSON, output_stream, prettyPrint_);
     }
+
+    void write_tin_relief_feature(
+        std::ostream& output_stream, const std::string& id,
+        const std::vector<std::vector<LinearRing>>& components,
+        const AttributeMapRow& attributes) override {
+      nlohmann::json outputJSON;
+      outputJSON["type"] = "CityJSONFeature";
+      outputJSON["id"] = id;
+      outputJSON["CityObjects"] = nlohmann::json::object();
+
+      std::map<arr3d, size_t> vertex_map;
+      std::set<arr3d> vertex_set;
+      std::vector<arr3d> vertex_vec;
+      TBox<double> terrain_bbox;
+      auto geometries = nlohmann::json::array();
+      for (const auto& component : components) {
+        std::vector<std::vector<std::vector<size_t>>> boundaries;
+        boundaries.reserve(component.size());
+        for (const auto& triangle : component) {
+          terrain_bbox.add(add_vertices_polygon(vertex_map, vertex_vec,
+                                                vertex_set, triangle));
+          boundaries.push_back(LinearRing2jboundary(vertex_map, triangle));
+        }
+
+        auto geometry = nlohmann::json::object();
+        geometry["type"] = "CompositeSurface";
+        geometry["lod"] = "1";
+        geometry["boundaries"] = boundaries;
+        geometries.push_back(std::move(geometry));
+      }
+
+      auto terrain = nlohmann::json::object();
+      terrain["type"] = "TINRelief";
+      terrain["attributes"] = attributes2json(attributes);
+      terrain["geographicalExtent"] = compute_geographical_extent(terrain_bbox);
+      terrain["geometry"] = std::move(geometries);
+      outputJSON["CityObjects"][id] = terrain;
+
+      std::vector<std::array<int, 3>> vertices_int;
+      vertices_int.reserve(vertex_vec.size());
+      for (const auto& vertex : vertex_vec) {
+        vertices_int.push_back({int((vertex[0] - translate_x_) / scale_x_),
+                                int((vertex[1] - translate_y_) / scale_y_),
+                                int((vertex[2] - translate_z_) / scale_z_)});
+      }
+      outputJSON["vertices"] = vertices_int;
+      write_to_stream(outputJSON, output_stream, prettyPrint_);
+    }
   };
 
   std::unique_ptr<CityJsonWriterInterface> createCityJsonWriter(