feat(hydrology): add flow breaching for drainage basin (tri and cell-based). Add flow breaching for TerrainTriMesh based on Disjktra shortest path

Author: otto-link

HighMap/include/highmap/geometry/cloud.hpp

@@ -505,6 +505,16 @@ Array cloud_sdf_to_array(const Cloud &cloud,
                          const Array *p_noise_x = nullptr,
                          const Array *p_noise_y = nullptr);
 
+/**
+ * @brief Checks whether the point cloud contains duplicate points.
+ *
+ * @param cloud Input point cloud.
+ * @param eps Distance tolerance for considering two points identical.
+ * @param xy_only If true, compares only the X and Y coordinates.
+ * @return True if duplicate points are found, false otherwise.
+ */
+bool has_duplicates(const Cloud &cloud, float eps = 1e-9f, bool xy_only = true);
+
 /**
  * @brief Interpolate values from an array at the points' `(x, y)` locations.
  *

HighMap/include/highmap/hydrology/drainage_basin.hpp

@@ -64,6 +64,10 @@ class DrainageBasin
       const std::vector<float> &erodibility,
       float                     m_exp) const;
 
+  void flow_breach();
+
+  std::vector<std::vector<glm::vec3>> flow_breach_paths();
+
   float update_elevations(const std::vector<float> &response_times,
                           float                     uplift_rate,
                           const std::vector<float> &max_slope);

HighMap/include/highmap/hydrology/drainage_basin_cell_based.hpp

@@ -33,6 +33,7 @@ class DrainageBasinCellBased
   void compute_receivers(unsigned int seed = 0, float noise_strength = 0.f);
   void compute_receivers_priority_flood();
   void update_stream_tree(unsigned int seed, float noise_strength);
+  void update_stream_tree();
   void update_traversals();
 
   std::vector<glm::ivec2> get_outlets() const;
@@ -59,6 +60,8 @@ class DrainageBasinCellBased
 
   void accumulate_area_by_outlet(Array &acc) const;
 
+  void flow_breach();
+
   // --- Members ---
 
   Array z;

HighMap/include/highmap/terrain_tri_mesh.hpp

@@ -74,6 +74,12 @@ class TerrainTriMesh
     }
   };
 
+  struct ShortestPathResult
+  {
+    std::vector<float>  distance;
+    std::vector<size_t> parent;
+  };
+
 public:
   TerrainTriMesh() = default;
   TerrainTriMesh(const std::vector<glm::vec3> &ref_points);
@@ -112,6 +118,8 @@ class TerrainTriMesh
 
   void subdivise();
 
+  void flow_breach(float epsilon);
+
   // --- Triangle walk and interp ---
 
   bool barycentric(const glm::vec2 &p,
@@ -151,6 +159,24 @@ class TerrainTriMesh
 
   size_t size() const;
 
+  // --- Shortest path ---
+
+  ShortestPathResult compute_shortest_paths_to_hull(
+      bool  use_delta_z = false,
+      float elevation_weight = 1.f) const;
+
+  ShortestPathResult compute_shortest_paths(size_t start,
+                                            bool   use_delta_z = false,
+                                            float elevation_weight = 1.f) const;
+
+  std::vector<size_t> shortest_path(size_t start,
+                                    size_t end,
+                                    bool   use_delta_z = false,
+                                    float  elevation_weight = 1.f) const;
+
+  std::vector<size_t> path_to_hull(size_t                    start,
+                                   const ShortestPathResult &result) const;
+
   // --- Accessors ---
 
   const std::vector<glm::vec3> &get_points() const;
@@ -190,4 +216,9 @@ TerrainTriMesh generate_terrain_tri_mesh_from_heightmap(const Array &z,
                                                         int max_triangles = 0,
                                                         int max_points = 0);
 
+TerrainTriMesh generate_terrain_tri_mesh_from_heightmap_random(
+    const Array  &z,
+    int           control_points_count,
+    std::uint32_t seed);
+
 } // namespace hmap

HighMap/src/geometry/cloud_functions.cpp

@@ -64,6 +64,34 @@ Array cloud_sdf_to_array(const Cloud &cloud,
   return array;
 }
 
+bool has_duplicates(const Cloud &cloud, float eps, bool xy_only)
+{
+  std::vector<glm::vec3> pts = cloud.to_vec3();
+
+  std::sort(pts.begin(),
+            pts.end(),
+            [](const auto &a, const auto &b)
+            { return a.x < b.x || (a.x == b.x && a.y < b.y); });
+
+  if (xy_only)
+  {
+    for (size_t i = 1; i < pts.size(); ++i)
+    {
+      glm::vec2 p0 = {pts[i].x, pts[i].y};
+      glm::vec2 p1 = {pts[i - 1].x, pts[i - 1].y};
+
+      if (glm::distance(p0, p1) < eps) return true;
+    }
+  }
+  else
+  {
+    for (size_t i = 1; i < pts.size(); ++i)
+      if (glm::distance(pts[i], pts[i - 1]) < eps) return true;
+  }
+
+  return false;
+}
+
 std::vector<float> interpolate_values_from_array(const Cloud &cloud,
                                                  const Array &array,
                                                  glm::vec4    bbox)

HighMap/src/hydrology/drainage_basin.cpp

@@ -320,6 +320,82 @@ std::pair<std::vector<size_t>, bool> DrainageBasin::find_subroots()
   return {subroot, has_lake};
 }
 
+void DrainageBasin::flow_breach()
+{
+  auto &pts = this->mesh.get_points();
+
+  for (const auto &[outlet, traversal] : this->traversals)
+  {
+    if (traversal.empty()) continue;
+
+    // traversal is upstream -> downstream
+    size_t p = traversal.front();
+    float  zc = pts[p].z; // current
+
+    // follow receivers downstream to outlet
+    while (true)
+    {
+      const size_t &r = this->receivers[p];
+
+      if (r == p) break; // reached outlet
+
+      if (pts[p].z < zc)
+        zc = pts[p].z;
+      else
+        pts[p].z = zc;
+
+      p = r;
+    }
+  }
+}
+
+std::vector<std::vector<glm::vec3>> DrainageBasin::flow_breach_paths()
+{
+  auto &pts = this->mesh.get_points();
+
+  std::vector<std::vector<glm::vec3>> paths;
+  paths.reserve(this->traversals.size());
+
+  for (const auto &[outlet, traversal] : this->traversals)
+  {
+    if (traversal.empty())
+    {
+      paths.emplace_back();
+      continue;
+    }
+
+    std::vector<glm::vec3> path;
+    path.reserve(traversal.size());
+
+    // traversal is upstream -> downstream
+    size_t p = traversal.front();
+    float  zc = pts[p].z; // current
+
+    path.push_back(pts[p]);
+
+    // follow receivers downstream to outlet
+    while (true)
+    {
+      const size_t &r = this->receivers[p];
+
+      if (r == p) break; // reached outlet
+
+      if (pts[p].z < zc)
+        zc = pts[p].z;
+      else
+        pts[p].z = zc;
+
+      path.push_back(pts[p]);
+
+      p = r;
+    }
+
+    paths.push_back(std::move(path));
+  }
+
+  return paths;
+}
+
 const std::vector<size_t> &DrainageBasin::for_each_upstream(size_t outlet) const
 {
   return traversals.at(outlet);
@@ -341,7 +417,7 @@ std::vector<std::vector<size_t>> DrainageBasin::get_main_channels() const
     std::vector<size_t> channel;
     channel.reserve(traversal.size());
 
-    // traversal is upstream -> downstream
+    // traversal is downstream -> upstream
     size_t p = traversal.front();
     channel.push_back(p);
 

HighMap/src/hydrology/drainage_basin_cell_based.cpp

@@ -293,6 +293,33 @@ std::pair<Mat<glm::ivec2>, bool> DrainageBasinCellBased::find_subroots()
   return {subroot, has_lake};
 }
 
+void DrainageBasinCellBased::flow_breach()
+{
+  for (const auto &[outlet, traversal] : traversals)
+  {
+    if (traversal.empty()) continue;
+
+    // traversal is upstream -> downstream
+    glm::ivec2 p = traversal.front();
+    float      zc = z(p);
+
+    // follow receivers downstream to outlet
+    while (true)
+    {
+      const glm::ivec2 &r = receivers(p.x, p.y);
+
+      if (r == p) break; // reached outlet
+
+      if (z(p) < zc)
+        zc = z(p);
+      else
+        z(p) = zc;
+
+      p = r;
+    }
+  }
+}
+
 std::vector<std::vector<glm::ivec2>> DrainageBasinCellBased::get_main_channels()
     const
 {
@@ -517,6 +544,11 @@ void DrainageBasinCellBased::update_stream_tree(unsigned int seed,
   this->update_traversals();
 }
 
+void DrainageBasinCellBased::update_stream_tree()
+{
+  this->update_stream_tree(0, 0.f);
+}
+
 void DrainageBasinCellBased::update_traversals()
 {
   this->traversals.clear();

HighMap/src/terrain_tri_mesh/flow_breach.cpp

@@ -0,0 +1,33 @@
+/* Copyright (c) 2026 Otto Link. Distributed under the terms of the GNU General
+ * Public License. The full license is in the file LICENSE, distributed with
+ * this software. */
+#include <limits>
+#include <queue>
+
+#include "highmap/terrain_tri_mesh.hpp"
+
+namespace hmap
+{
+
+void TerrainTriMesh::flow_breach(float epsilon)
+{
+  auto paths = compute_shortest_paths_to_hull(true);
+
+  std::vector<bool>   visited(points.size(), false);
+  std::vector<size_t> path;
+
+  for (size_t start = 0; start < points.size(); ++start)
+  {
+    std::vector<size_t> path = this->path_to_hull(start, paths);
+
+    for (size_t i = 0; i < path.size() - 1; ++i)
+    {
+      size_t child = path[i + 1];
+      size_t parent = path[i];
+      float  required = points[parent].z - epsilon;
+      if (points[child].z > required) points[child].z = required;
+    }
+  }
+}
+
+} // namespace hmap