← 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/hierarchical_clustering_index.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/hierarchical_clustering_index.h
 
1
/***********************************************************************
 
2
 * Software License Agreement (BSD License)
 
3
 *
 
4
 * Copyright 2008-2011  Marius Muja (mariusm@cs.ubc.ca). All rights reserved.
 
5
 * Copyright 2008-2011  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_HIERARCHICAL_CLUSTERING_INDEX_H_
 
32
#define OPENCV_FLANN_HIERARCHICAL_CLUSTERING_INDEX_H_
 
33
 
 
34
#include <algorithm>
 
35
#include <map>
 
36
#include <cassert>
 
37
#include <limits>
 
38
#include <cmath>
 
39
 
 
40
#include "general.h"
 
41
#include "nn_index.h"
 
42
#include "dist.h"
 
43
#include "matrix.h"
 
44
#include "result_set.h"
 
45
#include "heap.h"
 
46
#include "allocator.h"
 
47
#include "random.h"
 
48
#include "saving.h"
 
49
 
 
50
 
 
51
namespace cvflann
 
52
{
 
53
 
 
54
struct HierarchicalClusteringIndexParams : public IndexParams
 
55
{
 
56
    HierarchicalClusteringIndexParams(int branching = 32,
 
57
                                      flann_centers_init_t centers_init = FLANN_CENTERS_RANDOM,
 
58
                                      int trees = 4, int leaf_size = 100)
 
59
    {
 
60
        (*this)["algorithm"] = FLANN_INDEX_HIERARCHICAL;
 
61
        // The branching factor used in the hierarchical clustering
 
62
        (*this)["branching"] = branching;
 
63
        // Algorithm used for picking the initial cluster centers
 
64
        (*this)["centers_init"] = centers_init;
 
65
        // number of parallel trees to build
 
66
        (*this)["trees"] = trees;
 
67
        // maximum leaf size
 
68
        (*this)["leaf_size"] = leaf_size;
 
69
    }
 
70
};
 
71
 
 
72
 
 
73
/**
 
74
 * Hierarchical index
 
75
 *
 
76
 * Contains a tree constructed through a hierarchical clustering
 
77
 * and other information for indexing a set of points for nearest-neighbour matching.
 
78
 */
 
79
template <typename Distance>
 
80
class HierarchicalClusteringIndex : public NNIndex<Distance>
 
81
{
 
82
public:
 
83
    typedef typename Distance::ElementType ElementType;
 
84
    typedef typename Distance::ResultType DistanceType;
 
85
 
 
86
private:
 
87
 
 
88
 
 
89
    typedef void (HierarchicalClusteringIndex::* centersAlgFunction)(int, int*, int, int*, int&);
 
90
 
 
91
    /**
 
92
     * The function used for choosing the cluster centers.
 
93
     */
 
94
    centersAlgFunction chooseCenters;
 
95
 
 
96
 
 
97
 
 
98
    /**
 
99
     * Chooses the initial centers in the k-means clustering in a random manner.
 
100
     *
 
101
     * Params:
 
102
     *     k = number of centers
 
103
     *     vecs = the dataset of points
 
104
     *     indices = indices in the dataset
 
105
     *     indices_length = length of indices vector
 
106
     *
 
107
     */
 
108
    void chooseCentersRandom(int k, int* dsindices, int indices_length, int* centers, int& centers_length)
 
109
    {
 
110
        UniqueRandom r(indices_length);
 
111
 
 
112
        int index;
 
113
        for (index=0; index<k; ++index) {
 
114
            bool duplicate = true;
 
115
            int rnd;
 
116
            while (duplicate) {
 
117
                duplicate = false;
 
118
                rnd = r.next();
 
119
                if (rnd<0) {
 
120
                    centers_length = index;
 
121
                    return;
 
122
                }
 
123
 
 
124
                centers[index] = dsindices[rnd];
 
125
 
 
126
                for (int j=0; j<index; ++j) {
 
127
                    DistanceType sq = distance(dataset[centers[index]], dataset[centers[j]], dataset.cols);
 
128
                    if (sq<1e-16) {
 
129
                        duplicate = true;
 
130
                    }
 
131
                }
 
132
            }
 
133
        }
 
134
 
 
135
        centers_length = index;
 
136
    }
 
137
 
 
138
 
 
139
    /**
 
140
     * Chooses the initial centers in the k-means using Gonzales' algorithm
 
141
     * so that the centers are spaced apart from each other.
 
142
     *
 
143
     * Params:
 
144
     *     k = number of centers
 
145
     *     vecs = the dataset of points
 
146
     *     indices = indices in the dataset
 
147
     * Returns:
 
148
     */
 
149
    void chooseCentersGonzales(int k, int* dsindices, int indices_length, int* centers, int& centers_length)
 
150
    {
 
151
        int n = indices_length;
 
152
 
 
153
        int rnd = rand_int(n);
 
154
        assert(rnd >=0 && rnd < n);
 
155
 
 
156
        centers[0] = dsindices[rnd];
 
157
 
 
158
        int index;
 
159
        for (index=1; index<k; ++index) {
 
160
 
 
161
            int best_index = -1;
 
162
            DistanceType best_val = 0;
 
163
            for (int j=0; j<n; ++j) {
 
164
                DistanceType dist = distance(dataset[centers[0]],dataset[dsindices[j]],dataset.cols);
 
165
                for (int i=1; i<index; ++i) {
 
166
                    DistanceType tmp_dist = distance(dataset[centers[i]],dataset[dsindices[j]],dataset.cols);
 
167
                    if (tmp_dist<dist) {
 
168
                        dist = tmp_dist;
 
169
                    }
 
170
                }
 
171
                if (dist>best_val) {
 
172
                    best_val = dist;
 
173
                    best_index = j;
 
174
                }
 
175
            }
 
176
            if (best_index!=-1) {
 
177
                centers[index] = dsindices[best_index];
 
178
            }
 
179
            else {
 
180
                break;
 
181
            }
 
182
        }
 
183
        centers_length = index;
 
184
    }
 
185
 
 
186
 
 
187
    /**
 
188
     * Chooses the initial centers in the k-means using the algorithm
 
189
     * proposed in the KMeans++ paper:
 
190
     * Arthur, David; Vassilvitskii, Sergei - k-means++: The Advantages of Careful Seeding
 
191
     *
 
192
     * Implementation of this function was converted from the one provided in Arthur's code.
 
193
     *
 
194
     * Params:
 
195
     *     k = number of centers
 
196
     *     vecs = the dataset of points
 
197
     *     indices = indices in the dataset
 
198
     * Returns:
 
199
     */
 
200
    void chooseCentersKMeanspp(int k, int* dsindices, int indices_length, int* centers, int& centers_length)
 
201
    {
 
202
        int n = indices_length;
 
203
 
 
204
        double currentPot = 0;
 
205
        DistanceType* closestDistSq = new DistanceType[n];
 
206
 
 
207
        // Choose one random center and set the closestDistSq values
 
208
        int index = rand_int(n);
 
209
        assert(index >=0 && index < n);
 
210
        centers[0] = dsindices[index];
 
211
 
 
212
        // Computing distance^2 will have the advantage of even higher probability further to pick new centers
 
213
        // far from previous centers (and this complies to "k-means++: the advantages of careful seeding" article)
 
214
        for (int i = 0; i < n; i++) {
 
215
            closestDistSq[i] = distance(dataset[dsindices[i]], dataset[dsindices[index]], dataset.cols);
 
216
            closestDistSq[i] = ensureSquareDistance<Distance>( closestDistSq[i] );
 
217
            currentPot += closestDistSq[i];
 
218
        }
 
219
 
 
220
 
 
221
        const int numLocalTries = 1;
 
222
 
 
223
        // Choose each center
 
224
        int centerCount;
 
225
        for (centerCount = 1; centerCount < k; centerCount++) {
 
226
 
 
227
            // Repeat several trials
 
228
            double bestNewPot = -1;
 
229
            int bestNewIndex = 0;
 
230
            for (int localTrial = 0; localTrial < numLocalTries; localTrial++) {
 
231
 
 
232
                // Choose our center - have to be slightly careful to return a valid answer even accounting
 
233
                // for possible rounding errors
 
234
                double randVal = rand_double(currentPot);
 
235
                for (index = 0; index < n-1; index++) {
 
236
                    if (randVal <= closestDistSq[index]) break;
 
237
                    else randVal -= closestDistSq[index];
 
238
                }
 
239
 
 
240
                // Compute the new potential
 
241
                double newPot = 0;
 
242
                for (int i = 0; i < n; i++) {
 
243
                    DistanceType dist = distance(dataset[dsindices[i]], dataset[dsindices[index]], dataset.cols);
 
244
                    newPot += std::min( ensureSquareDistance<Distance>(dist), closestDistSq[i] );
 
245
                }
 
246
 
 
247
                // Store the best result
 
248
                if ((bestNewPot < 0)||(newPot < bestNewPot)) {
 
249
                    bestNewPot = newPot;
 
250
                    bestNewIndex = index;
 
251
                }
 
252
            }
 
253
 
 
254
            // Add the appropriate center
 
255
            centers[centerCount] = dsindices[bestNewIndex];
 
256
            currentPot = bestNewPot;
 
257
            for (int i = 0; i < n; i++) {
 
258
                DistanceType dist = distance(dataset[dsindices[i]], dataset[dsindices[bestNewIndex]], dataset.cols);
 
259
                closestDistSq[i] = std::min( ensureSquareDistance<Distance>(dist), closestDistSq[i] );
 
260
            }
 
261
        }
 
262
 
 
263
        centers_length = centerCount;
 
264
 
 
265
        delete[] closestDistSq;
 
266
    }
 
267
 
 
268
 
 
269
    /**
 
270
     * Chooses the initial centers in a way inspired by Gonzales (by Pierre-Emmanuel Viel):
 
271
     * select the first point of the list as a candidate, then parse the points list. If another
 
272
     * point is further than current candidate from the other centers, test if it is a good center
 
273
     * of a local aggregation. If it is, replace current candidate by this point. And so on...
 
274
     *
 
275
     * Used with KMeansIndex that computes centers coordinates by averaging positions of clusters points,
 
276
     * this doesn't make a real difference with previous methods. But used with HierarchicalClusteringIndex
 
277
     * class that pick centers among existing points instead of computing the barycenters, there is a real
 
278
     * improvement.
 
279
     *
 
280
     * Params:
 
281
     *     k = number of centers
 
282
     *     vecs = the dataset of points
 
283
     *     indices = indices in the dataset
 
284
     * Returns:
 
285
     */
 
286
    void GroupWiseCenterChooser(int k, int* dsindices, int indices_length, int* centers, int& centers_length)
 
287
    {
 
288
        const float kSpeedUpFactor = 1.3f;
 
289
 
 
290
        int n = indices_length;
 
291
 
 
292
        DistanceType* closestDistSq = new DistanceType[n];
 
293
 
 
294
        // Choose one random center and set the closestDistSq values
 
295
        int index = rand_int(n);
 
296
        assert(index >=0 && index < n);
 
297
        centers[0] = dsindices[index];
 
298
 
 
299
        for (int i = 0; i < n; i++) {
 
300
            closestDistSq[i] = distance(dataset[dsindices[i]], dataset[dsindices[index]], dataset.cols);
 
301
        }
 
302
 
 
303
 
 
304
        // Choose each center
 
305
        int centerCount;
 
306
        for (centerCount = 1; centerCount < k; centerCount++) {
 
307
 
 
308
            // Repeat several trials
 
309
            double bestNewPot = -1;
 
310
            int bestNewIndex = 0;
 
311
            DistanceType furthest = 0;
 
312
            for (index = 0; index < n; index++) {
 
313
 
 
314
                // We will test only the potential of the points further than current candidate
 
315
                if( closestDistSq[index] > kSpeedUpFactor * (float)furthest ) {
 
316
 
 
317
                    // Compute the new potential
 
318
                    double newPot = 0;
 
319
                    for (int i = 0; i < n; i++) {
 
320
                        newPot += std::min( distance(dataset[dsindices[i]], dataset[dsindices[index]], dataset.cols)
 
321
                                            , closestDistSq[i] );
 
322
                    }
 
323
 
 
324
                    // Store the best result
 
325
                    if ((bestNewPot < 0)||(newPot <= bestNewPot)) {
 
326
                        bestNewPot = newPot;
 
327
                        bestNewIndex = index;
 
328
                        furthest = closestDistSq[index];
 
329
                    }
 
330
                }
 
331
            }
 
332
 
 
333
            // Add the appropriate center
 
334
            centers[centerCount] = dsindices[bestNewIndex];
 
335
            for (int i = 0; i < n; i++) {
 
336
                closestDistSq[i] = std::min( distance(dataset[dsindices[i]], dataset[dsindices[bestNewIndex]], dataset.cols)
 
337
                                             , closestDistSq[i] );
 
338
            }
 
339
        }
 
340
 
 
341
        centers_length = centerCount;
 
342
 
 
343
        delete[] closestDistSq;
 
344
    }
 
345
 
 
346
 
 
347
public:
 
348
 
 
349
 
 
350
    /**
 
351
     * Index constructor
 
352
     *
 
353
     * Params:
 
354
     *          inputData = dataset with the input features
 
355
     *          params = parameters passed to the hierarchical k-means algorithm
 
356
     */
 
357
    HierarchicalClusteringIndex(const Matrix<ElementType>& inputData, const IndexParams& index_params = HierarchicalClusteringIndexParams(),
 
358
                                Distance d = Distance())
 
359
        : dataset(inputData), params(index_params), root(NULL), indices(NULL), distance(d)
 
360
    {
 
361
        memoryCounter = 0;
 
362
 
 
363
        size_ = dataset.rows;
 
364
        veclen_ = dataset.cols;
 
365
 
 
366
        branching_ = get_param(params,"branching",32);
 
367
        centers_init_ = get_param(params,"centers_init", FLANN_CENTERS_RANDOM);
 
368
        trees_ = get_param(params,"trees",4);
 
369
        leaf_size_ = get_param(params,"leaf_size",100);
 
370
 
 
371
        if (centers_init_==FLANN_CENTERS_RANDOM) {
 
372
            chooseCenters = &HierarchicalClusteringIndex::chooseCentersRandom;
 
373
        }
 
374
        else if (centers_init_==FLANN_CENTERS_GONZALES) {
 
375
            chooseCenters = &HierarchicalClusteringIndex::chooseCentersGonzales;
 
376
        }
 
377
        else if (centers_init_==FLANN_CENTERS_KMEANSPP) {
 
378
            chooseCenters = &HierarchicalClusteringIndex::chooseCentersKMeanspp;
 
379
        }
 
380
        else if (centers_init_==FLANN_CENTERS_GROUPWISE) {
 
381
            chooseCenters = &HierarchicalClusteringIndex::GroupWiseCenterChooser;
 
382
        }
 
383
        else {
 
384
            throw FLANNException("Unknown algorithm for choosing initial centers.");
 
385
        }
 
386
 
 
387
        trees_ = get_param(params,"trees",4);
 
388
        root = new NodePtr[trees_];
 
389
        indices = new int*[trees_];
 
390
 
 
391
        for (int i=0; i<trees_; ++i) {
 
392
            root[i] = NULL;
 
393
            indices[i] = NULL;
 
394
        }
 
395
    }
 
396
 
 
397
    HierarchicalClusteringIndex(const HierarchicalClusteringIndex&);
 
398
    HierarchicalClusteringIndex& operator=(const HierarchicalClusteringIndex&);
 
399
 
 
400
    /**
 
401
     * Index destructor.
 
402
     *
 
403
     * Release the memory used by the index.
 
404
     */
 
405
    virtual ~HierarchicalClusteringIndex()
 
406
    {
 
407
        free_elements();
 
408
 
 
409
        if (root!=NULL) {
 
410
            delete[] root;
 
411
        }
 
412
 
 
413
        if (indices!=NULL) {
 
414
            delete[] indices;
 
415
        }
 
416
    }
 
417
 
 
418
 
 
419
    /**
 
420
     * Release the inner elements of indices[]
 
421
     */
 
422
    void free_elements()
 
423
    {
 
424
        if (indices!=NULL) {
 
425
            for(int i=0; i<trees_; ++i) {
 
426
                if (indices[i]!=NULL) {
 
427
                    delete[] indices[i];
 
428
                    indices[i] = NULL;
 
429
                }
 
430
            }
 
431
        }
 
432
    }
 
433
 
 
434
 
 
435
    /**
 
436
     *  Returns size of index.
 
437
     */
 
438
    size_t size() const
 
439
    {
 
440
        return size_;
 
441
    }
 
442
 
 
443
    /**
 
444
     * Returns the length of an index feature.
 
445
     */
 
446
    size_t veclen() const
 
447
    {
 
448
        return veclen_;
 
449
    }
 
450
 
 
451
 
 
452
    /**
 
453
     * Computes the inde memory usage
 
454
     * Returns: memory used by the index
 
455
     */
 
456
    int usedMemory() const
 
457
    {
 
458
        return pool.usedMemory+pool.wastedMemory+memoryCounter;
 
459
    }
 
460
 
 
461
    /**
 
462
     * Builds the index
 
463
     */
 
464
    void buildIndex()
 
465
    {
 
466
        if (branching_<2) {
 
467
            throw FLANNException("Branching factor must be at least 2");
 
468
        }
 
469
 
 
470
        free_elements();
 
471
 
 
472
        for (int i=0; i<trees_; ++i) {
 
473
            indices[i] = new int[size_];
 
474
            for (size_t j=0; j<size_; ++j) {
 
475
                indices[i][j] = (int)j;
 
476
            }
 
477
            root[i] = pool.allocate<Node>();
 
478
            computeClustering(root[i], indices[i], (int)size_, branching_,0);
 
479
        }
 
480
    }
 
481
 
 
482
 
 
483
    flann_algorithm_t getType() const
 
484
    {
 
485
        return FLANN_INDEX_HIERARCHICAL;
 
486
    }
 
487
 
 
488
 
 
489
    void saveIndex(FILE* stream)
 
490
    {
 
491
        save_value(stream, branching_);
 
492
        save_value(stream, trees_);
 
493
        save_value(stream, centers_init_);
 
494
        save_value(stream, leaf_size_);
 
495
        save_value(stream, memoryCounter);
 
496
        for (int i=0; i<trees_; ++i) {
 
497
            save_value(stream, *indices[i], size_);
 
498
            save_tree(stream, root[i], i);
 
499
        }
 
500
 
 
501
    }
 
502
 
 
503
 
 
504
    void loadIndex(FILE* stream)
 
505
    {
 
506
        free_elements();
 
507
 
 
508
        if (root!=NULL) {
 
509
            delete[] root;
 
510
        }
 
511
 
 
512
        if (indices!=NULL) {
 
513
            delete[] indices;
 
514
        }
 
515
 
 
516
        load_value(stream, branching_);
 
517
        load_value(stream, trees_);
 
518
        load_value(stream, centers_init_);
 
519
        load_value(stream, leaf_size_);
 
520
        load_value(stream, memoryCounter);
 
521
 
 
522
        indices = new int*[trees_];
 
523
        root = new NodePtr[trees_];
 
524
        for (int i=0; i<trees_; ++i) {
 
525
            indices[i] = new int[size_];
 
526
            load_value(stream, *indices[i], size_);
 
527
            load_tree(stream, root[i], i);
 
528
        }
 
529
 
 
530
        params["algorithm"] = getType();
 
531
        params["branching"] = branching_;
 
532
        params["trees"] = trees_;
 
533
        params["centers_init"] = centers_init_;
 
534
        params["leaf_size"] = leaf_size_;
 
535
    }
 
536
 
 
537
 
 
538
    /**
 
539
     * Find set of nearest neighbors to vec. Their indices are stored inside
 
540
     * the result object.
 
541
     *
 
542
     * Params:
 
543
     *     result = the result object in which the indices of the nearest-neighbors are stored
 
544
     *     vec = the vector for which to search the nearest neighbors
 
545
     *     searchParams = parameters that influence the search algorithm (checks)
 
546
     */
 
547
    void findNeighbors(ResultSet<DistanceType>& result, const ElementType* vec, const SearchParams& searchParams)
 
548
    {
 
549
 
 
550
        int maxChecks = get_param(searchParams,"checks",32);
 
551
 
 
552
        // Priority queue storing intermediate branches in the best-bin-first search
 
553
        Heap<BranchSt>* heap = new Heap<BranchSt>((int)size_);
 
554
 
 
555
        std::vector<bool> checked(size_,false);
 
556
        int checks = 0;
 
557
        for (int i=0; i<trees_; ++i) {
 
558
            findNN(root[i], result, vec, checks, maxChecks, heap, checked);
 
559
        }
 
560
 
 
561
        BranchSt branch;
 
562
        while (heap->popMin(branch) && (checks<maxChecks || !result.full())) {
 
563
            NodePtr node = branch.node;
 
564
            findNN(node, result, vec, checks, maxChecks, heap, checked);
 
565
        }
 
566
        assert(result.full());
 
567
 
 
568
        delete heap;
 
569
 
 
570
    }
 
571
 
 
572
    IndexParams getParameters() const
 
573
    {
 
574
        return params;
 
575
    }
 
576
 
 
577
 
 
578
private:
 
579
 
 
580
    /**
 
581
     * Struture representing a node in the hierarchical k-means tree.
 
582
     */
 
583
    struct Node
 
584
    {
 
585
        /**
 
586
         * The cluster center index
 
587
         */
 
588
        int pivot;
 
589
        /**
 
590
         * The cluster size (number of points in the cluster)
 
591
         */
 
592
        int size;
 
593
        /**
 
594
         * Child nodes (only for non-terminal nodes)
 
595
         */
 
596
        Node** childs;
 
597
        /**
 
598
         * Node points (only for terminal nodes)
 
599
         */
 
600
        int* indices;
 
601
        /**
 
602
         * Level
 
603
         */
 
604
        int level;
 
605
    };
 
606
    typedef Node* NodePtr;
 
607
 
 
608
 
 
609
 
 
610
    /**
 
611
     * Alias definition for a nicer syntax.
 
612
     */
 
613
    typedef BranchStruct<NodePtr, DistanceType> BranchSt;
 
614
 
 
615
 
 
616
 
 
617
    void save_tree(FILE* stream, NodePtr node, int num)
 
618
    {
 
619
        save_value(stream, *node);
 
620
        if (node->childs==NULL) {
 
621
            int indices_offset = (int)(node->indices - indices[num]);
 
622
            save_value(stream, indices_offset);
 
623
        }
 
624
        else {
 
625
            for(int i=0; i<branching_; ++i) {
 
626
                save_tree(stream, node->childs[i], num);
 
627
            }
 
628
        }
 
629
    }
 
630
 
 
631
 
 
632
    void load_tree(FILE* stream, NodePtr& node, int num)
 
633
    {
 
634
        node = pool.allocate<Node>();
 
635
        load_value(stream, *node);
 
636
        if (node->childs==NULL) {
 
637
            int indices_offset;
 
638
            load_value(stream, indices_offset);
 
639
            node->indices = indices[num] + indices_offset;
 
640
        }
 
641
        else {
 
642
            node->childs = pool.allocate<NodePtr>(branching_);
 
643
            for(int i=0; i<branching_; ++i) {
 
644
                load_tree(stream, node->childs[i], num);
 
645
            }
 
646
        }
 
647
    }
 
648
 
 
649
 
 
650
 
 
651
 
 
652
    void computeLabels(int* dsindices, int indices_length,  int* centers, int centers_length, int* labels, DistanceType& cost)
 
653
    {
 
654
        cost = 0;
 
655
        for (int i=0; i<indices_length; ++i) {
 
656
            ElementType* point = dataset[dsindices[i]];
 
657
            DistanceType dist = distance(point, dataset[centers[0]], veclen_);
 
658
            labels[i] = 0;
 
659
            for (int j=1; j<centers_length; ++j) {
 
660
                DistanceType new_dist = distance(point, dataset[centers[j]], veclen_);
 
661
                if (dist>new_dist) {
 
662
                    labels[i] = j;
 
663
                    dist = new_dist;
 
664
                }
 
665
            }
 
666
            cost += dist;
 
667
        }
 
668
    }
 
669
 
 
670
    /**
 
671
     * The method responsible with actually doing the recursive hierarchical
 
672
     * clustering
 
673
     *
 
674
     * Params:
 
675
     *     node = the node to cluster
 
676
     *     indices = indices of the points belonging to the current node
 
677
     *     branching = the branching factor to use in the clustering
 
678
     *
 
679
     * TODO: for 1-sized clusters don't store a cluster center (it's the same as the single cluster point)
 
680
     */
 
681
    void computeClustering(NodePtr node, int* dsindices, int indices_length, int branching, int level)
 
682
    {
 
683
        node->size = indices_length;
 
684
        node->level = level;
 
685
 
 
686
        if (indices_length < leaf_size_) { // leaf node
 
687
            node->indices = dsindices;
 
688
            std::sort(node->indices,node->indices+indices_length);
 
689
            node->childs = NULL;
 
690
            return;
 
691
        }
 
692
 
 
693
        std::vector<int> centers(branching);
 
694
        std::vector<int> labels(indices_length);
 
695
 
 
696
        int centers_length;
 
697
        (this->*chooseCenters)(branching, dsindices, indices_length, &centers[0], centers_length);
 
698
 
 
699
        if (centers_length<branching) {
 
700
            node->indices = dsindices;
 
701
            std::sort(node->indices,node->indices+indices_length);
 
702
            node->childs = NULL;
 
703
            return;
 
704
        }
 
705
 
 
706
 
 
707
        //      assign points to clusters
 
708
        DistanceType cost;
 
709
        computeLabels(dsindices, indices_length, &centers[0], centers_length, &labels[0], cost);
 
710
 
 
711
        node->childs = pool.allocate<NodePtr>(branching);
 
712
        int start = 0;
 
713
        int end = start;
 
714
        for (int i=0; i<branching; ++i) {
 
715
            for (int j=0; j<indices_length; ++j) {
 
716
                if (labels[j]==i) {
 
717
                    std::swap(dsindices[j],dsindices[end]);
 
718
                    std::swap(labels[j],labels[end]);
 
719
                    end++;
 
720
                }
 
721
            }
 
722
 
 
723
            node->childs[i] = pool.allocate<Node>();
 
724
            node->childs[i]->pivot = centers[i];
 
725
            node->childs[i]->indices = NULL;
 
726
            computeClustering(node->childs[i],dsindices+start, end-start, branching, level+1);
 
727
            start=end;
 
728
        }
 
729
    }
 
730
 
 
731
 
 
732
 
 
733
    /**
 
734
     * Performs one descent in the hierarchical k-means tree. The branches not
 
735
     * visited are stored in a priority queue.
 
736
     *
 
737
     * Params:
 
738
     *      node = node to explore
 
739
     *      result = container for the k-nearest neighbors found
 
740
     *      vec = query points
 
741
     *      checks = how many points in the dataset have been checked so far
 
742
     *      maxChecks = maximum dataset points to checks
 
743
     */
 
744
 
 
745
 
 
746
    void findNN(NodePtr node, ResultSet<DistanceType>& result, const ElementType* vec, int& checks, int maxChecks,
 
747
                Heap<BranchSt>* heap, std::vector<bool>& checked)
 
748
    {
 
749
        if (node->childs==NULL) {
 
750
            if (checks>=maxChecks) {
 
751
                if (result.full()) return;
 
752
            }
 
753
            for (int i=0; i<node->size; ++i) {
 
754
                int index = node->indices[i];
 
755
                if (!checked[index]) {
 
756
                    DistanceType dist = distance(dataset[index], vec, veclen_);
 
757
                    result.addPoint(dist, index);
 
758
                    checked[index] = true;
 
759
                    ++checks;
 
760
                }
 
761
            }
 
762
        }
 
763
        else {
 
764
            DistanceType* domain_distances = new DistanceType[branching_];
 
765
            int best_index = 0;
 
766
            domain_distances[best_index] = distance(vec, dataset[node->childs[best_index]->pivot], veclen_);
 
767
            for (int i=1; i<branching_; ++i) {
 
768
                domain_distances[i] = distance(vec, dataset[node->childs[i]->pivot], veclen_);
 
769
                if (domain_distances[i]<domain_distances[best_index]) {
 
770
                    best_index = i;
 
771
                }
 
772
            }
 
773
            for (int i=0; i<branching_; ++i) {
 
774
                if (i!=best_index) {
 
775
                    heap->insert(BranchSt(node->childs[i],domain_distances[i]));
 
776
                }
 
777
            }
 
778
            delete[] domain_distances;
 
779
            findNN(node->childs[best_index],result,vec, checks, maxChecks, heap, checked);
 
780
        }
 
781
    }
 
782
 
 
783
private:
 
784
 
 
785
 
 
786
    /**
 
787
     * The dataset used by this index
 
788
     */
 
789
    const Matrix<ElementType> dataset;
 
790
 
 
791
    /**
 
792
     * Parameters used by this index
 
793
     */
 
794
    IndexParams params;
 
795
 
 
796
 
 
797
    /**
 
798
     * Number of features in the dataset.
 
799
     */
 
800
    size_t size_;
 
801
 
 
802
    /**
 
803
     * Length of each feature.
 
804
     */
 
805
    size_t veclen_;
 
806
 
 
807
    /**
 
808
     * The root node in the tree.
 
809
     */
 
810
    NodePtr* root;
 
811
 
 
812
    /**
 
813
     *  Array of indices to vectors in the dataset.
 
814
     */
 
815
    int** indices;
 
816
 
 
817
 
 
818
    /**
 
819
     * The distance
 
820
     */
 
821
    Distance distance;
 
822
 
 
823
    /**
 
824
     * Pooled memory allocator.
 
825
     *
 
826
     * Using a pooled memory allocator is more efficient
 
827
     * than allocating memory directly when there is a large
 
828
     * number small of memory allocations.
 
829
     */
 
830
    PooledAllocator pool;
 
831
 
 
832
    /**
 
833
     * Memory occupied by the index.
 
834
     */
 
835
    int memoryCounter;
 
836
 
 
837
    /** index parameters */
 
838
    int branching_;
 
839
    int trees_;
 
840
    flann_centers_init_t centers_init_;
 
841
    int leaf_size_;
 
842
 
 
843
 
 
844
};
 
845
 
 
846
}
 
847
 
 
848
#endif /* OPENCV_FLANN_HIERARCHICAL_CLUSTERING_INDEX_H_ */