Add accumulation_curvature #166

HighMap/include/highmap/features.hpp

@@ -17,6 +17,31 @@
 namespace hmap
 {
 
+/**
+ * @brief This specific type of curvature reflects how the shape of the
+ * heightmap influences the accumulation of water. Positive accumulation
+ * curvature indicates converging flow, where water tends to accumulate and
+ * concentrate, often leading to the formation of channels or gullies. Negative
+ * accumulation curvature suggests diverging flow, where water is dispersed over
+ * a broader area, which is typical of ridges or hilltops.
+ * @param z Input array.
+ * @param ir Pre-filter radius.
+ * @return Output array.
+ *
+ *
+ * **Example**
+ * @include ex_curvature_gaussian.cpp
+ *
+ * **Result**
+ * @image html ex_curvature_gaussian0.png
+ * @image html ex_curvature_gaussian1.png
+ * @image html ex_curvature_gaussian2.png
+ * @image html ex_curvature_gaussian3.png
+ * @image html ex_curvature_gaussian4.png
+ * @image html ex_curvature_gaussian5.png
+ */
+Array accumulation_curvature(const Array &z, int ir);
+
 /**
  * @brief Return the connected-component labelling of the array.
  *
@@ -54,6 +79,7 @@ Array connected_components(const Array &array,
  * @image html ex_curvature_gaussian2.png
  * @image html ex_curvature_gaussian3.png
  * @image html ex_curvature_gaussian4.png
+ * @image html ex_curvature_gaussian5.png
  */
 Array curvature_gaussian(const Array &z);
 
@@ -72,6 +98,7 @@ Array curvature_gaussian(const Array &z);
  * @image html ex_curvature_gaussian2.png
  * @image html ex_curvature_gaussian3.png
  * @image html ex_curvature_gaussian4.png
+ * @image html ex_curvature_gaussian5.png
  */
 Array curvature_mean(const Array &z);
 
@@ -177,6 +204,7 @@ Array rugosity(const Array &z, int ir);
  * @image html ex_curvature_gaussian2.png
  * @image html ex_curvature_gaussian3.png
  * @image html ex_curvature_gaussian4.png
+ * @image html ex_curvature_gaussian5.png
  */
 Array shape_index(const Array &z, int ir);
 
@@ -197,6 +225,7 @@ Array shape_index(const Array &z, int ir);
  * @image html ex_curvature_gaussian2.png
  * @image html ex_curvature_gaussian3.png
  * @image html ex_curvature_gaussian4.png
+ * @image html ex_curvature_gaussian5.png
  */
 Array unsphericity(const Array &z, int ir);
 

HighMap/src/features/curvature.cpp

@@ -33,6 +33,36 @@ Array compute_k(Array &zx, Array &zy, Array &zxx, Array &zxy, Array &zyy)
 
 // functions
 
+Array accumulation_curvature(const Array &z, int ir)
+{
+  Array ac = z;
+  if (ir > 0)
+    smooth_cpulse(ac, ir);
+
+  Array zx = gradient_x(ac);
+  Array zy = gradient_y(ac);
+  Array zxx = gradient_x(zx);
+  Array zxy = gradient_y(zx);
+  Array zyy = gradient_y(zy);
+
+  Array zx2 = zx * zx;
+  Array zy2 = zy * zy;
+  Array p = zx2 + zy2;
+  Array n = zy2 * zxx - 2.f * zxy * zx * zy + zx2 * zyy;
+
+  Array horizontal_curvature = p * pow(p + 1.f, 0.5f);
+  Array vertical_curvature = p * pow(p + 1.f, 1.5f);
+
+  ac = n * n / (horizontal_curvature * vertical_curvature + 1e-30);
+
+  if (ir > 0)
+    extrapolate_borders(ac, ir + 1, 0.1f);
+  else
+    extrapolate_borders(ac);
+
+  return ac;
+}
+
 Array curvature_gaussian(const Array &z)
 {
   Array k = Array(z.shape); // output

examples/ex_curvature_gaussian/ex_curvature_gaussian.cpp

@@ -19,10 +19,13 @@ int main(void)
 
   hmap::Array un = hmap::unsphericity(z, ir);
 
+  hmap::Array ac = hmap::accumulation_curvature(z, ir);
+
   z.to_png("ex_curvature_gaussian0.png", hmap::cmap::terrain);
   k.to_png("ex_curvature_gaussian1.png", hmap::cmap::jet);
   h.to_png("ex_curvature_gaussian2.png", hmap::cmap::jet);
 
   si.to_png("ex_curvature_gaussian3.png", hmap::cmap::jet);
   un.to_png("ex_curvature_gaussian4.png", hmap::cmap::jet);
+  ac.to_png("ex_curvature_gaussian5.png", hmap::cmap::jet);
 }