r.li.daemon: fix sampling frame and sample area logic, consider current region

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

raster/r.li/r.li.daemon/daemon.c

@@ -195,6 +195,7 @@ int parseSetup(char *path, struct list *l, struct g_area *g, char *raster)
     double rel_x, rel_y, rel_rl, rel_cl;
     double sf_n, sf_s, sf_e, sf_w;
     int sf_x, sf_y, sf_rl, sf_cl;
+    int sa_x, sa_y, sa_rl, sa_cl;
     int size;
 
     if (stat(path, &s) != 0)
@@ -220,14 +221,28 @@ int parseSetup(char *path, struct list *l, struct g_area *g, char *raster)
     rel_rl = atof(strtok(NULL, "|"));
     rel_cl = atof(strtok(NULL, "\n"));
 
-    /* find raster map */
-    Rast_get_cellhd(raster, "", &cellhd);
+    /* use current region ! */
+    Rast_get_window(&cellhd);
 
     /* calculate absolute sampling frame definition */
     sf_x = (int)rint(cellhd.cols * rel_x);
     sf_y = (int)rint(cellhd.rows * rel_y);
     sf_rl = (int)rint(cellhd.rows * rel_rl);
     sf_cl = (int)rint(cellhd.cols * rel_cl);
+    
+    /* sanity check */
+    if (sf_x < 0)
+	sf_x = 0;
+    if (sf_y < 0)
+	sf_y = 0;
+    if (sf_x > cellhd.cols)
+	sf_x = cellhd.cols;
+    if (sf_y > cellhd.rows)
+	sf_y = cellhd.rows;
+    if (sf_rl > cellhd.rows - sf_y)
+	sf_rl = cellhd.rows - sf_y;
+    if (sf_cl > cellhd.cols - sf_x)
+	sf_cl = cellhd.cols - sf_x;
 
     /* calculate sample frame boundaries */
     sf_n = cellhd.north - (cellhd.ns_res * sf_y);
@@ -250,19 +265,28 @@ int parseSetup(char *path, struct list *l, struct g_area *g, char *raster)
 
 	    if (rel_sa_x == -1.0 && rel_sa_y == -1.0) {
 		/* runtime disposition */
-		int sa_rl, sa_cl;
 
 		sa_rl = (int)rint(cellhd.rows * rel_sa_rl);
 		sa_cl = (int)rint(cellhd.cols * rel_sa_cl);
+
+		/* sanity check */
+		if (sa_rl > cellhd.rows - sf_y)
+		    sa_rl = cellhd.rows - sf_y;
+		if (sa_cl > cellhd.cols - sf_x)
+		    sa_cl = cellhd.cols - sf_x;
+
+		/* total sample area */
 		g->rows = sf_rl;
 		g->cols = sf_cl;
+		g->x = sf_x;
+		g->y = sf_y;
+		/* current sample area (subset of total sample area) */
 		g->rl = sa_rl;
 		g->cl = sa_cl;
-		g->count = 1;
 		g->sf_x = sf_x;
 		g->sf_y = sf_y;
-		g->x = sf_x;
-		g->y = sf_y;
+
+		g->count = 1;
 		g->maskname = NULL;
 
 		return disposeAreas(l, g, strtok(NULL, "\n"));
@@ -273,10 +297,30 @@ int parseSetup(char *path, struct list *l, struct g_area *g, char *raster)
 		toReturn = NORMAL;
 		/*read file and create list */
 		m.type = AREA;
-		m.f.f_a.x = (int)rint(cellhd.cols * rel_sa_x);
-		m.f.f_a.y = (int)rint(cellhd.rows * rel_sa_y);
-		m.f.f_a.rl = (int)rint(cellhd.rows * rel_sa_rl);
-		m.f.f_a.cl = (int)rint(cellhd.cols * rel_sa_cl);
+		/* current sample area (subset of total sample area) */
+		sa_x = (int)rint(cellhd.cols * rel_sa_x);
+		sa_y = (int)rint(cellhd.rows * rel_sa_y);
+		sa_rl = (int)rint(cellhd.rows * rel_sa_rl);
+		sa_cl = (int)rint(cellhd.cols * rel_sa_cl);
+
+		/* sanity check */
+		if (sa_x < 0)
+		    sa_x = 0;
+		if (sa_y < 0)
+		    sa_y = 0;
+		if (sa_x > cellhd.cols)
+		    sa_x = cellhd.cols;
+		if (sa_y > cellhd.rows)
+		    sa_y = cellhd.rows;
+		if (sa_rl > cellhd.rows - sa_y)
+		    sa_rl = cellhd.rows - sa_y;
+		if (sa_cl > cellhd.cols - sa_x)
+		    sa_cl = cellhd.cols - sa_x;
+
+		m.f.f_a.x = sa_x;
+		m.f.f_a.y = sa_y;
+		m.f.f_a.rl = sa_rl;
+		m.f.f_a.cl = sa_cl;
 		m.f.f_a.aid = aid;
 		aid++;
 		insertNode(l, m);
@@ -302,17 +346,25 @@ int parseSetup(char *path, struct list *l, struct g_area *g, char *raster)
 
 	    if (rel_sa_x == -1 && rel_sa_y == -1) {
 		/* runtime disposition */
-		int sa_rl, sa_cl;
 
 		sa_rl = (int)rint(cellhd.rows * rel_sa_rl);
 		sa_cl = (int)rint(cellhd.cols * rel_sa_cl);
+
+		/* sanity check */
+		if (sa_rl > cellhd.rows - sf_y)
+		    sa_rl = cellhd.rows - sf_y;
+		if (sa_cl > cellhd.cols - sf_x)
+		    sa_cl = cellhd.cols - sf_x;
+
+		/* total sample area */
 		g->rows = sf_rl;
 		g->cols = sf_cl;
+		g->x = sf_x;
+		g->y = sf_y;
+		/* current sample area (subset of total sample area) */
 		g->rl = sa_rl;
 		g->cl = sa_cl;
 		g->count = 1;
-		g->x = sf_x;
-		g->y = sf_y;
 		g->maskname = maskname;
 		return disposeAreas(l, g, strtok(NULL, "\n"));
 	    }
@@ -321,10 +373,30 @@ int parseSetup(char *path, struct list *l, struct g_area *g, char *raster)
 		msg m;
 
 		m.type = MASKEDAREA;
-		m.f.f_ma.x = (int)rint(cellhd.cols * rel_sa_x);
-		m.f.f_ma.y = (int)rint(cellhd.rows * rel_sa_y);
-		m.f.f_ma.rl = (int)rint(cellhd.rows * rel_sa_rl);
-		m.f.f_ma.cl = (int)rint(cellhd.cols * rel_sa_cl);
+		/* current sample area (subset of total sample area) */
+		sa_x = (int)rint(cellhd.cols * rel_sa_x);
+		sa_y = (int)rint(cellhd.rows * rel_sa_y);
+		sa_rl = (int)rint(cellhd.rows * rel_sa_rl);
+		sa_cl = (int)rint(cellhd.cols * rel_sa_cl);
+
+		/* sanity check */
+		if (sa_x < 0)
+		    sa_x = 0;
+		if (sa_y < 0)
+		    sa_y = 0;
+		if (sa_x > cellhd.cols)
+		    sa_x = cellhd.cols;
+		if (sa_y > cellhd.rows)
+		    sa_y = cellhd.rows;
+		if (sa_rl > cellhd.rows - sa_y)
+		    sa_rl = cellhd.rows - sa_y;
+		if (sa_cl > cellhd.cols - sa_x)
+		    sa_cl = cellhd.cols - sa_x;
+
+		m.f.f_ma.x = sa_x;
+		m.f.f_ma.y = sa_y;
+		m.f.f_ma.rl = sa_rl;
+		m.f.f_ma.cl = sa_cl;
 		m.f.f_ma.aid = aid;
 		strcpy(m.f.f_ma.mask, maskname);
 		aid++;
@@ -342,12 +414,11 @@ int parseSetup(char *path, struct list *l, struct g_area *g, char *raster)
 	double sa_n, sa_s, sa_w, sa_e;
 	int aid = 1;
 	char maskname[GNAME_MAX] = {'\0'};
-	struct Cell_head window;
 	msg m;
 
 	/* Get the window setting. g.region rast=<input raster> */
 	/*   ? same as cellhd above ? */
-	G_get_window(&window);
+	/* no. the current window might be different */
 
 	do {
 	    strcpy(maskname, strtok(NULL, "|"));
@@ -362,16 +433,36 @@ int parseSetup(char *path, struct list *l, struct g_area *g, char *raster)
 		raster via v.to.rast. See r.li.setup/sample_area_vector.sh.
 		This is used only for reading the region (NS, EW). */
 
+	    /* current sample area (subset of total sample area) */
+
 	    /* Get start x and y position of masked overlay raster with
-		respect to input raster region from window.
-		sa_n, sa_e are read from configuration file. */
-	    m.f.f_ma.x = (int)Rast_easting_to_col(sa_w, &window);
-	    m.f.f_ma.y = (int)Rast_northing_to_row(sa_n, &window);
+		respect to current region.
+		sa_n, sa_w are read from configuration file. */
+	    sa_x = (int)rint((sa_w - cellhd.west) / cellhd.ew_res);
+	    sa_y = (int)rint((cellhd.north - sa_n) / cellhd.ns_res);
 
 	    /* Get row count and column count of overlay raster */
-	    m.f.f_ma.rl = (int)rint((sa_n - sa_s) / cellhd.ns_res);
-	    m.f.f_ma.cl = (int)rint((sa_e - sa_w) / cellhd.ew_res);
-
+	    sa_rl = (int)rint((sa_n - sa_s) / cellhd.ns_res);
+	    sa_cl = (int)rint((sa_e - sa_w) / cellhd.ew_res);
+
+	    /* sanity check */
+	    if (sa_x < 0)
+		sa_x = 0;
+	    if (sa_y < 0)
+		sa_y = 0;
+	    if (sa_x > cellhd.cols)
+		sa_x = cellhd.cols;
+	    if (sa_y > cellhd.rows)
+		sa_y = cellhd.rows;
+	    if (sa_rl > cellhd.rows - sa_y)
+		sa_rl = cellhd.rows - sa_y;
+	    if (sa_cl > cellhd.cols - sa_x)
+		sa_cl = cellhd.cols - sa_x;
+
+	    m.f.f_ma.x = sa_x;
+	    m.f.f_ma.y = sa_y;
+	    m.f.f_ma.rl = sa_rl;
+	    m.f.f_ma.cl = sa_cl;
 	    m.f.f_ma.aid = aid;
 	    strcpy(m.f.f_ma.mask, maskname);
 	    aid++;

raster/r.li/r.li.daemon/daemon.h

@@ -96,10 +96,10 @@ typedef struct fcell_memory_entry *fcell_manager;
 
  /**
   * \brief fields of an area descriptor
-  * \member x the x coordinate of upper left corner
-  * \member y the y coordinate of upper left corner
-  * \member rl area length in rows
-  * \member cl area length in columns
+  * \member x column offset = start of sample area
+  * \member y row offset = start of sample area
+  * \member rl sample area length in rows
+  * \member cl sample area length in columns
   * \member rc number of rows in the cache
   * \member mask file descriptor of mask raster file (-1 if there is no mask)
  */
@@ -212,7 +212,7 @@ void worker_end(void);
   * \param cl the lenght in cols of sample area
   * \return the name of mask raster file to use
   */
-char *mask_preprocessing(char *mask, char *raster, int rl, int cl);
+char *mask_preprocessing(char *mask, char *raster, struct area_entry *ad);
 
  /**
   * \brief writes the output for a raster file

raster/r.li/r.li.daemon/list.h

@@ -68,13 +68,15 @@ void removeNode(struct list *l);
   * \param dist inter-area distance
   * \param add_row cell to add in rows
   * \param add_col cell to add in columns
-  * \param x x coordinate of next area
-  * \param y y coordinate of next area
-  * \param rl area length in rows
-  * \param cl areal length in columns
+  * \param rows area length in rows
+  * \param cols area length in columns
+  * \param x column offset of next area
+  * \param y row offset of next area
+  * \param rl sample area length in rows
+  * \param cl sample area length in columns
   * \param count identifier of next area
-  * \param sf_x x coordinate of sample frame
-  * \param sf_y y coordinate of sample frame
+  * \param sf_x column offset of sample frame
+  * \param sf_y row offset of sample frame
   * \param maskname name of mask for the area
   */
 struct g_area

raster/r.li/r.li.daemon/r.li.daemon.html

@@ -88,7 +88,7 @@ Claudio Porta and Lucio Davide Spano, students of Computer Science
 University of Pisa (Italy).<br>
 Commission from Faunalia Pontedera (PI)<br>
 Rewritten from "r.le.setup" by William L. Baker<br>
-Rewritten to produce correct results by Markus Metz<br>
+Various bug fixes by Markus Metz<br>
 This program is free software under the GPL (>=v2)
 Read the COPYING file that comes with GRASS for details.
 

raster/r.li/r.li.daemon/worker.c

@@ -64,7 +64,9 @@ void worker_init(char *r, int f(int, char **, struct area_entry *, double *), ch
 
     /* open raster map */
     fd = Rast_open_old(raster, "");
-    Rast_get_cellhd(raster, "", &hd);
+
+    /* get current window */
+    Rast_get_window(&hd);
 
     /* read data type to allocate cache */
     data_type = Rast_map_type(raster, "");
@@ -127,8 +129,8 @@ void worker_process(msg * ret, msg * m)
 	ad->raster = raster;
 
 	/* mask preprocessing */
-	ad->mask_name = mask_preprocessing(m->f.f_ma.mask,
-					   raster, ad->rl, ad->cl);
+	ad->mask_name = mask_preprocessing(m->f.f_ma.mask, raster, ad);
+
 	if (ad->mask_name == NULL) {
 	    G_message(_("unable to open <%s> mask ... continuing without!"),
 		      m->f.f_ma.mask);
@@ -145,6 +147,12 @@ void worker_process(msg * ret, msg * m)
 	G_fatal_error("Program error, worker() type=%d", m->type);
 	break;
     }
+    
+    /* sanity check on the sample area ? */
+    /* 0 <= ad->x < hd.cols */
+    /* 0 <= ad->y < hd.rows */
+    /* ad->rl + ad->y <= hd.rows */
+    /* ad->cl + ad->x <= hd.cols */
 
     /* memory menagement */
     if (ad->rc > used) {
@@ -208,45 +216,33 @@ void worker_end(void)
     Rast_close(fd);
 }
 
-char *mask_preprocessing(char *mask, char *raster, int rl, int cl)
+char *mask_preprocessing(char *mask, char *raster, struct area_entry *ad)
 {
     const char *tmp_file;
-    struct Cell_head cell, oldcell;
     int mask_fd, old_fd, *buf, i, j;
     CELL *old;
-    double add_row, add_col;
 
-    buf = G_malloc(cl * sizeof(int));
+    buf = G_malloc(ad->cl * sizeof(int));
 
     G_debug(3, "daemon mask preproc: raster=[%s] mask=[%s]  rl=%d cl=%d",
-	    raster, mask, rl, cl);
-
-    /* open raster */
-    Rast_get_cellhd(raster, "", &cell);
-
-    /* open raster */
-    Rast_get_cellhd(mask, "", &oldcell);
-
-    add_row = 1.0 * oldcell.rows / rl;
-    add_col = 1.0 * oldcell.cols / cl;
+	    raster, mask, ad->rl, ad->cl);
 
     tmp_file = G_tempfile();
     mask_fd = open(tmp_file, O_RDWR | O_CREAT, 0755);
     old_fd = Rast_open_old(mask, "");
     old = Rast_allocate_c_buf();
 
-    for (i = 0; i < rl; i++) {
-	int riga;
+    /* write out sample area size: ad->rl rows and ad->cl columns */
+
+    for (i = 0; i < ad->rl; i++) {
 
-	riga = (int)rint(i * add_row);
-	Rast_get_c_row_nomask(old_fd, old, riga);
-	for (j = 0; j < cl; j++) {
-	    int colonna;
+	Rast_get_c_row_nomask(old_fd, old, i + ad->y);
+	for (j = 0; j < ad->cl; j++) {
 
-	    colonna = (int)rint(j * add_col);
-	    buf[j] = old[colonna];
+	    /* NULL -> 0, else 1 */
+	    buf[j] = !Rast_is_c_null_value(&old[j + ad->x]);
 	}
-	if (write(mask_fd, buf, cl * sizeof(int)) < 0)
+	if (write(mask_fd, buf, ad->cl * sizeof(int)) < 0)
 	    return NULL;
     }