← Back to branch summary

~paparazzi-uav/paparazzi/v5.0-manual

« back to all changes in this revision

Viewing changes to sw/ext/opencv_bebop/opencv/modules/flann/include/opencv2/flann/dist.h

  • Committer: Paparazzi buildbot
  • Date: 2016-05-18 15:00:29 UTC
  • Revision ID: felix.ruess+docbot@gmail.com-20160518150029-e8lgzi5kvb4p7un9
Manual import commit 4b8bbb730080dac23cf816b98908dacfabe2a8ec from v5.0 branch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
sw/ext/opencv_bebop/opencv/modules/flann/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, size_t 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 (size_t i = 0; i < size; i++) {
 
404
            result += popCountTable[a[i] ^ b[i]];
 
405
        }
 
406
        return result;
 
407
    }
 
408
};
 
409
 
 
410
/**
 
411
 * Hamming distance functor (pop count between two binary vectors, i.e. xor them and count the number of bits set)
 
412
 * That code was taken from brief.cpp in OpenCV
 
413
 */
 
414
template<class T>
 
415
struct Hamming
 
416
{
 
417
    typedef False is_kdtree_distance;
 
418
    typedef False is_vector_space_distance;
 
419
 
 
420
 
 
421
    typedef T ElementType;
 
422
    typedef int ResultType;
 
423
 
 
424
    template<typename Iterator1, typename Iterator2>
 
425
    ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
 
426
    {
 
427
        ResultType result = 0;
 
428
#ifdef __ARM_NEON__
 
429
        {
 
430
            uint32x4_t bits = vmovq_n_u32(0);
 
431
            for (size_t i = 0; i < size; i += 16) {
 
432
                uint8x16_t A_vec = vld1q_u8 (a + i);
 
433
                uint8x16_t B_vec = vld1q_u8 (b + i);
 
434
                uint8x16_t AxorB = veorq_u8 (A_vec, B_vec);
 
435
                uint8x16_t bitsSet = vcntq_u8 (AxorB);
 
436
                uint16x8_t bitSet8 = vpaddlq_u8 (bitsSet);
 
437
                uint32x4_t bitSet4 = vpaddlq_u16 (bitSet8);
 
438
                bits = vaddq_u32(bits, bitSet4);
 
439
            }
 
440
            uint64x2_t bitSet2 = vpaddlq_u32 (bits);
 
441
            result = vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),0);
 
442
            result += vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),2);
 
443
        }
 
444
#elif __GNUC__
 
445
        {
 
446
            //for portability just use unsigned long -- and use the __builtin_popcountll (see docs for __builtin_popcountll)
 
447
            typedef unsigned long long pop_t;
 
448
            const size_t modulo = size % sizeof(pop_t);
 
449
            const pop_t* a2 = reinterpret_cast<const pop_t*> (a);
 
450
            const pop_t* b2 = reinterpret_cast<const pop_t*> (b);
 
451
            const pop_t* a2_end = a2 + (size / sizeof(pop_t));
 
452
 
 
453
            for (; a2 != a2_end; ++a2, ++b2) result += __builtin_popcountll((*a2) ^ (*b2));
 
454
 
 
455
            if (modulo) {
 
456
                //in the case where size is not dividable by sizeof(size_t)
 
457
                //need to mask off the bits at the end
 
458
                pop_t a_final = 0, b_final = 0;
 
459
                memcpy(&a_final, a2, modulo);
 
460
                memcpy(&b_final, b2, modulo);
 
461
                result += __builtin_popcountll(a_final ^ b_final);
 
462
            }
 
463
        }
 
464
#else // NO NEON and NOT GNUC
 
465
        typedef unsigned long long pop_t;
 
466
        HammingLUT lut;
 
467
        result = lut(reinterpret_cast<const unsigned char*> (a),
 
468
                     reinterpret_cast<const unsigned char*> (b), size * sizeof(pop_t));
 
469
#endif
 
470
        return result;
 
471
    }
 
472
};
 
473
 
 
474
template<typename T>
 
475
struct Hamming2
 
476
{
 
477
    typedef False is_kdtree_distance;
 
478
    typedef False is_vector_space_distance;
 
479
 
 
480
    typedef T ElementType;
 
481
    typedef int ResultType;
 
482
 
 
483
    /** This is popcount_3() from:
 
484
     * http://en.wikipedia.org/wiki/Hamming_weight */
 
485
    unsigned int popcnt32(uint32_t n) const
 
486
    {
 
487
        n -= ((n >> 1) & 0x55555555);
 
488
        n = (n & 0x33333333) + ((n >> 2) & 0x33333333);
 
489
        return (((n + (n >> 4))& 0xF0F0F0F)* 0x1010101) >> 24;
 
490
    }
 
491
 
 
492
#ifdef FLANN_PLATFORM_64_BIT
 
493
    unsigned int popcnt64(uint64_t n) const
 
494
    {
 
495
        n -= ((n >> 1) & 0x5555555555555555);
 
496
        n = (n & 0x3333333333333333) + ((n >> 2) & 0x3333333333333333);
 
497
        return (((n + (n >> 4))& 0x0f0f0f0f0f0f0f0f)* 0x0101010101010101) >> 56;
 
498
    }
 
499
#endif
 
500
 
 
501
    template <typename Iterator1, typename Iterator2>
 
502
    ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
 
503
    {
 
504
#ifdef FLANN_PLATFORM_64_BIT
 
505
        const uint64_t* pa = reinterpret_cast<const uint64_t*>(a);
 
506
        const uint64_t* pb = reinterpret_cast<const uint64_t*>(b);
 
507
        ResultType result = 0;
 
508
        size /= (sizeof(uint64_t)/sizeof(unsigned char));
 
509
        for(size_t i = 0; i < size; ++i ) {
 
510
            result += popcnt64(*pa ^ *pb);
 
511
            ++pa;
 
512
            ++pb;
 
513
        }
 
514
#else
 
515
        const uint32_t* pa = reinterpret_cast<const uint32_t*>(a);
 
516
        const uint32_t* pb = reinterpret_cast<const uint32_t*>(b);
 
517
        ResultType result = 0;
 
518
        size /= (sizeof(uint32_t)/sizeof(unsigned char));
 
519
        for(size_t i = 0; i < size; ++i ) {
 
520
            result += popcnt32(*pa ^ *pb);
 
521
            ++pa;
 
522
            ++pb;
 
523
        }
 
524
#endif
 
525
        return result;
 
526
    }
 
527
};
 
528
 
 
529
 
 
530
 
 
531
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
532
 
 
533
template<class T>
 
534
struct HistIntersectionDistance
 
535
{
 
536
    typedef True is_kdtree_distance;
 
537
    typedef True is_vector_space_distance;
 
538
 
 
539
    typedef T ElementType;
 
540
    typedef typename Accumulator<T>::Type ResultType;
 
541
 
 
542
    /**
 
543
     *  Compute the histogram intersection distance
 
544
     */
 
545
    template <typename Iterator1, typename Iterator2>
 
546
    ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
 
547
    {
 
548
        ResultType result = ResultType();
 
549
        ResultType min0, min1, min2, min3;
 
550
        Iterator1 last = a + size;
 
551
        Iterator1 lastgroup = last - 3;
 
552
 
 
553
        /* Process 4 items with each loop for efficiency. */
 
554
        while (a < lastgroup) {
 
555
            min0 = (ResultType)(a[0] < b[0] ? a[0] : b[0]);
 
556
            min1 = (ResultType)(a[1] < b[1] ? a[1] : b[1]);
 
557
            min2 = (ResultType)(a[2] < b[2] ? a[2] : b[2]);
 
558
            min3 = (ResultType)(a[3] < b[3] ? a[3] : b[3]);
 
559
            result += min0 + min1 + min2 + min3;
 
560
            a += 4;
 
561
            b += 4;
 
562
            if ((worst_dist>0)&&(result>worst_dist)) {
 
563
                return result;
 
564
            }
 
565
        }
 
566
        /* Process last 0-3 pixels.  Not needed for standard vector lengths. */
 
567
        while (a < last) {
 
568
            min0 = (ResultType)(*a < *b ? *a : *b);
 
569
            result += min0;
 
570
            ++a;
 
571
            ++b;
 
572
        }
 
573
        return result;
 
574
    }
 
575
 
 
576
    /**
 
577
     * Partial distance, used by the kd-tree.
 
578
     */
 
579
    template <typename U, typename V>
 
580
    inline ResultType accum_dist(const U& a, const V& b, int) const
 
581
    {
 
582
        return a<b ? a : b;
 
583
    }
 
584
};
 
585
 
 
586
 
 
587
 
 
588
template<class T>
 
589
struct HellingerDistance
 
590
{
 
591
    typedef True is_kdtree_distance;
 
592
    typedef True is_vector_space_distance;
 
593
 
 
594
    typedef T ElementType;
 
595
    typedef typename Accumulator<T>::Type ResultType;
 
596
 
 
597
    /**
 
598
     *  Compute the Hellinger distance
 
599
     */
 
600
    template <typename Iterator1, typename Iterator2>
 
601
    ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
 
602
    {
 
603
        ResultType result = ResultType();
 
604
        ResultType diff0, diff1, diff2, diff3;
 
605
        Iterator1 last = a + size;
 
606
        Iterator1 lastgroup = last - 3;
 
607
 
 
608
        /* Process 4 items with each loop for efficiency. */
 
609
        while (a < lastgroup) {
 
610
            diff0 = sqrt(static_cast<ResultType>(a[0])) - sqrt(static_cast<ResultType>(b[0]));
 
611
            diff1 = sqrt(static_cast<ResultType>(a[1])) - sqrt(static_cast<ResultType>(b[1]));
 
612
            diff2 = sqrt(static_cast<ResultType>(a[2])) - sqrt(static_cast<ResultType>(b[2]));
 
613
            diff3 = sqrt(static_cast<ResultType>(a[3])) - sqrt(static_cast<ResultType>(b[3]));
 
614
            result += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3;
 
615
            a += 4;
 
616
            b += 4;
 
617
        }
 
618
        while (a < last) {
 
619
            diff0 = sqrt(static_cast<ResultType>(*a++)) - sqrt(static_cast<ResultType>(*b++));
 
620
            result += diff0 * diff0;
 
621
        }
 
622
        return result;
 
623
    }
 
624
 
 
625
    /**
 
626
     * Partial distance, used by the kd-tree.
 
627
     */
 
628
    template <typename U, typename V>
 
629
    inline ResultType accum_dist(const U& a, const V& b, int) const
 
630
    {
 
631
        ResultType diff = sqrt(static_cast<ResultType>(a)) - sqrt(static_cast<ResultType>(b));
 
632
        return diff * diff;
 
633
    }
 
634
};
 
635
 
 
636
 
 
637
template<class T>
 
638
struct ChiSquareDistance
 
639
{
 
640
    typedef True is_kdtree_distance;
 
641
    typedef True is_vector_space_distance;
 
642
 
 
643
    typedef T ElementType;
 
644
    typedef typename Accumulator<T>::Type ResultType;
 
645
 
 
646
    /**
 
647
     *  Compute the chi-square distance
 
648
     */
 
649
    template <typename Iterator1, typename Iterator2>
 
650
    ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
 
651
    {
 
652
        ResultType result = ResultType();
 
653
        ResultType sum, diff;
 
654
        Iterator1 last = a + size;
 
655
 
 
656
        while (a < last) {
 
657
            sum = (ResultType)(*a + *b);
 
658
            if (sum>0) {
 
659
                diff = (ResultType)(*a - *b);
 
660
                result += diff*diff/sum;
 
661
            }
 
662
            ++a;
 
663
            ++b;
 
664
 
 
665
            if ((worst_dist>0)&&(result>worst_dist)) {
 
666
                return result;
 
667
            }
 
668
        }
 
669
        return result;
 
670
    }
 
671
 
 
672
    /**
 
673
     * Partial distance, used by the kd-tree.
 
674
     */
 
675
    template <typename U, typename V>
 
676
    inline ResultType accum_dist(const U& a, const V& b, int) const
 
677
    {
 
678
        ResultType result = ResultType();
 
679
        ResultType sum, diff;
 
680
 
 
681
        sum = (ResultType)(a+b);
 
682
        if (sum>0) {
 
683
            diff = (ResultType)(a-b);
 
684
            result = diff*diff/sum;
 
685
        }
 
686
        return result;
 
687
    }
 
688
};
 
689
 
 
690
 
 
691
template<class T>
 
692
struct KL_Divergence
 
693
{
 
694
    typedef True is_kdtree_distance;
 
695
    typedef True is_vector_space_distance;
 
696
 
 
697
    typedef T ElementType;
 
698
    typedef typename Accumulator<T>::Type ResultType;
 
699
 
 
700
    /**
 
701
     *  Compute the Kullback–Leibler divergence
 
702
     */
 
703
    template <typename Iterator1, typename Iterator2>
 
704
    ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType worst_dist = -1) const
 
705
    {
 
706
        ResultType result = ResultType();
 
707
        Iterator1 last = a + size;
 
708
 
 
709
        while (a < last) {
 
710
            if (* b != 0) {
 
711
                ResultType ratio = (ResultType)(*a / *b);
 
712
                if (ratio>0) {
 
713
                    result += *a * log(ratio);
 
714
                }
 
715
            }
 
716
            ++a;
 
717
            ++b;
 
718
 
 
719
            if ((worst_dist>0)&&(result>worst_dist)) {
 
720
                return result;
 
721
            }
 
722
        }
 
723
        return result;
 
724
    }
 
725
 
 
726
    /**
 
727
     * Partial distance, used by the kd-tree.
 
728
     */
 
729
    template <typename U, typename V>
 
730
    inline ResultType accum_dist(const U& a, const V& b, int) const
 
731
    {
 
732
        ResultType result = ResultType();
 
733
        if( *b != 0 ) {
 
734
            ResultType ratio = (ResultType)(a / b);
 
735
            if (ratio>0) {
 
736
                result = a * log(ratio);
 
737
            }
 
738
        }
 
739
        return result;
 
740
    }
 
741
};
 
742
 
 
743
 
 
744
 
 
745
/*
 
746
 * This is a "zero iterator". It basically behaves like a zero filled
 
747
 * array to all algorithms that use arrays as iterators (STL style).
 
748
 * It's useful when there's a need to compute the distance between feature
 
749
 * and origin it and allows for better compiler optimisation than using a
 
750
 * zero-filled array.
 
751
 */
 
752
template <typename T>
 
753
struct ZeroIterator
 
754
{
 
755
 
 
756
    T operator*()
 
757
    {
 
758
        return 0;
 
759
    }
 
760
 
 
761
    T operator[](int)
 
762
    {
 
763
        return 0;
 
764
    }
 
765
 
 
766
    const ZeroIterator<T>& operator ++()
 
767
    {
 
768
        return *this;
 
769
    }
 
770
 
 
771
    ZeroIterator<T> operator ++(int)
 
772
    {
 
773
        return *this;
 
774
    }
 
775
 
 
776
    ZeroIterator<T>& operator+=(int)
 
777
    {
 
778
        return *this;
 
779
    }
 
780
 
 
781
};
 
782
 
 
783
 
 
784
/*
 
785
 * Depending on processed distances, some of them are already squared (e.g. L2)
 
786
 * and some are not (e.g.Hamming). In KMeans++ for instance we want to be sure
 
787
 * we are working on ^2 distances, thus following templates to ensure that.
 
788
 */
 
789
template <typename Distance, typename ElementType>
 
790
struct squareDistance
 
791
{
 
792
    typedef typename Distance::ResultType ResultType;
 
793
    ResultType operator()( ResultType dist ) { return dist*dist; }
 
794
};
 
795
 
 
796
 
 
797
template <typename ElementType>
 
798
struct squareDistance<L2_Simple<ElementType>, ElementType>
 
799
{
 
800
    typedef typename L2_Simple<ElementType>::ResultType ResultType;
 
801
    ResultType operator()( ResultType dist ) { return dist; }
 
802
};
 
803
 
 
804
template <typename ElementType>
 
805
struct squareDistance<L2<ElementType>, ElementType>
 
806
{
 
807
    typedef typename L2<ElementType>::ResultType ResultType;
 
808
    ResultType operator()( ResultType dist ) { return dist; }
 
809
};
 
810
 
 
811
 
 
812
template <typename ElementType>
 
813
struct squareDistance<MinkowskiDistance<ElementType>, ElementType>
 
814
{
 
815
    typedef typename MinkowskiDistance<ElementType>::ResultType ResultType;
 
816
    ResultType operator()( ResultType dist ) { return dist; }
 
817
};
 
818
 
 
819
template <typename ElementType>
 
820
struct squareDistance<HellingerDistance<ElementType>, ElementType>
 
821
{
 
822
    typedef typename HellingerDistance<ElementType>::ResultType ResultType;
 
823
    ResultType operator()( ResultType dist ) { return dist; }
 
824
};
 
825
 
 
826
template <typename ElementType>
 
827
struct squareDistance<ChiSquareDistance<ElementType>, ElementType>
 
828
{
 
829
    typedef typename ChiSquareDistance<ElementType>::ResultType ResultType;
 
830
    ResultType operator()( ResultType dist ) { return dist; }
 
831
};
 
832
 
 
833
 
 
834
template <typename Distance>
 
835
typename Distance::ResultType ensureSquareDistance( typename Distance::ResultType dist )
 
836
{
 
837
    typedef typename Distance::ElementType ElementType;
 
838
 
 
839
    squareDistance<Distance, ElementType> dummy;
 
840
    return dummy( dist );
 
841
}
 
842
 
 
843
 
 
844
/*
 
845
 * ...and a template to ensure the user that he will process the normal distance,
 
846
 * and not squared distance, without loosing processing time calling sqrt(ensureSquareDistance)
 
847
 * that will result in doing actually sqrt(dist*dist) for L1 distance for instance.
 
848
 */
 
849
template <typename Distance, typename ElementType>
 
850
struct simpleDistance
 
851
{
 
852
    typedef typename Distance::ResultType ResultType;
 
853
    ResultType operator()( ResultType dist ) { return dist; }
 
854
};
 
855
 
 
856
 
 
857
template <typename ElementType>
 
858
struct simpleDistance<L2_Simple<ElementType>, ElementType>
 
859
{
 
860
    typedef typename L2_Simple<ElementType>::ResultType ResultType;
 
861
    ResultType operator()( ResultType dist ) { return sqrt(dist); }
 
862
};
 
863
 
 
864
template <typename ElementType>
 
865
struct simpleDistance<L2<ElementType>, ElementType>
 
866
{
 
867
    typedef typename L2<ElementType>::ResultType ResultType;
 
868
    ResultType operator()( ResultType dist ) { return sqrt(dist); }
 
869
};
 
870
 
 
871
 
 
872
template <typename ElementType>
 
873
struct simpleDistance<MinkowskiDistance<ElementType>, ElementType>
 
874
{
 
875
    typedef typename MinkowskiDistance<ElementType>::ResultType ResultType;
 
876
    ResultType operator()( ResultType dist ) { return sqrt(dist); }
 
877
};
 
878
 
 
879
template <typename ElementType>
 
880
struct simpleDistance<HellingerDistance<ElementType>, ElementType>
 
881
{
 
882
    typedef typename HellingerDistance<ElementType>::ResultType ResultType;
 
883
    ResultType operator()( ResultType dist ) { return sqrt(dist); }
 
884
};
 
885
 
 
886
template <typename ElementType>
 
887
struct simpleDistance<ChiSquareDistance<ElementType>, ElementType>
 
888
{
 
889
    typedef typename ChiSquareDistance<ElementType>::ResultType ResultType;
 
890
    ResultType operator()( ResultType dist ) { return sqrt(dist); }
 
891
};
 
892
 
 
893
 
 
894
template <typename Distance>
 
895
typename Distance::ResultType ensureSimpleDistance( typename Distance::ResultType dist )
 
896
{
 
897
    typedef typename Distance::ElementType ElementType;
 
898
 
 
899
    simpleDistance<Distance, ElementType> dummy;
 
900
    return dummy( dist );
 
901
}
 
902
 
 
903
}
 
904
 
 
905
#endif //OPENCV_FLANN_DIST_H_