diff --git a/apps/roofer-app/config.hpp b/apps/roofer-app/config.hpp index 0c36cfcc..048fd1c4 100644 --- a/apps/roofer-app/config.hpp +++ b/apps/roofer-app/config.hpp @@ -58,10 +58,12 @@ struct InputPointcloud { roofer::vec1f pt_densities; roofer::vec1b is_glass_roof; roofer::vec1b lod11_forced; + roofer::vec1b pointcloud_insufficient; std::vector nodata_circles; std::vector building_clouds; std::vector building_rasters; roofer::vec1f ground_elevations; + roofer::vec1f roof_elevations; roofer::vec1i acquisition_years; std::unique_ptr rtree; @@ -84,7 +86,7 @@ struct RooferConfig { int lod11_fallback_area = 69000; float lod11_fallback_density = 5; roofer::arr2f tilesize = {1000, 1000}; - bool clear_if_insufficient = false; + bool clear_if_insufficient = true; bool write_crop_outputs = false; bool output_all = false; @@ -115,12 +117,14 @@ struct RooferConfig { std::string output_path; // reconstruct + int lod11_fallback_planes = 900; + int lod11_fallback_time = 1800000; roofer::ReconstructionConfig rec; // output attribute names std::unordered_map n = { - {"status", "rf_status"}, - {"reconstruction_time", "rf_reconstruction_time"}, + {"success", "rf_success"}, + {"reconstruction_time", "rf_t_run"}, {"val3dity_lod12", "rf_val3dity_lod12"}, {"val3dity_lod13", "rf_val3dity_lod13"}, {"val3dity_lod22", "rf_val3dity_lod22"}, @@ -138,7 +142,7 @@ struct RooferConfig { {"h_roof_70p", "rf_h_roof_70p"}, {"h_roof_min", "rf_h_roof_min"}, {"h_roof_max", "rf_h_roof_max"}, - {"roof_n_planes", "rf_roof_n_planes"}, + {"roof_n_planes", "rf_roof_planes"}, {"rmse_lod12", "rf_rmse_lod12"}, {"rmse_lod13", "rf_rmse_lod13"}, {"rmse_lod22", "rf_rmse_lod22"}, @@ -148,6 +152,8 @@ struct RooferConfig { {"h_ground", "rf_h_ground"}, {"slope", "rf_slope"}, {"azimuth", "rf_azimuth"}, + {"extrusion_mode", "rf_extrusion_mode"}, + {"pointcloud_unusable", "rf_pointcloud_unusable"}, }; }; @@ -642,7 +648,8 @@ struct RooferConfigHandler { _cfg.cellsize, {roofer::v::HigherThan(0)}); add("lod11-fallback-area", "lod11 fallback area", _cfg.lod11_fallback_area, {roofer::v::HigherThan(0)}); - add("skip-insufficient", "skip buildings with insufficient pointcloud data", + add("reconstruct-insufficient", + "reconstruct buildings with insufficient pointcloud data", _cfg.clear_if_insufficient, {}); // add("lod11-fallback-density", "lod11 fallback density", // _cfg.lod11_fallback_density, {roofer::v::HigherThan(0)}}); @@ -682,6 +689,13 @@ struct RooferConfigHandler { _cfg.cj_scale, {}); add("cj-translate", "Translation applied to CityJSON output vertices", _cfg.cj_translate, {}); + add("lod11-fallback-plane-cnt", + "Fallback to LoD11 if number of detected planes exceeds this value.", + _cfg.lod11_fallback_planes, {roofer::v::HigherThan(0)}); + add("lod11-fallback-time", + "Fallback to LoD11 if time spent on detecting planes exceeds this " + "value. In milliseconds.", + _cfg.lod11_fallback_time, {roofer::v::HigherThan(0)}); addr("plane-detect-k", "plane detect k", _cfg.rec.plane_detect_k, {roofer::v::HigherThan(0)}); addr("plane-detect-min-points", "plane detect min points", diff --git a/apps/roofer-app/crop_tile.hpp b/apps/roofer-app/crop_tile.hpp index ec32ab3b..fe12c944 100644 --- a/apps/roofer-app/crop_tile.hpp +++ b/apps/roofer-app/crop_tile.hpp @@ -94,7 +94,8 @@ bool crop_tile(const roofer::TBox& tile, PointCloudCropper->process( lasfiles, footprints, buffered_footprints, ipc.building_clouds, - ipc.ground_elevations, ipc.acquisition_years, polygon_extent, + ipc.ground_elevations, ipc.acquisition_years, + ipc.pointcloud_insufficient, polygon_extent, {.ground_class = ipc.grnd_class, .building_class = ipc.bld_class, .clear_if_insufficient = cfg.clear_if_insufficient, @@ -130,7 +131,9 @@ bool crop_tile(const roofer::TBox& tile, ipc.nodata_fractions.resize(N_fp); ipc.pt_densities.resize(N_fp); ipc.is_glass_roof.reserve(N_fp); + ipc.roof_elevations.reserve(N_fp); ipc.lod11_forced.reserve(N_fp); + ipc.pointcloud_insufficient.reserve(N_fp); if (cfg.write_index) ipc.nodata_circles.resize(N_fp); // auto& r_nodata = attributes.insert_vec("r_nodata_"+ipc.name); @@ -145,6 +148,8 @@ bool crop_tile(const roofer::TBox& tile, roofer::misc::computePointDensity(ipc.building_rasters[i]); ipc.is_glass_roof[i] = roofer::misc::testForGlassRoof(ipc.building_rasters[i]); + ipc.roof_elevations[i] = + roofer::misc::computeRoofElevation(ipc.building_rasters[i], 0.7); auto target_density = cfg.ceil_point_density; bool do_force_lod11 = @@ -187,27 +192,19 @@ bool crop_tile(const roofer::TBox& tile, // add raster stats attributes from PointCloudCropper to footprint attributes for (auto& ipc : input_pointclouds) { - auto& h_ground = - attributes.insert_vec(cfg.n.at("h_ground") + "_" + ipc.name); auto& nodata_r = attributes.insert_vec(cfg.n.at("nodata_r") + "_" + ipc.name); auto& nodata_frac = attributes.insert_vec(cfg.n.at("nodata_frac") + "_" + ipc.name); auto& pt_density = attributes.insert_vec(cfg.n.at("pt_density") + "_" + ipc.name); - auto& is_glass_roof = - attributes.insert_vec(cfg.n.at("is_glass_roof") + "_" + ipc.name); - h_ground.reserve(N_fp); nodata_r.reserve(N_fp); nodata_frac.reserve(N_fp); 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] + (*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]); - is_glass_roof.push_back(ipc.is_glass_roof[i]); } } @@ -346,7 +343,18 @@ bool crop_tile(const roofer::TBox& tile, building.footprint = footprints[i]; building.h_ground = input_pointclouds[selected->index].ground_elevations[i]; + building.h_roof_70p_rough = + input_pointclouds[selected->index].roof_elevations[i]; building.force_lod11 = input_pointclouds[selected->index].lod11_forced[i]; + building.pointcloud_insufficient = + input_pointclouds[selected->index].pointcloud_insufficient[i]; + building.is_glass_roof = + input_pointclouds[selected->index].is_glass_roof[i]; + + if (input_pointclouds[selected->index].lod11_forced[i]) { + building.extrusion_mode = ExtrusionMode::LOD11_FALLBACK; + force_lod11_vec[i] = input_pointclouds[selected->index].lod11_forced[i]; + } output_building_tile.attributes = attributes; building.jsonl_path = fmt::format( @@ -354,9 +362,6 @@ bool crop_tile(const roofer::TBox& tile, fmt::arg("pc_name", input_pointclouds[selected->index].name), fmt::arg("path", cfg.output_path)); } - if (input_pointclouds[selected->index].lod11_forced[i]) { - force_lod11_vec[i] = input_pointclouds[selected->index].lod11_forced[i]; - } if (cfg.write_crop_outputs) { { diff --git a/apps/roofer-app/reconstruct_building.hpp b/apps/roofer-app/reconstruct_building.hpp index 7c47c098..5091504b 100644 --- a/apps/roofer-app/reconstruct_building.hpp +++ b/apps/roofer-app/reconstruct_building.hpp @@ -19,18 +19,70 @@ // Author(s): // Ravi Peters -enum LOD { LOD12 = 12, LOD13 = 13, LOD22 = 22 }; +enum LOD { LOD11 = 11, LOD12 = 12, LOD13 = 13, LOD22 = 22 }; + +void multisolid_post_process(BuildingObject& building, RooferConfig* rfcfg, + LOD lod, + std::unordered_map& multisolid, + std::optional& rmse, + std::optional& volume, + std::optional& attr_val3dity) { + auto MeshPropertyCalculator = roofer::misc::createMeshPropertyCalculator(); + for (auto& [label, mesh] : multisolid) { + mesh.get_attributes().resize(mesh.get_polygons().size()); + MeshPropertyCalculator->compute_roof_height( + mesh, {.z_offset = building.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"]}); + if (lod == LOD22) { + MeshPropertyCalculator->compute_roof_orientation( + mesh, {.slope = rfcfg->n["slope"], .azimuth = rfcfg->n["azimuth"]}); + } + } + + auto MeshTriangulator = + roofer::reconstruction::createMeshTriangulatorLegacy(); + MeshTriangulator->compute(multisolid); + volume = MeshTriangulator->volumes.at(0); + // logger.debug("Completed MeshTriangulator"); +#ifdef RF_USE_RERUN + rec.log(worldname + "MeshTriangulator", + rerun::Mesh3D(MeshTriangulator->triangles) + .with_vertex_normals(MeshTriangulator->normals) + .with_class_ids(MeshTriangulator->ring_ids)); +#endif + + auto PC2MeshDistCalculator = roofer::misc::createPC2MeshDistCalculator(); + PC2MeshDistCalculator->compute(building.pointcloud_building, + MeshTriangulator->multitrianglecol, + MeshTriangulator->ring_ids); + rmse = PC2MeshDistCalculator->rms_error; + // logger.debug("Completed PC2MeshDistCalculator. RMSE={}", + // PC2MeshDistCalculator->rms_error); +#ifdef RF_USE_RERUN +// rec.log(worldname+"PC2MeshDistCalculator", +// rerun::Mesh3D(PC2MeshDistCalculator->triangles).with_vertex_normals(MeshTriangulator->normals).with_class_ids(MeshTriangulator->ring_ids)); +#endif -std::unordered_map extrude( - roofer::Arrangement_2 arrangement, const float& floor_elevation, - const float& z_offset, #ifdef RF_USE_VAL3DITY - std::string& attr_val3dity, + if (multisolid.size() > 0) { + auto Val3dator = roofer::misc::createVal3dator(); + Val3dator->compute(multisolid); + attr_val3dity = Val3dator->errors.front(); + } + // logger.debug("Completed Val3dator. Errors={}", Val3dator->errors.front()); #endif - float& rmse, float& volume, +} + +std::unordered_map extrude_lod22( + roofer::Arrangement_2 arrangement, BuildingObject& building, + RooferConfig* rfcfg, roofer::reconstruction::SegmentRasteriserInterface* SegmentRasteriser, - roofer::reconstruction::PlaneDetectorInterface* PlaneDetector, LOD lod, - RooferConfig* rfcfg) { + LOD lod, std::optional& rmse, std::optional& volume, + std::optional& attr_val3dity) { auto* cfg = &(rfcfg->rec); bool dissolve_step_edges = false; bool dissolve_all_interior = false; @@ -71,7 +123,7 @@ std::unordered_map extrude( auto ArrangementExtruder = roofer::reconstruction::createArrangementExtruder(); - ArrangementExtruder->compute(arrangement, floor_elevation, + ArrangementExtruder->compute(arrangement, building.h_ground, {.LoD2 = extrude_LoD2}); // logger.debug("Completed ArrangementExtruder"); #ifdef RF_USE_RERUN @@ -80,65 +132,41 @@ std::unordered_map 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"]}); - if (lod == LOD22) { - 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", - rerun::Mesh3D(MeshTriangulator->triangles) - .with_vertex_normals(MeshTriangulator->normals) - .with_class_ids(MeshTriangulator->ring_ids)); -#endif + multisolid_post_process(building, rfcfg, lod, ArrangementExtruder->multisolid, + rmse, volume, attr_val3dity); - auto PC2MeshDistCalculator = roofer::misc::createPC2MeshDistCalculator(); - PC2MeshDistCalculator->compute(PlaneDetector->pts_per_roofplane, - MeshTriangulator->multitrianglecol, - MeshTriangulator->ring_ids); - rmse = PC2MeshDistCalculator->rms_error; - // logger.debug("Completed PC2MeshDistCalculator. RMSE={}", - // PC2MeshDistCalculator->rms_error); -#ifdef RF_USE_RERUN -// rec.log(worldname+"PC2MeshDistCalculator", -// rerun::Mesh3D(PC2MeshDistCalculator->triangles).with_vertex_normals(MeshTriangulator->normals).with_class_ids(MeshTriangulator->ring_ids)); -#endif + return ArrangementExtruder->multisolid; +} +void extrude_lod11(BuildingObject& building, RooferConfig* rfcfg) { + auto SimplePolygonExtruder = + roofer::reconstruction::createSimplePolygonExtruder(); + SimplePolygonExtruder->compute(building.footprint, building.h_ground, + building.h_roof_70p_rough); + // std::vector> multisolidvec; + building.multisolids_lod12 = SimplePolygonExtruder->multisolid; + building.multisolids_lod13 = SimplePolygonExtruder->multisolid; + building.multisolids_lod22 = SimplePolygonExtruder->multisolid; + + multisolid_post_process( + building, rfcfg, LOD11, SimplePolygonExtruder->multisolid, + building.rmse_lod12, building.volume_lod12, building.val3dity_lod12); + building.rmse_lod13 = building.rmse_lod12; + building.rmse_lod22 = building.rmse_lod12; + building.volume_lod13 = building.volume_lod12; + building.volume_lod22 = building.volume_lod12; #ifdef RF_USE_VAL3DITY - if (ArrangementExtruder->multisolid.size() > 0) { - auto Val3dator = roofer::misc::createVal3dator(); - Val3dator->compute(ArrangementExtruder->multisolid); - attr_val3dity = Val3dator->errors.front(); - } - // logger.debug("Completed Val3dator. Errors={}", Val3dator->errors.front()); + building.val3dity_lod13 = building.val3dity_lod12; + building.val3dity_lod22 = building.val3dity_lod12; #endif - return ArrangementExtruder->multisolid; + building.extrusion_mode = LOD11_FALLBACK; + building.roof_elevation_70p = building.h_roof_70p_rough; } void reconstruct_building(BuildingObject& building, RooferConfig* rfcfg) { auto* cfg = &(rfcfg->rec); auto& logger = roofer::logger::Logger::get_logger(); - // split into ground and roof points - - // logger.debug("{} ground points and {} roof points", points_ground.size(), - // points_roof.size()); #ifdef RF_USE_RERUN // Create a new `RecordingStream` which sends data over TCP to the viewer @@ -165,73 +193,66 @@ void reconstruct_building(BuildingObject& building, RooferConfig* rfcfg) { std::unordered_map> timings; - auto t0 = std::chrono::high_resolution_clock::now(); - auto PlaneDetector = roofer::reconstruction::createPlaneDetector(); - PlaneDetector->detect( - building.pointcloud_building, - { - .metrics_plane_k = cfg->plane_detect_k, - .metrics_plane_min_points = cfg->plane_detect_min_points, - .metrics_plane_epsilon = cfg->plane_detect_epsilon, - .metrics_plane_normal_threshold = cfg->plane_detect_normal_angle, - }); - timings["PlaneDetector"] = std::chrono::high_resolution_clock::now() - t0; - // logger.debug("Completed PlaneDetector (roof), found {} roofplanes", - // PlaneDetector->pts_per_roofplane.size()); - - building.roof_type = PlaneDetector->roof_type; - building.roof_elevation_50p = PlaneDetector->roof_elevation_50p; - building.roof_elevation_70p = PlaneDetector->roof_elevation_70p; - building.roof_elevation_min = PlaneDetector->roof_elevation_min; - building.roof_elevation_max = PlaneDetector->roof_elevation_max; - building.roof_n_planes = PlaneDetector->pts_per_roofplane.size(); - - bool pointcloud_insufficient = PlaneDetector->roof_type == "no points" || - PlaneDetector->roof_type == "no planes"; - if (pointcloud_insufficient) { - // building.was_skipped = true; - // TODO: return building status pointcloud_insufficient - // CityJsonWriter->write(path_output_jsonl, footprints, nullptr, nullptr, - // nullptr, attributes); - return; + if (building.pointcloud_insufficient) { + building.extrusion_mode = SKIP; } - t0 = std::chrono::high_resolution_clock::now(); - auto PlaneDetector_ground = roofer::reconstruction::createPlaneDetector(); - PlaneDetector_ground->detect(building.pointcloud_ground); - timings["PlaneDetector_ground"] = - std::chrono::high_resolution_clock::now() - t0; - // logger.debug("Completed PlaneDetector (ground), found {} groundplanes", - // PlaneDetector_ground->pts_per_roofplane.size()); - -#ifdef RF_USE_RERUN - rec.log("world/segmented_points", - rerun::Collection{rerun::components::AnnotationContext{ - rerun::datatypes::AnnotationInfo( - 0, "no plane", rerun::datatypes::Rgba32(30, 30, 30))}}); - rec.log("world/segmented_points", - rerun::Points3D(points_roof).with_class_ids(PlaneDetector->plane_id)); -#endif - - // check skip_attribute - // logger.debug("force LoD1.1 = {}", building.force_lod11); - - if (building.force_lod11) { - auto SimplePolygonExtruder = - roofer::reconstruction::createSimplePolygonExtruder(); - SimplePolygonExtruder->compute(building.footprint, building.h_ground, - PlaneDetector->roof_elevation_70p); - std::vector> multisolidvec; - building.multisolids_lod12 = SimplePolygonExtruder->multisolid; - building.multisolids_lod13 = SimplePolygonExtruder->multisolid; - building.multisolids_lod22 = SimplePolygonExtruder->multisolid; - // TODO: return building status + if (building.extrusion_mode == SKIP) { + return; + } else if (building.extrusion_mode == LOD11_FALLBACK) { + extrude_lod11(building, rfcfg); return; - // multisolidvec.push_back(SimplePolygonExtruder->multisolid); - // CityJsonWriter->write(path_output_jsonl, footprints, &multisolidvec, - // &multisolidvec, &multisolidvec, attributes); - // logger.debug("Completed CityJsonWriter to {}", path_output_jsonl); } else { + auto t0 = std::chrono::high_resolution_clock::now(); + auto PlaneDetector = roofer::reconstruction::createPlaneDetector(); + auto PlaneDetector_ground = roofer::reconstruction::createPlaneDetector(); + try { + auto plane_detector_cfg = roofer::reconstruction::PlaneDetectorConfig{ + .metrics_plane_k = cfg->plane_detect_k, + .metrics_plane_min_points = cfg->plane_detect_min_points, + .metrics_plane_epsilon = cfg->plane_detect_epsilon, + .metrics_plane_normal_threshold = cfg->plane_detect_normal_angle, + .with_limits = true, + .limit_n_regions = rfcfg->lod11_fallback_planes, + .limit_n_milliseconds = rfcfg->lod11_fallback_time, + }; + PlaneDetector->detect(building.pointcloud_building, plane_detector_cfg); + timings["PlaneDetector"] = std::chrono::high_resolution_clock::now() - t0; + t0 = std::chrono::high_resolution_clock::now(); + PlaneDetector_ground->detect(building.pointcloud_ground, + plane_detector_cfg); + timings["PlaneDetector_ground"] = + std::chrono::high_resolution_clock::now() - t0; + + building.roof_type = PlaneDetector->roof_type; + building.roof_elevation_50p = PlaneDetector->roof_elevation_50p; + building.roof_elevation_70p = PlaneDetector->roof_elevation_70p; + building.roof_elevation_min = PlaneDetector->roof_elevation_min; + building.roof_elevation_max = PlaneDetector->roof_elevation_max; + building.roof_n_planes = PlaneDetector->pts_per_roofplane.size(); + + bool pointcloud_insufficient = PlaneDetector->roof_type == "no points" || + PlaneDetector->roof_type == "no planes"; + if (pointcloud_insufficient) { + building.extrusion_mode = SKIP; + return; + } + } catch (const std::runtime_error& e) { + extrude_lod11(building, rfcfg); + logger.warning("[reconstructor] {}, LoD1.1 fallback: {}", + building.jsonl_path.string(), e.what()); + return; + } + // #ifdef RF_USE_RERUN + // rec.log("world/segmented_points", + // rerun::Collection{rerun::components::AnnotationContext{ + // rerun::datatypes::AnnotationInfo( + // 0, "no plane", rerun::datatypes::Rgba32(30, 30, + // 30))}}); + // rec.log( + // "world/segmented_points", + // rerun::Points3D(points_roof).with_class_ids(PlaneDetector->plane_id)); + // #endif t0 = std::chrono::high_resolution_clock::now(); auto AlphaShaper = roofer::reconstruction::createAlphaShaper(); AlphaShaper->compute(PlaneDetector->pts_per_roofplane, @@ -354,33 +375,21 @@ void reconstruct_building(BuildingObject& building, RooferConfig* rfcfg) { // logger.debug("LoD={}", cfg->lod); t0 = std::chrono::high_resolution_clock::now(); if (cfg->lod == 0 || cfg->lod == 12) { - building.multisolids_lod12 = - extrude(arrangement, building.h_ground, building.z_offset, -#ifdef RF_USE_VAL3DITY - building.val3dity_lod12, -#endif - building.rmse_lod12, building.volume_lod12, - SegmentRasteriser.get(), PlaneDetector.get(), LOD12, rfcfg); + building.multisolids_lod12 = extrude_lod22( + arrangement, building, rfcfg, SegmentRasteriser.get(), LOD12, + building.rmse_lod12, building.volume_lod12, building.val3dity_lod12); } if (cfg->lod == 0 || cfg->lod == 13) { - building.multisolids_lod13 = - extrude(arrangement, building.h_ground, building.z_offset, -#ifdef RF_USE_VAL3DITY - building.val3dity_lod13, -#endif - building.rmse_lod13, building.volume_lod13, - SegmentRasteriser.get(), PlaneDetector.get(), LOD13, rfcfg); + building.multisolids_lod13 = extrude_lod22( + arrangement, building, rfcfg, SegmentRasteriser.get(), LOD13, + building.rmse_lod13, building.volume_lod13, building.val3dity_lod13); } if (cfg->lod == 0 || cfg->lod == 22) { - building.multisolids_lod22 = - extrude(arrangement, building.h_ground, building.z_offset, -#ifdef RF_USE_VAL3DITY - building.val3dity_lod22, -#endif - building.rmse_lod22, building.volume_lod22, - SegmentRasteriser.get(), PlaneDetector.get(), LOD22, rfcfg); + building.multisolids_lod22 = extrude_lod22( + arrangement, building, rfcfg, SegmentRasteriser.get(), LOD22, + building.rmse_lod22, building.volume_lod22, building.val3dity_lod22); } timings["extrude"] = std::chrono::high_resolution_clock::now() - t0; diff --git a/apps/roofer-app/roofer-app.cpp b/apps/roofer-app/roofer-app.cpp index e05d94cf..ba5f8619 100644 --- a/apps/roofer-app/roofer-app.cpp +++ b/apps/roofer-app/roofer-app.cpp @@ -85,6 +85,8 @@ namespace fs = std::filesystem; #include "config.hpp" +enum ExtrusionMode { STANDARD, LOD11_FALLBACK, SKIP }; + /** * @brief A single building object * @@ -108,24 +110,29 @@ struct BuildingObject { // set in crop fs::path jsonl_path; float h_ground; + float h_roof_70p_rough; bool force_lod11; // force_lod11 / fallback_lod11 + bool pointcloud_insufficient; + bool is_glass_roof; + ExtrusionMode extrusion_mode = STANDARD; // set in reconstruction - std::string roof_type; - float roof_elevation_50p; - float roof_elevation_70p; - float roof_elevation_min; - float roof_elevation_max; - int roof_n_planes; - 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; + // optionals may not get assigned a valid value + std::string roof_type = "unknown"; + std::optional roof_elevation_50p; + std::optional roof_elevation_70p; + std::optional roof_elevation_min; + std::optional roof_elevation_max; + std::optional roof_n_planes; + std::optional rmse_lod12; + std::optional rmse_lod13; + std::optional rmse_lod22; + std::optional volume_lod12; + std::optional volume_lod13; + std::optional volume_lod22; + std::optional val3dity_lod12; + std::optional val3dity_lod13; + std::optional val3dity_lod22; // bool was_skipped; // b3_reconstructie_onvolledig; }; @@ -980,15 +987,13 @@ int main(int argc, const char* argv[]) { logger.debug("[serializer] Serializing tile {}", building_tile); // create status attribute - auto& attr_status = building_tile.attributes.insert_vec( - roofer_cfg.n["status"]); + auto& attr_success = building_tile.attributes.insert_vec( + roofer_cfg.n["success"]); for (auto& progress : building_tile.buildings_progresses) { - if (progress == RECONSTRUCTION_FAILED) { - attr_status.push_back("reconstruction_failed"); - } else if (progress == RECONSTRUCTION_SUCCEEDED) { - attr_status.push_back("reconstruction_succeeded"); + if (progress == RECONSTRUCTION_SUCCEEDED) { + attr_success.push_back(true); } else { - attr_status.push_back("unknown"); + attr_success.push_back(false); } } // create time attribute @@ -1066,49 +1071,58 @@ int main(int argc, const char* argv[]) { building.attribute_index); attrow.insert(roofer_cfg.n["h_ground"], building.h_ground); + attrow.insert(roofer_cfg.n["is_glass_roof"], + building.is_glass_roof); + attrow.insert(roofer_cfg.n["pointcloud_unusable"], + building.pointcloud_insufficient); attrow.insert(roofer_cfg.n["roof_type"], building.roof_type); - attrow.insert(roofer_cfg.n["h_roof_50p"], - building.roof_elevation_50p); - attrow.insert(roofer_cfg.n["h_roof_70p"], - building.roof_elevation_70p); - attrow.insert(roofer_cfg.n["h_roof_min"], - building.roof_elevation_min); - attrow.insert(roofer_cfg.n["h_roof_max"], - building.roof_elevation_max); - attrow.insert(roofer_cfg.n["roof_n_planes"], - building.roof_n_planes); + attrow.insert_optional(roofer_cfg.n["h_roof_50p"], + building.roof_elevation_50p); + attrow.insert_optional(roofer_cfg.n["h_roof_70p"], + building.roof_elevation_70p); + attrow.insert_optional(roofer_cfg.n["h_roof_min"], + building.roof_elevation_min); + attrow.insert_optional(roofer_cfg.n["h_roof_max"], + building.roof_elevation_max); + attrow.insert_optional(roofer_cfg.n["roof_n_planes"], + building.roof_n_planes); + attrow.insert(roofer_cfg.n["extrusion_mode"], + building.extrusion_mode); std::unordered_map* ms12 = nullptr; std::unordered_map* ms13 = nullptr; std::unordered_map* ms22 = nullptr; 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); + attrow.insert_optional(roofer_cfg.n["rmse_lod12"], + building.rmse_lod12); + attrow.insert_optional(roofer_cfg.n["volume_lod12"], + building.volume_lod12); #if RF_USE_VAL3DITY - attrow.insert(roofer_cfg.n["val3dity_lod12"], - building.val3dity_lod12); + attrow.insert_optional(roofer_cfg.n["val3dity_lod12"], + building.val3dity_lod12); #endif } 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); + attrow.insert_optional(roofer_cfg.n["rmse_lod13"], + building.rmse_lod13); + attrow.insert_optional(roofer_cfg.n["volume_lod13"], + building.volume_lod13); #if RF_USE_VAL3DITY - attrow.insert(roofer_cfg.n["val3dity_lod13"], - building.val3dity_lod13); + attrow.insert_optional(roofer_cfg.n["val3dity_lod13"], + building.val3dity_lod13); #endif } 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); + attrow.insert_optional(roofer_cfg.n["rmse_lod22"], + building.rmse_lod22); + attrow.insert_optional(roofer_cfg.n["volume_lod22"], + building.volume_lod22); #if RF_USE_VAL3DITY - attrow.insert(roofer_cfg.n["val3dity_lod22"], - building.val3dity_lod22); + attrow.insert_optional(roofer_cfg.n["val3dity_lod22"], + building.val3dity_lod22); #endif } CityJsonWriter->write_feature(ofs, building.footprint, ms12, ms13, diff --git a/include/roofer/common/common.hpp b/include/roofer/common/common.hpp index baf944cf..ed3be0d6 100644 --- a/include/roofer/common/common.hpp +++ b/include/roofer/common/common.hpp @@ -139,6 +139,13 @@ namespace roofer { void insert(const std::string& name, T value) { _attributes[name] = value; }; + template + void insert_optional(const std::string& name, std::optional opt) { + if (opt.has_value()) + _attributes[name] = opt.value(); + else + _attributes[name] = std::monostate(); + }; void set_null(const std::string& name); bool is_null(const std::string& name) const; diff --git a/include/roofer/io/StreamCropper.hpp b/include/roofer/io/StreamCropper.hpp index a3fded74..2be3aaa9 100644 --- a/include/roofer/io/StreamCropper.hpp +++ b/include/roofer/io/StreamCropper.hpp @@ -50,7 +50,8 @@ namespace roofer::io { std::vector& polygons, std::vector& buf_polygons, std::vector& point_clouds, vec1f& ground_elevations, - vec1i& acquisition_years, const Box& polygon_extent, + vec1i& acquisition_years, vec1b& pointcloud_insufficient, + const Box& polygon_extent, PointCloudCropperConfig cfg = PointCloudCropperConfig{}) = 0; // virtual void process( // std::string filepaths, std::vector& polygons, diff --git a/include/roofer/misc/PointcloudRasteriser.hpp b/include/roofer/misc/PointcloudRasteriser.hpp index 3340520e..3c0b18ce 100644 --- a/include/roofer/misc/PointcloudRasteriser.hpp +++ b/include/roofer/misc/PointcloudRasteriser.hpp @@ -65,6 +65,8 @@ namespace roofer::misc { float computePointDensity(const ImageMap& image_bundle); + float computeRoofElevation(const ImageMap& image_bundle, float percentile); + // Determine if the two point clouds describe the same object. bool isMutated(const ImageMap& a, const ImageMap& b, const float& threshold_mutation_fraction, diff --git a/include/roofer/reconstruction/PlaneDetector.hpp b/include/roofer/reconstruction/PlaneDetector.hpp index 10678baf..ccc1bed7 100644 --- a/include/roofer/reconstruction/PlaneDetector.hpp +++ b/include/roofer/reconstruction/PlaneDetector.hpp @@ -49,6 +49,11 @@ namespace roofer::reconstruction { bool regularize_orthogonality_ = false; bool regularize_coplanarity_ = false; bool regularize_axis_symmetry_ = false; + + // limits + bool with_limits = false; + int limit_n_regions = 100; + int limit_n_milliseconds = 10000; }; struct PlaneDetectorInterface { diff --git a/include/roofer/reconstruction/RegionGrower.hpp b/include/roofer/reconstruction/RegionGrower.hpp index cb2ef9d1..399d85e5 100644 --- a/include/roofer/reconstruction/RegionGrower.hpp +++ b/include/roofer/reconstruction/RegionGrower.hpp @@ -24,6 +24,7 @@ #include #include #include +#include namespace roofer { @@ -113,6 +114,40 @@ namespace roofer { } } }; + template + void grow_regions_with_limits(candidateDS& cds, Tester& tester, + size_t limit_n_regions, + size_t limit_n_milliseconds) { + std::vector new_regions; + std::deque seeds = cds.get_seeds(); + + region_ids.resize(cds.size, 0); + // first region means unsegmented + regions.push_back(regionType(0)); + + // region growing from seed points + auto t_start = std::chrono::high_resolution_clock::now(); + while (seeds.size() > 0) { + auto idx = seeds.front(); + seeds.erase(seeds.begin()); + if (region_ids[idx] == 0) { + grow_one_region(cds, tester, idx); + ++cur_region_id; + if (regions.size() >= limit_n_regions || + std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - t_start) + .count() >= limit_n_milliseconds) { + throw std::runtime_error( + "Region growing limit reached. Time = " + + std::to_string( + std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - t_start) + .count()) + + "ms, regioncount = " + std::to_string(regions.size())); + } + } + } + }; }; } // namespace regiongrower diff --git a/src/core/reconstruction/PlaneDetector.cpp b/src/core/reconstruction/PlaneDetector.cpp index 2c0dbd0a..ebcbc6ad 100644 --- a/src/core/reconstruction/PlaneDetector.cpp +++ b/src/core/reconstruction/PlaneDetector.cpp @@ -269,7 +269,12 @@ namespace roofer { R; R.min_segment_count = cfg.metrics_plane_min_points; if (points.size() > cfg.metrics_plane_min_points) - R.grow_regions(PDS, DNTester); + if (cfg.with_limits) { + R.grow_regions_with_limits(PDS, DNTester, cfg.limit_n_regions, + cfg.limit_n_milliseconds); + } else { + R.grow_regions(PDS, DNTester); + } total_plane_cnt = R.regions.size(); // classify horizontal/vertical planes using plane normals diff --git a/src/extra/io/StreamCropper.cpp b/src/extra/io/StreamCropper.cpp index 383ce892..a2823932 100644 --- a/src/extra/io/StreamCropper.cpp +++ b/src/extra/io/StreamCropper.cpp @@ -38,6 +38,7 @@ namespace roofer::io { std::vector& point_clouds; vec1f& ground_elevations; vec1i& acquisition_years; + vec1b& pointcloud_insufficient; // ground elevations std::vector> ground_buffer_points; @@ -60,6 +61,7 @@ namespace roofer::io { std::vector& point_clouds, vec1f& ground_elevations, vec1i& acquisition_years, + vec1b& pointcloud_insufficient, const Box& completearea_bb, float& cellsize, float& buffer, int ground_class = 2, int building_class = 6, @@ -71,6 +73,7 @@ namespace roofer::io { ground_class(ground_class), building_class(building_class), acquisition_years(acquisition_years), + pointcloud_insufficient(pointcloud_insufficient), handle_overlap_points(handle_overlap_points) { // point_clouds_ground.resize(polygons.size()); point_clouds.resize(polygons.size()); @@ -204,8 +207,7 @@ namespace roofer::io { void do_post_process(float& ground_percentile, float& max_density_delta, float& coverage_threshold, bool& clear_if_insufficient, vec1f& poly_areas, vec1i& poly_pt_counts_bld, - vec1i& poly_pt_counts_grd, - vec1s& poly_ptcoverage_class, vec1f& poly_densities) { + vec1i& poly_pt_counts_grd, vec1f& poly_densities) { // compute poly properties struct PolyInfo { size_t pt_count_bld; @@ -346,15 +348,8 @@ namespace roofer::io { for (size_t poly_i = 0; poly_i < polygons.size(); ++poly_i) { auto& info = poly_info[poly_i]; - if ((info.pt_count_bld / info.area) < cov_thres) { - if (clear_if_insufficient) { - auto& point_cloud = point_clouds.at(poly_i); - point_cloud.clear(); - } - poly_ptcoverage_class.push_back(std::string("insufficient")); - } else { - poly_ptcoverage_class.push_back(std::string("sufficient")); - } + pointcloud_insufficient.push_back((info.pt_count_bld / info.area) < + cov_thres); // info.pt_count = point_cloud.size(); poly_areas.push_back(float(info.area)); poly_pt_counts_bld.push_back(int(info.pt_count_bld)); @@ -444,27 +439,21 @@ namespace roofer::io { std::vector& buf_polygons, std::vector& point_clouds, vec1f& ground_elevations, vec1i& acquisition_years, - const Box& polygon_extent, PointCloudCropperConfig cfg) { + vec1b& pointcloud_insufficient, const Box& polygon_extent, + PointCloudCropperConfig cfg) { // vec1f ground_elevations; vec1f poly_areas; vec1i poly_pt_counts_bld; vec1i poly_pt_counts_grd; - vec1s poly_ptcoverage_class; vec1f poly_densities; auto& logger = logger::Logger::get_logger(); - PointsInPolygonsCollector pip_collector{polygons, - buf_polygons, - point_clouds, - ground_elevations, - acquisition_years, - polygon_extent, - cfg.cellsize, - cfg.buffer, - cfg.ground_class, - cfg.building_class, - cfg.handle_overlap_points}; + PointsInPolygonsCollector pip_collector{ + polygons, buf_polygons, point_clouds, + ground_elevations, acquisition_years, pointcloud_insufficient, + polygon_extent, cfg.cellsize, cfg.buffer, + cfg.ground_class, cfg.building_class, cfg.handle_overlap_points}; for (auto lasfile : lasfiles) { LASreadOpener lasreadopener; @@ -535,7 +524,7 @@ namespace roofer::io { pip_collector.do_post_process( cfg.ground_percentile, cfg.max_density_delta, cfg.coverage_threshold, cfg.clear_if_insufficient, poly_areas, poly_pt_counts_bld, - poly_pt_counts_grd, poly_ptcoverage_class, poly_densities); + poly_pt_counts_grd, poly_densities); } // void process(std::string source, std::vector& polygons, diff --git a/src/extra/misc/PointcloudRasteriser.cpp b/src/extra/misc/PointcloudRasteriser.cpp index 3f403013..e450884d 100644 --- a/src/extra/misc/PointcloudRasteriser.cpp +++ b/src/extra/misc/PointcloudRasteriser.cpp @@ -223,6 +223,27 @@ namespace roofer::misc { } } + float computeRoofElevation(const ImageMap& pc, float percentile) { + const auto& fp = pc.at("fp").array; + const auto& h_max = pc.at("max").array; + const auto& nodata = pc.at("max").nodataval; + // auto& cellsize = pc.at("fp").cellsize; + std::vector h_max_in_fp; + h_max_in_fp.reserve(fp.size()); + for (size_t i = 0; i < fp.size(); ++i) { + if (fp[i] != 0 && h_max[i] != nodata) { + h_max_in_fp.push_back(h_max[i]); + } + } + if (h_max_in_fp.size() == 0) { + return 0; + } + // sort the values + std::sort(h_max_in_fp.begin(), h_max_in_fp.end()); + // get the 70th percentile + return h_max_in_fp[h_max_in_fp.size() * percentile]; + } + bool testForGlassRoof(const ImageMap& pc, float threshold_glass_roof) { auto& grp = pc.at("grp").array; auto& cellsize = pc.at("grp").cellsize;