[GRASS-SVN] r65505 - in grass-addons/grass7/vector/v.kriging: . backup_version_with_PCL
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jun 20 05:01:04 PDT 2015
Author: evas
Date: 2015-06-20 05:01:04 -0700 (Sat, 20 Jun 2015)
New Revision: 65505
Added:
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/Makefile
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/geostat.cpp
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/getval.cpp
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/local_proto.h
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/main.cpp
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/quantile.cpp
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/stat.cpp
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils.cpp
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_kriging.cpp
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_raster.cpp
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_write.cpp
grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/v.kriging.html
Log:
v.kriging: backup of PCL version
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/Makefile
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/Makefile 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,22 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.kriging
+
+# other PCL version:
+#-lpcl_apps -lpcl_common -lpcl_features -lpcl_filters -lpcl_io -lpcl_kdtree -lpcl_keypoints -lpcl_octree -lpcl_registration -lpcl_sample_consensus -lpcl_search -lpcl_segmentation -lpcl_surface -lpcl_visualization
+LIBES = $(RASTER3DLIB) $(RASTERLIB) $(VECTORLIB) $(DBMILIB) $(GMATHLIB) $(GISLIB) $(IOSTREAMLIB) $(RTREELIB) -lpcl_common -lpcl_features -lpcl_filters -lpcl_io_ply -lpcl_io -lpcl_kdtree -lpcl_keypoints -lpcl_octree -lpcl_outofcore -lpcl_people -lpcl_recognition -lpcl_registration -lpcl_sample_consensus -lpcl_search -lpcl_segmentation -lpcl_surface -lpcl_tracking -lpcl_visualization
+
+DEPENDENCIES= $(RASTER3DDEP) $(RASTERDEP) $(VECTORDEP) $(DBMIDEP) $(GMATHDEP) $(GISDEP) $(IOSTREAMDEP) $(RTREEDEP)
+EXTRA_INC = $(VECT_INC) -I/usr/include/pcl-1.7 -I/usr/include/eigen3 $(PCL_INCLUDE_DIRS) $(PROJINC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+EXTRA_CFLAGS=-fopenmp
+EXTRA_LIBS=$(GISLIB) -lgomp $(MATHLIB) -lboost_system
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+LINK = $(CXX)
+
+ifneq ($(strip $(CXX)),)
+default: cmd
+endif
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/geostat.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/geostat.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/geostat.cpp 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,720 @@
+#include "local_proto.h"
+
+/* experimental variogram
+ * based on 2D variogram (alghalandis.com/?page_id=463)) */
+void E_variogram(int type, struct int_par *xD, struct points *pnts, struct pcl_utils *pclp, struct var_par *pars)
+{
+ // Variogram type
+ struct parameters *var_par, *var_par_vert;
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts;
+ switch (type) {
+ case 0:
+ var_par = &pars->hz;
+ pcl_pnts = pclp->pnts_hz;
+ break;
+ case 1:
+ var_par = &pars->vert;
+ pcl_pnts = pclp->pnts_vert;
+ break;
+ default:
+ var_par = &pars->fin;
+ pcl_pnts = pclp->pnts;
+ break;
+ }
+
+ // Local variables
+ int n = pnts->n; // # of input points
+ double *r; // xyz coordinates of input points
+ double *vals = pnts->invals; // values to be used for interpolation
+ int phase = xD->phase;
+ int function = var_par->function; // type of theoretical variogram
+
+ int nLag = var_par->nLag; // # of horizontal bins
+ int nLag_vert;
+ double *vert;
+ double dir = var_par->dir; // azimuth
+ double td = var_par->td; // angle tolerance
+ double lag = var_par->lag; // size of horizontal bin
+ double lag_vert;
+
+ /*switch (type) {
+ case 1:
+ td = var_par->tz;
+ break;
+ case 2:
+ td = var_par->horizontal->td;
+ tz = var_par->vertical->td;
+ break;
+ default:
+ td = var_par->td;
+ break;
+ }*/
+
+ // depend on variogram type:
+ int k; // index of variogram to be computed
+ int nrows = nLag;
+
+ double max_dist = var_par->max_dist; // max distance of interpolated points
+ double radius = SQUARE(max_dist); // radius of interpolation in the horizontal plane
+ struct write *report = &xD->report;
+
+ if (type == 2) {
+ nLag_vert = var_par->nLag_vert; // # of horizontal bins
+ lag_vert = var_par->lag_vert; // size of horizontal bin
+ }
+
+ int ncols = type == 2 ? nLag_vert : 1;
+
+ // Variogram processing variables
+ int s; // index of horizontal segment
+ int b; // index of verical segment (bivariate only)
+ int i, j; // indices of the input points
+ int err0=0; // number of invalid values
+
+ double *dr; // coordinate differences of each couple
+ double *h; // distance of boundary of the horizontal segment
+ double tv; // bearing of the line between the couple of the input points
+ double ddir; // the azimuth of computing variogram
+ double rv; // radius of the couple of points
+ double drv; // distance between the points
+ double rvh; // difference between point distance and the horizontal segment boundary
+ double dv; // difference of the values to be interpolated that are located on the couple of points
+
+ double *i_vals, *j_vals; // values located on the point couples
+
+ int n1 = n+1; // number of points + 1
+ pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree (new pcl::KdTreeFLANN<pcl::PointXYZ>);
+
+ pcl::PointXYZ *searchPt; // search point (each input point)
+ std::vector<int> ind; // vector of indices
+ std::vector<float> sqDist;
+ int n_ind; // number of elements (including 1st irrelevant)
+ int *indices; // save indices into matrix to use them later
+ int *idx; // index to be set (output)
+ int *ii, *i0; // difference of indices between rellevant input points
+ int *ip; // index to be set into the index matrix (input)
+
+ double gamma_sum; // sum of dissimilarities in one bin
+ double cpls; // # of dissimilarities in one bin
+ double gamma_E; // average of dissimilarities in one bin (element of gamma matrix)
+ mat_struct *gamma_temp; // gamma matrix (hz, vert or bivar)
+ mat_struct *c_temp; // matrix of # of dissimilarities
+ double *gamma; // pointer to gamma matrix
+ double *c; // pointer to c matrix
+
+ unsigned int percents = 50;
+ int counter;
+ int count;
+
+ double gamma_stat; // sum of gamma elements (non-nan)
+ int gamma_n; // # of gamma elements (non-nan)
+ double gamma_sill; // sill
+
+ double nugget, sum_nugget = 0.;
+ int n_nugget = 0;
+
+ /* Allocated vertices and matrices:
+ * --------------------------------
+ * dr - triple of coordinate differences (e.g. dx, dy, dz)
+ * h - vector of length pieces [nL x 1]
+ * vert - vector of vertical pieces [nZ x 1] (3D only)
+ * c - matrix containing number of points in each sector defined by h (and vert) [nL x nZ; in 2D nZ = 1]
+ * gama - matrix of values of experimental variogram
+ */
+
+ // allocate:
+ dr = (double *) G_malloc(3 * sizeof(double)); // vector of coordinate differences
+ var_par->h = (double *) G_malloc(nLag * sizeof(double)); // vector of bins
+
+ if (type == 2) {
+ var_par->vert = (double *) G_malloc(nLag_vert * sizeof(double)); // vector of bins
+ }
+
+ // control initialization:
+ if (dr == NULL || var_par->h == NULL || (type == 2 && var_par->vert == NULL)) {
+ if (xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Memory allocation failed..."));
+ }
+
+ r = &pnts->r[0]; // set up pointer to input coordinates
+ indices = (int *) G_malloc(n1*n * sizeof(int)); // matrix of indices of neighbouring points
+ if (indices == NULL) {
+ if (xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Memory allocation failed..."));
+ }
+
+ kd_tree->setInputCloud(pcl_pnts);
+ searchPt = &pcl_pnts->points[0];
+
+ // initialize:
+ c_temp = G_matrix_init(nrows, ncols, nrows); // temporal matrix of counts
+ gamma_temp = G_matrix_init(nrows, ncols, nrows); // temporal matrix (vector) of gammas
+
+ if (c_temp == NULL || gamma_temp == NULL) {
+ if (xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Memory allocation failed..."));
+ }
+
+ // set up pointers
+ c = &c_temp->vals[0];
+ gamma = &gamma_temp->vals[0];
+
+ // set up starting values
+ gamma_stat = 0.; //
+ gamma_n = 0; //
+
+ if (percents) {
+ G_percent_reset();
+ }
+ counter = nrows * ncols;
+ count = 0;
+
+ if (type == 2) {
+ vert = &var_par->vert[0];
+ }
+
+ for (b = 0; b < ncols; b++) {
+ if (type == 2) {
+ *vert = (b+1) * lag_vert;
+ }
+
+ // set up vector of distances...
+ h = &var_par->h[0]; // ... horizontal
+
+ for (s = 0; s < nrows; s++) { // process of the horizontal piece (isotrophy!!!)
+ *h = (s+1) * lag; // distance of horizontal lag
+
+ // for every bs cycle ...
+ i_vals = &vals[0]; // ... i-vals start at the beginning
+ gamma_sum = 0.; // gamma in dir direction and h distance
+ cpls = 0.; // count of couples in dir direction and h distance
+
+ /* Compute variogram for points in relevant neighbourhood */
+ for (i = 0; i < n-1; i++) { // for every input point...
+ if (b == 0 && s == 0) { // ... create vector of points in neighborhood
+ if (kd_tree->radiusSearch(*searchPt, max_dist, ind, sqDist) > 0) {
+ n_ind = ind.size (); // size of indices to 1st row
+ ip = &ind.data ()[1]; // values of indices to the next rows
+ idx = &indices[n1*i]; // begin new column
+ *idx = n_ind-1; // size of index vector as 1st element of row
+ for (j = 1; j < n_ind; j++) {
+ idx++; // next row
+ *idx = *ip; // set indices to the rest of the rows
+ ip++; // next index
+ } // end setting elements
+ } // end radiusSearch
+ else {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("There are no points in the surrounding of point no. %d..."), i);
+ } // end error
+ searchPt++; // next search point
+ } // end set up indices
+
+ idx = &indices[n1*i]; // ... new column of matrix is attached
+ ii = &indices[n1*i+1];
+ i0 = &indices[n1*i+1];
+ n_ind = *idx; // number of neighbouring points
+
+ j_vals = &vals[*ii];
+ for (j = 1; j < n_ind; j++) { // points within searchRadius
+ if (*ii > i) {
+ coord_diff(i, *ii, r, dr); // coordinate differences
+
+ /* Variogram processing */
+ if (type == 1) { // in case of vertical variogram...
+ tv = 0.5 * PI - atan2(dr[2], sqrt(SQUARE(dr[0]) + SQUARE(dr[1]))); // ... compute zenith angle instead of bearing
+ td = 0.5 * PI; // 0.25 usually
+ }
+ else
+ tv = atan2(*(dr+1), *dr); // bearing
+
+ ddir = tv - dir; // difference between bearing and azimuth
+ if (fabsf(ddir) <= td) { // is the difference acceptable according to tolerance?
+ // test distance
+ rv = type == 1 ? 0. : SQUARE(dr[0]) + SQUARE(dr[1]); // horizontal distance
+ if (type == 1 || type == 3) { // anisotropic only
+ rv += SQUARE(dr[2]); // vertical or anisotropic distance
+ }
+ rvh = sqrt(rv) - *h; // difference between distance and lag
+ if (rv <= radius && 0. <= rvh && rvh <= lag) { // distance test
+ if (type == 2) { // vertical test for bivariate variogram
+ rvh = *(dr+2) - *vert;
+ if (fabsf(rvh) <= lag_vert) { // elevation test
+ goto delta_V;
+ }
+ else {
+ goto end;
+ }
+ }
+ delta_V:
+ dv = *j_vals - *i_vals; // difference of values located on pair of points i, j
+ gamma_sum += dv * dv; // sum of squared differences
+ cpls += 1.; // count of elements
+ } // end distance test
+ } // end bearing test
+ //} // end dr2 test
+ } // end point selection
+ end:
+ *i0 = *ii;
+ ii++;
+ j_vals += *ii - *i0;
+ } // end j
+ i_vals++;
+ } // end i
+
+ if (isnan(gamma_sum) || cpls == 0.0) {
+ err0++;
+ goto gamma_isnan;
+ }
+
+ gamma_E = 0.5 * gamma_sum / cpls; // element of gamma matrix
+ *c = cpls; // # of values that have been used to compute gamma_e
+ *gamma = gamma_E;
+
+ gamma_stat += gamma_E; // sum of gamma elements (non-nan)
+ gamma_n++; // # of gamma elements (non-nan)
+
+ gamma_isnan:
+ h++;
+ c++;
+ gamma++;
+ count++;
+ } // end s
+
+ if (type == 2) {
+ vert++;
+ }
+ } // end b
+
+ free(indices);
+
+ if (err0 == nLag) { // todo> kedy nie je riesitelny teoreticky variogram?
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to compute experimental variogram..."));
+ } // end error
+
+ var_par->gamma = G_matrix_copy(gamma_temp);
+ var_par->c = G_matrix_copy(c_temp);
+ plot_experimental_variogram(xD, var_par->h, gamma_temp, var_par);
+
+ if (report->name) {
+ write2file_variogram_E(xD, var_par);
+ }
+
+ // Compute sill
+ if (phase < 2) {
+ switch (type) {
+ case 2:
+ int nh, nv;
+ double c_h, c_v, gamma_h, gamma_v;
+ char type[12];
+ nh = nv = 0;
+ gamma_h = gamma_v = 0.;
+
+ for (i=0; i<nrows; i++) {
+ c_h = G_matrix_get_element(c_temp, i, 0);
+ if (c_h != 0.) {
+ gamma_h += G_matrix_get_element(gamma_temp, i, 0);
+ nh++;
+ }
+ }
+
+ for (j=0; j<ncols; j++) {
+ c_v = G_matrix_get_element(c_temp, 0, j);
+ if (c_v != 0.) {
+ gamma_v += G_matrix_get_element(gamma_temp, 0, j);
+ nv++;
+ }
+ }
+
+ var_par->horizontal.sill = gamma_h / nh;
+ var_par->vertical.sill = gamma_v / nv;
+ break;
+
+ default:
+ var_par->sill = gamma_stat / gamma_n; // mean of gamma elements
+
+ variogram_type(var_par->type, type);
+ G_message(_("Default sill of %s variogram: %f"), type, var_par->sill);
+ break;
+ }
+ } // end compute sill
+
+ write_temporary2file(xD, var_par, gamma_temp);
+
+ G_matrix_free(c_temp);
+ G_matrix_free(gamma_temp);
+}
+
+/* theoretical variogram */
+void T_variogram(int type, struct opts opt, struct parameters *var_par, struct write *report)
+{
+ // report
+ if (report->name) {
+ time(&report->now);
+ if (type != 1) {
+ fprintf(report->fp, "\nComputation of theoretical variogram started on %s\n", ctime(&report->now));
+ }
+ }
+
+ // set up:
+ var_par->type = type;
+
+ char *variogram;
+ switch (type) {
+ case 0: // horizontal
+ var_par->nugget = atof(opt.nugget_hz->answer);
+ var_par->h_range = atof(opt.range_hz->answer);
+ if (opt.sill_hz->answer) {
+ var_par->sill = atof(opt.sill_hz->answer);
+ }
+ variogram = opt.function_var_hz->answer;
+ if (report->name) {
+ fprintf(report->fp, "Parameters of horizontal variogram:\n");
+ fprintf(report->fp, "Nugget effect: %f\n", var_par->nugget);
+ fprintf(report->fp, "Sill: %f\n", var_par->sill);
+ fprintf(report->fp, "Range: %f\n", var_par->h_range);
+ }
+ break;
+ case 1: // vertical
+ var_par->nugget = atof(opt.nugget_vert->answer);
+ var_par->h_range = atof(opt.range_vert->answer);
+ if (opt.sill_vert->answer) {
+ var_par->sill = atof(opt.sill_vert->answer);
+ }
+ variogram = opt.function_var_vert->answer;
+ if (report->name) {
+ fprintf(report->fp, "Parameters of vertical variogram:\n");
+ fprintf(report->fp, "Nugget effect: %f\n", var_par->nugget);
+ fprintf(report->fp, "Sill: %f\n", var_par->sill);
+ fprintf(report->fp, "Range: %f\n", var_par->h_range);
+ }
+ break;
+ case 2: // bivariate variogram
+ if (opt.function_var_hz->answer == NULL && opt.function_var_vert->answer == NULL) { // linear function
+ int nZ, nL, nr, i, j, nc;
+ double *h, *vert, *gamma, *c;
+ mat_struct *gR;
+
+ var_par->function = 5;
+
+ nL = var_par->gamma->rows; // # of cols (horizontal bins)
+ nZ = var_par->gamma->cols; // # of rows (vertical bins)
+ gamma = &var_par->gamma->vals[0]; // pointer to experimental variogram
+ c = &var_par->c->vals[0];
+
+ // Test length of design matrix
+ nr = 0; // # of rows of design matrix A - depend on # of non-nan gammas
+ for (i = 0; i < nZ; i++) {
+ for (j = 0; j < nL; j++) {
+ if (*c != 0.) { // todo: upgrade: !isnan(*c) L434 too
+ nr++;
+ }
+ gamma++;
+ c++;
+ } // end j
+ } // end i
+
+ // Number of columns of design matrix A
+ nc = var_par->function == 5 ? 3 : 2;
+
+ var_par->A = G_matrix_init(nr, nc, nr); // initialise design matrix
+ gR = G_matrix_init(nr, 1, nr); // initialise vector of observations
+
+ // Set design matrix A
+ nr = 0; // index of matrix row
+ gamma = &var_par->gamma->vals[0]; // pointer to experimental variogram
+ c = &var_par->c->vals[0];
+
+ if (var_par->function == 5) { // in case of bivariate variogram
+ vert = &var_par->vert[0]; // use separate variable for vertical direction
+ }
+
+ for (i = 0; i < nZ; i++) {
+ h = &var_par->h[0];
+ for (j = 0; j < nL; j++) {
+ if (*c != 0.) { // write to matrix - each valid element of gamma
+ switch (var_par->function) { // function of theoretical variogram
+ case 5: // bivariate planar
+ G_matrix_set_element(var_par->A, nr, 0, *h);
+ G_matrix_set_element(var_par->A, nr, 1, *vert);
+ G_matrix_set_element(var_par->A, nr, 2, 1.);
+ G_matrix_set_element(gR, nr, 0, *gamma);
+ break;
+ default:
+ G_matrix_set_element(var_par->A, nr, 0, *h);
+ G_matrix_set_element(var_par->A, nr, 1, 1.);
+ G_matrix_set_element(gR, nr, 0, *gamma);
+ break;
+ } // end switch variogram fuction
+ nr++; // length of vector of valid elements (not null)
+ } // end test if !isnan(*gamma)
+ h++;
+ gamma++;
+ c++;
+ } // end j
+ if (var_par->function == 5) {
+ vert++;
+ }
+ } // end i
+ end_loop:
+
+ // Estimate theoretical variogram coefficients
+ var_par->T = LSM(var_par->A, gR); // Least Square Method
+
+ // Test of theoretical variogram estimation
+ if (var_par->T->vals == NULL) { // NULL values of theoretical variogram
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error("Unable to compute 3D theoretical variogram...");
+ } // error
+
+ // constant raster
+ G_debug(0,"%f %f %f", var_par->T->vals[0], var_par->T->vals[1], var_par->T->vals[2]);
+ if (var_par->T->vals[0] == 0. && var_par->T->vals[1] == 0.) { //to do: otestuj pre 2d
+ var_par->const_val = 1;
+ G_message(_("Input values to be interpolated are constant."));
+ } // todo: cnstant raster for exponential etc.
+ else {
+ var_par->const_val = 0;
+ }
+
+ // coefficients of theoretical variogram (linear)
+ if (report->name) {
+ fprintf(report->fp, "Parameters of bivariate variogram:\n");
+ fprintf(report->fp, "Function: linear\n\n");
+ fprintf(report->fp, "gamma(h, vert) = %f * h + %f * vert + c\n", var_par->T->vals[0], var_par->T->vals[1], var_par->T->vals[2]);
+ }
+ }
+
+ else {
+ var_par->horizontal.nugget = atof(opt.nugget_hz->answer);
+ var_par->horizontal.h_range = atof(opt.range_hz->answer);
+ //var_par->horizontal.variogram = opt.variogram_hz->answer;
+ if (opt.sill_hz->answer) {
+ var_par->horizontal.sill = atof(opt.sill_hz->answer);
+ }
+
+ var_par->vertical.nugget = atof(opt.nugget_vert->answer);
+ var_par->vertical.h_range = atof(opt.range_vert->answer);
+ //var_par->vertical.variogram = opt.variogram_vert->answer;
+ if (opt.sill_hz->answer) {
+ var_par->vertical.sill = atof(opt.sill_vert->answer);
+ }
+
+ var_par->horizontal.function = set_function(opt.function_var_hz->answer, var_par, report);
+ var_par->vertical.function = set_function(opt.function_var_vert->answer, var_par, report);
+
+ if (report->name) {
+ fprintf(report->fp, "Parameters of bivariate variogram:\n");
+ fprintf(report->fp, "Nugget effect (hz): %f\n", var_par->horizontal.nugget);
+ fprintf(report->fp, "Sill (hz): %f\n", var_par->horizontal.sill);
+ fprintf(report->fp, "Range (hz): %f\n", var_par->horizontal.h_range);
+ fprintf(report->fp, "Function: %s\n\n", opt.function_var_hz->answer);
+ fprintf(report->fp, "Nugget effect (vert): %f\n", var_par->vertical.nugget);
+ fprintf(report->fp, "Sill (vert): %f\n", var_par->vertical.sill);
+ fprintf(report->fp, "Range (vert): %f\n", var_par->vertical.h_range);
+ fprintf(report->fp, "Function: %s\n", opt.function_var_vert->answer);
+ }
+ }
+
+ plot_var(TRUE, var_par->h, var_par->gamma, var_par); // Plot variogram using gnuplot
+ break;
+ case 3: // anisotropic
+ var_par->nugget = atof(opt.nugget_final->answer);
+ var_par->h_range = atof(opt.range_final->answer);
+ if (opt.sill_final->answer) {
+ var_par->sill = atof(opt.sill_final->answer);
+ }
+ variogram = opt.function_var_final->answer;
+ if (report->name) {
+ fprintf(report->fp, "Parameters of anisotropic variogram:\n");
+ fprintf(report->fp, "Nugget effect: %f\n", var_par->nugget);
+ fprintf(report->fp, "Sill: %f\n", var_par->sill);
+ fprintf(report->fp, "Range: %f\n", var_par->h_range);
+ }
+ break;
+ }
+
+ if (type != 2) {
+ var_par->function = set_function(variogram, var_par, report);
+ plot_var(FALSE, var_par->h, var_par->gamma, var_par); // Plot variogram using gnuplot
+ }
+}
+
+void ordinary_kriging(struct int_par *xD, struct reg_par *reg, struct points *pnts, struct pcl_utils *pclp, struct var_par *pars, struct output *out)
+{
+ // Local variables
+ int n = pnts->n;
+ double *r = pnts->r; // xyz coordinates of input poinst
+ double *vals = pnts->invals; // values to be used for interpolation
+ struct write *report = &xD->report;
+ struct write *crossvalid = &xD->crossvalid;
+ struct parameters *var_par;
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts;
+ double ratio;
+
+ var_par = &pars->fin;
+ pcl_pnts = pclp->pnts;
+ ratio = var_par->type == 3 ? xD->aniso_ratio : 1.;
+
+ unsigned int passed=0; // number of successfully interpolated valuesy
+ unsigned int percents=50; // counter
+ unsigned int nrcd; // number of cells/voxels
+ unsigned int row, col, dep; // indices of cells/voxels
+ double rslt_OK; // interpolated value located on r0
+
+ int i, j;
+ unsigned int n1;
+ double *r0; // xyz coordinates of cell/voxel centre
+ mat_struct *GM;
+ mat_struct *GM_sub; // submatrix of the points that are rellevant to the interpolation due the distance
+ mat_struct *GM_Inv; // inverted GM (GM_sub) matrix
+ mat_struct *g0; // vector of diffences between known and unknown values estimated using distances and the theoretical variogram
+ mat_struct *w0; // weights of values located on the input points
+
+ // Cell/voxel center coords (location of interpolated value)
+ r0 = (double *) G_malloc(3 * sizeof(double));
+ pcl::KdTreeFLANN<pcl::PointXYZ> kd_tree;
+ kd_tree.setInputCloud(pcl_pnts); // Set up kd-tree
+ pcl::PointXYZ cellCent;
+ std::vector<int> ind;
+ std::vector<float> sqDist;
+
+ double radius;
+ FILE *fp;
+
+ if (report->name) {
+ report->write2file = TRUE;
+ report->fp = fopen(report->name, "a");
+ time(&report->now);
+ fprintf(report->fp, "Interpolating values started on %s\n\n", ctime(&report->now));
+ }
+
+ G_message(_("Interpolating unknown values..."));
+ if (percents) {
+ G_percent_reset();
+ }
+
+ open_layer(xD, reg, out); // open 2D/3D raster
+
+ if (var_par->const_val == 1) {
+ goto constant_voxel_centre;
+ }
+
+ radius = var_par->max_dist;
+
+ GM = set_up_G(pnts, var_par, &xD->report);
+ var_par->GM = G_matrix_copy(GM);
+ if (n < 1000) {
+ GM_Inv = G_matrix_inverse(GM);
+ }
+
+ if (crossvalid->name) {
+ crossvalidation(xD, pnts, pcl_pnts, var_par);
+ }
+
+ constant_voxel_centre:
+ for (dep=0; dep < reg->ndeps; dep++) {
+ if (xD->i3 == TRUE) {
+ if (percents) {
+ G_percent(dep, reg->ndeps, 1);
+ }
+ }
+ for (row=0; row < reg->nrows; row++) {
+ if (xD->i3 == FALSE) {
+ if (percents) {
+ G_percent(row, reg->nrows, 1);
+ }
+ }
+ //#pragma omp parallel for private(col, r0, cellCent, ind, sqDist, n1, GM, GM_Inv, g0, w0, rslt_OK)
+ for (col=0; col < reg->ncols; col++) {
+
+ if (var_par->const_val == 1) {
+ goto constant_voxel_val;
+ }
+
+ cell_centre(col, row, dep, xD, reg, r0, var_par); // Output point
+ if (n < 1000)
+ goto small_sample;
+
+ cellCent.x = *r0;
+ cellCent.y = *(r0+1);
+ cellCent.z = *(r0+2); // if 2D: 0.
+
+ if (kd_tree.radiusSearch(cellCent, radius, ind, sqDist) > 0 ) {
+ GM_sub = submatrix(ind, GM, report);
+ GM_Inv = G_matrix_inverse(GM_sub);
+ G_matrix_free(GM_sub);
+
+ small_sample:
+ g0 = set_up_g0(xD, ind, pnts, r0, 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
+
+ if (n >= 1000) {
+ G_matrix_free(GM_Inv);
+ G_matrix_free(g0);
+ }
+
+ rslt_OK = result(ind, pnts, w0); // Estimated cell/voxel value rslt_OK = w x inputs
+ //if (pnts->trend != NULL) {
+ /* normal gravity: rslt_OK += pnts->trend->vals[0] * r0[0] + pnts->trend->vals[1] * r0[0]*r0[1]*r0[2]/ratio + pnts->trend->vals[2] * r0[0]*r0[1] + pnts->trend->vals[3] * r0[0]*r0[2]/ratio + pnts->trend->vals[4] * r0[1]*r0[2]/ratio + pnts->trend->vals[5] * r0[2]/ratio + pnts->trend->vals[6];*/
+ //}
+
+ // Create output
+ constant_voxel_val:
+ if (xD->const_val == 1)
+ rslt_OK = (double) *vals;
+ if (write2layer(xD, reg, out, col, row, dep, rslt_OK) == 0) {
+ if (report->name) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Error writing result into output layer..."));
+ }
+ } // end if searchRadius
+ else {
+ if (report->name) { // close report file
+ 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 col
+ } // end row
+ } // end dep
+
+ if (report->name) {
+ fprintf(report->fp, "\n************************************************\n\n");
+ time(&report->now);
+ fprintf(report->fp, "v.kriging completed on %s", ctime(&report->now));
+ fclose(report->fp);
+ }
+
+ switch (xD->i3) {
+ case TRUE:
+ Rast3d_close(out->fd_3d); // Close 3D raster map
+ break;
+ case FALSE:
+ Rast_close(out->fd_2d); // Close 2D raster map
+ break;
+ }
+}
+
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/getval.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/getval.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/getval.cpp 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,551 @@
+#include "local_proto.h"
+
+#ifndef MAX
+#define MIN(a,b) ((a<b) ? a : b)
+#define MAX(a,b) ((a>b) ? a : b)
+#endif
+
+/* get array of values from attribute column (function based on part of v.buffer2 (Radim Blazek, Rosen Matev)) */
+double *get_col_values(struct Map_info *map, struct int_par *xD, struct points *pnts, int field, const char *column, int detrend)
+{
+ struct select *in_reg = &pnts->in_reg;
+
+ int i, n, nrec, type, ctype;
+ struct field_info *Fi;
+
+ dbCatValArray cvarr;
+ dbDriver *Driver;
+
+ int *index;
+ double *values, *vals;
+
+ int save;
+ if (xD == NULL)
+ save = 0;
+ else
+ save = 1;
+
+ G_message(_("Reading values from attribute column <%s>..."), column);
+ db_CatValArray_init(&cvarr); /* array of categories and values initialised */
+
+ Fi = Vect_get_field(map, field); /* info about call of DB */
+ if (Fi == NULL) {
+ if (save == 1 && xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Database connection not defined for layer %d"), field);
+ }
+ Driver = db_start_driver_open_database(Fi->driver, Fi->database); /* started connection to DB */
+ if (Driver == NULL) {
+ if (save == 1 && xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
+ Fi->database, Fi->driver);
+ }
+
+ /* Note do not check if the column exists in the table because it may be expression */
+
+ /* TODO: only select values we need instead of all in column */
+
+ /* Select pairs key/value to array, values are sorted by key (must be integer) */
+ nrec = db_select_CatValArray(Driver, Fi->table, Fi->key, column, NULL, &cvarr);
+ if (nrec < 0) {
+ if (save == 1 && xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Unable to select data from table <%s>"), Fi->table);
+ }
+
+ ctype = cvarr.ctype;
+ if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE) {
+ if (save == 1 && xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Column must be numeric"));
+ }
+
+ db_close_database_shutdown_driver(Driver);
+
+ /* Output cats/values list */
+ switch (in_reg->out) {
+ case 0:
+ index = NULL;
+ n = cvarr.n_values;
+ break;
+ default:
+ index = &in_reg->indices[0];
+ n = in_reg->n;
+ break;
+ }
+
+ values = (double *) G_malloc(n * sizeof(double));
+ vals = &values[0];
+
+ for (i = 0; i < n; i++) {
+ /* TODO: Only for point data:
+ * - store indexes of skipped points in read_points and compare them with database
+ */
+ if (index == NULL || (index != NULL && *index == i)) {
+ if (ctype == DB_C_TYPE_INT)
+ *vals = (double) cvarr.value[i].val.i;
+ else if (ctype == DB_C_TYPE_DOUBLE) {
+ *vals = cvarr.value[i].val.d;
+ //if (i<78)
+ // G_debug(0, "%d %f", i, *vals);
+ }
+ vals++;
+ if (index != NULL)
+ index++;
+ }
+ }
+
+ if (detrend == 1) {
+ double *x, *y, *z;
+ x = (double *) G_malloc(n * sizeof(double));
+ y = (double *) G_malloc(n * sizeof(double));
+ z = (double *) G_malloc(n * sizeof(double));
+
+ mat_struct *A, *gR, *T, *res;
+ x = &pnts->r[0];
+ y = &pnts->r[1];
+ z = &pnts->r[2];
+ vals = &values[0];
+ // Number of columns of design matrix A
+ A = G_matrix_init(n, 2, n); // initialise design matrix, normal grav: 7
+ gR = G_matrix_init(n, 1, n); // initialise vector of observations
+
+ for (i = 0; i < n; i++) {
+ /* normal gravity:
+ G_matrix_set_element(A, i, 0, *x);
+ G_matrix_set_element(A, i, 1, *x * *y * *z);
+ G_matrix_set_element(A, i, 2, *x * *y);
+ G_matrix_set_element(A, i, 3, *x * *z);
+ G_matrix_set_element(A, i, 4, *z * *y);
+ G_matrix_set_element(A, i, 5, *z);
+ G_matrix_set_element(A, i, 6, 1.);
+ */
+ G_matrix_set_element(A, i, 0, *x);
+ G_matrix_set_element(A, i, 1, *y);
+ G_matrix_set_element(gR, i, 0, *vals);
+ x+=3;
+ y+=3;
+ z+=3;
+ vals++;
+ }
+ T = LSM(A, gR); // Least Square Method
+ res = G_matrix_product(A,T);
+
+ FILE *fp;
+ x = &pnts->r[0];
+ y = &pnts->r[1];
+ z = &pnts->r[2];
+ fp = fopen("trend.txt","w");
+
+ vals = &values[0];
+ double *resid;
+ resid = &res->vals[0];
+ for (i = 0; i < n; i++) {
+ *vals = *vals - *resid;
+ fprintf(fp,"%f %f %f %f\n", *x, *y, *z, *vals);
+ x+=3;
+ y+=3;
+ z+=3;
+ vals++;
+ resid++;
+ }
+ fclose(fp);
+ //G_debug(0,"a=%f b=%f c=%f d=%f", T->vals[0], T->vals[1], T->vals[2], T->vals[3]);
+ pnts->trend = T;
+ }
+
+ if (xD->phase == 0 && xD->report.name) {
+ write2file_values(&xD->report, column);
+ test_normality(n, values, &xD->report);
+ }
+
+ return values;
+}
+
+/* get coordinates of input points */
+void read_points(struct Map_info *map, struct reg_par *reg, struct points *point, struct pcl_utils *pclp, struct int_par xD, const char *zcol, int field, struct write *report)
+{
+ int type; // check elements of vector layer
+ int n; // # of all elements in vector layer
+ int nskipped; // # of skipped elements (all except points)
+ int n_in_reg = 0; // # of points within a region
+ int out_reg = 0; // # of points out of region
+ double *z_attr; // vertical coordinate from attribute value
+
+ struct line_pnts *Points; // structures to hold line *Points (map)
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts (new pcl::PointCloud<pcl::PointXYZ>);
+ pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree (new pcl::KdTreeFLANN<pcl::PointXYZ>);
+
+ double *rx, *ry, *rz, *r; // pointers to coordinates
+ pcl::PointXYZ *pclr; // pointer to PCL structure
+ int ind = 0; // index of the point
+ int *indices; // indices of the points within the region
+ int *index; // pointer to the vector of indices
+
+ G_message(_("Reading coordinates..."));
+ Points = Vect_new_line_struct(); // Make point structure
+ n = Vect_get_num_primitives(map, GV_POINT); // Number of points - topology required
+ r = (double *) G_malloc(n * 3 * sizeof(double)); // x0 y0 z0 ... xn yn zn
+ indices = (int *) G_malloc(n * sizeof(int));
+ point->r_min = (double *) G_malloc(3 * sizeof(double)); // minimum
+ point->r_max = (double *) G_malloc(3 * sizeof(double)); // maximum
+
+ rx = &r[0]; // set up pointer to x coordinate
+ ry = &r[1]; // set up pointer to y coordinate
+ rz = &r[2]; // set up pointer to z coordinate
+ index = &indices[0]; // set up pointer to indices
+
+ pcl_pnts->points.resize (n); // set up size of PCL structure
+ pclr = &pcl_pnts->points[0]; // set up pointer to PCL structure
+
+ // Get 3rd coordinate of 2D points from attribute column -> 3D interpolation
+ if (xD.v3 == FALSE && zcol != NULL)
+ z_attr = get_col_values(map, NULL, point, field, zcol, FALSE);
+
+ nskipped = 0; // # of skipped elements (lines, areas)
+
+ // Read points and set them into the structures
+ while (TRUE) {
+ type = Vect_read_next_line(map, Points, NULL);
+ if (type == -1) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to read vector map"));
+ }
+ if (type == -2)
+ break;
+
+ // skip everything except points
+ if (type != GV_POINT) {
+ nskipped++;
+ continue;
+ }
+
+ if (isnan(Points->x[0]) || isnan(Points->y[0]) || isnan(Points->z[0])) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Input coordinates must be a number. Check point no. %d..."), ind);
+ }
+
+ // take every point in region...
+ if ((reg->west <= Points->x[0] && Points->x[0] <= reg->east) && (reg->south <= Points->y[0] && Points->y[0] <= reg->north)) {
+ *rx = Points->x[0]; // set up x coordinate
+ *ry = Points->y[0]; // set up y coordinate
+
+ pclr->x = Points->x[0]; // set up x coordinate to PCL structure
+ pclr->y = Points->y[0]; // set up y coordinate to PCL structure
+
+ if (xD.i3 == TRUE) { // 3D interpolation
+ if (reg->bot <= Points->z[0] && Points->z[0] <= reg->top) {
+ if (zcol == NULL) { // 3D layer
+ *rz = Points->z[0]; // set up z coordinate
+ //if (n_in_reg < 79)
+ // G_debug(0,"%d %f %f %f", n_in_reg, *rx, *ry, *rz);
+ pclr->z = Points->z[0]; // set up z coordinate to PCL structure
+ if (*rz != Points->z[0]) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Error reading input coordinates z..."));
+ }
+ }
+ else { // 2D layer with z-column
+ *rz = *z_attr; // set up x coordinate
+ pclr->z = *z_attr; // set up x coordinate to PCL structure
+ } // end else (2,5D layer)
+ } // end if rb, rt
+ else
+ goto out_of_region;
+ } // end if 3D interpolation
+
+ else { // 2D interpolation
+ *rz = 0.; // set up z coordinate
+ pclr->z = 0.; // set up z coordinate to PCL structure
+ }
+
+ // Find extends
+ if (n_in_reg == 0)
+ point->r_min = point->r_max = triple(*rx, *ry, *rz);
+
+ else {
+ point->r_min = triple(MIN(*point->r_min, *rx), MIN(*(point->r_min+1), *ry), MIN(*(point->r_min+2), *rz));
+ point->r_max = triple(MAX(*point->r_max, *rx), MAX(*(point->r_max+1), *ry), MAX(*(point->r_max+2), *rz));
+ } // end else (extent test)
+
+ n_in_reg++; // # of points in region
+ *index = ind; // store index of point within the region
+ index++; // go to next index
+
+ rx += 3; // go to new x coordinate
+ ry += 3; // go to new y coordinate
+ rz += 3; // go to new z coordinate
+ pclr++; // go to new PCL point
+ goto finish;
+ } // end in region
+
+ out_of_region:
+ out_reg++;
+ continue;
+
+ finish:
+ if (zcol != NULL) // todo: skontroluj, co robi pre 3d
+ z_attr++;
+ ind++;
+ } // end while (TRUE)
+
+ point->in_reg.total = n;
+ point->in_reg.n = n_in_reg;
+ point->in_reg.out = out_reg;
+ point->n = n_in_reg;
+
+ point->r = (double *) G_malloc(3 * point->n * sizeof(double));
+ memcpy(point->r, r, 3 * n_in_reg * sizeof(double));
+
+ point->in_reg.indices = (int *) G_malloc(n_in_reg * sizeof(int));
+ memcpy(point->in_reg.indices, indices, n_in_reg * sizeof(int));
+
+ pcl_pnts->points.resize (n_in_reg);
+
+ pclp->pnts = pcl_pnts;
+ kd_tree->setInputCloud(pclp->pnts);
+ pclp->kd_tree = kd_tree;
+
+ if (nskipped > 0)
+ G_warning(_("%d features skipped, only points are accepted"),
+ nskipped);
+
+ Vect_destroy_line_struct(Points);
+
+ if (n_in_reg == 0) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to read coordinates of point layer"));
+ }
+
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts_hz (new pcl::PointCloud<pcl::PointXYZ>);
+ pcl_pnts_hz->points.resize (n);
+
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts_vert (new pcl::PointCloud<pcl::PointXYZ>);
+ pcl_pnts_vert->points.resize (n);
+
+ for (int i=0; i<n_in_reg; i++) {
+ pcl_pnts_hz->points[i].x = pcl_pnts->points[i].x;
+ pcl_pnts_hz->points[i].y = pcl_pnts->points[i].y;
+ pcl_pnts_hz->points[i].z = 0.;
+
+ pcl_pnts_vert->points[i].x = 0.;
+ pcl_pnts_vert->points[i].y = 0.;
+ pcl_pnts_vert->points[i].z = pcl_pnts->points[i].z;
+ } // end switch var_par->function
+
+ pclp->pnts_hz = pcl_pnts_hz;
+ pclp->pnts_vert = pcl_pnts_vert;
+
+ //write2file_vector(&xD, point);
+}
+
+extern "C" {
+ void get_region_pars(struct int_par *xD, struct reg_par *reg)
+ {
+ /* Get extent and resolution of output raster 2D/3D */
+ if (xD->i3 == FALSE) { /* 2D region */
+ /* Get database region parameters */
+ Rast_get_window(®->reg_2d); /* stores the current default window in region */
+ reg->west = reg->reg_2d.west;
+ reg->east = reg->reg_2d.east;
+ reg->north = reg->reg_2d.north;
+ reg->south = reg->reg_2d.south;
+ reg->ew_res = reg->reg_2d.ew_res;
+ reg->ns_res = reg->reg_2d.ns_res;
+ reg->nrows = reg->reg_2d.rows;
+ reg->ncols = reg->reg_2d.cols;
+ reg->ndeps = 1;
+ }
+ if (xD->i3 == TRUE) { /* 3D region */
+ Rast3d_init_defaults();
+ /* Get database region parameters */
+ Rast3d_get_window(®->reg_3d); /* stores the current default window in region */
+ if (reg->reg_3d.bottom == 0 && (reg->reg_3d.tb_res == 0 || reg->reg_3d.depths == 0)) {
+ if (xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error("To process 3D interpolation, please set 3D region settings.");
+ }
+ reg->west = reg->reg_3d.west;
+ reg->east = reg->reg_3d.east;
+ reg->north = reg->reg_3d.north;
+ reg->south = reg->reg_3d.south;
+ reg->bot = reg->reg_3d.bottom; /* output 3D raster */
+ reg->top = reg->reg_3d.top; /* output 3D raster */
+ reg->ew_res = reg->reg_3d.ew_res;
+ reg->ns_res = reg->reg_3d.ns_res;
+ reg->bt_res = reg->reg_3d.tb_res;
+ reg->nrows = reg->reg_3d.rows;
+ reg->ncols = reg->reg_3d.cols;
+ reg->ndeps = xD->i3 == TRUE ? reg->reg_3d.depths : 1; /* 2D/3D interpolation */
+ }
+ reg->nrcd = reg->nrows * reg->ncols * reg->ndeps; /* number of cells */
+ }
+}
+
+void read_tmp_vals(const char *file_name, struct parameters *var_par, struct int_par *xD)
+{
+ FILE *fp;
+
+ fp = fopen(file_name, "r");
+ if (fp == NULL)
+ G_fatal_error(_("File is missing..."));
+
+ else {
+ int i, type;
+ int nLag;
+ double lag, max_dist, td_hz, sill;
+ double *h, *h_elm, *c, *gamma;
+ int function;
+ int file, file_length;
+
+ for (i=0; i<2; i++) {
+ if (fscanf(fp, "%d", &file_length) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+ if (file_length > 3) {
+ if (fscanf(fp, "%d", &file) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+ if (file == 9) {
+ xD->report.name = (char *) G_malloc(file_length * sizeof(char));
+ if (fscanf(fp, "%s", xD->report.name) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+ continue;
+ }
+ else if (file == 8) {
+ xD->crossvalid.name = (char *) G_malloc(file_length * sizeof(char));
+ if (fscanf(fp, "%s", xD->crossvalid.name) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+ continue;
+ }
+ }
+ else {
+ type = file_length;
+ goto no_file;
+ }
+ } // todo: test without report and crossval
+
+ no_file:
+
+ switch (type) {
+ case 2:
+ int j, nLag_vert;
+ double lag_vert, max_dist_vert;
+ double *v_elm, *g_elm, *c_elm;
+ double sill_hz, sill_vert;
+
+ if (fscanf(fp, "%d %lf %lf %d %lf %lf %lf", &nLag_vert, &lag_vert, &max_dist_vert, &nLag, &lag, &max_dist, &td_hz) < 7)
+ G_fatal_error(_("Nothing to scan..."));
+
+
+ var_par->h = (double *) G_malloc(nLag * sizeof(double));
+ h_elm = &var_par->h[0];
+ for (i=0; i<nLag; i++) {
+ if (fscanf(fp, "%lf", h_elm) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ h_elm++;
+ }
+
+ var_par->vert = (double *) G_malloc(nLag_vert * sizeof(double));
+ v_elm = &var_par->vert[0];
+ for (i=0; i<nLag_vert; i++) {
+ if (fscanf(fp, "%lf", v_elm) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ v_elm++;
+ }
+
+ var_par->c = G_matrix_init(nLag, nLag_vert, nLag);
+ c_elm = &var_par->c->vals[0];
+ for (i=0; i<nLag_vert; i++) {
+ for (j=0; j<nLag; j++) {
+ if (fscanf(fp, "%lf", c_elm) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+ c_elm++;
+ }
+ }
+
+ var_par->gamma = G_matrix_init(nLag, nLag_vert, nLag);
+ g_elm = &var_par->gamma->vals[0];
+ for (i=0; i<nLag_vert; i++) {
+ for (j=0; j<nLag; j++) {
+ if (fscanf(fp, "%lf", g_elm) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+ g_elm++;
+ }
+ }
+
+ if (fscanf(fp, "%lf %lf", &sill_hz, &sill_vert) < 2)
+ G_fatal_error(_("Nothing to scan..."));
+
+ var_par->nLag_vert = nLag_vert;
+ var_par->lag_vert = lag_vert;
+ var_par->max_dist_vert = max_dist_vert;
+ var_par->horizontal.sill = sill_hz;
+ var_par->vertical.sill = sill_vert;
+ break;
+ default:
+ if (type == 3) {
+ if (fscanf(fp, "%lf", &xD->aniso_ratio) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ if (fscanf(fp, "%d %lf %lf", &nLag, &lag, &max_dist) < 3) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ if (type != 1)
+ if (fscanf(fp, "%lf", &td_hz) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+
+ var_par->h = (double *) G_malloc(nLag * sizeof(double));
+ h_elm = &var_par->h[0];
+
+ var_par->c = G_matrix_init(nLag, 1, nLag);
+ c = &var_par->c->vals[0];
+
+ var_par->gamma = G_matrix_init(nLag, 1, nLag);
+ gamma = &var_par->gamma->vals[0];
+
+ for (i=0; i<nLag; i++) {
+ if (fscanf(fp, "%lf %lf %lf", h_elm, c, gamma) < 3)
+ G_fatal_error(_("Nothing to scan..."));
+ h_elm++;
+ c++;
+ gamma++;
+ }
+
+ if (fscanf(fp, "%lf", &sill) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+ var_par->sill = sill;
+ break;
+ }
+
+ var_par->type = type;
+ var_par->nLag = nLag;
+ var_par->lag = lag;
+ var_par->max_dist = max_dist;
+ var_par->td = td_hz;
+ }
+ fclose(fp);
+}
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/local_proto.h
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/local_proto.h (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/local_proto.h 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,264 @@
+#ifndef __LOCAL_PROTO__
+#define __LOCAL_PROTO__
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/config.h>
+#include <grass/gis.h>
+#include <grass/gmath.h>
+#include <grass/la.h>
+#include <grass/raster3d.h>
+#include <grass/glocale.h>
+#include <grass/rtree.h>
+
+#ifndef PI
+#define PI M_PI
+#define DEG2RAD(ang) (ang / 180. * PI)
+#define RAD2DEG(ang) (ang / PI * 180.)
+#endif
+
+#ifndef SQUARE
+#define SQUARE(a) (a*a)
+#define POW3(a) (a*a*a)
+#define POW4(a) (a*a*a*a)
+#endif
+
+struct opts
+{
+ struct Option *input, *output, *phase, *report, *crossvalid, *function_var_hz, *function_var_vert, *function_var_final, *function_var_final_vert, *form_file, *field, *intpl, *zcol, *var_dir_hz, *var_dir_vert, *nL, *nZ, *td_hz, *td_vert, *nugget_hz, *nugget_vert, *nugget_final, *nugget_final_vert, *sill_hz, *sill_vert, *sill_final, *sill_final_vert, *range_hz, *range_vert, *range_final, *range_final_vert;
+};
+
+struct flgs
+{
+ struct Flag *d23; /* 2D/3D interpolation */
+ struct Flag *bivariate;
+ struct Flag *univariate;
+ struct Flag *detrend;
+};
+
+struct select
+{
+ int n; // # of selected
+ int out; // # of the others
+ int total; // total number
+ int *indices; // indices of selected
+};
+
+struct points // inputs
+{
+ int n; // number of points
+ double *r; // triples of coordinates (e.g. x0 y0 z0... xn yn zn)
+ struct RTree *R_tree; // spatial index
+ struct RTree *Rtree_hz; // spatial index
+ struct RTree *Rtree_vert; // spatial index
+ double *r_min; // min coords
+ double *r_max; // max coords
+ double max_dist; // maximum distance
+ double center[3]; // center
+ double *invals; // values to be interpolated
+ struct select in_reg; // points in region
+ mat_struct *trend;
+};
+
+struct bivar
+{
+ int vert;
+ char *variogram;
+
+ double *h;
+ double max_dist;
+ int nLag;
+ double lag;
+
+ int function;
+ double sill;
+ double nugget;
+ double h_range;
+};
+
+struct parameters
+{
+ int function; // variogram function: lin, exp, spher, Gauss, bivar
+ int type; // variogram type: hz / vert / aniso / bivar
+ int const_val;
+ double dir; // azimuth for variogram computing
+ double td; // maximum azimuth
+
+ double radius; // radius (squared maximum distance)
+ double max_dist; // maximum distance - hz / aniso
+ double max_dist_vert; // maximum distance - vert
+
+ int nLag; // number of lags - hz / aniso
+ double lag; // lag size
+
+ int nLag_vert; // number of lags - vert
+ double lag_vert; // lag size
+
+ double *h; // lag distance from search point - hz / aniso
+ double *vert; // lag distance from search point - vert
+
+ int gamma_n; // # of dissimilarities between input points
+ mat_struct *gamma; // experimental variogram matrix
+ double gamma_sum; // sum of gamma values
+
+ double nugget;
+ double sill;
+ double part_sill;
+ double h_range;
+
+ struct bivar horizontal; // horizontal variogram properties
+ struct bivar vertical; // vertical variogram properties
+
+ mat_struct *A; // plan matrix
+ mat_struct *T; // coefficients of theoretical variogram
+ mat_struct *GM; // GM = theor_var(dist: input, output points)
+
+ char *name; // name of input vector layer
+ char term[12]; // output format - gnuplot terminal
+ char ext[4]; // output format - file format
+};
+
+struct var_par // parameters of experimental variogram
+{
+ struct parameters hz;
+ struct parameters vert;
+ struct parameters fin;
+
+ //double Lmax; // maximal length in horizontal direction
+ //double dzmax; // maximal elevation difference
+ //double radius; // radius for kd-tree
+};
+
+struct write
+{
+ int write2file; // write to file or not
+ char *name; // filename
+ FILE *fp;
+ time_t now;
+};
+
+struct int_par // Interpolation settings
+{
+ int i3; // TRUE = 3D interpolation, FALSE = 2D interpolation (user sets by flag)
+ char v3; // TRUE = 3D layer, FALSE = 2D layer (according to Vect_is_3d())
+ int phase;
+ int bivar;
+ int univar;
+ double aniso_ratio;
+ int const_val;
+ struct write report;
+ struct write crossvalid;
+};
+
+struct reg_par // Region settings -> output extent and resolution
+{
+ struct Cell_head reg_2d; // region for 2D interpolation
+ RASTER3D_Region reg_3d; // region for 3D interpolation
+ double west; // region.west
+ double east; // region.east
+ double north; // region.north
+ double south; // region.south
+ double bot; // region.bottom
+ double top; // region.top
+ double ew_res; // east-west resolution
+ double ns_res; // north-south resolution
+ double bt_res; // bottom-top resolution
+ int nrows; // # of rows
+ int ncols; // # of cols
+ int ndeps; // # of deps
+ int nrcd; // # of cells
+};
+
+struct output
+{
+ DCELL *dcell;
+ int fd_2d;
+ RASTER3D_Map *fd_3d;
+ char *name; // name of output 2D/3D raster
+};
+
+double *get_col_values(struct Map_info *, struct int_par *, struct points *, int, const char *, int);
+void test_normality(int, double *, struct write *);
+void read_points(struct Map_info *, struct reg_par *, struct points *, struct int_par *, const char *, int, struct write *);
+double min(double *, struct points *);
+double max(double *, struct points *);
+
+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_n_NNs(int, int, struct points *, int);
+double sum_NN(int, int, struct ilist *, struct points *);
+
+void correct_indices(int, struct ilist *, double *, struct points *, struct parameters *);
+int cmpVals(const void *, const void *);
+void coord_diff(int, int, double *, double *);
+double distance_diff(double *);
+double radius_hz_diff(double *);
+double zenith_angle(double *);
+void triple(double, double, double, double *);
+double lag_size(int, struct int_par *, struct points *, struct parameters *, struct write *);
+int lag_number(double, double *);
+void variogram_restricts(struct int_par *, struct points *, struct parameters *);
+void geometric_anisotropy(struct int_par *, struct points *);
+double find_intersect_x(double *, double *, double *, double *, struct write *);
+double find_intersect_y(double *, double *, double *, double *, double , struct write *);
+mat_struct *LSM(mat_struct *, mat_struct *);
+mat_struct *nonlin_LMS(int , double *, double *);
+
+void E_variogram(int, struct int_par *, struct points *, struct reg_par *, struct var_par *);
+void T_variogram(int, int, struct opts, struct parameters *, struct write *);
+void ordinary_kriging(struct int_par *, struct reg_par *, struct points *, struct var_par *, struct output *);
+
+void linear_variogram(struct parameters *, struct write *);
+double bivar_sill(int, mat_struct *);
+void sill(struct parameters *);
+int sill_compare(struct int_par *, struct flgs *, struct var_par *, struct points *);
+int set_function(char *, struct write *);
+double RBF(double *);
+double linear(double, double, double);
+double exponential(double, double, double, double);
+double spherical(double, double, double, double);
+double gaussian(double, double, double, double);
+double variogram_fction(struct parameters *, double *);
+void set_gnuplot(char *, struct parameters *);
+void plot_experimental_variogram(struct int_par *, struct parameters *);
+void plot_var(int, int, struct parameters *);
+void variogram_type(int, char *);
+void write2file_basics(struct int_par *, struct opts *);
+void write2file_vector(struct int_par *, struct points *);
+void write2file_values(struct write *, const char *);
+void write2file_varSetsIntro(int, struct write *);
+void write2file_varSets(struct write *, struct parameters *);
+void write2file_variogram_E(struct int_par *, struct parameters *, mat_struct *);
+void write2file_variogram_T(struct write *);
+void write_temporary2file(struct int_par *, struct parameters *);
+void read_tmp_vals(const char *, struct parameters *, struct int_par *);
+
+mat_struct *set_up_G(struct points *, struct parameters *, struct write *);
+mat_struct *set_up_g0(struct int_par *, struct points *, struct ilist *, double *, struct parameters *);
+mat_struct *submatrix(struct ilist *, mat_struct *, struct write *);
+double result(struct points *, struct ilist *, mat_struct *);
+
+void crossvalidation(struct int_par *, struct points *, struct parameters *);
+void cell_centre(unsigned int, unsigned int, unsigned int, struct int_par *, struct reg_par *, double *, struct parameters *);
+
+void get_region_pars(struct int_par *, struct reg_par *);
+void open_layer(struct int_par *, struct reg_par *, struct output *);
+int write2layer(struct int_par *, struct reg_par *, struct output *, unsigned int, unsigned int, unsigned int, double);
+
+static inline double get_quantile(int);
+static void get_slot_counts(int, double *);
+static void initialize_bins(void);
+static void fill_bins(int, double *);
+static int compare(const void *, const void *);
+static void sort_bins(void);
+static void compute_quantiles(int, double *, struct write *);
+double quantile(double, int, double *, struct write *);
+#endif
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/main.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/main.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/main.cpp 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,496 @@
+
+/****************************************************************
+ *
+ * MODULE: v.kriging
+ * AUTHOR: Eva Stopková
+ * in case of functions taken from another modules,
+ * they are cited above the function
+ * or at the beginning of the file (e.g. quantile.cpp
+ * that uses slightly modified functions
+ * taken from the module r.quantile (Clemens, G.))
+ * PURPOSE: Module interpolates the values to two- or three-dimensional grid using input values
+ * located on 2D/3D point vector layer. Interpolation method
+ * Ordinary kriging has been extended for 3D points (v = f(x,y) -> v = f(x,y,z)).
+ *
+ * COPYRIGHT: (C) 2012-2014 Eva Stopková and by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that
+ * comes with GRASS for details.
+ *
+ **************************************************************/
+#include "local_proto.h"
+
+
+int main(int argc, char *argv[])
+{
+ // Vector layer and module
+ struct Map_info map; // Input vector map
+ struct GModule *module; // Module
+
+ struct reg_par reg; // Region parameters
+ struct points pnts; // Input points (coordinates, extent, values, etc.)
+ struct pcl_utils pclp; // PCL structure to store coordinates of input points
+
+ // Geostatistical parameters
+ struct int_par xD; // 2D/3D interpolation for 2D/3D vector layer
+ struct var_par var_par; // Variogram (experimental and theoretical)
+
+ // Outputs
+ struct output out; // Output layer properties
+ FILE *fp;
+
+ // Settings
+ int field;
+ struct opts opt;
+ struct flgs flg;
+
+ /* ------- Module creation ------- */
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("3D raster"));
+ G_add_keyword(_("ordinary kriging - for 2D and 3D data"));
+ module->description =
+ _("Interpolates 2D or 3D raster based on input values located on 2D or 3D point vector layer (method ordinary kriging extended to 3D).");
+
+ // Setting options
+ opt.input = G_define_standard_option(G_OPT_V_INPUT); // Vector input layer
+ opt.input->label = _("Name of input vector points map");
+
+ flg.d23 = G_define_flag(); // Option to process 2D or 3D interpolation
+ flg.d23->key = '2';
+ flg.d23->description = _("Force 2D interpolation even if input is 3D");
+ flg.d23->guisection = _("3D");
+
+ opt.field = G_define_standard_option(G_OPT_V_FIELD);
+
+ opt.phase = G_define_option();
+ opt.phase->key = "phase";
+ opt.phase->options = "initial, middle, final";
+ opt.phase->description = _("Phase of interpolation. In the initial phase, there is empirical variogram computed. In the middle phase, function of theoretical variogram is chosen by the user and its coefficients are estimated empirically. In the final phase, unknown values are interpolated using theoretical variogram from previous phase.");
+ opt.phase->required = YES;
+
+ opt.output = G_define_option(); // Output layer
+ opt.output->key = "output";
+ opt.output->description =
+ _("Name for output 2D/3D raster map");
+ opt.output->guisection = _("Final");
+
+ opt.report = G_define_standard_option(G_OPT_F_OUTPUT); // Report file
+ opt.report->key = "report";
+ opt.report->description = _("Path and name of a file where the report should be written (initial phase only)");
+ opt.report->required = NO;
+ opt.report->guisection = _("Initial");
+
+ opt.crossvalid = G_define_standard_option(G_OPT_F_OUTPUT); // Report file
+ opt.crossvalid->key = "crossvalid";
+ opt.crossvalid->description = _("Path and name of a file where the results of cross validation should be written (final phase only)");
+ opt.crossvalid->required = NO;
+ opt.crossvalid->guisection = _("Final");
+
+ flg.bivariate = G_define_flag();
+ flg.bivariate->key = 'b';
+ flg.bivariate->description = _("Compute bivariate variogram (just in case of 3D interpolation)");
+ flg.bivariate->guisection = _("Middle");
+
+ flg.univariate = G_define_flag();
+ flg.univariate->key = 'u';
+ flg.univariate->description = _("Compute univariate variogram (just in case of 3D interpolation)");
+ flg.univariate->guisection = _("Middle");
+
+ opt.function_var_hz = G_define_option(); // Variogram type
+ opt.function_var_hz->key = "hz_function";
+ opt.function_var_hz->options = "linear, exponential, spherical, gaussian, bivariate";
+ opt.function_var_hz->description = _("Type of horizontal variogram function (just in middle phase)");
+ opt.function_var_hz->guisection = _("Middle");
+
+ opt.function_var_vert = G_define_option(); // Variogram type
+ opt.function_var_vert->key = "vert_function";
+ opt.function_var_vert->options = "linear, exponential, spherical, gaussian, bivariate";
+ opt.function_var_vert->description = _("Type of vertical variogram function (just in middle phase)");
+ opt.function_var_vert->guisection = _("Middle");
+
+ opt.function_var_final = G_define_option(); // Variogram type
+ opt.function_var_final->key = "final_function";
+ opt.function_var_final->options = "linear, exponential, spherical, gaussian, bivariate";
+ opt.function_var_final->description = _("Type of final variogram function (just in final phase)");
+ opt.function_var_final->guisection = _("Final");
+
+ flg.detrend = G_define_flag();
+ flg.detrend->key = 't';
+ flg.detrend->description = _("Eliminate trend if variogram is parabolic (just in initial phase)");
+ flg.detrend->guisection = _("Initial");
+
+ opt.form_file = G_define_option(); // Variogram plot - various output formats
+ opt.form_file->key = "fileformat";
+ opt.form_file->options = "cdr,dxf,eps,tex,pdf,png,svg";
+ opt.form_file->description = _("Variogram plot (set up in the middle phase) can be saved in several file formats or just presented in Gnuplot terminal (if fileformat is not set)");
+ opt.form_file->guisection = _("Middle");
+
+ opt.intpl = G_define_standard_option(G_OPT_DB_COLUMN); // Input values for interpolation
+ opt.intpl->key = "icolumn";
+ opt.intpl->description =
+ _("Name of the attribute column containing input values for interpolation");
+ opt.intpl->required = YES;
+
+ opt.zcol = G_define_standard_option(G_OPT_DB_COLUMN); // Column with z coord (2D points)
+ opt.zcol->key = "zcolumn";
+ opt.zcol->description =
+ _("Name of the attribute column containing z coordinates - set only if you want to perform 3D interpolation based on 2D point layer");
+ opt.zcol->required = NO;
+ opt.zcol->guisection = _("3D");
+
+ opt.var_dir_hz = G_define_option();
+ opt.var_dir_hz->key = "azimuth";
+ opt.var_dir_hz->type = TYPE_DOUBLE;
+ opt.var_dir_hz->required = NO;
+ opt.var_dir_hz->answer = "0.0";
+ opt.var_dir_hz->description =
+ _("Azimuth of variogram computing (isotrophic), initial phase only ");
+ opt.var_dir_hz->guisection = _("Initial");
+
+ opt.var_dir_vert = G_define_option();
+ opt.var_dir_vert->key = "zenith_angle";
+ opt.var_dir_vert->type = TYPE_DOUBLE;
+ opt.var_dir_vert->required = NO;
+ opt.var_dir_vert->answer = "0.0";
+ opt.var_dir_vert->description =
+ _("Zenith angle of variogram computing (isotrophic), initial phase only");
+ opt.var_dir_vert->guisection = _("Initial");
+
+ opt.nL = G_define_option();
+ opt.nL->key = "lpieces";
+ opt.nL->type = TYPE_INTEGER;
+ opt.nL->required = NO;
+ opt.nL->description = _("Count of length pieces, initial phase only");
+ opt.nL->guisection = _("Initial");
+
+ opt.nZ = G_define_option();
+ opt.nZ->key = "vpieces";
+ opt.nZ->type = TYPE_INTEGER;
+ opt.nZ->required = NO;
+ opt.nZ->description =
+ _("Count of vertical pieces (set only for computing 3D variogram), initial phase only");
+ opt.nZ->guisection = _("Initial");
+
+ opt.td_hz = G_define_option();
+ opt.td_hz->key = "td";
+ opt.td_hz->type = TYPE_DOUBLE;
+ opt.td_hz->answer = "90.0";
+ opt.td_hz->description = _("Angle of variogram processing, initial phase only");
+ opt.td_hz->required = NO;
+ opt.td_hz->guisection = _("Initial");
+
+ opt.nugget_hz = G_define_option();
+ opt.nugget_hz->key = "hz_nugget";
+ opt.nugget_hz->type = TYPE_DOUBLE;
+ opt.nugget_hz->description = _("Nugget effect of horizontal variogram, middle phase only");
+ opt.nugget_hz->required = NO;
+ opt.nugget_hz->guisection = _("Middle");
+
+ opt.nugget_vert = G_define_option();
+ opt.nugget_vert->key = "vert_nugget";
+ opt.nugget_vert->type = TYPE_DOUBLE;
+ opt.nugget_vert->description = _("Nugget effect of vertical variogram, middle phase only");
+ opt.nugget_vert->required = NO;
+ opt.nugget_vert->guisection = _("Middle");
+
+ opt.nugget_final = G_define_option();
+ opt.nugget_final->key = "final_nugget";
+ opt.nugget_final->type = TYPE_DOUBLE;
+ opt.nugget_final->description = _("Nugget effect of final variogram, final phase only");
+ opt.nugget_final->required = NO;
+ opt.nugget_final->guisection = _("Final");
+
+ opt.sill_hz = G_define_option();
+ opt.sill_hz->key = "hz_sill";
+ opt.sill_hz->type = TYPE_DOUBLE;
+ opt.sill_hz->description = _("Sill of horizontal variogram, middle phase only");
+ opt.sill_hz->required = NO;
+ opt.sill_hz->guisection = _("Middle");
+
+ opt.sill_vert = G_define_option();
+ opt.sill_vert->key = "vert_sill";
+ opt.sill_vert->type = TYPE_DOUBLE;
+ opt.sill_vert->description = _("Sill of vertical variogram, middle phase only");
+ opt.sill_vert->required = NO;
+ opt.sill_vert->guisection = _("Middle");
+
+ opt.sill_final = G_define_option();
+ opt.sill_final->key = "final_sill";
+ opt.sill_final->type = TYPE_DOUBLE;
+ opt.sill_final->description = _("Sill of final variogram, final phase only");
+ opt.sill_final->required = NO;
+ opt.sill_final->guisection = _("Final");
+
+ opt.range_hz = G_define_option();
+ opt.range_hz->key = "hz_range";
+ opt.range_hz->type = TYPE_DOUBLE;
+ opt.range_hz->description = _("Range of horizontal variogram, middle phase only");
+ opt.range_hz->required = NO;
+ opt.range_hz->guisection = _("Middle");
+
+ opt.range_vert = G_define_option();
+ opt.range_vert->key = "vert_range";
+ opt.range_vert->type = TYPE_DOUBLE;
+ opt.range_vert->description = _("Range of vertical variogram, middle phase only");
+ opt.range_vert->required = NO;
+ opt.range_vert->guisection = _("Middle");
+
+ opt.range_final = G_define_option();
+ opt.range_final->key = "final_range";
+ opt.range_final->type = TYPE_DOUBLE;
+ opt.range_final->description = _("Range of final variogram, final phase only");
+ opt.range_final->required = NO;
+ opt.range_final->guisection = _("Final");
+ /* --------------------------------------------------------- */
+
+ G_gisinit(argv[0]);
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* Get parameters from the parser */
+ if (strcmp(opt.phase->answer, "initial") == 0)
+ xD.phase = 0; // estimate range
+ else if (strcmp(opt.phase->answer, "middle") == 0)
+ xD.phase = 1; // estimate anisotropic variogram
+ else if (strcmp(opt.phase->answer, "final") == 0)
+ xD.phase = 2; // compute kriging
+
+ // Open report file if desired
+ if (opt.report->answer) {
+ xD.report.write2file = TRUE;
+ xD.report.name = opt.report->answer;
+ xD.report.fp = fopen(xD.report.name, "w");
+ time(&xD.report.now);
+ fprintf(xD.report.fp, "v.kriging started on %s\n\n", ctime(&xD.report.now));
+ G_message(_("Report is being written to %s..."), xD.report.name);
+ }
+ else
+ xD.report.write2file = FALSE;
+
+ if (opt.crossvalid->answer) {
+ xD.crossvalid.write2file = TRUE;
+ xD.crossvalid.name = opt.crossvalid->answer;
+ xD.crossvalid.fp = fopen(xD.crossvalid.name, "w");
+ }
+ else
+ xD.crossvalid.write2file = FALSE;
+
+ var_par.hz.name = var_par.vert.name = var_par.fin.name = opt.input->answer; // set name of variogram
+ var_par.hz.dir = opt.var_dir_hz->answer ? DEG2RAD(atof(opt.var_dir_hz->answer)) : 0.; // Azimuth
+ var_par.vert.dir = opt.var_dir_vert->answer ? DEG2RAD(atof(opt.var_dir_vert->answer)) : 0.; // Zenith angle
+ var_par.fin.dir = opt.var_dir_hz->answer ? DEG2RAD(atof(opt.var_dir_hz->answer)) : 0.; // Azimuth
+
+ var_par.hz.td = DEG2RAD(atof(opt.td_hz->answer)); // Angle of variogram processing
+ //var_par.vert.td = DEG2RAD(atof(opt.td_vert->answer)); // Angle of variogram processing
+
+ /*if (phase > 0) {
+ //var_par.sill = atof(opt.sill->answer);
+ var_par.nugget = atof(opt.nugget->answer);
+ var_par.h_range = atof(opt.range->answer);
+ if (/*(var_par.sill != -1. && var_par.sill <= 0.) || *//*(var_par.nugget != -1. && var_par.nugget < 0.) || (var_par.h_range != -1. && var_par.h_range <= 0.)) {
+ if (xD.report.write2file == TRUE) {
+ fclose(xD.report.fp);
+ remove(xD.report.name);
+ }
+ G_fatal_error(_("Variogram parameters should be positive..."));
+ } // error
+ } */
+
+ /*if (opt.nL->answer) { // Test if nL have been set up (optional)
+ if (var_par.nL < 1) // Invalid value
+ G_message(_("Number of horizontal pieces must be at least 1. Default value will be used..."));
+ else
+ var_par.nL = atof(opt.nL->answer); // todo osetrit aby sa neprepisovalo
+ }
+
+ if (opt.variogram->answer) // Function of theoretical variogram
+ set_function(opt.variogram->answer, &var_par, &xD.report);*/
+
+ if (opt.form_file->answer) { // Plotting variogram
+ set_gnuplot(opt.form_file->answer, &var_par.hz);
+ set_gnuplot(opt.form_file->answer, &var_par.vert);
+ set_gnuplot(opt.form_file->answer, &var_par.fin);
+ }
+ else {
+ strcpy(var_par.hz.term, ""); /* TO DO: ked pouzivatel zada hlupost */
+ strcpy(var_par.vert.term, "");
+ strcpy(var_par.fin.term, "");
+ }
+
+ xD.bivar = flg.bivariate->answer == TRUE ? TRUE : FALSE;
+ xD.univar = flg.univariate->answer == TRUE ? TRUE : FALSE;
+ if (xD.bivar == TRUE && xD.univar == TRUE) {
+ if (xD.report.write2file == TRUE) {
+ fclose(xD.report.fp);
+ remove(xD.report.name);
+ }
+ G_fatal_error(_("You should mark either univariate, or bivariate variogram, not both of them..."));
+ } // error
+
+ /* ---------------------------------------------------------- */
+ Vect_set_open_level(2); // Open input vector map
+
+ if (0 > Vect_open_old2(&map, opt.input->answer, "", opt.field->answer)) {
+ if (xD.report.write2file == TRUE) {
+ fclose(xD.report.fp);
+ remove(xD.report.name);
+ }
+ G_fatal_error(_("Unable to open vector map <%s>"), opt.input->answer);
+ } // error
+ Vect_set_error_handler_io(&map, NULL);
+ /* ---------------------------------------------------------- */
+
+ /* Perform 2D or 3D interpolation ? */
+ xD.i3 = flg.d23->answer ? FALSE : TRUE; // 2D/3D interpolation
+ xD.v3 = Vect_is_3d(&map); // 2D/3D layer
+
+ /* What could happen:
+ *-------------------
+ * 3D interpolation + 3D points = 3D GIS
+ * 3D interpolation + 2D points = 2,5D -> 3D GIS (needs attribute column with z and 3D region)
+ * 2D interpolation + 3D points = 3D -> 2,5D GIS
+ * 2D interpolation + 2D points = 2,5D GIS */
+
+ // 3D interpolation
+ if (xD.i3 == TRUE) {
+ if (xD.v3 == FALSE) { // 2D input
+ if (!opt.zcol->answer) {
+ if (xD.report.write2file == TRUE) { // close report file
+ fprintf(xD.report.fp, "Error (see standard output). Process killed...");
+ fclose(xD.report.fp);
+ }
+ G_fatal_error(_("To process 3D interpolation based on 2D input, please set attribute column containing z coordinates or switch to 2D interpolation."));
+ }
+ } // end if zcol == NULL
+ // 3D or 2,5D input
+ if (opt.nZ->answer) { // Test if nZ have been set up (optional)
+ if (var_par.vert.nLag < 1) // Invalid value
+ G_message(_("Number of vertical pieces must be at least 1. Default value will be used..."));
+ else
+ var_par.vert.nLag = atof(opt.nZ->answer);
+ } // end if nZ->answer
+ /* TO DO: chyba, ak nie je cele cislo */
+ } // end if 3D interpolation
+
+ /* 2D interpolation */
+ else {
+ G_fatal_error(_("Access to 2D kriging is denied temporarily. Please feel free to use 3D kriging..."));
+ var_par.vert.nLag = -1; // abs will be used in next steps
+ }
+
+ field = Vect_get_field_number(&map, opt.field->answer);
+ if (xD.report.write2file == TRUE)
+ write2file_basics(&xD, &opt);
+ /* ---------------------------------------------------------- */
+
+ /* Get... */
+ get_region_pars(&xD, ®); // ... region parameters
+ read_points(&map, ®, &pnts, &pclp, xD, opt.zcol->answer, field, &xD.report); // ... coordinates of points
+ pnts.invals = get_col_values(&map, &xD, &pnts, field, opt.intpl->answer, flg.detrend->answer); // ... values for interpolation
+
+ /* Estimate 2D/3D variogram */
+ switch (xD.phase) {
+ case 0:
+ /* Determine maximal horizontal (and vertical) distance + lags */
+ var_par.hz.type = 0;
+ if (xD.i3 == TRUE) {
+ var_par.vert.type = 1;
+ }
+
+ variogram_restricts(&xD, &pnts, pclp.pnts_hz, &var_par.hz);
+ if (xD.i3 == TRUE) {
+ variogram_restricts(&xD, &pnts, pclp.pnts_vert, &var_par.vert);
+ }
+
+ if (xD.i3 == TRUE && var_par.vert.nLag > var_par.hz.nLag) {
+ var_par.vert.nLag = var_par.hz.nLag;
+ var_par.vert.lag = var_par.vert.max_dist / var_par.vert.nLag;
+ write2file_varSets(&xD.report, &var_par.vert);
+ }
+
+ E_variogram(0, &xD, &pnts, &pclp, &var_par);
+ if (xD.i3 == TRUE) {
+ E_variogram(1, &xD, &pnts, &pclp, &var_par);
+ G_message(_("You may continue to computing theoretical variograms (middle phase)..."));
+ }
+ else {
+ G_message(_("You may continue to computing theoretical variograms (final phase)..."));
+ }
+ goto end;
+ case 1:
+ read_tmp_vals("variogram_hz_tmp.txt", &var_par.hz, &xD);
+ read_tmp_vals("variogram_vert_tmp.txt", &var_par.vert, &xD);
+
+ if (xD.report.name) {
+ xD.report.write2file = TRUE;
+ xD.report.fp = fopen(xD.report.name, "a");
+ }
+
+ T_variogram(0, opt, &var_par.hz, &xD.report); // Theoretical variogram
+ T_variogram(1, opt, &var_par.vert, &xD.report); // Theoretical variogram
+
+
+ /* difference between ranges
+ - if greater than 5% - bivariate
+ - if smaller than 5% - anisotropic
+ */
+ double sill95, diff_sill;
+ sill95 = var_par.hz.sill > var_par.vert.sill ? 0.05 * var_par.hz.sill : 0.05 * var_par.vert.sill;
+ diff_sill = fabsf(var_par.hz.sill - var_par.vert.sill);
+ if (xD.bivar == TRUE || (!flg.univariate->answer && diff_sill > sill95)) { // zonal anisotropy
+ var_par.fin.type = 2;
+ var_par.fin.td = var_par.hz.td;
+ var_par.fin.max_dist = var_par.hz.max_dist;
+ var_par.fin.max_dist_vert = var_par.vert.max_dist;
+ variogram_restricts(&xD, &pnts, pclp.pnts, &var_par.fin);
+ E_variogram(2, &xD, &pnts, &pclp, &var_par);
+ }
+
+ else if (xD.univar == TRUE || (!flg.bivariate->answer && diff_sill <= sill95)) { // geometric anisotropy
+ var_par.fin.type = 3;
+ var_par.fin.max_dist = var_par.hz.max_dist;
+ var_par.fin.td = var_par.hz.td;
+
+ xD.aniso_ratio = var_par.hz.h_range / var_par.vert.h_range;
+ geometric_anisotropy(&xD, &pnts, pclp.pnts);
+
+ variogram_restricts(&xD, &pnts, pclp.pnts, &var_par.fin);
+
+ E_variogram(3, &xD, &pnts, &pclp, &var_par);
+ }
+ G_message(_("You may continue to interpolate values (final phase)..."));
+ goto end;
+ case 2:
+ if (!opt.output->answer)
+ G_fatal_error(_("Please set up name of output layer..."));
+ out.name = opt.output->answer; // Output layer name
+
+ read_tmp_vals("variogram_final_tmp.txt", &var_par.fin, &xD);
+
+ if (xD.report.name) {
+ xD.report.write2file = TRUE;
+ xD.report.fp = fopen(xD.report.name, "a");
+ if (xD.report.fp == NULL)
+ G_fatal_error(_("Cannot open the file..."));
+ }
+
+ if (var_par.fin.type == 3)
+ geometric_anisotropy(&xD, &pnts, pclp.pnts);
+
+ T_variogram(var_par.fin.type, opt, &var_par.fin, &xD.report); // Theoretical variogram
+ break;
+ }
+
+ /* Ordinary kriging (including 2D/3D raster creation) */
+ ordinary_kriging(&xD, ®, &pnts, &pclp, &var_par, &out);
+ /* ---------------------------------------------------------- */
+
+ end:
+ Vect_close(&map); // Close vector map
+
+ exit(EXIT_SUCCESS);
+}
+
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/quantile.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/quantile.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/quantile.cpp 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,223 @@
+#include "local_proto.h"
+
+/* These functions are taken from the module r.quantile (Clemens, G.) */
+
+struct bin
+{
+ unsigned long origin;
+ double minimum, maximum;
+ int base, count;
+};
+
+static double minimum, maximum;
+static int num_quants;
+static double *quants;
+static int num_slots;
+static unsigned int *slots;
+static double slot_size;
+static unsigned long total;
+static int num_values;
+static unsigned short *slot_bins;
+static int num_bins;
+static struct bin *bins;
+static double *values;
+
+static inline int get_slot(double c)
+{
+ int i = (int) floor((c - minimum) / slot_size);
+
+ if (i < 0)
+ i = 0;
+ if (i > num_slots - 1)
+ i = num_slots - 1;
+ return i;
+}
+
+static inline double get_quantile(int n)
+{
+ return (double)total * quants[n]; // # of lower values
+}
+
+static void get_slot_counts(int n, double *data)
+{
+ int i;
+
+ total = 0;
+ for (i = 0; i < n; i++) {
+ int j;
+
+ // todo nan
+ j = get_slot(data[i]);
+
+ slots[j]++; // rise value of j-th slot
+ total++;
+ }
+ //G_percent(i, n, 2);
+}
+
+static void initialize_bins(void)
+{
+ int slot; // index of slot
+ double next; // percentile
+ int bin = 0; // adjacent bin
+ unsigned long accum = 0;
+ int quant = 0; // index of percentile
+
+ num_values = 0;
+ next = get_quantile(quant); // # of lower values
+
+ for (slot = 0; slot < num_slots; slot++) {
+ unsigned int count = slots[slot]; // assign # of elements in each slots to the count
+ unsigned long accum2 = accum + count; // total # of elements
+
+ if (accum2 > next) { // # of values is greater than percentile
+ struct bin *b = &bins[bin]; // make bin
+ slot_bins[slot] = ++bin;
+
+ b->origin = accum; // origin of bin = total # of elements
+ b->base = num_values;
+ b->count = 0;
+ b->minimum = minimum + slot_size * slot;
+ b->maximum = minimum + slot_size * (slot + 1);
+
+ while (accum2 > next)
+ next = get_quantile(++quant);
+
+ num_values += count;
+ }
+
+ accum = accum2;
+ }
+
+ num_bins = bin;
+}
+
+static void fill_bins(int n, double *data)
+{
+ int i;
+
+ for (i = 0; i < n; i++) {
+ int j, bin;
+ struct bin *b;
+
+ // todo nan
+
+ j = get_slot(data[i]);
+
+ if (!slot_bins[j])
+ continue;
+
+ bin = slot_bins[j] - 1;
+ b = &bins[bin];
+
+ values[b->base + b->count++] = data[i];
+ }
+
+ //G_percent(i, n, 2);
+}
+
+static int compare(const void *aa, const void *bb)
+{
+ double a = *(const double *)aa;
+ double b = *(const double *)bb;
+
+ if (a < b)
+ return -1;
+ if (a > b)
+ return 1;
+ return 0;
+}
+
+static void sort_bins(void)
+{
+ int bin;
+
+ for (bin = 0; bin < num_bins; bin++) {
+ struct bin *b = &bins[bin];
+
+ qsort(&values[b->base], b->count, sizeof(double), compare);
+ //G_percent(bin, num_bins, 2);
+ }
+ //G_percent(bin, num_bins, 2);
+}
+
+static void compute_quantiles(int recode, double *quantile, struct write *report)
+{
+ int bin = 0;
+ double prev_v = minimum;
+ int quant;
+
+ for (quant = 0; quant < num_quants; quant++) {
+ struct bin *b;
+ double next = get_quantile(quant);
+ double k, v;
+ int i0, i1;
+
+ for (; bin < num_bins; bin++) {
+ b = &bins[bin];
+ if (b->origin + b->count >= next)
+ break;
+ }
+
+ if (bin < num_bins) {
+ k = next - b->origin;
+ i0 = (int)floor(k);
+ i1 = (int)ceil(k);
+
+ if (i1 > b->count - 1)
+ i1 = b->count - 1;
+
+ v = (i0 == i1)
+ ? values[b->base + i0]
+ : values[b->base + i0] * (i1 - k) +
+ values[b->base + i1] * (k - i0);
+ }
+ else
+ v = maximum;
+
+ *quantile = v;
+ prev_v = v;
+ if (report != NULL && report->write2file == TRUE)
+ fprintf(report->fp, "%f:%f:%f\n", prev_v, v, quants[quant]);
+ }
+}
+
+double quantile(double q, int n, double *data, struct write *report)
+{
+ int recode=FALSE;
+ double quantile;
+
+ num_slots = 1000000; // # of slots
+
+ num_quants = 1;
+ quants = (double *) G_malloc(num_quants * sizeof(double));
+ quants[0] = q; // compute quantile for 0.05
+
+ // initialize list of slots
+ slots = (unsigned int *) G_calloc(num_slots, sizeof(unsigned int));
+ slot_bins = (unsigned short *) G_calloc(num_slots, sizeof(unsigned short));
+
+ // get min and max of the data
+ minimum = data[0];
+ maximum = data[0];
+ for (int i=0; i<n; i++) {
+ minimum = MIN(data[i], minimum);
+ maximum = MAX(data[i], maximum);
+ }
+
+ slot_size = (maximum - minimum) / num_slots;
+ get_slot_counts(n, data); // # of data in each slot
+
+ bins = (struct bin *) G_calloc(num_quants, sizeof(struct bin));
+ initialize_bins();
+ G_free(slots);
+
+ values = (double *) G_calloc(num_values, sizeof(double));
+ fill_bins(n, data);
+ G_free(slot_bins);
+
+ sort_bins();
+ compute_quantiles(recode, &quantile, report);
+
+ return quantile;
+}
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/stat.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/stat.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/stat.cpp 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,52 @@
+#include "local_proto.h"
+
+void test_normality(int n, double *data, struct write *report)
+{
+ int i;
+ double *elm;
+ double sum_data = 0.;
+ double res;
+ double sum_res2, sum_res3, sum_res4;
+ double tmp;
+ double mean, skewness, kurtosis;
+ double limit5, limit1;
+
+ // compute mean
+ elm = &data[0];
+ for (i=0; i<n; i++) {
+ sum_data += *elm;
+ elm++;
+ }
+ mean = sum_data / n;
+
+ // compute skewness and curtosis
+ elm = &data[0];
+ sum_res2 = sum_res3 = sum_res4 = 0.;
+ for (i=0; i<n; i++) {
+ res = *elm - mean;
+ sum_res2 += SQUARE(res);
+ sum_res3 += POW3(res);
+ sum_res4 += POW4(res);
+ elm++;
+ }
+
+ tmp = sum_res2 / (n-1);
+ skewness = (sum_res3 / (n-1)) / sqrt(POW3(tmp));
+ kurtosis = (sum_res4 / (n-1)) / SQUARE(tmp) - 3.;
+
+ // expected RMSE
+ limit5 = sqrt(6./n);
+ limit1 = sqrt(24./n);
+
+ fprintf(report->fp,"\n");
+ fprintf(report->fp,"Test of normality\n");
+ fprintf(report->fp,"-----------------\n");
+ fprintf(report->fp,"Test statistics | Value\n");
+ fprintf(report->fp,"skewness | %f\n", skewness);
+ fprintf(report->fp,"kurtosis | %f\n", kurtosis);
+ fprintf(report->fp,"\n");
+ fprintf(report->fp,"Level of confidence | Critical value\n");
+ fprintf(report->fp," 0.05 | %f\n", limit5);
+ fprintf(report->fp," 0.01 | %f\n", limit1);
+ fprintf(report->fp,"-----------------\n\n");
+}
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils.cpp 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,790 @@
+#include "local_proto.h"
+
+#define GNUPLOT "gnuplot -persist"
+
+// make coordinate triples xyz
+double *triple(double x, double y, double z)
+{
+ double *t;
+ t = (double *) G_malloc(3 * sizeof(double));
+
+ *t = x;
+ *(t+1) = y;
+ *(t+2) = z;
+
+ return t;
+}
+
+// compute coordinate differences
+void coord_diff(int i, int j, double *r, double *dr)
+{
+ int k=3*i,l=3*j;
+ double *rk, *rl, zi, zj, *drt;
+ rk = &r[k];
+ rl = &r[l];
+ drt = &dr[0];
+
+ if (*rk == 0. && *(rk+1) == 0. && *(rk+2) == 0.)
+ G_fatal_error(_("Coordinates of point no. %d are zeros."),i);
+
+ *drt = *rl - *rk; // dx
+ *(drt+1) = *(rl+1) - *(rk+1); // dy
+ *(drt+2) = *(rl+2) - *(rk+2); // dz
+}
+
+// compute horizontal radius from coordinate differences
+double radius_hz_diff(double *dr)
+{
+ double rds;
+ rds = SQUARE(dr[0]) + SQUARE(dr[1]);
+ return rds;
+}
+
+// compute size of the lag
+double lag_distance(int direction, struct points *pnts, pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts, struct parameters *var_par, struct write *report)
+{
+ int n = pnts->n; // # of input points
+ double *r = pnts->r; // xyz of input points
+ int fction = var_par->function;
+ double max_dist = (var_par->type == 2 && direction == 3) ? var_par->max_dist_vert : var_par->max_dist; // maximum horizontal distance
+
+ int i, j, k;
+ pcl::PointXYZ *searchPt;
+ searchPt = &pcl_pnts->points[0];
+
+ pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree (new pcl::KdTreeFLANN<pcl::PointXYZ>);
+ kd_tree->setInputCloud(pcl_pnts);
+
+ std::vector<int> ind; // indices of NN
+ std::vector<float> sqDist; // sqDistances of NN
+
+ double dr[3]; // coordinate differences
+ double quant_dist2 = var_par->radius; // quantile
+ double lagNN; // squared quantile - lag
+ double *dist2, *d; // radius
+ dist2 = (double *) G_malloc(n * sizeof(double));
+ d = &dist2[0];
+
+ int perc5, zeros;
+ int nN;
+ perc5 = (int) round(0.05 * n);
+ double dist0;
+ int j0 = 0;
+
+ for (i=0; i < n; i++) {
+ zeros = 0; // # of zero distances
+ dist0 = 0.; // initial value of 5% distance
+ while (dist0 != 0.) { // while 5% distance is zero
+ perc5 += zeros; // 5% plus points in zero distance
+ // find 5% nearest points
+ if (kd_tree->nearestKSearch(*searchPt, perc5, ind, sqDist) > 0) {
+ j = ind.size () - 1; // store index of 5%-th point
+ dist0 = sqDist.data ()[j]; // 5% distance
+
+ zeros = 0; // count zero distances
+ for (k=j0; k <= j; k++) {
+ if (sqDist.data ()[k] == 0.)
+ zeros++;
+ }
+
+ if (zeros > 0) {
+ dist0 = 0.;
+ j0 = j+1;
+ }
+ } // end nearestKSearch
+
+ else {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Error searching nearest neighbours of point %d..."), i);
+ } // end error
+ } // end while dist0 == 0
+
+ coord_diff(i, j, r, dr); // coordinate differences
+
+ // squared distance
+ switch (direction) {
+ case 12: // horizontal
+ *d = SQUARE(dr[0]) + SQUARE(dr[1]);
+ break;
+ case 3: // vertical
+ *d = SQUARE(dr[2]);
+ break;
+ case 0: // all
+ *d = SQUARE(dr[0]) + SQUARE(dr[1]) + SQUARE(dr[2]);
+ break;
+ }
+ if (*d != 0.)
+ quant_dist2 = MIN(quant_dist2, *d);
+ d++;
+ searchPt++;
+ }
+
+ lagNN = sqrt(quant_dist2);
+
+ return lagNN;
+}
+
+double linear(double x, double a, double b)
+{
+ double y;
+ y = a * x + b;
+
+ return y;
+}
+
+double exponential(double x, double nugget, double part_sill, double h_range)
+{
+ double y;
+ y = nugget + part_sill * (1. - exp(- 3. * x/h_range)); // practical
+
+ return y;
+}
+
+double spherical(double x, double a, double b, double c)
+{
+ double y;
+ double ratio;
+ if (x < c) {
+ ratio = x / c;
+ y = a + b * (1.5 * ratio - 0.5 * POW3(ratio));
+ }
+ else
+ y = a + b;
+
+ return y;
+}
+
+double gaussian(double x, double a, double b, double c)
+{
+ double y;
+ double c2 = SQUARE(c);
+ double ratio;
+ ratio = x / c2;
+ y = a + b * (1 - exp(-3. * ratio));
+
+ return y;
+}
+
+// compute # of lags
+int lag_number(double lag, double *varpar_max)
+{
+ int n; // number of lags
+ n = (int) round(*varpar_max / lag); // compute number of lags
+ *varpar_max = n * lag; // set up modified maximum distance (to not have empty lags)
+
+ return n;
+}
+
+// maximal horizontal and vertical distance to restrict variogram computing
+void variogram_restricts(struct int_par *xD, struct points *pnts, pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts, struct parameters *var_par)
+{
+ int n = pnts->n; // # of points
+ double *r = pnts->r; // xyz of points
+ struct write *report = &xD->report;
+
+ int i;
+ double *min, *max; // extend
+ double dr[3]; // coordinate differences
+ double h_maxG; // modify maximum distance (to not have empty lags)
+
+ char type[12];
+ variogram_type(var_par->type, type);
+
+ // Find extent
+ G_message(_("Computing %s variogram properties..."), type);
+ min = &pnts->r_min[0]; // set up pointer to minimum xyz
+ max = &pnts->r_max[0]; // set up pointer to maximum xyz
+
+ dr[0] = *max - *min; // x range
+ dr[1] = *(max+1) - *(min+1); // y range
+ dr[2] = *(max+2) - *(min+2); // z range
+
+ switch (var_par->type) {
+ case 1:
+ var_par->max_dist = dr[2]; // zmax - zmin
+ var_par->radius = SQUARE(var_par->max_dist); // anisotropic distance (todo: try also dz)
+ break;
+ default:
+ var_par->radius = radius_hz_diff(dr) / 9.; // horizontal radius (1/9)
+ var_par->max_dist = sqrt(var_par->radius); // 1/3 max horizontal dist (Surfer, Golden SW)
+ break;
+ }
+
+ if (report->name)
+ write2file_varSetsIntro(var_par->type, report);
+
+ int dimension;
+ switch (var_par->type) {
+ case 0:
+ dimension = 12;
+ break;
+ case 1:
+ dimension = 3;
+ break;
+ case 3:
+ dimension = 0;
+ break;
+ }
+
+ if (var_par->type == 2) {
+ var_par->lag = lag_distance(12, pnts, pcl_pnts, var_par, report); // the shortest distance between NN in horizontal direction
+ var_par->nLag = lag_number(var_par->lag, &var_par->max_dist);
+ var_par->max_dist = var_par->nLag * var_par->lag;
+
+ var_par->lag_vert = lag_distance(3, pnts, pcl_pnts, var_par, report); // the shortest distance between NN in horizontal direction
+ var_par->nLag_vert = lag_number(var_par->lag_vert, &var_par->max_dist_vert);
+ var_par->max_dist_vert = var_par->nLag_vert * var_par->lag_vert;
+ }
+
+ else {
+ var_par->lag = lag_distance(dimension, pnts, pcl_pnts, var_par, report); // the shortest distance between NN in horizontal direction
+ var_par->nLag = lag_number(var_par->lag, &var_par->max_dist);
+ }
+
+ write2file_varSets(report, var_par);
+}
+
+void geometric_anisotropy(struct int_par *xD, struct points *pnts, pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts)
+{
+ int i;
+ double ratio = xD->aniso_ratio;
+ double *z;
+ z = &pnts->r[2];
+
+ for (i=0; i<pnts->n; i++) {
+ *z = pcl_pnts->points[i].z = ratio * *z;
+ z += 3;
+ }
+}
+
+// Least Squares Method
+mat_struct *LSM(mat_struct *A, mat_struct *x)
+{
+ int i, nr, nc;
+ mat_struct *AT, *ATA, *ATA_Inv, *ATx, *T;
+
+ /* A[nr x nc] */
+ nr = A->rows;
+ nc = A->cols;
+
+ /* LMS */
+ AT = G_matrix_transpose(A); /* Transposed design matrix */
+ ATA = G_matrix_product(AT, A); /* AT*A */
+ ATA_Inv = G_matrix_inverse(ATA); /* (AT*A)^-1 */
+ ATx = G_matrix_product(AT, x); /* AT*x */
+ T = G_matrix_product(ATA_Inv, ATx); /* theta = (AT*A)^(-1)*AT*x */
+
+ return T;
+}
+
+mat_struct *nonlin_LMS(int n, double *dist, double *gamma)
+{
+ int i=0, j;
+ double ctrl=1.;
+ double a, b, c, mean_dy = 1., sum_dy, *h, *g, *add;
+ mat_struct *y, *param, *param0, *J, *yf, *JT, *JTJ, *JT_Inv, *JTdy, *dy, *delta;
+ y = G_matrix_init(n, 1, n);
+ param = G_matrix_init(3, 1, 3);
+ param0 = G_matrix_init(3, 1, 3);
+ J = G_matrix_init(n, 3, n);
+ yf = G_matrix_init(n, 1, n);
+
+ // matrix of input elements
+ g = &gamma[0];
+ for (j=0; j<n; j++) {
+ G_matrix_set_element(y, j, 0, *g);
+ g++;
+ }
+
+ while (fabsf(mean_dy) >= 0.0001) {
+ if (i==0) { // set initial parameters
+ G_matrix_set_element(param, 0, 0, 0.);
+ G_matrix_set_element(param, 1, 0, 45.);
+ G_matrix_set_element(param, 2, 0, 45.);
+ }
+ a = param->vals[0];
+ b = param->vals[1];
+ c = param->vals[2];
+ G_debug(0,"%d a=%f b=%f c=%f", i, a, b, c);
+
+ h = &dist[0];
+ for (j=0; j<n; j++) {
+ // Jacobi matrix
+ // todo: rozdelit na variogramy
+ G_matrix_set_element(J, j, 0, 1.);
+ G_matrix_set_element(J, j, 1, 1. - exp(-3. * *h / c));
+ G_matrix_set_element(J, j, 2, -3. * b * *h/SQUARE(c) * exp(-3. * *h*c));
+
+ // f(a,b,c)
+ G_matrix_set_element(yf, j, 0, b * (1. - exp(-3. * *h/c)));
+ G_debug(0, "%f %f %f", G_matrix_get_element(J, j, 0), G_matrix_get_element(J, j, 1), G_matrix_get_element(J, j, 2));
+ h++;
+ }
+
+ // differences
+ dy = G_matrix_subtract(y, yf);
+ sum_dy = 0.;
+ add = &dy->vals[0];
+ for (j=0; j<n; j++) {
+ sum_dy += SQUARE(*add);
+ add++;
+ }
+ mean_dy = sum_dy / n;
+ G_debug(0,"%d mean=%f", i, mean_dy);
+
+ JT = G_matrix_transpose(J);
+ JTJ = G_matrix_product(JT, J);
+ JT_Inv = G_matrix_inverse(JTJ);
+ JTdy = G_matrix_product(JT, dy);
+ delta = G_matrix_product(JT_Inv, JTdy);
+ G_debug(0,"a=%f b=%f c=%f", delta->vals[0], delta->vals[1], delta->vals[2]);
+
+ param0 = G_matrix_copy(param);
+ param = G_matrix_subtract(param0, delta);
+ i++;
+ if (i>3)
+ G_fatal_error(_("stoj"));
+ }
+ G_debug(0,"nugget=%f sill=%f range=%f", param->vals[0], param->vals[1], param->vals[2]);
+ return param;
+}
+
+// set up type of function according to GUI (user)
+int set_function(char *variogram, struct parameters *var_par, struct write *report)
+{
+ int function;
+ if (strcmp(variogram, "linear") == 0) {
+ function = 0;
+ // var_par->bivar = FALSE;
+ }
+ else if (strcmp(variogram, "parabolic") == 0) {
+ function = 1;
+ //var_par->bivar = FALSE;
+ }
+ else if (strcmp(variogram, "exponential") == 0) {
+ function = 2;
+ // var_par->bivar = FALSE;
+ }
+ else if (strcmp(variogram, "spherical") == 0) {
+ function = 3;
+ //var_par->bivar = FALSE;
+ }
+ else if (strcmp(variogram, "gaussian") == 0) {
+ function = 4;
+ //var_par->bivar = FALSE;
+ }
+ else if (strcmp(variogram, "bivariate") == 0) {
+ function = 5;
+ //var_par->bivar = TRUE;
+ }
+ else {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Set up correct name of variogram function..."));
+ }
+
+ return function;
+}
+
+// set up terminal and format of output plot
+void set_gnuplot(char *fileformat, struct parameters *var_par)
+{
+ if (strcmp(fileformat,"cdr") == 0) {
+ strcpy(var_par->term, "corel");
+ strcpy(var_par->ext, "cdr");
+ }
+ if (strcmp(fileformat,"dxf") == 0) {
+ strcpy(var_par->term, "dxf");
+ strcpy(var_par->ext, "dxf");
+ }
+ if (strcmp(fileformat,"eps") == 0) {
+ strcpy(var_par->term, "postscript");
+ strcpy(var_par->ext, "eps");
+ }
+ if (strcmp(fileformat,"pdf") == 0) {
+ strcpy(var_par->term, "pdfcairo");
+ strcpy(var_par->ext, "pdf");
+ }
+ if (strcmp(fileformat,"png") == 0) {
+ strcpy(var_par->term, "png");
+ strcpy(var_par->ext, "png");
+ }
+ if (strcmp(fileformat,"svg") == 0) {
+ strcpy(var_par->term, "svg");
+ strcpy(var_par->ext, "svg");
+ }
+}
+
+// plot experimental variogram
+void plot_experimental_variogram(struct int_par *xD, double *h_exp, mat_struct *gamma_exp, struct parameters *var_par)
+{
+ int bivar = var_par->type == 2 ? TRUE : FALSE;
+
+ int i, j; // indices
+ int nr, nc; // # of rows, cols
+ double *h; // pointer to horizontal or anisotropic bins
+ double *vert;
+ double *gamma_var; // pointer to gamma matrix
+ double *c; // pointer to number of elements
+ FILE *gp; // pointer to file
+
+ nr = gamma_exp->rows; // # of rows of gamma matrix
+ nc = gamma_exp->cols; // # of cols of gamma matrix
+
+ gamma_var = &gamma_exp->vals[0]; // values of gamma matrix
+ c = &var_par->c->vals[0];
+
+ gp = fopen("dataE.dat","w"); // open file to write experimental variogram
+ if (gp == NULL)
+ G_fatal_error(_("Error writing file"));
+
+ /* 3D experimental variogram */
+ if (bivar == TRUE) {
+ vert = &var_par->vert[0];
+ for (i=0; i < nc+1; i++) { // for each row (nZ)
+ h = &var_par->h[0];
+ for (j=0; j < nr+1; j++) { // for each col (nL)
+ if (i == 0 && j == 0)
+ fprintf(gp, "%d ", nr); // write to file count of columns
+ else if (i == 0 && j != 0) { // 1st row
+ fprintf(gp, "%f", *h); // write h to 1st row of the file
+ h++;
+ if (j < nr)
+ fprintf(gp, " "); // spaces between elements
+ } // end 1st row setting
+
+ else {
+ if (j == 0 && i != 0) // 1st column
+ fprintf(gp, "%f ", *vert);// write vert to 1st column of the file
+ else {
+ if (*c == 0.)
+ fprintf(gp, "NaN"); // write other elements
+ else
+ fprintf(gp, "%f", *gamma_var); // write other elements
+
+ if (j != nr)
+ fprintf(gp, " "); // spaces between elements
+ gamma_var++;
+ c++;
+ }
+ } // end columns settings
+ } // end j
+ fprintf(gp, "\n");
+ if (i!=0) // vert must not increase in the 1st loop
+ vert++;
+ } // end i
+ } // end 3D
+
+ /* 2D experimental variogram */
+ else {
+ h = &h_exp[0];
+ for (i=0; i < nr; ++i) {
+ if (*c == 0.)
+ fprintf(gp, "%f NaN\n", *h);
+ else
+ fprintf(gp, "%f %f\n", *h, *gamma_var);
+ c++;
+ h++;
+ gamma_var++;
+ }
+ }
+ fclose(gp);
+
+ gp = popen(GNUPLOT, "w"); /* open file to create plots */
+ if (gp == NULL)
+ G_message("Unable to plot variogram");
+
+ fprintf(gp, "set terminal wxt %d size 800,450\n", var_par->type);
+
+ if (bivar == TRUE) { // bivariate variogram
+ fprintf(gp, "set title \"Experimental variogram (3D) of <%s>\"\n", var_par->name);
+ fprintf(gp, "set xlabel \"h interval No\"\n");
+ fprintf(gp, "set ylabel \"dz interval No\"\n");
+ fprintf(gp, "set zlabel \"gamma(h,dz) [units^2]\" rotate by 90 left\n");
+ fprintf(gp, "set key below\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp, "set dgrid3d %d,%d\n", nc, nr);
+ fprintf(gp, "splot 'dataE.dat' every ::1:1 matrix title \"experimental variogram\"\n");
+ }
+
+ else { // univariate variogram
+ char dim[6], dist[2];
+ if (xD->i3 == TRUE) { // 3D
+ if (var_par->type == 0)
+ strcpy(dim,"hz");
+ if (var_par->type == 1)
+ strcpy(dim, "vert");
+ if (var_par->type == 3)
+ strcpy(dim, "aniso");
+ fprintf(gp, "set title \"Experimental variogram (3D %s) of <%s>\"\n", dim, var_par->name);
+ }
+ else // 2D
+ fprintf(gp, "set title \"Experimental variogram (2D) of <%s>\"\n", var_par->name);
+
+ fprintf(gp, "set xlabel \"h [m]\"\n");
+ fprintf(gp, "set ylabel \"gamma(h) [units^2]\"\n");
+ fprintf(gp, "set key bottom right\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp, "plot 'dataE.dat' using 1:2 title \"experimental variogram\" pointtype 5\n");
+ }
+ fclose(gp);
+
+ //remove("dataE.dat");
+}
+
+// plot experimental and theoretical variogram
+void plot_var(int bivar, double *h_exp, mat_struct *gamma_exp, struct parameters *var_par)
+{
+ int function, func_h, func_v;
+ double nugget, nugget_h, nugget_v;
+ double sill, sill_h, sill_v;
+ double h_range, h_range_h, h_range_v;
+
+ switch (bivar) {
+ case FALSE:
+ function = var_par->function;
+ nugget = var_par->nugget;
+ sill = var_par->sill - nugget;
+ h_range = var_par->h_range;
+ break;
+ case TRUE:
+ if (var_par->function != 5) {
+ func_h = var_par->horizontal.function;
+ func_v = var_par->vertical.function;
+
+ nugget_h = var_par->horizontal.nugget;
+ sill_h = var_par->horizontal.sill - nugget_h;
+ h_range_h = var_par->horizontal.h_range;
+
+ nugget_v = var_par->vertical.nugget;
+ sill_v = var_par->vertical.sill - nugget_v;
+ h_range_v = var_par->vertical.h_range;
+ }
+ break;
+ }
+
+ int i, j; // indices
+ int nr, nc; // # of rows, cols
+ double *h; // pointer to horizontal or anisotropic bins
+ double *vert; // pointer to vertical bins
+ //double *T_cf; // coefficients of theoretical variogram - vals
+ double h_ratio;
+ double g_teor;
+ double *c;
+ double *gamma_var; // pointer to gamma matrix
+ double *T; // coefficients of theoretical variogram - matrix
+ FILE *gp; // pointer to file
+
+ nr = gamma_exp->rows; // # of rows of gamma matrix
+ nc = gamma_exp->cols; // # of cols of gamma matrix
+
+ c = &var_par->c->vals[0];
+ gamma_var = &gamma_exp->vals[0]; // values of gamma matrix
+
+ if (var_par->type != 2) { // univariate variogram
+ gp = fopen("dataE.dat","w"); // open file to write experimental variogram
+ if (gp == NULL)
+ G_fatal_error(_("Error writing file"));
+
+ h = &h_exp[0];
+ for (i=0; i < nr; i++) {
+ if (*c == 0.)
+ fprintf(gp, "%f NaN\n", *h); // write other elements
+ else
+ fprintf(gp, "%f %f\n", *h, *gamma_var);
+ h++;
+ c++;
+ gamma_var++;
+ }
+ }
+ else {
+ gp = fopen("dataE.dat","r"); // open file to write experimental variogram
+ if (gp == NULL)
+ G_fatal_error(_("You have probably deleted dataE.dat - process middle phase again, please."));
+ }
+
+ if (fclose(gp) != 0)
+ G_fatal_error(_("Error closing file..."));
+
+ /* Theoretical variogram */
+ gp = fopen("dataT.dat","w"); // open file to write theoretical variogram
+ if (gp == NULL)
+ G_fatal_error(_("Error opening file..."));
+
+ double dd;
+ double hh;
+ double hv[2];
+ h = &h_exp[0];
+
+
+
+ switch (bivar) {
+ case 0: // Univariate variogram
+ for (i=0; i <= nr; i++) {
+ hh = i==0 ? 0. : *h;
+
+ switch (function) {
+ case 0: // linear
+ dd = *T * hh + *(T+1);
+ break;
+ case 1: // parabolic
+ dd = *T * SQUARE(hh) + *(T+1);
+ break;
+ case 2: // exponential
+ dd = nugget + sill * (1 - exp(- 3. * hh/h_range)); // practical
+ break;
+ case 3: // spherical
+ if (hh < h_range) {
+ h_ratio = hh / h_range;
+ dd = nugget + sill * (1.5 * h_ratio - 0.5 * POW3(h_ratio));
+ }
+ else
+ dd = sill + nugget;
+ break;
+ case 4: // Gaussian
+ h_ratio = SQUARE(hh) / SQUARE(h_range);
+ dd = nugget + sill * (1 - exp(-3. * h_ratio));
+ break;
+ }
+ fprintf(gp, "%f %f\n", hh, dd);
+
+ if (i > 0)
+ h++;
+ } // end i for cycle
+ break;
+
+ case 1: // bivariate (3D)
+ // Coefficients of theoretical variogram
+ T = &var_par->T->vals[0]; // values
+ vert = &var_par->vert[0];
+ for (i=0; i < nc+1; i++) { // for each row
+ h = &var_par->h[0];
+ for (j=0; j < nr+1; j++) { // for each col
+ if (i == 0 && j == 0) // the 0-th element...
+ fprintf(gp, "%d ", nr); // ... is number of columns
+ else if (i == 0 && j != 0) { // for the other elements in 1st row...
+ fprintf(gp, "%f", *h); // ... write h
+ if (j < nr)
+ fprintf(gp, " "); // ... write a space
+ h++; // go to next h
+ }
+ else { // other rows
+ if (j == 0 && i != 0) // write vert to 1st column
+ fprintf(gp, "%f ", *vert);
+ else { // fill the other columns with g_teor
+ hv[0] = *h;
+ hv[1] = *vert;
+ g_teor = variogram_fction(var_par, hv);
+ fprintf(gp, "%f", g_teor);
+ if (j != nr)
+ fprintf(gp, " ");
+ h++;
+ } // end else: fill the other columns
+ } // end fill the other rows
+ } //end j
+ if (i != 0 && j != 0)
+ vert++;
+ fprintf(gp, "\n");
+ } //end i
+ break;
+ }
+
+ if (fclose(gp) == EOF)
+ G_fatal_error(_("Error closing file..."));
+
+ gp = popen(GNUPLOT, "w"); /* open file to create plots */
+ if (gp == NULL)
+ G_message(_("Unable to plot variogram"));
+
+ if (strcmp(var_par->term,"") != 0) {
+ //fprintf(gp, "set terminal %s size 750,450\n", var_par->term);
+ fprintf(gp, "set terminal wxt size 750,450\n", var_par->type);
+ switch (var_par->type) {
+ case 0:
+ if (bivar == TRUE)
+ fprintf(gp, "set output \"variogram_bivar.%s\"\n", var_par->ext);
+ else
+ fprintf(gp, "set output \"variogram_horizontal.%s\"\n", var_par->ext);
+ break;
+ case 1:
+ fprintf(gp, "set output \"variogram_vertical.%s\"\n", var_par->ext);
+ break;
+ case 2:
+ if (bivar == TRUE)
+ fprintf(gp, "set output \"variogram_bivariate.%s\"\n", var_par->ext);
+ else
+ fprintf(gp, "set output \"variogram_anisotropic.%s\"\n", var_par->ext);
+ break;
+ }
+ }
+ else {
+ G_warning("\nVariogram output> If you wish some special file format, please set it at the beginning.\n");
+ fprintf(gp, "set terminal wxt %d size 800,450\n", var_par->type);
+ }
+
+ if (bivar == TRUE) { // bivariate variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D) of <%s>\"\n", var_par->name);
+
+ fprintf(gp, "set xlabel \"h interval No\"\n");
+ fprintf(gp, "set ylabel \"dz interval No\"\n");
+ fprintf(gp, "set zlabel \"gamma(h,dz) [units^2]\" rotate by 90 left\n");
+ fprintf(gp, "set key below\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp, "set dgrid3d %d,%d\n", nc, nr);
+ fprintf(gp, "splot 'dataE.dat' every ::1:1 matrix title \"experimental variogram\", 'dataT.dat' every ::1:1 matrix with lines title \"theoretical variogram\" palette\n");
+ }
+
+ else { // univariate variogram*/
+ //if (xD->i3 == TRUE) { // 3D
+ if (var_par->type == 0) { // horizontal variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D hz) of <%s>\"\n", var_par->name);
+ //fprintf(gp, "set label \"linear: gamma(h) = %f*h + %f\" at screen 0.30,0.90\n", var_par->T->vals[0], var_par->T->vals[1]);
+ }
+ else if (var_par->type == 1) { // vertical variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D vert) of <%s>\"\n", var_par->name);
+ //fprintf(gp, "set label \"linear: gamma(h) = %f*h + %f\" at screen 0.30,0.90\n", var_par->T->vals[0], var_par->T->vals[1]);
+ }
+ else if (var_par->type == 3) // anisotropic variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D) of <%s>\"\n", var_par->name);
+ //}
+
+/*else // 2D
+ fprintf(gp, "set title \"Experimental and theoretical variogram (2D) of <%s>\"\n", var_par->name);*/
+
+ //if (var_par->type == 2/* || xD->i3 == FALSE*/) {
+ switch (var_par->function) {
+ case 0: // linear
+ fprintf(gp, "set label \"linear: gamma(h) = %f*h + %f\" at screen 0.30,0.90\n", var_par->T->vals[0], var_par->T->vals[1]);
+ break;
+ case 1: // parabolic
+ fprintf(gp, "set label \"parabolic: gamma(h) = %f*h^2 + %f\" at screen 0.25,0.90\n", var_par->T->vals[0], var_par->T->vals[1]);
+ break;
+ case 2: // exponential
+ fprintf(gp, "set rmargin 5\n");
+ fprintf(gp, "set label \"exponential: gamma(h) = %f + %f * (1 - exp(-3*h / %f))\" at screen 0.10,0.90\n", var_par->nugget, var_par->sill - var_par->nugget, var_par->h_range);
+ break;
+ case 3: // spherical
+ fprintf(gp, "set rmargin 5\n");
+ fprintf(gp, "set label \"spherical: gamma(h) = %f+%f*(1.5*h/%f-0.5*(h/%f)^3)\" at screen 0.05,0.90\n", var_par->nugget, var_par->sill - var_par->nugget, var_par->h_range, var_par->h_range);
+ break;
+ case 4: // gaussian
+ fprintf(gp, "set label \"gaussian: gamma(h) = %f + %f * (1 - exp(-3*(h / %f)^2))\" at screen 0.10,0.90\n", var_par->nugget, var_par->sill - var_par->nugget, var_par->h_range);
+ break;
+ }
+ //}
+ fprintf(gp, "set xlabel \"h [m]\"\n");
+ fprintf(gp, "set ylabel \"gamma(h) [units^2]\"\n");
+ fprintf(gp, "set key bottom right\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp, "plot 'dataE.dat' using 1:2 title \"experimental variogram\" pointtype 5, 'dataT.dat' using 1:2 title \"theoretical variogram\" with linespoints\n");
+ }
+ fclose(gp);
+
+ //remove("dataE.dat");
+ //remove("dataT.dat");
+}
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_kriging.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_kriging.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_kriging.cpp 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,473 @@
+#include "local_proto.h"
+
+ __inline double square(double x)
+ {
+ return x*x;
+ }
+
+// formulation of variogram functions
+double variogram_fction(struct parameters *var_par, double *dr)
+{
+ // Local variables
+ int variogram = var_par->function; // theoretical variogram
+ int type = var_par->type;
+ double nugget;
+ double part_sill;
+ double h_range;
+ double *T;
+
+ if (type == 2) {
+ T = &var_par->T->vals[0]; // coefficients of theoretical variogram
+ }
+
+ double radius; // square distance of the point couple
+ double h;
+ double vert;
+ double h_ratio;
+ double teor_var, result = 0.; // estimated value of the theoretical variogram
+
+ int n_cycles = (type == 2 && var_par->function != 5) ? 2 : 1; // # of cycles (bi- or univariate)
+
+ for (int i=0; i < n_cycles; i++) {
+ if (n_cycles == 2) { // bivariate variogram
+ variogram = i==0 ? var_par->horizontal.function : var_par->vertical.function;
+ h = i==0 ? dr[0] : dr[1];
+ radius = SQUARE(h);
+ nugget = i==0 ? var_par->horizontal.nugget : var_par->vertical.nugget;
+ part_sill = i==0 ? var_par->horizontal.sill - nugget : var_par->vertical.sill - nugget;
+ h_range = i==0 ? var_par->horizontal.h_range : var_par->vertical.h_range;
+ }
+ else { // univariate variogram or linear bivariate
+ variogram = var_par->function;
+ if (variogram == 5) {
+ h = dr[0];
+ vert = dr[1];
+ }
+ else {
+ radius = SQUARE(dr[0]) + SQUARE(dr[1]) + SQUARE(dr[2]); // weighted dz
+ h = sqrt(radius);
+
+ nugget = var_par->nugget;
+ part_sill = var_par->sill - nugget;
+ h_range = var_par->h_range;
+ }
+ }
+
+ switch (variogram) {
+ case 0: // linear function
+ teor_var = linear(h, *T, *(T+1));
+ break;
+ case 1: // parabolic function TODO: delete
+ teor_var = *T * radius + *(T+1);
+ break;
+ case 2: // exponential function
+ teor_var = exponential(h, nugget, part_sill, h_range);
+ break;
+ case 3: // spherical function
+ teor_var = spherical(h, nugget, part_sill, h_range);
+ break;
+ case 4: // Gaussian function
+ teor_var = gaussian(radius, nugget, part_sill, h_range);
+ break;
+ case 5:
+ teor_var = *T * h + *(T+1) * vert + *(T+2);
+ break;
+ } // end switch (variogram)
+
+ result += teor_var;
+ }
+
+ return result;
+}
+
+// compute coordinates of reference point - cell centre
+void cell_centre(unsigned int col, unsigned int row, unsigned int dep, struct int_par *xD, struct reg_par *reg, double *r0, struct parameters *var_par)
+{
+ // Local variables
+ int i3 = xD->i3;
+
+ r0[0] = reg->west + (col + 0.5) * reg->ew_res; // x0
+ r0[1] = reg->north - (row + 0.5) * reg->ns_res; // y0
+ if (i3 == TRUE)
+ switch (var_par->function) {
+ case 5:
+ r0[2] = reg->bot + (dep + 0.5) * reg->bt_res; // z0
+ break;
+ default:
+ r0[2] = xD->aniso_ratio * (reg->bot + (dep + 0.5) * reg->bt_res); // z0
+ break;
+ }
+ else
+ r0[2] = 0.;
+}
+
+// set up elements of G matrix
+mat_struct *set_up_G(struct points *pnts, struct parameters *var_par, struct write *report)
+{
+ // Local variables
+ int n = pnts->n; // number of input points
+ double *r = pnts->r; // xyz coordinates of input points
+ int type = var_par->type;
+
+ int i, j; // indices of matrix rows/cols
+ int n1 = n+1; // number of matrix rows/cols
+ double teor_var; // estimated value based on theoretical variogram and distance of the input points
+ double *dr; // coordinate differences dx, dy, dz between couples of points
+ mat_struct *GM; // GM matrix
+
+ dr = (double *) G_malloc(3 * sizeof(double));
+ GM = G_matrix_init(n1, n1, n1); // G[n1][n1] matrix
+
+ doublereal *md, *mu, *ml, *dbu, *dbl, *m1r, *m1c;
+ dbu = &GM->vals[0]; // upper matrix elements
+ dbl = &GM->vals[0]; // lower matrix elements
+ md = &GM->vals[0]; // diagonal matrix elements
+ m1r = &GM->vals[n1-1]; // GM - last matrix row
+ m1c = &GM->vals[n1*n]; // GM - last matrix col
+
+ // G[n;n]
+ for (i = 0; i < n1; i++) { // for elements in each row
+ dbu += n1; // new row of the U
+ dbl++; // new col of the L
+ mu = dbu; // 1st element in the new U row
+ ml = dbl; // 1st element in the new L col
+ for (j = i; j < n1; j++) { // for elements in each col of upper matrix
+ if (i != j && i != n && j != n) { // non-diagonal elements except last row/col
+ coord_diff(i, j, r, dr); // compute coordinate differences
+ if (type == 2) {
+ dr[0] = sqrt(SQUARE(dr[0]) + SQUARE(dr[1]));
+ dr[1] = dr[2];
+ }
+ teor_var = variogram_fction(var_par, dr);
+
+ if (isnan(teor_var)) {
+ if (report->name) {
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Theoretical variogram is NAN..."));
+ }
+
+ *mu = *ml = (doublereal) teor_var; // set the value to the matrix
+ mu += n1; // go to next element in the U row
+ ml++; // go to next element in the L col
+ } // end non-diagonal elements condition
+ } // end j loop
+
+ // go to the diagonal element in the next row
+ dbu++; // U
+ dbl += n1; // L
+ *md = 0.0; // set diagonal
+ md += n1+1;
+ // Set ones
+ if (i<n1-1) { // all elements except the last one...
+ *m1r = *m1c = 1.0; // ... shall be 1
+ m1r += n1; // go to next col in last row
+ m1c++; // go to next row in last col
+ } // end "last 1" condition
+ } // end i loop
+
+ free(dr);
+ return GM;
+}
+
+// make G submatrix for rellevant points
+mat_struct *submatrix(std::vector<int> index, mat_struct *GM_all, struct write *report)
+{
+ // Local variables
+ int n = index.size ();
+ int *lind;
+ lind = &index.data ()[0];
+ mat_struct *GM = GM_all;
+
+ int i, j, k, N1 = GM->rows, n1 = n+1, *dinR0, *dinR, *dini, *dini0, *dinj0, *dinj;
+ doublereal *dbo, *dbx, *dbu, *dbl, *md, *mu, *ml, *m1r, *m1c;
+
+ mat_struct *sub;
+ sub = G_matrix_init(n1,n1,n1);
+ if (sub == NULL) {
+ if (report->name) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to initialize G-submatrix..."));
+ }
+
+ // Initialize indexes: GM[0;1]
+ // - cols
+ dini = &lind[0]; // i - set processing column
+ dini0 = &lind[0]; // i-1 - previous column
+ // initialize new col
+ dinR0 = &lind[1]; // return to first processing row - lower GM
+ dinR = &lind[1]; // return to first processing row - lower GM
+
+ dbo = &GM->vals[*dini * N1 + *dinR]; // 1st value in GM_all
+ md = &sub->vals[0]; // GM - diagonal
+ dbu = &sub->vals[n1]; // GM - upper
+ dbl = &sub->vals[1]; // GM - lower
+ m1r = &sub->vals[n1-1]; // GM - last row
+ m1c = &sub->vals[n1*n]; // GM - last col
+
+ for (i=0; i<=n; i++) { // set up cols of submatrix
+ // original matrix
+ dbx = dbo; // new col in GM_all
+ dinj = dinR; // orig: inicialize element (row) in (col) - upper matrix
+ dinj0 = dinR; // orig: inicialize element (row - 1) in (col)
+ // submatrix
+ mu = dbu; // sub: start new column
+ ml = dbl; // sub: start new row
+ for (j=i+1; j<n; j++) { // for rows
+ //submatrix
+ *mu = *ml = *dbx; // general elements
+ //G_debug(0,"GM=%f", *dbx);
+ mu += n1; // go to element in next column
+ ml++; // go to element in next row
+ // Original matrix
+ dinj0 = dinj; // save current ind
+ dinj++; // set next ind
+ if (j < n-1)
+ dbx += *dinj - *dinj0;
+ } // end j
+ // Original matrix
+ dini0 = dini; // save current ind
+ dini++; // set next ind
+ if (i < n-1) {
+ dbu += n1+1; // new col in GM
+ dbl += n1+1; // new row in GM
+ }
+ dinR0 = dinR;
+ dinR++;
+ dbo += (*dini - *dini0) * N1 + (*dinR - *dinR0);
+ // Set ones
+ *m1r = *m1c = 1.0;
+ m1r += n1;
+ m1c++;
+ // Set diagonals
+ *md = 0.0;
+ md += n1+1;
+ } // end i
+
+ return sub;
+}
+
+// set up g0 vector
+mat_struct *set_up_g0(struct int_par *xD, std::vector<int> ind, struct points *pnts, double *r0, struct parameters *var_par)
+{
+ // Local variables
+ int i3 = xD->i3; // 2D or 3D interpolation
+ int bivar = xD->bivar;
+ int type = var_par->type;
+ int n;
+ double *r; // xyz coordinates of input points
+
+ int *lind, *lind0;
+
+ if (ind.empty()) {
+ n = pnts->n;
+ r = &pnts->r[0];
+ }
+ else {
+ n = ind.size (); // number of input points
+ lind = &ind[0];
+ lind0 = &ind[0];
+ r = &pnts->r[*lind * 3];
+ }
+ int n1 = n + 1;
+
+ int i; // index of elements and length of g0 vector
+ double teor_var; // estimated value based on theoretical variogram and distance of the input points
+ double *dr; // coordinate differences dx, dy, dz between couples of points
+ mat_struct *g0; // vector of estimated differences between known and unknown values
+
+ dr = (double *) G_malloc(3 * sizeof(double)); // Coordinate differences
+ g0 = G_matrix_init(n1,1,n1);
+
+ double *ro;
+ doublereal *g;
+ g = &g0->vals[0];
+
+ for (i=0; i<n; i++) { // count of input points
+ // Coord diffs (input points and cell/voxel center)
+ dr[0] = *r - r0[0]; // dx
+ dr[1] = *(r+1) - r0[1]; // dy
+ switch (i3) {
+ case FALSE:
+ // Cell value diffs estimation using linear variogram
+ dr[2] = 0.0; // !!! optimalize - not to be necessary to use this line
+ break;
+ case TRUE:
+ dr[2] = *(r+2) - r0[2]; // dz
+ break;
+ }
+ if (type == 2) {
+ dr[0] = sqrt(SQUARE(dr[0]) + SQUARE(dr[1]));
+ dr[1] = dr[2];
+ }
+
+ *g = variogram_fction(var_par, dr);
+ g++;
+
+ if (ind.empty())
+ r += 3;
+ else {
+ *lind0 = *lind;
+ lind++;
+ r += 3 * (*lind - *lind0);
+ }
+ } // end i loop
+
+ // condition: sum of weights = 1
+ *g = 1.0; // g0 [n1x1]
+
+ return g0;
+}
+
+// compute result to be located on reference point
+double result(std::vector<int> ind, struct points *pnts, mat_struct *w0)
+{
+ int n;
+ double *vo;
+ mat_struct *loc_w = w0;
+
+ int *indo, *indo0, *lind;
+
+ if (ind.empty())
+ n = pnts->n;
+ else {
+ n = ind.size ();
+ lind = &ind[0];
+ }
+
+ int i;
+ mat_struct *ins, *w, *rslt_OK;
+ doublereal *vt, *wo, *wt;
+ ins = G_matrix_init(n, 1, n);
+ w = G_matrix_init(1, n, 1);
+
+ if (ind.empty())
+ vo = &pnts->invals[0];// original input values
+ else {
+ indo = &lind[0]; // original indexes
+ indo0 = &lind[0];
+ vo = &pnts->invals[*indo];// original input values
+ }
+
+ vt = &ins->vals[0]; // selected inputs values
+ wo = &loc_w->vals[0]; // selected inputs values
+ wt = &w->vals[0]; // weights without last one
+
+ for (i=0; i<n; i++) {
+ *vt = *vo;
+ *wt = *wo;
+ if (i<n-1) {
+ indo0 = indo;
+ indo++;
+ if (ind.empty())
+ vo++;
+ else
+ vo += *indo - *indo0;
+ vt++;
+ wo++;
+ wt++;
+ }
+ }
+ rslt_OK = G_matrix_product(w,ins);
+
+ G_matrix_free(w);
+ G_matrix_free(ins);
+
+ return rslt_OK->vals[0];
+}
+
+// validation
+void crossvalidation(struct int_par *xD, struct points *pnts, pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts, struct parameters *var_par)
+{
+ int n = pnts->n; // # of input points
+ double *r = pnts->r; // coordinates of points
+ double *vals = pnts->invals;
+ double ratio;
+ mat_struct *GM = var_par->GM;
+
+ ratio = var_par->type == 3 ? xD->aniso_ratio : 1.;
+
+ int i, j;
+ int n1;
+ double radius;
+ double r0[3];
+ mat_struct *GM_sub;
+ mat_struct *GM_Inv, *g0, *w0;
+ double rslt_OK;
+ pcl::KdTreeFLANN<pcl::PointXYZ> kd_tree;
+ kd_tree.setInputCloud(pcl_pnts); // Set up kd-tree
+ pcl::PointXYZ cellCent;
+ std::vector<int> ind;
+ std::vector<float> sqDist;
+ struct write *crossvalid = &xD->crossvalid;
+ struct write *report = &xD->report;
+ double *normal, *absval, *norm, *av;
+ normal = (double *) G_malloc(n * sizeof(double));
+ absval = (double *) G_malloc(n * sizeof(double));
+ norm = &normal[0];
+ av = &absval[0];
+
+ FILE *fp;
+ if (crossvalid->name)
+ fp = fopen(crossvalid->name, "w");
+
+ radius = var_par->max_dist;
+
+ for (i=0; i<n; i++) {
+ // set up processing point
+ *r0 = cellCent.x = *r;
+ *(r0+1) = cellCent.y = *(r+1);
+ *(r0+2) = cellCent.z = *(r+2); // if 2D: 0.
+
+ if (kd_tree.radiusSearch(cellCent, radius, ind, sqDist) > 0 ) {
+ n1 = ind.size () + 1; // number of rellevant points + weights
+ GM_sub = submatrix(ind, GM, &xD->report);
+ GM_Inv = G_matrix_inverse(GM_sub);
+ G_matrix_free(GM_sub);
+
+ g0 = set_up_g0(xD, ind, pnts, r0, 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);
+ rslt_OK = result(ind, pnts, w0); // Estimated cell/voxel value rslt_OK = w x inputs
+ G_matrix_free(w0);
+
+ //Create output
+ *norm = *vals - rslt_OK;
+ *av = fabsf(*norm);
+ if (xD->i3 == TRUE)
+ fprintf(fp,"%d %.3f %.3f %.2f %f %f %f\n", i, *r, *(r+1), *(r+2) / ratio, pnts->invals[i], rslt_OK, *norm);
+ else
+ fprintf(fp,"%d %.3f %.3f %f %f %f\n", i, *r, *(r+1), *vals, rslt_OK, *norm);
+ }
+ else
+ fprintf(fp,"%d. point does not have neighbours in given radius\n", i);
+ r += 3;
+ vals++;
+ norm++;
+ av++;
+ }
+ fclose(fp);
+ G_message(_("Cross validation results have been written into <%s>"), crossvalid->name);
+
+ if (report->name) {
+ double quant05, quant10, quant25, quant50, quant75, quant90, quant95;
+ fprintf(report->fp, "\n************************************************\n");
+ fprintf(report->fp, "*** Cross validation results ***\n");
+
+ test_normality(n, normal, report);
+
+ fprintf(report->fp, "Quantile of absolute values\n");
+ quant95 = quantile(0.95, n, absval, report);
+ //quant90 = quantile(0.90, n, absval, report);
+ //quant75 = quantile(0.75, n, absval, report);
+ //quant50 = quantile(0.50, n, absval, report);
+ //quant25 = quantile(0.25, n, absval, report);
+ //quant10 = quantile(0.10, n, absval, report);
+ //quant05 = quantile(0.05, n, absval, report);
+ }
+}
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_raster.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_raster.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_raster.cpp 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,61 @@
+#include "local_proto.h"
+
+extern "C" {
+ /* open raster layer */
+ void open_layer(struct int_par *xD, struct reg_par *reg, struct output *out)
+ {
+ struct write *report = &xD->report;
+ /* 2D Raster layer */
+ if (xD->i3 == FALSE) {
+ out->dcell = (DCELL *) Rast_allocate_buf(DCELL_TYPE);
+ out->fd_2d = Rast_open_new(out->name, DCELL_TYPE); /* open output raster */
+ if (out->fd_2d < 0) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to create 2D raster <%s>"), out->name);
+ }
+ }
+ /* 3D Raster layer */
+ if (xD->i3 == TRUE) {
+ out->fd_3d = (RASTER3D_Map *) Rast3d_open_cell_new(out->name, DCELL_TYPE, RASTER3D_USE_CACHE_XYZ, ®->reg_3d); // initialize pointer to 3D region
+ if (out->fd_3d == NULL) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to create 3D raster <%s>"), out->name);
+ }
+ }
+ }
+
+ int write2layer(struct int_par *xD, struct reg_par *reg, struct output *out, unsigned int col, unsigned int row, unsigned int dep, double rslt_OK)
+ {
+ // Local variables
+ int i3 = xD->i3;
+ struct write *report = &xD->report;
+
+ int pass = 0; /* Number of processed cells */
+ switch (i3) {
+ case FALSE: /* set value to cell (2D) */
+ out->dcell[col] = (DCELL) (rslt_OK);
+ pass++;
+ if (col == reg->ncols-1)
+ Rast_put_row(out->fd_2d, out->dcell, DCELL_TYPE);
+ break;
+
+ case TRUE: /* set value to voxel (based on part of r3.gwflow (Soeren Gebbert)) */
+ if (Rast3d_put_double(out->fd_3d, col, row, dep, rslt_OK) == 0) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Error writing cell (%d,%d,%d) with value %f, nrows = %d"), row, col, dep, rslt_OK, reg->nrows);
+ }
+ pass++;
+ break;
+ }
+ return pass;
+ }
+}
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_write.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_write.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/utils_write.cpp 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,305 @@
+#include "local_proto.h"
+
+void variogram_type(int code, char *type)
+{
+ switch (code) {
+ case 0:
+ strcpy(type, "horizontal");
+ break;
+ case 1:
+ strcpy(type, "vertical");
+ break;
+ case 2:
+ strcpy(type, "bivariate");
+ break;
+ case 3:
+ strcpy(type, "anisotropic");
+ break;
+ }
+}
+
+void write2file_basics(struct int_par *xD, struct opts *opt)
+{
+ struct write *report = &xD->report;
+
+ fprintf(report->fp, "************************************************\n");
+ fprintf(report->fp, "*** Input ***\n\n");
+ fprintf(report->fp, "- vector layer: %s\n", opt->input->answer);
+ switch (xD->v3) {
+ case TRUE:
+ fprintf(report->fp, "- is 3D: yes\n");
+ break;
+ case FALSE:
+ fprintf(report->fp, "- is 3D: no\n");
+ break;
+ }
+ // todo: move to end
+ fprintf(report->fp, "\n");
+ fprintf(report->fp, "*** Output *** \n\n");
+ switch (xD->i3) {
+ case FALSE:
+ fprintf(report->fp, "- raster layer: %s\n", opt->output->answer);
+ fprintf(report->fp, "- is 3D: no\n");
+ break;
+ case TRUE:
+ fprintf(report->fp, "- volume layer: %s\n", opt->output->answer);
+ fprintf(report->fp, "- is 3D: yes\n");
+ break;
+ }
+}
+
+void write2file_vector( struct int_par *xD, struct points *pnts)
+{
+ struct select in_reg = pnts->in_reg;
+ struct write *report = &xD->report;
+ fprintf(report->fp, "\n");
+ fprintf(report->fp, "************************************************\n");
+ fprintf(report->fp, "*** Input vector layer properties ***\n\n");
+ fprintf(report->fp, "# of points (total): %d\n", in_reg.total);
+ fprintf(report->fp, "# of points (in region): %d\n", in_reg.n);
+ fprintf(report->fp, "# of points (out of region): %d\n\n", in_reg.out);
+ fprintf(report->fp, "- extent:\n");
+ switch (xD->i3) {
+ case TRUE:
+ fprintf(report->fp, "xmin=%f ymin=%f zmin=%f\n", pnts->r_min[0], pnts->r_min[1], pnts->r_min[2]);
+ fprintf(report->fp, "xmax=%f ymax=%f zmax=%f\n", pnts->r_max[0], pnts->r_max[1], pnts->r_max[2]);
+ break;
+ case FALSE:
+ fprintf(report->fp, "xmin=%f ymin=%f\n", pnts->r_min[0], pnts->r_min[1]);
+ fprintf(report->fp, "xmax=%f ymax=%f\n", pnts->r_max[0], pnts->r_max[1]);
+ break;
+ }
+}
+
+void write2file_values(struct write *report, const char *column)
+{
+ fprintf(report->fp, "\n");
+ fprintf(report->fp, "************************************************\n");
+ fprintf(report->fp, "*** Values to be interpolated ***\n\n");
+ fprintf(report->fp, "- attribute column: %s\n", column);
+}
+
+void write2file_varSetsIntro(int code, struct write *report)
+{
+ char type[12];
+ variogram_type(code, type);
+
+ fprintf(report->fp, "\n************************************************\n");
+ fprintf(report->fp, "*** Variogram settings - %s ***\n\n", type);
+}
+
+void write2file_varSets(struct write *report, struct parameters *var_par)
+{
+ char type[12];
+ variogram_type(var_par->type, type);
+
+ fprintf(report->fp, "- number of lags in %s direction: %d\n", type, var_par->nLag);
+ fprintf(report->fp, "- lag distance (%s): %f\n", type, var_par->lag);
+
+ if (var_par->function == 5) {
+ fprintf(report->fp, "- number of lags in vertical direction: %d\n", var_par->nLag_vert);
+ fprintf(report->fp, "- lag distance (vertical): %f\n", var_par->lag);
+ }
+
+ fprintf(report->fp, "- azimuth: %f°\n", RAD2DEG(var_par->dir));
+ fprintf(report->fp, "- angular tolerance: %f°\n", RAD2DEG(var_par->td));
+ //fprintf(report->fp, "- variogram plotted: %s\n", );
+
+ switch (var_par->function) {
+ case 5:
+ fprintf(report->fp, "- maximum distance in horizontal direction (1/3): %f\n", var_par->max_dist);
+ break;
+ default:
+ fprintf(report->fp, "- maximum distance (1/3): %f\n", var_par->max_dist);
+ break;
+ }
+}
+
+void write2file_variogram_E(struct int_par *xD, struct parameters *var_par)
+{
+ int i, j;
+ int fction = var_par->function;
+ double *h, *vert;
+ double *c = var_par->c->vals;
+ double *gamma = var_par->gamma->vals;
+ struct write *report = &xD->report;
+
+ char type[12];
+ variogram_type(var_par->type, type);
+
+ fprintf(report->fp, "\n************************************************\n");
+ fprintf(report->fp, "*** Experimental variogram - %s ***\n\n", type);
+
+ if (fction == 5) { // if variogram is bivariate
+ vert = &var_par->vert[0]; // use vertical variable
+ for (i=0; i<var_par->nLag_vert; i++) { // write header - verts
+ if (i == 0)
+ fprintf(report->fp, " lagV ||"); // column for h
+ fprintf(report->fp, " %f ||", *vert);
+ //fprintf(report->fp, " # of pairs | average ||");
+ //fprintf(report->fp, "------------------------");
+ vert++;
+ }
+ fprintf(report->fp, "\n");
+
+ for (i=0; i<var_par->nLag_vert; i++) { // write header - h c gamma
+ if (i == 0)
+ fprintf(report->fp, " lagHZ ||"); // column for h
+ fprintf(report->fp, " c | ave ||", *vert);
+ }
+ fprintf(report->fp, "\n");
+
+ for (i=0; i<var_par->nLag_vert; i++) { // write header - h c gamma
+ if (i == 0)
+ fprintf(report->fp, "--------||"); // column for h
+ fprintf(report->fp, "-----------", *vert);
+ }
+ fprintf(report->fp, "\n");
+ }
+ else {
+ fprintf(report->fp, " lag || # of pairs | average\n");
+ fprintf(report->fp, "------------------------------------\n");
+ }
+
+ // Write values
+ h = &var_par->h[0];
+ for (i=0; i<var_par->nLag; i++) {
+ fprintf(report->fp, "%f ||", *h);
+ if (fction == 5) {
+ //for (j=0; j<var_par->nZ; j++) {
+ fprintf(report->fp, " %d | %f ||", (int) *c, *gamma);
+ c++;
+ gamma++;
+ //} // end for j
+ fprintf(report->fp, "\n");
+ //for (j=0; j<var_par->nZ; j++)
+ fprintf(report->fp, "-----------", (int) *c, *gamma);
+ fprintf(report->fp, "\n");
+ }
+ else {
+ fprintf(report->fp, " %d | %f\n", (int) *c, *gamma);
+ c++;
+ gamma++;
+ }
+ h++;
+ }
+ fprintf(report->fp, "------------------------------------\n");
+}
+
+void write2file_variogram_T(struct write *report)
+{
+ fprintf(report->fp, "\n");
+ fprintf(report->fp, "************************************************\n");
+ fprintf(report->fp, "*** Theoretical variogram ***\n\n");
+}
+
+void write_temporary2file(struct int_par *xD, struct parameters *var_par, mat_struct *gammaMT)
+{
+ // local variables
+ int type = var_par->type;
+ int nLag = var_par->nLag;
+ int nLag_vert;
+ double *h = var_par->h;
+ double *vert;
+ double *c = var_par->c->vals;
+ double *gamma = gammaMT->vals;
+ double sill = var_par->sill;
+
+ int i, j; // index
+ FILE *fp;
+ int file_length;
+
+ switch (type) {
+ case 0: // horizontal variogram
+ fp = fopen("variogram_hz_tmp.txt", "w");
+ if (xD->report.name) { // write name of report file
+ file_length = strlen(xD->report.name);
+ if (file_length < 4) // 4 types of variogram
+ G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ }
+ fprintf(fp,"%d\n", var_par->type); // write # of lags
+ break;
+ case 1: // vertical variogram
+ fp = fopen("variogram_vert_tmp.txt", "w");
+ if (xD->report.name) { // write name of report file
+ file_length = strlen(xD->report.name);
+ if (file_length < 3)
+ G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ }
+ fprintf(fp,"%d\n", var_par->type); // write type
+ break;
+ case 2:
+ fp = fopen("variogram_final_tmp.txt", "w");
+ if (xD->report.name) { // write name of report file
+ file_length = strlen(xD->report.name);
+ if (file_length < 4) // 4 types of variogram
+ G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ }
+
+ fprintf(fp,"%d\n", var_par->type); // write # of lags
+ fprintf(fp,"%d\n", var_par->nLag_vert); // write # of lags
+ fprintf(fp,"%f\n", var_par->lag_vert); // write size of lag
+ fprintf(fp,"%f\n", var_par->max_dist_vert); // write maximum distance
+ break;
+ case 3: // anisotropic variogram
+ fp = fopen("variogram_final_tmp.txt", "w");
+ if (xD->report.name) { // write name of report file
+ file_length = strlen(xD->report.name);
+ if (file_length < 4) // 4 types of variogram
+ G_fatal_error(_("File name must contain more than 4 characters...")); // todo: error
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ }
+ fprintf(fp,"%d\n", var_par->type); // write type
+ fprintf(fp,"%f\n", xD->aniso_ratio); // write ratio of anisotropy
+ break;
+ }
+
+ fprintf(fp,"%d\n", nLag); // write # of lags
+ fprintf(fp,"%f\n", var_par->lag); // write size of lag
+ fprintf(fp,"%f\n", var_par->max_dist); // write maximum distance
+ if (type != 1)
+ fprintf(fp,"%f\n", var_par->td); // write maximum distance
+
+ switch (type) {
+ case 2:
+ nLag_vert = var_par->nLag_vert;
+ // write h
+ for (i=0; i<nLag; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *h);
+ h++;
+ }
+ // write vert
+ vert = &var_par->vert[0];
+ for (i=0; i<nLag_vert; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *vert);
+ vert++;
+ }
+ // write c
+ for (i=0; i<nLag * nLag_vert; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *c);
+ c++;
+ }
+ // write gamma
+ for (i=0; i<nLag * nLag_vert; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *gamma);
+ gamma++;
+ }
+ fprintf(fp,"%f\n", var_par->horizontal.sill);// write sill
+ fprintf(fp,"%f\n", var_par->vertical.sill);// write sill
+ break;
+ default:
+ for (i=0; i<nLag; i++) { // write experimental variogram
+ fprintf(fp,"%f %f %f\n", *h, *c, *gamma);
+ h++;
+ c++;
+ gamma++;
+ }
+ fprintf(fp,"%f", sill);// write sill
+ break;
+ }
+
+ fclose(fp);
+}
Added: grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/v.kriging.html
===================================================================
--- grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/v.kriging.html (rev 0)
+++ grass-addons/grass7/vector/v.kriging/backup_version_with_PCL/v.kriging.html 2015-06-20 12:01:04 UTC (rev 65505)
@@ -0,0 +1,60 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.kriging</em> interpolates unknown values using method <i>ordinary
+kriging</i>. Output can be 2D or 3D.
+
+<h2>EXAMPLES</h2> Input layer should contain 3D coordinates (xyz) and
+values to be interpolated (in attribute table). The commands can look
+like this:
+
+<div class="code"><pre>
+v.kriging phase=initial in=input_layer icol=name report=report_file.txt file=png
+</pre></div>
+<div class="code"><pre>
+v.kriging in=input_layer phase=middle hz_fun=exponential vert_fun=exponential \
+ ic=name file=png hz_nugget=0 vert_nugget=0 hz_range=double vert_range=double -b
+</pre></div>
+<div class="code"><pre>
+v.kriging in=input_layer phase=final final_fun=bivariate icol=name file=png \
+ out=name crossval=crossval_file.txt
+</pre></div>
+
+<h2>TODO</h2>
+<ul>
+<li><b>anisotropy</b> in horizontal direction missing
+<li>current version is suitable just for <b>metric coordinate systems</b>
+<li>enable <b>mask usage</b>
+<li><b>bivariate variogram</b> needs to be rebuilt (theory)
+<li><b>2D interpolation from 3D input layer</b> needs to be rebuilt (especially in case that there are too many points located on identical horizontal coordinates with different elevation)
+</ul>
+
+<h2>BUGS</h2>
+<ul>
+<li>2D variogram gives quite inaccurate results for <b>sparse (or spatially heterogeneous) data</b>
+<li><b>bivariate variogram</b> - too inaccurate results
+</ul>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="v.vol.rst.html">v.vol.rst</a><br>
+<a href="v.krige.html">v.krige</a>
+</em>
+
+<h2>REQUIREMENTS</h2>
+<ul>
+<li><b>Gnuplot</b> graphing utility, <a href="http://www.gnuplot.info/">more</a><br>
+<li><b>LAPACK / BLAS</b> (libraries for numerical computing) for
+ GMATH library (GRASS Numerical Library)<br>
+ <a href="http://www.netlib.org/lapack">http://www.netlib.org/lapack</a>
+ (usually available on Linux distros)
+</ul>
+
+<h2>AUTHOR</h2>
+
+Eva Stopkova<br>
+functions taken from another modules are cited above the function or at
+the beginning of the file (e.g. <i>quantile.cpp</i> that uses slightly
+modified functions taken from the module <i>r.quantile</i> (Clemens, G.))
+<p>
+<i>Last changed: $Date: 2015-06-20 13:23:47 +0200 (Sat, 20 Jun 2015) $</i>
More information about the grass-commit
mailing list