kdtree.hpp

/* $Id$ */

#ifndef SAMS_KDTREE_HPP
#define SAMS_KDTREE_HPP

//
// (c) Matthew B. Kennel, Institute for Nonlinear Science, UCSD (2004)
//
// Licensed under the Academic Free License version 1.1 found at the end of this file

//
// Modified by Maxime Cottret 
//

//
// Implement a kd tree for fast searching of points in a fixed data base
// in k-dimensional Euclidean space.
//
// 

#include <vector>
#include <algorithm>
#include "jmath/jblas.hpp"

using namespace std;

namespace jafar {
  namespace sams {

    typedef struct {
      double lower, upper;
    } interval;


    // let the compiler know that this is a names of classes. 
    class Kdtree_node; 
    class searchrecord;

    //
    // struct KDTREE_RESULT
    // class  KDTREE_RESULT
    //

    struct Kdtree_result {
      // 
      // the search routines return a (wrapped) vector
      // of these. 
      //
    public:
      double dis;  
      int idx;    
      
      
      bool operator<(const Kdtree_result& e1) const{
  return (dis < e1.dis);
      }

    }; 
    

    class Kdtree_result_vector : public vector<Kdtree_result> {
      // inherit a vector<Kdtree_result>
      // but, optionally maintain it in heap form as a priority
      // queue.
    public:

    
      void push_element_and_heapify(Kdtree_result& e);
      double replace_maxpri_elt_return_new_maxpri(Kdtree_result& e);


     double max_value(); 
      
    };


    class Kdtree {
    public: 
      const jblas::mat& the_data;   

      const int N;   
      int dim; 
      bool sort_results;  
      const bool rearrange;  

    public:
      
      Kdtree(jblas::mat& data_in,bool sort_out = true, bool rearrange_in = true,int dim_in=-1);

      ~Kdtree();
  

    public:
      // search routines

      void n_nearest_brute_force(jblas::vec& qv, int nn, Kdtree_result_vector& result);
      
      void n_nearest(jblas::vec& qv, int nn, Kdtree_result_vector& result);
      
      void n_nearest_around_point(int idxin, int correltime, int nn,
          Kdtree_result_vector& result);
       
      void r_nearest(jblas::vec& qv, double r2,Kdtree_result_vector& result); 
       
      void r_nearest_around_point(int idxin, int correltime, double r2,
          Kdtree_result_vector& result);
     
      int r_count(jblas::vec& qv, double r2);
      
      int r_count_around_point(int idxin, int correltime, double r2);
     

      friend class Kdtree_node;
      friend class searchrecord;
    private:
      // private data members

      Kdtree_node* root; // the root pointer

      const jblas::mat* data;
      // pointing either to the_data or an internal
      // rearranged data as necessary

      vector<int> ind; 
      // the index for the tree leaves.  Data in a leaf with bounds [l,u] are
      // in  'the_data[ind[l],*] to the_data[ind[u],*]

      jblas::mat rearranged_data;  
      // if rearrange is true then this is the rearranged data storage. 


      static const int bucketsize = 12;  // global constant. 

    private:
      void set_data(jblas::mat& din); 
      void build_tree(); // builds the tree.  Used upon construction. 
      Kdtree_node* build_tree_for_range(int l, int u, Kdtree_node* parent);
      void select_on_coordinate(int c, int k, int l, int u); 
      int select_on_coordinate_value(int c, double alpha, int l, int u); 
      void spread_in_coordinate(int c, int l, int u, interval& interv);
    };

#ifndef SWIG

    //
    // search record substructure
    //
    // one of these is created for each search.
    // this holds useful information  to be used
    // during the search



    static const double infinity = 1.0e38;

    class searchrecord {

    private:
      friend class Kdtree;
      friend class Kdtree_node;

      jblas::vec& qv; 
      int dim;
      bool rearrange;
      unsigned int nn; // , nfound;
      double ballsize;
      int centeridx, correltime;

      Kdtree_result_vector& result;  // results
      const jblas::mat* data; 
      const vector<int>& ind; 
      // constructor

    public:
      searchrecord(jblas::vec& qv_in, Kdtree& tree_in,
       Kdtree_result_vector& result_in) :  
  qv(qv_in),
  result(result_in),
  data(tree_in.data),
  ind(tree_in.ind) 
      {
  dim = tree_in.dim;
  rearrange = tree_in.rearrange;
  ballsize = infinity; 
  nn = 0; 
      };

    };



    //
    // class KDTREE_NODE
    // 
    // a node in the tree.  Many are created per tree dynamically.. 
    //


    class Kdtree_node {
    public:
      // constructor
      Kdtree_node(int dim);
      //, int cut_dim_in,
      //         double cut_val_in, double cut_val_left_in, 
      //         double cut_val_right_in);
      // destructor
      ~Kdtree_node();

    private:
      // visible to self and Kdtree.
      friend class Kdtree;  // allow Kdtree to access private 

      int cut_dim;                                 // dimension to cut; 
      double cut_val, cut_val_left, cut_val_right;  //cut value
      int l,u;  // extents in index array for searching

      vector<interval> box; // [min,max] of the box enclosing all points
  
      Kdtree_node *left, *right;  // pointers to left and right nodes. 

      void search(searchrecord& sr); 
      // recursive innermost core routine for searching.. 

      bool box_in_search_range(searchrecord& sr);
      // return true if the bounding box for this node is within the
      // search range given by the searchvector and maximum ballsize in 'sr'. 

      void check_query_in_bound(searchrecord& sr); // debugging only

      // for processing final buckets. 
      void process_terminal_node(searchrecord& sr);
      void process_terminal_node_fixedball(searchrecord& sr);


    };

#endif // SWIG

  } // namespace sams
} // namespace jafar

#endif // SAMS_KDTREE_HPP


// The KDTREE2 software is licensed under the terms of the Academic Free
// Software License, listed herein.  In addition, users of this software
// must give appropriate citation in relevant technical documentation or
// journal paper to the author, Matthew B. Kennel, Institute For
// Nonlinear Science, preferably via a reference to the www.arxiv.org
// repository of this document, {\tt www.arxiv.org e-print:
// physics/0408067}.  This requirement will be deemed to be advisory and
// not mandatory as is necessary to permit the free inclusion of the
// present software with any software licensed under the terms of any
// version of the GNU General Public License, or GNU Library General
// Public License.


// Academic Free License
// Version 1.1

// This Academic Free License applies to any original work of authorship
// (the "Original Work") whose owner (the "Licensor") has placed the
// following notice immediately following the copyright notice for the
// Original Work: "Licensed under the Academic Free License version 1.1."

// Grant of License. Licensor hereby grants to any person obtaining a
// copy of the Original Work ("You") a world-wide, royalty-free,
// non-exclusive, perpetual, non-sublicenseable license (1) to use, copy,
// modify, merge, publish, perform, distribute and/or sell copies of the
// Original Work and derivative works thereof, and (2) under patent
// claims owned or controlled by the Licensor that are embodied in the
// Original Work as furnished by the Licensor, to make, use, sell and
// offer for sale the Original Work and derivative works thereof, subject
// to the following conditions.

// Right of Attribution. Redistributions of the Original Work must
// reproduce all copyright notices in the Original Work as furnished by
// the Licensor, both in the Original Work itself and in any
// documentation and/or other materials provided with the distribution of
// the Original Work in executable form.

// Exclusions from License Grant. Neither the names of Licensor, nor the
// names of any contributors to the Original Work, nor any of their
// trademarks or service marks, may be used to endorse or promote
// products derived from this Original Work without express prior written
// permission of the Licensor.

// WARRANTY AND DISCLAIMERS. LICENSOR WARRANTS THAT THE COPYRIGHT IN AND
// TO THE ORIGINAL WORK IS OWNED BY THE LICENSOR OR THAT THE ORIGINAL
// WORK IS DISTRIBUTED BY LICENSOR UNDER A VALID CURRENT LICENSE FROM THE
// COPYRIGHT OWNER. EXCEPT AS EXPRESSLY STATED IN THE IMMEDIATELY
// PRECEEDING SENTENCE, THE ORIGINAL WORK IS PROVIDED UNDER THIS LICENSE
// ON AN "AS IS" BASIS, WITHOUT WARRANTY, EITHER EXPRESS OR IMPLIED,
// INCLUDING, WITHOUT LIMITATION, THE WARRANTY OF NON-INFRINGEMENT AND
// WARRANTIES THAT THE ORIGINAL WORK IS MERCHANTABLE OR FIT FOR A
// PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY OF THE ORIGINAL
// WORK IS WITH YOU. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL
// PART OF THIS LICENSE. NO LICENSE TO ORIGINAL WORK IS GRANTED HEREUNDER
// EXCEPT UNDER THIS DISCLAIMER.

// LIMITATION OF LIABILITY. UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL
// THEORY, WHETHER TORT (INCLUDING NEGLIGENCE), CONTRACT, OR OTHERWISE,
// SHALL THE LICENSOR BE LIABLE TO ANY PERSON FOR ANY DIRECT, INDIRECT,
// SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER ARISING
// AS A RESULT OF THIS LICENSE OR THE USE OF THE ORIGINAL WORK INCLUDING,
// WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK STOPPAGE,
// COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL OTHER COMMERCIAL
// DAMAGES OR LOSSES, EVEN IF SUCH PERSON SHALL HAVE BEEN INFORMED OF THE
// POSSIBILITY OF SUCH DAMAGES. THIS LIMITATION OF LIABILITY SHALL NOT
// APPLY TO LIABILITY FOR DEATH OR PERSONAL INJURY RESULTING FROM SUCH
// PARTY'S NEGLIGENCE TO THE EXTENT APPLICABLE LAW PROHIBITS SUCH
// LIMITATION. SOME JURISDICTIONS DO NOT ALLOW THE EXCLUSION OR
// LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THIS EXCLUSION
// AND LIMITATION MAY NOT APPLY TO YOU.

// License to Source Code. The term "Source Code" means the preferred
// form of the Original Work for making modifications to it and all
// available documentation describing how to access and modify the
// Original Work. Licensor hereby agrees to provide a machine-readable
// copy of the Source Code of the Original Work along with each copy of
// the Original Work that Licensor distributes. Licensor reserves the
// right to satisfy this obligation by placing a machine-readable copy of
// the Source Code in an information repository reasonably calculated to
// permit inexpensive and convenient access by You for as long as
// Licensor continues to distribute the Original Work, and by publishing
// the address of that information repository in a notice immediately
// following the copyright notice that applies to the Original Work.

// Mutual Termination for Patent Action. This License shall terminate
// automatically and You may no longer exercise any of the rights granted
// to You by this License if You file a lawsuit in any court alleging
// that any OSI Certified open source software that is licensed under any
// license containing this "Mutual Termination for Patent Action" clause
// infringes any patent claims that are essential to use that software.

// This license is Copyright (C) 2002 Lawrence E. Rosen. All rights
// reserved. Permission is hereby granted to copy and distribute this
// license without modification. This license may not be modified without
// the express written permission of its copyright owner.