r.geomorphon: Add two more comparison modes (#1096)

The geomorphon computation includes a procedure done for each of the
eight cardinal directions. The procedure traces through the elevation
data to find the angles of zenith and nadir line of sight and to decide
if the cardinal point is above (1), below (-1) or level (0) with the
central point.

The existing nadir/zenith comparison defaults to 0 when zenith and nadir
angles are exactly equal and at least one angle is over its threshold.
For example, this would be the case if the zenith point was 20 raster
cells away at height +20m and the nadir point was 1 raster cell away at
height -1m. This is much more likely to happen when the elevation data
is low-resolution. Specifically, at 30m raster resolution and 1m
vertical resolution this happened every 1 out of 25 cases on average.

In that edge case the computed distance to the cardinal point remains
0, which results in incorrect derived cartesian coordinates and related
results, most notably azimuth and elongation. Besides that, the
comparison does not specifically check that when it comes to a 1 or -1
result, the corresponding angle is above its threshold.

Address these issues by introducing a new "comparison" command-line
option, in which "anglev1" stands for the exact existing behaviour and
"anglev2" and "anglev2_distance" enable more sophisticated comparison
methods.

raster/r.geomorphon/local_proto.h

@@ -82,6 +82,13 @@ GLOBAL double search_distance, flat_distance;
 GLOBAL double flat_threshold, flat_threshold_height;
 GLOBAL struct Cell_head window;
 GLOBAL int cell_step;
+/* Zenith/nadir comparison modes. */
+GLOBAL enum
+{
+    ANGLEV1,
+    ANGLEV2,
+    ANGLEV2_DISTANCE
+} compmode;
 
 /* memory */
 int open_map(MAPS * rast);

raster/r.geomorphon/main.c

@@ -84,6 +84,7 @@ int main(int argc, char **argv)
         *par_skip_radius,
         *par_flat_threshold,
         *par_flat_distance,
+        *par_comparison,
         *par_multi_prefix, *par_multi_step, *par_multi_start;
     struct Flag *flag_units, *flag_extended;
 
@@ -151,6 +152,15 @@ int main(int argc, char **argv)
         par_flat_distance->description =
             _("Flatness distance, zero for none");
 
+        par_comparison = G_define_option();
+        par_comparison->key = "comparison";
+        par_comparison->type = TYPE_STRING;
+        par_comparison->options = "anglev1,anglev2,anglev2_distance";
+        par_comparison->answer = "anglev1";
+        par_comparison->required = NO;
+        par_comparison->description =
+            _("Comparison mode for zenith/nadir line-of-sight search");
+
         par_multi_prefix = G_define_option();
         par_multi_prefix->key = "prefix";
         par_multi_prefix->type = TYPE_STRING;
@@ -193,6 +203,14 @@ int main(int argc, char **argv)
         double ns_resolution;
 
         multires = (par_multi_prefix->answer) ? 1 : 0;
+        if (!strcmp(par_comparison->answer, "anglev1"))
+            compmode = ANGLEV1;
+        else if (!strcmp(par_comparison->answer, "anglev2"))
+            compmode = ANGLEV2;
+        else if (!strcmp(par_comparison->answer, "anglev2_distance"))
+            compmode = ANGLEV2_DISTANCE;
+        else
+            G_fatal_error(_("Failed parsing <%s>"), par_comparison->answer);
         for (i = o_forms; i < o_size; ++i)      /* check for outputs */
             if (opt_output[i]->answer) {
                 if (G_legal_filename(opt_output[i]->answer) < 0)

raster/r.geomorphon/pattern.c

@@ -7,6 +7,65 @@
 static int nextr[NUM_DIRS] = { -1, -1, -1, 0, 1, 1, 1, 0 };
 static int nextc[NUM_DIRS] = { 1, 0, -1, -1, -1, 0, 1, 1 };
 
+/*
+ * A more thorough comparison using a few factors of different priority
+ * similar to the BGP best path selection algorithm. When the distances are
+ * equal, it becomes an improved version of the original r.geomorphon
+ * comparison (which is different from the original paper), in that it applies
+ * each threshold to its respective angle and does not default to 0 when there
+ * is a tie. Each angle must be non-negative.
+ */
+static int compare_multi(const double nadir_angle, const double zenith_angle,
+                         const double nadir_threshold,
+                         const double zenith_threshold,
+                         const double nadir_distance,
+                         const double zenith_distance)
+{
+    const unsigned char
+        nadir_over = nadir_angle > nadir_threshold,
+        zenith_over = zenith_angle > zenith_threshold;
+
+    /*
+     * If neither angle exceeds its threshold, consider the elevation profile
+     * flat enough.
+     */
+    if (!nadir_over && !zenith_over)
+        return 0;
+    /*
+     * If exactly one angle exceeds its threshold, that angle represents the
+     * elevation profile.
+     */
+    if (!nadir_over && zenith_over)
+        return 1;
+    if (nadir_over && !zenith_over)
+        return -1;
+    /*
+     * If both angles exceed their thresholds, the greater angle (if any)
+     * represents the elevation profile better.
+     */
+    if (nadir_angle < zenith_angle)
+        return 1;
+    if (nadir_angle > zenith_angle)
+        return -1;
+    /*
+     * Here each angle is above its threshold and the angles are exactly equal
+     * (which happens quite often when the elevation values are integer rather
+     * than real). Consider the angle computed over the greater distance to
+     * represent the elevation profile better since it is based on a greater
+     * number of cells.
+     */
+    if (nadir_distance < zenith_distance)
+        return 1;
+    if (nadir_distance > zenith_distance)
+        return -1;
+    /*
+     * If there is still a tie, 0 would not be a valid result because both
+     * angles exceed their thresholds hence the elevation profile definitely
+     * is not flat. Resolve this with a preferred constant value.
+     */
+    return 1;
+}
+
 int calc_pattern(PATTERN * pattern, int row, int cur_row, int col)
 {
     /* calculate parameters of geomorphons and store it in the struct pattern */
@@ -121,6 +180,48 @@ int calc_pattern(PATTERN * pattern, int row, int cur_row, int col)
         if (nadir_angle < -nadir_threshold)
             pattern->negatives += i;
 
+        if (compmode != ANGLEV1) {
+            /*
+             * One of the differences from ANGLEV1 is that even if there is a
+             * tie, the second switch block will eventually use one of the
+             * distances instead of 0 to set the cardinal point distance.
+             */
+            switch (compmode) {
+            case ANGLEV2:
+                pattern->pattern[i] =
+                    compare_multi(fabs(nadir_angle), fabs(zenith_angle),
+                                  nadir_threshold, zenith_threshold, 0, 0);
+                break;
+            case ANGLEV2_DISTANCE:
+                pattern->pattern[i] =
+                    compare_multi(fabs(nadir_angle), fabs(zenith_angle),
+                                  nadir_threshold, zenith_threshold,
+                                  nadir_distance, zenith_distance);
+                break;
+            default:
+                G_fatal_error(_("Internal error in %s()"), __func__);
+            }
+
+            switch (pattern->pattern[i]) {
+            case 1:
+                pattern->elevation[i] = zenith_height;
+                pattern->distance[i] = zenith_distance;
+                pattern->num_positives++;
+                break;
+            case -1:
+                pattern->elevation[i] = nadir_height;
+                pattern->distance[i] = nadir_distance;
+                pattern->num_negatives++;
+                break;
+            case 0:
+                pattern->distance[i] = search_distance;
+                break;
+            }
+
+            continue;
+        }                       /* if (compmode != ANGLEV1) */
+
+        /* ANGLEV1 */
         if (fabs(zenith_angle) > zenith_threshold ||
             fabs(nadir_angle) > nadir_threshold) {
             if (fabs(nadir_angle) < fabs(zenith_angle)) {
@@ -135,6 +236,10 @@ int calc_pattern(PATTERN * pattern, int row, int cur_row, int col)
                 pattern->distance[i] = nadir_distance;
                 pattern->num_negatives++;
             }
+            /*
+             * If the angles are exactly equal, the cardinal direction search
+             * results are the values set at the beginning of the for() loop.
+             */
         }
         else {
             pattern->distance[i] = search_distance;

raster/r.geomorphon/r.geomorphon.html

@@ -67,6 +67,13 @@ <h2>OPTIONS</h2>
 <DD>The difference (in degrees) between zenith and nadir line-of-sight which indicate flat direction. If higher threshold produce more flat maps. If resolution of the map is low (more than 1 km per cell) threshold should be very small (much smaller than 1 degree) because on such distance 1 degree of difference means several meters of high difference.</DD>
 <DT><b>dist</b></DT>
 <DD>>Flat distance. This is additional parameter defining the distance above which the threshold starts to decrease to avoid problems with pseudo-flat line-of-sights if real elevation difference appears on the distance where its value is higher (TO BE CORRECTED).</DD>
+<DT><b>comparison</b></DT>
+<DD>Comparison mode for zenith/nadir line-of-sight search. "anglev1" is
+the original r.geomorphon comparison mode. "anglev2" is an improved
+mode, which better handles angle thresholds and zenith/nadir angles
+that are exactly equal. "anglev2_distance" in addition to that takes
+the zenith/nadir distances into account when the angles are exactly
+equal.</DD>
 <DT><b>forms</b></DT>
 <DD>Returns geomorphic map with 10 most popular terrestrial forms. Legend for forms, its definition by the number of <em>+</em> and <em>-</em> and its idealized  visualisation are presented at the image.
 <center>