Smooth hexbin-generated boundaries.

include/pdal/Geometry.hpp

@@ -0,0 +1,131 @@
+/******************************************************************************
+* Copyright (c) 2015, Hobu Inc. (info@hobu.co)
+*
+* All rights reserved.
+*
+* Redistribution and use in source and binary forms, with or without
+* modification, are permitted provided that the following
+* conditions are met:
+*
+*     * Redistributions of source code must retain the above copyright
+*       notice, this list of conditions and the following disclaimer.
+*     * Redistributions in binary form must reproduce the above copyright
+*       notice, this list of conditions and the following disclaimer in
+*       the documentation and/or other materials provided
+*       with the distribution.
+*     * Neither the name of Hobu, Inc. nor the
+*       names of its contributors may be used to endorse or promote
+*       products derived from this software without specific prior
+*       written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
+* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
+* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+* OF SUCH DAMAGE.
+****************************************************************************/
+
+#include <pdal/pdal_types.hpp>
+
+#include <cstdarg>
+
+#ifdef PDAL_HAVE_GEOS
+#include <geos_c.h>
+#endif
+
+#pragma once
+
+namespace pdal
+{
+
+#ifdef PDAL_HAVE_GEOS
+namespace
+{
+
+static GEOSContextHandle_t s_environment(NULL);
+static int s_contextCnt(0);
+
+static void GEOSErrorHandler(const char *fmt, ...)
+{
+    va_list args;
+
+    va_start(args, fmt);
+    char buf[1024];
+
+    vsnprintf(buf, sizeof(buf), fmt, args);
+    std::cout << "GEOS error: " << buf << std::endl;
+
+    va_end(args);
+}
+
+static void GEOSWarningHandler(const char *fmt, ...)
+{
+    va_list args;
+
+    char buf[1024];
+    vsnprintf(buf, sizeof(buf), fmt, args);
+    std::cout << "GEOS warning: " << buf << std::endl;
+
+    va_end(args);
+}
+
+} // Unnamed namespace
+
+namespace Geometry
+{
+
+static GEOSContextHandle_t init()
+{
+    if (s_contextCnt == 0)
+        s_environment = initGEOS_r(GEOSWarningHandler, GEOSErrorHandler);
+    s_contextCnt++;
+    return s_environment;
+}
+
+static void finish()
+{
+    s_contextCnt--;
+    if (s_contextCnt == 0)
+        finishGEOS_r(s_environment);
+}
+
+static std::string smoothPolygon(const std::string& wkt, double tolerance)
+{
+    GEOSContextHandle_t env = init();
+
+    GEOSGeometry *geom = GEOSGeomFromWKT_r(env, wkt.c_str());
+    if (!geom)
+        return "";
+
+    GEOSGeometry *smoothed = GEOSTopologyPreserveSimplify_r(env, geom,
+        tolerance);
+    if (!smoothed)
+        return "";
+    char *smoothedWkt = GEOSGeomToWKT_r(env, smoothed);
+    std::string outWkt(smoothedWkt);
+    GEOSFree_r(env, smoothedWkt);
+    GEOSGeom_destroy_r(env, geom);
+    GEOSGeom_destroy_r(env, smoothed);
+    finish();
+    return outWkt;
+}
+#else
+
+static std::string smoothPolygon(const std::string& wkt, double tolerance)
+{
+    throw pdal_error("Can't call smoothPolygon.  PDAL not built with GEOS.");
+}
+
+#endif
+
+} // namespace Geometry
+
+} // namespace pdal
+

kernels/info/InfoKernel.cpp

@@ -397,6 +397,7 @@ void InfoKernel::dump(MetadataNode& root)
         PointViewSet viewSet = m_manager->views();
         assert(viewSet.size() == 1);
         root.add(m_hexbinStage->getMetadata().clone("boundary"));
+        root.add(m_hexbinStage->getMetadata().clone("smoothed_boundary"));
     }
 }
 

kernels/tindex/TIndexKernel.cpp

@@ -86,6 +86,7 @@ TIndexKernel::TIndexKernel()
     , m_dataset(NULL)
     , m_layer(NULL)
     , m_fastBoundary(false)
+    , m_smoothBoundary(false)
 
 {
     m_log.setLeader("pdal tindex");
@@ -104,9 +105,12 @@ void TIndexKernel::addSwitches()
             "OGR-readable/writeable tile index output")
         ("filespec", po::value<std::string>(&m_filespec),
             "Build: Pattern of files to index. Merge: Output filename")
+        ("smooth_boundary", po::value<bool>(&m_smoothBoundary)->
+            zero_tokens()->implicit_value(false),
+            "use smoothed version of full boundary")
         ("fast_boundary", po::value<bool>(&m_fastBoundary)->
             zero_tokens()->implicit_value(true),
-            "use extend instead of exact boundary")
+            "use extent instead of exact boundary")
         ("lyr_name", po::value<std::string>(&m_layerName),
             "OGR layer name to write into datasource")
         ("tindex_name", po::value<std::string>(&m_tileIndexColumnName)->
@@ -180,6 +184,9 @@ void TIndexKernel::validateSwitches()
             throw pdal_error("--bounds option not supported when building "
                 "index.");
     }
+    if (m_smoothBoundary && m_fastBoundary)
+        throw pdal_error("Can't request both smooth_boundary and "
+            "fast_boundary.");
 }
 
 
@@ -415,7 +422,6 @@ void TIndexKernel::mergeFile()
             Options reproOptions;
             reproOptions.add("out_srs", m_tgtSrsString);
             reproOptions.add("in_srs", f.m_srs);
-    std::cerr << "Repro = " << m_tgtSrsString << "/" << f.m_srs << "!\n";
             repro->setOptions(reproOptions);
             premerge = repro;
         }
@@ -589,8 +595,12 @@ TIndexKernel::FileInfo TIndexKernel::getFileInfo(KernelFactory& factory,
         hexer->prepare(table);
         hexer->execute(table);
 
-        fileInfo.m_boundary =
-            table.metadata().findChild("filters.hexbin:boundary").value();
+        MetadataNode m = table.metadata();
+        m = m_smoothBoundary ?
+            m.findChild("filters.hexbin:smooth_boundary") :
+            m.findChild("filters.hexbin:boundary");
+        fileInfo.m_boundary = m.value();
+
         if (!table.spatialRef().empty())
             fileInfo.m_srs = table.spatialRef().getWKT();
     }

kernels/tindex/TIndexKernel.hpp

@@ -113,6 +113,7 @@ class PDAL_DLL TIndexKernel : public Kernel
     std::string m_tgtSrsString;
     std::string m_assignSrsString;
     bool m_fastBoundary;
+    bool m_smoothBoundary;
 };
 
 } // namespace pdal

plugins/hexbin/filters/HexBin.cpp

@@ -35,6 +35,7 @@
 #include "HexBin.hpp"
 
 #include <hexer/HexIter.hpp>
+#include <pdal/Geometry.hpp>
 #include <pdal/StageFactory.hpp>
 
 using namespace hexer;
@@ -141,6 +142,15 @@ void HexBin::done(PointTableRef table)
     m_grid->toWKT(polygon);
     m_metadata.add("boundary", polygon.str(),
         "Boundary MULTIPOLYGON of domain");
+
+    const double SQRT_3(1.732050808);
+    // This is just a little longer than the prependicular distance of the
+    // "humps."
+//    double tolerance = 1.1 * m_grid->height() / (2 * SQRT_3);
+    double tolerance = 1.1 * m_grid->height() / 2;
+    m_metadata.add("smoothed_boundary",
+        Geometry::smoothPolygon(polygon.str(), tolerance),
+        "Smoothed boundary MULTIPOLYGON of domain");
 }
 
 } // namespace pdal

plugins/pcl/dartsample/dart_sample.hpp

@@ -80,7 +80,7 @@ pcl::DartSample<PointT>::applyFilter (std::vector<int> &indices)
 
   // Create octree, seeded with the first index
   pcl::octree::OctreePointCloudSearch<PointT> tree (radius_ / std::sqrt (3));
-  pcl::PointCloud<PointT>::Ptr cloud_t (new pcl::PointCloud<PointT>);
+  typename pcl::PointCloud<PointT>::Ptr cloud_t (new pcl::PointCloud<PointT>);
   tree.setInputCloud (cloud_t);
   tree.addPointToCloud (input_->points[shuffled_indices[0]], cloud_t);
 

plugins/pcl/filters/HeightFilter.cpp

@@ -158,7 +158,7 @@ void HeightFilter::filter(PointView& view)
     ground_tree.reset(new pcl::search::KdTree<pcl::PointXYZ> (false));
     ground_tree->setInputCloud(cloud_ground_projected);
 
-    for (int i = 0; i < cloud_nonground_projected->size(); ++i)
+    for (size_t i = 0; i < cloud_nonground_projected->size(); ++i)
     {
         pcl::PointXYZ nonground_query = cloud_nonground_projected->points[i];
         std::vector<int> neighbors(1);

src/CMakeLists.txt

@@ -43,6 +43,7 @@ set(PDAL_BASE_HPP
   "${PDAL_HEADERS_DIR}/Filter.hpp"
   "${PDAL_HEADERS_DIR}/FlexWriter.hpp"
   "${PDAL_HEADERS_DIR}/GDALUtils.hpp"
+  "${PDAL_HEADERS_DIR}/Geometry.hpp"
   "${PDAL_HEADERS_DIR}/GlobalEnvironment.hpp"
   "${PDAL_HEADERS_DIR}/gitsha.h"
   "${PDAL_HEADERS_DIR}/KDIndex.hpp"