computing nominal point density for tessellated polygon in filters.he…

…xbin

include/pdal/Geometry.hpp

@@ -116,6 +116,21 @@ static std::string smoothPolygon(const std::string& wkt, double tolerance)
     finish();
     return outWkt;
 }
+
+static double computeArea(const std::string& wkt)
+{
+    GEOSContextHandle_t env = init();
+
+    GEOSGeometry *geom = GEOSGeomFromWKT_r(env, wkt.c_str());
+    if (!geom)
+        return 0.0;
+
+    double output(0.0);
+    int er = GEOSArea_r(env, geom, &output);
+    GEOSGeom_destroy_r(env, geom);
+    finish();
+    return output;
+}
 #else
 
 static std::string smoothPolygon(const std::string& wkt, double tolerance)

plugins/hexbin/filters/HexBin.cpp

@@ -83,6 +83,7 @@ void HexBin::filter(PointView& view)
         double x = view.getFieldAs<double>(pdal::Dimension::Id::X, idx);
         double y = view.getFieldAs<double>(pdal::Dimension::Id::Y, idx);
         m_grid->addPoint(x, y);
+        m_count++;
     }
 }
 
@@ -151,6 +152,16 @@ void HexBin::done(PointTableRef table)
     m_metadata.add("smoothed_boundary",
         Geometry::smoothPolygon(polygon.str(), tolerance),
         "Smoothed boundary MULTIPOLYGON of domain");
+    double area = Geometry::computeArea(polygon.str());
+
+    double density = (double) m_count / area ;
+    m_metadata.add("density",
+            density,
+        "Number of points per square unit");
+    m_metadata.add("area",
+            area,
+        "Area in square units of tessellated polygon");
+
 }
 
 } // namespace pdal

plugins/hexbin/filters/HexBin.hpp

@@ -62,6 +62,7 @@ class PDAL_DLL HexBin : public Filter
     int32_t m_density;
     double m_edgeLength;
     bool m_outputTesselation;
+    point_count_t m_count;
 
     virtual void processOptions(const Options& options);
     virtual void ready(PointTableRef table);