Add unsphericity curvature feature

HighMap/include/highmap/features.hpp

@@ -53,6 +53,7 @@ Array connected_components(const Array &array,
  * @image html ex_curvature_gaussian1.png
  * @image html ex_curvature_gaussian2.png
  * @image html ex_curvature_gaussian3.png
+ * @image html ex_curvature_gaussian4.png
  */
 Array curvature_gaussian(const Array &z);
 
@@ -70,6 +71,7 @@ Array curvature_gaussian(const Array &z);
  * @image html ex_curvature_gaussian1.png
  * @image html ex_curvature_gaussian2.png
  * @image html ex_curvature_gaussian3.png
+ * @image html ex_curvature_gaussian4.png
  */
 Array curvature_mean(const Array &z);
 
@@ -174,9 +176,30 @@ Array rugosity(const Array &z, int ir);
  * @image html ex_curvature_gaussian1.png
  * @image html ex_curvature_gaussian2.png
  * @image html ex_curvature_gaussian3.png
+ * @image html ex_curvature_gaussian4.png
  */
 Array shape_index(const Array &z, int ir);
 
+/**
+ * @brief Unsesphericity refers to the degree to which an object or surface
+ * deviates from being perfectly spherical or symmetrical.
+ * @param z Input array.
+ * @param ir Pre-filter radius.
+ * @return Resulting array.
+ * @return Resulting array (> 0.5 for convex and < 0.5 for concave).
+ *
+ * **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
+ */
+Array unsphericity(const Array &z, int ir);
+
 /**
  * @brief Return the "valley width", corresponding to the distance to the
  * concave region frontier (in this concave frontier).

HighMap/src/features/curvature.cpp

@@ -72,13 +72,32 @@ Array shape_index(const Array &z, int ir)
   Array h = compute_h(zx, zy, zxx, zxy, zyy); // mean
 
   Array d = pow(h * h - k, 0.5f);
-  si = 2.f / M_PI * atan(h / (d + 1e-30));
+  si = -2.f / M_PI * atan(h / (d + 1e-30));
   si *= 0.5f;
   si += 0.5f;
 
   return si;
 }
 
+Array unsphericity(const Array &z, int ir)
+{
+  Array si = z;
+  if (ir > 0)
+    smooth_cpulse(si, ir);
+
+  // compute curvature criteria
+  Array zx = gradient_x(si);
+  Array zy = gradient_y(si);
+  Array zxx = gradient_x(zx);
+  Array zxy = gradient_y(zx);
+  Array zyy = gradient_y(zy);
+
+  Array k = compute_k(zx, zy, zxx, zxy, zyy); // gaussian
+  Array h = compute_h(zx, zy, zxx, zxy, zyy); // mean
+
+  return pow(h * h - k, 0.5f);
+}
+
 Array valley_width(const Array &z, int ir)
 {
   Array vw = z;

examples/ex_curvature_gaussian/ex_curvature_gaussian.cpp

@@ -2,18 +2,27 @@
 
 int main(void)
 {
-  hmap::Vec2<int> shape = {256, 256};
+  hmap::Vec2<int>   shape = {256, 256};
+  hmap::Vec2<float> kw = {4.f, 4.f};
+  int               seed = 1;
+
+  hmap::Array z = hmap::noise_fbm(hmap::NoiseType::n_perlin_half,
+                                  shape,
+                                  kw,
+                                  seed);
 
-  hmap::Array z = hmap::gabor(shape, 4.f, 30.f);
   hmap::Array k = hmap::curvature_gaussian(z);
   hmap::Array h = hmap::curvature_mean(z);
 
-  int         ir = 4;
+  int         ir = 16;
   hmap::Array si = hmap::shape_index(z, ir);
 
+  hmap::Array un = hmap::unsphericity(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);
 }