[GRASS-SVN] r41062 - in grass/branches/develbranch_6/vector/lidar:
lidarlib v.lidar.correction v.lidar.edgedetection
v.lidar.growing v.outlier v.surf.bspline
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Feb 17 12:56:53 EST 2010
Author: mmetz
Date: 2010-02-17 12:56:52 -0500 (Wed, 17 Feb 2010)
New Revision: 41062
Modified:
grass/branches/develbranch_6/vector/lidar/lidarlib/PolimiFunct.h
grass/branches/develbranch_6/vector/lidar/lidarlib/raster.c
grass/branches/develbranch_6/vector/lidar/lidarlib/zones.c
grass/branches/develbranch_6/vector/lidar/v.lidar.correction/correction.c
grass/branches/develbranch_6/vector/lidar/v.lidar.correction/correction.h
grass/branches/develbranch_6/vector/lidar/v.lidar.correction/main.c
grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/edgedetection.c
grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/edgedetection.h
grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/main.c
grass/branches/develbranch_6/vector/lidar/v.lidar.growing/main.c
grass/branches/develbranch_6/vector/lidar/v.outlier/main.c
grass/branches/develbranch_6/vector/lidar/v.outlier/outlier.c
grass/branches/develbranch_6/vector/lidar/v.outlier/outlier.h
grass/branches/develbranch_6/vector/lidar/v.surf.bspline/crosscorr.c
grass/branches/develbranch_6/vector/lidar/v.surf.bspline/main.c
Log:
fix lidar tools, backport from trunk r41061
Modified: grass/branches/develbranch_6/vector/lidar/lidarlib/PolimiFunct.h
===================================================================
--- grass/branches/develbranch_6/vector/lidar/lidarlib/PolimiFunct.h 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/lidarlib/PolimiFunct.h 2010-02-17 17:56:52 UTC (rev 41062)
@@ -28,36 +28,36 @@
#include <grass/gmath.h>
/*----------------------------------------------------------------------------------------------------------*/
-/*CONSTANS DECLARATION */
+/*CONSTANTS DECLARATION */
-#define NSPLX_MAX 200 /* Maximum number of splines along East direction used in the subregions interpolation */
-#define NSPLY_MAX 200 /* Maximum number of splines along North direction used in the subregions interpolation */
-#define OVERLAP_SIZE 20 /* Subregions overlapping size. */
-#define LATO 2000 /* Side's size for v.to.qrast command. */
-#define CONTOUR 15 /**/
+#define NSPLX_MAX 150 /* Maximum number of splines along East direction used in the subregions interpolation */
+#define NSPLY_MAX 150 /* Maximum number of splines along North direction used in the subregions interpolation */
+#define OVERLAP_SIZE 10 /* Subregions overlapping size. */
+#define LATO 1000 /* Side's size for v.lidar.growing. */
+#define CONTOUR 15 /**/
#define GENERAL_ROW 0
-#define GENERAL_COLUMN 1
+#define GENERAL_COLUMN 1
#define FIRST_ROW 2
#define LAST_ROW 3
#define FIRST_COLUMN 4
#define LAST_COLUMN 5
/* FIELDS ID */
#define F_EDGE_DETECTION_CLASS 1
-#define F_CLASSIFICATION 2
+#define F_CLASSIFICATION 2
#define F_INTERPOLATION 3
-#define F_COUNTER_OBJ 4
+#define F_COUNTER_OBJ 4
/* PRE-CLASSIFICATION */
#define PRE_TERRAIN 1
#define PRE_EDGE 2
#define PRE_UNKNOWN 3
/* FINAL CLASSIFICATION */
-#define TERRAIN_SINGLE 1
-#define TERRAIN_DOUBLE 2
+#define TERRAIN_SINGLE 1
+#define TERRAIN_DOUBLE 2
#define OBJECT_DOUBLE 3
#define OBJECT_SINGLE 4
/* SINGLE OR DOUBLE PULSE */
-#define SINGLE_PULSE 1
-#define DOUBLE_PULSE 2
+#define SINGLE_PULSE 1
+#define DOUBLE_PULSE 2
/* INTERPOLATOR */
#define P_BILINEAR 1
#define P_BICUBIC 0
@@ -96,6 +96,7 @@
/*FUNCTIONS DECLARATION */
/*zones */
void P_zero_dim(struct Reg_dimens * /**/);
+int P_set_dim(struct Reg_dimens *, double, double, int *, int *);
int P_set_regions(struct Cell_head *, /**/
BOUND_BOX *, /**/
@@ -105,6 +106,8 @@
int P_get_BandWidth(int, /**/ int /**/);
+double P_estimate_splinestep(struct Map_info *, double *, double *);
+
struct Point *P_Read_Vector_Region_Map(struct Map_info *, /**/
struct Cell_head *, /**/
int *, /**/ int, /**/ int /**/);
@@ -142,8 +145,10 @@
int, /**/ int, /**/ int, /**/ int, /**/ int /**/);
/*----------------------------------------------------------------------------------------------------------*/
-int P_Create_Aux_Table(dbDriver *, /**/ char * /**/);
+int P_Create_Aux2_Table(dbDriver *, /**/ char * /**/);
+int P_Create_Aux4_Table(dbDriver *, /**/ char * /**/);
+
int P_Drop_Aux_Table(dbDriver *, /**/ char * /**/);
/*----------------------------------------------------------------------------------------------------------*/
Modified: grass/branches/develbranch_6/vector/lidar/lidarlib/raster.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/lidarlib/raster.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/lidarlib/raster.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -21,16 +21,15 @@
{
int i;
char buf[1024];
+ dbString sql;
double interpolation, csi, eta, weight;
struct line_pnts *point;
- dbString sql;
-
- /*CREARE LA TABELLA */
-
point = Vect_new_line_struct();
+ db_begin_transaction(driver);
+
for (i = 0; i < num_points; i++) {
if (Vect_point_in_box(obs[i][0], obs[i][1], mean, &General)) { /*Here mean is just for asking if obs point is in box */
@@ -52,10 +51,8 @@
if (Vect_point_in_box(obs[i][0], obs[i][1], interpolation, &Overlap)) { /*(5) */
Vect_write_line(Out, GV_POINT, point, categories);
-
}
else {
-
db_init_string(&sql);
sprintf(buf, "INSERT INTO %s (ID, X, Y, Interp)", tab_name);
@@ -67,12 +64,11 @@
obs[i][1]);
db_append_string(&sql, buf);
- if ((point->x[0] > Overlap.E)) {
-
- if ((*point->y > Overlap.N)) { /*(3) */
- csi = (*point->x - Overlap.E) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - csi) * (1 - eta);
+ if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = csi * eta;
*point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
@@ -81,14 +77,13 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to access table <%s>"),
buf);
-
}
- else if ((*point->y < Overlap.S)) { /*(1) */
- csi = (*point->x - Overlap.E) / overlap;
- eta = (Overlap.S - *point->y) / overlap;
- weight = (1 - csi) * eta;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (*point->y - General.S) / overlap;
+ weight = csi * eta;
*point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
@@ -97,13 +92,12 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to access table <%s>"),
buf);
-
}
- else { /*(1) */
- weight = (*point->x - Overlap.E) / overlap;
- *point->z = (1 - weight) * interpolation;
+ else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(1) */
+ weight = (General.E - *point->x) / overlap;
+ *point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
db_append_string(&sql, buf);
@@ -111,16 +105,15 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to access table <%s>"),
buf);
}
-
}
- else if ((point->x[0] < Overlap.W)) {
- if ((*point->y > Overlap.N)) { /*(4) */
- csi = (Overlap.W - *point->x) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - eta) * csi;
+ else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(4) */
+ csi = (*point->x - General.W) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = eta * csi;
*point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
@@ -129,13 +122,12 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to access table <%s>"),
buf);
-
}
- else if ((*point->y < Overlap.S)) { /*(2) */
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(2) */
csi = (*point->x - General.W) / overlap;
- eta = (Overlap.S - *point->y) / overlap;
+ eta = (*point->y - General.S) / overlap;
weight = csi * eta;
*point->z = weight * interpolation;
@@ -145,13 +137,12 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to access table <%s>"),
buf);
-
}
- else { /*(2) */
- weight = (Overlap.W - *point->x) / overlap;
- *point->z = (1 - weight) * interpolation;
+ else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) { /*(2) */
+ weight = (*point->x - General.W) / overlap;
+ *point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
db_append_string(&sql, buf);
@@ -159,15 +150,14 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to access table <%s>"),
buf);
}
-
}
- else {
- if ((point->y[0] > Overlap.N)) { /*(3) */
- weight = (*point->y - Overlap.N) / overlap;
- *point->z = (1 - weight) * interpolation;
+ else if ((*point->x >= Overlap.W) && (*point->x <= Overlap.E)){
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ weight = (General.N - *point->y) / overlap;
+ *point->z = weight * interpolation;
sprintf(buf, "%lf", *point->z);
db_append_string(&sql, buf);
@@ -175,11 +165,11 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to access table <%s>"),
buf);
}
- else { /*(1) */
- weight = (Overlap.S - *point->y) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ weight = (*point->y - General.S) / overlap;
*point->z = (1 - weight) * interpolation;
sprintf(buf, "%lf", *point->z);
@@ -188,14 +178,17 @@
db_append_string(&sql, buf);
if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("Unable to create table: %s"),
+ G_fatal_error(_("Unable to access table <%s>"),
buf);
}
}
}
- }
- /*IF*/}
- /*FOR*/ return;
+ } /*IF*/
+ } /*FOR*/
+
+ db_commit_transaction(driver);
+
+ return;
}
@@ -207,13 +200,36 @@
int ncols, int bilin)
{
- int col, row;
+ int col, row, startcol, endcol, startrow, endrow;
double X, Y, interpolation, weight, csi, eta;
struct Cell_head Original;
G_get_window(&Original);
- for (row = 0; row < nrows; row++) {
- for (col = 0; col < ncols; col++) {
+ if (Original.north > General.N)
+ startrow = (Original.north - General.N) / Original.ns_res -1;
+ else
+ startrow = 0;
+ if (Original.north > General.S) {
+ endrow = (Original.north - General.S) / Original.ns_res + 1;
+ if (endrow > nrows)
+ endrow = nrows;
+ }
+ else
+ endrow = nrows;
+ if (General.W > Original.west)
+ startcol = (General.W - Original.west) / Original.ew_res - 1;
+ else
+ startcol = 0;
+ if (General.E > Original.west) {
+ endcol = (General.E - Original.west) / Original.ew_res + 1;
+ if (endcol > ncols)
+ endcol = ncols;
+ }
+ else
+ endcol = ncols;
+
+ for (row = startrow; row < endrow; row++) {
+ for (col = startcol; col < endcol; col++) {
X = G_col_to_easting((double)(col) + 0.5, &Original);
Y = G_row_to_northing((double)(row) + 0.5, &Original);
@@ -234,70 +250,59 @@
if (Vect_point_in_box(X, Y, interpolation, &Overlap)) { /* (5) */
matrix[row][col] = interpolation;
-
}
else {
-
- if ((X > Overlap.E)) {
-
- if ((Y > Overlap.N)) { /* (3) */
- csi = (X - Overlap.E) / overlap;
- eta = (Y - Overlap.N) / overlap;
- weight = (1 - csi) * (1 - eta);
+ if ((X > Overlap.E) && (X < General.E)) {
+ if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
+ csi = (General.E - X) / overlap;
+ eta = (General.N - Y) / overlap;
+ weight = csi * eta;
interpolation *= weight;
matrix[row][col] += interpolation;
-
}
- else if ((Y < Overlap.S)) { /* (1) */
- csi = (X - Overlap.E) / overlap;
+ else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
+ csi = (General.E - X) / overlap;
eta = (Y - General.S) / overlap;
- weight = (1 - csi) * eta;
+ weight = csi * eta;
interpolation *= weight;
matrix[row][col] = interpolation;
-
}
- else { /* (1) */
- weight = (X - Overlap.E) / overlap;
- interpolation *= 1 - weight;
+ else if ((Y >= Overlap.S) && (Y <= Overlap.N)) { /* (1) */
+ weight = (General.E - X ) / overlap;
+ interpolation *= weight;
matrix[row][col] = interpolation;
}
-
}
- else if ((X < Overlap.W)) {
-
- if ((Y > Overlap.N)) { /* (4) */
+ else if ((X < Overlap.W) && (X > General.W)) {
+ if ((Y > Overlap.N) && (Y < General.N)) { /* (4) */
csi = (X - General.W) / overlap;
- eta = (Y - Overlap.N) / overlap;
- weight = (1 - eta) * csi;
+ eta = (General.N - Y) / overlap;
+ weight = eta * csi;
interpolation *= weight;
matrix[row][col] += interpolation;
-
}
- else if ((Y < Overlap.S)) { /* (2) */
+ else if ((Y < Overlap.S) && (Y > General.S)) { /* (2) */
csi = (X - General.W) / overlap;
eta = (Y - General.S) / overlap;
weight = csi * eta;
interpolation *= weight;
matrix[row][col] += interpolation;
-
}
- else { /* (2) */
- weight = (Overlap.W - X) / overlap;
- interpolation *= 1 - weight;
+ else if ((Y >= Overlap.S) && (Y <= Overlap.N)) { /* (2) */
+ weight = (X - General.W) / overlap;
+ interpolation *= weight;
matrix[row][col] += interpolation;
}
-
}
- else {
- if ((Y > Overlap.N)) { /* (3) */
- weight = (Y - Overlap.N) / overlap;
- interpolation *= 1 - weight;
+ else if ((X >= Overlap.W) && (X <= Overlap.E)) {
+ if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
+ weight = (General.N - Y) / overlap;
+ interpolation *= weight;
matrix[row][col] += interpolation;
-
}
- else { /* (1) */
- weight = (Overlap.S - Y) / overlap;
- interpolation *= 1 - weight;
+ else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
+ weight = (Y - General.S) / overlap;
+ interpolation *= weight;
matrix[row][col] = interpolation;
}
}
Modified: grass/branches/develbranch_6/vector/lidar/lidarlib/zones.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/lidarlib/zones.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/lidarlib/zones.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -22,11 +22,44 @@
}
/*----------------------------------------------------------------------------------------*/
+/*
+ --------------------------------------------
+ | Elaboration region |
+ | ------------------------------------ |
+ | | General region | |
+ | | ---------------------------- | |
+ | | | | | |
+ | | | | | |
+ | | | | | |
+ | | | Overlap region | | |
+ | | | | | |
+ | | | | | |
+ | | | | | |
+ | | ---------------------------- | |
+ | | | |
+ | ------------------------------------ |
+ | |
+ --------------------------------------------
+
+ The terminology is misleading:
+ The Overlap region does NOT overlap with neighbouring segments,
+ but the Elaboration and General region do overlap
+
+ Elaboration is used for interpolation
+ Interpolated points in Elaboration but outside General are discarded
+ Interpolated points in General but outside Overlap are weighed by
+ their distance to Overlap and summed up
+ Interpolated points in Overlap are taken as they are
+
+ The buffer zones Elaboration - General and General - Overlap must be
+ large enough to avoid artefacts
+ */
+
int
P_set_regions(struct Cell_head *Elaboration, BOUND_BOX * General,
BOUND_BOX * Overlap, struct Reg_dimens dim, int type)
{
- /* Set the Elaborationoration region limits-> Also set the limits of the orlo and overlapping regions->
+ /* Set the Elaboration, General, and Overlap region limits
* Returns 0 on success; -1 on failure*/
struct Cell_head orig;
@@ -53,34 +86,34 @@
Overlap->E = General->E - dim.overlap;
return 0;
- case FIRST_ROW: /* It is just started with first row */
- Elaboration->north = orig.north;
+ case FIRST_ROW: /* Just started with first row */
+ Elaboration->north = orig.north + 2 * dim.orlo_h;
Elaboration->south = Elaboration->north - dim.latoN;
- General->N = Elaboration->north;
+ General->N = Elaboration->north - 2 * dim.orlo_h;
General->S = Elaboration->south + dim.orlo_h;
- Overlap->N = Elaboration->north;
+ Overlap->N = General->N;
Overlap->S = General->S + dim.overlap;
return 0;
- case LAST_ROW: /* It is reached last row */
- Elaboration->south = orig.south;
- Overlap->S = Elaboration->south;
- General->S = Elaboration->south;
+ case LAST_ROW: /* Reached last row */
+ Elaboration->south = orig.south - 2 * dim.orlo_h;
+ General->S = Elaboration->south + 2 * dim.orlo_h;
+ Overlap->S = General->S;
return 0;
- case FIRST_COLUMN: /* It is just started with first column */
- Elaboration->west = orig.west;
+ case FIRST_COLUMN: /* Just started with first column */
+ Elaboration->west = orig.west - 2 * dim.orlo_v;
Elaboration->east = Elaboration->west + dim.latoE;
- General->W = Elaboration->west;
+ General->W = Elaboration->west + 2 * dim.orlo_v;
General->E = Elaboration->east - dim.orlo_v;
- Overlap->W = Elaboration->west;
+ Overlap->W = General->W;
Overlap->E = General->E - dim.overlap;
return 0;
- case LAST_COLUMN: /* It is reached last column */
- Elaboration->east = orig.east;
- Overlap->E = Elaboration->east;
- General->E = Elaboration->east;
+ case LAST_COLUMN: /* Reached last column */
+ Elaboration->east = orig.east + 2 * dim.orlo_v;
+ General->E = Elaboration->east - 2 * dim.orlo_v;
+ Overlap->E = General->E;
return 0;
}
@@ -88,22 +121,104 @@
}
/*----------------------------------------------------------------------------------------*/
+int P_set_dim(struct Reg_dimens *dim, double pe, double pn, int *nsplx, int *nsply)
+{
+ int total_splines, orlo_splines, n_windows, minsplines;
+ double E_extension, N_extension, orloE, orloN;
+ struct Cell_head orig;
+ int ret = 0;
+
+ G_get_window(&orig);
+
+ E_extension = orig.east - orig.west;
+ N_extension = orig.north - orig.south;
+ dim->latoE = *nsplx * pe;
+ dim->latoN = *nsply * pn;
+ orloE = dim->latoE - dim->overlap - 2 * dim->orlo_v;
+ orloN = dim->latoN - dim->overlap - 2 * dim->orlo_h;
+
+ /* number of moving windows: E_extension / orloE */
+ /* remaining steps: total steps - (floor(E_extension / orloE) * E_extension) / passoE */
+ /* remaining steps must be larger than orlo_v + overlap + half of overlap window */
+ total_splines = ceil(E_extension / pe);
+ orlo_splines = orloE / pe;
+ n_windows = floor(E_extension / orloE); /* without last one */
+ if (n_windows > 0) {
+ minsplines = ceil((double)(dim->latoE / 2.0 - dim->orlo_v - dim->overlap) / pe);
+ while (total_splines - orlo_splines * n_windows < minsplines) {
+ *nsplx -= 1;
+ dim->latoE = *nsplx * pe;
+ orloE = dim->latoE - dim->overlap - 2 * dim->orlo_v;
+
+ orlo_splines = orloE / pe;
+ n_windows = floor(E_extension / orloE); /* without last one */
+ minsplines = ceil((double)(dim->latoE / 2.0 - dim->orlo_v - dim->overlap) / pe);
+ if (ret == 0)
+ ret = 1;
+ }
+ while (total_splines - orlo_splines * n_windows < minsplines * 2 && minsplines > 30) {
+ *nsplx -= 1;
+ dim->latoE = *nsplx * pe;
+ orloE = dim->latoE - dim->overlap - 2 * dim->orlo_v;
+
+ orlo_splines = orloE / pe;
+ n_windows = floor(E_extension / orloE); /* without last one */
+ minsplines = ceil((double)(dim->latoE / 2.0 - dim->orlo_v - dim->overlap) / pe);
+ if (ret == 0)
+ ret = 1;
+ }
+ }
+
+ total_splines = ceil(N_extension / pn);
+ orlo_splines = orloN / pn;
+ n_windows = floor(N_extension / orloN); /* without last one */
+ if (n_windows > 0) {
+ minsplines = ceil((double)(dim->latoN / 2.0 - dim->orlo_h - dim->overlap) / pn);
+ while (total_splines - orlo_splines * n_windows < minsplines) {
+ *nsply -= 1;
+ dim->latoN = *nsply * pn;
+ orloN = dim->latoN - dim->overlap - 2 * dim->orlo_h;
+
+ orlo_splines = orloN / pn;
+ n_windows = floor(N_extension / orloN); /* without last one */
+ minsplines = ceil((double)(dim->latoN / 2.0 - dim->orlo_h - dim->overlap) / pn);
+ if (ret < 2)
+ ret += 2;
+ }
+ while (total_splines - orlo_splines * n_windows < minsplines * 2 && minsplines > 30) {
+ *nsply -= 1;
+ dim->latoN = *nsply * pn;
+ orloN = dim->latoN - dim->overlap - 2 * dim->orlo_h;
+
+ orlo_splines = orloN / pn;
+ n_windows = floor(N_extension / orloN); /* without last one */
+ minsplines = ceil((double)(dim->latoN / 2.0 - dim->orlo_h - dim->overlap) / pn);
+ if (ret < 2)
+ ret += 2;
+ }
+ }
+ return 0;
+}
+
+/*----------------------------------------------------------------------------------------*/
int P_get_orlo(int interpolator, struct Reg_dimens *dim, double pe, double pn)
{
- /* Set the orlo regions dimension->
- * Returns 1 on success of bilinear; 2 on success of bicubic, 0 on failure-> */
- if (interpolator == 1) { /* the interpolator's bilinear */
- dim->orlo_v = 30 * pe; /*4 */
- dim->orlo_h = 30 * pn;
+ /* Set the orlo regions dimension
+ * Returns 1 on success of bilinear; 2 on success of bicubic, 0 on failure */
+ if (interpolator == P_BILINEAR) {
+ /* in case of edge artefacts, increase as multiples of 3 */
+ dim->orlo_v = 9 * pe;
+ dim->orlo_h = 9 * pn;
return 1;
}
- else if (interpolator == 0) { /* the interpolator's bicubic */
- dim->orlo_v = 40 * pe; /*3 */
- dim->orlo_h = 40 * pn;
+ else if (interpolator == P_BICUBIC) {
+ /* in case of edge artefacts, increase as multiples of 4 */
+ dim->orlo_v = 12 * pe; /*3 */
+ dim->orlo_h = 12 * pn;
return 2;
}
else
- return 0; /* The interpolator it's not bilinear nor bicubic!! */
+ return 0; /* The interpolator is neither bilinear nor bicubic!! */
}
/*----------------------------------------------------------------------------------------*/
@@ -147,12 +262,73 @@
return mean;
}
+
+double P_estimate_splinestep(struct Map_info *Map, double *dens, double *dist)
+{
+ int type, npoints = 0;
+ double xmin = 0, xmax = 0, ymin = 0, ymax = 0;
+ double x, y, z;
+ struct line_pnts *points;
+ struct line_cats *categories;
+ BOUND_BOX region_box;
+ struct Cell_head orig;
+
+ G_get_set_window(&orig);
+ Vect_region_box(&orig, ®ion_box);
+
+ points = Vect_new_line_struct();
+ categories = Vect_new_cats_struct();
+
+ Vect_rewind(Map);
+ while ((type = Vect_read_next_line(Map, points, categories)) > 0) {
+ if (!(type & GV_POINT))
+ continue;
+
+ x = points->x[0];
+ y = points->y[0];
+ if (points->z != NULL)
+ z = points->z[0];
+ else
+ z = 0.0;
+
+ /* only use points in current region */
+ if (Vect_point_in_box(x, y, z, ®ion_box)) {
+ npoints++;
+
+ if (npoints > 1) {
+ if (xmin > x)
+ xmin = x;
+ else if (xmax < x)
+ xmax = x;
+ if (ymin > y)
+ ymin = y;
+ else if (ymax < y)
+ ymax = y;
+ }
+ else {
+ xmin = xmax = x;
+ ymin = ymax = y;
+ }
+ }
+ }
+ if (npoints > 0) {
+ /* estimated average distance between points in map units */
+ *dist = sqrt(((xmax - xmin) * (ymax - ymin)) / npoints);
+ /* estimated point density as number of points per square map unit */
+ *dens = npoints / ((xmax - xmin) * (ymax - ymin));
+ return 0;
+ }
+ else {
+ return -1;
+ }
+}
+
struct Point *P_Read_Vector_Region_Map(struct Map_info *Map,
struct Cell_head *Elaboration,
int *num_points, int dim_vect,
int layer)
{
- int line_num, pippo, npoints, cat;
+ int line_num, pippo, npoints, cat, type;
double x, y, z;
struct Point *obs;
struct line_pnts *points;
@@ -165,15 +341,18 @@
points = Vect_new_line_struct();
categories = Vect_new_cats_struct();
- /* Reading the elaboration zone points */
+ /* Reading points inside elaboration zone */
Vect_region_box(Elaboration, &elaboration_box);
npoints = 0;
line_num = 0;
Vect_rewind(Map);
- while (Vect_read_next_line(Map, points, categories) > 0) {
+ while ((type = Vect_read_next_line(Map, points, categories)) > 0) {
+ if (!(type & GV_POINT))
+ continue;
+
line_num++;
x = points->x[0];
@@ -183,7 +362,7 @@
else
z = 0.0;
- /* Reading and storing points only if it is in elaboration_reg */
+ /* Reading and storing points only if in elaboration_reg */
if (Vect_point_in_box(x, y, z, &elaboration_box)) {
npoints++;
if (npoints >= pippo) {
@@ -212,42 +391,69 @@
}
/*------------------------------------------------------------------------------------------------*/
-int P_Create_Aux_Table(dbDriver * driver, char *tab_name)
+int P_Create_Aux2_Table(dbDriver * driver, char *tab_name)
{
dbTable *auxiliar_tab;
dbColumn *column;
- int created = FALSE;
- auxiliar_tab = db_alloc_table(4);
+ auxiliar_tab = db_alloc_table(2);
db_set_table_name(auxiliar_tab, tab_name);
db_set_table_description(auxiliar_tab,
"Intermediate interpolated values");
- column = db_get_table_column(auxiliar_tab, 2);
- db_set_column_name(column, "Y");
- db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+ column = db_get_table_column(auxiliar_tab, 0);
+ db_set_column_name(column, "ID");
+ db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
column = db_get_table_column(auxiliar_tab, 1);
- db_set_column_name(column, "X");
- db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+ db_set_column_name(column, "Interp");
+ db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
+ if (db_create_table(driver, auxiliar_tab) == DB_OK) {
+ G_debug(1, _("<%s> created in database."), tab_name);
+ return TRUE;
+ }
+ else
+ G_warning(_("<%s> has not been created in database."), tab_name);
+
+ return FALSE;
+}
+
+/*------------------------------------------------------------------------------------------------*/
+int P_Create_Aux4_Table(dbDriver * driver, char *tab_name)
+{
+ dbTable *auxiliar_tab;
+ dbColumn *column;
+
+ auxiliar_tab = db_alloc_table(4);
+ db_set_table_name(auxiliar_tab, tab_name);
+ db_set_table_description(auxiliar_tab,
+ "Intermediate interpolated values");
+
column = db_get_table_column(auxiliar_tab, 0);
db_set_column_name(column, "ID");
db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
- column = db_get_table_column(auxiliar_tab, 3);
+ column = db_get_table_column(auxiliar_tab, 1);
db_set_column_name(column, "Interp");
db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
+ column = db_get_table_column(auxiliar_tab, 2);
+ db_set_column_name(column, "X");
+ db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+
+ column = db_get_table_column(auxiliar_tab, 3);
+ db_set_column_name(column, "Y");
+ db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
+
if (db_create_table(driver, auxiliar_tab) == DB_OK) {
G_debug(1, _("<%s> created in database."), tab_name);
- created = TRUE;
- return created;
+ return TRUE;
}
else
- G_fatal_error(_("<%s> has not been created in database."), tab_name);
+ G_warning(_("<%s> has not been created in database."), tab_name);
- return created;
+ return FALSE;
}
/*------------------------------------------------------------------------------------------------*/
@@ -269,7 +475,6 @@
struct Cell_head Original;
G_get_window(&Original);
- G_set_window(&Original);
nrows = G_window_rows();
ncols = G_window_cols();
@@ -286,6 +491,7 @@
}
G_put_d_raster_row(fd, raster);
}
+ G_percent(row, nrows, 2);
}
/*------------------------------------------------------------------------------------------------*/
@@ -294,7 +500,7 @@
char *tab_name)
{
- int more, ltype, line_num, type, count = 0;
+ int more, line_num, type, count = 0;
double coordX, coordY, coordZ;
struct line_pnts *point;
@@ -337,7 +543,7 @@
value = db_get_column_value(column);
else
continue;
- coordX = db_get_value_double(value);
+ coordZ = db_get_value_double(value);
column = db_get_table_column(table, 2);
type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
@@ -345,7 +551,7 @@
value = db_get_column_value(column);
else
continue;
- coordY = db_get_value_double(value);
+ coordX = db_get_value_double(value);
column = db_get_table_column(table, 3);
type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
@@ -353,7 +559,7 @@
value = db_get_column_value(column);
else
continue;
- coordZ = db_get_value_double(value);
+ coordY = db_get_value_double(value);
Vect_copy_xyz_to_pnts(point, &coordX, &coordY, &coordZ, 1);
Vect_reset_cats(cat);
@@ -372,7 +578,6 @@
G_get_window(&Original);
- G_set_window(&Original);
nrows = G_window_rows();
ncols = G_window_cols();
@@ -384,10 +589,30 @@
}
#endif
-/*------------------------------------------------------------------------------------------------*/
-#ifdef notdef
-/* Subzones definition */
------------------------|4 | 3 | 3 | ||-----------------------|||||||2 | 5 | 1
- | ||||||||-----------------------|2 | 1 | 1 |
- ||-----------------------||||||||||||||||||-----------------------||||||-----------------------
-#endif
+/*! DEFINITION OF THE SUBZONES
+
+ 5: inside Overlap region
+ all others: inside General region but outside Overlap region
+
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | |4| 3 |3| | |
+ ---------------------------------
+ | | | | | | | |
+ | | |2| 5 |1| | |
+ | | | | | | | |
+ ---------------------------------
+ | | |2| 1 |1| | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ */
Modified: grass/branches/develbranch_6/vector/lidar/v.lidar.correction/correction.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.lidar.correction/correction.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.lidar.correction/correction.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -20,6 +20,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include <grass/gis.h>
#include <grass/Vect.h>
#include <grass/dbmi.h>
@@ -29,183 +30,198 @@
void
P_Sparse_Correction(struct Map_info *In, struct Map_info *Out,
struct Map_info *Terrain, struct Cell_head *Elaboration,
- BOUND_BOX General, BOUND_BOX Overlap, double **obs,
- double *param, int *line_num, double passoN,
- double passoE, double overlap, double HighThresh,
- double LowThresh, int nsplx, int nsply, int num_points,
- dbDriver * driver, double mean)
+ BOUND_BOX General, BOUND_BOX Overlap,
+ double **obs, struct lidar_cat *lcat, double *param,
+ int *line_num, double passoN, double passoE,
+ double overlap, double HighThresh, double LowThresh,
+ int nsplx, int nsply, int num_points,
+ dbDriver * driver, double mean, char *tab_name)
{
int i = 0, class;
double interpolation, csi, eta, weight;
- struct line_pnts *points;
+ BOUND_BOX elaboration_box;
+ struct line_pnts *point;
struct line_cats *cats;
- points = Vect_new_line_struct();
+ point = Vect_new_line_struct();
cats = Vect_new_cats_struct();
- Vect_rewind(In);
- while (Vect_read_next_line(In, points, cats) > 0) {
- i++;
- if (Vect_point_in_box(*points->x, *points->y, *points->z, &General)) {
+ Vect_region_box(Elaboration, &elaboration_box);
+
+ db_begin_transaction(driver);
+
+ for (i = 0; i < num_points; i++) { /* Sparse points */
+ G_percent(i, num_points, 2);
+ Vect_reset_line(point);
+ Vect_reset_cats(cats);
+
+ if (Vect_point_in_box(obs[i][0], obs[i][1], mean, &General)) {
interpolation =
- dataInterpolateBilin(*points->x, *points->y, passoE, passoN,
+ dataInterpolateBilin(obs[i][0], obs[i][1], passoE, passoN,
nsplx, nsply, Elaboration->west,
Elaboration->south, param);
interpolation += mean;
- if (Vect_point_in_box(*points->x, *points->y, *points->z, &Overlap)) { /*(5) */
+ Vect_copy_xyz_to_pnts(point, &obs[i][0], &obs[i][1], &obs[i][2],
+ 1);
+ *point->z += mean;
- Vect_cat_get(cats, F_CLASSIFICATION, &class);
+ if (Vect_point_in_box(*point->x, *point->y, *point->z, &Overlap)) { /*(5) */
+
+ Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+ Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+ Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+ class = lcat[i].cat_class;
class =
- correction(class, *points->z, interpolation, HighThresh,
+ correction(class, *point->z, interpolation, HighThresh,
LowThresh);
- Vect_cat_del(cats, F_CLASSIFICATION);
Vect_cat_set(cats, F_CLASSIFICATION, class);
if (class == TERRAIN_SINGLE || class == TERRAIN_DOUBLE)
- Vect_write_line(Terrain, GV_POINT, points, cats);
+ Vect_write_line(Terrain, GV_POINT, point, cats);
- Vect_write_line(Out, GV_POINT, points, cats);
-
+ Vect_write_line(Out, GV_POINT, point, cats);
}
else {
-
- if ((*points->x > Overlap.E)) {
-
- if ((*points->y > Overlap.N)) { /*(3) */
- csi = (*points->x - Overlap.E) / overlap;
- eta = (*points->y - Overlap.N) / overlap;
- weight = (1 - csi) * (1 - eta);
+ if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = csi * eta;
interpolation *= weight;
if (Select_Correction
- (&interpolation, line_num[i], driver) != DB_OK)
+ (&interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
if (UpDate_Correction
- (interpolation, line_num[i], driver) != DB_OK)
+ (interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to update the database"));
-
}
- else if ((*points->y < Overlap.S)) { /*(1) */
- csi = (*points->x - Overlap.E) / overlap;
- eta = (*points->y - General.S) / overlap;
- weight = (1 - csi) * eta;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (*point->y - General.S) / overlap;
+ weight = csi * eta;
if (Insert_Correction
(interpolation * weight, line_num[i],
- driver) != DB_OK)
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to write in the database"));
-
}
- else { /*(1) */
- weight = (*points->x - Overlap.E) / overlap;
+ else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) { /*(1) */
+ weight = (General.E - *point->x) / overlap;
if (Insert_Correction
(interpolation * weight, line_num[i],
- driver) != DB_OK)
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to write in the database"));
}
}
- else if ((*points->x < Overlap.W)) {
- if ((*points->y > Overlap.N)) { /*(4) */
- csi = (*points->x - General.W) / overlap;
- eta = (*points->y - Overlap.N) / overlap;
- weight = (1 - eta) * csi;
+ else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(4) */
+ csi = (*point->x - General.W) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = eta * csi;
+ interpolation *= weight;
- interpolation *= weight;
if (Select_Correction
- (&interpolation, line_num[i], driver) != DB_OK)
+ (&interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
- Vect_cat_get(cats, F_CLASSIFICATION, &class);
+ Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+ Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+ Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+ class = lcat[i].cat_class;
class =
- correction(class, *points->z, interpolation,
+ correction(class, *point->z, interpolation,
HighThresh, LowThresh);
Vect_cat_set(cats, F_CLASSIFICATION, class);
- Vect_write_line(Out, GV_POINT, points, cats);
+ Vect_write_line(Out, GV_POINT, point, cats);
if (class == TERRAIN_SINGLE ||
class == TERRAIN_DOUBLE)
- Vect_write_line(Terrain, GV_POINT, points, cats);
-
+ Vect_write_line(Terrain, GV_POINT, point, cats);
}
- else if ((*points->y < Overlap.S)) { /*(2) */
- csi = (*points->x - General.W) / overlap;
- eta = (*points->y - General.S) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(2) */
+ csi = (*point->x - General.W) / overlap;
+ eta = (*point->y - General.S) / overlap;
weight = csi * eta;
interpolation *= weight;
if (Select_Correction
- (&interpolation, line_num[i], driver) != DB_OK)
+ (&interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
if (UpDate_Correction
- (interpolation, line_num[i], driver) != DB_OK)
+ (interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to update the database"));
}
- else { /*(2) */
- weight = (Overlap.W - *points->x) / overlap;
+ else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) { /*(2) */
+ weight = (*point->x - General.W) / overlap;
interpolation *= weight;
if (Select_Correction
- (&interpolation, line_num[i], driver) != DB_OK)
+ (&interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
- Vect_cat_get(cats, F_CLASSIFICATION, &class);
+ Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+ Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+ Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+ class = lcat[i].cat_class;
class =
- correction(class, *points->z, interpolation,
+ correction(class, *point->z, interpolation,
HighThresh, LowThresh);
Vect_cat_set(cats, F_CLASSIFICATION, class);
- Vect_write_line(Out, GV_POINT, points, cats);
+ Vect_write_line(Out, GV_POINT, point, cats);
if (class == TERRAIN_SINGLE ||
class == TERRAIN_DOUBLE)
- Vect_write_line(Terrain, GV_POINT, points, cats);
-
+ Vect_write_line(Terrain, GV_POINT, point, cats);
}
-
}
- else {
- if ((*points->y > Overlap.N)) { /*(3) */
- weight = (*points->y - Overlap.N) / overlap;
+ else if ((*point->x >= Overlap.W) && (*point->x <= Overlap.E)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ weight = (General.N - *point->y) / overlap;
interpolation *= weight;
if (Select_Correction
- (&interpolation, line_num[i], driver) != DB_OK)
+ (&interpolation, line_num[i], driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
- Vect_cat_get(cats, F_CLASSIFICATION, &class);
+ Vect_cat_set(cats, F_EDGE_DETECTION_CLASS, lcat[i].cat_edge);
+ Vect_cat_set(cats, F_INTERPOLATION, lcat[i].cat_interp);
+ Vect_cat_set(cats, F_COUNTER_OBJ, lcat[i].cat_obj);
+ class = lcat[i].cat_class;
class =
- correction(class, *points->z, interpolation,
+ correction(class, *point->z, interpolation,
HighThresh, LowThresh);
Vect_cat_set(cats, F_CLASSIFICATION, class);
- Vect_write_line(Out, GV_POINT, points, cats);
+ Vect_write_line(Out, GV_POINT, point, cats);
if (class == TERRAIN_SINGLE ||
class == TERRAIN_DOUBLE)
- Vect_write_line(Terrain, GV_POINT, points, cats);
+ Vect_write_line(Terrain, GV_POINT, point, cats);
}
- else { /*(1) */
- weight = (Overlap.S - *points->y) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ weight = (*point->y - General.S) / overlap;
if (Insert_Correction
(interpolation * weight, line_num[i],
- driver) != DB_OK)
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to write in the database"));
}
}
}
- }
- /*IF*/ Vect_reset_line(points);
- Vect_reset_cats(cats);
-
- }
- /*FOR*/ Vect_destroy_line_struct(points);
+ } /* if in General box */
+ } /* while */
+ G_percent(num_points, num_points, 2);
+ Vect_destroy_line_struct(point);
Vect_destroy_cats_struct(cats);
+ db_commit_transaction(driver);
+
return;
}
@@ -228,7 +244,7 @@
return class;
}
-int Select_Correction(double *Interp, int line_num, dbDriver * driver)
+int Select_Correction(double *Interp, int line_num, dbDriver * driver, char *tab_name)
{
int more;
char buf[1024];
@@ -239,10 +255,8 @@
dbValue *Interp_value;
db_init_string(&sql);
- db_zero_string(&sql);
- sprintf(buf, "SELECT Interp FROM Auxiliar_correction_table WHERE ID=%d",
- line_num);
+ sprintf(buf, "SELECT ID, Interp FROM %s WHERE ID=%d", tab_name, line_num);
db_append_string(&sql, buf);
if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
@@ -260,55 +274,61 @@
*Interp += db_get_value_double(Interp_value);
}
+ db_close_cursor(&cursor);
+ db_free_string(&sql);
return DB_OK;
}
-int Insert_Correction(double Interp, int line_num, dbDriver * driver)
+int Insert_Correction(double Interp, int line_num, dbDriver * driver, char *tab_name)
{
char buf[1024];
dbString sql;
+ int ret;
db_init_string(&sql);
- sprintf(buf, "INSERT INTO Auxiliar_correction_table (ID, Interp)");
+ sprintf(buf, "INSERT INTO %s (ID, Interp)", tab_name);
db_append_string(&sql, buf);
sprintf(buf, " VALUES (%d, %lf)", line_num, Interp);
db_append_string(&sql, buf);
- return db_execute_immediate(driver, &sql);
+ ret = db_execute_immediate(driver, &sql);
+ db_free_string(&sql);
+ return ret;
}
-int UpDate_Correction(double Interp, int line_num, dbDriver * driver)
+int UpDate_Correction(double Interp, int line_num, dbDriver * driver, char *tab_name)
{
char buf[1024];
dbString sql;
+ int ret;
db_init_string(&sql);
sprintf(buf,
- "UPDATE Auxiliar_correction_table SET Interp=%lf WHERE ID=%d",
- Interp, line_num);
+ "UPDATE %s SET Interp=%lf WHERE ID=%d", tab_name, Interp, line_num);
db_append_string(&sql, buf);
- return db_execute_immediate(driver, &sql);
+ ret = db_execute_immediate(driver, &sql);
+ db_free_string(&sql);
+ return ret;
}
struct Point *P_Read_Vector_Correction(struct Map_info *Map,
struct Cell_head *Elaboration,
int *num_points, int *num_terrain,
- int dim_vect)
+ int dim_vect, struct lidar_cat **lcat)
{
- int line_num, npoints, nterrain, pippo, cat_edge;
+ int line_num, npoints, nterrain, pippo, cat;
double x, y, z;
struct Point *obs;
+ struct lidar_cat *lc;
struct line_pnts *points;
struct line_cats *categories;
BOUND_BOX elaboration_box;
pippo = dim_vect;
- /*if (first_it) */
obs = (struct Point *)G_calloc(pippo, sizeof(struct Point));
- /*else */
- /*obs = (struct Point*) G_realloc ((void *) obs, pippo * sizeof(struct Point)); */
+ lc = (struct lidar_cat *)G_calloc(pippo, sizeof(struct lidar_cat));
points = Vect_new_line_struct();
categories = Vect_new_cats_struct();
@@ -338,12 +358,10 @@
npoints++;
if (npoints >= pippo) {
pippo += dim_vect;
- obs =
- (struct Point *)G_realloc((void *)obs,
- (signed int)(pippo *
- sizeof(struct
- Point)));
- /*lineID = (int *) G_realloc ((void *) lineID, pippo * sizeof(int)); */
+ obs = (struct Point *)G_realloc((void *)obs,
+ (signed int)(pippo * sizeof(struct Point)));
+ lc = (struct lidar_cat *)G_realloc((void *)lc,
+ (signed int)(pippo * sizeof(struct lidar_cat)));
}
/* Storing observation vector */
@@ -351,27 +369,55 @@
obs[npoints - 1].coordY = y;
obs[npoints - 1].coordZ = z;
obs[npoints - 1].lineID = line_num; /* Storing also the line's number */
- Vect_cat_get(categories, F_EDGE_DETECTION_CLASS, &cat_edge);
- obs[npoints - 1].cat = cat_edge;
+ Vect_cat_get(categories, F_EDGE_DETECTION_CLASS, &cat);
+ obs[npoints - 1].cat = cat;
+ lc[npoints - 1].cat_edge = cat;
+
+ /* Only terrain points */
+ if (cat == TERRAIN_SINGLE)
+ nterrain++;
+
+ Vect_cat_get(categories, F_CLASSIFICATION, &cat);
+ lc[npoints - 1].cat_class = cat;
+ Vect_cat_get(categories, F_INTERPOLATION, &cat);
+ lc[npoints - 1].cat_interp = cat;
+ Vect_cat_get(categories, F_COUNTER_OBJ, &cat);
+ lc[npoints - 1].cat_obj = cat;
}
-
- /* Only terrain points */
- if (cat_edge == TERRAIN_SINGLE)
- nterrain++;
}
Vect_destroy_line_struct(points);
Vect_destroy_cats_struct(categories);
*num_points = npoints;
*num_terrain = nterrain;
+ *lcat = lc;
return obs;
}
+/*! DEFINITION OF THE SUBZONES
+ 5: inside Overlap region
+ all others: inside General region but outside Overlap region
-#ifdef notdef
-/* SUBZONES DEFINITION */
------------------------|4 | 3 | 3 | ||-----------------------|||||||2 | 5 | 1
- | ||||||||-----------------------|2 | 1 | 1 |
- ||-----------------------||||||||||||||||||-----------------------||||||-----------------------
-#endif
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | |4| 3 |3| | |
+ ---------------------------------
+ | | | | | | | |
+ | | |2| 5 |1| | |
+ | | | | | | | |
+ ---------------------------------
+ | | |2| 1 |1| | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ */
Modified: grass/branches/develbranch_6/vector/lidar/v.lidar.correction/correction.h
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.lidar.correction/correction.h 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.lidar.correction/correction.h 2010-02-17 17:56:52 UTC (rev 41062)
@@ -7,6 +7,14 @@
#include <grass/PolimiFunct.h>
+struct lidar_cat
+{
+ int cat_edge; /* category in layer F_EDGE_DETECTION_CLASS */
+ int cat_class; /* category in layer F_CLASSIFICATION */
+ int cat_interp; /* category in layer F_INTERPOLATION */
+ int cat_obj; /* category in layer F_COUNTER_OBJ */
+};
+
void
P_Sparse_Correction(struct Map_info *, /**/
struct Map_info *, /**/
@@ -15,6 +23,7 @@
BOUND_BOX, /**/
BOUND_BOX, /**/
double **, /**/
+ struct lidar_cat *,
double *, /**/
int *, /**/
double, /**/
@@ -23,19 +32,21 @@
double, /**/
double, /**/
int, /**/
- int, /**/ int, /**/ dbDriver *, /**/ double /**/);
+ int, /**/
+ int, /**/
+ dbDriver *, /**/
+ double, /**/
+ char *);
-int Insert_Correction(double, /**/ int, /**/ dbDriver * /**/);
+int Insert_Correction(double, /**/ int, /**/ dbDriver *, /**/ char *);
-int UpDate_Correction(double, /**/ int, /**/ dbDriver * /**/);
+int UpDate_Correction(double, /**/ int, /**/ dbDriver *, /**/ char *);
-int Select_Correction(double *, /**/ int, /**/ dbDriver * /**/);
+int Select_Correction(double *, /**/ int, /**/ dbDriver *, /**/ char *);
-int Create_AuxEdge_Table(dbDriver *);
-
struct Point *P_Read_Vector_Correction(struct Map_info *, /**/
struct Cell_head *, /**/
- int *, /**/ int *, /**/ int /**/);
+ int *, /**/ int *, /**/ int /**/,
+ struct lidar_cat **);
-void P_Aux_to_Raster(double **, int);
int correction(int, double, double, double, double);
Modified: grass/branches/develbranch_6/vector/lidar/v.lidar.correction/main.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.lidar.correction/main.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.lidar.correction/main.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -20,6 +20,7 @@
/*INCLUDES*/
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include <grass/config.h>
#include <grass/gis.h>
#include <grass/Vect.h>
@@ -32,10 +33,16 @@
int main(int argc, char *argv[])
{
/* Declarations */
- int dim_vect, nparameters, BW, npoints, nrows, ncols, nsply, nsplx;
- char *dvr, *db, *mapset, table_name[1024];
+ int dim_vect, nparameters, BW, npoints, nrows, ncols;
+ int nsply, nsplx, nsplx_adj, nsply_adj;
+ int nsubregion_col, nsubregion_row;
+ int subregion = 0, nsubregions = 0;
+ const char *dvr, *db, *mapset;
+ char table_name[GNAME_MAX];
+ char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
double lambda, ew_resol, ns_resol, mean, passoN, passoE, HighThresh,
LowThresh;
+ double N_extension, E_extension, orloE, orloN;
int i, nterrain, count_terrain;
@@ -43,11 +50,12 @@
int *lineVect;
double *TN, *Q, *parVect; /* Interpolating and least-square vectors */
- double **N, **obsVect; /* Interpolation and least-square matrix */
+ double **N, **obsVect, **obsVect_all; /* Interpolation and least-square matrix */
struct Map_info In, Out, Terrain;
struct Option *in_opt, *out_opt, *out_terrain_opt, *passoE_opt,
*passoN_opt, *lambda_f_opt, *Thresh_A_opt, *Thresh_B_opt;
+ struct Flag *spline_step_flag;
struct GModule *module;
struct Cell_head elaboration_reg, original_reg;
@@ -55,6 +63,7 @@
BOUND_BOX general_box, overlap_box;
struct Point *observ;
+ struct lidar_cat *lcat;
dbDriver *driver;
@@ -65,6 +74,12 @@
module->description =
_("Correction of the v.lidar.growing output. It is the last of the three algorithms for LIDAR filtering.");
+ spline_step_flag = G_define_flag();
+ spline_step_flag->key = 'e';
+ spline_step_flag->label = _("Estimate point density and distance");
+ spline_step_flag->description =
+ _("Estimate point density and distance for the input vector points within the current region extends and quit");
+
in_opt = G_define_standard_option(G_OPT_V_INPUT);
in_opt->description =
_("Input observation vector map name (v.lidar.growing output)");
@@ -132,12 +147,32 @@
lambda = atof(lambda_f_opt->answer);
HighThresh = atof(Thresh_A_opt->answer);
LowThresh = atof(Thresh_B_opt->answer);
- dvr = G__getenv2("DB_DRIVER", G_VAR_MAPSET);
- db = G__getenv2("DB_DATABASE", G_VAR_MAPSET);
+ if (!(db = G__getenv2("DB_DATABASE", G_VAR_MAPSET)))
+ G_fatal_error(_("Unable to read name of database"));
+
+ if (!(dvr = G__getenv2("DB_DRIVER", G_VAR_MAPSET)))
+ G_fatal_error(_("Unable to read name of driver"));
+
/* Setting auxiliar table's name */
- sprintf(table_name, "%s_aux", out_opt->answer);
+ if (G__name_is_fully_qualified(out_opt->answer, xname, xmapset)) {
+ sprintf(table_name, "%s_aux", xname);
+ }
+ else
+ sprintf(table_name, "%s_aux", out_opt->answer);
+ /* Something went wrong in a previous v.lidar.correction execution */
+ if (db_table_exists(dvr, db, table_name)) {
+ /* Start driver and open db */
+ driver = db_start_driver_open_database(dvr, db);
+ if (driver == NULL)
+ G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
+ dvr);
+ if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
+ G_fatal_error(_("Old auxiliar table could not be dropped"));
+ db_close_database_shutdown_driver(driver);
+ }
+
/* Checking vector names */
Vect_check_input_output_name(in_opt->answer, out_opt->answer,
GV_FATAL_EXIT);
@@ -146,10 +181,28 @@
if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL)
G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
- Vect_set_open_level(1); /*without topology */
+ Vect_set_open_level(1); /* without topology */
if (1 > Vect_open_old(&In, in_opt->answer, mapset))
G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
+ /* Input vector must be 3D */
+ if (!Vect_is_3d(&In))
+ G_fatal_error(_("Input vector map <%s> is not 3D!"), in_opt->answer);
+
+ /* Estimate point density and mean distance for current region */
+ if (spline_step_flag->answer) {
+ double dens, dist;
+ if (P_estimate_splinestep(&In, &dens, &dist) == 0) {
+ G_message("Estimated point density: %.4g", dens);
+ G_message("Estimated mean distance between points: %.4g", dist);
+ }
+ else
+ G_warning(_("No points in current region!"));
+
+ Vect_close(&In);
+ exit(EXIT_SUCCESS);
+ }
+
/* Open output vector */
if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z)) {
Vect_close(&In);
@@ -176,32 +229,75 @@
G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
dvr);
+ /* Create auxiliar table */
+ if ((flag_auxiliar =
+ P_Create_Aux2_Table(driver, table_name)) == FALSE) {
+ Vect_close(&In);
+ Vect_close(&Out);
+ Vect_close(&Terrain);
+ exit(EXIT_FAILURE);
+ }
+
+ db_create_index2(driver, table_name, "ID");
+ /* sqlite likes that */
+ db_close_database_shutdown_driver(driver);
+ driver = db_start_driver_open_database(dvr, db);
+
/* Setting regions and boxes */
G_get_set_window(&original_reg);
G_get_set_window(&elaboration_reg);
Vect_region_box(&elaboration_reg, &overlap_box);
Vect_region_box(&elaboration_reg, &general_box);
- /* Fixxing parameters of the elaboration region */
- /*! Each original_region will be divided into several subregions. These
- * subregion will be overlaped by its neibourgh subregions. This overlapping
- * is calculated as OVERLAP_PASS times the east-west resolution. */
-
nrows = G_window_rows();
ncols = G_window_cols();
ew_resol = original_reg.ew_res;
ns_resol = original_reg.ns_res;
+ /*------------------------------------------------------------------
+ | Subdividing and working with tiles:
+ | Each original region will be divided into several subregions.
+ | Each one will be overlaped by its neighbouring subregions.
+ | The overlapping is calculated as a fixed OVERLAP_SIZE times
+ | the largest spline step plus 2 * orlo
+ ----------------------------------------------------------------*/
+
+ /* Fixing parameters of the elaboration region */
P_zero_dim(&dims);
- dims.latoE = NSPLX_MAX * passoE;
- dims.latoN = NSPLY_MAX * passoN;
- dims.overlap = OVERLAP_SIZE * ew_resol;
+
+ nsplx_adj = NSPLX_MAX;
+ nsply_adj = NSPLY_MAX;
+ if (passoN > passoE)
+ dims.overlap = OVERLAP_SIZE * passoN;
+ else
+ dims.overlap = OVERLAP_SIZE * passoE;
P_get_orlo(P_BILINEAR, &dims, passoE, passoN);
+ P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
- /* Subdividing and working with tiles */
+ G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
+ G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
+
+ /* calculate number of subregions */
+ orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
+ orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+
+ N_extension = original_reg.north - original_reg.south;
+ E_extension = original_reg.east - original_reg.west;
+
+ nsubregion_col = ceil(E_extension / orloE) + 0.5;
+ nsubregion_row = ceil(N_extension / orloN) + 0.5;
+
+ if (nsubregion_col < 0)
+ nsubregion_col = 0;
+ if (nsubregion_row < 0)
+ nsubregion_row = 0;
+
+ nsubregions = nsubregion_row * nsubregion_col;
+
elaboration_reg.south = original_reg.north;
last_row = FALSE;
+
while (last_row == FALSE) { /* For each row */
P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
@@ -219,11 +315,13 @@
}
nsply =
- ceil((elaboration_reg.north - elaboration_reg.south) / passoN) +
- 1;
+ ceil((elaboration_reg.north -
+ elaboration_reg.south) / passoN) + 0.5;
+ /*
if (nsply > NSPLY_MAX) {
nsply = NSPLY_MAX;
}
+ */
G_debug(1, _("nsply = %d"), nsply);
elaboration_reg.east = original_reg.west;
@@ -231,6 +329,10 @@
while (last_column == FALSE) { /* For each column */
+ subregion++;
+ if (nsubregions > 1)
+ G_message(_("subregion %d of %d"), subregion, nsubregions);
+
P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
GENERAL_COLUMN);
@@ -247,33 +349,36 @@
nsplx =
ceil((elaboration_reg.east - elaboration_reg.west) / passoE) +
- 1;
+ 0.5;
+ /*
if (nsplx > NSPLX_MAX) {
nsplx = NSPLX_MAX;
}
+ */
G_debug(1, _("nsplx = %d"), nsplx);
dim_vect = nsplx * nsply;
G_debug(1, _("read vector region map"));
observ =
P_Read_Vector_Correction(&In, &elaboration_reg, &npoints,
- &nterrain, dim_vect);
+ &nterrain, dim_vect, &lcat);
G_debug(5, _("npoints = %d, nterrain = %d"), npoints, nterrain);
if (npoints > 0) { /* If there is any point falling into elaboration_reg. */
count_terrain = 0;
nparameters = nsplx * nsply;
- /* Mean's calculation */
- G_debug(3, _("Mean's calculation"));
+ /* Mean calculation */
+ G_debug(3, _("Mean calculation"));
mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
/*Least Squares system */
BW = P_get_BandWidth(P_BILINEAR, nsply); /* Bilinear interpolation */
N = G_alloc_matrix(nparameters, BW); /* Normal matrix */
TN = G_alloc_vector(nparameters); /* vector */
- parVect = G_alloc_vector(nparameters); /* Bicubic parameters vector */
- obsVect = G_alloc_matrix(nterrain + 1, 3); /* Observation vector */
+ parVect = G_alloc_vector(nparameters); /* Bilinear parameters vector */
+ obsVect = G_alloc_matrix(nterrain + 1, 3); /* Observation vector with terrain points */
+ obsVect_all = G_alloc_matrix(npoints + 1, 3); /* Observation vector with all points */
Q = G_alloc_vector(nterrain + 1); /* "a priori" var-cov matrix */
lineVect = G_alloc_ivector(npoints + 1);
@@ -288,9 +393,14 @@
count_terrain++;
}
lineVect[i] = observ[i].lineID;
+ obsVect_all[i][0] = observ[i].coordX;
+ obsVect_all[i][1] = observ[i].coordY;
+ obsVect_all[i][2] = observ[i].coordZ - mean;
}
- G_debug(3, _("M.Q. solution"));
+ G_free(observ);
+
+ G_verbose_message(_("Bilinear interpolation"));
normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
nsply, elaboration_reg.west,
elaboration_reg.south, nterrain, nparameters,
@@ -301,30 +411,25 @@
G_free_matrix(N);
G_free_vector(TN);
G_free_vector(Q);
+ G_free_matrix(obsVect);
- if (flag_auxiliar == FALSE) {
- if ((flag_auxiliar =
- P_Create_Aux_Table(driver, table_name)) == FALSE) {
- Vect_close(&In);
- Vect_close(&Out);
- Vect_close(&Terrain);
- exit(EXIT_FAILURE);
- }
- }
-
- G_debug(3, _("Correction and creation of terrain vector"));
+ G_verbose_message( _("Correction and creation of terrain vector"));
P_Sparse_Correction(&In, &Out, &Terrain, &elaboration_reg,
- general_box, overlap_box, obsVect,
+ general_box, overlap_box, obsVect_all, lcat,
parVect, lineVect, passoN, passoE,
dims.overlap, HighThresh, LowThresh,
- nsplx, nsply, npoints, driver, mean);
+ nsplx, nsply, npoints, driver, mean, table_name);
- G_free_matrix(obsVect);
+ G_free_vector(parVect);
+ G_free_matrix(obsVect_all);
G_free_ivector(lineVect);
-
- G_free_vector(parVect);
}
- G_free(observ);
+ else {
+ G_free(observ);
+ G_warning(_("No data within this subregion. "
+ "Consider changing the spline step."));
+ }
+ G_free(lcat);
} /*! END WHILE; last_column = TRUE */
} /*! END WHILE; last_row = TRUE */
Modified: grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/edgedetection.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/edgedetection.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/edgedetection.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -26,11 +26,13 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include <grass/gis.h>
#include <grass/Vect.h>
#include <grass/dbmi.h>
/* #include <grass/PolimiFunct.h> */
#include "edgedetection.h"
+
int edge_detection(struct Cell_head elaboration_reg, BOUND_BOX Overlap_Box,
double *parBilin, double obsX, double obsY,
double *partial, double alpha, double residual,
@@ -41,14 +43,14 @@
/* 3 = PRE_UNKNOWN */
int c1, c2;
- double g[9][2], *gradient, gradPto, dirPto;
+ double g[9][2], gradient[2], gradPto, dirPto;
extern double passoE, passoN;
static struct Cell_head Elaboration;
g[0][0] = partial[0];
g[0][1] = partial[1];
- gradPto = sqrt(g[0][0] * g[0][0] + g[0][1] * g[0][1]);
+ gradPto = g[0][0] * g[0][0] + g[0][1] * g[0][1];
dirPto = atan(g[0][1] / g[0][0]) + M_PI / 2; /* radiants */
Elaboration = elaboration_reg;
@@ -59,72 +61,64 @@
else if ((gradPto > gradLow) && (residual > 0)) { /* Soft condition for 'edge' points */
if (Vect_point_in_box(obsX, obsY, 0.0, &Overlap_Box)) {
- gradient =
- Get_Gradient(Elaboration, obsX + passoE * cos(dirPto),
- obsY + passoN * sin(dirPto), parBilin);
+ Get_Gradient(Elaboration, obsX + passoE * cos(dirPto),
+ obsY + passoN * sin(dirPto), parBilin, gradient);
g[2][0] = gradient[0];
g[2][1] = gradient[1];
- gradient =
- Get_Gradient(Elaboration, obsX + passoE * cos(dirPto + M_PI),
- obsY + passoN * sin(dirPto + M_PI), parBilin);
+ Get_Gradient(Elaboration, obsX + passoE * cos(dirPto + M_PI),
+ obsY + passoN * sin(dirPto + M_PI), parBilin, gradient);
g[7][0] = gradient[0];
g[7][1] = gradient[1];
if ((fabs(atan(g[2][1] / g[2][0]) + M_PI / 2 - dirPto) < alpha) &&
(fabs(atan(g[7][1] / g[7][0]) + M_PI / 2 - dirPto) < alpha)) {
- gradient =
- Get_Gradient(Elaboration,
+ Get_Gradient(Elaboration,
obsX + passoE * cos(dirPto + M_PI / 4),
obsY + passoN * sin(dirPto + M_PI / 4),
- parBilin);
+ parBilin, gradient);
g[1][0] = gradient[0];
g[1][1] = gradient[1];
- gradient =
- Get_Gradient(Elaboration,
+ Get_Gradient(Elaboration,
obsX + passoE * cos(dirPto - M_PI / 4),
obsY + passoN * sin(dirPto - M_PI / 4),
- parBilin);
+ parBilin, gradient);
g[3][0] = gradient[0];
g[3][1] = gradient[1];
- gradient =
- Get_Gradient(Elaboration,
+ Get_Gradient(Elaboration,
obsX + passoE * cos(dirPto + M_PI / 2),
obsY + passoN * sin(dirPto + M_PI / 2),
- parBilin);
+ parBilin, gradient);
g[4][0] = gradient[0];
g[4][1] = gradient[1];
- gradient =
- Get_Gradient(Elaboration,
+ Get_Gradient(Elaboration,
obsX + passoE * cos(dirPto - M_PI / 2),
obsY + passoN * sin(dirPto - M_PI / 2),
- parBilin);
+ parBilin, gradient);
g[5][0] = gradient[0];
g[5][1] = gradient[1];
- gradient =
- Get_Gradient(Elaboration,
+ Get_Gradient(Elaboration,
obsX + passoE * cos(dirPto + M_PI * 3 / 4),
obsY + passoN * sin(dirPto + M_PI * 3 / 4),
- parBilin);
+ parBilin, gradient);
g[6][0] = gradient[0];
g[6][1] = gradient[1];
- gradient =
- Get_Gradient(Elaboration,
+ Get_Gradient(Elaboration,
obsX + passoE * cos(dirPto - M_PI * 3 / 4),
obsY + passoN * sin(dirPto - M_PI * 3 / 4),
- parBilin);
+ parBilin, gradient);
g[8][0] = gradient[0];
g[8][1] = gradient[1];
c2 = 0;
for (c1 = 0; c1 < 9; c1++)
- if (sqrt(g[c1][0] * g[c1][0] + g[c1][1] * g[c1][1]) >
+ if (g[c1][0] * g[c1][0] + g[c1][1] * g[c1][1] >
gradHigh)
c2++;
@@ -137,23 +131,21 @@
return PRE_TERRAIN;
}
else
- return UNKNOWN;
+ return PRE_UNKNOWN;
} /* END ELSE IF */
else
return PRE_TERRAIN;
}
-double *Get_Gradient(struct Cell_head Elaboration, double X, double Y,
- double *parVect)
+int Get_Gradient(struct Cell_head Elaboration, double X, double Y,
+ double *parVect, double *grad)
{
int row, col, N;
- double csi, eta, d, b, a, c, *grad;
+ double csi, eta, d, b, a, c;
extern int nsply;
extern double passoN, passoE;
- grad = (double *)G_calloc(2, sizeof(double));
-
row = (int)((Y - Elaboration.south) / passoN);
col = (int)((X - Elaboration.west) / passoE);
N = nsply * col + row;
@@ -165,7 +157,7 @@
c = parVect[N + 1 + nsply] - a - b - d;
grad[0] = (a + c * eta);
grad[1] = (b + c * csi);
- return grad;
+ return 0;
}
void classification(struct Map_info *Out, struct Cell_head Elaboration,
@@ -173,10 +165,10 @@
double *parBilin, double *parBicub, double mean,
double alpha, double gradHigh, double gradLow,
double overlap, int *line_num, int num_points,
- dbDriver * driver, char *vect_name)
+ dbDriver * driver, char *tabint_name, char *tab_name)
{
int i, edge;
- double interpolation, weight, residual, eta, csi, *gradient;
+ double interpolation, weight, residual, eta, csi, gradient[2];
extern int nsplx, nsply, line_out_counter;
extern double passoN, passoE;
@@ -187,8 +179,10 @@
point = Vect_new_line_struct();
categories = Vect_new_cats_struct();
+ db_begin_transaction(driver);
+
for (i = 0; i < num_points; i++) { /* Sparse points */
- G_percent(i, num_points, 6);
+ G_percent(i, num_points, 2);
Vect_reset_line(point);
Vect_reset_cats(categories);
@@ -200,15 +194,14 @@
Elaboration.south, parBicub);
interpolation += mean;
- Vect_copy_xyz_to_pnts(point, &obs[i][0], &obs[i][1], &obs[i][2],
- 1);
- gradient =
- Get_Gradient(Elaboration, obs[i][0], obs[i][1], parBilin);
+ Vect_copy_xyz_to_pnts(point, &obs[i][0], &obs[i][1], &obs[i][2], 1);
+
+ Get_Gradient(Elaboration, obs[i][0], obs[i][1], parBilin, gradient);
+
*point->z += mean;
/*Vect_cat_set (categories, F_INTERPOLATION, line_out_counter); */
if (Vect_point_in_box(obs[i][0], obs[i][1], interpolation, &Overlap)) { /*(5) */
-
residual = *point->z - interpolation;
edge =
edge_detection(Elaboration, Overlap, parBilin, *point->x,
@@ -219,71 +212,70 @@
Vect_cat_set(categories, F_INTERPOLATION, line_out_counter);
Vect_write_line(Out, GV_POINT, point, categories);
Insert_Interpolation(interpolation, line_out_counter, driver,
- vect_name);
+ tabint_name);
line_out_counter++;
-
}
else {
- if ((*point->x > Overlap.E) && (*point->x != General.E)) {
+ if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = csi * eta;
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(3) */
- csi = (*point->x - Overlap.E) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - csi) * (1 - eta);
-
gradient[0] *= weight;
gradient[1] *= weight;
interpolation *= weight;
- if (Select
- (&gradient[0], &gradient[1], &interpolation,
- line_num[i], driver) != DB_OK)
- G_fatal_error(_("Impossible to read the database"));
+ if (Select(&gradient[0], &gradient[1],
+ &interpolation, line_num[i], driver,
+ tab_name) != DB_OK)
+ G_fatal_error(_("Impossible to read from aux table"));
- if (UpDate
- (gradient[0], gradient[1], interpolation,
- line_num[i], driver) != DB_OK)
- G_fatal_error(_("Impossible to update the database"));
-
+ if (UpDate(gradient[0], gradient[1],
+ interpolation, line_num[i], driver,
+ tab_name) != DB_OK)
+ G_fatal_error(_("Impossible to update aux table"));
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(1) */
- csi = (*point->x - Overlap.E) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ csi = (General.E - *point->x) / overlap;
eta = (*point->y - General.S) / overlap;
- weight = (1 - csi) * eta;
+ weight = csi * eta;
- if (Insert
- (gradient[0] * weight, gradient[1] * weight,
- interpolation * weight, line_num[i],
- driver) != DB_OK)
- G_fatal_error(_("Impossible to write in the database"));
+ gradient[0] *= weight;
+ gradient[1] *= weight;
+ interpolation *= weight;
+ if (Insert(gradient[0], gradient[1], interpolation,
+ line_num[i], driver, tab_name) != DB_OK)
+ G_fatal_error(_("Impossible to write to aux table"));
+
}
else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(1) */
- weight = (*point->x - Overlap.E) / overlap;
+ weight = (General.E - *point->x) / overlap;
- if (Insert
- (gradient[0] * weight, gradient[1] * weight,
- interpolation * weight, line_num[i],
- driver) != DB_OK)
- G_fatal_error(_("Impossible to write in the database"));
- }
+ gradient[0] *= weight;
+ gradient[1] *= weight;
+ interpolation *= weight;
+ if (Insert(gradient[0], gradient[1], interpolation,
+ line_num[i], driver, tab_name) != DB_OK)
+ G_fatal_error(_("Impossible to write to aux table"));
+ }
}
- else if ((*point->x < Overlap.W) && (*point->x != General.W)) {
-
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(4) */
+ else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(4) */
csi = (*point->x - General.W) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - eta) * csi;
+ eta = (General.N - *point->y) / overlap;
+ weight = eta * csi;
gradient[0] *= weight;
gradient[1] *= weight;
interpolation *= weight;
- if (Select
- (&gradient[0], &gradient[1], &interpolation,
- line_num[i], driver) != DB_OK)
- G_fatal_error(_("Impossible to read the database"));
+ if (Select(&gradient[0], &gradient[1],
+ &interpolation, line_num[i], driver,
+ tab_name) != DB_OK)
+ G_fatal_error(_("Impossible to read from aux table"));
residual = *point->z - interpolation;
edge =
@@ -298,11 +290,10 @@
line_out_counter);
Vect_write_line(Out, GV_POINT, point, categories);
Insert_Interpolation(interpolation, line_out_counter,
- driver, vect_name);
+ driver, tabint_name);
line_out_counter++;
-
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(2) */
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(2) */
csi = (*point->x - General.W) / overlap;
eta = (*point->y - General.S) / overlap;
weight = csi * eta;
@@ -311,28 +302,28 @@
gradient[1] *= weight;
interpolation *= weight;
- if (Select
- (&gradient[0], &gradient[1], &interpolation,
- line_num[i], driver) != DB_OK)
- G_fatal_error(_("Impossible to read the database"));
+ if (Select(&gradient[0], &gradient[1],
+ &interpolation, line_num[i], driver,
+ tab_name) != DB_OK)
+ G_fatal_error(_("Impossible to read from aux table"));
- if (UpDate
- (gradient[0], gradient[1], interpolation,
- line_num[i], driver) != DB_OK)
- G_fatal_error(_("Impossible to update the database"));
+ if (UpDate(gradient[0], gradient[1],
+ interpolation, line_num[i], driver,
+ tab_name) != DB_OK)
+ G_fatal_error(_("Impossible to update aux table"));
}
else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(2) */
- weight = (Overlap.W - *point->x) / overlap;
+ weight = (*point->x - General.W) / overlap;
gradient[0] *= weight;
gradient[1] *= weight;
interpolation *= weight;
- if (Select
- (&gradient[0], &gradient[1], &interpolation,
- line_num[i], driver) != DB_OK)
- G_fatal_error(_("Impossible to read the database"));
+ if (Select(&gradient[0], &gradient[1],
+ &interpolation, line_num[i], driver,
+ tab_name) != DB_OK)
+ G_fatal_error(_("Impossible to read from aux table"));
residual = *point->z - interpolation;
edge =
@@ -347,23 +338,22 @@
line_out_counter);
Vect_write_line(Out, GV_POINT, point, categories);
Insert_Interpolation(interpolation, line_out_counter,
- driver, vect_name);
+ driver, tabint_name);
line_out_counter++;
}
-
}
else if ((*point->x <= Overlap.E) && (*point->x >= Overlap.W)) {
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(3) */
- weight = (*point->y - Overlap.N) / overlap;
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ weight = (General.N - *point->y) / overlap;
gradient[0] *= weight;
gradient[1] *= weight;
interpolation *= weight;
- if (Select
- (&gradient[0], &gradient[1], &interpolation,
- line_num[i], driver) != DB_OK)
- G_fatal_error(_("Impossible to read the database"));
+ if (Select(&gradient[0], &gradient[1],
+ &interpolation, line_num[i], driver,
+ tab_name) != DB_OK)
+ G_fatal_error(_("Impossible to read from aux table"));
residual = *point->z - interpolation;
edge =
@@ -378,84 +368,96 @@
line_out_counter);
Vect_write_line(Out, GV_POINT, point, categories);
Insert_Interpolation(interpolation, line_out_counter,
- driver, vect_name);
+ driver, tabint_name);
line_out_counter++;
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(1) */
- weight = (Overlap.S - *point->y) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ weight = (*point->y - General.S) / overlap;
- if (Insert
- (gradient[0] * weight, gradient[1] * weight,
- interpolation * weight, line_num[i], driver)
- != DB_OK)
- G_fatal_error(_("Impossible to write in the database"));
+ gradient[0] *= weight;
+ gradient[1] *= weight;
+ interpolation *= weight;
+ if (Insert(gradient[0], gradient[1], interpolation,
+ line_num[i], driver, tab_name) != DB_OK)
+ G_fatal_error(_("Impossible to write to aux table"));
+
} /*else (1) */
} /*else */
}
} /*end if obs */
} /*end for */
- G_percent(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)
+int Insert(double partialX, double partialY, double Interp,
+ int line_num, dbDriver * driver, char *tab_name)
{
-
char buf[1024];
dbString sql;
+ int ret;
db_init_string(&sql);
sprintf(buf,
- "INSERT INTO Auxiliar_edge_table (ID, Interp, partialX, partialY)");
+ "INSERT INTO %s (ID, Interp, X, Y)", tab_name);
db_append_string(&sql, buf);
- sprintf(buf, " VALUES (%d, %lf, %lf, %lf)", line_num, Interp, partialX,
- partialY);
+ sprintf(buf, " VALUES (%d, %lf, %lf, %lf)", line_num, Interp,
+ partialX, partialY);
db_append_string(&sql, buf);
- return db_execute_immediate(driver, &sql);
+ ret = db_execute_immediate(driver, &sql);
+ db_free_string(&sql);
+
+ return ret;
}
int Insert_Interpolation(double Interp, int line_num, dbDriver * driver,
- char *name)
+ char *tab_name)
{
-
char buf[1024];
dbString sql;
+ int ret;
db_init_string(&sql);
- sprintf(buf, "INSERT INTO %s_edge_Interpolation (ID, Interp)", name);
+ sprintf(buf, "INSERT INTO %s (ID, Interp)", tab_name);
db_append_string(&sql, buf);
sprintf(buf, " VALUES (%d, %lf)", line_num, Interp);
db_append_string(&sql, buf);
- return db_execute_immediate(driver, &sql);
+ ret = db_execute_immediate(driver, &sql);
+ db_free_string(&sql);
+
+ return ret;
}
-int UpDate(double partialX, double partialY, double Interp, int line_num,
- dbDriver * driver)
+int UpDate(double partialX, double partialY, double Interp,
+ int line_num, dbDriver * driver, char *tab_name)
{
-
char buf[1024];
dbString sql;
+ int ret;
db_init_string(&sql);
sprintf(buf,
- "UPDATE Auxiliar_edge_table SET Interp=%lf, PartialX=%lf, PartialY=%lf WHERE ID=%d",
- Interp, partialX, partialY, line_num);
+ "UPDATE %s SET Interp=%lf, X=%lf, Y=%lf WHERE ID=%d",
+ tab_name, Interp, partialX, partialY, line_num);
db_append_string(&sql, buf);
- return db_execute_immediate(driver, &sql);
+ ret = db_execute_immediate(driver, &sql);
+ db_free_string(&sql);
+
+ return ret;
}
-int Select(double *PartialX, double *PartialY, double *Interp, int line_num,
- dbDriver * driver)
+int Select(double *PartialX, double *PartialY, double *Interp,
+ int line_num, dbDriver * driver, char *tab_name)
{
-
int more;
char buf[1024];
dbString sql;
@@ -466,15 +468,16 @@
db_init_string(&sql);
sprintf(buf,
- "SELECT ID, Interp, partialX, partialY FROM Auxiliar_edge_table WHERE ID=%d",
- line_num);
+ "SELECT ID, Interp, X, Y FROM %s WHERE ID=%d",
+ tab_name, line_num);
db_append_string(&sql, buf);
if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
return -1;
+ table = db_get_cursor_table(&cursor);
+
while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
- table = db_get_cursor_table(&cursor);
Interp_col = db_get_table_column(table, 1);
PartialX_col = db_get_table_column(table, 2);
@@ -503,107 +506,34 @@
*PartialY += db_get_value_double(PartialY_value);
}
db_close_cursor(&cursor);
+ db_free_string(&sql);
return DB_OK;
}
-
-int Create_AuxEdge_Table(dbDriver * driver)
-{
- dbTable *table;
- dbColumn *ID_col, *Interp_col, *PartialX_col, *PartialY_col;
- int created;
-
- table = db_alloc_table(4);
- db_set_table_name(table, "Auxiliar_edge_table");
- db_set_table_description(table,
- "It is used for the intermediate interpolated and gradient values");
-
- ID_col = db_get_table_column(table, 0);
- db_set_column_name(ID_col, "ID");
- db_set_column_sqltype(ID_col, DB_SQL_TYPE_INTEGER);
-
- Interp_col = db_get_table_column(table, 1);
- db_set_column_name(Interp_col, "Interp");
- db_set_column_sqltype(Interp_col, DB_SQL_TYPE_REAL);
-
- PartialX_col = db_get_table_column(table, 2);
- db_set_column_name(PartialX_col, "PartialX");
- db_set_column_sqltype(PartialX_col, DB_SQL_TYPE_REAL);
-
- PartialY_col = db_get_table_column(table, 3);
- db_set_column_name(PartialY_col, "PartialY");
- db_set_column_sqltype(PartialY_col, DB_SQL_TYPE_REAL);
-
- if (db_create_table(driver, table) == DB_OK) {
- G_debug(3, _("<Auxiliar_edge_table> created in database."));
- created = TRUE;
- }
- else
- return FALSE;
-
- return created;
-}
-
-int Drop_Aux_Table(dbDriver * driver)
-{
- dbString drop;
-
- db_init_string(&drop);
- db_append_string(&drop, "drop table ");
- db_append_string(&drop, "Auxiliar_edge_table");
- return db_execute_immediate(driver, &drop);
-
-}
-
-
-int Create_Interpolation_Table(char *vect_name, dbDriver * driver)
-{
- char table_name[1024];
-
- dbTable *table;
- dbColumn *ID_col, *Interp_col;
-
- sprintf(table_name, "%s_edge_Interpolation", vect_name);
-
- table = db_alloc_table(2);
- db_set_table_name(table, table_name);
- db_set_table_description(table,
- "This table is the bicubic interpolation of the input vector");
-
- ID_col = db_get_table_column(table, 0);
- db_set_column_name(ID_col, "ID");
- db_set_column_sqltype(ID_col, DB_SQL_TYPE_INTEGER);
-
- Interp_col = db_get_table_column(table, 1);
- db_set_column_name(Interp_col, "Interp");
- db_set_column_sqltype(Interp_col, DB_SQL_TYPE_REAL);
-
- if (db_create_table(driver, table) == DB_OK) {
- G_debug(3, _("<%s> created in database."), db_get_table_name(table));
- return DB_OK;
- }
- else
- return !DB_OK;
-}
-
-
-#ifdef notdef
/*! DEFINITION OF THE SUBZONES
- -----------------------
- |4| 3 |3| | |
- -----------------------
- | | | | | |
- |2| 5 |1| | |
- | | | | | |
- -----------------------
- |2| 1 |1| | |
- -----------------------
- | | | | | |
- | | | | | |
- | | | | | |
- -----------------------
- | | | | | |
- -----------------------
+ 5: inside Overlap region
+ all others: inside General region but outside Overlap region
+
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | |4| 3 |3| | |
+ ---------------------------------
+ | | | | | | | |
+ | | |2| 5 |1| | |
+ | | | | | | | |
+ ---------------------------------
+ | | |2| 1 |1| | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
*/
-#endif
Modified: grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/edgedetection.h
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/edgedetection.h 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/edgedetection.h 2010-02-17 17:56:52 UTC (rev 41062)
@@ -38,8 +38,8 @@
double *, /**/
double, /**/ double, /**/ double, /**/ double /**/);
-double *Get_Gradient(struct Cell_head, /**/
- double, /**/ double, /**/ double * /**/);
+int Get_Gradient(struct Cell_head, /**/
+ double, /**/ double, /**/ double *, /**/ double *);
void classification(struct Map_info *, /**/
struct Cell_head, /**/
@@ -53,16 +53,12 @@
double, /**/
double, /**/
double, /**/
- int *, /**/ int, /**/ dbDriver *, /**/ char * /**/);
+ int *, /**/ int, /**/ dbDriver *, /**/ char *, /**/ char *);
-int Insert(double, /**/ double, /**/ double, /**/ int, /**/ dbDriver * /**/);
+int Insert(double, /**/ double, /**/ double, /**/ int, /**/ dbDriver *, /**/ char *);
-int UpDate(double, /**/ double, /**/ double, /**/ int, /**/ dbDriver * /**/);
+int UpDate(double, /**/ double, /**/ double, /**/ int, /**/ dbDriver *, /**/ char *);
-int Select(double *, /**/
- double *, /**/ double *, /**/ int, /**/ dbDriver * /**/);
+int Select(double *, /**/ double *, /**/ double *, /**/ int, /**/ dbDriver *, /**/ char *);
-int Create_AuxEdge_Table(dbDriver *);
-int Create_Interpolation_Table(char *, dbDriver *);
int Insert_Interpolation(double, int, dbDriver *, char *);
-int Drop_Aux_Table(dbDriver *);
Modified: grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/main.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/main.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.lidar.edgedetection/main.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -25,6 +25,7 @@
/*INCLUDES*/
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include <grass/gis.h>
#include <grass/Vect.h>
#include <grass/dbmi.h>
@@ -32,21 +33,25 @@
#include <grass/config.h>
#include <grass/PolimiFunct.h>
#include "edgedetection.h"
-int nsply, nsplx, line_out_counter, first_it;
+
+int nsply, nsplx, line_out_counter;
double passoN, passoE;
-
/**************************************************************************************
**************************************************************************************/
int main(int argc, char *argv[])
{
/* Variables' declarations */
+ int nsplx_adj, nsply_adj;
+ int nsubregion_col, nsubregion_row, subregion = 0, nsubregions = 0;
+ double N_extension, E_extension, orloE, orloN;
int dim_vect, nparameters, BW, npoints;
- double lambda_B, lambda_F, grad_H, grad_L, alpha, ew_resol, ns_resol,
- mean;
- char *dvr, *db, *mapset, table_interpolation[1024], table_name[1024];
+ double lambda_B, lambda_F, grad_H, grad_L, alpha, mean;
+ const char *dvr, *db, *mapset;
+ char table_interpolation[GNAME_MAX], table_name[GNAME_MAX];
+ char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
- int last_row, last_column, flag_auxiliar = FALSE, pass_num;
+ int last_row, last_column, flag_auxiliar = FALSE;
int *lineVect;
double *TN, *Q, *parVect_bilin, *parVect_bicub; /* Interpolating and least-square vectors */
@@ -56,6 +61,7 @@
struct Map_info In, Out;
struct Option *in_opt, *out_opt, *passoE_opt, *passoN_opt,
*lambdaF_opt, *lambdaB_opt, *gradH_opt, *gradL_opt, *alfa_opt;
+ struct Flag *spline_step_flag;
struct GModule *module;
struct Cell_head elaboration_reg, original_reg;
@@ -73,6 +79,12 @@
module->description =
_("Detects the object's edges from a LIDAR data set.");
+ spline_step_flag = G_define_flag();
+ spline_step_flag->key = 'e';
+ spline_step_flag->label = _("Estimate point density and distance");
+ spline_step_flag->description =
+ _("Estimate point density and distance for the input vector points within the current region extends and quit");
+
in_opt = G_define_standard_option(G_OPT_V_INPUT);
out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
@@ -153,6 +165,10 @@
grad_H = atof(gradH_opt->answer);
grad_L = atof(gradL_opt->answer);
alpha = atof(alfa_opt->answer);
+
+ grad_L = grad_L * grad_L;
+ grad_H = grad_H * grad_H;
+
if (!(db = G__getenv2("DB_DATABASE", G_VAR_MAPSET)))
G_fatal_error(_("Unable to read name of database"));
@@ -160,8 +176,39 @@
G_fatal_error(_("Unable to read name of driver"));
/* Setting auxiliar table's name */
- sprintf(table_name, "%s_aux", out_opt->answer);
+ if (G__name_is_fully_qualified(out_opt->answer, xname, xmapset)) {
+ sprintf(table_name, "%s_aux", xname);
+ sprintf(table_interpolation, "%s_edge_Interpolation", xname);
+ }
+ else {
+ sprintf(table_name, "%s_aux", out_opt->answer);
+ sprintf(table_interpolation, "%s_edge_Interpolation", out_opt->answer);
+ }
+ /* Something went wrong in a previous v.lidar.edgedetection execution */
+ if (db_table_exists(dvr, db, table_name)) {
+ /* Start driver and open db */
+ driver = db_start_driver_open_database(dvr, db);
+ if (driver == NULL)
+ G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
+ dvr);
+ if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
+ G_fatal_error(_("Old auxiliar table could not be dropped"));
+ db_close_database_shutdown_driver(driver);
+ }
+
+ /* Something went wrong in a previous v.lidar.edgedetection execution */
+ if (db_table_exists(dvr, db, table_interpolation)) {
+ /* Start driver and open db */
+ driver = db_start_driver_open_database(dvr, db);
+ if (driver == NULL)
+ G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
+ dvr);
+ if (P_Drop_Aux_Table(driver, table_interpolation) != DB_OK)
+ G_fatal_error(_("Old auxiliar table could not be dropped"));
+ db_close_database_shutdown_driver(driver);
+ }
+
/* Checking vector names */
Vect_check_input_output_name(in_opt->answer, out_opt->answer,
GV_FATAL_EXIT);
@@ -170,15 +217,33 @@
G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
}
- /* Open output vector */
- if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z))
- G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
-
Vect_set_open_level(1);
/* Open input vector */
if (1 > Vect_open_old(&In, in_opt->answer, mapset))
G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
+ /* Input vector must be 3D */
+ if (!Vect_is_3d(&In))
+ G_fatal_error(_("Input vector map <%s> is not 3D!"), in_opt->answer);
+
+ /* Estimate point density and mean distance for current region */
+ if (spline_step_flag->answer) {
+ double dens, dist;
+ if (P_estimate_splinestep(&In, &dens, &dist) == 0) {
+ G_message("Estimated point density: %.4g", dens);
+ G_message("Estimated mean distance between points: %.4g", dist);
+ }
+ else
+ G_warning(_("No points in current region!"));
+
+ Vect_close(&In);
+ exit(EXIT_SUCCESS);
+ }
+
+ /* Open output vector */
+ if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z))
+ G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
+
/* Copy vector Head File */
Vect_copy_head_data(&In, &Out);
Vect_hist_copy(&In, &Out);
@@ -190,36 +255,68 @@
G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
dvr);
- if (Create_Interpolation_Table(out_opt->answer, driver) != DB_OK)
+ /* Create auxiliar and interpolation table */
+ if ((flag_auxiliar = P_Create_Aux4_Table(driver, table_name)) == FALSE)
+ G_fatal_error(_("It was impossible to create <%s>."), table_name);
+
+ if (P_Create_Aux2_Table(driver, table_interpolation) == FALSE)
G_fatal_error(_("It was impossible to create <%s> interpolation table in database."),
out_opt->answer);
+ db_create_index2(driver, table_name, "ID");
+ db_create_index2(driver, table_interpolation, "ID");
+ /* sqlite likes that */
+ db_close_database_shutdown_driver(driver);
+ driver = db_start_driver_open_database(dvr, db);
+
/* Setting regions and boxes */
G_get_set_window(&original_reg);
G_get_set_window(&elaboration_reg);
Vect_region_box(&elaboration_reg, &overlap_box);
Vect_region_box(&elaboration_reg, &general_box);
+ /*------------------------------------------------------------------
+ | 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_PASS times the east-west resolution. */
-
- ew_resol = original_reg.ew_res;
- ns_resol = original_reg.ns_res;
-
P_zero_dim(&dims);
- dims.latoE = NSPLX_MAX * passoE;
- dims.latoN = NSPLY_MAX * passoN;
- dims.overlap = OVERLAP_SIZE * ew_resol;
+ nsplx_adj = NSPLX_MAX;
+ nsply_adj = NSPLY_MAX;
+ if (passoN > passoE)
+ dims.overlap = OVERLAP_SIZE * passoN;
+ else
+ dims.overlap = OVERLAP_SIZE * passoE;
P_get_orlo(P_BICUBIC, &dims, passoE, passoN);
+ P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
+ G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
+ G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
+
+ /* calculate number of subregions */
+ orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
+ orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+
+ N_extension = original_reg.north - original_reg.south;
+ E_extension = original_reg.east - original_reg.west;
+
+ nsubregion_col = ceil(E_extension / orloE) + 0.5;
+ nsubregion_row = ceil(N_extension / orloN) + 0.5;
+
+ if (nsubregion_col < 0)
+ nsubregion_col = 0;
+ if (nsubregion_row < 0)
+ nsubregion_row = 0;
+
+ nsubregions = nsubregion_row * nsubregion_col;
+
elaboration_reg.south = original_reg.north;
-
last_row = FALSE;
- first_it = TRUE;
- pass_num = 0;
while (last_row == FALSE) { /* For each row */
@@ -239,10 +336,12 @@
nsply =
ceil((elaboration_reg.north - elaboration_reg.south) / passoN) +
- 1;
+ 0.5;
+ /*
if (nsply > NSPLY_MAX) {
nsply = NSPLY_MAX;
}
+ */
G_debug(1, "nsply = %d", nsply);
elaboration_reg.east = original_reg.west;
@@ -250,6 +349,10 @@
while (last_column == FALSE) { /* For each column */
+ subregion++;
+ if (nsubregions > 1)
+ G_message(_("subregion %d of %d"), subregion, nsubregions);
+
P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
GENERAL_COLUMN);
@@ -266,10 +369,12 @@
nsplx =
ceil((elaboration_reg.east - elaboration_reg.west) / passoE) +
- 1;
+ 0.5;
+ /*
if (nsplx > NSPLX_MAX) {
nsplx = NSPLX_MAX;
}
+ */
G_debug(1, "nsplx = %d", nsplx);
/*Setting the active region */
@@ -282,9 +387,6 @@
if (npoints > 0) { /* If there is any point falling into elaboration_reg... */
int i, tn;
- pass_num++;
- G_verbose_message(_("Pass %d:"), pass_num);
-
nparameters = nsplx * nsply;
/* Mean's calculation */
@@ -295,7 +397,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 */
@@ -310,7 +412,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,
@@ -341,36 +443,30 @@
G_free_vector(TN);
G_free_vector(Q);
- if (flag_auxiliar == FALSE) {
- G_debug(1,
- _("Creating auxiliar table for archiving overlapping zones"));
- if ((flag_auxiliar =
- Create_AuxEdge_Table(driver)) == FALSE)
- G_fatal_error(_("It was impossible to create <%s>."),
- table_name);
- }
G_verbose_message(_("Point classification"));
classification(&Out, elaboration_reg, general_box,
overlap_box, obsVect, parVect_bilin,
parVect_bicub, mean, alpha, grad_H, grad_L,
dims.overlap, lineVect, npoints, driver,
- out_opt->answer);
+ table_interpolation, table_name);
G_free_vector(parVect_bilin);
G_free_vector(parVect_bicub);
G_free_matrix(obsVect);
G_free_ivector(lineVect);
} /* IF */
-
- G_free(observ);
- first_it = FALSE;
+ else {
+ G_free(observ);
+ G_warning(_("No data within this subregion. "
+ "Consider changing the spline step."));
+ }
} /*! END WHILE; last_column = TRUE */
} /*! END WHILE; last_row = TRUE */
/* Dropping auxiliar table */
if (npoints > 0) {
G_debug(1, _("Dropping <%s>"), table_name);
- if (Drop_Aux_Table(driver) != DB_OK)
+ if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
G_warning(_("Auxiliar table could not be dropped"));
}
@@ -378,7 +474,6 @@
Vect_close(&In);
- sprintf(table_interpolation, "%s_edge_Interpolation", out_opt->answer);
Vect_map_add_dblink(&Out, F_INTERPOLATION, NULL, table_interpolation,
"id", db, dvr);
Modified: grass/branches/develbranch_6/vector/lidar/v.lidar.growing/main.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.lidar.growing/main.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.lidar.growing/main.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -22,6 +22,7 @@
/*INCLUDES*/
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include <grass/config.h>
#include <grass/gis.h>
#include <grass/Vect.h>
@@ -38,6 +39,8 @@
/* Variables' declarations */
int row, nrows, col, ncols, MaxPoints;
+ int nsubregion_col, nsubregion_row;
+ int subregion = 0, nsubregions = 0;
int last_row, last_column;
int nlines, nlines_first, line_num;
int more;
@@ -45,7 +48,8 @@
double Z_interp;
double Thres_j, Thres_d, ew_resol, ns_resol;
double minNS, minEW, maxNS, maxEW;
- char *mapset, buf[1024];
+ const char *mapset;
+ char buf[1024], table_name[GNAME_MAX];
char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
int colorBordo, ripieno, conta, lungPunti, lungHull, xi, c1, c2;
@@ -68,8 +72,6 @@
dbString sql;
dbTable *table;
dbCursor cursor;
- const char *p;
- int found;
/*------------------------------------------------------------------------------------------*/
/* Options' declaration */ ;
@@ -80,6 +82,8 @@
"algorithm for determining the building inside");
in_opt = G_define_standard_option(G_OPT_V_INPUT);
+ in_opt->description =
+ _("Input vector (v.lidar.edgedetection output");
out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
@@ -121,27 +125,22 @@
Vect_check_input_output_name(in_opt->answer, out_opt->answer,
GV_FATAL_EXIT);
- sprintf(xname, "%s", in_opt->answer);
- found = 0;
- for (p = in_opt->answer; *p; p++)
- if (*p == '@') {
- found = 1;
- break;
- }
-
- if (found) {
- if (G__name_is_fully_qualified(in_opt->answer, xname, xmapset) < 0) /* strip off mapset from input for SQL */
- G_fatal_error(_("Vector map <%s> not found"), xname);
- }
if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL) {
G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
}
- /*Vect_set_open_level (2); WITH TOPOLOGY */
+ /* Setting auxiliar table's name */
+ if (G__name_is_fully_qualified(in_opt->answer, xname, xmapset)) {
+ sprintf(table_name, "%s_edge_Interpolation", xname);
+ }
+ else
+ sprintf(table_name, "%s_edge_Interpolation", in_opt->answer);
+
Vect_set_open_level(1); /* WITHOUT TOPOLOGY */
if (Vect_open_old(&In, in_opt->answer, mapset) < 1)
G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
+ Vect_set_open_level(1); /* WITHOUT TOPOLOGY */
if (Vect_open_old(&First, first_opt->answer, mapset) < 1)
G_fatal_error(_("Unable to open vector map <%s>"), first_opt->answer);
@@ -167,40 +166,65 @@
G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
field->driver);
- /* Setting regions and boxes */
- G_debug(1, _("Setting regions and boxes"));
- G_get_set_window(&original_reg);
- G_get_set_window(&elaboration_reg);
-
- /* Fixxing parameters of the elaboration region */
- /*! The original_region will be divided into several subregions */
- ew_resol = original_reg.ew_res;
- ns_resol = original_reg.ns_res;
-
- /* Subdividing and working with tiles */
- elaboration_reg.south = original_reg.north;
- last_row = FALSE;
-
+ /* is this the right place to open the cursor ??? */
+
db_init_string(&sql);
db_zero_string(&sql);
- sprintf(buf, "SELECT Interp,ID FROM %s_edge_Interpolation", xname);
+ sprintf(buf, "SELECT Interp,ID FROM %s", table_name);
G_debug(1, "buf: %s", buf);
db_append_string(&sql, buf);
if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
- G_fatal_error(_("Unable to open table <%s_edge_Interpolation>"), xname);
+ G_fatal_error(_("Unable to open table <%s>"), table_name);
count_obj = 1;
- nlines = Vect_get_num_lines(&In);
+ /* no topology, get number of lines in input vector */
+ nlines = 0;
points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
+ Vect_rewind(&In);
+ while (Vect_read_next_line(&In, points, Cats) > 0) {
+ nlines++;
+ }
+ Vect_rewind(&In);
- nlines_first = Vect_get_num_lines(&First);
+ /* no topology, get number of lines in first pulse input vector */
+ nlines_first = 0;
points_first = Vect_new_line_struct();
Cats_first = Vect_new_cats_struct();
+ Vect_rewind(&First);
+ while (Vect_read_next_line(&First, points_first, Cats_first) > 0) {
+ nlines_first++;
+ }
+ Vect_rewind(&First);
+ /* Setting regions and boxes */
+ G_debug(1, _("Setting regions and boxes"));
+ G_get_set_window(&original_reg);
+ G_get_set_window(&elaboration_reg);
+
+ /* Fixing parameters of the elaboration region */
+ /*! The original_region will be divided into subregions */
+ ew_resol = original_reg.ew_res;
+ ns_resol = original_reg.ns_res;
+
+ /* calculate number of subregions */
+ nsubregion_col = ceil((original_reg.east - original_reg.west) / (LATO * ew_resol)) + 0.5;
+ nsubregion_row = ceil((original_reg.north - original_reg.south) / (LATO * ns_resol)) + 0.5;
+
+ if (nsubregion_col < 0)
+ nsubregion_col = 0;
+ if (nsubregion_row < 0)
+ nsubregion_row = 0;
+
+ nsubregions = nsubregion_row * nsubregion_col;
+
+ /* Subdividing and working with tiles */
+ elaboration_reg.south = original_reg.north;
+ last_row = FALSE;
+
while (last_row == FALSE) { /* For each strip of LATO rows */
elaboration_reg.north = elaboration_reg.south;
@@ -220,6 +244,10 @@
while (last_column == FALSE) { /* For each strip of LATO columns */
BOUND_BOX elaboration_box;
+ subregion++;
+ if (nsubregions > 1)
+ G_message(_("subregion %d of %d"), subregion, nsubregions);
+
elaboration_reg.west = elaboration_reg.east;
if (elaboration_reg.west < original_reg.west) /* First column */
elaboration_reg.west = original_reg.west;
@@ -231,10 +259,13 @@
last_column = TRUE;
}
- /*Setting the active region */
- G_set_window(&elaboration_reg);
- nrows = G_window_rows();
- ncols = G_window_cols();
+ /* Setting the active region */
+ elaboration_reg.ns_res = ns_resol;
+ elaboration_reg.ew_res = ew_resol;
+ nrows = (elaboration_reg.north - elaboration_reg.south) / ns_resol + 0.1;
+ ncols = (elaboration_reg.east - elaboration_reg.west) / ew_resol + 0.1;
+ elaboration_reg.rows = nrows;
+ elaboration_reg.cols = ncols;
G_debug(1, _("Rows = %d"), nrows);
G_debug(1, _("Columns = %d"), ncols);
@@ -257,13 +288,13 @@
}
}
+ G_verbose_message(_("read points in input vector"));
+ Vect_region_box(&elaboration_reg, &elaboration_box);
line_num = 0;
Vect_rewind(&In);
while (Vect_read_next_line(&In, points, Cats) > 0) {
line_num++;
- Vect_region_box(&elaboration_reg, &elaboration_box);
-
if ((Vect_point_in_box
(points->x[0], points->y[0], points->z[0],
&elaboration_box)) &&
@@ -279,6 +310,37 @@
(int)(G_easting_to_col
(points->x[0], &elaboration_reg));
+ Z_interp = 0;
+ /* TODO: make sure the current db_fetch() usage works */
+ /* why not: */
+ /*
+ db_init_string(&sql);
+ sprintf(buf, "SELECT Interp,ID FROM %s WHERE ID=%d", table_name, line_num);
+ db_append_string(&sql, buf);
+
+ if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
+ G_fatal_error(_("Unable to open table <%s>"), table_name);
+
+ while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
+ dbColumn *Z_Interp_col;
+ dbValue *Z_Interp_value;
+ table = db_get_cursor_table(&cursor);
+
+ Z_Interp_col = db_get_table_column(table, 1);
+
+ if (db_sqltype_to_Ctype(db_get_column_sqltype(Z_Interp_col)) ==
+ DB_C_TYPE_DOUBLE)
+ Z_Interp_value = db_get_column_value(Z_Interp_col);
+ else
+ continue;
+
+ Z_interp = db_get_value_double(Z_Interp_value);
+ break;
+ }
+ db_close_cursor(&cursor);
+ db_free_string(&sql);
+ */
+ /* instead of */
while (1) {
if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK ||
!more)
@@ -353,13 +415,9 @@
(points_first->x[0], points_first->y[0],
points_first->z[0], &elaboration_box)) &&
((points->x[0] != elaboration_reg.west) ||
- (points->x[0] == original_reg.west)) && ((points->y[0]
- !=
- elaboration_reg.
- north) ||
- (points->y[0] ==
- original_reg.
- north))) {
+ (points->x[0] == original_reg.west)) &&
+ ((points->y[0] != elaboration_reg.north) ||
+ (points->y[0] == original_reg.north))) {
row =
(int)(G_northing_to_row
@@ -379,7 +437,7 @@
/* REGION GROWING */
if (region == TRUE) {
- G_debug(1, "Region Growing");
+ G_verbose_message(_("Region Growing"));
punti_bordo = G_alloc_matrix(MaxPoints, 3);
P = Pvector(0, MaxPoints);
@@ -388,6 +446,7 @@
ripieno = 6;
for (row = 0; row <= nrows; row++) {
+ G_percent(row, nrows, 2);
for (col = 0; col <= ncols; col++) {
if ((raster_matrix[row][col].clas >= Thres_j) &&
@@ -416,7 +475,7 @@
punti_bordo, &lungPunti, row, col,
colorBordo, Thres_j, MaxPoints);
- /*CONVEX-HULL COMPUTATION */
+ /* CONVEX-HULL COMPUTATION */
lungHull = ch2d(P, lungPunti);
cvxHull = G_alloc_matrix(lungHull, 3);
Modified: grass/branches/develbranch_6/vector/lidar/v.outlier/main.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.outlier/main.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.outlier/main.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -20,6 +20,7 @@
/*INCLUDES*/
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include <grass/config.h>
#include <grass/gis.h>
#include <grass/Vect.h>
@@ -28,16 +29,22 @@
#include <grass/PolimiFunct.h>
#include "outlier.h"
/* GLOBAL VARIABLES DEFINITIONS */
-int nsply, nsplx, first_it;
+int nsply, nsplx;
double passoN, passoE, Thres_Outlier;
/*--------------------------------------------------------------------------------------*/
int main(int argc, char *argv[])
{
/* Variables' declarations */
+ int nsplx_adj, nsply_adj;
+ int nsubregion_col, nsubregion_row;
+ int subregion = 0, nsubregions = 0;
+ double N_extension, E_extension, orloE, orloN;
int dim_vect, nparameters, BW, npoints;
- double ew_resol, ns_resol, mean, lambda;
- char *dvr, *db, *mapset, table_name[1024];
+ double mean, lambda;
+ const char *dvr, *db, *mapset;
+ char table_name[GNAME_MAX];
+ char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
int last_row, last_column, flag_auxiliar = FALSE;
@@ -49,7 +56,7 @@
struct Map_info In, Out, Outlier, Qgis;
struct Option *in_opt, *out_opt, *outlier_opt, *qgis_opt, *passoE_opt,
*passoN_opt, *lambda_f_opt, *Thres_O_opt;
- /*struct Flag *qgis_flag; */
+ struct Flag *spline_step_flag;
struct GModule *module;
struct Cell_head elaboration_reg, original_reg;
@@ -66,6 +73,12 @@
module->keywords = _("vector, statistics");
module->description = _("Removes outliers from vector point data.");
+ spline_step_flag = G_define_flag();
+ spline_step_flag->key = 'e';
+ spline_step_flag->label = _("Estimate point density and distance");
+ spline_step_flag->description =
+ _("Estimate point density and distance for the input vector points within the current region extends and quit");
+
in_opt = G_define_standard_option(G_OPT_V_INPUT);
out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
@@ -133,9 +146,7 @@
lambda = atof(lambda_f_opt->answer);
Thres_Outlier = atof(Thres_O_opt->answer);
- /* Setting auxiliar table's name */
- /* sprintf(table_name, "%s_aux", out_opt->answer); */
- sprintf(table_name, "Auxiliar_outlier_table");
+ flag_auxiliar = FALSE;
/* Checking vector names */
Vect_check_input_output_name(in_opt->answer, out_opt->answer,
@@ -145,6 +156,49 @@
G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
}
+ /* Setting auxiliar table's name */
+ if (G__name_is_fully_qualified(out_opt->answer, xname, xmapset)) {
+ sprintf(table_name, "%s_aux", xname);
+ }
+ else
+ sprintf(table_name, "%s_aux", out_opt->answer);
+
+ /* Something went wrong in a previous v.outlier execution */
+ if (db_table_exists(dvr, db, table_name)) {
+ /* Start driver and open db */
+ driver = db_start_driver_open_database(dvr, db);
+ if (driver == NULL)
+ G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
+ dvr);
+ if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
+ G_fatal_error(_("Old auxiliar table could not be dropped"));
+ db_close_database_shutdown_driver(driver);
+ }
+
+ Vect_set_open_level(1);
+ /* Open input vector */
+ if (1 > Vect_open_old(&In, in_opt->answer, mapset))
+ G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
+ in_opt->answer);
+
+ /* Input vector must be 3D */
+ if (!Vect_is_3d(&In))
+ G_fatal_error(_("Input vector map <%s> is not 3D!"), in_opt->answer);
+
+ /* Estimate point density and mean distance for current region */
+ if (spline_step_flag->answer) {
+ double dens, dist;
+ if (P_estimate_splinestep(&In, &dens, &dist) == 0) {
+ G_message("Estimated point density: %.4g", dens);
+ G_message("Estimated mean distance between points: %.4g", dist);
+ }
+ else
+ G_warning(_("No points in current region!"));
+
+ Vect_close(&In);
+ exit(EXIT_SUCCESS);
+ }
+
/* Open output vector */
if (qgis_opt->answer)
if (0 > Vect_open_new(&Qgis, qgis_opt->answer, WITHOUT_Z))
@@ -162,12 +216,6 @@
G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
}
- Vect_set_open_level(1);
- /* Open input vector */
- if (1 > Vect_open_old(&In, in_opt->answer, mapset))
- G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
- in_opt->answer);
-
/* Copy vector Head File */
Vect_copy_head_data(&In, &Out);
Vect_hist_copy(&In, &Out);
@@ -189,37 +237,64 @@
G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
dvr);
- /* Something went wrong in a previous v.outlier execution */
- if (db_table_exists(dvr, db, table_name)) {
- if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
- G_fatal_error(_("Old auxiliar table could not be dropped"));
- }
+ /* Create auxiliar table */
+ if ((flag_auxiliar =
+ P_Create_Aux2_Table(driver, table_name)) == FALSE)
+ G_fatal_error(_("It was impossible to create <%s> table."), table_name);
+ db_create_index2(driver, table_name, "ID");
+ /* sqlite likes that */
+ db_close_database_shutdown_driver(driver);
+ driver = db_start_driver_open_database(dvr, db);
+
/* Setting regions and boxes */
G_get_set_window(&original_reg);
G_get_set_window(&elaboration_reg);
Vect_region_box(&elaboration_reg, &overlap_box);
Vect_region_box(&elaboration_reg, &general_box);
+ /*------------------------------------------------------------------
+ | Subdividing and working with tiles:
+ | Each original region will be divided into several subregions.
+ | Each one will be overlaped by its neighbouring subregions.
+ | The overlapping is calculated as a fixed OVERLAP_SIZE times
+ | the largest spline step plus 2 * orlo
+ ----------------------------------------------------------------*/
+
/* Fixing parameters of the elaboration region */
- /*! Each original_region will be divided into several subregions. These
- * subregion will be overlapped by its neibourgh subregions. This overlapping
- * is calculated as OVERLAP_PASS times the east-west resolution. */
+ P_zero_dim(&dims);
- ew_resol = original_reg.ew_res;
- ns_resol = original_reg.ns_res;
+ nsplx_adj = NSPLX_MAX;
+ nsply_adj = NSPLY_MAX;
+ if (passoN > passoE)
+ dims.overlap = OVERLAP_SIZE * passoN;
+ else
+ dims.overlap = OVERLAP_SIZE * passoE;
+ P_get_orlo(P_BILINEAR, &dims, passoE, passoN);
+ P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
- P_zero_dim(&dims);
+ G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
+ G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
- dims.latoE = NSPLX_MAX * passoE;
- dims.latoN = NSPLY_MAX * passoN;
- dims.overlap = OVERLAP_SIZE * ew_resol;
- P_get_orlo(P_BICUBIC, &dims, passoE, passoN);
+ /* calculate number of subregions */
+ orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
+ orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+ N_extension = original_reg.north - original_reg.south;
+ E_extension = original_reg.east - original_reg.west;
+
+ nsubregion_col = ceil(E_extension / orloE) + 0.5;
+ nsubregion_row = ceil(N_extension / orloN) + 0.5;
+
+ if (nsubregion_col < 0)
+ nsubregion_col = 0;
+ if (nsubregion_row < 0)
+ nsubregion_row = 0;
+
+ nsubregions = nsubregion_row * nsubregion_col;
+
elaboration_reg.south = original_reg.north;
-
last_row = FALSE;
- first_it = TRUE;
while (last_row == FALSE) { /* For each row */
@@ -239,10 +314,12 @@
nsply =
ceil((elaboration_reg.north - elaboration_reg.south) / passoN) +
- 1;
+ 0.5;
+ /*
if (nsply > NSPLY_MAX) {
nsply = NSPLY_MAX;
}
+ */
G_debug(1, "nsply = %d", nsply);
elaboration_reg.east = original_reg.west;
@@ -250,6 +327,10 @@
while (last_column == FALSE) { /* For each column */
+ subregion++;
+ if (nsubregions > 1)
+ G_message(_("subregion %d of %d"), subregion, nsubregions);
+
P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
GENERAL_COLUMN);
@@ -266,10 +347,12 @@
nsplx =
ceil((elaboration_reg.east - elaboration_reg.west) / passoE) +
- 1;
+ 0.5;
+ /*
if (nsplx > NSPLX_MAX) {
nsplx = NSPLX_MAX;
}
+ */
G_debug(1, "nsplx = %d", nsplx);
/*Setting the active region */
@@ -305,7 +388,9 @@
Q[i] = 1; /* Q=I */
}
- G_debug(1, "Bilinear interpolation");
+ G_free(observ);
+
+ G_verbose_message(_("Bilinear interpolation"));
normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
nsply, elaboration_reg.west,
elaboration_reg.south, npoints, nparameters,
@@ -317,22 +402,15 @@
G_free_vector(TN);
G_free_vector(Q);
- if (flag_auxiliar == FALSE) {
- G_debug(1,
- "Creating auxiliar table for archiving overlapping zones");
- if ((flag_auxiliar =
- P_Create_Aux_Table(driver, table_name)) == FALSE)
- G_fatal_error(_("It was impossible to create <Auxiliar_outlier_table>."));
- }
-
+ G_verbose_message(_("Outlier detection"));
if (qgis_opt->answer)
P_Outlier(&Out, &Outlier, &Qgis, elaboration_reg,
general_box, overlap_box, obsVect, parVect,
- mean, dims.overlap, lineVect, npoints, driver);
+ mean, dims.overlap, lineVect, npoints, driver, table_name);
else
P_Outlier(&Out, &Outlier, NULL, elaboration_reg,
general_box, overlap_box, obsVect, parVect,
- mean, dims.overlap, lineVect, npoints, driver);
+ mean, dims.overlap, lineVect, npoints, driver, table_name);
G_free_vector(parVect);
@@ -340,8 +418,11 @@
G_free_ivector(lineVect);
} /*! END IF; npoints > 0 */
- G_free(observ);
- first_it = FALSE;
+ else {
+ G_free(observ);
+ G_warning(_("No data within this subregion. "
+ "Consider changing the spline step."));
+ }
} /*! END WHILE; last_column = TRUE */
} /*! END WHILE; last_row = TRUE */
Modified: grass/branches/develbranch_6/vector/lidar/v.outlier/outlier.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.outlier/outlier.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.outlier/outlier.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -1,31 +1,36 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include <grass/gis.h>
#include <grass/Vect.h>
#include <grass/dbmi.h>
#include "outlier.h"
+extern double Thres_Outlier;
+
void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
struct Map_info *Qgis, struct Cell_head Elaboration,
- BOUND_BOX General, BOUND_BOX Overlap, double **obs,
- double *parBilin, double mean, double overlap, int *line_num,
- int num_points, dbDriver * driver)
+ BOUND_BOX General, BOUND_BOX Overlap,
+ double **obs, double *parBilin, double mean, double overlap,
+ int *line_num, int num_points, dbDriver * driver,
+ char *tab_name)
{
int i;
- double interpolation, weight, residual, eta, csi, gradient[2];
-
+ double interpolation, weight, residual, eta, csi;
extern int nsplx, nsply;
extern double passoN, passoE;
-
struct line_pnts *point;
struct line_cats *categories;
point = Vect_new_line_struct();
categories = Vect_new_cats_struct();
+ db_begin_transaction(driver);
+
for (i = 0; i < num_points; i++) { /* Sparse points */
+ G_percent(i, num_points, 2);
Vect_reset_line(point);
Vect_reset_cats(categories);
@@ -57,60 +62,53 @@
}
else {
- if ((*point->x > Overlap.E) && (*point->x != General.E)) {
+ if ((*point->x > Overlap.E) && (*point->x < General.E)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ csi = (General.E - *point->x) / overlap;
+ eta = (General.N - *point->y) / overlap;
+ weight = csi * eta;
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(3) */
- csi = (*point->x - Overlap.E) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - csi) * (1 - eta);
-
- gradient[0] *= weight;
- gradient[1] *= weight;
interpolation *= weight;
- if (Select_Outlier(&gradient[0], line_num[i], driver)
- != DB_OK)
+ if (Select_Outlier(&interpolation, line_num[i],
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
- if (UpDate_Outlier(gradient[0], line_num[i], driver)
- != DB_OK)
+ if (UpDate_Outlier(interpolation, line_num[i],
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to update the database"));
-
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(1) */
- csi = (*point->x - Overlap.E) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ csi = (General.E - *point->x) / overlap;
eta = (*point->y - General.S) / overlap;
- weight = (1 - csi) * eta;
+ weight = csi * eta;
- if (Insert_Outlier
- (interpolation * weight, line_num[i],
- driver) != DB_OK)
+ interpolation *= weight;
+
+ if (Insert_Outlier(interpolation, line_num[i],
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to write in the database"));
-
}
else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(1) */
- weight = (*point->x - Overlap.E) / overlap;
+ weight = (General.E - *point->x) / overlap;
- if (Insert_Outlier
- (interpolation * weight, line_num[i],
- driver) != DB_OK)
+ interpolation *= weight;
+
+ if (Insert_Outlier(interpolation, line_num[i],
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to write in the database"));
}
-
}
- else if ((*point->x < Overlap.W) && (*point->x != General.W)) {
-
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(4) */
+ else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(4) */
csi = (*point->x - General.W) / overlap;
- eta = (*point->y - Overlap.N) / overlap;
- weight = (1 - eta) * csi;
+ eta = (General.N - *point->y) / overlap;
+ weight = eta * csi;
- gradient[0] *= weight;
- gradient[1] *= weight;
interpolation *= weight;
- if (Select_Outlier(&gradient[0], line_num[i], driver)
- != DB_OK)
+ if (Select_Outlier(&interpolation, line_num[i],
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
residual = *point->z - interpolation;
@@ -122,38 +120,34 @@
Vect_write_line(Qgis, GV_POINT, point,
categories);
}
- else
+ else {
Vect_write_line(Outlier, GV_POINT, point,
categories);
-
+ G_message("here we are");
+ }
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(2) */
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(2) */
csi = (*point->x - General.W) / overlap;
eta = (*point->y - General.S) / overlap;
weight = csi * eta;
- gradient[0] *= weight;
- gradient[1] *= weight;
interpolation *= weight;
- if (Select_Outlier(&gradient[0], line_num[i], driver)
- != DB_OK)
+ if (Select_Outlier(&interpolation, line_num[i],
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
- if (UpDate_Outlier(gradient[0], line_num[i], driver)
- != DB_OK)
+ if (UpDate_Outlier(interpolation, line_num[i],
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to update the database"));
-
}
else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(2) */
- weight = (Overlap.W - *point->x) / overlap;
+ weight = (*point->x - General.W) / overlap;
- gradient[0] *= weight;
- gradient[1] *= weight;
interpolation *= weight;
- if (Select_Outlier(&gradient[0], line_num[i], driver)
- != DB_OK)
+ if (Select_Outlier(&interpolation, line_num[i],
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
residual = *point->z - interpolation;
@@ -169,18 +163,15 @@
Vect_write_line(Outlier, GV_POINT, point,
categories);
}
-
}
else if ((*point->x <= Overlap.E) && (*point->x >= Overlap.W)) {
- if ((*point->y > Overlap.N) && (*point->y != General.N)) { /*(3) */
- weight = (*point->y - Overlap.N) / overlap;
+ if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
+ weight = (General.N - *point->y) / overlap;
- gradient[0] *= weight;
- gradient[1] *= weight;
interpolation *= weight;
- if (Select_Outlier(&gradient[0], line_num[i], driver)
- != DB_OK)
+ if (Select_Outlier(&interpolation, line_num[i],
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to read the database"));
residual = *point->z - interpolation;
@@ -195,14 +186,14 @@
else
Vect_write_line(Outlier, GV_POINT, point,
categories);
-
}
- else if ((*point->y < Overlap.S) && (*point->y != General.S)) { /*(1) */
- weight = (Overlap.S - *point->y) / overlap;
+ else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
+ weight = (*point->y - General.S) / overlap;
- if (Insert_Outlier
- (interpolation * weight, line_num[i], driver)
- != DB_OK)
+ interpolation *= weight;
+
+ if (Insert_Outlier(interpolation, line_num[i],
+ driver, tab_name) != DB_OK)
G_fatal_error(_("Impossible to write in the database"));
} /*else (1) */
@@ -210,42 +201,56 @@
}
} /*end if obs */
} /*end for */
+
+ G_percent(num_points, num_points, 2);
+ G_debug(2, "P_outlier: done");
+
+ db_commit_transaction(driver);
+
Vect_destroy_line_struct(point);
Vect_destroy_cats_struct(categories);
} /*end puntisparsi_select */
-int Insert_Outlier(double Interp, int line_num, dbDriver * driver)
+int Insert_Outlier(double Interp, int line_num, dbDriver * driver,
+ char *tab_name)
{
-
char buf[1024];
dbString sql;
+ int ret;
db_init_string(&sql);
- sprintf(buf, "INSERT INTO Auxiliar_outlier_table (ID, Interp)");
+ sprintf(buf, "INSERT INTO %s (ID, Interp)", tab_name);
db_append_string(&sql, buf);
sprintf(buf, " VALUES (%d, %lf)", line_num, Interp);
db_append_string(&sql, buf);
- return db_execute_immediate(driver, &sql);
+ ret = db_execute_immediate(driver, &sql);
+ db_free_string(&sql);
+
+ return ret;
}
-int UpDate_Outlier(double Interp, int line_num, dbDriver * driver)
+int UpDate_Outlier(double Interp, int line_num, dbDriver * driver,
+ char *tab_name)
{
-
char buf[1024];
dbString sql;
+ int ret;
db_init_string(&sql);
- sprintf(buf, "UPDATE Auxiliar_outlier_table SET Interp=%lf WHERE ID=%d",
- Interp, line_num);
+ sprintf(buf, "UPDATE %s SET Interp=%lf WHERE ID=%d",
+ tab_name, Interp, line_num);
db_append_string(&sql, buf);
- return db_execute_immediate(driver, &sql);
+ ret = db_execute_immediate(driver, &sql);
+ db_free_string(&sql);
+
+ return ret;
}
-int Select_Outlier(double *Interp, int line_num, dbDriver * driver)
+int Select_Outlier(double *Interp, int line_num, dbDriver * driver,
+ char *tab_name)
{
-
int more;
char buf[1024];
dbString sql;
@@ -255,8 +260,7 @@
dbValue *Interp_value;
db_init_string(&sql);
- sprintf(buf, "SELECT ID, Interp FROM Auxiliar_outlier_table WHERE ID=%d",
- line_num);
+ sprintf(buf, "SELECT ID, Interp FROM %s WHERE ID=%d", tab_name, line_num);
db_append_string(&sql, buf);
if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
@@ -276,36 +280,42 @@
*Interp += db_get_value_double(Interp_value);
}
db_close_cursor(&cursor);
+ db_free_string(&sql);
return DB_OK;
}
int P_is_outlier(double pippo)
{
- extern double Thres_Outlier;
-
if (fabs(pippo) < Thres_Outlier)
return FALSE;
return TRUE;
}
-#ifdef notdef
/*! DEFINITION OF THE SUBZONES
- -----------------------
- |4| 3 |3| | |
- -----------------------
- | | | | | |
- |2| 5 |1| | |
- | | | | | |
- -----------------------
- |2| 1 |1| | |
- -----------------------
- | | | | | |
- | | | | | |
- | | | | | |
- -----------------------
- | | | | | |
- -----------------------
+ 5: inside Overlap region
+ all others: inside General region but outside Overlap region
+
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | |4| 3 |3| | |
+ ---------------------------------
+ | | | | | | | |
+ | | |2| 5 |1| | |
+ | | | | | | | |
+ ---------------------------------
+ | | |2| 1 |1| | |
+ ---------------------------------
+ | | | | | | | |
+ | | | | | | | |
+ | | | | | | | |
+ ---------------------------------
+ | | | | | | | |
+ ---------------------------------
*/
-#endif
Modified: grass/branches/develbranch_6/vector/lidar/v.outlier/outlier.h
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.outlier/outlier.h 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.outlier/outlier.h 2010-02-17 17:56:52 UTC (rev 41062)
@@ -16,12 +16,12 @@
BOUND_BOX, /**/
double **, /**/
double *, /**/
- double, /**/ double, /**/ int *, /**/ int, /**/ dbDriver * /**/);
+ double, /**/ double, /**/ int *, /**/ int, /**/ dbDriver *, /**/ char *);
-int Insert_Outlier(double, /**/ int, /**/ dbDriver * /**/);
+int Insert_Outlier(double, /**/ int, /**/ dbDriver *, /**/ char *);
-int UpDate_Outlier(double, /**/ int, /**/ dbDriver * /**/);
+int UpDate_Outlier(double, /**/ int, /**/ dbDriver *, /**/ char *);
-int Select_Outlier(double *, /**/ int, /**/ dbDriver * /**/);
+int Select_Outlier(double *, /**/ int, /**/ dbDriver *, /**/ char *);
int P_is_outlier(double);
Modified: grass/branches/develbranch_6/vector/lidar/v.surf.bspline/crosscorr.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.surf.bspline/crosscorr.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.surf.bspline/crosscorr.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -20,6 +20,7 @@
/*INCLUDES*/
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include <grass/gis.h>
#include <grass/Vect.h>
#include <grass/dbmi.h>
@@ -48,7 +49,8 @@
int nsplx, nsply, nparam_spl, ndata;
double *mean, *rms, *stdev, rms_min, stdev_min;
- double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; /* Fixed values (by the moment) */
+ /* double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; */ /* Fixed values (by the moment) */
+ double lambda[PARAM_LAMBDA] = { 0.01, 0.05, 0.1, 0.2, 0.3, 0.4 }; /* Fixed values (by the moment) */
double *TN, *Q, *parVect; /* Interpolation and least-square vectors */
double **N, **obsVect; /* Interpolation and least-square matrix */
@@ -151,7 +153,7 @@
nsplx, nsply);
BW = P_get_BandWidth(bilin, nsply);
- /**/
+ /**/
/*Least Squares system */
N = G_alloc_matrix(nparam_spl, BW); /* Normal matrix */
TN = G_alloc_vector(nparam_spl); /* vector */
Modified: grass/branches/develbranch_6/vector/lidar/v.surf.bspline/main.c
===================================================================
--- grass/branches/develbranch_6/vector/lidar/v.surf.bspline/main.c 2010-02-17 17:05:36 UTC (rev 41061)
+++ grass/branches/develbranch_6/vector/lidar/v.surf.bspline/main.c 2010-02-17 17:56:52 UTC (rev 41062)
@@ -1,5 +1,5 @@
-/***********************************************************************
+/**********************************************************************
*
* MODULE: v.surf.bspline
*
@@ -15,11 +15,12 @@
* Read the file COPYING that comes with GRASS
* for details.
*
- **************************************************************************/
+ **********************************************************************/
- /*INCLUDES*/
+/* INCLUDES */
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include <grass/gis.h>
#include <grass/Vect.h>
#include <grass/dbmi.h>
@@ -27,22 +28,25 @@
#include <grass/config.h>
#include <grass/PolimiFunct.h>
#include "bspline.h"
- /*GLOBAL VARIABLES */
+
+/* GLOBAL VARIABLES */
int bspline_field;
char *bspline_column;
-/*-------------------------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------*/
int main(int argc, char *argv[])
{
-
- /* Variables' declarations */
- int nsply, nsplx, nlines, nrows, ncols;
+ /* Variable declarations */
+ int nsply, nsplx, nrows, ncols, nsplx_adj, nsply_adj;
int nsubregion_col, nsubregion_row, subregion_row, subregion_col;
+ int subregion = 0, nsubregions = 0;
int last_row, last_column, grid, bilin, ext, flag_auxiliar, cross; /* booleans */
double passoN, passoE, lambda, mean;
double N_extension, E_extension, orloE, orloN;
- char *mapset, *dvr, *db, *vector, *map, table_name[1024], title[64];
+ const char *mapset, *dvr, *db, *vector, *map;
+ char table_name[GNAME_MAX], title[64];
+ char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
int dim_vect, nparameters, BW;
int *lineVect; /* Vector restoring primitive's ID */
@@ -50,15 +54,15 @@
double *TN, *Q, *parVect; /* Interpolating and least-square vectors */
double **N, **obsVect; /* Interpolation and least-square matrix */
- /* Structs' declarations */
+ /* Structs declarations */
int raster;
struct Map_info In, In_ext, Out;
struct History history;
struct GModule *module;
struct Option *in_opt, *in_ext_opt, *out_opt, *out_map_opt, *passoE_opt,
- *passoN_opt, *lambda_f_opt, *type, *dfield_opt, *col_opt;
- struct Flag *cross_corr_flag;
+ *passoN_opt, *lambda_f_opt, *type_opt, *dfield_opt, *col_opt;
+ struct Flag *cross_corr_flag, *spline_step_flag;
struct Reg_dimens dims;
struct Cell_head elaboration_reg, original_reg;
@@ -73,28 +77,31 @@
struct field_info *Fi;
dbDriver *driver, *driver_cats;
- /*-------------------------------------------------------------------------------------------*/
- /* Options' declaration */
- module = G_define_module(); {
- module->keywords = _("vector, interpolation");
- module->description =
- _("Bicubic or bilinear spline interpolation with Tykhonov regularization.");
- }
+ /*----------------------------------------------------------------*/
+ /* Options declarations */
+ module = G_define_module();
+ module->keywords = _("vector, interpolation");
+ module->description =
+ _("Bicubic or bilinear spline interpolation with Tykhonov regularization.");
- cross_corr_flag = G_define_flag(); {
- cross_corr_flag->key = 'c';
- cross_corr_flag->description =
- _("Find the best Tykhonov regularizing parameter using a \"leave-one-out\" cross validation method");
- }
+ cross_corr_flag = G_define_flag();
+ cross_corr_flag->key = 'c';
+ cross_corr_flag->description =
+ _("Find the best Tykhonov regularizing parameter using a \"leave-one-out\" cross validation method");
+ spline_step_flag = G_define_flag();
+ spline_step_flag->key = 'e';
+ spline_step_flag->label = _("Estimate point density and distance");
+ spline_step_flag->description =
+ _("Estimate point density and distance for the input vector points within the current region extends and quit");
+
in_opt = G_define_standard_option(G_OPT_V_INPUT);
- in_ext_opt = G_define_standard_option(G_OPT_V_INPUT); {
- in_ext_opt->key = "sparse";
- in_ext_opt->required = NO;
- in_ext_opt->description =
- _("Name of input vector map of sparse points");
- }
+ in_ext_opt = G_define_standard_option(G_OPT_V_INPUT);
+ in_ext_opt->key = "sparse";
+ in_ext_opt->required = NO;
+ in_ext_opt->description =
+ _("Name of input vector map of sparse points");
out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
out_opt->required = NO;
@@ -103,69 +110,71 @@
out_map_opt->key = "raster";
out_map_opt->required = NO;
- passoE_opt = G_define_option(); {
- passoE_opt->key = "sie";
- passoE_opt->type = TYPE_DOUBLE;
- passoE_opt->required = NO;
- passoE_opt->answer = "4";
- passoE_opt->description =
- _("Length of each spline step in the east-west direction");
- passoE_opt->guisection = _("Settings");
- }
+ passoE_opt = G_define_option();
+ passoE_opt->key = "sie";
+ passoE_opt->type = TYPE_DOUBLE;
+ passoE_opt->required = NO;
+ passoE_opt->answer = "4";
+ passoE_opt->description =
+ _("Length of each spline step in the east-west direction");
+ passoE_opt->guisection = _("Settings");
- passoN_opt = G_define_option(); {
- passoN_opt->key = "sin";
- passoN_opt->type = TYPE_DOUBLE;
- passoN_opt->required = NO;
- passoN_opt->answer = "4";
- passoN_opt->description =
- _("Length of each spline step in the north-south direction");
- passoN_opt->guisection = _("Settings");
- }
+ passoN_opt = G_define_option();
+ passoN_opt->key = "sin";
+ passoN_opt->type = TYPE_DOUBLE;
+ passoN_opt->required = NO;
+ passoN_opt->answer = "4";
+ passoN_opt->description =
+ _("Length of each spline step in the north-south direction");
+ passoN_opt->guisection = _("Settings");
- type = G_define_option(); {
- type->key = "type";
- type->type = TYPE_STRING;
- type->required = NO;
- type->description = _("Spline interpolation algorithm");
- type->options = "bilinear,bicubic";
- type->answer = "bilinear";
- type->guisection = _("Settings");
- }
+ type_opt = G_define_option();
+ type_opt->key = "method";
+ type_opt->type = TYPE_STRING;
+ type_opt->required = NO;
+ type_opt->description = _("Spline interpolation algorithm");
+ type_opt->options = "bilinear,bicubic";
+ type_opt->answer = "bilinear";
+ type_opt->guisection = _("Settings");
- lambda_f_opt = G_define_option(); {
- lambda_f_opt->key = "lambda_i";
- lambda_f_opt->type = TYPE_DOUBLE;
- lambda_f_opt->required = NO;
- lambda_f_opt->description = _("Tykhonov regularization parameter (affects smoothing)");
- lambda_f_opt->answer = "1";
- lambda_f_opt->guisection = _("Settings");
- }
+ lambda_f_opt = G_define_option();
+ lambda_f_opt->key = "lambda_i";
+ lambda_f_opt->type = TYPE_DOUBLE;
+ lambda_f_opt->required = NO;
+ lambda_f_opt->description = _("Tykhonov regularization parameter (affects smoothing)");
+ lambda_f_opt->answer = "1";
+ lambda_f_opt->guisection = _("Settings");
- dfield_opt = G_define_standard_option(G_OPT_V_FIELD); {
- dfield_opt->description =
- _("If set to 0, z coordinates are used. (3D vector only)");
- dfield_opt->answer = "0";
- dfield_opt->guisection = _("Settings");
- dfield_opt->gisprompt = "old_layer,layer,layer_zero";
- }
+ dfield_opt = G_define_standard_option(G_OPT_V_FIELD);
+ dfield_opt->description =
+ _("If set to -1, z coordinates are used. (3D vector only)");
+ dfield_opt->answer = "-1";
+ dfield_opt->guisection = _("Settings");
- col_opt = G_define_option(); {
- col_opt->key = "column";
- col_opt->type = TYPE_STRING;
- col_opt->required = NO;
- col_opt->description =
- _("Attribute table column with values to interpolate (if layer>0)");
- col_opt->guisection = _("Settings");
- }
+ col_opt = G_define_option();
+ col_opt->key = "column";
+ col_opt->type = TYPE_STRING;
+ col_opt->required = NO;
+ col_opt->description =
+ _("Attribute table column with values to interpolate (if layer>0)");
+ col_opt->guisection = _("Settings");
-/*-------------------------------------------------------------------------------------------*/
+ /*----------------------------------------------------------------*/
/* Parsing */
G_gisinit(argv[0]);
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
- if (!strcmp(type->answer, "bilinear"))
+ vector = out_opt->answer;
+ map = out_map_opt->answer;
+
+ if (vector && map)
+ G_fatal_error(_("Choose either vector or raster output, not both"));
+
+ if (!vector && !map && !cross_corr_flag->answer)
+ G_fatal_error(_("No raster or vector or cross-validation output"));
+
+ if (!strcmp(type_opt->answer, "bilinear"))
bilin = P_BILINEAR;
else
bilin = P_BICUBIC;
@@ -175,12 +184,8 @@
lambda = atof(lambda_f_opt->answer);
bspline_field = atoi(dfield_opt->answer);
bspline_column = col_opt->answer;
-
+
flag_auxiliar = FALSE;
- if (cross_corr_flag->answer) {
- }
- vector = out_opt->answer;
- map = out_map_opt->answer;
if (!(db = G__getenv2("DB_DATABASE", G_VAR_MAPSET)))
G_fatal_error(_("Unable to read name of database"));
@@ -189,39 +194,26 @@
G_fatal_error(_("Unable to read name of driver"));
/* Setting auxiliar table's name */
- if (vector)
- sprintf(table_name, "%s_aux", out_opt->answer);
+ if (vector) {
+ if (G__name_is_fully_qualified(out_opt->answer, xname, xmapset)) {
+ sprintf(table_name, "%s_aux", xname);
+ }
+ else
+ sprintf(table_name, "%s_aux", out_opt->answer);
+ }
- /* Open driver and database */
- if (db_table_exists(dvr, db, table_name)) { /*Something went wrong in a
- previous v.surf.bspline execution */
- dbString sql;
- char buf[1024];
-
-
+ /* Something went wrong in a previous v.surf.bspline execution */
+ if (db_table_exists(dvr, db, table_name)) {
+ /* Start driver and open db */
driver = db_start_driver_open_database(dvr, db);
if (driver == NULL)
- G_fatal_error(_("No database connection for driver <%s> is defined. "
- "Run db.connect."), dvr);
-
- db_init_string(&sql);
- db_zero_string(&sql);
- sprintf(buf, "drop table %s", table_name);
- db_append_string(&sql, buf);
- if (db_execute_immediate(driver, &sql) != DB_OK)
- G_fatal_error(_("It was not possible to drop <%s> table. "
- "Nothing will be done. Try to drop it manually."),
- table_name);
-
+ G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
+ dvr);
+ if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
+ G_fatal_error(_("Old auxiliar table could not be dropped"));
db_close_database_shutdown_driver(driver);
}
- if (vector && map)
- G_fatal_error(_("Choose either vector or raster output, not both"));
-#ifdef notdef
- if (!vector && !map && !cross_corr_flag->answer)
- G_fatal_error(_("No raster nor vector output"));
-#endif
/* Open input vector */
if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL)
G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
@@ -231,6 +223,45 @@
G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
in_opt->answer);
+ /* check availability of z values
+ * column option overrrides 3D z coordinates */
+ if (!Vect_is_3d(&In) && (bspline_field <= 0 || bspline_column == NULL))
+ G_fatal_error(_("Need either 3D vector or layer and column with z values"));
+ if (bspline_field > 0 && bspline_column == NULL)
+ G_fatal_error(_("Layer but not column with z values given"));
+
+ /* Estimate point density and mean distance for current region */
+ if (spline_step_flag->answer) {
+ double dens, dist;
+ if (P_estimate_splinestep(&In, &dens, &dist) == 0) {
+ G_message("Estimated point density: %.4g", dens);
+ G_message("Estimated mean distance between points: %.4g", dist);
+ }
+ else
+ G_warning(_("No points in current region!"));
+
+ Vect_close(&In);
+ exit(EXIT_SUCCESS);
+ }
+
+ /*----------------------------------------------------------------*/
+ /* Cross-correlation begins */
+ if (cross_corr_flag->answer) {
+ G_debug(1, "CrossCorrelation()");
+ cross = cross_correlation(&In, passoE, passoN);
+
+ if (cross != TRUE)
+ G_fatal_error(_("Cross validation didn't finish correctly"));
+ else {
+ G_debug(1, "Cross validation finished correctly");
+
+ Vect_close(&In);
+
+ G_done_msg(_("Cross validation finished for sie = %f and sin = %f"), passoE, passoN);
+ exit(EXIT_SUCCESS);
+ }
+ }
+
/* Open input ext vector */
if (!in_ext_opt->answer) {
ext = FALSE;
@@ -284,16 +315,12 @@
G_set_fp_type(DCELL_TYPE);
if (!vector && map) {
grid = TRUE;
- /*
- if (G_find_cell (out_map_opt->answer, G_mapset()) != NULL)
- G_fatal_error (_("Raster <%s> already exist."), out_map_opt->answer);
- */
-
if ((raster = G_open_fp_cell_new(out_map_opt->answer)) < 0)
G_fatal_error(_("Unable to create raster map <%s>"),
out_map_opt->answer);
}
+ /* read z values from attribute table */
if (bspline_field > 0) {
db_CatValArray_init(&cvarr);
Fi = Vect_get_field(&In, bspline_field);
@@ -324,54 +351,30 @@
db_close_database_shutdown_driver(driver_cats);
}
-/*-------------------------------------------------------------------------------------------*/
- /*Cross-correlation begins */
- if (cross_corr_flag->answer) { /* CROSS-CORRELATION WILL BE DONE */
- G_debug(1, "CrossCorrelation()");
- /*cross = cross_correlation (&In, &passoE, &passoN, &lambda); */
- cross = cross_correlation(&In, passoE, passoN);
-
- if (cross != TRUE)
- G_fatal_error(_("Cross validation didn't finish correctly"));
- else {
- G_debug(1, "Cross validation finished correctly");
-
- Vect_close(&In);
- if (ext != FALSE)
- Vect_close(&In_ext);
- if (vector)
- Vect_close(&Out);
-
- if (map)
- G_close_cell(raster);
-
- G_done_msg(_("Cross Validation was success!"));
- exit(EXIT_SUCCESS);
- }
-
- /*G_debug (0, _("passoE = %f, passoN=%f, lambda_i=%f"), passoE, passoN, lambda); */
- }
-
-/*-------------------------------------------------------------------------------------------*/
+ /*----------------------------------------------------------------*/
/* Interpolation begins */
G_debug(1, "Interpolation()");
- if (map) {
- nrows = G_window_rows();
- ncols = G_window_cols();
-
- if (!G_alloc_matrix(nrows, ncols))
- G_fatal_error(_("Interpolation: The region resolution is too high: "
- "%d cells. Consider to change it."),
- nrows * ncols);
- }
-
/* Open driver and database */
driver = db_start_driver_open_database(dvr, db);
if (driver == NULL)
G_fatal_error(_("No database connection for driver <%s> is defined. "
"Run db.connect."), dvr);
+ /* Create auxiliar table */
+ if (vector) {
+ if ((flag_auxiliar = P_Create_Aux4_Table(driver, table_name)) == FALSE) {
+ P_Drop_Aux_Table(driver, table_name);
+ G_fatal_error(_("Interpolation: Creating table: "
+ "It was impossible to create table <%s>."),
+ table_name);
+ }
+ db_create_index2(driver, table_name, "ID");
+ /* sqlite likes that */
+ db_close_database_shutdown_driver(driver);
+ driver = db_start_driver_open_database(dvr, db);
+ }
+
/* Setting regions and boxes */
G_debug(1, "Interpolation: Setting regions and boxes");
G_get_window(&original_reg);
@@ -379,49 +382,66 @@
Vect_region_box(&elaboration_reg, &overlap_box);
Vect_region_box(&elaboration_reg, &general_box);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+
/* Alloc raster matrix */
- if (map)
+ if (grid == TRUE) {
if (!(raster_matrix = G_alloc_matrix(nrows, ncols)))
G_fatal_error(_("Cannot allocate memory for auxiliar matrix."
"Consider changing region resolution"));
+ }
- /* Fixxing parameters of the elaboration region */
- P_zero_dim(&dims); /* Set to zero the dim struct */
- dims.latoE = NSPLX_MAX * passoE;
- dims.latoN = NSPLY_MAX * passoN;
- dims.overlap = OVERLAP_SIZE * passoE;
- P_get_orlo(bilin, &dims, passoE, passoN); /* Set the last two dim elements */
+ /*------------------------------------------------------------------
+ | Subdividing and working with tiles:
+ | Each original region will be divided into several subregions.
+ | Each one will be overlaped by its neighbouring subregions.
+ | The overlapping is calculated as a fixed OVERLAP_SIZE times
+ | the largest spline step plus 2 * orlo
+ ----------------------------------------------------------------*/
- N_extension = original_reg.ns_res * original_reg.rows;
- E_extension = original_reg.ew_res * original_reg.cols;
- orloE = dims.latoE - dims.overlap - 2 * dims.orlo_h;
- orloN = dims.latoN - dims.overlap - 2 * dims.orlo_v;
- nsubregion_col = 2 + ((E_extension - dims.latoE) / orloE);
- nsubregion_row = 2 + ((N_extension - dims.latoN) / orloN);
+ /* Fixing parameters of the elaboration region */
+ P_zero_dim(&dims); /* Set dim struct to zero */
+ nsplx_adj = NSPLX_MAX;
+ nsply_adj = NSPLY_MAX;
+ if (passoN > passoE)
+ dims.overlap = OVERLAP_SIZE * passoN;
+ else
+ dims.overlap = OVERLAP_SIZE * passoE;
+ P_get_orlo(bilin, &dims, passoE, passoN);
+ P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
+
+ G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
+ G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
+
+ /* calculate number of subregions */
+ orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
+ orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+
+ N_extension = original_reg.north - original_reg.south;
+ E_extension = original_reg.east - original_reg.west;
+
+ nsubregion_col = ceil(E_extension / orloE) + 0.5;
+ nsubregion_row = ceil(N_extension / orloN) + 0.5;
+
if (nsubregion_col < 0)
nsubregion_col = 0;
if (nsubregion_row < 0)
nsubregion_row = 0;
+ nsubregions = nsubregion_row * nsubregion_col;
+
/* Creating line and categories structs */
points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
- nlines = Vect_get_num_lines(&In);
+ Vect_cat_set(Cats, 1, 0);
- /*------------------------------------------------------------------------------------------
- | Subdividing and working with tiles:
- | Each original_region will be divided into several subregions.
- | Each one will be overlaped by its neibourgh subregions.
- | The overlapping was calculated as a fixed OVERLAP_SIZE times the east-west resolution
- ------------------------------------------------------------------------------------------*/
-
- elaboration_reg.south = original_reg.north;
-
- G_percent(0, 1, 10);
subregion_row = 0;
+ elaboration_reg.south = original_reg.north;
last_row = FALSE;
- while (last_row == FALSE) { /* For each row */
+
+ while (last_row == FALSE) { /* For each subregion row */
subregion_row++;
P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
GENERAL_ROW);
@@ -430,36 +450,35 @@
P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
FIRST_ROW);
- nsply =
- ceil((elaboration_reg.north -
- elaboration_reg.south) / passoN);
- G_debug(1, "Interpolation: nsply = %d", nsply);
- if (nsply > NSPLY_MAX)
- nsply = NSPLY_MAX;
}
if (elaboration_reg.south <= original_reg.south) { /* Last row */
P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
LAST_ROW);
- nsply =
- ceil((elaboration_reg.north -
- elaboration_reg.south) / passoN);
last_row = TRUE;
- G_debug(1, "Interpolation: nsply = %d", nsply);
- if (nsply > NSPLY_MAX)
- nsply = NSPLY_MAX;
}
+ nsply =
+ ceil((elaboration_reg.north -
+ elaboration_reg.south) / passoN) + 0.5;
+ G_debug(1, "Interpolation: nsply = %d", nsply);
+ /*
+ if (nsply > NSPLY_MAX)
+ nsply = NSPLY_MAX;
+ */
elaboration_reg.east = original_reg.west;
last_column = FALSE;
-
subregion_col = 0;
- while (last_column == FALSE) { /* For each column */
+
+ while (last_column == FALSE) { /* For each subregion column */
int npoints = 0;
- int nsubzones = 0, subzone = 0;
subregion_col++;
+ subregion++;
+ if (nsubregions > 1)
+ G_message(_("subregion %d of %d"), subregion, nsubregions);
+
P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
GENERAL_COLUMN);
@@ -467,12 +486,6 @@
P_set_regions(&elaboration_reg, &general_box, &overlap_box,
dims, FIRST_COLUMN);
- nsplx =
- ceil((elaboration_reg.east -
- elaboration_reg.west) / passoE);
- G_debug(1, "Interpolation: nsply = %d", nsply);
- if (nsplx > NSPLX_MAX)
- nsplx = NSPLX_MAX;
}
if (elaboration_reg.east >= original_reg.east) { /* Last column */
@@ -480,13 +493,15 @@
P_set_regions(&elaboration_reg, &general_box, &overlap_box,
dims, LAST_COLUMN);
last_column = TRUE;
- nsplx =
- ceil((elaboration_reg.east -
- elaboration_reg.west) / passoE);
- G_debug(1, "Interpolation: nsply = %d", nsply);
- if (nsplx > NSPLX_MAX)
- nsplx = NSPLX_MAX;
}
+ nsplx =
+ ceil((elaboration_reg.east -
+ elaboration_reg.west) / passoE) + 0.5;
+ G_debug(1, "Interpolation: nsplx = %d", nsplx);
+ /*
+ if (nsplx > NSPLX_MAX)
+ nsplx = NSPLX_MAX;
+ */
G_debug(1, "Interpolation: (%d,%d): subregion bounds",
subregion_row, subregion_col);
G_debug(1, "Interpolation: \t\tNORTH:%.2f\t",
@@ -496,7 +511,7 @@
G_debug(1, "Interpolation: \t\tSOUTH:%.2f",
elaboration_reg.south);
- /*Setting the active region */
+ /* reading points in interpolation region */
dim_vect = nsplx * nsply;
observ =
P_Read_Vector_Region_Map(&In, &elaboration_reg, &npoints,
@@ -507,23 +522,18 @@
if (npoints > 0) { /* */
int i;
- double *obs_mean;
nparameters = nsplx * nsply;
- BW = P_get_BandWidth(bilin, nsply);
+ BW = P_get_BandWidth(P_BILINEAR, nsply) + 2 * (bilin == P_BICUBIC);
- /*Least Squares system */
+ /* Least Squares system */
N = G_alloc_matrix(nparameters, BW); /* Normal matrix */
TN = G_alloc_vector(nparameters); /* vector */
parVect = G_alloc_vector(nparameters); /* Parameters vector */
obsVect = G_alloc_matrix(npoints, 3); /* Observation vector */
Q = G_alloc_vector(npoints); /* "a priori" var-cov matrix */
lineVect = G_alloc_ivector(npoints); /* */
- obs_mean = G_alloc_vector(npoints);
- if (bspline_field <= 0)
- mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
-
for (i = 0; i < npoints; i++) { /* Setting obsVect vector & Q matrix */
double dval;
@@ -532,11 +542,10 @@
obsVect[i][0] = observ[i].coordX;
obsVect[i][1] = observ[i].coordY;
+ /* read z coordinates from attribute table */
if (bspline_field > 0) {
- int cat, ival, ret, type;
+ int cat, ival, ret;
- if (!(type & GV_POINTS))
- continue;
cat = observ[i].cat;
if (cat < 0)
continue;
@@ -546,14 +555,14 @@
db_CatValArray_get_value_int(&cvarr, cat,
&ival);
obsVect[i][2] = ival;
- obs_mean[i] = ival;
+ observ[i].coordZ = ival;
}
else { /* DB_C_TYPE_DOUBLE */
ret =
db_CatValArray_get_value_double(&cvarr, cat,
&dval);
obsVect[i][2] = dval;
- obs_mean[i] = dval;
+ observ[i].coordZ = dval;
}
if (ret != DB_OK) {
G_warning(_("Interpolation: (%d,%d): No record for point (cat = %d)"),
@@ -561,26 +570,25 @@
continue;
}
}
-
+ /* use z coordinates of 3D vector */
else {
obsVect[i][2] = observ[i].coordZ;
- obs_mean[i] = observ[i].coordZ;
- } /*obsVect[i][2] = observ[i].coordZ - mean; */
+ }
}
/* Mean calculation for every point */
- if (bspline_field > 0)
- mean = calc_mean(obs_mean, npoints);
+ mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
G_debug(1, "Interpolation: (%d,%d): mean=%lf",
subregion_row, subregion_col, mean);
+ G_free(observ);
+
for (i = 0; i < npoints; i++)
obsVect[i][2] -= mean;
- G_free(observ);
-
- if (bilin) { /* Bilinear interpolation */
+ /* Bilinear interpolation */
+ if (bilin) {
G_debug(1,
"Interpolation: (%d,%d): Bilinear interpolation...",
subregion_row, subregion_col);
@@ -590,6 +598,7 @@
nparameters, BW);
nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
}
+ /* Bicubic interpolation */
else {
G_debug(1,
"Interpolation: (%d,%d): Bicubic interpolation...",
@@ -607,22 +616,16 @@
G_free_vector(TN);
G_free_vector(Q);
- if (grid == FALSE) { /*OBSERVATION POINTS INTERPOLATION */
- /* Auxiliar table creation */
- if (flag_auxiliar == FALSE) {
- G_debug(1,
- "Interpolation: Creating auxiliar table for archiving "
- "overlapping zones");
- if ((flag_auxiliar =
- P_Create_Aux_Table(driver,
- table_name)) == FALSE) {
- P_Drop_Aux_Table(driver, table_name);
- G_fatal_error(_("Interpolation: Creating table: "
- "It was impossible to create table <%s>."),
- table_name);
- }
- }
-
+ if (grid == TRUE) { /* GRID INTERPOLATION ==> INTERPOLATION INTO A RASTER */
+ G_debug(1, "Interpolation: (%d,%d): Regular_Points...",
+ subregion_row, subregion_col);
+ raster_matrix =
+ P_Regular_Points(&elaboration_reg, general_box,
+ overlap_box, raster_matrix, parVect,
+ passoN, passoE, dims.overlap, mean,
+ nsplx, nsply, nrows, ncols, bilin);
+ }
+ else { /* OBSERVATION POINTS INTERPOLATION */
if (ext == FALSE) {
G_debug(1, "Interpolation: (%d,%d): Sparse_Points...",
subregion_row, subregion_col);
@@ -632,11 +635,7 @@
dims.overlap, nsplx, nsply, npoints,
bilin, Cats, driver, mean,
table_name);
-
- G_free_matrix(obsVect);
- G_free_vector(parVect);
}
-
else { /* FLAG_EXT == TRUE */
int npoints_ext, *lineVect_ext = NULL;
double **obsVect_ext; /*, mean_ext = .0; */
@@ -670,53 +669,38 @@
mean, table_name);
G_free_matrix(obsVect_ext);
- G_free_vector(parVect);
G_free_ivector(lineVect_ext);
} /* END FLAG_EXT == TRUE */
- } /* END IF GRID == FLASE */
-
- else { /*GRID INTERPOLATION ==> INTERPOLATION INTO A RASTER */
- G_free_matrix(obsVect);
- flag_auxiliar = TRUE;
- G_debug(1, "Interpolation: (%d,%d): Regular_Points...",
- subregion_row, subregion_col);
- raster_matrix =
- P_Regular_Points(&elaboration_reg, general_box,
- overlap_box, raster_matrix, parVect,
- passoN, passoE, dims.overlap, mean,
- nsplx, nsply, nrows, ncols, bilin);
- G_free_vector(parVect);
- }
+ } /* END GRID == FALSE */
+ G_free_vector(parVect);
+ G_free_matrix(obsVect);
G_free_ivector(lineVect);
}
- else
- G_warning(_("No data within this subzone. "
+ else {
+ G_free(observ);
+ G_warning(_("No data within this subregion. "
"Consider changing the spline step."));
-
- subzone = (subregion_row - 1) * nsubregion_col + subregion_col;
- nsubzones = nsubregion_row * nsubregion_col;
- G_percent(subzone, nsubzones, 10);
-
+ }
} /*! END WHILE; last_column = TRUE */
} /*! END WHILE; last_row = TRUE */
- /* Writing into the output vector map the points from the overlapping zones */
- if (flag_auxiliar == TRUE) {
- if (grid == FALSE) {
- if (ext == FALSE)
- P_Aux_to_Vector(&In, &Out, driver, table_name);
- else
- P_Aux_to_Vector(&In_ext, &Out, driver, table_name);
+ G_verbose_message(_("Writing output..."));
+ /* Writing the output raster map */
+ if (grid == TRUE) {
+ P_Aux_to_Raster(raster_matrix, raster);
+ G_free_matrix(raster_matrix);
+ }
+ /* Writing to the output vector map the points from the overlapping zones */
+ else if (flag_auxiliar == TRUE) {
+ if (ext == FALSE)
+ P_Aux_to_Vector(&In, &Out, driver, table_name);
+ else
+ P_Aux_to_Vector(&In_ext, &Out, driver, table_name);
- /* Dropping auxiliar table */
- G_debug(1, "%s: Dropping <%s>", argv[0], table_name);
- if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
- G_fatal_error(_("Auxiliar table could not be dropped"));
- }
- else {
- P_Aux_to_Raster(raster_matrix, raster);
- G_free_matrix(raster_matrix);
- }
+ /* Dropping auxiliar table */
+ G_debug(1, "%s: Dropping <%s>", argv[0], table_name);
+ if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
+ G_fatal_error(_("Auxiliar table could not be dropped"));
}
db_close_database_shutdown_driver(driver);
@@ -732,7 +716,7 @@
/* set map title */
sprintf(title, "%s interpolation with Tykhonov regularization",
- type->answer);
+ type_opt->answer);
G_put_cell_title(out_map_opt->answer, title);
/* write map history */
G_short_history(out_map_opt->answer, "raster", &history);
More information about the grass-commit
mailing list