[GRASS-SVN] r40831 - in grass/trunk/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
Fri Feb 5 13:51:51 EST 2010


Author: mmetz
Date: 2010-02-05 13:51:50 -0500 (Fri, 05 Feb 2010)
New Revision: 40831

Modified:
   grass/trunk/vector/lidar/lidarlib/PolimiFunct.h
   grass/trunk/vector/lidar/lidarlib/raster.c
   grass/trunk/vector/lidar/lidarlib/zones.c
   grass/trunk/vector/lidar/v.lidar.correction/correction.c
   grass/trunk/vector/lidar/v.lidar.correction/correction.h
   grass/trunk/vector/lidar/v.lidar.correction/main.c
   grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.c
   grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.h
   grass/trunk/vector/lidar/v.lidar.edgedetection/main.c
   grass/trunk/vector/lidar/v.lidar.growing/main.c
   grass/trunk/vector/lidar/v.outlier/main.c
   grass/trunk/vector/lidar/v.outlier/outlier.c
   grass/trunk/vector/lidar/v.outlier/outlier.h
   grass/trunk/vector/lidar/v.surf.bspline/main.c
Log:
lidar tools overhaul

Modified: grass/trunk/vector/lidar/lidarlib/PolimiFunct.h
===================================================================
--- grass/trunk/vector/lidar/lidarlib/PolimiFunct.h	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/lidarlib/PolimiFunct.h	2010-02-05 18:51:50 UTC (rev 40831)
@@ -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 		     2000	/* Side's size for v.to.qrast command. */
+#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 *, /**/
 		  struct 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 /**/);

Modified: grass/trunk/vector/lidar/lidarlib/raster.c
===================================================================
--- grass/trunk/vector/lidar/lidarlib/raster.c	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/lidarlib/raster.c	2010-02-05 18:51:50 UTC (rev 40831)
@@ -22,14 +22,11 @@
 {
     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();
 
     for (i = 0; i < num_points; i++) {
@@ -208,13 +205,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 = Rast_col_to_easting((double)(col) + 0.5, &Original);
 	    Y = Rast_row_to_northing((double)(row) + 0.5, &Original);
 
@@ -238,7 +258,6 @@
 
 		}
 		else {
-
 		    if ((X > Overlap.E)) {
 
 			if ((Y > Overlap.N)) {	/* (3) */

Modified: grass/trunk/vector/lidar/lidarlib/zones.c
===================================================================
--- grass/trunk/vector/lidar/lidarlib/zones.c	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/lidarlib/zones.c	2010-02-05 18:51:50 UTC (rev 40831)
@@ -27,7 +27,7 @@
 P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
 	      struct 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 region limits-> Also set the limits of the orlo and overlapping regions->
      * Returns 0 on success; -1 on failure*/
     struct Cell_head orig;
 
@@ -55,33 +55,33 @@
 	return 0;
 
     case FIRST_ROW:		/* It is just started with first row */
-	Elaboration->north = orig.north;
+	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;
+	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;
+	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;
+	Elaboration->east = orig.east + 2 * dim.orlo_v;
+	General->E = Elaboration->east - 2 * dim.orlo_v;
+	Overlap->E = General->E;
 	return 0;
     }
 
@@ -89,18 +89,98 @@
 }
 
 /*----------------------------------------------------------------------------------------*/
+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;
+	dim->orlo_v = 9 * pe;	/*4 */
+	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;
+	dim->orlo_v = 12 * pe;	/*3 */
+	dim->orlo_h = 12 * pn;
 	return 2;
     }
     else
@@ -148,12 +228,71 @@
 
     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;
+    struct bound_box region_box;
+    struct Cell_head orig;
+
+    G_get_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;
+
+	/* Reading and storing points only if it is in elaboration_reg */
+	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) {
+	*dist = sqrt(((xmax - xmin) * (ymax - ymin)) / npoints);
+	*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;
@@ -173,8 +312,11 @@
     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];
@@ -287,6 +429,7 @@
 	}
 	Rast_put_d_row(fd, raster);
     }
+    G_percent(row, nrows, 2);
 }
 
 /*------------------------------------------------------------------------------------------------*/
@@ -295,7 +438,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;

Modified: grass/trunk/vector/lidar/v.lidar.correction/correction.c
===================================================================
--- grass/trunk/vector/lidar/v.lidar.correction/correction.c	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.lidar.correction/correction.c	2010-02-05 18:51:50 UTC (rev 40831)
@@ -34,7 +34,7 @@
 		    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)
+		    dbDriver * driver, double mean, char *tab_name)
 {
     int i = 0, class;
     double interpolation, csi, eta, weight;
@@ -81,11 +81,11 @@
 			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"));
 
 		    }
@@ -96,7 +96,7 @@
 
 			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"));
 
 		    }
@@ -105,7 +105,7 @@
 
 			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"));
 		    }
 		}
@@ -117,7 +117,7 @@
 
 			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);
@@ -139,11 +139,11 @@
 			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) */
@@ -151,7 +151,7 @@
 			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);
@@ -174,7 +174,7 @@
 			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);
@@ -194,7 +194,7 @@
 
 			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"));
 		    }
 		}
@@ -229,7 +229,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];
@@ -242,8 +242,7 @@
     db_init_string(&sql);
     db_zero_string(&sql);
 
-    sprintf(buf, "SELECT Interp FROM Auxiliar_correction_table WHERE ID=%d",
-	    line_num);
+    sprintf(buf, "SELECT 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)
@@ -261,35 +260,41 @@
 
 	*Interp += db_get_value_double(Interp_value);
     }
+    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,

Modified: grass/trunk/vector/lidar/v.lidar.correction/correction.h
===================================================================
--- grass/trunk/vector/lidar/v.lidar.correction/correction.h	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.lidar.correction/correction.h	2010-02-05 18:51:50 UTC (rev 40831)
@@ -23,19 +23,20 @@
 		    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 /**/);
 
-void P_Aux_to_Raster(double **, int);
 int correction(int, double, double, double, double);

Modified: grass/trunk/vector/lidar/v.lidar.correction/main.c
===================================================================
--- grass/trunk/vector/lidar/v.lidar.correction/main.c	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.lidar.correction/main.c	2010-02-05 18:51:50 UTC (rev 40831)
@@ -34,11 +34,16 @@
 int main(int argc, char *argv[])
 {
     /* Declarations */
-    int dim_vect, nparameters, BW, npoints, nrows, ncols, nsply, nsplx;
+    int dim_vect, nparameters, BW, npoints, nrows, ncols;
+    int nsply, nsplx, nsplx_adj, nsply_adj;
+    int nsubregion_col, nsubregion_row;
+    int subzone = 0, nsubzones = 0;
     const char *dvr, *db, *mapset;
-    char table_name[1024];
+    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;
 
@@ -51,6 +56,7 @@
     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;
@@ -69,6 +75,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 spline step value");
+    spline_step_flag->description =
+	_("Estimate a good spline step value 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)");
@@ -136,12 +148,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);
@@ -154,6 +186,20 @@
     if (1 > Vect_open_old(&In, in_opt->answer, mapset))
 	G_fatal_error(_("Unable to open vector map <%s>"), 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);
@@ -180,29 +226,58 @@
 	G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
 		      dvr);
 
+    /* Create auxiliar table */
+    if ((flag_auxiliar =
+	 P_Create_Aux_Table(driver, table_name)) == FALSE) {
+	Vect_close(&In);
+	Vect_close(&Out);
+	Vect_close(&Terrain);
+	exit(EXIT_FAILURE);
+    }
+
     /* 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 = Rast_window_rows();
     ncols = Rast_window_cols();
 
     ew_resol = original_reg.ew_res;
     ns_resol = original_reg.ns_res;
 
+    /* Fixing parameters of the elaboration region */
+    /*! Each original_region will be divided into several subregions. These
+     *  subregions will overlap with their neibourghing subregions. This overlapping
+     *  is calculated as OVERLAP_SIZE times the east-west spline step. */
+
     P_zero_dim(&dims);
-    dims.latoE = NSPLX_MAX * passoE;
-    dims.latoN = NSPLY_MAX * passoN;
-    dims.overlap = OVERLAP_SIZE * ew_resol;
+    N_extension = original_reg.north - original_reg.south;
+    E_extension = original_reg.east - original_reg.west;
+
+    nsplx_adj = NSPLX_MAX;
+    nsply_adj = NSPLY_MAX;
+    dims.overlap = OVERLAP_SIZE * passoE;
     P_get_orlo(P_BILINEAR, &dims, passoE, passoN);
+    P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
 
+    G_verbose_message("adjusted EW splines %d", nsplx_adj);
+    G_verbose_message("adjusted NS splines %d", nsply_adj);
+
+    orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
+    orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+
+    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;
+
+    nsubzones = nsubregion_row * nsubregion_col;
+
     /* Subdividing and working with tiles */
     elaboration_reg.south = original_reg.north;
     last_row = FALSE;
@@ -223,11 +298,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;
@@ -235,6 +312,10 @@
 
 	while (last_column == FALSE) {	/* For each column */
 
+	    subzone++;
+	    if (nsubzones > 1)
+		G_message(_("subzone %d of %d"), subzone, nsubzones);
+
 	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 			  GENERAL_COLUMN);
 
@@ -251,10 +332,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);
 
 	    dim_vect = nsplx * nsply;
@@ -306,27 +389,16 @@
 		G_free_vector(TN);
 		G_free_vector(Q);
 
-		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"));
 		P_Sparse_Correction(&In, &Out, &Terrain, &elaboration_reg,
 				    general_box, overlap_box, obsVect,
 				    parVect, lineVect, passoN, passoE,
 				    dims.overlap, HighThresh, LowThresh,
-				    nsplx, nsply, npoints, driver, mean);
+				    nsplx, nsply, npoints, driver, mean, table_name);
 
+		G_free_vector(parVect);
 		G_free_matrix(obsVect);
 		G_free_ivector(lineVect);
-
-		G_free_vector(parVect);
 	    }
 	    G_free(observ);
 	}			/*! END WHILE; last_column = TRUE */

Modified: grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.c
===================================================================
--- grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.c	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.c	2010-02-05 18:51:50 UTC (rev 40831)
@@ -33,6 +33,7 @@
 #include <grass/dbmi.h>
     /* #include <grass/PolimiFunct.h> */
 #include "edgedetection.h"
+
 int edge_detection(struct Cell_head elaboration_reg, struct bound_box Overlap_Box,
 		   double *parBilin, double obsX, double obsY,
 		   double *partial, double alpha, double residual,
@@ -43,14 +44,14 @@
     /* 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]);
+    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;
@@ -61,72 +62,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 + M_PI),
-			     obsY + passoN * sin(dirPto + M_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]) + M_PI / 2 - dirPto) < alpha) &&
 		(fabs(atan(g[7][1] / g[7][0]) + M_PI / 2 - dirPto) < alpha)) {
 
-		gradient =
-		    Get_Gradient(Elaboration,
+		Get_Gradient(Elaboration,
 				 obsX + passoE * cos(dirPto + M_PI / 4),
 				 obsY + passoN * sin(dirPto + M_PI / 4),
-				 parBilin);
+				 parBilin, gradient);
 		g[1][0] = gradient[0];
 		g[1][1] = gradient[1];
 
-		gradient =
-		    Get_Gradient(Elaboration,
+		Get_Gradient(Elaboration,
 				 obsX + passoE * cos(dirPto - M_PI / 4),
 				 obsY + passoN * sin(dirPto - M_PI / 4),
-				 parBilin);
+				 parBilin, gradient);
 		g[3][0] = gradient[0];
 		g[3][1] = gradient[1];
 
-		gradient =
-		    Get_Gradient(Elaboration,
+		Get_Gradient(Elaboration,
 				 obsX + passoE * cos(dirPto + M_PI / 2),
 				 obsY + passoN * sin(dirPto + M_PI / 2),
-				 parBilin);
+				 parBilin, gradient);
 		g[4][0] = gradient[0];
 		g[4][1] = gradient[1];
 
-		gradient =
-		    Get_Gradient(Elaboration,
+		Get_Gradient(Elaboration,
 				 obsX + passoE * cos(dirPto - M_PI / 2),
 				 obsY + passoN * sin(dirPto - M_PI / 2),
-				 parBilin);
+				 parBilin, gradient);
 		g[5][0] = gradient[0];
 		g[5][1] = gradient[1];
 
-		gradient =
-		    Get_Gradient(Elaboration,
+		Get_Gradient(Elaboration,
 				 obsX + passoE * cos(dirPto + M_PI * 3 / 4),
 				 obsY + passoN * sin(dirPto + M_PI * 3 / 4),
-				 parBilin);
+				 parBilin, gradient);
 		g[6][0] = gradient[0];
 		g[6][1] = gradient[1];
 
-		gradient =
-		    Get_Gradient(Elaboration,
+		Get_Gradient(Elaboration,
 				 obsX + passoE * cos(dirPto - M_PI * 3 / 4),
 				 obsY + passoN * sin(dirPto - M_PI * 3 / 4),
-				 parBilin);
+				 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++;
 
@@ -139,23 +132,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;
@@ -167,7 +158,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,
@@ -175,10 +166,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;
@@ -190,7 +181,7 @@
     categories = Vect_new_cats_struct();
 
     for (i = 0; i < num_points; i++) {	/* Sparse points */
-	G_percent(i, num_points, 6);
+	G_percent(i, num_points, 2);
 
 	Vect_reset_line(point);
 	Vect_reset_cats(categories);
@@ -202,10 +193,10 @@
 				       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); */
 
@@ -221,7 +212,7 @@
 		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++;
 
 	    }
@@ -237,15 +228,15 @@
 			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) */
@@ -253,21 +244,25 @@
 			eta = (*point->y - General.S) / overlap;
 			weight = (1 - 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;
 
-			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"));
 		    }
 
 		}
@@ -282,10 +277,10 @@
 			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 =
@@ -300,7 +295,7 @@
 				     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++;
 
 		    }
@@ -313,15 +308,15 @@
 			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) */
@@ -331,10 +326,10 @@
 			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 =
@@ -349,7 +344,7 @@
 				     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++;
 		    }
 
@@ -362,10 +357,10 @@
 			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 =
@@ -380,18 +375,20 @@
 				     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;
 
-			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 */
 	    }
@@ -403,61 +400,69 @@
     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, partialX, partialY)", 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, PartialX=%lf, PartialY=%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;
@@ -468,15 +473,16 @@
 
     db_init_string(&sql);
     sprintf(buf,
-	    "SELECT ID, Interp, partialX, partialY FROM Auxiliar_edge_table WHERE ID=%d",
-	    line_num);
+	    "SELECT ID, Interp, partialX, partialY 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);
@@ -505,18 +511,18 @@
 	*PartialY += db_get_value_double(PartialY_value);
     }
     db_close_cursor(&cursor);
+    db_free_string(&sql);
     return DB_OK;
 }
 
-
-int Create_AuxEdge_Table(dbDriver * driver)
+int Create_AuxEdge_Table(dbDriver * driver, char *tab_name)
 {
     dbTable *table;
     dbColumn *ID_col, *Interp_col, *PartialX_col, *PartialY_col;
     int created;
 
     table = db_alloc_table(4);
-    db_set_table_name(table, "Auxiliar_edge_table");
+    db_set_table_name(table, tab_name);
     db_set_table_description(table,
 			     "It is used for the intermediate interpolated and gradient values");
 
@@ -537,38 +543,24 @@
     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."));
+	G_debug(1, "<%s> table created in database.", db_get_table_name(table));
 	created = TRUE;
     }
     else
 	return FALSE;
 
+    db_create_index2(driver, tab_name, "ID");
+
     return created;
 }
 
-int Drop_Aux_Table(dbDriver * driver)
+int Create_Interpolation_Table(dbDriver * driver, char *tab_name)
 {
-    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_name(table, tab_name);
     db_set_table_description(table,
 			     "This table is the bicubic interpolation of the input vector");
 
@@ -581,11 +573,11 @@
     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));
+	G_debug(1, _("<%s> table created in database."), db_get_table_name(table));
 	return DB_OK;
     }
     else
-	return !DB_OK;
+	return DB_FAILED;
 }
 
 

Modified: grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.h
===================================================================
--- grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.h	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.h	2010-02-05 18:51:50 UTC (rev 40831)
@@ -38,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, /**/
@@ -53,16 +53,14 @@
 		    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 Create_AuxEdge_Table(dbDriver *, char *);
+int Create_Interpolation_Table(dbDriver *, char *);
 int Insert_Interpolation(double, int, dbDriver *, char *);
-int Drop_Aux_Table(dbDriver *);

Modified: grass/trunk/vector/lidar/v.lidar.edgedetection/main.c
===================================================================
--- grass/trunk/vector/lidar/v.lidar.edgedetection/main.c	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.lidar.edgedetection/main.c	2010-02-05 18:51:50 UTC (rev 40831)
@@ -33,22 +33,25 @@
 #include <grass/config.h>
 #include <grass/PolimiFunct.h>
 #include "edgedetection.h"
+
 int nsply, nsplx, line_out_counter, first_it;
 double passoN, passoE;
 
-
 /**************************************************************************************
 **************************************************************************************/
 int main(int argc, char *argv[])
 {
     /* Variables' declarations */
+    int nsplx_adj, nsply_adj;
+    int nsubregion_col, nsubregion_row, subzone = 0, nsubzones = 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;
+    double lambda_B, lambda_F, grad_H, grad_L, alpha, mean;
     const char *dvr, *db, *mapset;
-    char table_interpolation[1024], table_name[1024];
+    char table_interpolation[GNAME_MAX], table_name[GNAME_MAX];
+    char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
-    int last_row, last_column, flag_auxiliar = FALSE, pass_num;
+    int last_row, last_column, flag_auxiliar = FALSE;
 
     int *lineVect;
     double *TN, *Q, *parVect_bilin, *parVect_bicub;	/* Interpolating and least-square vectors */
@@ -58,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;
@@ -77,6 +81,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 spline step value");
+    spline_step_flag->description =
+	_("Estimate a good spline step value 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);
@@ -157,6 +167,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"));
 
@@ -164,8 +178,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);
@@ -174,15 +219,29 @@
 	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);
 
+    /* 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);
@@ -194,7 +253,11 @@
 	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 = Create_AuxEdge_Table(driver, table_name)) == FALSE)
+	G_fatal_error(_("It was impossible to create <%s>."), table_name);
+
+    if (Create_Interpolation_Table(driver, table_interpolation) != DB_OK)
 	G_fatal_error(_("It was impossible to create <%s> interpolation table in database."),
 		      out_opt->answer);
 
@@ -207,23 +270,39 @@
     /* Fixing parameters of the elaboration region */
     /*! Each original_region will be divided into several subregions. These
      *  subregion will be overlapped by its neighboring subregions. This overlapping
-     *  is calculated as OVERLAP_PASS times the east-west resolution. */
+     *  is calculated as OVERLAP_SIZE times the east-west spline step. */
 
-    ew_resol = original_reg.ew_res;
-    ns_resol = original_reg.ns_res;
-
     P_zero_dim(&dims);
 
-    dims.latoE = NSPLX_MAX * passoE;
-    dims.latoN = NSPLY_MAX * passoN;
-    dims.overlap = OVERLAP_SIZE * ew_resol;
-    P_get_orlo(P_BICUBIC, &dims, passoE, passoN);
+    N_extension = original_reg.north - original_reg.south;
+    E_extension = original_reg.east - original_reg.west;
 
+    nsplx_adj = NSPLX_MAX;
+    nsply_adj = NSPLY_MAX;
+    dims.overlap = OVERLAP_SIZE * passoE;
+    P_get_orlo(P_BICUBIC, &dims, passoE, passoN);	/* Set orlo_h|v */
+    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);
+
+    orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
+    orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+
+    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;
+
+    nsubzones = nsubregion_row * nsubregion_col;
+
     elaboration_reg.south = original_reg.north;
 
     last_row = FALSE;
     first_it = TRUE;
-    pass_num = 0;
 
     while (last_row == FALSE) {	/* For each row */
 
@@ -243,10 +322,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;
@@ -254,6 +335,10 @@
 
 	while (last_column == FALSE) {	/* For each column */
 
+	    subzone++;
+	    if (nsubzones > 1)
+		G_message(_("subzone %d of %d"), subzone, nsubzones);
+
 	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 			  GENERAL_COLUMN);
 
@@ -270,10 +355,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 */
@@ -286,9 +373,6 @@
 	    if (npoints > 0) {	/* If there is any point falling into elaboration_reg... */
 		int i, tn;
 
-		pass_num++;
-		G_verbose_message(_("Pass %d:"), pass_num);
-
 		nparameters = nsplx * nsply;
 
 		/* Mean's calculation */
@@ -345,20 +429,12 @@
 		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_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);
@@ -374,7 +450,7 @@
     /* 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"));
     }
 
@@ -382,7 +458,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/trunk/vector/lidar/v.lidar.growing/main.c
===================================================================
--- grass/trunk/vector/lidar/v.lidar.growing/main.c	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.lidar.growing/main.c	2010-02-05 18:51:50 UTC (rev 40831)
@@ -25,8 +25,8 @@
 #include <math.h>
 #include <grass/config.h>
 #include <grass/gis.h>
+#include <grass/vector.h>
 #include <grass/raster.h>
-#include <grass/vector.h>
 #include <grass/dbmi.h>
 #include <grass/glocale.h>
 #include "growing.h"
@@ -40,6 +40,8 @@
 
     /* Variables' declarations */
     int row, nrows, col, ncols, MaxPoints;
+    int nsubregion_col, nsubregion_row;
+    int subzone = 0, nsubzones = 0;
     int last_row, last_column;
     int nlines, nlines_first, line_num;
     int more;
@@ -48,7 +50,7 @@
     double Thres_j, Thres_d, ew_resol, ns_resol;
     double minNS, minEW, maxNS, maxEW;
     const char *mapset;
-    char buf[1024];
+    char buf[1024], table_name[GNAME_MAX];
     char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
     int colorBordo, ripieno, conta, lungPunti, lungHull, xi, c1, c2;
@@ -71,8 +73,6 @@
     dbString sql;
     dbTable *table;
     dbCursor cursor;
-    const char *p;
-    int found;
 
 /*------------------------------------------------------------------------------------------*/
     /* Options' declaration */ ;
@@ -124,27 +124,23 @@
     /* Open input vector */
     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);
 
@@ -175,11 +171,21 @@
     G_get_set_window(&original_reg);
     G_get_set_window(&elaboration_reg);
 
-    /*  Fixxing parameters of the elaboration region */
+    /*  Fixing parameters of the elaboration region */
     /*! The original_region will be divided into several subregions */
     ew_resol = original_reg.ew_res;
     ns_resol = original_reg.ns_res;
 
+    nsubregion_col = ceil((original_reg.east - original_reg.west) / (LATO * ew_resol)) + 0.5;
+    nsubregion_row = ceil((original_reg.north - original_reg.south) / (LATO * ns_resol)) + 0.5;
+
+    if (nsubregion_col < 0)
+	nsubregion_col = 0;
+    if (nsubregion_row < 0)
+	nsubregion_row = 0;
+
+    nsubzones = nsubregion_row * nsubregion_col;
+
     /* Subdividing and working with tiles */
     elaboration_reg.south = original_reg.north;
     last_row = FALSE;
@@ -187,22 +193,32 @@
     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);
+    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);
+    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);
 
     while (last_row == FALSE) {	/* For each strip of LATO rows */
 
@@ -223,6 +239,10 @@
 	while (last_column == FALSE) {	/* For each strip of LATO columns */
 	    struct bound_box elaboration_box;
 
+	    subzone++;
+	    if (nsubzones > 1)
+		G_message(_("subzone %d of %d"), subzone, nsubzones);
+
 	    elaboration_reg.west = elaboration_reg.east;
 	    if (elaboration_reg.west < original_reg.west)	/* First column */
 		elaboration_reg.west = original_reg.west;
@@ -234,10 +254,11 @@
 		last_column = TRUE;
 	    }
 
-	    /*Setting the active region */
-	    Rast_set_window(&elaboration_reg);
-	    nrows = Rast_window_rows();
-	    ncols = Rast_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;
 
 	    G_debug(1, _("Rows = %d"), nrows);
 	    G_debug(1, _("Columns = %d"), ncols);
@@ -282,6 +303,7 @@
 			(int)(Rast_easting_to_col
 			      (points->x[0], &elaboration_reg));
 
+		    Z_interp = 0;
 		    while (1) {
 			if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK ||
 			    !more)
@@ -356,13 +378,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)(Rast_northing_to_row

Modified: grass/trunk/vector/lidar/v.outlier/main.c
===================================================================
--- grass/trunk/vector/lidar/v.outlier/main.c	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.outlier/main.c	2010-02-05 18:51:50 UTC (rev 40831)
@@ -36,10 +36,15 @@
 int main(int argc, char *argv[])
 {
     /* Variables' declarations */
+    int nsplx_adj, nsply_adj;
+    int nsubregion_col, nsubregion_row;
+    int subzone = 0, nsubzones = 0;
+    double N_extension, E_extension, orloE, orloN;
     int dim_vect, nparameters, BW, npoints;
     double ew_resol, ns_resol, mean, lambda;
     const char *dvr, *db, *mapset;
-    char table_name[1024];
+    char table_name[GNAME_MAX];
+    char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
     int last_row, last_column, flag_auxiliar = FALSE;
 
@@ -51,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;
@@ -69,6 +74,12 @@
     G_add_keyword(_("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 spline step value");
+    spline_step_flag->description =
+	_("Estimate a good spline step value 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);
@@ -136,10 +147,6 @@
     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");
-
     /* Checking vector names */
     Vect_check_input_output_name(in_opt->answer, out_opt->answer,
 				 GV_FATAL_EXIT);
@@ -148,6 +155,45 @@
 	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);
+
+    /* 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))
@@ -165,12 +211,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);
@@ -192,12 +232,6 @@
 	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"));
-    }
-
     /* Setting regions and boxes */
     G_get_set_window(&original_reg);
     G_get_set_window(&elaboration_reg);
@@ -207,18 +241,38 @@
     /* 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. */
+     *  is calculated as OVERLAP_SIZE times the east-west spline step. */
 
     ew_resol = original_reg.ew_res;
     ns_resol = original_reg.ns_res;
 
     P_zero_dim(&dims);
 
-    dims.latoE = NSPLX_MAX * passoE;
-    dims.latoN = NSPLY_MAX * passoN;
-    dims.overlap = OVERLAP_SIZE * ew_resol;
-    P_get_orlo(P_BICUBIC, &dims, passoE, passoN);
+    N_extension = original_reg.north - original_reg.south;
+    E_extension = original_reg.east - original_reg.west;
 
+    nsplx_adj = NSPLX_MAX;
+    nsply_adj = NSPLY_MAX;
+    dims.overlap = OVERLAP_SIZE * passoE;
+    P_get_orlo(P_BILINEAR, &dims, passoE, passoN);
+    P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
+
+    G_verbose_message("adjusted EW splines %d", nsplx_adj);
+    G_verbose_message("adjusted NS splines %d", nsply_adj);
+
+    orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
+    orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+
+    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;
+
+    nsubzones = nsubregion_row * nsubregion_col;
+
     elaboration_reg.south = original_reg.north;
 
     last_row = FALSE;
@@ -242,10 +296,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;
@@ -253,6 +309,10 @@
 
 	while (last_column == FALSE) {	/* For each column */
 
+	    subzone++;
+	    if (nsubzones > 1)
+		G_message(_("subzone %d of %d"), subzone, nsubzones);
+
 	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 			  GENERAL_COLUMN);
 
@@ -269,10 +329,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 */
@@ -308,7 +370,7 @@
 		    Q[i] = 1;	/* Q=I */
 		}
 
-		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,
@@ -328,14 +390,15 @@
 			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);

Modified: grass/trunk/vector/lidar/v.outlier/outlier.c
===================================================================
--- grass/trunk/vector/lidar/v.outlier/outlier.c	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.outlier/outlier.c	2010-02-05 18:51:50 UTC (rev 40831)
@@ -8,18 +8,19 @@
 
 #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,
-	       struct bound_box General, struct bound_box Overlap, double **obs,
-	       double *parBilin, double mean, double overlap, int *line_num,
-	       int num_points, dbDriver * driver)
+	       struct bound_box General, struct 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;
-
+    double interpolation, weight, residual, eta, csi;
     extern int nsplx, nsply;
     extern double passoN, passoE;
-
     struct line_pnts *point;
     struct line_cats *categories;
 
@@ -27,6 +28,7 @@
     categories = Vect_new_cats_struct();
 
     for (i = 0; i < num_points; i++) {	/* Sparse points */
+	G_percent(i, num_points, 2);
 	Vect_reset_line(point);
 	Vect_reset_cats(categories);
 
@@ -65,16 +67,14 @@
 			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"));
 
 		    }
@@ -83,18 +83,20 @@
 			eta = (*point->y - General.S) / overlap;
 			weight = (1 - 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;
 
-			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"));
 		    }
 
@@ -106,12 +108,10 @@
 			eta = (*point->y - Overlap.N) / overlap;
 			weight = (1 - 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;
@@ -133,28 +133,24 @@
 			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;
 
-			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;
@@ -176,12 +172,10 @@
 		    if ((*point->y > Overlap.N) && (*point->y != General.N)) {	/*(3) */
 			weight = (*point->y - Overlap.N) / 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;
@@ -201,9 +195,10 @@
 		    else if ((*point->y < Overlap.S) && (*point->y != General.S)) {	/*(1) */
 			weight = (Overlap.S - *point->y) / 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) */
@@ -211,42 +206,54 @@
 	    }
 	}			/*end if obs */
     }				/*end for */
+
+    G_percent(num_points, num_points, 2);
+    G_debug(2, "P_outlier: done");
+
     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;
@@ -256,8 +263,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)
@@ -277,13 +283,12 @@
 	*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;
 

Modified: grass/trunk/vector/lidar/v.outlier/outlier.h
===================================================================
--- grass/trunk/vector/lidar/v.outlier/outlier.h	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.outlier/outlier.h	2010-02-05 18:51:50 UTC (rev 40831)
@@ -16,12 +16,12 @@
 	  struct 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/trunk/vector/lidar/v.surf.bspline/main.c
===================================================================
--- grass/trunk/vector/lidar/v.surf.bspline/main.c	2010-02-05 10:03:42 UTC (rev 40830)
+++ grass/trunk/vector/lidar/v.surf.bspline/main.c	2010-02-05 18:51:50 UTC (rev 40831)
@@ -1,5 +1,5 @@
 
-/***********************************************************************
+/**********************************************************************
  *
  * MODULE:       v.surf.bspline
  *
@@ -15,9 +15,9 @@
  *               Read the file COPYING that comes with GRASS
  *               for details.
  *
- **************************************************************************/
+ **********************************************************************/
 
- /*INCLUDES*/
+/* INCLUDES */
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
@@ -29,23 +29,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, nlines, nrows, ncols, nsplx_adj, nsply_adj;
     int nsubregion_col, nsubregion_row, subregion_row, subregion_col;
+    int subzone = 0, nsubzones = 0;
     int last_row, last_column, grid, bilin, ext, flag_auxiliar, cross;	/* booleans */
     double passoN, passoE, lambda, mean;
     double N_extension, E_extension, orloE, orloN;
 
     const char *mapset, *dvr, *db, *vector, *map;
-    char table_name[1024], title[64];
+    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 */
@@ -53,15 +55,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;
@@ -76,29 +78,32 @@
     struct field_info *Fi;
     dbDriver *driver, *driver_cats;
 
-    /*-------------------------------------------------------------------------------------------*/
-    /* Options' declaration */
-    module = G_define_module(); {
-G_add_keyword(_("vector"));
-G_add_keyword(_("interpolation"));
-	module->description =
-	    _("Bicubic or bilinear spline interpolation with Tykhonov regularization.");
-    }
+    /*----------------------------------------------------------------*/
+    /* Options declarations */
+    module = G_define_module();
+    G_add_keyword(_("vector"));
+    G_add_keyword(_("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;
@@ -107,68 +112,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 = "type";
+    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 = 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");
 
-    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;
@@ -180,10 +188,6 @@
     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"));
@@ -192,39 +196,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);
@@ -234,6 +225,41 @@
 	G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
 		      in_opt->answer);
 
+    if (!Vect_is_3d(&In) && bspline_field <= 0 && bspline_column)
+	G_fatal_error(_("Need either 3D vector or layer and column with z values"));
+
+    /* 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;
@@ -287,11 +313,6 @@
     Rast_set_fp_type(DCELL_TYPE);
     if (!vector && map) {
 	grid = TRUE;
-	/*
-	   if (G_find_raster (out_map_opt->answer, G_mapset()) != NULL) 
-	   G_fatal_error (_("Raster <%s> already exist."), out_map_opt->answer);
-	 */
-
 	raster = Rast_open_fp_new(out_map_opt->answer);
     }
 
@@ -325,54 +346,26 @@
 	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)
-		Rast_close(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 = Rast_window_rows();
-	ncols = Rast_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);
 
+    /* Auxiliar table creation */
+    if (vector) {
+	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);
+	}
+    }
+
     /* Setting regions and boxes */
     G_debug(1, "Interpolation: Setting regions and boxes");
     G_get_window(&original_reg);
@@ -380,49 +373,64 @@
     Vect_region_box(&elaboration_reg, &overlap_box);
     Vect_region_box(&elaboration_reg, &general_box);
 
+    nrows = Rast_window_rows();
+    ncols = Rast_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 */
+    /* Fixing 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;
+
+    N_extension = original_reg.north - original_reg.south;
+    E_extension = original_reg.east - original_reg.west;
+
+    nsplx_adj = NSPLX_MAX;
+    nsply_adj = NSPLY_MAX;
     dims.overlap = OVERLAP_SIZE * passoE;
-    P_get_orlo(bilin, &dims, passoE, passoN);	/* Set the last two dim elements */
+    P_get_orlo(bilin, &dims, passoE, passoN);	/* Set orlo_h|v */
+    P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
 
-    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);
+    G_verbose_message("adjusted EW splines %d", nsplx_adj);
+    G_verbose_message("adjusted NS splines %d", nsply_adj);
 
+    orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
+    orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+
+    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;
 
+    nsubzones = nsubregion_row * nsubregion_col;
+
     /* Creating line and categories structs */
     points = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
+    Vect_cat_set(Cats, 1, 0);
     nlines = Vect_get_num_lines(&In);
 
-    /*------------------------------------------------------------------------------------------
+    /*------------------------------------------------------------------
       | Subdividing and working with tiles: 									
       | Each original_region will be divided into several subregions. 
       | Each one will be overlaped by its neibourgh subregions. 
-      | The overlapping was calculated as a fixed OVERLAP_SIZE times the east-west resolution
-      ------------------------------------------------------------------------------------------*/
+      | The overlapping was calculated as a fixed OVERLAP_SIZE times
+      | the east-west spline step
+      ----------------------------------------------------------------*/
 
     elaboration_reg.south = original_reg.north;
 
     G_percent(0, 1, 10);
     subregion_row = 0;
     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);
@@ -431,36 +439,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++;
+	    subzone++;
+	    if (nsubzones > 1)
+		G_message(_("subzone %d of %d"), subzone, nsubzones);
+
 	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 			  GENERAL_COLUMN);
 
@@ -468,12 +475,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 */
@@ -481,13 +482,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",
@@ -497,7 +500,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,
@@ -511,9 +514,9 @@
 		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 */
@@ -522,9 +525,6 @@
 		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;
 
@@ -534,10 +534,8 @@
 		    obsVect[i][1] = observ[i].coordY;
 
 		    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;
@@ -566,13 +564,16 @@
 		    else {
 			obsVect[i][2] = observ[i].coordZ;
 			obs_mean[i] = observ[i].coordZ;
-		    }		/*obsVect[i][2] = observ[i].coordZ - mean; */
+		    }		/* obsVect[i][2] = observ[i].coordZ - mean; */
 		}
 
 		/* Mean calculation for every point */
 		if (bspline_field > 0)
 		    mean = calc_mean(obs_mean, npoints);
+		else
+		    mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
 
+
 		G_debug(1, "Interpolation: (%d,%d): mean=%lf",
 			subregion_row, subregion_col, mean);
 
@@ -608,22 +609,17 @@
 		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 */
+		    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);
+		}
+		else {		/* OBSERVATION POINTS INTERPOLATION */
 		    if (ext == FALSE) {
 			G_debug(1, "Interpolation: (%d,%d): Sparse_Points...",
 				subregion_row, subregion_col);
@@ -633,11 +629,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; */
@@ -671,53 +663,37 @@
 					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);
+		G_free_vector(obs_mean);
 	    }
 	    else
 		G_warning(_("No data within this subzone. "
 			    "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);
@@ -733,7 +709,7 @@
 
 	/* set map title */
 	sprintf(title, "%s interpolation with Tykhonov regularization",
-		type->answer);
+		type_opt->answer);
 	Rast_put_cell_title(out_map_opt->answer, title);
 	/* write map history */
 	Rast_short_history(out_map_opt->answer, "raster", &history);



More information about the grass-commit mailing list