Skip to content

Commit

Permalink
Working filters.poisson
Browse files Browse the repository at this point in the history
  • Loading branch information
@abellgithub
abellgithub committed Jun 27, 2017
1 parent 28fe672 commit 0b58e05
Show file tree
Hide file tree
Showing 56 changed files with 21,396 additions and 18 deletions.
4 changes: 3 additions & 1 deletion CMakeLists.txt
View file
Expand Up @@ -58,7 +58,7 @@ endif()
set(PDAL_UTIL_LIB_NAME pdal_util)
set(PDAL_PLANG_LIB_NAME pdal_plang)
set(PDAL_BOOST_LIB_NAME pdal_boost)
#set(PDAL_ARBITER_LIB_NAME pdal_arbiter)
set(PDAL_KAZHDAN_LIB_NAME pdal_kazhdan)

set(CMAKE_INCLUDE_DIRECTORIES_PROJECT_BEFORE ON)

Expand Down Expand Up @@ -172,6 +172,7 @@ endif()
add_subdirectory(dimbuilder)
add_subdirectory(vendor/pdalboost)
add_subdirectory(vendor/arbiter)
add_subdirectory(vendor/kazhdan)

if (NOT PDAL_HAVE_JSONCPP)
add_subdirectory(vendor/jsoncpp/dist)
Expand Down Expand Up @@ -251,6 +252,7 @@ target_link_libraries(${PDAL_BASE_LIB_NAME}
${PDAL_REEXPORT}
${PDAL_UTIL_LIB_NAME}
${PDAL_ARBITER_LIB_NAME}
${PDAL_KAZHDAN_LIB_NAME}
${JSON_CPP_LINK_TYPE}
${PDAL_JSONCPP_LIB_NAME}
INTERFACE
Expand Down
38 changes: 26 additions & 12 deletions doc/stages/filters.poisson.rst
View file
Expand Up @@ -4,14 +4,26 @@
filters.poisson
===============================================================================

The Poisson filter passes data through the Point Cloud Library (`PCL`_) Poisson
surface reconstruction algorithm.

Poisson is an implementation of the method described in [Kazhdan2006]_.
The poisson filter passes data Mischa Kazhdan's poisson surface reconstruction
algorithm. [Kazhdan2006]_ It creates a watertight surface from the original
point set by creating an entirely new point set representing the imputed
isosurface. The algorithm requires normal vectors to each point in order
to run. If the x, y and z normal dimensions are present in the input point
set, they will be used by the algorithm. If they don't exist, the poisson
filter will invoke the PDAL normal filter to create them before running.

The poisson algorithm will usually create a larger output point set
than the input point set. Because the algorithm constructs new points, data
associated with the original points set will be lost, as the algorithm has
limited ability to impute associated data. However, if color dimensions
(red, green and blue) are present in the input, colors will be reconstruced
in the output point set.

.. [Kazhdan2006] Kazhdan, Michael, Matthew Bolitho, and Hugues Hoppe. "Poisson surface reconstruction." Proceedings of the fourth Eurographics symposium on Geometry processing. Vol. 7. 2006.
.. _`PCL`: http://www.pointclouds.org
This integration of the algorithm with PDAL only supports a limited set of
the options available to the implementation. If you need support for further
options, please let us know.

Example
-------------------------------------------------------------------------------
Expand All @@ -23,12 +35,10 @@ Example
"dense.las",
{
"type":"filters.poisson",
"depth":"8",
"point_weight":"4"
},
{
"type":"writers.las",
"filename":"thinned.las",
"type":"writers.ply",
"filename":"isosurface.ply",
}
]
}
Expand All @@ -37,8 +47,12 @@ Example
Options
-------------------------------------------------------------------------------

density
Write an estimate of neighborhood density for each point in the output
set.

depth
Maximum depth of the tree used for reconstruction. [Default: **8**]
Maximum depth of the tree used for reconstruction. The output is sentsitve
to this parameter. Increase if the results appear unsatisfactory.
[Default: **8**]

point_weight
Importance of interpolation of point samples in the screened Poisson equation. [Default: **4.0**]
230 changes: 227 additions & 3 deletions filters/PoissonFilter.cpp
View file
Expand Up @@ -33,26 +33,250 @@
****************************************************************************/

#include "PoissonFilter.hpp"
#include "NormalFilter.hpp"

#include <pdal/pdal_macros.hpp>

#include <kazhdan/PoissonRecon.h>
#include <kazhdan/point_source/PointSource.h>

namespace pdal
{

class PointViewSource : public PointSource
{
public:
PointViewSource(PointView& view) : m_view(view), m_current(0)
{}

virtual void reset()
{ m_current = 0; }
virtual bool nextPoint(Point& point)
{
using namespace Dimension;

if (m_current >= m_view.size())
return false;
point.p.coords[0] = m_view.getFieldAs<double>(Id::X, m_current);
point.p.coords[1] = m_view.getFieldAs<double>(Id::Y, m_current);
point.p.coords[2] = m_view.getFieldAs<double>(Id::Z, m_current);
point.n.coords[0] = m_view.getFieldAs<double>(Id::NormalX, m_current);
point.n.coords[1] = m_view.getFieldAs<double>(Id::NormalY, m_current);
point.n.coords[2] = m_view.getFieldAs<double>(Id::NormalZ, m_current);
m_current++;
return true;
}

private:
PointView& m_view;
PointId m_current;
};

class ColorPointViewSource : public ColorPointSource
{
public:
ColorPointViewSource(PointView& view) : m_view(view), m_current(0)
{}

virtual void reset()
{ m_current = 0; }
virtual bool nextPoint(Point& point, Point3D<double>& color)
{
using namespace Dimension;

if (m_current >= m_view.size())
return false;
point.p.coords[0] = m_view.getFieldAs<double>(Id::X, m_current);
point.p.coords[1] = m_view.getFieldAs<double>(Id::Y, m_current);
point.p.coords[2] = m_view.getFieldAs<double>(Id::Z, m_current);
point.n.coords[0] = m_view.getFieldAs<double>(Id::NormalX, m_current);
point.n.coords[1] = m_view.getFieldAs<double>(Id::NormalY, m_current);
point.n.coords[2] = m_view.getFieldAs<double>(Id::NormalZ, m_current);
color.coords[0] = m_view.getFieldAs<double>(Id::Red, m_current);
color.coords[1] = m_view.getFieldAs<double>(Id::Green, m_current);
color.coords[2] = m_view.getFieldAs<double>(Id::Blue, m_current);
m_current++;
return true;
}

private:
PointView& m_view;
PointId m_current;
};

class PointViewMesh : public Kazhdan::Mesh
{
public:
PointViewMesh(PointView& view, bool color) : m_view(view),
m_mesh(*(m_view.createMesh("poisson"))), m_doColor(color)
{ resetIterator(); }

virtual int pointCount() const
{ return m_view.size(); }
virtual int polygonCount() const
{ return m_mesh.size(); }
virtual int newPoint(const std::array<double, 3>& position)
{
PointId cnt = m_view.size();
m_view.setField(Dimension::Id::X, cnt, position[0]);
m_view.setField(Dimension::Id::Y, cnt, position[1]);
m_view.setField(Dimension::Id::Z, cnt, position[2]);
return cnt;
}

virtual int newPoint(const std::array<double, 3>& position, double density)
{
PointId cnt = m_view.size();
m_view.setField(Dimension::Id::X, cnt, position[0]);
m_view.setField(Dimension::Id::Y, cnt, position[1]);
m_view.setField(Dimension::Id::Z, cnt, position[2]);
m_view.setField(Dimension::Id::Density, cnt, density);
return cnt;
}

virtual int newPoint(const std::array<double, 3>& position,
const std::array<uint8_t, 3>& color, double density)
{
PointId cnt = m_view.size();
m_view.setField(Dimension::Id::X, cnt, position[0]);
m_view.setField(Dimension::Id::Y, cnt, position[1]);
m_view.setField(Dimension::Id::Z, cnt, position[2]);
m_view.setField(Dimension::Id::Red, cnt, color[0]);
m_view.setField(Dimension::Id::Green, cnt, color[1]);
m_view.setField(Dimension::Id::Blue, cnt, color[2]);
m_view.setField(Dimension::Id::Density, cnt, density);
return cnt;
}

virtual void newPolygon(std::vector<int>& poly)
{
assert(poly.size() == 3);
m_mesh.add(poly[0], poly[1], poly[2]);
}

virtual bool hasDensity() const
{ return m_view.hasDim(Dimension::Id::Density); }

virtual void resetIterator()
{
m_polyIdx = 0;
m_pointIdx = 0;
}

virtual bool nextPolygon(Kazhdan::Polygon& poly)
{
if (m_polyIdx >= m_mesh.size())
return false;

const Triangle& t = m_mesh[m_polyIdx];
poly.insert(poly.end(), { (int)t.m_a, (int)t.m_b, (int)t.m_c } );
m_polyIdx++;
return true;
}

virtual bool nextPoint(Kazhdan::Point& point)
{
if (m_pointIdx > m_view.size())
return false;
point.m_position[0] = m_view.getFieldAs<double>(Dimension::Id::X,
m_pointIdx);
point.m_position[1] = m_view.getFieldAs<double>(Dimension::Id::Y,
m_pointIdx);
point.m_position[2] = m_view.getFieldAs<double>(Dimension::Id::Z,
m_pointIdx);
point.m_density = m_view.getFieldAs<double>(Dimension::Id::Density,
m_pointIdx);
if (m_doColor)
{
point.m_color[0] = m_view.getFieldAs<uint8_t>(Dimension::Id::Red,
m_pointIdx);
point.m_color[1] = m_view.getFieldAs<uint8_t>(Dimension::Id::Green,
m_pointIdx);
point.m_color[2] = m_view.getFieldAs<uint8_t>(Dimension::Id::Blue,
m_pointIdx);
}
m_pointIdx++;
return true;
}

private:
PointView& m_view;
TriangularMesh& m_mesh;
size_t m_polyIdx;
size_t m_pointIdx;
bool m_doColor;
};

static PluginInfo const s_info =
PluginInfo("filters.poisson", "Poisson Surface Reconstruction Filter",
"http://pdal.io/stages/filters.poisson.html");

CREATE_STATIC_PLUGIN(1, 0, PoissonFilter, Filter, s_info)

std::string PoissonFilter::getName() const
std::string PoissonFilter::getName() const { return s_info.name; }

void PoissonFilter::addDimensions(PointLayoutPtr layout)
{
return s_info.name;
if (layout->hasDim(Dimension::Id::Red) &&
layout->hasDim(Dimension::Id::Green) &&
layout->hasDim(Dimension::Id::Blue))
m_doColor = true;

if (layout->hasDim(Dimension::Id::NormalX))
{
if ((!layout->hasDim(Dimension::Id::NormalY)) ||
(!layout->hasDim(Dimension::Id::NormalZ)))
throwError("If normals are provided as part of the input "
"dataset, all of X, Y and Z must be provided.");
m_normalsProvided = true;
}
else
layout->registerDims( {Dimension::Id::NormalX, Dimension::Id::NormalY,
Dimension::Id::NormalZ} );
}


void PoissonFilter::filter(PointView& view)
void PoissonFilter::addArgs(ProgramArgs& args)
{
args.add("density", "Output density estimates", m_density);
args.add("depth", "Maximum depth of the octree used for reconstruction",
m_depth, 8);
}


PointViewSet PoissonFilter::run(PointViewPtr view)
{
if (!m_normalsProvided)
NormalFilter().doFilter(*view);

std::unique_ptr<PointSource> source;
if (m_doColor)
source.reset(new ColorPointViewSource(*view));
else
source.reset(new PointViewSource(*view));

PoissonOpts<double> opts;

opts.m_depth = m_depth;
opts.m_density = m_density;
opts.m_solveDepth = m_depth;
opts.m_kernelDepth = m_depth - 2;
if (m_doColor)
{
opts.m_color = 16;
opts.m_hasColor = true;
}

PoissonRecon<double> recon(opts, *source);
recon.execute();
recon.evaluate();

PointViewSet s;
PointViewPtr outView = view->makeNew();
s.insert(outView);
PointViewMesh mesh(*outView, m_doColor);
recon.extractMesh(mesh);
return s;
}

} // namespace pdal
10 changes: 8 additions & 2 deletions filters/PoissonFilter.hpp
View file
Expand Up @@ -46,7 +46,7 @@ namespace pdal
class PDAL_DLL PoissonFilter : public Filter
{
public:
PoissonFilter() : Filter()
PoissonFilter() : Filter(), m_normalsProvided(false)
{}
PoissonFilter& operator=(const PoissonFilter&) = delete;
PoissonFilter(const PoissonFilter&) = delete;
Expand All @@ -56,8 +56,14 @@ class PDAL_DLL PoissonFilter : public Filter
std::string getName() const;

private:
virtual void filter(PointView& view);
bool m_density;
int m_depth;
bool m_normalsProvided;
bool m_doColor;

virtual void addDimensions(PointLayoutPtr layout);
virtual PointViewSet run(PointViewPtr view);
virtual void addArgs(ProgramArgs& args);
};

} // namespace pdal
9 changes: 9 additions & 0 deletions pdal/Dimension.json
View file
Expand Up @@ -340,22 +340,31 @@
},
{
"name": "NormalX",
"alt_names": "nx",
"type": "double",
"description": "X component of a vector normal to surface at this point"
},
{
"name": "NormalY",
"alt_names": "ny",
"type": "double",
"description": "Y component of a vector normal to surface at this point"
},
{
"name": "NormalZ",
"alt_names": "nz",
"type": "double",
"description": "Z component of a vector normal to surface at this point"
},
{
"name": "Curvature",
"type": "double",
"description": "Curvature of surface at this point"
},
{
"name": "Density",
"type": "double",
"description": "Estimate of point density"
}

] }

0 comments on commit 0b58e05

Please sign in to comment.