Merge remote-tracking branch 'origin/master' into issue-3525

cmake/win32_compiler_options.cmake

@@ -6,7 +6,8 @@ function(pdal_target_compile_settings target)
     set_property(TARGET ${target} PROPERTY CXX_STANDARD 11)
     set_property(TARGET ${target} PROPERTY CXX_STANDARD_REQUIRED TRUE)
     target_compile_definitions(${target} PRIVATE
-        -DWIN32_LEAN_AND_MEAN)
+        -DWIN32_LEAN_AND_MEAN
+        -D_USE_MATH_DEFINES)
     if (MSVC)
         # check for MSVC 8+
         if (NOT (MSVC_VERSION VERSION_LESS 1400))

doc/stages/readers.ept.rst

@@ -103,12 +103,19 @@ origin
 polygon
   The clipping polygon, expressed in a well-known text string,
   eg: "POLYGON((0 0, 5000 10000, 10000 0, 0 0))".  This option can be
-  specified more than once by placing values in an array.
+  specified more than once by placing values in an array, in which case all of
+  them will be unioned together, acting as a single multipolygon.
 
 .. note::
 
     When using ``pdal info --summary``, using the ``polygon`` option will cause the resulting bounds to be clipped to the maximal extents of all provided polygons, and the resulting number of points to be an upper bound for this polygon selection.
 
+.. note::
+
+    When a ``bounds`` option is specified alongside the ``polygon`` option, only
+    the points that fall within *both* the bounds and the polygon(s) will be
+    returned.
+
 ogr
   A JSON object representing an OGR query to fetch polygons to use for filtering. The polygons
   fetched from the query are treated exactly like those specified in the ``polygon`` option.

filters/GreedyProjection.hpp

@@ -39,10 +39,7 @@
 
 #pragma once
 
-// This is for M_PI on Windows.
-#define _USE_MATH_DEFINES
-#include <math.h>
-
+#include <cmath>
 #include <fstream>
 #include <iostream>
 

io/EptReader.cpp

@@ -556,7 +556,7 @@ bool EptReader::passesSpatialFilter(const BOX3D& tileBounds) const
     auto boxOverlaps = [this, &reproject, &tileBounds]() -> bool
     {
         if (!m_p->bounds.box.valid())
-            return false;
+            return true;
 
         // If the reprojected source bounds doesn't overlap our query bounds, we're done.
         return reproject(tileBounds, m_p->bounds.xform).overlaps(m_p->bounds.box);
@@ -566,6 +566,9 @@ bool EptReader::passesSpatialFilter(const BOX3D& tileBounds) const
     // we can skip
     auto polysOverlap = [this, &reproject, &tileBounds]() -> bool
     {
+        if (m_p->polys.empty())
+            return true;
+
         for (auto& ps : m_p->polys)
             if (!ps.poly.disjoint(reproject(tileBounds, ps.xform)))
                 return true;
@@ -580,7 +583,7 @@ bool EptReader::passesSpatialFilter(const BOX3D& tileBounds) const
     // at the same time, it seems to get corrupted. There may be other instances
     // that need to be locked.
     std::lock_guard<std::mutex> lock(m_p->mutex);
-    return boxOverlaps() || polysOverlap();
+    return boxOverlaps() && polysOverlap();
 }
 
 
@@ -675,13 +678,16 @@ bool EptReader::processPoint(PointRef& dst, const TileContents& tile)
     auto passesBoundsFilter = [this](double x, double y, double z)
     {
         if (!m_p->bounds.box.valid())
-            return false;
+            return true;
         m_p->bounds.xform.transform(x, y, z);
         return m_p->bounds.box.contains(x, y, z);
     };
 
     auto passesPolyFilter = [this](double xo, double yo, double zo)
     {
+        if (m_p->polys.empty())
+            return true;
+
         for (PolyXform& ps : m_p->polys)
         {
             double x = xo;
@@ -701,7 +707,7 @@ bool EptReader::processPoint(PointRef& dst, const TileContents& tile)
 
     // If there is a spatial filter, make sure it passes.
     if (hasSpatialFilter())
-        if (!passesBoundsFilter(x, y, z) && !passesPolyFilter(x, y, z))
+        if (!passesBoundsFilter(x, y, z) || !passesPolyFilter(x, y, z))
             return false;
 
     for (auto& el : m_p->info->dims())

io/OptechReader.cpp

@@ -32,7 +32,6 @@
 * OF SUCH DAMAGE.
 ****************************************************************************/
 
-#define _USE_MATH_DEFINES
 #include "OptechReader.hpp"
 
 #include <cmath>

io/SbetReader.cpp

@@ -38,8 +38,7 @@
 #include <pdal/PointRef.hpp>
 #include <pdal/util/FileUtils.hpp>
 
-#define _USE_MATH_DEFINES
-#include <math.h>
+#include <cmath>
 
 namespace pdal
 {

io/SbetWriter.cpp

@@ -37,8 +37,7 @@
 #include <pdal/PointView.hpp>
 #include <pdal/util/ProgramArgs.hpp>
 
-#define _USE_MATH_DEFINES
-#include <math.h>
+#include <cmath>
 
 namespace pdal
 {

test/unit/GeoreferenceTest.cpp

@@ -32,13 +32,11 @@
 * OF SUCH DAMAGE.
 ****************************************************************************/
 
-#define _USE_MATH_DEFINES
 #include <pdal/pdal_test_main.hpp>
 #include <pdal/util/Georeference.hpp>
 
 #include <cmath>
 
-
 namespace pdal
 {
 namespace georeference

test/unit/io/EptReaderTest.cpp

@@ -629,6 +629,86 @@ TEST(EptReaderTest, boundedCrop)
     EXPECT_EQ(sourceNp, 47u);
 }
 
+TEST(EptReaderTest, polygonAndBoundsCrop)
+{
+    std::string wkt = FileUtils::readFileIntoString(
+        Support::datapath("autzen/autzen-selection.wkt"));
+
+    // This box is approximately the bounding box of the WKT above, with the
+    // eastmost 25% of the bounds omitted.  So this should shrink our query
+    // results from the "boundedCrop" test above since we are further limiting
+    // our spatial selection.
+    std::string boxstring = "([636577.1, 637297.4225], [850571.42, 851489.34])";
+    BOX2D box;
+    Utils::fromString(boxstring, box);
+
+    // First we'll query the EptReader for these bounds.
+    EptReader reader;
+    {
+        Options options;
+        options.add("filename", eptAutzenPath);
+        Option polygon("polygon", wkt + "/ EPSG:3644");
+        options.add(polygon);
+        Option bounds("bounds", boxstring);
+        options.add(bounds);
+        reader.setOptions(options);
+    }
+
+    PointTable eptTable;
+    reader.prepare(eptTable);
+
+    uint64_t eptNp(0);
+    for (const PointViewPtr& view : reader.execute(eptTable))
+    {
+        eptNp += view->size();
+    }
+
+    // Now we'll check the result against a crop filter of the source file with
+    // the same bounds.
+    LasReader source;
+    {
+        Options options;
+        options.add("filename", Support::datapath("las/1.2-with-color.las"));
+        source.setOptions(options);
+    }
+    CropFilter boundsCrop;
+    {
+        Options options;
+        Option bounds("bounds", boxstring);
+        options.add(bounds);
+        boundsCrop.setOptions(options);
+        boundsCrop.setInput(source);
+    }
+    CropFilter polygonCrop;
+    {
+        Options options;
+        Option polygon("polygon", wkt + "/ EPSG:3644");
+        options.add(polygon);
+        polygonCrop.setOptions(options);
+        polygonCrop.setInput(boundsCrop);
+    }
+    PointTable sourceTable;
+    polygonCrop.prepare(sourceTable);
+    uint64_t sourceNp(0);
+
+    BOX2D got;
+    for (const PointViewPtr& view : polygonCrop.execute(sourceTable))
+    {
+        sourceNp += view->size();
+        for (std::size_t i = 0; i < view->size(); ++i) {
+            EXPECT_TRUE(
+                box.contains(
+                    view->getFieldAs<double>(pdal::Dimension::Id::X, i),
+                    view->getFieldAs<double>(pdal::Dimension::Id::Y, i)));
+        }
+    }
+
+    EXPECT_EQ(eptNp, sourceNp);
+    EXPECT_EQ(eptNp, 38u);
+    EXPECT_EQ(sourceNp, 38u);
+}
+
+
 TEST(EptReaderTest, boundedCropReprojection)
 {
     std::string selection = FileUtils::readFileIntoString(