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. #ifdef __ARM_NEON__
  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. #ifdef __ARM_NEON__
  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. typedef unsigned long long pop_t;
  403. HammingLUT lut;
  404. result = lut(reinterpret_cast<const unsigned char*> (a),
  405. reinterpret_cast<const unsigned char*> (b), size * sizeof(pop_t));
  406. #endif
  407. return result;
  408. }
  409. };
  410. template<typename T>
  411. struct Hamming2
  412. {
  413. typedef False is_kdtree_distance;
  414. typedef False is_vector_space_distance;
  415. typedef T ElementType;
  416. typedef int ResultType;
  417. /** This is popcount_3() from:
  418. * http://en.wikipedia.org/wiki/Hamming_weight */
  419. unsigned int popcnt32(uint32_t n) const
  420. {
  421. n -= ((n >> 1) & 0x55555555);
  422. n = (n & 0x33333333) + ((n >> 2) & 0x33333333);
  423. return (((n + (n >> 4))& 0xF0F0F0F)* 0x1010101) >> 24;
  424. }
  425. #ifdef FLANN_PLATFORM_64_BIT
  426. unsigned int popcnt64(uint64_t n) const
  427. {
  428. n -= ((n >> 1) & 0x5555555555555555);
  429. n = (n & 0x3333333333333333) + ((n >> 2) & 0x3333333333333333);
  430. return (((n + (n >> 4))& 0x0f0f0f0f0f0f0f0f)* 0x0101010101010101) >> 56;
  431. }
  432. #endif
  433. template <typename Iterator1, typename Iterator2>
  434. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
  435. {
  436. #ifdef FLANN_PLATFORM_64_BIT
  437. const uint64_t* pa = reinterpret_cast<const uint64_t*>(a);
  438. const uint64_t* pb = reinterpret_cast<const uint64_t*>(b);
  439. ResultType result = 0;
  440. size /= (sizeof(uint64_t)/sizeof(unsigned char));
  441. for(size_t i = 0; i < size; ++i ) {
  442. result += popcnt64(*pa ^ *pb);
  443. ++pa;
  444. ++pb;
  445. }
  446. #else
  447. const uint32_t* pa = reinterpret_cast<const uint32_t*>(a);
  448. const uint32_t* pb = reinterpret_cast<const uint32_t*>(b);
  449. ResultType result = 0;
  450. size /= (sizeof(uint32_t)/sizeof(unsigned char));
  451. for(size_t i = 0; i < size; ++i ) {
  452. result += popcnt32(*pa ^ *pb);
  453. ++pa;
  454. ++pb;
  455. }
  456. #endif
  457. return result;
  458. }
  459. };
  460. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  461. template<class T>
  462. struct HistIntersectionDistance
  463. {
  464. typedef True is_kdtree_distance;
  465. typedef True is_vector_space_distance;
  466. typedef T ElementType;
  467. typedef typename Accumulator<T>::Type ResultType;
  468. /**
  469. * Compute the histogram intersection distance
  470. */
  471. template <typename Iterator1, typename Iterator2>
  472. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
  473. {
  474. ResultType result = ResultType();
  475. ResultType min0, min1, min2, min3;
  476. Iterator1 last = a + size;
  477. Iterator1 lastgroup = last - 3;
  478. /* Process 4 items with each loop for efficiency. */
  479. while (a < lastgroup) {
  480. min0 = (ResultType)(a[0] < b[0] ? a[0] : b[0]);
  481. min1 = (ResultType)(a[1] < b[1] ? a[1] : b[1]);
  482. min2 = (ResultType)(a[2] < b[2] ? a[2] : b[2]);
  483. min3 = (ResultType)(a[3] < b[3] ? a[3] : b[3]);
  484. result += min0 + min1 + min2 + min3;
  485. a += 4;
  486. b += 4;
  487. if ((worst_dist>0)&&(result>worst_dist)) {
  488. return result;
  489. }
  490. }
  491. /* Process last 0-3 pixels. Not needed for standard vector lengths. */
  492. while (a < last) {
  493. min0 = (ResultType)(*a < *b ? *a : *b);
  494. result += min0;
  495. ++a;
  496. ++b;
  497. }
  498. return result;
  499. }
  500. /**
  501. * Partial distance, used by the kd-tree.
  502. */
  503. template <typename U, typename V>
  504. inline ResultType accum_dist(const U& a, const V& b, int) const
  505. {
  506. return a<b ? a : b;
  507. }
  508. };
  509. template<class T>
  510. struct HellingerDistance
  511. {
  512. typedef True is_kdtree_distance;
  513. typedef True is_vector_space_distance;
  514. typedef T ElementType;
  515. typedef typename Accumulator<T>::Type ResultType;
  516. /**
  517. * Compute the Hellinger distance
  518. */
  519. template <typename Iterator1, typename Iterator2>
  520. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
  521. {
  522. ResultType result = ResultType();
  523. ResultType diff0, diff1, diff2, diff3;
  524. Iterator1 last = a + size;
  525. Iterator1 lastgroup = last - 3;
  526. /* Process 4 items with each loop for efficiency. */
  527. while (a < lastgroup) {
  528. diff0 = sqrt(static_cast<ResultType>(a[0])) - sqrt(static_cast<ResultType>(b[0]));
  529. diff1 = sqrt(static_cast<ResultType>(a[1])) - sqrt(static_cast<ResultType>(b[1]));
  530. diff2 = sqrt(static_cast<ResultType>(a[2])) - sqrt(static_cast<ResultType>(b[2]));
  531. diff3 = sqrt(static_cast<ResultType>(a[3])) - sqrt(static_cast<ResultType>(b[3]));
  532. result += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3;
  533. a += 4;
  534. b += 4;
  535. }
  536. while (a < last) {
  537. diff0 = sqrt(static_cast<ResultType>(*a++)) - sqrt(static_cast<ResultType>(*b++));
  538. result += diff0 * diff0;
  539. }
  540. return result;
  541. }
  542. /**
  543. * Partial distance, used by the kd-tree.
  544. */
  545. template <typename U, typename V>
  546. inline ResultType accum_dist(const U& a, const V& b, int) const
  547. {
  548. ResultType diff = sqrt(static_cast<ResultType>(a)) - sqrt(static_cast<ResultType>(b));
  549. return diff * diff;
  550. }
  551. };
  552. template<class T>
  553. struct ChiSquareDistance
  554. {
  555. typedef True is_kdtree_distance;
  556. typedef True is_vector_space_distance;
  557. typedef T ElementType;
  558. typedef typename Accumulator<T>::Type ResultType;
  559. /**
  560. * Compute the chi-square distance
  561. */
  562. template <typename Iterator1, typename Iterator2>
  563. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
  564. {
  565. ResultType result = ResultType();
  566. ResultType sum, diff;
  567. Iterator1 last = a + size;
  568. while (a < last) {
  569. sum = (ResultType)(*a + *b);
  570. if (sum>0) {
  571. diff = (ResultType)(*a - *b);
  572. result += diff*diff/sum;
  573. }
  574. ++a;
  575. ++b;
  576. if ((worst_dist>0)&&(result>worst_dist)) {
  577. return result;
  578. }
  579. }
  580. return result;
  581. }
  582. /**
  583. * Partial distance, used by the kd-tree.
  584. */
  585. template <typename U, typename V>
  586. inline ResultType accum_dist(const U& a, const V& b, int) const
  587. {
  588. ResultType result = ResultType();
  589. ResultType sum, diff;
  590. sum = (ResultType)(a+b);
  591. if (sum>0) {
  592. diff = (ResultType)(a-b);
  593. result = diff*diff/sum;
  594. }
  595. return result;
  596. }
  597. };
  598. template<class T>
  599. struct KL_Divergence
  600. {
  601. typedef True is_kdtree_distance;
  602. typedef True is_vector_space_distance;
  603. typedef T ElementType;
  604. typedef typename Accumulator<T>::Type ResultType;
  605. /**
  606. * Compute the Kullback-Leibler divergence
  607. */
  608. template <typename Iterator1, typename Iterator2>
  609. ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
  610. {
  611. ResultType result = ResultType();
  612. Iterator1 last = a + size;
  613. while (a < last) {
  614. if (* b != 0) {
  615. ResultType ratio = (ResultType)(*a / *b);
  616. if (ratio>0) {
  617. result += *a * log(ratio);
  618. }
  619. }
  620. ++a;
  621. ++b;
  622. if ((worst_dist>0)&&(result>worst_dist)) {
  623. return result;
  624. }
  625. }
  626. return result;
  627. }
  628. /**
  629. * Partial distance, used by the kd-tree.
  630. */
  631. template <typename U, typename V>
  632. inline ResultType accum_dist(const U& a, const V& b, int) const
  633. {
  634. ResultType result = ResultType();
  635. if( *b != 0 ) {
  636. ResultType ratio = (ResultType)(a / b);
  637. if (ratio>0) {
  638. result = a * log(ratio);
  639. }
  640. }
  641. return result;
  642. }
  643. };
  644. /*
  645. * This is a "zero iterator". It basically behaves like a zero filled
  646. * array to all algorithms that use arrays as iterators (STL style).
  647. * It's useful when there's a need to compute the distance between feature
  648. * and origin it and allows for better compiler optimisation than using a
  649. * zero-filled array.
  650. */
  651. template <typename T>
  652. struct ZeroIterator
  653. {
  654. T operator*()
  655. {
  656. return 0;
  657. }
  658. T operator[](int)
  659. {
  660. return 0;
  661. }
  662. const ZeroIterator<T>& operator ++()
  663. {
  664. return *this;
  665. }
  666. ZeroIterator<T> operator ++(int)
  667. {
  668. return *this;
  669. }
  670. ZeroIterator<T>& operator+=(int)
  671. {
  672. return *this;
  673. }
  674. };
  675. /*
  676. * Depending on processed distances, some of them are already squared (e.g. L2)
  677. * and some are not (e.g.Hamming). In KMeans++ for instance we want to be sure
  678. * we are working on ^2 distances, thus following templates to ensure that.
  679. */
  680. template <typename Distance, typename ElementType>
  681. struct squareDistance
  682. {
  683. typedef typename Distance::ResultType ResultType;
  684. ResultType operator()( ResultType dist ) { return dist*dist; }
  685. };
  686. template <typename ElementType>
  687. struct squareDistance<L2_Simple<ElementType>, ElementType>
  688. {
  689. typedef typename L2_Simple<ElementType>::ResultType ResultType;
  690. ResultType operator()( ResultType dist ) { return dist; }
  691. };
  692. template <typename ElementType>
  693. struct squareDistance<L2<ElementType>, ElementType>
  694. {
  695. typedef typename L2<ElementType>::ResultType ResultType;
  696. ResultType operator()( ResultType dist ) { return dist; }
  697. };
  698. template <typename ElementType>
  699. struct squareDistance<MinkowskiDistance<ElementType>, ElementType>
  700. {
  701. typedef typename MinkowskiDistance<ElementType>::ResultType ResultType;
  702. ResultType operator()( ResultType dist ) { return dist; }
  703. };
  704. template <typename ElementType>
  705. struct squareDistance<HellingerDistance<ElementType>, ElementType>
  706. {
  707. typedef typename HellingerDistance<ElementType>::ResultType ResultType;
  708. ResultType operator()( ResultType dist ) { return dist; }
  709. };
  710. template <typename ElementType>
  711. struct squareDistance<ChiSquareDistance<ElementType>, ElementType>
  712. {
  713. typedef typename ChiSquareDistance<ElementType>::ResultType ResultType;
  714. ResultType operator()( ResultType dist ) { return dist; }
  715. };
  716. template <typename Distance>
  717. typename Distance::ResultType ensureSquareDistance( typename Distance::ResultType dist )
  718. {
  719. typedef typename Distance::ElementType ElementType;
  720. squareDistance<Distance, ElementType> dummy;
  721. return dummy( dist );
  722. }
  723. /*
  724. * ...and a template to ensure the user that he will process the normal distance,
  725. * and not squared distance, without losing processing time calling sqrt(ensureSquareDistance)
  726. * that will result in doing actually sqrt(dist*dist) for L1 distance for instance.
  727. */
  728. template <typename Distance, typename ElementType>
  729. struct simpleDistance
  730. {
  731. typedef typename Distance::ResultType ResultType;
  732. ResultType operator()( ResultType dist ) { return dist; }
  733. };
  734. template <typename ElementType>
  735. struct simpleDistance<L2_Simple<ElementType>, ElementType>
  736. {
  737. typedef typename L2_Simple<ElementType>::ResultType ResultType;
  738. ResultType operator()( ResultType dist ) { return sqrt(dist); }
  739. };
  740. template <typename ElementType>
  741. struct simpleDistance<L2<ElementType>, ElementType>
  742. {
  743. typedef typename L2<ElementType>::ResultType ResultType;
  744. ResultType operator()( ResultType dist ) { return sqrt(dist); }
  745. };
  746. template <typename ElementType>
  747. struct simpleDistance<MinkowskiDistance<ElementType>, ElementType>
  748. {
  749. typedef typename MinkowskiDistance<ElementType>::ResultType ResultType;
  750. ResultType operator()( ResultType dist ) { return sqrt(dist); }
  751. };
  752. template <typename ElementType>
  753. struct simpleDistance<HellingerDistance<ElementType>, ElementType>
  754. {
  755. typedef typename HellingerDistance<ElementType>::ResultType ResultType;
  756. ResultType operator()( ResultType dist ) { return sqrt(dist); }
  757. };
  758. template <typename ElementType>
  759. struct simpleDistance<ChiSquareDistance<ElementType>, ElementType>
  760. {
  761. typedef typename ChiSquareDistance<ElementType>::ResultType ResultType;
  762. ResultType operator()( ResultType dist ) { return sqrt(dist); }
  763. };
  764. template <typename Distance>
  765. typename Distance::ResultType ensureSimpleDistance( typename Distance::ResultType dist )
  766. {
  767. typedef typename Distance::ElementType ElementType;
  768. simpleDistance<Distance, ElementType> dummy;
  769. return dummy( dist );
  770. }
  771. }
  772. #endif //OPENCV_FLANN_DIST_H_