r.topidx: fix bug with steep slopes (#193)

fixed bug which showed up if there are negative TWI values (happens in very steep slope)

raster/r.topidx/global.h

@@ -1,3 +1,5 @@
+#include <float.h>
+
 #include <grass/gis.h>
 #include <grass/raster.h>
 
@@ -6,8 +8,13 @@
 #define	atbv(i,j)		atb[i][j]
 #define	is_cv_null(i,j)		Rast_is_d_null_value(&cv(i,j))
 #define	is_atbv_null(i,j)	Rast_is_d_null_value(&atbv(i,j))
+#define is_atbv_unprocessed(i,j) (atbv(i,j) == UNPROCESSED)
 
+#ifndef DBL_MAX
+#define DBL_MAX 1.797693E308  /* DBL_MAX approximation */
+#endif
 #define	ZERO			0.0000001
+#define	UNPROCESSED		-DBL_MAX
 
 #ifdef _MAIN_C_
 #define GLOBAL

raster/r.topidx/topidx.c

@@ -19,7 +19,7 @@ void initialize(void)
 		Rast_set_d_null_value(&atbv(i, j), 1);
 	    }
 	    else {
-		atbv(i, j) = -10.0;
+		atbv(i, j) = UNPROCESSED;
 	    }
 	}
     }
@@ -55,7 +55,7 @@ void calculate_atanb(void)
 		    continue;
 
 		/* skip squares already done */
-		if (is_atbv_null(i, j) || atbv(i, j) >= ZERO)
+		if (is_atbv_null(i, j) || !is_atbv_unprocessed(i,j))
 		    continue;
 
 		/* check the 8 possible flow directions for
@@ -66,47 +66,47 @@ void calculate_atanb(void)
 			(is_cv_null(i - 1, j - 1) ||
 			 cv(i - 1, j - 1) > cv(i, j)) &&
 			!is_atbv_null(i - 1, j - 1) &&
-			atbv(i - 1, j - 1) < ZERO)
+			is_atbv_unprocessed(i - 1, j - 1))
 			continue;
 
 		    if ((is_cv_null(i - 1, j) ||
 			 cv(i - 1, j) > cv(i, j)) &&
-			!is_atbv_null(i - 1, j) && atbv(i - 1, j) < ZERO)
+			!is_atbv_null(i - 1, j) && is_atbv_unprocessed(i - 1, j))
 			continue;
 
 		    if (j + 1 < window.cols &&
 			(is_cv_null(i - 1, j + 1) ||
 			 cv(i - 1, j + 1) > cv(i, j)) &&
 			!is_atbv_null(i - 1, j + 1) &&
-			atbv(i - 1, j + 1) < ZERO)
+			is_atbv_unprocessed(i - 1, j + 1))
 			continue;
 		}
 		if (j > 0 &&
 		    (is_cv_null(i, j - 1) ||
 		     cv(i, j - 1) > cv(i, j)) &&
-		    !is_atbv_null(i, j - 1) && atbv(i, j - 1) < ZERO)
+		    !is_atbv_null(i, j - 1) && is_atbv_unprocessed(i, j - 1))
 		    continue;
 		if (j + 1 < window.cols &&
 		    (is_cv_null(i, j + 1) ||
 		     cv(i, j + 1) > cv(i, j)) &&
-		    !is_atbv_null(i, j + 1) && atbv(i, j + 1) < ZERO)
+		    !is_atbv_null(i, j + 1) && is_atbv_unprocessed(i, j + 1))
 		    continue;
 		if (i + 1 < window.rows) {
 		    if (j > 0 &&
 			(is_cv_null(i + 1, j - 1) ||
 			 cv(i + 1, j - 1) > cv(i, j)) &&
 			!is_atbv_null(i + 1, j - 1) &&
-			atbv(i + 1, j - 1) < ZERO)
+			is_atbv_unprocessed(i + 1, j - 1))
 			continue;
 		    if ((is_cv_null(i + 1, j) ||
 			 cv(i + 1, j) > cv(i, j)) &&
-			!is_atbv_null(i + 1, j) && atbv(i + 1, j) < ZERO)
+			!is_atbv_null(i + 1, j) && is_atbv_unprocessed(i + 1, j))
 			continue;
 		    if (j + 1 < window.cols &&
 			(is_cv_null(i + 1, j + 1) ||
 			 cv(i + 1, j + 1) > cv(i, j)) &&
 			!is_atbv_null(i + 1, j + 1) &&
-			atbv(i + 1, j + 1) < ZERO)
+			is_atbv_unprocessed(i + 1, j + 1))
 			continue;
 		}
 		/* find the outflow directions and calculate