beautify flat areas: more tuning, clean up

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@42269 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.watershed/ram/Gwater.h

@@ -101,7 +101,7 @@ double get_slope(int, int, int, int, CELL, CELL);
 int replace(int, int, int, int);
 
 /* do_flatarea.c */
-int do_flatarea(int, CELL);
+int do_flatarea(int, CELL, CELL *, CELL *);
 
 /* do_cum.c */
 int do_cum(void);

raster/r.watershed/ram/do_astar.c

@@ -104,10 +104,11 @@ int do_astar(void)
 		skip_diag = 0;
 		
 		alt_nbr[ct_dir] = alt[index_up];
-		if (flat_flag && !is_worked) {
-		    alt_val = alt[index_doer];
+		if (flat_flag && !is_in_list && !is_worked) {
+		    alt_val = alt_bak[index_doer];
+		    alt_nbr[ct_dir] = alt_bak[index_up];
 		    if (!flat_is_done && alt_nbr[ct_dir] == alt_val) {
-			do_flatarea(index_doer, alt_val);
+			do_flatarea(index_doer, alt_val, alt_bak, alt);
 			alt_nbr[ct_dir] = alt[index_up];
 			flat_is_done = 1;
 			nbr_flat_is_done = 1;

raster/r.watershed/ram/do_flatarea.c

@@ -49,26 +49,21 @@ struct pq *pq_create(void)
 int pq_add(int idx, struct pq *q)
 {
     assert(q->last);
-    if (q->last->idx != -1)
-	G_fatal_error("idx is %d", q->last->idx);
+    assert(q->last->idx == -1);
 
     q->last->idx = idx;
-    if (q->last->next == NULL) {
-	struct pq_node *n = (struct pq_node *) G_malloc(sizeof(struct pq_node));
-	n->next = NULL;
-	n->idx = -1;
-	q->last->next = n;
-	q->last = q->last->next;
-    }
-    else {
-	q->last = q->last->next;
-	assert(NULL);
+    if (q->last->next != NULL) {
+	G_fatal_error(_("Beautify flat areas: priority queue error"));
     }
 
+    struct pq_node *n = (struct pq_node *) G_malloc(sizeof(struct pq_node));
+    n->next = NULL;
+    n->idx = -1;
+    q->last->next = n;
+    q->last = q->last->next;
+
     assert(q->last != q->last->next);
     assert(q->first != q->last);
-    assert(q->last->idx == -1);
-    assert(q->last->next == NULL);
     q->size++;
 
     return 0;
@@ -79,8 +74,6 @@ int pq_drop(struct pq *q)
     int idx = q->first->idx;
     struct pq_node *n = q->first;
 
-    assert(q->first != q->first->next);
-    assert(q->size);
     q->size--;
 
     q->first = q->first->next;
@@ -119,8 +112,11 @@ int cmp_orders(const void *a, const void *b)
     
     return (oa->index < ob->index ? -1 : (oa->index > ob->index));
 }
-
-int do_flatarea(int index, CELL ele)
+/*
+ * return 0 if nothing was modidied
+ * return 1 if elevation was modified
+ */
+int do_flatarea(int index, CELL ele, CELL *alt_org, CELL *alt_new)
 {
     int upr, upc, r, c, ct_dir;
     CELL is_in_list, is_worked, this_in_list;
@@ -137,18 +133,20 @@ int do_flatarea(int index, CELL ele)
     struct RB_TREE *order_tree = rbtree_create(cmp_orders, sizeof(struct orders));
 
     pq_add(index, down_pq);
+    pq_add(index, up_pq);
     inc_order.downhill = -1;
     inc_order.uphill = 0;
     inc_order.index = index;
-    inc_order.flag = 1;
+    inc_order.flag = 0;
     rbtree_insert(order_tree, &inc_order);
 
-    n_flat_cells++;
+    n_flat_cells = 1;
 
     min_ele_diff = INT_MAX;
     max_uphill_order = max_downhill_order = 0;
 
     /* get uphill start points */
+    G_debug(2, "get uphill start points");
     counter = 0;
     while (down_pq->size) {
 	if ((index_doer = pq_drop(down_pq)) == -1)
@@ -156,6 +154,8 @@ int do_flatarea(int index, CELL ele)
 
 	seg_index_rc(alt_seg, index_doer, &r, &c);
 
+	FLAG_SET(flat_done, r, c);
+
 	/* check all neighbours, breadth first search */
 	for (ct_dir = 0; ct_dir < sides; ct_dir++) {
 	    /* get r, c (upr, upc) for this neighbour */
@@ -166,7 +166,7 @@ int do_flatarea(int index, CELL ele)
 		index_up = SEG_INDEX(alt_seg, upr, upc);
 		is_in_list = FLAG_GET(in_list, upr, upc);
 		is_worked = FLAG_GET(worked, upr, upc);
-		ele_nbr = alt[index_up];
+		ele_nbr = alt_org[index_up];
 
 		if (ele_nbr == ele && !is_worked) {
 
@@ -186,7 +186,6 @@ int do_flatarea(int index, CELL ele)
 			    pq_add(index_up, up_pq);
 			    /* set uphill order to 0 */
 			    inc_order.uphill = 0;
-			    inc_order.flag = 1;
 			    counter++;
 			}
 			rbtree_insert(order_tree, &inc_order);
@@ -195,22 +194,24 @@ int do_flatarea(int index, CELL ele)
 	    }
 	}
     }
-
+    /* flat area too small, not worth the effort */
     if (n_flat_cells < 5) {
 	/* clean up */
 	pq_destroy(up_pq);
 	pq_destroy(down_pq);
 	rbtree_destroy(order_tree);
 
-	return 1;
+	return 0;
     }
     
-    G_debug(2, "%d flat cells, %d cells in tree, %d start cells", n_flat_cells, (int)order_tree->count, counter);
+    G_debug(2, "%d flat cells, %d cells in tree, %d start cells",
+               n_flat_cells, (int)order_tree->count, counter);
 
     pq_destroy(down_pq);
     down_pq = pq_create();
 
     /* got uphill start points, do uphill correction */
+    G_debug(2, "got uphill start points, do uphill correction");
     counter = 0;
     uphill_order = 1;
     while (up_pq->size) {
@@ -231,7 +232,8 @@ int do_flatarea(int index, CELL ele)
 	uphill_order = order_found->uphill;
 
 	if (last_order > uphill_order)
-	    G_warning("queue error: last uphill order %d > current uphill order %d", last_order, uphill_order);
+	    G_warning(_("queue error: last uphill order %d > current uphill order %d"),
+	              last_order, uphill_order);
 
 	/* debug */
 	if (uphill_order == -1)
@@ -253,9 +255,10 @@ int do_flatarea(int index, CELL ele)
 		index_up = SEG_INDEX(alt_seg, upr, upc);
 		is_in_list = FLAG_GET(in_list, upr, upc);
 		is_worked = FLAG_GET(worked, upr, upc);
-		ele_nbr = alt[index_up];
+		ele_nbr = alt_org[index_up];
 
-		/* all cells that are in_list should have been added previously as uphill start points */
+		/* all cells that are in_list should have been added
+		 * previously as uphill start points */
 		if (ele_nbr == ele && !is_worked) {
 
 		    inc_order.index = index_up;
@@ -264,14 +267,13 @@ int do_flatarea(int index, CELL ele)
 		    }
 
 		    /* not yet added to queue */
-		    if (nbr_order_found->flag != 1) {
+		    if (nbr_order_found->uphill == -1) {
 			if (is_in_list)
 			    G_warning("cell should be in queue");
 			/* add to up queue */
 			pq_add(index_up, up_pq);
 			/* set nbr uphill order = current uphill order + 1 */
 			nbr_order_found->uphill = uphill_order;
-			nbr_order_found->flag = 1;
 		    }
 		}
 		/* add focus cell to down queue */
@@ -287,18 +289,20 @@ int do_flatarea(int index, CELL ele)
 	    }
 	}
     }
-    /* debug: all flags should be set to 1 */
+    /* debug: all flags should be set to 0 */
 
     pq_destroy(up_pq);
     up_pq = pq_create();
 
     /* got downhill start points, do downhill correction */
+    G_debug(2, "got downhill start points, do downhill correction");
     downhill_order = 1;
     while (down_pq->size) {
 	if ((index_doer = pq_drop(down_pq)) == -1)
-	    G_fatal_error("downhill order: no more points in down queue");
+	    G_fatal_error(_("downhill order: no more points in down queue"));
 
 	seg_index_rc(alt_seg, index_doer, &r, &c);
+	this_in_list = FLAG_GET(in_list, r, c);
 
 	/* get downhill order for this point */
 	inc_order.index = index_doer;
@@ -309,12 +313,16 @@ int do_flatarea(int index, CELL ele)
 	downhill_order = order_found->downhill;
 
 	if (last_order > downhill_order)
-	    G_warning("queue error: last downhill order %d > current downhill order %d", last_order, downhill_order);
+	    G_warning(_("queue error: last downhill order %d > current downhill order %d"),
+	              last_order, downhill_order);
 
 	/* debug */
 	if (downhill_order == -1)
 	    G_fatal_error(_("downhill order: downhill order not set"));
 
+	if (max_downhill_order < downhill_order)
+	    max_downhill_order = downhill_order;
+
 	downhill_order++;
 
 	/* check all neighbours, breadth first search */
@@ -327,7 +335,7 @@ int do_flatarea(int index, CELL ele)
 		index_up = SEG_INDEX(alt_seg, upr, upc);
 		is_in_list = FLAG_GET(in_list, upr, upc);
 		is_worked = FLAG_GET(worked, upr, upc);
-		ele_nbr = alt[index_up];
+		ele_nbr = alt_org[index_up];
 
 		if (ele_nbr == ele && !is_worked) {
 
@@ -343,15 +351,11 @@ int do_flatarea(int index, CELL ele)
 			/* set nbr downhill order = current downhill order + 1 */
 			nbr_order_found->downhill = downhill_order;
 
-			if (max_downhill_order < downhill_order)
-			    max_downhill_order = downhill_order;
-
-
 			/* add to up queue */
 			if (is_in_list) {
 			    pq_add(index_up, up_pq);
-			    /* unset flag */
-			    nbr_order_found->flag = 0;
+			    /* set flag */
+			    nbr_order_found->flag = 1;
 			}
 		    }
 		}
@@ -363,8 +367,8 @@ int do_flatarea(int index, CELL ele)
 
     /* increment: ele += uphill_order +  max_downhill_order - downhill_order */
     /* decrement: ele += uphill_order - max_uphill_order - downhill_order */
-    //max_downhill_order++;
 
+    G_debug(2, "adjust ele");
     while (up_pq->size) {
 	if ((index_doer = pq_drop(up_pq)) == -1)
 	    G_fatal_error("no more points in up queue");
@@ -372,8 +376,6 @@ int do_flatarea(int index, CELL ele)
 	seg_index_rc(alt_seg, index_doer, &r, &c);
 	this_in_list = FLAG_GET(in_list, r, c);
 
-	FLAG_SET(flat_done, r, c);
-
 	/* get uphill and downhill order for this point */
 	inc_order.index = index_doer;
 	if ((order_found = rbtree_find(order_tree, &inc_order)) == NULL)
@@ -385,7 +387,7 @@ int do_flatarea(int index, CELL ele)
 	/* debug */
 	if (uphill_order == -1)
 	    G_fatal_error(_("adjustment: uphill order not set"));
-	if (downhill_order == -1)
+	if (!this_in_list && downhill_order == -1)
 	    G_fatal_error(_("adjustment: downhill order not set"));
 
 	/* increment */
@@ -393,7 +395,8 @@ int do_flatarea(int index, CELL ele)
 	    downhill_order = max_downhill_order;
 	    uphill_order = 0;
 	}
-	alt[index_doer] += uphill_order + (double)(max_downhill_order - downhill_order) / 2.0 + 0.5;
+	alt_new[index_doer] +=
+	    (uphill_order + (double)(max_downhill_order - downhill_order) / 2.0 + 0.5) / 2.0 + 0.5;
 
 	/* check all neighbours, breadth first search */
 	for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -405,7 +408,7 @@ int do_flatarea(int index, CELL ele)
 		index_up = SEG_INDEX(alt_seg, upr, upc);
 		is_in_list = FLAG_GET(in_list, upr, upc);
 		is_worked = FLAG_GET(worked, upr, upc);
-		ele_nbr = alt[index_up];
+		ele_nbr = alt_org[index_up];
 
 		if (ele_nbr == ele && !is_worked) {
 
@@ -414,12 +417,12 @@ int do_flatarea(int index, CELL ele)
 			G_fatal_error(_("flat cell escaped in adjustment"));
 
 		    /* not yet added to queue */
-		    if (nbr_order_found->flag != 0) {
+		    if (nbr_order_found->flag == 0) {
 			if (is_in_list)
 			    G_warning("adjustment: in_list cell should be in queue");
 			/* add to up queue */
 			pq_add(index_up, up_pq);
-			nbr_order_found->flag = 0;
+			nbr_order_found->flag = 1;
 		    }
 		}
 	    }
@@ -431,5 +434,5 @@ int do_flatarea(int index, CELL ele)
     pq_destroy(down_pq);
     rbtree_destroy(order_tree);
     
-    return 0;
+    return 1;
 }