code clean-up and more documentation

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

vector/lidar/lidarlib/PolimiFunct.h

@@ -33,7 +33,7 @@
 #define NSPLX_MAX 	      150	/* Maximum number of splines along East direction used in the subregions interpolation */
 #define NSPLY_MAX	      150	/* Maximum number of splines along North direction used in the subregions interpolation */
 #define OVERLAP_SIZE 	       10	/* Subregions overlapping size. */
-#define LATO 		     2000	/* Side's size for v.to.qrast command. */
+#define LATO 		     1000	/* Side's size for v.lidar.growing. */
 #define CONTOUR		       15 	/**/
 #define GENERAL_ROW 		0
 #define GENERAL_COLUMN 	        1
@@ -145,7 +145,9 @@ double **P_Regular_Points(struct Cell_head *, /**/
 			  int, /**/ int, /**/ int, /**/ int, /**/ int /**/);
 
 /*----------------------------------------------------------------------------------------------------------*/
-int P_Create_Aux_Table(dbDriver *, /**/ char * /**/);
+int P_Create_Aux2_Table(dbDriver *, /**/ char * /**/);
+
+int P_Create_Aux4_Table(dbDriver *, /**/ char * /**/);
 
 int P_Drop_Aux_Table(dbDriver *, /**/ char * /**/);
 

vector/lidar/lidarlib/raster.c

@@ -29,6 +29,8 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
     point = Vect_new_line_struct();
 
+    db_begin_transaction(driver);
+    
     for (i = 0; i < num_points; i++) {
 
 	if (Vect_point_in_box(obs[i][0], obs[i][1], mean, &General)) {	/*Here mean is just for asking if obs point is in box */
@@ -50,10 +52,8 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
 	    if (Vect_point_in_box(obs[i][0], obs[i][1], interpolation, &Overlap)) {	/*(5) */
 		Vect_write_line(Out, GV_POINT, point, categories);
-
 	    }
 	    else {
-
 		db_init_string(&sql);
 
 		sprintf(buf, "INSERT INTO %s (ID, X, Y, Interp)", tab_name);
@@ -65,12 +65,11 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 			obs[i][1]);
 		db_append_string(&sql, buf);
 
-		if ((point->x[0] > Overlap.E)) {
-
-		    if ((*point->y > Overlap.N)) {	/*(3) */
-			csi = (*point->x - Overlap.E) / overlap;
-			eta = (*point->y - Overlap.N) / overlap;
-			weight = (1 - csi) * (1 - eta);
+		if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(3) */
+			csi = (General.E - *point->x) / overlap;
+			eta = (General.N - *point->y) / overlap;
+			weight = csi * eta;
 			*point->z = weight * interpolation;
 
 			sprintf(buf, "%lf", *point->z);
@@ -79,14 +78,13 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 			db_append_string(&sql, buf);
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
-			    G_fatal_error(_("Unable to create table: %s"),
+			    G_fatal_error(_("Unable to write to table: %s"),
 					  buf);
-
 		    }
-		    else if ((*point->y < Overlap.S)) {	/*(1) */
-			csi = (*point->x - Overlap.E) / overlap;
-			eta = (Overlap.S - *point->y) / overlap;
-			weight = (1 - csi) * eta;
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(1) */
+			csi = (General.E - *point->x) / overlap;
+			eta = (*point->y - General.S) / overlap;
+			weight = csi * eta;
 			*point->z = weight * interpolation;
 
 			sprintf(buf, "%lf", *point->z);
@@ -95,13 +93,12 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 			db_append_string(&sql, buf);
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
-			    G_fatal_error(_("Unable to create table: %s"),
+			    G_fatal_error(_("Unable to write to table: %s"),
 					  buf);
-
 		    }
-		    else {	/*(1) */
-			weight = (*point->x - Overlap.E) / overlap;
-			*point->z = (1 - weight) * interpolation;
+		    else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) {	/*(1) */
+			weight = (General.E - *point->x) / overlap;
+			*point->z = weight * interpolation;
 
 			sprintf(buf, "%lf", *point->z);
 			db_append_string(&sql, buf);
@@ -109,16 +106,15 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 			db_append_string(&sql, buf);
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
-			    G_fatal_error(_("Unable to create table: %s"),
+			    G_fatal_error(_("Unable to write to table: %s"),
 					  buf);
 		    }
-
 		}
-		else if ((point->x[0] < Overlap.W)) {
-		    if ((*point->y > Overlap.N)) {	/*(4) */
-			csi = (Overlap.W - *point->x) / overlap;
-			eta = (*point->y - Overlap.N) / overlap;
-			weight = (1 - eta) * csi;
+		else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(4) */
+			csi = (*point->x - General.W) / overlap;
+			eta = (General.N - *point->y) / overlap;
+			weight = eta * csi;
 			*point->z = weight * interpolation;
 
 			sprintf(buf, "%lf", *point->z);
@@ -127,13 +123,12 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 			db_append_string(&sql, buf);
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
-			    G_fatal_error(_("Unable to create table: %s"),
+			    G_fatal_error(_("Unable to write to table: %s"),
 					  buf);
-
 		    }
-		    else if ((*point->y < Overlap.S)) {	/*(2) */
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(2) */
 			csi = (*point->x - General.W) / overlap;
-			eta = (Overlap.S - *point->y) / overlap;
+			eta = (*point->y - General.S) / overlap;
 			weight = csi * eta;
 			*point->z = weight * interpolation;
 
@@ -143,13 +138,12 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 			db_append_string(&sql, buf);
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
-			    G_fatal_error(_("Unable to create table: %s"),
+			    G_fatal_error(_("Unable to write to table: %s"),
 					  buf);
-
 		    }
-		    else {	/*(2) */
-			weight = (Overlap.W - *point->x) / overlap;
-			*point->z = (1 - weight) * interpolation;
+		    else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) {	/*(2) */
+			weight = (*point->x - General.W) / overlap;
+			*point->z = weight * interpolation;
 
 			sprintf(buf, "%lf", *point->z);
 			db_append_string(&sql, buf);
@@ -157,15 +151,14 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 			db_append_string(&sql, buf);
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
-			    G_fatal_error(_("Unable to create table: %s"),
+			    G_fatal_error(_("Unable to write to table: %s"),
 					  buf);
 		    }
-
 		}
-		else {
-		    if ((point->y[0] > Overlap.N)) {	/*(3) */
-			weight = (*point->y - Overlap.N) / overlap;
-			*point->z = (1 - weight) * interpolation;
+		else if ((*point->x >= Overlap.W) && (*point->x <= Overlap.E)){
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(3) */
+			weight = (General.N - *point->y) / overlap;
+			*point->z = weight * interpolation;
 
 			sprintf(buf, "%lf", *point->z);
 			db_append_string(&sql, buf);
@@ -173,11 +166,11 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 			db_append_string(&sql, buf);
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
-			    G_fatal_error(_("Unable to create table: %s"),
+			    G_fatal_error(_("Unable to write to table: %s"),
 					  buf);
 		    }
-		    else {	/*(1) */
-			weight = (Overlap.S - *point->y) / overlap;
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(1) */
+			weight = (*point->y - General.S) / overlap;
 			*point->z = (1 - weight) * interpolation;
 
 			sprintf(buf, "%lf", *point->z);
@@ -186,14 +179,17 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 			db_append_string(&sql, buf);
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
-			    G_fatal_error(_("Unable to create table: %s"),
+			    G_fatal_error(_("Unable to write to table: %s"),
 					  buf);
 		    }
 		}
 	    }
-	}
-    /*IF*/}
-    /*FOR*/ return;
+	}  /*IF*/
+    }  /*FOR*/
+
+    db_commit_transaction(driver);
+
+    return;
 }
 
 
@@ -255,69 +251,59 @@ double **P_Regular_Points(struct Cell_head *Elaboration, struct bound_box Genera
 
 		if (Vect_point_in_box(X, Y, interpolation, &Overlap)) {	/* (5) */
 		    matrix[row][col] = interpolation;
-
 		}
 		else {
-		    if ((X > Overlap.E)) {
-
-			if ((Y > Overlap.N)) {	/* (3) */
-			    csi = (X - Overlap.E) / overlap;
-			    eta = (Y - Overlap.N) / overlap;
-			    weight = (1 - csi) * (1 - eta);
+		    if ((X > Overlap.E) && (X < General.E)) {
+			if ((Y > Overlap.N) && (Y < General.N)) {	/* (3) */
+			    csi = (General.E - X) / overlap;
+			    eta = (General.N - Y) / overlap;
+			    weight = csi * eta;
 			    interpolation *= weight;
 			    matrix[row][col] += interpolation;
-
 			}
-			else if ((Y < Overlap.S)) {	/* (1) */
-			    csi = (X - Overlap.E) / overlap;
+			else if ((Y < Overlap.S) && (Y > General.S)) {	/* (1) */
+			    csi = (General.E - X) / overlap;
 			    eta = (Y - General.S) / overlap;
-			    weight = (1 - csi) * eta;
+			    weight = csi * eta;
 			    interpolation *= weight;
 			    matrix[row][col] = interpolation;
-
 			}
-			else {	/* (1) */
-			    weight = (X - Overlap.E) / overlap;
-			    interpolation *= 1 - weight;
+			else if ((Y >= Overlap.S) && (Y <= Overlap.N)) {	/* (1) */
+			    weight = (General.E - X ) / overlap;
+			    interpolation *= weight;
 			    matrix[row][col] = interpolation;
 			}
-
 		    }
-		    else if ((X < Overlap.W)) {
-
-			if ((Y > Overlap.N)) {	/* (4) */
+		    else if ((X < Overlap.W) && (X > General.W)) {
+			if ((Y > Overlap.N) && (Y < General.N)) {	/* (4) */
 			    csi = (X - General.W) / overlap;
-			    eta = (Y - Overlap.N) / overlap;
-			    weight = (1 - eta) * csi;
+			    eta = (General.N - Y) / overlap;
+			    weight = eta * csi;
 			    interpolation *= weight;
 			    matrix[row][col] += interpolation;
-
 			}
-			else if ((Y < Overlap.S)) {	/* (2) */
+			else if ((Y < Overlap.S) && (Y > General.S)) {	/* (2) */
 			    csi = (X - General.W) / overlap;
 			    eta = (Y - General.S) / overlap;
 			    weight = csi * eta;
 			    interpolation *= weight;
 			    matrix[row][col] += interpolation;
-
 			}
-			else {	/* (2) */
-			    weight = (Overlap.W - X) / overlap;
-			    interpolation *= 1 - weight;
+			else if ((Y >= Overlap.S) && (Y <= Overlap.N)) {	/* (2) */
+			    weight = (X - General.W) / overlap;
+			    interpolation *= weight;
 			    matrix[row][col] += interpolation;
 			}
-
 		    }
-		    else {
-			if ((Y > Overlap.N)) {	/* (3) */
-			    weight = (Y - Overlap.N) / overlap;
-			    interpolation *= 1 - weight;
+		    else if ((X >= Overlap.W) && (X <= Overlap.E)) {
+			if ((Y > Overlap.N) && (Y < General.N)) {	/* (3) */
+			    weight = (General.N - Y) / overlap;
+			    interpolation *= weight;
 			    matrix[row][col] += interpolation;
-
 			}
-			else {	/* (1) */
-			    weight = (Overlap.S - Y) / overlap;
-			    interpolation *= 1 - weight;
+			else if ((Y < Overlap.S) && (Y > General.S)) {	/* (1) */
+			    weight = (Y - General.S) / overlap;
+			    interpolation *= weight;
 			    matrix[row][col] = interpolation;
 			}
 		    }

vector/lidar/lidarlib/zones.c

@@ -23,11 +23,44 @@ void P_zero_dim(struct Reg_dimens *dim)
 }
 
 /*----------------------------------------------------------------------------------------*/
+/*
+ --------------------------------------------
+ |            Elaboration region            |
+ |   ------------------------------------   |
+ |   |          General region          |   |
+ |   |   ----------------------------   |   |
+ |   |   |                          |   |   |
+ |   |   |                          |   |   |
+ |   |   |                          |   |   |
+ |   |   |      Overlap region      |   |   |
+ |   |   |                          |   |   |
+ |   |   |                          |   |   |
+ |   |   |                          |   |   |
+ |   |   ----------------------------   |   |
+ |   |                                  |   |
+ |   ------------------------------------   |
+ |                                          |
+ --------------------------------------------
+
+ The terminology is misleading:
+ The Overlap region does NOT overlap with neighbouring segments,
+ but the Elaboration and General region do overlap
+
+ Elaboration is used for interpolation
+ Interpolated points in Elaboration but outside General are discarded
+ Interpolated points in General but outside Overlap are weighed by
+ their distance to Overlap and summed up
+ Interpolated points in Overlap are taken as they are
+
+ The buffer zones Elaboration - General and General - Overlap must be
+ large enough to avoid artefacts 
+ */
+
 int
 P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
 	      struct bound_box * Overlap, struct Reg_dimens dim, int type)
 {
-    /* Set the Elaboration region limits-> Also set the limits of the orlo and overlapping regions->
+    /* Set the Elaboration, General, and Overlap region limits
      * Returns 0 on success; -1 on failure*/
     struct Cell_head orig;
 
@@ -54,7 +87,7 @@ P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
 	Overlap->E = General->E - dim.overlap;
 	return 0;
 
-    case FIRST_ROW:		/* It is just started with first row */
+    case FIRST_ROW:		/* Just started with first row */
 	Elaboration->north = orig.north + 2 * dim.orlo_h;
 	Elaboration->south = Elaboration->north - dim.latoN;
 	General->N = Elaboration->north - 2 * dim.orlo_h;
@@ -63,13 +96,13 @@ P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
 	Overlap->S = General->S + dim.overlap;
 	return 0;
 
-    case LAST_ROW:		/* It is reached last row */
+    case LAST_ROW:		/* Reached last row */
 	Elaboration->south = orig.south - 2 * dim.orlo_h;
 	General->S = Elaboration->south + 2 * dim.orlo_h;
 	Overlap->S = General->S;
 	return 0;
 
-    case FIRST_COLUMN:		/* It is just started with first column */
+    case FIRST_COLUMN:		/* Just started with first column */
 	Elaboration->west = orig.west - 2 * dim.orlo_v;
 	Elaboration->east = Elaboration->west + dim.latoE;
 	General->W = Elaboration->west + 2 * dim.orlo_v;
@@ -78,7 +111,7 @@ P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
 	Overlap->E = General->E - dim.overlap;
 	return 0;
 
-    case LAST_COLUMN:		/* It is reached last column */
+    case LAST_COLUMN:		/* Reached last column */
 	Elaboration->east = orig.east + 2 * dim.orlo_v;
 	General->E = Elaboration->east - 2 * dim.orlo_v;
 	Overlap->E = General->E;
@@ -171,20 +204,22 @@ int P_set_dim(struct Reg_dimens *dim, double pe, double pn, int *nsplx, int *nsp
 /*----------------------------------------------------------------------------------------*/
 int P_get_orlo(int interpolator, struct Reg_dimens *dim, double pe, double pn)
 {
-    /* Set the orlo regions dimension->
-     * Returns 1 on success of bilinear; 2 on success of bicubic, 0 on failure-> */
-    if (interpolator == 1) {	/* the interpolator's bilinear */
-	dim->orlo_v = 9 * pe;	/*4 */
+    /* Set the orlo regions dimension
+     * Returns 1 on success of bilinear; 2 on success of bicubic, 0 on failure */
+    if (interpolator == P_BILINEAR) {
+       	/* in case of edge artefacts, increase as multiples of 3 */
+	dim->orlo_v = 9 * pe;
 	dim->orlo_h = 9 * pn;
 	return 1;
     }
-    else if (interpolator == 0) {	/* the interpolator's bicubic */
+    else if (interpolator == P_BICUBIC) {
+       	/* in case of edge artefacts, increase as multiples of 4 */
 	dim->orlo_v = 12 * pe;	/*3 */
 	dim->orlo_h = 12 * pn;
 	return 2;
     }
     else
-	return 0;		/* The interpolator it's not bilinear nor bicubic!! */
+	return 0;		/* The interpolator is neither bilinear nor bicubic!! */
 }
 
 /*----------------------------------------------------------------------------------------*/
@@ -239,7 +274,7 @@ double P_estimate_splinestep(struct Map_info *Map, double *dens, double *dist)
     struct bound_box region_box;
     struct Cell_head orig;
 
-    G_get_window(&orig);
+    G_get_set_window(&orig);
     Vect_region_box(&orig, &region_box);
 
     points = Vect_new_line_struct();
@@ -257,7 +292,7 @@ double P_estimate_splinestep(struct Map_info *Map, double *dens, double *dist)
 	else
 	    z = 0.0;
 
-	/* Reading and storing points only if it is in elaboration_reg */
+	/* only use points in current region */
 	if (Vect_point_in_box(x, y, z, &region_box)) {
 	    npoints++;
 
@@ -278,7 +313,9 @@ double P_estimate_splinestep(struct Map_info *Map, double *dens, double *dist)
 	}
     }
     if (npoints > 0) {
+	/* estimated average distance between points in map units */
 	*dist = sqrt(((xmax - xmin) * (ymax - ymin)) / npoints);
+	/* estimated point density as number of points per square map unit */
 	*dens = npoints / ((xmax - xmin) * (ymax - ymin));
 	return 0;
     }
@@ -305,7 +342,7 @@ struct Point *P_Read_Vector_Region_Map(struct Map_info *Map,
     points = Vect_new_line_struct();
     categories = Vect_new_cats_struct();
 
-    /* Reading the elaboration zone points */
+    /* Reading points inside elaboration zone */
     Vect_region_box(Elaboration, &elaboration_box);
 
     npoints = 0;
@@ -326,7 +363,7 @@ struct Point *P_Read_Vector_Region_Map(struct Map_info *Map,
 	else
 	    z = 0.0;
 
-	/* Reading and storing points only if it is in elaboration_reg */
+	/* Reading and storing points only if in elaboration_reg */
 	if (Vect_point_in_box(x, y, z, &elaboration_box)) {
 	    npoints++;
 	    if (npoints >= pippo) {
@@ -355,42 +392,69 @@ struct Point *P_Read_Vector_Region_Map(struct Map_info *Map,
 }
 
 /*------------------------------------------------------------------------------------------------*/
-int P_Create_Aux_Table(dbDriver * driver, char *tab_name)
+int P_Create_Aux2_Table(dbDriver * driver, char *tab_name)
 {
     dbTable *auxiliar_tab;
     dbColumn *column;
-    int created = FALSE;
 
-    auxiliar_tab = db_alloc_table(4);
+    auxiliar_tab = db_alloc_table(2);
     db_set_table_name(auxiliar_tab, tab_name);
     db_set_table_description(auxiliar_tab,
 			     "Intermediate interpolated values");
 
-    column = db_get_table_column(auxiliar_tab, 2);
-    db_set_column_name(column, "Y");
-    db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+    column = db_get_table_column(auxiliar_tab, 0);
+    db_set_column_name(column, "ID");
+    db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
 
     column = db_get_table_column(auxiliar_tab, 1);
-    db_set_column_name(column, "X");
-    db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+    db_set_column_name(column, "Interp");
+    db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
+
+    if (db_create_table(driver, auxiliar_tab) == DB_OK) {
+	G_debug(1, _("<%s> created in database."), tab_name);
+	return TRUE;
+    }
+    else
+	G_warning(_("<%s> has not been created in database."), tab_name);
+
+    return FALSE;
+}
+
+/*------------------------------------------------------------------------------------------------*/
+int P_Create_Aux4_Table(dbDriver * driver, char *tab_name)
+{
+    dbTable *auxiliar_tab;
+    dbColumn *column;
+
+    auxiliar_tab = db_alloc_table(4);
+    db_set_table_name(auxiliar_tab, tab_name);
+    db_set_table_description(auxiliar_tab,
+			     "Intermediate interpolated values");
 
     column = db_get_table_column(auxiliar_tab, 0);
     db_set_column_name(column, "ID");
     db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
 
-    column = db_get_table_column(auxiliar_tab, 3);
+    column = db_get_table_column(auxiliar_tab, 1);
     db_set_column_name(column, "Interp");
     db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
 
+    column = db_get_table_column(auxiliar_tab, 2);
+    db_set_column_name(column, "X");
+    db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+
+    column = db_get_table_column(auxiliar_tab, 3);
+    db_set_column_name(column, "Y");
+    db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+
     if (db_create_table(driver, auxiliar_tab) == DB_OK) {
 	G_debug(1, _("<%s> created in database."), tab_name);
-	created = TRUE;
-	return created;
+	return TRUE;
     }
     else
-	G_fatal_error(_("<%s> has not been created in database."), tab_name);
+	G_warning(_("<%s> has not been created in database."), tab_name);
 
-    return created;
+    return FALSE;
 }
 
 /*------------------------------------------------------------------------------------------------*/
@@ -481,7 +545,7 @@ P_Aux_to_Vector(struct Map_info *Map, struct Map_info *Out, dbDriver * driver,
 	    value = db_get_column_value(column);
 	else
 	    continue;
-	coordX = db_get_value_double(value);
+	coordZ = db_get_value_double(value);
 
 	column = db_get_table_column(table, 2);
 	type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
@@ -489,7 +553,7 @@ P_Aux_to_Vector(struct Map_info *Map, struct Map_info *Out, dbDriver * driver,
 	    value = db_get_column_value(column);
 	else
 	    continue;
-	coordY = db_get_value_double(value);
+	coordX = db_get_value_double(value);
 
 	column = db_get_table_column(table, 3);
 	type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
@@ -497,7 +561,7 @@ P_Aux_to_Vector(struct Map_info *Map, struct Map_info *Out, dbDriver * driver,
 	    value = db_get_column_value(column);
 	else
 	    continue;
-	coordZ = db_get_value_double(value);
+	coordY = db_get_value_double(value);
 
 	Vect_copy_xyz_to_pnts(point, &coordX, &coordY, &coordZ, 1);
 	Vect_reset_cats(cat);
@@ -528,10 +592,30 @@ double **P_Null_Matrix(double **matrix)
 }
 #endif
 
-/*------------------------------------------------------------------------------------------------*/
-#ifdef notdef
-/* Subzones definition */
------------------------|4 | 3 | 3 | ||-----------------------|||||||2 | 5 | 1
-    | ||||||||-----------------------|2 | 1 | 1 |
-    ||-----------------------||||||||||||||||||-----------------------||||||-----------------------
-#endif
+/*! DEFINITION OF THE SUBZONES
+
+  5: inside Overlap region
+  all others: inside General region but outside Overlap region
+
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       |4|   3   |3|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       |2|   5   |1|       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       |2|   1   |1|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
+ */

vector/lidar/v.lidar.correction/correction.c

@@ -30,54 +30,65 @@
 void
 P_Sparse_Correction(struct Map_info *In, struct Map_info *Out,
 		    struct Map_info *Terrain, struct Cell_head *Elaboration,
-		    struct bound_box General, struct bound_box Overlap, double **obs,
-		    double *param, int *line_num, double passoN,
-		    double passoE, double overlap, double HighThresh,
-		    double LowThresh, int nsplx, int nsply, int num_points,
+		    struct bound_box General, struct bound_box Overlap,
+		    double **obs, struct lidar_cat *lcat, double *param,
+		    int *line_num, double passoN, double passoE,
+		    double overlap, double HighThresh, double LowThresh,
+		    int nsplx, int nsply, int num_points,
 		    dbDriver * driver, double mean, char *tab_name)
 {
     int i = 0, class;
     double interpolation, csi, eta, weight;
 
-    struct line_pnts *points;
+    struct bound_box elaboration_box;
+    struct line_pnts *point;
     struct line_cats *cats;
 
-    points = Vect_new_line_struct();
+    point = Vect_new_line_struct();
     cats = Vect_new_cats_struct();
 
-    Vect_rewind(In);
-    while (Vect_read_next_line(In, points, cats) > 0) {
-	i++;
-	if (Vect_point_in_box(*points->x, *points->y, *points->z, &General)) {
+    Vect_region_box(Elaboration, &elaboration_box);
+
+    db_begin_transaction(driver);
+    
+    for (i = 0; i < num_points; i++) {	/* Sparse points */
+	G_percent(i, num_points, 2);
+	Vect_reset_line(point);
+	Vect_reset_cats(cats);
+	
+	if (Vect_point_in_box(obs[i][0], obs[i][1], mean, &General)) {
 	    interpolation =
-		dataInterpolateBilin(*points->x, *points->y, passoE, passoN,
+		dataInterpolateBilin(obs[i][0], obs[i][1], passoE, passoN,
 				     nsplx, nsply, Elaboration->west,
 				     Elaboration->south, param);
 	    interpolation += mean;
 
-	    if (Vect_point_in_box(*points->x, *points->y, *points->z, &Overlap)) {	/*(5) */
+	    Vect_copy_xyz_to_pnts(point, &obs[i][0], &obs[i][1], &obs[i][2],
+				  1);
+	    *point->z += mean;
+
+	    if (Vect_point_in_box(*point->x, *point->y, *point->z, &Overlap)) {	/*(5) */
 
-		Vect_cat_get(cats, F_CLASSIFICATION, &class);
+		Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+		Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+		Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+		class = lcat[i].cat_class;
 		class =
-		    correction(class, *points->z, interpolation, HighThresh,
+		    correction(class, *point->z, interpolation, HighThresh,
 			       LowThresh);
-		Vect_cat_del(cats, F_CLASSIFICATION);
 		Vect_cat_set(cats, F_CLASSIFICATION, class);
 
 		if (class == TERRAIN_SINGLE || class == TERRAIN_DOUBLE)
-		    Vect_write_line(Terrain, GV_POINT, points, cats);
-
-		Vect_write_line(Out, GV_POINT, points, cats);
+		    Vect_write_line(Terrain, GV_POINT, point, cats);
 
+		Vect_write_line(Out, GV_POINT, point, cats);
 	    }
 	    else {
-
-		if ((*points->x > Overlap.E)) {
-
-		    if ((*points->y > Overlap.N)) {	/*(3) */
-			csi = (*points->x - Overlap.E) / overlap;
-			eta = (*points->y - Overlap.N) / overlap;
-			weight = (1 - csi) * (1 - eta);
+		if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(3) */
+			csi = (General.E - *point->x) / overlap;
+			eta = (General.N - *point->y) / overlap;
+			weight = csi * eta;
 			interpolation *= weight;
 
 			if (Select_Correction
@@ -87,21 +98,19 @@ P_Sparse_Correction(struct Map_info *In, struct Map_info *Out,
 			if (UpDate_Correction
 			    (interpolation, line_num[i], driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to update the database"));
-
 		    }
-		    else if ((*points->y < Overlap.S)) {	/*(1) */
-			csi = (*points->x - Overlap.E) / overlap;
-			eta = (*points->y - General.S) / overlap;
-			weight = (1 - csi) * eta;
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(1) */
+			csi = (General.E - *point->x) / overlap;
+			eta = (*point->y - General.S) / overlap;
+			weight = csi * eta;
 
 			if (Insert_Correction
 			    (interpolation * weight, line_num[i],
 			     driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to write in the database"));
-
 		    }
-		    else {	/*(1) */
-			weight = (*points->x - Overlap.E) / overlap;
+		    else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) {	/*(1) */
+			weight = (General.E - *point->x) / overlap;
 
 			if (Insert_Correction
 			    (interpolation * weight, line_num[i],
@@ -109,32 +118,34 @@ P_Sparse_Correction(struct Map_info *In, struct Map_info *Out,
 			    G_fatal_error(_("Impossible to write in the database"));
 		    }
 		}
-		else if ((*points->x < Overlap.W)) {
-		    if ((*points->y > Overlap.N)) {	/*(4) */
-			csi = (*points->x - General.W) / overlap;
-			eta = (*points->y - Overlap.N) / overlap;
-			weight = (1 - eta) * csi;
-
+		else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(4) */
+			csi = (*point->x - General.W) / overlap;
+			eta = (General.N - *point->y) / overlap;
+			weight = eta * csi;
 			interpolation *= weight;
+
 			if (Select_Correction
 			    (&interpolation, line_num[i], driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to read the database"));
 
-			Vect_cat_get(cats, F_CLASSIFICATION, &class);
+			Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+			Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+			Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+			class = lcat[i].cat_class;
 			class =
-			    correction(class, *points->z, interpolation,
+			    correction(class, *point->z, interpolation,
 				       HighThresh, LowThresh);
 			Vect_cat_set(cats, F_CLASSIFICATION, class);
 
-			Vect_write_line(Out, GV_POINT, points, cats);
+			Vect_write_line(Out, GV_POINT, point, cats);
 			if (class == TERRAIN_SINGLE ||
 			    class == TERRAIN_DOUBLE)
-			    Vect_write_line(Terrain, GV_POINT, points, cats);
-
+			    Vect_write_line(Terrain, GV_POINT, point, cats);
 		    }
-		    else if ((*points->y < Overlap.S)) {	/*(2) */
-			csi = (*points->x - General.W) / overlap;
-			eta = (*points->y - General.S) / overlap;
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(2) */
+			csi = (*point->x - General.W) / overlap;
+			eta = (*point->y - General.S) / overlap;
 			weight = csi * eta;
 			interpolation *= weight;
 
@@ -146,51 +157,55 @@ P_Sparse_Correction(struct Map_info *In, struct Map_info *Out,
 			    (interpolation, line_num[i], driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to update the database"));
 		    }
-		    else {	/*(2) */
-			weight = (Overlap.W - *points->x) / overlap;
+		    else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) {	/*(2) */
+			weight = (*point->x - General.W) / overlap;
 			interpolation *= weight;
 
 			if (Select_Correction
 			    (&interpolation, line_num[i], driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to read the database"));
 
-			Vect_cat_get(cats, F_CLASSIFICATION, &class);
+			Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+			Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+			Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+			class = lcat[i].cat_class;
 			class =
-			    correction(class, *points->z, interpolation,
+			    correction(class, *point->z, interpolation,
 				       HighThresh, LowThresh);
 			Vect_cat_set(cats, F_CLASSIFICATION, class);
 
-			Vect_write_line(Out, GV_POINT, points, cats);
+			Vect_write_line(Out, GV_POINT, point, cats);
 			if (class == TERRAIN_SINGLE ||
 			    class == TERRAIN_DOUBLE)
-			    Vect_write_line(Terrain, GV_POINT, points, cats);
-
+			    Vect_write_line(Terrain, GV_POINT, point, cats);
 		    }
-
 		}
-		else {
-		    if ((*points->y > Overlap.N)) {	/*(3) */
-			weight = (*points->y - Overlap.N) / overlap;
+		else if ((*point->x >= Overlap.W) && (*point->x <= Overlap.E)) {
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(3) */
+			weight = (General.N - *point->y) / overlap;
 			interpolation *= weight;
 
 			if (Select_Correction
 			    (&interpolation, line_num[i], driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to read the database"));
 
-			Vect_cat_get(cats, F_CLASSIFICATION, &class);
+			Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+			Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+			Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+			class = lcat[i].cat_class;
 			class =
-			    correction(class, *points->z, interpolation,
+			    correction(class, *point->z, interpolation,
 				       HighThresh, LowThresh);
 			Vect_cat_set(cats, F_CLASSIFICATION, class);
 
-			Vect_write_line(Out, GV_POINT, points, cats);
+			Vect_write_line(Out, GV_POINT, point, cats);
 			if (class == TERRAIN_SINGLE ||
 			    class == TERRAIN_DOUBLE)
-			    Vect_write_line(Terrain, GV_POINT, points, cats);
+			    Vect_write_line(Terrain, GV_POINT, point, cats);
 
 		    }
-		    else {	/*(1) */
-			weight = (Overlap.S - *points->y) / overlap;
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(1) */
+			weight = (*point->y - General.S) / overlap;
 
 			if (Insert_Correction
 			    (interpolation * weight, line_num[i],
@@ -199,14 +214,14 @@ P_Sparse_Correction(struct Map_info *In, struct Map_info *Out,
 		    }
 		}
 	    }
-	}
-	/*IF*/ Vect_reset_line(points);
-	Vect_reset_cats(cats);
-
-    }
-    /*FOR*/ Vect_destroy_line_struct(points);
+	}  /* if in General box */
+    }  /* while */
+    G_percent(num_points, num_points, 2);
+    Vect_destroy_line_struct(point);
     Vect_destroy_cats_struct(cats);
 
+    db_commit_transaction(driver);
+
     return;
 }
 
@@ -240,9 +255,8 @@ int Select_Correction(double *Interp, int line_num, dbDriver * driver, char *tab
     dbValue *Interp_value;
 
     db_init_string(&sql);
-    db_zero_string(&sql);
 
-    sprintf(buf, "SELECT Interp FROM %s WHERE ID=%d", tab_name, line_num);
+    sprintf(buf, "SELECT ID, Interp FROM %s WHERE ID=%d", tab_name, line_num);
     db_append_string(&sql, buf);
 
     if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
@@ -260,6 +274,7 @@ int Select_Correction(double *Interp, int line_num, dbDriver * driver, char *tab
 
 	*Interp += db_get_value_double(Interp_value);
     }
+    db_close_cursor(&cursor);
     db_free_string(&sql);
     return DB_OK;
 }
@@ -300,21 +315,20 @@ int UpDate_Correction(double Interp, int line_num, dbDriver * driver, char *tab_
 struct Point *P_Read_Vector_Correction(struct Map_info *Map,
 				       struct Cell_head *Elaboration,
 				       int *num_points, int *num_terrain,
-				       int dim_vect)
+				       int dim_vect, struct lidar_cat **lcat)
 {
-    int line_num, npoints, nterrain, pippo, cat_edge;
+    int line_num, npoints, nterrain, pippo, cat;
     double x, y, z;
     struct Point *obs;
+    struct lidar_cat *lc;
     struct line_pnts *points;
     struct line_cats *categories;
     struct bound_box elaboration_box;
 
     pippo = dim_vect;
 
-    /*if (first_it) */
     obs = (struct Point *)G_calloc(pippo, sizeof(struct Point));
-    /*else */
-    /*obs = (struct Point*) G_realloc ((void *) obs, pippo * sizeof(struct Point)); */
+    lc = (struct lidar_cat *)G_calloc(pippo, sizeof(struct lidar_cat));
 
     points = Vect_new_line_struct();
     categories = Vect_new_cats_struct();
@@ -344,12 +358,10 @@ struct Point *P_Read_Vector_Correction(struct Map_info *Map,
 	    npoints++;
 	    if (npoints >= pippo) {
 		pippo += dim_vect;
-		obs =
-		    (struct Point *)G_realloc((void *)obs,
-					      (signed int)(pippo *
-							   sizeof(struct
-								  Point)));
-		/*lineID = (int *) G_realloc ((void *) lineID, pippo * sizeof(int)); */
+		obs = (struct Point *)G_realloc((void *)obs,
+			    (signed int)(pippo * sizeof(struct Point)));
+		lc = (struct lidar_cat *)G_realloc((void *)lc,
+			    (signed int)(pippo * sizeof(struct lidar_cat)));
 	    }
 
 	    /* Storing observation vector */
@@ -357,27 +369,55 @@ struct Point *P_Read_Vector_Correction(struct Map_info *Map,
 	    obs[npoints - 1].coordY = y;
 	    obs[npoints - 1].coordZ = z;
 	    obs[npoints - 1].lineID = line_num;	/* Storing also the line's number */
-	    Vect_cat_get(categories, F_EDGE_DETECTION_CLASS, &cat_edge);
-	    obs[npoints - 1].cat = cat_edge;
+	    Vect_cat_get(categories, F_EDGE_DETECTION_CLASS, &cat);
+	    obs[npoints - 1].cat = cat;
+	    lc[npoints - 1].cat_edge = cat;
+
+	    /* Only terrain points */
+	    if (cat == TERRAIN_SINGLE)
+		nterrain++;
+
+	    Vect_cat_get(categories, F_CLASSIFICATION, &cat);
+	    lc[npoints - 1].cat_class = cat;
+	    Vect_cat_get(categories, F_INTERPOLATION, &cat);
+	    lc[npoints - 1].cat_interp = cat;
+	    Vect_cat_get(categories, F_COUNTER_OBJ, &cat);
+	    lc[npoints - 1].cat_obj = cat;
 	}
-
-	/* Only terrain points */
-	if (cat_edge == TERRAIN_SINGLE)
-	    nterrain++;
     }
     Vect_destroy_line_struct(points);
     Vect_destroy_cats_struct(categories);
 
     *num_points = npoints;
     *num_terrain = nterrain;
+    *lcat = lc;
     return obs;
 }
 
-
-
-#ifdef notdef
-/* SUBZONES DEFINITION */
------------------------|4 | 3 | 3 | ||-----------------------|||||||2 | 5 | 1
-    | ||||||||-----------------------|2 | 1 | 1 |
-    ||-----------------------||||||||||||||||||-----------------------||||||-----------------------
-#endif
+/*! DEFINITION OF THE SUBZONES 
+
+  5: inside Overlap region
+  all others: inside General region but outside Overlap region
+
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       |4|   3   |3|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       |2|   5   |1|       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       |2|   1   |1|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
+ */

vector/lidar/v.lidar.correction/correction.h

@@ -7,6 +7,14 @@
 
 #include <grass/PolimiFunct.h>
 
+struct lidar_cat
+{
+    int cat_edge;     /* category in layer F_EDGE_DETECTION_CLASS */
+    int cat_class;    /* category in layer F_CLASSIFICATION */
+    int cat_interp;   /* category in layer F_INTERPOLATION */
+    int cat_obj;      /* category in layer F_COUNTER_OBJ */
+};
+
 void
 P_Sparse_Correction(struct Map_info *, /**/
 		    struct Map_info *, /**/
@@ -15,6 +23,7 @@ P_Sparse_Correction(struct Map_info *, /**/
 		    struct bound_box, /**/
 		    struct bound_box, /**/
 		    double **, /**/
+		    struct lidar_cat *,
 		    double *, /**/
 		    int *, /**/
 		    double, /**/
@@ -37,6 +46,7 @@ int Select_Correction(double *, /**/ int, /**/ dbDriver *, /**/ char *);
 
 struct Point *P_Read_Vector_Correction(struct Map_info *, /**/
 				       struct Cell_head *, /**/
-				       int *, /**/ int *, /**/ int /**/);
+				       int *, /**/ int *, /**/ int /**/,
+				       struct lidar_cat **);
 
 int correction(int, double, double, double, double);

vector/lidar/v.lidar.correction/main.c

@@ -51,7 +51,7 @@ int main(int argc, char *argv[])
 
     int *lineVect;
     double *TN, *Q, *parVect;	/* Interpolating and least-square vectors */
-    double **N, **obsVect;	/* Interpolation and least-square matrix */
+    double **N, **obsVect, **obsVect_all;	/* Interpolation and least-square matrix */
 
     struct Map_info In, Out, Terrain;
     struct Option *in_opt, *out_opt, *out_terrain_opt, *passoE_opt,
@@ -64,6 +64,7 @@ int main(int argc, char *argv[])
     struct bound_box general_box, overlap_box;
 
     struct Point *observ;
+    struct lidar_cat *lcat;
 
     dbDriver *driver;
 
@@ -77,9 +78,9 @@ int main(int argc, char *argv[])
 
     spline_step_flag = G_define_flag();
     spline_step_flag->key = 'e';
-    spline_step_flag->label = _("Estimate spline step value");
+    spline_step_flag->label = _("Estimate point density and distance");
     spline_step_flag->description =
-	_("Estimate a good spline step value for the input vector points within the current region extends and quit");
+	_("Estimate point density and distance for the input vector points within the current region extends and quit");
 
     in_opt = G_define_standard_option(G_OPT_V_INPUT);
     in_opt->description =
@@ -182,10 +183,14 @@ int main(int argc, char *argv[])
     if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL)
 	G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
 
-    Vect_set_open_level(1);	/*without topology */
+    Vect_set_open_level(1);	/* without topology */
     if (1 > Vect_open_old(&In, in_opt->answer, mapset))
 	G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
 
+    /* Input vector must be 3D */
+    if (!Vect_is_3d(&In))
+	G_fatal_error(_("Input vector map <%s> is not 3D!"), in_opt->answer);
+
     /* Estimate point density and mean distance for current region */
     if (spline_step_flag->answer) {
 	double dens, dist;
@@ -228,13 +233,18 @@ int main(int argc, char *argv[])
 
     /* Create auxiliar table */
     if ((flag_auxiliar =
-	 P_Create_Aux_Table(driver, table_name)) == FALSE) {
+	 P_Create_Aux2_Table(driver, table_name)) == FALSE) {
 	Vect_close(&In);
 	Vect_close(&Out);
 	Vect_close(&Terrain);
 	exit(EXIT_FAILURE);
     }
 
+    db_create_index2(driver, table_name, "ID");
+    /* sqlite likes that */
+    db_close_database_shutdown_driver(driver);
+    driver = db_start_driver_open_database(dvr, db);
+
     /* Setting regions and boxes */
     G_get_set_window(&original_reg);
     G_get_set_window(&elaboration_reg);
@@ -247,27 +257,36 @@ int main(int argc, char *argv[])
     ew_resol = original_reg.ew_res;
     ns_resol = original_reg.ns_res;
 
-    /* Fixing parameters of the elaboration region */
-    /*! Each original_region will be divided into several subregions. These
-     *  subregions will overlap with their neibourghing subregions. This overlapping
-     *  is calculated as OVERLAP_SIZE times the east-west spline step. */
+    /*------------------------------------------------------------------
+      | Subdividing and working with tiles: 									
+      | Each original region will be divided into several subregions. 
+      | Each one will be overlaped by its neighbouring subregions. 
+      | The overlapping is calculated as a fixed OVERLAP_SIZE times
+      | the largest spline step plus 2 * orlo
+      ----------------------------------------------------------------*/
 
+    /* Fixing parameters of the elaboration region */
     P_zero_dim(&dims);
-    N_extension = original_reg.north - original_reg.south;
-    E_extension = original_reg.east - original_reg.west;
 
     nsplx_adj = NSPLX_MAX;
     nsply_adj = NSPLY_MAX;
-    dims.overlap = OVERLAP_SIZE * passoE;
+    if (passoN > passoE)
+	dims.overlap = OVERLAP_SIZE * passoN;
+    else
+	dims.overlap = OVERLAP_SIZE * passoE;
     P_get_orlo(P_BILINEAR, &dims, passoE, passoN);
     P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
 
     G_verbose_message("adjusted EW splines %d", nsplx_adj);
     G_verbose_message("adjusted NS splines %d", nsply_adj);
 
+    /* calculate number of subzones */
     orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
     orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
 
+    N_extension = original_reg.north - original_reg.south;
+    E_extension = original_reg.east - original_reg.west;
+
     nsubregion_col = ceil(E_extension / orloE) + 0.5;
     nsubregion_row = ceil(N_extension / orloN) + 0.5;
 
@@ -278,9 +297,9 @@ int main(int argc, char *argv[])
 
     nsubzones = nsubregion_row * nsubregion_col;
 
-    /* Subdividing and working with tiles */
     elaboration_reg.south = original_reg.north;
     last_row = FALSE;
+
     while (last_row == FALSE) {	/* For each row */
 
 	P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
@@ -344,23 +363,24 @@ int main(int argc, char *argv[])
 	    G_debug(1, _("read vector region map"));
 	    observ =
 		P_Read_Vector_Correction(&In, &elaboration_reg, &npoints,
-					 &nterrain, dim_vect);
+					 &nterrain, dim_vect, &lcat);
 
 	    G_debug(5, _("npoints = %d, nterrain = %d"), npoints, nterrain);
 	    if (npoints > 0) {	/* If there is any point falling into elaboration_reg. */
 		count_terrain = 0;
 		nparameters = nsplx * nsply;
 
-		/* Mean's calculation */
-		G_debug(3, _("Mean's calculation"));
+		/* Mean calculation */
+		G_debug(3, _("Mean calculation"));
 		mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
 
 		/*Least Squares system */
 		BW = P_get_BandWidth(P_BILINEAR, nsply);	/* Bilinear interpolation */
 		N = G_alloc_matrix(nparameters, BW);	/* Normal matrix */
 		TN = G_alloc_vector(nparameters);	/* vector */
-		parVect = G_alloc_vector(nparameters);	/* Bicubic parameters vector */
-		obsVect = G_alloc_matrix(nterrain + 1, 3);	/* Observation vector */
+		parVect = G_alloc_vector(nparameters);	/* Bilinear parameters vector */
+		obsVect = G_alloc_matrix(nterrain + 1, 3);	/* Observation vector with terrain points */
+		obsVect_all = G_alloc_matrix(npoints + 1, 3);	/* Observation vector with all points */
 		Q = G_alloc_vector(nterrain + 1);	/* "a priori" var-cov matrix */
 		lineVect = G_alloc_ivector(npoints + 1);
 
@@ -375,9 +395,14 @@ int main(int argc, char *argv[])
 			count_terrain++;
 		    }
 		    lineVect[i] = observ[i].lineID;
+		    obsVect_all[i][0] = observ[i].coordX;
+		    obsVect_all[i][1] = observ[i].coordY;
+		    obsVect_all[i][2] = observ[i].coordZ - mean;
 		}
 
-		G_debug(3, _("M.Q. solution"));
+		G_free(observ);
+
+		G_verbose_message(_("Bilinear interpolation"));
 		normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
 			       nsply, elaboration_reg.west,
 			       elaboration_reg.south, nterrain, nparameters,
@@ -388,19 +413,25 @@ int main(int argc, char *argv[])
 		G_free_matrix(N);
 		G_free_vector(TN);
 		G_free_vector(Q);
+		G_free_matrix(obsVect);
 
-		G_debug(3, _("Correction and creation of terrain vector"));
+		G_verbose_message( _("Correction and creation of terrain vector"));
 		P_Sparse_Correction(&In, &Out, &Terrain, &elaboration_reg,
-				    general_box, overlap_box, obsVect,
+				    general_box, overlap_box, obsVect_all, lcat,
 				    parVect, lineVect, passoN, passoE,
 				    dims.overlap, HighThresh, LowThresh,
 				    nsplx, nsply, npoints, driver, mean, table_name);
 
 		G_free_vector(parVect);
-		G_free_matrix(obsVect);
+		G_free_matrix(obsVect_all);
 		G_free_ivector(lineVect);
 	    }
-	    G_free(observ);
+	    else {
+		G_free(observ);
+		G_warning(_("No data within this subzone. "
+			    "Consider changing the spline step."));
+	    }
+	    G_free(lcat);
 	}			/*! END WHILE; last_column = TRUE */
     }				/*! END WHILE; last_row = TRUE */
 

vector/lidar/v.lidar.edgedetection/edgedetection.c

@@ -180,6 +180,8 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
     point = Vect_new_line_struct();
     categories = Vect_new_cats_struct();
 
+    db_begin_transaction(driver);
+    
     for (i = 0; i < num_points; i++) {	/* Sparse points */
 	G_percent(i, num_points, 2);
 
@@ -201,7 +203,6 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 	    /*Vect_cat_set (categories, F_INTERPOLATION, line_out_counter); */
 
 	    if (Vect_point_in_box(obs[i][0], obs[i][1], interpolation, &Overlap)) {	/*(5) */
-
 		residual = *point->z - interpolation;
 		edge =
 		    edge_detection(Elaboration, Overlap, parBilin, *point->x,
@@ -214,15 +215,13 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 		Insert_Interpolation(interpolation, line_out_counter, driver,
 				     tabint_name);
 		line_out_counter++;
-
 	    }
 	    else {
-		if ((*point->x > Overlap.E) && (*point->x != General.E)) {
-
-		    if ((*point->y > Overlap.N) && (*point->y != General.N)) {	/*(3) */
-			csi = (*point->x - Overlap.E) / overlap;
-			eta = (*point->y - Overlap.N) / overlap;
-			weight = (1 - csi) * (1 - eta);
+		if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(3) */
+			csi = (General.E - *point->x) / overlap;
+			eta = (General.N - *point->y) / overlap;
+			weight = csi * eta;
 
 			gradient[0] *= weight;
 			gradient[1] *= weight;
@@ -237,12 +236,11 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 			           interpolation, line_num[i], driver,
 			           tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to update aux table"));
-
 		    }
-		    else if ((*point->y < Overlap.S) && (*point->y != General.S)) {	/*(1) */
-			csi = (*point->x - Overlap.E) / overlap;
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(1) */
+			csi = (General.E - *point->x) / overlap;
 			eta = (*point->y - General.S) / overlap;
-			weight = (1 - csi) * eta;
+			weight = csi * eta;
 
 			gradient[0] *= weight;
 			gradient[1] *= weight;
@@ -254,7 +252,7 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 
 		    }
 		    else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) {	/*(1) */
-			weight = (*point->x - Overlap.E) / overlap;
+			weight = (General.E - *point->x) / overlap;
 
 			gradient[0] *= weight;
 			gradient[1] *= weight;
@@ -264,14 +262,12 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 			           line_num[i], driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to write to aux table"));
 		    }
-
 		}
-		else if ((*point->x < Overlap.W) && (*point->x != General.W)) {
-
-		    if ((*point->y > Overlap.N) && (*point->y != General.N)) {	/*(4) */
+		else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(4) */
 			csi = (*point->x - General.W) / overlap;
-			eta = (*point->y - Overlap.N) / overlap;
-			weight = (1 - eta) * csi;
+			eta = (General.N - *point->y) / overlap;
+			weight = eta * csi;
 
 			gradient[0] *= weight;
 			gradient[1] *= weight;
@@ -297,9 +293,8 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 			Insert_Interpolation(interpolation, line_out_counter,
 					     driver, tabint_name);
 			line_out_counter++;
-
 		    }
-		    else if ((*point->y < Overlap.S) && (*point->y != General.S)) {	/*(2) */
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(2) */
 			csi = (*point->x - General.W) / overlap;
 			eta = (*point->y - General.S) / overlap;
 			weight = csi * eta;
@@ -320,7 +315,7 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 
 		    }
 		    else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) {	/*(2) */
-			weight = (Overlap.W - *point->x) / overlap;
+			weight = (*point->x - General.W) / overlap;
 
 			gradient[0] *= weight;
 			gradient[1] *= weight;
@@ -347,11 +342,10 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 					     driver, tabint_name);
 			line_out_counter++;
 		    }
-
 		}
 		else if ((*point->x <= Overlap.E) && (*point->x >= Overlap.W)) {
-		    if ((*point->y > Overlap.N) && (*point->y != General.N)) {	/*(3) */
-			weight = (*point->y - Overlap.N) / overlap;
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(3) */
+			weight = (General.N - *point->y) / overlap;
 
 			gradient[0] *= weight;
 			gradient[1] *= weight;
@@ -378,8 +372,8 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 					     driver, tabint_name);
 			line_out_counter++;
 		    }
-		    else if ((*point->y < Overlap.S) && (*point->y != General.S)) {	/*(1) */
-			weight = (Overlap.S - *point->y) / overlap;
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(1) */
+			weight = (*point->y - General.S) / overlap;
 
 			gradient[0] *= weight;
 			gradient[1] *= weight;
@@ -394,11 +388,13 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 	    }
 	}			/*end if obs */
     }				/*end for */
-    G_percent(i, num_points, 6); /* finish it */
+    G_percent(num_points, num_points, 2); /* finish it */
+
+    db_commit_transaction(driver);
 
     Vect_destroy_line_struct(point);
     Vect_destroy_cats_struct(categories);
-}				/*end puntisparsi_select */
+}
 
 int Insert(double partialX, double partialY, double Interp,
            int line_num, dbDriver * driver, char *tab_name)
@@ -409,7 +405,7 @@ int Insert(double partialX, double partialY, double Interp,
 
     db_init_string(&sql);
     sprintf(buf,
-	    "INSERT INTO %s (ID, Interp, partialX, partialY)", tab_name);
+	    "INSERT INTO %s (ID, Interp, X, Y)", tab_name);
     db_append_string(&sql, buf);
     sprintf(buf, " VALUES (%d, %lf, %lf, %lf)", line_num, Interp,
             partialX, partialY);
@@ -450,7 +446,7 @@ int UpDate(double partialX, double partialY, double Interp,
 
     db_init_string(&sql);
     sprintf(buf,
-	    "UPDATE %s SET Interp=%lf, PartialX=%lf, PartialY=%lf WHERE ID=%d",
+	    "UPDATE %s SET Interp=%lf, X=%lf, Y=%lf WHERE ID=%d",
 	    tab_name, Interp, partialX, partialY, line_num);
     db_append_string(&sql, buf);
 
@@ -473,7 +469,7 @@ int Select(double *PartialX, double *PartialY, double *Interp,
 
     db_init_string(&sql);
     sprintf(buf,
-	    "SELECT ID, Interp, partialX, partialY FROM %s WHERE ID=%d",
+	    "SELECT ID, Interp, X, Y FROM %s WHERE ID=%d",
 	    tab_name, line_num);
     db_append_string(&sql, buf);
 
@@ -515,89 +511,30 @@ int Select(double *PartialX, double *PartialY, double *Interp,
     return DB_OK;
 }
 
-int Create_AuxEdge_Table(dbDriver * driver, char *tab_name)
-{
-    dbTable *table;
-    dbColumn *ID_col, *Interp_col, *PartialX_col, *PartialY_col;
-    int created;
-
-    table = db_alloc_table(4);
-    db_set_table_name(table, tab_name);
-    db_set_table_description(table,
-			     "It is used for the intermediate interpolated and gradient values");
-
-    ID_col = db_get_table_column(table, 0);
-    db_set_column_name(ID_col, "ID");
-    db_set_column_sqltype(ID_col, DB_SQL_TYPE_INTEGER);
-
-    Interp_col = db_get_table_column(table, 1);
-    db_set_column_name(Interp_col, "Interp");
-    db_set_column_sqltype(Interp_col, DB_SQL_TYPE_REAL);
-
-    PartialX_col = db_get_table_column(table, 2);
-    db_set_column_name(PartialX_col, "PartialX");
-    db_set_column_sqltype(PartialX_col, DB_SQL_TYPE_REAL);
-
-    PartialY_col = db_get_table_column(table, 3);
-    db_set_column_name(PartialY_col, "PartialY");
-    db_set_column_sqltype(PartialY_col, DB_SQL_TYPE_REAL);
-
-    if (db_create_table(driver, table) == DB_OK) {
-	G_debug(1, "<%s> table created in database.", db_get_table_name(table));
-	created = TRUE;
-    }
-    else
-	return FALSE;
-
-    db_create_index2(driver, tab_name, "ID");
-
-    return created;
-}
-
-int Create_Interpolation_Table(dbDriver * driver, char *tab_name)
-{
-    dbTable *table;
-    dbColumn *ID_col, *Interp_col;
-
-    table = db_alloc_table(2);
-    db_set_table_name(table, tab_name);
-    db_set_table_description(table,
-			     "This table is the bicubic interpolation of the input vector");
-
-    ID_col = db_get_table_column(table, 0);
-    db_set_column_name(ID_col, "ID");
-    db_set_column_sqltype(ID_col, DB_SQL_TYPE_INTEGER);
-
-    Interp_col = db_get_table_column(table, 1);
-    db_set_column_name(Interp_col, "Interp");
-    db_set_column_sqltype(Interp_col, DB_SQL_TYPE_REAL);
-
-    if (db_create_table(driver, table) == DB_OK) {
-	G_debug(1, _("<%s> table created in database."), db_get_table_name(table));
-	return DB_OK;
-    }
-    else
-	return DB_FAILED;
-}
-
-
-#ifdef notdef
 /*! DEFINITION OF THE SUBZONES 
 
-   -----------------------
-   |4|   3   |3|       | |
-   -----------------------
-   | |       | |       | |
-   |2|   5   |1|       | |
-   | |       | |       | |
-   -----------------------
-   |2|   1   |1|       | |
-   -----------------------
-   | |       | |       | |
-   | |       | |       | |
-   | |       | |       | |
-   -----------------------
-   | |       | |       | |
-   -----------------------
+  5: inside Overlap region
+  all others: inside General region but outside Overlap region
+
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       |4|   3   |3|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       |2|   5   |1|       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       |2|   1   |1|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
  */
-#endif

vector/lidar/v.lidar.edgedetection/edgedetection.h

@@ -61,6 +61,4 @@ int UpDate(double, /**/ double, /**/ double, /**/ int, /**/ dbDriver *, /**/ cha
 
 int Select(double *, /**/ double *, /**/ double *, /**/ int, /**/ dbDriver *, /**/ char *);
 
-int Create_AuxEdge_Table(dbDriver *, char *);
-int Create_Interpolation_Table(dbDriver *, char *);
 int Insert_Interpolation(double, int, dbDriver *, char *);

vector/lidar/v.lidar.edgedetection/main.c

@@ -34,7 +34,7 @@
 #include <grass/PolimiFunct.h>
 #include "edgedetection.h"
 
-int nsply, nsplx, line_out_counter, first_it;
+int nsply, nsplx, line_out_counter;
 double passoN, passoE;
 
 /**************************************************************************************
@@ -83,9 +83,9 @@ int main(int argc, char *argv[])
 
     spline_step_flag = G_define_flag();
     spline_step_flag->key = 'e';
-    spline_step_flag->label = _("Estimate spline step value");
+    spline_step_flag->label = _("Estimate point density and distance");
     spline_step_flag->description =
-	_("Estimate a good spline step value for the input vector points within the current region extends and quit");
+	_("Estimate point density and distance for the input vector points within the current region extends and quit");
 
     in_opt = G_define_standard_option(G_OPT_V_INPUT);
 
@@ -224,6 +224,10 @@ int main(int argc, char *argv[])
     if (1 > Vect_open_old(&In, in_opt->answer, mapset))
 	G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
 
+    /* Input vector must be 3D */
+    if (!Vect_is_3d(&In))
+	G_fatal_error(_("Input vector map <%s> is not 3D!"), in_opt->answer);
+
     /* Estimate point density and mean distance for current region */
     if (spline_step_flag->answer) {
 	double dens, dist;
@@ -254,41 +258,55 @@ int main(int argc, char *argv[])
 		      dvr);
 
     /* Create auxiliar and interpolation table */
-    if ((flag_auxiliar = Create_AuxEdge_Table(driver, table_name)) == FALSE)
+    if ((flag_auxiliar = P_Create_Aux4_Table(driver, table_name)) == FALSE)
 	G_fatal_error(_("It was impossible to create <%s>."), table_name);
 
-    if (Create_Interpolation_Table(driver, table_interpolation) != DB_OK)
+    if (P_Create_Aux2_Table(driver, table_interpolation) == FALSE)
 	G_fatal_error(_("It was impossible to create <%s> interpolation table in database."),
 		      out_opt->answer);
 
+    db_create_index2(driver, table_name, "ID");
+    db_create_index2(driver, table_interpolation, "ID");
+    /* sqlite likes that */
+    db_close_database_shutdown_driver(driver);
+    driver = db_start_driver_open_database(dvr, db);
+
     /* Setting regions and boxes */
     G_get_set_window(&original_reg);
     G_get_set_window(&elaboration_reg);
     Vect_region_box(&elaboration_reg, &overlap_box);
     Vect_region_box(&elaboration_reg, &general_box);
 
-    /* Fixing parameters of the elaboration region */
-    /*! Each original_region will be divided into several subregions. These
-     *  subregion will be overlapped by its neighboring subregions. This overlapping
-     *  is calculated as OVERLAP_SIZE times the east-west spline step. */
+    /*------------------------------------------------------------------
+      | Subdividing and working with tiles: 									
+      | Each original region will be divided into several subregions. 
+      | Each one will be overlaped by its neighbouring subregions. 
+      | The overlapping is calculated as a fixed OVERLAP_SIZE times
+      | the largest spline step plus 2 * orlo
+      ----------------------------------------------------------------*/
 
+    /* Fixing parameters of the elaboration region */
     P_zero_dim(&dims);
 
-    N_extension = original_reg.north - original_reg.south;
-    E_extension = original_reg.east - original_reg.west;
-
     nsplx_adj = NSPLX_MAX;
     nsply_adj = NSPLY_MAX;
-    dims.overlap = OVERLAP_SIZE * passoE;
-    P_get_orlo(P_BICUBIC, &dims, passoE, passoN);	/* Set orlo_h|v */
+    if (passoN > passoE)
+	dims.overlap = OVERLAP_SIZE * passoN;
+    else
+	dims.overlap = OVERLAP_SIZE * passoE;
+    P_get_orlo(P_BICUBIC, &dims, passoE, passoN);
     P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
 
     G_verbose_message("adjusted EW splines %d", nsplx_adj);
     G_verbose_message("adjusted NS splines %d", nsply_adj);
 
+    /* calculate number of subzones */
     orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
     orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
 
+    N_extension = original_reg.north - original_reg.south;
+    E_extension = original_reg.east - original_reg.west;
+
     nsubregion_col = ceil(E_extension / orloE) + 0.5;
     nsubregion_row = ceil(N_extension / orloN) + 0.5;
 
@@ -300,9 +318,7 @@ int main(int argc, char *argv[])
     nsubzones = nsubregion_row * nsubregion_col;
 
     elaboration_reg.south = original_reg.north;
-
     last_row = FALSE;
-    first_it = TRUE;
 
     while (last_row == FALSE) {	/* For each row */
 
@@ -383,7 +399,7 @@ int main(int argc, char *argv[])
 		BW = P_get_BandWidth(P_BILINEAR, nsply);	/* Bilinear interpolation */
 		N = G_alloc_matrix(nparameters, BW);	/* Normal matrix */
 		TN = G_alloc_vector(nparameters);	/* vector */
-		parVect_bilin = G_alloc_vector(nparameters);	/* Bicubic parameters vector */
+		parVect_bilin = G_alloc_vector(nparameters);	/* Bilinear parameters vector */
 		obsVect = G_alloc_matrix(npoints + 1, 3);	/* Observation vector */
 		Q = G_alloc_vector(npoints + 1);	/* "a priori" var-cov matrix */
 
@@ -398,7 +414,7 @@ int main(int argc, char *argv[])
 		    Q[i] = 1;	/* Q=I */
 		}
 
-		/*G_free (observ); */
+		G_free(observ);
 
 		G_verbose_message(_("Bilinear interpolation"));
 		normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
@@ -441,9 +457,11 @@ int main(int argc, char *argv[])
 		G_free_matrix(obsVect);
 		G_free_ivector(lineVect);
 	    }			/* IF */
-
-	    G_free(observ);
-	    first_it = FALSE;
+	    else {
+		G_free(observ);
+		G_warning(_("No data within this subzone. "
+			    "Consider changing the spline step."));
+	    }
 	}			/*! END WHILE; last_column = TRUE */
     }				/*! END WHILE; last_row = TRUE */
 

vector/lidar/v.lidar.growing/main.c

@@ -84,6 +84,8 @@ int main(int argc, char *argv[])
 	  "algorithm for determining the building inside");
 
     in_opt = G_define_standard_option(G_OPT_V_INPUT);
+    in_opt->description =
+	_("Input vector (v.lidar.edgedetection output");
 
     out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
 
@@ -166,30 +168,8 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
 		      field->driver);
 
-    /* Setting regions and boxes */
-    G_debug(1, _("Setting regions and boxes"));
-    G_get_set_window(&original_reg);
-    G_get_set_window(&elaboration_reg);
-
-    /*  Fixing parameters of the elaboration region */
-    /*! The original_region will be divided into several subregions */
-    ew_resol = original_reg.ew_res;
-    ns_resol = original_reg.ns_res;
-
-    nsubregion_col = ceil((original_reg.east - original_reg.west) / (LATO * ew_resol)) + 0.5;
-    nsubregion_row = ceil((original_reg.north - original_reg.south) / (LATO * ns_resol)) + 0.5;
-
-    if (nsubregion_col < 0)
-	nsubregion_col = 0;
-    if (nsubregion_row < 0)
-	nsubregion_row = 0;
-
-    nsubzones = nsubregion_row * nsubregion_col;
-
-    /* Subdividing and working with tiles */
-    elaboration_reg.south = original_reg.north;
-    last_row = FALSE;
-
+    /* is this the right place to open the cursor ??? */
+    
     db_init_string(&sql);
     db_zero_string(&sql);
 
@@ -202,6 +182,7 @@ int main(int argc, char *argv[])
 
     count_obj = 1;
 
+    /* no topology, get number of lines in input vector */
     nlines = 0;
     points = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
@@ -211,6 +192,7 @@ int main(int argc, char *argv[])
     }
     Vect_rewind(&In);
 
+    /* no topology, get number of lines in first pulse input vector */
     nlines_first = 0;
     points_first = Vect_new_line_struct();
     Cats_first = Vect_new_cats_struct();
@@ -220,6 +202,31 @@ int main(int argc, char *argv[])
     }
     Vect_rewind(&First);
 
+    /* Setting regions and boxes */
+    G_debug(1, _("Setting regions and boxes"));
+    G_get_set_window(&original_reg);
+    G_get_set_window(&elaboration_reg);
+
+    /*  Fixing parameters of the elaboration region */
+    /*! The original_region will be divided into subregions */
+    ew_resol = original_reg.ew_res;
+    ns_resol = original_reg.ns_res;
+
+    /* calculate number of subzones */
+    nsubregion_col = ceil((original_reg.east - original_reg.west) / (LATO * ew_resol)) + 0.5;
+    nsubregion_row = ceil((original_reg.north - original_reg.south) / (LATO * ns_resol)) + 0.5;
+
+    if (nsubregion_col < 0)
+	nsubregion_col = 0;
+    if (nsubregion_row < 0)
+	nsubregion_row = 0;
+
+    nsubzones = nsubregion_row * nsubregion_col;
+
+    /* Subdividing and working with tiles */
+    elaboration_reg.south = original_reg.north;
+    last_row = FALSE;
+
     while (last_row == FALSE) {	/* For each strip of LATO rows */
 
 	elaboration_reg.north = elaboration_reg.south;
@@ -259,6 +266,8 @@ int main(int argc, char *argv[])
 	    elaboration_reg.ew_res = ew_resol;
 	    nrows = (elaboration_reg.north - elaboration_reg.south) / ns_resol + 0.1;
 	    ncols = (elaboration_reg.east - elaboration_reg.west) / ew_resol + 0.1;
+	    elaboration_reg.rows = nrows;
+	    elaboration_reg.cols = ncols;
 
 	    G_debug(1, _("Rows = %d"), nrows);
 	    G_debug(1, _("Columns = %d"), ncols);
@@ -281,13 +290,13 @@ int main(int argc, char *argv[])
 		}
 	    }
 
+	    G_verbose_message(_("read points in input vector"));
+	    Vect_region_box(&elaboration_reg, &elaboration_box);
 	    line_num = 0;
 	    Vect_rewind(&In);
 	    while (Vect_read_next_line(&In, points, Cats) > 0) {
 		line_num++;
 
-		Vect_region_box(&elaboration_reg, &elaboration_box);
-
 		if ((Vect_point_in_box
 		     (points->x[0], points->y[0], points->z[0],
 		      &elaboration_box)) &&
@@ -304,6 +313,36 @@ int main(int argc, char *argv[])
 			      (points->x[0], &elaboration_reg));
 
 		    Z_interp = 0;
+		    /* TODO: make sure the current db_fetch() usage works */
+		    /* why not: */
+		    /*
+		    db_init_string(&sql);
+		    sprintf(buf, "SELECT Interp,ID FROM %s WHERE ID=%d", table_name, line_num);
+		    db_append_string(&sql, buf);
+
+		    if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
+			G_fatal_error(_("Unable to open table <%s>"), table_name);
+
+		    while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
+			dbColumn *Z_Interp_col;
+			dbValue *Z_Interp_value;
+			table = db_get_cursor_table(&cursor);
+
+			Z_Interp_col = db_get_table_column(table, 1);
+
+			if (db_sqltype_to_Ctype(db_get_column_sqltype(Z_Interp_col)) ==
+			    DB_C_TYPE_DOUBLE)
+			    Z_Interp_value = db_get_column_value(Z_Interp_col);
+			else
+			    continue;
+
+			Z_interp = db_get_value_double(Z_Interp_value);
+			break;
+		    }
+		    db_close_cursor(&cursor);
+		    db_free_string(&sql);
+		    */
+		    /* instead of */
 		    while (1) {
 			if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK ||
 			    !more)
@@ -400,7 +439,7 @@ int main(int argc, char *argv[])
 
 	    /* REGION GROWING */
 	    if (region == TRUE) {
-		G_debug(1, "Region Growing");
+		G_verbose_message(_("Region Growing"));
 
 		punti_bordo = G_alloc_matrix(MaxPoints, 3);
 		P = Pvector(0, MaxPoints);
@@ -409,6 +448,7 @@ int main(int argc, char *argv[])
 		ripieno = 6;
 
 		for (row = 0; row <= nrows; row++) {
+		    G_percent(row, nrows, 2);
 		    for (col = 0; col <= ncols; col++) {
 
 			if ((raster_matrix[row][col].clas >= Thres_j) &&
@@ -437,7 +477,7 @@ int main(int argc, char *argv[])
 				     punti_bordo, &lungPunti, row, col,
 				     colorBordo, Thres_j, MaxPoints);
 
-			    /*CONVEX-HULL COMPUTATION */
+			    /* CONVEX-HULL COMPUTATION */
 			    lungHull = ch2d(P, lungPunti);
 			    cvxHull = G_alloc_matrix(lungHull, 3);
 

vector/lidar/v.outlier/main.c

@@ -29,7 +29,7 @@
 #include <grass/PolimiFunct.h>
 #include "outlier.h"
     /* GLOBAL VARIABLES DEFINITIONS */
-int nsply, nsplx, first_it;
+int nsply, nsplx;
 double passoN, passoE, Thres_Outlier;
 
 /*--------------------------------------------------------------------------------------*/
@@ -41,7 +41,7 @@ int main(int argc, char *argv[])
     int subzone = 0, nsubzones = 0;
     double N_extension, E_extension, orloE, orloN;
     int dim_vect, nparameters, BW, npoints;
-    double ew_resol, ns_resol, mean, lambda;
+    double mean, lambda;
     const char *dvr, *db, *mapset;
     char table_name[GNAME_MAX];
     char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
@@ -76,9 +76,9 @@ int main(int argc, char *argv[])
 
     spline_step_flag = G_define_flag();
     spline_step_flag->key = 'e';
-    spline_step_flag->label = _("Estimate spline step value");
+    spline_step_flag->label = _("Estimate point density and distance");
     spline_step_flag->description =
-	_("Estimate a good spline step value for the input vector points within the current region extends and quit");
+	_("Estimate point density and distance for the input vector points within the current region extends and quit");
 
     in_opt = G_define_standard_option(G_OPT_V_INPUT);
 
@@ -147,6 +147,8 @@ int main(int argc, char *argv[])
     lambda = atof(lambda_f_opt->answer);
     Thres_Outlier = atof(Thres_O_opt->answer);
 
+    flag_auxiliar = FALSE;
+
     /* Checking vector names */
     Vect_check_input_output_name(in_opt->answer, out_opt->answer,
 				 GV_FATAL_EXIT);
@@ -180,6 +182,10 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
 		      in_opt->answer);
 
+    /* Input vector must be 3D */
+    if (!Vect_is_3d(&In))
+	G_fatal_error(_("Input vector map <%s> is not 3D!"), in_opt->answer);
+
     /* Estimate point density and mean distance for current region */
     if (spline_step_flag->answer) {
 	double dens, dist;
@@ -232,37 +238,52 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
 		      dvr);
 
+    /* Create auxiliar table */
+    if ((flag_auxiliar =
+	 P_Create_Aux2_Table(driver, table_name)) == FALSE)
+	G_fatal_error(_("It was impossible to create <%s> table."), table_name);
+
+    db_create_index2(driver, table_name, "ID");
+    /* sqlite likes that */
+    db_close_database_shutdown_driver(driver);
+    driver = db_start_driver_open_database(dvr, db);
+
     /* Setting regions and boxes */
     G_get_set_window(&original_reg);
     G_get_set_window(&elaboration_reg);
     Vect_region_box(&elaboration_reg, &overlap_box);
     Vect_region_box(&elaboration_reg, &general_box);
 
-    /* Fixing parameters of the elaboration region */
-    /*! Each original_region will be divided into several subregions. These
-     *  subregion will be overlapped by its neibourgh subregions. This overlapping
-     *  is calculated as OVERLAP_SIZE times the east-west spline step. */
-
-    ew_resol = original_reg.ew_res;
-    ns_resol = original_reg.ns_res;
+    /*------------------------------------------------------------------
+      | Subdividing and working with tiles: 									
+      | Each original region will be divided into several subregions. 
+      | Each one will be overlaped by its neighbouring subregions. 
+      | The overlapping is calculated as a fixed OVERLAP_SIZE times
+      | the largest spline step plus 2 * orlo
+      ----------------------------------------------------------------*/
 
+    /* Fixing parameters of the elaboration region */
     P_zero_dim(&dims);
 
-    N_extension = original_reg.north - original_reg.south;
-    E_extension = original_reg.east - original_reg.west;
-
     nsplx_adj = NSPLX_MAX;
     nsply_adj = NSPLY_MAX;
-    dims.overlap = OVERLAP_SIZE * passoE;
+    if (passoN > passoE)
+	dims.overlap = OVERLAP_SIZE * passoN;
+    else
+	dims.overlap = OVERLAP_SIZE * passoE;
     P_get_orlo(P_BILINEAR, &dims, passoE, passoN);
     P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
 
     G_verbose_message("adjusted EW splines %d", nsplx_adj);
     G_verbose_message("adjusted NS splines %d", nsply_adj);
 
+    /* calculate number of subzones */
     orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
     orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
 
+    N_extension = original_reg.north - original_reg.south;
+    E_extension = original_reg.east - original_reg.west;
+
     nsubregion_col = ceil(E_extension / orloE) + 0.5;
     nsubregion_row = ceil(N_extension / orloN) + 0.5;
 
@@ -274,9 +295,7 @@ int main(int argc, char *argv[])
     nsubzones = nsubregion_row * nsubregion_col;
 
     elaboration_reg.south = original_reg.north;
-
     last_row = FALSE;
-    first_it = TRUE;
 
     while (last_row == FALSE) {	/* For each row */
 
@@ -370,6 +389,8 @@ int main(int argc, char *argv[])
 		    Q[i] = 1;	/* Q=I */
 		}
 
+		G_free(observ);
+
 		G_verbose_message(_("Bilinear interpolation"));
 		normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
 			       nsply, elaboration_reg.west,
@@ -382,14 +403,6 @@ int main(int argc, char *argv[])
 		G_free_vector(TN);
 		G_free_vector(Q);
 
-		if (flag_auxiliar == FALSE) {
-		    G_debug(1,
-			    "Creating auxiliar table for archiving overlapping zones");
-		    if ((flag_auxiliar =
-			 P_Create_Aux_Table(driver, table_name)) == FALSE)
-			G_fatal_error(_("It was impossible to create <Auxiliar_outlier_table>."));
-		}
-
 		G_verbose_message(_("Outlier detection"));
 		if (qgis_opt->answer)
 		    P_Outlier(&Out, &Outlier, &Qgis, elaboration_reg,
@@ -406,8 +419,11 @@ int main(int argc, char *argv[])
 		G_free_ivector(lineVect);
 
 	    }			/*! END IF; npoints > 0 */
-	    G_free(observ);
-	    first_it = FALSE;
+	    else {
+		G_free(observ);
+		G_warning(_("No data within this subzone. "
+			    "Consider changing the spline step."));
+	    }
 	}			/*! END WHILE; last_column = TRUE */
     }				/*! END WHILE; last_row = TRUE */
 

vector/lidar/v.outlier/outlier.c

@@ -27,6 +27,8 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
     point = Vect_new_line_struct();
     categories = Vect_new_cats_struct();
 
+    db_begin_transaction(driver);
+    
     for (i = 0; i < num_points; i++) {	/* Sparse points */
 	G_percent(i, num_points, 2);
 	Vect_reset_line(point);
@@ -60,12 +62,11 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 
 	    }
 	    else {
-		if ((*point->x > Overlap.E) && (*point->x != General.E)) {
-
-		    if ((*point->y > Overlap.N) && (*point->y != General.N)) {	/*(3) */
-			csi = (*point->x - Overlap.E) / overlap;
-			eta = (*point->y - Overlap.N) / overlap;
-			weight = (1 - csi) * (1 - eta);
+		if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(3) */
+			csi = (General.E - *point->x) / overlap;
+			eta = (General.N - *point->y) / overlap;
+			weight = csi * eta;
 
 			interpolation *= weight;
 
@@ -76,22 +77,20 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 			if (UpDate_Outlier(interpolation, line_num[i],
 			                   driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to update the database"));
-
 		    }
-		    else if ((*point->y < Overlap.S) && (*point->y != General.S)) {	/*(1) */
-			csi = (*point->x - Overlap.E) / overlap;
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(1) */
+			csi = (General.E - *point->x) / overlap;
 			eta = (*point->y - General.S) / overlap;
-			weight = (1 - csi) * eta;
+			weight = csi * eta;
 
 			interpolation *= weight;
 
 			if (Insert_Outlier(interpolation, line_num[i],
 					   driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to write in the database"));
-
 		    }
 		    else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) {	/*(1) */
-			weight = (*point->x - Overlap.E) / overlap;
+			weight = (General.E - *point->x) / overlap;
 
 			interpolation *= weight;
 
@@ -99,14 +98,12 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 					   driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to write in the database"));
 		    }
-
 		}
-		else if ((*point->x < Overlap.W) && (*point->x != General.W)) {
-
-		    if ((*point->y > Overlap.N) && (*point->y != General.N)) {	/*(4) */
+		else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(4) */
 			csi = (*point->x - General.W) / overlap;
-			eta = (*point->y - Overlap.N) / overlap;
-			weight = (1 - eta) * csi;
+			eta = (General.N - *point->y) / overlap;
+			weight = eta * csi;
 
 			interpolation *= weight;
 
@@ -123,12 +120,13 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 				Vect_write_line(Qgis, GV_POINT, point,
 						categories);
 			}
-			else
+			else {
 			    Vect_write_line(Outlier, GV_POINT, point,
 					    categories);
-
+			    G_message("here we are");
+			}
 		    }
-		    else if ((*point->y < Overlap.S) && (*point->y != General.S)) {	/*(2) */
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(2) */
 			csi = (*point->x - General.W) / overlap;
 			eta = (*point->y - General.S) / overlap;
 			weight = csi * eta;
@@ -142,10 +140,9 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 			if (UpDate_Outlier(interpolation, line_num[i],
 			                   driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to update the database"));
-
 		    }
 		    else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) {	/*(2) */
-			weight = (Overlap.W - *point->x) / overlap;
+			weight = (*point->x - General.W) / overlap;
 
 			interpolation *= weight;
 
@@ -166,11 +163,10 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 			    Vect_write_line(Outlier, GV_POINT, point,
 					    categories);
 		    }
-
 		}
 		else if ((*point->x <= Overlap.E) && (*point->x >= Overlap.W)) {
-		    if ((*point->y > Overlap.N) && (*point->y != General.N)) {	/*(3) */
-			weight = (*point->y - Overlap.N) / overlap;
+		    if ((*point->y > Overlap.N) && (*point->y < General.N)) {	/*(3) */
+			weight = (General.N - *point->y) / overlap;
 
 			interpolation *= weight;
 
@@ -190,10 +186,9 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 			else
 			    Vect_write_line(Outlier, GV_POINT, point,
 					    categories);
-
 		    }
-		    else if ((*point->y < Overlap.S) && (*point->y != General.S)) {	/*(1) */
-			weight = (Overlap.S - *point->y) / overlap;
+		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(1) */
+			weight = (*point->y - General.S) / overlap;
 
 			interpolation *= weight;
 
@@ -210,6 +205,8 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
     G_percent(num_points, num_points, 2);
     G_debug(2, "P_outlier: done");
 
+    db_commit_transaction(driver);
+
     Vect_destroy_line_struct(point);
     Vect_destroy_cats_struct(categories);
 }				/*end puntisparsi_select */
@@ -295,23 +292,30 @@ int P_is_outlier(double pippo)
     return TRUE;
 }
 
-#ifdef notdef
 /*! DEFINITION OF THE SUBZONES 
 
-   -----------------------
-   |4|   3   |3|       | |
-   -----------------------
-   | |       | |       | |
-   |2|   5   |1|       | |
-   | |       | |       | |
-   -----------------------
-   |2|   1   |1|       | |
-   -----------------------
-   | |       | |       | |
-   | |       | |       | |
-   | |       | |       | |
-   -----------------------
-   | |       | |       | |
-   -----------------------
+  5: inside Overlap region
+  all others: inside General region but outside Overlap region
+
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       |4|   3   |3|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       |2|   5   |1|       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       |2|   1   |1|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
  */
-#endif

vector/lidar/v.surf.bspline/crosscorr.c

@@ -49,7 +49,8 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
     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.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 *TN, *Q, *parVect;	/* Interpolation and least-square vectors */
     double **N, **obsVect;	/* Interpolation and least-square matrix */

vector/lidar/v.surf.bspline/main.c

@@ -38,7 +38,7 @@ char *bspline_column;
 int main(int argc, char *argv[])
 {
     /* Variable declarations */
-    int nsply, nsplx, nlines, nrows, ncols, nsplx_adj, nsply_adj;
+    int nsply, nsplx, nrows, ncols, nsplx_adj, nsply_adj;
     int nsubregion_col, nsubregion_row, subregion_row, subregion_col;
     int subzone = 0, nsubzones = 0;
     int last_row, last_column, grid, bilin, ext, flag_auxiliar, cross;	/* booleans */
@@ -225,8 +225,12 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
 		      in_opt->answer);
 
-    if (!Vect_is_3d(&In) && bspline_field <= 0 && bspline_column)
+    /* check availability of z values
+     * column option overrrides 3D z coordinates */
+    if (!Vect_is_3d(&In) && bspline_field <= 0 && bspline_column == NULL)
 	G_fatal_error(_("Need either 3D vector or layer and column with z values"));
+    if (bspline_field > 0 && bspline_column == NULL)
+	G_fatal_error(_("Layer but not column with z values given"));
 
     /* Estimate point density and mean distance for current region */
     if (spline_step_flag->answer) {
@@ -316,6 +320,7 @@ int main(int argc, char *argv[])
 	raster = Rast_open_fp_new(out_map_opt->answer);
     }
 
+    /* read z values from attribute table */
     if (bspline_field > 0) {
 	db_CatValArray_init(&cvarr);
 	Fi = Vect_get_field(&In, bspline_field);
@@ -356,14 +361,18 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("No database connection for driver <%s> is defined. "
 			"Run db.connect."), dvr);
 
-    /* Auxiliar table creation */
+    /* Create auxiliar table */
     if (vector) {
-	if ((flag_auxiliar = P_Create_Aux_Table(driver, table_name)) == FALSE) {
+	if ((flag_auxiliar = P_Create_Aux4_Table(driver, table_name)) == FALSE) {
 	    P_Drop_Aux_Table(driver, table_name);
 	    G_fatal_error(_("Interpolation: Creating table: "
 			    "It was impossible to create table <%s>."),
 			  table_name);
 	}
+	db_create_index2(driver, table_name, "ID");
+	/* sqlite likes that */
+	db_close_database_shutdown_driver(driver);
+	driver = db_start_driver_open_database(dvr, db);
     }
 
     /* Setting regions and boxes */
@@ -383,24 +392,36 @@ int main(int argc, char *argv[])
 			    "Consider changing region resolution"));
     }
 
-    /* Fixing parameters of the elaboration region */
-    P_zero_dim(&dims);		/* Set to zero the dim struct */
+    /*------------------------------------------------------------------
+      | Subdividing and working with tiles: 									
+      | Each original region will be divided into several subregions. 
+      | Each one will be overlaped by its neighbouring subregions. 
+      | The overlapping is calculated as a fixed OVERLAP_SIZE times
+      | the largest spline step plus 2 * orlo
+      ----------------------------------------------------------------*/
 
-    N_extension = original_reg.north - original_reg.south;
-    E_extension = original_reg.east - original_reg.west;
+    /* Fixing parameters of the elaboration region */
+    P_zero_dim(&dims);		/* Set dim struct to zero */
 
     nsplx_adj = NSPLX_MAX;
     nsply_adj = NSPLY_MAX;
-    dims.overlap = OVERLAP_SIZE * passoE;
-    P_get_orlo(bilin, &dims, passoE, passoN);	/* Set orlo_h|v */
+    if (passoN > passoE)
+	dims.overlap = OVERLAP_SIZE * passoN;
+    else
+	dims.overlap = OVERLAP_SIZE * passoE;
+    P_get_orlo(bilin, &dims, passoE, passoN);
     P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
 
     G_verbose_message("adjusted EW splines %d", nsplx_adj);
     G_verbose_message("adjusted NS splines %d", nsply_adj);
 
+    /* calculate number of subzones */
     orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
     orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
 
+    N_extension = original_reg.north - original_reg.south;
+    E_extension = original_reg.east - original_reg.west;
+
     nsubregion_col = ceil(E_extension / orloE) + 0.5;
     nsubregion_row = ceil(N_extension / orloN) + 0.5;
 
@@ -415,21 +436,11 @@ int main(int argc, char *argv[])
     points = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
     Vect_cat_set(Cats, 1, 0);
-    nlines = Vect_get_num_lines(&In);
-
-    /*------------------------------------------------------------------
-      | Subdividing and working with tiles: 									
-      | Each original_region will be divided into several subregions. 
-      | Each one will be overlaped by its neibourgh subregions. 
-      | The overlapping was calculated as a fixed OVERLAP_SIZE times
-      | the east-west spline step
-      ----------------------------------------------------------------*/
-
-    elaboration_reg.south = original_reg.north;
 
-    G_percent(0, 1, 10);
     subregion_row = 0;
+    elaboration_reg.south = original_reg.north;
     last_row = FALSE;
+
     while (last_row == FALSE) {	/* For each subregion row */
 	subregion_row++;
 	P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
@@ -511,7 +522,6 @@ int main(int argc, char *argv[])
 
 	    if (npoints > 0) {	/*  */
 		int i;
-		double *obs_mean;
 
 		nparameters = nsplx * nsply;
 		BW = P_get_BandWidth(P_BILINEAR, nsply) + 2 * (bilin == P_BICUBIC);
@@ -523,7 +533,6 @@ int main(int argc, char *argv[])
 		obsVect = G_alloc_matrix(npoints, 3);	/* Observation vector */
 		Q = G_alloc_vector(npoints);	/* "a priori" var-cov matrix */
 		lineVect = G_alloc_ivector(npoints);	/*  */
-		obs_mean = G_alloc_vector(npoints);
 
 		for (i = 0; i < npoints; i++) {	/* Setting obsVect vector & Q matrix */
 		    double dval;
@@ -533,6 +542,7 @@ int main(int argc, char *argv[])
 		    obsVect[i][0] = observ[i].coordX;
 		    obsVect[i][1] = observ[i].coordY;
 
+		    /* read z coordinates from attribute table */
 		    if (bspline_field > 0) {
 			int cat, ival, ret;
 
@@ -545,14 +555,14 @@ int main(int argc, char *argv[])
 				db_CatValArray_get_value_int(&cvarr, cat,
 							     &ival);
 			    obsVect[i][2] = ival;
-			    obs_mean[i] = ival;
+			    observ[i].coordZ = ival;
 			}
 			else {	/* DB_C_TYPE_DOUBLE */
 			    ret =
 				db_CatValArray_get_value_double(&cvarr, cat,
 								&dval);
 			    obsVect[i][2] = dval;
-			    obs_mean[i] = dval;
+			    observ[i].coordZ = dval;
 			}
 			if (ret != DB_OK) {
 			    G_warning(_("Interpolation: (%d,%d): No record for point (cat = %d)"),
@@ -560,29 +570,25 @@ int main(int argc, char *argv[])
 			    continue;
 			}
 		    }
-
+		    /* use z coordinates of 3D vector */
 		    else {
 			obsVect[i][2] = observ[i].coordZ;
-			obs_mean[i] = observ[i].coordZ;
-		    }		/* obsVect[i][2] = observ[i].coordZ - mean; */
+		    }
 		}
 
 		/* Mean calculation for every point */
-		if (bspline_field > 0)
-		    mean = calc_mean(obs_mean, npoints);
-		else
-		    mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
-
+		mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
 
 		G_debug(1, "Interpolation: (%d,%d): mean=%lf",
 			subregion_row, subregion_col, mean);
 
+		G_free(observ);
+
 		for (i = 0; i < npoints; i++)
 		    obsVect[i][2] -= mean;
 
-		G_free(observ);
-
-		if (bilin) {	/* Bilinear interpolation */
+		/* Bilinear interpolation */
+		if (bilin) {
 		    G_debug(1,
 			    "Interpolation: (%d,%d): Bilinear interpolation...",
 			    subregion_row, subregion_col);
@@ -592,6 +598,7 @@ int main(int argc, char *argv[])
 				   nparameters, BW);
 		    nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
 		}
+		/* Bicubic interpolation */
 		else {
 		    G_debug(1,
 			    "Interpolation: (%d,%d): Bicubic interpolation...",
@@ -610,7 +617,6 @@ int main(int argc, char *argv[])
 		G_free_vector(Q);
 
 		if (grid == TRUE) {	/* GRID INTERPOLATION ==> INTERPOLATION INTO A RASTER */
-		    flag_auxiliar = TRUE;
 		    G_debug(1, "Interpolation: (%d,%d): Regular_Points...",
 			    subregion_row, subregion_col);
 		    raster_matrix =
@@ -669,11 +675,12 @@ int main(int argc, char *argv[])
 		G_free_vector(parVect);
 		G_free_matrix(obsVect);
 		G_free_ivector(lineVect);
-		G_free_vector(obs_mean);
 	    }
-	    else
+	    else {
+		G_free(observ);
 		G_warning(_("No data within this subzone. "
 			    "Consider changing the spline step."));
+	    }
 	}			/*! END WHILE; last_column = TRUE */
     }				/*! END WHILE; last_row = TRUE */