lib/arraystats: re-enable "Discont" algorithm

Addresses heap buffer overflow in previously disabled 'Discont'
algorithm.

The immediate cause was too small allocation of variable `num` in
`AS_class_discont()`, which was fixed by an increase it's size (broadly
mirroring the original Fortran code).

This fix was complemented by some code related improvements:

- Code clean up and restructure of `AS_class_discont()`
- Re-enable "dis" algorithm in v.class and d.vect.thematic
- Pass `classbreaks` array by pointer to pointer to enable in- and
  output.

See: https://lists.osgeo.org/pipermail/grass-dev/2008-July/038951.html

Fixes: OSGeo#5486

Author: nilason

display/d.vect.thematic/main.c

@@ -119,11 +119,10 @@ int main(int argc, char **argv)
     algo_opt->options = "int,std,qua,equ,dis";
     algo_opt->description = _("Algorithm to use for classification");
     desc = NULL;
-    G_asprintf(&desc, "int;%s;std;%s;qua;%s;equ;%s", _("simple intervals"),
-               _("standard deviations"), _("quantiles"),
-               _("equiprobable (normal distribution)"));
+    G_asprintf(&desc, "int;%s;std;%s;qua;%s;equ;%s;dis;%s",
+               _("simple intervals"), _("standard deviations"), _("quantiles"),
+               _("equiprobable (normal distribution)"), _("discontinuities"));
     algo_opt->descriptions = desc;
-    /*currently disabled because of bugs       "dis;discontinuities"); */
     algo_opt->guisection = _("Classes");
 
     nbclass_opt = G_define_option();
@@ -415,7 +414,7 @@ int main(int argc, char **argv)
              * algorithms */
             class_info =
                 AS_class_apply_algorithm(AS_option_to_algorithm(algo_opt), data,
-                                         nrec, &nbreaks, breakpoints);
+                                         nrec, &nbreaks, &breakpoints);
         }
         else {
 

include/grass/defs/arraystats.h

@@ -2,17 +2,17 @@
 #define GRASS_ARRAYSTATSDEFS_H
 
 /* basic.c */
-void AS_eqdrt(double[], double[], int, int, double *);
+void AS_eqdrt(double[], double[], int, int, double *, double *, double *);
 void AS_basic_stats(double *, int, struct GASTATS *);
 
 /* class.c */
 int AS_option_to_algorithm(const struct Option *);
-double AS_class_apply_algorithm(int, double *, int, int *, double *);
-int AS_class_interval(double *, int, int, double *);
-int AS_class_quant(double *, int, int, double *);
-double AS_class_discont(double *, int, int, double *);
-double AS_class_stdev(double *, int, int, double *);
-int AS_class_equiprob(double *, int, int *, double *);
+double AS_class_apply_algorithm(int, double *, int, int *, double **);
+int AS_class_interval(double *, int, int, double **);
+int AS_class_quant(double *, int, int, double **);
+double AS_class_discont(double *, int, int, double **);
+double AS_class_stdev(double *, int, int, double **);
+int AS_class_equiprob(double *, int, int *, double **);
 int AS_class_frequencies(double *, int, int, double *, int *);
 
 #endif

lib/arraystats/basic.c

@@ -1,5 +1,7 @@
 #include <math.h>
+
 #include <grass/arraystats.h>
+#include <grass/gis.h>
 
 /*provides basic univar stats */
 void AS_basic_stats(double *data, int count, struct GASTATS *stats)
@@ -34,31 +36,43 @@ void AS_basic_stats(double *data, int count, struct GASTATS *stats)
     return;
 }
 
-void AS_eqdrt(double vectx[], double vecty[], int i1, int i2, double *vabc)
+void AS_eqdrt(double vectx[], double vecty[], int i1, int i2, double *a,
+              double *b, double *c)
 {
-    double bn = 0, bd = 0, x1 = 0, y1 = 0;
+    double bn = 0.0, bd = 0.0, x1 = 0.0, y1 = 0.0;
+    double x2 = 0.0, y2 = 0.0;
+
+    *a = 0.0;
+    *b = 0.0;
+    *c = 0.0;
+
+    x1 = vectx[i1];
+    y1 = vecty[i1];
+    x2 = vectx[i2];
+    y2 = vecty[i2];
 
-    vabc[0] = 0;
-    vabc[1] = 0;
-    vabc[2] = 0;
     if (i1 == 0) {
-        x1 = 0;
-        y1 = 0;
+        x1 = 0.0;
+        y1 = 0.0;
     }
     else {
         x1 = vectx[i1];
         y1 = vecty[i1];
     }
-    bn = y1 - vecty[i2];
-    bd = x1 - vectx[i2];
-    if (bd != 0) {
-        vabc[1] = bn / bd;
-        vabc[0] = y1 - vabc[1] * x1;
-        return;
+
+    bn = y1 - y2;
+    bd = x1 - x2;
+
+    if (fabs(bd) < GRASS_EPSILON) {
+        if (fabs(bn) < GRASS_EPSILON) {
+            printf("Points are equal\n");
+        }
+        else {
+            *c = x1;
+        }
+    }
+    else {
+        *b = bn / bd;
+        *a = y1 - *b * x1;
     }
-    if (bn != 0)
-        vabc[2] = x1;
-    else
-        G_debug(3, "Points are equal\n");
-    return;
 }