[GRASS-SVN] r69874 - grass/branches/releasebranch_7_0/raster/r.texture

[svn_grass at osgeo.org]

Author: mmetz
Date: 2016-11-22 14:12:05 -0800 (Tue, 22 Nov 2016)
New Revision: 69874

Modified:
   grass/branches/releasebranch_7_0/raster/r.texture/h_measure.c
Log:
r.texture: handle possible division by zero (backport from trunk r69872)

Modified: grass/branches/releasebranch_7_0/raster/r.texture/h_measure.c
===================================================================
--- grass/branches/releasebranch_7_0/raster/r.texture/h_measure.c	2016-11-22 22:11:34 UTC (rev 69873)
+++ grass/branches/releasebranch_7_0/raster/r.texture/h_measure.c	2016-11-22 22:12:05 UTC (rev 69874)
@@ -407,21 +407,24 @@
 float f2_contrast(void)
 {
     int i, j /*, n */;
-    float sum, bigsum = 0;
+    float /* sum,*/ bigsum = 0;
     float **P = P_matrix;
 
+    /* the three-loop version does not work 
+     * when gray tones that do not occur in the current window 
+     * have been removed in tone and P* */
     /*
     for (n = 0; n < Ng; n++) {
 	sum = 0;
 	for (i = 0; i < Ng; i++) {
 	    for (j = 0; j < Ng; j++) {
-		if ((tone[i] - tone[j]) == tone[n] ||
-		    (tone[j] - tone[i]) == tone[n]) {
+		if ((i - j) == n ||
+		    (j - i) == n) {
 		    sum += P[i][j];
 		}
 	    }
 	}
-	bigsum += tone[n] * tone[n] * sum;
+	bigsum += n * n * sum;
     }
     */
 
@@ -462,6 +465,9 @@
 	    tmp += tone[i] * tone[j] * P[i][j];
     }
     stddev = sqrt(sum_sqr - (mean * mean));
+    
+    if (stddev == 0)	/* stddev < GRASS_EPSILON or similar ? */
+	return 0;
 
     return (tmp - mean * mean) / (stddev * stddev);
 }
@@ -594,8 +600,10 @@
 	sum_sqr += P[i] * P[i];
     }
     /* tmp = Ng * Ng; */
-    tmp = (tone[Ng - 1] - tone[0]) * (tone[Ng - 1] - tone[0]);
-    var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp);
+    if (Ng > 1) {
+	tmp = (tone[Ng - 1] - tone[0]) * (tone[Ng - 1] - tone[0]);
+	var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp);
+    }
 
     return var;
 }
@@ -639,6 +647,9 @@
     }
 
     /* fprintf(stderr,"hxy1=%f\thxy=%f\thx=%f\thy=%f\n",hxy1,hxy,hx,hy); */
+    if (hx == 0 && hy == 0)
+	return 0;
+
     return ((hxy - hxy1) / (hx > hy ? hx : hy));
 }