vcs.maemo.org Git - opencv/blob - src/cv/cvlsh.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2009, Xavier Delacour, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41
  42 // 2009-01-12, Xavier Delacour <xavier.delacour@gmail.com>
  43
  44
  45 // * hash perf could be improved
  46 // * in particular, implement integer only (converted fixed from float input)
  47
  48 // * number of hash functions could be reduced (andoni paper)
  49
  50 // * redundant distance computations could be suppressed
  51
  52 // * rework CvLSHOperations interface-- move some of the loops into it to
  53 // * allow efficient async storage
  54
  55
  56 // Datar, M., Immorlica, N., Indyk, P., and Mirrokni, V. S. 2004. Locality-sensitive hashing
  57 // scheme based on p-stable distributions. In Proceedings of the Twentieth Annual Symposium on
  58 // Computational Geometry (Brooklyn, New York, USA, June 08 - 11, 2004). SCG '04. ACM, New York,
  59 // NY, 253-262. DOI= http://doi.acm.org/10.1145/997817.997857
  60
  61 #include "_cv.h"
  62 #include "cv.hpp"
  63 #include "cxmisc.h"
  64 #include <math.h>
  65 #include <vector>
  66 #include <algorithm>
  67 #include <limits>
  68
  69 template <class T>
  70 class memory_hash_ops : public CvLSHOperations {
  71   int d;
  72   std::vector<T> data;
  73   std::vector<int> free_data;
  74   struct node {
  75     int i, h2, next;
  76   };
  77   std::vector<node> nodes;
  78   std::vector<int> free_nodes;
  79   std::vector<int> bins;
  80
  81 public:
  82   memory_hash_ops(int _d, int n) : d(_d) {
  83     bins.resize(n, -1);
  84   }
  85
  86   virtual int vector_add(const void* _p) {
  87     const T* p = (const T*)_p;
  88     int i;
  89     if (free_data.empty()) {
  90       i = (int)data.size();
  91       data.insert(data.end(), d, 0);
  92     } else {
  93       i = free_data.end()[-1];
  94       free_data.pop_back();
  95     }
  96     std::copy(p, p + d, data.begin() + i);
  97     return i / d;
  98   }
  99   virtual void vector_remove(int i) {
 100     free_data.push_back(i * d);
 101   }
 102   virtual const void* vector_lookup(int i) {
 103     return &data[i * d];
 104   }
 105   virtual void vector_reserve(int n) {
 106     data.reserve(n * d);
 107   }
 108   virtual unsigned int vector_count() {
 109     return (unsigned)(data.size() / d - free_data.size());
 110   }
 111
 112   virtual void hash_insert(lsh_hash h, int /*l*/, int i) {
 113     int ii;
 114     if (free_nodes.empty()) {
 115       ii = (int)nodes.size();
 116       nodes.push_back(node());
 117     } else {
 118       ii = free_nodes.end()[-1];
 119       free_nodes.pop_back();
 120     }
 121     node& n = nodes[ii];
 122     int h1 = h.h1 % bins.size();
 123     n.i = i;
 124     n.h2 = h.h2;
 125     n.next = bins[h1];
 126     bins[h1] = ii;
 127   }
 128   virtual void hash_remove(lsh_hash h, int /*l*/, int i) {
 129     int h1 = h.h1 % bins.size();
 130     for (int ii = bins[h1], iin, iip = -1; ii != -1; iip = ii, ii = iin) {
 131       iin = nodes[ii].next;
 132       if (nodes[ii].h2 == h.h2 && nodes[ii].i == i) {
 133         free_nodes.push_back(ii);
 134         if (iip == -1)
 135           bins[h1] = iin;
 136         else
 137           nodes[iip].next = iin;
 138       }
 139     }
 140   }
 141   virtual int hash_lookup(lsh_hash h, int /*l*/, int* ret_i, int ret_i_max) {
 142     int h1 = h.h1 % bins.size();
 143     int k = 0;
 144     for (int ii = bins[h1]; ii != -1 && k < ret_i_max; ii = nodes[ii].next)
 145       if (nodes[ii].h2 == h.h2)
 146         ret_i[k++] = nodes[ii].i;
 147     return k;
 148   }
 149 };
 150
 151 template <class T,int cvtype>
 152 class pstable_l2_func {
 153   CvMat *a, *b, *r1, *r2;
 154   int d, k;
 155   double r;
 156   pstable_l2_func(const pstable_l2_func& x);
 157   pstable_l2_func& operator= (const pstable_l2_func& rhs);
 158 public:
 159   typedef T scalar_type;
 160   typedef T accum_type;
 161   pstable_l2_func(int _d, int _k, double _r, CvRNG& rng)
 162     : d(_d), k(_k), r(_r) {
 163     assert(sizeof(T) == CV_ELEM_SIZE1(cvtype));
 164     a = cvCreateMat(k, d, cvtype);
 165     b = cvCreateMat(k, 1, cvtype);
 166     r1 = cvCreateMat(k, 1, CV_32SC1);
 167     r2 = cvCreateMat(k, 1, CV_32SC1);
 168     cvRandArr(&rng, a, CV_RAND_NORMAL, cvScalar(0), cvScalar(1));
 169     cvRandArr(&rng, b, CV_RAND_UNI, cvScalar(0), cvScalar(r));
 170     cvRandArr(&rng, r1, CV_RAND_UNI,
 171               cvScalar(std::numeric_limits<int>::min()),
 172               cvScalar(std::numeric_limits<int>::max()));
 173     cvRandArr(&rng, r2, CV_RAND_UNI,
 174               cvScalar(std::numeric_limits<int>::min()),
 175               cvScalar(std::numeric_limits<int>::max()));
 176   }
 177   ~pstable_l2_func() {
 178     cvReleaseMat(&a);
 179     cvReleaseMat(&b);
 180     cvReleaseMat(&r1);
 181     cvReleaseMat(&r2);
 182   }
 183
 184   // * factor all L functions into this (reduces number of matrices to 4 total;
 185   // * simpler syntax in lsh_table). give parameter l here that tells us which
 186   // * row to use etc.
 187
 188   lsh_hash operator() (const T* x) const {
 189     const T* aj = (const T*)a->data.ptr;
 190     const T* bj = (const T*)b->data.ptr;
 191
 192     lsh_hash h;
 193     h.h1 = h.h2 = 0;
 194     for (int j = 0; j < k; ++j) {
 195       accum_type s = 0;
 196       for (int jj = 0; jj < d; ++jj)
 197         s += aj[jj] * x[jj];
 198       s += *bj;
 199       s = accum_type(s/r);
 200       int si = int(s);
 201       h.h1 += r1->data.i[j] * si;
 202       h.h2 += r2->data.i[j] * si;
 203
 204       aj += d;
 205       bj++;
 206     }
 207     return h;
 208   }
 209   accum_type distance(const T* p, const T* q) const {
 210     accum_type s = 0;
 211     for (int j = 0; j < d; ++j) {
 212       accum_type d1 = p[j] - q[j];
 213       s += d1 * d1;
 214     }
 215     return s;
 216   }
 217 };
 218
 219 template <class H>
 220 class lsh_table {
 221 public:
 222   typedef typename H::scalar_type scalar_type;
 223   typedef typename H::accum_type accum_type;
 224 private:
 225   std::vector<H*> g;
 226   CvLSHOperations* ops;
 227   int d, L, k;
 228   double r;
 229
 230   static accum_type comp_dist(const std::pair<int,accum_type>& x,
 231                               const std::pair<int,accum_type>& y) {
 232     return x.second < y.second;
 233   }
 234
 235   lsh_table(const lsh_table& x);
 236   lsh_table& operator= (const lsh_table& rhs);
 237 public:
 238   lsh_table(CvLSHOperations* _ops, int _d, int _L, int _k, double _r, CvRNG& rng)
 239     : ops(_ops), d(_d), L(_L), k(_k), r(_r) {
 240     g.resize(L);
 241     for (int j = 0; j < L; ++j)
 242       g[j] = new H(d, k, r, rng);
 243   }
 244   ~lsh_table() {
 245     for (int j = 0; j < L; ++j)
 246       delete g[j];
 247     delete ops;
 248   }
 249
 250   int dims() const {
 251     return d;
 252   }
 253   unsigned int size() const {
 254     return ops->vector_count();
 255   }
 256
 257   void add(const scalar_type* data, int n, int* ret_indices = 0) {
 258     for (int j=0;j<n;++j) {
 259       const scalar_type* x = data+j*d;
 260       int i = ops->vector_add(x);
 261       if (ret_indices)
 262         ret_indices[j] = i;
 263
 264       for (int l = 0; l < L; ++l) {
 265         lsh_hash h = (*g[l])(x);
 266         ops->hash_insert(h, l, i);
 267       }
 268     }
 269   }
 270   void remove(const int* indices, int n) {
 271     for (int j = 0; j < n; ++j) {
 272       int i = indices[n];
 273       const scalar_type* x = (const scalar_type*)ops->vector_lookup(i);
 274
 275       for (int l = 0; l < L; ++l) {
 276         lsh_hash h = (*g[l])(x);
 277         ops->hash_remove(h, l, i);
 278       }
 279       ops->vector_remove(i);
 280     }
 281   }
 282   void query(const scalar_type* q, int k0, int emax, double* dist, int* results) {
 283     int* tmp = (int*)cvStackAlloc(sizeof(int) * emax);
 284     typedef std::pair<int, accum_type> dr_type; // * swap int and accum_type here, for naming consistency
 285     dr_type* dr = (dr_type*)cvStackAlloc(sizeof(dr_type) * k0);
 286     int k1 = 0;
 287
 288     // * handle k0 >= emax, in which case don't track max distance
 289
 290     for (int l = 0; l < L && emax > 0; ++l) {
 291       lsh_hash h = (*g[l])(q);
 292       int m = ops->hash_lookup(h, l, tmp, emax);
 293       for (int j = 0; j < m && emax > 0; ++j, --emax) {
 294         int i = tmp[j];
 295         const scalar_type* p = (const scalar_type*)ops->vector_lookup(i);
 296         accum_type pd = (*g[l]).distance(p, q);
 297         if (k1 < k0) {
 298           dr[k1++] = std::make_pair(i, pd);
 299           std::push_heap(dr, dr + k1, comp_dist);
 300         } else if (pd < dr[0].second) {
 301           std::pop_heap(dr, dr + k0, comp_dist);
 302           dr[k0 - 1] = std::make_pair(i, pd);
 303           std::push_heap(dr, dr + k0, comp_dist);
 304         }
 305       }
 306     }
 307
 308     for (int j = 0; j < k1; ++j)
 309       dist[j] = dr[j].second, results[j] = dr[j].first;
 310     std::fill(dist + k1, dist + k0, 0);
 311     std::fill(results + k1, results + k0, -1);
 312   }
 313   void query(const scalar_type* data, int n, int k0, int emax, double* dist, int* results) {
 314     for (int j = 0; j < n; ++j) {
 315       query(data, k0, emax, dist, results);
 316       data += d; // * this may not agree with step for some scalar_types
 317       dist += k0;
 318       results += k0;
 319     }
 320   }
 321 };
 322
 323 typedef lsh_table<pstable_l2_func<float, CV_32FC1> > lsh_pstable_l2_32f;
 324 typedef lsh_table<pstable_l2_func<double, CV_64FC1> > lsh_pstable_l2_64f;
 325
 326 struct CvLSH {
 327   int type;
 328   union {
 329     lsh_pstable_l2_32f* lsh_32f;
 330     lsh_pstable_l2_64f* lsh_64f;
 331   } u;
 332 };
 333
 334 CvLSH* cvCreateLSH(CvLSHOperations* ops, int d, int L, int k, int type, double r, int64 seed) {
 335   CvLSH* lsh = 0;
 336   CvRNG rng = cvRNG(seed);
 337
 338   __BEGIN__;
 339   CV_FUNCNAME("cvCreateLSH");
 340
 341   if (type != CV_32FC1 && type != CV_64FC1)
 342     CV_ERROR(CV_StsUnsupportedFormat, "vectors must be either CV_32FC1 or CV_64FC1");
 343   lsh = new CvLSH;
 344   lsh->type = type;
 345   switch (type) {
 346   case CV_32FC1: lsh->u.lsh_32f = new lsh_pstable_l2_32f(ops, d, L, k, r, rng); break;
 347   case CV_64FC1: lsh->u.lsh_64f = new lsh_pstable_l2_64f(ops, d, L, k, r, rng); break;
 348   }
 349
 350   __END__;
 351   return lsh;
 352 }
 353
 354 CvLSH* cvCreateMemoryLSH(int d, int n, int L, int k, int type, double r, int64 seed) {
 355   CvLSHOperations* ops = 0;
 356
 357   switch (type) {
 358   case CV_32FC1: ops = new memory_hash_ops<float>(d,n); break;
 359   case CV_64FC1: ops = new memory_hash_ops<double>(d,n); break;
 360   }
 361   return cvCreateLSH(ops, d, L, k, type, r, seed);
 362 }
 363
 364 void cvReleaseLSH(CvLSH** lsh) {
 365   switch ((*lsh)->type) {
 366   case CV_32FC1: delete (*lsh)->u.lsh_32f; break;
 367   case CV_64FC1: delete (*lsh)->u.lsh_64f; break;
 368   default: assert(0);
 369   }
 370   delete *lsh;
 371   *lsh = 0;
 372 }
 373
 374 unsigned int LSHSize(CvLSH* lsh) {
 375   switch (lsh->type) {
 376   case CV_32FC1: return lsh->u.lsh_32f->size(); break;
 377   case CV_64FC1: return lsh->u.lsh_64f->size(); break;
 378   default: assert(0);
 379   }
 380   return 0;
 381 }
 382
 383
 384 void cvLSHAdd(CvLSH* lsh, const CvMat* data, CvMat* indices) {
 385   int dims, n;
 386   int* ret_indices = 0;
 387
 388   __BEGIN__;
 389   CV_FUNCNAME("cvLSHAdd");
 390
 391   switch (lsh->type) {
 392   case CV_32FC1: dims = lsh->u.lsh_32f->dims(); break;
 393   case CV_64FC1: dims = lsh->u.lsh_64f->dims(); break;
 394   default: assert(0); return;
 395   }
 396
 397   n = data->rows;
 398
 399   if (dims != data->cols)
 400     CV_ERROR(CV_StsBadSize, "data must be n x d, where d is what was used to construct LSH");
 401
 402   if (CV_MAT_TYPE(data->type) != lsh->type)
 403     CV_ERROR(CV_StsUnsupportedFormat, "type of data and constructed LSH must agree");
 404   if (indices) {
 405     if (CV_MAT_TYPE(indices->type) != CV_32SC1)
 406       CV_ERROR(CV_StsUnsupportedFormat, "indices must be CV_32SC1");
 407     if (indices->rows * indices->cols != n)
 408       CV_ERROR(CV_StsBadSize, "indices must be n x 1 or 1 x n for n x d data");
 409     ret_indices = indices->data.i;
 410   }
 411
 412   switch (lsh->type) {
 413   case CV_32FC1: lsh->u.lsh_32f->add(data->data.fl, n, ret_indices); break;
 414   case CV_64FC1: lsh->u.lsh_64f->add(data->data.db, n, ret_indices); break;
 415   default: assert(0); return;
 416   }
 417   __END__;
 418 }
 419
 420 void cvLSHRemove(CvLSH* lsh, const CvMat* indices) {
 421   int n;
 422
 423   __BEGIN__;
 424   CV_FUNCNAME("cvLSHRemove");
 425
 426   if (CV_MAT_TYPE(indices->type) != CV_32SC1)
 427     CV_ERROR(CV_StsUnsupportedFormat, "indices must be CV_32SC1");
 428   n = indices->rows * indices->cols;
 429   switch (lsh->type) {
 430   case CV_32FC1: lsh->u.lsh_32f->remove(indices->data.i, n); break;
 431   case CV_64FC1: lsh->u.lsh_64f->remove(indices->data.i, n); break;
 432   default: assert(0); return;
 433   }
 434   __END__;
 435 }
 436
 437 void cvLSHQuery(CvLSH* lsh, const CvMat* data, CvMat* indices, CvMat* dist, int k, int emax) {
 438   int dims;
 439
 440   __BEGIN__;
 441   CV_FUNCNAME("cvLSHQuery");
 442
 443   switch (lsh->type) {
 444   case CV_32FC1: dims = lsh->u.lsh_32f->dims(); break;
 445   case CV_64FC1: dims = lsh->u.lsh_64f->dims(); break;
 446   default: assert(0); return;
 447   }
 448
 449   if (k<1)
 450     CV_ERROR(CV_StsOutOfRange, "k must be positive");
 451   if (CV_MAT_TYPE(data->type) != lsh->type)
 452     CV_ERROR(CV_StsUnsupportedFormat, "type of data and constructed LSH must agree");
 453   if (dims != data->cols)
 454     CV_ERROR(CV_StsBadSize, "data must be n x d, where d is what was used to construct LSH");
 455   if (dist->rows != data->rows || dist->cols != k)
 456     CV_ERROR(CV_StsBadSize, "dist must be n x k for n x d data");
 457   if (dist->rows != indices->rows || dist->cols != indices->cols)
 458     CV_ERROR(CV_StsBadSize, "dist and indices must be same size");
 459   if (CV_MAT_TYPE(dist->type) != CV_64FC1)
 460     CV_ERROR(CV_StsUnsupportedFormat, "dist must be CV_64FC1");
 461   if (CV_MAT_TYPE(indices->type) != CV_32SC1)
 462     CV_ERROR(CV_StsUnsupportedFormat, "indices must be CV_32SC1");
 463
 464   switch (lsh->type) {
 465   case CV_32FC1: lsh->u.lsh_32f->query(data->data.fl, data->rows,
 466                                        k, emax, dist->data.db, indices->data.i); break;
 467   case CV_64FC1: lsh->u.lsh_64f->query(data->data.db, data->rows,
 468                                        k, emax, dist->data.db, indices->data.i); break;
 469   default: assert(0); return;
 470   }
 471   __END__;
 472 }