[GRASS-SVN] r41063 - in grass/branches/releasebranch_6_4/vector/lidar: lidarlib v.lidar.correction v.lidar.edgedetection v.lidar.growing v.outlier v.surf.bspline

svn_grass at osgeo.org svn_grass at osgeo.org
Wed Feb 17 12:59:30 EST 2010


Author: mmetz
Date: 2010-02-17 12:59:30 -0500 (Wed, 17 Feb 2010)
New Revision: 41063

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/raster.c
   grass/branches/releasebranch_6_4/vector/lidar/lidarlib/zones.c
   grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/correction.c
   grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/correction.h
   grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/main.c
   grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/edgedetection.c
   grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/edgedetection.h
   grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/main.c
   grass/branches/releasebranch_6_4/vector/lidar/v.lidar.growing/main.c
   grass/branches/releasebranch_6_4/vector/lidar/v.outlier/main.c
   grass/branches/releasebranch_6_4/vector/lidar/v.outlier/outlier.c
   grass/branches/releasebranch_6_4/vector/lidar/v.outlier/outlier.h
   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:
fix lidar tools, backport from trunk r41061

Modified: grass/branches/releasebranch_6_4/vector/lidar/lidarlib/PolimiFunct.h
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/lidarlib/PolimiFunct.h	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/lidarlib/PolimiFunct.h	2010-02-17 17:59:30 UTC (rev 41063)
@@ -28,36 +28,36 @@
 #include <grass/gmath.h>
 
 /*----------------------------------------------------------------------------------------------------------*/
-/*CONSTANS DECLARATION */
+/*CONSTANTS DECLARATION */
 
-#define NSPLX_MAX 		200	/* Maximum number of splines along East direction used in the subregions interpolation */
-#define NSPLY_MAX		200	/* Maximum number of splines along North direction used in the subregions interpolation */
-#define OVERLAP_SIZE 		20	/* Subregions overlapping size. */
-#define LATO 			2000	/* Side's size for v.to.qrast command. */
-#define CONTOUR		15 /**/
+#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 		     1000	/* Side's size for v.lidar.growing. */
+#define CONTOUR		       15 	/**/
 #define GENERAL_ROW 		0
-#define GENERAL_COLUMN 	1
+#define GENERAL_COLUMN 	        1
 #define FIRST_ROW 		2
 #define LAST_ROW 		3
 #define FIRST_COLUMN	 	4
 #define LAST_COLUMN 		5
     /* FIELDS ID */
 #define F_EDGE_DETECTION_CLASS	1
-#define F_CLASSIFICATION		2
+#define F_CLASSIFICATION	2
 #define F_INTERPOLATION		3
-#define F_COUNTER_OBJ			4
+#define F_COUNTER_OBJ		4
     /* PRE-CLASSIFICATION */
 #define PRE_TERRAIN 		1
 #define PRE_EDGE		2
 #define PRE_UNKNOWN		3
     /* FINAL CLASSIFICATION */
-#define TERRAIN_SINGLE 	1
-#define TERRAIN_DOUBLE	2
+#define TERRAIN_SINGLE 	        1
+#define TERRAIN_DOUBLE	        2
 #define OBJECT_DOUBLE		3
 #define OBJECT_SINGLE		4
     /* SINGLE OR DOUBLE PULSE */
-#define SINGLE_PULSE	1
-#define DOUBLE_PULSE	2
+#define SINGLE_PULSE	        1
+#define DOUBLE_PULSE	        2
     /* INTERPOLATOR */
 #define P_BILINEAR 		1
 #define P_BICUBIC 		0
@@ -96,6 +96,7 @@
 /*FUNCTIONS DECLARATION */
 /*zones */
 void P_zero_dim(struct Reg_dimens * /**/);
+int P_set_dim(struct Reg_dimens *, double, double, int *, int *);
 
 int P_set_regions(struct Cell_head *, /**/
 		  BOUND_BOX *, /**/
@@ -105,6 +106,8 @@
 
 int P_get_BandWidth(int, /**/ int /**/);
 
+double P_estimate_splinestep(struct Map_info *, double *, double *);
+
 struct Point *P_Read_Vector_Region_Map(struct Map_info *, /**/
 				       struct Cell_head *, /**/
 				       int *, /**/ int, /**/ int /**/);
@@ -142,8 +145,10 @@
 			  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 * /**/);
 
 /*----------------------------------------------------------------------------------------------------------*/

Modified: grass/branches/releasebranch_6_4/vector/lidar/lidarlib/TcholBand.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/lidarlib/TcholBand.c	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/lidarlib/TcholBand.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -9,15 +9,20 @@
 
 void tcholDec(double **N, double **T, int n, int BW)
 {
-    int i, j, k;
+    int i, j, k, end;
     double somma;
 
+    G_debug(2, "tcholDec(): n=%d  BW=%d", n, BW);
+
     for (i = 0; i < n; i++) {
+	G_percent(i, n, 2);
 	for (j = 0; j < BW; j++) {
 	    somma = N[i][j];
-	    for (k = 1; k < BW; k++)
-		if (((i - k) >= 0) && ((j + k) < BW))
-		    somma -= T[i - k][k] * T[i - k][j + k];
+	    /* start = 1 */
+	    /* end = BW - j or i + 1 */
+	    end = ((BW - j) < (i + 1) ? (BW - j) : (i + 1));
+	    for (k = 1; k < end; k++)
+		somma -= T[i - k][k] * T[i - k][j + k];
 	    if (j == 0) {
 		if (somma <= 0.0)
 		    G_fatal_error(_("Decomposition failed"));
@@ -28,6 +33,7 @@
 	}
     }
 
+    G_percent(i, n, 2);
     return;
 }
 
@@ -38,7 +44,7 @@
 {
 
     double **T;
-    int i, j;
+    int i, j, start, end;
 
 	/*--------------------------------------*/
     T = G_alloc_matrix(n, BW);
@@ -50,18 +56,22 @@
     parVect[0] = TN[0] / T[0][0];
     for (i = 1; i < n; i++) {
 	parVect[i] = TN[i];
-	for (j = 0; j < i; j++)
-	    if ((i - j) < BW)
-		parVect[i] -= T[j][i - j] * parVect[j];
+	/* start = 0 or i - BW + 1 */
+	start = ((i - BW + 1) < 0 ? 0 : (i - BW + 1));
+	/* end = i */
+	for (j = start; j < i; j++)
+	    parVect[i] -= T[j][i - j] * parVect[j];
 	parVect[i] = parVect[i] / T[i][0];
     }
 
     /* Backward substitution */
     parVect[n - 1] = parVect[n - 1] / T[n - 1][0];
     for (i = n - 2; i >= 0; i--) {
-	for (j = i + 1; j < n; j++)
-	    if ((j - i) < BW)
-		parVect[i] -= T[i][j - i] * parVect[j];
+	/* start = i + 1 */
+	/* end = n or BW + i */
+	end = (n < (BW + i) ? n : (BW + i));
+	for (j = i + 1; j < end; j++)
+	    parVect[i] -= T[i][j - i] * parVect[j];
 	parVect[i] = parVect[i] / T[i][0];
     }
 
@@ -80,24 +90,28 @@
 		 int BW)
 {
 
-    int i, j;
+    int i, j, start, end;
 
     /* Forward substitution */
     parVect[0] = TN[0] / T[0][0];
     for (i = 1; i < n; i++) {
 	parVect[i] = TN[i];
-	for (j = 0; j < i; j++)
-	    if ((i - j) < BW)
-		parVect[i] -= T[j][i - j] * parVect[j];
+	/* start = 0 or i - BW + 1 */
+	start = ((i - BW + 1) < 0 ? 0 : (i - BW + 1));
+	/* end = i */
+	for (j = start; j < i; j++)
+	    parVect[i] -= T[j][i - j] * parVect[j];
 	parVect[i] = parVect[i] / T[i][0];
     }
 
     /* Backward substitution */
     parVect[n - 1] = parVect[n - 1] / T[n - 1][0];
     for (i = n - 2; i >= 0; i--) {
-	for (j = i + 1; j < n; j++)
-	    if ((j - i) < BW)
-		parVect[i] -= T[i][j - i] * parVect[j];
+	/* start = i + 1 */
+	/* end = n or BW + i */
+	end = (n < (BW + i) ? n : (BW + i));
+	for (j = i + 1; j < end; j++)
+	    parVect[i] -= T[i][j - i] * parVect[j];
 	parVect[i] = parVect[i] / T[i][0];
     }
 
@@ -111,7 +125,7 @@
 {
     double **T = NULL;
     double *vect = NULL;
-    int i, j, k;
+    int i, j, k, start;
     double somma;
 
 	/*--------------------------------------*/
@@ -132,9 +146,11 @@
 	invNdiag[i] = vect[0] * vect[0];
 	for (j = i + 1; j < n; j++) {
 	    somma = 0.0;
-	    for (k = i; k < j; k++) {
-		if ((j - k) < BW)
-		    somma -= vect[k - i] * T[k][j - k];
+	    /* start = i or j - BW + 1 */
+	    start = ((j - BW + 1) < i ? i : (j - BW + 1));
+	    /* end = j */
+	    for (k = start; k < j; k++) {
+		somma -= vect[k - i] * T[k][j - k];
 	    }
 	    vect[j - i] = somma * T[j][0];
 	    invNdiag[i] += vect[j - i] * vect[j - i];
@@ -157,7 +173,7 @@
 
     double **T = NULL;
     double *vect = NULL;
-    int i, j, k;
+    int i, j, k, start, end;
     double somma;
 
 	/*--------------------------------------*/
@@ -171,18 +187,22 @@
     parVect[0] = TN[0] / T[0][0];
     for (i = 1; i < n; i++) {
 	parVect[i] = TN[i];
-	for (j = 0; j < i; j++)
-	    if ((i - j) < BW)
-		parVect[i] -= T[j][i - j] * parVect[j];
+	/* start = 0 or i - BW + 1 */
+	start = ((i - BW + 1) < 0 ? 0 : (i - BW + 1));
+	/* end = i */
+	for (j = start; j < i; j++)
+	    parVect[i] -= T[j][i - j] * parVect[j];
 	parVect[i] = parVect[i] / T[i][0];
     }
 
     /* Backward substitution */
     parVect[n - 1] = parVect[n - 1] / T[n - 1][0];
     for (i = n - 2; i >= 0; i--) {
-	for (j = i + 1; j < n; j++)
-	    if ((j - i) < BW)
-		parVect[i] -= T[i][j - i] * parVect[j];
+	/* start = i + 1 */
+	/* end = n or BW + i */
+	end = (n < (BW + i) ? n : (BW + i));
+	for (j = i + 1; j < end; j++)
+	    parVect[i] -= T[i][j - i] * parVect[j];
 	parVect[i] = parVect[i] / T[i][0];
     }
 
@@ -197,9 +217,11 @@
 	invNdiag[i] = vect[0] * vect[0];
 	for (j = i + 1; j < n; j++) {
 	    somma = 0.0;
-	    for (k = i; k < j; k++) {
-		if ((j - k) < BW)
-		    somma -= vect[k - i] * T[k][j - k];
+	    /* start = i or j - BW + 1 */
+	    start = ((j - BW + 1) < i ? i : (j - BW + 1));
+	    /* end = j */
+	    for (k = start; k < j; k++) {
+		somma -= vect[k - i] * T[k][j - k];
 	    }
 	    vect[j - i] = somma * T[j][0];
 	    invNdiag[i] += vect[j - i] * vect[j - i];

Modified: grass/branches/releasebranch_6_4/vector/lidar/lidarlib/raster.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/lidarlib/raster.c	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/lidarlib/raster.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -21,16 +21,15 @@
 {
     int i;
     char buf[1024];
+    dbString sql;
 
     double interpolation, csi, eta, weight;
     struct line_pnts *point;
-    dbString sql;
 
-
-    /*CREARE LA TABELLA */
-
     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 */
@@ -52,10 +51,8 @@
 
 	    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);
@@ -67,12 +64,11 @@
 			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);
@@ -81,14 +77,13 @@
 			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 access 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);
@@ -97,13 +92,12 @@
 			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 access 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);
@@ -111,16 +105,15 @@
 			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 access 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);
@@ -129,13 +122,12 @@
 			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 access 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;
 
@@ -145,13 +137,12 @@
 			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 access 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);
@@ -159,15 +150,14 @@
 			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 access 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);
@@ -175,11 +165,11 @@
 			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 access 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);
@@ -188,14 +178,17 @@
 			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 access table <%s>"),
 					  buf);
 		    }
 		}
 	    }
-	}
-     /*IF*/}
-     /*FOR*/ return;
+	}  /*IF*/
+    }  /*FOR*/
+
+    db_commit_transaction(driver);
+
+    return;
 }
 
 
@@ -207,13 +200,36 @@
 			  int ncols, int bilin)
 {
 
-    int col, row;
+    int col, row, startcol, endcol, startrow, endrow;
     double X, Y, interpolation, weight, csi, eta;
     struct Cell_head Original;
 
     G_get_window(&Original);
-    for (row = 0; row < nrows; row++) {
-	for (col = 0; col < ncols; col++) {
+    if (Original.north > General.N)
+	startrow = (Original.north - General.N) / Original.ns_res -1;
+    else
+	startrow = 0;
+    if (Original.north > General.S) {
+	endrow = (Original.north - General.S) / Original.ns_res + 1;
+	if (endrow > nrows)
+	    endrow = nrows;
+    }
+    else
+	endrow = nrows;
+    if (General.W > Original.west)
+	startcol = (General.W - Original.west) / Original.ew_res - 1;
+    else
+	startcol = 0;
+    if (General.E > Original.west) {
+	endcol = (General.E - Original.west) / Original.ew_res + 1;
+	if (endcol > ncols)
+	    endcol = ncols;
+    }
+    else
+	endcol = ncols;
+
+    for (row = startrow; row < endrow; row++) {
+	for (col = startcol; col < endcol; col++) {
 	    X = G_col_to_easting((double)(col) + 0.5, &Original);
 	    Y = G_row_to_northing((double)(row) + 0.5, &Original);
 
@@ -234,70 +250,59 @@
 
 		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;
 			}
 		    }

Modified: grass/branches/releasebranch_6_4/vector/lidar/lidarlib/zones.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/lidarlib/zones.c	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/lidarlib/zones.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -22,11 +22,44 @@
 }
 
 /*----------------------------------------------------------------------------------------*/
+/*
+ --------------------------------------------
+ |            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, BOUND_BOX * General,
 	      BOUND_BOX * Overlap, struct Reg_dimens dim, int type)
 {
-    /* Set the Elaborationoration 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;
 
@@ -53,34 +86,34 @@
 	Overlap->E = General->E - dim.overlap;
 	return 0;
 
-    case FIRST_ROW:		/* It is just started with first row */
-	Elaboration->north = orig.north;
+    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;
+	General->N = Elaboration->north - 2 * dim.orlo_h;
 	General->S = Elaboration->south + dim.orlo_h;
-	Overlap->N = Elaboration->north;
+	Overlap->N = General->N;
 	Overlap->S = General->S + dim.overlap;
 	return 0;
 
-    case LAST_ROW:		/* It is reached last row */
-	Elaboration->south = orig.south;
-	Overlap->S = Elaboration->south;
-	General->S = Elaboration->south;
+    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 */
-	Elaboration->west = orig.west;
+    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;
+	General->W = Elaboration->west + 2 * dim.orlo_v;
 	General->E = Elaboration->east - dim.orlo_v;
-	Overlap->W = Elaboration->west;
+	Overlap->W = General->W;
 	Overlap->E = General->E - dim.overlap;
 	return 0;
 
-    case LAST_COLUMN:		/* It is reached last column */
-	Elaboration->east = orig.east;
-	Overlap->E = Elaboration->east;
-	General->E = Elaboration->east;
+    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;
 	return 0;
     }
 
@@ -88,22 +121,104 @@
 }
 
 /*----------------------------------------------------------------------------------------*/
+int P_set_dim(struct Reg_dimens *dim, double pe, double pn, int *nsplx, int *nsply)
+{
+    int total_splines, orlo_splines, n_windows, minsplines;
+    double E_extension, N_extension, orloE, orloN;
+    struct Cell_head orig;
+    int ret = 0;
+
+    G_get_window(&orig);
+
+    E_extension = orig.east - orig.west;
+    N_extension = orig.north - orig.south;
+    dim->latoE = *nsplx * pe;
+    dim->latoN = *nsply * pn;
+    orloE = dim->latoE - dim->overlap - 2 * dim->orlo_v;
+    orloN = dim->latoN - dim->overlap - 2 * dim->orlo_h;
+
+    /* number of moving windows: E_extension / orloE */
+    /* remaining steps: total steps - (floor(E_extension / orloE) * E_extension) / passoE */
+    /* remaining steps must be larger than orlo_v + overlap + half of overlap window */
+    total_splines = ceil(E_extension / pe);
+    orlo_splines = orloE / pe;
+    n_windows = floor(E_extension / orloE); /* without last one */
+    if (n_windows > 0) {
+	minsplines = ceil((double)(dim->latoE / 2.0 - dim->orlo_v - dim->overlap) / pe);
+	while (total_splines - orlo_splines * n_windows < minsplines) {
+	    *nsplx -= 1;
+	    dim->latoE = *nsplx * pe;
+	    orloE = dim->latoE - dim->overlap - 2 * dim->orlo_v;
+
+	    orlo_splines = orloE / pe;
+	    n_windows = floor(E_extension / orloE); /* without last one */
+	    minsplines = ceil((double)(dim->latoE / 2.0 - dim->orlo_v - dim->overlap) / pe);
+	    if (ret == 0)
+		ret = 1;
+	}
+	while (total_splines - orlo_splines * n_windows < minsplines * 2 && minsplines > 30) {
+	    *nsplx -= 1;
+	    dim->latoE = *nsplx * pe;
+	    orloE = dim->latoE - dim->overlap - 2 * dim->orlo_v;
+
+	    orlo_splines = orloE / pe;
+	    n_windows = floor(E_extension / orloE); /* without last one */
+	    minsplines = ceil((double)(dim->latoE / 2.0 - dim->orlo_v - dim->overlap) / pe);
+	    if (ret == 0)
+		ret = 1;
+	}
+    }
+
+    total_splines = ceil(N_extension / pn);
+    orlo_splines = orloN / pn;
+    n_windows = floor(N_extension / orloN); /* without last one */
+    if (n_windows > 0) {
+	minsplines = ceil((double)(dim->latoN / 2.0 - dim->orlo_h - dim->overlap) / pn);
+	while (total_splines - orlo_splines * n_windows < minsplines) {
+	    *nsply -= 1;
+	    dim->latoN = *nsply * pn;
+	    orloN = dim->latoN - dim->overlap - 2 * dim->orlo_h;
+
+	    orlo_splines = orloN / pn;
+	    n_windows = floor(N_extension / orloN); /* without last one */
+	    minsplines = ceil((double)(dim->latoN / 2.0 - dim->orlo_h - dim->overlap) / pn);
+	    if (ret < 2)
+		ret += 2;
+	}
+	while (total_splines - orlo_splines * n_windows < minsplines * 2 && minsplines > 30) {
+	    *nsply -= 1;
+	    dim->latoN = *nsply * pn;
+	    orloN = dim->latoN - dim->overlap - 2 * dim->orlo_h;
+
+	    orlo_splines = orloN / pn;
+	    n_windows = floor(N_extension / orloN); /* without last one */
+	    minsplines = ceil((double)(dim->latoN / 2.0 - dim->orlo_h - dim->overlap) / pn);
+	    if (ret < 2)
+		ret += 2;
+	}
+    }
+    return 0;
+}
+
+/*----------------------------------------------------------------------------------------*/
 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 = 30 * pe;	/*4 */
-	dim->orlo_h = 30 * pn;
+    /* 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 */
-	dim->orlo_v = 40 * pe;	/*3 */
-	dim->orlo_h = 40 * pn;
+    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!! */
 }
 
 /*----------------------------------------------------------------------------------------*/
@@ -147,12 +262,73 @@
 
     return mean;
 }
+
+double P_estimate_splinestep(struct Map_info *Map, double *dens, double *dist)
+{
+    int type, npoints = 0;
+    double xmin = 0, xmax = 0, ymin = 0, ymax = 0;
+    double x, y, z;
+    struct line_pnts *points;
+    struct line_cats *categories;
+    BOUND_BOX region_box;
+    struct Cell_head orig;
+
+    G_get_set_window(&orig);
+    Vect_region_box(&orig, &region_box);
+
+    points = Vect_new_line_struct();
+    categories = Vect_new_cats_struct();
+
+    Vect_rewind(Map);
+    while ((type = Vect_read_next_line(Map, points, categories)) > 0) {
+	if (!(type & GV_POINT))
+	    continue;
+
+	x = points->x[0];
+	y = points->y[0];
+	if (points->z != NULL)
+	    z = points->z[0];
+	else
+	    z = 0.0;
+
+	/* only use points in current region */
+	if (Vect_point_in_box(x, y, z, &region_box)) {
+	    npoints++;
+
+	    if (npoints > 1) {
+		if (xmin > x)
+		    xmin = x;
+		else if (xmax < x)
+		    xmax = x;
+		if (ymin > y)
+		    ymin = y;
+		else if (ymax < y)
+		    ymax = y;
+	    }
+	    else {
+		xmin = xmax = x;
+		ymin = ymax = y;
+	    }
+	}
+    }
+    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;
+    }
+    else {
+	return -1;
+    }
+}
+
 struct Point *P_Read_Vector_Region_Map(struct Map_info *Map,
 				       struct Cell_head *Elaboration,
 				       int *num_points, int dim_vect,
 				       int layer)
 {
-    int line_num, pippo, npoints, cat;
+    int line_num, pippo, npoints, cat, type;
     double x, y, z;
     struct Point *obs;
     struct line_pnts *points;
@@ -165,15 +341,18 @@
     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;
     line_num = 0;
 
     Vect_rewind(Map);
-    while (Vect_read_next_line(Map, points, categories) > 0) {
+    while ((type = Vect_read_next_line(Map, points, categories)) > 0) {
 
+	if (!(type & GV_POINT))
+	    continue;
+
 	line_num++;
 
 	x = points->x[0];
@@ -183,7 +362,7 @@
 	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) {
@@ -212,42 +391,69 @@
 }
 
 /*------------------------------------------------------------------------------------------------*/
-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;
 }
 
 /*------------------------------------------------------------------------------------------------*/
@@ -269,7 +475,6 @@
     struct Cell_head Original;
 
     G_get_window(&Original);
-    G_set_window(&Original);
     nrows = G_window_rows();
     ncols = G_window_cols();
 
@@ -286,6 +491,7 @@
 	}
 	G_put_d_raster_row(fd, raster);
     }
+    G_percent(row, nrows, 2);
 }
 
 /*------------------------------------------------------------------------------------------------*/
@@ -294,7 +500,7 @@
 		char *tab_name)
 {
 
-    int more, ltype, line_num, type, count = 0;
+    int more, line_num, type, count = 0;
     double coordX, coordY, coordZ;
 
     struct line_pnts *point;
@@ -337,7 +543,7 @@
 	    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));
@@ -345,7 +551,7 @@
 	    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));
@@ -353,7 +559,7 @@
 	    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);
@@ -372,7 +578,6 @@
 
 
     G_get_window(&Original);
-    G_set_window(&Original);
     nrows = G_window_rows();
     ncols = G_window_cols();
 
@@ -384,10 +589,30 @@
 }
 #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|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
+ */

Modified: grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/correction.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/correction.c	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/correction.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -20,6 +20,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/Vect.h>
 #include <grass/dbmi.h>
@@ -29,183 +30,198 @@
 void
 P_Sparse_Correction(struct Map_info *In, struct Map_info *Out,
 		    struct Map_info *Terrain, struct Cell_head *Elaboration,
-		    BOUND_BOX General, 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,
-		    dbDriver * driver, double mean)
+		    BOUND_BOX General, 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;
+    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;
 
-		Vect_cat_get(cats, F_CLASSIFICATION, &class);
+	    if (Vect_point_in_box(*point->x, *point->y, *point->z, &Overlap)) {	/*(5) */
+
+		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(Terrain, GV_POINT, point, cats);
 
-		Vect_write_line(Out, GV_POINT, points, 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
-			    (&interpolation, line_num[i], driver) != DB_OK)
+			    (&interpolation, line_num[i], driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to read the database"));
 
 			if (UpDate_Correction
-			    (interpolation, line_num[i], driver) != DB_OK)
+			    (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) != DB_OK)
+			     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],
-			     driver) != DB_OK)
+			     driver, tab_name) != DB_OK)
 			    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;
 
-			interpolation *= weight;
 			if (Select_Correction
-			    (&interpolation, line_num[i], driver) != DB_OK)
+			    (&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;
 
 			if (Select_Correction
-			    (&interpolation, line_num[i], driver) != DB_OK)
+			    (&interpolation, line_num[i], driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to read the database"));
 
 			if (UpDate_Correction
-			    (interpolation, line_num[i], driver) != DB_OK)
+			    (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) != DB_OK)
+			    (&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) != DB_OK)
+			    (&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],
-			     driver) != DB_OK)
+			     driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to write in the database"));
 		    }
 		}
 	    }
-	}
-	 /*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;
 }
 
@@ -228,7 +244,7 @@
     return class;
 }
 
-int Select_Correction(double *Interp, int line_num, dbDriver * driver)
+int Select_Correction(double *Interp, int line_num, dbDriver * driver, char *tab_name)
 {
     int more;
     char buf[1024];
@@ -239,10 +255,8 @@
     dbValue *Interp_value;
 
     db_init_string(&sql);
-    db_zero_string(&sql);
 
-    sprintf(buf, "SELECT Interp FROM Auxiliar_correction_table WHERE ID=%d",
-	    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,55 +274,61 @@
 
 	*Interp += db_get_value_double(Interp_value);
     }
+    db_close_cursor(&cursor);
+    db_free_string(&sql);
     return DB_OK;
 }
 
-int Insert_Correction(double Interp, int line_num, dbDriver * driver)
+int Insert_Correction(double Interp, int line_num, dbDriver * driver, char *tab_name)
 {
     char buf[1024];
     dbString sql;
+    int ret;
 
     db_init_string(&sql);
-    sprintf(buf, "INSERT INTO Auxiliar_correction_table (ID, Interp)");
+    sprintf(buf, "INSERT INTO %s (ID, Interp)", tab_name);
     db_append_string(&sql, buf);
     sprintf(buf, " VALUES (%d, %lf)", line_num, Interp);
     db_append_string(&sql, buf);
 
-    return db_execute_immediate(driver, &sql);
+    ret = db_execute_immediate(driver, &sql);
+    db_free_string(&sql);
+    return ret;
 }
 
-int UpDate_Correction(double Interp, int line_num, dbDriver * driver)
+int UpDate_Correction(double Interp, int line_num, dbDriver * driver, char *tab_name)
 {
     char buf[1024];
     dbString sql;
+    int ret;
 
     db_init_string(&sql);
     sprintf(buf,
-	    "UPDATE Auxiliar_correction_table SET Interp=%lf WHERE ID=%d",
-	    Interp, line_num);
+	    "UPDATE %s SET Interp=%lf WHERE ID=%d", tab_name, Interp, line_num);
     db_append_string(&sql, buf);
 
-    return db_execute_immediate(driver, &sql);
+    ret = db_execute_immediate(driver, &sql);
+    db_free_string(&sql);
+    return ret;
 }
 
 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;
     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();
@@ -338,12 +358,10 @@
 	    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 */
@@ -351,27 +369,55 @@
 	    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;
 }
 
+/*! DEFINITION OF THE SUBZONES 
 
+  5: inside Overlap region
+  all others: inside General region but outside Overlap region
 
-#ifdef notdef
-/* SUBZONES DEFINITION */
------------------------|4 | 3 | 3 | ||-----------------------|||||||2 | 5 | 1
-    | ||||||||-----------------------|2 | 1 | 1 |
-    ||-----------------------||||||||||||||||||-----------------------||||||-----------------------
-#endif
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       |4|   3   |3|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       |2|   5   |1|       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       |2|   1   |1|       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   | |       | |       | |       | |
+   ---------------------------------
+   | |       | |       | |       | |
+   ---------------------------------
+ */

Modified: grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/correction.h
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/correction.h	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/correction.h	2010-02-17 17:59:30 UTC (rev 41063)
@@ -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 @@
 		    BOUND_BOX, /**/
 		    BOUND_BOX, /**/
 		    double **, /**/
+		    struct lidar_cat *,
 		    double *, /**/
 		    int *, /**/
 		    double, /**/
@@ -23,19 +32,21 @@
 		    double, /**/
 		    double, /**/
 		    int, /**/
-		    int, /**/ int, /**/ dbDriver *, /**/ double /**/);
+		    int, /**/
+		    int, /**/
+		    dbDriver *, /**/
+		    double, /**/
+		    char *);
 
-int Insert_Correction(double, /**/ int, /**/ dbDriver * /**/);
+int Insert_Correction(double, /**/ int, /**/ dbDriver *, /**/ char *);
 
-int UpDate_Correction(double, /**/ int, /**/ dbDriver * /**/);
+int UpDate_Correction(double, /**/ int, /**/ dbDriver *, /**/ char *);
 
-int Select_Correction(double *, /**/ int, /**/ dbDriver * /**/);
+int Select_Correction(double *, /**/ int, /**/ dbDriver *, /**/ char *);
 
-int Create_AuxEdge_Table(dbDriver *);
-
 struct Point *P_Read_Vector_Correction(struct Map_info *, /**/
 				       struct Cell_head *, /**/
-				       int *, /**/ int *, /**/ int /**/);
+				       int *, /**/ int *, /**/ int /**/,
+				       struct lidar_cat **);
 
-void P_Aux_to_Raster(double **, int);
 int correction(int, double, double, double, double);

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	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.lidar.correction/main.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -20,6 +20,7 @@
  /*INCLUDES*/
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #include <grass/config.h>
 #include <grass/gis.h>
 #include <grass/Vect.h>
@@ -32,10 +33,16 @@
 int main(int argc, char *argv[])
 {
     /* Declarations */
-    int dim_vect, nparameters, BW, npoints, nrows, ncols, nsply, nsplx;
-    char *dvr, *db, *mapset, table_name[1024];
+    int dim_vect, nparameters, BW, npoints, nrows, ncols;
+    int nsply, nsplx, nsplx_adj, nsply_adj;
+    int nsubregion_col, nsubregion_row;
+    int subregion = 0, nsubregions = 0;
+    const char *dvr, *db, *mapset;
+    char table_name[GNAME_MAX];
+    char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
     double lambda, ew_resol, ns_resol, mean, passoN, passoE, HighThresh,
 	LowThresh;
+    double N_extension, E_extension, orloE, orloN;
 
     int i, nterrain, count_terrain;
 
@@ -43,11 +50,12 @@
 
     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,
 	*passoN_opt, *lambda_f_opt, *Thresh_A_opt, *Thresh_B_opt;
+    struct Flag *spline_step_flag;
     struct GModule *module;
 
     struct Cell_head elaboration_reg, original_reg;
@@ -55,6 +63,7 @@
     BOUND_BOX general_box, overlap_box;
 
     struct Point *observ;
+    struct lidar_cat *lcat;
 
     dbDriver *driver;
 
@@ -65,6 +74,12 @@
     module->description =
 	_("Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.");
 
+    spline_step_flag = G_define_flag();
+    spline_step_flag->key = 'e';
+    spline_step_flag->label = _("Estimate point density and distance");
+    spline_step_flag->description =
+	_("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 =
 	_("Input observation vector map name (v.lidar.growing output)");
@@ -132,12 +147,32 @@
     lambda = atof(lambda_f_opt->answer);
     HighThresh = atof(Thresh_A_opt->answer);
     LowThresh = atof(Thresh_B_opt->answer);
-    dvr = G__getenv2("DB_DRIVER", G_VAR_MAPSET);
-    db = G__getenv2("DB_DATABASE", G_VAR_MAPSET);
 
+    if (!(db = G__getenv2("DB_DATABASE", G_VAR_MAPSET)))
+	G_fatal_error(_("Unable to read name of database"));
+
+    if (!(dvr = G__getenv2("DB_DRIVER", G_VAR_MAPSET)))
+	G_fatal_error(_("Unable to read name of driver"));
+
     /* Setting auxiliar table's name */
-    sprintf(table_name, "%s_aux", out_opt->answer);
+    if (G__name_is_fully_qualified(out_opt->answer, xname, xmapset)) {
+	sprintf(table_name, "%s_aux", xname);
+    }
+    else
+	sprintf(table_name, "%s_aux", out_opt->answer);
 
+    /* Something went wrong in a previous v.lidar.correction execution */
+    if (db_table_exists(dvr, db, table_name)) {
+	/* Start driver and open db */
+	driver = db_start_driver_open_database(dvr, db);
+	if (driver == NULL)
+	    G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
+			  dvr);
+	if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
+	    G_fatal_error(_("Old auxiliar table could not be dropped"));
+	db_close_database_shutdown_driver(driver);
+    }
+
     /* Checking vector names */
     Vect_check_input_output_name(in_opt->answer, out_opt->answer,
 				 GV_FATAL_EXIT);
@@ -146,10 +181,28 @@
     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;
+	if (P_estimate_splinestep(&In, &dens, &dist) == 0) {
+	    G_message("Estimated point density: %.4g", dens);
+	    G_message("Estimated mean distance between points: %.4g", dist);
+	}
+	else
+	    G_warning(_("No points in current region!"));
+	
+	Vect_close(&In);
+	exit(EXIT_SUCCESS);
+    }
+
     /* Open output vector */
     if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z)) {
 	Vect_close(&In);
@@ -176,32 +229,75 @@
 	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) {
+	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);
     Vect_region_box(&elaboration_reg, &overlap_box);
     Vect_region_box(&elaboration_reg, &general_box);
 
-    /* Fixxing parameters of the elaboration region */
-    /*! Each original_region will be divided into several subregions. These
-     *  subregion will be overlaped by its neibourgh subregions. This overlapping
-     *  is calculated as OVERLAP_PASS times the east-west resolution. */
-
     nrows = G_window_rows();
     ncols = G_window_cols();
 
     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);
-    dims.latoE = NSPLX_MAX * passoE;
-    dims.latoN = NSPLY_MAX * passoN;
-    dims.overlap = OVERLAP_SIZE * ew_resol;
+
+    nsplx_adj = NSPLX_MAX;
+    nsply_adj = NSPLY_MAX;
+    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);
 
-    /* Subdividing and working with tiles */
+    G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
+    G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
+
+    /* calculate number of subregions */
+    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;
+
+    if (nsubregion_col < 0)
+	nsubregion_col = 0;
+    if (nsubregion_row < 0)
+	nsubregion_row = 0;
+
+    nsubregions = nsubregion_row * nsubregion_col;
+
     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,
@@ -219,11 +315,13 @@
 	}
 
 	nsply =
-	    ceil((elaboration_reg.north - elaboration_reg.south) / passoN) +
-	    1;
+	    ceil((elaboration_reg.north -
+		  elaboration_reg.south) / passoN) + 0.5;
+	/*
 	if (nsply > NSPLY_MAX) {
 	    nsply = NSPLY_MAX;
 	}
+	*/
 	G_debug(1, _("nsply = %d"), nsply);
 
 	elaboration_reg.east = original_reg.west;
@@ -231,6 +329,10 @@
 
 	while (last_column == FALSE) {	/* For each column */
 
+	    subregion++;
+	    if (nsubregions > 1)
+		G_message(_("subregion %d of %d"), subregion, nsubregions);
+
 	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 			  GENERAL_COLUMN);
 
@@ -247,33 +349,36 @@
 
 	    nsplx =
 		ceil((elaboration_reg.east - elaboration_reg.west) / passoE) +
-		1;
+		0.5;
+	    /*
 	    if (nsplx > NSPLX_MAX) {
 		nsplx = NSPLX_MAX;
 	    }
+	    */
 	    G_debug(1, _("nsplx = %d"), nsplx);
 
 	    dim_vect = nsplx * nsply;
 	    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);
 
@@ -288,9 +393,14 @@
 			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,
@@ -301,30 +411,25 @@
 		G_free_matrix(N);
 		G_free_vector(TN);
 		G_free_vector(Q);
+		G_free_matrix(obsVect);
 
-		if (flag_auxiliar == FALSE) {
-		    if ((flag_auxiliar =
-			 P_Create_Aux_Table(driver, table_name)) == FALSE) {
-			Vect_close(&In);
-			Vect_close(&Out);
-			Vect_close(&Terrain);
-			exit(EXIT_FAILURE);
-		    }
-		}
-
-		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);
+				    nsplx, nsply, npoints, driver, mean, table_name);
 
-		G_free_matrix(obsVect);
+		G_free_vector(parVect);
+		G_free_matrix(obsVect_all);
 		G_free_ivector(lineVect);
-
-		G_free_vector(parVect);
 	    }
-	    G_free(observ);
+	    else {
+		G_free(observ);
+		G_warning(_("No data within this subregion. "
+			    "Consider changing the spline step."));
+	    }
+	    G_free(lcat);
 	}			/*! END WHILE; last_column = TRUE */
     }				/*! END WHILE; last_row = TRUE */
 

Modified: grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/edgedetection.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/edgedetection.c	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/edgedetection.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -26,11 +26,13 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/Vect.h>
 #include <grass/dbmi.h>
     /* #include <grass/PolimiFunct.h> */
 #include "edgedetection.h"
+
 int edge_detection(struct Cell_head elaboration_reg, BOUND_BOX Overlap_Box,
 		   double *parBilin, double obsX, double obsY,
 		   double *partial, double alpha, double residual,
@@ -41,15 +43,15 @@
     /* 3 = PRE_UNKNOWN */
 
     int c1, c2;
-    double g[9][2], *gradient, gradPto, dirPto;
+    double g[9][2], gradient[2], gradPto, dirPto;
     extern double passoE, passoN;
     static struct Cell_head Elaboration;
 
     g[0][0] = partial[0];
     g[0][1] = partial[1];
 
-    gradPto = sqrt(g[0][0] * g[0][0] + g[0][1] * g[0][1]);
-    dirPto = atan(g[0][1] / g[0][0]) + PI / 2;	/* radiants */
+    gradPto = g[0][0] * g[0][0] + g[0][1] * g[0][1];
+    dirPto = atan(g[0][1] / g[0][0]) + M_PI / 2;	/* radiants */
 
     Elaboration = elaboration_reg;
 
@@ -59,72 +61,64 @@
     else if ((gradPto > gradLow) && (residual > 0)) {	/* Soft condition for 'edge' points */
 
 	if (Vect_point_in_box(obsX, obsY, 0.0, &Overlap_Box)) {
-	    gradient =
-		Get_Gradient(Elaboration, obsX + passoE * cos(dirPto),
-			     obsY + passoN * sin(dirPto), parBilin);
+	    Get_Gradient(Elaboration, obsX + passoE * cos(dirPto),
+			     obsY + passoN * sin(dirPto), parBilin, gradient);
 	    g[2][0] = gradient[0];
 	    g[2][1] = gradient[1];
 
-	    gradient =
-		Get_Gradient(Elaboration, obsX + passoE * cos(dirPto + PI),
-			     obsY + passoN * sin(dirPto + PI), parBilin);
+	    Get_Gradient(Elaboration, obsX + passoE * cos(dirPto + M_PI),
+			     obsY + passoN * sin(dirPto + M_PI), parBilin, gradient);
 	    g[7][0] = gradient[0];
 	    g[7][1] = gradient[1];
 
-	    if ((fabs(atan(g[2][1] / g[2][0]) + PI / 2 - dirPto) < alpha) &&
-		(fabs(atan(g[7][1] / g[7][0]) + PI / 2 - dirPto) < alpha)) {
+	    if ((fabs(atan(g[2][1] / g[2][0]) + M_PI / 2 - dirPto) < alpha) &&
+		(fabs(atan(g[7][1] / g[7][0]) + M_PI / 2 - dirPto) < alpha)) {
 
-		gradient =
-		    Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto + PI / 4),
-				 obsY + passoN * sin(dirPto + PI / 4),
-				 parBilin);
+		Get_Gradient(Elaboration,
+				 obsX + passoE * cos(dirPto + M_PI / 4),
+				 obsY + passoN * sin(dirPto + M_PI / 4),
+				 parBilin, gradient);
 		g[1][0] = gradient[0];
 		g[1][1] = gradient[1];
 
-		gradient =
-		    Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto - PI / 4),
-				 obsY + passoN * sin(dirPto - PI / 4),
-				 parBilin);
+		Get_Gradient(Elaboration,
+				 obsX + passoE * cos(dirPto - M_PI / 4),
+				 obsY + passoN * sin(dirPto - M_PI / 4),
+				 parBilin, gradient);
 		g[3][0] = gradient[0];
 		g[3][1] = gradient[1];
 
-		gradient =
-		    Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto + PI / 2),
-				 obsY + passoN * sin(dirPto + PI / 2),
-				 parBilin);
+		Get_Gradient(Elaboration,
+				 obsX + passoE * cos(dirPto + M_PI / 2),
+				 obsY + passoN * sin(dirPto + M_PI / 2),
+				 parBilin, gradient);
 		g[4][0] = gradient[0];
 		g[4][1] = gradient[1];
 
-		gradient =
-		    Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto - PI / 2),
-				 obsY + passoN * sin(dirPto - PI / 2),
-				 parBilin);
+		Get_Gradient(Elaboration,
+				 obsX + passoE * cos(dirPto - M_PI / 2),
+				 obsY + passoN * sin(dirPto - M_PI / 2),
+				 parBilin, gradient);
 		g[5][0] = gradient[0];
 		g[5][1] = gradient[1];
 
-		gradient =
-		    Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto + PI * 3 / 4),
-				 obsY + passoN * sin(dirPto + PI * 3 / 4),
-				 parBilin);
+		Get_Gradient(Elaboration,
+				 obsX + passoE * cos(dirPto + M_PI * 3 / 4),
+				 obsY + passoN * sin(dirPto + M_PI * 3 / 4),
+				 parBilin, gradient);
 		g[6][0] = gradient[0];
 		g[6][1] = gradient[1];
 
-		gradient =
-		    Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto - PI * 3 / 4),
-				 obsY + passoN * sin(dirPto - PI * 3 / 4),
-				 parBilin);
+		Get_Gradient(Elaboration,
+				 obsX + passoE * cos(dirPto - M_PI * 3 / 4),
+				 obsY + passoN * sin(dirPto - M_PI * 3 / 4),
+				 parBilin, gradient);
 		g[8][0] = gradient[0];
 		g[8][1] = gradient[1];
 
 		c2 = 0;
 		for (c1 = 0; c1 < 9; c1++)
-		    if (sqrt(g[c1][0] * g[c1][0] + g[c1][1] * g[c1][1]) >
+		    if (g[c1][0] * g[c1][0] + g[c1][1] * g[c1][1] >
 			gradHigh)
 			c2++;
 
@@ -137,23 +131,21 @@
 		return PRE_TERRAIN;
 	}
 	else
-	    return UNKNOWN;
+	    return PRE_UNKNOWN;
     }				/* END ELSE IF */
     else
 	return PRE_TERRAIN;
 }
 
-double *Get_Gradient(struct Cell_head Elaboration, double X, double Y,
-		     double *parVect)
+int Get_Gradient(struct Cell_head Elaboration, double X, double Y,
+		     double *parVect, double *grad)
 {
     int row, col, N;
-    double csi, eta, d, b, a, c, *grad;
+    double csi, eta, d, b, a, c;
 
     extern int nsply;
     extern double passoN, passoE;
 
-    grad = (double *)G_calloc(2, sizeof(double));
-
     row = (int)((Y - Elaboration.south) / passoN);
     col = (int)((X - Elaboration.west) / passoE);
     N = nsply * col + row;
@@ -165,7 +157,7 @@
     c = parVect[N + 1 + nsply] - a - b - d;
     grad[0] = (a + c * eta);
     grad[1] = (b + c * csi);
-    return grad;
+    return 0;
 }
 
 void classification(struct Map_info *Out, struct Cell_head Elaboration,
@@ -173,10 +165,10 @@
 		    double *parBilin, double *parBicub, double mean,
 		    double alpha, double gradHigh, double gradLow,
 		    double overlap, int *line_num, int num_points,
-		    dbDriver * driver, char *vect_name)
+		    dbDriver * driver, char *tabint_name, char *tab_name)
 {
     int i, edge;
-    double interpolation, weight, residual, eta, csi, *gradient;
+    double interpolation, weight, residual, eta, csi, gradient[2];
 
     extern int nsplx, nsply, line_out_counter;
     extern double passoN, passoE;
@@ -187,7 +179,10 @@
     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);
 	Vect_reset_cats(categories);
@@ -199,15 +194,14 @@
 				       Elaboration.south, parBicub);
 	    interpolation += mean;
 
-	    Vect_copy_xyz_to_pnts(point, &obs[i][0], &obs[i][1], &obs[i][2],
-				  1);
-	    gradient =
-		Get_Gradient(Elaboration, obs[i][0], obs[i][1], parBilin);
+	    Vect_copy_xyz_to_pnts(point, &obs[i][0], &obs[i][1], &obs[i][2], 1);
+
+	    Get_Gradient(Elaboration, obs[i][0], obs[i][1], parBilin, gradient);
+
 	    *point->z += mean;
 	    /*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,
@@ -218,71 +212,70 @@
 		Vect_cat_set(categories, F_INTERPOLATION, line_out_counter);
 		Vect_write_line(Out, GV_POINT, point, categories);
 		Insert_Interpolation(interpolation, line_out_counter, driver,
-				     vect_name);
+				     tabint_name);
 		line_out_counter++;
-
 	    }
 	    else {
-		if ((*point->x > Overlap.E) && (*point->x != General.E)) {
+		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;
 
-		    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);
-
 			gradient[0] *= weight;
 			gradient[1] *= weight;
 			interpolation *= weight;
 
-			if (Select
-			    (&gradient[0], &gradient[1], &interpolation,
-			     line_num[i], driver) != DB_OK)
-			    G_fatal_error(_("Impossible to read the database"));
+			if (Select(&gradient[0], &gradient[1],
+			           &interpolation, line_num[i], driver,
+				   tab_name) != DB_OK)
+			    G_fatal_error(_("Impossible to read from aux table"));
 
-			if (UpDate
-			    (gradient[0], gradient[1], interpolation,
-			     line_num[i], driver) != DB_OK)
-			    G_fatal_error(_("Impossible to update the database"));
-
+			if (UpDate(gradient[0], gradient[1],
+			           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;
 
-			if (Insert
-			    (gradient[0] * weight, gradient[1] * weight,
-			     interpolation * weight, line_num[i],
-			     driver) != DB_OK)
-			    G_fatal_error(_("Impossible to write in the database"));
+			gradient[0] *= weight;
+			gradient[1] *= weight;
+			interpolation *= weight;
 
+			if (Insert(gradient[0], gradient[1], interpolation,
+			           line_num[i], driver, tab_name) != DB_OK)
+			    G_fatal_error(_("Impossible to write to aux table"));
+
 		    }
 		    else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) {	/*(1) */
-			weight = (*point->x - Overlap.E) / overlap;
+			weight = (General.E - *point->x) / overlap;
 
-			if (Insert
-			    (gradient[0] * weight, gradient[1] * weight,
-			     interpolation * weight, line_num[i],
-			     driver) != DB_OK)
-			    G_fatal_error(_("Impossible to write in the database"));
-		    }
+			gradient[0] *= weight;
+			gradient[1] *= weight;
+			interpolation *= weight;
 
+			if (Insert(gradient[0], gradient[1], interpolation,
+			           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;
 			interpolation *= weight;
 
-			if (Select
-			    (&gradient[0], &gradient[1], &interpolation,
-			     line_num[i], driver) != DB_OK)
-			    G_fatal_error(_("Impossible to read the database"));
+			if (Select(&gradient[0], &gradient[1],
+			           &interpolation, line_num[i], driver,
+			           tab_name) != DB_OK)
+			    G_fatal_error(_("Impossible to read from aux table"));
 
 			residual = *point->z - interpolation;
 			edge =
@@ -297,11 +290,10 @@
 				     line_out_counter);
 			Vect_write_line(Out, GV_POINT, point, categories);
 			Insert_Interpolation(interpolation, line_out_counter,
-					     driver, vect_name);
+					     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;
@@ -310,28 +302,28 @@
 			gradient[1] *= weight;
 			interpolation *= weight;
 
-			if (Select
-			    (&gradient[0], &gradient[1], &interpolation,
-			     line_num[i], driver) != DB_OK)
-			    G_fatal_error(_("Impossible to read the database"));
+			if (Select(&gradient[0], &gradient[1],
+			           &interpolation, line_num[i], driver,
+			           tab_name) != DB_OK)
+			    G_fatal_error(_("Impossible to read from aux table"));
 
-			if (UpDate
-			    (gradient[0], gradient[1], interpolation,
-			     line_num[i], driver) != DB_OK)
-			    G_fatal_error(_("Impossible to update the database"));
+			if (UpDate(gradient[0], gradient[1],
+			           interpolation, line_num[i], driver,
+			    tab_name) != DB_OK)
+			    G_fatal_error(_("Impossible to update aux table"));
 
 		    }
 		    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;
 			interpolation *= weight;
 
-			if (Select
-			    (&gradient[0], &gradient[1], &interpolation,
-			     line_num[i], driver) != DB_OK)
-			    G_fatal_error(_("Impossible to read the database"));
+			if (Select(&gradient[0], &gradient[1],
+			           &interpolation, line_num[i], driver,
+			           tab_name) != DB_OK)
+			    G_fatal_error(_("Impossible to read from aux table"));
 
 			residual = *point->z - interpolation;
 			edge =
@@ -346,23 +338,22 @@
 				     line_out_counter);
 			Vect_write_line(Out, GV_POINT, point, categories);
 			Insert_Interpolation(interpolation, line_out_counter,
-					     driver, vect_name);
+					     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;
 			interpolation *= weight;
 
-			if (Select
-			    (&gradient[0], &gradient[1], &interpolation,
-			     line_num[i], driver) != DB_OK)
-			    G_fatal_error(_("Impossible to read the database"));
+			if (Select(&gradient[0], &gradient[1],
+			           &interpolation, line_num[i], driver,
+			           tab_name) != DB_OK)
+			    G_fatal_error(_("Impossible to read from aux table"));
 
 			residual = *point->z - interpolation;
 			edge =
@@ -377,82 +368,96 @@
 				     line_out_counter);
 			Vect_write_line(Out, GV_POINT, point, categories);
 			Insert_Interpolation(interpolation, line_out_counter,
-					     driver, vect_name);
+					     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;
 
-			if (Insert
-			    (gradient[0] * weight, gradient[1] * weight,
-			     interpolation * weight, line_num[i], driver)
-			    != DB_OK)
-			    G_fatal_error(_("Impossible to write in the database"));
+			gradient[0] *= weight;
+			gradient[1] *= weight;
+			interpolation *= weight;
 
+			if (Insert(gradient[0], gradient[1], interpolation,
+			           line_num[i], driver, tab_name) != DB_OK)
+			    G_fatal_error(_("Impossible to write to aux table"));
+
 		    }		/*else (1) */
 		}		/*else */
 	    }
 	}			/*end if obs */
     }				/*end for */
+    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)
+int Insert(double partialX, double partialY, double Interp,
+           int line_num, dbDriver * driver, char *tab_name)
 {
-
     char buf[1024];
     dbString sql;
+    int ret;
 
     db_init_string(&sql);
     sprintf(buf,
-	    "INSERT INTO Auxiliar_edge_table (ID, Interp, partialX, partialY)");
+	    "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);
+    sprintf(buf, " VALUES (%d, %lf, %lf, %lf)", line_num, Interp,
+            partialX, partialY);
     db_append_string(&sql, buf);
 
-    return db_execute_immediate(driver, &sql);
+    ret = db_execute_immediate(driver, &sql);
+    db_free_string(&sql);
+
+    return ret;
 }
 
 int Insert_Interpolation(double Interp, int line_num, dbDriver * driver,
-			 char *name)
+			 char *tab_name)
 {
-
     char buf[1024];
     dbString sql;
+    int ret;
 
     db_init_string(&sql);
-    sprintf(buf, "INSERT INTO %s_edge_Interpolation (ID, Interp)", name);
+    sprintf(buf, "INSERT INTO %s (ID, Interp)", tab_name);
     db_append_string(&sql, buf);
     sprintf(buf, " VALUES (%d, %lf)", line_num, Interp);
     db_append_string(&sql, buf);
 
-    return db_execute_immediate(driver, &sql);
+    ret = db_execute_immediate(driver, &sql);
+    db_free_string(&sql);
+
+    return ret;
 }
 
 
-int UpDate(double partialX, double partialY, double Interp, int line_num,
-	   dbDriver * driver)
+int UpDate(double partialX, double partialY, double Interp,
+	   int line_num, dbDriver * driver, char *tab_name)
 {
-
     char buf[1024];
     dbString sql;
+    int ret;
 
     db_init_string(&sql);
     sprintf(buf,
-	    "UPDATE Auxiliar_edge_table SET Interp=%lf, PartialX=%lf, PartialY=%lf WHERE ID=%d",
-	    Interp, partialX, partialY, line_num);
+	    "UPDATE %s SET Interp=%lf, X=%lf, Y=%lf WHERE ID=%d",
+	    tab_name, Interp, partialX, partialY, line_num);
     db_append_string(&sql, buf);
 
-    return db_execute_immediate(driver, &sql);
+    ret = db_execute_immediate(driver, &sql);
+    db_free_string(&sql);
+
+    return ret;
 }
 
-int Select(double *PartialX, double *PartialY, double *Interp, int line_num,
-	   dbDriver * driver)
+int Select(double *PartialX, double *PartialY, double *Interp,
+           int line_num, dbDriver * driver, char *tab_name)
 {
-
     int more;
     char buf[1024];
     dbString sql;
@@ -463,15 +468,16 @@
 
     db_init_string(&sql);
     sprintf(buf,
-	    "SELECT ID, Interp, partialX, partialY FROM Auxiliar_edge_table WHERE ID=%d",
-	    line_num);
+	    "SELECT ID, Interp, X, Y 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)
 	return -1;
 
+    table = db_get_cursor_table(&cursor);
+
     while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
-	table = db_get_cursor_table(&cursor);
 
 	Interp_col = db_get_table_column(table, 1);
 	PartialX_col = db_get_table_column(table, 2);
@@ -500,107 +506,34 @@
 	*PartialY += db_get_value_double(PartialY_value);
     }
     db_close_cursor(&cursor);
+    db_free_string(&sql);
     return DB_OK;
 }
 
-
-int Create_AuxEdge_Table(dbDriver * driver)
-{
-    dbTable *table;
-    dbColumn *ID_col, *Interp_col, *PartialX_col, *PartialY_col;
-    int created;
-
-    table = db_alloc_table(4);
-    db_set_table_name(table, "Auxiliar_edge_table");
-    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(3, _("<Auxiliar_edge_table> created in database."));
-	created = TRUE;
-    }
-    else
-	return FALSE;
-
-    return created;
-}
-
-int Drop_Aux_Table(dbDriver * driver)
-{
-    dbString drop;
-
-    db_init_string(&drop);
-    db_append_string(&drop, "drop table ");
-    db_append_string(&drop, "Auxiliar_edge_table");
-    return db_execute_immediate(driver, &drop);
-
-}
-
-
-int Create_Interpolation_Table(char *vect_name, dbDriver * driver)
-{
-    char table_name[1024];
-
-    dbTable *table;
-    dbColumn *ID_col, *Interp_col;
-
-    sprintf(table_name, "%s_edge_Interpolation", vect_name);
-
-    table = db_alloc_table(2);
-    db_set_table_name(table, table_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(3, _("<%s> created in database."), db_get_table_name(table));
-	return DB_OK;
-    }
-    else
-	return !DB_OK;
-}
-
-
-#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

Modified: grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/edgedetection.h
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/edgedetection.h	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/edgedetection.h	2010-02-17 17:59:30 UTC (rev 41063)
@@ -28,7 +28,6 @@
 #include <grass/dbmi.h>
 #include <grass/glocale.h>
 #include <grass/PolimiFunct.h>
-#define PI 3.141592
 
 /*---------------------------------------------------------------------------------------*/
 int edge_detection(struct Cell_head, /**/
@@ -39,8 +38,8 @@
 		   double *, /**/
 		   double, /**/ double, /**/ double, /**/ double /**/);
 
-double *Get_Gradient(struct Cell_head, /**/
-		     double, /**/ double, /**/ double * /**/);
+int Get_Gradient(struct Cell_head, /**/
+		     double, /**/ double, /**/ double *, /**/ double *);
 
 void classification(struct Map_info *, /**/
 		    struct Cell_head, /**/
@@ -54,16 +53,12 @@
 		    double, /**/
 		    double, /**/
 		    double, /**/
-		    int *, /**/ int, /**/ dbDriver *, /**/ char * /**/);
+		    int *, /**/ int, /**/ dbDriver *, /**/ char *, /**/ char *);
 
-int Insert(double, /**/ double, /**/ double, /**/ int, /**/ dbDriver * /**/);
+int Insert(double, /**/ double, /**/ double, /**/ int, /**/ dbDriver *, /**/ char *);
 
-int UpDate(double, /**/ double, /**/ double, /**/ int, /**/ dbDriver * /**/);
+int UpDate(double, /**/ double, /**/ double, /**/ int, /**/ dbDriver *, /**/ char *);
 
-int Select(double *, /**/
-	   double *, /**/ double *, /**/ int, /**/ dbDriver * /**/);
+int Select(double *, /**/ double *, /**/ double *, /**/ int, /**/ dbDriver *, /**/ char *);
 
-int Create_AuxEdge_Table(dbDriver *);
-int Create_Interpolation_Table(char *, dbDriver *);
 int Insert_Interpolation(double, int, dbDriver *, char *);
-int Drop_Aux_Table(dbDriver *);

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	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.lidar.edgedetection/main.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -25,6 +25,7 @@
  /*INCLUDES*/
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/Vect.h>
 #include <grass/dbmi.h>
@@ -32,19 +33,23 @@
 #include <grass/config.h>
 #include <grass/PolimiFunct.h>
 #include "edgedetection.h"
-int nsply, nsplx, line_out_counter, first_it;
+
+int nsply, nsplx, line_out_counter;
 double passoN, passoE;
 
-
 /**************************************************************************************
 **************************************************************************************/
 int main(int argc, char *argv[])
 {
     /* Variables' declarations */
+    int nsplx_adj, nsply_adj;
+    int nsubregion_col, nsubregion_row, subregion = 0, nsubregions = 0;
+    double N_extension, E_extension, orloE, orloN;
     int dim_vect, nparameters, BW, npoints;
-    double lambda_B, lambda_F, grad_H, grad_L, alpha, ew_resol, ns_resol,
-	mean;
-    char *dvr, *db, *mapset, table_interpolation[1024], table_name[1024];
+    double lambda_B, lambda_F, grad_H, grad_L, alpha, mean;
+    const char *dvr, *db, *mapset;
+    char table_interpolation[GNAME_MAX], table_name[GNAME_MAX];
+    char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
     int last_row, last_column, flag_auxiliar = FALSE;
 
@@ -56,6 +61,7 @@
     struct Map_info In, Out;
     struct Option *in_opt, *out_opt, *passoE_opt, *passoN_opt,
 	*lambdaF_opt, *lambdaB_opt, *gradH_opt, *gradL_opt, *alfa_opt;
+    struct Flag *spline_step_flag;
     struct GModule *module;
 
     struct Cell_head elaboration_reg, original_reg;
@@ -73,6 +79,12 @@
     module->description =
 	_("Detects the object's edges from a LIDAR data set.");
 
+    spline_step_flag = G_define_flag();
+    spline_step_flag->key = 'e';
+    spline_step_flag->label = _("Estimate point density and distance");
+    spline_step_flag->description =
+	_("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);
 
     out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
@@ -153,6 +165,10 @@
     grad_H = atof(gradH_opt->answer);
     grad_L = atof(gradL_opt->answer);
     alpha = atof(alfa_opt->answer);
+
+    grad_L = grad_L * grad_L;
+    grad_H = grad_H * grad_H;
+
     if (!(db = G__getenv2("DB_DATABASE", G_VAR_MAPSET)))
 	G_fatal_error(_("Unable to read name of database"));
 
@@ -160,8 +176,39 @@
 	G_fatal_error(_("Unable to read name of driver"));
 
     /* Setting auxiliar table's name */
-    sprintf(table_name, "%s_aux", out_opt->answer);
+    if (G__name_is_fully_qualified(out_opt->answer, xname, xmapset)) {
+	sprintf(table_name, "%s_aux", xname);
+	sprintf(table_interpolation, "%s_edge_Interpolation", xname);
+    }
+    else {
+	sprintf(table_name, "%s_aux", out_opt->answer);
+	sprintf(table_interpolation, "%s_edge_Interpolation", out_opt->answer);
+    }
 
+    /* Something went wrong in a previous v.lidar.edgedetection execution */
+    if (db_table_exists(dvr, db, table_name)) {
+	/* Start driver and open db */
+	driver = db_start_driver_open_database(dvr, db);
+	if (driver == NULL)
+	    G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
+			  dvr);
+	if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
+	    G_fatal_error(_("Old auxiliar table could not be dropped"));
+	db_close_database_shutdown_driver(driver);
+    }
+
+    /* Something went wrong in a previous v.lidar.edgedetection execution */
+    if (db_table_exists(dvr, db, table_interpolation)) {
+	/* Start driver and open db */
+	driver = db_start_driver_open_database(dvr, db);
+	if (driver == NULL)
+	    G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
+			  dvr);
+	if (P_Drop_Aux_Table(driver, table_interpolation) != DB_OK)
+	    G_fatal_error(_("Old auxiliar table could not be dropped"));
+	db_close_database_shutdown_driver(driver);
+    }
+
     /* Checking vector names */
     Vect_check_input_output_name(in_opt->answer, out_opt->answer,
 				 GV_FATAL_EXIT);
@@ -170,15 +217,33 @@
 	G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
     }
 
-    /* Open output vector */
-    if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z))
-	G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
-
     Vect_set_open_level(1);
     /* Open input vector */
     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;
+	if (P_estimate_splinestep(&In, &dens, &dist) == 0) {
+	    G_message("Estimated point density: %.4g", dens);
+	    G_message("Estimated mean distance between points: %.4g", dist);
+	}
+	else
+	    G_warning(_("No points in current region!"));
+	
+	Vect_close(&In);
+	exit(EXIT_SUCCESS);
+    }
+
+    /* Open output vector */
+    if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z))
+	G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
+
     /* Copy vector Head File */
     Vect_copy_head_data(&In, &Out);
     Vect_hist_copy(&In, &Out);
@@ -190,35 +255,68 @@
 	G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
 		      dvr);
 
-    if (Create_Interpolation_Table(out_opt->answer, driver) != DB_OK)
+    /* Create auxiliar and interpolation table */
+    if ((flag_auxiliar = P_Create_Aux4_Table(driver, table_name)) == FALSE)
+	G_fatal_error(_("It was impossible to create <%s>."), table_name);
+
+    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);
 
-    /* Fixxing 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_PASS times the east-west resolution. */
+    /*------------------------------------------------------------------
+      | 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
+      ----------------------------------------------------------------*/
 
-    ew_resol = original_reg.ew_res;
-    ns_resol = original_reg.ns_res;
-
+    /* Fixing parameters of the elaboration region */
     P_zero_dim(&dims);
 
-    dims.latoE = NSPLX_MAX * passoE;
-    dims.latoN = NSPLY_MAX * passoN;
-    dims.overlap = OVERLAP_SIZE * ew_resol;
+    nsplx_adj = NSPLX_MAX;
+    nsply_adj = NSPLY_MAX;
+    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 subregions */
+    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;
+
+    if (nsubregion_col < 0)
+	nsubregion_col = 0;
+    if (nsubregion_row < 0)
+	nsubregion_row = 0;
+
+    nsubregions = nsubregion_row * nsubregion_col;
+
     elaboration_reg.south = original_reg.north;
-
     last_row = FALSE;
-    first_it = TRUE;
 
     while (last_row == FALSE) {	/* For each row */
 
@@ -238,17 +336,23 @@
 
 	nsply =
 	    ceil((elaboration_reg.north - elaboration_reg.south) / passoN) +
-	    1;
+	    0.5;
+	/*
 	if (nsply > NSPLY_MAX) {
 	    nsply = NSPLY_MAX;
 	}
-	G_debug(1, _("nsply = %d"), nsply);
+	*/
+	G_debug(1, "nsply = %d", nsply);
 
 	elaboration_reg.east = original_reg.west;
 	last_column = FALSE;
 
 	while (last_column == FALSE) {	/* For each column */
 
+	    subregion++;
+	    if (nsubregions > 1)
+		G_message(_("subregion %d of %d"), subregion, nsubregions);
+
 	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 			  GENERAL_COLUMN);
 
@@ -265,15 +369,17 @@
 
 	    nsplx =
 		ceil((elaboration_reg.east - elaboration_reg.west) / passoE) +
-		1;
+		0.5;
+	    /*
 	    if (nsplx > NSPLX_MAX) {
 		nsplx = NSPLX_MAX;
 	    }
-	    G_debug(1, _("nsplx = %d"), nsplx);
+	    */
+	    G_debug(1, "nsplx = %d", nsplx);
 
 	    /*Setting the active region */
 	    dim_vect = nsplx * nsply;
-	    G_debug(1, _("read vector region map"));
+	    G_debug(1, "read vector region map");
 	    observ =
 		P_Read_Vector_Region_Map(&In, &elaboration_reg, &npoints,
 					 dim_vect, 1);
@@ -287,11 +393,11 @@
 		mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
 
 		/* Least Squares system */
-		G_debug(1, _("Allocation memory for bilinear interpolation"));
+		G_debug(1, _("Allocating memory for bilinear interpolation"));
 		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 */
 
@@ -306,9 +412,9 @@
 		    Q[i] = 1;	/* Q=I */
 		}
 
-		/*G_free (observ); */
+		G_free(observ);
 
-		G_debug(1, _("Bilinear interpolation"));
+		G_verbose_message(_("Bilinear interpolation"));
 		normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
 			       nsply, elaboration_reg.west,
 			       elaboration_reg.south, npoints, nparameters,
@@ -320,12 +426,12 @@
 		for (tn = 0; tn < nparameters; tn++)
 		    TN[tn] = 0;
 
-		G_debug(1, _("Allocation memory for bicubic interpolation"));
+		G_debug(1, _("Allocating memory for bicubic interpolation"));
 		BW = P_get_BandWidth(P_BICUBIC, nsply);
 		N = G_alloc_matrix(nparameters, BW);	/* Normal matrix */
 		parVect_bicub = G_alloc_vector(nparameters);	/* Bicubic parameters vector */
 
-		G_debug(1, _("Bicubic interpolation"));
+		G_verbose_message(_("Bicubic interpolation"));
 		normalDefBicubic(N, TN, Q, obsVect, passoE, passoN, nsplx,
 				 nsply, elaboration_reg.west,
 				 elaboration_reg.south, npoints, nparameters,
@@ -337,36 +443,30 @@
 		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 =
-			 Create_AuxEdge_Table(driver)) == FALSE)
-			G_fatal_error(_("It was impossible to create <%s>."),
-				      table_name);
-		}
-		G_debug(1, _("Point classification"));
+		G_verbose_message(_("Point classification"));
 		classification(&Out, elaboration_reg, general_box,
 			       overlap_box, obsVect, parVect_bilin,
 			       parVect_bicub, mean, alpha, grad_H, grad_L,
 			       dims.overlap, lineVect, npoints, driver,
-			       out_opt->answer);
+			       table_interpolation, table_name);
 
 		G_free_vector(parVect_bilin);
 		G_free_vector(parVect_bicub);
 		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 subregion. "
+			    "Consider changing the spline step."));
+	    }
 	}			/*! END WHILE; last_column = TRUE */
     }				/*! END WHILE; last_row = TRUE */
 
     /* Dropping auxiliar table */
     if (npoints > 0) {
 	G_debug(1, _("Dropping <%s>"), table_name);
-	if (Drop_Aux_Table(driver) != DB_OK)
+	if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
 	    G_warning(_("Auxiliar table could not be dropped"));
     }
 
@@ -374,7 +474,6 @@
 
     Vect_close(&In);
 
-    sprintf(table_interpolation, "%s_edge_Interpolation", out_opt->answer);
     Vect_map_add_dblink(&Out, F_INTERPOLATION, NULL, table_interpolation,
 			"id", db, dvr);
 

Modified: grass/branches/releasebranch_6_4/vector/lidar/v.lidar.growing/main.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.lidar.growing/main.c	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.lidar.growing/main.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -22,6 +22,7 @@
  /*INCLUDES*/
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #include <grass/config.h>
 #include <grass/gis.h>
 #include <grass/Vect.h>
@@ -38,6 +39,8 @@
 
     /* Variables' declarations */
     int row, nrows, col, ncols, MaxPoints;
+    int nsubregion_col, nsubregion_row;
+    int subregion = 0, nsubregions = 0;
     int last_row, last_column;
     int nlines, nlines_first, line_num;
     int more;
@@ -45,7 +48,8 @@
     double Z_interp;
     double Thres_j, Thres_d, ew_resol, ns_resol;
     double minNS, minEW, maxNS, maxEW;
-    char *mapset, buf[1024];
+    const char *mapset;
+    char buf[1024], table_name[GNAME_MAX];
     char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
     int colorBordo, ripieno, conta, lungPunti, lungHull, xi, c1, c2;
@@ -68,8 +72,6 @@
     dbString sql;
     dbTable *table;
     dbCursor cursor;
-    const char *p;
-    int found;
 
 /*------------------------------------------------------------------------------------------*/
     /* Options' declaration */ ;
@@ -80,6 +82,8 @@
 	  "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);
 
@@ -121,27 +125,22 @@
     Vect_check_input_output_name(in_opt->answer, out_opt->answer,
 				 GV_FATAL_EXIT);
 
-    sprintf(xname, "%s", in_opt->answer);
-    found = 0;
-    for (p = in_opt->answer; *p; p++)
-	if (*p == '@') {
-	    found = 1;
-	    break;
-	}
-
-    if (found) {
-	if (G__name_is_fully_qualified(in_opt->answer, xname, xmapset) < 0)	/* strip off mapset from input for SQL */
-	    G_fatal_error(_("Vector map <%s> not found"), xname);
-    }
     if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL) {
 	G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
     }
 
-    /*Vect_set_open_level (2);          WITH TOPOLOGY */
+    /* Setting auxiliar table's name */
+    if (G__name_is_fully_qualified(in_opt->answer, xname, xmapset)) {
+	sprintf(table_name, "%s_edge_Interpolation", xname);
+    }
+    else
+	sprintf(table_name, "%s_edge_Interpolation", in_opt->answer);
+
     Vect_set_open_level(1);	/* WITHOUT TOPOLOGY */
     if (Vect_open_old(&In, in_opt->answer, mapset) < 1)
 	G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
 
+    Vect_set_open_level(1);	/* WITHOUT TOPOLOGY */
     if (Vect_open_old(&First, first_opt->answer, mapset) < 1)
 	G_fatal_error(_("Unable to open vector map <%s>"), first_opt->answer);
 
@@ -167,40 +166,65 @@
 	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);
-
-    /*  Fixxing 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;
-
-    /* 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);
 
-    sprintf(buf, "SELECT Interp,ID FROM %s_edge_Interpolation", xname);
+    sprintf(buf, "SELECT Interp,ID FROM %s", table_name);
     G_debug(1, "buf: %s", buf);
     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_edge_Interpolation>"), xname);
+	G_fatal_error(_("Unable to open table <%s>"), table_name);
 
     count_obj = 1;
 
-    nlines = Vect_get_num_lines(&In);
+    /* no topology, get number of lines in input vector */
+    nlines = 0;
     points = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
+    Vect_rewind(&In);
+    while (Vect_read_next_line(&In, points, Cats) > 0) {
+	nlines++;
+    }
+    Vect_rewind(&In);
 
-    nlines_first = Vect_get_num_lines(&First);
+    /* 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();
+    Vect_rewind(&First);
+    while (Vect_read_next_line(&First, points_first, Cats_first) > 0) {
+	nlines_first++;
+    }
+    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 subregions */
+    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;
+
+    nsubregions = 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;
@@ -220,6 +244,10 @@
 	while (last_column == FALSE) {	/* For each strip of LATO columns */
 	    BOUND_BOX elaboration_box;
 
+	    subregion++;
+	    if (nsubregions > 1)
+		G_message(_("subregion %d of %d"), subregion, nsubregions);
+
 	    elaboration_reg.west = elaboration_reg.east;
 	    if (elaboration_reg.west < original_reg.west)	/* First column */
 		elaboration_reg.west = original_reg.west;
@@ -231,10 +259,13 @@
 		last_column = TRUE;
 	    }
 
-	    /*Setting the active region */
-	    G_set_window(&elaboration_reg);
-	    nrows = G_window_rows();
-	    ncols = G_window_cols();
+	    /* Setting the active region */
+	    elaboration_reg.ns_res = ns_resol;
+	    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);
@@ -257,13 +288,13 @@
 		}
 	    }
 
+	    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)) &&
@@ -279,6 +310,37 @@
 			(int)(G_easting_to_col
 			      (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)
@@ -353,13 +415,9 @@
 		     (points_first->x[0], points_first->y[0],
 		      points_first->z[0], &elaboration_box)) &&
 		    ((points->x[0] != elaboration_reg.west) ||
-		     (points->x[0] == original_reg.west)) && ((points->y[0]
-							       !=
-							       elaboration_reg.
-							       north) ||
-							      (points->y[0] ==
-							       original_reg.
-							       north))) {
+		     (points->x[0] == original_reg.west)) &&
+		    ((points->y[0] != elaboration_reg.north) ||
+		     (points->y[0] == original_reg.north))) {
 
 		    row =
 			(int)(G_northing_to_row
@@ -379,7 +437,7 @@
 
 	    /* 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);
@@ -388,6 +446,7 @@
 		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) &&
@@ -416,7 +475,7 @@
 				     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);
 

Modified: grass/branches/releasebranch_6_4/vector/lidar/v.outlier/main.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.outlier/main.c	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.outlier/main.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -20,6 +20,7 @@
  /*INCLUDES*/
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #include <grass/config.h>
 #include <grass/gis.h>
 #include <grass/Vect.h>
@@ -28,16 +29,22 @@
 #include <grass/PolimiFunct.h>
 #include "outlier.h"
     /* GLOBAL VARIABLES DEFINITIONS */
-int nsply, nsplx, first_it;
+int nsply, nsplx;
 double passoN, passoE, Thres_Outlier;
 
 /*--------------------------------------------------------------------------------------*/
 int main(int argc, char *argv[])
 {
     /* Variables' declarations */
+    int nsplx_adj, nsply_adj;
+    int nsubregion_col, nsubregion_row;
+    int subregion = 0, nsubregions = 0;
+    double N_extension, E_extension, orloE, orloN;
     int dim_vect, nparameters, BW, npoints;
-    double ew_resol, ns_resol, mean, lambda;
-    char *dvr, *db, *mapset, table_name[1024];
+    double mean, lambda;
+    const char *dvr, *db, *mapset;
+    char table_name[GNAME_MAX];
+    char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
     int last_row, last_column, flag_auxiliar = FALSE;
 
@@ -49,7 +56,7 @@
     struct Map_info In, Out, Outlier, Qgis;
     struct Option *in_opt, *out_opt, *outlier_opt, *qgis_opt, *passoE_opt,
 	*passoN_opt, *lambda_f_opt, *Thres_O_opt;
-    /*struct Flag *qgis_flag; */
+    struct Flag *spline_step_flag;
     struct GModule *module;
 
     struct Cell_head elaboration_reg, original_reg;
@@ -66,6 +73,12 @@
     module->keywords = _("vector, statistics");
     module->description = _("Removes outliers from vector point data.");
 
+    spline_step_flag = G_define_flag();
+    spline_step_flag->key = 'e';
+    spline_step_flag->label = _("Estimate point density and distance");
+    spline_step_flag->description =
+	_("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);
 
     out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
@@ -133,9 +146,7 @@
     lambda = atof(lambda_f_opt->answer);
     Thres_Outlier = atof(Thres_O_opt->answer);
 
-    /* Setting auxiliar table's name */
-    /* sprintf(table_name, "%s_aux", out_opt->answer); */
-    sprintf(table_name, "Auxiliar_outlier_table");
+    flag_auxiliar = FALSE;
 
     /* Checking vector names */
     Vect_check_input_output_name(in_opt->answer, out_opt->answer,
@@ -145,6 +156,49 @@
 	G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
     }
 
+    /* Setting auxiliar table's name */
+    if (G__name_is_fully_qualified(out_opt->answer, xname, xmapset)) {
+	sprintf(table_name, "%s_aux", xname);
+    }
+    else
+	sprintf(table_name, "%s_aux", out_opt->answer);
+
+    /* Something went wrong in a previous v.outlier execution */
+    if (db_table_exists(dvr, db, table_name)) {
+	/* Start driver and open db */
+	driver = db_start_driver_open_database(dvr, db);
+	if (driver == NULL)
+	    G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
+			  dvr);
+	if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
+	    G_fatal_error(_("Old auxiliar table could not be dropped"));
+	db_close_database_shutdown_driver(driver);
+    }
+
+    Vect_set_open_level(1);
+    /* Open input vector */
+    if (1 > Vect_open_old(&In, in_opt->answer, mapset))
+	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;
+	if (P_estimate_splinestep(&In, &dens, &dist) == 0) {
+	    G_message("Estimated point density: %.4g", dens);
+	    G_message("Estimated mean distance between points: %.4g", dist);
+	}
+	else
+	    G_warning(_("No points in current region!"));
+	
+	Vect_close(&In);
+	exit(EXIT_SUCCESS);
+    }
+
     /* Open output vector */
     if (qgis_opt->answer)
 	if (0 > Vect_open_new(&Qgis, qgis_opt->answer, WITHOUT_Z))
@@ -162,12 +216,6 @@
 	G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
     }
 
-    Vect_set_open_level(1);
-    /* Open input vector */
-    if (1 > Vect_open_old(&In, in_opt->answer, mapset))
-	G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
-		      in_opt->answer);
-
     /* Copy vector Head File */
     Vect_copy_head_data(&In, &Out);
     Vect_hist_copy(&In, &Out);
@@ -189,37 +237,64 @@
 	G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
 		      dvr);
 
-    /* Something went wrong in a previous v.outlier execution */
-    if (db_table_exists(dvr, db, table_name)) {
-	if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
-	    G_fatal_error(_("Old auxiliar table could not be dropped"));
-    }
+    /* 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);
 
+    /*------------------------------------------------------------------
+      | 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 */
-    /*! Each original_region will be divided into several subregions. These
-     *  subregion will be overlapped by its neibourgh subregions. This overlapping
-     *  is calculated as OVERLAP_PASS times the east-west resolution. */
+    P_zero_dim(&dims);
 
-    ew_resol = original_reg.ew_res;
-    ns_resol = original_reg.ns_res;
+    nsplx_adj = NSPLX_MAX;
+    nsply_adj = NSPLY_MAX;
+    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);
 
-    P_zero_dim(&dims);
+    G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
+    G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
 
-    dims.latoE = NSPLX_MAX * passoE;
-    dims.latoN = NSPLY_MAX * passoN;
-    dims.overlap = OVERLAP_SIZE * ew_resol;
-    P_get_orlo(P_BICUBIC, &dims, passoE, passoN);
+    /* calculate number of subregions */
+    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;
+
+    if (nsubregion_col < 0)
+	nsubregion_col = 0;
+    if (nsubregion_row < 0)
+	nsubregion_row = 0;
+
+    nsubregions = nsubregion_row * nsubregion_col;
+
     elaboration_reg.south = original_reg.north;
-
     last_row = FALSE;
-    first_it = TRUE;
 
     while (last_row == FALSE) {	/* For each row */
 
@@ -239,10 +314,12 @@
 
 	nsply =
 	    ceil((elaboration_reg.north - elaboration_reg.south) / passoN) +
-	    1;
+	    0.5;
+	/*
 	if (nsply > NSPLY_MAX) {
 	    nsply = NSPLY_MAX;
 	}
+	*/
 	G_debug(1, "nsply = %d", nsply);
 
 	elaboration_reg.east = original_reg.west;
@@ -250,6 +327,10 @@
 
 	while (last_column == FALSE) {	/* For each column */
 
+	    subregion++;
+	    if (nsubregions > 1)
+		G_message(_("subregion %d of %d"), subregion, nsubregions);
+
 	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 			  GENERAL_COLUMN);
 
@@ -266,10 +347,12 @@
 
 	    nsplx =
 		ceil((elaboration_reg.east - elaboration_reg.west) / passoE) +
-		1;
+		0.5;
+	    /*
 	    if (nsplx > NSPLX_MAX) {
 		nsplx = NSPLX_MAX;
 	    }
+	    */
 	    G_debug(1, "nsplx = %d", nsplx);
 
 	    /*Setting the active region */
@@ -305,7 +388,9 @@
 		    Q[i] = 1;	/* Q=I */
 		}
 
-		G_debug(1, "Bilinear interpolation");
+		G_free(observ);
+
+		G_verbose_message(_("Bilinear interpolation"));
 		normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
 			       nsply, elaboration_reg.west,
 			       elaboration_reg.south, npoints, nparameters,
@@ -317,22 +402,15 @@
 		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,
 			      general_box, overlap_box, obsVect, parVect,
-			      mean, dims.overlap, lineVect, npoints, driver);
+			      mean, dims.overlap, lineVect, npoints, driver, table_name);
 		else
 		    P_Outlier(&Out, &Outlier, NULL, elaboration_reg,
 			      general_box, overlap_box, obsVect, parVect,
-			      mean, dims.overlap, lineVect, npoints, driver);
+			      mean, dims.overlap, lineVect, npoints, driver, table_name);
 
 
 		G_free_vector(parVect);
@@ -340,8 +418,11 @@
 		G_free_ivector(lineVect);
 
 	    }			/*! END IF; npoints > 0 */
-	    G_free(observ);
-	    first_it = FALSE;
+	    else {
+		G_free(observ);
+		G_warning(_("No data within this subregion. "
+			    "Consider changing the spline step."));
+	    }
 	}			/*! END WHILE; last_column = TRUE */
     }				/*! END WHILE; last_row = TRUE */
 

Modified: grass/branches/releasebranch_6_4/vector/lidar/v.outlier/outlier.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.outlier/outlier.c	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.outlier/outlier.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -1,31 +1,36 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/Vect.h>
 #include <grass/dbmi.h>
 
 #include "outlier.h"
 
+extern double Thres_Outlier;
+
 void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 	       struct Map_info *Qgis, struct Cell_head Elaboration,
-	       BOUND_BOX General, BOUND_BOX Overlap, double **obs,
-	       double *parBilin, double mean, double overlap, int *line_num,
-	       int num_points, dbDriver * driver)
+	       BOUND_BOX General, BOUND_BOX Overlap,
+	       double **obs, double *parBilin, double mean, double overlap,
+	       int *line_num, int num_points, dbDriver * driver,
+	       char *tab_name)
 {
     int i;
-    double interpolation, weight, residual, eta, csi, gradient[2];
-
+    double interpolation, weight, residual, eta, csi;
     extern int nsplx, nsply;
     extern double passoN, passoE;
-
     struct line_pnts *point;
     struct line_cats *categories;
 
     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);
 	Vect_reset_cats(categories);
 
@@ -57,60 +62,53 @@
 
 	    }
 	    else {
-		if ((*point->x > Overlap.E) && (*point->x != General.E)) {
+		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;
 
-		    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);
-
-			gradient[0] *= weight;
-			gradient[1] *= weight;
 			interpolation *= weight;
 
-			if (Select_Outlier(&gradient[0], line_num[i], driver)
-			    != DB_OK)
+			if (Select_Outlier(&interpolation, line_num[i],
+					   driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to read the database"));
 
-			if (UpDate_Outlier(gradient[0], line_num[i], driver)
-			    != DB_OK)
+			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;
 
-			if (Insert_Outlier
-			    (interpolation * weight, line_num[i],
-			     driver) != DB_OK)
+			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;
 
-			if (Insert_Outlier
-			    (interpolation * weight, line_num[i],
-			     driver) != DB_OK)
+			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->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;
 			interpolation *= weight;
 
-			if (Select_Outlier(&gradient[0], line_num[i], driver)
-			    != DB_OK)
+			if (Select_Outlier(&interpolation, line_num[i],
+			                   driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to read the database"));
 
 			residual = *point->z - interpolation;
@@ -122,38 +120,34 @@
 				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;
 
-			gradient[0] *= weight;
-			gradient[1] *= weight;
 			interpolation *= weight;
 
-			if (Select_Outlier(&gradient[0], line_num[i], driver)
-			    != DB_OK)
+			if (Select_Outlier(&interpolation, line_num[i],
+			                   driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to read the database"));
 
-			if (UpDate_Outlier(gradient[0], line_num[i], driver)
-			    != DB_OK)
+			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;
 
-			gradient[0] *= weight;
-			gradient[1] *= weight;
 			interpolation *= weight;
 
-			if (Select_Outlier(&gradient[0], line_num[i], driver)
-			    != DB_OK)
+			if (Select_Outlier(&interpolation, line_num[i],
+			                   driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to read the database"));
 
 			residual = *point->z - interpolation;
@@ -169,18 +163,15 @@
 			    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;
 
-			gradient[0] *= weight;
-			gradient[1] *= weight;
 			interpolation *= weight;
 
-			if (Select_Outlier(&gradient[0], line_num[i], driver)
-			    != DB_OK)
+			if (Select_Outlier(&interpolation, line_num[i],
+			                   driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to read the database"));
 
 			residual = *point->z - interpolation;
@@ -195,14 +186,14 @@
 			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;
 
-			if (Insert_Outlier
-			    (interpolation * weight, line_num[i], driver)
-			    != DB_OK)
+			interpolation *= weight;
+
+			if (Insert_Outlier(interpolation, line_num[i],
+					   driver, tab_name) != DB_OK)
 			    G_fatal_error(_("Impossible to write in the database"));
 
 		    }		/*else (1) */
@@ -210,42 +201,56 @@
 	    }
 	}			/*end if obs */
     }				/*end for */
+
+    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 */
 
-int Insert_Outlier(double Interp, int line_num, dbDriver * driver)
+int Insert_Outlier(double Interp, int line_num, dbDriver * driver,
+		   char *tab_name)
 {
-
     char buf[1024];
     dbString sql;
+    int ret;
 
     db_init_string(&sql);
-    sprintf(buf, "INSERT INTO Auxiliar_outlier_table (ID, Interp)");
+    sprintf(buf, "INSERT INTO %s (ID, Interp)", tab_name);
     db_append_string(&sql, buf);
     sprintf(buf, " VALUES (%d, %lf)", line_num, Interp);
     db_append_string(&sql, buf);
 
-    return db_execute_immediate(driver, &sql);
+    ret = db_execute_immediate(driver, &sql);
+    db_free_string(&sql);
+
+    return ret;
 }
 
-int UpDate_Outlier(double Interp, int line_num, dbDriver * driver)
+int UpDate_Outlier(double Interp, int line_num, dbDriver * driver,
+		   char *tab_name)
 {
-
     char buf[1024];
     dbString sql;
+    int ret;
 
     db_init_string(&sql);
-    sprintf(buf, "UPDATE Auxiliar_outlier_table SET Interp=%lf WHERE ID=%d",
-	    Interp, line_num);
+    sprintf(buf, "UPDATE %s SET Interp=%lf WHERE ID=%d",
+	    tab_name, Interp, line_num);
     db_append_string(&sql, buf);
 
-    return db_execute_immediate(driver, &sql);
+    ret = db_execute_immediate(driver, &sql);
+    db_free_string(&sql);
+
+    return ret;
 }
 
-int Select_Outlier(double *Interp, int line_num, dbDriver * driver)
+int Select_Outlier(double *Interp, int line_num, dbDriver * driver,
+		   char *tab_name)
 {
-
     int more;
     char buf[1024];
     dbString sql;
@@ -255,8 +260,7 @@
     dbValue *Interp_value;
 
     db_init_string(&sql);
-    sprintf(buf, "SELECT ID, Interp FROM Auxiliar_outlier_table WHERE ID=%d",
-	    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)
@@ -276,36 +280,42 @@
 	*Interp += db_get_value_double(Interp_value);
     }
     db_close_cursor(&cursor);
+    db_free_string(&sql);
     return DB_OK;
 }
 
 int P_is_outlier(double pippo)
 {
-    extern double Thres_Outlier;
-
     if (fabs(pippo) < Thres_Outlier)
 	return FALSE;
 
     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

Modified: grass/branches/releasebranch_6_4/vector/lidar/v.outlier/outlier.h
===================================================================
--- grass/branches/releasebranch_6_4/vector/lidar/v.outlier/outlier.h	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.outlier/outlier.h	2010-02-17 17:59:30 UTC (rev 41063)
@@ -16,12 +16,12 @@
 	  BOUND_BOX, /**/
 	  double **, /**/
 	  double *, /**/
-	  double, /**/ double, /**/ int *, /**/ int, /**/ dbDriver * /**/);
+	  double, /**/ double, /**/ int *, /**/ int, /**/ dbDriver *,  /**/ char *);
 
-int Insert_Outlier(double, /**/ int, /**/ dbDriver * /**/);
+int Insert_Outlier(double, /**/ int, /**/ dbDriver *, /**/ char *);
 
-int UpDate_Outlier(double, /**/ int, /**/ dbDriver * /**/);
+int UpDate_Outlier(double, /**/ int, /**/ dbDriver *, /**/ char *);
 
-int Select_Outlier(double *, /**/ int, /**/ dbDriver * /**/);
+int Select_Outlier(double *, /**/ int, /**/ dbDriver *, /**/ char *);
 
 int P_is_outlier(double);

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	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.surf.bspline/crosscorr.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -20,6 +20,7 @@
  /*INCLUDES*/
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/Vect.h>
 #include <grass/dbmi.h>
@@ -48,7 +49,8 @@
     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 */
@@ -151,7 +153,7 @@
 			  nsplx, nsply);
 
 	BW = P_get_BandWidth(bilin, nsply);
-	 /**/
+	/**/
 	    /*Least Squares system */
 	    N = G_alloc_matrix(nparam_spl, BW);	/* Normal matrix */
 	TN = G_alloc_vector(nparam_spl);	/* vector */

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	2010-02-17 17:56:52 UTC (rev 41062)
+++ grass/branches/releasebranch_6_4/vector/lidar/v.surf.bspline/main.c	2010-02-17 17:59:30 UTC (rev 41063)
@@ -1,5 +1,5 @@
 
-/***********************************************************************
+/**********************************************************************
  *
  * MODULE:       v.surf.bspline
  *
@@ -15,11 +15,12 @@
  *               Read the file COPYING that comes with GRASS
  *               for details.
  *
- **************************************************************************/
+ **********************************************************************/
 
- /*INCLUDES*/
+/* INCLUDES */
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/Vect.h>
 #include <grass/dbmi.h>
@@ -27,22 +28,25 @@
 #include <grass/config.h>
 #include <grass/PolimiFunct.h>
 #include "bspline.h"
-    /*GLOBAL VARIABLES */
+
+/* GLOBAL VARIABLES */
 int bspline_field;
 char *bspline_column;
 
-/*-------------------------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------*/
 int main(int argc, char *argv[])
 {
-
-    /* Variables' declarations */
-    int nsply, nsplx, nlines, nrows, ncols;
+    /* Variable declarations */
+    int nsply, nsplx, nrows, ncols, nsplx_adj, nsply_adj;
     int nsubregion_col, nsubregion_row, subregion_row, subregion_col;
+    int subregion = 0, nsubregions = 0;
     int last_row, last_column, grid, bilin, ext, flag_auxiliar, cross;	/* booleans */
     double passoN, passoE, lambda, mean;
     double N_extension, E_extension, orloE, orloN;
 
-    char *mapset, *dvr, *db, *vector, *map, table_name[1024], title[64];
+    const char *mapset, *dvr, *db, *vector, *map;
+    char table_name[GNAME_MAX], title[64];
+    char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
     int dim_vect, nparameters, BW;
     int *lineVect;		/* Vector restoring primitive's ID */
@@ -50,15 +54,15 @@
     double *TN, *Q, *parVect;	/* Interpolating and least-square vectors */
     double **N, **obsVect;	/* Interpolation and least-square matrix */
 
-    /* Structs' declarations */
+    /* Structs declarations */
     int raster;
     struct Map_info In, In_ext, Out;
     struct History history;
 
     struct GModule *module;
     struct Option *in_opt, *in_ext_opt, *out_opt, *out_map_opt, *passoE_opt,
-	*passoN_opt, *lambda_f_opt, *type, *dfield_opt, *col_opt;
-    struct Flag *cross_corr_flag;
+	*passoN_opt, *lambda_f_opt, *type_opt, *dfield_opt, *col_opt;
+    struct Flag *cross_corr_flag, *spline_step_flag;
 
     struct Reg_dimens dims;
     struct Cell_head elaboration_reg, original_reg;
@@ -73,28 +77,31 @@
     struct field_info *Fi;
     dbDriver *driver, *driver_cats;
 
-    /*-------------------------------------------------------------------------------------------*/
-    /* Options' declaration */
-    module = G_define_module(); {
-	module->keywords = _("vector, interpolation");
-	module->description =
-	    _("Bicubic or bilinear spline interpolation with Tykhonov regularization.");
-    }
+    /*----------------------------------------------------------------*/
+    /* Options declarations */
+    module = G_define_module();
+    module->keywords = _("vector, interpolation");
+    module->description =
+	_("Bicubic or bilinear spline interpolation with Tykhonov regularization.");
 
-    cross_corr_flag = G_define_flag(); {
-	cross_corr_flag->key = 'c';
-	cross_corr_flag->description =
-	    _("Find the best Tykhonov regularizing parameter using a \"leave-one-out\" cross validation method");
-    }
+    cross_corr_flag = G_define_flag();
+    cross_corr_flag->key = 'c';
+    cross_corr_flag->description =
+	_("Find the best Tykhonov regularizing parameter using a \"leave-one-out\" cross validation method");
 
+    spline_step_flag = G_define_flag();
+    spline_step_flag->key = 'e';
+    spline_step_flag->label = _("Estimate point density and distance");
+    spline_step_flag->description =
+	_("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_ext_opt = G_define_standard_option(G_OPT_V_INPUT); {
-	in_ext_opt->key = "sparse";
-	in_ext_opt->required = NO;
-	in_ext_opt->description =
-	    _("Name of input vector map of sparse points");
-    }
+    in_ext_opt = G_define_standard_option(G_OPT_V_INPUT);
+    in_ext_opt->key = "sparse";
+    in_ext_opt->required = NO;
+    in_ext_opt->description =
+	_("Name of input vector map of sparse points");
 
     out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
     out_opt->required = NO;
@@ -103,69 +110,71 @@
     out_map_opt->key = "raster";
     out_map_opt->required = NO;
 
-    passoE_opt = G_define_option(); {
-	passoE_opt->key = "sie";
-	passoE_opt->type = TYPE_DOUBLE;
-	passoE_opt->required = NO;
-	passoE_opt->answer = "4";
-	passoE_opt->description =
-	    _("Length of each spline step in the east-west direction");
-	passoE_opt->guisection = _("Settings");
-    }
+    passoE_opt = G_define_option();
+    passoE_opt->key = "sie";
+    passoE_opt->type = TYPE_DOUBLE;
+    passoE_opt->required = NO;
+    passoE_opt->answer = "4";
+    passoE_opt->description =
+	_("Length of each spline step in the east-west direction");
+    passoE_opt->guisection = _("Settings");
 
-    passoN_opt = G_define_option(); {
-	passoN_opt->key = "sin";
-	passoN_opt->type = TYPE_DOUBLE;
-	passoN_opt->required = NO;
-	passoN_opt->answer = "4";
-	passoN_opt->description =
-	    _("Length of each spline step in the north-south direction");
-	passoN_opt->guisection = _("Settings");
-    }
+    passoN_opt = G_define_option();
+    passoN_opt->key = "sin";
+    passoN_opt->type = TYPE_DOUBLE;
+    passoN_opt->required = NO;
+    passoN_opt->answer = "4";
+    passoN_opt->description =
+	_("Length of each spline step in the north-south direction");
+    passoN_opt->guisection = _("Settings");
 
-    type = G_define_option(); {
-	type->key = "type";
-	type->type = TYPE_STRING;
-	type->required = NO;
-	type->description = _("Spline interpolation algorithm");
-	type->options = "bilinear,bicubic";
-	type->answer = "bilinear";
-	type->guisection = _("Settings");
-    }
+    type_opt = G_define_option();
+    type_opt->key = "method";
+    type_opt->type = TYPE_STRING;
+    type_opt->required = NO;
+    type_opt->description = _("Spline interpolation algorithm");
+    type_opt->options = "bilinear,bicubic";
+    type_opt->answer = "bilinear";
+    type_opt->guisection = _("Settings");
 
-    lambda_f_opt = G_define_option(); {
-	lambda_f_opt->key = "lambda_i";
-	lambda_f_opt->type = TYPE_DOUBLE;
-	lambda_f_opt->required = NO;
-	lambda_f_opt->description = _("Tykhonov regularization parameter (affects smoothing)");
-	lambda_f_opt->answer = "1";
-	lambda_f_opt->guisection = _("Settings");
-    }
+    lambda_f_opt = G_define_option();
+    lambda_f_opt->key = "lambda_i";
+    lambda_f_opt->type = TYPE_DOUBLE;
+    lambda_f_opt->required = NO;
+    lambda_f_opt->description = _("Tykhonov regularization parameter (affects smoothing)");
+    lambda_f_opt->answer = "1";
+    lambda_f_opt->guisection = _("Settings");
 
-    dfield_opt = G_define_standard_option(G_OPT_V_FIELD); {
-	dfield_opt->description =
-	    _("If set to 0, z coordinates are used. (3D vector only)");
-	dfield_opt->answer = "0";
-	dfield_opt->guisection = _("Settings");
-	dfield_opt->gisprompt = "old_layer,layer,layer_zero";
-    }
+    dfield_opt = G_define_standard_option(G_OPT_V_FIELD);
+    dfield_opt->description =
+	_("If set to -1, z coordinates are used. (3D vector only)");
+    dfield_opt->answer = "-1";
+    dfield_opt->guisection = _("Settings");
 
-    col_opt = G_define_option(); {
-	col_opt->key = "column";
-	col_opt->type = TYPE_STRING;
-	col_opt->required = NO;
-	col_opt->description =
-	    _("Attribute table column with values to interpolate (if layer>0)");
-	col_opt->guisection = _("Settings");
-    }
+    col_opt = G_define_option();
+    col_opt->key = "column";
+    col_opt->type = TYPE_STRING;
+    col_opt->required = NO;
+    col_opt->description =
+	_("Attribute table column with values to interpolate (if layer>0)");
+    col_opt->guisection = _("Settings");
 
-/*-------------------------------------------------------------------------------------------*/
+    /*----------------------------------------------------------------*/
     /* Parsing */
     G_gisinit(argv[0]);
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
-    if (!strcmp(type->answer, "bilinear"))
+    vector = out_opt->answer;
+    map = out_map_opt->answer;
+
+    if (vector && map)
+	G_fatal_error(_("Choose either vector or raster output, not both"));
+
+    if (!vector && !map && !cross_corr_flag->answer)
+	G_fatal_error(_("No raster or vector or cross-validation output"));
+
+    if (!strcmp(type_opt->answer, "bilinear"))
 	bilin = P_BILINEAR;
     else
 	bilin = P_BICUBIC;
@@ -175,12 +184,8 @@
     lambda = atof(lambda_f_opt->answer);
     bspline_field = atoi(dfield_opt->answer);
     bspline_column = col_opt->answer;
-
+    
     flag_auxiliar = FALSE;
-    if (cross_corr_flag->answer) {
-    }
-    vector = out_opt->answer;
-    map = out_map_opt->answer;
 
     if (!(db = G__getenv2("DB_DATABASE", G_VAR_MAPSET)))
 	G_fatal_error(_("Unable to read name of database"));
@@ -189,39 +194,26 @@
 	G_fatal_error(_("Unable to read name of driver"));
 
     /* Setting auxiliar table's name */
-    if (vector)
-	sprintf(table_name, "%s_aux", out_opt->answer);
+    if (vector) {
+	if (G__name_is_fully_qualified(out_opt->answer, xname, xmapset)) {
+	    sprintf(table_name, "%s_aux", xname);
+	}
+	else
+	    sprintf(table_name, "%s_aux", out_opt->answer);
+    }
 
-    /* Open driver and database */
-    if (db_table_exists(dvr, db, table_name)) {	/*Something went wrong in a 
-						   previous v.surf.bspline execution */
-	dbString sql;
-	char buf[1024];
-
-
+    /* Something went wrong in a previous v.surf.bspline execution */
+    if (db_table_exists(dvr, db, table_name)) {
+	/* Start driver and open db */
 	driver = db_start_driver_open_database(dvr, db);
 	if (driver == NULL)
-	    G_fatal_error(_("No database connection for driver <%s> is defined. "
-			   "Run db.connect."), dvr);
-
-	db_init_string(&sql);
-	db_zero_string(&sql);
-	sprintf(buf, "drop table %s", table_name);
-	db_append_string(&sql, buf);
-	if (db_execute_immediate(driver, &sql) != DB_OK)
-	    G_fatal_error(_("It was not possible to drop <%s> table. "
-			    "Nothing will be done. Try to drop it manually."),
-			  table_name);
-
+	    G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
+			  dvr);
+	if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
+	    G_fatal_error(_("Old auxiliar table could not be dropped"));
 	db_close_database_shutdown_driver(driver);
     }
 
-    if (vector && map)
-	G_fatal_error(_("Choose either vector or raster output, not both"));
-#ifdef notdef
-    if (!vector && !map && !cross_corr_flag->answer)
-	G_fatal_error(_("No raster nor vector output"));
-#endif
     /* Open input vector */
     if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL)
 	G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
@@ -231,6 +223,45 @@
 	G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
 		      in_opt->answer);
 
+    /* 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) {
+	double dens, dist;
+	if (P_estimate_splinestep(&In, &dens, &dist) == 0) {
+	    G_message("Estimated point density: %.4g", dens);
+	    G_message("Estimated mean distance between points: %.4g", dist);
+	}
+	else
+	    G_warning(_("No points in current region!"));
+	
+	Vect_close(&In);
+	exit(EXIT_SUCCESS);
+    }
+
+    /*----------------------------------------------------------------*/
+    /* Cross-correlation begins */
+    if (cross_corr_flag->answer) {
+	G_debug(1, "CrossCorrelation()");
+	cross = cross_correlation(&In, passoE, passoN);
+
+	if (cross != TRUE)
+	    G_fatal_error(_("Cross validation didn't finish correctly"));
+	else {
+	    G_debug(1, "Cross validation finished correctly");
+
+	    Vect_close(&In);
+
+	    G_done_msg(_("Cross validation finished for sie = %f and sin = %f"), passoE, passoN);
+	    exit(EXIT_SUCCESS);
+	}
+    }
+
     /* Open input ext vector */
     if (!in_ext_opt->answer) {
 	ext = FALSE;
@@ -284,16 +315,12 @@
     G_set_fp_type(DCELL_TYPE);
     if (!vector && map) {
 	grid = TRUE;
-	/*
-	   if (G_find_cell (out_map_opt->answer, G_mapset()) != NULL) 
-	   G_fatal_error (_("Raster <%s> already exist."), out_map_opt->answer);
-	 */
-
 	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);
     }
 
+    /* read z values from attribute table */
     if (bspline_field > 0) {
 	db_CatValArray_init(&cvarr);
 	Fi = Vect_get_field(&In, bspline_field);
@@ -324,54 +351,30 @@
 	db_close_database_shutdown_driver(driver_cats);
     }
 
-/*-------------------------------------------------------------------------------------------*/
-    /*Cross-correlation begins */
-    if (cross_corr_flag->answer) {	/* CROSS-CORRELATION WILL BE DONE */
-	G_debug(1, "CrossCorrelation()");
-	/*cross = cross_correlation (&In, &passoE, &passoN, &lambda); */
-	cross = cross_correlation(&In, passoE, passoN);
-
-	if (cross != TRUE)
-	    G_fatal_error(_("Cross validation didn't finish correctly"));
-	else {
-	    G_debug(1, "Cross validation finished correctly");
-
-	    Vect_close(&In);
-	    if (ext != FALSE)
-		Vect_close(&In_ext);
-	    if (vector)
-		Vect_close(&Out);
-
-	    if (map)
-		G_close_cell(raster);
-
-	    G_done_msg(_("Cross Validation was success!"));
-	    exit(EXIT_SUCCESS);
-	}
-
-	/*G_debug (0, _("passoE = %f, passoN=%f, lambda_i=%f"), passoE, passoN, lambda); */
-    }
-
-/*-------------------------------------------------------------------------------------------*/
+    /*----------------------------------------------------------------*/
     /* Interpolation begins */
     G_debug(1, "Interpolation()");
 
-    if (map) {
-	nrows = G_window_rows();
-	ncols = G_window_cols();
-
-	if (!G_alloc_matrix(nrows, ncols))
-	    G_fatal_error(_("Interpolation: The region resolution is too high: "
-			   "%d cells. Consider to change it."),
-			  nrows * ncols);
-    }
-
     /* Open driver and database */
     driver = db_start_driver_open_database(dvr, db);
     if (driver == NULL)
 	G_fatal_error(_("No database connection for driver <%s> is defined. "
 			"Run db.connect."), dvr);
 
+    /* Create auxiliar table */
+    if (vector) {
+	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 */
     G_debug(1, "Interpolation: Setting regions and boxes");
     G_get_window(&original_reg);
@@ -379,49 +382,66 @@
     Vect_region_box(&elaboration_reg, &overlap_box);
     Vect_region_box(&elaboration_reg, &general_box);
 
+    nrows = G_window_rows();
+    ncols = G_window_cols();
+
     /* Alloc raster matrix */
-    if (map)
+    if (grid == TRUE) {
 	if (!(raster_matrix = G_alloc_matrix(nrows, ncols)))
 	    G_fatal_error(_("Cannot allocate memory for auxiliar matrix."
 			    "Consider changing region resolution"));
+    }
 
-    /* Fixxing parameters of the elaboration region */
-    P_zero_dim(&dims);		/* Set to zero the dim struct */
-    dims.latoE = NSPLX_MAX * passoE;
-    dims.latoN = NSPLY_MAX * passoN;
-    dims.overlap = OVERLAP_SIZE * passoE;
-    P_get_orlo(bilin, &dims, passoE, passoN);	/* Set the last two dim elements */
+    /*------------------------------------------------------------------
+      | 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.ns_res * original_reg.rows;
-    E_extension = original_reg.ew_res * original_reg.cols;
-    orloE = dims.latoE - dims.overlap - 2 * dims.orlo_h;
-    orloN = dims.latoN - dims.overlap - 2 * dims.orlo_v;
-    nsubregion_col = 2 + ((E_extension - dims.latoE) / orloE);
-    nsubregion_row = 2 + ((N_extension - dims.latoN) / orloN);
+    /* Fixing parameters of the elaboration region */
+    P_zero_dim(&dims);		/* Set dim struct to zero */
 
+    nsplx_adj = NSPLX_MAX;
+    nsply_adj = NSPLY_MAX;
+    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 subregions */
+    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;
+
     if (nsubregion_col < 0)
 	nsubregion_col = 0;
     if (nsubregion_row < 0)
 	nsubregion_row = 0;
 
+    nsubregions = nsubregion_row * nsubregion_col;
+
     /* Creating line and categories structs */
     points = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
-    nlines = Vect_get_num_lines(&In);
+    Vect_cat_set(Cats, 1, 0);
 
-    /*------------------------------------------------------------------------------------------
-      | 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 resolution
-      ------------------------------------------------------------------------------------------*/
-
-    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 row */
+
+    while (last_row == FALSE) {	/* For each subregion row */
 	subregion_row++;
 	P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 		      GENERAL_ROW);
@@ -430,36 +450,35 @@
 
 	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 			  FIRST_ROW);
-	    nsply =
-		ceil((elaboration_reg.north -
-		      elaboration_reg.south) / passoN);
-	    G_debug(1, "Interpolation: nsply = %d", nsply);
-	    if (nsply > NSPLY_MAX)
-		nsply = NSPLY_MAX;
 	}
 
 	if (elaboration_reg.south <= original_reg.south) {	/* Last row */
 
 	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 			  LAST_ROW);
-	    nsply =
-		ceil((elaboration_reg.north -
-		      elaboration_reg.south) / passoN);
 	    last_row = TRUE;
-	    G_debug(1, "Interpolation: nsply = %d", nsply);
-	    if (nsply > NSPLY_MAX)
-		nsply = NSPLY_MAX;
 	}
 
+	nsply =
+	    ceil((elaboration_reg.north -
+		  elaboration_reg.south) / passoN) + 0.5;
+	G_debug(1, "Interpolation: nsply = %d", nsply);
+	/*
+	if (nsply > NSPLY_MAX)
+	    nsply = NSPLY_MAX;
+	*/
 	elaboration_reg.east = original_reg.west;
 	last_column = FALSE;
-
 	subregion_col = 0;
-	while (last_column == FALSE) {	/* For each column */
+
+	while (last_column == FALSE) {	/* For each subregion column */
 	    int npoints = 0;
-	    int nsubzones = 0, subzone = 0;
 
 	    subregion_col++;
+	    subregion++;
+	    if (nsubregions > 1)
+		G_message(_("subregion %d of %d"), subregion, nsubregions);
+
 	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 			  GENERAL_COLUMN);
 
@@ -467,12 +486,6 @@
 
 		P_set_regions(&elaboration_reg, &general_box, &overlap_box,
 			      dims, FIRST_COLUMN);
-		nsplx =
-		    ceil((elaboration_reg.east -
-			  elaboration_reg.west) / passoE);
-		G_debug(1, "Interpolation: nsply = %d", nsply);
-		if (nsplx > NSPLX_MAX)
-		    nsplx = NSPLX_MAX;
 	    }
 
 	    if (elaboration_reg.east >= original_reg.east) {	/* Last column */
@@ -480,13 +493,15 @@
 		P_set_regions(&elaboration_reg, &general_box, &overlap_box,
 			      dims, LAST_COLUMN);
 		last_column = TRUE;
-		nsplx =
-		    ceil((elaboration_reg.east -
-			  elaboration_reg.west) / passoE);
-		G_debug(1, "Interpolation: nsply = %d", nsply);
-		if (nsplx > NSPLX_MAX)
-		    nsplx = NSPLX_MAX;
 	    }
+	    nsplx =
+		ceil((elaboration_reg.east -
+		      elaboration_reg.west) / passoE) + 0.5;
+	    G_debug(1, "Interpolation: nsplx = %d", nsplx);
+	    /*
+	    if (nsplx > NSPLX_MAX)
+		nsplx = NSPLX_MAX;
+	    */
 	    G_debug(1, "Interpolation: (%d,%d): subregion bounds",
 		    subregion_row, subregion_col);
 	    G_debug(1, "Interpolation: \t\tNORTH:%.2f\t",
@@ -496,7 +511,7 @@
 	    G_debug(1, "Interpolation: \t\tSOUTH:%.2f",
 		    elaboration_reg.south);
 
-	    /*Setting the active region */
+	    /* reading points in interpolation region */
 	    dim_vect = nsplx * nsply;
 	    observ =
 		P_Read_Vector_Region_Map(&In, &elaboration_reg, &npoints,
@@ -507,23 +522,18 @@
 
 	    if (npoints > 0) {	/*  */
 		int i;
-		double *obs_mean;
 
 		nparameters = nsplx * nsply;
-		BW = P_get_BandWidth(bilin, nsply);
+		BW = P_get_BandWidth(P_BILINEAR, nsply) + 2 * (bilin == P_BICUBIC);
 
-		/*Least Squares system */
+		/* Least Squares system */
 		N = G_alloc_matrix(nparameters, BW);	/* Normal matrix */
 		TN = G_alloc_vector(nparameters);	/* vector */
 		parVect = G_alloc_vector(nparameters);	/* Parameters vector */
 		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);
 
-		if (bspline_field <= 0)
-		    mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
-
 		for (i = 0; i < npoints; i++) {	/* Setting obsVect vector & Q matrix */
 		    double dval;
 
@@ -532,11 +542,10 @@
 		    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, type;
+			int cat, ival, ret;
 
-			if (!(type & GV_POINTS))
-			    continue;
 			cat = observ[i].cat;
 			if (cat < 0)
 			    continue;
@@ -546,14 +555,14 @@
 				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)"),
@@ -561,26 +570,25 @@
 			    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);
+		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);
@@ -590,6 +598,7 @@
 				   nparameters, BW);
 		    nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
 		}
+		/* Bicubic interpolation */
 		else {
 		    G_debug(1,
 			    "Interpolation: (%d,%d): Bicubic interpolation...",
@@ -607,22 +616,16 @@
 		G_free_vector(TN);
 		G_free_vector(Q);
 
-		if (grid == FALSE) {	/*OBSERVATION POINTS INTERPOLATION */
-		    /* Auxiliar table creation */
-		    if (flag_auxiliar == FALSE) {
-			G_debug(1,
-				"Interpolation: Creating auxiliar table for archiving "
-				"overlapping zones");
-			if ((flag_auxiliar =
-			     P_Create_Aux_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);
-			}
-		    }
-
+		if (grid == TRUE) {	/* GRID INTERPOLATION ==> INTERPOLATION INTO A RASTER */
+		    G_debug(1, "Interpolation: (%d,%d): Regular_Points...",
+			    subregion_row, subregion_col);
+		    raster_matrix =
+			P_Regular_Points(&elaboration_reg, general_box,
+					 overlap_box, raster_matrix, parVect,
+					 passoN, passoE, dims.overlap, mean,
+					 nsplx, nsply, nrows, ncols, bilin);
+		}
+		else {		/* OBSERVATION POINTS INTERPOLATION */
 		    if (ext == FALSE) {
 			G_debug(1, "Interpolation: (%d,%d): Sparse_Points...",
 				subregion_row, subregion_col);
@@ -632,11 +635,7 @@
 					dims.overlap, nsplx, nsply, npoints,
 					bilin, Cats, driver, mean,
 					table_name);
-
-			G_free_matrix(obsVect);
-			G_free_vector(parVect);
 		    }
-
 		    else {	/* FLAG_EXT == TRUE */
 			int npoints_ext, *lineVect_ext = NULL;
 			double **obsVect_ext;	/*, mean_ext = .0; */
@@ -670,53 +669,38 @@
 					mean, table_name);
 
 			G_free_matrix(obsVect_ext);
-			G_free_vector(parVect);
 			G_free_ivector(lineVect_ext);
 		    }		/* END FLAG_EXT == TRUE */
-		}		/* END IF GRID == FLASE */
-
-		else {		/*GRID INTERPOLATION ==> INTERPOLATION INTO A RASTER */
-		    G_free_matrix(obsVect);
-		    flag_auxiliar = TRUE;
-		    G_debug(1, "Interpolation: (%d,%d): Regular_Points...",
-			    subregion_row, subregion_col);
-		    raster_matrix =
-			P_Regular_Points(&elaboration_reg, general_box,
-					 overlap_box, raster_matrix, parVect,
-					 passoN, passoE, dims.overlap, mean,
-					 nsplx, nsply, nrows, ncols, bilin);
-		    G_free_vector(parVect);
-		}
+		}		/* END GRID == FALSE */
+		G_free_vector(parVect);
+		G_free_matrix(obsVect);
 		G_free_ivector(lineVect);
 	    }
-	    else
-		G_warning(_("No data within this subzone. "
+	    else {
+		G_free(observ);
+		G_warning(_("No data within this subregion. "
 			    "Consider changing the spline step."));
-
-	    subzone = (subregion_row - 1) * nsubregion_col + subregion_col;
-	    nsubzones = nsubregion_row * nsubregion_col;
-	    G_percent(subzone, nsubzones, 10);
-
+	    }
 	}			/*! END WHILE; last_column = TRUE */
     }				/*! END WHILE; last_row = TRUE */
 
-    /* Writing into the output vector map the points from the overlapping zones */
-    if (flag_auxiliar == TRUE) {
-	if (grid == FALSE) {
-	    if (ext == FALSE)
-		P_Aux_to_Vector(&In, &Out, driver, table_name);
-	    else
-		P_Aux_to_Vector(&In_ext, &Out, driver, table_name);
+    G_verbose_message(_("Writing output..."));
+    /* Writing the output raster map */
+    if (grid == TRUE) {
+	P_Aux_to_Raster(raster_matrix, raster);
+	G_free_matrix(raster_matrix);
+    }
+    /* Writing to the output vector map the points from the overlapping zones */
+    else if (flag_auxiliar == TRUE) {
+	if (ext == FALSE)
+	    P_Aux_to_Vector(&In, &Out, driver, table_name);
+	else
+	    P_Aux_to_Vector(&In_ext, &Out, driver, table_name);
 
-	    /* Dropping auxiliar table */
-	    G_debug(1, "%s: Dropping <%s>", argv[0], table_name);
-	    if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
-		G_fatal_error(_("Auxiliar table could not be dropped"));
-	}
-	else {
-	    P_Aux_to_Raster(raster_matrix, raster);
-	    G_free_matrix(raster_matrix);
-	}
+	/* Dropping auxiliar table */
+	G_debug(1, "%s: Dropping <%s>", argv[0], table_name);
+	if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
+	    G_fatal_error(_("Auxiliar table could not be dropped"));
     }
 
     db_close_database_shutdown_driver(driver);
@@ -732,7 +716,7 @@
 
 	/* set map title */
 	sprintf(title, "%s interpolation with Tykhonov regularization",
-		type->answer);
+		type_opt->answer);
 	G_put_cell_title(out_map_opt->answer, title);
 	/* write map history */
 	G_short_history(out_map_opt->answer, "raster", &history);



More information about the grass-commit mailing list