[GRASS-SVN] r53175 - in grass/branches/releasebranch_6_4/vector/lidar: lidarlib v.lidar.correction v.lidar.edgedetection v.outlier v.surf.bspline
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Sep 12 12:11:17 PDT 2012
Author: neteler
Date: 2012-09-12 12:11:17 -0700 (Wed, 12 Sep 2012)
New Revision: 53175
Modified:
grass/branches/releasebranch_6_4/vector/lidar/lidarlib/PolimiFunct.h
grass/branches/releasebranch_6_4/vector/lidar/lidarlib/TcholBand.c
grass/branches/releasebranch_6_4/vector/lidar/lidarlib/zones.c
grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/main.c
grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/main.c
grass/branches/releasebranch_6_4/vector/lidar/v.outlier/main.c
grass/branches/releasebranch_6_4/vector/lidar/v.surf.bspline/crosscorr.c
grass/branches/releasebranch_6_4/vector/lidar/v.surf.bspline/main.c
Log:
Backport of improvements from G6.5.svn; revert r53173 as no longer needed
Modified: grass/branches/releasebranch_6_4/vector/lidar/lidarlib/PolimiFunct.h
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/lidarlib/PolimiFunct.h 2012-09-12 19:03:34 UTC (rev 53174)
+++ grass/branches/releasebranch_6_4/vector/lidar/lidarlib/PolimiFunct.h 2012-09-12 19:11:17 UTC (rev 53175)
@@ -161,13 +161,13 @@
/*----------------------------------------------------------------------------------------------------------*/
/*tcholBand */
-void tcholDec(double **N, double **T, int n, int BW, int CV);
-void tcholSolve(double **N, double *TN, double *parVect, int n, int BW, int CV);
+void tcholDec(double **N, double **T, int n, int BW);
+void tcholSolve(double **N, double *TN, double *parVect, int n, int BW);
void tcholSolve2(double **N, double *TN, double **T, double *parVect, int n,
int BW);
-void tcholInv(double **N, double *invNdiag, int n, int BW, int CV);
+void tcholInv(double **N, double *invNdiag, int n, int BW);
void tcholSolveInv(double **N, double *TN, double *invNdiag, double *parVect,
- int n, int BW, int CV);
+ int n, int BW);
/*---------------------------------------------------------------------------------------*/
/*interpSpline */
Modified: grass/branches/releasebranch_6_4/vector/lidar/lidarlib/TcholBand.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/lidarlib/TcholBand.c 2012-09-12 19:03:34 UTC (rev 53174)
+++ grass/branches/releasebranch_6_4/vector/lidar/lidarlib/TcholBand.c 2012-09-12 19:11:17 UTC (rev 53175)
@@ -7,7 +7,7 @@
/*--------------------------------------------------------------------------------------*/
/* Cholesky decomposition -> T= Lower Triangular Matrix */
-void tcholDec(double **N, double **T, int n, int BW, int CV)
+void tcholDec(double **N, double **T, int n, int BW)
{
int i, j, k, end;
double somma;
@@ -15,7 +15,7 @@
G_debug(2, "tcholDec(): n=%d BW=%d", n, BW);
for (i = 0; i < n; i++) {
- if (CV == 0) G_percent(i, n, 2);
+ G_percent(i, n, 2);
for (j = 0; j < BW; j++) {
somma = N[i][j];
/* start = 1 */
@@ -33,14 +33,14 @@
}
}
- if (CV == 0) G_percent(i, n, 2);
+ G_percent(i, n, 2);
return;
}
/*--------------------------------------------------------------------------------------*/
/* Cholesky matrix solution */
-void tcholSolve(double **N, double *TN, double *parVect, int n, int BW, int CV)
+void tcholSolve(double **N, double *TN, double *parVect, int n, int BW)
{
double **T;
@@ -50,7 +50,7 @@
T = G_alloc_matrix(n, BW);
/*--------------------------------------*/
- tcholDec(N, T, n, BW, CV); /* T computation */
+ tcholDec(N, T, n, BW); /* T computation */
/* Forward substitution */
parVect[0] = TN[0] / T[0][0];
@@ -121,7 +121,7 @@
/*--------------------------------------------------------------------------------------*/
/* Cholesky matrix inversion */
-void tcholInv(double **N, double *invNdiag, int n, int BW, int CV)
+void tcholInv(double **N, double *invNdiag, int n, int BW)
{
double **T = NULL;
double *vect = NULL;
@@ -133,7 +133,7 @@
vect = G_alloc_vector(n);
/* T computation */
- tcholDec(N, T, n, BW, CV);
+ tcholDec(N, T, n, BW);
/* T Diagonal invertion */
for (i = 0; i < n; i++) {
@@ -168,7 +168,7 @@
/* Cholesky matrix solution and inversion */
void tcholSolveInv(double **N, double *TN, double *invNdiag, double *parVect,
- int n, int BW, int CV)
+ int n, int BW)
{
double **T = NULL;
@@ -181,7 +181,7 @@
vect = G_alloc_vector(n);
/* T computation */
- tcholDec(N, T, n, BW, CV);
+ tcholDec(N, T, n, BW);
/* Forward substitution */
parVect[0] = TN[0] / T[0][0];
Modified: grass/branches/releasebranch_6_4/vector/lidar/lidarlib/zones.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/lidarlib/zones.c 2012-09-12 19:03:34 UTC (rev 53174)
+++ grass/branches/releasebranch_6_4/vector/lidar/lidarlib/zones.c 2012-09-12 19:11:17 UTC (rev 53175)
@@ -226,7 +226,7 @@
{
/* Returns the interpolation matrixes BandWidth dimension */
- if (interpolator == 1) {
+ if (interpolator == P_BILINEAR) {
return (2 * nsplines + 1);
}
else {
Modified: grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/main.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/main.c 2012-09-12 19:03:34 UTC (rev 53174)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/main.c 2012-09-12 19:11:17 UTC (rev 53175)
@@ -406,7 +406,7 @@
elaboration_reg.south, nterrain, nparameters,
BW);
nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
- tcholSolve(N, TN, parVect, nparameters, BW, 0);
+ tcholSolve(N, TN, parVect, nparameters, BW);
G_free_matrix(N);
G_free_vector(TN);
Modified: grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/main.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/main.c 2012-09-12 19:03:34 UTC (rev 53174)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/main.c 2012-09-12 19:11:17 UTC (rev 53175)
@@ -420,7 +420,7 @@
elaboration_reg.south, npoints, nparameters,
BW);
nCorrectGrad(N, lambda_B, nsplx, nsply, passoE, passoN);
- tcholSolve(N, TN, parVect_bilin, nparameters, BW, 0);
+ tcholSolve(N, TN, parVect_bilin, nparameters, BW);
G_free_matrix(N);
for (tn = 0; tn < nparameters; tn++)
@@ -437,7 +437,7 @@
elaboration_reg.south, npoints, nparameters,
BW);
nCorrectLapl(N, lambda_F, nsplx, nsply, passoE, passoN);
- tcholSolve(N, TN, parVect_bicub, nparameters, BW, 0);
+ tcholSolve(N, TN, parVect_bicub, nparameters, BW);
G_free_matrix(N);
G_free_vector(TN);
Modified: grass/branches/releasebranch_6_4/vector/lidar/v.outlier/main.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.outlier/main.c 2012-09-12 19:03:34 UTC (rev 53174)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.outlier/main.c 2012-09-12 19:11:17 UTC (rev 53175)
@@ -396,7 +396,7 @@
elaboration_reg.south, npoints, nparameters,
BW);
nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
- tcholSolve(N, TN, parVect, nparameters, BW, 0);
+ tcholSolve(N, TN, parVect, nparameters, BW);
G_free_matrix(N);
G_free_vector(TN);
Modified: grass/branches/releasebranch_6_4/vector/lidar/v.surf.bspline/crosscorr.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.surf.bspline/crosscorr.c 2012-09-12 19:03:34 UTC (rev 53174)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.surf.bspline/crosscorr.c 2012-09-12 19:11:17 UTC (rev 53175)
@@ -36,7 +36,7 @@
/*-------------------------------------------------------------------------------------------*/
int cross_correlation(struct Map_info *Map, double passWE, double passNS)
/*
- Map: Map in which cross crorrelation will take values
+ Map: Vector map from which cross-crorrelation will take values
passWE: spline step in West-East direction
passNS: spline step in North-South direction
@@ -49,8 +49,11 @@
int nsplx, nsply, nparam_spl, ndata;
double *mean, *rms, *stdev, rms_min, stdev_min;
- /* double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; */ /* Fixed values (by the moment) */
- double lambda[PARAM_LAMBDA] = { 0.01, 0.05, 0.1, 0.2, 0.3, 0.4 }; /* Fixed values (by the moment) */
+ /* double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; */ /* Fixed values (at the moment) */
+ double lambda[PARAM_LAMBDA] = { 0.01, 0.05, 0.1, 0.2, 0.3, 0.4 }; /* Fixed values (at the moment) */
+ /* a more exhaustive search:
+ #define PARAM_LAMBDA 11
+ double lambda[PARAM_LAMBDA] = { 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1.0, 5.0, 10.0 }; */
double *TN, *Q, *parVect; /* Interpolation and least-square vectors */
double **N, **obsVect; /* Interpolation and least-square matrix */
@@ -97,7 +100,7 @@
int BW, lbd_min; /* lbd_min: index where minimun is found */
double mean_reg, *obs_mean;
- int nrec, ctype = 0;
+ int nrec, ctype = 0, verbosity;
struct field_info *Fi;
dbDriver *driver_cats;
@@ -105,6 +108,8 @@
rms = G_alloc_vector(PARAM_LAMBDA); /* number of parameter used used for cross validation */
stdev = G_alloc_vector(PARAM_LAMBDA);
+ verbosity = G_verbose(); /* store for later reset */
+
/* Working with attributes */
if (bspline_field > 0) {
db_CatValArray_init(&cvarr);
@@ -159,7 +164,7 @@
TN = G_alloc_vector(nparam_spl); /* vector */
parVect = G_alloc_vector(nparam_spl); /* Parameters vector */
obsVect = G_alloc_matrix(ndata, 3); /* Observation vector */
- Q = G_alloc_vector(ndata); /* "a priori" var-cov matrix */
+ Q = G_alloc_vector(ndata); /* "a priori" var-cov matrix */
obs_mean = G_alloc_vector(ndata);
stat_vect = alloc_Stats(ndata);
@@ -167,22 +172,23 @@
for (lbd = 0; lbd < PARAM_LAMBDA; lbd++) { /* For each lambda value */
G_message(_("Beginning cross validation with "
- "lambda_i=%.4f..."), lambda[lbd]);
-
+ "lambda_i=%.4f ... (%d of %d)"), lambda[lbd],
+ lbd+1, PARAM_LAMBDA);
+
/*
- How cross correlation algorithm is done:
- For each cicle, only the first ndata-1 "observ" elements are considered for the
+ How the cross correlation algorithm is done:
+ For each cycle, only the first ndata-1 "observ" elements are considered for the
interpolation. Within every interpolation mean is calculated to lowering border
- errors. The point let out will be used for an estimation. The error between the
+ errors. The point left out will be used for an estimation. The error between the
estimation and the observation is recorded for further statistics.
- At the end of the cicle, the last point, that is, the ndata-1 index, and the point
+ At the end of the cycle, the last point, that is, the ndata-1 index, and the point
with j index are swapped.
*/
for (j = 0; j < ndata; j++) { /* Cross Correlation will use all ndata points */
- double out_x, out_y, out_z; /* This point is let out */
+ double out_x, out_y, out_z; /* This point is left out */
for (i = 0; i < ndata; i++) { /* Each time, only the first ndata-1 points */
- double dval; /* are considered in the interpolation */
+ double dval; /* are considered in the interpolation */
/* Setting obsVect vector & Q matrix */
Q[i] = 1; /* Q=I */
@@ -233,7 +239,7 @@
for (i = 0; i < ndata; i++)
obsVect[i][2] -= mean_reg;
- /* This is let out */
+ /* This is left out */
out_x = observ[ndata - 1].coordX;
out_y = observ[ndata - 1].coordY;
out_z = obsVect[ndata - 1][2];
@@ -257,7 +263,9 @@
if (bilin) interpolation (&interp, P_BILINEAR);
else interpolation (&interp, P_BICUBIC);
*/
- tcholSolve(N, TN, parVect, nparam_spl, BW, 1);
+ G_set_verbose(G_verbose_min());
+ tcholSolve(N, TN, parVect, nparam_spl, BW);
+ G_set_verbose(verbosity);
/* Estimation of j-point */
if (bilin)
@@ -272,12 +280,15 @@
nsplx, nsply, region.west,
region.south, parVect);
- /*Difference between estimated and observated i-point */
+ /* Difference between estimated and observated i-point */
stat_vect.error[j] = out_z - stat_vect.estima[j];
G_debug(1, "CrossCorrelation: stat_vect.error[%d] = %lf", j,
stat_vect.error[j]);
- observ = swap(observ, j, ndata - 1); /* Once the last value is let out, it is swap with j-value */
+ /* Once the last value is left out, it is swaped with j-value */
+ observ = swap(observ, j, ndata - 1);
+
+ G_percent(j, ndata, 2);
}
mean[lbd] = calc_mean(stat_vect.error, stat_vect.n_points);
@@ -287,8 +298,9 @@
calc_standard_deviation(stat_vect.error, stat_vect.n_points);
G_message(_("Mean = %.5lf"), mean[lbd]);
- G_message(_("Root Means Square (RMS) = %.5lf"),
+ G_message(_("Root Mean Square (RMS) = %.5lf"),
rms[lbd]);
+ G_message("---");
} /* ENDFOR each lambda value */
G_free_matrix(N);
@@ -298,7 +310,7 @@
G_free_vector(parVect);
#ifdef nodef
/*TODO: if the minimum lambda is wanted, the function declaration must be changed */
- /*At this moment, consider rms only */
+ /* At this moment, consider rms only */
rms_min = find_minimum(rms, &lbd_min);
stdev_min = find_minimum(stdev, &lbd_min);
@@ -313,11 +325,11 @@
*lambda_min = lambda[lbd_min];
#endif
- G_message(_("Results into a table:"));
- G_message(_(" lambda | mean | rms |"));
+ G_message(_("Table of results:"));
+ fprintf(stdout, _(" lambda | mean | rms |\n"));
for (lbd = 0; lbd < PARAM_LAMBDA; lbd++) {
- G_message(_(" %-10.5f| %-12.4f| %-12.4f|"), lambda[lbd],
- mean[lbd], rms[lbd]);
+ fprintf(stdout, " %9.5f | %10.4f | %10.4f |\n", lambda[lbd],
+ mean[lbd], rms[lbd]);
}
G_free_vector(mean);
@@ -433,8 +445,8 @@
struct Point *swap(struct Point *point, int a, int b)
{
-
- SWAP(point[a].coordX, point[b].coordX); /* Once the last value is let out, it is swap with j-value */
+ /* Once the last value is left out, it is swaped with j-value */
+ SWAP(point[a].coordX, point[b].coordX);
SWAP(point[a].coordY, point[b].coordY);
SWAP(point[a].coordZ, point[b].coordZ);
SWAP(point[a].cat, point[b].cat);
Modified: grass/branches/releasebranch_6_4/vector/lidar/v.surf.bspline/main.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.surf.bspline/main.c 2012-09-12 19:03:34 UTC (rev 53174)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.surf.bspline/main.c 2012-09-12 19:11:17 UTC (rev 53175)
@@ -263,13 +263,8 @@
}
/* Open input ext vector */
- if (!in_ext_opt->answer) {
- ext = FALSE;
- G_message(_("No vector map of sparse points to interpolate was specified. "
- "Interpolation will be done with <%s> vector map"),
- in_opt->answer);
- }
- else {
+ ext = FALSE;
+ if (in_ext_opt->answer) {
ext = TRUE;
G_message(_("Vector map <%s> of sparse points will be interpolated"),
in_ext_opt->answer);
@@ -308,6 +303,9 @@
Vect_hist_copy(&In_ext, &Out);
}
Vect_hist_command(&Out);
+
+ G_message(_("Points in input vector map <%s> will be interpolated"),
+ vector);
}
/* raster output */
@@ -318,6 +316,9 @@
if ((raster = G_open_fp_cell_new(out_map_opt->answer)) < 0)
G_fatal_error(_("Unable to create raster map <%s>"),
out_map_opt->answer);
+
+ G_message(_("Cells for raster map <%s> will be interpolated"),
+ map);
}
/* read z values from attribute table */
@@ -473,6 +474,10 @@
while (last_column == FALSE) { /* For each subregion column */
int npoints = 0;
+ /* needed for sparse points interpolation */
+ int npoints_ext, *lineVect_ext = NULL;
+ double **obsVect_ext; /*, mean_ext = .0; */
+ struct Point *observ_ext;
subregion_col++;
subregion++;
@@ -513,14 +518,31 @@
/* reading points in interpolation region */
dim_vect = nsplx * nsply;
- observ =
- P_Read_Vector_Region_Map(&In, &elaboration_reg, &npoints,
- dim_vect, bspline_field);
+ if (grid == FALSE && ext == TRUE) {
+ observ_ext =
+ P_Read_Vector_Region_Map(&In_ext,
+ &elaboration_reg,
+ &npoints_ext, dim_vect,
+ 1);
+ }
+ else
+ npoints_ext = 1;
+
+ observ = NULL;
+ if (npoints_ext > 0) {
+ observ =
+ P_Read_Vector_Region_Map(&In, &elaboration_reg, &npoints,
+ dim_vect, bspline_field);
+ }
+ else
+ npoints = 0;
+
G_debug(1,
"Interpolation: (%d,%d): Number of points in <elaboration_box> is %d",
subregion_row, subregion_col, npoints);
- if (npoints > 0) { /* */
+ /* only interpolate if there are any points in current subregion */
+ if (npoints > 0 && npoints_ext > 0) {
int i;
nparameters = nsplx * nsply;
@@ -610,7 +632,7 @@
nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
}
- tcholSolve(N, TN, parVect, nparameters, BW, cross_corr_flag->answer);
+ tcholSolve(N, TN, parVect, nparameters, BW);
G_free_matrix(N);
G_free_vector(TN);
@@ -637,8 +659,11 @@
table_name);
}
else { /* FLAG_EXT == TRUE */
+
+ /* done that earlier */
+ /*
int npoints_ext, *lineVect_ext = NULL;
- double **obsVect_ext; /*, mean_ext = .0; */
+ double **obsVect_ext;
struct Point *observ_ext;
observ_ext =
@@ -646,6 +671,7 @@
&elaboration_reg,
&npoints_ext, dim_vect,
1);
+ */
obsVect_ext = G_alloc_matrix(npoints_ext, 3); /* Observation vector_ext */
lineVect_ext = G_alloc_ivector(npoints_ext);
@@ -677,9 +703,11 @@
G_free_ivector(lineVect);
}
else {
- G_free(observ);
- G_warning(_("No data within this subregion. "
- "Consider changing the spline step."));
+ if (observ)
+ G_free(observ);
+ if (npoints == 0)
+ G_warning(_("No data within this subregion. "
+ "Consider increasing spline step values."));
}
} /*! END WHILE; last_column = TRUE */
} /*! END WHILE; last_row = TRUE */
More information about the grass-commit
mailing list