[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