Listing 1: Template for the nearest neighbor search

#if !defined(TNEAR_H_INCLUDED)
#define TNEAR_H_INCLUDED


#include <limits.h>
#include <float.h>
#include <math.h>
#include <vector>

//==============================================================
template <typename T> class CNearTree
{
    T * m_ptLeft; // left object stored in this node
    T * m_ptRight;// right object stored in this node
    double m_dMaxLeft; //longest distance from the left object
                       // to anything below it in the tree
    double m_dMaxRight;// longest distance from the right object
                       // to anything below it in the tree
    CNearTree * m_pLeftBranch;  //tree descending from the left
    CNearTree * m_pRightBranch; //tree descending from the right

public:

//==============================================================
   CNearTree(void)  // constructor
   {
      m_ptLeft       = 0;
      m_ptRight      = 0;
      m_pLeftBranch  = 0;
      m_pRightBranch = 0;
      m_dMaxLeft     = DBL_MIN;
      m_dMaxRight    = DBL_MIN;
   }  //  CNearTree constructor

//==============================================================
   ~CNearTree(void)  // destructor
   {
      delete m_pLeftBranch  ;  m_pLeftBranch  =0;
      delete m_pRightBranch ;  m_pRightBranch =0;
      delete m_ptLeft       ;  m_ptLeft       =0;
      delete m_ptRight      ;  m_ptRight      =0;

      m_dMaxLeft     = DBL_MIN;
      m_dMaxRight    = DBL_MIN;
   }  //  ~CNearTree

//==============================================================
   void m_fnInsert( const T& t )
   {
      // do a bit of precomputing if possible so that we can
      // reduce the number of calls to operator 'double' as much
      // as possible; 'double' might use square roots
      double dTempRight =  0;
      double dTempLeft  =  0;
      if ( m_ptRight  != 0 )
      {
         dTempRight  = fabs( double( t - *m_ptRight ));
         dTempLeft   = fabs( double( t - *m_ptLeft  ));
      }
      if ( m_ptLeft == 0 )
      {
         m_ptLeft = new T( t );
      }
      else if ( m_ptRight == 0 )
      {
         m_ptRight   = new T( t );
      }
      else if ( dTempLeft > dTempRight )
      {
         if (m_pRightBranch==0) m_pRightBranch=new CNearTree;
         // note that the next line assumes that m_dMaxRight
         // is negative for a new node
         if (m_dMaxRight<dTempRight) m_dMaxRight=dTempRight;
         m_pRightBranch->m_fnInsert( t );
      }
      else
      {
         if (m_pLeftBranch ==0) m_pLeftBranch=new CNearTree;
         // note that the next line assumes that m_dMaxLeft
         // is negative for a new node
         if (m_dMaxLeft<dTempLeft) m_dMaxLeft=dTempLeft;
         m_pLeftBranch->m_fnInsert( t );
      }
   }  //  m_fnInsert

//==============================================================
   bool m_bfnNearestNeighbor ( const double& dRadius,
                             T& tClosest,   const T& t ) const
   {
      double dSearchRadius = dRadius;
      return ( m_bfnNearest ( dSearchRadius, tClosest, t ) );
   }  //  m_bfnNearestNeighbor

//==============================================================
   bool m_bfnFarthestNeighbor ( T& tFarthest, const T& t ) const
   {
      double dSearchRadius = DBL_MIN;
      return ( m_bfnFindFarthest (dSearchRadius, tFarthest, t));
   }  //  m_bfnFarthestNeighbor

//==============================================================
   long m_lfnFindInSphere ( const double& dRadius,
               std::vector<  T >& tClosest,   const T& t ) const
   {   // t is the probe point, tClosest is a vector of contacts
      // clear the contents of the return vector so that
      // things don't accidentally accumulate
      tClosest.clear( );
      return ( m_lfnInSphere( dRadius, tClosest, t ) );
   }  //  m_lfnFindInSphere

private:

//==============================================================
   bool m_bfnNearest ( double& dRadius,
                       T& tClosest,   const T& t ) const
   {
      double   dTempRadius;
      bool  bRet = false;
      // first test each of the left and right positions to see
      // if one holds a point nearer than the nearest so far.
      if (( m_ptLeft!=0 ) &&
        ((dTempRadius = fabs(double(t-*m_ptLeft))) <= dRadius))
      {
         dRadius  = dTempRadius;
         tClosest = *m_ptLeft;
         bRet     = true;
      }
      if (( m_ptRight!=0 ) &&
        (( dTempRadius = fabs(double(t-*m_ptRight)))<=dRadius))
      {
         dRadius  = dTempRadius;
         tClosest = *m_ptRight;
         bRet     = true;
      }
      // Now we test to see if the branches below might hold an
      // object nearer than the best so far found. The triangle
      // rule is used to test whether it's even necessary to
      // descend.
      if (( m_pLeftBranch  != 0 )  &&
        ((dRadius+m_dMaxLeft) >= fabs(double(t-*m_ptLeft))))
      {
         bRet|=m_pLeftBranch->m_bfnNearest(dRadius,tClosest,t);
      }

      if (( m_pRightBranch != 0 )  &&
        ((dRadius+m_dMaxRight) >= fabs(double(t-*m_ptRight))))
      {
         bRet|=m_pRightBranch->m_bfnNearest(dRadius,tClosest,t);
      }
      return ( bRet );
   };   // m_bfnNearest

//==============================================================
   long m_lfnInSphere ( const double& dRadius,
               std::vector<  T >& tClosest,   const T& t ) const
   {   // t is the probe point, tClosest is a vector of contacts
      long lReturn = 0;
      // first test each of the left and right positions to see
      // if one holds a point nearer than the search radius.
      if ((m_ptLeft!=0) && (fabs(double(t-*m_ptLeft))<=dRadius))
      {
         tClosest.push_back( *m_ptLeft ); // It's a keeper
         lReturn++ ;
      }
      if ((m_ptRight!=0)&&(fabs(double(t-*m_ptRight))<=dRadius))
      {
         tClosest.push_back( *m_ptRight ); // It's a keeper
         lReturn++ ;
      }
      //
      // Now we test to see if the branches below might hold an
      // object nearer than the search radius. The triangle rule
      // is used to test whether it's even necessary to descend.
      //
      if (( m_pLeftBranch  != 0 )  &&
         ((dRadius+m_dMaxLeft) >= fabs(double(t-*m_ptLeft))))
      {
         lReturn +=
           m_pLeftBranch->m_lfnInSphere( dRadius, tClosest, t );
      }
      if (( m_pRightBranch != 0 )  &&
         ((dRadius+m_dMaxRight) >= fabs(double(t-*m_ptRight))))
      {
         lReturn +=
           m_pRightBranch->m_lfnInSphere( dRadius, tClosest, t );
      }
      return ( lReturn );
   }  //  m_lfnInSphere

//==============================================================
   bool m_bfnFindFarthest ( double& dRad,
                        T& tFarthest,   const T& t ) const
   {
//   deleted from the journal listing since it is quite similar
//   to nearest
        return ( false );
   };   // m_bfnFindFarthest


}; // template class TNear

#endif // !defined(TNEAR_H_INCLUDED)
— End of Listing —