[GRASS-SVN] r65520 - grass-addons/grass7/vector/v.kriging

svn_grass at osgeo.org svn_grass at osgeo.org
Thu Jun 25 08:05:04 PDT 2015


Author: evas
Date: 2015-06-25 08:05:04 -0700 (Thu, 25 Jun 2015)
New Revision: 65520

Modified:
   grass-addons/grass7/vector/v.kriging/Makefile
   grass-addons/grass7/vector/v.kriging/geostat.c
   grass-addons/grass7/vector/v.kriging/local_proto.h
   grass-addons/grass7/vector/v.kriging/spatial_index.c
   grass-addons/grass7/vector/v.kriging/utils.c
   grass-addons/grass7/vector/v.kriging/utils_kriging.c
Log:
v.kriging: fixed bug in search point definition

Modified: grass-addons/grass7/vector/v.kriging/Makefile
===================================================================
--- grass-addons/grass7/vector/v.kriging/Makefile	2015-06-24 22:23:40 UTC (rev 65519)
+++ grass-addons/grass7/vector/v.kriging/Makefile	2015-06-25 15:05:04 UTC (rev 65520)
@@ -1,4 +1,4 @@
-MODULE_TOPDIR = ../..
+MODULE_TOPDIR = /home/evergreen/grass7
 
 PGM = v.kriging
 

Modified: grass-addons/grass7/vector/v.kriging/geostat.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/geostat.c	2015-06-24 22:23:40 UTC (rev 65519)
+++ grass-addons/grass7/vector/v.kriging/geostat.c	2015-06-25 15:05:04 UTC (rev 65520)
@@ -39,6 +39,7 @@
   int n = pnts->n;                  // # of input points
   double *pnts_r = pnts->r;         // xyz coordinates of input points
   double *r;                        // pointer to xyz coordinates
+  double *search;                   // pointer to search point coordinates
   double *vals = pnts->invals;      // values to be used for interpolation
   int i3 = xD->i3;
   int phase = xD->phase;            // phase: initial / middle / final
@@ -124,6 +125,8 @@
 
   // allocate:
   dr = (double *) G_malloc(3 * sizeof(double));   // vector of coordinate differences
+  search = (double *) G_malloc(3 * sizeof(double));   // vector of search point coordinates
+
   if (type != 2) {
     var_pars->h = (double *) G_malloc(nLag * sizeof(double)); // vector of bins
   }
@@ -134,7 +137,7 @@
   }
 
   // control initialization:
-  if (dr == NULL || var_pars->h == NULL || (type == 2 && var_pars->vert == NULL)) {
+  if (dr == NULL || search == NULL || var_pars->h == NULL || (type == 2 && var_pars->vert == NULL)) {
     if (xD->report.write2file == TRUE) { // close report file
       fprintf(xD->report.fp, "Error (see standard output). Process killed...");
       fclose(xD->report.fp);
@@ -189,15 +192,16 @@
 
       /* Compute variogram for points in relevant neighbourhood */
       for (i = 0; i < n-1; i++) {    // for each input point...	
+	search = &pnts->r[3 * i];
        	switch (type) { // find NNs according to variogram type
 	case 0: // horizontal variogram
-	  list = find_NNs_within(2, i, pnts, max_dist, -1);
+	  list = find_NNs_within(2, search, pnts, max_dist, -1);
 	  break;
 	case 1: // vertical variogram
-	  list = find_NNs_within(1, i, pnts, -1, max_dist);
+	  list = find_NNs_within(1, search, pnts, -1, max_dist);
 	  break;
 	default: // anisotropic or bivariate variogram    
-	  list = find_NNs_within(3, i, pnts, max_dist, max_dist_vert);
+	  list = find_NNs_within(3, search, pnts, max_dist, max_dist_vert);
 	  break;
 	}
 	
@@ -467,6 +471,8 @@
 
   int type = var_par->type;
   double max_dist = type == 2 ? var_par->horizontal.max_dist : var_par->max_dist;
+  //max_dist = sqrt(0.5 * SQUARE(max_dist));
+
   double max_dist_vert = type == 2 ? var_par->vertical.max_dist : var_par->max_dist;
   double ratio = var_par->type == 3 ? xD->aniso_ratio : 1.;
        
@@ -547,16 +553,14 @@
 	list = G_new_ilist();            // create list of overlapping rectangles
        	
 	if (i3 == TRUE) { // 3D kriging:
-	  list = find_NNs_within(3, i, pnts, max_dist, max_dist_vert);
+	  list = find_NNs_within(3, r0, pnts, max_dist, max_dist_vert);
 	}
 	else { // 2D kriging:
-	  list = find_NNs_within(2, i, pnts, max_dist, max_dist_vert);
+	  list = find_NNs_within(2, r0, pnts, 3. * max_dist, max_dist_vert);
 	}
-	 
-	n_vals = list->n_values;         // # of selected points
-	if (n_vals > 0) { // positive # of selected points: 
+
+	if (list->n_values > 1) { // positive # of selected points: 
 	  correct_indices(i3, list, r0, pnts, var_par);
-
 	  GM_sub = submatrix(list, GM, report); // make submatrix for selected points
 	  GM_Inv = G_matrix_inverse(GM_sub);    // invert submatrix
 	  G_matrix_free(GM_sub);
@@ -568,31 +572,34 @@
 	  G_matrix_free(g0);
 
 	  rslt_OK = result(pnts, list, w0); // Estimated cell/voxel value rslt_OK = w x inputs
-
-  	  // Create output
-	constant_voxel_val: 
-	  if (var_par->const_val == 1) { // constant input values:
-	    rslt_OK = (double) *vals; // setup input as output
-	  }
-
-	  // write output to the (3D) raster layer
-	  if (write2layer(xD, reg, out, col, row, dep, rslt_OK) == 0) {
-	    if (report->name) { // report file available
-	      fprintf(report->fp, "Error (see standard output). Process killed...");
-	      fclose(report->fp); // close report file
-	    }
-	    G_fatal_error(_("Error writing result into output layer..."));
-	  }
-	  G_free_ilist(list);    // free list memory
-	  
-	} // end if searchRadius
-	else { // no selected points:
+	}
+	else if (list->n_values == 1) {
+	  rslt_OK = vals[list->value[0] - 1]; // Estimated cell/voxel value rslt_OK = w x inputs
+	}
+	else if (list->n_values == 0) {
 	  if (report->name) { // report file available:
 	    fprintf(report->fp, "Error (see standard output). Process killed...");
 	    fclose(report->fp);
 	  }
 	  G_fatal_error(_("This point does not have neighbours in given radius..."));
 	} // end else: error
+
+  	  // Create output
+      constant_voxel_val: 
+	if (var_par->const_val == 1) { // constant input values:
+	  rslt_OK = (double) *vals; // setup input as output
+	}
+
+	// write output to the (3D) raster layer
+	if (write2layer(xD, reg, out, col, row, dep, rslt_OK) == 0) {
+	  if (report->name) { // report file available
+	    fprintf(report->fp, "Error (see standard output). Process killed...");
+	    fclose(report->fp); // close report file
+	  }
+	  G_fatal_error(_("Error writing result into output layer..."));
+	}
+
+	G_free_ilist(list);    // free list memory  
       } // end col
     } // end row 
   } // end dep

Modified: grass-addons/grass7/vector/v.kriging/local_proto.h
===================================================================
--- grass-addons/grass7/vector/v.kriging/local_proto.h	2015-06-24 22:23:40 UTC (rev 65519)
+++ grass-addons/grass7/vector/v.kriging/local_proto.h	2015-06-25 15:05:04 UTC (rev 65520)
@@ -192,7 +192,7 @@
 
 struct RTree *create_spatial_index(struct int_par *);
 void insert_rectangle(int, int, struct points *);
-struct ilist *find_NNs_within(int, int, struct points *, double, double);
+struct ilist *find_NNs_within(int, double *, struct points *, double, double);
 struct ilist *find_n_NNs(int, int, struct points *, int);
 double sum_NN(int, int, struct ilist *, struct points *);
 

Modified: grass-addons/grass7/vector/v.kriging/spatial_index.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/spatial_index.c	2015-06-24 22:23:40 UTC (rev 65519)
+++ grass-addons/grass7/vector/v.kriging/spatial_index.c	2015-06-25 15:05:04 UTC (rev 65520)
@@ -34,42 +34,51 @@
 }
 
 /* find neighbors in cubic (square) surroundings */
-struct ilist *find_NNs_within(int dim, int i, struct points *pnts, double max_dist, double max_dist_vert)
+struct ilist *find_NNs_within(int dim, double *search_pt, struct points *pnts, double max_dist, double max_dist_vert)
 {
   // local variables
-  double *r;                        // pointer to vector of the coordinates
-  r = &pnts->r[3 * i];
+  double *r;                         // pointer to vector of the coordinates
+  r = &search_pt[0];
 
-  struct RTree_Rect *search;        // search rectangle
+  int NN_sum = 0;                    // # of NN
+  double dist_step = max_dist; // distance iteration in case of empty closest surrounding of search point
+  double dist_step_vert = 0.1 * max_dist_vert; // vertical distance iteration
+  struct RTree_Rect *search;         // search rectangle
   struct ilist *list;
   list = G_new_ilist();
-  
-  switch (dim) {
-  case 3: // 3D:
-    search = RTreeAllocRect(pnts->R_tree); // allocate new rectangle
-    RTreeSetRect3D(search, pnts->R_tree, 
+
+  while (NN_sum < 2) { // TODO: this might be tricky - what # (density) is optimal?
+    switch (dim) {
+    case 3: // 3D:
+      search = RTreeAllocRect(pnts->R_tree); // allocate new rectangle
+      RTreeSetRect3D(search, pnts->R_tree, 
 		   *r - max_dist, *r + max_dist, 
 		   *(r+1) - max_dist, *(r+1) + max_dist, 
-		   *(r+2) - max_dist_vert, *(r+2) + max_dist_vert); // set up searching rectangle
-    RTreeSearch2(pnts->R_tree, search, list);   // search the nearest rectangle
-    break;
-  case 2: // 2D:
-    search = RTreeAllocRect(pnts->Rtree_hz); // allocate new rectangle
-    RTreeSetRect2D(search, pnts->Rtree_hz, 
+		     *(r+2) - max_dist_vert, *(r+2) + max_dist_vert); // set up searching rectangle
+      RTreeSearch2(pnts->R_tree, search, list);   // search the nearest rectangle
+      break;
+    case 2: // 2D:
+      search = RTreeAllocRect(pnts->Rtree_hz); // allocate new rectangle
+      RTreeSetRect2D(search, pnts->Rtree_hz, 
 		   *r - max_dist, *r + max_dist, 
 		   *(r+1) - max_dist, *(r+1) + max_dist);           // set up searching rectangle
-    RTreeSearch2(pnts->Rtree_hz, search, list);   // search the nearest rectangle
-    break;
-  case 1: // 1D: 
-    //G_debug(0, "%d %f %f", pnts->Rtree_vert->n_nodes, *(r+2), max_dist_vert);   
-    search = RTreeAllocRect(pnts->Rtree_vert); // allocate new rectangle
-    RTreeSetRect1D(search, pnts->Rtree_vert, 
+      RTreeSearch2(pnts->Rtree_hz, search, list);   // search the nearest rectangle
+      break;
+    case 1: // 1D:   
+      search = RTreeAllocRect(pnts->Rtree_vert); // allocate new rectangle
+      RTreeSetRect1D(search, pnts->Rtree_vert, 
 		   *(r+2) - max_dist_vert, *(r+2) + max_dist_vert); // set up searching rectangle
-    RTreeSearch2(pnts->Rtree_vert, search, list);   // search the nearest rectangle
-    break;
-  }
-  
-  RTreeFreeRect(search);
+
+      RTreeSearch2(pnts->Rtree_vert, search, list);   // search the nearest rectangle
+      break;
+    }
+    RTreeFreeRect(search);
+
+    NN_sum = list->n_values;
+    max_dist += dist_step;
+    max_dist_vert += dist_step_vert;
+  } // end while: nonzero # of NNs
+
   return list;
 }
 
@@ -83,9 +92,11 @@
   int n_vals = 0;              // # of the nearest neighbours (NN)
   int iter = 1;                // iteration number (searching NN)
   double max_dist = max_dist0; // iterated search radius
+  double *search;
+  search = &pnts->r[3 * i];
 
   while (n_vals < n) { // until some NN is found (2 points: identical and NN)
-    list = find_NNs_within(dim, i, pnts, max_dist, -1);
+    list = find_NNs_within(dim, search, pnts, max_dist, -1);
     n_vals = list->n_values;                 // set up control variable
     
     if (n_vals < n) {
@@ -95,7 +106,7 @@
   }
 
   // check spherical (circular) surrounding
-  list = find_NNs_within(dim, i, pnts, sqrt(2.) * max_dist, -1);
+  list = find_NNs_within(dim, search, pnts, sqrt(2.) * max_dist, -1);
 
   return list;
 }

Modified: grass-addons/grass7/vector/v.kriging/utils.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils.c	2015-06-24 22:23:40 UTC (rev 65519)
+++ grass-addons/grass7/vector/v.kriging/utils.c	2015-06-25 15:05:04 UTC (rev 65520)
@@ -21,34 +21,43 @@
   int i;
   int n = list->n_values;
   int *vals = list->value;
-  double *pnt = pnts->r;
+  double dx, dy, dz;
   int type = var_pars->type;
   double sqDist = type == 2 ? SQUARE(var_pars->horizontal.max_dist) : SQUARE(var_pars->max_dist);
   
-  int n_new = 0;
-  int *newvals;
+  int n_new = 0;  // # of effective indices (not identical points, nor too far)
+  int *newvals, *save;   // effective indices
   newvals = (int *) G_malloc(n * sizeof(int));
+  save = &newvals[0];
 
-  double *r;
+  double *r, sqDist_i;
 
   for (i=0; i<n; i++) {
     *vals -= 1;
-    r = &pnt[3 * *vals];
+    r = &pnts->r[3 * *vals];
 
-    if (type != 2 && SQUARE(*r - *r0) + SQUARE(*(r+1) - *(r0+1)) + SQUARE(*(r+2) - *(r0+2)) <= sqDist) { 
-      *newvals = *vals;
+    dx = *r - *r0;
+    dy = *(r+1) - *(r0+1);
+    dz = *(r+2) - *(r0+2);
+    sqDist_i = SQUARE(dx) + SQUARE(dy) + SQUARE(dz);
+    
+    if (type != 2 && (n == 1 || (sqDist_i != 0.))) { 
+      *save = *vals;
       n_new++;
-      newvals++;
+      save++;
     }
     vals++;
   }
-  
+
   if (type != 2 && n_new < n) {
     list->n_values = n_new;
-    G_realloc(list->value, n_new);
-    G_realloc(newvals, n_new);
+    list->value = (int *) G_realloc(list->value, n_new * sizeof(int));
     memcpy(list->value, newvals, n_new * sizeof(int)); 
   }
+
+  G_free(newvals);
+
+  return;
 }
 
 // make coordinate triples xyz
@@ -60,6 +69,8 @@
   *t = x;
   *(t+1) = y;
   *(t+2) = z;
+
+  return;
 }
 
 // compute coordinate differences
@@ -78,6 +89,8 @@
   *drt = *rl - *rk; // dx
   *(drt+1) = *(rl+1) - *(rk+1); // dy
   *(drt+2) = *(rl+2) - *(rk+2); // dz
+
+  return;
 }
 
 // compute horizontal radius from coordinate differences
@@ -85,6 +98,7 @@
 {
   double rds;
   rds = SQUARE(dr[0]) + SQUARE(dr[1]);
+
   return rds;
 }
 
@@ -119,22 +133,26 @@
   double lagNN;                    // square root of d_min -> lag
   double dist2;                    // square of the distances between search point and NN
 
+  double *search;                  // search point
+  search = (double *) G_malloc(3 * sizeof(double));
+
   int perc5;                       // 5% from total number of NN
   int add_ident;                   // number of identical points to add to perc5
   
   for (i=0; i < n; i++) { // for each input point...
   
     add_ident = 0;        // number of points identical to search point equals to zero
+    search = &pnts->r[3 * i];
 
     switch (direction) {
     case 12: // horizontal variogram
-      list = find_NNs_within(2, i, pnts, max_dist, -1); 
+      list = find_NNs_within(2, search, pnts, max_dist, -1); 
       break;
     case 3: // vertical variogram
-      list = find_NNs_within(1, i, pnts, max_dist, max_dist_vert);    
+      list = find_NNs_within(1, search, pnts, max_dist, max_dist_vert);    
       break;
     default: // anisotropic and bivariate variogram  
-      list = find_NNs_within(3, i, pnts, max_dist, max_dist_vert);   
+      list = find_NNs_within(3, search, pnts, max_dist, max_dist_vert);   
       break;
     }
     n_vals = list->n_values;                 // number of overlapping rectangles
@@ -196,9 +214,10 @@
     r += 3; // go to next search point  
   } // end i for loop: distance of i-th search pt and 5%-th NN 
 
-  G_free_ilist(list); // free list memory 
   lagNN = sqrt(d_min);
 
+  G_free_ilist(list); // free list memory
+ 
   return lagNN;
 }
 

Modified: grass-addons/grass7/vector/v.kriging/utils_kriging.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_kriging.c	2015-06-24 22:23:40 UTC (rev 65519)
+++ grass-addons/grass7/vector/v.kriging/utils_kriging.c	2015-06-25 15:05:04 UTC (rev 65520)
@@ -592,6 +592,7 @@
   double max_dist = type == 2 ? var_par->horizontal.max_dist : var_par->max_dist;
   double max_dist_vert = type == 2 ? var_par->vertical.max_dist : var_par->max_dist;
   
+  double *search;               // coordinates of the search point
   struct ilist *list;
   struct RTree_Rect *r_cell;
   mat_struct *GM_sub;
@@ -602,6 +603,7 @@
   struct write *report = &xD->report;
   double *normal, *absval, *norm, *av;
 
+  search = (double *) G_malloc(3 * sizeof(double));
   normal = (double *) G_malloc(n * sizeof(double));
   absval = (double *) G_malloc(n * sizeof(double));
   norm = &normal[0];
@@ -615,11 +617,12 @@
   for (i=0; i<n; i++) {   // for each input point [r0]:
     list = G_new_ilist();                  // create list of overlapping rectangles
 
+    search = &pnts->r[3 * i];
     if (i3 == TRUE) {
-      list = find_NNs_within(3, i, pnts, max_dist, max_dist_vert);
+      list = find_NNs_within(3, search, pnts, max_dist, max_dist_vert);
     }
     else {
-      list = find_NNs_within(2, i, pnts, max_dist, max_dist_vert);
+      list = find_NNs_within(2, search, pnts, max_dist, max_dist_vert);
     }
 	 
     n_vals = list->n_values;               // # of overlapping rectangles
@@ -633,7 +636,7 @@
       
       g0 = set_up_g0(xD, pnts, list, r, var_par); // Diffs inputs - unknowns (incl. cond. 1)) 
       w0 = G_matrix_product(GM_Inv, g0); // Vector of weights, condition SUM(w) = 1 in last row
-       
+      
       G_matrix_free(g0);
       G_matrix_free(GM_Inv);
 



More information about the grass-commit mailing list