Skip to content

Commit

Permalink
Streaming for crop filter.
Browse files Browse the repository at this point in the history
  • Loading branch information
@abellgithub
abellgithub committed Dec 14, 2015
1 parent 736993b commit 16c3063
Show file tree
Hide file tree
Showing 4 changed files with 220 additions and 35 deletions.
89 changes: 68 additions & 21 deletions filters/crop/CropFilter.cpp
View file
Expand Up @@ -237,6 +237,32 @@ Options CropFilter::getDefaultOptions()
}


void CropFilter::ready(PointTableRef table)
{
// If we're using view-based SRSs, handle geometry in the view.
if (table.supportsView())
return;

for (auto& geom : m_geoms)
preparePolygon(geom, table.anySpatialReference());
}


bool CropFilter::processOne(PointRef& point)
{
#ifdef PDAL_HAVE_GEOS
for (auto& geom : m_geoms)
if (!crop(point, geom))
return false;
#endif
for (auto& box : m_bounds)
if (!crop(point, box))
return false;

return true;
}


PointViewSet CropFilter::run(PointViewPtr view)
{
PointViewSet viewSet;
Expand Down Expand Up @@ -276,15 +302,23 @@ PointViewSet CropFilter::run(PointViewPtr view)
}


bool CropFilter::crop(PointRef& point, const BOX2D& box)
{
double x = point.getFieldAs<double>(Dimension::Id::X);
double y = point.getFieldAs<double>(Dimension::Id::Y);

// Return true if we're keeping a point.
return (m_cropOutside != box.contains(x, y));
}


void CropFilter::crop(const BOX2D& box, PointView& input, PointView& output)
{
PointRef point = input.point(0);
for (PointId idx = 0; idx < input.size(); ++idx)
{
double x = input.getFieldAs<double>(Dimension::Id::X, idx);
double y = input.getFieldAs<double>(Dimension::Id::Y, idx);

bool contained = box.contains(x, y);
if (m_cropOutside != box.contains(x, y))
point.setPointId(idx);
if (crop(point, box))
output.appendPoint(input, idx);
}
}
Expand All @@ -311,33 +345,46 @@ GEOSGeometry *CropFilter::createPoint(double x, double y, double z)
}


void CropFilter::crop(const GeomPkg& g, PointView& input, PointView& output)
bool CropFilter::crop(PointRef& point, const GeomPkg& g)
{
bool logOutput = (log()->getLevel() > LogLevel::Debug4);

//ABELL - This makes no sense. Why are we setting this here but never
// using it, as it gets reset to 10 below.
if (logOutput)
log()->floatPrecision(8);

for (PointId idx = 0; idx < input.size(); ++idx)
double x = point.getFieldAs<double>(Dimension::Id::X);
double y = point.getFieldAs<double>(Dimension::Id::Y);
double z = point.getFieldAs<double>(Dimension::Id::Z);

if (logOutput)
{
double x = input.getFieldAs<double>(Dimension::Id::X, idx);
double y = input.getFieldAs<double>(Dimension::Id::Y, idx);
double z = input.getFieldAs<double>(Dimension::Id::Z, idx);
log()->floatPrecision(10);
log()->get(LogLevel::Debug5) << "input: " << x << " y: " << y <<
" z: " << z << std::endl;
}

if (logOutput)
{
log()->floatPrecision(10);
log()->get(LogLevel::Debug5) << "input: " << x << " y: " << y <<
" z: " << z << std::endl;
}
GEOSGeometry *p = createPoint(x, y, z);
bool covers = (bool)(GEOSPreparedCovers_r(m_geosEnvironment,
g.m_prepGeom, p));
bool keep = (m_cropOutside != covers);
GEOSGeom_destroy_r(m_geosEnvironment, p);
return keep;
}

GEOSGeometry *p = createPoint(x, y, z);
bool covers = (bool)(GEOSPreparedCovers_r(m_geosEnvironment,
g.m_prepGeom, p));
if (m_cropOutside != covers)

void CropFilter::crop(const GeomPkg& g, PointView& input, PointView& output)
{
PointRef point = input.point(0);
for (PointId idx = 0; idx < input.size(); ++idx)
{
point.setPointId(idx);
if (crop(point, g))
output.appendPoint(input, idx);
GEOSGeom_destroy_r(m_geosEnvironment, p);
}
}

#endif


Expand Down
4 changes: 4 additions & 0 deletions filters/crop/CropFilter.hpp
View file
Expand Up @@ -86,9 +86,13 @@ class PDAL_DLL CropFilter : public Filter
std::vector<GeomPkg> m_geoms;

virtual void processOptions(const Options& options);
virtual void ready(PointTableRef table);
virtual bool processOne(PointRef& point);
virtual PointViewSet run(PointViewPtr view);
virtual void done(PointTableRef table);
bool crop(PointRef& point, const BOX2D& box);
void crop(const BOX2D& box, PointView& input, PointView& output);
bool crop(PointRef& point, const GeomPkg& g);
void crop(const GeomPkg& g, PointView& input, PointView& output);
#ifdef PDAL_HAVE_GEOS
GEOSGeometry *validatePolygon(const std::string& poly);
Expand Down
2 changes: 1 addition & 1 deletion src/Stage.cpp
View file
Expand Up @@ -198,7 +198,7 @@ void Stage::execute(StreamPointTable& table)
bool finished = false;
while (!finished)
{
// Clear the spatial
// Clear the spatial reference when processing starts.
table.clearSpatialReferences();
PointId idx = 0;
PointRef point(&table, idx);
Expand Down
160 changes: 147 additions & 13 deletions test/unit/filters/CropFilterTest.cpp
View file
Expand Up @@ -37,11 +37,13 @@
#include <pdal/util/FileUtils.hpp>
#include <pdal/PointView.hpp>
#include <pdal/StageFactory.hpp>
#include <pdal/BufferReader.hpp>
#include <CropFilter.hpp>
#include <FauxReader.hpp>
#include <LasReader.hpp>
#include <ReprojectionFilter.hpp>
#include <StatsFilter.hpp>
#include <StreamCallbackFilter.hpp>
#include "Support.hpp"

using namespace pdal;
Expand Down Expand Up @@ -200,26 +202,158 @@ TEST(CropFilterTest, test_crop_polygon_reprojection)
#endif
}

/**
TEST(CropFilterTest, multibounds)
{
using namespace Dimension;

PointTable table;
table.layout->registerDim(Id::X);
table.layout->registerDim(Id::Y);
table.layout->registerDim(Id::Z);
table.layout()->registerDim(Id::X);
table.layout()->registerDim(Id::Y);
table.layout()->registerDim(Id::Z);

PointView view(table);
view.setField(Id::X, 0, 1);
view.setField(Id::Y, 0, 1);
PointViewPtr view(new PointView(table));
view->setField(Id::X, 0, 2);
view->setField(Id::Y, 0, 2);

view->setField(Id::X, 1, 4);
view->setField(Id::Y, 1, 2);

view->setField(Id::X, 2, 6);
view->setField(Id::Y, 2, 2);

view->setField(Id::X, 3, 8);
view->setField(Id::Y, 3, 2);

view->setField(Id::X, 4, 10);
view->setField(Id::Y, 4, 2);

view->setField(Id::X, 5, 12);
view->setField(Id::Y, 5, 2);

BufferReader r;
r.addView(view);

CropFilter crop;
Options o;
o.add("bounds", "([1, 3], [1, 3])");
o.add("bounds", "([5, 7], [1, 3])");
o.add("polygon", "POLYGON ((9 1, 11 1, 11 3, 9 3, 9 1))");
crop.setInput(r);
crop.setOptions(o);

crop.prepare(table);
PointViewSet s = crop.execute(table);
// Make sure we get three views, one with the point 2, 2, one with the
// point 6, 2 and one with 10,2.
EXPECT_EQ(s.size(), 3u);
static int total_cnt = 0;
for (auto v : s)
{
int cnt = 0;
EXPECT_EQ(v->size(), 1u);
double x = v->getFieldAs<double>(Dimension::Id::X, 0);
double y = v->getFieldAs<double>(Dimension::Id::Y, 0);
if (x == 2 && y == 2)
cnt = 1;
if (x == 6 && y == 2)
cnt = 2;
if (x == 10 && y == 2)
cnt = 4;
EXPECT_TRUE(cnt > 0);
total_cnt += cnt;
}
EXPECT_EQ(total_cnt, 7);
}

view.setField(Id::X, 1, 2);
view.setField(Id::Y, 1, 6);
TEST(CropFilterTest, stream)
{
using namespace Dimension;

view.setField(Id::X, 2, 4);
view.setField(Id::Y, 2, 4);
FixedPointTable table(2);
table.layout()->registerDim(Id::X);
table.layout()->registerDim(Id::Y);
table.layout()->registerDim(Id::Z);

class StreamReader : public Reader
{
public:
std::string getName() const
{ return "readers.stream"; }
bool processOne(PointRef& point)
{
static int i = 0;

if (i == 0)
{
point.setField(Id::X, 2);
point.setField(Id::Y, 2);
}
else if (i == 1)
{
point.setField(Id::X, 6);
point.setField(Id::Y, 2);
}
else if (i == 2)
{
point.setField(Id::X, 8);
point.setField(Id::Y, 2);
}
else if (i == 3)
{
point.setField(Id::X, 10);
point.setField(Id::Y, 2);
}
else if (i == 4)
{
point.setField(Id::X, 12);
point.setField(Id::Y, 2);
}
else
return false;
i++;
return true;
}
};

StreamReader r;

BOX3D p
CropFilter crop;
Options o;
o.add("bounds", "([1, 3], [1, 3])");
o.add("bounds", "([5, 7], [1, 3])");
o.add("polygon", "POLYGON ((9 1, 11 1, 11 3, 9 3, 9 1))");
crop.setInput(r);
crop.setOptions(o);

auto cb = [](PointRef& point)
{
static int i = 0;
if (i == 0)
{
EXPECT_EQ(point.getFieldAs<int>(Id::X), 2);
EXPECT_EQ(point.getFieldAs<int>(Id::Y), 2);
}
if (i == 1)
{
EXPECT_EQ(point.getFieldAs<int>(Id::X), 6);
EXPECT_EQ(point.getFieldAs<int>(Id::Y), 2);
}
if (i == 2)
{
EXPECT_EQ(point.getFieldAs<int>(Id::X), 10);
EXPECT_EQ(point.getFieldAs<int>(Id::Y), 2);
}
else
return false;
i++;
return true;
};

StreamCallbackFilter f;
f.setCallback(cb);
f.setInput(crop);

f.prepare(table);
f.execute(table);
}
**/

0 comments on commit 16c3063

Please sign in to comment.