[GRASS-SVN] r40860 - 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
Mon Feb 8 02:51:27 EST 2010
Author: mmetz
Date: 2010-02-08 02:51:25 -0500 (Mon, 08 Feb 2010)
New Revision: 40860
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.surf.bspline/crosscorr.c
grass/trunk/vector/lidar/v.surf.bspline/main.c
Log:
code clean-up and more documentation
Modified: grass/trunk/vector/lidar/lidarlib/PolimiFunct.h
===================================================================
--- grass/trunk/vector/lidar/lidarlib/PolimiFunct.h 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/lidarlib/PolimiFunct.h 2010-02-08 07:51:25 UTC (rev 40860)
@@ -33,7 +33,7 @@
#define NSPLX_MAX 150 /* Maximum number of splines along East direction used in the subregions interpolation */
#define NSPLY_MAX 150 /* Maximum number of splines along North direction used in the subregions interpolation */
#define OVERLAP_SIZE 10 /* Subregions overlapping size. */
-#define LATO 2000 /* Side's size for v.to.qrast command. */
+#define LATO 1000 /* Side's size for v.lidar.growing. */
#define CONTOUR 15 /**/
#define GENERAL_ROW 0
#define GENERAL_COLUMN 1
@@ -145,8 +145,10 @@
int, /**/ int, /**/ int, /**/ int, /**/ int /**/);
/*----------------------------------------------------------------------------------------------------------*/
-int P_Create_Aux_Table(dbDriver *, /**/ char * /**/);
+int P_Create_Aux2_Table(dbDriver *, /**/ char * /**/);
+int P_Create_Aux4_Table(dbDriver *, /**/ char * /**/);
+
int P_Drop_Aux_Table(dbDriver *, /**/ char * /**/);
/*----------------------------------------------------------------------------------------------------------*/
Modified: grass/trunk/vector/lidar/lidarlib/raster.c
===================================================================
--- grass/trunk/vector/lidar/lidarlib/raster.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/lidarlib/raster.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -29,6 +29,8 @@
point = Vect_new_line_struct();
+ db_begin_transaction(driver);
+
for (i = 0; i < num_points; i++) {
if (Vect_point_in_box(obs[i][0], obs[i][1], mean, &General)) { /*Here mean is just for asking if obs point is in box */
@@ -50,10 +52,8 @@
if (Vect_point_in_box(obs[i][0], obs[i][1], interpolation, &Overlap)) { /*(5) */
Vect_write_line(Out, GV_POINT, point, categories);
-
}
else {
-
db_init_string(&sql);
sprintf(buf, "INSERT INTO %s (ID, X, Y, Interp)", tab_name);
@@ -65,12 +65,11 @@
obs[i][1]);
db_append_string(&sql, buf);
- if ((point->x[0] > Overlap.E)) {
-
- if ((*point->y > Overlap.N)) { /*(3) */
- csi = (*point->x - Overlap.E) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - csi) * (1 - eta);
+ if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = csi * eta;
*point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
@@ -79,14 +78,13 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to write to table: %s"),
buf);
-
}
- else if ((*point->y < Overlap.S)) { /*(1) */
- csi = (*point->x - Overlap.E) / overlap;
- eta = (Overlap.S - *point->y) / overlap;
- weight = (1 - csi) * eta;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (*point->y - General.S) / overlap;
+ weight = csi * eta;
*point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
@@ -95,13 +93,12 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to write to table: %s"),
buf);
-
}
- else { /*(1) */
- weight = (*point->x - Overlap.E) / overlap;
- *point->z = (1 - weight) * interpolation;
+ else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(1) */
+ weight = (General.E - *point->x) / overlap;
+ *point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
db_append_string(&sql, buf);
@@ -109,16 +106,15 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to write to table: %s"),
buf);
}
-
}
- else if ((point->x[0] < Overlap.W)) {
- if ((*point->y > Overlap.N)) { /*(4) */
- csi = (Overlap.W - *point->x) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - eta) * csi;
+ else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(4) */
+ csi = (*point->x - General.W) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = eta * csi;
*point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
@@ -127,13 +123,12 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to write to table: %s"),
buf);
-
}
- else if ((*point->y < Overlap.S)) { /*(2) */
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(2) */
csi = (*point->x - General.W) / overlap;
- eta = (Overlap.S - *point->y) / overlap;
+ eta = (*point->y - General.S) / overlap;
weight = csi * eta;
*point->z = weight * interpolation;
@@ -143,13 +138,12 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to write to table: %s"),
buf);
-
}
- else { /*(2) */
- weight = (Overlap.W - *point->x) / overlap;
- *point->z = (1 - weight) * interpolation;
+ else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) { /*(2) */
+ weight = (*point->x - General.W) / overlap;
+ *point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
db_append_string(&sql, buf);
@@ -157,15 +151,14 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to write to table: %s"),
buf);
}
-
}
- else {
- if ((point->y[0] > Overlap.N)) { /*(3) */
- weight = (*point->y - Overlap.N) / overlap;
- *point->z = (1 - weight) * interpolation;
+ else if ((*point->x >= Overlap.W) && (*point->x <= Overlap.E)){
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ weight = (General.N - *point->y) / overlap;
+ *point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
db_append_string(&sql, buf);
@@ -173,11 +166,11 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to write to table: %s"),
buf);
}
- else { /*(1) */
- weight = (Overlap.S - *point->y) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ weight = (*point->y - General.S) / overlap;
*point->z = (1 - weight) * interpolation;
sprintf(buf, "%lf", *point->z);
@@ -186,14 +179,17 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to write to table: %s"),
buf);
}
}
}
- }
- /*IF*/}
- /*FOR*/ return;
+ } /*IF*/
+ } /*FOR*/
+
+ db_commit_transaction(driver);
+
+ return;
}
@@ -255,69 +251,59 @@
if (Vect_point_in_box(X, Y, interpolation, &Overlap)) { /* (5) */
matrix[row][col] = interpolation;
-
}
else {
- if ((X > Overlap.E)) {
-
- if ((Y > Overlap.N)) { /* (3) */
- csi = (X - Overlap.E) / overlap;
- eta = (Y - Overlap.N) / overlap;
- weight = (1 - csi) * (1 - eta);
+ if ((X > Overlap.E) && (X < General.E)) {
+ if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
+ csi = (General.E - X) / overlap;
+ eta = (General.N - Y) / overlap;
+ weight = csi * eta;
interpolation *= weight;
matrix[row][col] += interpolation;
-
}
- else if ((Y < Overlap.S)) { /* (1) */
- csi = (X - Overlap.E) / overlap;
+ else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
+ csi = (General.E - X) / overlap;
eta = (Y - General.S) / overlap;
- weight = (1 - csi) * eta;
+ weight = csi * eta;
interpolation *= weight;
matrix[row][col] = interpolation;
-
}
- else { /* (1) */
- weight = (X - Overlap.E) / overlap;
- interpolation *= 1 - weight;
+ else if ((Y >= Overlap.S) && (Y <= Overlap.N)) { /* (1) */
+ weight = (General.E - X ) / overlap;
+ interpolation *= weight;
matrix[row][col] = interpolation;
}
-
}
- else if ((X < Overlap.W)) {
-
- if ((Y > Overlap.N)) { /* (4) */
+ else if ((X < Overlap.W) && (X > General.W)) {
+ if ((Y > Overlap.N) && (Y < General.N)) { /* (4) */
csi = (X - General.W) / overlap;
- eta = (Y - Overlap.N) / overlap;
- weight = (1 - eta) * csi;
+ eta = (General.N - Y) / overlap;
+ weight = eta * csi;
interpolation *= weight;
matrix[row][col] += interpolation;
-
}
- else if ((Y < Overlap.S)) { /* (2) */
+ else if ((Y < Overlap.S) && (Y > General.S)) { /* (2) */
csi = (X - General.W) / overlap;
eta = (Y - General.S) / overlap;
weight = csi * eta;
interpolation *= weight;
matrix[row][col] += interpolation;
-
}
- else { /* (2) */
- weight = (Overlap.W - X) / overlap;
- interpolation *= 1 - weight;
+ else if ((Y >= Overlap.S) && (Y <= Overlap.N)) { /* (2) */
+ weight = (X - General.W) / overlap;
+ interpolation *= weight;
matrix[row][col] += interpolation;
}
-
}
- else {
- if ((Y > Overlap.N)) { /* (3) */
- weight = (Y - Overlap.N) / overlap;
- interpolation *= 1 - weight;
+ else if ((X >= Overlap.W) && (X <= Overlap.E)) {
+ if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
+ weight = (General.N - Y) / overlap;
+ interpolation *= weight;
matrix[row][col] += interpolation;
-
}
- else { /* (1) */
- weight = (Overlap.S - Y) / overlap;
- interpolation *= 1 - weight;
+ else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
+ weight = (Y - General.S) / overlap;
+ interpolation *= weight;
matrix[row][col] = interpolation;
}
}
Modified: grass/trunk/vector/lidar/lidarlib/zones.c
===================================================================
--- grass/trunk/vector/lidar/lidarlib/zones.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/lidarlib/zones.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -23,11 +23,44 @@
}
/*----------------------------------------------------------------------------------------*/
+/*
+ --------------------------------------------
+ | Elaboration region |
+ | ------------------------------------ |
+ | | General region | |
+ | | ---------------------------- | |
+ | | | | | |
+ | | | | | |
+ | | | | | |
+ | | | Overlap region | | |
+ | | | | | |
+ | | | | | |
+ | | | | | |
+ | | ---------------------------- | |
+ | | | |
+ | ------------------------------------ |
+ | |
+ --------------------------------------------
+
+ The terminology is misleading:
+ The Overlap region does NOT overlap with neighbouring segments,
+ but the Elaboration and General region do overlap
+
+ Elaboration is used for interpolation
+ Interpolated points in Elaboration but outside General are discarded
+ Interpolated points in General but outside Overlap are weighed by
+ their distance to Overlap and summed up
+ Interpolated points in Overlap are taken as they are
+
+ The buffer zones Elaboration - General and General - Overlap must be
+ large enough to avoid artefacts
+ */
+
int
P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
struct bound_box * Overlap, struct Reg_dimens dim, int type)
{
- /* Set the Elaboration region limits-> Also set the limits of the orlo and overlapping regions->
+ /* Set the Elaboration, General, and Overlap region limits
* Returns 0 on success; -1 on failure*/
struct Cell_head orig;
@@ -54,7 +87,7 @@
Overlap->E = General->E - dim.overlap;
return 0;
- case FIRST_ROW: /* It is just started with first row */
+ case FIRST_ROW: /* Just started with first row */
Elaboration->north = orig.north + 2 * dim.orlo_h;
Elaboration->south = Elaboration->north - dim.latoN;
General->N = Elaboration->north - 2 * dim.orlo_h;
@@ -63,13 +96,13 @@
Overlap->S = General->S + dim.overlap;
return 0;
- case LAST_ROW: /* It is reached last row */
+ case LAST_ROW: /* Reached last row */
Elaboration->south = orig.south - 2 * dim.orlo_h;
General->S = Elaboration->south + 2 * dim.orlo_h;
Overlap->S = General->S;
return 0;
- case FIRST_COLUMN: /* It is just started with first column */
+ case FIRST_COLUMN: /* Just started with first column */
Elaboration->west = orig.west - 2 * dim.orlo_v;
Elaboration->east = Elaboration->west + dim.latoE;
General->W = Elaboration->west + 2 * dim.orlo_v;
@@ -78,7 +111,7 @@
Overlap->E = General->E - dim.overlap;
return 0;
- case LAST_COLUMN: /* It is reached last column */
+ case LAST_COLUMN: /* Reached last column */
Elaboration->east = orig.east + 2 * dim.orlo_v;
General->E = Elaboration->east - 2 * dim.orlo_v;
Overlap->E = General->E;
@@ -171,20 +204,22 @@
/*----------------------------------------------------------------------------------------*/
int P_get_orlo(int interpolator, struct Reg_dimens *dim, double pe, double pn)
{
- /* Set the orlo regions dimension->
- * Returns 1 on success of bilinear; 2 on success of bicubic, 0 on failure-> */
- if (interpolator == 1) { /* the interpolator's bilinear */
- dim->orlo_v = 9 * pe; /*4 */
+ /* Set the orlo regions dimension
+ * Returns 1 on success of bilinear; 2 on success of bicubic, 0 on failure */
+ if (interpolator == P_BILINEAR) {
+ /* in case of edge artefacts, increase as multiples of 3 */
+ dim->orlo_v = 9 * pe;
dim->orlo_h = 9 * pn;
return 1;
}
- else if (interpolator == 0) { /* the interpolator's bicubic */
+ else if (interpolator == P_BICUBIC) {
+ /* in case of edge artefacts, increase as multiples of 4 */
dim->orlo_v = 12 * pe; /*3 */
dim->orlo_h = 12 * pn;
return 2;
}
else
- return 0; /* The interpolator it's not bilinear nor bicubic!! */
+ return 0; /* The interpolator is neither bilinear nor bicubic!! */
}
/*----------------------------------------------------------------------------------------*/
@@ -239,7 +274,7 @@
struct bound_box region_box;
struct Cell_head orig;
- G_get_window(&orig);
+ G_get_set_window(&orig);
Vect_region_box(&orig, ®ion_box);
points = Vect_new_line_struct();
@@ -257,7 +292,7 @@
else
z = 0.0;
- /* Reading and storing points only if it is in elaboration_reg */
+ /* only use points in current region */
if (Vect_point_in_box(x, y, z, ®ion_box)) {
npoints++;
@@ -278,7 +313,9 @@
}
}
if (npoints > 0) {
+ /* estimated average distance between points in map units */
*dist = sqrt(((xmax - xmin) * (ymax - ymin)) / npoints);
+ /* estimated point density as number of points per square map unit */
*dens = npoints / ((xmax - xmin) * (ymax - ymin));
return 0;
}
@@ -305,7 +342,7 @@
points = Vect_new_line_struct();
categories = Vect_new_cats_struct();
- /* Reading the elaboration zone points */
+ /* Reading points inside elaboration zone */
Vect_region_box(Elaboration, &elaboration_box);
npoints = 0;
@@ -326,7 +363,7 @@
else
z = 0.0;
- /* Reading and storing points only if it is in elaboration_reg */
+ /* Reading and storing points only if in elaboration_reg */
if (Vect_point_in_box(x, y, z, &elaboration_box)) {
npoints++;
if (npoints >= pippo) {
@@ -355,42 +392,69 @@
}
/*------------------------------------------------------------------------------------------------*/
-int P_Create_Aux_Table(dbDriver * driver, char *tab_name)
+int P_Create_Aux2_Table(dbDriver * driver, char *tab_name)
{
dbTable *auxiliar_tab;
dbColumn *column;
- int created = FALSE;
- auxiliar_tab = db_alloc_table(4);
+ auxiliar_tab = db_alloc_table(2);
db_set_table_name(auxiliar_tab, tab_name);
db_set_table_description(auxiliar_tab,
"Intermediate interpolated values");
- column = db_get_table_column(auxiliar_tab, 2);
- db_set_column_name(column, "Y");
- db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+ column = db_get_table_column(auxiliar_tab, 0);
+ db_set_column_name(column, "ID");
+ db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
column = db_get_table_column(auxiliar_tab, 1);
- db_set_column_name(column, "X");
- db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+ db_set_column_name(column, "Interp");
+ db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
+ if (db_create_table(driver, auxiliar_tab) == DB_OK) {
+ G_debug(1, _("<%s> created in database."), tab_name);
+ return TRUE;
+ }
+ else
+ G_warning(_("<%s> has not been created in database."), tab_name);
+
+ return FALSE;
+}
+
+/*------------------------------------------------------------------------------------------------*/
+int P_Create_Aux4_Table(dbDriver * driver, char *tab_name)
+{
+ dbTable *auxiliar_tab;
+ dbColumn *column;
+
+ auxiliar_tab = db_alloc_table(4);
+ db_set_table_name(auxiliar_tab, tab_name);
+ db_set_table_description(auxiliar_tab,
+ "Intermediate interpolated values");
+
column = db_get_table_column(auxiliar_tab, 0);
db_set_column_name(column, "ID");
db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
- column = db_get_table_column(auxiliar_tab, 3);
+ column = db_get_table_column(auxiliar_tab, 1);
db_set_column_name(column, "Interp");
db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
+ column = db_get_table_column(auxiliar_tab, 2);
+ db_set_column_name(column, "X");
+ db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+
+ column = db_get_table_column(auxiliar_tab, 3);
+ db_set_column_name(column, "Y");
+ db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+
if (db_create_table(driver, auxiliar_tab) == DB_OK) {
G_debug(1, _("<%s> created in database."), tab_name);
- created = TRUE;
- return created;
+ return TRUE;
}
else
- G_fatal_error(_("<%s> has not been created in database."), tab_name);
+ G_warning(_("<%s> has not been created in database."), tab_name);
- return created;
+ return FALSE;
}
/*------------------------------------------------------------------------------------------------*/
@@ -481,7 +545,7 @@
value = db_get_column_value(column);
else
continue;
- coordX = db_get_value_double(value);
+ coordZ = db_get_value_double(value);
column = db_get_table_column(table, 2);
type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
@@ -489,7 +553,7 @@
value = db_get_column_value(column);
else
continue;
- coordY = db_get_value_double(value);
+ coordX = db_get_value_double(value);
column = db_get_table_column(table, 3);
type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
@@ -497,7 +561,7 @@
value = db_get_column_value(column);
else
continue;
- coordZ = db_get_value_double(value);
+ coordY = db_get_value_double(value);
Vect_copy_xyz_to_pnts(point, &coordX, &coordY, &coordZ, 1);
Vect_reset_cats(cat);
@@ -528,10 +592,30 @@
}
#endif
-/*------------------------------------------------------------------------------------------------*/
-#ifdef notdef
-/* Subzones definition */
------------------------|4 | 3 | 3 | ||-----------------------|||||||2 | 5 | 1
- | ||||||||-----------------------|2 | 1 | 1 |
- ||-----------------------||||||||||||||||||-----------------------||||||-----------------------
-#endif
+/*! DEFINITION OF THE SUBZONES
+
+ 5: inside Overlap region
+ all others: inside General region but outside Overlap region
+
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | |4| 3 |3| | |
+ ---------------------------------
+ | | | | | | | |
+ | | |2| 5 |1| | |
+ | | | | | | | |
+ ---------------------------------
+ | | |2| 1 |1| | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ */
Modified: grass/trunk/vector/lidar/v.lidar.correction/correction.c
===================================================================
--- grass/trunk/vector/lidar/v.lidar.correction/correction.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.lidar.correction/correction.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -30,54 +30,65 @@
void
P_Sparse_Correction(struct Map_info *In, struct Map_info *Out,
struct Map_info *Terrain, struct Cell_head *Elaboration,
- struct bound_box General, struct bound_box Overlap, double **obs,
- double *param, int *line_num, double passoN,
- double passoE, double overlap, double HighThresh,
- double LowThresh, int nsplx, int nsply, int num_points,
+ struct bound_box General, struct bound_box Overlap,
+ double **obs, struct lidar_cat *lcat, double *param,
+ int *line_num, double passoN, double passoE,
+ double overlap, double HighThresh, double LowThresh,
+ int nsplx, int nsply, int num_points,
dbDriver * driver, double mean, char *tab_name)
{
int i = 0, class;
double interpolation, csi, eta, weight;
- struct line_pnts *points;
+ struct bound_box elaboration_box;
+ struct line_pnts *point;
struct line_cats *cats;
- points = Vect_new_line_struct();
+ point = Vect_new_line_struct();
cats = Vect_new_cats_struct();
- Vect_rewind(In);
- while (Vect_read_next_line(In, points, cats) > 0) {
- i++;
- if (Vect_point_in_box(*points->x, *points->y, *points->z, &General)) {
+ Vect_region_box(Elaboration, &elaboration_box);
+
+ db_begin_transaction(driver);
+
+ for (i = 0; i < num_points; i++) { /* Sparse points */
+ G_percent(i, num_points, 2);
+ Vect_reset_line(point);
+ Vect_reset_cats(cats);
+
+ if (Vect_point_in_box(obs[i][0], obs[i][1], mean, &General)) {
interpolation =
- dataInterpolateBilin(*points->x, *points->y, passoE, passoN,
+ dataInterpolateBilin(obs[i][0], obs[i][1], passoE, passoN,
nsplx, nsply, Elaboration->west,
Elaboration->south, param);
interpolation += mean;
- if (Vect_point_in_box(*points->x, *points->y, *points->z, &Overlap)) { /*(5) */
+ Vect_copy_xyz_to_pnts(point, &obs[i][0], &obs[i][1], &obs[i][2],
+ 1);
+ *point->z += mean;
- Vect_cat_get(cats, F_CLASSIFICATION, &class);
+ if (Vect_point_in_box(*point->x, *point->y, *point->z, &Overlap)) { /*(5) */
+
+ Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+ Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+ Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+ class = lcat[i].cat_class;
class =
- correction(class, *points->z, interpolation, HighThresh,
+ correction(class, *point->z, interpolation, HighThresh,
LowThresh);
- Vect_cat_del(cats, F_CLASSIFICATION);
Vect_cat_set(cats, F_CLASSIFICATION, class);
if (class == TERRAIN_SINGLE || class == TERRAIN_DOUBLE)
- Vect_write_line(Terrain, GV_POINT, points, cats);
+ Vect_write_line(Terrain, GV_POINT, point, cats);
- Vect_write_line(Out, GV_POINT, points, cats);
-
+ Vect_write_line(Out, GV_POINT, point, cats);
}
else {
-
- if ((*points->x > Overlap.E)) {
-
- if ((*points->y > Overlap.N)) { /*(3) */
- csi = (*points->x - Overlap.E) / overlap;
- eta = (*points->y - Overlap.N) / overlap;
- weight = (1 - csi) * (1 - eta);
+ if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = csi * eta;
interpolation *= weight;
if (Select_Correction
@@ -87,21 +98,19 @@
if (UpDate_Correction
(interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to update the database"));
-
}
- else if ((*points->y < Overlap.S)) { /*(1) */
- csi = (*points->x - Overlap.E) / overlap;
- eta = (*points->y - General.S) / overlap;
- weight = (1 - csi) * eta;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (*point->y - General.S) / overlap;
+ weight = csi * eta;
if (Insert_Correction
(interpolation * weight, line_num[i],
driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to write in the database"));
-
}
- else { /*(1) */
- weight = (*points->x - Overlap.E) / overlap;
+ else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) { /*(1) */
+ weight = (General.E - *point->x) / overlap;
if (Insert_Correction
(interpolation * weight, line_num[i],
@@ -109,32 +118,34 @@
G_fatal_error(_("Impossible to write in the database"));
}
}
- else if ((*points->x < Overlap.W)) {
- if ((*points->y > Overlap.N)) { /*(4) */
- csi = (*points->x - General.W) / overlap;
- eta = (*points->y - Overlap.N) / overlap;
- weight = (1 - eta) * csi;
+ else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(4) */
+ csi = (*point->x - General.W) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = eta * csi;
+ interpolation *= weight;
- interpolation *= weight;
if (Select_Correction
(&interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
- Vect_cat_get(cats, F_CLASSIFICATION, &class);
+ Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+ Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+ Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+ class = lcat[i].cat_class;
class =
- correction(class, *points->z, interpolation,
+ correction(class, *point->z, interpolation,
HighThresh, LowThresh);
Vect_cat_set(cats, F_CLASSIFICATION, class);
- Vect_write_line(Out, GV_POINT, points, cats);
+ Vect_write_line(Out, GV_POINT, point, cats);
if (class == TERRAIN_SINGLE ||
class == TERRAIN_DOUBLE)
- Vect_write_line(Terrain, GV_POINT, points, cats);
-
+ Vect_write_line(Terrain, GV_POINT, point, cats);
}
- else if ((*points->y < Overlap.S)) { /*(2) */
- csi = (*points->x - General.W) / overlap;
- eta = (*points->y - General.S) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(2) */
+ csi = (*point->x - General.W) / overlap;
+ eta = (*point->y - General.S) / overlap;
weight = csi * eta;
interpolation *= weight;
@@ -146,51 +157,55 @@
(interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to update the database"));
}
- else { /*(2) */
- weight = (Overlap.W - *points->x) / overlap;
+ else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) { /*(2) */
+ weight = (*point->x - General.W) / overlap;
interpolation *= weight;
if (Select_Correction
(&interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
- Vect_cat_get(cats, F_CLASSIFICATION, &class);
+ Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+ Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+ Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+ class = lcat[i].cat_class;
class =
- correction(class, *points->z, interpolation,
+ correction(class, *point->z, interpolation,
HighThresh, LowThresh);
Vect_cat_set(cats, F_CLASSIFICATION, class);
- Vect_write_line(Out, GV_POINT, points, cats);
+ Vect_write_line(Out, GV_POINT, point, cats);
if (class == TERRAIN_SINGLE ||
class == TERRAIN_DOUBLE)
- Vect_write_line(Terrain, GV_POINT, points, cats);
-
+ Vect_write_line(Terrain, GV_POINT, point, cats);
}
-
}
- else {
- if ((*points->y > Overlap.N)) { /*(3) */
- weight = (*points->y - Overlap.N) / overlap;
+ else if ((*point->x >= Overlap.W) && (*point->x <= Overlap.E)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ weight = (General.N - *point->y) / overlap;
interpolation *= weight;
if (Select_Correction
(&interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
- Vect_cat_get(cats, F_CLASSIFICATION, &class);
+ Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+ Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+ Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+ class = lcat[i].cat_class;
class =
- correction(class, *points->z, interpolation,
+ correction(class, *point->z, interpolation,
HighThresh, LowThresh);
Vect_cat_set(cats, F_CLASSIFICATION, class);
- Vect_write_line(Out, GV_POINT, points, cats);
+ Vect_write_line(Out, GV_POINT, point, cats);
if (class == TERRAIN_SINGLE ||
class == TERRAIN_DOUBLE)
- Vect_write_line(Terrain, GV_POINT, points, cats);
+ Vect_write_line(Terrain, GV_POINT, point, cats);
}
- else { /*(1) */
- weight = (Overlap.S - *points->y) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ weight = (*point->y - General.S) / overlap;
if (Insert_Correction
(interpolation * weight, line_num[i],
@@ -199,14 +214,14 @@
}
}
}
- }
- /*IF*/ Vect_reset_line(points);
- Vect_reset_cats(cats);
-
- }
- /*FOR*/ Vect_destroy_line_struct(points);
+ } /* if in General box */
+ } /* while */
+ G_percent(num_points, num_points, 2);
+ Vect_destroy_line_struct(point);
Vect_destroy_cats_struct(cats);
+ db_commit_transaction(driver);
+
return;
}
@@ -240,9 +255,8 @@
dbValue *Interp_value;
db_init_string(&sql);
- db_zero_string(&sql);
- sprintf(buf, "SELECT Interp FROM %s WHERE ID=%d", tab_name, line_num);
+ sprintf(buf, "SELECT ID, Interp FROM %s WHERE ID=%d", tab_name, line_num);
db_append_string(&sql, buf);
if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
@@ -260,6 +274,7 @@
*Interp += db_get_value_double(Interp_value);
}
+ db_close_cursor(&cursor);
db_free_string(&sql);
return DB_OK;
}
@@ -300,21 +315,20 @@
struct Point *P_Read_Vector_Correction(struct Map_info *Map,
struct Cell_head *Elaboration,
int *num_points, int *num_terrain,
- int dim_vect)
+ int dim_vect, struct lidar_cat **lcat)
{
- int line_num, npoints, nterrain, pippo, cat_edge;
+ int line_num, npoints, nterrain, pippo, cat;
double x, y, z;
struct Point *obs;
+ struct lidar_cat *lc;
struct line_pnts *points;
struct line_cats *categories;
struct bound_box elaboration_box;
pippo = dim_vect;
- /*if (first_it) */
obs = (struct Point *)G_calloc(pippo, sizeof(struct Point));
- /*else */
- /*obs = (struct Point*) G_realloc ((void *) obs, pippo * sizeof(struct Point)); */
+ lc = (struct lidar_cat *)G_calloc(pippo, sizeof(struct lidar_cat));
points = Vect_new_line_struct();
categories = Vect_new_cats_struct();
@@ -344,12 +358,10 @@
npoints++;
if (npoints >= pippo) {
pippo += dim_vect;
- obs =
- (struct Point *)G_realloc((void *)obs,
- (signed int)(pippo *
- sizeof(struct
- Point)));
- /*lineID = (int *) G_realloc ((void *) lineID, pippo * sizeof(int)); */
+ obs = (struct Point *)G_realloc((void *)obs,
+ (signed int)(pippo * sizeof(struct Point)));
+ lc = (struct lidar_cat *)G_realloc((void *)lc,
+ (signed int)(pippo * sizeof(struct lidar_cat)));
}
/* Storing observation vector */
@@ -357,27 +369,55 @@
obs[npoints - 1].coordY = y;
obs[npoints - 1].coordZ = z;
obs[npoints - 1].lineID = line_num; /* Storing also the line's number */
- Vect_cat_get(categories, F_EDGE_DETECTION_CLASS, &cat_edge);
- obs[npoints - 1].cat = cat_edge;
+ Vect_cat_get(categories, F_EDGE_DETECTION_CLASS, &cat);
+ obs[npoints - 1].cat = cat;
+ lc[npoints - 1].cat_edge = cat;
+
+ /* Only terrain points */
+ if (cat == TERRAIN_SINGLE)
+ nterrain++;
+
+ Vect_cat_get(categories, F_CLASSIFICATION, &cat);
+ lc[npoints - 1].cat_class = cat;
+ Vect_cat_get(categories, F_INTERPOLATION, &cat);
+ lc[npoints - 1].cat_interp = cat;
+ Vect_cat_get(categories, F_COUNTER_OBJ, &cat);
+ lc[npoints - 1].cat_obj = cat;
}
-
- /* Only terrain points */
- if (cat_edge == TERRAIN_SINGLE)
- nterrain++;
}
Vect_destroy_line_struct(points);
Vect_destroy_cats_struct(categories);
*num_points = npoints;
*num_terrain = nterrain;
+ *lcat = lc;
return obs;
}
+/*! DEFINITION OF THE SUBZONES
+ 5: inside Overlap region
+ all others: inside General region but outside Overlap region
-#ifdef notdef
-/* SUBZONES DEFINITION */
------------------------|4 | 3 | 3 | ||-----------------------|||||||2 | 5 | 1
- | ||||||||-----------------------|2 | 1 | 1 |
- ||-----------------------||||||||||||||||||-----------------------||||||-----------------------
-#endif
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | |4| 3 |3| | |
+ ---------------------------------
+ | | | | | | | |
+ | | |2| 5 |1| | |
+ | | | | | | | |
+ ---------------------------------
+ | | |2| 1 |1| | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ */
Modified: grass/trunk/vector/lidar/v.lidar.correction/correction.h
===================================================================
--- grass/trunk/vector/lidar/v.lidar.correction/correction.h 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.lidar.correction/correction.h 2010-02-08 07:51:25 UTC (rev 40860)
@@ -7,6 +7,14 @@
#include <grass/PolimiFunct.h>
+struct lidar_cat
+{
+ int cat_edge; /* category in layer F_EDGE_DETECTION_CLASS */
+ int cat_class; /* category in layer F_CLASSIFICATION */
+ int cat_interp; /* category in layer F_INTERPOLATION */
+ int cat_obj; /* category in layer F_COUNTER_OBJ */
+};
+
void
P_Sparse_Correction(struct Map_info *, /**/
struct Map_info *, /**/
@@ -15,6 +23,7 @@
struct bound_box, /**/
struct bound_box, /**/
double **, /**/
+ struct lidar_cat *,
double *, /**/
int *, /**/
double, /**/
@@ -37,6 +46,7 @@
struct Point *P_Read_Vector_Correction(struct Map_info *, /**/
struct Cell_head *, /**/
- int *, /**/ int *, /**/ int /**/);
+ int *, /**/ int *, /**/ int /**/,
+ struct lidar_cat **);
int correction(int, double, double, double, double);
Modified: grass/trunk/vector/lidar/v.lidar.correction/main.c
===================================================================
--- grass/trunk/vector/lidar/v.lidar.correction/main.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.lidar.correction/main.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -51,7 +51,7 @@
int *lineVect;
double *TN, *Q, *parVect; /* Interpolating and least-square vectors */
- double **N, **obsVect; /* Interpolation and least-square matrix */
+ double **N, **obsVect, **obsVect_all; /* Interpolation and least-square matrix */
struct Map_info In, Out, Terrain;
struct Option *in_opt, *out_opt, *out_terrain_opt, *passoE_opt,
@@ -64,6 +64,7 @@
struct bound_box general_box, overlap_box;
struct Point *observ;
+ struct lidar_cat *lcat;
dbDriver *driver;
@@ -77,9 +78,9 @@
spline_step_flag = G_define_flag();
spline_step_flag->key = 'e';
- spline_step_flag->label = _("Estimate spline step value");
+ spline_step_flag->label = _("Estimate point density and distance");
spline_step_flag->description =
- _("Estimate a good spline step value for the input vector points within the current region extends and quit");
+ _("Estimate point density and distance for the input vector points within the current region extends and quit");
in_opt = G_define_standard_option(G_OPT_V_INPUT);
in_opt->description =
@@ -182,10 +183,14 @@
if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL)
G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
- Vect_set_open_level(1); /*without topology */
+ Vect_set_open_level(1); /* without topology */
if (1 > Vect_open_old(&In, in_opt->answer, mapset))
G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
+ /* Input vector must be 3D */
+ if (!Vect_is_3d(&In))
+ G_fatal_error(_("Input vector map <%s> is not 3D!"), in_opt->answer);
+
/* Estimate point density and mean distance for current region */
if (spline_step_flag->answer) {
double dens, dist;
@@ -228,13 +233,18 @@
/* Create auxiliar table */
if ((flag_auxiliar =
- P_Create_Aux_Table(driver, table_name)) == FALSE) {
+ P_Create_Aux2_Table(driver, table_name)) == FALSE) {
Vect_close(&In);
Vect_close(&Out);
Vect_close(&Terrain);
exit(EXIT_FAILURE);
}
+ db_create_index2(driver, table_name, "ID");
+ /* sqlite likes that */
+ db_close_database_shutdown_driver(driver);
+ driver = db_start_driver_open_database(dvr, db);
+
/* Setting regions and boxes */
G_get_set_window(&original_reg);
G_get_set_window(&elaboration_reg);
@@ -247,27 +257,36 @@
ew_resol = original_reg.ew_res;
ns_resol = original_reg.ns_res;
+ /*------------------------------------------------------------------
+ | Subdividing and working with tiles:
+ | Each original region will be divided into several subregions.
+ | Each one will be overlaped by its neighbouring subregions.
+ | The overlapping is calculated as a fixed OVERLAP_SIZE times
+ | the largest spline step plus 2 * orlo
+ ----------------------------------------------------------------*/
+
/* Fixing parameters of the elaboration region */
- /*! 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);
- N_extension = original_reg.north - original_reg.south;
- E_extension = original_reg.east - original_reg.west;
nsplx_adj = NSPLX_MAX;
nsply_adj = NSPLY_MAX;
- dims.overlap = OVERLAP_SIZE * passoE;
+ if (passoN > passoE)
+ dims.overlap = OVERLAP_SIZE * passoN;
+ else
+ dims.overlap = OVERLAP_SIZE * passoE;
P_get_orlo(P_BILINEAR, &dims, passoE, passoN);
P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
G_verbose_message("adjusted EW splines %d", nsplx_adj);
G_verbose_message("adjusted NS splines %d", nsply_adj);
+ /* calculate number of subzones */
orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+ N_extension = original_reg.north - original_reg.south;
+ E_extension = original_reg.east - original_reg.west;
+
nsubregion_col = ceil(E_extension / orloE) + 0.5;
nsubregion_row = ceil(N_extension / orloN) + 0.5;
@@ -278,9 +297,9 @@
nsubzones = nsubregion_row * nsubregion_col;
- /* Subdividing and working with tiles */
elaboration_reg.south = original_reg.north;
last_row = FALSE;
+
while (last_row == FALSE) { /* For each row */
P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
@@ -344,23 +363,24 @@
G_debug(1, _("read vector region map"));
observ =
P_Read_Vector_Correction(&In, &elaboration_reg, &npoints,
- &nterrain, dim_vect);
+ &nterrain, dim_vect, &lcat);
G_debug(5, _("npoints = %d, nterrain = %d"), npoints, nterrain);
if (npoints > 0) { /* If there is any point falling into elaboration_reg. */
count_terrain = 0;
nparameters = nsplx * nsply;
- /* Mean's calculation */
- G_debug(3, _("Mean's calculation"));
+ /* Mean calculation */
+ G_debug(3, _("Mean calculation"));
mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
/*Least Squares system */
BW = P_get_BandWidth(P_BILINEAR, nsply); /* Bilinear interpolation */
N = G_alloc_matrix(nparameters, BW); /* Normal matrix */
TN = G_alloc_vector(nparameters); /* vector */
- parVect = G_alloc_vector(nparameters); /* Bicubic parameters vector */
- obsVect = G_alloc_matrix(nterrain + 1, 3); /* Observation vector */
+ parVect = G_alloc_vector(nparameters); /* Bilinear parameters vector */
+ obsVect = G_alloc_matrix(nterrain + 1, 3); /* Observation vector with terrain points */
+ obsVect_all = G_alloc_matrix(npoints + 1, 3); /* Observation vector with all points */
Q = G_alloc_vector(nterrain + 1); /* "a priori" var-cov matrix */
lineVect = G_alloc_ivector(npoints + 1);
@@ -375,9 +395,14 @@
count_terrain++;
}
lineVect[i] = observ[i].lineID;
+ obsVect_all[i][0] = observ[i].coordX;
+ obsVect_all[i][1] = observ[i].coordY;
+ obsVect_all[i][2] = observ[i].coordZ - mean;
}
- G_debug(3, _("M.Q. solution"));
+ G_free(observ);
+
+ G_verbose_message(_("Bilinear interpolation"));
normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
nsply, elaboration_reg.west,
elaboration_reg.south, nterrain, nparameters,
@@ -388,19 +413,25 @@
G_free_matrix(N);
G_free_vector(TN);
G_free_vector(Q);
+ G_free_matrix(obsVect);
- G_debug(3, _("Correction and creation of terrain vector"));
+ G_verbose_message( _("Correction and creation of terrain vector"));
P_Sparse_Correction(&In, &Out, &Terrain, &elaboration_reg,
- general_box, overlap_box, obsVect,
+ general_box, overlap_box, obsVect_all, lcat,
parVect, lineVect, passoN, passoE,
dims.overlap, HighThresh, LowThresh,
nsplx, nsply, npoints, driver, mean, table_name);
G_free_vector(parVect);
- G_free_matrix(obsVect);
+ G_free_matrix(obsVect_all);
G_free_ivector(lineVect);
}
- G_free(observ);
+ else {
+ G_free(observ);
+ G_warning(_("No data within this subzone. "
+ "Consider changing the spline step."));
+ }
+ G_free(lcat);
} /*! END WHILE; last_column = TRUE */
} /*! END WHILE; last_row = TRUE */
Modified: grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.c
===================================================================
--- grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -180,6 +180,8 @@
point = Vect_new_line_struct();
categories = Vect_new_cats_struct();
+ db_begin_transaction(driver);
+
for (i = 0; i < num_points; i++) { /* Sparse points */
G_percent(i, num_points, 2);
@@ -201,7 +203,6 @@
/*Vect_cat_set (categories, F_INTERPOLATION, line_out_counter); */
if (Vect_point_in_box(obs[i][0], obs[i][1], interpolation, &Overlap)) { /*(5) */
-
residual = *point->z - interpolation;
edge =
edge_detection(Elaboration, Overlap, parBilin, *point->x,
@@ -214,16 +215,14 @@
Insert_Interpolation(interpolation, line_out_counter, driver,
tabint_name);
line_out_counter++;
-
}
else {
- if ((*point->x > Overlap.E) && (*point->x != General.E)) {
+ if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = csi * eta;
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(3) */
- csi = (*point->x - Overlap.E) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - csi) * (1 - eta);
-
gradient[0] *= weight;
gradient[1] *= weight;
interpolation *= weight;
@@ -237,12 +236,11 @@
interpolation, line_num[i], driver,
tab_name) != DB_OK)
G_fatal_error(_("Impossible to update aux table"));
-
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(1) */
- csi = (*point->x - Overlap.E) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ csi = (General.E - *point->x) / overlap;
eta = (*point->y - General.S) / overlap;
- weight = (1 - csi) * eta;
+ weight = csi * eta;
gradient[0] *= weight;
gradient[1] *= weight;
@@ -254,7 +252,7 @@
}
else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(1) */
- weight = (*point->x - Overlap.E) / overlap;
+ weight = (General.E - *point->x) / overlap;
gradient[0] *= weight;
gradient[1] *= weight;
@@ -264,14 +262,12 @@
line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to write to aux table"));
}
-
}
- else if ((*point->x < Overlap.W) && (*point->x != General.W)) {
-
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(4) */
+ else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(4) */
csi = (*point->x - General.W) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - eta) * csi;
+ eta = (General.N - *point->y) / overlap;
+ weight = eta * csi;
gradient[0] *= weight;
gradient[1] *= weight;
@@ -297,9 +293,8 @@
Insert_Interpolation(interpolation, line_out_counter,
driver, tabint_name);
line_out_counter++;
-
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(2) */
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(2) */
csi = (*point->x - General.W) / overlap;
eta = (*point->y - General.S) / overlap;
weight = csi * eta;
@@ -320,7 +315,7 @@
}
else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(2) */
- weight = (Overlap.W - *point->x) / overlap;
+ weight = (*point->x - General.W) / overlap;
gradient[0] *= weight;
gradient[1] *= weight;
@@ -347,11 +342,10 @@
driver, tabint_name);
line_out_counter++;
}
-
}
else if ((*point->x <= Overlap.E) && (*point->x >= Overlap.W)) {
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(3) */
- weight = (*point->y - Overlap.N) / overlap;
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ weight = (General.N - *point->y) / overlap;
gradient[0] *= weight;
gradient[1] *= weight;
@@ -378,8 +372,8 @@
driver, tabint_name);
line_out_counter++;
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(1) */
- weight = (Overlap.S - *point->y) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ weight = (*point->y - General.S) / overlap;
gradient[0] *= weight;
gradient[1] *= weight;
@@ -394,11 +388,13 @@
}
} /*end if obs */
} /*end for */
- G_percent(i, num_points, 6); /* finish it */
+ G_percent(num_points, num_points, 2); /* finish it */
+ db_commit_transaction(driver);
+
Vect_destroy_line_struct(point);
Vect_destroy_cats_struct(categories);
-} /*end puntisparsi_select */
+}
int Insert(double partialX, double partialY, double Interp,
int line_num, dbDriver * driver, char *tab_name)
@@ -409,7 +405,7 @@
db_init_string(&sql);
sprintf(buf,
- "INSERT INTO %s (ID, Interp, partialX, partialY)", tab_name);
+ "INSERT INTO %s (ID, Interp, X, Y)", tab_name);
db_append_string(&sql, buf);
sprintf(buf, " VALUES (%d, %lf, %lf, %lf)", line_num, Interp,
partialX, partialY);
@@ -450,7 +446,7 @@
db_init_string(&sql);
sprintf(buf,
- "UPDATE %s SET Interp=%lf, PartialX=%lf, PartialY=%lf WHERE ID=%d",
+ "UPDATE %s SET Interp=%lf, X=%lf, Y=%lf WHERE ID=%d",
tab_name, Interp, partialX, partialY, line_num);
db_append_string(&sql, buf);
@@ -473,7 +469,7 @@
db_init_string(&sql);
sprintf(buf,
- "SELECT ID, Interp, partialX, partialY FROM %s WHERE ID=%d",
+ "SELECT ID, Interp, X, Y FROM %s WHERE ID=%d",
tab_name, line_num);
db_append_string(&sql, buf);
@@ -515,89 +511,30 @@
return DB_OK;
}
-int Create_AuxEdge_Table(dbDriver * driver, char *tab_name)
-{
- dbTable *table;
- dbColumn *ID_col, *Interp_col, *PartialX_col, *PartialY_col;
- int created;
-
- table = db_alloc_table(4);
- db_set_table_name(table, tab_name);
- db_set_table_description(table,
- "It is used for the intermediate interpolated and gradient values");
-
- ID_col = db_get_table_column(table, 0);
- db_set_column_name(ID_col, "ID");
- db_set_column_sqltype(ID_col, DB_SQL_TYPE_INTEGER);
-
- Interp_col = db_get_table_column(table, 1);
- db_set_column_name(Interp_col, "Interp");
- db_set_column_sqltype(Interp_col, DB_SQL_TYPE_REAL);
-
- PartialX_col = db_get_table_column(table, 2);
- db_set_column_name(PartialX_col, "PartialX");
- db_set_column_sqltype(PartialX_col, DB_SQL_TYPE_REAL);
-
- PartialY_col = db_get_table_column(table, 3);
- db_set_column_name(PartialY_col, "PartialY");
- db_set_column_sqltype(PartialY_col, DB_SQL_TYPE_REAL);
-
- if (db_create_table(driver, table) == DB_OK) {
- G_debug(1, "<%s> table created in database.", db_get_table_name(table));
- created = TRUE;
- }
- else
- return FALSE;
-
- db_create_index2(driver, tab_name, "ID");
-
- return created;
-}
-
-int Create_Interpolation_Table(dbDriver * driver, char *tab_name)
-{
- dbTable *table;
- dbColumn *ID_col, *Interp_col;
-
- table = db_alloc_table(2);
- db_set_table_name(table, tab_name);
- db_set_table_description(table,
- "This table is the bicubic interpolation of the input vector");
-
- ID_col = db_get_table_column(table, 0);
- db_set_column_name(ID_col, "ID");
- db_set_column_sqltype(ID_col, DB_SQL_TYPE_INTEGER);
-
- Interp_col = db_get_table_column(table, 1);
- db_set_column_name(Interp_col, "Interp");
- db_set_column_sqltype(Interp_col, DB_SQL_TYPE_REAL);
-
- if (db_create_table(driver, table) == DB_OK) {
- G_debug(1, _("<%s> table created in database."), db_get_table_name(table));
- return DB_OK;
- }
- else
- return DB_FAILED;
-}
-
-
-#ifdef notdef
/*! DEFINITION OF THE SUBZONES
- -----------------------
- |4| 3 |3| | |
- -----------------------
- | | | | | |
- |2| 5 |1| | |
- | | | | | |
- -----------------------
- |2| 1 |1| | |
- -----------------------
- | | | | | |
- | | | | | |
- | | | | | |
- -----------------------
- | | | | | |
- -----------------------
+ 5: inside Overlap region
+ all others: inside General region but outside Overlap region
+
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | |4| 3 |3| | |
+ ---------------------------------
+ | | | | | | | |
+ | | |2| 5 |1| | |
+ | | | | | | | |
+ ---------------------------------
+ | | |2| 1 |1| | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
*/
-#endif
Modified: grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.h
===================================================================
--- grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.h 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.lidar.edgedetection/edgedetection.h 2010-02-08 07:51:25 UTC (rev 40860)
@@ -61,6 +61,4 @@
int Select(double *, /**/ double *, /**/ double *, /**/ int, /**/ dbDriver *, /**/ char *);
-int Create_AuxEdge_Table(dbDriver *, char *);
-int Create_Interpolation_Table(dbDriver *, char *);
int Insert_Interpolation(double, int, dbDriver *, char *);
Modified: grass/trunk/vector/lidar/v.lidar.edgedetection/main.c
===================================================================
--- grass/trunk/vector/lidar/v.lidar.edgedetection/main.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.lidar.edgedetection/main.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -34,7 +34,7 @@
#include <grass/PolimiFunct.h>
#include "edgedetection.h"
-int nsply, nsplx, line_out_counter, first_it;
+int nsply, nsplx, line_out_counter;
double passoN, passoE;
/**************************************************************************************
@@ -83,9 +83,9 @@
spline_step_flag = G_define_flag();
spline_step_flag->key = 'e';
- spline_step_flag->label = _("Estimate spline step value");
+ spline_step_flag->label = _("Estimate point density and distance");
spline_step_flag->description =
- _("Estimate a good spline step value for the input vector points within the current region extends and quit");
+ _("Estimate point density and distance for the input vector points within the current region extends and quit");
in_opt = G_define_standard_option(G_OPT_V_INPUT);
@@ -224,6 +224,10 @@
if (1 > Vect_open_old(&In, in_opt->answer, mapset))
G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
+ /* Input vector must be 3D */
+ if (!Vect_is_3d(&In))
+ G_fatal_error(_("Input vector map <%s> is not 3D!"), in_opt->answer);
+
/* Estimate point density and mean distance for current region */
if (spline_step_flag->answer) {
double dens, dist;
@@ -254,41 +258,55 @@
dvr);
/* Create auxiliar and interpolation table */
- if ((flag_auxiliar = Create_AuxEdge_Table(driver, table_name)) == FALSE)
+ if ((flag_auxiliar = P_Create_Aux4_Table(driver, table_name)) == FALSE)
G_fatal_error(_("It was impossible to create <%s>."), table_name);
- if (Create_Interpolation_Table(driver, table_interpolation) != DB_OK)
+ if (P_Create_Aux2_Table(driver, table_interpolation) == FALSE)
G_fatal_error(_("It was impossible to create <%s> interpolation table in database."),
out_opt->answer);
+ db_create_index2(driver, table_name, "ID");
+ db_create_index2(driver, table_interpolation, "ID");
+ /* sqlite likes that */
+ db_close_database_shutdown_driver(driver);
+ driver = db_start_driver_open_database(dvr, db);
+
/* Setting regions and boxes */
G_get_set_window(&original_reg);
G_get_set_window(&elaboration_reg);
Vect_region_box(&elaboration_reg, &overlap_box);
Vect_region_box(&elaboration_reg, &general_box);
+ /*------------------------------------------------------------------
+ | Subdividing and working with tiles:
+ | Each original region will be divided into several subregions.
+ | Each one will be overlaped by its neighbouring subregions.
+ | The overlapping is calculated as a fixed OVERLAP_SIZE times
+ | the largest spline step plus 2 * orlo
+ ----------------------------------------------------------------*/
+
/* Fixing parameters of the elaboration region */
- /*! Each original_region will be divided into several subregions. These
- * subregion will be overlapped by its neighboring subregions. This overlapping
- * is calculated as OVERLAP_SIZE times the east-west spline step. */
-
P_zero_dim(&dims);
- N_extension = original_reg.north - original_reg.south;
- E_extension = original_reg.east - original_reg.west;
-
nsplx_adj = NSPLX_MAX;
nsply_adj = NSPLY_MAX;
- dims.overlap = OVERLAP_SIZE * passoE;
- P_get_orlo(P_BICUBIC, &dims, passoE, passoN); /* Set orlo_h|v */
+ if (passoN > passoE)
+ dims.overlap = OVERLAP_SIZE * passoN;
+ else
+ dims.overlap = OVERLAP_SIZE * passoE;
+ P_get_orlo(P_BICUBIC, &dims, passoE, passoN);
P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
G_verbose_message("adjusted EW splines %d", nsplx_adj);
G_verbose_message("adjusted NS splines %d", nsply_adj);
+ /* calculate number of subzones */
orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+ N_extension = original_reg.north - original_reg.south;
+ E_extension = original_reg.east - original_reg.west;
+
nsubregion_col = ceil(E_extension / orloE) + 0.5;
nsubregion_row = ceil(N_extension / orloN) + 0.5;
@@ -300,9 +318,7 @@
nsubzones = nsubregion_row * nsubregion_col;
elaboration_reg.south = original_reg.north;
-
last_row = FALSE;
- first_it = TRUE;
while (last_row == FALSE) { /* For each row */
@@ -383,7 +399,7 @@
BW = P_get_BandWidth(P_BILINEAR, nsply); /* Bilinear interpolation */
N = G_alloc_matrix(nparameters, BW); /* Normal matrix */
TN = G_alloc_vector(nparameters); /* vector */
- parVect_bilin = G_alloc_vector(nparameters); /* Bicubic parameters vector */
+ parVect_bilin = G_alloc_vector(nparameters); /* Bilinear parameters vector */
obsVect = G_alloc_matrix(npoints + 1, 3); /* Observation vector */
Q = G_alloc_vector(npoints + 1); /* "a priori" var-cov matrix */
@@ -398,7 +414,7 @@
Q[i] = 1; /* Q=I */
}
- /*G_free (observ); */
+ G_free(observ);
G_verbose_message(_("Bilinear interpolation"));
normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
@@ -441,9 +457,11 @@
G_free_matrix(obsVect);
G_free_ivector(lineVect);
} /* IF */
-
- G_free(observ);
- first_it = FALSE;
+ else {
+ G_free(observ);
+ G_warning(_("No data within this subzone. "
+ "Consider changing the spline step."));
+ }
} /*! END WHILE; last_column = TRUE */
} /*! END WHILE; last_row = TRUE */
Modified: grass/trunk/vector/lidar/v.lidar.growing/main.c
===================================================================
--- grass/trunk/vector/lidar/v.lidar.growing/main.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.lidar.growing/main.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -84,6 +84,8 @@
"algorithm for determining the building inside");
in_opt = G_define_standard_option(G_OPT_V_INPUT);
+ in_opt->description =
+ _("Input vector (v.lidar.edgedetection output");
out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
@@ -166,30 +168,8 @@
G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
field->driver);
- /* Setting regions and boxes */
- G_debug(1, _("Setting regions and boxes"));
- G_get_set_window(&original_reg);
- G_get_set_window(&elaboration_reg);
-
- /* Fixing parameters of the elaboration region */
- /*! The original_region will be divided into several subregions */
- ew_resol = original_reg.ew_res;
- ns_resol = original_reg.ns_res;
-
- nsubregion_col = ceil((original_reg.east - original_reg.west) / (LATO * ew_resol)) + 0.5;
- nsubregion_row = ceil((original_reg.north - original_reg.south) / (LATO * ns_resol)) + 0.5;
-
- if (nsubregion_col < 0)
- nsubregion_col = 0;
- if (nsubregion_row < 0)
- nsubregion_row = 0;
-
- nsubzones = nsubregion_row * nsubregion_col;
-
- /* Subdividing and working with tiles */
- elaboration_reg.south = original_reg.north;
- last_row = FALSE;
-
+ /* is this the right place to open the cursor ??? */
+
db_init_string(&sql);
db_zero_string(&sql);
@@ -202,6 +182,7 @@
count_obj = 1;
+ /* no topology, get number of lines in input vector */
nlines = 0;
points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
@@ -211,6 +192,7 @@
}
Vect_rewind(&In);
+ /* no topology, get number of lines in first pulse input vector */
nlines_first = 0;
points_first = Vect_new_line_struct();
Cats_first = Vect_new_cats_struct();
@@ -220,6 +202,31 @@
}
Vect_rewind(&First);
+ /* Setting regions and boxes */
+ G_debug(1, _("Setting regions and boxes"));
+ G_get_set_window(&original_reg);
+ G_get_set_window(&elaboration_reg);
+
+ /* Fixing parameters of the elaboration region */
+ /*! The original_region will be divided into subregions */
+ ew_resol = original_reg.ew_res;
+ ns_resol = original_reg.ns_res;
+
+ /* calculate number of subzones */
+ nsubregion_col = ceil((original_reg.east - original_reg.west) / (LATO * ew_resol)) + 0.5;
+ nsubregion_row = ceil((original_reg.north - original_reg.south) / (LATO * ns_resol)) + 0.5;
+
+ if (nsubregion_col < 0)
+ nsubregion_col = 0;
+ if (nsubregion_row < 0)
+ nsubregion_row = 0;
+
+ nsubzones = nsubregion_row * nsubregion_col;
+
+ /* Subdividing and working with tiles */
+ elaboration_reg.south = original_reg.north;
+ last_row = FALSE;
+
while (last_row == FALSE) { /* For each strip of LATO rows */
elaboration_reg.north = elaboration_reg.south;
@@ -259,6 +266,8 @@
elaboration_reg.ew_res = ew_resol;
nrows = (elaboration_reg.north - elaboration_reg.south) / ns_resol + 0.1;
ncols = (elaboration_reg.east - elaboration_reg.west) / ew_resol + 0.1;
+ elaboration_reg.rows = nrows;
+ elaboration_reg.cols = ncols;
G_debug(1, _("Rows = %d"), nrows);
G_debug(1, _("Columns = %d"), ncols);
@@ -281,13 +290,13 @@
}
}
+ G_verbose_message(_("read points in input vector"));
+ Vect_region_box(&elaboration_reg, &elaboration_box);
line_num = 0;
Vect_rewind(&In);
while (Vect_read_next_line(&In, points, Cats) > 0) {
line_num++;
- Vect_region_box(&elaboration_reg, &elaboration_box);
-
if ((Vect_point_in_box
(points->x[0], points->y[0], points->z[0],
&elaboration_box)) &&
@@ -304,6 +313,36 @@
(points->x[0], &elaboration_reg));
Z_interp = 0;
+ /* TODO: make sure the current db_fetch() usage works */
+ /* why not: */
+ /*
+ db_init_string(&sql);
+ sprintf(buf, "SELECT Interp,ID FROM %s WHERE ID=%d", table_name, line_num);
+ db_append_string(&sql, buf);
+
+ if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
+ G_fatal_error(_("Unable to open table <%s>"), table_name);
+
+ while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
+ dbColumn *Z_Interp_col;
+ dbValue *Z_Interp_value;
+ table = db_get_cursor_table(&cursor);
+
+ Z_Interp_col = db_get_table_column(table, 1);
+
+ if (db_sqltype_to_Ctype(db_get_column_sqltype(Z_Interp_col)) ==
+ DB_C_TYPE_DOUBLE)
+ Z_Interp_value = db_get_column_value(Z_Interp_col);
+ else
+ continue;
+
+ Z_interp = db_get_value_double(Z_Interp_value);
+ break;
+ }
+ db_close_cursor(&cursor);
+ db_free_string(&sql);
+ */
+ /* instead of */
while (1) {
if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK ||
!more)
@@ -400,7 +439,7 @@
/* REGION GROWING */
if (region == TRUE) {
- G_debug(1, "Region Growing");
+ G_verbose_message(_("Region Growing"));
punti_bordo = G_alloc_matrix(MaxPoints, 3);
P = Pvector(0, MaxPoints);
@@ -409,6 +448,7 @@
ripieno = 6;
for (row = 0; row <= nrows; row++) {
+ G_percent(row, nrows, 2);
for (col = 0; col <= ncols; col++) {
if ((raster_matrix[row][col].clas >= Thres_j) &&
@@ -437,7 +477,7 @@
punti_bordo, &lungPunti, row, col,
colorBordo, Thres_j, MaxPoints);
- /*CONVEX-HULL COMPUTATION */
+ /* CONVEX-HULL COMPUTATION */
lungHull = ch2d(P, lungPunti);
cvxHull = G_alloc_matrix(lungHull, 3);
Modified: grass/trunk/vector/lidar/v.outlier/main.c
===================================================================
--- grass/trunk/vector/lidar/v.outlier/main.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.outlier/main.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -29,7 +29,7 @@
#include <grass/PolimiFunct.h>
#include "outlier.h"
/* GLOBAL VARIABLES DEFINITIONS */
-int nsply, nsplx, first_it;
+int nsply, nsplx;
double passoN, passoE, Thres_Outlier;
/*--------------------------------------------------------------------------------------*/
@@ -41,7 +41,7 @@
int subzone = 0, nsubzones = 0;
double N_extension, E_extension, orloE, orloN;
int dim_vect, nparameters, BW, npoints;
- double ew_resol, ns_resol, mean, lambda;
+ double mean, lambda;
const char *dvr, *db, *mapset;
char table_name[GNAME_MAX];
char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
@@ -76,9 +76,9 @@
spline_step_flag = G_define_flag();
spline_step_flag->key = 'e';
- spline_step_flag->label = _("Estimate spline step value");
+ spline_step_flag->label = _("Estimate point density and distance");
spline_step_flag->description =
- _("Estimate a good spline step value for the input vector points within the current region extends and quit");
+ _("Estimate point density and distance for the input vector points within the current region extends and quit");
in_opt = G_define_standard_option(G_OPT_V_INPUT);
@@ -147,6 +147,8 @@
lambda = atof(lambda_f_opt->answer);
Thres_Outlier = atof(Thres_O_opt->answer);
+ flag_auxiliar = FALSE;
+
/* Checking vector names */
Vect_check_input_output_name(in_opt->answer, out_opt->answer,
GV_FATAL_EXIT);
@@ -180,6 +182,10 @@
G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
in_opt->answer);
+ /* Input vector must be 3D */
+ if (!Vect_is_3d(&In))
+ G_fatal_error(_("Input vector map <%s> is not 3D!"), in_opt->answer);
+
/* Estimate point density and mean distance for current region */
if (spline_step_flag->answer) {
double dens, dist;
@@ -232,37 +238,52 @@
G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
dvr);
+ /* Create auxiliar table */
+ if ((flag_auxiliar =
+ P_Create_Aux2_Table(driver, table_name)) == FALSE)
+ G_fatal_error(_("It was impossible to create <%s> table."), table_name);
+
+ db_create_index2(driver, table_name, "ID");
+ /* sqlite likes that */
+ db_close_database_shutdown_driver(driver);
+ driver = db_start_driver_open_database(dvr, db);
+
/* Setting regions and boxes */
G_get_set_window(&original_reg);
G_get_set_window(&elaboration_reg);
Vect_region_box(&elaboration_reg, &overlap_box);
Vect_region_box(&elaboration_reg, &general_box);
+ /*------------------------------------------------------------------
+ | Subdividing and working with tiles:
+ | Each original region will be divided into several subregions.
+ | Each one will be overlaped by its neighbouring subregions.
+ | The overlapping is calculated as a fixed OVERLAP_SIZE times
+ | the largest spline step plus 2 * orlo
+ ----------------------------------------------------------------*/
+
/* Fixing parameters of the elaboration region */
- /*! Each original_region will be divided into several subregions. These
- * subregion will be overlapped by its neibourgh subregions. This overlapping
- * is calculated as OVERLAP_SIZE times the east-west spline step. */
-
- ew_resol = original_reg.ew_res;
- ns_resol = original_reg.ns_res;
-
P_zero_dim(&dims);
- N_extension = original_reg.north - original_reg.south;
- E_extension = original_reg.east - original_reg.west;
-
nsplx_adj = NSPLX_MAX;
nsply_adj = NSPLY_MAX;
- dims.overlap = OVERLAP_SIZE * passoE;
+ if (passoN > passoE)
+ dims.overlap = OVERLAP_SIZE * passoN;
+ else
+ dims.overlap = OVERLAP_SIZE * passoE;
P_get_orlo(P_BILINEAR, &dims, passoE, passoN);
P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
G_verbose_message("adjusted EW splines %d", nsplx_adj);
G_verbose_message("adjusted NS splines %d", nsply_adj);
+ /* calculate number of subzones */
orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+ N_extension = original_reg.north - original_reg.south;
+ E_extension = original_reg.east - original_reg.west;
+
nsubregion_col = ceil(E_extension / orloE) + 0.5;
nsubregion_row = ceil(N_extension / orloN) + 0.5;
@@ -274,9 +295,7 @@
nsubzones = nsubregion_row * nsubregion_col;
elaboration_reg.south = original_reg.north;
-
last_row = FALSE;
- first_it = TRUE;
while (last_row == FALSE) { /* For each row */
@@ -370,6 +389,8 @@
Q[i] = 1; /* Q=I */
}
+ G_free(observ);
+
G_verbose_message(_("Bilinear interpolation"));
normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
nsply, elaboration_reg.west,
@@ -382,14 +403,6 @@
G_free_vector(TN);
G_free_vector(Q);
- if (flag_auxiliar == FALSE) {
- G_debug(1,
- "Creating auxiliar table for archiving overlapping zones");
- if ((flag_auxiliar =
- P_Create_Aux_Table(driver, table_name)) == FALSE)
- G_fatal_error(_("It was impossible to create <Auxiliar_outlier_table>."));
- }
-
G_verbose_message(_("Outlier detection"));
if (qgis_opt->answer)
P_Outlier(&Out, &Outlier, &Qgis, elaboration_reg,
@@ -406,8 +419,11 @@
G_free_ivector(lineVect);
} /*! END IF; npoints > 0 */
- G_free(observ);
- first_it = FALSE;
+ else {
+ G_free(observ);
+ G_warning(_("No data within this subzone. "
+ "Consider changing the spline step."));
+ }
} /*! END WHILE; last_column = TRUE */
} /*! END WHILE; last_row = TRUE */
Modified: grass/trunk/vector/lidar/v.outlier/outlier.c
===================================================================
--- grass/trunk/vector/lidar/v.outlier/outlier.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.outlier/outlier.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -27,6 +27,8 @@
point = Vect_new_line_struct();
categories = Vect_new_cats_struct();
+ db_begin_transaction(driver);
+
for (i = 0; i < num_points; i++) { /* Sparse points */
G_percent(i, num_points, 2);
Vect_reset_line(point);
@@ -60,13 +62,12 @@
}
else {
- if ((*point->x > Overlap.E) && (*point->x != General.E)) {
+ if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = csi * eta;
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(3) */
- csi = (*point->x - Overlap.E) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - csi) * (1 - eta);
-
interpolation *= weight;
if (Select_Outlier(&interpolation, line_num[i],
@@ -76,22 +77,20 @@
if (UpDate_Outlier(interpolation, line_num[i],
driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to update the database"));
-
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(1) */
- csi = (*point->x - Overlap.E) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ csi = (General.E - *point->x) / overlap;
eta = (*point->y - General.S) / overlap;
- weight = (1 - csi) * eta;
+ weight = csi * eta;
interpolation *= weight;
if (Insert_Outlier(interpolation, line_num[i],
driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to write in the database"));
-
}
else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(1) */
- weight = (*point->x - Overlap.E) / overlap;
+ weight = (General.E - *point->x) / overlap;
interpolation *= weight;
@@ -99,14 +98,12 @@
driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to write in the database"));
}
-
}
- else if ((*point->x < Overlap.W) && (*point->x != General.W)) {
-
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(4) */
+ else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(4) */
csi = (*point->x - General.W) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - eta) * csi;
+ eta = (General.N - *point->y) / overlap;
+ weight = eta * csi;
interpolation *= weight;
@@ -123,12 +120,13 @@
Vect_write_line(Qgis, GV_POINT, point,
categories);
}
- else
+ else {
Vect_write_line(Outlier, GV_POINT, point,
categories);
-
+ G_message("here we are");
+ }
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(2) */
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(2) */
csi = (*point->x - General.W) / overlap;
eta = (*point->y - General.S) / overlap;
weight = csi * eta;
@@ -142,10 +140,9 @@
if (UpDate_Outlier(interpolation, line_num[i],
driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to update the database"));
-
}
else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(2) */
- weight = (Overlap.W - *point->x) / overlap;
+ weight = (*point->x - General.W) / overlap;
interpolation *= weight;
@@ -166,11 +163,10 @@
Vect_write_line(Outlier, GV_POINT, point,
categories);
}
-
}
else if ((*point->x <= Overlap.E) && (*point->x >= Overlap.W)) {
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(3) */
- weight = (*point->y - Overlap.N) / overlap;
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ weight = (General.N - *point->y) / overlap;
interpolation *= weight;
@@ -190,10 +186,9 @@
else
Vect_write_line(Outlier, GV_POINT, point,
categories);
-
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(1) */
- weight = (Overlap.S - *point->y) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ weight = (*point->y - General.S) / overlap;
interpolation *= weight;
@@ -210,6 +205,8 @@
G_percent(num_points, num_points, 2);
G_debug(2, "P_outlier: done");
+ db_commit_transaction(driver);
+
Vect_destroy_line_struct(point);
Vect_destroy_cats_struct(categories);
} /*end puntisparsi_select */
@@ -295,23 +292,30 @@
return TRUE;
}
-#ifdef notdef
/*! DEFINITION OF THE SUBZONES
- -----------------------
- |4| 3 |3| | |
- -----------------------
- | | | | | |
- |2| 5 |1| | |
- | | | | | |
- -----------------------
- |2| 1 |1| | |
- -----------------------
- | | | | | |
- | | | | | |
- | | | | | |
- -----------------------
- | | | | | |
- -----------------------
+ 5: inside Overlap region
+ all others: inside General region but outside Overlap region
+
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | |4| 3 |3| | |
+ ---------------------------------
+ | | | | | | | |
+ | | |2| 5 |1| | |
+ | | | | | | | |
+ ---------------------------------
+ | | |2| 1 |1| | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
*/
-#endif
Modified: grass/trunk/vector/lidar/v.surf.bspline/crosscorr.c
===================================================================
--- grass/trunk/vector/lidar/v.surf.bspline/crosscorr.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.surf.bspline/crosscorr.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -49,7 +49,8 @@
int nsplx, nsply, nparam_spl, ndata;
double *mean, *rms, *stdev, rms_min, stdev_min;
- double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; /* Fixed values (by the moment) */
+ /* double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; */ /* Fixed values (by the moment) */
+ double lambda[PARAM_LAMBDA] = { 0.01, 0.05, 0.1, 0.2, 0.3, 0.4 }; /* Fixed values (by the moment) */
double *TN, *Q, *parVect; /* Interpolation and least-square vectors */
double **N, **obsVect; /* Interpolation and least-square matrix */
Modified: grass/trunk/vector/lidar/v.surf.bspline/main.c
===================================================================
--- grass/trunk/vector/lidar/v.surf.bspline/main.c 2010-02-08 06:18:45 UTC (rev 40859)
+++ grass/trunk/vector/lidar/v.surf.bspline/main.c 2010-02-08 07:51:25 UTC (rev 40860)
@@ -38,7 +38,7 @@
int main(int argc, char *argv[])
{
/* Variable declarations */
- int nsply, nsplx, nlines, nrows, ncols, nsplx_adj, nsply_adj;
+ int nsply, nsplx, nrows, ncols, nsplx_adj, nsply_adj;
int nsubregion_col, nsubregion_row, subregion_row, subregion_col;
int subzone = 0, nsubzones = 0;
int last_row, last_column, grid, bilin, ext, flag_auxiliar, cross; /* booleans */
@@ -225,8 +225,12 @@
G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
in_opt->answer);
- if (!Vect_is_3d(&In) && bspline_field <= 0 && bspline_column)
+ /* check availability of z values
+ * column option overrrides 3D z coordinates */
+ if (!Vect_is_3d(&In) && bspline_field <= 0 && bspline_column == NULL)
G_fatal_error(_("Need either 3D vector or layer and column with z values"));
+ if (bspline_field > 0 && bspline_column == NULL)
+ G_fatal_error(_("Layer but not column with z values given"));
/* Estimate point density and mean distance for current region */
if (spline_step_flag->answer) {
@@ -316,6 +320,7 @@
raster = Rast_open_fp_new(out_map_opt->answer);
}
+ /* read z values from attribute table */
if (bspline_field > 0) {
db_CatValArray_init(&cvarr);
Fi = Vect_get_field(&In, bspline_field);
@@ -356,14 +361,18 @@
G_fatal_error(_("No database connection for driver <%s> is defined. "
"Run db.connect."), dvr);
- /* Auxiliar table creation */
+ /* Create auxiliar table */
if (vector) {
- if ((flag_auxiliar = P_Create_Aux_Table(driver, table_name)) == FALSE) {
+ if ((flag_auxiliar = P_Create_Aux4_Table(driver, table_name)) == FALSE) {
P_Drop_Aux_Table(driver, table_name);
G_fatal_error(_("Interpolation: Creating table: "
"It was impossible to create table <%s>."),
table_name);
}
+ db_create_index2(driver, table_name, "ID");
+ /* sqlite likes that */
+ db_close_database_shutdown_driver(driver);
+ driver = db_start_driver_open_database(dvr, db);
}
/* Setting regions and boxes */
@@ -383,24 +392,36 @@
"Consider changing region resolution"));
}
+ /*------------------------------------------------------------------
+ | Subdividing and working with tiles:
+ | Each original region will be divided into several subregions.
+ | Each one will be overlaped by its neighbouring subregions.
+ | The overlapping is calculated as a fixed OVERLAP_SIZE times
+ | the largest spline step plus 2 * orlo
+ ----------------------------------------------------------------*/
+
/* Fixing parameters of the elaboration region */
- P_zero_dim(&dims); /* Set to zero the dim struct */
+ P_zero_dim(&dims); /* Set dim struct to zero */
- 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 orlo_h|v */
+ if (passoN > passoE)
+ dims.overlap = OVERLAP_SIZE * passoN;
+ else
+ dims.overlap = OVERLAP_SIZE * passoE;
+ P_get_orlo(bilin, &dims, passoE, passoN);
P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
G_verbose_message("adjusted EW splines %d", nsplx_adj);
G_verbose_message("adjusted NS splines %d", nsply_adj);
+ /* calculate number of subzones */
orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+ N_extension = original_reg.north - original_reg.south;
+ E_extension = original_reg.east - original_reg.west;
+
nsubregion_col = ceil(E_extension / orloE) + 0.5;
nsubregion_row = ceil(N_extension / orloN) + 0.5;
@@ -415,21 +436,11 @@
points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
Vect_cat_set(Cats, 1, 0);
- nlines = Vect_get_num_lines(&In);
- /*------------------------------------------------------------------
- | Subdividing and working with tiles:
- | Each original_region will be divided into several subregions.
- | Each one will be overlaped by its neibourgh subregions.
- | The overlapping was calculated as a fixed OVERLAP_SIZE times
- | the east-west spline step
- ----------------------------------------------------------------*/
-
+ subregion_row = 0;
elaboration_reg.south = original_reg.north;
+ last_row = FALSE;
- G_percent(0, 1, 10);
- subregion_row = 0;
- last_row = FALSE;
while (last_row == FALSE) { /* For each subregion row */
subregion_row++;
P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
@@ -511,7 +522,6 @@
if (npoints > 0) { /* */
int i;
- double *obs_mean;
nparameters = nsplx * nsply;
BW = P_get_BandWidth(P_BILINEAR, nsply) + 2 * (bilin == P_BICUBIC);
@@ -523,7 +533,6 @@
obsVect = G_alloc_matrix(npoints, 3); /* Observation vector */
Q = G_alloc_vector(npoints); /* "a priori" var-cov matrix */
lineVect = G_alloc_ivector(npoints); /* */
- obs_mean = G_alloc_vector(npoints);
for (i = 0; i < npoints; i++) { /* Setting obsVect vector & Q matrix */
double dval;
@@ -533,6 +542,7 @@
obsVect[i][0] = observ[i].coordX;
obsVect[i][1] = observ[i].coordY;
+ /* read z coordinates from attribute table */
if (bspline_field > 0) {
int cat, ival, ret;
@@ -545,14 +555,14 @@
db_CatValArray_get_value_int(&cvarr, cat,
&ival);
obsVect[i][2] = ival;
- obs_mean[i] = ival;
+ observ[i].coordZ = ival;
}
else { /* DB_C_TYPE_DOUBLE */
ret =
db_CatValArray_get_value_double(&cvarr, cat,
&dval);
obsVect[i][2] = dval;
- obs_mean[i] = dval;
+ observ[i].coordZ = dval;
}
if (ret != DB_OK) {
G_warning(_("Interpolation: (%d,%d): No record for point (cat = %d)"),
@@ -560,29 +570,25 @@
continue;
}
}
-
+ /* use z coordinates of 3D vector */
else {
obsVect[i][2] = observ[i].coordZ;
- obs_mean[i] = observ[i].coordZ;
- } /* obsVect[i][2] = observ[i].coordZ - mean; */
+ }
}
/* Mean calculation for every point */
- if (bspline_field > 0)
- mean = calc_mean(obs_mean, npoints);
- else
- mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
+ mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
-
G_debug(1, "Interpolation: (%d,%d): mean=%lf",
subregion_row, subregion_col, mean);
+ G_free(observ);
+
for (i = 0; i < npoints; i++)
obsVect[i][2] -= mean;
- G_free(observ);
-
- if (bilin) { /* Bilinear interpolation */
+ /* Bilinear interpolation */
+ if (bilin) {
G_debug(1,
"Interpolation: (%d,%d): Bilinear interpolation...",
subregion_row, subregion_col);
@@ -592,6 +598,7 @@
nparameters, BW);
nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
}
+ /* Bicubic interpolation */
else {
G_debug(1,
"Interpolation: (%d,%d): Bicubic interpolation...",
@@ -610,7 +617,6 @@
G_free_vector(Q);
if (grid == TRUE) { /* GRID INTERPOLATION ==> INTERPOLATION INTO A RASTER */
- flag_auxiliar = TRUE;
G_debug(1, "Interpolation: (%d,%d): Regular_Points...",
subregion_row, subregion_col);
raster_matrix =
@@ -669,11 +675,12 @@
G_free_vector(parVect);
G_free_matrix(obsVect);
G_free_ivector(lineVect);
- G_free_vector(obs_mean);
}
- else
+ else {
+ G_free(observ);
G_warning(_("No data within this subzone. "
"Consider changing the spline step."));
+ }
} /*! END WHILE; last_column = TRUE */
} /*! END WHILE; last_row = TRUE */
More information about the grass-commit
mailing list