lsh_table.h 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492
  1. /***********************************************************************
  2. * Software License Agreement (BSD License)
  3. *
  4. * Copyright 2008-2009 Marius Muja (mariusm@cs.ubc.ca). All rights reserved.
  5. * Copyright 2008-2009 David G. Lowe (lowe@cs.ubc.ca). All rights reserved.
  6. *
  7. * THE BSD LICENSE
  8. *
  9. * Redistribution and use in source and binary forms, with or without
  10. * modification, are permitted provided that the following conditions
  11. * are met:
  12. *
  13. * 1. Redistributions of source code must retain the above copyright
  14. * notice, this list of conditions and the following disclaimer.
  15. * 2. Redistributions in binary form must reproduce the above copyright
  16. * notice, this list of conditions and the following disclaimer in the
  17. * documentation and/or other materials provided with the distribution.
  18. *
  19. * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  20. * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
  21. * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  22. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
  23. * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  24. * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  25. * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  26. * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  27. * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  28. * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  29. *************************************************************************/
  30. /***********************************************************************
  31. * Author: Vincent Rabaud
  32. *************************************************************************/
  33. #ifndef OPENCV_FLANN_LSH_TABLE_H_
  34. #define OPENCV_FLANN_LSH_TABLE_H_
  35. #include <algorithm>
  36. #include <iostream>
  37. #include <iomanip>
  38. #include <limits.h>
  39. // TODO as soon as we use C++0x, use the code in USE_UNORDERED_MAP
  40. #ifdef __GXX_EXPERIMENTAL_CXX0X__
  41. # define USE_UNORDERED_MAP 1
  42. #else
  43. # define USE_UNORDERED_MAP 0
  44. #endif
  45. #if USE_UNORDERED_MAP
  46. #include <unordered_map>
  47. #else
  48. #include <map>
  49. #endif
  50. #include <math.h>
  51. #include <stddef.h>
  52. #include "dynamic_bitset.h"
  53. #include "matrix.h"
  54. namespace cvflann
  55. {
  56. namespace lsh
  57. {
  58. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  59. /** What is stored in an LSH bucket
  60. */
  61. typedef uint32_t FeatureIndex;
  62. /** The id from which we can get a bucket back in an LSH table
  63. */
  64. typedef unsigned int BucketKey;
  65. /** A bucket in an LSH table
  66. */
  67. typedef std::vector<FeatureIndex> Bucket;
  68. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  69. /** POD for stats about an LSH table
  70. */
  71. struct LshStats
  72. {
  73. std::vector<unsigned int> bucket_sizes_;
  74. size_t n_buckets_;
  75. size_t bucket_size_mean_;
  76. size_t bucket_size_median_;
  77. size_t bucket_size_min_;
  78. size_t bucket_size_max_;
  79. size_t bucket_size_std_dev;
  80. /** Each contained vector contains three value: beginning/end for interval, number of elements in the bin
  81. */
  82. std::vector<std::vector<unsigned int> > size_histogram_;
  83. };
  84. /** Overload the << operator for LshStats
  85. * @param out the streams
  86. * @param stats the stats to display
  87. * @return the streams
  88. */
  89. inline std::ostream& operator <<(std::ostream& out, const LshStats& stats)
  90. {
  91. int w = 20;
  92. out << "Lsh Table Stats:\n" << std::setw(w) << std::setiosflags(std::ios::right) << "N buckets : "
  93. << stats.n_buckets_ << "\n" << std::setw(w) << std::setiosflags(std::ios::right) << "mean size : "
  94. << std::setiosflags(std::ios::left) << stats.bucket_size_mean_ << "\n" << std::setw(w)
  95. << std::setiosflags(std::ios::right) << "median size : " << stats.bucket_size_median_ << "\n" << std::setw(w)
  96. << std::setiosflags(std::ios::right) << "min size : " << std::setiosflags(std::ios::left)
  97. << stats.bucket_size_min_ << "\n" << std::setw(w) << std::setiosflags(std::ios::right) << "max size : "
  98. << std::setiosflags(std::ios::left) << stats.bucket_size_max_;
  99. // Display the histogram
  100. out << std::endl << std::setw(w) << std::setiosflags(std::ios::right) << "histogram : "
  101. << std::setiosflags(std::ios::left);
  102. for (std::vector<std::vector<unsigned int> >::const_iterator iterator = stats.size_histogram_.begin(), end =
  103. stats.size_histogram_.end(); iterator != end; ++iterator) out << (*iterator)[0] << "-" << (*iterator)[1] << ": " << (*iterator)[2] << ", ";
  104. return out;
  105. }
  106. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  107. /** Lsh hash table. As its key is a sub-feature, and as usually
  108. * the size of it is pretty small, we keep it as a continuous memory array.
  109. * The value is an index in the corpus of features (we keep it as an unsigned
  110. * int for pure memory reasons, it could be a size_t)
  111. */
  112. template<typename ElementType>
  113. class LshTable
  114. {
  115. public:
  116. /** A container of all the feature indices. Optimized for space
  117. */
  118. #if USE_UNORDERED_MAP
  119. typedef std::unordered_map<BucketKey, Bucket> BucketsSpace;
  120. #else
  121. typedef std::map<BucketKey, Bucket> BucketsSpace;
  122. #endif
  123. /** A container of all the feature indices. Optimized for speed
  124. */
  125. typedef std::vector<Bucket> BucketsSpeed;
  126. /** Default constructor
  127. */
  128. LshTable()
  129. {
  130. }
  131. /** Default constructor
  132. * Create the mask and allocate the memory
  133. * @param feature_size is the size of the feature (considered as a ElementType[])
  134. * @param key_size is the number of bits that are turned on in the feature
  135. */
  136. LshTable(unsigned int feature_size, unsigned int key_size)
  137. {
  138. (void)feature_size;
  139. (void)key_size;
  140. std::cerr << "LSH is not implemented for that type" << std::endl;
  141. assert(0);
  142. }
  143. /** Add a feature to the table
  144. * @param value the value to store for that feature
  145. * @param feature the feature itself
  146. */
  147. void add(unsigned int value, const ElementType* feature)
  148. {
  149. // Add the value to the corresponding bucket
  150. BucketKey key = (lsh::BucketKey)getKey(feature);
  151. switch (speed_level_) {
  152. case kArray:
  153. // That means we get the buckets from an array
  154. buckets_speed_[key].push_back(value);
  155. break;
  156. case kBitsetHash:
  157. // That means we can check the bitset for the presence of a key
  158. key_bitset_.set(key);
  159. buckets_space_[key].push_back(value);
  160. break;
  161. case kHash:
  162. {
  163. // That means we have to check for the hash table for the presence of a key
  164. buckets_space_[key].push_back(value);
  165. break;
  166. }
  167. }
  168. }
  169. /** Add a set of features to the table
  170. * @param dataset the values to store
  171. */
  172. void add(Matrix<ElementType> dataset)
  173. {
  174. #if USE_UNORDERED_MAP
  175. buckets_space_.rehash((buckets_space_.size() + dataset.rows) * 1.2);
  176. #endif
  177. // Add the features to the table
  178. for (unsigned int i = 0; i < dataset.rows; ++i) add(i, dataset[i]);
  179. // Now that the table is full, optimize it for speed/space
  180. optimize();
  181. }
  182. /** Get a bucket given the key
  183. * @param key
  184. * @return
  185. */
  186. inline const Bucket* getBucketFromKey(BucketKey key) const
  187. {
  188. // Generate other buckets
  189. switch (speed_level_) {
  190. case kArray:
  191. // That means we get the buckets from an array
  192. return &buckets_speed_[key];
  193. break;
  194. case kBitsetHash:
  195. // That means we can check the bitset for the presence of a key
  196. if (key_bitset_.test(key)) return &buckets_space_.find(key)->second;
  197. else return 0;
  198. break;
  199. case kHash:
  200. {
  201. // That means we have to check for the hash table for the presence of a key
  202. BucketsSpace::const_iterator bucket_it, bucket_end = buckets_space_.end();
  203. bucket_it = buckets_space_.find(key);
  204. // Stop here if that bucket does not exist
  205. if (bucket_it == bucket_end) return 0;
  206. else return &bucket_it->second;
  207. break;
  208. }
  209. }
  210. return 0;
  211. }
  212. /** Compute the sub-signature of a feature
  213. */
  214. size_t getKey(const ElementType* /*feature*/) const
  215. {
  216. std::cerr << "LSH is not implemented for that type" << std::endl;
  217. assert(0);
  218. return 1;
  219. }
  220. /** Get statistics about the table
  221. * @return
  222. */
  223. LshStats getStats() const;
  224. private:
  225. /** defines the speed fo the implementation
  226. * kArray uses a vector for storing data
  227. * kBitsetHash uses a hash map but checks for the validity of a key with a bitset
  228. * kHash uses a hash map only
  229. */
  230. enum SpeedLevel
  231. {
  232. kArray, kBitsetHash, kHash
  233. };
  234. /** Initialize some variables
  235. */
  236. void initialize(size_t key_size)
  237. {
  238. const size_t key_size_lower_bound = 1;
  239. //a value (size_t(1) << key_size) must fit the size_t type so key_size has to be strictly less than size of size_t
  240. const size_t key_size_upper_bound = std::min(sizeof(BucketKey) * CHAR_BIT + 1, sizeof(size_t) * CHAR_BIT);
  241. if (key_size < key_size_lower_bound || key_size >= key_size_upper_bound)
  242. {
  243. CV_Error(cv::Error::StsBadArg, cv::format("Invalid key_size (=%d). Valid values for your system are %d <= key_size < %d.", (int)key_size, (int)key_size_lower_bound, (int)key_size_upper_bound));
  244. }
  245. speed_level_ = kHash;
  246. key_size_ = (unsigned)key_size;
  247. }
  248. /** Optimize the table for speed/space
  249. */
  250. void optimize()
  251. {
  252. // If we are already using the fast storage, no need to do anything
  253. if (speed_level_ == kArray) return;
  254. // Use an array if it will be more than half full
  255. if (buckets_space_.size() > ((size_t(1) << key_size_) / 2)) {
  256. speed_level_ = kArray;
  257. // Fill the array version of it
  258. buckets_speed_.resize(size_t(1) << key_size_);
  259. for (BucketsSpace::const_iterator key_bucket = buckets_space_.begin(); key_bucket != buckets_space_.end(); ++key_bucket) buckets_speed_[key_bucket->first] = key_bucket->second;
  260. // Empty the hash table
  261. buckets_space_.clear();
  262. return;
  263. }
  264. // If the bitset is going to use less than 10% of the RAM of the hash map (at least 1 size_t for the key and two
  265. // for the vector) or less than 512MB (key_size_ <= 30)
  266. if (((std::max(buckets_space_.size(), buckets_speed_.size()) * CHAR_BIT * 3 * sizeof(BucketKey)) / 10
  267. >= (size_t(1) << key_size_)) || (key_size_ <= 32)) {
  268. speed_level_ = kBitsetHash;
  269. key_bitset_.resize(size_t(1) << key_size_);
  270. key_bitset_.reset();
  271. // Try with the BucketsSpace
  272. for (BucketsSpace::const_iterator key_bucket = buckets_space_.begin(); key_bucket != buckets_space_.end(); ++key_bucket) key_bitset_.set(key_bucket->first);
  273. }
  274. else {
  275. speed_level_ = kHash;
  276. key_bitset_.clear();
  277. }
  278. }
  279. /** The vector of all the buckets if they are held for speed
  280. */
  281. BucketsSpeed buckets_speed_;
  282. /** The hash table of all the buckets in case we cannot use the speed version
  283. */
  284. BucketsSpace buckets_space_;
  285. /** What is used to store the data */
  286. SpeedLevel speed_level_;
  287. /** If the subkey is small enough, it will keep track of which subkeys are set through that bitset
  288. * That is just a speedup so that we don't look in the hash table (which can be mush slower that checking a bitset)
  289. */
  290. DynamicBitset key_bitset_;
  291. /** The size of the sub-signature in bits
  292. */
  293. unsigned int key_size_;
  294. // Members only used for the unsigned char specialization
  295. /** The mask to apply to a feature to get the hash key
  296. * Only used in the unsigned char case
  297. */
  298. std::vector<size_t> mask_;
  299. };
  300. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  301. // Specialization for unsigned char
  302. template<>
  303. inline LshTable<unsigned char>::LshTable(unsigned int feature_size, unsigned int subsignature_size)
  304. {
  305. initialize(subsignature_size);
  306. // Allocate the mask
  307. mask_ = std::vector<size_t>((size_t)ceil((float)(feature_size * sizeof(char)) / (float)sizeof(size_t)), 0);
  308. // A bit brutal but fast to code
  309. std::vector<size_t> indices(feature_size * CHAR_BIT);
  310. for (size_t i = 0; i < feature_size * CHAR_BIT; ++i) indices[i] = i;
  311. std::random_shuffle(indices.begin(), indices.end());
  312. // Generate a random set of order of subsignature_size_ bits
  313. for (unsigned int i = 0; i < key_size_; ++i) {
  314. size_t index = indices[i];
  315. // Set that bit in the mask
  316. size_t divisor = CHAR_BIT * sizeof(size_t);
  317. size_t idx = index / divisor; //pick the right size_t index
  318. mask_[idx] |= size_t(1) << (index % divisor); //use modulo to find the bit offset
  319. }
  320. // Set to 1 if you want to display the mask for debug
  321. #if 0
  322. {
  323. size_t bcount = 0;
  324. BOOST_FOREACH(size_t mask_block, mask_){
  325. out << std::setw(sizeof(size_t) * CHAR_BIT / 4) << std::setfill('0') << std::hex << mask_block
  326. << std::endl;
  327. bcount += __builtin_popcountll(mask_block);
  328. }
  329. out << "bit count : " << std::dec << bcount << std::endl;
  330. out << "mask size : " << mask_.size() << std::endl;
  331. return out;
  332. }
  333. #endif
  334. }
  335. /** Return the Subsignature of a feature
  336. * @param feature the feature to analyze
  337. */
  338. template<>
  339. inline size_t LshTable<unsigned char>::getKey(const unsigned char* feature) const
  340. {
  341. // no need to check if T is dividable by sizeof(size_t) like in the Hamming
  342. // distance computation as we have a mask
  343. const size_t* feature_block_ptr = reinterpret_cast<const size_t*> ((const void*)feature);
  344. // Figure out the subsignature of the feature
  345. // Given the feature ABCDEF, and the mask 001011, the output will be
  346. // 000CEF
  347. size_t subsignature = 0;
  348. size_t bit_index = 1;
  349. for (std::vector<size_t>::const_iterator pmask_block = mask_.begin(); pmask_block != mask_.end(); ++pmask_block) {
  350. // get the mask and signature blocks
  351. size_t feature_block = *feature_block_ptr;
  352. size_t mask_block = *pmask_block;
  353. while (mask_block) {
  354. // Get the lowest set bit in the mask block
  355. size_t lowest_bit = mask_block & (-(ptrdiff_t)mask_block);
  356. // Add it to the current subsignature if necessary
  357. subsignature += (feature_block & lowest_bit) ? bit_index : 0;
  358. // Reset the bit in the mask block
  359. mask_block ^= lowest_bit;
  360. // increment the bit index for the subsignature
  361. bit_index <<= 1;
  362. }
  363. // Check the next feature block
  364. ++feature_block_ptr;
  365. }
  366. return subsignature;
  367. }
  368. template<>
  369. inline LshStats LshTable<unsigned char>::getStats() const
  370. {
  371. LshStats stats;
  372. stats.bucket_size_mean_ = 0;
  373. if ((buckets_speed_.empty()) && (buckets_space_.empty())) {
  374. stats.n_buckets_ = 0;
  375. stats.bucket_size_median_ = 0;
  376. stats.bucket_size_min_ = 0;
  377. stats.bucket_size_max_ = 0;
  378. return stats;
  379. }
  380. if (!buckets_speed_.empty()) {
  381. for (BucketsSpeed::const_iterator pbucket = buckets_speed_.begin(); pbucket != buckets_speed_.end(); ++pbucket) {
  382. stats.bucket_sizes_.push_back((lsh::FeatureIndex)pbucket->size());
  383. stats.bucket_size_mean_ += pbucket->size();
  384. }
  385. stats.bucket_size_mean_ /= buckets_speed_.size();
  386. stats.n_buckets_ = buckets_speed_.size();
  387. }
  388. else {
  389. for (BucketsSpace::const_iterator x = buckets_space_.begin(); x != buckets_space_.end(); ++x) {
  390. stats.bucket_sizes_.push_back((lsh::FeatureIndex)x->second.size());
  391. stats.bucket_size_mean_ += x->second.size();
  392. }
  393. stats.bucket_size_mean_ /= buckets_space_.size();
  394. stats.n_buckets_ = buckets_space_.size();
  395. }
  396. std::sort(stats.bucket_sizes_.begin(), stats.bucket_sizes_.end());
  397. // BOOST_FOREACH(int size, stats.bucket_sizes_)
  398. // std::cout << size << " ";
  399. // std::cout << std::endl;
  400. stats.bucket_size_median_ = stats.bucket_sizes_[stats.bucket_sizes_.size() / 2];
  401. stats.bucket_size_min_ = stats.bucket_sizes_.front();
  402. stats.bucket_size_max_ = stats.bucket_sizes_.back();
  403. // TODO compute mean and std
  404. /*float mean, stddev;
  405. stats.bucket_size_mean_ = mean;
  406. stats.bucket_size_std_dev = stddev;*/
  407. // Include a histogram of the buckets
  408. unsigned int bin_start = 0;
  409. unsigned int bin_end = 20;
  410. bool is_new_bin = true;
  411. for (std::vector<unsigned int>::iterator iterator = stats.bucket_sizes_.begin(), end = stats.bucket_sizes_.end(); iterator
  412. != end; )
  413. if (*iterator < bin_end) {
  414. if (is_new_bin) {
  415. stats.size_histogram_.push_back(std::vector<unsigned int>(3, 0));
  416. stats.size_histogram_.back()[0] = bin_start;
  417. stats.size_histogram_.back()[1] = bin_end - 1;
  418. is_new_bin = false;
  419. }
  420. ++stats.size_histogram_.back()[2];
  421. ++iterator;
  422. }
  423. else {
  424. bin_start += 20;
  425. bin_end += 20;
  426. is_new_bin = true;
  427. }
  428. return stats;
  429. }
  430. // End the two namespaces
  431. }
  432. }
  433. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  434. #endif /* OPENCV_FLANN_LSH_TABLE_H_ */