Skip to content

Commit

Permalink
Minor updates to several filters
Browse files Browse the repository at this point in the history
Initial focus is on updating static filters to always register official
Dimensions.

* LiTreeFilter now registers existing ClusterID instead of custom
  TreeID.

* LOFFilter now registers existing NNDistance instead of custom
  KDistance.

* OptimalNeighborhoodFilter now registers existing OptimalRadius and
  OptimalKNN.

* Register the following as "official" Dimensions: Coplanar,
  Eigenvalue0, Eigenvalue1, Eigenvalue2, LocalOutlierFactor,
  LocalReachabilityDistance, Miniball, PlaneFit, RadialDensity, Rank,
  Reciprocity

While updating Dimensions, also took the opportunity to:

* Be more explicit about datatypes, instead of relying on auto, to
  improve readability.

* Use PointView range-based for loops where it make sense. Generally
  makes code a little more succint.

* Do not pass KDIndex as frequently, when it can typically be grabbed
  off the view at no cost.
  • Loading branch information
chambbj committed Oct 21, 2020
1 parent c4c5885 commit 043e2f3
Show file tree
Hide file tree
Showing 23 changed files with 195 additions and 169 deletions.
9 changes: 7 additions & 2 deletions doc/stages/filters.litree.rst
Original file line number Diff line number Diff line change
Expand Up @@ -6,15 +6,20 @@ filters.litree

The purpose of the Li tree filter is to segment individual trees from an input
``PointView``. In the output ``PointView`` points that are deemed to be part of
a tree are labeled with a ``TreeID``. Tree IDs start at 1, with non-tree points
given a ``TreeID`` of 0.
a tree are labeled with a ``ClusterID``. Tree IDs start at 1, with non-tree points
given a ``ClusterID`` of 0.

.. note::

The filter differs only slightly from the paper in the addition of a few
conditions on size of tree, minimum height above ground for tree seeding, and
flexible radius for non-tree seed insertion.

.. note::

In earlier PDAL releases, ``ClusterID`` was stored in the ``TreeID``
Dimemsion.

.. embed::

Example
Expand Down
11 changes: 8 additions & 3 deletions doc/stages/filters.lof.rst
Original file line number Diff line number Diff line change
Expand Up @@ -8,15 +8,15 @@ of determining the degree to which an object is an outlier. This filter
is an implementation of the method
described in [Breunig2000]_.

The filter creates three new dimensions, ``KDistance``,
The filter creates three new dimensions, ``NNDistance``,
``LocalReachabilityDistance`` and ``LocalOutlierFactor``, all of which are
double-precision floating values. The ``KDistance`` dimension records the
double-precision floating values. The ``NNDistance`` dimension records the
Euclidean distance between a point and it's k-th nearest neighbor (the number
of k neighbors is set with the minpts_ option). The
``LocalReachabilityDistance`` is the inverse of the mean
of all reachability distances for a neighborhood of points. This reachability
distance is defined as the max of the Euclidean distance to a neighboring point
and that neighbor's own previously computed ``KDistance``. Finally, each point
and that neighbor's own previously computed ``NNDistance``. Finally, each point
has a ``LocalOutlierFactor`` which is the mean of all
``LocalReachabilityDistance`` values for the neighborhood. In each case, the
neighborhood is the set of k nearest neighbors.
Expand All @@ -32,6 +32,11 @@ users of this filter should find instructive.
To inspect the newly created, non-standard dimensions, be sure to write to an
output format that can support arbitrary dimensions, such as BPF.

.. note::

In earlier PDAL releases, ``NNDistance`` was stored in the ``KDistance``
Dimemsion.

.. embed::

Example
Expand Down
28 changes: 14 additions & 14 deletions filters/ApproximateCoplanarFilter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,13 @@
#include <Eigen/Dense>

#include <string>
#include <vector>

namespace pdal
{

using namespace Dimension;
using namespace Eigen;

static StaticPluginInfo const s_info
{
"filters.approximatecoplanar",
Expand All @@ -71,35 +73,33 @@ void ApproximateCoplanarFilter::addArgs(ProgramArgs& args)

void ApproximateCoplanarFilter::addDimensions(PointLayoutPtr layout)
{
m_coplanar = layout->registerOrAssignDim("Coplanar",
Dimension::Type::Unsigned8);
layout->registerDim(Id::Coplanar);
}


void ApproximateCoplanarFilter::filter(PointView& view)
{
using namespace Eigen;

KD3Index& kdi = view.build3dIndex();
const KD3Index& kdi = view.build3dIndex();

for (PointId i = 0; i < view.size(); ++i)
for (PointRef p : view)
{
// find the k-nearest neighbors
auto ids = kdi.neighbors(i, m_knn);
PointIdList ids = kdi.neighbors(p, m_knn);

// compute covariance of the neighborhood
auto B = math::computeCovariance(view, ids);
Matrix3d B = math::computeCovariance(view, ids);

// perform the eigen decomposition
SelfAdjointEigenSolver<Matrix3d> solver(B);
if (solver.info() != Success)
Eigen::SelfAdjointEigenSolver<Matrix3d> solver(B);
if (solver.info() != Eigen::Success)
throwError("Cannot perform eigen decomposition.");
auto ev = solver.eigenvalues();
Vector3d ev = solver.eigenvalues();

// test eigenvalues to label points that are approximately coplanar
if ((ev[1] > m_thresh1 * ev[0]) && (m_thresh2 * ev[1] > ev[2]))
view.setField(m_coplanar, i, 1u);
p.setField(Id::Coplanar, 1u);
else
view.setField(m_coplanar, i, 0u);
p.setField(Id::Coplanar, 0u);
}
}

Expand Down
4 changes: 0 additions & 4 deletions filters/ApproximateCoplanarFilter.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -36,14 +36,11 @@

#include <pdal/Filter.hpp>

#include <cstdint>
#include <memory>
#include <string>

namespace pdal
{

class Options;
class PointLayout;
class PointView;

Expand All @@ -62,7 +59,6 @@ class PDAL_DLL ApproximateCoplanarFilter : public Filter
int m_knn;
double m_thresh1;
double m_thresh2;
Dimension::Id m_coplanar;

virtual void addDimensions(PointLayoutPtr layout);
virtual void addArgs(ProgramArgs& args);
Expand Down
34 changes: 17 additions & 17 deletions filters/EigenvaluesFilter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,13 @@
#include <Eigen/Dense>

#include <string>
#include <vector>

namespace pdal
{

using namespace Dimension;
using namespace Eigen;

static StaticPluginInfo const s_info
{
"filters.eigenvalues",
Expand Down Expand Up @@ -86,9 +88,9 @@ void EigenvaluesFilter::addArgs(ProgramArgs& args)

void EigenvaluesFilter::addDimensions(PointLayoutPtr layout)
{
m_e0 = layout->registerOrAssignDim("Eigenvalue0", Dimension::Type::Double);
m_e1 = layout->registerOrAssignDim("Eigenvalue1", Dimension::Type::Double);
m_e2 = layout->registerOrAssignDim("Eigenvalue2", Dimension::Type::Double);
layout->registerDim(Id::Eigenvalue0);
layout->registerDim(Id::Eigenvalue1);
layout->registerDim(Id::Eigenvalue2);
}

void EigenvaluesFilter::prepared(PointTableRef table)
Expand All @@ -112,17 +114,15 @@ void EigenvaluesFilter::prepared(PointTableRef table)

void EigenvaluesFilter::filter(PointView& view)
{
using namespace Eigen;

KD3Index& kdi = view.build3dIndex();
const KD3Index& kdi = view.build3dIndex();

for (PointId i = 0; i < view.size(); ++i)
for (PointRef p : view)
{
// find neighbors, either by radius or k nearest neighbors
PointIdList ids;
if (m_args->m_radiusArg->set())
{
ids = kdi.radius(i, m_args->m_radius);
ids = kdi.radius(p, m_args->m_radius);

// if insufficient number of neighbors, eigen solver will fail
// anyway, it may be okay to silently return without setting any of
Expand All @@ -132,27 +132,27 @@ void EigenvaluesFilter::filter(PointView& view)
}
else
{
ids = kdi.neighbors(i, m_args->m_knn + 1, m_args->m_stride);
ids = kdi.neighbors(p, m_args->m_knn + 1, m_args->m_stride);
}

// compute covariance of the neighborhood
auto B = math::computeCovariance(view, ids);
Matrix3d B = math::computeCovariance(view, ids);

// perform the eigen decomposition
SelfAdjointEigenSolver<Matrix3d> solver(B);
if (solver.info() != Success)
Eigen::SelfAdjointEigenSolver<Matrix3d> solver(B);
if (solver.info() != Eigen::Success)
throwError("Cannot perform eigen decomposition.");
auto ev = solver.eigenvalues();
Vector3d ev = solver.eigenvalues();

if (m_args->m_normalize)
{
double sum = ev[0] + ev[1] + ev[2];
ev /= sum;
}

view.setField(m_e0, i, ev[0]);
view.setField(m_e1, i, ev[1]);
view.setField(m_e2, i, ev[2]);
p.setField(Id::Eigenvalue0, ev[0]);
p.setField(Id::Eigenvalue1, ev[1]);
p.setField(Id::Eigenvalue2, ev[2]);
}
}

Expand Down
4 changes: 0 additions & 4 deletions filters/EigenvaluesFilter.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -36,14 +36,11 @@

#include <pdal/Filter.hpp>

#include <cstdint>
#include <memory>
#include <string>

namespace pdal
{

class Options;
class PointLayout;
class PointView;
struct EigenvalueArgs;
Expand All @@ -59,7 +56,6 @@ class PDAL_DLL EigenvaluesFilter : public Filter

private:
std::unique_ptr<EigenvalueArgs> m_args;
Dimension::Id m_e0, m_e1, m_e2;

virtual void addDimensions(PointLayoutPtr layout);
virtual void addArgs(ProgramArgs& args);
Expand Down
15 changes: 7 additions & 8 deletions filters/EstimateRankFilter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,6 @@
#include "EstimateRankFilter.hpp"

#include <string>
#include <vector>

#include <pdal/KDIndex.hpp>
#include <pdal/util/ProgramArgs.hpp>
Expand All @@ -44,6 +43,8 @@
namespace pdal
{

using namespace Dimension;

static StaticPluginInfo const s_info
{
"filters.estimaterank",
Expand All @@ -68,19 +69,17 @@ void EstimateRankFilter::addArgs(ProgramArgs& args)

void EstimateRankFilter::addDimensions(PointLayoutPtr layout)
{
m_rank = layout->registerOrAssignDim("Rank", Dimension::Type::Unsigned8);
layout->registerDim(Id::Rank);
}

void EstimateRankFilter::filter(PointView& view)
{
KD3Index& kdi = view.build3dIndex();
const KD3Index& kdi = view.build3dIndex();

for (PointId i = 0; i < view.size(); ++i)
for (PointRef p : view)
{
// find the k-nearest neighbors
auto ids = kdi.neighbors(i, m_knn);

view.setField(m_rank, i, math::computeRank(view, ids, m_thresh));
PointIdList ids = kdi.neighbors(p, m_knn);
p.setField(Id::Rank, math::computeRank(view, ids, m_thresh));
}
}

Expand Down
4 changes: 0 additions & 4 deletions filters/EstimateRankFilter.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -36,14 +36,11 @@

#include <pdal/Filter.hpp>

#include <cstdint>
#include <memory>
#include <string>

namespace pdal
{

class Options;
class PointLayout;
class PointView;

Expand All @@ -60,7 +57,6 @@ class PDAL_DLL EstimateRankFilter : public Filter
private:
int m_knn;
double m_thresh;
Dimension::Id m_rank;

virtual void addDimensions(PointLayoutPtr layout);
virtual void addArgs(ProgramArgs& args);
Expand Down
28 changes: 15 additions & 13 deletions filters/LOFFilter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,8 @@
namespace pdal
{

using namespace Dimension;

static StaticPluginInfo const s_info
{
"filters.lof",
Expand All @@ -63,17 +65,14 @@ void LOFFilter::addArgs(ProgramArgs& args)

void LOFFilter::addDimensions(PointLayoutPtr layout)
{
using namespace Dimension;
m_kdist = layout->registerOrAssignDim("KDistance", Type::Double);
m_lrd = layout->registerOrAssignDim("LocalReachabilityDistance", Type::Double);
m_lof = layout->registerOrAssignDim("LocalOutlierFactor", Type::Double);
layout->registerDim(Id::NNDistance);
layout->registerDim(Id::LocalReachabilityDistance);
layout->registerDim(Id::LocalOutlierFactor);
}

void LOFFilter::filter(PointView& view)
{
using namespace Dimension;

KD3Index& index = view.build3dIndex();
const KD3Index& index = view.build3dIndex();

// Increment the minimum number of points, as knnSearch will be returning
// the neighbors along with the query point.
Expand All @@ -97,7 +96,7 @@ void LOFFilter::filter(PointView& view)
for (size_t j = 0; j < m_minpts; ++j)
mat[std::make_pair(i, j)] = std::make_pair(indices[j], std::sqrt(sqr_dists[j]));

view.setField(m_kdist, i, std::sqrt(sqr_dists[m_minpts - 1]));
view.setField(Id::NNDistance, i, std::sqrt(sqr_dists[m_minpts - 1]));
}

// Second pass: Compute the local reachability distance for each point.
Expand All @@ -113,27 +112,30 @@ void LOFFilter::filter(PointView& view)
for (size_t j = 0; j < m_minpts; ++j)
{
auto val = mat[std::make_pair(i, j)];
double k = view.getFieldAs<double>(m_kdist, val.first);
double k = view.getFieldAs<double>(Id::NNDistance, val.first);
double reachdist = (std::max)(k, val.second);
M1 += (reachdist - M1) / ++n;
}
view.setField(m_lrd, i, 1.0 / M1);
view.setField(Id::LocalReachabilityDistance, i, 1.0 / M1);
}

// Third pass: Compute the local outlier factor for each point.
// The LOF is the average of the lrd's for a neighborhood of points.
log()->get(LogLevel::Debug) << "Computing LOF...\n";
for (PointId i = 0; i < view.size(); ++i)
{
double lrdp = view.getFieldAs<double>(m_lrd, i);
double lrdp = view.getFieldAs<double>(Id::LocalReachabilityDistance, i);
double M1 = 0.0;
point_count_t n = 0;
for (size_t j = 0; j < m_minpts; ++j)
{
auto val = mat[std::make_pair(i, j)];
M1 += (view.getFieldAs<double>(m_lrd, val.first) / lrdp - M1) / ++n;
PointRef p = view.point(val.first);
double ratio =
p.getFieldAs<double>(Id::LocalReachabilityDistance) / lrdp;
M1 += (ratio - M1) / ++n;
}
view.setField(m_lof, i, M1);
view.setField(Id::LocalOutlierFactor, i, M1);
}
}

Expand Down

0 comments on commit 043e2f3

Please sign in to comment.