Skip to content

Commit

Permalink
Merge pull request #2455 from helix-re/filter/localfeature
Browse files Browse the repository at this point in the history
Local feature filter
  • Loading branch information
chambbj committed Apr 16, 2019
2 parents dcd59b8 + e3d380a commit 9ee9211
Show file tree
Hide file tree
Showing 7 changed files with 487 additions and 12 deletions.
4 changes: 4 additions & 0 deletions doc/references.rst
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,12 @@ Reference
.. [Cook1986] Cook, Robert L. "Stochastic sampling in computer graphics." *ACM Transactions on Graphics (TOG)* 5.1 (1986): 51-72.
.. [Demantke2011] Demantké J., Mallet C., David N., Vallet, B. "Dimensionality Based Scale Selection in 3d LIDAR Point Clouds." Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci, XXXVIII-5/W12, 97-102, 2011
.. [Dippe1985] Dippé, Mark AZ, and Erling Henry Wold. "Antialiasing through stochastic sampling." *ACM Siggraph Computer Graphics* 19.3 (1985): 69-78.
.. [Guinard2017] Guinard S., Landrieu L. "Weakly Supervised Segmented-Aided Classification of Urban Scenes From 3D LIDAR Point Clouds." Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci., XLII-1/W1, 151-157, 2017
.. [Kazhdan2006] Kazhdan, Michael, Matthew Bolitho, and Hugues Hoppe. "Poisson surface reconstruction." Proceedings of the fourth Eurographics symposium on Geometry processing. Vol. 7. 2006.
.. [Li2010] Li, Ruosi, et al. "Polygonizing extremal surfaces with manifold guarantees." Proceedings of the 14th ACM Symposium on Solid and Physical Modeling. ACM, 2010.
Expand Down
65 changes: 65 additions & 0 deletions doc/stages/filters.covariancefeatures.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
.. _filters.covariancefeatures:

===============================================================================
filters.covariancefeatures
===============================================================================

This filter implements various local feature descriptors introduced that are based on the covariance matrix of a
point's neighborhood. The user can pick a set of feature descriptors by
setting the ``feature_set`` option. Currently, the only supported feature is the dimensionality_
set of feature descriptors introduced below.

Example
-------------------------------------------------------------------------------

.. code-block:: json
[
"input.las",
{
"type":"filters.covariancefeatures",
"knn":8,
"threads": 2,
"feature_set": "Dimensionality"
},
{
"type":"writers.bpf",
"filename":"output.las",
"output_dims":"X,Y,Z,Linearity,Planarity,Scattering,Verticality"
}
]
Options
-------------------------------------------------------------------------------

knn
The number of k nearest neighbors used for calculating the covariance matrix. [Default: 10]

threads
The number of threads used for computing the feature descriptors. [Default: 1]

feature_set
The features to be computed. Currently only supports ``Dimensionality``. [Default: "Dimensionality"]

.. _dimensionality:

Dimensionality feature set
................................................................................
The features introduced in [Demantke2011]_ describe the shape
of the neighborhood, indicating whether
the local geometry is more linear (1D), planar (2D) or volumetric (3D) while the one introduced in
[Guinard2017]_ adds the idea of a structure being vertical.


The dimensionality filter introduces the following four descriptors that are computed from the covariance matrix of the ``knn`` neighbors:

* linearity - higher for long thin strips
* planarity - higher for planar surfaces
* scattering - higher for complex 3d neighbourhoods
* verticality - higher for vertical structures, highest for thin vertical strips

It introduces four new dimensions that hold each one of these values: ``Linearity`` ``Planarity`` ``Scattering``
and ``Verticality``.



155 changes: 155 additions & 0 deletions filters/CovarianceFeaturesFilter.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,155 @@
/******************************************************************************
* Copyright (c) 2019, Helix Re Inc. nicolas@helix.re
*
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following
* conditions are met:
*
* * Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* * Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in
* the documentation and/or other materials provided
* with the distribution.
* * Neither the name of Helix Re Inc. nor the
* names of its contributors may be used to endorse or promote
* products derived from this software without specific prior
* written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
* OF SUCH DAMAGE.
****************************************************************************/

// This is an implementation of the local feature descriptors introduced in
// WEAKLY SUPERVISED SEGMENTATION-AIDED CLASSIFICATION OF URBANSCENES FROM 3D LIDAR POINT CLOUDS
// Stéphane Guinard, Loïc Landrieu, 2017

#include "CovarianceFeaturesFilter.hpp"

#include <pdal/EigenUtils.hpp>
#include <pdal/KDIndex.hpp>
#include <pdal/util/ProgramArgs.hpp>

#include <Eigen/Dense>

#include <string>
#include <vector>
#include <cmath>

namespace pdal
{

static StaticPluginInfo const s_info
{
"filters.covariancefeatures",
"Filter that calculates local features based on the covariance matrix of a point's neighborhood.",
"http://pdal.io/stages/filters.covariancefeatures.html"
};

CREATE_STATIC_STAGE(CovarianceFeaturesFilter, s_info)

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

void CovarianceFeaturesFilter::addArgs(ProgramArgs& args)
{
args.add("knn", "k-Nearest neighbors", m_knn, 10);
args.add("threads", "Number of threads used to run this filter", m_threads, 1);
args.add("feature_set", "Set of features to be computed", m_featureSet, "Dimensionality");
}

void CovarianceFeaturesFilter::addDimensions(PointLayoutPtr layout)
{
if (m_featureSet == "Dimensionality")
{
for (auto dim: {"Linearity", "Planarity", "Scattering", "Verticality"})
m_extraDims[dim] = layout->registerOrAssignDim(dim, Dimension::Type::Float);
}
}

void CovarianceFeaturesFilter::filter(PointView& view)
{

KD3Index& kdi = view.build3dIndex();

point_count_t nloops = view.size();
std::vector<std::thread> threadPool(m_threads);
for(int t = 0;t<m_threads;t++)
{
threadPool[t] = std::thread(std::bind(
[&](const PointId start, const PointId end, const PointId t)
{
for(PointId i = start;i<end;i++)
setDimensionality(view, i, kdi);
},
t*nloops/m_threads,(t+1)==m_threads?nloops:(t+1)*nloops/m_threads,t));
}
for (auto &t: threadPool)
t.join();
}

void CovarianceFeaturesFilter::setDimensionality(PointView &view, const PointId &id, const KD3Index &kid)
{
using namespace Eigen;

// find the k-nearest neighbors
auto ids = kid.neighbors(id, m_knn + 1);

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

// perform the eigen decomposition
SelfAdjointEigenSolver<Matrix3f> solver(B);
if (solver.info() != Success)
throwError("Cannot perform eigen decomposition.");

// Extract eigenvalues and eigenvectors in decreasing order (largest eigenvalue first)
auto ev = solver.eigenvalues();
std::vector<float> lambda = {(std::max(ev[2],0.f)),
(std::max(ev[1],0.f)),
(std::max(ev[0],0.f))};

if (lambda[0] == 0)
throwError("Eigenvalues are all 0. Can't compute local features.");

auto eigenVectors = solver.eigenvectors();
std::vector<float> v1(3), v2(3), v3(3);
for (int i=0; i < 3; i++)
{
v1[i] = eigenVectors.col(2)(i);
v2[i] = eigenVectors.col(1)(i);
v3[i] = eigenVectors.col(0)(i);
}

float linearity = (sqrtf(lambda[0]) - sqrtf(lambda[1])) / sqrtf(lambda[0]);
float planarity = (sqrtf(lambda[1]) - sqrtf(lambda[2])) / sqrtf(lambda[0]);
float scattering = sqrtf(lambda[2]) / sqrtf(lambda[0]);
view.setField(m_extraDims["Linearity"], id, linearity);
view.setField(m_extraDims["Planarity"], id, planarity);
view.setField(m_extraDims["Scattering"], id, scattering);

std::vector<float> unary_vector(3);
float norm = 0;
for (int i=0; i <3 ; i++)
{
unary_vector[i] = lambda[0] * fabsf(v1[i]) + lambda[1] * fabsf(v2[i]) + lambda[2] * fabsf(v3[i]);
norm += unary_vector[i] * unary_vector[i];
}
norm = sqrtf(norm);
view.setField(m_extraDims["Verticality"], id, unary_vector[2] / norm);
}
}
66 changes: 66 additions & 0 deletions filters/CovarianceFeaturesFilter.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
/******************************************************************************
* Copyright (c) 2019, Helix Re Inc. nicolas@helix.re
*
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following
* conditions are met:
*
* * Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* * Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in
* the documentation and/or other materials provided
* with the distribution.
* * Neither the name of Helix Re Inc. nor the
* names of its contributors may be used to endorse or promote
* products derived from this software without specific prior
* written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
* OF SUCH DAMAGE.
****************************************************************************/

#pragma once

#include <thread>

#include <pdal/Filter.hpp>

namespace pdal {

class PDAL_DLL CovarianceFeaturesFilter: public Filter
{
public:
CovarianceFeaturesFilter() : Filter() {}
CovarianceFeaturesFilter &operator=(const CovarianceFeaturesFilter &) = delete;
CovarianceFeaturesFilter(const CovarianceFeaturesFilter &) = delete;

std::string getName() const;

private:

int m_knn;
int m_threads;
std::string m_featureSet;
std::map<std::string,Dimension::Id> m_extraDims;

virtual void addDimensions(PointLayoutPtr layout);
virtual void addArgs(ProgramArgs &args);
virtual void filter(PointView &view);

void setDimensionality(PointView &view, const PointId &id, const KD3Index &kid);
};
}

24 changes: 12 additions & 12 deletions pdal/KDIndex.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -101,25 +101,25 @@ class PDAL_DLL KD2Index : public KDIndex<2>
throw pdal_error("KD2Index: point view missing 'Y' dimension.");
}

PointId neighbor(double x, double y)
PointId neighbor(double x, double y) const
{
std::vector<PointId> ids = neighbors(x, y, 1);
return (ids.size() ? ids[0] : 0);
}

PointId neighbor(PointId idx)
PointId neighbor(PointId idx) const
{
std::vector<PointId> ids = neighbors(idx, 1);
return (ids.size() ? ids[0] : 0);
}

PointId neighbor(PointRef &point)
PointId neighbor(PointRef &point) const
{
std::vector<PointId> ids = neighbors(point, 1);
return (ids.size() ? ids[0] : 0);
}

std::vector<PointId> neighbors(double x, double y, point_count_t k)
std::vector<PointId> neighbors(double x, double y, point_count_t k) const
{
k = (std::min)(m_buf.size(), k);
std::vector<PointId> output(k);
Expand All @@ -135,15 +135,15 @@ class PDAL_DLL KD2Index : public KDIndex<2>
return output;
}

std::vector<PointId> neighbors(PointId idx, point_count_t k)
std::vector<PointId> neighbors(PointId idx, point_count_t k) const
{
double x = m_buf.getFieldAs<double>(Dimension::Id::X, idx);
double y = m_buf.getFieldAs<double>(Dimension::Id::Y, idx);

return neighbors(x, y, k);
}

std::vector<PointId> neighbors(PointRef &point, point_count_t k)
std::vector<PointId> neighbors(PointRef &point, point_count_t k) const
{
double x = point.getFieldAs<double>(Dimension::Id::X);
double y = point.getFieldAs<double>(Dimension::Id::Y);
Expand Down Expand Up @@ -226,26 +226,26 @@ class PDAL_DLL KD3Index : public KDIndex<3>
throw pdal_error("KD3Index: point view missing 'Z' dimension.");
}

PointId neighbor(double x, double y, double z)
PointId neighbor(double x, double y, double z) const
{
std::vector<PointId> ids = neighbors(x, y, z, 1);
return (ids.size() ? ids[0] : 0);
}

PointId neighbor(PointId idx)
PointId neighbor(PointId idx) const
{
std::vector<PointId> ids = neighbors(idx, 1);
return (ids.size() ? ids[0] : 0);
}

PointId neighbor(PointRef &point)
PointId neighbor(PointRef &point) const
{
std::vector<PointId> ids = neighbors(point, 1);
return (ids.size() ? ids[0] : 0);
}

std::vector<PointId> neighbors(double x, double y, double z,
point_count_t k)
point_count_t k) const
{
k = (std::min)(m_buf.size(), k);
std::vector<PointId> output(k);
Expand All @@ -262,7 +262,7 @@ class PDAL_DLL KD3Index : public KDIndex<3>
return output;
}

std::vector<PointId> neighbors(PointId idx, point_count_t k)
std::vector<PointId> neighbors(PointId idx, point_count_t k) const
{
double x = m_buf.getFieldAs<double>(Dimension::Id::X, idx);
double y = m_buf.getFieldAs<double>(Dimension::Id::Y, idx);
Expand All @@ -271,7 +271,7 @@ class PDAL_DLL KD3Index : public KDIndex<3>
return neighbors(x, y, z, k);
}

std::vector<PointId> neighbors(PointRef &point, point_count_t k)
std::vector<PointId> neighbors(PointRef &point, point_count_t k) const
{
double x = point.getFieldAs<double>(Dimension::Id::X);
double y = point.getFieldAs<double>(Dimension::Id::Y);
Expand Down

0 comments on commit 9ee9211

Please sign in to comment.