[GRASS-SVN] r42140 - in grass/trunk/vector/lidar: . lidarlib
v.surf.bspline
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri May 7 04:04:46 EDT 2010
Author: mmetz
Date: 2010-05-07 04:04:42 -0400 (Fri, 07 May 2010)
New Revision: 42140
Modified:
grass/trunk/vector/lidar/Makefile
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.surf.bspline/crosscorr.c
grass/trunk/vector/lidar/v.surf.bspline/main.c
Log:
lidar tools update
Modified: grass/trunk/vector/lidar/Makefile
===================================================================
--- grass/trunk/vector/lidar/Makefile 2010-05-07 08:03:03 UTC (rev 42139)
+++ grass/trunk/vector/lidar/Makefile 2010-05-07 08:04:42 UTC (rev 42140)
@@ -5,7 +5,8 @@
v.outlier \
v.lidar.correction \
v.lidar.edgedetection \
- v.lidar.growing
+ v.lidar.growing \
+ r.resamp.bspline
SUBDIRS = lidarlib $(SUBDIRS1)
Modified: grass/trunk/vector/lidar/lidarlib/PolimiFunct.h
===================================================================
--- grass/trunk/vector/lidar/lidarlib/PolimiFunct.h 2010-05-07 08:03:03 UTC (rev 42139)
+++ grass/trunk/vector/lidar/lidarlib/PolimiFunct.h 2010-05-07 08:04:42 UTC (rev 42140)
@@ -112,6 +112,11 @@
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 @@
dbDriver *, /**/ double, /**/ char * /**/);
double **P_Regular_Points(struct Cell_head *, /**/
+ struct Cell_head *, /**/
struct bound_box, /**/
struct bound_box, /**/
double **, /**/
Modified: grass/trunk/vector/lidar/lidarlib/raster.c
===================================================================
--- grass/trunk/vector/lidar/lidarlib/raster.c 2010-05-07 08:03:03 UTC (rev 42139)
+++ grass/trunk/vector/lidar/lidarlib/raster.c 2010-05-07 08:04:42 UTC (rev 42140)
@@ -79,7 +79,7 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
if (db_execute_immediate(driver, &sql) != DB_OK)
G_fatal_error(_("Unable to access table <%s>"),
- buf);
+ tab_name);
}
}
}
@@ -194,35 +194,35 @@
/*------------------------------------------------------------------------------------------------*/
-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 @@
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 */
Modified: grass/trunk/vector/lidar/lidarlib/zones.c
===================================================================
--- grass/trunk/vector/lidar/lidarlib/zones.c 2010-05-07 08:03:03 UTC (rev 42139)
+++ grass/trunk/vector/lidar/lidarlib/zones.c 2010-05-07 08:04:42 UTC (rev 42140)
@@ -90,7 +90,7 @@
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 @@
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 @@
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 @@
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)
{
Modified: grass/trunk/vector/lidar/v.surf.bspline/crosscorr.c
===================================================================
--- grass/trunk/vector/lidar/v.surf.bspline/crosscorr.c 2010-05-07 08:03:03 UTC (rev 42139)
+++ grass/trunk/vector/lidar/v.surf.bspline/crosscorr.c 2010-05-07 08:04:42 UTC (rev 42140)
@@ -50,7 +50,7 @@
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 */
Modified: grass/trunk/vector/lidar/v.surf.bspline/main.c
===================================================================
--- grass/trunk/vector/lidar/v.surf.bspline/main.c 2010-05-07 08:03:03 UTC (rev 42139)
+++ grass/trunk/vector/lidar/v.surf.bspline/main.c 2010-05-07 08:04:42 UTC (rev 42140)
@@ -378,7 +378,7 @@
"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 @@
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 @@
/* 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 @@
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 @@
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 @@
&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 @@
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 */
More information about the grass-commit
mailing list