modules/legacy/src/lsh.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 "precomp.hpp"
  62 #include <math.h>
  63 #include <vector>
  64 #include <algorithm>
  65 #include <limits>
  66
  67 template <class T>
  68 class memory_hash_ops : public CvLSHOperations {
  69   int d;
  70   std::vector<T> data;
  71   std::vector<int> free_data;
  72   struct node {
  73     int i, h2, next;
  74   };
  75   std::vector<node> nodes;
  76   std::vector<int> free_nodes;
  77   std::vector<int> bins;
  78
  79 public:
  80   memory_hash_ops(int _d, int n) : d(_d) {
  81     bins.resize(n, -1);
  82   }
  83
  84   virtual int vector_add(const void* _p) {
  85     const T* p = (const T*)_p;
  86     int i;
  87     if (free_data.empty()) {
  88       i = (int)data.size();
  89       data.insert(data.end(), d, 0);
  90     } else {
  91       i = free_data.end()[-1];
  92       free_data.pop_back();
  93     }
  94     std::copy(p, p + d, data.begin() + i);
  95     return i / d;
  96   }
  97   virtual void vector_remove(int i) {
  98     free_data.push_back(i * d);
  99   }
 100   virtual const void* vector_lookup(int i) {
 101     return &data[i * d];
 102   }
 103   virtual void vector_reserve(int n) {
 104     data.reserve(n * d);
 105   }
 106   virtual unsigned int vector_count() {
 107     return (unsigned)(data.size() / d - free_data.size());
 108   }
 109
 110   virtual void hash_insert(lsh_hash h, int /*l*/, int i) {
 111     int ii;
 112     if (free_nodes.empty()) {
 113       ii = (int)nodes.size();
 114       nodes.push_back(node());
 115     } else {
 116       ii = free_nodes.end()[-1];
 117       free_nodes.pop_back();
 118     }
 119     node& n = nodes[ii];
 120     int h1 = (int)(h.h1 % bins.size());
 121     n.i = i;
 122     n.h2 = h.h2;
 123     n.next = bins[h1];
 124     bins[h1] = ii;
 125   }
 126   virtual void hash_remove(lsh_hash h, int /*l*/, int i) {
 127     int h1 = (int)(h.h1 % bins.size());
 128     for (int ii = bins[h1], iin, iip = -1; ii != -1; iip = ii, ii = iin) {
 129       iin = nodes[ii].next;
 130       if (nodes[ii].h2 == h.h2 && nodes[ii].i == i) {
 131     free_nodes.push_back(ii);
 132     if (iip == -1)
 133       bins[h1] = iin;
 134     else
 135       nodes[iip].next = iin;
 136       }
 137     }
 138   }
 139   virtual int hash_lookup(lsh_hash h, int /*l*/, int* ret_i, int ret_i_max) {
 140     int h1 = (int)(h.h1 % bins.size());
 141     int k = 0;
 142     for (int ii = bins[h1]; ii != -1 && k < ret_i_max; ii = nodes[ii].next)
 143       if (nodes[ii].h2 == h.h2)
 144     ret_i[k++] = nodes[ii].i;
 145     return k;
 146   }
 147 };
 148
 149 template <class T,int cvtype>
 150 class pstable_l2_func {
 151   CvMat *a, *b, *r1, *r2;
 152   int d, k;
 153   double r;
 154   pstable_l2_func(const pstable_l2_func& x);
 155   pstable_l2_func& operator= (const pstable_l2_func& rhs);
 156 public:
 157   typedef T scalar_type;
 158   typedef T accum_type;
 159   pstable_l2_func(int _d, int _k, double _r, CvRNG& rng)
 160     : d(_d), k(_k), r(_r) {
 161     assert(sizeof(T) == CV_ELEM_SIZE1(cvtype));
 162     a = cvCreateMat(k, d, cvtype);
 163     b = cvCreateMat(k, 1, cvtype);
 164     r1 = cvCreateMat(k, 1, CV_32SC1);
 165     r2 = cvCreateMat(k, 1, CV_32SC1);
 166     cvRandArr(&rng, a, CV_RAND_NORMAL, cvScalar(0), cvScalar(1));
 167     cvRandArr(&rng, b, CV_RAND_UNI, cvScalar(0), cvScalar(r));
 168     cvRandArr(&rng, r1, CV_RAND_UNI,
 169           cvScalar(std::numeric_limits<int>::min()),
 170           cvScalar(std::numeric_limits<int>::max()));
 171     cvRandArr(&rng, r2, CV_RAND_UNI,
 172           cvScalar(std::numeric_limits<int>::min()),
 173           cvScalar(std::numeric_limits<int>::max()));
 174   }
 175   ~pstable_l2_func() {
 176     cvReleaseMat(&a);
 177     cvReleaseMat(&b);
 178     cvReleaseMat(&r1);
 179     cvReleaseMat(&r2);
 180   }
 181
 182   // * factor all L functions into this (reduces number of matrices to 4 total;
 183   // * simpler syntax in lsh_table). give parameter l here that tells us which
 184   // * row to use etc.
 185
 186   lsh_hash operator() (const T* x) const {
 187     const T* aj = (const T*)a->data.ptr;
 188     const T* bj = (const T*)b->data.ptr;
 189
 190     lsh_hash h;
 191     h.h1 = h.h2 = 0;
 192     for (int j = 0; j < k; ++j) {
 193       accum_type s = 0;
 194       for (int jj = 0; jj < d; ++jj)
 195     s += aj[jj] * x[jj];
 196       s += *bj;
 197       s = accum_type(s/r);
 198       int si = int(s);
 199       h.h1 += r1->data.i[j] * si;
 200       h.h2 += r2->data.i[j] * si;
 201
 202       aj += d;
 203       bj++;
 204     }
 205     return h;
 206   }
 207   accum_type distance(const T* p, const T* q) const {
 208     accum_type s = 0;
 209     for (int j = 0; j < d; ++j) {
 210       accum_type d1 = p[j] - q[j];
 211       s += d1 * d1;
 212     }
 213     return s;
 214   }
 215 };
 216
 217 template <class H>
 218 class lsh_table {
 219 public:
 220   typedef typename H::scalar_type scalar_type;
 221   typedef typename H::accum_type accum_type;
 222 private:
 223   std::vector<H*> g;
 224   CvLSHOperations* ops;
 225   int d, L, k;
 226   double r;
 227
 228   static accum_type comp_dist(const std::pair<int,accum_type>& x,
 229                   const std::pair<int,accum_type>& y) {
 230     return x.second < y.second;
 231   }
 232
 233   lsh_table(const lsh_table& x);
 234   lsh_table& operator= (const lsh_table& rhs);
 235 public:
 236   lsh_table(CvLSHOperations* _ops, int _d, int Lval, int _k, double _r, CvRNG& rng)
 237     : ops(_ops), d(_d), L(Lval), k(_k), r(_r) {
 238     g.resize(L);
 239     for (int j = 0; j < L; ++j)
 240       g[j] = new H(d, k, r, rng);
 241   }
 242   ~lsh_table() {
 243     for (int j = 0; j < L; ++j)
 244       delete g[j];
 245     delete ops;
 246   }
 247
 248   int dims() const {
 249     return d;
 250   }
 251   unsigned int size() const {
 252     return ops->vector_count();
 253   }
 254
 255   void add(const scalar_type* data, int n, int* ret_indices = 0) {
 256     for (int j=0;j<n;++j) {
 257       const scalar_type* x = data+j*d;
 258       int i = ops->vector_add(x);
 259       if (ret_indices)
 260     ret_indices[j] = i;
 261
 262       for (int l = 0; l < L; ++l) {
 263     lsh_hash h = (*g[l])(x);
 264     ops->hash_insert(h, l, i);
 265       }
 266     }
 267   }
 268   void remove(const int* indices, int n) {
 269     for (int j = 0; j < n; ++j) {
 270       int i = indices[n];
 271       const scalar_type* x = (const scalar_type*)ops->vector_lookup(i);
 272
 273       for (int l = 0; l < L; ++l) {
 274     lsh_hash h = (*g[l])(x);
 275     ops->hash_remove(h, l, i);
 276       }
 277       ops->vector_remove(i);
 278     }
 279   }
 280   void query(const scalar_type* q, int k0, int emax, double* dist, int* results) {
 281     cv::AutoBuffer<int> tmp(emax);
 282     typedef std::pair<int, accum_type> dr_type; // * swap int and accum_type here, for naming consistency
 283     cv::AutoBuffer<dr_type> dr(k0);
 284     int k1 = 0;
 285
 286     // * handle k0 >= emax, in which case don't track max distance
 287
 288     for (int l = 0; l < L && emax > 0; ++l) {
 289       lsh_hash h = (*g[l])(q);
 290       int m = ops->hash_lookup(h, l, tmp, emax);
 291       for (int j = 0; j < m && emax > 0; ++j, --emax) {
 292     int i = tmp[j];
 293     const scalar_type* p = (const scalar_type*)ops->vector_lookup(i);
 294     accum_type pd = (*g[l]).distance(p, q);
 295     if (k1 < k0) {
 296       dr[k1++] = std::make_pair(i, pd);
 297       std::push_heap(&dr[0], &dr[k1], comp_dist);
 298     } else if (pd < dr[0].second) {
 299       std::pop_heap(&dr[0], &dr[k0], comp_dist);
 300       dr[k0 - 1] = std::make_pair(i, pd);
 301       std::push_heap(&dr[0], &dr[k0], comp_dist);
 302     }
 303       }
 304     }
 305
 306     for (int j = 0; j < k1; ++j)
 307       dist[j] = dr[j].second, results[j] = dr[j].first;
 308     std::fill(dist + k1, dist + k0, 0);
 309     std::fill(results + k1, results + k0, -1);
 310   }
 311   void query(const scalar_type* data, int n, int k0, int emax, double* dist, int* results) {
 312     for (int j = 0; j < n; ++j) {
 313       query(data, k0, emax, dist, results);
 314       data += d; // * this may not agree with step for some scalar_types
 315       dist += k0;
 316       results += k0;
 317     }
 318   }
 319 };
 320
 321 typedef lsh_table<pstable_l2_func<float, CV_32FC1> > lsh_pstable_l2_32f;
 322 typedef lsh_table<pstable_l2_func<double, CV_64FC1> > lsh_pstable_l2_64f;
 323
 324 struct CvLSH {
 325   int type;
 326   union {
 327     lsh_pstable_l2_32f* lsh_32f;
 328     lsh_pstable_l2_64f* lsh_64f;
 329   } u;
 330 };
 331
 332 CvLSH* cvCreateLSH(CvLSHOperations* ops, int d, int L, int k, int type, double r, int64 seed) {
 333   CvLSH* lsh = 0;
 334   CvRNG rng = cvRNG(seed);
 335
 336   if (type != CV_32FC1 && type != CV_64FC1)
 337     CV_Error(CV_StsUnsupportedFormat, "vectors must be either CV_32FC1 or CV_64FC1");
 338   lsh = new CvLSH;
 339   try
 340   {
 341       lsh->type = type;
 342       switch (type) {
 343       case CV_32FC1: lsh->u.lsh_32f = new lsh_pstable_l2_32f(ops, d, L, k, r, rng); break;
 344       case CV_64FC1: lsh->u.lsh_64f = new lsh_pstable_l2_64f(ops, d, L, k, r, rng); break;
 345       }
 346   }
 347   catch(...)
 348   {
 349       delete lsh;
 350       throw;
 351   }
 352
 353   return lsh;
 354 }
 355
 356 CvLSH* cvCreateMemoryLSH(int d, int n, int L, int k, int type, double r, int64 seed) {
 357   CvLSHOperations* ops = 0;
 358
 359   switch (type) {
 360   case CV_32FC1: ops = new memory_hash_ops<float>(d,n); break;
 361   case CV_64FC1: ops = new memory_hash_ops<double>(d,n); break;
 362   }
 363   return cvCreateLSH(ops, d, L, k, type, r, seed);
 364 }
 365
 366 void cvReleaseLSH(CvLSH** lsh) {
 367   switch ((*lsh)->type) {
 368   case CV_32FC1: delete (*lsh)->u.lsh_32f; break;
 369   case CV_64FC1: delete (*lsh)->u.lsh_64f; break;
 370   default: assert(0);
 371   }
 372   delete *lsh;
 373   *lsh = 0;
 374 }
 375
 376 unsigned int LSHSize(CvLSH* lsh) {
 377   switch (lsh->type) {
 378   case CV_32FC1: return lsh->u.lsh_32f->size(); break;
 379   case CV_64FC1: return lsh->u.lsh_64f->size(); break;
 380   default: assert(0);
 381   }
 382   return 0;
 383 }
 384
 385
 386 void cvLSHAdd(CvLSH* lsh, const CvMat* data, CvMat* indices) {
 387   int dims, n;
 388   int* ret_indices = 0;
 389
 390   switch (lsh->type) {
 391   case CV_32FC1: dims = lsh->u.lsh_32f->dims(); break;
 392   case CV_64FC1: dims = lsh->u.lsh_64f->dims(); break;
 393   default: assert(0); return;
 394   }
 395
 396   n = data->rows;
 397
 398   if (dims != data->cols)
 399     CV_Error(CV_StsBadSize, "data must be n x d, where d is what was used to construct LSH");
 400
 401   if (CV_MAT_TYPE(data->type) != lsh->type)
 402     CV_Error(CV_StsUnsupportedFormat, "type of data and constructed LSH must agree");
 403   if (indices) {
 404     if (CV_MAT_TYPE(indices->type) != CV_32SC1)
 405       CV_Error(CV_StsUnsupportedFormat, "indices must be CV_32SC1");
 406     if (indices->rows * indices->cols != n)
 407       CV_Error(CV_StsBadSize, "indices must be n x 1 or 1 x n for n x d data");
 408     ret_indices = indices->data.i;
 409   }
 410
 411   switch (lsh->type) {
 412   case CV_32FC1: lsh->u.lsh_32f->add(data->data.fl, n, ret_indices); break;
 413   case CV_64FC1: lsh->u.lsh_64f->add(data->data.db, n, ret_indices); break;
 414   default: assert(0); return;
 415   }
 416 }
 417
 418 void cvLSHRemove(CvLSH* lsh, const CvMat* indices) {
 419   int n;
 420
 421   if (CV_MAT_TYPE(indices->type) != CV_32SC1)
 422     CV_Error(CV_StsUnsupportedFormat, "indices must be CV_32SC1");
 423   n = indices->rows * indices->cols;
 424   switch (lsh->type) {
 425   case CV_32FC1: lsh->u.lsh_32f->remove(indices->data.i, n); break;
 426   case CV_64FC1: lsh->u.lsh_64f->remove(indices->data.i, n); break;
 427   default: assert(0); return;
 428   }
 429 }
 430
 431 void cvLSHQuery(CvLSH* lsh, const CvMat* data, CvMat* indices, CvMat* dist, int k, int emax) {
 432   int dims;
 433
 434   switch (lsh->type) {
 435   case CV_32FC1: dims = lsh->u.lsh_32f->dims(); break;
 436   case CV_64FC1: dims = lsh->u.lsh_64f->dims(); break;
 437   default: assert(0); return;
 438   }
 439
 440   if (k<1)
 441     CV_Error(CV_StsOutOfRange, "k must be positive");
 442   if (CV_MAT_TYPE(data->type) != lsh->type)
 443     CV_Error(CV_StsUnsupportedFormat, "type of data and constructed LSH must agree");
 444   if (dims != data->cols)
 445     CV_Error(CV_StsBadSize, "data must be n x d, where d is what was used to construct LSH");
 446   if (dist->rows != data->rows || dist->cols != k)
 447     CV_Error(CV_StsBadSize, "dist must be n x k for n x d data");
 448   if (dist->rows != indices->rows || dist->cols != indices->cols)
 449     CV_Error(CV_StsBadSize, "dist and indices must be same size");
 450   if (CV_MAT_TYPE(dist->type) != CV_64FC1)
 451     CV_Error(CV_StsUnsupportedFormat, "dist must be CV_64FC1");
 452   if (CV_MAT_TYPE(indices->type) != CV_32SC1)
 453     CV_Error(CV_StsUnsupportedFormat, "indices must be CV_32SC1");
 454
 455   switch (lsh->type) {
 456   case CV_32FC1: lsh->u.lsh_32f->query(data->data.fl, data->rows,
 457                        k, emax, dist->data.db, indices->data.i); break;
 458   case CV_64FC1: lsh->u.lsh_64f->query(data->data.db, data->rows,
 459                        k, emax, dist->data.db, indices->data.i); break;
 460   default: assert(0); return;
 461   }
 462 }