Merge pull request #101 from poncateam/add_knngraphs

Add knngraphs to Ponca

CHANGELOG

@@ -9,6 +9,7 @@ improved doc.
 
 - API
     - [spatialPartitioning] Add kd-tree traits type (#80)
+    - [spatialPartitioning] Add knn-graphs (#101)
     - [fitting] Add Primitive::getNumNeighbors (#86)
     - [fitting] Add Primitive::getWeightFunc (#99)
     - [fitting] Change naming convention and rational for principal curvatures in CurvatureEstimatorBase (#94)

Ponca/SpatialPartitioning

@@ -13,3 +13,5 @@
 #include "src/SpatialPartitioning/query.h"
 #include "src/SpatialPartitioning/KdTree/kdTree.h"
 #include "src/SpatialPartitioning/KdTree/kdTreeTraits.h"
+#include "src/SpatialPartitioning/KnnGraph/knnGraph.h"
+#include "src/SpatialPartitioning/KnnGraph/knnGraphTraits.h"

Ponca/src/SpatialPartitioning/KnnGraph/Iterator/knnGraphRangeIterator.h

@@ -0,0 +1,41 @@
+/*
+ This Source Code Form is subject to the terms of the Mozilla Public
+ License, v. 2.0. If a copy of the MPL was not distributed with this
+ file, You can obtain one at http://mozilla.org/MPL/2.0/.
+*/
+
+#pragma once
+
+namespace Ponca {
+
+template <typename Traits>
+class KnnGraphRangeQuery;
+
+template <typename Traits>
+class KnnGraphRangeIterator
+{
+protected:
+    friend class KnnGraphRangeQuery<Traits>;
+
+public:
+    inline KnnGraphRangeIterator(KnnGraphRangeQuery<Traits>* query, int index = -1) : m_query(query), m_index(index) {}
+
+public:
+    bool operator != (const KnnGraphRangeIterator& other) const{
+        return m_index != other.m_index;
+    }
+
+    void operator ++ (){
+        m_query->advance(*this);
+    }
+
+    int  operator *  () const{
+        return m_index;
+    }
+
+protected:
+    KnnGraphRangeQuery<Traits>* m_query {nullptr};
+    int m_index {-1};
+};
+
+} // namespace Ponca

Ponca/src/SpatialPartitioning/KnnGraph/Query/knnGraphKNearestQuery.h

@@ -0,0 +1,54 @@
+/*
+ This Source Code Form is subject to the terms of the Mozilla Public
+ License, v. 2.0. If a copy of the MPL was not distributed with this
+ file, You can obtain one at http://mozilla.org/MPL/2.0/.
+*/
+
+#pragma once
+
+#include "../Iterator/knnGraphRangeIterator.h"
+#include <vector>
+
+namespace Ponca {
+
+template <typename Traits>class KnnGraphBase; // Need forward declaration to avoid mutual inclusion
+
+
+#ifndef PARSED_WITH_DOXYGEN
+struct KnnGraphQueryOutputType : public QueryOutputBase{
+    using OutputParameter = typename QueryOutputBase::DummyOutputParameter;
+};
+#endif
+
+template <typename Traits>class KnnGraphKNearestQuery
+#ifdef PARSED_WITH_DOXYGEN
+: public KNearestIndexQuery<typename Traits::IndexType, typename Traits::DataPoint::Scalar>
+#else
+    // we skip output because we don't need it: k is static, and already stored in the index array
+: public Query<QueryInputIsIndex<typename Traits::IndexType>,KnnGraphQueryOutputType>
+#endif
+{
+public:
+    using Iterator = typename Traits::IndexContainer::const_iterator;
+#ifdef PARSED_WITH_DOXYGEN
+    using QueryType = KNearestIndexQuery<typename Traits::IndexType, typename Traits::DataPoint::Scalar>;
+#else
+    using QueryType = Query<QueryInputIsIndex<typename Traits::IndexType>,KnnGraphQueryOutputType>;
+#endif
+
+public:
+    inline KnnGraphKNearestQuery(const KnnGraphBase<Traits>* graph, int index)
+        : m_graph(graph), QueryType(index){}
+
+    inline Iterator begin() const{
+        return m_graph->index_data().begin() + QueryType::input() * m_graph->k();
+    }
+    inline Iterator end() const{
+        return m_graph->index_data().begin() + (QueryType::input()+1) * m_graph->k();
+    }
+
+protected:
+    const KnnGraphBase<Traits>* m_graph {nullptr};
+};
+
+} // namespace Ponca

Ponca/src/SpatialPartitioning/KnnGraph/Query/knnGraphRangeQuery.h

@@ -0,0 +1,103 @@
+/*
+ This Source Code Form is subject to the terms of the Mozilla Public
+ License, v. 2.0. If a copy of the MPL was not distributed with this
+ file, You can obtain one at http://mozilla.org/MPL/2.0/.
+*/
+
+#pragma once
+
+#include "../../query.h"
+#include "../Iterator/knnGraphRangeIterator.h"
+
+#include <vector>
+#include <stack>
+
+namespace Ponca {
+template <typename Traits> class KnnGraphBase;
+
+template <typename Traits>
+class KnnGraphRangeQuery : public RangeIndexQuery<typename Traits::IndexType, typename Traits::DataPoint::Scalar>
+{
+protected:
+    using QueryType = RangeIndexQuery<typename Traits::IndexType, typename Traits::DataPoint::Scalar>;
+    friend class KnnGraphRangeIterator<Traits>; // This type must be equal to KnnGraphRangeQuery::Iterator
+
+public:
+    using DataPoint  = typename Traits::DataPoint;
+    using IndexType  = typename Traits::IndexType;
+    using Scalar     = typename DataPoint::Scalar;
+    using VectorType = typename DataPoint::VectorType;
+    using Iterator   = KnnGraphRangeIterator<Traits>;
+
+public:
+    inline KnnGraphRangeQuery(const KnnGraphBase<Traits>* graph, Scalar radius, int index):
+            QueryType(radius, index),
+            m_graph(graph),
+            m_flag(graph->size()),
+            m_stack() {}
+
+public:
+    inline Iterator begin(){
+        KnnGraphRangeIterator it(this);
+        this->initialize(it);
+        this->advance(it);
+        return it;
+    }
+
+    inline Iterator end(){
+        return KnnGraphRangeIterator(this, m_graph->size());
+    }
+
+protected:
+    inline void initialize(Iterator& iterator){
+        m_flag.resize(m_graph->size());
+        std::fill(m_flag.begin(), m_flag.end(), false);
+
+        PONCA_DEBUG_ASSERT(m_stack.empty());
+        m_stack.push(QueryType::input());
+        m_flag[QueryType::input()] = true;
+
+        iterator.m_index = -1;
+    }
+
+    inline void advance(Iterator& iterator){
+        const auto& points  = m_graph->m_kdTreePoints;
+        const auto& point   = points[QueryType::input()].pos();
+
+        if(! (iterator != end())) return;
+
+        if(m_stack.empty())
+        {
+            iterator = end();
+        }
+        else
+        {
+            int idx_current = m_stack.top();
+            m_stack.pop();
+
+            PONCA_DEBUG_ASSERT((point - points[idx_current]).squaredNorm() < m_squared_radius);
+
+            iterator.m_index = idx_current;
+
+            for(int idx_nei : m_graph->k_nearest_neighbors(idx_current))
+            {
+                PONCA_DEBUG_ASSERT(idx_nei>0);
+                Scalar d  = (point - points[idx_nei].pos()).squaredNorm();
+                Scalar th = QueryType::descentDistanceThreshold();
+                if(!m_flag[idx_nei] && (point - points[idx_nei].pos()).squaredNorm() < QueryType::descentDistanceThreshold())
+                {
+                    m_flag[idx_nei] = true;
+                    m_stack.push(idx_nei);
+                }
+            }
+            if (iterator.m_index == QueryType::input()) advance(iterator); // query is not included in returned set
+        }
+    }
+
+protected:
+    const KnnGraphBase<Traits>*   m_graph {nullptr};
+    std::vector<bool> m_flag;  ///< hold ids status (ids range from 0 to point cloud size)
+    std::stack<int>   m_stack; ///< hold ids (ids range from 0 to point cloud size)
+};
+
+} // namespace Ponca

Ponca/src/SpatialPartitioning/KnnGraph/knnGraph.h

@@ -0,0 +1,134 @@
+/*
+ This Source Code Form is subject to the terms of the Mozilla Public
+ License, v. 2.0. If a copy of the MPL was not distributed with this
+ file, You can obtain one at http://mozilla.org/MPL/2.0/.
+*/
+
+#pragma once
+
+#include "./knnGraphTraits.h"
+
+#include "Query/knnGraphKNearestQuery.h"
+#include "Query/knnGraphRangeQuery.h"
+
+#include "../KdTree/kdTree.h"
+
+#include <memory>
+
+namespace Ponca {
+
+template <typename Traits> class KnnGraphBase;
+
+/*!
+ * \brief Public interface for KnnGraph datastructure.
+ *
+ * Provides default implementation of the KnnGraph
+ *
+ * \see KnnGraphDefaultTraits for the default trait interface documentation.
+ * \see KnnGraphBase for complete API
+ */
+template <typename DataPoint>
+using KnnGraph = KnnGraphBase<KnnGraphDefaultTraits<DataPoint>>;
+
+/*!
+ * \brief Customizable base class for KnnGraph datastructure
+ *
+ * \see Ponca::KnGraph
+ *
+ * \tparam Traits Traits type providing the types and constants used by the KnnGraph. Must have the
+ * same interface as the default traits type.
+ *
+ * \see KnnGraphDefaultTraits for the trait interface documentation.
+ *
+ */
+template <typename Traits> class KnnGraphBase
+{
+public:
+    using DataPoint  = typename Traits::DataPoint; ///< DataPoint given by user via Traits
+    using Scalar     = typename DataPoint::Scalar; ///< Scalar given by user via DataPoint
+    using VectorType = typename DataPoint::VectorType; ///< VectorType given by user via DataPoint
+
+    using IndexType      = typename Traits::IndexType;
+    using PointContainer = typename Traits::PointContainer; ///< Container for DataPoint used inside the KdTree
+    using IndexContainer = typename Traits::IndexContainer; ///< Container for indices used inside the KdTree
+
+    using KNearestIndexQuery = KnnGraphKNearestQuery<Traits>;
+    using RangeIndexQuery    = KnnGraphRangeQuery<Traits>;
+
+    friend class KnnGraphKNearestQuery<Traits>; // This type must be equal to KnnGraphBase::KNearestIndexQuery
+    friend class KnnGraphRangeQuery<Traits>;    // This type must be equal to KnnGraphBase::RangeIndexQuery
+
+    // knnGraph ----------------------------------------------------------------
+public:
+    /// \brief Build a KnnGraph from a KdTree
+    ///
+    /// \warning In the current version, the graph does not support kdtree with subsampling
+    /// \param k Number of requested neighbors. Might be reduced if k is larger than the kdtree size - 1
+    ///          (query point is not included in query output, thus -1)
+    ///
+    /// \warning Stores a const reference to kdtree.point_data()
+    /// \warning KdTreeTraits compatibility is checked with static assertion
+    template<typename KdTreeTraits>
+    inline KnnGraphBase(const KdTreeBase<KdTreeTraits>& kdtree, int k = 6)
+            : m_k(std::min(k,kdtree.index_count()-1)),
+              m_kdTreePoints(kdtree.point_data())
+    {
+        static_assert( std::is_same_v<typename Traits::DataPoint, typename KdTreeTraits::DataPoint>,
+                       "KdTreeTraits::DataPoint is not equal to Traits::DataPoint" );
+        static_assert( std::is_same_v<typename Traits::PointContainer, typename KdTreeTraits::PointContainer>,
+                       "KdTreeTraits::PointContainer is not equal to Traits::PointContainer" );
+        static_assert( std::is_same_v<typename Traits::IndexContainer, typename KdTreeTraits::IndexContainer>,
+                       "KdTreeTraits::IndexContainer is not equal to Traits::IndexContainer" );
+
+        // We need to account for the entire point set, irrespectively of the sampling. This is because the kdtree
+        // (k_nearest_neighbors) return ids of the entire point set, not it sub-sampled list of ids.
+        // \fixme Update API to properly handle kdtree subsampling
+        const int cloudSize   = kdtree.point_count();
+        {
+            const int samplesSize = kdtree.index_count();
+            eigen_assert(cloudSize == samplesSize);
+        }
+
+        m_indices.resize(cloudSize * m_k, -1);
+
+#pragma omp parallel for shared(kdtree, cloudSize) default(none)
+        for(int i=0; i<cloudSize; ++i)
+        {
+            int j = 0;
+            for(auto n : kdtree.k_nearest_neighbors(typename KdTreeTraits::IndexType(i),
+                                                   typename KdTreeTraits::IndexType(m_k)))
+            {
+                m_indices[i * m_k + j] = n;
+                ++j;
+            }
+        }
+    }
+
+    // Query -------------------------------------------------------------------
+public:
+    inline KNearestIndexQuery k_nearest_neighbors(int index) const{
+        return KnnGraphKNearestQuery(this, index);
+    }
+
+    inline RangeIndexQuery    range_neighbors(int index, Scalar r) const{
+        return RangeIndexQuery(this, r, index);
+    }
+
+    // Accessors ---------------------------------------------------------------
+public:
+    /// \brief Number of neighbor per vertex
+    inline int k() const { return m_k; }
+    /// \brief Number of vertices in the neighborhood graph
+    inline int size() const { return m_indices.size()/m_k; }
+
+    // Data --------------------------------------------------------------------
+private:
+    const int m_k;
+    IndexContainer m_indices; ///< \brief Stores neighborhood relations
+
+protected: // for friends relations
+    const PointContainer& m_kdTreePoints;
+    inline const IndexContainer& index_data() const { return m_indices; };
+};
+
+} // namespace Ponca

Ponca/src/SpatialPartitioning/KnnGraph/knnGraphTraits.h

@@ -0,0 +1,48 @@
+/*
+ This Source Code Form is subject to the terms of the Mozilla Public
+ License, v. 2.0. If a copy of the MPL was not distributed with this
+ file, You can obtain one at http://mozilla.org/MPL/2.0/.
+*/
+
+
+#pragma once
+
+#include <Eigen/Geometry>
+
+namespace Ponca {
+
+/*!
+ * \brief The default traits type used by the kd-tree.
+ */
+template <typename _DataPoint>
+struct KnnGraphDefaultTraits
+{
+    /*!
+     * \brief The type used to store point data.
+     *
+     * Must provide `Scalar` and `VectorType` typedefs.
+     *
+     * `VectorType` must provide a `squaredNorm()` function returning a `Scalar`, as well as a
+     * `maxCoeff(int*)` function returning the dimension index of its largest scalar in its output
+     * parameter (e.g. 0 for *x*, 1 for *y*, etc.).
+     */
+    using DataPoint = _DataPoint;
+
+private:
+    using Scalar     = typename DataPoint::Scalar;
+    using VectorType = typename DataPoint::VectorType;
+
+public:
+    /*!
+     * \brief The type used to calculate node bounding boxes.
+     *
+     * Must provide `min()`, `max()`, and `center()` functions, all returning a `VectorType`.
+     */
+    using AabbType = Eigen::AlignedBox<Scalar, DataPoint::Dim>;
+
+    // Containers
+    using IndexType      = int;
+    using PointContainer = std::vector<DataPoint>;
+    using IndexContainer = std::vector<IndexType>;
+};
+} // namespace Ponca