r.li.shape: calculate shape index

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

raster/r.li/r.li.shape/main.c

@@ -18,12 +18,24 @@
 
 #include <stdlib.h>
 #include <fcntl.h>
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
 #include <grass/glocale.h>
 #include "../r.li.daemon/daemon.h"
+#include "../r.li.daemon/GenericCell.h"
+
+/* type, cell count and edge count of each patch */
+struct pst {
+    generic_cell type;
+    long cells;
+    long edges;
+};
 
 rli_func shape_index;
+int calculate(int fd, struct area_entry *ad, double *result);
+int calculateD(int fd, struct area_entry *ad, double *result);
+int calculateF(int fd, struct area_entry *ad, double *result);
 
 int main(int argc, char *argv[])
 {
@@ -59,55 +71,754 @@ int main(int argc, char *argv[])
 
 int shape_index(int fd, char **par, struct area_entry *ad, double *result)
 {
-    double area;
-    struct Cell_head hd;
-    CELL complete_value;
-    double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
-    int mask_fd = -1, null_count = 0;
-    int i = 0, k = 0;
-    int *mask_buf = NULL;
+    int ris = RLI_OK;
+    double indice = 0;
+
+    switch (ad->data_type) {
+    case CELL_TYPE:
+	{
+	    ris = calculate(fd, ad, &indice);
+	    break;
+	}
+    case DCELL_TYPE:
+	{
+	    ris = calculateD(fd, ad, &indice);
+	    break;
+	}
+    case FCELL_TYPE:
+	{
+	    ris = calculateF(fd, ad, &indice);
+	    break;
+	}
+    default:
+	{
+	    G_fatal_error("data type unknown");
+	    return RLI_ERRORE;
+	}
+    }
+
+    if (ris != RLI_OK) {
+	return RLI_ERRORE;
+    }
 
-    Rast_set_c_null_value(&complete_value, 1);
+    *result = indice;
+
+    return RLI_OK;
+}
 
-    Rast_get_window(&hd);
+
+int calculate(int fd, struct area_entry *ad, double *result)
+{
+    CELL *buf, *buf_sup, *buf_null;
+    CELL corrCell, precCell, supCell;
+    long npatch, area; 
+    long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
+    struct pst *pst;
+    long nalloc, incr;
+    int i, j, k;
+    int connected;
+    int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
+
+    buf_null = Rast_allocate_c_buf();
+    Rast_set_c_null_value(buf_null, Rast_window_cols());
+    buf_sup = buf_null;
+
+    /* initialize patch ids */
+    pid_corr = G_malloc(ad->cl * sizeof(long));
+    pid_sup = G_malloc(ad->cl * sizeof(long));
+
+    for (j = 0; j < ad->cl; j++) {
+	pid_corr[j] = 0;
+	pid_sup[j] = 0;
+    }
 
     /* open mask if needed */
+    mask_fd = -1;
+    mask_buf = mask_sup = NULL;
+    masked = FALSE;
     if (ad->mask == 1) {
 	if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
-	    return 0;
+	    return RLI_ERRORE;
+	mask_buf = G_malloc(ad->cl * sizeof(int));
+	if (mask_buf == NULL) {
+	    G_fatal_error("malloc mask_buf failed");
+	    return RLI_ERRORE;
+	}
+	mask_sup = G_malloc(ad->cl * sizeof(int));
+	if (mask_sup == NULL) {
+	    G_fatal_error("malloc mask_buf failed");
+	    return RLI_ERRORE;
+	}
+	for (j = 0; j < ad->cl; j++)
+	    mask_buf[j] = 0;
+
+	masked = TRUE;
+    }
+
+    /* calculate number of patches */
+    npatch = 0;
+    area = 0;
+    pid = 0;
+
+    /* patch size and type */
+    incr = 1024;
+    if (incr > ad->rl)
+	incr = ad->rl;
+    if (incr > ad->cl)
+	incr = ad->cl;
+    if (incr < 2)
+	incr = 2;
+    nalloc = incr;
+    pst = G_malloc(nalloc * sizeof(struct pst));
+    for (k = 0; k < nalloc; k++) {
+	pst[k].cells = 0;
+	pst[k].edges = 0;
+    }
+
+    for (i = 0; i < ad->rl; i++) {
+	buf = RLI_get_cell_raster_row(fd, i + ad->y, ad);
+	if (i > 0) {
+	    buf_sup = RLI_get_cell_raster_row(fd, i - 1 + ad->y, ad);
+	}
+
+	if (masked) {
+	    mask_tmp = mask_sup;
+	    mask_sup = mask_buf;
+	    mask_buf = mask_tmp;
+	    if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
+		return 0;
+	}
+
+	ltmp = pid_sup;
+	pid_sup = pid_corr;
+	pid_corr = ltmp;
+
+	Rast_set_c_null_value(&precCell, 1);
+
+	connected = 0;
+	for (j = 0; j < ad->cl; j++) {
+	    pid_corr[j] = 0;
+	    
+	    corrCell = buf[j + ad->x];
+	    if (masked && (mask_buf[j] == 0)) {
+		Rast_set_c_null_value(&corrCell, 1);
+	    }
+	    else {
+		/* total sample area */
+		area++;
+	    }
+
+	    supCell = buf_sup[j + ad->x];
+	    if (masked && (mask_sup[j] == 0)) {
+		Rast_set_c_null_value(&supCell, 1);
+	    }
+
+	    if (Rast_is_c_null_value(&corrCell)) {
+		if (!Rast_is_c_null_value(&precCell))
+		    pst[pid_corr[j - 1]].edges++;
+		if (!Rast_is_c_null_value(&supCell))
+		    pst[pid_sup[j]].edges++;
+		connected = 0;
+		precCell = corrCell;
+		continue;
+	    }
+
+	    if (!Rast_is_c_null_value(&precCell) && corrCell == precCell) {
+		pid_corr[j] = pid_corr[j - 1];
+		connected = 1;
+		pst[pid_corr[j]].cells++;
+	    }
+	    else {
+		connected = 0;
+	    }
+
+	    if (!Rast_is_c_null_value(&supCell) && corrCell == supCell) {
+
+		if (pid_corr[j] != pid_sup[j]) {
+		    /* connect or merge */
+		    /* after r.clump */
+		    if (connected) {
+			npatch--;
 
+			if (npatch == 0) {
+			    G_fatal_error("npatch == 0 at row %d, col %d", i, j);
+			}
+		    }
+
+		    old_pid = pid_corr[j];
+		    new_pid = pid_sup[j];
+		    pid_corr[j] = new_pid;
+		    if (old_pid > 0) {
+			/* merge */
+			/* update left side of the current row */
+			for (k = 0; k < j; k++) {
+			    if (pid_corr[k] == old_pid)
+				pid_corr[k] = new_pid;
+			}
+			/* update right side of the previous row */
+			for (k = j + 1; k < ad->cl; k++) {
+			    if (pid_sup[k] == old_pid)
+				pid_sup[k] = new_pid;
+			}
+			pst[new_pid].cells += pst[old_pid].cells;
+			pst[old_pid].cells = 0;
+			pst[new_pid].edges += pst[old_pid].edges;
+			pst[old_pid].edges = 0;
+			
+			if (old_pid == pid)
+			    pid--;
+		    }
+		    else {
+			pst[new_pid].cells++;
+		    }
+		}
+		connected = 1;
+	    }
+
+	    if (!connected) {
+		/* start new patch */
+		npatch++;
+		pid++;
+		pid_corr[j] = pid;
+
+		if (pid >= nalloc) {
+		    pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
+
+		    for (k = nalloc; k < pid + incr; k++) {
+			pst[k].cells = 0;
+			pst[k].edges = 0;
+		    }
+			
+		    nalloc = pid + incr;
+		}
+
+		pst[pid].cells = 1;
+		pst[pid].type.t = CELL_TYPE;
+		pst[pid].type.val.c = corrCell;
+	    }
+	    /* update edge count for corr */
+	    if (Rast_is_c_null_value(&precCell) || precCell != corrCell)
+		pst[pid_corr[j]].edges++;
+	    if (Rast_is_c_null_value(&supCell) || supCell != corrCell)
+		pst[pid_corr[j]].edges++;
+	    if (i == ad->rl - 1)
+		pst[pid_corr[j]].edges++;
+	    if (j == ad->cl - 1)
+		pst[pid_corr[j]].edges++;
+	    /* update edge count for prec */
+	    if (!Rast_is_c_null_value(&precCell) && precCell != corrCell)
+		pst[pid_corr[j - 1]].edges++;
+	    /* update edge count for sup */
+	    if (!Rast_is_c_null_value(&supCell) && supCell != corrCell)
+		pst[pid_sup[j]].edges++;
+
+	    precCell = corrCell;
+	}
+    }
+
+    if (npatch > 0) {
+	double edges, cells;
+
+	edges = cells = 0;
+	
+	for (i = 1; i <= pid; i++) {
+	    cells += pst[i].cells;
+	    edges += pst[i].edges;
+	}
+
+	*result = 0.25 * edges / sqrt(cells);
+    }
+    else {
+	Rast_set_d_null_value(result, 1);
+    }
+
+    if (masked) {
+	close(mask_fd);
+	G_free(mask_buf);
+	G_free(mask_sup);
+    }
+    G_free(buf_null);
+    G_free(pid_corr);
+    G_free(pid_sup);
+    G_free(pst);
+
+    return RLI_OK;
+}
+
+
+int calculateD(int fd, struct area_entry *ad, double *result)
+{
+    DCELL *buf, *buf_sup, *buf_null;
+    DCELL corrCell, precCell, supCell;
+    long npatch, area; 
+    long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
+    struct pst *pst;
+    long nalloc, incr;
+    int i, j, k;
+    int connected;
+    int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
+
+    buf_null = Rast_allocate_d_buf();
+    Rast_set_d_null_value(buf_null, Rast_window_cols());
+    buf_sup = buf_null;
+
+    /* initialize patch ids */
+    pid_corr = G_malloc(ad->cl * sizeof(long));
+    pid_sup = G_malloc(ad->cl * sizeof(long));
+
+    for (j = 0; j < ad->cl; j++) {
+	pid_corr[j] = 0;
+	pid_sup[j] = 0;
+    }
+
+    /* open mask if needed */
+    mask_fd = -1;
+    mask_buf = mask_sup = NULL;
+    masked = FALSE;
+    if (ad->mask == 1) {
+	if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
+	    return RLI_ERRORE;
 	mask_buf = G_malloc(ad->cl * sizeof(int));
+	if (mask_buf == NULL) {
+	    G_fatal_error("malloc mask_buf failed");
+	    return RLI_ERRORE;
+	}
+	mask_sup = G_malloc(ad->cl * sizeof(int));
+	if (mask_sup == NULL) {
+	    G_fatal_error("malloc mask_buf failed");
+	    return RLI_ERRORE;
+	}
+	for (j = 0; j < ad->cl; j++)
+	    mask_buf[j] = 0;
+
+	masked = TRUE;
+    }
+
+    /* calculate number of patches */
+    npatch = 0;
+    area = 0;
+    pid = 0;
 
-	for (i = 0; i < ad->rl; i++) {
+    /* patch size and type */
+    incr = 1024;
+    if (incr > ad->rl)
+	incr = ad->rl;
+    if (incr > ad->cl)
+	incr = ad->cl;
+    if (incr < 2)
+	incr = 2;
+    nalloc = incr;
+    pst = G_malloc(nalloc * sizeof(struct pst));
+    for (k = 0; k < nalloc; k++) {
+	pst[k].cells = 0;
+	pst[k].edges = 0;
+    }
+
+    for (i = 0; i < ad->rl; i++) {
+	buf = RLI_get_dcell_raster_row(fd, i + ad->y, ad);
+	if (i > 0) {
+	    buf_sup = RLI_get_dcell_raster_row(fd, i - 1 + ad->y, ad);
+	}
+
+	if (masked) {
+	    mask_tmp = mask_sup;
+	    mask_sup = mask_buf;
+	    mask_buf = mask_tmp;
 	    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) {
-		    null_count++;
+	ltmp = pid_sup;
+	pid_sup = pid_corr;
+	pid_corr = ltmp;
+
+	Rast_set_d_null_value(&precCell, 1);
+
+	connected = 0;
+	for (j = 0; j < ad->cl; j++) {
+	    pid_corr[j] = 0;
+	    
+	    corrCell = buf[j + ad->x];
+	    if (masked && (mask_buf[j] == 0)) {
+		Rast_set_d_null_value(&corrCell, 1);
+	    }
+	    else {
+		/* total sample area */
+		area++;
+	    }
+
+	    supCell = buf_sup[j + ad->x];
+	    if (masked && (mask_sup[j] == 0)) {
+		Rast_set_d_null_value(&supCell, 1);
+	    }
+
+	    if (Rast_is_d_null_value(&corrCell)) {
+		if (!Rast_is_d_null_value(&precCell))
+		    pst[pid_corr[j - 1]].edges++;
+		if (!Rast_is_d_null_value(&supCell))
+		    pst[pid_sup[j]].edges++;
+		connected = 0;
+		precCell = corrCell;
+		continue;
+	    }
+
+	    if (!Rast_is_d_null_value(&precCell) && corrCell == precCell) {
+		pid_corr[j] = pid_corr[j - 1];
+		connected = 1;
+		pst[pid_corr[j]].cells++;
+	    }
+	    else {
+		connected = 0;
+	    }
+
+	    if (!Rast_is_d_null_value(&supCell) && corrCell == supCell) {
+
+		if (pid_corr[j] != pid_sup[j]) {
+		    /* connect or merge */
+		    /* after r.clump */
+		    if (connected) {
+			npatch--;
+
+			if (npatch == 0) {
+			    G_fatal_error("npatch == 0 at row %d, col %d", i, j);
+			}
+		    }
+
+		    old_pid = pid_corr[j];
+		    new_pid = pid_sup[j];
+		    pid_corr[j] = new_pid;
+		    if (old_pid > 0) {
+			/* merge */
+			/* update left side of the current row */
+			for (k = 0; k < j; k++) {
+			    if (pid_corr[k] == old_pid)
+				pid_corr[k] = new_pid;
+			}
+			/* update right side of the previous row */
+			for (k = j + 1; k < ad->cl; k++) {
+			    if (pid_sup[k] == old_pid)
+				pid_sup[k] = new_pid;
+			}
+			pst[new_pid].cells += pst[old_pid].cells;
+			pst[old_pid].cells = 0;
+			pst[new_pid].edges += pst[old_pid].edges;
+			pst[old_pid].edges = 0;
+			
+			if (old_pid == pid)
+			    pid--;
+		    }
+		    else {
+			pst[new_pid].cells++;
+		    }
 		}
+		connected = 1;
 	    }
+
+	    if (!connected) {
+		/* start new patch */
+		npatch++;
+		pid++;
+		pid_corr[j] = pid;
+
+		if (pid >= nalloc) {
+		    pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
+
+		    for (k = nalloc; k < pid + incr; k++) {
+			pst[k].cells = 0;
+			pst[k].edges = 0;
+		    }
+			
+		    nalloc = pid + incr;
+		}
+
+		pst[pid].cells = 1;
+		pst[pid].type.t = CELL_TYPE;
+		pst[pid].type.val.c = corrCell;
+	    }
+	    /* update edge count for corr */
+	    if (Rast_is_d_null_value(&precCell) || precCell != corrCell)
+		pst[pid_corr[j]].edges++;
+	    if (Rast_is_d_null_value(&supCell) || supCell != corrCell)
+		pst[pid_corr[j]].edges++;
+	    if (i == ad->rl - 1)
+		pst[pid_corr[j]].edges++;
+	    if (j == ad->cl - 1)
+		pst[pid_corr[j]].edges++;
+	    /* update edge count for prec */
+	    if (!Rast_is_d_null_value(&precCell) && precCell != corrCell)
+		pst[pid_corr[j - 1]].edges++;
+	    /* update edge count for sup */
+	    if (!Rast_is_d_null_value(&supCell) && supCell != corrCell)
+		pst[pid_sup[j]].edges++;
+
+	    precCell = corrCell;
 	}
+    }
+
+    if (npatch > 0) {
+	double edges, cells;
+
+	edges = cells = 0;
+	
+	for (i = 1; i <= pid; i++) {
+	    cells += pst[i].cells;
+	    edges += pst[i].edges;
+	}
+
+	*result = 0.25 * edges / sqrt(cells);
+    }
+    else {
+	Rast_set_d_null_value(result, 1);
+    }
+
+    if (masked) {
+	close(mask_fd);
 	G_free(mask_buf);
+	G_free(mask_sup);
     }
+    G_free(buf_null);
+    G_free(pid_corr);
+    G_free(pid_sup);
+    G_free(pst);
+
+    return RLI_OK;
+}
+
+
+int calculateF(int fd, struct area_entry *ad, double *result)
+{
+    FCELL *buf, *buf_sup, *buf_null;
+    FCELL corrCell, precCell, supCell;
+    long npatch, area; 
+    long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
+    struct pst *pst;
+    long nalloc, incr;
+    int i, j, k;
+    int connected;
+    int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
+
+    buf_null = Rast_allocate_f_buf();
+    Rast_set_f_null_value(buf_null, Rast_window_cols());
+    buf_sup = buf_null;
 
-    /* calculate distance */
-    G_begin_distance_calculations();
-    /* EW Dist at North edge */
-    EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
-    /* EW Dist at South Edge */
-    EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
-    /* NS Dist at East edge */
-    NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
-    /* NS Dist at West edge */
-    NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
+    /* initialize patch ids */
+    pid_corr = G_malloc(ad->cl * sizeof(long));
+    pid_sup = G_malloc(ad->cl * sizeof(long));
 
+    for (j = 0; j < ad->cl; j++) {
+	pid_corr[j] = 0;
+	pid_sup[j] = 0;
+    }
 
-    area = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
-	(((NS_DIST1 + NS_DIST2) / 2) / hd.rows) *
-	(ad->rl * ad->cl - null_count);
+    /* open mask if needed */
+    mask_fd = -1;
+    mask_buf = mask_sup = NULL;
+    masked = FALSE;
+    if (ad->mask == 1) {
+	if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
+	    return RLI_ERRORE;
+	mask_buf = G_malloc(ad->cl * sizeof(int));
+	if (mask_buf == NULL) {
+	    G_fatal_error("malloc mask_buf failed");
+	    return RLI_ERRORE;
+	}
+	mask_sup = G_malloc(ad->cl * sizeof(int));
+	if (mask_sup == NULL) {
+	    G_fatal_error("malloc mask_buf failed");
+	    return RLI_ERRORE;
+	}
+	for (j = 0; j < ad->cl; j++)
+	    mask_buf[j] = 0;
+
+	masked = TRUE;
+    }
 
-    *result = area;
+    /* calculate number of patches */
+    npatch = 0;
+    area = 0;
+    pid = 0;
+
+    /* patch size and type */
+    incr = 1024;
+    if (incr > ad->rl)
+	incr = ad->rl;
+    if (incr > ad->cl)
+	incr = ad->cl;
+    if (incr < 2)
+	incr = 2;
+    nalloc = incr;
+    pst = G_malloc(nalloc * sizeof(struct pst));
+    for (k = 0; k < nalloc; k++) {
+	pst[k].cells = 0;
+	pst[k].edges = 0;
+    }
+
+    for (i = 0; i < ad->rl; i++) {
+	buf = RLI_get_fcell_raster_row(fd, i + ad->y, ad);
+	if (i > 0) {
+	    buf_sup = RLI_get_fcell_raster_row(fd, i - 1 + ad->y, ad);
+	}
+
+	if (masked) {
+	    mask_tmp = mask_sup;
+	    mask_sup = mask_buf;
+	    mask_buf = mask_tmp;
+	    if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
+		return 0;
+	}
+
+	ltmp = pid_sup;
+	pid_sup = pid_corr;
+	pid_corr = ltmp;
+
+	Rast_set_f_null_value(&precCell, 1);
+
+	connected = 0;
+	for (j = 0; j < ad->cl; j++) {
+	    pid_corr[j] = 0;
+	    
+	    corrCell = buf[j + ad->x];
+	    if (masked && (mask_buf[j] == 0)) {
+		Rast_set_f_null_value(&corrCell, 1);
+	    }
+	    else {
+		/* total sample area */
+		area++;
+	    }
+
+	    supCell = buf_sup[j + ad->x];
+	    if (masked && (mask_sup[j] == 0)) {
+		Rast_set_f_null_value(&supCell, 1);
+	    }
+
+	    if (Rast_is_f_null_value(&corrCell)) {
+		if (!Rast_is_f_null_value(&precCell))
+		    pst[pid_corr[j - 1]].edges++;
+		if (!Rast_is_f_null_value(&supCell))
+		    pst[pid_sup[j]].edges++;
+		connected = 0;
+		precCell = corrCell;
+		continue;
+	    }
+
+	    if (!Rast_is_f_null_value(&precCell) && corrCell == precCell) {
+		pid_corr[j] = pid_corr[j - 1];
+		connected = 1;
+		pst[pid_corr[j]].cells++;
+	    }
+	    else {
+		connected = 0;
+	    }
+
+	    if (!Rast_is_f_null_value(&supCell) && corrCell == supCell) {
+
+		if (pid_corr[j] != pid_sup[j]) {
+		    /* connect or merge */
+		    /* after r.clump */
+		    if (connected) {
+			npatch--;
+
+			if (npatch == 0) {
+			    G_fatal_error("npatch == 0 at row %d, col %d", i, j);
+			}
+		    }
+
+		    old_pid = pid_corr[j];
+		    new_pid = pid_sup[j];
+		    pid_corr[j] = new_pid;
+		    if (old_pid > 0) {
+			/* merge */
+			/* update left side of the current row */
+			for (k = 0; k < j; k++) {
+			    if (pid_corr[k] == old_pid)
+				pid_corr[k] = new_pid;
+			}
+			/* update right side of the previous row */
+			for (k = j + 1; k < ad->cl; k++) {
+			    if (pid_sup[k] == old_pid)
+				pid_sup[k] = new_pid;
+			}
+			pst[new_pid].cells += pst[old_pid].cells;
+			pst[old_pid].cells = 0;
+			pst[new_pid].edges += pst[old_pid].edges;
+			pst[old_pid].edges = 0;
+			
+			if (old_pid == pid)
+			    pid--;
+		    }
+		    else {
+			pst[new_pid].cells++;
+		    }
+		}
+		connected = 1;
+	    }
+
+	    if (!connected) {
+		/* start new patch */
+		npatch++;
+		pid++;
+		pid_corr[j] = pid;
+
+		if (pid >= nalloc) {
+		    pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
+
+		    for (k = nalloc; k < pid + incr; k++) {
+			pst[k].cells = 0;
+			pst[k].edges = 0;
+		    }
+			
+		    nalloc = pid + incr;
+		}
+
+		pst[pid].cells = 1;
+		pst[pid].type.t = CELL_TYPE;
+		pst[pid].type.val.c = corrCell;
+	    }
+	    /* update edge count for corr */
+	    if (Rast_is_f_null_value(&precCell) || precCell != corrCell)
+		pst[pid_corr[j]].edges++;
+	    if (Rast_is_f_null_value(&supCell) || supCell != corrCell)
+		pst[pid_corr[j]].edges++;
+	    if (i == ad->rl - 1)
+		pst[pid_corr[j]].edges++;
+	    if (j == ad->cl - 1)
+		pst[pid_corr[j]].edges++;
+	    /* update edge count for prec */
+	    if (!Rast_is_f_null_value(&precCell) && precCell != corrCell)
+		pst[pid_corr[j - 1]].edges++;
+	    /* update edge count for sup */
+	    if (!Rast_is_f_null_value(&supCell) && supCell != corrCell)
+		pst[pid_sup[j]].edges++;
+
+	    precCell = corrCell;
+	}
+    }
+
+    if (npatch > 0) {
+	double edges, cells;
+
+	edges = cells = 0;
+	
+	for (i = 1; i <= pid; i++) {
+	    cells += pst[i].cells;
+	    edges += pst[i].edges;
+	}
+
+	*result = 0.25 * edges / sqrt(cells);
+    }
+    else {
+	Rast_set_d_null_value(result, 1);
+    }
+
+    if (masked) {
+	close(mask_fd);
+	G_free(mask_buf);
+	G_free(mask_sup);
+    }
+    G_free(buf_null);
+    G_free(pid_corr);
+    G_free(pid_sup);
+    G_free(pst);
 
-    return 1;
+    return RLI_OK;
 }

raster/r.li/r.li.shape/r.li.shape.html

@@ -1,7 +1,15 @@
 <h2>DESCRIPTION</h2>
 
-<em>r.li.shape</em> calculates the "shape index" as:<br>
-<I> f(sample_area)= Area </I><br>
+<em>r.li.shape</em> calculates the landscape shape index as:
+<div class="code"><pre>
+LSI = 0.25 * E / sqrt(A)
+</pre></div>
+
+with:
+<ul>
+  <li><b>E</b>: sum of all edges</li>
+  <li><b>A</b>: sum of all patch areas</li>
+</ul>
 
 <h2>NOTES</h2>
 
@@ -12,13 +20,7 @@ the <tt>$HOME/.grass7/r.li/output/</tt> folder.
 <p>
 <!-- TODO: verify next: -->
 If the input raster map contains only NULL values then <em>r.li.shape</em>
-returns -1.<br>
-<!-- does that mean the program exit code or raster values?? -->
-If you want to change these -1 values to NULL, run subsequently on the resulting map:
-<div class="code"><pre>
-r.null setnull=-1 input=my_map
-</pre></div>
-after index calculation.
+returns NULL.
 
 <h2>EXAMPLES</h2>
 To calculate the shape index on map <em>my_map</em>, using