dist.h 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905
  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. #ifndef OPENCV_FLANN_DIST_H_
  31. #define OPENCV_FLANN_DIST_H_
  32. #include <cmath>
  33. #include <cstdlib>
  34. #include <string.h>
  35. #ifdef _MSC_VER
  36. typedef unsigned __int32 uint32_t;
  37. typedef unsigned __int64 uint64_t;
  38. #else
  39. #include <stdint.h>
  40. #endif
  41. #include "defines.h"
  42. #if defined _WIN32 && defined(_M_ARM)
  43. # include <Intrin.h>
  44. #endif
  45. #if defined(__ARM_NEON__) && !defined(__CUDACC__)
  46. # include "arm_neon.h"
  47. #endif
  48. namespace cvflann
  49. {
  50. template<typename T>
  51. inline T abs(T x) { return (x<0) ? -x : x; }
  52. template<>
  53. inline int abs<int>(int x) { return ::abs(x); }
  54. template<>
  55. inline float abs<float>(float x) { return fabsf(x); }
  56. template<>
  57. inline double abs<double>(double x) { return fabs(x); }
  58. template<typename T>
  59. struct Accumulator { typedef T Type; };
  60. template<>
  61. struct Accumulator<unsigned char> { typedef float Type; };
  62. template<>
  63. struct Accumulator<unsigned short> { typedef float Type; };
  64. template<>
  65. struct Accumulator<unsigned int> { typedef float Type; };
  66. template<>
  67. struct Accumulator<char> { typedef float Type; };
  68. template<>
  69. struct Accumulator<short> { typedef float Type; };
  70. template<>
  71. struct Accumulator<int> { typedef float Type; };
  72. #undef True
  73. #undef False
  74. class True
  75. {
  76. };
  77. class False
  78. {
  79. };
  80. /**
  81. * Squared Euclidean distance functor.
  82. *
  83. * This is the simpler, unrolled version. This is preferable for
  84. * very low dimensionality data (eg 3D points)
  85. */
  86. template<class T>
  87. struct L2_Simple
  88. {
  89. typedef True is_kdtree_distance;
  90. typedef True is_vector_space_distance;
  91. typedef T ElementType;
  92. typedef typename Accumulator<T>::Type ResultType;
  93. template <typename Iterator1, typename Iterator2>
  94. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
  95. {
  96. ResultType result = ResultType();
  97. ResultType diff;
  98. for(size_t i = 0; i < size; ++i ) {
  99. diff = *a++ - *b++;
  100. result += diff*diff;
  101. }
  102. return result;
  103. }
  104. template <typename U, typename V>
  105. inline ResultType accum_dist(const U& a, const V& b, int) const
  106. {
  107. return (a-b)*(a-b);
  108. }
  109. };
  110. /**
  111. * Squared Euclidean distance functor, optimized version
  112. */
  113. template<class T>
  114. struct L2
  115. {
  116. typedef True is_kdtree_distance;
  117. typedef True is_vector_space_distance;
  118. typedef T ElementType;
  119. typedef typename Accumulator<T>::Type ResultType;
  120. /**
  121. * Compute the squared Euclidean distance between two vectors.
  122. *
  123. * This is highly optimised, with loop unrolling, as it is one
  124. * of the most expensive inner loops.
  125. *
  126. * The computation of squared root at the end is omitted for
  127. * efficiency.
  128. */
  129. template <typename Iterator1, typename Iterator2>
  130. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
  131. {
  132. ResultType result = ResultType();
  133. ResultType diff0, diff1, diff2, diff3;
  134. Iterator1 last = a + size;
  135. Iterator1 lastgroup = last - 3;
  136. /* Process 4 items with each loop for efficiency. */
  137. while (a < lastgroup) {
  138. diff0 = (ResultType)(a[0] - b[0]);
  139. diff1 = (ResultType)(a[1] - b[1]);
  140. diff2 = (ResultType)(a[2] - b[2]);
  141. diff3 = (ResultType)(a[3] - b[3]);
  142. result += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3;
  143. a += 4;
  144. b += 4;
  145. if ((worst_dist>0)&&(result>worst_dist)) {
  146. return result;
  147. }
  148. }
  149. /* Process last 0-3 pixels. Not needed for standard vector lengths. */
  150. while (a < last) {
  151. diff0 = (ResultType)(*a++ - *b++);
  152. result += diff0 * diff0;
  153. }
  154. return result;
  155. }
  156. /**
  157. * Partial euclidean distance, using just one dimension. This is used by the
  158. * kd-tree when computing partial distances while traversing the tree.
  159. *
  160. * Squared root is omitted for efficiency.
  161. */
  162. template <typename U, typename V>
  163. inline ResultType accum_dist(const U& a, const V& b, int) const
  164. {
  165. return (a-b)*(a-b);
  166. }
  167. };
  168. /*
  169. * Manhattan distance functor, optimized version
  170. */
  171. template<class T>
  172. struct L1
  173. {
  174. typedef True is_kdtree_distance;
  175. typedef True is_vector_space_distance;
  176. typedef T ElementType;
  177. typedef typename Accumulator<T>::Type ResultType;
  178. /**
  179. * Compute the Manhattan (L_1) distance between two vectors.
  180. *
  181. * This is highly optimised, with loop unrolling, as it is one
  182. * of the most expensive inner loops.
  183. */
  184. template <typename Iterator1, typename Iterator2>
  185. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
  186. {
  187. ResultType result = ResultType();
  188. ResultType diff0, diff1, diff2, diff3;
  189. Iterator1 last = a + size;
  190. Iterator1 lastgroup = last - 3;
  191. /* Process 4 items with each loop for efficiency. */
  192. while (a < lastgroup) {
  193. diff0 = (ResultType)abs(a[0] - b[0]);
  194. diff1 = (ResultType)abs(a[1] - b[1]);
  195. diff2 = (ResultType)abs(a[2] - b[2]);
  196. diff3 = (ResultType)abs(a[3] - b[3]);
  197. result += diff0 + diff1 + diff2 + diff3;
  198. a += 4;
  199. b += 4;
  200. if ((worst_dist>0)&&(result>worst_dist)) {
  201. return result;
  202. }
  203. }
  204. /* Process last 0-3 pixels. Not needed for standard vector lengths. */
  205. while (a < last) {
  206. diff0 = (ResultType)abs(*a++ - *b++);
  207. result += diff0;
  208. }
  209. return result;
  210. }
  211. /**
  212. * Partial distance, used by the kd-tree.
  213. */
  214. template <typename U, typename V>
  215. inline ResultType accum_dist(const U& a, const V& b, int) const
  216. {
  217. return abs(a-b);
  218. }
  219. };
  220. template<class T>
  221. struct MinkowskiDistance
  222. {
  223. typedef True is_kdtree_distance;
  224. typedef True is_vector_space_distance;
  225. typedef T ElementType;
  226. typedef typename Accumulator<T>::Type ResultType;
  227. int order;
  228. MinkowskiDistance(int order_) : order(order_) {}
  229. /**
  230. * Compute the Minkowsky (L_p) distance between two vectors.
  231. *
  232. * This is highly optimised, with loop unrolling, as it is one
  233. * of the most expensive inner loops.
  234. *
  235. * The computation of squared root at the end is omitted for
  236. * efficiency.
  237. */
  238. template <typename Iterator1, typename Iterator2>
  239. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
  240. {
  241. ResultType result = ResultType();
  242. ResultType diff0, diff1, diff2, diff3;
  243. Iterator1 last = a + size;
  244. Iterator1 lastgroup = last - 3;
  245. /* Process 4 items with each loop for efficiency. */
  246. while (a < lastgroup) {
  247. diff0 = (ResultType)abs(a[0] - b[0]);
  248. diff1 = (ResultType)abs(a[1] - b[1]);
  249. diff2 = (ResultType)abs(a[2] - b[2]);
  250. diff3 = (ResultType)abs(a[3] - b[3]);
  251. result += pow(diff0,order) + pow(diff1,order) + pow(diff2,order) + pow(diff3,order);
  252. a += 4;
  253. b += 4;
  254. if ((worst_dist>0)&&(result>worst_dist)) {
  255. return result;
  256. }
  257. }
  258. /* Process last 0-3 pixels. Not needed for standard vector lengths. */
  259. while (a < last) {
  260. diff0 = (ResultType)abs(*a++ - *b++);
  261. result += pow(diff0,order);
  262. }
  263. return result;
  264. }
  265. /**
  266. * Partial distance, used by the kd-tree.
  267. */
  268. template <typename U, typename V>
  269. inline ResultType accum_dist(const U& a, const V& b, int) const
  270. {
  271. return pow(static_cast<ResultType>(abs(a-b)),order);
  272. }
  273. };
  274. template<class T>
  275. struct MaxDistance
  276. {
  277. typedef False is_kdtree_distance;
  278. typedef True is_vector_space_distance;
  279. typedef T ElementType;
  280. typedef typename Accumulator<T>::Type ResultType;
  281. /**
  282. * Compute the max distance (L_infinity) between two vectors.
  283. *
  284. * This distance is not a valid kdtree distance, it's not dimensionwise additive.
  285. */
  286. template <typename Iterator1, typename Iterator2>
  287. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
  288. {
  289. ResultType result = ResultType();
  290. ResultType diff0, diff1, diff2, diff3;
  291. Iterator1 last = a + size;
  292. Iterator1 lastgroup = last - 3;
  293. /* Process 4 items with each loop for efficiency. */
  294. while (a < lastgroup) {
  295. diff0 = abs(a[0] - b[0]);
  296. diff1 = abs(a[1] - b[1]);
  297. diff2 = abs(a[2] - b[2]);
  298. diff3 = abs(a[3] - b[3]);
  299. if (diff0>result) {result = diff0; }
  300. if (diff1>result) {result = diff1; }
  301. if (diff2>result) {result = diff2; }
  302. if (diff3>result) {result = diff3; }
  303. a += 4;
  304. b += 4;
  305. if ((worst_dist>0)&&(result>worst_dist)) {
  306. return result;
  307. }
  308. }
  309. /* Process last 0-3 pixels. Not needed for standard vector lengths. */
  310. while (a < last) {
  311. diff0 = abs(*a++ - *b++);
  312. result = (diff0>result) ? diff0 : result;
  313. }
  314. return result;
  315. }
  316. /* This distance functor is not dimension-wise additive, which
  317. * makes it an invalid kd-tree distance, not implementing the accum_dist method */
  318. };
  319. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  320. /**
  321. * Hamming distance functor - counts the bit differences between two strings - useful for the Brief descriptor
  322. * bit count of A exclusive XOR'ed with B
  323. */
  324. struct HammingLUT
  325. {
  326. typedef False is_kdtree_distance;
  327. typedef False is_vector_space_distance;
  328. typedef unsigned char ElementType;
  329. typedef int ResultType;
  330. /** this will count the bits in a ^ b
  331. */
  332. ResultType operator()(const unsigned char* a, const unsigned char* b, size_t size) const
  333. {
  334. static const uchar popCountTable[] =
  335. {
  336. 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  337. 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  338. 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  339. 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  340. 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  341. 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  342. 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  343. 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
  344. };
  345. ResultType result = 0;
  346. for (size_t i = 0; i < size; i++) {
  347. result += popCountTable[a[i] ^ b[i]];
  348. }
  349. return result;
  350. }
  351. };
  352. /**
  353. * Hamming distance functor (pop count between two binary vectors, i.e. xor them and count the number of bits set)
  354. * That code was taken from brief.cpp in OpenCV
  355. */
  356. template<class T>
  357. struct Hamming
  358. {
  359. typedef False is_kdtree_distance;
  360. typedef False is_vector_space_distance;
  361. typedef T ElementType;
  362. typedef int ResultType;
  363. template<typename Iterator1, typename Iterator2>
  364. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
  365. {
  366. ResultType result = 0;
  367. #if defined(__ARM_NEON__) && !defined(__CUDACC__)
  368. {
  369. uint32x4_t bits = vmovq_n_u32(0);
  370. for (size_t i = 0; i < size; i += 16) {
  371. uint8x16_t A_vec = vld1q_u8 (a + i);
  372. uint8x16_t B_vec = vld1q_u8 (b + i);
  373. uint8x16_t AxorB = veorq_u8 (A_vec, B_vec);
  374. uint8x16_t bitsSet = vcntq_u8 (AxorB);
  375. uint16x8_t bitSet8 = vpaddlq_u8 (bitsSet);
  376. uint32x4_t bitSet4 = vpaddlq_u16 (bitSet8);
  377. bits = vaddq_u32(bits, bitSet4);
  378. }
  379. uint64x2_t bitSet2 = vpaddlq_u32 (bits);
  380. result = vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),0);
  381. result += vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),2);
  382. }
  383. #elif __GNUC__
  384. {
  385. //for portability just use unsigned long -- and use the __builtin_popcountll (see docs for __builtin_popcountll)
  386. typedef unsigned long long pop_t;
  387. const size_t modulo = size % sizeof(pop_t);
  388. const pop_t* a2 = reinterpret_cast<const pop_t*> (a);
  389. const pop_t* b2 = reinterpret_cast<const pop_t*> (b);
  390. const pop_t* a2_end = a2 + (size / sizeof(pop_t));
  391. for (; a2 != a2_end; ++a2, ++b2) result += __builtin_popcountll((*a2) ^ (*b2));
  392. if (modulo) {
  393. //in the case where size is not dividable by sizeof(size_t)
  394. //need to mask off the bits at the end
  395. pop_t a_final = 0, b_final = 0;
  396. memcpy(&a_final, a2, modulo);
  397. memcpy(&b_final, b2, modulo);
  398. result += __builtin_popcountll(a_final ^ b_final);
  399. }
  400. }
  401. #else // NO NEON and NOT GNUC
  402. HammingLUT lut;
  403. result = lut(reinterpret_cast<const unsigned char*> (a),
  404. reinterpret_cast<const unsigned char*> (b), size);
  405. #endif
  406. return result;
  407. }
  408. };
  409. template<typename T>
  410. struct Hamming2
  411. {
  412. typedef False is_kdtree_distance;
  413. typedef False is_vector_space_distance;
  414. typedef T ElementType;
  415. typedef int ResultType;
  416. /** This is popcount_3() from:
  417. * http://en.wikipedia.org/wiki/Hamming_weight */
  418. unsigned int popcnt32(uint32_t n) const
  419. {
  420. n -= ((n >> 1) & 0x55555555);
  421. n = (n & 0x33333333) + ((n >> 2) & 0x33333333);
  422. return (((n + (n >> 4))& 0xF0F0F0F)* 0x1010101) >> 24;
  423. }
  424. #ifdef FLANN_PLATFORM_64_BIT
  425. unsigned int popcnt64(uint64_t n) const
  426. {
  427. n -= ((n >> 1) & 0x5555555555555555);
  428. n = (n & 0x3333333333333333) + ((n >> 2) & 0x3333333333333333);
  429. return (((n + (n >> 4))& 0x0f0f0f0f0f0f0f0f)* 0x0101010101010101) >> 56;
  430. }
  431. #endif
  432. template <typename Iterator1, typename Iterator2>
  433. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
  434. {
  435. #ifdef FLANN_PLATFORM_64_BIT
  436. const uint64_t* pa = reinterpret_cast<const uint64_t*>(a);
  437. const uint64_t* pb = reinterpret_cast<const uint64_t*>(b);
  438. ResultType result = 0;
  439. size /= (sizeof(uint64_t)/sizeof(unsigned char));
  440. for(size_t i = 0; i < size; ++i ) {
  441. result += popcnt64(*pa ^ *pb);
  442. ++pa;
  443. ++pb;
  444. }
  445. #else
  446. const uint32_t* pa = reinterpret_cast<const uint32_t*>(a);
  447. const uint32_t* pb = reinterpret_cast<const uint32_t*>(b);
  448. ResultType result = 0;
  449. size /= (sizeof(uint32_t)/sizeof(unsigned char));
  450. for(size_t i = 0; i < size; ++i ) {
  451. result += popcnt32(*pa ^ *pb);
  452. ++pa;
  453. ++pb;
  454. }
  455. #endif
  456. return result;
  457. }
  458. };
  459. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  460. template<class T>
  461. struct HistIntersectionDistance
  462. {
  463. typedef True is_kdtree_distance;
  464. typedef True is_vector_space_distance;
  465. typedef T ElementType;
  466. typedef typename Accumulator<T>::Type ResultType;
  467. /**
  468. * Compute the histogram intersection distance
  469. */
  470. template <typename Iterator1, typename Iterator2>
  471. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
  472. {
  473. ResultType result = ResultType();
  474. ResultType min0, min1, min2, min3;
  475. Iterator1 last = a + size;
  476. Iterator1 lastgroup = last - 3;
  477. /* Process 4 items with each loop for efficiency. */
  478. while (a < lastgroup) {
  479. min0 = (ResultType)(a[0] < b[0] ? a[0] : b[0]);
  480. min1 = (ResultType)(a[1] < b[1] ? a[1] : b[1]);
  481. min2 = (ResultType)(a[2] < b[2] ? a[2] : b[2]);
  482. min3 = (ResultType)(a[3] < b[3] ? a[3] : b[3]);
  483. result += min0 + min1 + min2 + min3;
  484. a += 4;
  485. b += 4;
  486. if ((worst_dist>0)&&(result>worst_dist)) {
  487. return result;
  488. }
  489. }
  490. /* Process last 0-3 pixels. Not needed for standard vector lengths. */
  491. while (a < last) {
  492. min0 = (ResultType)(*a < *b ? *a : *b);
  493. result += min0;
  494. ++a;
  495. ++b;
  496. }
  497. return result;
  498. }
  499. /**
  500. * Partial distance, used by the kd-tree.
  501. */
  502. template <typename U, typename V>
  503. inline ResultType accum_dist(const U& a, const V& b, int) const
  504. {
  505. return a<b ? a : b;
  506. }
  507. };
  508. template<class T>
  509. struct HellingerDistance
  510. {
  511. typedef True is_kdtree_distance;
  512. typedef True is_vector_space_distance;
  513. typedef T ElementType;
  514. typedef typename Accumulator<T>::Type ResultType;
  515. /**
  516. * Compute the Hellinger distance
  517. */
  518. template <typename Iterator1, typename Iterator2>
  519. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
  520. {
  521. ResultType result = ResultType();
  522. ResultType diff0, diff1, diff2, diff3;
  523. Iterator1 last = a + size;
  524. Iterator1 lastgroup = last - 3;
  525. /* Process 4 items with each loop for efficiency. */
  526. while (a < lastgroup) {
  527. diff0 = sqrt(static_cast<ResultType>(a[0])) - sqrt(static_cast<ResultType>(b[0]));
  528. diff1 = sqrt(static_cast<ResultType>(a[1])) - sqrt(static_cast<ResultType>(b[1]));
  529. diff2 = sqrt(static_cast<ResultType>(a[2])) - sqrt(static_cast<ResultType>(b[2]));
  530. diff3 = sqrt(static_cast<ResultType>(a[3])) - sqrt(static_cast<ResultType>(b[3]));
  531. result += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3;
  532. a += 4;
  533. b += 4;
  534. }
  535. while (a < last) {
  536. diff0 = sqrt(static_cast<ResultType>(*a++)) - sqrt(static_cast<ResultType>(*b++));
  537. result += diff0 * diff0;
  538. }
  539. return result;
  540. }
  541. /**
  542. * Partial distance, used by the kd-tree.
  543. */
  544. template <typename U, typename V>
  545. inline ResultType accum_dist(const U& a, const V& b, int) const
  546. {
  547. ResultType diff = sqrt(static_cast<ResultType>(a)) - sqrt(static_cast<ResultType>(b));
  548. return diff * diff;
  549. }
  550. };
  551. template<class T>
  552. struct ChiSquareDistance
  553. {
  554. typedef True is_kdtree_distance;
  555. typedef True is_vector_space_distance;
  556. typedef T ElementType;
  557. typedef typename Accumulator<T>::Type ResultType;
  558. /**
  559. * Compute the chi-square distance
  560. */
  561. template <typename Iterator1, typename Iterator2>
  562. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
  563. {
  564. ResultType result = ResultType();
  565. ResultType sum, diff;
  566. Iterator1 last = a + size;
  567. while (a < last) {
  568. sum = (ResultType)(*a + *b);
  569. if (sum>0) {
  570. diff = (ResultType)(*a - *b);
  571. result += diff*diff/sum;
  572. }
  573. ++a;
  574. ++b;
  575. if ((worst_dist>0)&&(result>worst_dist)) {
  576. return result;
  577. }
  578. }
  579. return result;
  580. }
  581. /**
  582. * Partial distance, used by the kd-tree.
  583. */
  584. template <typename U, typename V>
  585. inline ResultType accum_dist(const U& a, const V& b, int) const
  586. {
  587. ResultType result = ResultType();
  588. ResultType sum, diff;
  589. sum = (ResultType)(a+b);
  590. if (sum>0) {
  591. diff = (ResultType)(a-b);
  592. result = diff*diff/sum;
  593. }
  594. return result;
  595. }
  596. };
  597. template<class T>
  598. struct KL_Divergence
  599. {
  600. typedef True is_kdtree_distance;
  601. typedef True is_vector_space_distance;
  602. typedef T ElementType;
  603. typedef typename Accumulator<T>::Type ResultType;
  604. /**
  605. * Compute the Kullback-Leibler divergence
  606. */
  607. template <typename Iterator1, typename Iterator2>
  608. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
  609. {
  610. ResultType result = ResultType();
  611. Iterator1 last = a + size;
  612. while (a < last) {
  613. if (* b != 0) {
  614. ResultType ratio = (ResultType)(*a / *b);
  615. if (ratio>0) {
  616. result += *a * log(ratio);
  617. }
  618. }
  619. ++a;
  620. ++b;
  621. if ((worst_dist>0)&&(result>worst_dist)) {
  622. return result;
  623. }
  624. }
  625. return result;
  626. }
  627. /**
  628. * Partial distance, used by the kd-tree.
  629. */
  630. template <typename U, typename V>
  631. inline ResultType accum_dist(const U& a, const V& b, int) const
  632. {
  633. ResultType result = ResultType();
  634. if( *b != 0 ) {
  635. ResultType ratio = (ResultType)(a / b);
  636. if (ratio>0) {
  637. result = a * log(ratio);
  638. }
  639. }
  640. return result;
  641. }
  642. };
  643. /*
  644. * This is a "zero iterator". It basically behaves like a zero filled
  645. * array to all algorithms that use arrays as iterators (STL style).
  646. * It's useful when there's a need to compute the distance between feature
  647. * and origin it and allows for better compiler optimisation than using a
  648. * zero-filled array.
  649. */
  650. template <typename T>
  651. struct ZeroIterator
  652. {
  653. T operator*()
  654. {
  655. return 0;
  656. }
  657. T operator[](int)
  658. {
  659. return 0;
  660. }
  661. const ZeroIterator<T>& operator ++()
  662. {
  663. return *this;
  664. }
  665. ZeroIterator<T> operator ++(int)
  666. {
  667. return *this;
  668. }
  669. ZeroIterator<T>& operator+=(int)
  670. {
  671. return *this;
  672. }
  673. };
  674. /*
  675. * Depending on processed distances, some of them are already squared (e.g. L2)
  676. * and some are not (e.g.Hamming). In KMeans++ for instance we want to be sure
  677. * we are working on ^2 distances, thus following templates to ensure that.
  678. */
  679. template <typename Distance, typename ElementType>
  680. struct squareDistance
  681. {
  682. typedef typename Distance::ResultType ResultType;
  683. ResultType operator()( ResultType dist ) { return dist*dist; }
  684. };
  685. template <typename ElementType>
  686. struct squareDistance<L2_Simple<ElementType>, ElementType>
  687. {
  688. typedef typename L2_Simple<ElementType>::ResultType ResultType;
  689. ResultType operator()( ResultType dist ) { return dist; }
  690. };
  691. template <typename ElementType>
  692. struct squareDistance<L2<ElementType>, ElementType>
  693. {
  694. typedef typename L2<ElementType>::ResultType ResultType;
  695. ResultType operator()( ResultType dist ) { return dist; }
  696. };
  697. template <typename ElementType>
  698. struct squareDistance<MinkowskiDistance<ElementType>, ElementType>
  699. {
  700. typedef typename MinkowskiDistance<ElementType>::ResultType ResultType;
  701. ResultType operator()( ResultType dist ) { return dist; }
  702. };
  703. template <typename ElementType>
  704. struct squareDistance<HellingerDistance<ElementType>, ElementType>
  705. {
  706. typedef typename HellingerDistance<ElementType>::ResultType ResultType;
  707. ResultType operator()( ResultType dist ) { return dist; }
  708. };
  709. template <typename ElementType>
  710. struct squareDistance<ChiSquareDistance<ElementType>, ElementType>
  711. {
  712. typedef typename ChiSquareDistance<ElementType>::ResultType ResultType;
  713. ResultType operator()( ResultType dist ) { return dist; }
  714. };
  715. template <typename Distance>
  716. typename Distance::ResultType ensureSquareDistance( typename Distance::ResultType dist )
  717. {
  718. typedef typename Distance::ElementType ElementType;
  719. squareDistance<Distance, ElementType> dummy;
  720. return dummy( dist );
  721. }
  722. /*
  723. * ...and a template to ensure the user that he will process the normal distance,
  724. * and not squared distance, without losing processing time calling sqrt(ensureSquareDistance)
  725. * that will result in doing actually sqrt(dist*dist) for L1 distance for instance.
  726. */
  727. template <typename Distance, typename ElementType>
  728. struct simpleDistance
  729. {
  730. typedef typename Distance::ResultType ResultType;
  731. ResultType operator()( ResultType dist ) { return dist; }
  732. };
  733. template <typename ElementType>
  734. struct simpleDistance<L2_Simple<ElementType>, ElementType>
  735. {
  736. typedef typename L2_Simple<ElementType>::ResultType ResultType;
  737. ResultType operator()( ResultType dist ) { return sqrt(dist); }
  738. };
  739. template <typename ElementType>
  740. struct simpleDistance<L2<ElementType>, ElementType>
  741. {
  742. typedef typename L2<ElementType>::ResultType ResultType;
  743. ResultType operator()( ResultType dist ) { return sqrt(dist); }
  744. };
  745. template <typename ElementType>
  746. struct simpleDistance<MinkowskiDistance<ElementType>, ElementType>
  747. {
  748. typedef typename MinkowskiDistance<ElementType>::ResultType ResultType;
  749. ResultType operator()( ResultType dist ) { return sqrt(dist); }
  750. };
  751. template <typename ElementType>
  752. struct simpleDistance<HellingerDistance<ElementType>, ElementType>
  753. {
  754. typedef typename HellingerDistance<ElementType>::ResultType ResultType;
  755. ResultType operator()( ResultType dist ) { return sqrt(dist); }
  756. };
  757. template <typename ElementType>
  758. struct simpleDistance<ChiSquareDistance<ElementType>, ElementType>
  759. {
  760. typedef typename ChiSquareDistance<ElementType>::ResultType ResultType;
  761. ResultType operator()( ResultType dist ) { return sqrt(dist); }
  762. };
  763. template <typename Distance>
  764. typename Distance::ResultType ensureSimpleDistance( typename Distance::ResultType dist )
  765. {
  766. typedef typename Distance::ElementType ElementType;
  767. simpleDistance<Distance, ElementType> dummy;
  768. return dummy( dist );
  769. }
  770. }
  771. #endif //OPENCV_FLANN_DIST_H_