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