KDTree.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/KDTree.h - Released 2024-06-18T15:48:54Z
// ----------------------------------------------------------------------------
// This file is part of the PixInsight Class Library (PCL).
// PCL is a multiplatform C++ framework for development of PixInsight modules.
//
// Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
//
// Redistribution and use in both source and binary forms, with or without
// modification, is permitted provided that the following conditions are met:
//
// 1. All redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. All 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.
//
// 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
//    of their contributors, may be used to endorse or promote products derived
//    from this software without specific prior written permission. For written
//    permission, please contact info@pixinsight.com.
//
// 4. All products derived from this software, in any form whatsoever, must
//    reproduce the following acknowledgment in the end-user documentation
//    and/or other materials provided with the product:
//
//    "This product is based on software from the PixInsight project, developed
//    by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
//
//    Alternatively, if that is where third-party acknowledgments normally
//    appear, this acknowledgment must be reproduced in the product itself.
//
// THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS 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 PLEIADES ASTROPHOTO OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
// INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
// DATA OR PROFITS) 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.
// ----------------------------------------------------------------------------
 
#ifndef __PCL_KDTree_h
#define __PCL_KDTree_h
 
 
#include <pcl/Defs.h>
 
#include <pcl/Array.h>
#include <pcl/Vector.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
template <class T>
class PCL_CLASS KDTree
{
public:
 
   using point = T;
 
   using component = typename point::component;
 
   using component_vector = GenericVector<component>;
 
   using point_list = Array<point>;
 
   KDTree() = default;
 
   KDTree( const point_list& points, int bucketCapacity = 16 )
   {
      Build( points, bucketCapacity );
   }
 
   KDTree( const point_list& points, int dimension, int bucketCapacity )
   {
      Build( points, dimension, bucketCapacity );
   }
 
   KDTree( const KDTree& ) = delete;
 
   KDTree& operator =( const KDTree& ) = delete;
 
   KDTree( KDTree&& x )
      : m_root( x.m_root )
      , m_dimension( x.m_dimension )
      , m_bucketCapacity( x.m_bucketCapacity )
      , m_length( x.m_length )
   {
      x.m_root = nullptr;
      x.m_length = 0;
   }
 
   KDTree& operator =( KDTree&& x )
   {
      if ( &x != this )
      {
         DestroyTree( m_root );
         m_root = x.m_root;
         m_dimension = x.m_dimension;
         m_bucketCapacity = x.m_bucketCapacity;
         m_length = x.m_length;
         x.m_root = nullptr;
         x.m_length = 0;
      }
      return *this;
   }
 
   ~KDTree()
   {
      Clear();
   }
 
   void Clear()
   {
      DestroyTree( m_root );
      m_root = nullptr;
      m_length = 0;
   }
 
   void Build( const point_list& points, int bucketCapacity = 16 )
   {
      Clear();
      m_bucketCapacity = Max( 1, bucketCapacity );
      if ( !points.IsEmpty() )
      {
         m_dimension = points[0].Length();
         if ( m_dimension < 1 )
            throw Error( "KDTree::Build(): Invalid point space dimension." );
         m_root = BuildTree( points, 0 );
      }
   }
 
   void Build( const point_list& points, int dimension, int bucketCapacity )
   {
      Clear();
      m_bucketCapacity = Max( 1, bucketCapacity );
      if ( (m_dimension = dimension) < 1 )
         throw Error( "KDTree::Build(): Invalid point space dimension." );
      m_root = BuildTree( points, 0 );
   }
 
   point_list Search( const point& pt, component epsilon ) const
   {
      component_vector p0( m_dimension );
      component_vector p1( m_dimension );
      for ( int i = 0; i < m_dimension; ++i )
      {
         p0[i] = pt[i] - epsilon;
         p1[i] = pt[i] + epsilon;
      }
      point_list found;
      SearchTree( found, p0, p1, m_root, 0 );
      return found;
   }
 
   template <class F>
   void Search( const point& pt, component epsilon, F callback, void* data ) const
   {
      component_vector p0( m_dimension );
      component_vector p1( m_dimension );
      for ( int i = 0; i < m_dimension; ++i )
      {
         p0[i] = pt[i] - epsilon;
         p1[i] = pt[i] + epsilon;
      }
      SearchTree( p0, p1, callback, data, m_root, 0 );
   }
 
   int Dimension() const
   {
      return m_dimension;
   }
 
   size_type Length() const
   {
      return m_length;
   }
 
   bool IsEmpty()
   {
      return m_root == nullptr;
   }
 
   friend void Swap( KDTree& x1, KDTree& x2 )
   {
      pcl::Swap( x1.m_root,           x2.m_root );
      pcl::Swap( x1.m_dimension,      x2.m_dimension );
      pcl::Swap( x1.m_bucketCapacity, x2.m_bucketCapacity );
      pcl::Swap( x1.m_length,         x2.m_length );
   }
 
private:
 
   struct Node
   {
      double split = 0;       // position of this node's splitting hyperplane
      Node*  left = nullptr;  // child points at coordinates <= split
      Node*  right = nullptr; // child points at coordinates > split
 
      Node( double s = 0 )
         : split( s )
      {
      }
 
      bool IsLeaf() const
      {
         return left == nullptr && right == nullptr;
      }
   };
 
   struct LeafNode : public Node
   {
      point_list points;
 
      LeafNode( const point_list& p )
         : points( p )
      {
      }
   };
 
   Node*     m_root = nullptr;
   int       m_dimension = 0;
   int       m_bucketCapacity = 0;
   size_type m_length = 0;
 
   LeafNode* NewLeafNode( const point_list& points )
   {
      m_length += points.Length();
      return new LeafNode( points );
   }
 
   Node* BuildTree( const point_list& points, int depth )
   {
      if ( points.IsEmpty() )
         return nullptr;
 
      if ( points.Length() <= size_type( m_bucketCapacity ) )
         return NewLeafNode( points );
 
      int index = depth % m_dimension;
 
      Node* node = new Node( SplitValue( points, index ) );
 
      point_list left, right;
      for ( const point& p : points )
         if ( p[index] <= node->split )
            left.Add( p );
         else
            right.Add( p );
 
      // If we are about to build a degenerate subtree, abort this branch of
      // recursion right now.
      if ( left.IsEmpty() || right.IsEmpty() )
      {
         delete node;
         return NewLeafNode( points );
      }
 
      node->left  = BuildTree( left, depth+1 );
      node->right = BuildTree( right, depth+1 );
 
      // Further degeneracies cannot happen in theory, but let's prevent them
      // for extra safety.
      if ( node->IsLeaf() )
      {
         delete node;
         return NewLeafNode( points );
      }
 
      return node;
   }
 
   void SearchTree( point_list& found, const component_vector& p0, const component_vector& p1, const Node* node, int depth ) const
   {
      if ( node != nullptr )
         if ( node->IsLeaf() )
         {
            const LeafNode* leaf = static_cast<const LeafNode*>( node );
            for ( const point& p : leaf->points )
               for ( int j = 0; ; )
               {
                  component x = p[j];
                  if ( x < p0[j] || p1[j] < x )
                     break;
                  if ( ++j == m_dimension )
                  {
                     found.Add( p );
                     break;
                  }
               }
         }
         else
         {
            int index = depth % m_dimension;
            if ( p0[index] <= node->split )
               SearchTree( found, p0, p1, node->left, depth+1 );
            if ( p1[index] > node->split )
               SearchTree( found, p0, p1, node->right, depth+1 );
         }
   }
 
   template <class F>
   void SearchTree( const component_vector& p0, const component_vector& p1, F callback, void* data, const Node* node, int depth ) const
   {
      if ( node != nullptr )
         if ( node->IsLeaf() )
         {
            const LeafNode* leaf = static_cast<const LeafNode*>( node );
            for ( const point& p : leaf->points )
               for ( int j = 0; ; )
               {
                  component x = p[j];
                  if ( x < p0[j] || p1[j] < x )
                     break;
                  if ( ++j == m_dimension )
                  {
                     callback( p, data );
                     break;
                  }
               }
         }
         else
         {
            int index = depth % m_dimension;
            if ( p0[index] <= node->split )
               SearchTree( p0, p1, callback, data, node->left, depth+1 );
            if ( p1[index] > node->split )
               SearchTree( p0, p1, callback, data, node->right, depth+1 );
         }
   }
 
   void DestroyTree( Node* node )
   {
      if ( node != nullptr )
         if ( node->IsLeaf() )
            delete static_cast<LeafNode*>( node );
         else
         {
            DestroyTree( node->left );
            DestroyTree( node->right );
            delete node;
         }
   }
 
   double SplitValue( const point_list& points, int index )
   {
      component_vector v( points.Length() );
      for ( int i = 0; i < v.Length(); ++i )
         v[i] = points[i][index];
      return v.Median();
   }
};
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif   // __PCL_KDTree_h
 
// ----------------------------------------------------------------------------
// EOF pcl/KDTree.h - Released 2024-06-18T15:48:54Z