Added SmoothVectors tool

Added SmoothVectors tool

README.md

@@ -298,6 +298,7 @@ The library currently contains more than 370 tools, which are each grouped based
 - ***ReclassFromFile***: Reclassifies the values in a raster image using reclass ranges in a text file.
 - ***RelatedCircumscribingCircle***: Calculates the related circumscribing circle of vector polygons.
 - ***ShapeComplexityIndex***: Calculates overall polygon shape complexity or irregularity.
+- ***SmoothVectors***: Smooths a vector coverage of either a POLYLINE or POLYGON base ShapeType.
 - ***SumOverlay***: Calculates the sum for each grid cell from a group of raster images.
 - ***TINGridding***: Creates a raster grid based on a triangular irregular network (TIN) fitted to vector points.
 - ***VectorHexBinning***: Hex-bins a set of vector points.

readme.txt

@@ -60,24 +60,31 @@ Version 0.11.0 (XX-XX-2018)
 - The following tools were added to the project:
     AddPointCoordinatesToTable
     CentroidVector
-    CompactnessRatio and PerimeterAreaRatio
+    CompactnessRatio
+    ConstructVectorTIN
     ElongationRatio
     ExtendVectorLines
     HoleProportion
     LayerFootprint
-    LidarConstructVectorTIN and ConstructVectorTIN
-    LidarTINGridding and TINGridding
+    LidarConstructVectorTIN
+    LidarTINGridding
     LinesToPolygons
     Medoid
-    MinimumBoundingCircle and MinimumBoundingEnvelope
-    MultiPartToSinglePart and SinglePartToMultiPart
-    PolygonArea and PolygonPerimeter
+    MinimumBoundingCircle
+    MinimumBoundingEnvelope
+    MultiPartToSinglePart
+    PerimeterAreaRatio
+    PolygonArea
+    PolygonPerimeter
     RasterStreamsToVector
     RasterToVectorPoints
     RelatedCircumscribingCircle
     RemovePolygonHoles
     ShapeComplexityIndex
+    SinglePartToMultiPart
+    SmoothVectors
     SumOverlay
+    TINGridding
 
 - Added a minimum number of neighbours criteria in the neighbourhood search of the
   LidarGroundPointFilter tool. In this way, if the fixed-radius search yields fewer

src/algorithms/mod.rs

@@ -10,8 +10,8 @@ mod convex_hull;
 mod delaunay_triangulation;
 mod is_clockwise_order;
 mod minimum_bounding_box;
-mod point_in_poly;
 mod poly_area;
+mod poly_ops;
 mod poly_perimeter;
 mod smallest_enclosing_circle;
 
@@ -20,7 +20,7 @@ pub use self::convex_hull::convex_hull;
 pub use self::delaunay_triangulation::{triangulate, Triangulation};
 pub use self::is_clockwise_order::is_clockwise_order;
 pub use self::minimum_bounding_box::{minimum_bounding_box, MinimizationCriterion};
-pub use self::point_in_poly::{point_in_poly, poly_in_poly};
 pub use self::poly_area::polygon_area;
+pub use self::poly_ops::{point_in_poly, poly_in_poly, poly_is_convex, winding_number};
 pub use self::poly_perimeter::polygon_perimeter;
 pub use self::smallest_enclosing_circle::smallest_enclosing_circle;

src/algorithms/poly_ops.rs

@@ -0,0 +1,201 @@
+/* 
+This tool is part of the WhiteboxTools geospatial analysis library.
+Authors: Dr. John Lindsay
+Created: 30/08/2018
+Last Modified: 30/08/2018
+License: MIT
+*/
+
+use structures::Point2D;
+
+/// Tests if a point is Left|On|Right of an infinite line,
+/// based on http://geomalgorithms.com/a03-_inclusion.html.
+///
+/// Input:  three points p0, p1, and p2
+///
+/// Return: > 0 for p2 left of the line through p0 and p1
+///         = 0 for p2 on the line through p0 and p1
+///         < 0 for p2 right of the line through p0 and p1
+fn is_left(p0: &Point2D, p1: &Point2D, p2: &Point2D) -> f64 {
+    (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y)
+}
+
+/// Tests whether a point is within in a polygon using the winding number (wn).
+/// The point falls within the test polygon if the winding number is non-zero.
+/// Notice that points on the edge of the poly will be deemed outside.
+/// Input:   p = a point,
+///          poly[] = vertex points of a polygon v[n+1] with v[n]=v[0]
+pub fn point_in_poly(p: &Point2D, poly: &[Point2D]) -> bool {
+    // winding_number(&p, &poly) != 0i32
+    winding_number(&p, &poly) % 2 == 1i32
+}
+
+/// Calculates the Winding number (wn) test for a point in a polygon operation
+/// in order to determine whether the point is within the polygon. The
+/// point falls within the test polygon if the winding number is non-zero.
+///
+/// Input:   p = a point,
+///          poly[] = vertex points of a polygon poly[n+1] with poly[n]=poly[0]
+pub fn winding_number(p: &Point2D, poly: &[Point2D]) -> i32 {
+    if poly[0] != poly[poly.len() - 1] {
+        panic!("Error (from poly_ops::winding_num): point squence does not form a closed polygon.");
+    }
+    let mut wn = 0i32;
+    // loop through all edges of the polygon
+    for i in 0..poly.len() - 1 {
+        // edge from poly[i] to poly[i+1]
+        if poly[i].y <= p.y {
+            // start y <= p.y
+            if poly[i + 1].y > p.y {
+                // an upward crossing
+                if is_left(&poly[i], &poly[i + 1], &p) > 0f64 {
+                    // p left of edge
+                    wn += 1i32; // have a valid up intersect
+                }
+            }
+        } else {
+            // start y > p.y (no test needed)
+            if poly[i + 1].y <= p.y {
+                // a downward crossing
+                if is_left(&poly[i], &poly[i + 1], &p) < 0f64 {
+                    // p right of edge
+                    wn -= 1i32; // have a valid down intersect
+                }
+            }
+        }
+    }
+    wn
+}
+
+/// Tests whether one polygon is contained within another polygon. Notice that
+/// for polygons that are not contained within the test poly, failure occurs
+/// very quickly. In the case of disjoint (non-overlapping) polys, the function
+/// returns from the first tested vertex.
+pub fn poly_in_poly(contained_poly: &[Point2D], containing_poly: &[Point2D]) -> bool {
+    for p in contained_poly {
+        if !point_in_poly(p, containing_poly) {
+            return false;
+        }
+    }
+    true
+}
+
+/// Return true if the polygon is convex.
+pub fn poly_is_convex(poly: &[Point2D]) -> bool {
+    // For each set of three adjacent points A, B, C,
+    // find the cross product AB · BC. If the sign of
+    // all the cross products is the same, the angles
+    // are all positive or negative (depending on the
+    // order in which we visit them) so the polygon
+    // is convex.
+    let mut got_negative = false;
+    let mut got_positive = false;
+    let num_points = poly.len();
+    let (mut b, mut c): (usize, usize);
+    let mut cross_product: f64;
+    for a in 0..num_points {
+        b = (a + 1) % num_points;
+        c = (b + 1) % num_points;
+
+        cross_product = (poly[a].x - poly[b].x) * (poly[c].y - poly[b].y)
+            - (poly[a].y - poly[b].y) * (poly[c].x - poly[b].x);
+        if cross_product < 0f64 {
+            got_negative = true;
+        } else if cross_product > 0f64 {
+            got_positive = true;
+        }
+        if got_negative && got_positive {
+            return false;
+        }
+    }
+
+    // If we got this far, the polygon is convex.
+    return true;
+}
+
+// pub fn line_poly_intersections(line: &[Point2D], poly: &[Point2D]) -> bool {
+//     for p in line {
+//         if !point_in_poly(p, containing_poly) {
+//             return false;
+//         }
+//     }
+// }
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use structures::Point2D;
+    #[test]
+    fn test_point_in_poly() {
+        let poly = [
+            Point2D::new(0.0, 0.0),
+            Point2D::new(5.0, 0.0),
+            Point2D::new(5.0, 5.0),
+            Point2D::new(0.0, 0.0),
+        ];
+        // point inside rectangle
+        assert!(point_in_poly(&Point2D::new(2.0, 2.0), &poly));
+        // point outside rectangle
+        assert_eq!(point_in_poly(&Point2D::new(12.0, 12.0), &poly), false);
+    }
+
+    #[test]
+    fn test_winding_number() {
+        let poly = [
+            Point2D::new(0.0, 0.0),
+            Point2D::new(5.0, 0.0),
+            Point2D::new(5.0, 5.0),
+            Point2D::new(0.0, 0.0),
+        ];
+        // point on rectangle
+        assert_eq!(winding_number(&Point2D::new(5.0, 2.0), &poly), 0i32);
+        assert_eq!(winding_number(&Point2D::new(4.0, 2.0), &poly), 1i32);
+        assert_eq!(winding_number(&Point2D::new(6.0, 2.0), &poly), 0i32);
+    }
+
+    #[test]
+    fn test_poly_in_poly() {
+        let poly1 = [
+            Point2D::new(0.0, 0.0),
+            Point2D::new(5.0, 0.0),
+            Point2D::new(5.0, 5.0),
+            Point2D::new(0.0, 0.0),
+        ];
+
+        let poly2 = [
+            Point2D::new(-1.0, -1.0),
+            Point2D::new(6.0, -1.0),
+            Point2D::new(6.0, 6.0),
+            Point2D::new(-1.0, -1.0),
+        ];
+
+        assert!(poly_in_poly(&poly1, &poly2));
+
+        assert_eq!(poly_in_poly(&poly2, &poly1), false);
+    }
+
+    #[test]
+    fn test_poly_is_convex() {
+        let poly = [
+            Point2D::new(0.0, 0.0),
+            Point2D::new(5.0, 0.0),
+            Point2D::new(5.0, 5.0),
+            Point2D::new(0.0, 5.0),
+            Point2D::new(0.0, 0.0),
+        ];
+        // point on rectangle
+        assert!(poly_is_convex(&poly));
+
+        let poly = [
+            Point2D::new(0.0, 0.0),
+            Point2D::new(5.0, 0.0),
+            Point2D::new(5.0, 5.0),
+            Point2D::new(2.5, 3.0),
+            Point2D::new(0.0, 5.0),
+            Point2D::new(0.0, 0.0),
+        ];
+        // point on rectangle
+        assert_eq!(poly_is_convex(&poly), false);
+    }
+
+}

src/tools/gis_analysis/mod.rs

@@ -56,6 +56,7 @@ mod reclass_equal_interval;
 mod reclass_from_file;
 mod related_circumscribing_circle;
 mod shape_complexity_index;
+mod smooth_vectors;
 mod sum_overlay;
 mod tin_gridding;
 mod vector_hex_bin;
@@ -120,6 +121,7 @@ pub use self::reclass_equal_interval::ReclassEqualInterval;
 pub use self::reclass_from_file::ReclassFromFile;
 pub use self::related_circumscribing_circle::RelatedCircumscribingCircle;
 pub use self::shape_complexity_index::ShapeComplexityIndex;
+pub use self::smooth_vectors::SmoothVectors;
 pub use self::sum_overlay::SumOverlay;
 pub use self::tin_gridding::TINGridding;
 pub use self::vector_hex_bin::VectorHexBinning;