[GRASS-SVN] r42949 - in grass-addons/raster/r.pi: . r.pi.index

svn_grass at osgeo.org svn_grass at osgeo.org
Thu Jul 29 12:21:45 EDT 2010


Author: wegmann
Date: 2010-07-29 16:21:45 +0000 (Thu, 29 Jul 2010)
New Revision: 42949

Added:
   grass-addons/raster/r.pi/r.pi.index/
   grass-addons/raster/r.pi/r.pi.index/Makefile
   grass-addons/raster/r.pi/r.pi.index/description.html
   grass-addons/raster/r.pi/r.pi.index/frag.c
   grass-addons/raster/r.pi/r.pi.index/func.c
   grass-addons/raster/r.pi/r.pi.index/local_proto.h
   grass-addons/raster/r.pi/r.pi.index/main.c
Log:
initial commit: r.pi.index provides basic patch based spatial indices like area, shape, euclidean distance to nearest neighbour etc.

Added: grass-addons/raster/r.pi/r.pi.index/Makefile
===================================================================
--- grass-addons/raster/r.pi/r.pi.index/Makefile	                        (rev 0)
+++ grass-addons/raster/r.pi/r.pi.index/Makefile	2010-07-29 16:21:45 UTC (rev 42949)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.pi.index
+
+LIBES = $(STATSLIB) $(GISLIB)
+DEPENDENCIES = $(STATSDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd


Property changes on: grass-addons/raster/r.pi/r.pi.index/Makefile
___________________________________________________________________
Added: svn:executable
   + *

Added: grass-addons/raster/r.pi/r.pi.index/description.html
===================================================================
--- grass-addons/raster/r.pi/r.pi.index/description.html	                        (rev 0)
+++ grass-addons/raster/r.pi/r.pi.index/description.html	2010-07-29 16:21:45 UTC (rev 42949)
@@ -0,0 +1,40 @@
+<H2>DESCRIPTION</H2>
+
+Available options for the index to be computed for patches within a certain class are: area (area), perimeter (perim), SHAPE (shape), Border-Index (bor), Compactness (comp), Asymmetry (asym), area-perimeter ratio (apr), fractal dimension (fract), distance to euclidean nearest neighbour (ENN)</DD>
+
+<P>
+
+The program will be run non-interactively if the user specifies program arguments (see OPTIONS) on the command
+line.  Alternately, the user can simply type <B>r.pi.index</B> on the command line, without program
+arguments.  In this case, the user will be prompted for flag settings and parameter values.
+
+
+<H2>NOTES</H2>
+
+The <EM>Nearest Neighbour Index</EM> (ENN) analyse the Euclidean Nearest 
+Neighbour to the first neighbouring patch. The output value is in pixel and can
+be converted to a distance values using g.region resolution information.
+<EM>r.pi.ENN</EM> and <EM>r.pi.FNN</EM> provide the same analysis concerning the
+ first nearest neighbour (NN), but are extended to the n-th NN. However due to
+code construction does the <EM>r.pi.index</EM> distance analysis to first ENN perform faster.
+
+
+<P>
+
+<H2>BUGS</H2>
+Landscapes with more than 10 000 patches might cause a memory
+allocation error depending on your system.
+
+
+<H2>SEE ALSO</H2>
+
+<EM><A HREF="r.pi.ENN.html">r.pi.ENN</A></EM><br>
+<EM><A HREF="r.pi.FNN.html">r.pi.FNN</A></EM><br>
+<EM><A HREF="r.pi.dist.html">r.pi.dist</A></EM><br>
+<EM><A HREF="r.li.setup.html">r.li.setup</A></EM><br>
+
+<H2>AUTHOR</H2>
+Programming: Elshad Shirinov<br>
+Scientific concept: Dr. Martin Wegmann <br>
+Department of Remote Sensing <br>Remote Sensing and Biodiversity Unit<br> University of Wuerzburg, Germany
+


Property changes on: grass-addons/raster/r.pi/r.pi.index/description.html
___________________________________________________________________
Added: svn:executable
   + *

Added: grass-addons/raster/r.pi/r.pi.index/frag.c
===================================================================
--- grass-addons/raster/r.pi/r.pi.index/frag.c	                        (rev 0)
+++ grass-addons/raster/r.pi/r.pi.index/frag.c	2010-07-29 16:21:45 UTC (rev 42949)
@@ -0,0 +1,108 @@
+#include <string.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/stats.h>
+#include <math.h>
+#include "local_proto.h"
+
+typedef struct {
+		int x, y;
+	} Position;
+
+int getNeighbors(Position* res, int x, int y, int nx, int ny, int nbr_cnt);	
+	
+int getNeighbors(Position* res, int x, int y, int nx, int ny, int nbr_cnt) {
+	int left, right, top, bottom;
+	int i, j;
+	int cnt = 0;
+	
+	switch(nbr_cnt) {
+		case 4: /* von Neumann neighborhood */
+			if (x > 0 && flagbuf[y * nx + x - 1] == 1) {
+				res[cnt].x = x - 1; res[cnt].y = y; cnt++;				
+			}
+			if (y > 0 && flagbuf[(y - 1) * nx + x] == 1) {
+				res[cnt].x = x; res[cnt].y = y - 1; cnt++;
+			}
+			if (x < nx - 1 && flagbuf[y * nx + x + 1] == 1) {
+				res[cnt].x = x + 1; res[cnt].y = y; cnt++;
+			}
+			if (y < ny - 1 && flagbuf[(y + 1) * nx + x] == 1) {
+				res[cnt].x = x; res[cnt].y = y + 1; cnt++;
+			}
+			break;
+		case 8: /* Moore neighborhood */
+			left = x > 0 ? x - 1 : 0;
+			top = y > 0 ? y - 1 : 0;
+			right = x < nx - 1 ? x + 1 : nx - 1;
+			bottom = y < ny - 1 ? y + 1 : ny - 1;			
+			for (i = left; i <= right; i++) {
+				for (j = top; j <= bottom; j++) {
+					if(!(i == x && j == y) && flagbuf[j * nx + i] == 1) {
+						res[cnt].x = i; res[cnt].y = j; cnt++;
+					}
+				}
+			}
+			break;
+	}
+	
+	return cnt;
+}	
+	
+void writeFrag(int row, int col, int nbr_cnt) {
+	int x, y, i;
+	Position* list = (Position *) G_malloc(nrows * ncols * sizeof (Position));
+	Position* first = list;
+	Position* last = list;
+	Position* nbr_list = (Position *) G_malloc(8 * sizeof (Position));
+	
+	/* count neighbors */
+	int neighbors = 0;
+	if (col > 0 && flagbuf[row * ncols + col - 1] != 0) neighbors++;
+	if (row > 0 && flagbuf[(row - 1) * ncols + col] != 0) neighbors++;
+	if (col < ncols - 1 && flagbuf[row * ncols + col + 1] != 0) neighbors++;
+	if (row < nrows - 1 && flagbuf[(row + 1) * ncols + col] != 0) neighbors++;
+	
+	/* write first cell */
+	actpos->x = col; actpos->y = row; actpos->neighbors = neighbors; actpos++;
+	flagbuf[row * ncols + col] = -1;
+
+	/* push position on fifo-list */
+	last->x = col; last->y = row; last++;	
+
+	while (first < last) {
+		/* get position from fifo-list */
+		int r = first->y; int c = first->x; first++;
+		
+		int left = c > 0 ? c - 1 : 0;
+		int top = r > 0 ? r - 1 : 0;
+		int right = c < ncols - 1 ? c + 1 : ncols - 1;
+		int bottom = r < nrows - 1 ? r + 1 : nrows - 1;
+		
+		/* add neighbors to fifo-list */
+		int cnt = getNeighbors(nbr_list, c, r, ncols, nrows, nbr_cnt);
+				
+		for(i = 0; i < cnt; i++) {
+			x = nbr_list[i].x; y = nbr_list[i].y;
+			
+			/* add position to fifo-list */
+			last->x = x; last->y = y; last++;
+				
+			/* count neighbors */
+			neighbors = 0;
+			if (x > 0 && flagbuf[y * ncols + x - 1] != 0) neighbors++;
+			if (y > 0 && flagbuf[(y - 1) * ncols + x] != 0) neighbors++;
+			if (x < ncols - 1 && flagbuf[y * ncols + x + 1] != 0) neighbors++;
+			if (y < nrows - 1 && flagbuf[(y + 1) * ncols + x] != 0) neighbors++;
+	
+			/* set values */
+			actpos->x = x; actpos->y = y; actpos->neighbors = neighbors; actpos++;
+			flagbuf[y * ncols + x] = -1;		
+		}
+	}
+	
+	G_free(list);
+	G_free(nbr_list);
+}

Added: grass-addons/raster/r.pi/r.pi.index/func.c
===================================================================
--- grass-addons/raster/r.pi/r.pi.index/func.c	                        (rev 0)
+++ grass-addons/raster/r.pi/r.pi.index/func.c	2010-07-29 16:21:45 UTC (rev 42949)
@@ -0,0 +1,264 @@
+#include <string.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/stats.h>
+#include <math.h>
+#include "local_proto.h"
+
+int f_area(DCELL *vals, Coords **frags, int count) {
+	int i;
+	/* for all patches */
+	for(i = 0; i < count; i++) {
+		// display progress
+		if (verbose)
+			G_percent (i, count, 2);
+		
+		vals[i] = (DCELL) (frags[i+1] - frags[i]);
+	}
+
+	return 0;
+}
+
+int f_perim(DCELL *vals, Coords **frags, int count) {
+	Coords *p;
+	int i;
+	/* for all patches */
+	for(i = 0; i < count; i++) {
+		// display progress
+		if (verbose)
+			G_percent (i, count, 2);
+		
+		int border = 0;
+		/* for all cells in a patch */
+		for(p = frags[i]; p < frags[i+1]; p++) {
+			border+= 4 - p->neighbors;
+		}
+		vals[i] = (DCELL) border;
+	}
+	
+	return 0;
+}
+
+int f_shapeindex(DCELL *vals, Coords **frags, int count) {
+	Coords *p;
+	int i;
+	/* for all patches */
+	for(i = 0; i < count; i++) {
+		// display progress
+		if (verbose)
+			G_percent (i, count, 2);
+		
+		int border = 0;
+		int area = (DCELL) (frags[i+1] - frags[i]);
+		/* for all cells in a patch */
+		for(p = frags[i]; p < frags[i+1]; p++) {
+			border+= 4 - p->neighbors;
+		}
+		vals[i] = (DCELL) border / ( 4 * sqrt((DCELL) area) );
+	}
+	
+	return 0;
+}
+
+int f_borderindex(DCELL *vals, Coords **frags, int count) {
+	Coords *p;
+	int i;
+	/* for all patches */
+	for(i = 0; i < count; i++) {
+		// display progress
+		if (verbose)
+			G_percent (i, count, 2);
+		
+		int border = 0;
+		int maxx, maxy, minx, miny;
+		int l = 0;
+		int w = 0;
+		maxx = minx = frags[i]->x;
+		maxy = miny = frags[i]->y;
+		/* for all cells in a patch */
+		for(p = frags[i]; p < frags[i+1]; p++) {
+			border+= 4 - p->neighbors;
+			maxx = p->x > maxx ? p->x : maxx;
+			minx = p->x < minx ? p->x : minx;
+			maxy = p->y > maxy ? p->y : maxy;
+			miny = p->y < miny ? p->y : miny;
+		}
+		l = maxx - minx + 1; w = maxy - miny + 1;
+		vals[i] = (DCELL) border / ( 2 * ( (DCELL)l + (DCELL)w ) );
+	}
+	
+	return 0;
+}
+
+int f_compactness(DCELL *vals, Coords **frags, int count) {
+	Coords *p;
+	int i;
+	/* for all patches */
+	for(i = 0; i < count; i++) {
+		// display progress
+		if (verbose)
+			G_percent (i, count, 2);
+		
+		int area = 0;
+		int maxx, maxy, minx, miny;
+		int l = 0;
+		int w = 0;
+		maxx = minx = frags[i]->x;
+		maxy = miny = frags[i]->y;
+		/* for all cells in a patch */
+		for(p = frags[i]; p < frags[i+1]; p++) {
+			area++;
+			maxx = p->x > maxx ? p->x : maxx;
+			minx = p->x < minx ? p->x : minx;
+			maxy = p->y > maxy ? p->y : maxy;
+			miny = p->y < miny ? p->y : miny;
+		}
+		l = maxx - minx + 1; w = maxy - miny + 1;
+		vals[i] = (DCELL)l * (DCELL)w / (DCELL) area;
+	}
+	
+	return 0;
+}
+
+int f_asymmetry(DCELL *vals, Coords **frags, int count) {
+	Coords *p;
+	int i;
+	/* for all patches */
+	for(i = 0; i < count; i++) {
+		// display progress
+		if (verbose)
+			G_percent (i, count, 2);
+		
+		/* compute variance for x,y and xy in the patch */
+		/* formula: a(x) = sum(x_i), b(x) = sum(x_i²), var(x) = (b(x) - a(x)² / n) / n */
+		/* covar(x*y) = (a(x * y) - a(x) * a(y) / n) /n */
+		int ax, ay, axy;
+		int bx, by;
+		DCELL vx, vy, vxy, vsum, invn;
+		int n = 0;
+		ax = ay = axy = 0; bx = by = 0;
+		/* for all cells in a patch */
+//		fprintf(stderr, "\npatch %d: ", i);
+		for(p = frags[i]; p < frags[i+1]; p++, n++) {
+			int x = p->x; int y = p->y; int xy = p->x * p->y;
+			ax += x; ay += y; axy += xy;
+			bx += x * x; by += y * y;
+//			fprintf(stderr, "x_%d = %d, y_%d = %d; ", n, x, n, y);
+		}
+		invn = 1.0 / (DCELL)n;
+		vx = ( (DCELL)bx - (DCELL)ax * (DCELL)ax * invn ) * invn;
+		vy = ( (DCELL)by - (DCELL)ay * (DCELL)ay * invn ) * invn;
+		vxy = ( (DCELL)axy - (DCELL)ax * (DCELL)ay * invn ) * invn;
+//		fprintf(stderr, " axy = %d, ax = %d, ay = %d, n = %d", axy, ax, ay, n);
+		vsum = vx + vy;
+		vals[i] = 2 * sqrt(0.25 * vsum * vsum + vxy * vxy - vx * vy) / vsum;
+	}
+}
+
+int f_area_perim_ratio(DCELL *vals, Coords **frags, int count) {
+	Coords *p;
+	int i;
+	/* for all patches */
+	for(i = 0; i < count; i++) {
+		// display progress
+		if (verbose)
+			G_percent (i, count, 2);
+		
+		int border = 0;
+		int area = (DCELL) (frags[i+1] - frags[i]);
+		/* for all cells in a patch */
+		for(p = frags[i]; p < frags[i+1]; p++) {
+			border+= 4 - p->neighbors;
+		}
+		vals[i] = (DCELL) area / (DCELL) border;
+	}
+	
+	return 0;
+}
+
+int f_frac_dim(DCELL *vals, Coords **frags, int count) {
+	Coords *p;
+	int i;
+	/* for all patches */
+	for(i = 0; i < count; i++) {
+		// display progress
+		if (verbose)
+			G_percent (i, count, 2);
+		
+		int border = 0;
+		int area = (DCELL) (frags[i+1] - frags[i]);
+		/* for all cells in a patch */
+		for(p = frags[i]; p < frags[i+1]; p++) {
+			border+= 4 - p->neighbors;
+		}
+		vals[i] = 2 * log(0.25 * (DCELL) border) / log((DCELL) area);
+	}
+	
+	return 0;
+}
+
+DCELL dist(Coords *p1, Coords *p2) {
+	int x1 = p1->x; int y1 = p1->y;
+	int x2 = p2->x; int y2 = p2->y;
+	int dx = x2 - x1; int dy = y2 - y1;
+	return sqrt(dx * dx + dy * dy);
+}
+
+DCELL min_dist(Coords **frags, int n1, int n2) {
+	Coords *p1, *p2;
+	DCELL min = 1000000.0;
+	// for all cells in the first patch
+	for (p1 = frags[n1]; p1 < frags[n1+1]; p1++) {	
+		// if cell at the border
+		if (p1->neighbors < 4) {
+			// for all cells in the second patch
+			for (p2 = frags[n2]; p2 < frags[n2+1]; p2++) {
+				// if cell at the border
+				if (p2->neighbors < 4) {
+					DCELL d = dist(p1, p2);
+					if (d < min) {
+						min = d;
+					}
+				}
+			}
+		}
+	}
+	return min;
+}
+
+int f_nearest_dist(DCELL *vals, Coords **frags, int count) {
+	int i, j;
+	/* for all patches */
+	for (i = 0; i < count; i++) {
+		// display progress
+		if (verbose)
+			G_percent (i, count, 2);
+		
+		DCELL min = 1000000.0;
+		for (j = 0; j < count; j++) {
+			if (i != j) {
+				DCELL d = min_dist(frags, i, j);
+				if (d < min) {
+					min = d;
+				}
+			}
+		}
+		vals[i] = min;
+	}
+
+	return 0;
+}
+
+int f_proximity(DCELL *vals, Coords **frags, int count, DCELL param) {
+	int i, j;
+	for (i = 0; i < count; i++) {
+		for (j = 0; j < count; j++) {
+			if (i != j) {
+			}
+		}
+	}
+
+	return 0;
+}

Added: grass-addons/raster/r.pi/r.pi.index/local_proto.h
===================================================================
--- grass-addons/raster/r.pi/r.pi.index/local_proto.h	                        (rev 0)
+++ grass-addons/raster/r.pi/r.pi.index/local_proto.h	2010-07-29 16:21:45 UTC (rev 42949)
@@ -0,0 +1,30 @@
+#ifndef LOCAL_PROTO_H
+#define LOCAL_PROTO_H
+
+#define MAX_DOUBLE 1000000.0
+
+typedef struct {
+		int x, y;
+		int neighbors;
+	} Coords;
+
+extern void writeFrag(int row, int col, int nbr_cnt);
+
+extern int f_area(DCELL *vals, Coords **frags, int);
+extern int f_perim(DCELL *vals, Coords **frags, int);
+extern int f_shapeindex(DCELL*, Coords**, int);
+extern int f_borderindex(DCELL*, Coords**, int);
+extern int f_compactness(DCELL*, Coords**, int);
+extern int f_asymmetry(DCELL*, Coords**, int);
+extern int f_area_perim_ratio(DCELL*, Coords**, int);
+extern int f_frac_dim(DCELL*, Coords**, int);
+extern int f_nearest_dist(DCELL*, Coords**, int);
+
+/* global variables */
+int nrows, ncols;
+Coords ** fragments;
+int *flagbuf;
+Coords * actpos;
+int verbose;
+
+#endif /* LOCAL_PROTO_H */

Added: grass-addons/raster/r.pi/r.pi.index/main.c
===================================================================
--- grass-addons/raster/r.pi/r.pi.index/main.c	                        (rev 0)
+++ grass-addons/raster/r.pi/r.pi.index/main.c	2010-07-29 16:21:45 UTC (rev 42949)
@@ -0,0 +1,306 @@
+#include <string.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/stats.h>
+#include <math.h>
+#include "local_proto.h"
+
+/*
+	Fragmentation analysis - basic indices in r.pi.index
+	
+	by Elshad Shirinov.
+*/
+
+typedef int (*f_func)(DCELL*, Coords**, int);
+
+struct menu
+{
+	f_func method;		/* routine to compute new value */
+	char *name;		/* method name */
+	char *text;		/* menu display - full description */
+};
+
+static struct menu menu[] =
+{
+	{f_area, "area", "area of specified patches"},
+	{f_perim, "perimeter", "perimeter of specified patches"},
+	{f_shapeindex, "shape", "shape index of specified patches"},
+	{f_borderindex, "border", "border index of specified patches"},
+	{f_compactness, "compactness", "compactness of specified patches"},
+	{f_asymmetry, "asymmetry", "asymmetry of specified patches"},
+	{f_area_perim_ratio, "area-perimeter", "area-perimeter ratio of specified patches"},
+	{f_frac_dim, "fractal", "fractal dimension of specified patches"},
+	{f_nearest_dist, "ENN", "euclidean distance to nearest patch"},
+	{0,0,0}
+};
+
+int main (int argc, char *argv[])
+{
+	/* input */
+	char *newname, *oldname, *newmapset, *oldmapset;
+	char title[1024];
+
+	/* in and out file pointers */
+	int in_fd;
+	int out_fd;
+	DCELL *result, res[30];
+	
+	/* map_type and categories */
+	RASTER_MAP_TYPE map_type;
+	struct Categories cats;
+	
+	int method;
+	f_func compute_values;
+	
+	/* neighbors count */
+	int neighb_count;
+
+	char *p;
+
+	int row, col, i, j;
+	int readrow;
+	int keyval;
+	DCELL range = MAX_DOUBLE;
+	
+	int n;
+	int copycolr;
+	struct Colors colr;
+	struct GModule *module;
+	struct
+	{
+		struct Option *input, *output;
+		struct Option *keyval, *method;
+		struct Option *title;
+	} parm;
+	struct
+	{
+		struct Flag *adjacent, *quiet;
+	} flag;
+
+	DCELL *values;
+	Coords *cells;
+	int fragcount = 0;
+	
+	struct Cell_head ch, window;
+
+	G_gisinit (argv[0]);
+
+	module = G_define_module();
+	module->keywords = _("raster");
+	module->description =
+		_("Provides basic patch based indices, like area, SHAPE or Nearest Neighbour distance.");
+
+	parm.input = G_define_option() ;
+	parm.input->key        = "input" ;
+	parm.input->type       = TYPE_STRING ;
+	parm.input->required   = YES ;
+	parm.input->gisprompt  = "old,cell,raster" ;
+	parm.input->description= _("Raster map containing categories") ;
+
+	parm.output = G_define_option() ;
+	parm.output->key        = "output" ;
+	parm.output->type       = TYPE_STRING ;
+	parm.output->required   = YES ;
+	parm.output->gisprompt  = "new,cell,raster" ;
+	parm.output->description= _("Output patch-based result as raster map") ;
+
+	parm.keyval = G_define_option() ;
+	parm.keyval->key        = "keyval" ;
+	parm.keyval->type       = TYPE_INTEGER ;
+	parm.keyval->required   = YES ;
+	parm.keyval->description= _("The value of the class to be analysed") ;
+
+	parm.method = G_define_option() ;
+	parm.method->key        = "method" ;
+	parm.method->type       = TYPE_STRING ;
+	parm.method->required   = YES ;
+	p = parm.method->options  = G_malloc(1024);
+	for (n = 0; menu[n].name; n++)
+	{
+		if (n)
+			strcat (p, ",");
+		else
+			*p = 0;
+		strcat (p, menu[n].name);
+	}
+	parm.method->description= _("Operation to perform on fragments") ;
+
+	parm.title = G_define_option() ;
+	parm.title->key        = "title" ;
+	parm.title->key_desc   = "\"phrase\"" ;
+	parm.title->type       = TYPE_STRING ;
+	parm.title->required   = NO ;
+	parm.title->description= _("Title of the output raster file") ;
+
+	flag.adjacent = G_define_flag() ;
+	flag.adjacent->key        = 'a' ;
+	flag.adjacent->description= _("Set for 8 cell-neighbors. 4 cell-neighbors are default.") ;
+	
+	flag.quiet = G_define_flag();
+	flag.quiet->key = 'q';
+	flag.quiet->description = _("Run quietly");
+
+	if (G_parser(argc,argv))
+		exit(1);
+
+	/* get names of input files */
+	oldname = parm.input->answer;
+
+	/* test input files existance */
+	if(NULL == (oldmapset = G_find_cell2(oldname,"")))
+	{
+		fprintf (stderr, "%s: <%s> raster file not found\n",
+			 G_program_name(), oldname);
+		exit(1);
+	}
+
+	/* check if the new file name is correct */
+	newname = parm.output->answer;
+	if (G_legal_filename(newname) < 0)
+	{
+		fprintf (stderr, "%s: <%s> illegal file name\n",
+			 G_program_name(), newname);
+		exit(1);
+	}
+	newmapset = G_mapset();
+
+	/* get size */
+	nrows = G_window_rows();
+	ncols = G_window_cols();
+	
+	/* open cell files */
+	if ((in_fd = G_open_cell_old (oldname, oldmapset)) < 0)
+	{
+		char msg[200];
+		sprintf(msg,"can't open cell file <%s> in mapset %s\n",
+			oldname, oldmapset);
+		G_fatal_error (msg);
+		exit(-1);
+	}
+
+	/* get map type */
+	map_type = DCELL_TYPE;//G_raster_map_type(oldname, oldmapset);
+
+	/* copy color table */
+	copycolr = (G_read_colors (oldname, oldmapset, &colr) > 0);
+	
+	/* get key value */
+	sscanf(parm.keyval->answer, "%d", &keyval);
+
+	/* get the method */
+	for (method = 0; (p = menu[method].name); method++)
+		if ((strcmp(p, parm.method->answer) == 0))
+			break;
+	if (!p)
+	{
+		G_warning (_("<%s=%s> unknown %s"),
+			 parm.method->key, parm.method->answer, parm.method->key);
+		G_usage();
+		exit(EXIT_FAILURE);
+	}
+
+	/* establish the newvalue routine */
+	compute_values = menu[method].method;
+	
+	/* get number of cell-neighbors */
+	neighb_count = flag.adjacent->answer ? 8 : 4;
+		
+	/* allocate the cell buffers */
+	cells = (Coords *) G_malloc ( nrows * ncols * sizeof(Coords));
+	actpos = cells;
+	fragments = (Coords **) G_malloc (nrows * ncols * sizeof (Coords*));
+	fragments[0] = cells;
+	flagbuf = (int *) G_malloc (nrows * ncols * sizeof (int));
+	result = G_allocate_d_raster_buf();
+
+	/* get title, initialize the category and stat info */
+	if (parm.title->answer)
+		strcpy (title, parm.title->answer);
+	else
+		sprintf (title,"Fragmentation of file: %s", oldname);
+
+	/* open the new cellfile  */
+	out_fd = G_open_raster_new (newname, map_type);
+	if (out_fd < 0) {
+		char msg[200];
+		sprintf(msg,"can't create new cell file <%s> in mapset %s\n",
+			newname, newmapset);
+		G_fatal_error (msg);
+		exit(1);
+	}
+
+	if (verbose = !flag.quiet->answer)
+		fprintf (stderr, "Loading Patches ... ");
+
+	/* find fragments */
+	for (row = 0; row < nrows; row++) {
+		G_get_d_raster_row (in_fd, result, row);
+		for (col = 0; col < ncols; col++) {
+			if(result[col] == keyval)
+				flagbuf[row * ncols + col] = 1;
+		}
+		
+		if (verbose)
+			G_percent (row, nrows, 2);
+	}
+
+	for(row = 0; row < nrows; row++) {
+		for(col = 0; col < ncols; col++) {
+			if(flagbuf[row * ncols + col] == 1) {
+				fragcount++;
+				writeFrag(row, col, neighb_count);
+				fragments[fragcount] = actpos;
+			}
+		}
+	}
+	if (verbose)
+		G_percent (nrows, nrows, 2);	
+
+	/* perform actual function on the patches */
+	if (verbose = !flag.quiet->answer)
+		fprintf (stderr, "Performing operation ... ");
+	values = (DCELL *) G_malloc ( fragcount * sizeof(DCELL));
+	compute_values(values, fragments, fragcount);
+	if (verbose)
+		G_percent (fragcount, fragcount, 2);
+
+	/* write the output file */
+	if (verbose = !flag.quiet->answer)
+		fprintf (stderr, "Writing output ... ");
+	for(row = 0; row < nrows; row++) {
+		G_set_d_null_value(result, ncols);
+	
+		for(i = 0; i < fragcount; i++) {
+			for(actpos = fragments[i]; actpos < fragments[i+1]; actpos++) {
+				if(actpos->y == row) {
+					result[actpos->x] = values[i];
+				}
+			}
+		}
+		
+		G_put_d_raster_row(out_fd, result);
+		
+		if (verbose)
+			G_percent (row, nrows, 2);
+	}
+	
+	if (verbose)
+		G_percent (nrows, nrows, 2);
+		
+	G_close_cell (out_fd);
+	G_close_cell (in_fd);
+
+	G_free(cells);
+	G_free(fragments);
+	G_free(flagbuf);
+
+	G_init_cats (0, title, &cats);
+	G_write_cats (newname, &cats);
+
+	if(copycolr)
+		G_write_colors (newname, newmapset, &colr);
+
+	exit(0);
+}


Property changes on: grass-addons/raster/r.pi/r.pi.index/main.c
___________________________________________________________________
Added: svn:executable
   + *



More information about the grass-commit mailing list