lidar tools update

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

vector/lidar/Makefile

@@ -5,7 +5,8 @@ SUBDIRS1 = \
 	v.outlier \
 	v.lidar.correction \
 	v.lidar.edgedetection \
-	v.lidar.growing
+	v.lidar.growing \
+        r.resamp.bspline
 
 SUBDIRS = lidarlib $(SUBDIRS1)
 

vector/lidar/lidarlib/PolimiFunct.h

@@ -112,6 +112,11 @@ struct Point *P_Read_Vector_Region_Map(struct Map_info *, /**/
 				       struct Cell_head *, /**/
 				       int *, /**/ int, /**/ int /**/);
 
+struct Point *P_Read_Raster_Region_Map(double **, /**/
+				       struct Cell_head *, /**/
+				       struct Cell_head *, /**/
+				       int *, /**/ int *, /**/ int /**/);
+
 double P_Mean_Calc(struct Cell_head *, /**/ struct Point *, /**/ int /**/);
 
 /*----------------------------------------------------------------------------------------------------------*/
@@ -134,6 +139,7 @@ P_Sparse_Points(struct Map_info *, /**/
 		dbDriver *, /**/ double, /**/ char * /**/);
 
 double **P_Regular_Points(struct Cell_head *, /**/
+                          struct Cell_head *, /**/
 			  struct bound_box, /**/
 			  struct bound_box, /**/
 			  double **, /**/

vector/lidar/lidarlib/raster.c

@@ -79,7 +79,7 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
 			    G_fatal_error(_("Unable to access table <%s>"),
-					  buf);
+					  tab_name);
 		    }
 		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(1) */
 			csi = (General.E - *point->x) / overlap;
@@ -94,7 +94,7 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
 			    G_fatal_error(_("Unable to access table <%s>"),
-					  buf);
+					  tab_name);
 		    }
 		    else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) {	/*(1) */
 			weight = (General.E - *point->x) / overlap;
@@ -107,7 +107,7 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
 			    G_fatal_error(_("Unable to access table <%s>"),
-					  buf);
+					  tab_name);
 		    }
 		}
 		else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
@@ -124,7 +124,7 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
 			    G_fatal_error(_("Unable to access table <%s>"),
-					  buf);
+					  tab_name);
 		    }
 		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(2) */
 			csi = (*point->x - General.W) / overlap;
@@ -139,7 +139,7 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
 			    G_fatal_error(_("Unable to access table <%s>"),
-					  buf);
+					  tab_name);
 		    }
 		    else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) {	/*(2) */
 			weight = (*point->x - General.W) / overlap;
@@ -152,7 +152,7 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
 			    G_fatal_error(_("Unable to access table <%s>"),
-					  buf);
+					  tab_name);
 		    }
 		}
 		else if ((*point->x >= Overlap.W) && (*point->x <= Overlap.E)){
@@ -167,7 +167,7 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
 			    G_fatal_error(_("Unable to access table <%s>"),
-					  buf);
+					  tab_name);
 		    }
 		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(1) */
 			weight = (*point->y - General.S) / overlap;
@@ -180,7 +180,7 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
 			if (db_execute_immediate(driver, &sql) != DB_OK)
 			    G_fatal_error(_("Unable to access table <%s>"),
-					  buf);
+					  tab_name);
 		    }
 		}
 	    }
@@ -194,35 +194,35 @@ P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
 
 
 /*------------------------------------------------------------------------------------------------*/
-double **P_Regular_Points(struct Cell_head *Elaboration, struct bound_box General,
-			  struct bound_box Overlap, double **matrix, double *param,
+double **P_Regular_Points(struct Cell_head *Elaboration, struct Cell_head *Original,
+                          struct bound_box General, struct bound_box Overlap,
+			  double **matrix, double *param,
 			  double passoN, double passoE, double overlap,
-			  double mean, int nsplx, int nsply, int nrows,
-			  int ncols, int bilin)
+			  double mean, int nsplx, int nsply,
+			  int nrows, int ncols, int bilin)
 {
 
     int col, row, startcol, endcol, startrow, endrow;
     double X, Y, interpolation, weight, csi, eta;
-    struct Cell_head Original;
 
-    G_get_window(&Original);
-    if (Original.north > General.N)
-	startrow = (Original.north - General.N) / Original.ns_res -1;
+    /* G_get_window(&Original); */
+    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 (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;
+    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 (General.E > Original->west) {
+	endcol = (General.E - Original->west) / Original->ew_res + 1;
 	if (endcol > ncols)
 	    endcol = ncols;
     }
@@ -231,8 +231,8 @@ double **P_Regular_Points(struct Cell_head *Elaboration, struct bound_box Genera
 
     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);
+	    X = Rast_col_to_easting((double)(col) + 0.5, Original);
+	    Y = Rast_row_to_northing((double)(row) + 0.5, Original);
 
 	    if (Vect_point_in_box(X, Y, mean, &General)) {	/* Here, mean is just for asking if obs point is in box */
 

vector/lidar/lidarlib/zones.c

@@ -90,7 +90,7 @@ P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
     case FIRST_ROW:		/* Just started with first row */
 	Elaboration->north = orig.north + 2 * dim.edge_h;
 	Elaboration->south = Elaboration->north - dim.sn_size;
-	General->N = Elaboration->north - 2 * dim.edge_h;
+	General->N = orig.north;
 	General->S = Elaboration->south + dim.edge_h;
 	Overlap->N = General->N;
 	Overlap->S = General->S + dim.overlap;
@@ -98,14 +98,14 @@ P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
 
     case LAST_ROW:		/* Reached last row */
 	Elaboration->south = orig.south - 2 * dim.edge_h;
-	General->S = Elaboration->south + 2 * dim.edge_h;
+	General->S = orig.south;
 	Overlap->S = General->S;
 	return 0;
 
     case FIRST_COLUMN:		/* Just started with first column */
 	Elaboration->west = orig.west - 2 * dim.edge_v;
 	Elaboration->east = Elaboration->west + dim.ew_size;
-	General->W = Elaboration->west + 2 * dim.edge_v;
+	General->W = orig.west;
 	General->E = Elaboration->east - dim.edge_v;
 	Overlap->W = General->W;
 	Overlap->E = General->E - dim.overlap;
@@ -113,7 +113,7 @@ P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
 
     case LAST_COLUMN:		/* Reached last column */
 	Elaboration->east = orig.east + 2 * dim.edge_v;
-	General->E = Elaboration->east - 2 * dim.edge_v;
+	General->E = orig.east;
 	Overlap->E = General->E;
 	return 0;
     }
@@ -391,6 +391,87 @@ struct Point *P_Read_Vector_Region_Map(struct Map_info *Map,
     return obs;
 }
 
+struct Point *P_Read_Raster_Region_Map(double **matrix,
+				       struct Cell_head *Elaboration,
+				       struct Cell_head *Original,
+				       int *num_points, int *num_nulls, int dim_vect)
+{
+    int col, row, startcol, endcol, startrow, endrow, nrows, ncols;
+    int pippo, npoints, nnulls;
+    double x, y, z;
+    struct Point *obs;
+    struct bound_box elaboration_box;
+
+    pippo = dim_vect;
+    obs = (struct Point *)G_calloc(pippo, sizeof(struct Point));
+
+    /* Reading points inside elaboration zone */
+    Vect_region_box(Elaboration, &elaboration_box);
+
+    npoints = nnulls = 0;
+    nrows = Original->rows;
+    ncols = Original->cols;
+
+    if (Original->north > Elaboration->north)
+	startrow = (Original->north - Elaboration->north) / Original->ns_res - 1;
+    else
+	startrow = 0;
+    if (Original->north > Elaboration->south) {
+	endrow = (Original->north - Elaboration->south) / Original->ns_res + 1;
+	if (endrow > nrows)
+	    endrow = nrows;
+    }
+    else
+	endrow = nrows;
+    if (Elaboration->west > Original->west)
+	startcol = (Elaboration->west - Original->west) / Original->ew_res - 1;
+    else
+	startcol = 0;
+    if (Elaboration->east > Original->west) {
+	endcol = (Elaboration->east - 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++) {
+
+	    z = matrix[row][col];
+
+	    if (!Rast_is_d_null_value(&z)) {
+		x = Rast_col_to_easting((double)(col) + 0.5, Original);
+		y = Rast_row_to_northing((double)(row) + 0.5, Original);
+
+		if (Vect_point_in_box(x, y, 0, &elaboration_box)) {
+		    npoints++;
+		    if (npoints >= pippo) {
+			pippo += dim_vect;
+			obs =
+			    (struct Point *)G_realloc((void *)obs,
+						      (signed int)pippo *
+						      sizeof(struct Point));
+		    }
+
+		    /* Storing observation vector */
+		    obs[npoints - 1].coordX = x;
+		    obs[npoints - 1].coordY = y;
+		    obs[npoints - 1].coordZ = z;
+		}
+	    }
+	    else {
+		/* if point in output region */
+	        nnulls++;
+	    }
+	}
+    }
+
+    *num_points = npoints;
+    *num_nulls = nnulls;
+    return obs;
+}
+
 /*------------------------------------------------------------------------------------------------*/
 int P_Create_Aux2_Table(dbDriver * driver, char *tab_name)
 {

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

@@ -50,7 +50,7 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
     double *mean, *rms, *stdev;
 
     /* double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; */	/* Fixed values (by the moment) */
-    double lambda[PARAM_LAMBDA] = { 0.01, 0.05, 0.1, 0.2, 0.3, 0.4 };	/* Fixed values (by the moment) */
+    double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.005, 0.01, 0.02, 0.05 };	/* Fixed values (by the moment) */
 
     double *TN, *Q, *parVect;	/* Interpolation and least-square vectors */
     double **N, **obsVect;	/* Interpolation and least-square matrix */

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

@@ -378,7 +378,7 @@ int main(int argc, char *argv[])
 			    "It was impossible to create table <%s>."),
 			  table_name);
 	}
-	db_create_index2(driver, table_name, "ID");
+	/* db_create_index2(driver, table_name, "ID"); */
 	/* sqlite likes that */
 	db_close_database_shutdown_driver(driver);
 	driver = db_start_driver_open_database(dvr, db);
@@ -482,6 +482,10 @@ int main(int argc, char *argv[])
 
 	while (last_column == FALSE) {	/* For each subregion column */
 	    int npoints = 0;
+	    /* needed for sparse points interpolation */
+	    int npoints_ext, *lineVect_ext = NULL;
+	    double **obsVect_ext;	/*, mean_ext = .0; */
+	    struct Point *observ_ext;
 
 	    subregion_col++;
 	    subregion++;
@@ -522,14 +526,30 @@ int main(int argc, char *argv[])
 
 	    /* reading points in interpolation region */
 	    dim_vect = nsplx * nsply;
-	    observ =
-		P_Read_Vector_Region_Map(&In, &elaboration_reg, &npoints,
-					 dim_vect, bspline_field);
+	    if (grid == FALSE && ext == TRUE) {
+		observ_ext =
+		    P_Read_Vector_Region_Map(&In_ext,
+					     &elaboration_reg,
+					     &npoints_ext, dim_vect,
+					     1);
+	    }
+	    else
+		npoints_ext = 1;
+
+	    observ = NULL;
+	    if (npoints_ext > 0) {
+		observ =
+		    P_Read_Vector_Region_Map(&In, &elaboration_reg, &npoints,
+					     dim_vect, bspline_field);
+	    }
+	    else
+		npoints = 0;
+
 	    G_debug(1,
 		    "Interpolation: (%d,%d): Number of points in <elaboration_box> is %d",
 		    subregion_row, subregion_col, npoints);
 
-	    if (npoints > 0) {	/*  */
+	    if (npoints > 0 && npoints_ext > 0) {	/*  */
 		int i;
 
 		nparameters = nsplx * nsply;
@@ -629,7 +649,7 @@ int main(int argc, char *argv[])
 		    G_debug(1, "Interpolation: (%d,%d): Regular_Points...",
 			    subregion_row, subregion_col);
 		    raster_matrix =
-			P_Regular_Points(&elaboration_reg, general_box,
+			P_Regular_Points(&elaboration_reg, &original_reg, general_box,
 					 overlap_box, raster_matrix, parVect,
 					 stepN, stepE, dims.overlap, mean,
 					 nsplx, nsply, nrows, ncols, bilin);
@@ -646,8 +666,11 @@ int main(int argc, char *argv[])
 					table_name);
 		    }
 		    else {	/* FLAG_EXT == TRUE */
+
+			/* done that earlier */
+			/*
 			int npoints_ext, *lineVect_ext = NULL;
-			double **obsVect_ext;	/*, mean_ext = .0; */
+			double **obsVect_ext;
 			struct Point *observ_ext;
 
 			observ_ext =
@@ -655,6 +678,7 @@ int main(int argc, char *argv[])
 						     &elaboration_reg,
 						     &npoints_ext, dim_vect,
 						     1);
+			*/
 
 			obsVect_ext = G_alloc_matrix(npoints_ext, 3);	/* Observation vector_ext */
 			lineVect_ext = G_alloc_ivector(npoints_ext);
@@ -686,9 +710,11 @@ int main(int argc, char *argv[])
 		G_free_ivector(lineVect);
 	    }
 	    else {
-		G_free(observ);
-		G_warning(_("No data within this subregion. "
-			    "Consider changing the spline step."));
+		if (observ)
+		    G_free(observ);
+		if (npoints == 0)
+		    G_warning(_("No data within this subregion. "
+				"Consider changing the spline step."));
 	    }
 	}			/*! END WHILE; last_column = TRUE */
     }				/*! END WHILE; last_row = TRUE */