r.li.patchdensity: fix number of patches

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

raster/r.li/r.li.patchdensity/main.c

@@ -59,30 +59,36 @@ int main(int argc, char *argv[])
 
 int patch_density(int fd, char **par, struct area_entry *ad, double *result)
 {
-    CELL *buf, *sup;
-    int count = 0, i, j, connected = 0, complete_line = 1, other_above = 0;
-    double area;
+    CELL *buf, *sup, *cnull;
+    CELL pid, old_pid, *pid_curr, *pid_sup, *ctmp;
+    int count, i, j, k, connected, other_above;
     struct Cell_head hd;
-    CELL complete_value;
+    int mask_fd, *mask_buf, null_count;
+    double area;
     double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
-    int mask_fd = -1, *mask_buf, *mask_sup, null_count = 0;
-
-    G_debug(1, "begin patch_density() index");
 
-    Rast_set_c_null_value(&complete_value, 1);
     Rast_get_cellhd(ad->raster, "", &hd);
 
+    cnull = Rast_allocate_c_buf();
+    Rast_set_c_null_value(cnull, Rast_window_cols());
+    sup = cnull;
+
+    /* initialize patch ids */
+    pid_curr = Rast_allocate_c_buf();
+    Rast_set_c_null_value(pid_curr, Rast_window_cols());
+    pid_sup = Rast_allocate_c_buf();
+    Rast_set_c_null_value(pid_sup, Rast_window_cols());
+
     /* open mask if needed */
+    mask_fd = -1;
+    mask_buf = NULL;
+    null_count = 0;
     if (ad->mask == 1) {
 	if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
 	    return 0;
-
-	mask_buf = malloc(ad->cl * sizeof(int));
-	mask_sup = malloc(ad->cl * sizeof(int));
+	mask_buf = G_malloc(ad->cl * sizeof(int));
     }
 
-    sup = Rast_allocate_c_buf();
-
     /* calculate distance */
     G_begin_distance_calculations();
     /* EW Dist at North edge */
@@ -94,103 +100,96 @@ int patch_density(int fd, char **par, struct area_entry *ad, double *result)
     /* NS Dist at West edge */
     NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
 
-
-    /* calculate number of patch */
+    /* calculate number of patches */
+    count = 0;
+    connected = 0;
+    other_above = 0;
+    pid = 0;
 
     for (i = 0; i < ad->rl; i++) {
 	buf = RLI_get_cell_raster_row(fd, i + ad->y, ad);
-
 	if (i > 0) {
 	    sup = RLI_get_cell_raster_row(fd, i - 1 + ad->y, ad);
 	}
-
 	/* mask values */
 	if (ad->mask == 1) {
-	    int k;
-
-	    if (i > 0) {
-		int *tmp;
-
-		tmp = mask_sup;
-		mask_buf = mask_sup;
-	    }
 	    if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
 		return 0;
 	    for (k = 0; k < ad->cl; k++) {
 		if (mask_buf[k] == 0) {
-		    Rast_set_c_null_value(mask_buf + k, 1);
+		    Rast_set_c_null_value(&buf[k + ad->x], 1);
 		    null_count++;
 		}
 	    }
-
 	}
+	
+	ctmp = pid_sup;
+	pid_sup = pid_curr;
+	pid_curr = ctmp;
+	Rast_set_c_null_value(pid_curr, Rast_window_cols());
+
+	connected = 0;
+	other_above = 1;
+	for (j = 0; j < ad->cl; j++) {
+	    
+	    if (Rast_is_null_value(&(buf[j + ad->x]), CELL_TYPE)) {
+		connected = 0;
+		other_above = 1;
+		continue;
+	    }
 
-
-	if (complete_line) {
-	    if (!Rast_is_null_value(&(buf[ad->x]), CELL_TYPE) &&
-		buf[ad->x] != complete_value)
-		count++;
-
-	    for (j = 0; j < ad->cl - 1; j++) {
-
-		if (buf[j + ad->x] != buf[j + 1 + ad->x]) {
-		    complete_line = 0;
-		    if (!Rast_is_null_value(&(buf[j + 1 + ad->x]), CELL_TYPE)
-			&& buf[j + 1 + ad->x] != complete_value)
-			count++;
+	    if (sup[j + ad->x] == buf[j + ad->x]) {
+		
+		if (!connected) {
+		    pid_curr[j + ad->x] = pid_sup[j + ad->x];
 		}
-
-	    }
-	    if (complete_line) {
-		complete_value = buf[ad->x];
-	    }
-	}
-	else {
-	    complete_line = 1;
-	    connected = 0;
-	    other_above = 0;
-	    for (j = 0; j < ad->cl; j++) {
-		if (sup[j + ad->x] == buf[j + ad->x]) {
-		    connected = 1;
-		    if (other_above) {
-			other_above = 0;
+		
+		if (pid_curr[j + ad->x] != pid_sup[j + ad->x]) {
+		    if (connected) {
 			count--;
 		    }
-		}
-		else {
-		    if (connected &&
-			!Rast_is_null_value(&(buf[j + ad->x]), CELL_TYPE))
-			other_above = 1;
-		}
-		if (j < ad->cl - 1 && buf[j + ad->x] != buf[j + 1 + ad->x]) {
-		    complete_line = 0;
-		    if (!connected &&
-			!Rast_is_null_value(&(buf[j + ad->x]), CELL_TYPE)) {
-
-			count++;
-			connected = 0;
-			other_above = 0;
+		    if (other_above) {
+			pid_curr[j + ad->x] = pid_sup[j + ad->x];
+			for (k = j + ad->x - 1; k >= ad->x; k--) {
+			    if (buf[k] != buf[j + ad->x])
+				break;
+			    pid_curr[k] = pid_sup[j + ad->x];
+			}
 		    }
 		    else {
-			connected = 0;
-			other_above = 0;
+			old_pid = pid_sup[j + ad->x];
+			pid_sup[j + ad->x] = pid_curr[j + ad->x];
+			
+			for (k = j + 1; k < ad->cl; k++) {
+			    if (pid_sup[k + ad->x] == old_pid) {
+				pid_sup[k + ad->x] = pid_curr[j + ad->x];
+			    }
+			}
 		    }
 		}
-	    }
-	    if (!connected &&
-		sup[ad->cl - 1 + ad->x] != buf[ad->cl - 1 + ad->x]) {
-		if (!Rast_is_null_value
-		    (&(buf[ad->cl - 1 + ad->x]), CELL_TYPE)) {
-		    count++;
-		    complete_line = 0;
-		}
+
+		other_above = 0;
+		connected = 1;
 	    }
 
-	    if (complete_line)
-		complete_value = buf[ad->x];
+	    if (!connected) {
+		count++;
+		pid++;
+		pid_curr[j + ad->x] = pid;
+	    }
 
+	    if (j < ad->cl - 1) {
+		if (buf[j + ad->x] == buf[j + 1 + ad->x]) {
+		    
+		    connected = 1;
+		    pid_curr[j + 1 + ad->x] = pid_curr[j + ad->x];
+		}
+		else {
+		    other_above = 1;
+		    connected = 0;
+		}
+	    }
 	}
-
     }
 
     area = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
@@ -200,9 +199,15 @@ int patch_density(int fd, char **par, struct area_entry *ad, double *result)
     if (area != 0)
 	*result = (count / area) * 1000000;
     else
-	*result = -1;
+	Rast_set_d_null_value(result, 1);
+
+    if (ad->mask == 1) {
+	G_free(mask_buf);
+    }
 
-    G_free(sup);
+    G_free(cnull);
+    G_free(pid_curr);
+    G_free(pid_sup);
 
     return RLI_OK;
 }

raster/r.li/r.li.patchdensity/r.li.patchdensity.html

@@ -5,8 +5,10 @@
 f(sample_area) = (Patch_Number/Area) * 1000000
 </pre></div>
 
-that is 1000000 by number of patch for area unit.
-This index is calculated using a 4 neighbour algorithm.
+that is the number of patches per 1,000,000 area units.
+<p>
+This index is calculated using a 4 neighbour algorithm, diagonal cells 
+are ignored when tracing a patch.
 
 <h2>NOTES</h2>
 
@@ -17,11 +19,12 @@ the <tt>$HOME/.grass7/r.li/output/</tt> folder.
 <p>
 <!-- TODO: verify next: -->
 A map of NULL values is considered to have zero patches. <br>
-If raster area is 0, <em>r.li.patchdensity</em> returns -1. This is only
-possible if the raster is masked.<br>
-If you want to change these -1 values to NULL, run subsequently on the resulting map:
+If the size of the sample area is 0, <em>r.li.patchdensity</em> returns 
+NULL. This is only possible if the raster is masked.<br>
+If you want to change these NULL values to some number, run subsequently 
+on the resulting map:
 <div class="code"><pre>
-r.null setnull=-1 input=my_out
+r.null null=-1 input=my_out
 </pre></div>
 after index calculation.