contrib/mul/mbl/mbl_clusters.txx
Go to the documentation of this file.
00001 // This is mul/mbl/mbl_clusters.txx
00002 #ifndef mbl_clusters_txx_
00003 #define mbl_clusters_txx_
00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00005 #pragma interface
00006 #endif
00007 //:
00008 // \file
00009 // \brief  Class to record clusters of data, for faster neighbour finding
00010 // \author Tim Cootes
00011 
00012 #include "mbl_clusters.h"
00013 #include <vcl_iostream.h>
00014 #include <vcl_cassert.h>
00015 #include <vcl_algorithm.h>
00016 #include <vsl/vsl_binary_io.h>
00017 #include <vsl/vsl_vector_io.h>
00018 
00019 //: Default constructor
00020 template<class T, class D>
00021 mbl_clusters<T,D>::mbl_clusters()
00022   : data_(0)
00023 {
00024 }
00025 
00026 //: Define maximum radius for each cluster
00027 template<class T, class D>
00028 void mbl_clusters<T,D>::set_max_r(double r)
00029 {
00030   max_r_ = r;
00031 }
00032 
00033 //: Empty clusters
00034 template<class T, class D>
00035 void mbl_clusters<T,D>::empty()
00036 {
00037   p_.resize(0);
00038   index_.resize(0);
00039   r_.resize(0);
00040 }
00041 
00042 //: Define external data array (pointer retained)
00043 //  Empty existing clusters, then process every element of data
00044 //  to create clusters, by calling add_object()
00045 template<class T, class D>
00046 void mbl_clusters<T,D>::set_data(const vcl_vector<T>& data)
00047 {
00048   empty();
00049   data_ = &data;
00050   unsigned n = data.size();
00051   for (unsigned i=0;i<n;++i)  add_object(i);
00052 }
00053 
00054 //: Define external data array (pointer retained)
00055 //  Use carefully! This sets the internal pointer to
00056 //  point to data.  Really only to be used after loading
00057 //  internals using b_read(bfs).
00058 template<class T, class D>
00059 void mbl_clusters<T,D>::set_data_ptr(const vcl_vector<T>& data)
00060 {
00061   data_=&data;
00062 }
00063 
00064 //: Return index of nearest object in data() to t
00065 //  Nearest object in data() to t is given by data()[nearest(t,d)];
00066 //  The distance to the point is d
00067 template<class T, class D>
00068 unsigned mbl_clusters<T,D>::nearest(const T& t, double& d) const
00069 {
00070   assert(data_!=0);
00071 
00072   // Initialise with first in data
00073   unsigned best_i = 0;
00074   d = D::d(data()[0],t);
00075 
00076   const T* data_ptr = &data()[0];
00077 
00078   // Try each cluster in turn
00079   for (unsigned j=0;j<p_.size();++j)
00080   {
00081     double dj = D::d(t,p_[j]);
00082     if (dj-r_[j]<d)
00083     {
00084       // There may be a point in the cluster closer than the current best
00085       const vcl_vector<unsigned>& ind = index_[j];
00086       for (unsigned i=0;i<ind.size();++i)
00087       {
00088         double di=D::d(data_ptr[ind[i]],t);
00089         if (di<d) { d=di; best_i=ind[i]; }
00090       }
00091     }
00092   }
00093   return best_i;
00094 }
00095 
00096 //: Return index of nearest object in data() to t
00097 //  Consider only objects in clusters given in c_list
00098 //  Nearest object in data() to t is given by data()[nearest(t,d)];
00099 //  The distance to the point is d
00100 template<class T, class D>
00101 unsigned mbl_clusters<T,D>::nearest(const T& t, double& d,
00102                    const vcl_vector<unsigned>& c_list) const
00103 {
00104   assert(data_!=0);
00105 
00106   // Initialise with first in set for c_list[0]
00107   unsigned best_i = 0;
00108   d = D::d(data()[index_[c_list[0]][0]],t);
00109 
00110   const T* data_ptr = &data()[0];
00111 
00112   // Try each cluster in turn
00113   for (unsigned k=0;k<c_list.size();++k)
00114   {
00115     unsigned j=c_list[k];
00116     double dj = D::d(t,p_[j]);
00117     if (dj-r_[j]<d)
00118     {
00119       // There may be a point in the cluster closer than the current best
00120       const vcl_vector<unsigned>& ind = index_[j];
00121       for (unsigned i=0;i<ind.size();++i)
00122       {
00123         double di=D::d(data_ptr[ind[i]],t);
00124         if (di<d) { d=di; best_i=ind[i]; }
00125       }
00126     }
00127   }
00128   return best_i;
00129 }
00130 
00131 
00132 //: Return index of nearest cluster in data() to t
00133 //  Finds nearest cluster key point to t
00134 //  The distance to the point is d
00135 template<class T, class D>
00136 unsigned mbl_clusters<T,D>::nearest_cluster(const T& t, double& d) const
00137 {
00138   assert(p_.size()>0);
00139 
00140   d = D::d(p_[0],t);
00141   unsigned best_j = 0;
00142 
00143   // Try each cluster in turn
00144   for (unsigned j=1;j<p_.size();++j)
00145   {
00146     double dj = D::d(t,p_[j]);
00147     if (dj<d) { d=dj; best_j=j; }
00148   }
00149 
00150   return best_j;
00151 }
00152 
00153 //: Return indices of clusters which may contain nearest point to t
00154 //  Searches through all the clusters
00155 template<class T, class D>
00156 void mbl_clusters<T,D>::nearest_clusters(const T& t, double& max_d,
00157                            vcl_vector<unsigned>& near_c) const
00158 {
00159   assert(p_.size()>0);
00160 
00161   vcl_vector<unsigned> c1;
00162   vcl_vector<double> d1;
00163 
00164   double d = D::d(p_[0],t);
00165   c1.push_back(0);
00166   d1.push_back(d);
00167   max_d = d+r_[0];
00168 
00169   // Try each cluster in turn, recording any that might include closest
00170   for (unsigned j=1;j<p_.size();++j)
00171   {
00172     double dj = D::d(t,p_[j]);
00173     if (dj-r_[j]<=max_d)
00174     {
00175       c1.push_back(j);
00176       d1.push_back(dj);
00177       max_d=vcl_min(max_d,dj+r_[j]);
00178     }
00179   }
00180 
00181   // Pass through the data again to prune out far clusters
00182   near_c.resize(0);
00183   for (unsigned i=0;i<c1.size();++i)
00184   {
00185     if (d1[i]-r_[c1[i]]<=max_d) near_c.push_back(c1[i]);
00186   }
00187 }
00188 
00189 //: Return indices of clusters which may contain nearest point to t
00190 //  Searches through clusters listed in c_list.
00191 //  On input, max_d gives initial limit on distance.
00192 //  On exit, max_d gives the revised limit on the distance
00193 template<class T, class D>
00194 void mbl_clusters<T,D>::nearest_clusters(const T& t, double& max_d,
00195                            const vcl_vector<unsigned>& c_list,
00196                            vcl_vector<unsigned>& near_c) const
00197 {
00198   assert(p_.size()>0);
00199 
00200   // Storage for first pass
00201   vcl_vector<unsigned> c1;
00202   vcl_vector<double> d1;
00203 
00204   // Try each cluster in turn, recording any that might include closest
00205   for (unsigned k=0;k<c_list.size();++k)
00206   {
00207     unsigned j=c_list[k];
00208     double dj = D::d(t,p_[j]);
00209     if (dj-r_[j]<=max_d)
00210     {
00211       c1.push_back(j);
00212       d1.push_back(dj);
00213       max_d=vcl_min(max_d,dj+r_[j]);
00214     }
00215   }
00216 
00217   // Pass through the data again to prune out far clusters
00218   near_c.resize(0);
00219   for (unsigned i=0;i<c1.size();++i)
00220   {
00221     if (d1[i]-r_[c1[i]]<=max_d) near_c.push_back(c1[i]);
00222   }
00223 }
00224 
00225 //: Create a new cluster around point index i
00226 //  Return index of cluster
00227 template<class T, class D>
00228 unsigned mbl_clusters<T,D>::create_cluster(unsigned new_i, double r)
00229 {
00230   // Create a new cluster using this as a key point
00231   p_.push_back(data()[new_i]);
00232   r_.push_back(r);
00233   vcl_vector<unsigned> ind(1);
00234   ind[0]=new_i;
00235   index_.push_back(ind);
00236   return p_.size()-1;
00237 }
00238 
00239 //: Append new object with index i and assign to a cluster
00240 //  Assumes that new object data()[i] is available.
00241 //  Deduce which cluster it belongs to and add it.
00242 //  Create new cluster if further than max_r() from any.
00243 //  Return index of cluster it is assigned to
00244 template<class T, class D>
00245 unsigned mbl_clusters<T,D>::add_object(unsigned new_i, double r)
00246 {
00247   assert(new_i<data_->size());
00248 
00249   // If initially empty, create one cluster
00250   if (p_.size()==0) return create_cluster(new_i,r);
00251 
00252   double d;
00253   unsigned j=nearest_cluster(data()[new_i],d);
00254   d+=r;  // Allow for new_i being key point of another cluster, radius r
00255   if (d<max_r_)
00256   {
00257     // Add it to cluster j
00258     index_[j].push_back(new_i);
00259     if (d>r_[j]) r_[j]=d;  // Update max radius of cluster
00260   }
00261   else
00262     j=create_cluster(new_i,r);
00263 
00264   return j;
00265 }
00266 
00267 //: Assign object data()[i] to cluster ci, knowing distance r.
00268 //  r is the distance D::d(data()[i],p()[ci])
00269 template<class T, class D>
00270 void mbl_clusters<T,D>::assign_to_cluster(unsigned i, unsigned ci,
00271                                           double r)
00272 {
00273   index_[ci].push_back(i);
00274   if (r>r_[ci]) r_[ci]=r;
00275 }
00276 
00277 
00278 //: Append new object with index i (data()[i]), creating a new cluster
00279 //  Return index of resulting cluster, which is initialised with
00280 //  given radius.
00281 unsigned add_cluster(unsigned i, double r=0.0);
00282 
00283 
00284 //: Finds list of clusters whose keypoint is within d of t
00285 //  Returns number of such clusters. If >0, then nearest_c
00286 //  gives index of cluster with centre nearest to t
00287 template<class T, class D>
00288 unsigned mbl_clusters<T,D>::clusters_within_d(const T& t, double d,
00289                                               vcl_vector<unsigned>& c_list,
00290                                               unsigned& nearest_c,
00291                                               double& min_d)
00292 {
00293   c_list.resize(0);
00294   nearest_c=0;
00295   min_d = d+1;
00296   for (unsigned i=0;i<p_.size();++i)
00297   {
00298     double di=D::d(t,p_[i]);
00299     if (di<=d)
00300     {
00301       c_list.push_back(i);
00302       if (di<min_d) { nearest_c=i; min_d=di; }
00303     }
00304   }
00305   return c_list.size();
00306 }
00307 
00308 //: Finds list of clusters whose keypoint is within d of t
00309 //  Returns number of such clusters. If >0, then nearest_c
00310 //  gives index of cluster with centre nearest to t
00311 template<class T, class D>
00312 unsigned mbl_clusters<T,D>::clusters_within_d(const T& t, double d,
00313                                               const vcl_vector<unsigned>& in_list,
00314                                               vcl_vector<unsigned>& c_list,
00315                                               unsigned& nearest_c,
00316                                               double& min_d)
00317 {
00318   c_list.resize(0);
00319   nearest_c=0;
00320   min_d = d+1;
00321   for (unsigned j=0;j<in_list.size();++j)
00322   {
00323     unsigned i=in_list[j];
00324     double di=D::d(t,p_[i]);
00325     if (di<=d)
00326     {
00327       c_list.push_back(i);
00328       if (di<min_d) { nearest_c=i; min_d=di; }
00329     }
00330   }
00331   return c_list.size();
00332 }
00333 
00334 
00335 //: Finds list of clusters whose keypoint is within max_r() of t
00336 //  Returns number of such clusters. If >0, then nearest_c
00337 //  gives index of cluster with centre nearest to t
00338 template<class T, class D>
00339 unsigned mbl_clusters<T,D>::clusters_within_max_r(const T& t,
00340                                                   vcl_vector<unsigned>& c_list,
00341                                                   unsigned& nearest_c,
00342                                                   double& min_d)
00343 {
00344   return clusters_within_d(t,max_r(),c_list,nearest_c,min_d);
00345 }
00346 
00347 //: Finds list of clusters whose keypoint is within max_r() of t
00348 template<class T, class D>
00349 unsigned mbl_clusters<T,D>::clusters_within_max_r(const T& t,
00350                                                   const vcl_vector<unsigned>& in_list,
00351                                                   vcl_vector<unsigned>& c_list,
00352                                                   unsigned& nearest_c,
00353                                                   double& min_d)
00354 {
00355   return clusters_within_d(t,max_r(),in_list,c_list,nearest_c,min_d);
00356 }
00357 
00358 //: Create list of object indices in listed clusters
00359 //  Concatenates lists of indices for each cluster in c_list
00360 template<class T, class D>
00361 void mbl_clusters<T,D>::in_clusters(const vcl_vector<unsigned>& c_list,
00362                                     vcl_vector<unsigned>& o_list) const
00363 {
00364   o_list.resize(0);
00365   for (unsigned i=0;i<c_list.size();++i)
00366   {
00367     const vcl_vector<unsigned>& ind = index()[c_list[i]];
00368     for (unsigned j=0;j<ind.size();++j) o_list.push_back(ind[j]);
00369   }
00370 }
00371 
00372 
00373 //: Write out list of elements in each cluster
00374 template<class T, class D>
00375 void mbl_clusters<T,D>::print_cluster_sets(vcl_ostream& os) const
00376 {
00377   for (unsigned i=0;i<index_.size();++i)
00378   {
00379     os << i << ") ";
00380     for (unsigned j=0;j<index_[i].size();++j)
00381       os<<index_[i][j] << ' ';
00382     os<<'\n';
00383   }
00384 }
00385 
00386 //: Write out list of elements in each cluster
00387 template<class T, class D>
00388 void mbl_clusters<T,D>::print_summary(vcl_ostream& os) const
00389 {
00390   os << " max_r: " << max_r_ << " n_clusters: " << p_.size();
00391 }
00392 
00393 template<class T, class D>
00394 short mbl_clusters<T,D>::version_no() const
00395 {
00396     return 1;
00397 }
00398 
00399 template<class T, class D>
00400 void mbl_clusters<T,D>::b_write(vsl_b_ostream& bfs) const
00401 {
00402   vsl_b_write(bfs,version_no());
00403   vsl_b_write(bfs,p_);
00404   vsl_b_write(bfs,r_);
00405   vsl_b_write(bfs,max_r_);
00406   vsl_b_write(bfs,index_);
00407 }
00408 
00409 template<class T, class D>
00410 void mbl_clusters<T,D>::b_read(vsl_b_istream& bfs)
00411 {
00412   short version;
00413   vsl_b_read(bfs,version);
00414   switch (version)
00415   {
00416     case 1:
00417       vsl_b_read(bfs,p_);
00418       vsl_b_read(bfs,r_);
00419       vsl_b_read(bfs,max_r_);
00420       vsl_b_read(bfs,index_);
00421     break;
00422 
00423   default:
00424     vcl_cerr << "mbl_clusters<T,D>::b_read() "
00425       "Unexpected version number " << version << vcl_endl;
00426     bfs.is().clear(vcl_ios::badbit); // Set an unrecoverable IO error on stream
00427     return;
00428   }
00429 
00430   data_=0;
00431 }
00432 
00433 //: Binary file stream output operator for class reference
00434 template<class T, class D>
00435 void vsl_b_write(vsl_b_ostream& bfs, const mbl_clusters<T,D>& c)
00436 {
00437   c.b_write(bfs);
00438 }
00439 
00440 //: Binary file stream input operator for class reference
00441 template<class T, class D>
00442 void vsl_b_read(vsl_b_istream& bfs, mbl_clusters<T,D>& c)
00443 {
00444   c.b_read(bfs);
00445 }
00446 
00447 //: Stream output operator for class reference
00448 template<class T, class D>
00449 vcl_ostream& operator<<(vcl_ostream& os,const mbl_clusters<T,D>& c)
00450 {
00451   c.print_summary(os);
00452   return os;
00453 }
00454 
00455 
00456 #define MBL_CLUSTERS_INSTANTIATE(T,D) \
00457 template class mbl_clusters< T,D >; \
00458 template void vsl_b_write(vsl_b_ostream& bfs, const mbl_clusters<T,D >& c); \
00459 template void vsl_b_read(vsl_b_istream& bfs, mbl_clusters<T,D >& c); \
00460 template vcl_ostream& operator<<(vcl_ostream& os,const mbl_clusters<T,D >& c)
00461 
00462 #endif // mbl_clusters_txx_