fix #2904 - support cropping with 3D bounds

filters/CropFilter.cpp

@@ -121,7 +121,7 @@ void CropFilter::initialize()
 
     m_boxes.clear();
     for (auto& bound : m_args->m_bounds)
-        m_boxes.push_back(bound.to2d());
+        m_boxes.push_back(bound);
 
     m_distance2 = m_args->m_distance * m_args->m_distance;
 }
@@ -151,8 +151,17 @@ bool CropFilter::processOne(PointRef& point)
                 return true;
 
     for (auto& box : m_boxes)
-        if (crop(point, box))
-            return true;
+        if (box.is3d())
+        {
+            if (crop(point, box.to3d()))
+                return true;
+        }
+        else
+        {
+            if (crop(point, box.to2d()))
+                return true;
+        }
+
 
     for (auto& center: m_args->m_centers)
         if (crop(point, center))
@@ -247,6 +256,15 @@ PointViewSet CropFilter::run(PointViewPtr view)
     return viewSet;
 }
 
+bool CropFilter::crop(const PointRef& point, const BOX3D& box)
+{
+    double x = point.getFieldAs<double>(Dimension::Id::X);
+    double y = point.getFieldAs<double>(Dimension::Id::Y);
+    double z = point.getFieldAs<double>(Dimension::Id::Z);
+
+    // Return true if we're keeping a point.
+    return (m_args->m_cropOutside != box.contains(x, y, z));
+}
 
 bool CropFilter::crop(const PointRef& point, const BOX2D& box)
 {
@@ -257,6 +275,29 @@ bool CropFilter::crop(const PointRef& point, const BOX2D& box)
     return (m_args->m_cropOutside != box.contains(x, y));
 }
 
+void CropFilter::crop(const Bounds& box, PointView& input, PointView& output)
+{
+    PointRef point = input.point(0);
+    bool is3d = box.is3d();
+    for (PointId idx = 0; idx < input.size(); ++idx)
+    {
+        point.setPointId(idx);
+        if (is3d)
+        {
+            if (crop(point, box.to3d()))
+            {
+                output.appendPoint(input, idx);
+            }
+        } else 
+        {
+            if (crop(point, box.to2d()))
+            {
+                output.appendPoint(input, idx);
+            }
+        }
+    }
+}
+
 
 void CropFilter::crop(const BOX2D& box, PointView& input, PointView& output)
 {

filters/CropFilter.hpp

@@ -76,7 +76,7 @@ class PDAL_DLL CropFilter : public Filter, public Streamable
     std::unique_ptr<CropArgs> m_args;
     double m_distance2;
     std::vector<ViewGeom> m_geoms;
-    std::vector<BOX2D> m_boxes;
+    std::vector<Bounds> m_boxes;
 
     void addArgs(ProgramArgs& args);
     virtual void initialize();
@@ -86,7 +86,9 @@ class PDAL_DLL CropFilter : public Filter, public Streamable
     virtual bool processOne(PointRef& point);
     virtual PointViewSet run(PointViewPtr view);
     bool crop(const PointRef& point, const BOX2D& box);
+    bool crop(const PointRef& point, const BOX3D& box);
     void crop(const BOX2D& box, PointView& input, PointView& output);
+    void crop(const Bounds& box, PointView& input, PointView& output);
     bool crop(const PointRef& point, GridPnp& g);
     void crop(const ViewGeom& g, PointView& input, PointView& output);
     bool crop(const PointRef& point, const filter::Point& center);

pdal/GDALUtils.cpp

@@ -195,6 +195,18 @@ bool reprojectBounds(BOX3D& box, const std::string& srcSrs,
 }
 
 
+bool reprojectBounds(Bounds& box, const std::string& srcSrs,
+    const std::string& dstSrs)
+{
+    SrsTransform transform(srcSrs, dstSrs);
+
+    BOX3D b = box.to3d();
+    bool ok = transform.transform(b.minx, b.miny, b.minz);
+    if (ok)
+        ok = transform.transform(b.maxx, b.maxy, b.maxz);
+    return ok;
+}
+
 /**
   Reproject a bounds box from a source projection to a destination.
   \param box  2D Bounds box to be reprojected in-place.

pdal/GDALUtils.hpp

@@ -69,6 +69,8 @@ using ITER_VAL = typename std::iterator_traits<ITER>::value_type;
 
 PDAL_DLL void registerDrivers();
 PDAL_DLL void unregisterDrivers();
+PDAL_DLL bool reprojectBounds(Bounds& box, const std::string& srcSrs,
+    const std::string& dstSrs);
 PDAL_DLL bool reprojectBounds(BOX3D& box, const std::string& srcSrs,
     const std::string& dstSrs);
 PDAL_DLL bool reprojectBounds(BOX2D& box, const std::string& srcSrs,

test/unit/filters/CropFilterTest.cpp

@@ -116,6 +116,66 @@ TEST(CropFilterTest, test_crop)
     EXPECT_NEAR(avgZ, 500.00000, delZ);
 }
 
+TEST(CropFilterTest, test_crop_3d)
+{
+    BOX3D srcBounds(0.0, 0.0, 0.0, 10.0, 100.0, 1000.0);
+    Options opts;
+    opts.add("bounds", srcBounds);
+    opts.add("count", 1000);
+    opts.add("mode", "ramp");
+    FauxReader reader;
+    reader.setOptions(opts);
+
+    // crop the window to 20% the size in each dimension
+    BOX3D dstBounds(2.0, 20.0, 200.0, 4.0, 40.0, 400.0);
+    Options cropOpts;
+    cropOpts.add("bounds", dstBounds);
+
+    CropFilter filter;
+    filter.setOptions(cropOpts);
+    filter.setInput(reader);
+
+    Options statOpts;
+
+    StatsFilter stats;
+    stats.setOptions(statOpts);
+    stats.setInput(filter);
+
+    PointTable table;
+    stats.prepare(table);
+    PointViewSet viewSet = stats.execute(table);
+    EXPECT_EQ(viewSet.size(), 1u);
+    PointViewPtr buf = *viewSet.begin();
+
+    const stats::Summary& statsX = stats.getStats(Dimension::Id::X);
+    const stats::Summary& statsY = stats.getStats(Dimension::Id::Y);
+    const stats::Summary& statsZ = stats.getStats(Dimension::Id::Z);
+    EXPECT_EQ(buf->size(), 200u);
+
+    const double minX = statsX.minimum();
+    const double minY = statsY.minimum();
+    const double minZ = statsZ.minimum();
+    const double maxX = statsX.maximum();
+    const double maxY = statsY.maximum();
+    const double maxZ = statsZ.maximum();
+    const double avgX = statsX.average();
+    const double avgY = statsY.average();
+    const double avgZ = statsZ.average();
+
+    const double delX = 10.0 / 999.0 * 100.0;
+    const double delY = 100.0 / 999.0 * 100.0;
+    const double delZ = 1000.0 / 999.0 * 100.0;
+
+    EXPECT_NEAR(minX, 2.0, delX);
+    EXPECT_NEAR(minY, 20.0, delY);
+    EXPECT_NEAR(minZ, 200.0, delZ);
+    EXPECT_NEAR(maxX, 4.0, delX);
+    EXPECT_NEAR(maxY, 40.0, delY);
+    EXPECT_NEAR(maxZ, 400.0, delZ);
+    EXPECT_NEAR(avgX, 3.00000, delX);
+    EXPECT_NEAR(avgY, 30.00000, delY);
+    EXPECT_NEAR(avgZ, 300.00000, delZ);
+}
 
 TEST(CropFilterTest, test_crop_polygon)
 {