v.outlier: add option to filter out only positive or negative outliers

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

vector/v.outlier/main.c

@@ -42,6 +42,7 @@ int main(int argc, char *argv[])
     char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
     int last_row, last_column, flag_auxiliar = FALSE;
+    int filter_mode;
 
     int *lineVect;
     double *TN, *Q, *parVect;	/* Interpolating and least-square vectors */
@@ -50,7 +51,7 @@ int main(int argc, char *argv[])
     /* Structs' declarations */
     struct Map_info In, Out, Outlier, Qgis;
     struct Option *in_opt, *out_opt, *outlier_opt, *qgis_opt, *stepE_opt,
-	*stepN_opt, *lambda_f_opt, *Thres_O_opt;
+	*stepN_opt, *lambda_f_opt, *Thres_O_opt, *filter_opt;
     struct Flag *spline_step_flag;
     struct GModule *module;
 
@@ -125,6 +126,14 @@ int main(int argc, char *argv[])
     Thres_O_opt->description = _("Threshold for the outliers");
     Thres_O_opt->answer = "50";
 
+    filter_opt = G_define_option();
+    filter_opt->key = "filter";
+    filter_opt->type = TYPE_STRING;
+    filter_opt->required = NO;
+    filter_opt->description = _("Filtering option");
+    filter_opt->options = "both,positive,negative";
+    filter_opt->answer = "both";
+
     /* Parsing */
     G_gisinit(argv[0]);
 
@@ -142,6 +151,13 @@ int main(int argc, char *argv[])
     lambda = atof(lambda_f_opt->answer);
     Thres_Outlier = atof(Thres_O_opt->answer);
 
+    filter_mode = 0;
+    if (strcmp(filter_opt->answer, "positive") == 0)
+	filter_mode = 1;
+    else if (strcmp(filter_opt->answer, "negative") == 0)
+	filter_mode = -1;
+    P_set_outlier_fn(filter_mode);
+
     flag_auxiliar = FALSE;
 
     /* Checking vector names */
@@ -402,11 +418,13 @@ int main(int argc, char *argv[])
 		if (qgis_opt->answer)
 		    P_Outlier(&Out, &Outlier, &Qgis, elaboration_reg,
 			      general_box, overlap_box, obsVect, parVect,
-			      mean, dims.overlap, lineVect, npoints, driver, table_name);
+			      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, table_name);
+			      mean, dims.overlap, lineVect, npoints,
+			      driver, table_name);
 
 
 		G_free_vector(parVect);

vector/v.outlier/outlier.c

@@ -4,6 +4,20 @@
 #include <math.h>
 #include "outlier.h"
 
+typedef int (*outlier_fn)(double);
+
+static outlier_fn is_outlier;
+
+void P_set_outlier_fn(int filter_mode)
+{
+    if (filter_mode < 0)
+	is_outlier = P_is_outlier_n;
+    else if (filter_mode > 0)
+	is_outlier = P_is_outlier_p;
+    else
+	is_outlier = P_is_outlier;
+}
+
 extern double Thres_Outlier;
 
 void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
@@ -20,6 +34,7 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
     struct line_pnts *point;
     struct line_cats *categories;
 
+
     point = Vect_new_line_struct();
     categories = Vect_new_cats_struct();
 
@@ -46,7 +61,7 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 
 		residual = *point->z - interpolation;
 
-		if (FALSE == P_is_outlier(residual)) {
+		if (FALSE == is_outlier(residual)) {
 		    Vect_write_line(Out, GV_POINT, point, categories);
 		    Vect_cat_set(categories, 1, (int)*point->z);
 		    if (Qgis)
@@ -109,7 +124,7 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 
 			residual = *point->z - interpolation;
 
-			if (FALSE == P_is_outlier(residual)) {
+			if (FALSE == is_outlier(residual)) {
 			    Vect_write_line(Out, GV_POINT, point, categories);
 			    Vect_cat_set(categories, 1, (int)*point->z);
 			    if (Qgis)
@@ -147,7 +162,7 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 
 			residual = *point->z - interpolation;
 
-			if (FALSE == P_is_outlier(residual)) {
+			if (FALSE == is_outlier(residual)) {
 			    Vect_write_line(Out, GV_POINT, point, categories);
 			    Vect_cat_set(categories, 1, (int)*point->z);
 			    if (Qgis)
@@ -171,7 +186,7 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 
 			residual = *point->z - interpolation;
 
-			if (FALSE == P_is_outlier(residual)) {
+			if (FALSE == is_outlier(residual)) {
 			    Vect_write_line(Out, GV_POINT, point, categories);
 			    Vect_cat_set(categories, 1, (int)*point->z);
 			    if (Qgis)
@@ -287,6 +302,22 @@ int P_is_outlier(double pippo)
     return TRUE;
 }
 
+int P_is_outlier_p(double pippo)
+{
+    if (pippo < Thres_Outlier)
+	return FALSE;
+
+    return TRUE;
+}
+
+int P_is_outlier_n(double pippo)
+{
+    if (pippo > Thres_Outlier)
+	return FALSE;
+
+    return TRUE;
+}
+
 /*! DEFINITION OF THE SUBZONES 
 
   5: inside Overlap region

vector/v.outlier/outlier.h

@@ -15,7 +15,7 @@ P_Outlier(struct Map_info *, /**/
 	  struct bound_box, /**/
 	  double **, /**/
 	  double *, /**/
-	  double, /**/ double, /**/ int *, /**/ int, /**/ dbDriver *,  /**/ char *);
+	  double, /**/ double, /**/ int *, /**/ int, /**/ dbDriver *, /**/ char *);
 
 int Insert_Outlier(double, /**/ int, /**/ dbDriver *, /**/ char *);
 
@@ -23,4 +23,7 @@ int UpDate_Outlier(double, /**/ int, /**/ dbDriver *, /**/ char *);
 
 int Select_Outlier(double *, /**/ int, /**/ dbDriver *, /**/ char *);
 
+void P_set_outlier_fn(int);
 int P_is_outlier(double);
+int P_is_outlier_p(double);
+int P_is_outlier_n(double);

vector/v.outlier/v.outlier.html

@@ -8,6 +8,15 @@ an absolute value more than the given threshold from a fixed value,
 reckoned from its surroundings by the interpolation, are considered as
 an outlier, and hence are removed.
 
+<p>
+The <em>filter</em> option specifies if all outliers will be removed
+(default), or only positive or only negative outliers. Filtering out
+only positive outliers can be useful to filter out vegetation returns
+(e.g. from forest canopies) from LIDAR point clouds, in order to
+extract Digital Terrain Models. Filtering out only negative outliers
+can be useful to estimate vegetation height.
+
+<p>
 There is a flag to create a vector that can be visualizated by
 qgis. That means that topology is build and the z coordinate is
 considered as a category.