contrib/mul/mbl/mbl_cluster_tree.txx
Go to the documentation of this file.
00001 // This is mul/mbl/mbl_cluster_tree.txx
00002 #ifndef mbl_cluster_txx_
00003 #define mbl_cluster_txx_
00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00005 #pragma interface
00006 #endif
00007 //:
00008 // \file
00009 // \brief  Record trees of clusters of data, for faster neighbour finding
00010 // \author Tim Cootes
00011 
00012 #include "mbl_cluster_tree.h"
00013 #include <vcl_iostream.h>
00014 #include <vcl_cassert.h>
00015 #include <vsl/vsl_binary_io.h>
00016 #include <vsl/vsl_vector_io.h>
00017 
00018 //: Default constructor
00019 template<class T, class D>
00020 mbl_cluster_tree<T,D>::mbl_cluster_tree()
00021   : data_(0)
00022 {
00023 }
00024 
00025 //: Empty clusters
00026 template<class T, class D>
00027 void mbl_cluster_tree<T,D>::empty()
00028 {
00029   for (unsigned i=0;i<cluster_.size();++i)
00030   {
00031     cluster_[i].empty();
00032     parent_[i].resize(0);
00033   }
00034 }
00035 
00036 //: Define number of levels and max radius of clusters at each level
00037 template<class T, class D>
00038 void mbl_cluster_tree<T,D>::set_max_r(const vcl_vector<double>& r)
00039 {
00040   empty();
00041   unsigned nL = r.size();
00042   cluster_.resize(nL);
00043   for (unsigned i=0;i<nL;++i) cluster_[i].set_max_r(r[i]);
00044   parent_.resize( nL);
00045   for (unsigned i=0;i<nL;++i) parent_[i].resize(0);
00046 
00047   cluster_[0].set_data(data_);
00048   for (unsigned i=1;i<nL;++i)
00049     cluster_[i].set_data(cluster_[i-1].p());
00050 }
00051 
00052 
00053 //: Define external data array (pointer retained)
00054 //  Empty existing clusters, then process every element of data
00055 //  to create clusters, by calling add_object()
00056 template<class T, class D>
00057 void mbl_cluster_tree<T,D>::set_data(const vcl_vector<T>& data)
00058 {
00059   empty();
00060   data_ = data;
00061   parent_[0].resize(data.size());
00062   unsigned n = data.size();
00063   for (unsigned i=0;i<n;++i)  add_object(i);
00064 }
00065 
00066 //: Add an extra element to data()
00067 template<class T, class D>
00068 void mbl_cluster_tree<T,D>::push_back(const T& t)
00069 {
00070   data_.push_back(t);
00071   parent_[0].push_back(0);  // Create space for parent
00072   add_object(data().size()-1);
00073 }
00074 
00075 //: Return index of nearest object in data() to t
00076 //  Nearest object in data() to t is given by data()[nearest(t,d)];
00077 //  The distance to the point is d
00078 template<class T, class D>
00079 unsigned mbl_cluster_tree<T,D>::nearest(const T& t, double& d) const
00080 {
00081   assert(data().size()>0);
00082   assert(cluster_.size()>0);
00083 
00084   if (cluster_.size()==1) return cluster_[0].nearest(t,d);
00085 
00086   // Perform hierarchical search
00087   // Find possible clusters at top level
00088   unsigned L=cluster_.size()-1;
00089   vcl_vector<unsigned> near_c0, near_c1;
00090   double max_d;
00091   cluster_[L].nearest_clusters(t,max_d,near_c1);
00092 
00093   cluster_[L].in_clusters(near_c1,near_c0);  // Find objects in next L
00094   L--;
00095   while (L!=0)
00096   {
00097     cluster_[L].nearest_clusters(t,max_d,near_c0,near_c1);
00098     cluster_[L].in_clusters(near_c1,near_c0);   // Find objects in next L
00099     --L;
00100   }
00101   return cluster_[0].nearest(t,d,near_c0);
00102 }
00103 
00104 //: Append new object with index i and assign to clusters
00105 //  Assumes that new object data()[i] is available.
00106 //  Deduce which cluster it belongs to and add it.
00107 //  Create new cluster if further than max_r() from any.
00108 template<class T, class D>
00109 void mbl_cluster_tree<T,D>::add_object(unsigned new_i)
00110 {
00111   assert(new_i<data().size());
00112 
00113   if (cluster_.size()==1) { cluster_[0].add_object(new_i); return; }
00114 
00115   const T& t = data()[new_i];
00116 
00117   unsigned Lhi = cluster_.size()-1;
00118 
00119   // Find any clusters at top level which could hold t
00120   vcl_vector<unsigned> c_list0;
00121   vcl_vector<unsigned>  nearest_c(Lhi+1);
00122   vcl_vector<double>  nearest_d(Lhi+1);
00123 
00124   if (cluster_[Lhi].clusters_within_max_r(t,c_list0,nearest_c[Lhi],
00125                                                    nearest_d[Lhi])==0)
00126   {
00127     // No clusters at the top level can include t, so create new
00128     // cluster at each level.
00129     unsigned cL0=cluster_[0].create_cluster(new_i);
00130     parent_[0][new_i]=cL0;
00131     for (unsigned L=1;L<=Lhi;++L)
00132     {
00133       unsigned cL1=cluster_[L].create_cluster(cL0);
00134       parent_[L].push_back(cL1);
00135       cL0=cL1;  // Record position
00136     }
00137     return;
00138   }
00139 
00140   for (unsigned L0=Lhi;L0!=0;--L0)
00141   {
00142     unsigned L = L0 - 1;
00143     vcl_vector<unsigned> c_list1=c_list0;
00144 
00145      // Generate list of elements in each cluster to process at next level
00146     cluster_[L0].in_clusters(c_list1,c_list0);
00147 
00148     // Find clusters at level L which could contain t
00149     unsigned nc=cluster_[L].clusters_within_max_r(t,c_list0,c_list1,
00150                                        nearest_c[L],nearest_d[L]);
00151 
00152     if (nc==0)
00153     {
00154       // No clusters at this level can include t,
00155       // so create new cluster at levels 0..L
00156       unsigned cL=cluster_[0].create_cluster(new_i);
00157       parent_[0][new_i]=cL;
00158       if (parent_.size()>1)
00159         parent_[1].push_back(0);  // Make space for new one
00160 
00161       for (unsigned L1=1;L1<=L;++L1)
00162       {
00163         unsigned cL1=cluster_[L1].create_cluster(cL);
00164         if (L1<Lhi)
00165           parent_[L1+1].push_back(0);  // Create space for record of new
00166         parent_[L1][cL]=cL1;
00167         cL=cL1;  // Record index
00168       }
00169 
00170       // cL is cluster containing new object at level L
00171       // Assign to nearest cluster in level above (c)
00172       unsigned c=nearest_c[L0];
00173       cluster_[L0].assign_to_cluster(cL,c,nearest_d[L0]);
00174       parent_[L0][cL]=c;  // Record parent
00175 
00176         // Track back through ancestors of c, updating radii
00177       for (unsigned L1=L0+1;L1<=Lhi;++L1)
00178       {
00179         c = parent_[L1][c];
00180         double d=D::d(t,cluster_[L1].p()[c]);
00181         if (d>cluster_[L1].r()[c])
00182           cluster_[L1].set_r(c,d);
00183       }
00184       return;
00185     }
00186     vcl_swap(c_list0,c_list1);  // Set c_list0 to current valid list
00187   }
00188 
00189   // If reached here, then t is in range of a cluster at every level
00190   unsigned c = nearest_c[0];
00191   cluster_[0].assign_to_cluster(new_i,c,nearest_d[0]);
00192   parent_[0][new_i]=c;
00193   for (unsigned L=1;L<=Lhi;++L)
00194   {
00195     // Track back through parents of c, updating radii
00196     c = parent_[L][c];
00197     double d=D::d(t,cluster_[L].p()[c]);
00198     if (d>cluster_[L].r()[c])
00199       cluster_[L].set_r(c,d);
00200   }
00201 }
00202 
00203 //: Print ancestry of every element
00204 template<class T, class D>
00205 void mbl_cluster_tree<T,D>::print_tree(vcl_ostream& os) const
00206 {
00207   for (unsigned i=0;i<data().size();++i)
00208   {
00209     os << i;
00210     unsigned p=i;
00211     for (unsigned L=0;L<cluster_.size();++L)
00212     {
00213       p=parent_[L][p];
00214       os << " - "<<p;
00215     }
00216     os<<'\n';
00217   }
00218 }
00219 
00220 //: Print summary information
00221 template<class T, class D>
00222 void mbl_cluster_tree<T,D>::print_summary(vcl_ostream& os) const
00223 {
00224   for (unsigned i=0;i<cluster_.size();++i)
00225   {
00226     os << "Level "<<i<<") max_r: "<<cluster_[i].max_r()
00227        << " n_clusters: "<<cluster_[i].p().size()<<'\n';
00228   }
00229 }
00230 
00231 template<class T, class D>
00232 short mbl_cluster_tree<T,D>::version_no() const
00233 {
00234     return 1;
00235 }
00236 
00237 template<class T, class D>
00238 void mbl_cluster_tree<T,D>::b_write(vsl_b_ostream& bfs) const
00239 {
00240   vsl_b_write(bfs,version_no());
00241   vsl_b_write(bfs,data_);
00242   vsl_b_write(bfs,parent_);
00243   // Write out clusters explicitly to avoid creating another template
00244   vsl_b_write(bfs,cluster_.size());
00245   for (unsigned L=0;L<cluster_.size();++L)
00246     vsl_b_write(bfs,cluster_[L]);
00247 }
00248 
00249 template<class T, class D>
00250 void mbl_cluster_tree<T,D>::b_read(vsl_b_istream& bfs)
00251 {
00252   short version;
00253   vsl_b_read(bfs,version);
00254   unsigned nc=0;
00255   switch (version)
00256   {
00257     case 1:
00258       vsl_b_read(bfs,data_);
00259       vsl_b_read(bfs,parent_);
00260       vsl_b_read(bfs,nc);
00261       cluster_.resize(nc);
00262       for (unsigned L=0;L<nc;++L)
00263         vsl_b_read(bfs,cluster_[L]);
00264    break;
00265 
00266   default:
00267     vcl_cerr << "mbl_cluster_tree<T,D>::b_read() "
00268       "Unexpected version number " << version << vcl_endl;
00269     bfs.is().clear(vcl_ios::badbit); // Set an unrecoverable IO error on stream
00270     return;
00271   }
00272 
00273   // Connect each cluster to data below
00274   cluster_[0].set_data(data_);
00275   for (unsigned L=1;L<cluster_.size();++L)
00276     cluster_[L].set_data(cluster_[L-1].p());
00277 }
00278 
00279 //: Binary file stream output operator for class reference
00280 template<class T, class D>
00281 void vsl_b_write(vsl_b_ostream& bfs, const mbl_cluster_tree<T,D>& c)
00282 {
00283   c.b_write(bfs);
00284 }
00285 
00286 //: Binary file stream input operator for class reference
00287 template<class T, class D>
00288 void vsl_b_read(vsl_b_istream& bfs, mbl_cluster_tree<T,D>& c)
00289 {
00290   c.b_read(bfs);
00291 }
00292 
00293 //: Stream output operator for class reference
00294 template<class T, class D>
00295 vcl_ostream& operator<<(vcl_ostream& os,const mbl_cluster_tree<T,D>& c)
00296 {
00297   c.print_summary(os);
00298   return os;
00299 }
00300 
00301 
00302 #define MBL_CLUSTER_TREE_INSTANTIATE(T,D) \
00303 template class mbl_cluster_tree<T,D >; \
00304 template void vsl_b_write(vsl_b_ostream& bfs, const mbl_cluster_tree<T,D >& c); \
00305 template void vsl_b_read(vsl_b_istream& bfs, mbl_cluster_tree<T,D >& c); \
00306 template vcl_ostream& operator<<(vcl_ostream&os,const mbl_cluster_tree<T,D >& c)
00307 
00308 #endif // mbl_cluster_txx_