diff --git a/cmake/SetupCMakeProject.cmake b/cmake/SetupCMakeProject.cmake index 09245677..b5ce265e 100644 --- a/cmake/SetupCMakeProject.cmake +++ b/cmake/SetupCMakeProject.cmake @@ -37,7 +37,7 @@ function(_alp_build_and_install NAME SRC_DIR BUILD_DIR INSTALL_DIR) "-DCMAKE_EXE_LINKER_FLAGS=${CMAKE_EXE_LINKER_FLAGS} ${_alp_sanitizer_flags}" "-DCMAKE_MODULE_LINKER_FLAGS=${CMAKE_MODULE_LINKER_FLAGS} ${_alp_sanitizer_flags}" "-DCMAKE_SHARED_LINKER_FLAGS=${CMAKE_SHARED_LINKER_FLAGS} ${_alp_sanitizer_flags}" - -DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH} + "-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH}" -DCMAKE_INSTALL_PREFIX=${INSTALL_DIR} -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} ${ARGN} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e8d11f46..a6839538 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -16,7 +16,7 @@ find_package(TBB REQUIRED) include(../cmake/SetupGDAL.cmake) alp_setup_gdal(GDAL_VERSION v3.10.3 PROJ_VERSION 9.6.0) -alp_setup_cmake_project(cgal URL https://github.com/CGAL/cgal.git COMMITISH "v6.0.1" CMAKE_ARGUMENTS -DCGAL_HEADER_ONLY=ON -DBUILD_TESTING=OFF -DBUILD_EXAMPLES=OFF) +alp_setup_cmake_project(cgal URL https://github.com/CGAL/cgal.git COMMITISH "v6.0.1" CMAKE_ARGUMENTS -DBUILD_TESTING=OFF -DWITH_examples=OFF) find_package(CGAL REQUIRED) include(../cmake/SetupOpenCV.cmake) diff --git a/src/dag_builder/ContinuationMode.h b/src/dag_builder/ContinuationMode.h index c27f1f98..46a707ab 100644 --- a/src/dag_builder/ContinuationMode.h +++ b/src/dag_builder/ContinuationMode.h @@ -1,7 +1,7 @@ #pragma once enum class ContinuationMode { + Error, Resume, Overwrite, - Fail }; diff --git a/src/dag_builder/Partitioning.h b/src/dag_builder/Partitioning.h new file mode 100644 index 00000000..67c74893 --- /dev/null +++ b/src/dag_builder/Partitioning.h @@ -0,0 +1,10 @@ +#pragma once + +#include +#include + + +struct [[nodiscard]] Partitioning { + uint32_t partition_count = 0; + std::vector cluster_partitions; // original cluster index -> partition index +}; diff --git a/src/dag_builder/VertexInCluster.h b/src/dag_builder/VertexInCluster.h new file mode 100644 index 00000000..cafdccfa --- /dev/null +++ b/src/dag_builder/VertexInCluster.h @@ -0,0 +1,21 @@ +#pragma once + +#include + +#include "hash_utils.h" + +struct VertexInCluster { + uint32_t cluster_index; + uint32_t local_vertex_index; + + auto operator<=>(const VertexInCluster &other) const = default; +}; + +namespace std { + template <> + struct hash { + size_t operator()(const VertexInCluster &v) const noexcept { + return ::hash::combine(v.cluster_index, v.local_vertex_index); + } + }; +} diff --git a/src/dag_builder/VertexInClustering.h b/src/dag_builder/VertexInClustering.h new file mode 100644 index 00000000..4600141c --- /dev/null +++ b/src/dag_builder/VertexInClustering.h @@ -0,0 +1,21 @@ +#pragma once + +#include + +#include "hash_utils.h" + +struct VertexInClustering { + uint32_t clustering_index; + uint32_t global_vertex_index; + + auto operator<=>(const VertexInClustering &other) const = default; +}; + +namespace std { +template <> +struct hash { + size_t operator()(const VertexInClustering &v) const noexcept { + return ::hash::combine(v.clustering_index, v.global_vertex_index); + } +}; +} // namespace std diff --git a/src/dag_builder/atlas/Packer.cpp b/src/dag_builder/atlas/Packer.cpp index 72be5c87..4bfb7cde 100644 --- a/src/dag_builder/atlas/Packer.cpp +++ b/src/dag_builder/atlas/Packer.cpp @@ -4,6 +4,7 @@ #include #include +#include #include "log.h" #include "Packer.h" @@ -37,10 +38,13 @@ inline void normalize_uvs(std::span uvs) { struct Packer::Impl { std::unique_ptr atlas; + double total_effective_pixel_area = 0.0; Impl() : atlas(xatlas::Create(), xatlas::Destroy) {} void add_uv_mesh(const UvMesh& mesh) { + this->total_effective_pixel_area += glm::compMul(glm::dvec2(mesh.texture_size)); + std::vector scaled_uvs = transform_vector(mesh.uvs, [&](const glm::dvec2 &uv) { const glm::dvec2 scaled = uv * glm::dvec2(mesh.texture_size); return glm::vec2(scaled); @@ -87,9 +91,7 @@ struct Packer::Impl { normalize_uvs(uvs.flat()); - // TODO: return atlas->width and atlas->height or aspect? - - return Packing(uvs); + return Packing(uvs, this->atlas->utilization[0]); } }; diff --git a/src/dag_builder/atlas/Packer.h b/src/dag_builder/atlas/Packer.h index 450f8856..1e47a55a 100644 --- a/src/dag_builder/atlas/Packer.h +++ b/src/dag_builder/atlas/Packer.h @@ -1,9 +1,13 @@ #pragma once +#include +#include +#include #include #include +#include "number_utils.h" #include "SegmentedBuffer.h" namespace atlas { @@ -17,8 +21,10 @@ struct UvMesh { class Packing { public: Packing() = default; - Packing(std::vector uvs) : Packing(SegmentedBuffer(std::move(uvs))) {} - Packing(SegmentedBuffer uvs) : _uvs(std::move(uvs)) {} + Packing(std::vector uvs, const double effective_pixel_area, const float utilization = 1.0f) + : Packing(SegmentedBuffer(std::move(uvs)), effective_pixel_area, utilization) {} + Packing(SegmentedBuffer uvs, const double effective_pixel_area, const float utilization = 1.0f) + : _uvs(std::move(uvs)), _effective_pixel_area(effective_pixel_area), _utilization(utilization) {} std::span uvs_for_mesh(const uint32_t mesh_index) const { return this->_uvs.segment(mesh_index); @@ -31,8 +37,22 @@ class Packing { return this->_uvs.backing(); } + // Sum of effective input pixels (source-texture pixels actually covered by UVs) across all packed meshes. + double effective_pixel_area() const { + return this->_effective_pixel_area; + } + + // Fraction of the packed atlas actually covered by charts. 1.0 for + // packings that were not produced by the chart packer (e.g. a single + // mesh's own UVs), since there is no wasted atlas space to account for. + float utilization() const { + return this->_utilization; + } + private: SegmentedBuffer _uvs; + double _effective_pixel_area = 0.0; + float _utilization = 1.0f; }; class Packer { diff --git a/src/dag_builder/atlas/TextureBaker.h b/src/dag_builder/atlas/TextureBaker.h index e5f242fb..fd31d9c9 100644 --- a/src/dag_builder/atlas/TextureBaker.h +++ b/src/dag_builder/atlas/TextureBaker.h @@ -1,5 +1,6 @@ #pragma once +#include #include #include #include @@ -12,10 +13,12 @@ #include #include +#include #include "atlas/Packer.h" #include "atlas/pull_reproject_texture.h" #include "mesh/bounds.h" +#include "number_utils.h" #include "opencv_utils.h" #include "range_utils.h" #include "TextureSet.h" @@ -84,13 +87,25 @@ inline glm::uvec2 get_effective_texture_size(const TextureMap &map) { inline glm::uvec2 get_effective_texture_size(const TextureComposition &comp) { double area = 0; for (const auto &map : comp.maps) { - const auto size = get_effective_texture_size(map); - area += size.x * size.y; + area += glm::compMul(glm::dvec2(get_effective_texture_size(map))); } return glm::uvec2(std::ceil(std::sqrt(area))); } } // namespace detail +// Resolution needed to give a compression_ratio fraction of effective_pixel_area +// worth of resolution, after accounting for atlas packing waste (unused texels +// between packed charts), capped at max_texture_size. +inline glm::uvec2 compute_bake_texture_size(const double effective_pixel_area, const float utilization, const double compression_ratio, const uint32_t max_texture_size) { + const double target_area = effective_pixel_area * compression_ratio / std::max(utilization, 1e-3); + const uint32_t side = static_cast(std::ceil(std::sqrt(target_area))); + return glm::uvec2(std::min(next_power_of_two(side), max_texture_size)); +} + +inline glm::uvec2 compute_bake_texture_size(const atlas::Packing &packing, const double compression_ratio, const uint32_t max_texture_size) { + return compute_bake_texture_size(packing.effective_pixel_area(), packing.utilization(), compression_ratio, max_texture_size); +} + class TextureBaker { public: TextureMapId add_mesh(TexturedMesh mesh) { @@ -128,36 +143,72 @@ class TextureBaker { return id; } - TextureCake bake(const glm::uvec2 texture_size) { + // Exact output resolution. + TextureCake bake(const glm::uvec2 texture_size) const { + return this->bake(this->pack(), texture_size); + } + + // Runs atlas packing (chart computation for multi-mesh atlases; a cheap + // passthrough for a single entry). Callers that need packing.utilization() + // or packing.effective_pixel_area() to choose an output resolution should + // call this once and pass the result to bake(packing, texture_size) below, + // rather than letting bake(texture_size) compute it again. + atlas::Packing pack() const { if (this->_entries.size() == 1) { - return this->bake_single(this->_entries.front(), texture_size); - } else { - return this->bake_multi(texture_size); + return this->pack_single(this->_entries.front()); + } + + atlas::Packer packer; + + std::vector target_triangles; + for (const Entry& entry : this->_entries) { + visit_entry(entry, + [&](const TexturedMesh& mesh) { + packer.add_uv_mesh(mesh.triangles, mesh.map.uvs, detail::get_effective_texture_size(mesh.map)); + }, + [&](const TextureComposition& comp) { + target_triangles.clear(); + target_triangles.reserve(comp.triangles.size()); + for (const MappedTriangle &triangle : comp.triangles) { + target_triangles.push_back(triangle.target); + } + packer.add_uv_mesh( + target_triangles, + comp.target_uvs, + detail::get_effective_texture_size(comp) + ); + }); } + + return packer.pack(); } -private: - using Entry = std::variant; + // Bakes using an already-computed packing (e.g. from pack()), so no atlas repacking happens. + TextureCake bake(const atlas::Packing &packing, const glm::uvec2 texture_size) const { + if (this->_entries.size() == 1) { + return visit_entry(this->_entries.front(), + [&](const TexturedMesh &mesh) { + TextureCake cake; + cake._packing = packing; + cake._texture = rescale_texture(mesh.map.texture, texture_size); + return cake; + }, + [&](const TextureComposition &) { + TextureCake cake; + cake._packing = packing; + cake._texture = this->bake_texture(cake._packing, texture_size); + return cake; + }); + } - TextureCake bake_multi(const glm::uvec2 texture_size) const { TextureCake cake; - cake._packing = this->pack(); + cake._packing = packing; cake._texture = this->bake_texture(cake._packing, texture_size); return cake; } - TextureCake bake_single(const Entry &entry, const glm::uvec2 texture_size) const { - return visit_entry(entry, - [&](const TexturedMesh &mesh) { - TextureCake cake; - cake._packing = this->pack_single(mesh); - cake._texture = rescale_texture(mesh.map.texture, texture_size); - return cake; - }, - [&](const TextureComposition &) { - return this->bake_multi(texture_size); - }); - } +private: + using Entry = std::variant; Texture bake_texture(const atlas::Packing &packing, const glm::uvec2 texture_size) const { TextureReprojector reprojector(texture_size, CV_8UC3); @@ -233,43 +284,24 @@ class TextureBaker { } atlas::Packing pack_single(const Entry &entry) const { + const double area = effective_pixel_area(entry); return visit_entry(entry, [&](const TexturedMesh &mesh) { - return atlas::Packing(mesh.map.uvs); + return atlas::Packing(mesh.map.uvs, area); }, [&](const TextureComposition &comp) { - return atlas::Packing(comp.target_uvs); + return atlas::Packing(comp.target_uvs, area); }); } - atlas::Packing pack() const { - if (this->_entries.size() == 1) { - return this->pack_single(this->_entries.front()); - } - - atlas::Packer packer; - - std::vector target_triangles; - for (const Entry& entry : this->_entries) { - visit_entry(entry, - [&](const TexturedMesh& mesh) { - packer.add_uv_mesh(mesh.triangles, mesh.map.uvs, detail::get_effective_texture_size(mesh.map)); - }, - [&](const TextureComposition& comp) { - target_triangles.clear(); - target_triangles.reserve(comp.triangles.size()); - for (const MappedTriangle &triangle : comp.triangles) { - target_triangles.push_back(triangle.target); - } - packer.add_uv_mesh( - target_triangles, - comp.target_uvs, - detail::get_effective_texture_size(comp) - ); - }); - } - - return packer.pack(); + static double effective_pixel_area(const Entry &entry) { + return visit_entry(entry, + [](const TexturedMesh &mesh) { + return glm::compMul(glm::dvec2(detail::get_effective_texture_size(mesh.map))); + }, + [](const TextureComposition &comp) { + return glm::compMul(glm::dvec2(detail::get_effective_texture_size(comp))); + }); } template diff --git a/src/dag_builder/atlas/pull_reproject_texture.h b/src/dag_builder/atlas/pull_reproject_texture.h index a3ec847e..749ce8f2 100644 --- a/src/dag_builder/atlas/pull_reproject_texture.h +++ b/src/dag_builder/atlas/pull_reproject_texture.h @@ -47,11 +47,11 @@ glm::dvec2 clamp_uv(const glm::dvec2 &uv) { } FixedPoint sample_cell_center(const glm::uvec2 sample) { - return FixedPoint::from_real(glm::dvec2(sample) + 0.5); -} - -bool is_left_of_or_on_edge(const FixedPoint edge_start, const FixedPoint edge_end, const FixedPoint sample) { - return fp::orient(edge_start, edge_end, sample) >= FixedScalar::zero(); + // Fixed-point encoding of sample + 0.5; exactly representable, so no rounding happens. + const FixedCoord::repr_t half = FixedCoord::scale / 2; + return FixedPoint{{ + FixedCoord::from_integer(FixedCoord::repr_t{sample.x}).value + half, + FixedCoord::from_integer(FixedCoord::repr_t{sample.y}).value + half}}; } radix::geometry::Aabb2ui clamp_bounds( @@ -79,23 +79,51 @@ radix::geometry::Aabb2ui get_pixel_bounds( return clamp_bounds(pixel_bounds, grid_bounds); } +// Edge function of a directed edge, evaluated incrementally over the sample grid. +// The value is the raw cross product (edge_end - edge_start) x (sample - edge_start) +// without rescaling, so its sign is the exact orientation of the sample. +struct EdgeFunction { + FixedCoord::repr_t value = 0; + FixedCoord::repr_t step_x = 0; + FixedCoord::repr_t step_y = 0; + + EdgeFunction(const FixedPoint edge_start, const FixedPoint edge_end, const FixedPoint sample) { + const auto edge = edge_end.raw - edge_start.raw; + const auto offset = sample.raw - edge_start.raw; + this->value = edge.x * offset.y - edge.y * offset.x; + this->step_x = -edge.y * FixedCoord::scale; + this->step_y = edge.x * FixedCoord::scale; + } +}; + void rasterize_triangle( const FixedTriangle &triangle, const int32_t triangle_id, const glm::uvec2 grid_size, Vector2D &triangle_index_map) { const radix::geometry::Aabb2ui bounds = get_pixel_bounds(triangle, grid_size); + const FixedPoint start = sample_cell_center({bounds.min.x, bounds.min.y}); + + std::array edges = { + EdgeFunction(triangle[0], triangle[1], start), + EdgeFunction(triangle[1], triangle[2], start), + EdgeFunction(triangle[2], triangle[0], start)}; + for (uint32_t sample_y = bounds.min.y; sample_y <= bounds.max.y; sample_y++) { + std::array row = {edges[0].value, edges[1].value, edges[2].value}; + for (uint32_t sample_x = bounds.min.x; sample_x <= bounds.max.x; sample_x++) { - const FixedPoint sample = sample_cell_center({sample_x, sample_y}); + if (row[0] >= 0 && row[1] >= 0 && row[2] >= 0) { + triangle_index_map(sample_y, sample_x) = triangle_id; + } - if (!is_left_of_or_on_edge(triangle[0], triangle[1], sample) || - !is_left_of_or_on_edge(triangle[1], triangle[2], sample) || - !is_left_of_or_on_edge(triangle[2], triangle[0], sample)) { - continue; + for (uint32_t k = 0; k < 3; k++) { + row[k] += edges[k].step_x; } + } - triangle_index_map(sample_y, sample_x) = triangle_id; + for (uint32_t k = 0; k < 3; k++) { + edges[k].value += edges[k].step_y; } } } @@ -182,7 +210,7 @@ glm::dvec2 target_to_source(const PreparedTriangle &triangle, const FixedPoint t // Clamp-to-edge: out-of-bounds coordinates reuse the nearest border pixel. cv::Vec3f source_texel(const cv::Mat &image, const glm::ivec2 pos) { const glm::ivec2 clamped = glm::clamp(pos, glm::ivec2(0), glm::ivec2(image.cols - 1, image.rows - 1)); - return cv::Vec3f(image.at(clamped.y, clamped.x)) * (1.0f / 255.0f); + return cv::Vec3f(image.at(clamped.y, clamped.x)); } cv::Vec3f sample_source_linear( @@ -245,8 +273,51 @@ class TextureReprojector { const Vector2D triangle_index_map = build_triangle_index_map(prepared_triangles, grid_size); - cv::Mat output(this->_output_size.y, this->_output_size.x, CV_32FC3, cv::Scalar(0, 0, 0)); + cv::Mat output(this->_output_size.y, this->_output_size.x, this->_output_type, cv::Scalar::all(0)); + this->render_by_depth(output, prepared_triangles, triangle_index_map); + return output; + } + +private: + // Resolve the destination channel type once so the per-pixel write is a + // compile-time cv::saturate_cast, matching convertTo's rounding without its + // per-pixel dispatch or a full-image float accumulator. + void render_by_depth( + cv::Mat &output, + const std::vector &prepared_triangles, + const Vector2D &triangle_index_map) { + switch (CV_MAT_DEPTH(this->_output_type)) { + case CV_8U: + this->render_into(output, prepared_triangles, triangle_index_map); + break; + case CV_8S: + this->render_into(output, prepared_triangles, triangle_index_map); + break; + case CV_16U: + this->render_into(output, prepared_triangles, triangle_index_map); + break; + case CV_16S: + this->render_into(output, prepared_triangles, triangle_index_map); + break; + case CV_32S: + this->render_into(output, prepared_triangles, triangle_index_map); + break; + case CV_32F: + this->render_into(output, prepared_triangles, triangle_index_map); + break; + case CV_64F: + this->render_into(output, prepared_triangles, triangle_index_map); + break; + default: + throw std::invalid_argument("TextureReprojector: unsupported output depth"); + } + } + template + void render_into( + cv::Mat &output, + const std::vector &prepared_triangles, + const Vector2D &triangle_index_map) { for (uint32_t y = 0; y < this->_output_size.y; y++) { for (uint32_t x = 0; x < this->_output_size.x; x++) { cv::Vec3f color(0, 0, 0); @@ -279,17 +350,16 @@ class TextureReprojector { // Dividing by sample_count (not sample_scale^2) correctly averages partially covered pixels. if (sample_count > 0) { - output.at(y, x) = color / static_cast(sample_count); + const cv::Vec3f averaged = color / static_cast(sample_count); + output.at>(y, x) = cv::Vec( + cv::saturate_cast(averaged[0]), + cv::saturate_cast(averaged[1]), + cv::saturate_cast(averaged[2])); } } } - - cv::Mat converted; - output.convertTo(converted, this->_output_type); - return converted; } -private: glm::uvec2 _output_size; int _output_type; ReprojectionOptions _options; diff --git a/src/dag_builder/build.cpp b/src/dag_builder/build.cpp index bcb833b6..aa8f6b45 100644 --- a/src/dag_builder/build.cpp +++ b/src/dag_builder/build.cpp @@ -15,7 +15,6 @@ #include "clusterize.h" #include "compact.h" #include "encoded.h" -#include "error_bounds.h" #include "int_math.h" #include "log.h" #include "merge/clusterings.h" @@ -38,6 +37,7 @@ #include "utils.h" #include "vertex_lock.h" #include "parallel.h" +#include "ContinuationMode.h" namespace dag { @@ -74,11 +74,11 @@ struct LodResult { }; // Run the full LOD pipeline on a clustering: partition, simplify, re-clusterize, and build child map. -LodResult build_lod(const Clustering &input, const BuildOptions &options) { - const Partitioning partitioning = create_partitioning(input, PartitionOptions{ - .clusters_per_partition = options.clusters_per_partition}); +LodResult build_lod(const Clustering &clusters, const BuildOptions &options, const radix::geometry::Aabb3d &node_bounds) { + const Partitioning partitioning = create_partitioning(clusters, PartitionOptions{ + .clusters_per_partition = options.clusters_per_partition}); // Construct new clusters and generate new textures. - Clustering clustering = apply_partitioning(input, partitioning); + Clustering clustering = apply_partitioning(clusters, partitioning, options.uv_unwrap_algorithm); // Create partition to cluster map std::vector> partition_to_clusters(partitioning.partition_count); for (const auto [cluster_index, partition_index] : enumerate(partitioning.cluster_partitions)) { @@ -89,29 +89,36 @@ LodResult build_lod(const Clustering &input, const BuildOptions &options) { trim_textures_inplace(clustering); // Find vertices to lock - const std::vector vertex_lock = find_vertices_to_lock(clustering); - - // Simplify each cluster - clustering = simplify(clustering, SimplifyOptions{ - .target_ratio = options.target_ratio, - .vertex_lock = VertexLock::mask(vertex_lock)}); + const std::vector vertex_lock = find_vertices_to_lock(clustering, node_bounds); + + // Convert relative_target_error (a fraction of the node bounds) to an absolute + // error; both targets stay optional so simplify() can apply its own default. + std::optional absolute_target_error; + if (options.relative_target_error) { + const double node_extent = glm::compMax(node_bounds.size()); + absolute_target_error = static_cast(options.relative_target_error.value() * node_extent); + } + const SimplifyOptions simplify_options{ + .target_ratio = options.target_ratio, + .absolute_target_error = absolute_target_error, + .vertex_lock = VertexLock::mask(vertex_lock), + .error_mode = ErrorMode::Add, + .preserve_cluster_count = true + }; + clustering = simplify(clustering, simplify_options); remove_unused_vertices_inplace(clustering); // Split each cluster into roughly 4 parts auto result = clusterize(clustering); - clustering = result.clustering; - // Build child map - std::vector> child_map(clustering.cluster_count()); + // Build child map. Simplification preserves the cluster count, so each final + // cluster's backward mapping index is its partition index. + std::vector> child_map(result.clustering.cluster_count()); for (const auto [final_index, partition_index] : enumerate(result.backward_mapping)) { - auto &children = child_map[final_index]; - const auto &original_clusters = partition_to_clusters[partition_index]; - for (const uint32_t original_index : original_clusters) { - children.push_back(original_index); - } + child_map[final_index] = partition_to_clusters[partition_index]; } - return {clustering, child_map}; + return {std::move(result.clustering), std::move(child_map)}; } // Load input meshes, clusterize, and filter them to the target region. @@ -157,17 +164,21 @@ dag::ClusterBatch load_and_simplify_dag_nodes( const std::vector &dag_ids, const RegionFilter &filter, const double epsilon, + const radix::geometry::Aabb3d &node_bounds, const BuildContext &ctx) { std::vector cluster_sources; std::vector filtered; for (const octree::Id &id : dag_ids) { - const auto dag_node = ctx.output_storage.load(id); + // Load dag node clusters + auto dag_node = ctx.output_storage.load(id); if (!dag_node) { LOG_WARN("Failed to load DAG node {}, skipping", id); continue; } Clustering clustering = std::move(dag_node.value().clustering); + + // Assign canonical cluster ids for (auto &[cluster_index, cluster] : enumerate(clustering.clusters)) { cluster.id = cluster_sources.size(); cluster_sources.emplace_back(id, cluster_index); @@ -177,16 +188,23 @@ dag::ClusterBatch load_and_simplify_dag_nodes( } auto indices = find_clusters_matching(clustering, filter); + // Find relevant clusters from this node. if (indices.empty()) { continue; } filtered.push_back(slice_clusters(clustering, indices)); } + if (filtered.empty()) { + return {}; + } + // Merge remaining clusterings const Clustering merged = merge_clusterings(filtered, epsilon); - const auto [simplified, child_map] = build_lod(merged, ctx.options); + auto [simplified, child_map] = build_lod(merged, ctx.options, node_bounds); + // Run nanite-style grouping, simplification, splitting. + // Create map from cluster index to child cluster id. auto child_id_map = transform_vector(child_map, [&](const auto &children) { return transform_vector(children, [&](const uint32_t merged_index) { const uint32_t source_index = merged.clusters[merged_index].id; @@ -194,7 +212,7 @@ dag::ClusterBatch load_and_simplify_dag_nodes( }); }); - return {simplified, std::move(child_id_map)}; + return {std::move(simplified), std::move(child_id_map)}; } // Combine input clusters with inner clusters. @@ -239,7 +257,7 @@ dag::ClusterBatch combine_input_and_inner( // Compute epsilon value for merging clusters based on the size of the node bounds. double compute_epsilon(const radix::geometry::Aabb3d &bounds) { - return glm::compAdd(bounds.size()) / 30000.0; + return glm::compAdd(bounds.size()) / 10000; } // Build a single target node by loading input clusters (preserved as-is) and @@ -254,6 +272,8 @@ std::optional build_node( const auto node_bounds = ctx.shifted_space.get_node_bounds(target_id); const double epsilon = compute_epsilon(node_bounds); + // Prepare filter to only include clusters inside the target_id bounds. + // In CurrentAndCoarser mode, input regions represented by the relevant DAG nodes are excluded as well. RegionFilter input_filter; input_filter.include = {node_bounds}; if (ctx.options.include_mode == IncludeMode::CurrentAndCoarser) { @@ -265,8 +285,7 @@ std::optional build_node( RegionFilter dag_filter; dag_filter.include = {node_bounds}; - dag::ClusterBatch inner = load_and_simplify_dag_nodes( - dag_ids, dag_filter, epsilon, ctx); + dag::ClusterBatch inner = load_and_simplify_dag_nodes(dag_ids, dag_filter, epsilon, node_bounds, ctx); if (input_clusters.empty() && inner.clustering.is_empty()) { LOG_WARN("No valid clusters for node {}, skipping", target_id); @@ -313,7 +332,7 @@ std::vector find_relevant_input_nodes( } else if (include_mode == IncludeMode::CurrentAndCoarser) { // Include inputs at same or coarser levels whose bounds intersect target_id in the shifted tree. const radix::geometry::Aabb3d target_bounds = shifted_space.get_node_bounds(target_id); - for (uint32_t level = target_id.level(); level < input_by_level.size(); level++) { + for (uint32_t level = 0; level <= target_id.level(); level++) { const std::vector &input_nodes = input_by_level[level]; for (const octree::Id& input_id : input_nodes) { const radix::geometry::Aabb3d input_bounds = space.get_node_bounds(input_id); @@ -332,9 +351,9 @@ std::vector find_relevant_dag_nodes( const std::unordered_set &prev_level_built, const octree::OddLevelShifted &shifted_space) { std::vector result; - for (const octree::Id &parent : shifted_space.get_intersecting_nodes_on_level(target_id, target_id.level()+1)) { - if (prev_level_built.contains(parent)) { - result.push_back(parent); + for (const octree::Id &child : shifted_space.get_intersecting_nodes_on_level(target_id, target_id.level()+1)) { + if (prev_level_built.contains(child)) { + result.push_back(child); } } return result; @@ -351,7 +370,10 @@ std::unordered_set find_nodes_to_build_on_level( // Consider input nodes at this level const std::span level_input_ids = input_by_level[level]; for (const octree::Id &input_id : level_input_ids) { - target_set.insert(input_id); + const octree::IdRect shifted_nodes = ctx.shifted_space.find_intersecting_nodes_for_standard_id(input_id); + for (const octree::Id &shifted_id : shifted_nodes) { + target_set.insert(shifted_id); + } } // Consider parents of previously built nodes @@ -438,7 +460,7 @@ std::unordered_set build_level( LOG_INFO("Building level {} ({} targets)", level, targets.size()); std::unordered_set already_built; - if (ctx.options.resume) { + if (ctx.options.continuation_mode != ContinuationMode::Overwrite) { for (const octree::Id &target : targets) { if (ctx.output_storage.has(target)) { already_built.insert(target); @@ -446,6 +468,10 @@ std::unordered_set build_level( } } + if (ctx.options.continuation_mode == ContinuationMode::Error && already_built.empty()) { + LOG_ERROR_AND_EXIT("Found some of target nodes already in built in the output directory, use --resume or --overwrite"); + } + // Initialize debug storage if requested (contains .glb meshes) std::optional debug_storage; if (ctx.options.write_debug_meshes) { @@ -458,33 +484,28 @@ std::unordered_set build_level( tbb::concurrent_vector saved_ids; ProgressIndicator progress(targets.size()); - for (size_t i = 0; i < already_built.size(); i++) { - progress.task_finished(); - } auto progress_thread = progress.start_monitoring(); parallel_foreach(targets, [&](const octree::Id &target) { if (already_built.contains(target)) { + progress.task_finished(); return; } const auto dag_ids = find_value(inner_nodes, target).value_or(std::vector{}); const auto target_input_ids = find_value(input_sources, target).value_or(std::vector{}); - auto result = build_node( - target, - target_input_ids, - dag_ids, - ctx); - + auto result = build_node(target, target_input_ids, dag_ids, ctx); if (result) { const auto save_result = ctx.output_storage.save(target, *result); - DEBUG_ASSERT(save_result.has_value()); - ALP_UNUSED(save_result); - if (debug_storage) { - debug_storage->save(target, clustering_to_mesh(result->clustering)); + if (save_result) { + if (debug_storage) { + debug_storage->save(target, clustering_to_mesh(result->clustering)); + } + saved_ids.push_back(target); + } else { + LOG_ERROR("Failed to save node {}: {}", target, save_result.error()); } - saved_ids.push_back(target); } progress.task_finished(); }, ctx.options.parallelize); @@ -506,7 +527,7 @@ void build_levels( const octree::IndexedMeshStorage &input_storage, octree::IndexedDagStorage &output_storage, const BuildOptions &options, - const Range &level_range) { + const AnyRange &level_range) { const octree::OddLevelShifted shifted_space = octree::OddLevelShifted::earth(); const octree::Space space = octree::Space::earth(); const octree::Id root_node = options.root_node; @@ -522,16 +543,24 @@ void build_levels( const uint32_t max_input_level = max_input_level_opt.value(); const Range valid_range{root_node.level(), max_input_level + 1}; - const Range range = valid_range.intersection(level_range); - if (range.empty()) { + const Range range = valid_range.intersect(level_range).to_range(max_input_level + 1); + if (range.is_empty()) { LOG_WARN("Requested level range does not overlap with buildable levels {}-{}", root_node.level(), max_input_level); return; } + // Seed prev_level_built with any already-built nodes one level finer than the first level. + // TODO: use hierachical lookup here based on root_bounds and IdRect + std::unordered_set prev_level_built; + for (const auto &[id, status] : output_storage.index()) { + if (id.level() == range.end && status != octree::NodeStatus::Virtual) { + prev_level_built.insert(id); + } + } + BuildContext ctx{input_storage, ThreadSafeStorage(std::move(output_storage)), options, space, shifted_space, root_bounds}; - std::unordered_set prev_level_built; - for (uint32_t level = range.max; level-- > range.min;) { + for (uint32_t level = range.end; level-- > range.start;) { prev_level_built = build_level( level, input_by_level, @@ -547,7 +576,7 @@ void build_full( const octree::IndexedMeshStorage &input_storage, octree::IndexedDagStorage &output_storage, const BuildOptions &options) { - build_levels(input_storage, output_storage, options, full_range()); + build_levels(input_storage, output_storage, options, RangeFull{}); } } // namespace dag diff --git a/src/dag_builder/build.h b/src/dag_builder/build.h index dfd407e6..9d000e58 100644 --- a/src/dag_builder/build.h +++ b/src/dag_builder/build.h @@ -1,7 +1,9 @@ #pragma once #include +#include +#include "ContinuationMode.h" #include "Range.h" #include "build_config.h" #include "octree/storage/MeshStorage.h" @@ -18,13 +20,14 @@ enum class IncludeMode { struct BuildOptions { uint32_t clusters_per_partition; - float target_ratio; + std::optional target_ratio; + std::optional relative_target_error; uv::Algorithm uv_unwrap_algorithm; octree::Id root_node = octree::Id::root(); IncludeMode include_mode = IncludeMode::CurrentOnly; bool write_debug_meshes = IS_DEBUG_BUILD; bool parallelize = false; - bool resume = true; + ContinuationMode continuation_mode = ContinuationMode::Error; }; void build_full( @@ -36,6 +39,6 @@ void build_levels( const octree::IndexedMeshStorage &input_storage, octree::IndexedDagStorage &output_storage, const BuildOptions &options, - const Range &level_range); + const AnyRange &level_range); } // namespace dag diff --git a/src/dag_builder/cli.cpp b/src/dag_builder/cli.cpp index 34aa1eba..4d0030e2 100644 --- a/src/dag_builder/cli.cpp +++ b/src/dag_builder/cli.cpp @@ -7,9 +7,10 @@ #include #include #include -#include #include +#include #include +#include #include using namespace cli; @@ -26,6 +27,11 @@ const std::unordered_map log_level_names {"trace", spdlog::level::level_enum::trace}, }; +const std::unordered_map include_mode_names{ + {"current", dag::IncludeMode::CurrentOnly}, + {"current-and-coarser", dag::IncludeMode::CurrentAndCoarser}, +}; + const std::unordered_map uv_unwrap_algorithm_names{ {"AsRigidAsPossible", uv::Algorithm::AsRigidAsPossible}, {"arap", uv::Algorithm::AsRigidAsPossible}, @@ -45,19 +51,19 @@ const std::unordered_map uv_unwrap_algorithm_names{ {"TutteBarycentricMapping", uv::Algorithm::TutteBarycentricMapping}, {"tutte", uv::Algorithm::TutteBarycentricMapping}}; -Range make_level_range(const std::vector& input) { +AnyRange make_level_range(const std::vector& input) { switch (input.size()) { case 0: - return full_range(); + return RangeFull{}; case 1: - return Range(input[0]); + return RangeInclusive(input[0]); case 2: { const uint32_t min_level = input[0]; const uint32_t max_level = input[1]; if (min_level > max_level) { throw CLI::ValidationError("--levels minimum must be <= maximum"); } - return Range{min_level, max_level}; + return RangeInclusive{min_level, max_level}; } default: UNREACHABLE(); @@ -84,6 +90,34 @@ octree::Id make_octree_id(const std::vector &values) { throw CLI::ValidationError("--root-node expects either: or "); } +ContinuationMode make_continuation_mode(const bool resume, const bool overwrite) { + if (resume && overwrite) { + UNREACHABLE(); + } + if (resume) { + return ContinuationMode::Resume; + } else if (overwrite) { + return ContinuationMode::Overwrite; + } + + // TODO(ORIGINAL AUTHOR): !!! THIS FALLBACK NEEDS AN AUTHORITATIVE BEHAVIOUR DECISION !!! + // + // Neither --resume nor --overwrite is a normal command-line combination, but the intended + // continuation semantics are unclear. Do not replace this exception without first confirming + // the desired behaviour with the original author. + // + // Codex believes the correct implementation is: + // + // return ContinuationMode::Error; + // + // ContinuationMode::Error is already the default in Args and BuildOptions, and build_level() + // contains a dedicated branch for it. Returning Error would therefore preserve the apparent + // design: refuse to overwrite existing output unless the user explicitly selects --resume or + // --overwrite. The exception below intentionally keeps the unsupported path loud until the + // original author confirms or corrects that interpretation. + throw std::runtime_error("unsupported operation"); +} + } // namespace Args cli::parse(int argc, const char *const *argv) { @@ -97,12 +131,15 @@ Args cli::parse(int argc, const char *const *argv) { .input_path = {}, .output_path = {}, .root_node = octree::Id::root(), - .level_range = full_range(), + .level_range = RangeFull{}, .uv_unwrap_algorithm = {}, .clusters_per_partition = 8, - .target_ratio = 0.25, + .target_ratio = std::nullopt, + .target_error = std::nullopt, .write_debug_meshes = false, - .overwrite = false, + .parallelize = false, + .include_mode = dag::IncludeMode::CurrentOnly, + .continuation_mode = ContinuationMode::Error }; app.add_option("--input", args.input_path, "Path to input mesh dataset") @@ -133,10 +170,25 @@ Args cli::parse(int argc, const char *const *argv) { ->check(CLI::PositiveNumber); app.add_option("--target-ratio", args.target_ratio, "Simplification target ratio") - ->default_val(args.target_ratio) ->check(CLI::Range(0.0, 1.0)); - app.add_flag("--overwrite", args.overwrite, "Overwrite data already present in output"); + app.add_option("--target-error", args.target_error, "Simplification target error as a fraction of node bounds") + ->check(CLI::NonNegativeNumber); + + bool resume = false; + bool overwrite = false; + app.add_flag("--resume", resume, "Resume building the dag from then data in output") + ->excludes("--overwrite"); + app.add_flag("--overwrite", overwrite, "Overwrite data already present in output") + ->excludes("--resume"); + + app.add_flag("--write-debug-meshes", args.write_debug_meshes, "Write debug .glb meshes alongside the output"); + + app.add_flag("--parallelize", args.parallelize, "Build nodes within a level in parallel"); + + app.add_option("--include-mode", args.include_mode, "Which input nodes to include when building a level") + ->transform(CLI::CheckedTransformer(include_mode_names, CLI::ignore_case)) + ->default_val(args.include_mode); app.add_option("--verbosity", args.log_level, "Verbosity level of logging") ->transform(CLI::CheckedTransformer(log_level_names, CLI::ignore_case)) @@ -149,6 +201,7 @@ Args cli::parse(int argc, const char *const *argv) { if (!root_node_values.empty()) { args.root_node = make_octree_id(root_node_values); } + args.continuation_mode = make_continuation_mode(resume, overwrite); } catch (const CLI::ParseError &e) { std::exit(app.exit(e)); } diff --git a/src/dag_builder/cli.h b/src/dag_builder/cli.h index 5d78166b..3d88509e 100644 --- a/src/dag_builder/cli.h +++ b/src/dag_builder/cli.h @@ -7,6 +7,7 @@ #include +#include "build.h" #include "octree/Id.h" #include "Range.h" #include "uv/unwrap.h" @@ -19,14 +20,17 @@ struct Args { std::filesystem::path input_path; std::filesystem::path output_path; octree::Id root_node; - Range level_range; + AnyRange level_range; uv::Algorithm uv_unwrap_algorithm; uint32_t clusters_per_partition; - float target_ratio; + std::optional target_ratio; + std::optional target_error; bool write_debug_meshes; - bool overwrite; + bool parallelize; + dag::IncludeMode include_mode; + ContinuationMode continuation_mode; }; Args parse(int argc, const char *const *argv); diff --git a/src/dag_builder/cluster.h b/src/dag_builder/cluster.h index a673ed86..9ddca2da 100644 --- a/src/dag_builder/cluster.h +++ b/src/dag_builder/cluster.h @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -10,7 +11,7 @@ #include "TextureSet.h" struct Cluster { - uint32_t id = -1; + uint32_t id = std::numeric_limits::max(); std::vector vertex_indices; // indices into Clustering::positions std::vector local_triangles; // indices into this->vertex_indices std::vector uvs; // per local vertex diff --git a/src/dag_builder/clusterize.h b/src/dag_builder/clusterize.h index 43e5d377..df2e4a4c 100644 --- a/src/dag_builder/clusterize.h +++ b/src/dag_builder/clusterize.h @@ -27,7 +27,7 @@ struct ClusterOptions { uint32_t min_triangles = 256; uint32_t max_triangles = 256; float cone_weight = 0.5; - float split_factor = 0.0; + float split_factor = 2.0; }; inline std::vector clusterize( diff --git a/src/dag_builder/compact.h b/src/dag_builder/compact.h index b1b8343c..bda6630c 100644 --- a/src/dag_builder/compact.h +++ b/src/dag_builder/compact.h @@ -2,11 +2,14 @@ #include #include +#include #include #include "cluster.h" +#include "enumerate.h" #include "slice.h" +#include "validate.h" inline void compact_cluster_inplace(Cluster &cluster) { const size_t vertex_count = cluster.vertex_count(); @@ -23,7 +26,7 @@ inline void compact_cluster_inplace(Cluster &cluster) { } // Build remap - constexpr uint32_t invalid_remap = -1; + constexpr uint32_t invalid_remap = std::numeric_limits::max(); std::vector remap(vertex_count, invalid_remap); std::vector new_vertex_indices; std::vector new_uvs; @@ -57,7 +60,7 @@ inline void compact_cluster_inplace(Cluster &cluster) { inline void remove_unused_vertices_inplace(Clustering &clustering) { const uint32_t vertex_count = clustering.vertex_count(); - constexpr uint32_t invalid_vertex = -1; + constexpr uint32_t invalid_vertex = std::numeric_limits::max(); std::vector vertex_remap(vertex_count, invalid_vertex); uint32_t next_index = 0; for (Cluster &cluster : clustering.clusters) { diff --git a/src/dag_builder/dag_node.h b/src/dag_builder/dag_node.h index 13970d72..249f957b 100644 --- a/src/dag_builder/dag_node.h +++ b/src/dag_builder/dag_node.h @@ -21,6 +21,10 @@ struct ClusterBatch { static ClusterBatch make_leaves(Clustering clustering) { return {clustering, {}}; } + + bool is_leaves() const { + return this->child_map.empty(); + } }; /* diff --git a/src/dag_builder/error_bounds.h b/src/dag_builder/error_bounds.h deleted file mode 100644 index 127d3f23..00000000 --- a/src/dag_builder/error_bounds.h +++ /dev/null @@ -1,26 +0,0 @@ -#pragma once - -#include -#include - -constexpr double EARTH_CIRCUMFERENCE_M = 40075016.686; // Earth's circumference in metres -constexpr uint32_t MAX_LEVEL = 30; - -// Compute metres per pixel at equator for a given zoom level -constexpr double metres_per_pixel(uint32_t zoom) { - return EARTH_CIRCUMFERENCE_M / (256.0 * (1ULL << zoom)); -} - -constexpr std::array generate_meters_per_pixel_array() { - std::array arr = {}; - for (uint32_t zoom = 0; zoom < MAX_LEVEL; zoom++) { - arr[zoom] = metres_per_pixel(zoom); - } - return arr; -} - -constexpr std::array METERS_PER_PIXEL_AT_EQUATOR = generate_meters_per_pixel_array(); - -static_assert(std::abs(METERS_PER_PIXEL_AT_EQUATOR[0] - 156543) < 1, "level 0 mismatch"); -static_assert(std::abs(METERS_PER_PIXEL_AT_EQUATOR[10] - 152.874) < 0.001, "level 10 mismatch"); -static_assert(std::abs(METERS_PER_PIXEL_AT_EQUATOR[20] - 0.149) < 0.001, "level 20 mismatch"); diff --git a/src/dag_builder/main.cpp b/src/dag_builder/main.cpp index ed35633d..355a7923 100644 --- a/src/dag_builder/main.cpp +++ b/src/dag_builder/main.cpp @@ -9,6 +9,7 @@ #include "octree/storage/MeshStorage.h" #include "octree/storage/open.h" #include "storage.h" +#include "ContinuationMode.h" int main(int argc, char **argv) { const cli::Args args = cli::parse(argc, argv); @@ -18,15 +19,23 @@ int main(int argc, char **argv) { try { const octree::IndexedMeshStorage input_storage = octree::open_folder_indexed(args.input_path); octree::IndexedDagStorage output_storage = octree::open_folder_indexed(args.output_path); - output_storage.settings().allow_overwrite = args.overwrite; + output_storage.settings().allow_overwrite = args.continuation_mode == ContinuationMode::Overwrite; - const dag::BuildOptions options{ + dag::BuildOptions options{ .clusters_per_partition = args.clusters_per_partition, .target_ratio = args.target_ratio, + .relative_target_error = args.target_error, .uv_unwrap_algorithm = args.uv_unwrap_algorithm, .root_node = args.root_node, + .include_mode = args.include_mode, .write_debug_meshes = args.write_debug_meshes, + .parallelize = args.parallelize, + .continuation_mode = args.continuation_mode }; + // If the user specified neither target, fall back to a default error. + if (!args.target_ratio && !args.target_error) { + options.relative_target_error = 0.01f; + } dag::build_levels(input_storage, output_storage, options, args.level_range); output_storage.save_index(); diff --git a/src/dag_builder/merge/clusterings.cpp b/src/dag_builder/merge/clusterings.cpp index bf11223c..97ff5c91 100644 --- a/src/dag_builder/merge/clusterings.cpp +++ b/src/dag_builder/merge/clusterings.cpp @@ -1,118 +1,642 @@ +#include +#include +#include +#include #include +#include +#include +#include #include +#include +#include +#include +#include #include #include +#include +#include "OffsetTable.h" +#include "SegmentedBuffer.h" +#include "VertexInClustering.h" +#include "UnionFind.h" +#include "build_config.h" #include "cluster.h" -#include "validate.h" -#include "enumerate.h" #include "compact.h" +#include "enumerate.h" #include "merge/clusterings.h" +#include "mesh/boundary.h" #include "mesh/triangle_compare.h" +#include "range_utils.h" #include "spatial_lookup/Hashmap.h" +#include "validate.h" namespace detail { -/* -std::array hash_image(const cv::Mat &input_image) { - std::array hash_buffer; - cv::Mat hash_wrapper(1, 16, CV_8UC1, hash_buffer.data()); - cv::img_hash::averageHash(input_image, hash_wrapper); - return hash_buffer; + +uint32_t to_flat_index(const OffsetTable &base_offsets, const VertexInClustering &vertex) { + return base_offsets.offset(vertex.clustering_index, vertex.global_vertex_index); } -*/ + +std::vector find_boundary_vertices(const std::span clusterings, const OffsetTable &offsets) { + const size_t total_vertex_count = offsets.total_size(); + std::vector boundary_vertices; + boundary_vertices.reserve(total_vertex_count / 2); + + std::vector boundary_vertex_mask; + for (const auto &[clustering_index, clustering] : enumerate(clusterings)) { + for (const Cluster &cluster : clustering.clusters) { + mesh::build_boundary_vertex_mask(cluster.local_triangles, cluster.vertex_count(), boundary_vertex_mask); + + for (const auto &[local_vertex_index, on_boundary] : enumerate(boundary_vertex_mask)) { + if (!on_boundary) { + continue; + } + + const uint32_t global_vertex_index = cluster.vertex_indices[local_vertex_index]; + boundary_vertices.push_back(VertexInClustering{ + .clustering_index = static_cast(clustering_index), + .global_vertex_index = global_vertex_index + }); + } + } + } + + return boundary_vertices; } -Clustering merge_clusterings(const std::span clusterings, const double epsilon) { - if (clusterings.empty()) { +struct StrongDouble { + double value; + + std::strong_ordering operator<=>(const StrongDouble &other) const { + return std::strong_order(value, other.value); + } + + bool operator==(const StrongDouble &other) const { + return std::strong_order(value, other.value) == 0; + } +}; + +struct CandidateEdge { + uint32_t start_index; + uint32_t end_index; + double distance_sq; + + auto operator<=>(const CandidateEdge &other) const { + return std::tuple(StrongDouble{this->distance_sq}, start_index, end_index) <=> + std::tuple(StrongDouble{other.distance_sq}, other.start_index, other.end_index); + } + auto operator==(const CandidateEdge &other) const { + return std::tuple(StrongDouble{this->distance_sq}, start_index, end_index) == + std::tuple(StrongDouble{other.distance_sq}, other.start_index, other.end_index); + } +}; + +void validate_epsilon(const double epsilon) { + if (!(epsilon > 0.0) || !std::isfinite(epsilon)) { + throw std::invalid_argument("merge_clusterings: epsilon must be finite and positive"); + } +} + +OffsetTable build_offset_table(const std::span clusterings) { + OffsetTable offsets(clusterings.size()); + + for (const Clustering &clustering : clusterings) { + offsets.append_length(clustering.vertex_count()); + } + + return offsets; +} + +spatial_lookup::Hashmap3d build_spatial_map( + const std::span clusterings, + const std::span boundary_vertices, + const double epsilon) { + spatial_lookup::Hashmap3d spatial_map(epsilon * 5); + + for (const VertexInClustering &vertex : boundary_vertices) { + const glm::dvec3 &position = clusterings[vertex.clustering_index].positions[vertex.global_vertex_index]; + spatial_map.insert(position, vertex); + } + + return spatial_map; +} + +template +std::vector build_canonical_indices(UnionFind_ &union_find) { + const uint32_t total_vertex_count = union_find.size(); + constexpr const uint32_t INVALID_INDEX = std::numeric_limits::max(); + std::vector vertex_to_canonical(total_vertex_count, INVALID_INDEX); + uint32_t next_index = 0; + + for (uint32_t flat_index = 0; flat_index < total_vertex_count; flat_index++) { + const uint32_t repr_index = union_find.find(flat_index); + uint32_t &repr_canonical_vertex = vertex_to_canonical[repr_index]; + if (repr_canonical_vertex == INVALID_INDEX) { + repr_canonical_vertex = next_index; + next_index++; + } + vertex_to_canonical[flat_index] = repr_canonical_vertex; + } + + return vertex_to_canonical; +} + +std::vector build_average_positions( + const std::span clusterings, + const OffsetTable &base_offsets, + const std::vector &vertex_to_canonical +) { + if (vertex_to_canonical.empty()) { return {}; } - if (clusterings.size() == 1) { - return clusterings[0]; + + const uint32_t position_count = std::ranges::max(vertex_to_canonical) + 1; + std::vector new_positions(position_count); + std::vector counts(position_count, 0); + + for (const auto &[clustering_index, clustering] : enumerate(clusterings)) { + const uint32_t offset = base_offsets.get_begin(clustering_index); + + for (const auto &[global_vertex_index, position] : enumerate(clustering.positions)) { + const uint32_t flat_index = offset + global_vertex_index; + const uint32_t canonical_index = vertex_to_canonical[flat_index]; + + new_positions[canonical_index] += position; + counts[canonical_index]++; + } } - // Compute per-clustering base offsets for vertex_to_canonical lookup - const uint32_t total_vertex_count = sum(clusterings, [&](const auto &c) { return c.vertex_count(); }); - std::vector base_offsets(clusterings.size()); - { - uint32_t offset = 0; - for (size_t ci = 0; ci < clusterings.size(); ci++) { - base_offsets[ci] = offset; - offset += clusterings[ci].vertex_count(); - } - } - - // Build a position -> canonical index mapping using epsilon-ball deduplication - std::vector vertex_to_canonical(total_vertex_count); - spatial_lookup::Hashmap3d vertex_lookup(epsilon); - std::vector new_positions; - new_positions.reserve(total_vertex_count); - - std::vector matches; - for (size_t ci = 0; ci < clusterings.size(); ci++) { - for (uint32_t vi = 0; vi < clusterings[ci].vertex_count(); vi++) { - const glm::dvec3 &position = clusterings[ci].positions[vi]; - const uint32_t global_index = base_offsets[ci] + vi; - matches.clear(); - if (vertex_lookup.find_all_near(position, epsilon, matches)) { - vertex_to_canonical[global_index] = matches[0]; - } else { - const uint32_t next_index = static_cast(new_positions.size()); - vertex_lookup.insert(position, next_index); - new_positions.push_back(position); - vertex_to_canonical[global_index] = next_index; - } + for (const auto &[position_index, count] : enumerate(counts)) { + DEBUG_ASSERT(count > 0); + new_positions[position_index] /= static_cast(count); + } + + return new_positions; +} + +std::vector build_last_writer_positions( + const std::span clusterings, + const OffsetTable &base_offsets, + const std::vector &vertex_to_canonical +) { + if (vertex_to_canonical.empty()) { + return {}; + } + + const uint32_t position_count = std::ranges::max(vertex_to_canonical) + 1; + std::vector new_positions(position_count); + + for (const auto &[clustering_index, clustering] : enumerate(clusterings)) { + const uint32_t offset = base_offsets.get_begin(clustering_index); + + for (const auto &[global_vertex_index, position] : enumerate(clustering.positions)) { + const uint32_t flat_index = offset + global_vertex_index; + const uint32_t canonical_index = vertex_to_canonical[flat_index]; + new_positions[canonical_index] = position; } } - const size_t unique_count = new_positions.size(); - const size_t shared_count = total_vertex_count - unique_count; - LOG_DEBUG("Merging with {} shared and {} unique vertices", shared_count, unique_count); + return new_positions; +} +void filter_degenerate_triangle(Clustering& clustering) { + for (Cluster &cluster : clustering.clusters) { + const size_t removed = std::erase_if(cluster.local_triangles, [&](const glm::uvec3 &local_triangle) { + const glm::uvec3 global_triangle( + cluster.vertex_indices[local_triangle.x], + cluster.vertex_indices[local_triangle.y], + cluster.vertex_indices[local_triangle.z]); + + return mesh::is_degenerate(global_triangle); + }); + + // Compact vertices if we found any. + if (removed > 0) { + compact_cluster_inplace(cluster); + } + } +} + +Clustering rebuild_clustering( + const std::span clusterings, + const OffsetTable &base_offsets, + const std::vector &vertex_to_canonical, + const bool average_positions, + const bool remove_degenerate_triangles) { Clustering merged; - merged.positions = std::move(new_positions); - for (size_t ci = 0; ci < clusterings.size(); ci++) { - const Clustering &clustering = clusterings[ci]; - // Merge textures - std::vector new_texture_ids; - new_texture_ids.reserve(clustering.textures.size()); + // Build new position buffer + if (average_positions) { + merged.positions = build_average_positions(clusterings, base_offsets, vertex_to_canonical); + } else { + merged.positions = build_last_writer_positions(clusterings, base_offsets, vertex_to_canonical); + } + + std::vector texture_id_map; + for (const auto &[clustering_index, clustering] : enumerate(clusterings)) { + // Prepare mapping from old to new texture id + texture_id_map.clear(); + texture_id_map.reserve(clustering.textures.size()); for (const cv::Mat &texture : clustering.textures) { - new_texture_ids.push_back(merged.textures.add(texture)); + texture_id_map.push_back(merged.textures.add(texture)); } - // Merge clusters with adjusted indices + const uint32_t offset = base_offsets.get_begin(clustering_index); + for (const Cluster &cluster : clustering.clusters) { Cluster new_cluster; new_cluster.local_triangles = cluster.local_triangles; new_cluster.uvs = cluster.uvs; new_cluster.id = cluster.id; - new_cluster.texture_id = new_texture_ids[cluster.texture_id]; + new_cluster.texture_id = texture_id_map[cluster.texture_id]; + new_cluster.absolute_error = cluster.absolute_error; new_cluster.vertex_indices.reserve(cluster.vertex_count()); for (const uint32_t vertex_index : cluster.vertex_indices) { - new_cluster.vertex_indices.push_back(vertex_to_canonical[base_offsets[ci] + vertex_index]); + const uint32_t global_index = offset + vertex_index; + new_cluster.vertex_indices.push_back(vertex_to_canonical[global_index]); } merged.clusters.push_back(std::move(new_cluster)); } } - // Remove degenerate triangles - for (Cluster &cluster : merged.clusters) { - const size_t removed = std::erase_if(cluster.local_triangles, [&](const auto &local_triangle) { - const glm::uvec3 global_triangle( - cluster.vertex_indices[local_triangle.x], - cluster.vertex_indices[local_triangle.y], - cluster.vertex_indices[local_triangle.z]); - return mesh::is_degenerate(global_triangle); + if (remove_degenerate_triangles) { + filter_degenerate_triangle(merged); + } + + DEBUG_ASSERT(merged.cluster_count() == sum(clusterings, [](const Clustering &clustering) { return clustering.cluster_count(); })); + validate(merged); + return merged; +} + +struct BestMatch { + uint32_t flat_index; + double distance_sq; +}; + +UnionFindWithSizes build_epsilon_neighbourhoods( + const std::span clusterings, + const std::span boundary_vertices, + const OffsetTable &base_offsets, + const spatial_lookup::Hashmap3d &spatial_map, + const uint32_t total_vertex_count, + const double epsilon) { + UnionFindWithSizes union_find(total_vertex_count); + std::vector matches; + + for (const VertexInClustering &vertex : boundary_vertices) { + const glm::dvec3 &position = clusterings[vertex.clustering_index].positions[vertex.global_vertex_index]; + const uint32_t flat_index = to_flat_index(base_offsets, vertex); + + // Final all points within epsilon. + matches.clear(); + spatial_map.find_all_near(position, epsilon, matches); + + for (const VertexInClustering &match : matches) { + // Skip if from same source + if (match.clustering_index == vertex.clustering_index) { + continue; + } + + // Skip if already covered the reverse + const uint32_t match_flat_index = to_flat_index(base_offsets, match); + if (match_flat_index <= flat_index) { + continue; + } + + union_find.make_union(flat_index, match_flat_index); + } + } + + return union_find; +} +} // namespace detail + +// Greedily matches vertices within epsilon ball as long as they belong to a different clustering. +std::vector build_greedy_local_mapping( + const std::span clusterings, + const std::span boundary_vertices, + const OffsetTable &base_offsets, + const spatial_lookup::Hashmap3d &spatial_map, + const uint32_t total_vertex_count, + const double epsilon) { + + constexpr uint32_t INVALID_INDEX = std::numeric_limits::max(); + std::vector vertex_to_canonical(total_vertex_count, INVALID_INDEX); + uint32_t next_index = 0; + + // Preallocate for next loop + std::vector matches; + std::vector> best_by_cluster; + + for (const VertexInClustering &vertex : boundary_vertices) { + const glm::dvec3 &position = clusterings[vertex.clustering_index].positions[vertex.global_vertex_index]; + const uint32_t flat_index = detail::to_flat_index(base_offsets, vertex); + + uint32_t &canonical_index = vertex_to_canonical[flat_index]; + if (canonical_index != INVALID_INDEX) { + // Already mapped + continue; + } + + canonical_index = next_index; + next_index++; + + // Find all points within epsilon + matches.clear(); + spatial_map.find_all_near(position, epsilon, matches); + DEBUG_ASSERT(!matches.empty()); + + // Remove any matches that are in the same cluster as the current vertex + std::erase_if(matches, [&](const VertexInClustering &match) { + return match.clustering_index == vertex.clustering_index; }); - if (removed > 0) { - compact_cluster_inplace(cluster); + switch (matches.size()) { + case 0: + // Found no other vertex nearby -> unique + break; + + case 1: { + // Found a single other vertex nearby -> assign same canonical index + const VertexInClustering &other_vertex = matches[0]; + const uint32_t other_flat_index = detail::to_flat_index(base_offsets, other_vertex); + if (vertex_to_canonical[other_flat_index] == INVALID_INDEX) { + vertex_to_canonical[other_flat_index] = canonical_index; + } + break; + } + + default: { + // Found multiple nearby vertices -> merge with nearest per clustering + best_by_cluster.assign(clusterings.size(), std::nullopt); + + for (const VertexInClustering &match : matches) { + const uint32_t match_flat_index = detail::to_flat_index(base_offsets, match); + if (vertex_to_canonical[match_flat_index] != INVALID_INDEX) { + continue; + } + + const glm::dvec3 &match_position = clusterings[match.clustering_index].positions[match.global_vertex_index]; + const double distance_sq = glm::length2(match_position - position); + auto &best_opt = best_by_cluster[match.clustering_index]; + if (!best_opt.has_value() || distance_sq < best_opt->distance_sq) { + best_opt = detail::BestMatch{ + .flat_index = match_flat_index, + .distance_sq = distance_sq}; + } + } + + for (const auto &[cluster_index, best_opt] : enumerate(best_by_cluster)) { + if (best_opt.has_value()) { + vertex_to_canonical[best_opt->flat_index] = canonical_index; + } + } + + break; + } } } - validate(merged); - return merged; + // Assign new canonical index for non boundary vertices + for (uint32_t &canonical_index : vertex_to_canonical) { + if (canonical_index == INVALID_INDEX) { + canonical_index = next_index++; + } + } + + return vertex_to_canonical; +} + +// Merges vertices from different clusterings into the same output vertex if they are directly or transitively connected by epsilon-distance matches (thus can merge multiple from same clustering if there is an intermediate within epsilon distance of both). +std::vector build_connected_component_mapping( + const std::span clusterings, + const std::span boundary_vertices, + const OffsetTable &base_offsets, + const spatial_lookup::Hashmap3d &spatial_map, + const uint32_t total_vertex_count, + const double epsilon) { + UnionFindWithSizes neighbourhoods = detail::build_epsilon_neighbourhoods(clusterings, boundary_vertices, base_offsets, spatial_map, total_vertex_count, epsilon); + return detail::build_canonical_indices(neighbourhoods); +} + +// +std::vector build_multipartite_nearest_mapping( + const std::span clusterings, + const OffsetTable &base_offsets, + const std::span boundary_vertices, + const spatial_lookup::Hashmap3d &spatial_map, + const uint32_t total_vertex_count, + const double epsilon) { + UnionFindWithSizes neighbourhoods = detail::build_epsilon_neighbourhoods(clusterings, boundary_vertices, base_offsets, spatial_map, total_vertex_count, epsilon); + const auto sets = neighbourhoods.get_sets_sparse(); + + // Prepare output vector + constexpr const uint32_t INVALID_INDEX = std::numeric_limits::max(); + std::vector vertex_to_canonical(total_vertex_count, INVALID_INDEX); + uint32_t next_canonical_index = 0; + + // Preallocate structures used in the loop + // Keep in sync with the fallback check in merge_clusterings. + ASSERT(clusterings.size() <= 64); + std::vector> local_source_masks; + std::vector local_vertices; + UnionFind local_merges; + + std::vector edges; + for (const auto &[repr_index, set] : enumerate(sets.segments())) { + switch (set.size()) { + case 0: + // (0) Empty -> this vertex is part of another set and not its repr + continue; + case 1: + // (1) This vertex is unique -> add directly + vertex_to_canonical[set[0]] = next_canonical_index; + next_canonical_index++; + continue; + case 2: + // (2) Two vertices in the set -> already optimal and guaranteed from different sources + vertex_to_canonical[set[0]] = next_canonical_index; + vertex_to_canonical[set[1]] = next_canonical_index; + next_canonical_index++; + continue; + } + + // Precompute indices for local vertices + local_vertices.clear(); + local_vertices.reserve(set.size()); + for (const uint32_t flat_index : set) { + const uint32_t clustering_index = base_offsets.locate(flat_index).segment; + const uint32_t global_index = base_offsets.local_index(clustering_index, flat_index); + const VertexInClustering vertex{clustering_index, global_index}; + local_vertices.push_back(vertex); + } + + // Find candidate edges + edges.clear(); + if (set.size() > 16) { + LOG_WARN_BACKOFF("Merging clusterings with large epsilon ({}) encountered {} vertices within transitive epsilon neighbourhood.", epsilon, set.size()); + } + + for (uint32_t i = 0; i < set.size(); i++) { + const VertexInClustering vertex_a = local_vertices[i]; + const glm::dvec3 position_a = clusterings[vertex_a.clustering_index].positions[vertex_a.global_vertex_index]; + + for (uint32_t j = i + 1; j < set.size(); j++) { + const VertexInClustering vertex_b = local_vertices[j]; + if (vertex_a.clustering_index == vertex_b.clustering_index) { + // Skip if from same source + continue; + } + const glm::dvec3 position_b = clusterings[vertex_b.clustering_index].positions[vertex_b.global_vertex_index]; + + edges.push_back(detail::CandidateEdge{ + .start_index = i, + .end_index = j, + .distance_sq = glm::distance2(position_a, position_b) + }); + } + } + + // Sort the edges + std::ranges::sort(edges); + + // Find compatible local merges + local_merges.reset(set.size()); + local_source_masks.assign(set.size(), {}); + for (const auto &[vertex_id, vertex] : enumerate(local_vertices)) { + local_source_masks[vertex_id].set(vertex.clustering_index); + } + + for (const detail::CandidateEdge &edge : edges) { + const uint32_t start_repr = local_merges.find(edge.start_index); + const uint32_t end_repr = local_merges.find(edge.end_index); + + if (start_repr == end_repr) { + continue; + } + + const std::bitset<64> &start_sources = local_source_masks[start_repr]; + const std::bitset<64> &end_sources = local_source_masks[end_repr]; + if ((start_sources & end_sources).any()) { + continue; + } + + if constexpr (IS_DEBUG_BUILD) { + if (edge.distance_sq > epsilon * epsilon) { + LOG_WARN_BACKOFF("Transitively welding vertices {} apart with epsilon {}", std::sqrt(edge.distance_sq), epsilon); + } + } + + const uint32_t merged_repr = local_merges.make_union(start_repr, end_repr); + const std::bitset<64> merged_sources = local_source_masks[start_repr] | local_source_masks[end_repr]; + local_source_masks[merged_repr] = merged_sources; + } + + // Commit local merges to global merge map + const auto local_sets = local_merges.get_sets_compact(); + for (const auto &local_set : local_sets.segments()) { + if (local_set.size() < 2) { + continue; + } + + const uint32_t canonical_index = next_canonical_index; + next_canonical_index++; + + for (const uint32_t vertex_id : local_set) { + const VertexInClustering &vertex = local_vertices[vertex_id]; + const uint32_t flat_index = detail::to_flat_index(base_offsets, vertex); + vertex_to_canonical[flat_index] = canonical_index; + } + } + } + + // Assign new canonical index for non-boundary vertices + for (uint32_t &canonical_index : vertex_to_canonical) { + if (canonical_index == INVALID_INDEX) { + canonical_index = next_canonical_index; + next_canonical_index++; + } + } + + return vertex_to_canonical; +} + +Clustering merge_clusterings( + const std::span clusterings, + const double epsilon, + MergeMode merge_mode, + const bool average_positions) { + + if (clusterings.empty()) { + return {}; + } + + detail::validate_epsilon(epsilon); + + if (clusterings.size() == 1) { + return clusterings[0]; + } + + const uint32_t total_vertex_count = sum(clusterings, [](const auto &c) { return c.vertex_count(); }); + const OffsetTable base_offsets = detail::build_offset_table(clusterings); + const std::vector boundary_vertices = detail::find_boundary_vertices(clusterings, base_offsets); + const spatial_lookup::Hashmap3d spatial_map = detail::build_spatial_map(clusterings, boundary_vertices, epsilon); + + if (merge_mode == MergeMode::MultipartiteNearest && clusterings.size() > 64) { + LOG_WARN("Cannot merge {} clusterings with MergeMode::MultipartiteNearest (max 64 sources), falling back to MergeMode::GreedyLocal", clusterings.size()); + merge_mode = MergeMode::GreedyLocal; + } + + std::vector vertex_to_canonical; + + switch (merge_mode) { + case MergeMode::GreedyLocal: + vertex_to_canonical = build_greedy_local_mapping( + clusterings, + boundary_vertices, + base_offsets, + spatial_map, + total_vertex_count, + epsilon); + break; + + case MergeMode::ConnectedComponents: + vertex_to_canonical = build_connected_component_mapping( + clusterings, + boundary_vertices, + base_offsets, + spatial_map, + total_vertex_count, + epsilon); + break; + + case MergeMode::MultipartiteNearest: + vertex_to_canonical = build_multipartite_nearest_mapping( + clusterings, + base_offsets, + boundary_vertices, + spatial_map, + total_vertex_count, + epsilon); + break; + } + + const size_t unique_count = max(vertex_to_canonical) + 1; + + const size_t shared_count = static_cast(total_vertex_count) - unique_count; + + if (shared_count == 0) { + LOG_WARN("Merging clusterings with epsilon {} welded no vertices", epsilon); + } + + LOG_DEBUG("Merging with {} shared and {} unique vertices", shared_count, unique_count); + + return detail::rebuild_clustering( + clusterings, + base_offsets, + vertex_to_canonical, + average_positions, + merge_mode == MergeMode::ConnectedComponents); } diff --git a/src/dag_builder/merge/clusterings.h b/src/dag_builder/merge/clusterings.h index 64c44fb9..bee38f9a 100644 --- a/src/dag_builder/merge/clusterings.h +++ b/src/dag_builder/merge/clusterings.h @@ -4,4 +4,16 @@ #include "cluster.h" -Clustering merge_clusterings(const std::span clusterings, const double epsilon); +enum class MergeMode { + GreedyLocal, + ConnectedComponents, + MultipartiteNearest, +}; + +// Merges clusterings into a single vertex space, welding vertices of different clusterings +// that lie within epsilon of each other. Clusters are neither added, removed nor reordered. +Clustering merge_clusterings( + const std::span clusterings, + const double epsilon, + const MergeMode merge_mode = MergeMode::MultipartiteNearest, + const bool average_positions = true); diff --git a/src/dag_builder/merge/clusters.h b/src/dag_builder/merge/clusters.h index 07e8ce5b..13a832c5 100644 --- a/src/dag_builder/merge/clusters.h +++ b/src/dag_builder/merge/clusters.h @@ -25,12 +25,13 @@ #include "mesh/split.h" #include "mesh/topology.h" #include "opencv_utils.h" -#include "partition.h" #include "range_utils.h" #include "uv/unwrap.h" #include "vector_utils.h" #include "atlas/TextureBaker.h" #include "TinyVector.h" +#include "Partitioning.h" + namespace detail { inline bool check_all_same_texture(const Clustering &clustering, const std::span cluster_indices) { @@ -477,12 +478,16 @@ struct UvMap { std::vector uvs; }; +// Algorithm used to retry a UV unwrap when the requested algorithm fails. +inline constexpr uv::Algorithm fallback_algorithm = uv::Algorithm::TutteBarycentricMapping; + inline UvMap unwrap_merged_cluster( const Clustering &clustering, const Cluster &merged_cluster, HybridIndexPairMap &merged_to_original, - const std::span cluster_indices) { + const std::span cluster_indices, + const uv::Algorithm algorithm = uv::DEFAULT_ALGORITHM) { // Materialize cluster mesh const mesh::Simple merged_mesh = materialize_cluster(merged_cluster, clustering.positions); DEBUG_ASSERT(!merged_mesh.has_uvs()); @@ -501,6 +506,10 @@ inline UvMap unwrap_merged_cluster( // Prepare atlas for new texture TextureBaker baker; + // Texture map id each component is baked into, used to gather the packed + // UVs back into merged-vertex order. + std::vector component_map_ids(components.size()); + // Preallocate std::vector source_clusters; source_clusters.reserve(cluster_indices.size()); @@ -537,15 +546,14 @@ inline UvMap unwrap_merged_cluster( }); const uint32_t texture_id = clustering.clusters[cluster_index].texture_id; const cv::Mat texture = clustering.textures[texture_id]; - baker.add_mesh(TexturedMesh{component.triangles, TextureMap{uvs, texture}}); + component_map_ids[component_index] = baker.add_mesh(TexturedMesh{component.triangles, TextureMap{uvs, texture}}); } else { // If multiple clusters are relevant, we have to perform an unwrap. - /*auto result = uv::unwrap(component, uv::Algorithm::AsRigidAsPossible); - if (!result.has_value()) { - LOG_WARN("Failed to unwrap using ARAP: {}", result.error().description()); - result = uv::unwrap(component); - }*/ - auto result = uv::unwrap(component); + auto result = uv::unwrap(component, algorithm); + if (!result.has_value() && algorithm != fallback_algorithm) { + LOG_WARN("Failed to unwrap using requested algorithm: {}, falling back to TutteBarycentricMapping", result.error().description()); + result = uv::unwrap(component, fallback_algorithm); + } if (!result.has_value()) { LOG_ERROR_AND_EXIT("Failed to unwrap using default: {}", result.error().description()); } @@ -595,15 +603,34 @@ inline UvMap unwrap_merged_cluster( .target = target_triangle, }); } - baker.add_composition(comp); + component_map_ids[component_index] = baker.add_composition(comp); + } + } + + // Combine all textures together, compressing resolution in proportion to + // how many clusters are being merged so the DAG's texture budget shrinks + // going up the LOD hierarchy. + const atlas::Packing packing = baker.pack(); + const double compression_ratio = 1.0 / double(cluster_indices.size()); + const glm::uvec2 texture_size = compute_bake_texture_size(packing, compression_ratio, 4096); + const auto baked = baker.bake(packing, texture_size); + + // Gather the per-component packed UVs into merged-vertex order. Each + // component is baked using component-local vertex indices, so the flat + // buffer cannot be assigned to the merged cluster directly. + std::vector merged_uvs(merged_cluster.vertex_count()); + for (const auto &[component_index, component] : enumerate(components)) { + const std::span component_uvs = baked.uvs_for(component_map_ids[component_index]); + const std::vector &local_to_merged = component_to_merged[component_index]; + DEBUG_ASSERT(component_uvs.size() == local_to_merged.size()); + for (const auto [local_index, merged_index] : enumerate(local_to_merged)) { + merged_uvs[merged_index] = component_uvs[local_index]; } } - // Combine all textures together - const auto baked = baker.bake(glm::uvec2(2048)); return UvMap{ baked.texture(), - baked.uvs() + std::move(merged_uvs) }; } @@ -615,7 +642,8 @@ struct ClusterAndTexture { inline ClusterAndTexture merge_clusters_with_unwrap( const Clustering &clustering, const std::span cluster_indices, - const std::span vertex_remap) { + const std::span vertex_remap, + const uv::Algorithm algorithm = uv::DEFAULT_ALGORITHM) { // Merge raw geometry auto [merged_cluster, merged_to_original] = merge_cluster_geometry(clustering, cluster_indices, vertex_remap); @@ -633,14 +661,14 @@ inline ClusterAndTexture merge_clusters_with_unwrap( } // Unwrap cluster and generate new texture - auto [texture, uvs] = unwrap_merged_cluster(clustering, merged_cluster, merged_to_original, cluster_indices); + auto [texture, uvs] = unwrap_merged_cluster(clustering, merged_cluster, merged_to_original, cluster_indices, algorithm); merged_cluster.uvs = std::move(uvs); return ClusterAndTexture{merged_cluster, texture}; } } -inline Clustering merge_clusters(const Clustering &clustering, const Partitioning &partitioning) { +inline Clustering merge_clusters(const Clustering &clustering, const Partitioning &partitioning, const uv::Algorithm algorithm = uv::DEFAULT_ALGORITHM) { const uint32_t cluster_count = clustering.cluster_count(); const size_t partition_count = partitioning.partition_count; const std::vector &cluster_partitions = partitioning.cluster_partitions; @@ -664,9 +692,8 @@ inline Clustering merge_clusters(const Clustering &clustering, const Partitionin cluster_indices.push_back(i); } } - if (cluster_indices.empty()) { - continue; - } + // Empty partitions would break the cluster index == partition index mapping relied on by build_lod. + ASSERT(!cluster_indices.empty()); // Check if we need to perform a fresh uv unwrap due to different textures or inconsistent uvs const bool needs_unwrap = detail::check_merge_needs_unwrap(clustering, cluster_indices); @@ -675,7 +702,7 @@ inline Clustering merge_clusters(const Clustering &clustering, const Partitionin Texture texture; if (needs_unwrap) { // We need to perform a fresh uv unwrap and generate a new texture - const auto result = detail::merge_clusters_with_unwrap(clustering, cluster_indices, vertex_remap); + const auto result = detail::merge_clusters_with_unwrap(clustering, cluster_indices, vertex_remap, algorithm); merged_cluster = result.cluster; texture = result.texture; unwrap_count++; @@ -688,6 +715,14 @@ inline Clustering merge_clusters(const Clustering &clustering, const Partitionin reuse_count++; } merged_cluster.texture_id = textures.add(texture); + + // Carry the largest child error into the merged cluster. + double max_child_error = 0.0; + for (const uint32_t cluster_index : cluster_indices) { + max_child_error = std::max(max_child_error, clustering.clusters[cluster_index].absolute_error); + } + merged_cluster.absolute_error = max_child_error; + partitioned_clusters.push_back(std::move(merged_cluster)); } diff --git a/src/dag_builder/meshopt.h b/src/dag_builder/meshopt.h index 305fa834..54be4556 100644 --- a/src/dag_builder/meshopt.h +++ b/src/dag_builder/meshopt.h @@ -8,6 +8,7 @@ #include #include "glm_utils.h" +#include "range_utils.h" namespace meshopt { @@ -212,7 +213,7 @@ inline PartitionClustersResult partition_clusters( const size_t cluster_count = cluster_vertex_counts.size(); std::vector partitions(cluster_count); - const size_t total_index_count = std::accumulate(cluster_vertex_counts.begin(), cluster_vertex_counts.end(), 0u); + const size_t total_index_count = sum(cluster_vertex_counts, size_t{0}); const std::span positions_flat = flatten(vertex_positions); const size_t partition_count = meshopt_partitionClusters( diff --git a/src/dag_builder/optimize.h b/src/dag_builder/optimize.h index 15e68937..af76c8f6 100644 --- a/src/dag_builder/optimize.h +++ b/src/dag_builder/optimize.h @@ -13,7 +13,7 @@ inline Cluster optimize(Cluster cluster) { } inline void optimize_inplace(Clustering &clustering) { - for (auto &cluster : clustering) { + for (auto &cluster : clustering.clusters) { optimize_inplace(cluster); } } diff --git a/src/dag_builder/partition.h b/src/dag_builder/partition.h index 13565475..27c9e396 100644 --- a/src/dag_builder/partition.h +++ b/src/dag_builder/partition.h @@ -31,16 +31,14 @@ #include "vector_utils.h" #include "atlas/Packer.h" #include "mesh/igl/manifold.h" +#include "merge/clusters.h" +#include "Partitioning.h" + struct PartitionOptions { uint32_t clusters_per_partition = 4; }; -struct [[nodiscard]] Partitioning { - uint32_t partition_count = 0; - std::vector cluster_partitions; // original cluster index -> partition index -}; - inline Partitioning create_partitioning(const Clustering &clustering, const PartitionOptions options = {}) { const uint32_t cluster_count = clustering.clusters.size(); const uint32_t clusters_per_partition = options.clusters_per_partition; @@ -82,14 +80,12 @@ inline Partitioning create_partitioning(const Clustering &clustering, const Part std::move(partition_result.cluster_partitions)}; } -Clustering merge_clusters(const Clustering &clustering, const Partitioning &partitioning); -inline Clustering apply_partitioning(const Clustering &clustering, const Partitioning &partitioning) { - return merge_clusters(clustering, partitioning); +inline Clustering apply_partitioning(const Clustering &clustering, const Partitioning &partitioning, const uv::Algorithm algorithm = uv::DEFAULT_ALGORITHM) { + return merge_clusters(clustering, partitioning, algorithm); } -#include "merge/clusters.h" -inline Clustering partition(const Clustering &clustering, const PartitionOptions &options = {}) { +inline Clustering partition(const Clustering &clustering, const PartitionOptions &options = {}, const uv::Algorithm algorithm = uv::DEFAULT_ALGORITHM) { const Partitioning partitioning = create_partitioning(clustering, options); - return apply_partitioning(clustering, partitioning); + return apply_partitioning(clustering, partitioning, algorithm); } diff --git a/src/dag_builder/simplify.h b/src/dag_builder/simplify.h index 7ad1b1a3..a8cc5dd0 100644 --- a/src/dag_builder/simplify.h +++ b/src/dag_builder/simplify.h @@ -1,10 +1,12 @@ #pragma once +#include #include #include #include #include #include +#include #include #include @@ -52,14 +54,47 @@ struct VertexLock { explicit VertexLock(T t) : v(std::move(t)) {} }; +// Target ratio used when neither a ratio nor an error target is specified +inline constexpr float DEFAULT_TARGET_RATIO = 0.5f; + +// Controls how the input cluster's error is combined with the error introduced by this +// simplification when setting the simplified cluster's error. +enum class ErrorMode { + Overwrite, // keep only this simplification's error, discarding the input error + Add, // add this simplification's error to the input error + Max, // keep the larger of the input error and this simplification's error +}; + struct SimplifyOptions { - float target_ratio = 0.5; - float absolute_target_error = meshopt::NO_TARGET_ERROR; + std::optional target_ratio; + std::optional absolute_target_error; VertexLock vertex_lock = VertexLock::none(); float uv_weight = 0.5; + ErrorMode error_mode = ErrorMode::Overwrite; + // Emit exactly one output cluster per input cluster + bool preserve_cluster_count = false; + + // Default options apply a target ratio when no stop condition is specified. + static SimplifyOptions defaults() { + SimplifyOptions options; + options.target_ratio = DEFAULT_TARGET_RATIO; + return options; + } }; namespace detail { + inline double combine_error(const ErrorMode mode, const double input_error, const double current_error) { + switch (mode) { + case ErrorMode::Overwrite: + return current_error; + case ErrorMode::Add: + return input_error + current_error; + case ErrorMode::Max: + return std::max(input_error, current_error); + } + return current_error; + } + inline std::vector resolve_vertex_lock(const VertexLock& v, const Clustering& clustering, const Cluster& cluster) { constexpr const uint8_t UNLOCKED = VertexLock::UNLOCKED; constexpr const uint8_t LOCKED = VertexLock::LOCKED; @@ -99,7 +134,7 @@ namespace detail { [[nodiscard]] inline Clustering simplify( const Clustering& original_clustering, - const SimplifyOptions options = {} + const SimplifyOptions options = SimplifyOptions::defaults() ) { Clustering simplified_clustering; simplified_clustering.textures = original_clustering.textures; @@ -109,6 +144,10 @@ inline Clustering simplify( std::vector cluster_positions_f; std::vector cluster_uvs_f; std::vector vertex_remap; + + const float target_ratio = options.target_ratio.value_or(0.0f); + const float absolute_target_error = options.absolute_target_error.value_or(meshopt::NO_TARGET_ERROR); + for (const Cluster &original_cluster : original_clustering.clusters) { const size_t original_vertex_count = original_cluster.vertex_indices.size(); @@ -120,9 +159,16 @@ inline Clustering simplify( cluster_positions_f.clear(); cluster_positions_f.reserve(original_vertex_count); to_approximate_normalized(cluster_positions, cluster_positions_f, &bounds); - const float max_extents = glm::compMax(bounds.size()); - const float relative_target_error = options.absolute_target_error == meshopt::NO_TARGET_ERROR ? - meshopt::NO_TARGET_ERROR : options.absolute_target_error / (max_extents * 2); + const float max_extents = glm::compMax(bounds.size()) / 2.0f; + if (max_extents == 0.0f) { + // Empty or degenerate cluster + if (options.preserve_cluster_count) { + simplified_clustering.clusters.push_back(original_cluster); + } + continue; + } + const float relative_target_error = absolute_target_error == meshopt::NO_TARGET_ERROR ? + meshopt::NO_TARGET_ERROR : absolute_target_error / max_extents; // Prepare vertex attributes (uv) cluster_uvs_f.clear(); @@ -141,7 +187,11 @@ inline Clustering simplify( // Perform simplification const size_t original_triangle_count = original_cluster.local_triangles.size(); - const size_t target_triangle_count = static_cast(options.target_ratio * original_triangle_count); + size_t target_triangle_count = static_cast(target_ratio * original_triangle_count); + if (options.preserve_cluster_count) { + // Keep at least one triangle so the cluster cannot disappear + target_triangle_count = std::max(target_triangle_count, size_t{1}); + } meshopt::SimplifyResult result = meshopt::simplify_with_attributes( original_cluster.local_triangles, cluster_positions_f, @@ -153,7 +203,10 @@ inline Clustering simplify( relative_target_error, meshopt_SimplifyErrorAbsolute); if (result.triangles.empty()) { - // Simplification removed all triangles, go to next cluster + // Simplification removed all triangles + if (options.preserve_cluster_count) { + simplified_clustering.clusters.push_back(original_cluster); + } continue; } @@ -212,21 +265,26 @@ inline Clustering simplify( } } - // Make error absolute - const float absolute_error = result.relative_error * (max_extents * 2); + // Make error absolute and combine with the input cluster's error + const double absolute_error = result.relative_error * max_extents; + const double combined_error = detail::combine_error(options.error_mode, original_cluster.absolute_error, absolute_error); // Create new cluster Cluster simplified_cluster{ + .id = original_cluster.id, .vertex_indices = std::move(vertex_indices), .local_triangles = std::move(local_triangles), .uvs = std::move(uvs), .texture_id = original_cluster.texture_id, - .absolute_error = absolute_error + .absolute_error = combined_error }; validate(simplified_cluster, simplified_clustering.positions); simplified_clustering.clusters.push_back(std::move(simplified_cluster)); } + if (options.preserve_cluster_count) { + DEBUG_ASSERT(original_clustering.cluster_count() == simplified_clustering.cluster_count()); + } validate(simplified_clustering); return simplified_clustering; } diff --git a/src/dag_builder/slice.h b/src/dag_builder/slice.h index dd55c4fa..50c8a0b5 100644 --- a/src/dag_builder/slice.h +++ b/src/dag_builder/slice.h @@ -3,19 +3,38 @@ #include #include "cluster.h" +#include "enumerate.h" +#include "mesh/VertexMap.h" -inline Clustering slice_clusters(const Clustering &clustering, const std::span cluster_indices) { +struct ClusteringAndMap { + Clustering clustering; + VertexMap remap; +}; + +inline bool selects_all_clusters_in_order(const Clustering &clustering, const std::span cluster_indices) { + if (cluster_indices.size() != clustering.cluster_count()) { + return false; + } + for (const auto &[index, cluster_index] : enumerate(cluster_indices)) { + if (cluster_index != index) { + return false; + } + } + return true; +} + +inline ClusteringAndMap slice_clusters_with_map(const Clustering &clustering, const std::span cluster_indices) { const uint32_t cluster_count = clustering.cluster_count(); const uint32_t vertex_count = clustering.vertex_count(); - if (cluster_indices.size() == clustering.cluster_count()) { - return clustering; + if (selects_all_clusters_in_order(clustering, cluster_indices)) { + return {clustering, VertexMap::identity(vertex_count)}; } Clustering new_clustering; const uint32_t approximate_vertex_count = std::min(vertex_count, static_cast(vertex_count * static_cast(cluster_indices.size()) / static_cast(cluster_count) * 1.5)); new_clustering.positions.reserve(approximate_vertex_count); - constexpr uint32_t invalid_remap = -1; + constexpr uint32_t invalid_remap = VertexMap::invalid_index; std::vector vertex_remap(clustering.vertex_count(), invalid_remap); std::vector texture_remap(clustering.textures.size(), invalid_remap); @@ -50,5 +69,9 @@ inline Clustering slice_clusters(const Clustering &clustering, const std::span cluster_indices) { + return slice_clusters_with_map(clustering, cluster_indices).clustering; } diff --git a/src/dag_builder/split.h b/src/dag_builder/split.h index 066a35be..08f5cd0b 100644 --- a/src/dag_builder/split.h +++ b/src/dag_builder/split.h @@ -2,7 +2,6 @@ #include #include -#include #include #include @@ -17,33 +16,6 @@ #include "utils.h" #include "validate.h" -namespace { -template -inline constexpr auto int_div_ceil(A a, B b) { - static_assert(std::is_integral_v); - static_assert(std::is_integral_v); - using T = std::conditional_t<(sizeof(A) >= sizeof(B)), A, B>; - ASSERT(b != 0); - - const T a_t = static_cast(a); - const T b_t = static_cast(b); - const T q = a_t / b_t; - const T r = a_t % b_t; - - // Exact division -> already the ceiling - if (r == 0) { - return q; - } - - // Same sign - if ((a > 0 && b > 0) || (a < 0 && b < 0)) { - return q + 1; - } - - // Opposite signs -> truncation toward zero already gave the ceiling - return q; -} -} template inline Clustering split_each_into_equal_parts(const Clustering &input, const S num_parts) { @@ -66,10 +38,10 @@ inline Clustering split_each_into_equal_parts(const Clustering &input, const S n const bool has_uvs = !cluster.uvs.empty(); // Ensure we can meaningfully split triangles and vertices - ASSERT(triangle_count >= num_parts && "Cluster has fewer triangles than requested parts"); - ASSERT(vertex_count >= num_parts && "Cluster has fewer vertices than requested parts"); + ASSERT(triangle_count >= num_parts, "Cluster has fewer triangles than requested parts"); + ASSERT(vertex_count >= num_parts, "Cluster has fewer vertices than requested parts"); - // Initialize eptr and eind based on the triangle data< + // Initialize eptr and eind based on the triangle data std::vector eptr(triangle_count + 1); // Offsets in eind where faces start/end std::vector eind(triangle_count * 3); // Indices into vertex array @@ -137,7 +109,7 @@ inline Clustering split_each_into_equal_parts(const Clustering &input, const S n const uint32_t partition_index = epart[i]; Cluster &partition = partitions[partition_index]; glm::uvec3 triangle = cluster.local_triangles[i]; - for (uint8_t k=0; k<3; k++) { + for (uint8_t k = 0; k < 3; k++) { const uint32_t original_index = triangle[k]; uint32_t &new_index = remap[partition_index][original_index]; if (new_index == invalid_remap) { diff --git a/src/dag_builder/vertex_lock.h b/src/dag_builder/vertex_lock.h index 6c52cf4e..e334efa4 100644 --- a/src/dag_builder/vertex_lock.h +++ b/src/dag_builder/vertex_lock.h @@ -4,24 +4,29 @@ #include #include -#include +#include +#include #include +#include #include "TinyVector.h" #include "cluster.h" -#include "mesh/bounds.h" #include "mesh/triangle_compare.h" #include "simplify.h" #include "utils.h" +#include "VertexInCluster.h" -namespace detail { -struct VertexInCluster { - uint32_t cluster_index; - uint32_t local_vertex_index; -}; + +// Signed distance from a point to an axis-aligned box: negative inside (distance to the +// nearest face), positive outside (euclidean distance to the surface), zero on a face. +inline double signed_distance_to_bounds(const radix::geometry::Aabb3d &bounds, const glm::dvec3 &point) { + const glm::dvec3 gaps = glm::max(bounds.min - point, point - bounds.max); + const double outside = glm::length(glm::max(gaps, glm::dvec3(0.0))); + const double inside = glm::min(glm::compMax(gaps), 0.0); + return outside + inside; } -inline std::vector find_vertices_to_lock(const Clustering &clustering) { +inline std::vector find_vertices_to_lock(const Clustering &clustering, const radix::geometry::Aabb3d &node_bounds) { // Lock every triangle where any vertex is either on the border and near the bounds or shared between clusters. const uint32_t cluster_count = clustering.cluster_count(); const uint32_t vertex_count = clustering.vertex_count(); @@ -41,28 +46,22 @@ inline std::vector find_vertices_to_lock(const Clustering &clustering) } } - // Calculate safe bounds excluding the outer boundary of the mesh - const radix::geometry::Aabb3d bounds = calculate_bounds(clustering.positions); - const glm::dvec3 center = bounds.centre(); - const glm::dvec3 extents = bounds.size() / 2.0; - const glm::dvec3 unlocked_extents = extents * 0.99; - const radix::geometry::Aabb3d unlocked_bounds(center - unlocked_extents, center + unlocked_extents); - - // Lock vertices outside safe bounds that are on the boundary + // Lock boundary vertices on or beyond a node face, as these are shared with + // neighbouring nodes. + const double lock_margin = glm::compMax(node_bounds.size()) * 0.05; for (const uint32_t vertex_index : boundary_vertices) { - // Lock only vertices outside the unlocked bounds const glm::dvec3 &position = clustering.positions[vertex_index]; - if (!unlocked_bounds.contains(position)) { + if (signed_distance_to_bounds(node_bounds, position) >= -lock_margin) { vertices_to_lock.insert(vertex_index); } } // Find all vertices shared between at least 2 clusters - std::vector> cluster_membership(vertex_count); + std::vector> cluster_membership(vertex_count); for (uint32_t cluster_index = 0; cluster_index < cluster_count; cluster_index++) { const Cluster &cluster = clustering.clusters[cluster_index]; for (const auto [local_vertex_index, vertex_index] : enumerate(cluster.vertex_indices)) { - TinyVector &membership = cluster_membership[vertex_index]; + TinyVector &membership = cluster_membership[vertex_index]; membership.emplace_back(cluster_index, local_vertex_index); } } diff --git a/src/dag_convert_debug/main.cpp b/src/dag_convert_debug/main.cpp index dc46268b..2904c209 100644 --- a/src/dag_convert_debug/main.cpp +++ b/src/dag_convert_debug/main.cpp @@ -57,7 +57,7 @@ void export_storage_folder(const cli::Args &args) { continue; } - const mesh::Simple mesh = clustering_to_mesh(load_result.value().clustering, true); + const mesh::Simple mesh = clustering_to_mesh(load_result.value().clustering); const auto save_result = output_storage.save(id, mesh); if (!save_result.has_value()) { diff --git a/src/sf_merger/cut.h b/src/sf_merger/cut.h index 322359e2..b33d0a37 100644 --- a/src/sf_merger/cut.h +++ b/src/sf_merger/cut.h @@ -34,7 +34,7 @@ inline void cut_leaf_node( DEBUG_ASSERT(ctx.input.index().is(octree::NodeStatus::Leaf, id)); const SimpleMesh mesh = DEBUG_ASSERT_VAL(ctx.input.load(id)).value(); - LOG_TRACE("Cutting mesh at {} using mask with {} vertices and {} triangles", + LOG_TRACE("Cutting mesh at {} using mask with {} vertices and {} triangles", id, mask.mesh.vertex_count(), mask.mesh.face_count()); const Cow clipped = clip_on_mask(mesh, mask, ctx.keep_inside); if (clipped.is_ref()) { @@ -43,7 +43,7 @@ inline void cut_leaf_node( } else { const SimpleMesh &clipped_mesh = clipped; if (!clipped_mesh.is_empty()) { - LOG_TRACE("Mesh was clipped from {} vertices and {} triangles to {} vertices and {} triangles", + LOG_TRACE("Mesh was clipped from {} vertices and {} triangles to {} vertices and {} triangles", mesh.vertex_count(), mesh.face_count(), clipped_mesh.vertex_count(), clipped_mesh.face_count()); DEBUG_ASSERT_VAL(ctx.output.save(id, clipped_mesh)); } else { @@ -132,4 +132,3 @@ inline void cut_dataset( cut_dataset(input_dataset, mask, output_dataset, keep_inside); } - diff --git a/src/sf_merger/main.cpp b/src/sf_merger/main.cpp index 5a0bad24..bd20a16b 100644 --- a/src/sf_merger/main.cpp +++ b/src/sf_merger/main.cpp @@ -19,7 +19,7 @@ std::optional load_mask_from_path(const std::filesystem::path& path) { auto result = mask::load_from_path(path, radius_range); if (result.has_value()) { const auto mask = result.value(); - LOG_DEBUG("Loaded mask successfully ({} vertices, {} triangles)", + LOG_DEBUG("Loaded mask successfully ({} vertices, {} triangles)", mask.mesh.vertex_count(), mask.mesh.face_count()); return mask; } else { @@ -42,7 +42,7 @@ void run(const cli::MergeArgs& args) { LOG_TRACE("Creating output dataset at {}", args.output_path); std::filesystem::create_directories(args.output_path); octree::Storage output_dataset = octree::open_folder(args.output_path, false, octree::OpenOptions{.preferred_extension_with_dot = ".glb"}); - + std::optional mask = flatten(map(args.mask_path, load_mask_from_path)); return merge_datasets(base_dataset, new_dataset, output_dataset, mask); @@ -81,4 +81,3 @@ int main(int argc, char **argv) { // std::filesystem::copy_options::recursive | // std::filesystem::copy_options::overwrite_existing); } - diff --git a/src/terrainlib/MoreExact.h b/src/terrainlib/MoreExact.h deleted file mode 100644 index 266bb916..00000000 --- a/src/terrainlib/MoreExact.h +++ /dev/null @@ -1,20 +0,0 @@ -#pragma once - -#include - -// Generic helper to determine the "more exact" type, preferring floating point types -template -struct MoreExactType { -private: - static constexpr bool t1_is_fp = std::is_floating_point_v; - static constexpr bool t2_is_fp = std::is_floating_point_v; - -public: - using type = std::conditional_t< - (t1_is_fp && !t2_is_fp), T1, - std::conditional_t<(t2_is_fp && !t1_is_fp), T2, - std::conditional_t<(sizeof(T1) >= sizeof(T2)), T1, T2>>>; -}; - -template -using MoreExact = typename MoreExactType::type; \ No newline at end of file diff --git a/src/terrainlib/OffsetTable.h b/src/terrainlib/OffsetTable.h index 47aeb116..e476e2f5 100644 --- a/src/terrainlib/OffsetTable.h +++ b/src/terrainlib/OffsetTable.h @@ -6,23 +6,37 @@ #include #include -template -class OffsetTable { +template +class OffsetTable_ { public: using index_type = Index; struct range_type { index_type begin; // inclusive index_type end; // exclusive + + bool contains(index_type index) { + return this->begin <= index && index < end; + } + + bool empty() { + return this->begin >= this->end; + } }; struct locate_result { - index_type element; + index_type segment; range_type range; }; - OffsetTable() : _offsets{0} {} - + OffsetTable_() { + _offsets.push_back(0); + } + OffsetTable_(size_t size) { + _offsets.reserve(size + 1); + _offsets.push_back(0); + } + void clear() noexcept { this->_offsets.clear(); this->_offsets.push_back(0); @@ -44,13 +58,13 @@ class OffsetTable { return this->_offsets.back(); } - // Append an element given its length. + // Append an segment given its length. void append_length(const index_type length) { const index_type next_offset = this->_offsets.back() + length; this->_offsets.push_back(next_offset); } - // Append multiple elements given their lengths. + // Append multiple segments given their lengths. void append_lengths(const std::span lengths) { const size_t current_size = this->_offsets.size(); this->_offsets.reserve(current_size + lengths.size()); @@ -60,33 +74,33 @@ class OffsetTable { } } - // Modifies the start of an element, shrinking or extending the previous one accordingly. - void set_begin(const index_type element, const index_type new_begin) { - if (element == 0) { - throw std::out_of_range("cannot modify begin of first element"); + // Modifies the start of an segment, shrinking or extending the previous one accordingly. + void set_begin(const index_type segment, const index_type new_begin) { + if (segment == 0) { + throw std::out_of_range("cannot modify begin of first segment"); } - if (element >= this->size()) { - throw std::out_of_range("element out of range"); + if (segment >= this->size()) { + throw std::out_of_range("segment out of range"); } - const index_type prev_boundary = this->_offsets[element - 1]; - const index_type next_boundary = this->_offsets[element + 1]; + const index_type prev_boundary = this->_offsets[segment - 1]; + const index_type next_boundary = this->_offsets[segment + 1]; if (new_begin < prev_boundary || new_begin > next_boundary) { throw std::invalid_argument("new_begin out of valid adjacent boundaries"); } - this->_offsets[element] = new_begin; + this->_offsets[segment] = new_begin; } - // Modifies the end of an element, shrinking or extending the next one accordingly (if present). - void set_end(const index_type element, const index_type new_end) { - if (element >= this->size()) { - throw std::out_of_range("element out of range"); + // Modifies the end of an segment, shrinking or extending the next one accordingly (if present). + void set_end(const index_type segment, const index_type new_end) { + if (segment >= this->size()) { + throw std::out_of_range("segment out of range"); } - const index_type prev_boundary = this->_offsets[element]; - if (element + 1 < this->size()) { - const index_type next_boundary = this->_offsets[element + 2]; + const index_type prev_boundary = this->_offsets[segment]; + if (segment + 1 < this->size()) { + const index_type next_boundary = this->_offsets[segment + 2]; if (new_end < prev_boundary || new_end > next_boundary) { throw std::invalid_argument("new_end out of valid adjacent boundaries"); } @@ -96,27 +110,74 @@ class OffsetTable { } } - this->_offsets[element + 1] = new_end; + this->_offsets[segment + 1] = new_end; + } + + void remove_empty_segments() { + const auto new_end = std::unique(this->_offsets.begin(), this->_offsets.end()); + this->_offsets.erase(new_end, this->_offsets.end()); + } + + bool in_same_segment(const index_type first, const index_type second) const { + const range_type range = this->locate(first).range; + return range.contains(second); + } + + bool in_segment(const index_type segment, const index_type index) const { + const range_type range = this->segment_range(segment); + return range.contains(index); + } + + index_type get_begin(const index_type segment) const { + if (segment >= this->size()) { + throw std::out_of_range("segment out of range"); + } + + return this->_offsets[segment]; + } + + index_type get_end(const index_type segment) const { + if (segment >= this->size()) { + throw std::out_of_range("segment out of range"); + } + + return this->_offsets[segment + 1]; + } + + index_type offset(const index_type segment, const index_type local_index) const { + const range_type r = this->segment_range(segment); + + if (local_index >= r.end - r.begin) { + throw std::out_of_range("local index out of range"); + } + + return r.begin + local_index; + } + + // Inverse of offset(segment, local_index): buffer index -> local index within segment. + index_type local_index(const index_type segment, const index_type buffer_index) const { + const range_type r = this->segment_range(segment); + return buffer_index - r.begin; } // Element -> [begin, end) - range_type range(const index_type element) const { - if (element >= this->size()) { - throw std::out_of_range("element out of range"); + range_type segment_range(const index_type segment) const { + if (segment >= this->size()) { + throw std::out_of_range("segment out of range"); } - const index_type begin = this->_offsets[element]; - const index_type end = this->_offsets[element + 1]; + const index_type begin = this->_offsets[segment]; + const index_type end = this->_offsets[segment + 1]; return {begin, end}; } - index_type element_size(const index_type element) const { - const range_type r = this->range(element); + index_type segment_size(const index_type segment) const { + const range_type r = this->segment_range(segment); return r.end - r.begin; } - // Buffer index -> owning element (and its [begin,end)). + // Buffer index -> owning segment (and its [begin,end)). locate_result locate(const index_type buffer_index) const { if (buffer_index >= this->total_size()) { throw std::out_of_range("buffer index out of range"); @@ -128,11 +189,13 @@ class OffsetTable { this->_offsets.end(), buffer_index); - const index_type element = static_cast(std::distance(this->_offsets.begin(), it) - 1); + const index_type segment = static_cast(std::distance(this->_offsets.begin(), it) - 1); - return {element, this->range(element)}; + return {segment, this->segment_range(segment)}; } private: std::vector _offsets; }; + +using OffsetTable = OffsetTable_; diff --git a/src/terrainlib/Range.h b/src/terrainlib/Range.h index eb0000d6..128947f7 100644 --- a/src/terrainlib/Range.h +++ b/src/terrainlib/Range.h @@ -1,110 +1,11 @@ #pragma once -#include -#include -#include -#include -#include - -#include "number_utils.h" - -template -struct Range { - T min; - T max; - - constexpr Range() : Range(std::numeric_limits::max(), std::numeric_limits::lowest()) {} - constexpr Range(T value) : Range(Range::from_single(value)) {} - constexpr Range(T min_value, T max_value) : min(min_value), max(max_value) {} - - [[nodiscard]] constexpr bool empty() const noexcept { - return min >= max; - } - - [[nodiscard]] constexpr T size() const noexcept { - static_assert(std::is_arithmetic_v, "Range::size requires arithmetic T"); - return this->empty() ? T{0} : (max - min); - } - - [[nodiscard]] constexpr bool valid() const noexcept { - return this->min <= this->max; - } - - [[nodiscard]] constexpr bool contains(T value) const noexcept { - return value >= this->min && value < this->max; - } - - [[nodiscard]] constexpr bool contains(const Range &other) const noexcept { - if (other.empty()) { - return true; - } - return other.min >= this->min && other.max <= this->max; - } - - [[nodiscard]] constexpr bool overlaps(const Range &other) const noexcept { - if (other.empty()) { - return false; - } - return this->min < other.max && other.min < this->max; - } - - constexpr void expand(T value) noexcept { - this->min = std::min(min, value); - this->max = std::max(max, value + T{1}); - } - - constexpr void expand(const Range &other) noexcept { - if (other.empty()) { - return; - } - this->min = std::min(this->min, other.min); - this->max = std::max(this->max, other.max); - } - - constexpr void clamp(const Range &bounds) noexcept { - this->min = std::clamp(this->min, bounds.min, bounds.max); - this->max = std::clamp(this->max, bounds.min, bounds.max); - if (this->min > this->max) { - this->min = this->max; - } - } - - constexpr void translate(T offset) noexcept { - this->min += offset; - this->max += offset; - } - - [[nodiscard]] constexpr Range intersection(const Range &other) const noexcept { - return { - std::max(this->min, other.min), - std::min(this->max, other.max)}; - } - - [[nodiscard]] constexpr Range unite(const Range &other) const noexcept { - if (this->empty()) { - return other; - } - if (other.empty()) { - return *this; - } - return { - std::min(this->min, other.min), - std::max(this->max, other.max)}; - } - - [[nodiscard]] static constexpr Range from_single(T value) noexcept { - return {value, next_higher(value)}; - } - - [[nodiscard]] constexpr bool operator==(const Range &) const noexcept = default; -}; - -template -constexpr Range full_range() { - return Range(std::numeric_limits::lowest(), std::numeric_limits::max()); -} - -template -constexpr Range empty_range() { - return Range(std::numeric_limits::max(), std::numeric_limits::lowest()); -} +#include "range/AnyRange.h" +#include "range/Bound.h" +#include "range/Range.h" +#include "range/RangeBounds.h" +#include "range/RangeFrom.h" +#include "range/RangeFull.h" +#include "range/RangeInclusive.h" +#include "range/RangeTo.h" +#include "range/RangeToInclusive.h" diff --git a/src/terrainlib/SegmentedBuffer.h b/src/terrainlib/SegmentedBuffer.h index 292909af..84011fbc 100644 --- a/src/terrainlib/SegmentedBuffer.h +++ b/src/terrainlib/SegmentedBuffer.h @@ -18,7 +18,7 @@ class SegmentedBuffer { using value_type = TValue; using index_type = TValueIndex; using segment_index = TSegmentIndex; - using offset_range = typename OffsetTable::range_type; + using offset_range = typename OffsetTable_::range_type; SegmentedBuffer() { this->reset(); @@ -61,7 +61,7 @@ class SegmentedBuffer { } const segment_index last_segment = this->segment_count() - 1; - const index_type new_end = this->_offsets.range(last_segment).begin + size; + const index_type new_end = this->_offsets.segment_range(last_segment).begin + size; this->_offsets.set_end(last_segment, new_end); this->_data.resize(new_end, value); } @@ -88,7 +88,7 @@ class SegmentedBuffer { // Appends an item to the end of the last segment in the data buffer. void push_to_last_segment(const value_type &value) { const segment_index last_segment = this->segment_count() - 1; - const offset_range range = this->_offsets.range(last_segment); + const offset_range range = this->_offsets.segment_range(last_segment); this->_data.push_back(value); this->_offsets.set_end(last_segment, range.end + 1); @@ -97,7 +97,7 @@ class SegmentedBuffer { // Appends an item to the end of the last segment in the data buffer. void push_to_last_segment(value_type &&value) { const segment_index last_segment = this->segment_count() - 1; - const offset_range range = this->_offsets.range(last_segment); + const offset_range range = this->_offsets.segment_range(last_segment); this->_data.push_back(std::forward(value)); this->_offsets.set_end(last_segment, range.end + 1); @@ -114,6 +114,12 @@ class SegmentedBuffer { this->_offsets.append_length(0); } + // Removes all empty segments while preserving the order of non-empty segments. + // If all segments are empty, the buffer is reset to one empty segment. + void remove_empty_segments() { + this->_offsets.remove_empty_segments(); + } + // Accesses an element using segment-relative indexing. value_type &operator()(const segment_index segment_index, const index_type element_index) { return this->get_impl(*this, segment_index, element_index); @@ -131,7 +137,7 @@ class SegmentedBuffer { // Returns the number of elements contained within a specific segment. index_type segment_size(const segment_index segment_index) const noexcept { - return this->_offsets.element_size(segment_index); + return this->_offsets.segment_size(segment_index); } // Returns the total number of elements across all segments. @@ -140,32 +146,39 @@ class SegmentedBuffer { } // Returns a readonly view of one segment of the buffer. - std::span segment(const segment_index segment_index) const noexcept { + std::span segment(const segment_index segment_index) const & noexcept { return this->get_segment_impl(*this, segment_index); } // Returns a view of one segment of the buffer. - std::span segment(const segment_index segment_index) noexcept { + std::span segment(const segment_index segment_index) & noexcept { return this->get_segment_impl(*this, segment_index); } // Returns a readonly view of the last segment of the buffer. - std::span last_segment() const noexcept { + std::span last_segment() const & noexcept { return this->segment(this->segment_count() - 1); } // Returns a view of the last segment of the buffer. - std::span last_segment() noexcept { + std::span last_segment() & noexcept { return this->segment(this->segment_count() - 1); } // Returns a readonly flat view of the buffer. - std::span flat() const noexcept { + std::span flat() const & noexcept { return this->_data; } // Returns a flat view of the buffer. - std::span flat() noexcept { + std::span flat() & noexcept { return this->_data; } + auto segments() const & noexcept { + return this->segments_impl(*this); + } + auto segments() & noexcept { + return this->segments_impl(*this); + } + std::vector& backing() noexcept { return this->_data; } @@ -178,7 +191,7 @@ class SegmentedBuffer { static inline auto &get_impl(Self &self, const segment_index segment_index, const index_type element_index) { DEBUG_ASSERT(segment_index < self._offsets.size()); - const offset_range range = self._offsets.range(segment_index); + const offset_range range = self._offsets.segment_range(segment_index); const index_type flat_index = range.begin + element_index; DEBUG_ASSERT(flat_index < range.end); @@ -189,13 +202,21 @@ class SegmentedBuffer { static inline auto get_segment_impl(Self &self, const segment_index segment_index) { DEBUG_ASSERT(segment_index < self._offsets.size()); - const offset_range range = self._offsets.range(segment_index); + const offset_range range = self._offsets.segment_range(segment_index); const index_type length = range.end - range.begin; return std::span(self._data.data() + range.begin, length); } + template + static auto segments_impl(Self &self) noexcept { + return std::views::iota(segment_index{0}, self.segment_count()) + | std::views::transform([&self](const segment_index i) noexcept { + return self.segment(i); + }); + } + // Manages the boundaries of each segment. - OffsetTable _offsets; + OffsetTable_ _offsets; // The contiguous backing storage for all segments. std::vector _data; diff --git a/src/terrainlib/UnionFind.h b/src/terrainlib/UnionFind.h index fed28461..fa40155e 100644 --- a/src/terrainlib/UnionFind.h +++ b/src/terrainlib/UnionFind.h @@ -1,13 +1,15 @@ #pragma once #include -#include #include -#include #include +#include #include +#include "SegmentedBuffer.h" + + template class UnionFind_ { public: @@ -38,12 +40,12 @@ class UnionFind_ { return this->find_impl(*this, item); } - void make_union(const Index x, const Index y) noexcept { + Index make_union(const Index x, const Index y) noexcept { const Index x_rep = this->find(x); const Index y_rep = this->find(y); if (x_rep == y_rep) { - return; + return x_rep; } this->_parents[x_rep] = y_rep; @@ -51,7 +53,17 @@ class UnionFind_ { if constexpr (TrackSizes) { this->_sizes[y_rep] += this->_sizes[x_rep]; + this->_sizes[x_rep] = Size{0}; } + + return y_rep; + } + + bool is_same_set(const Index x, const Index y) const noexcept { + return this->is_same_set_impl(*this, x, y); + } + bool is_same_set(const Index x, const Index y) noexcept { + return this->is_same_set_impl(*this, x, y); } [[nodiscard]] Size size() const noexcept { @@ -75,12 +87,26 @@ class UnionFind_ { return !this->is_joint(); } - [[nodiscard]] std::unordered_map> get_sets() const { - return this->get_sets_impl(*this); + [[nodiscard]] std::unordered_map> get_sets_as_map() const { + return this->get_sets_as_map_impl(*this); } - [[nodiscard]] std::unordered_map> get_sets() { - return this->get_sets_impl(*this); + [[nodiscard]] std::unordered_map> get_sets_as_map() { + return this->get_sets_as_map_impl(*this); + } + + [[nodiscard]] SegmentedBuffer get_sets_compact() const { + return this->get_sets_compact_impl(*this); + } + [[nodiscard]] SegmentedBuffer get_sets_compact() { + return this->get_sets_compact_impl(*this); + } + + [[nodiscard]] SegmentedBuffer get_sets_sparse() const { + return this->get_sets_sparse_impl(*this); + } + [[nodiscard]] SegmentedBuffer get_sets_sparse() { + return this->get_sets_sparse_impl(*this); } private: @@ -106,7 +132,14 @@ class UnionFind_ { } template - static std::unordered_map> get_sets_impl(Self &self) { + static bool is_same_set_impl(Self &self, const Index x, const Index y) noexcept { + const Index x_rep = self.find(x); + const Index y_rep = self.find(y); + return x_rep == y_rep; + } + + template + static std::unordered_map> get_sets_as_map_impl(Self &self) { std::unordered_map> sets; const Size n = self.size(); @@ -119,6 +152,49 @@ class UnionFind_ { return sets; } + template + static SegmentedBuffer get_sets_compact_impl(Self &self) { + SegmentedBuffer sparse = get_sets_sparse_impl(self); + sparse.remove_empty_segments(); + return sparse; + } + + template + static SegmentedBuffer get_sets_sparse_impl(Self &self) { + const Size n = self.size(); + SegmentedBuffer sets; + + if constexpr (TrackSizes) { + sets.init(std::span(self._sizes)); + + std::vector cursors(n, Size{0}); + for (Index item = 0; item < n; item++) { + const Index rep = self.find(item); + std::span set = sets.segment(rep); + Size &cursor = cursors[rep]; + set[cursor] = item; + cursor++; + } + } else { + std::vector counts(n, Size{0}); + for (Index item = 0; item < n; item++) { + counts[self.find(item)]++; + } + sets.init(std::span(counts)); + + for (Index item = 0; item < n; item++) { + const Index rep = self.find(item); + std::span set = sets.segment(rep); + Size& remaining_count = counts[rep]; + Size cursor = set.size() - remaining_count; + set[cursor] = item; + remaining_count--; + } + } + + return sets; + } + private: std::vector _parents; std::conditional_t, std::monostate> _sizes; diff --git a/src/terrainlib/int_math.h b/src/terrainlib/int_math.h index 9ea61a46..a7dec2fe 100644 --- a/src/terrainlib/int_math.h +++ b/src/terrainlib/int_math.h @@ -3,6 +3,7 @@ #include #include #include +#include #include @@ -40,4 +41,34 @@ template template [[nodiscard]] constexpr bool is_odd(T value) noexcept { return !is_even(value); -} \ No newline at end of file +} + +template + requires(std::is_signed_v == std::is_signed_v) +[[nodiscard]] constexpr auto int_div_ceil(A a, B b) { + using T = std::common_type_t; + ASSERT(b != 0); + + const T a_t = static_cast(a); + const T b_t = static_cast(b); + + if constexpr (std::unsigned_integral) { + return (a_t + b_t - 1) / b_t; + } else { + const T q = a_t / b_t; + const T r = a_t % b_t; + + // Exact division -> already the ceiling + if (r == 0) { + return q; + } + + // Same sign + if ((a_t > 0 && b_t > 0) || (a_t < 0 && b_t < 0)) { + return q + 1; + } + + // Opposite signs -> truncation toward zero already gave the ceiling + return q; + } +} diff --git a/src/terrainlib/mesh/VertexMap.h b/src/terrainlib/mesh/VertexMap.h index 978bc064..d11f4679 100644 --- a/src/terrainlib/mesh/VertexMap.h +++ b/src/terrainlib/mesh/VertexMap.h @@ -25,7 +25,7 @@ class VertexMap { } static VertexMap identity(const std::span triangles) { const auto range = mesh::find_vertex_index_range(triangles); - return VertexMap::identity(range.size(), range.min); + return VertexMap::identity(range.size(), range.start); } static VertexMap from_forward(std::vector forward, const uint32_t offset = 0) { return from_forward(OffsetVector(std::move(forward), offset)); diff --git a/src/terrainlib/mesh/manifold.inl b/src/terrainlib/mesh/manifold.inl index 292300d6..0f4671d7 100644 --- a/src/terrainlib/mesh/manifold.inl +++ b/src/terrainlib/mesh/manifold.inl @@ -157,7 +157,7 @@ void duplicate_non_manifold_vertices( } // Collect connected face fans - std::unordered_map> face_fans = union_find.get_sets(); + std::unordered_map> face_fans = union_find.get_sets_as_map(); const auto largest_fan_it = std::max_element( face_fans.begin(), face_fans.end(), diff --git a/src/terrainlib/mesh/reindex.inl b/src/terrainlib/mesh/reindex.inl index 9eca4ded..658d901a 100644 --- a/src/terrainlib/mesh/reindex.inl +++ b/src/terrainlib/mesh/reindex.inl @@ -26,7 +26,7 @@ constexpr uint32_t invalid_index = VertexMap::invalid_index; inline VertexMap create_reindex_map(const std::span triangles) { const Range vertex_range = find_vertex_index_range(triangles); OffsetVector old_to_new; - old_to_new.offset = vertex_range.min; + old_to_new.offset = vertex_range.start; old_to_new.resize(vertex_range.size(), detail::invalid_index); uint32_t next_index = 0; diff --git a/src/terrainlib/octree/OddLevelShifted.cpp b/src/terrainlib/octree/OddLevelShifted.cpp index 3f2543d1..41ee1cd8 100644 --- a/src/terrainlib/octree/OddLevelShifted.cpp +++ b/src/terrainlib/octree/OddLevelShifted.cpp @@ -101,13 +101,17 @@ IdRect OddLevelShifted::find_intersecting_standard_nodes(const Id &id) const { } Bounds OddLevelShifted::get_node_bounds_with_children(const Id &id) const { + Bounds bounds = this->get_node_bounds(id); if (!id.has_children()) { - return this->get_node_bounds(id); + return bounds; } - const Bounds first_bounds = this->get_node_bounds(id.child(0).value()); - const Bounds last_bounds = this->get_node_bounds(id.child(7).value()); - return Bounds(first_bounds.min, last_bounds.max); + // The node and its children live in different shifted regions (children are + // one level finer, so an odd node has unshifted children and vice versa). + // Union both so the bounds cover the node regardless of its level parity. + bounds.expand_by(this->get_node_bounds(id.child(0).value())); + bounds.expand_by(this->get_node_bounds(id.child(7).value())); + return bounds; } Bounds OddLevelShifted::get_node_bounds(const Id &id) const { diff --git a/src/terrainlib/opencv_utils.h b/src/terrainlib/opencv_utils.h index 04128e0c..4b527935 100644 --- a/src/terrainlib/opencv_utils.h +++ b/src/terrainlib/opencv_utils.h @@ -3,6 +3,10 @@ #include #include +inline size_t mat_byte_size(const cv::Mat &mat) { + return mat.total() * mat.elemSize(); +} + inline void rescale_texture(const cv::Mat &source, cv::Mat &destination, const glm::uvec2 new_size) { cv::resize(source, destination, cv::Size(new_size.x, new_size.y), 0, 0); } diff --git a/src/terrainlib/range/AnyRange.h b/src/terrainlib/range/AnyRange.h new file mode 100644 index 00000000..91e94b40 --- /dev/null +++ b/src/terrainlib/range/AnyRange.h @@ -0,0 +1,49 @@ +#pragma once + +#include +#include + +#include "Bound.h" +#include "Range.h" +#include "RangeBounds.h" +#include "RangeFull.h" + +template +struct AnyRange : RangeBounds, T> { + Bound start; + Bound end; + + constexpr AnyRange(Bound start_bound, Bound end_bound) noexcept + : start(std::move(start_bound)), end(std::move(end_bound)) {} + + constexpr AnyRange(std::pair, Bound> bounds) noexcept + : start(std::move(bounds.first)), end(std::move(bounds.second)) {} + + template + constexpr AnyRange(const RangeBounds &range) noexcept + : start(static_cast(range).start_bound()), + end(static_cast(range).end_bound()) {} + + constexpr AnyRange(RangeFull range) noexcept + : start(range.template start_bound()), end(range.template end_bound()) {} + + [[nodiscard]] constexpr Bound start_bound() const noexcept { + return this->start; + } + + [[nodiscard]] constexpr Bound end_bound() const noexcept { + return this->end; + } + + // Resolves an Unbounded start to 0 and an Unbounded end to len, mirroring how Rust + // resolves RangeFull/RangeTo/RangeFrom against a slice's length. + [[nodiscard]] constexpr Range to_range(T len) const noexcept + requires std::unsigned_integral + { + const T range_start = this->start.kind == BoundKind::Unbounded ? T{0} : *this->start.value; + const T range_end = this->end.kind == BoundKind::Unbounded + ? len + : (this->end.kind == BoundKind::Included ? *this->end.value + 1 : *this->end.value); + return {range_start, range_end}; + } +}; diff --git a/src/terrainlib/range/Bound.h b/src/terrainlib/range/Bound.h new file mode 100644 index 00000000..0207925a --- /dev/null +++ b/src/terrainlib/range/Bound.h @@ -0,0 +1,30 @@ +#pragma once + +#include +#include + +enum class BoundKind { + Included, + Excluded, + Unbounded, +}; + +template +struct Bound { + BoundKind kind; + std::optional value; + + [[nodiscard]] static constexpr Bound included(T value) noexcept { + return {BoundKind::Included, std::move(value)}; + } + + [[nodiscard]] static constexpr Bound excluded(T value) noexcept { + return {BoundKind::Excluded, std::move(value)}; + } + + [[nodiscard]] static constexpr Bound unbounded() noexcept { + return {BoundKind::Unbounded, std::nullopt}; + } + + [[nodiscard]] constexpr bool operator==(const Bound &) const noexcept = default; +}; diff --git a/src/terrainlib/range/Range.h b/src/terrainlib/range/Range.h new file mode 100644 index 00000000..344561a8 --- /dev/null +++ b/src/terrainlib/range/Range.h @@ -0,0 +1,40 @@ +#pragma once + +#include + +#include "Bound.h" +#include "RangeBounds.h" +#include "number_utils.h" + +template +struct Range : RangeBounds, T> { + T start; + T end; + + constexpr Range() : Range({}, {}) {} + constexpr Range(T value) requires std::integral || std::floating_point + : Range(value, next_higher(value)) {} + constexpr Range(T start_value, T end_value) : start(start_value), end(end_value) {} + + [[nodiscard]] constexpr Bound start_bound() const noexcept { + return Bound::included(this->start); + } + + [[nodiscard]] constexpr Bound end_bound() const noexcept { + return Bound::excluded(this->end); + } + + [[nodiscard]] constexpr T size() const noexcept { + return this->end - this->start; + } + + [[nodiscard]] constexpr bool is_in_bounds(const T &len) const noexcept { + return this->start <= this->end && this->end <= len; + } + + [[nodiscard]] constexpr bool is_overlapping(const Range &other) const noexcept { + return this->start < other.end && other.start < this->end; + } + + [[nodiscard]] constexpr bool operator==(const Range &) const noexcept = default; +}; diff --git a/src/terrainlib/range/RangeBounds.h b/src/terrainlib/range/RangeBounds.h new file mode 100644 index 00000000..a9d1cb94 --- /dev/null +++ b/src/terrainlib/range/RangeBounds.h @@ -0,0 +1,82 @@ +#pragma once + +#include + +#include "Bound.h" + +struct RangeFull; + +template +struct AnyRange; + +template +[[nodiscard]] constexpr Bound tighter_start_bound(const Bound &a, const Bound &b) noexcept { + if (a.kind == BoundKind::Unbounded) { + return b; + } + if (b.kind == BoundKind::Unbounded) { + return a; + } + if (*a.value != *b.value) { + return *a.value > *b.value ? a : b; + } + return a.kind == BoundKind::Excluded ? a : b; +} + +template +[[nodiscard]] constexpr Bound tighter_end_bound(const Bound &a, const Bound &b) noexcept { + if (a.kind == BoundKind::Unbounded) { + return b; + } + if (b.kind == BoundKind::Unbounded) { + return a; + } + if (*a.value != *b.value) { + return *a.value < *b.value ? a : b; + } + return a.kind == BoundKind::Excluded ? a : b; +} + +template +struct RangeBounds { + [[nodiscard]] constexpr bool contains(const T &item) const noexcept { + const auto &self = static_cast(*this); + const Bound start = self.start_bound(); + const Bound end = self.end_bound(); + + const bool after_start = start.kind == BoundKind::Unbounded || + (start.kind == BoundKind::Included ? *start.value <= item : *start.value < item); + const bool before_end = end.kind == BoundKind::Unbounded || + (end.kind == BoundKind::Included ? item <= *end.value : item < *end.value); + + return after_start && before_end; + } + + [[nodiscard]] constexpr bool is_empty() const noexcept { + const auto &self = static_cast(*this); + const Bound start = self.start_bound(); + const Bound end = self.end_bound(); + + if (start.kind == BoundKind::Unbounded || end.kind == BoundKind::Unbounded) { + return false; + } + + return end.kind == BoundKind::Included ? !(*start.value <= *end.value) : !(*start.value < *end.value); + } + + template + [[nodiscard]] constexpr AnyRange intersect(const RangeBounds &other) const noexcept { + const auto &self = static_cast(*this); + const auto &other_self = static_cast(other); + return AnyRange( + tighter_start_bound(self.start_bound(), other_self.start_bound()), + tighter_end_bound(self.end_bound(), other_self.end_bound())); + } + + [[nodiscard]] constexpr AnyRange intersect(const RangeFull &) const noexcept { + const auto &self = static_cast(*this); + return AnyRange(self.start_bound(), self.end_bound()); + } + + [[nodiscard]] constexpr bool operator==(const RangeBounds &) const noexcept = default; +}; diff --git a/src/terrainlib/range/RangeFrom.h b/src/terrainlib/range/RangeFrom.h new file mode 100644 index 00000000..7bb1ca77 --- /dev/null +++ b/src/terrainlib/range/RangeFrom.h @@ -0,0 +1,19 @@ +#pragma once + +#include "Bound.h" +#include "RangeBounds.h" + +template +struct RangeFrom : RangeBounds, T> { + T start; + + [[nodiscard]] constexpr Bound start_bound() const noexcept { + return Bound::included(this->start); + } + + [[nodiscard]] constexpr Bound end_bound() const noexcept { + return Bound::unbounded(); + } + + [[nodiscard]] constexpr bool operator==(const RangeFrom &) const noexcept = default; +}; diff --git a/src/terrainlib/range/RangeFull.h b/src/terrainlib/range/RangeFull.h new file mode 100644 index 00000000..5da536dc --- /dev/null +++ b/src/terrainlib/range/RangeFull.h @@ -0,0 +1,27 @@ +#pragma once + +#include "Bound.h" + +struct RangeFull { + template + [[nodiscard]] constexpr Bound start_bound() const noexcept { + return Bound::unbounded(); + } + + template + [[nodiscard]] constexpr Bound end_bound() const noexcept { + return Bound::unbounded(); + } + + template + [[nodiscard]] constexpr bool contains(const T &) const noexcept { + return true; + } + + template + [[nodiscard]] constexpr bool is_empty() const noexcept { + return false; + } + + [[nodiscard]] constexpr bool operator==(const RangeFull &) const noexcept = default; +}; diff --git a/src/terrainlib/range/RangeInclusive.h b/src/terrainlib/range/RangeInclusive.h new file mode 100644 index 00000000..baa247a4 --- /dev/null +++ b/src/terrainlib/range/RangeInclusive.h @@ -0,0 +1,23 @@ +#pragma once + +#include "Bound.h" +#include "RangeBounds.h" + +template +struct RangeInclusive : RangeBounds, T> { + T start; + T end; + + constexpr RangeInclusive(T value) : RangeInclusive(value, value) {} + constexpr RangeInclusive(T start_value, T end_value) : start(start_value), end(end_value) {} + + [[nodiscard]] constexpr Bound start_bound() const noexcept { + return Bound::included(this->start); + } + + [[nodiscard]] constexpr Bound end_bound() const noexcept { + return Bound::included(this->end); + } + + [[nodiscard]] constexpr bool operator==(const RangeInclusive &) const noexcept = default; +}; diff --git a/src/terrainlib/range/RangeTo.h b/src/terrainlib/range/RangeTo.h new file mode 100644 index 00000000..ccdc7fc4 --- /dev/null +++ b/src/terrainlib/range/RangeTo.h @@ -0,0 +1,19 @@ +#pragma once + +#include "Bound.h" +#include "RangeBounds.h" + +template +struct RangeTo : RangeBounds, T> { + T end; + + [[nodiscard]] constexpr Bound start_bound() const noexcept { + return Bound::unbounded(); + } + + [[nodiscard]] constexpr Bound end_bound() const noexcept { + return Bound::excluded(this->end); + } + + [[nodiscard]] constexpr bool operator==(const RangeTo &) const noexcept = default; +}; diff --git a/src/terrainlib/range/RangeToInclusive.h b/src/terrainlib/range/RangeToInclusive.h new file mode 100644 index 00000000..71de400b --- /dev/null +++ b/src/terrainlib/range/RangeToInclusive.h @@ -0,0 +1,19 @@ +#pragma once + +#include "Bound.h" +#include "RangeBounds.h" + +template +struct RangeToInclusive : RangeBounds, T> { + T end; + + [[nodiscard]] constexpr Bound start_bound() const noexcept { + return Bound::unbounded(); + } + + [[nodiscard]] constexpr Bound end_bound() const noexcept { + return Bound::included(this->end); + } + + [[nodiscard]] constexpr bool operator==(const RangeToInclusive &) const noexcept = default; +}; diff --git a/src/terrainlib/range_utils.h b/src/terrainlib/range_utils.h index bf42de67..a5ff24e9 100644 --- a/src/terrainlib/range_utils.h +++ b/src/terrainlib/range_utils.h @@ -138,6 +138,7 @@ constexpr T sum(Range &&range, T init, F &&f) { } template + requires std::invocable> constexpr auto sum(Range &&range, F &&f) { using T = std::decay_t>>; return sum(std::forward(range), T{}, std::forward(f)); @@ -149,6 +150,26 @@ constexpr T sum(Range &&range, T init = T{}) { return sum(std::forward(range), init, std::identity{}); } +// Minimum element of a range, or the given default value if the range is empty. +template > + requires std::convertible_to, T> +constexpr T min(Range &&range, T default_value = T{}) { + if (std::ranges::empty(range)) { + return default_value; + } + return std::ranges::min(std::forward(range)); +} + +// Maximum element of a range, or the given default value if the range is empty. +template > + requires std::convertible_to, T> +constexpr T max(Range &&range, T default_value = T{}) { + if (std::ranges::empty(range)) { + return default_value; + } + return std::ranges::max(std::forward(range)); +} + template constexpr auto range(T begin, T end) { return std::views::iota(begin, end); diff --git a/src/terrainlib/spatial_lookup/CellBased.h b/src/terrainlib/spatial_lookup/CellBased.h index e29df1d4..e834ebfb 100644 --- a/src/terrainlib/spatial_lookup/CellBased.h +++ b/src/terrainlib/spatial_lookup/CellBased.h @@ -5,7 +5,6 @@ #include #include -#include "MoreExact.h" #include "log.h" #include "spatial_lookup/CellBasedStorage.h" #include "spatial_lookup/NDLoopHelper.h" @@ -17,7 +16,7 @@ namespace detail { // Calculate distance^2, works with all types template auto distance_sq(const glm::vec &a, const glm::vec &b) { - using T = MoreExact; + using T = std::common_type_t; using Vec = glm::vec; if constexpr (std::is_floating_point_v) { @@ -116,6 +115,16 @@ class CellBased { bool find_all_at(const Vec &point, Vector &out) { return find_all_at_impl(*this, point, out); } + + template + bool for_all_points(Func &&func) const { + return this->_storage.for_all_points(std::forward(func)); + } + + template + bool for_all_points(Func &&func) { + return this->_storage.for_all_points(std::forward(func)); + } private: Storage _storage; @@ -131,8 +140,8 @@ class CellBased { std::is_const_v, const Value &, Value &>; - // Promote to the more precise type to avoid e.g. float epsilon against a double grid - using E = MoreExact; + // Promote to the more precise type + using E = std::common_type_t; DEBUG_ASSERT(epsilon > 0); @@ -250,7 +259,8 @@ class CellBased { std::is_const_v, const Value &, Value &>; - using E = MoreExact; + // Promote to the more precise type + using E = std::common_type_t; out.clear(); diff --git a/src/terrainlib/spatial_lookup/CellBasedStorage.h b/src/terrainlib/spatial_lookup/CellBasedStorage.h index 60117ade..57ddd9ba 100644 --- a/src/terrainlib/spatial_lookup/CellBasedStorage.h +++ b/src/terrainlib/spatial_lookup/CellBasedStorage.h @@ -35,6 +35,14 @@ concept CellBasedStorage = requires( { const_storage.for_all_in_cell(index, [](const glm::vec&, const Value&) {}) } -> std::same_as; + + { + storage.for_all_points([](const glm::vec &, Value &) {}) + } -> std::same_as; + + { + const_storage.for_all_points([](const glm::vec &, const Value &) {}) + } -> std::same_as; }; } // namespace spatial_lookup diff --git a/src/terrainlib/spatial_lookup/Grid.h b/src/terrainlib/spatial_lookup/Grid.h index f0e5164a..9f86ca38 100644 --- a/src/terrainlib/spatial_lookup/Grid.h +++ b/src/terrainlib/spatial_lookup/Grid.h @@ -14,42 +14,3 @@ template using Grid3d = Grid<3, double, Value>; } - -/* -namespace { -radix::geometry::Aabb3d pad_bounds(const radix::geometry::Aabb3d &bounds, const double percentage) { - const glm::dvec3 bounds_padding = bounds.size() * percentage; - const radix::geometry::Aabb3d padded_bounds(bounds.min - bounds_padding, bounds.max + bounds_padding); - return padded_bounds; -} - -template -Grid _construct_grid_for_meshes(const radix::geometry::Aabb3d &bounds, const size_t vertex_count) { - const radix::geometry::Aabb3d padded_bounds = pad_bounds(bounds, 0.01); - - const double max_extends = max_component(padded_bounds.size()); - const glm::dvec3 relative_extends = padded_bounds.size() / max_extends; - const glm::uvec3 grid_divisions = glm::max(glm::uvec3(2 * std::cbrt(vertex_count) * relative_extends), glm::uvec3(1)); - Grid grid(padded_bounds.min, padded_bounds.size(), grid_divisions); - - return grid; -} -} // namespace - -template -inline Grid construct_grid_for_mesh(const SimpleMesh_ &mesh) { - const radix::geometry::Aabb3d bounds = calculate_bounds(mesh); - const size_t vertex_count = mesh.vertex_count(); - return _construct_grid_for_meshes(bounds, vertex_count); -} - -template -inline Grid construct_grid_for_meshes(const std::span> meshes) { - const radix::geometry::Aabb3d bounds = calculate_bounds(meshes); - const size_t maximal_merged_mesh_size = std::transform_reduce( - meshes.begin(), meshes.end(), 0, - [](const size_t a, const size_t b) { return a + b; }, - [](const SimpleMesh_ &mesh) { return mesh.vertex_count(); }); - return _construct_grid_for_meshes(bounds, maximal_merged_mesh_size); -} -*/ \ No newline at end of file diff --git a/src/terrainlib/spatial_lookup/GridStorage.h b/src/terrainlib/spatial_lookup/GridStorage.h index 23ba9c56..7df68f90 100644 --- a/src/terrainlib/spatial_lookup/GridStorage.h +++ b/src/terrainlib/spatial_lookup/GridStorage.h @@ -135,22 +135,22 @@ class GridStorage { template bool for_all_in_cell(const CellIndex index, Func &&func) const { - if (!this->is_valid_cell_index(index)) { - return false; - } + return for_all_in_cell_impl(*this, index, std::forward(func)); + } - const Cell &cell = this->_data[index]; - for (const Item &item : cell.items) { - func(item.point, item.value); - } + template + bool for_all_in_cell(const CellIndex index, Func &&func) { + return for_all_in_cell_impl(*this, index, std::forward(func)); + } - return !cell.items.empty(); + template + bool for_all_points(Func &&func) const { + return for_all_points_impl(*this, std::forward(func)); } + template - bool for_all_in_cell(const CellIndex index, Func &&func) { - return const_cast(this)->for_all_in_cell(index, [&](const Vec &vec, const Value &value) { - func(vec, const_cast(value)); - }); + bool for_all_points(Func &&func) { + return for_all_points_impl(*this, std::forward(func)); } [[nodiscard]] Vec cell_size() const { @@ -168,6 +168,33 @@ class GridStorage { Vec _size; glm::vec _divisions; + template + static bool for_all_in_cell_impl(SelfT &self, const CellIndex index, Func &&func) { + using ValueRef = std::conditional_t, const Value &, Value &>; + + if (!self.is_valid_cell_index(index)) { + return false; + } + + auto &cell = self._data[index]; + for (auto &item : cell.items) { + func(item.point, static_cast(item.value)); + } + return !cell.items.empty(); + } + + template + static bool for_all_points_impl(SelfT &self, Func &&func) { + using ValueRef = std::conditional_t, const Value &, Value &>; + + for (auto &cell : self._data) { + for (auto &item : cell.items) { + func(item.point, static_cast(item.value)); + } + } + return self._point_count > 0; + } + // static_assert(CellBasedStorage); }; diff --git a/src/terrainlib/spatial_lookup/HashmapPointMap.h b/src/terrainlib/spatial_lookup/HashmapPointMap.h deleted file mode 100644 index a576e2d0..00000000 --- a/src/terrainlib/spatial_lookup/HashmapPointMap.h +++ /dev/null @@ -1,61 +0,0 @@ -#pragma once - -#include -#include - -#include - -#include "spatial_lookup/PointMap.h" -#include "spatial_lookup/QuantizedHash.h" - -namespace spatial_lookup { - -// PointMap implementation backed by a quantized hashmap. -// At most one Value is stored per quantized cell; the first insert wins. -template -class HashmapPointMap : public PointMap { -public: - using Vec = glm::vec; - - explicit HashmapPointMap(Component epsilon) - : _epsilon(epsilon), _map(0, Hash(epsilon), Equal(epsilon)) {} - - Value find_or_insert(const Vec &point, Value value) override { - // try_emplace performs a single lookup; the first insert wins - const auto [it, inserted] = this->_map.try_emplace(point, std::move(value)); - return it->second; - } - - std::optional find(const Vec &point) const override { - const auto it = this->_map.find(point); - if (it == this->_map.end()) { - return std::nullopt; - } - return it->second; - } - - void clear() { - this->_map.clear(); - } - - [[nodiscard]] bool empty() const { - return this->_map.empty(); - } - - [[nodiscard]] size_t size() const { - return this->_map.size(); - } - - Component epsilon() const { - return this->_epsilon; - } - -private: - using Hash = detail::QuantizedVecHash; - using Equal = detail::QuantizedVecEqual; - - Component _epsilon; - std::unordered_map _map; -}; - -} // namespace spatial_lookup diff --git a/src/terrainlib/spatial_lookup/HashmapStorage.h b/src/terrainlib/spatial_lookup/HashmapStorage.h index 33af1576..c47a194a 100644 --- a/src/terrainlib/spatial_lookup/HashmapStorage.h +++ b/src/terrainlib/spatial_lookup/HashmapStorage.h @@ -68,18 +68,36 @@ class HashmapStorage { return for_all_in_cell_impl(*this, index, std::forward(func)); } + template + bool for_all_points(Func &&func) const { + return for_all_points_impl(*this, std::forward(func)); + } + + template + bool for_all_points(Func &&func) { + return for_all_points_impl(*this, std::forward(func)); + } + private: template static bool for_all_in_cell_impl(SelfT &self, const CellIndex &index, Func &&func) { auto [begin, end] = self._store.equal_range(index.quantized); bool any_found = false; - for (auto it = begin; it != end; ++it) { any_found = true; func(it->first, it->second); } + return any_found; + } + template + static bool for_all_points_impl(SelfT &self, Func &&func) { + bool any_found = false; + for (auto &[point, value] : self._store) { + any_found = true; + func(point, value); + } return any_found; } diff --git a/src/terrainlib/spatial_lookup/PointMap.h b/src/terrainlib/spatial_lookup/PointMap.h deleted file mode 100644 index 4f31ddc4..00000000 --- a/src/terrainlib/spatial_lookup/PointMap.h +++ /dev/null @@ -1,29 +0,0 @@ -#pragma once - -#include - -#include - -namespace spatial_lookup { - -// Maps spatial positions to a single canonical Value per quantized cell. -// The first value inserted for a given location wins; subsequent inserts of -// equivalent positions return that canonical value without overwriting it. -// "Equivalent" is defined by the concrete implementation (exact equality, -// epsilon radius, quantized bucket, etc.). -template -class PointMap { -public: - using Vec = glm::vec; - - virtual ~PointMap() = default; - - // If an equivalent point already exists, return its value. - // Otherwise insert (point -> value) and return value. - virtual Value find_or_insert(const Vec &point, Value value) = 0; - - // Returns the stored value for an equivalent point, or nullopt if none exists. - virtual std::optional find(const Vec &point) const = 0; -}; - -} // namespace spatial_lookup diff --git a/src/terrainlib/string_utils.h b/src/terrainlib/string_utils.h index 435eddbf..f2de1635 100644 --- a/src/terrainlib/string_utils.h +++ b/src/terrainlib/string_utils.h @@ -1,11 +1,28 @@ #pragma once +#include #include +#include #include #include #include #include +// Format a byte count as a human-readable string, e.g. "1.5 GiB". +inline std::string human_byte_size(const size_t bytes) { + constexpr std::array units = {"B", "KiB", "MiB", "GiB", "TiB"}; + double value = static_cast(bytes); + size_t unit = 0; + while (value >= 1024.0 && unit + 1 < units.size()) { + value /= 1024.0; + unit++; + } + std::array buffer; + const char *format = unit == 0 ? "%.0f %s" : "%.1f %s"; + std::snprintf(buffer.data(), buffer.size(), format, value, units[unit]); + return std::string(buffer.data()); +} + template std::optional from_chars(const std::basic_string_view sv) { static_assert(std::is_integral_v, "T must be an integral type"); diff --git a/src/tile_downloader/TileDownloader.h b/src/tile_downloader/TileDownloader.h index ceb5103b..4678376c 100644 --- a/src/tile_downloader/TileDownloader.h +++ b/src/tile_downloader/TileDownloader.h @@ -16,18 +16,33 @@ class TileDownloader { public: TileDownloader(const TileUrlBuilder &url_builder, std::string output_pattern, - bool early_skip, std::optional max_zoom_level) + bool early_skip, std::optional max_zoom_level, unsigned root_zoom_level) : _url_builder(url_builder), _output_pattern(std::move(output_pattern)), + _logger(root_zoom_level), _early_skip(early_skip), _max_zoom_level(max_zoom_level) {} void download_recursive(const radix::tile::Id &root_id) { - this->_logger.progress(root_id, "Connecting..."); + this->_logger.start(); + this->download_recursive_core(root_id); + this->_logger.finish(); + } + +private: + const TileUrlBuilder &_url_builder; + std::string _output_pattern; + HttpClient _http; + TileLogger _logger; + bool _early_skip; + std::optional _max_zoom_level; + + void download_recursive_core(const radix::tile::Id &root_id) { auto result = this->download_tile(root_id); - this->_logger.result(root_id, result); + this->_logger.report_error(root_id, result); if (is_failure(result)) { + this->_logger.missing(root_id); return; } @@ -45,18 +60,10 @@ class TileDownloader { continue; } - this->download_recursive(children[i]); + this->download_recursive_core(children[i]); } } -private: - const TileUrlBuilder &_url_builder; - std::string _output_pattern; - HttpClient _http; - TileLogger _logger; - bool _early_skip; - std::optional _max_zoom_level; - static bool is_failure(const TileResult::Status &result) { return std::holds_alternative(result) || std::holds_alternative(result) @@ -112,12 +119,9 @@ class TileDownloader { ensure_parent_dirs(path); const std::string url = this->_url_builder.build_url(tile); - auto progress_fn = [&](double fraction) { - this->_logger.progress(tile, fraction); - }; for (int attempt = 0; attempt < 100; attempt++) { - HttpResponse response = this->_http.get(url, progress_fn); + HttpResponse response = this->_http.get(url); if (response.curl_code == CURLE_OK && this->_http.is_image(response)) { write_file(path, response.body); diff --git a/src/tile_downloader/TileLogger.h b/src/tile_downloader/TileLogger.h index 754143bc..3919e74c 100644 --- a/src/tile_downloader/TileLogger.h +++ b/src/tile_downloader/TileLogger.h @@ -1,15 +1,20 @@ #pragma once +#include +#include +#include +#include +#include #include #include +#include #include #include #include #include -#include - +#include "ProgressIndicator.h" #include "log.h" struct TileResult { @@ -26,55 +31,86 @@ struct TileResult { class TileLogger { public: - void progress(const radix::tile::Id &tile, std::string_view status) const { - if (Log::get_logger()->level() <= spdlog::level::info) { - fmt::print("\33[2K\r{} ({})", this->format(tile), status); - std::fflush(nullptr); - } + explicit TileLogger(unsigned root_zoom_level) + : _root_zoom_level(root_zoom_level), _progress(PROGRESS_RESOLUTION) {} + + // Starts the progress display; pair with a call to finish(). + void start() { + this->_progress_thread = this->_progress.start_monitoring(); } - void progress(const radix::tile::Id &tile, double fraction) const { - if (fraction < 0) { - this->progress(tile, "Downloading..."); - } else { - this->progress(tile, fmt::format("Downloading {:.0f}%...", fraction * 100)); - } + // The tile turned out not to exist (or errored out), so its branch stops here. + void missing(const radix::tile::Id &tile) { + this->close_branch(tile); } - void skipped(const radix::tile::Id &tile) const { - this->final(tile, "Skipped"); + // We chose not to descend into this tile's branch (already on disk, or past --max-zoom-level). + void skipped(const radix::tile::Id &tile) { + this->close_branch(tile); } - void result(const radix::tile::Id &tile, const TileResult::Status &status) const { + void report_error(const radix::tile::Id &tile, const TileResult::Status &status) const { std::visit([&](const auto &r) { using T = std::decay_t; - if constexpr (std::is_same_v) { - this->final(tile, "Done"); - } else if constexpr (std::is_same_v) { - this->final(tile, "Skipped"); - } else if constexpr (std::is_same_v) { - this->final(tile, "Absent"); - } else if constexpr (std::is_same_v) { - this->final(tile, fmt::format("Error HTTP [{}]", r.status_code)); + if constexpr (std::is_same_v) { + LOG_WARN("{}: HTTP error [{}]", format(tile), r.status_code); } else if constexpr (std::is_same_v) { - this->final(tile, "Error, bad content type."); + LOG_WARN("{}: bad content type", format(tile)); } else if constexpr (std::is_same_v) { - this->final(tile, fmt::format("Error cURL [{}] {}", unsigned(r.code), curl_easy_strerror(r.code))); + LOG_WARN("{}: cURL error [{}] {}", format(tile), unsigned(r.code), curl_easy_strerror(r.code)); } else if constexpr (std::is_same_v) { - this->final(tile, "Failed after 100 attempts"); + LOG_WARN("{}: failed after 100 attempts", format(tile)); } }, status); } + // Fills in any steps still missing due to floating point rounding, so the + // progress indicator reaches its total, then joins the thread started by start(). + void finish() { + while (this->_steps_done < PROGRESS_RESOLUTION) { + this->_progress.task_finished(); + this->_steps_done++; + } + if (this->_progress_thread.joinable()) { + this->_progress_thread.join(); + } + } + private: + static constexpr size_t PROGRESS_RESOLUTION = 10'000; + + unsigned _root_zoom_level; + ProgressIndicator _progress; + std::jthread _progress_thread; + size_t _steps_done = 0; + std::map _level_counts; + static std::string format(const radix::tile::Id &tile) { return fmt::format("Tile[Zoom={}, X={}, Y={}]", tile.zoom_level, tile.coords.x, tile.coords.y); } - void final(const radix::tile::Id &tile, std::string_view status) const { - if (Log::get_logger()->level() <= spdlog::level::info) { - fmt::print("\33[2K\r{} ({})\n", this->format(tile), status); - std::fflush(nullptr); + // Every tile at `level` covers an equal 4^-(level - root) share of the whole + // download tree, so summing that share for every closed branch gives the + // fraction of the tree that is done. + double completed_fraction() const { + double fraction = 0.0; + for (const auto &[level, count] : this->_level_counts) { + fraction += static_cast(count) / std::pow(4.0, static_cast(level - this->_root_zoom_level)); + } + return fraction; + } + + void close_branch(const radix::tile::Id &tile) { + this->_level_counts[tile.zoom_level]++; + this->advance_progress(); + } + + void advance_progress() { + const auto target = std::min(PROGRESS_RESOLUTION, + static_cast(this->completed_fraction() * double(PROGRESS_RESOLUTION))); + while (this->_steps_done < target) { + this->_progress.task_finished(); + this->_steps_done++; } } }; diff --git a/src/tile_downloader/main.cpp b/src/tile_downloader/main.cpp index 79102590..09c69d64 100644 --- a/src/tile_downloader/main.cpp +++ b/src/tile_downloader/main.cpp @@ -27,7 +27,7 @@ int main(int argc, char *argv[]) { const radix::tile::Id root_id = {args.zoom, {args.x, args.y}, args.scheme}; - TileDownloader downloader(*url_builder, args.output, args.early_skip, args.max_zoom_level); + TileDownloader downloader(*url_builder, args.output, args.early_skip, args.max_zoom_level, root_id.zoom_level); downloader.download_recursive(root_id); return 0; diff --git a/unittests/CMakeLists.txt b/unittests/CMakeLists.txt index cbf59ea8..720a4648 100644 --- a/unittests/CMakeLists.txt +++ b/unittests/CMakeLists.txt @@ -51,6 +51,7 @@ add_executable(unittests_terrainlib terrainlib/polygon_triangulate.cpp terrainlib/progress_indicator_test.cpp terrainlib/quantize.cpp + terrainlib/range.cpp terrainlib/range_utils.cpp terrainlib/spatial_lookup.cpp terrainlib/string_utils.cpp @@ -89,6 +90,9 @@ endif() if(TARGET dagbuilderlib) add_executable(unittests_dagbuilder catch2_helpers.h + dag_builder/fixed_point.cpp + dag_builder/merge_clusterings.cpp + dag_builder/pull_reproject_texture.cpp dag_builder/slice.cpp ) target_link_libraries(unittests_dagbuilder PRIVATE dagbuilderlib Catch2::Catch2WithMain) diff --git a/unittests/dag_builder/fixed_point.cpp b/unittests/dag_builder/fixed_point.cpp new file mode 100644 index 00000000..0d7af622 --- /dev/null +++ b/unittests/dag_builder/fixed_point.cpp @@ -0,0 +1,109 @@ +#include + +#include "../catch2_helpers.h" + +#include "fixed.h" +#include "fixed_point.h" + +TEST_CASE("round_div rounds half away from zero", "[dag_builder][fixed_point]") { + using Q = fixed::Q<16>; + + // Positive half rounds up + CHECK(Q::round_div(1, 2) == 1); + CHECK(Q::round_div(3, 2) == 2); + + // Negative half rounds away from zero (down in absolute value) + CHECK(Q::round_div(-1, 2) == -1); + CHECK(Q::round_div(-3, 2) == -2); + + // Below half rounds down + CHECK(Q::round_div(7, 3) == 2); + CHECK(Q::round_div(-7, 3) == -2); + + // Above half rounds up + CHECK(Q::round_div(7, 4) == 2); + CHECK(Q::round_div(-7, 4) == -2); + + // Negative denominator normalizes correctly + CHECK(Q::round_div(7, -4) == Q::round_div(-7, 4)); + CHECK(Q::round_div(-7, -4) == Q::round_div(7, 4)); +} + +TEST_CASE("fixed::Q mul performs fixed-point multiplication with rounding", "[dag_builder][fixed_point]") { + using Q = fixed::Q<16>; + + // Small integer multiplication preserves scale + const auto a_val = Q::from_integer(2).value; + const auto b_val = Q::from_integer(3).value; + const auto expected = Q::from_integer(6).value; + CHECK(Q::mul(a_val, b_val) == expected); + + // Half-LSB case: multiplying fixed-point 0.5 (32768) by itself + const auto half = Q::from_real(0.5).value; + const auto quarter = Q::mul(half, half); + CHECK(quarter == Q::from_real(0.25).value); + + // Very small multiplication rounds to zero + CHECK(Q::mul(1, 1) == 0); + + // Inverse relationship with division for values at scale + const auto a = Q::scale; + const auto b = Q::scale * 3; + const auto product = Q::mul(a, b); + const auto divided = Q::div(product, a); + CHECK(divided == b); +} + +TEST_CASE("fixed::Q div performs fixed-point division with rounding", "[dag_builder][fixed_point]") { + using Q = fixed::Q<16>; + + // Division is inverse of multiplication for exactly representable values + const auto a_val = Q::from_integer(6).value; + const auto b_val = Q::from_integer(2).value; + const auto quotient = Q::div(a_val, b_val); + CHECK(quotient == Q::from_integer(3).value); + + // Dividing by same value gives 1 + CHECK(Q::div(Q::scale * 7, Q::scale * 7) == Q::scale); + + // Half division + const auto one = Q::scale; + const auto two = Q::from_integer(2).value; + const auto half = Q::div(one, two); + CHECK(half == Q::from_real(0.5).value); +} + +TEST_CASE("fp::orient detects triangle winding and collinearity", "[dag_builder][fixed_point]") { + using Q = fixed::Q<16>; + using FixedPoint = fp::Vec2; + using FixedScalar = fp::Scalar; + + // CCW triple gives positive orient + const auto a = FixedPoint::from_real(0.0, 0.0); + const auto b = FixedPoint::from_real(1.0, 0.0); + const auto p_ccw = FixedPoint::from_real(0.5, 0.5); + const auto orient_ccw = fp::orient(a, b, p_ccw); + CHECK(orient_ccw > FixedScalar::zero()); + + // CW triple gives negative orient + const auto p_cw = FixedPoint::from_real(0.5, -0.5); + const auto orient_cw = fp::orient(a, b, p_cw); + CHECK(orient_cw < FixedScalar::zero()); + + // Collinear points give zero orient + const auto p_collinear = FixedPoint::from_real(0.5, 0.0); + const auto orient_collinear = fp::orient(a, b, p_collinear); + CHECK(orient_collinear == FixedScalar::zero()); + + // Collinear with integer coordinates + const auto p1 = FixedPoint::from_real(0.0, 0.0); + const auto p2 = FixedPoint::from_real(2.0, 0.0); + const auto p3 = FixedPoint::from_real(1.0, 0.0); + CHECK(fp::orient(p1, p2, p3) == FixedScalar::zero()); + + // Half-coordinate collinearity + const auto q1 = FixedPoint::from_real(0.0, 0.0); + const auto q2 = FixedPoint::from_real(1.0, 0.0); + const auto q3 = FixedPoint::from_real(0.5, 0.0); + CHECK(fp::orient(q1, q2, q3) == FixedScalar::zero()); +} diff --git a/unittests/dag_builder/merge_clusterings.cpp b/unittests/dag_builder/merge_clusterings.cpp new file mode 100644 index 00000000..90abbb9c --- /dev/null +++ b/unittests/dag_builder/merge_clusterings.cpp @@ -0,0 +1,215 @@ +#include +#include +#include +#include + +#include +#include + +#include "../catch2_helpers.h" + +#include "cluster.h" +#include "merge/clusterings.h" + +namespace { + +cv::Mat make_texture(const uchar tag) { + cv::Mat texture(1, 1, CV_8UC3); + texture.at(0, 0) = cv::Vec3b(tag, tag, tag); + return texture; +} + +// A single, unconnected triangle. Since none of its edges have a twin, all +// three vertices are boundary vertices and are eligible for welding. +Clustering make_triangle_clustering(const glm::dvec3 &v0, const glm::dvec3 &v1, const glm::dvec3 &v2, const uchar texture_tag = 0) { + Clustering clustering; + clustering.positions = {v0, v1, v2}; + clustering.textures.add(make_texture(texture_tag)); + + Cluster cluster; + cluster.id = 0; + cluster.vertex_indices = {0, 1, 2}; + cluster.local_triangles = {glm::uvec3(0, 1, 2)}; + cluster.texture_id = 0; + clustering.clusters.push_back(std::move(cluster)); + + return clustering; +} + +// A closed fan of four triangles around a center vertex. Every edge touching +// the center vertex is shared by two triangles (and thus has a twin), so the +// center vertex is interior. The four outer vertices each touch one of the +// unmatched outer boundary edges, so they remain boundary vertices. +Clustering make_fan_clustering( + const glm::dvec3 ¢er, + const glm::dvec3 &outer0, + const glm::dvec3 &outer1, + const glm::dvec3 &outer2, + const glm::dvec3 &outer3, + const uchar texture_tag = 0) { + Clustering clustering; + clustering.positions = {outer0, outer1, outer2, outer3, center}; + clustering.textures.add(make_texture(texture_tag)); + + Cluster cluster; + cluster.id = 0; + cluster.vertex_indices = {0, 1, 2, 3, 4}; + cluster.local_triangles = { + glm::uvec3(4, 0, 1), + glm::uvec3(4, 1, 2), + glm::uvec3(4, 2, 3), + glm::uvec3(4, 3, 0), + }; + cluster.texture_id = 0; + clustering.clusters.push_back(std::move(cluster)); + + return clustering; +} + +bool contains_position_near(const std::vector &positions, const glm::dvec3 &target, const double tolerance) { + for (const glm::dvec3 &position : positions) { + if (glm::distance(position, target) <= tolerance) { + return true; + } + } + return false; +} + +} // namespace + +TEST_CASE("merge_clusterings returns an empty clustering for empty input", "[dag_builder][merge]") { + const std::vector clusterings; + + const Clustering merged = merge_clusterings(clusterings, 0.01); + + CHECK(merged.positions.empty()); + CHECK(merged.clusters.empty()); +} + +TEST_CASE("merge_clusterings returns a copy of the single input clustering unchanged", "[dag_builder][merge]") { + const Clustering clustering = make_triangle_clustering(glm::dvec3(0.0, 0.0, 0.0), glm::dvec3(1.0, 0.0, 0.0), glm::dvec3(0.0, 1.0, 0.0)); + const std::vector clusterings = {clustering}; + + const Clustering merged = merge_clusterings(clusterings, 0.01); + + REQUIRE(merged.positions.size() == clustering.positions.size()); + for (size_t i = 0; i < clustering.positions.size(); i++) { + CHECK(merged.positions[i] == clustering.positions[i]); + } + REQUIRE(merged.clusters.size() == clustering.clusters.size()); + CHECK(merged.clusters[0].vertex_indices == clustering.clusters[0].vertex_indices); + CHECK(merged.clusters[0].local_triangles == clustering.clusters[0].local_triangles); +} + +TEST_CASE("merge_clusterings throws std::invalid_argument for non-positive or non-finite epsilon", "[dag_builder][merge]") { + const std::vector clusterings = { + make_triangle_clustering(glm::dvec3(0.0, 0.0, 0.0), glm::dvec3(10.0, 0.0, 0.0), glm::dvec3(5.0, 10.0, 0.0)), + make_triangle_clustering(glm::dvec3(0.001, 0.0, 0.0), glm::dvec3(110.0, 0.0, 0.0), glm::dvec3(105.0, 10.0, 0.0)), + }; + + CHECK_THROWS_AS(merge_clusterings(clusterings, 0.0), std::invalid_argument); + CHECK_THROWS_AS(merge_clusterings(clusterings, -1.0), std::invalid_argument); + CHECK_THROWS_AS(merge_clusterings(clusterings, std::numeric_limits::quiet_NaN()), std::invalid_argument); + CHECK_THROWS_AS(merge_clusterings(clusterings, std::numeric_limits::infinity()), std::invalid_argument); +} + +TEST_CASE("merge_clusterings welds boundary vertices within epsilon and preserves triangles", "[dag_builder][merge]") { + const MergeMode merge_mode = GENERATE(MergeMode::GreedyLocal, MergeMode::ConnectedComponents, MergeMode::MultipartiteNearest); + + const Clustering a = make_triangle_clustering(glm::dvec3(0.0, 0.0, 0.0), glm::dvec3(10.0, 0.0, 0.0), glm::dvec3(5.0, 10.0, 0.0)); + const Clustering b = make_triangle_clustering(glm::dvec3(0.001, 0.0, 0.0), glm::dvec3(110.0, 0.0, 0.0), glm::dvec3(105.0, 10.0, 0.0)); + const std::vector clusterings = {a, b}; + + const Clustering merged = merge_clusterings(clusterings, 0.01, merge_mode); + + REQUIRE(merged.clusters.size() == 2); + CHECK(merged.clusters[0].local_triangles.size() == 1); + CHECK(merged.clusters[1].local_triangles.size() == 1); + REQUIRE(merged.vertex_count() == 5); + CHECK(contains_position_near(merged.positions, glm::dvec3(0.0005, 0.0, 0.0), 1e-9)); + CHECK(contains_position_near(merged.positions, glm::dvec3(10.0, 0.0, 0.0), 1e-9)); + CHECK(contains_position_near(merged.positions, glm::dvec3(110.0, 0.0, 0.0), 1e-9)); +} + +TEST_CASE("merge_clusterings does not weld vertices further apart than epsilon", "[dag_builder][merge]") { + const MergeMode merge_mode = GENERATE(MergeMode::GreedyLocal, MergeMode::ConnectedComponents, MergeMode::MultipartiteNearest); + + const Clustering a = make_triangle_clustering(glm::dvec3(0.0, 0.0, 0.0), glm::dvec3(10.0, 0.0, 0.0), glm::dvec3(5.0, 10.0, 0.0)); + const Clustering b = make_triangle_clustering(glm::dvec3(1.0, 0.0, 0.0), glm::dvec3(110.0, 0.0, 0.0), glm::dvec3(105.0, 10.0, 0.0)); + const std::vector clusterings = {a, b}; + + const Clustering merged = merge_clusterings(clusterings, 0.01, merge_mode); + + CHECK(merged.vertex_count() == 6); +} + +TEST_CASE("merge_clusterings welds mutually close boundary vertices across three clusterings", "[dag_builder][merge]") { + const MergeMode merge_mode = GENERATE(MergeMode::GreedyLocal, MergeMode::ConnectedComponents, MergeMode::MultipartiteNearest); + + const Clustering a = make_triangle_clustering(glm::dvec3(0.0, 0.0, 0.0), glm::dvec3(10.0, 0.0, 0.0), glm::dvec3(5.0, 10.0, 0.0)); + const Clustering b = make_triangle_clustering(glm::dvec3(0.001, 0.0, 0.0), glm::dvec3(110.0, 0.0, 0.0), glm::dvec3(105.0, 10.0, 0.0)); + const Clustering c = make_triangle_clustering(glm::dvec3(-0.001, 0.0, 0.0), glm::dvec3(210.0, 0.0, 0.0), glm::dvec3(205.0, 10.0, 0.0)); + const std::vector clusterings = {a, b, c}; + + const Clustering merged = merge_clusterings(clusterings, 0.01, merge_mode); + + REQUIRE(merged.clusters.size() == 3); + CHECK(merged.clusters[0].local_triangles.size() == 1); + CHECK(merged.clusters[1].local_triangles.size() == 1); + CHECK(merged.clusters[2].local_triangles.size() == 1); + REQUIRE(merged.vertex_count() == 7); +} + +TEST_CASE("merge_clusterings averages the position of a vertex welded across three clusterings", "[dag_builder][merge]") { + const MergeMode merge_mode = GENERATE(MergeMode::GreedyLocal, MergeMode::ConnectedComponents, MergeMode::MultipartiteNearest); + + const Clustering a = make_triangle_clustering(glm::dvec3(0.0, 0.0, 0.0), glm::dvec3(10.0, 0.0, 0.0), glm::dvec3(5.0, 10.0, 0.0)); + const Clustering b = make_triangle_clustering(glm::dvec3(0.001, 0.0, 0.0), glm::dvec3(110.0, 0.0, 0.0), glm::dvec3(105.0, 10.0, 0.0)); + const Clustering c = make_triangle_clustering(glm::dvec3(-0.001, 0.0, 0.0), glm::dvec3(210.0, 0.0, 0.0), glm::dvec3(205.0, 10.0, 0.0)); + const std::vector clusterings = {a, b, c}; + + const Clustering merged = merge_clusterings(clusterings, 0.01, merge_mode, true); + + REQUIRE(merged.vertex_count() == 7); + CHECK(contains_position_near(merged.positions, glm::dvec3(0.0, 0.0, 0.0), 1e-9)); +} + +TEST_CASE("merge_clusterings never welds interior vertices even when coincident", "[dag_builder][merge]") { + const MergeMode merge_mode = GENERATE(MergeMode::GreedyLocal, MergeMode::ConnectedComponents, MergeMode::MultipartiteNearest); + + const glm::dvec3 center(0.0, 0.0, 0.0); + const Clustering a = make_fan_clustering(center, glm::dvec3(10.0, 0.0, 0.0), glm::dvec3(0.0, 10.0, 0.0), glm::dvec3(-10.0, 0.0, 0.0), glm::dvec3(0.0, -10.0, 0.0)); + const Clustering b = make_fan_clustering(center, glm::dvec3(510.0, 0.0, 0.0), glm::dvec3(500.0, 10.0, 0.0), glm::dvec3(490.0, 0.0, 0.0), glm::dvec3(500.0, -10.0, 0.0)); + const std::vector clusterings = {a, b}; + + const Clustering merged = merge_clusterings(clusterings, 0.01, merge_mode); + + CHECK(merged.vertex_count() == 10); +} + +TEST_CASE("merge_clusterings averages positions of welded vertices when average_positions is true", "[dag_builder][merge]") { + const MergeMode merge_mode = GENERATE(MergeMode::GreedyLocal, MergeMode::ConnectedComponents, MergeMode::MultipartiteNearest); + + const Clustering a = make_triangle_clustering(glm::dvec3(0.0, 0.0, 0.0), glm::dvec3(10.0, 0.0, 0.0), glm::dvec3(5.0, 10.0, 0.0)); + const Clustering b = make_triangle_clustering(glm::dvec3(0.0, 0.0, 0.5), glm::dvec3(110.0, 0.0, 0.0), glm::dvec3(105.0, 10.0, 0.0)); + const std::vector clusterings = {a, b}; + + const Clustering merged = merge_clusterings(clusterings, 1.0, merge_mode, true); + + REQUIRE(merged.vertex_count() == 5); + CHECK(contains_position_near(merged.positions, glm::dvec3(0.0, 0.0, 0.25), 1e-9)); +} + +TEST_CASE("merge_clusterings keeps the last contributing position when average_positions is false", "[dag_builder][merge]") { + const MergeMode merge_mode = GENERATE(MergeMode::GreedyLocal, MergeMode::ConnectedComponents, MergeMode::MultipartiteNearest); + + const Clustering a = make_triangle_clustering(glm::dvec3(0.0, 0.0, 0.0), glm::dvec3(10.0, 0.0, 0.0), glm::dvec3(5.0, 10.0, 0.0)); + const Clustering b = make_triangle_clustering(glm::dvec3(0.0, 0.0, 0.5), glm::dvec3(110.0, 0.0, 0.0), glm::dvec3(105.0, 10.0, 0.0)); + const std::vector clusterings = {a, b}; + + const Clustering merged = merge_clusterings(clusterings, 1.0, merge_mode, false); + + REQUIRE(merged.vertex_count() == 5); + CHECK(contains_position_near(merged.positions, glm::dvec3(0.0, 0.0, 0.5), 1e-9)); + CHECK_FALSE(contains_position_near(merged.positions, glm::dvec3(0.0, 0.0, 0.0), 1e-9)); +} diff --git a/unittests/dag_builder/pull_reproject_texture.cpp b/unittests/dag_builder/pull_reproject_texture.cpp new file mode 100644 index 00000000..6cd920ee --- /dev/null +++ b/unittests/dag_builder/pull_reproject_texture.cpp @@ -0,0 +1,233 @@ +#include +#include + +#include "../catch2_helpers.h" + +#include "atlas/pull_reproject_texture.h" + +namespace { + +cv::Mat make_texture(const cv::Vec3b top_left, const cv::Vec3b top_right) { + cv::Mat texture(1, 2, CV_8UC3); + texture.at(0, 0) = top_left; + texture.at(0, 1) = top_right; + return texture; +} + +} // namespace + +TEST_CASE("TextureReprojector copies pixels under an identity mapping", "[dag_builder][pull_reproject_texture]") { + const cv::Mat source = make_texture(cv::Vec3b(10, 20, 30), cv::Vec3b(40, 50, 60)); + + std::array triangles; + triangles[0].source_image_index = 0; + triangles[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles[0].target_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles[1].source_image_index = 0; + triangles[1].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 1), glm::dvec2(0, 1)}; + triangles[1].target_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 1), glm::dvec2(0, 1)}; + + TextureReprojector reprojector(glm::uvec2(2, 1), CV_8UC3, ReprojectionOptions{1, cv::INTER_NEAREST}); + const cv::Mat output = reprojector.render(std::span(&source, 1), triangles); + + CHECK(output.at(0, 0) == cv::Vec3b(10, 20, 30)); + CHECK(output.at(0, 1) == cv::Vec3b(40, 50, 60)); +} + +TEST_CASE("TextureReprojector linearly blends between neighboring source pixels", "[dag_builder][pull_reproject_texture]") { + const cv::Mat source = make_texture(cv::Vec3b(0, 0, 0), cv::Vec3b(200, 200, 200)); + + std::array triangles; + triangles[0].source_image_index = 0; + triangles[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles[0].target_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + + TextureReprojector reprojector(glm::uvec2(1, 1), CV_8UC3, ReprojectionOptions{1, cv::INTER_LINEAR}); + const cv::Mat output = reprojector.render(std::span(&source, 1), triangles); + + // The single output sample lands halfway between the two source pixels. + CHECK(output.at(0, 0) == cv::Vec3b(100, 100, 100)); +} + +TEST_CASE("TextureReprojector leaves the background where triangles are degenerate", "[dag_builder][pull_reproject_texture]") { + const cv::Mat source = make_texture(cv::Vec3b(10, 20, 30), cv::Vec3b(40, 50, 60)); + + std::array triangles; + triangles[0].source_image_index = 0; + triangles[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles[0].target_uvs = {glm::dvec2(0.5, 0.5), glm::dvec2(0.5, 0.5), glm::dvec2(0.5, 0.5)}; + + TextureReprojector reprojector(glm::uvec2(1, 1), CV_8UC3, ReprojectionOptions{1, cv::INTER_NEAREST}); + const cv::Mat output = reprojector.render(std::span(&source, 1), triangles); + + CHECK(output.at(0, 0) == cv::Vec3b(0, 0, 0)); +} + +TEST_CASE("TextureReprojector supersamples and averages partial pixel coverage", "[dag_builder][pull_reproject_texture]") { + // Triangle covering only the lower-left sample at UV (0.25, 0.25) + cv::Mat source(1, 1, CV_8UC3, cv::Vec3b(100, 100, 100)); + + std::array triangles; + triangles[0].source_image_index = 0; + triangles[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + // Triangle covers lower-left quarter of output pixel + triangles[0].target_uvs = {glm::dvec2(0, 0), glm::dvec2(0.5, 0), glm::dvec2(0.5, 0.5)}; + + TextureReprojector reprojector(glm::uvec2(1, 1), CV_8UC3, ReprojectionOptions{2, cv::INTER_NEAREST}); + const cv::Mat output = reprojector.render(std::span(&source, 1), triangles); + + // Only 1 of 4 samples is covered, output is average of covered samples + CHECK(output.at(0, 0) == cv::Vec3b(100, 100, 100)); +} + +TEST_CASE("TextureReprojector handles non-uniform coverage with gradient source", "[dag_builder][pull_reproject_texture]") { + // 2x1 gradient source for testing coverage-aware averaging + cv::Mat source(1, 2, CV_8UC3); + source.at(0, 0) = cv::Vec3b(50, 50, 50); + source.at(0, 1) = cv::Vec3b(150, 150, 150); + + std::array triangles; + triangles[0].source_image_index = 0; + triangles[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles[0].target_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + + TextureReprojector reprojector(glm::uvec2(1, 1), CV_8UC3, ReprojectionOptions{2, cv::INTER_LINEAR}); + const cv::Mat output = reprojector.render(std::span(&source, 1), triangles); + + const auto pixel = output.at(0, 0); + // With linear interpolation across gradient, should be intermediate value + CHECK(pixel[0] > 50); + CHECK(pixel[0] < 150); +} + +TEST_CASE("TextureReprojector background stays zero when no samples covered", "[dag_builder][pull_reproject_texture]") { + cv::Mat source(1, 1, CV_8UC3, cv::Vec3b(255, 255, 255)); + + std::array triangles; + triangles[0].source_image_index = 0; + triangles[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + // Triangle outside output area + triangles[0].target_uvs = {glm::dvec2(2, 2), glm::dvec2(3, 2), glm::dvec2(3, 3)}; + + TextureReprojector reprojector(glm::uvec2(1, 1), CV_8UC3, ReprojectionOptions{2, cv::INTER_NEAREST}); + const cv::Mat output = reprojector.render(std::span(&source, 1), triangles); + + CHECK(output.at(0, 0) == cv::Vec3b(0, 0, 0)); +} + +TEST_CASE("TextureReprojector quad seam has no holes with supersampling", "[dag_builder][pull_reproject_texture]") { + // Full coverage with two triangles forming a quad and supersampling + cv::Mat source(1, 1, CV_8UC3, cv::Vec3b(128, 128, 128)); + + std::array triangles; + triangles[0].source_image_index = 0; + triangles[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles[0].target_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles[1].source_image_index = 0; + triangles[1].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 1), glm::dvec2(0, 1)}; + triangles[1].target_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 1), glm::dvec2(0, 1)}; + + TextureReprojector reprojector(glm::uvec2(4, 4), CV_8UC3, ReprojectionOptions{2, cv::INTER_NEAREST}); + const cv::Mat output = reprojector.render(std::span(&source, 1), triangles); + + // Every output pixel must equal source color (no holes on seams) + for (int y = 0; y < output.rows; y++) { + for (int x = 0; x < output.cols; x++) { + REQUIRE(output.at(y, x) == cv::Vec3b(128, 128, 128)); + } + } +} + +TEST_CASE("TextureReprojector renders identically regardless of winding order", "[dag_builder][pull_reproject_texture]") { + cv::Mat source(1, 1, CV_8UC3, cv::Vec3b(64, 64, 64)); + + // CCW triangle + std::array triangles_ccw; + triangles_ccw[0].source_image_index = 0; + triangles_ccw[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles_ccw[0].target_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + + // CW triangle (swap two vertices) + std::array triangles_cw; + triangles_cw[0].source_image_index = 0; + triangles_cw[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 1), glm::dvec2(1, 0)}; + triangles_cw[0].target_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 1), glm::dvec2(1, 0)}; + + TextureReprojector reprojector(glm::uvec2(2, 2), CV_8UC3, ReprojectionOptions{1, cv::INTER_NEAREST}); + + const cv::Mat output_ccw = reprojector.render(std::span(&source, 1), triangles_ccw); + const cv::Mat output_cw = reprojector.render(std::span(&source, 1), triangles_cw); + + // Winding repair should make both renders identical + for (int y = 0; y < output_ccw.rows; y++) { + for (int x = 0; x < output_ccw.cols; x++) { + CHECK(output_ccw.at(y, x) == output_cw.at(y, x)); + } + } +} + +TEST_CASE("TextureReprojector handles multiple source images with different colors", "[dag_builder][pull_reproject_texture]") { + // Two 1x1 sources with different colors + cv::Mat source0(1, 1, CV_8UC3, cv::Vec3b(255, 0, 0)); // red + cv::Mat source1(1, 1, CV_8UC3, cv::Vec3b(0, 255, 0)); // green + std::array sources = {source0, source1}; + + // Two triangles, each using a different source + std::array triangles; + // Left triangle uses source 0 (red) + triangles[0].source_image_index = 0; + triangles[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles[0].target_uvs = {glm::dvec2(0, 0), glm::dvec2(0.5, 0), glm::dvec2(0.5, 1)}; + // Right triangle uses source 1 (green) + triangles[1].source_image_index = 1; + triangles[1].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles[1].target_uvs = {glm::dvec2(0.5, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + + TextureReprojector reprojector(glm::uvec2(2, 1), CV_8UC3, ReprojectionOptions{1, cv::INTER_NEAREST}); + const cv::Mat output = reprojector.render(sources, triangles); + + CHECK(output.at(0, 0) == cv::Vec3b(255, 0, 0)); // red + CHECK(output.at(0, 1) == cv::Vec3b(0, 255, 0)); // green +} + +TEST_CASE("TextureReprojector clamps UVs and applies edge replication", "[dag_builder][pull_reproject_texture]") { + cv::Mat source(2, 2, CV_8UC3); + source.at(0, 0) = cv::Vec3b(255, 0, 0); // red + source.at(0, 1) = cv::Vec3b(0, 255, 0); // green + source.at(1, 0) = cv::Vec3b(0, 0, 255); // blue + source.at(1, 1) = cv::Vec3b(255, 255, 0); // yellow + + std::array triangles; + triangles[0].source_image_index = 0; + // All source UVs lie outside [0, 1] and get clamped to the (1, 1) corner. + triangles[0].source_uvs = {glm::dvec2(2, 2), glm::dvec2(3, 2), glm::dvec2(3, 3)}; + triangles[0].target_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + + TextureReprojector reprojector(glm::uvec2(1, 1), CV_8UC3, ReprojectionOptions{1, cv::INTER_NEAREST}); + const cv::Mat output = reprojector.render(std::span(&source, 1), triangles); + + // Every sample reads the bottom-right border texel (yellow). + CHECK(output.at(0, 0) == cv::Vec3b(255, 255, 0)); +} + +TEST_CASE("TextureReprojector pins rounding behavior for linear interpolation blends", "[dag_builder][pull_reproject_texture]") { + // Create a 2x1 source with black and white to test 127.5 blend + cv::Mat source(1, 2, CV_8UC3); + source.at(0, 0) = cv::Vec3b(0, 0, 0); // black + source.at(0, 1) = cv::Vec3b(255, 255, 255); // white + + std::array triangles; + triangles[0].source_image_index = 0; + triangles[0].source_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + triangles[0].target_uvs = {glm::dvec2(0, 0), glm::dvec2(1, 0), glm::dvec2(1, 1)}; + + TextureReprojector reprojector(glm::uvec2(1, 1), CV_8UC3, ReprojectionOptions{1, cv::INTER_LINEAR}); + const cv::Mat output = reprojector.render(std::span(&source, 1), triangles); + + const auto pixel = output.at(0, 0); + // Linear interpolation at center between 0 and 255 produces 127.5 + // OpenCV cvRound uses round-half-to-even: 127.5 rounds to 128 + CHECK(pixel[0] == 128); + CHECK(pixel[1] == 128); + CHECK(pixel[2] == 128); +} diff --git a/unittests/terrainlib/octree_odd_level_shifted.cpp b/unittests/terrainlib/octree_odd_level_shifted.cpp index 3559a59c..72e7fea0 100644 --- a/unittests/terrainlib/octree_odd_level_shifted.cpp +++ b/unittests/terrainlib/octree_odd_level_shifted.cpp @@ -380,15 +380,15 @@ TEST_CASE("OddLevelShifted::get_node_bounds all level-1 nodes together tile the CHECK(max_node.max.z == Catch::Approx(space.bounds().max.z)); } -TEST_CASE("OddLevelShifted::get_node_bounds_with_children non-leaf equals union of children bounds", "[octree::OddLevelShifted]") { +TEST_CASE("OddLevelShifted::get_node_bounds_with_children non-leaf equals union of node and children bounds", "[octree::OddLevelShifted]") { const OddLevelShifted space = OddLevelShifted::earth(); const Id id{2, {1, 1, 1}}; const Bounds actual = space.get_node_bounds_with_children(id); - const Bounds child_first = space.get_node_bounds(id.child(0).value()); - const Bounds child_last = space.get_node_bounds(id.child(7).value()); - const Bounds expected(child_first.min, child_last.max); + Bounds expected = space.get_node_bounds(id); + expected.expand_by(space.get_node_bounds(id.child(0).value())); + expected.expand_by(space.get_node_bounds(id.child(7).value())); CHECK(actual.min.x == Catch::Approx(expected.min.x)); CHECK(actual.min.y == Catch::Approx(expected.min.y)); @@ -398,7 +398,7 @@ TEST_CASE("OddLevelShifted::get_node_bounds_with_children non-leaf equals union CHECK(actual.max.z == Catch::Approx(expected.max.z)); } -TEST_CASE("OddLevelShifted::get_node_bounds_with_children even-level min extends and max shifts negatively", "[octree::OddLevelShifted]") { +TEST_CASE("OddLevelShifted::get_node_bounds_with_children even-level includes node and children", "[octree::OddLevelShifted]") { const OddLevelShifted space = OddLevelShifted::earth(); const Id id{2, {1, 1, 1}}; @@ -408,9 +408,9 @@ TEST_CASE("OddLevelShifted::get_node_bounds_with_children even-level min extends CHECK(with_children.min.x < node.min.x); CHECK(with_children.min.y < node.min.y); CHECK(with_children.min.z < node.min.z); - CHECK(with_children.max.x < node.max.x); - CHECK(with_children.max.y < node.max.y); - CHECK(with_children.max.z < node.max.z); + CHECK(with_children.max.x == Catch::Approx(node.max.x)); + CHECK(with_children.max.y == Catch::Approx(node.max.y)); + CHECK(with_children.max.z == Catch::Approx(node.max.z)); } TEST_CASE("OddLevelShifted::contains interior point returns true", "[octree::OddLevelShifted]") { diff --git a/unittests/terrainlib/offset_table.cpp b/unittests/terrainlib/offset_table.cpp index 1106613b..92aa1cf1 100644 --- a/unittests/terrainlib/offset_table.cpp +++ b/unittests/terrainlib/offset_table.cpp @@ -4,17 +4,15 @@ #include #include -using Table = OffsetTable; - TEST_CASE("OffsetTable default construction") { - Table table; + OffsetTable table; CHECK(table.empty()); CHECK(table.size() == 0); CHECK(table.total_size() == 0); } TEST_CASE("OffsetTable append_length") { - Table table; + OffsetTable table; table.append_length(3); table.append_length(5); table.append_length(2); @@ -22,91 +20,91 @@ TEST_CASE("OffsetTable append_length") { CHECK(table.size() == 3); CHECK(table.total_size() == 10); - auto r0 = table.range(0); + auto r0 = table.segment_range(0); CHECK(r0.begin == 0); CHECK(r0.end == 3); - auto r1 = table.range(1); + auto r1 = table.segment_range(1); CHECK(r1.begin == 3); CHECK(r1.end == 8); - auto r2 = table.range(2); + auto r2 = table.segment_range(2); CHECK(r2.begin == 8); CHECK(r2.end == 10); } TEST_CASE("OffsetTable append_lengths via span") { - Table table; + OffsetTable table; std::array lengths = {3, 5, 2}; table.append_lengths(lengths); CHECK(table.size() == 3); CHECK(table.total_size() == 10); - auto r0 = table.range(0); + auto r0 = table.segment_range(0); CHECK(r0.begin == 0); CHECK(r0.end == 3); - auto r1 = table.range(1); + auto r1 = table.segment_range(1); CHECK(r1.begin == 3); CHECK(r1.end == 8); - auto r2 = table.range(2); + auto r2 = table.segment_range(2); CHECK(r2.begin == 8); CHECK(r2.end == 10); } TEST_CASE("OffsetTable element_size") { - Table table; + OffsetTable table; table.append_length(3); table.append_length(5); table.append_length(2); - CHECK(table.element_size(0) == 3); - CHECK(table.element_size(1) == 5); - CHECK(table.element_size(2) == 2); + CHECK(table.segment_size(0) == 3); + CHECK(table.segment_size(1) == 5); + CHECK(table.segment_size(2) == 2); } TEST_CASE("OffsetTable locate") { - Table table; + OffsetTable table; table.append_length(3); table.append_length(5); table.append_length(2); { auto result = table.locate(0); - CHECK(result.element == 0); + CHECK(result.segment == 0); CHECK(result.range.begin == 0); CHECK(result.range.end == 3); } { auto result = table.locate(3); - CHECK(result.element == 1); + CHECK(result.segment == 1); CHECK(result.range.begin == 3); CHECK(result.range.end == 8); } { auto result = table.locate(7); - CHECK(result.element == 1); + CHECK(result.segment == 1); CHECK(result.range.begin == 3); CHECK(result.range.end == 8); } { auto result = table.locate(8); - CHECK(result.element == 2); + CHECK(result.segment == 2); CHECK(result.range.begin == 8); CHECK(result.range.end == 10); } { auto result = table.locate(9); - CHECK(result.element == 2); + CHECK(result.segment == 2); CHECK(result.range.begin == 8); CHECK(result.range.end == 10); } } TEST_CASE("OffsetTable locate out of range throws") { - Table table; + OffsetTable table; table.append_length(3); table.append_length(5); table.append_length(2); @@ -116,16 +114,16 @@ TEST_CASE("OffsetTable locate out of range throws") { } TEST_CASE("OffsetTable range out of range throws") { - Table table; + OffsetTable table; table.append_length(3); table.append_length(5); - CHECK_THROWS_AS(table.range(2), std::out_of_range); - CHECK_THROWS_AS(table.range(100), std::out_of_range); + CHECK_THROWS_AS(table.segment_range(2), std::out_of_range); + CHECK_THROWS_AS(table.segment_range(100), std::out_of_range); } TEST_CASE("OffsetTable set_end on last element to extend") { - Table table; + OffsetTable table; table.append_length(3); table.append_length(5); table.append_length(2); @@ -133,30 +131,30 @@ TEST_CASE("OffsetTable set_end on last element to extend") { table.set_end(2, 15); CHECK(table.total_size() == 15); - auto r2 = table.range(2); + auto r2 = table.segment_range(2); CHECK(r2.begin == 8); CHECK(r2.end == 15); } TEST_CASE("OffsetTable set_begin on non-first element") { - Table table; + OffsetTable table; table.append_length(3); table.append_length(5); table.append_length(2); table.set_begin(1, 4); - auto r0 = table.range(0); + auto r0 = table.segment_range(0); CHECK(r0.begin == 0); CHECK(r0.end == 4); - auto r1 = table.range(1); + auto r1 = table.segment_range(1); CHECK(r1.begin == 4); CHECK(r1.end == 8); } TEST_CASE("OffsetTable set_begin on element 0 throws") { - Table table; + OffsetTable table; table.append_length(3); table.append_length(5); @@ -164,7 +162,7 @@ TEST_CASE("OffsetTable set_begin on element 0 throws") { } TEST_CASE("OffsetTable clear then re-add elements") { - Table table; + OffsetTable table; table.append_length(3); table.append_length(5); @@ -183,11 +181,11 @@ TEST_CASE("OffsetTable clear then re-add elements") { CHECK(table.size() == 2); CHECK(table.total_size() == 10); - auto r0 = table.range(0); + auto r0 = table.segment_range(0); CHECK(r0.begin == 0); CHECK(r0.end == 4); - auto r1 = table.range(1); + auto r1 = table.segment_range(1); CHECK(r1.begin == 4); CHECK(r1.end == 10); } diff --git a/unittests/terrainlib/range.cpp b/unittests/terrainlib/range.cpp new file mode 100644 index 00000000..61466f04 --- /dev/null +++ b/unittests/terrainlib/range.cpp @@ -0,0 +1,172 @@ +#include "../catch2_helpers.h" +#include "Range.h" + +TEST_CASE("Range: contains value") { + Range r(2, 8); + + CHECK(r.contains(2)); + CHECK(r.contains(5)); + CHECK_FALSE(r.contains(8)); + CHECK_FALSE(r.contains(1)); + CHECK_FALSE(r.contains(9)); +} + +TEST_CASE("Range: is_empty") { + CHECK(Range(3, 3).is_empty()); + CHECK(Range(5, 2).is_empty()); + CHECK_FALSE(Range(2, 5).is_empty()); +} + +TEST_CASE("Range: size") { + CHECK(Range(2, 8).size() == 6); + CHECK(Range(5, 5).size() == 0); +} + +TEST_CASE("Range: single-value constructor") { + Range r(5); + + CHECK(r.start == 5); + CHECK(r.contains(5)); + CHECK_FALSE(r.contains(6)); +} + +TEST_CASE("Range: is_in_bounds") { + CHECK(Range(2, 8).is_in_bounds(8)); + CHECK(Range(2, 8).is_in_bounds(10)); + CHECK_FALSE(Range(2, 8).is_in_bounds(5)); + CHECK_FALSE(Range(5, 2).is_in_bounds(10)); +} + +TEST_CASE("Range: is_overlapping") { + CHECK(Range(2, 8).is_overlapping(Range(5, 12))); + CHECK_FALSE(Range(2, 4).is_overlapping(Range(6, 8))); + CHECK_FALSE(Range(2, 4).is_overlapping(Range(4, 8))); +} + +TEST_CASE("Range: intersect overlapping ranges") { + Range a(2, 8); + Range b(5, 12); + + AnyRange result = a.intersect(b); + + CHECK(result.start == Bound::included(5)); + CHECK(result.end == Bound::excluded(8)); +} + +TEST_CASE("Range: intersect disjoint ranges is empty") { + Range a(2, 4); + Range b(6, 8); + + CHECK(a.intersect(b).is_empty()); +} + +TEST_CASE("Range: intersect with RangeFull returns itself") { + Range a(2, 8); + + AnyRange result = a.intersect(RangeFull{}); + + CHECK(result.start == Bound::included(2)); + CHECK(result.end == Bound::excluded(8)); +} + +TEST_CASE("Range: equality") { + CHECK(Range(2, 8) == Range(2, 8)); + CHECK_FALSE(Range(2, 8) == Range(2, 9)); +} + +TEST_CASE("RangeInclusive: contains value") { + RangeInclusive r(2, 8); + + CHECK(r.contains(2)); + CHECK(r.contains(8)); + CHECK_FALSE(r.contains(1)); + CHECK_FALSE(r.contains(9)); +} + +TEST_CASE("RangeInclusive: is_empty") { + CHECK_FALSE(RangeInclusive(3, 3).is_empty()); + CHECK(RangeInclusive(5, 2).is_empty()); +} + +TEST_CASE("RangeFrom: contains value") { + RangeFrom r{{}, 5}; + + CHECK(r.contains(5)); + CHECK(r.contains(1000)); + CHECK_FALSE(r.contains(4)); +} + +TEST_CASE("RangeFrom: is_empty is always false") { + CHECK_FALSE((RangeFrom{{}, 5}.is_empty())); +} + +TEST_CASE("RangeTo: contains value") { + RangeTo r{{}, 5}; + + CHECK(r.contains(-1000)); + CHECK(r.contains(4)); + CHECK_FALSE(r.contains(5)); +} + +TEST_CASE("RangeToInclusive: contains value") { + RangeToInclusive r{{}, 5}; + + CHECK(r.contains(5)); + CHECK_FALSE(r.contains(6)); +} + +TEST_CASE("RangeFull: contains anything") { + RangeFull r; + + CHECK(r.contains(0)); + CHECK(r.contains(-1000000)); + CHECK(r.contains(1000000)); +} + +TEST_CASE("RangeFull: is_empty is always false") { + CHECK_FALSE(RangeFull{}.is_empty()); +} + +TEST_CASE("AnyRange: constructs from and behaves like each concrete range type") { + AnyRange from_range = Range(2, 8); + CHECK(from_range.contains(5)); + CHECK_FALSE(from_range.contains(8)); + + AnyRange from_inclusive = RangeInclusive(2, 8); + CHECK(from_inclusive.contains(8)); + + AnyRange from_from = RangeFrom{{}, 5}; + CHECK(from_from.contains(1000)); + CHECK_FALSE(from_from.contains(4)); + + AnyRange from_to = RangeTo{{}, 5}; + CHECK(from_to.contains(4)); + CHECK_FALSE(from_to.contains(5)); + + AnyRange from_full = RangeFull{}; + CHECK(from_full.contains(1000000)); + CHECK_FALSE(from_full.is_empty()); +} + +TEST_CASE("AnyRange: is_empty") { + CHECK(AnyRange(Range(3, 3)).is_empty()); + CHECK_FALSE(AnyRange(Range(2, 5)).is_empty()); + CHECK_FALSE(AnyRange(RangeFrom{{}, 5}).is_empty()); +} + +TEST_CASE("AnyRange: to_range resolves unbounded sides using the given length") { + const AnyRange full = RangeFull{}; + CHECK(full.to_range(10) == Range(0, 10)); + + const AnyRange from = RangeFrom{{}, 3}; + CHECK(from.to_range(10) == Range(3, 10)); + + const AnyRange to = RangeTo{{}, 3}; + CHECK(to.to_range(10) == Range(0, 3)); + + const AnyRange to_inclusive = RangeToInclusive{{}, 3}; + CHECK(to_inclusive.to_range(10) == Range(0, 4)); + + const AnyRange bounded = Range(2, 8); + CHECK(bounded.to_range(100) == Range(2, 8)); +} diff --git a/unittests/terrainlib/union_find.cpp b/unittests/terrainlib/union_find.cpp index b5678af1..a1c72732 100644 --- a/unittests/terrainlib/union_find.cpp +++ b/unittests/terrainlib/union_find.cpp @@ -59,7 +59,7 @@ TEST_CASE("UnionFind get_sets", "[UnionFind]") { uf.make_union(0, 2); uf.make_union(1, 3); - auto sets = uf.get_sets(); + auto sets = uf.get_sets_as_map(); REQUIRE(sets.size() == 3); // Each set should contain the right elements