Renesas / opencv-lib

openCV library for Renesas RZ/A

Dependents:   RZ_A2M_Mbed_samples

include/opencv2/flann/kmeans_index.h@0:0e0631af0305, 2021-01-29 (annotated)

Committer:
RyoheiHagimoto
Date:
Fri Jan 29 04:53:38 2021 +0000
Revision:
0:0e0631af0305
copied from https://github.com/d-kato/opencv-lib.

Who changed what in which revision?

UserRevisionLine numberNew contents of line
RyoheiHagimoto 0:0e0631af0305 1 /***********************************************************************
RyoheiHagimoto 0:0e0631af0305 2 * Software License Agreement (BSD License)
RyoheiHagimoto 0:0e0631af0305 3 *
RyoheiHagimoto 0:0e0631af0305 4 * Copyright 2008-2009 Marius Muja (mariusm@cs.ubc.ca). All rights reserved.
RyoheiHagimoto 0:0e0631af0305 5 * Copyright 2008-2009 David G. Lowe (lowe@cs.ubc.ca). All rights reserved.
RyoheiHagimoto 0:0e0631af0305 6 *
RyoheiHagimoto 0:0e0631af0305 7 * THE BSD LICENSE
RyoheiHagimoto 0:0e0631af0305 8 *
RyoheiHagimoto 0:0e0631af0305 9 * Redistribution and use in source and binary forms, with or without
RyoheiHagimoto 0:0e0631af0305 10 * modification, are permitted provided that the following conditions
RyoheiHagimoto 0:0e0631af0305 11 * are met:
RyoheiHagimoto 0:0e0631af0305 12 *
RyoheiHagimoto 0:0e0631af0305 13 * 1. Redistributions of source code must retain the above copyright
RyoheiHagimoto 0:0e0631af0305 14 * notice, this list of conditions and the following disclaimer.
RyoheiHagimoto 0:0e0631af0305 15 * 2. Redistributions in binary form must reproduce the above copyright
RyoheiHagimoto 0:0e0631af0305 16 * notice, this list of conditions and the following disclaimer in the
RyoheiHagimoto 0:0e0631af0305 17 * documentation and/or other materials provided with the distribution.
RyoheiHagimoto 0:0e0631af0305 18 *
RyoheiHagimoto 0:0e0631af0305 19 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
RyoheiHagimoto 0:0e0631af0305 20 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
RyoheiHagimoto 0:0e0631af0305 21 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
RyoheiHagimoto 0:0e0631af0305 22 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
RyoheiHagimoto 0:0e0631af0305 23 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
RyoheiHagimoto 0:0e0631af0305 24 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
RyoheiHagimoto 0:0e0631af0305 25 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
RyoheiHagimoto 0:0e0631af0305 26 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
RyoheiHagimoto 0:0e0631af0305 27 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
RyoheiHagimoto 0:0e0631af0305 28 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
RyoheiHagimoto 0:0e0631af0305 29 *************************************************************************/
RyoheiHagimoto 0:0e0631af0305 30
RyoheiHagimoto 0:0e0631af0305 31 #ifndef OPENCV_FLANN_KMEANS_INDEX_H_
RyoheiHagimoto 0:0e0631af0305 32 #define OPENCV_FLANN_KMEANS_INDEX_H_
RyoheiHagimoto 0:0e0631af0305 33
RyoheiHagimoto 0:0e0631af0305 34 #include <algorithm>
RyoheiHagimoto 0:0e0631af0305 35 #include <map>
RyoheiHagimoto 0:0e0631af0305 36 #include <cassert>
RyoheiHagimoto 0:0e0631af0305 37 #include <limits>
RyoheiHagimoto 0:0e0631af0305 38 #include <cmath>
RyoheiHagimoto 0:0e0631af0305 39
RyoheiHagimoto 0:0e0631af0305 40 #include "general.h"
RyoheiHagimoto 0:0e0631af0305 41 #include "nn_index.h"
RyoheiHagimoto 0:0e0631af0305 42 #include "dist.h"
RyoheiHagimoto 0:0e0631af0305 43 #include "matrix.h"
RyoheiHagimoto 0:0e0631af0305 44 #include "result_set.h"
RyoheiHagimoto 0:0e0631af0305 45 #include "heap.h"
RyoheiHagimoto 0:0e0631af0305 46 #include "allocator.h"
RyoheiHagimoto 0:0e0631af0305 47 #include "random.h"
RyoheiHagimoto 0:0e0631af0305 48 #include "saving.h"
RyoheiHagimoto 0:0e0631af0305 49 #include "logger.h"
RyoheiHagimoto 0:0e0631af0305 50
RyoheiHagimoto 0:0e0631af0305 51
RyoheiHagimoto 0:0e0631af0305 52 namespace cvflann
RyoheiHagimoto 0:0e0631af0305 53 {
RyoheiHagimoto 0:0e0631af0305 54
RyoheiHagimoto 0:0e0631af0305 55 struct KMeansIndexParams : public IndexParams
RyoheiHagimoto 0:0e0631af0305 56 {
RyoheiHagimoto 0:0e0631af0305 57 KMeansIndexParams(int branching = 32, int iterations = 11,
RyoheiHagimoto 0:0e0631af0305 58 flann_centers_init_t centers_init = FLANN_CENTERS_RANDOM, float cb_index = 0.2 )
RyoheiHagimoto 0:0e0631af0305 59 {
RyoheiHagimoto 0:0e0631af0305 60 (*this)["algorithm"] = FLANN_INDEX_KMEANS;
RyoheiHagimoto 0:0e0631af0305 61 // branching factor
RyoheiHagimoto 0:0e0631af0305 62 (*this)["branching"] = branching;
RyoheiHagimoto 0:0e0631af0305 63 // max iterations to perform in one kmeans clustering (kmeans tree)
RyoheiHagimoto 0:0e0631af0305 64 (*this)["iterations"] = iterations;
RyoheiHagimoto 0:0e0631af0305 65 // algorithm used for picking the initial cluster centers for kmeans tree
RyoheiHagimoto 0:0e0631af0305 66 (*this)["centers_init"] = centers_init;
RyoheiHagimoto 0:0e0631af0305 67 // cluster boundary index. Used when searching the kmeans tree
RyoheiHagimoto 0:0e0631af0305 68 (*this)["cb_index"] = cb_index;
RyoheiHagimoto 0:0e0631af0305 69 }
RyoheiHagimoto 0:0e0631af0305 70 };
RyoheiHagimoto 0:0e0631af0305 71
RyoheiHagimoto 0:0e0631af0305 72
RyoheiHagimoto 0:0e0631af0305 73 /**
RyoheiHagimoto 0:0e0631af0305 74 * Hierarchical kmeans index
RyoheiHagimoto 0:0e0631af0305 75 *
RyoheiHagimoto 0:0e0631af0305 76 * Contains a tree constructed through a hierarchical kmeans clustering
RyoheiHagimoto 0:0e0631af0305 77 * and other information for indexing a set of points for nearest-neighbour matching.
RyoheiHagimoto 0:0e0631af0305 78 */
RyoheiHagimoto 0:0e0631af0305 79 template <typename Distance>
RyoheiHagimoto 0:0e0631af0305 80 class KMeansIndex : public NNIndex<Distance>
RyoheiHagimoto 0:0e0631af0305 81 {
RyoheiHagimoto 0:0e0631af0305 82 public:
RyoheiHagimoto 0:0e0631af0305 83 typedef typename Distance::ElementType ElementType;
RyoheiHagimoto 0:0e0631af0305 84 typedef typename Distance::ResultType DistanceType;
RyoheiHagimoto 0:0e0631af0305 85
RyoheiHagimoto 0:0e0631af0305 86
RyoheiHagimoto 0:0e0631af0305 87
RyoheiHagimoto 0:0e0631af0305 88 typedef void (KMeansIndex::* centersAlgFunction)(int, int*, int, int*, int&);
RyoheiHagimoto 0:0e0631af0305 89
RyoheiHagimoto 0:0e0631af0305 90 /**
RyoheiHagimoto 0:0e0631af0305 91 * The function used for choosing the cluster centers.
RyoheiHagimoto 0:0e0631af0305 92 */
RyoheiHagimoto 0:0e0631af0305 93 centersAlgFunction chooseCenters;
RyoheiHagimoto 0:0e0631af0305 94
RyoheiHagimoto 0:0e0631af0305 95
RyoheiHagimoto 0:0e0631af0305 96
RyoheiHagimoto 0:0e0631af0305 97 /**
RyoheiHagimoto 0:0e0631af0305 98 * Chooses the initial centers in the k-means clustering in a random manner.
RyoheiHagimoto 0:0e0631af0305 99 *
RyoheiHagimoto 0:0e0631af0305 100 * Params:
RyoheiHagimoto 0:0e0631af0305 101 * k = number of centers
RyoheiHagimoto 0:0e0631af0305 102 * vecs = the dataset of points
RyoheiHagimoto 0:0e0631af0305 103 * indices = indices in the dataset
RyoheiHagimoto 0:0e0631af0305 104 * indices_length = length of indices vector
RyoheiHagimoto 0:0e0631af0305 105 *
RyoheiHagimoto 0:0e0631af0305 106 */
RyoheiHagimoto 0:0e0631af0305 107 void chooseCentersRandom(int k, int* indices, int indices_length, int* centers, int& centers_length)
RyoheiHagimoto 0:0e0631af0305 108 {
RyoheiHagimoto 0:0e0631af0305 109 UniqueRandom r(indices_length);
RyoheiHagimoto 0:0e0631af0305 110
RyoheiHagimoto 0:0e0631af0305 111 int index;
RyoheiHagimoto 0:0e0631af0305 112 for (index=0; index<k; ++index) {
RyoheiHagimoto 0:0e0631af0305 113 bool duplicate = true;
RyoheiHagimoto 0:0e0631af0305 114 int rnd;
RyoheiHagimoto 0:0e0631af0305 115 while (duplicate) {
RyoheiHagimoto 0:0e0631af0305 116 duplicate = false;
RyoheiHagimoto 0:0e0631af0305 117 rnd = r.next();
RyoheiHagimoto 0:0e0631af0305 118 if (rnd<0) {
RyoheiHagimoto 0:0e0631af0305 119 centers_length = index;
RyoheiHagimoto 0:0e0631af0305 120 return;
RyoheiHagimoto 0:0e0631af0305 121 }
RyoheiHagimoto 0:0e0631af0305 122
RyoheiHagimoto 0:0e0631af0305 123 centers[index] = indices[rnd];
RyoheiHagimoto 0:0e0631af0305 124
RyoheiHagimoto 0:0e0631af0305 125 for (int j=0; j<index; ++j) {
RyoheiHagimoto 0:0e0631af0305 126 DistanceType sq = distance_(dataset_[centers[index]], dataset_[centers[j]], dataset_.cols);
RyoheiHagimoto 0:0e0631af0305 127 if (sq<1e-16) {
RyoheiHagimoto 0:0e0631af0305 128 duplicate = true;
RyoheiHagimoto 0:0e0631af0305 129 }
RyoheiHagimoto 0:0e0631af0305 130 }
RyoheiHagimoto 0:0e0631af0305 131 }
RyoheiHagimoto 0:0e0631af0305 132 }
RyoheiHagimoto 0:0e0631af0305 133
RyoheiHagimoto 0:0e0631af0305 134 centers_length = index;
RyoheiHagimoto 0:0e0631af0305 135 }
RyoheiHagimoto 0:0e0631af0305 136
RyoheiHagimoto 0:0e0631af0305 137
RyoheiHagimoto 0:0e0631af0305 138 /**
RyoheiHagimoto 0:0e0631af0305 139 * Chooses the initial centers in the k-means using Gonzales' algorithm
RyoheiHagimoto 0:0e0631af0305 140 * so that the centers are spaced apart from each other.
RyoheiHagimoto 0:0e0631af0305 141 *
RyoheiHagimoto 0:0e0631af0305 142 * Params:
RyoheiHagimoto 0:0e0631af0305 143 * k = number of centers
RyoheiHagimoto 0:0e0631af0305 144 * vecs = the dataset of points
RyoheiHagimoto 0:0e0631af0305 145 * indices = indices in the dataset
RyoheiHagimoto 0:0e0631af0305 146 * Returns:
RyoheiHagimoto 0:0e0631af0305 147 */
RyoheiHagimoto 0:0e0631af0305 148 void chooseCentersGonzales(int k, int* indices, int indices_length, int* centers, int& centers_length)
RyoheiHagimoto 0:0e0631af0305 149 {
RyoheiHagimoto 0:0e0631af0305 150 int n = indices_length;
RyoheiHagimoto 0:0e0631af0305 151
RyoheiHagimoto 0:0e0631af0305 152 int rnd = rand_int(n);
RyoheiHagimoto 0:0e0631af0305 153 assert(rnd >=0 && rnd < n);
RyoheiHagimoto 0:0e0631af0305 154
RyoheiHagimoto 0:0e0631af0305 155 centers[0] = indices[rnd];
RyoheiHagimoto 0:0e0631af0305 156
RyoheiHagimoto 0:0e0631af0305 157 int index;
RyoheiHagimoto 0:0e0631af0305 158 for (index=1; index<k; ++index) {
RyoheiHagimoto 0:0e0631af0305 159
RyoheiHagimoto 0:0e0631af0305 160 int best_index = -1;
RyoheiHagimoto 0:0e0631af0305 161 DistanceType best_val = 0;
RyoheiHagimoto 0:0e0631af0305 162 for (int j=0; j<n; ++j) {
RyoheiHagimoto 0:0e0631af0305 163 DistanceType dist = distance_(dataset_[centers[0]],dataset_[indices[j]],dataset_.cols);
RyoheiHagimoto 0:0e0631af0305 164 for (int i=1; i<index; ++i) {
RyoheiHagimoto 0:0e0631af0305 165 DistanceType tmp_dist = distance_(dataset_[centers[i]],dataset_[indices[j]],dataset_.cols);
RyoheiHagimoto 0:0e0631af0305 166 if (tmp_dist<dist) {
RyoheiHagimoto 0:0e0631af0305 167 dist = tmp_dist;
RyoheiHagimoto 0:0e0631af0305 168 }
RyoheiHagimoto 0:0e0631af0305 169 }
RyoheiHagimoto 0:0e0631af0305 170 if (dist>best_val) {
RyoheiHagimoto 0:0e0631af0305 171 best_val = dist;
RyoheiHagimoto 0:0e0631af0305 172 best_index = j;
RyoheiHagimoto 0:0e0631af0305 173 }
RyoheiHagimoto 0:0e0631af0305 174 }
RyoheiHagimoto 0:0e0631af0305 175 if (best_index!=-1) {
RyoheiHagimoto 0:0e0631af0305 176 centers[index] = indices[best_index];
RyoheiHagimoto 0:0e0631af0305 177 }
RyoheiHagimoto 0:0e0631af0305 178 else {
RyoheiHagimoto 0:0e0631af0305 179 break;
RyoheiHagimoto 0:0e0631af0305 180 }
RyoheiHagimoto 0:0e0631af0305 181 }
RyoheiHagimoto 0:0e0631af0305 182 centers_length = index;
RyoheiHagimoto 0:0e0631af0305 183 }
RyoheiHagimoto 0:0e0631af0305 184
RyoheiHagimoto 0:0e0631af0305 185
RyoheiHagimoto 0:0e0631af0305 186 /**
RyoheiHagimoto 0:0e0631af0305 187 * Chooses the initial centers in the k-means using the algorithm
RyoheiHagimoto 0:0e0631af0305 188 * proposed in the KMeans++ paper:
RyoheiHagimoto 0:0e0631af0305 189 * Arthur, David; Vassilvitskii, Sergei - k-means++: The Advantages of Careful Seeding
RyoheiHagimoto 0:0e0631af0305 190 *
RyoheiHagimoto 0:0e0631af0305 191 * Implementation of this function was converted from the one provided in Arthur's code.
RyoheiHagimoto 0:0e0631af0305 192 *
RyoheiHagimoto 0:0e0631af0305 193 * Params:
RyoheiHagimoto 0:0e0631af0305 194 * k = number of centers
RyoheiHagimoto 0:0e0631af0305 195 * vecs = the dataset of points
RyoheiHagimoto 0:0e0631af0305 196 * indices = indices in the dataset
RyoheiHagimoto 0:0e0631af0305 197 * Returns:
RyoheiHagimoto 0:0e0631af0305 198 */
RyoheiHagimoto 0:0e0631af0305 199 void chooseCentersKMeanspp(int k, int* indices, int indices_length, int* centers, int& centers_length)
RyoheiHagimoto 0:0e0631af0305 200 {
RyoheiHagimoto 0:0e0631af0305 201 int n = indices_length;
RyoheiHagimoto 0:0e0631af0305 202
RyoheiHagimoto 0:0e0631af0305 203 double currentPot = 0;
RyoheiHagimoto 0:0e0631af0305 204 DistanceType* closestDistSq = new DistanceType[n];
RyoheiHagimoto 0:0e0631af0305 205
RyoheiHagimoto 0:0e0631af0305 206 // Choose one random center and set the closestDistSq values
RyoheiHagimoto 0:0e0631af0305 207 int index = rand_int(n);
RyoheiHagimoto 0:0e0631af0305 208 assert(index >=0 && index < n);
RyoheiHagimoto 0:0e0631af0305 209 centers[0] = indices[index];
RyoheiHagimoto 0:0e0631af0305 210
RyoheiHagimoto 0:0e0631af0305 211 for (int i = 0; i < n; i++) {
RyoheiHagimoto 0:0e0631af0305 212 closestDistSq[i] = distance_(dataset_[indices[i]], dataset_[indices[index]], dataset_.cols);
RyoheiHagimoto 0:0e0631af0305 213 closestDistSq[i] = ensureSquareDistance<Distance>( closestDistSq[i] );
RyoheiHagimoto 0:0e0631af0305 214 currentPot += closestDistSq[i];
RyoheiHagimoto 0:0e0631af0305 215 }
RyoheiHagimoto 0:0e0631af0305 216
RyoheiHagimoto 0:0e0631af0305 217
RyoheiHagimoto 0:0e0631af0305 218 const int numLocalTries = 1;
RyoheiHagimoto 0:0e0631af0305 219
RyoheiHagimoto 0:0e0631af0305 220 // Choose each center
RyoheiHagimoto 0:0e0631af0305 221 int centerCount;
RyoheiHagimoto 0:0e0631af0305 222 for (centerCount = 1; centerCount < k; centerCount++) {
RyoheiHagimoto 0:0e0631af0305 223
RyoheiHagimoto 0:0e0631af0305 224 // Repeat several trials
RyoheiHagimoto 0:0e0631af0305 225 double bestNewPot = -1;
RyoheiHagimoto 0:0e0631af0305 226 int bestNewIndex = -1;
RyoheiHagimoto 0:0e0631af0305 227 for (int localTrial = 0; localTrial < numLocalTries; localTrial++) {
RyoheiHagimoto 0:0e0631af0305 228
RyoheiHagimoto 0:0e0631af0305 229 // Choose our center - have to be slightly careful to return a valid answer even accounting
RyoheiHagimoto 0:0e0631af0305 230 // for possible rounding errors
RyoheiHagimoto 0:0e0631af0305 231 double randVal = rand_double(currentPot);
RyoheiHagimoto 0:0e0631af0305 232 for (index = 0; index < n-1; index++) {
RyoheiHagimoto 0:0e0631af0305 233 if (randVal <= closestDistSq[index]) break;
RyoheiHagimoto 0:0e0631af0305 234 else randVal -= closestDistSq[index];
RyoheiHagimoto 0:0e0631af0305 235 }
RyoheiHagimoto 0:0e0631af0305 236
RyoheiHagimoto 0:0e0631af0305 237 // Compute the new potential
RyoheiHagimoto 0:0e0631af0305 238 double newPot = 0;
RyoheiHagimoto 0:0e0631af0305 239 for (int i = 0; i < n; i++) {
RyoheiHagimoto 0:0e0631af0305 240 DistanceType dist = distance_(dataset_[indices[i]], dataset_[indices[index]], dataset_.cols);
RyoheiHagimoto 0:0e0631af0305 241 newPot += std::min( ensureSquareDistance<Distance>(dist), closestDistSq[i] );
RyoheiHagimoto 0:0e0631af0305 242 }
RyoheiHagimoto 0:0e0631af0305 243
RyoheiHagimoto 0:0e0631af0305 244 // Store the best result
RyoheiHagimoto 0:0e0631af0305 245 if ((bestNewPot < 0)||(newPot < bestNewPot)) {
RyoheiHagimoto 0:0e0631af0305 246 bestNewPot = newPot;
RyoheiHagimoto 0:0e0631af0305 247 bestNewIndex = index;
RyoheiHagimoto 0:0e0631af0305 248 }
RyoheiHagimoto 0:0e0631af0305 249 }
RyoheiHagimoto 0:0e0631af0305 250
RyoheiHagimoto 0:0e0631af0305 251 // Add the appropriate center
RyoheiHagimoto 0:0e0631af0305 252 centers[centerCount] = indices[bestNewIndex];
RyoheiHagimoto 0:0e0631af0305 253 currentPot = bestNewPot;
RyoheiHagimoto 0:0e0631af0305 254 for (int i = 0; i < n; i++) {
RyoheiHagimoto 0:0e0631af0305 255 DistanceType dist = distance_(dataset_[indices[i]], dataset_[indices[bestNewIndex]], dataset_.cols);
RyoheiHagimoto 0:0e0631af0305 256 closestDistSq[i] = std::min( ensureSquareDistance<Distance>(dist), closestDistSq[i] );
RyoheiHagimoto 0:0e0631af0305 257 }
RyoheiHagimoto 0:0e0631af0305 258 }
RyoheiHagimoto 0:0e0631af0305 259
RyoheiHagimoto 0:0e0631af0305 260 centers_length = centerCount;
RyoheiHagimoto 0:0e0631af0305 261
RyoheiHagimoto 0:0e0631af0305 262 delete[] closestDistSq;
RyoheiHagimoto 0:0e0631af0305 263 }
RyoheiHagimoto 0:0e0631af0305 264
RyoheiHagimoto 0:0e0631af0305 265
RyoheiHagimoto 0:0e0631af0305 266
RyoheiHagimoto 0:0e0631af0305 267 public:
RyoheiHagimoto 0:0e0631af0305 268
RyoheiHagimoto 0:0e0631af0305 269 flann_algorithm_t getType() const
RyoheiHagimoto 0:0e0631af0305 270 {
RyoheiHagimoto 0:0e0631af0305 271 return FLANN_INDEX_KMEANS;
RyoheiHagimoto 0:0e0631af0305 272 }
RyoheiHagimoto 0:0e0631af0305 273
RyoheiHagimoto 0:0e0631af0305 274 class KMeansDistanceComputer : public cv::ParallelLoopBody
RyoheiHagimoto 0:0e0631af0305 275 {
RyoheiHagimoto 0:0e0631af0305 276 public:
RyoheiHagimoto 0:0e0631af0305 277 KMeansDistanceComputer(Distance _distance, const Matrix<ElementType>& _dataset,
RyoheiHagimoto 0:0e0631af0305 278 const int _branching, const int* _indices, const Matrix<double>& _dcenters, const size_t _veclen,
RyoheiHagimoto 0:0e0631af0305 279 int* _count, int* _belongs_to, std::vector<DistanceType>& _radiuses, bool& _converged, cv::Mutex& _mtx)
RyoheiHagimoto 0:0e0631af0305 280 : distance(_distance)
RyoheiHagimoto 0:0e0631af0305 281 , dataset(_dataset)
RyoheiHagimoto 0:0e0631af0305 282 , branching(_branching)
RyoheiHagimoto 0:0e0631af0305 283 , indices(_indices)
RyoheiHagimoto 0:0e0631af0305 284 , dcenters(_dcenters)
RyoheiHagimoto 0:0e0631af0305 285 , veclen(_veclen)
RyoheiHagimoto 0:0e0631af0305 286 , count(_count)
RyoheiHagimoto 0:0e0631af0305 287 , belongs_to(_belongs_to)
RyoheiHagimoto 0:0e0631af0305 288 , radiuses(_radiuses)
RyoheiHagimoto 0:0e0631af0305 289 , converged(_converged)
RyoheiHagimoto 0:0e0631af0305 290 , mtx(_mtx)
RyoheiHagimoto 0:0e0631af0305 291 {
RyoheiHagimoto 0:0e0631af0305 292 }
RyoheiHagimoto 0:0e0631af0305 293
RyoheiHagimoto 0:0e0631af0305 294 void operator()(const cv::Range& range) const
RyoheiHagimoto 0:0e0631af0305 295 {
RyoheiHagimoto 0:0e0631af0305 296 const int begin = range.start;
RyoheiHagimoto 0:0e0631af0305 297 const int end = range.end;
RyoheiHagimoto 0:0e0631af0305 298
RyoheiHagimoto 0:0e0631af0305 299 for( int i = begin; i<end; ++i)
RyoheiHagimoto 0:0e0631af0305 300 {
RyoheiHagimoto 0:0e0631af0305 301 DistanceType sq_dist = distance(dataset[indices[i]], dcenters[0], veclen);
RyoheiHagimoto 0:0e0631af0305 302 int new_centroid = 0;
RyoheiHagimoto 0:0e0631af0305 303 for (int j=1; j<branching; ++j) {
RyoheiHagimoto 0:0e0631af0305 304 DistanceType new_sq_dist = distance(dataset[indices[i]], dcenters[j], veclen);
RyoheiHagimoto 0:0e0631af0305 305 if (sq_dist>new_sq_dist) {
RyoheiHagimoto 0:0e0631af0305 306 new_centroid = j;
RyoheiHagimoto 0:0e0631af0305 307 sq_dist = new_sq_dist;
RyoheiHagimoto 0:0e0631af0305 308 }
RyoheiHagimoto 0:0e0631af0305 309 }
RyoheiHagimoto 0:0e0631af0305 310 if (sq_dist > radiuses[new_centroid]) {
RyoheiHagimoto 0:0e0631af0305 311 radiuses[new_centroid] = sq_dist;
RyoheiHagimoto 0:0e0631af0305 312 }
RyoheiHagimoto 0:0e0631af0305 313 if (new_centroid != belongs_to[i]) {
RyoheiHagimoto 0:0e0631af0305 314 count[belongs_to[i]]--;
RyoheiHagimoto 0:0e0631af0305 315 count[new_centroid]++;
RyoheiHagimoto 0:0e0631af0305 316 belongs_to[i] = new_centroid;
RyoheiHagimoto 0:0e0631af0305 317 mtx.lock();
RyoheiHagimoto 0:0e0631af0305 318 converged = false;
RyoheiHagimoto 0:0e0631af0305 319 mtx.unlock();
RyoheiHagimoto 0:0e0631af0305 320 }
RyoheiHagimoto 0:0e0631af0305 321 }
RyoheiHagimoto 0:0e0631af0305 322 }
RyoheiHagimoto 0:0e0631af0305 323
RyoheiHagimoto 0:0e0631af0305 324 private:
RyoheiHagimoto 0:0e0631af0305 325 Distance distance;
RyoheiHagimoto 0:0e0631af0305 326 const Matrix<ElementType>& dataset;
RyoheiHagimoto 0:0e0631af0305 327 const int branching;
RyoheiHagimoto 0:0e0631af0305 328 const int* indices;
RyoheiHagimoto 0:0e0631af0305 329 const Matrix<double>& dcenters;
RyoheiHagimoto 0:0e0631af0305 330 const size_t veclen;
RyoheiHagimoto 0:0e0631af0305 331 int* count;
RyoheiHagimoto 0:0e0631af0305 332 int* belongs_to;
RyoheiHagimoto 0:0e0631af0305 333 std::vector<DistanceType>& radiuses;
RyoheiHagimoto 0:0e0631af0305 334 bool& converged;
RyoheiHagimoto 0:0e0631af0305 335 cv::Mutex& mtx;
RyoheiHagimoto 0:0e0631af0305 336 KMeansDistanceComputer& operator=( const KMeansDistanceComputer & ) { return *this; }
RyoheiHagimoto 0:0e0631af0305 337 };
RyoheiHagimoto 0:0e0631af0305 338
RyoheiHagimoto 0:0e0631af0305 339 /**
RyoheiHagimoto 0:0e0631af0305 340 * Index constructor
RyoheiHagimoto 0:0e0631af0305 341 *
RyoheiHagimoto 0:0e0631af0305 342 * Params:
RyoheiHagimoto 0:0e0631af0305 343 * inputData = dataset with the input features
RyoheiHagimoto 0:0e0631af0305 344 * params = parameters passed to the hierarchical k-means algorithm
RyoheiHagimoto 0:0e0631af0305 345 */
RyoheiHagimoto 0:0e0631af0305 346 KMeansIndex(const Matrix<ElementType>& inputData, const IndexParams& params = KMeansIndexParams(),
RyoheiHagimoto 0:0e0631af0305 347 Distance d = Distance())
RyoheiHagimoto 0:0e0631af0305 348 : dataset_(inputData), index_params_(params), root_(NULL), indices_(NULL), distance_(d)
RyoheiHagimoto 0:0e0631af0305 349 {
RyoheiHagimoto 0:0e0631af0305 350 memoryCounter_ = 0;
RyoheiHagimoto 0:0e0631af0305 351
RyoheiHagimoto 0:0e0631af0305 352 size_ = dataset_.rows;
RyoheiHagimoto 0:0e0631af0305 353 veclen_ = dataset_.cols;
RyoheiHagimoto 0:0e0631af0305 354
RyoheiHagimoto 0:0e0631af0305 355 branching_ = get_param(params,"branching",32);
RyoheiHagimoto 0:0e0631af0305 356 iterations_ = get_param(params,"iterations",11);
RyoheiHagimoto 0:0e0631af0305 357 if (iterations_<0) {
RyoheiHagimoto 0:0e0631af0305 358 iterations_ = (std::numeric_limits<int>::max)();
RyoheiHagimoto 0:0e0631af0305 359 }
RyoheiHagimoto 0:0e0631af0305 360 centers_init_ = get_param(params,"centers_init",FLANN_CENTERS_RANDOM);
RyoheiHagimoto 0:0e0631af0305 361
RyoheiHagimoto 0:0e0631af0305 362 if (centers_init_==FLANN_CENTERS_RANDOM) {
RyoheiHagimoto 0:0e0631af0305 363 chooseCenters = &KMeansIndex::chooseCentersRandom;
RyoheiHagimoto 0:0e0631af0305 364 }
RyoheiHagimoto 0:0e0631af0305 365 else if (centers_init_==FLANN_CENTERS_GONZALES) {
RyoheiHagimoto 0:0e0631af0305 366 chooseCenters = &KMeansIndex::chooseCentersGonzales;
RyoheiHagimoto 0:0e0631af0305 367 }
RyoheiHagimoto 0:0e0631af0305 368 else if (centers_init_==FLANN_CENTERS_KMEANSPP) {
RyoheiHagimoto 0:0e0631af0305 369 chooseCenters = &KMeansIndex::chooseCentersKMeanspp;
RyoheiHagimoto 0:0e0631af0305 370 }
RyoheiHagimoto 0:0e0631af0305 371 else {
RyoheiHagimoto 0:0e0631af0305 372 throw FLANNException("Unknown algorithm for choosing initial centers.");
RyoheiHagimoto 0:0e0631af0305 373 }
RyoheiHagimoto 0:0e0631af0305 374 cb_index_ = 0.4f;
RyoheiHagimoto 0:0e0631af0305 375
RyoheiHagimoto 0:0e0631af0305 376 }
RyoheiHagimoto 0:0e0631af0305 377
RyoheiHagimoto 0:0e0631af0305 378
RyoheiHagimoto 0:0e0631af0305 379 KMeansIndex(const KMeansIndex&);
RyoheiHagimoto 0:0e0631af0305 380 KMeansIndex& operator=(const KMeansIndex&);
RyoheiHagimoto 0:0e0631af0305 381
RyoheiHagimoto 0:0e0631af0305 382
RyoheiHagimoto 0:0e0631af0305 383 /**
RyoheiHagimoto 0:0e0631af0305 384 * Index destructor.
RyoheiHagimoto 0:0e0631af0305 385 *
RyoheiHagimoto 0:0e0631af0305 386 * Release the memory used by the index.
RyoheiHagimoto 0:0e0631af0305 387 */
RyoheiHagimoto 0:0e0631af0305 388 virtual ~KMeansIndex()
RyoheiHagimoto 0:0e0631af0305 389 {
RyoheiHagimoto 0:0e0631af0305 390 if (root_ != NULL) {
RyoheiHagimoto 0:0e0631af0305 391 free_centers(root_);
RyoheiHagimoto 0:0e0631af0305 392 }
RyoheiHagimoto 0:0e0631af0305 393 if (indices_!=NULL) {
RyoheiHagimoto 0:0e0631af0305 394 delete[] indices_;
RyoheiHagimoto 0:0e0631af0305 395 }
RyoheiHagimoto 0:0e0631af0305 396 }
RyoheiHagimoto 0:0e0631af0305 397
RyoheiHagimoto 0:0e0631af0305 398 /**
RyoheiHagimoto 0:0e0631af0305 399 * Returns size of index.
RyoheiHagimoto 0:0e0631af0305 400 */
RyoheiHagimoto 0:0e0631af0305 401 size_t size() const
RyoheiHagimoto 0:0e0631af0305 402 {
RyoheiHagimoto 0:0e0631af0305 403 return size_;
RyoheiHagimoto 0:0e0631af0305 404 }
RyoheiHagimoto 0:0e0631af0305 405
RyoheiHagimoto 0:0e0631af0305 406 /**
RyoheiHagimoto 0:0e0631af0305 407 * Returns the length of an index feature.
RyoheiHagimoto 0:0e0631af0305 408 */
RyoheiHagimoto 0:0e0631af0305 409 size_t veclen() const
RyoheiHagimoto 0:0e0631af0305 410 {
RyoheiHagimoto 0:0e0631af0305 411 return veclen_;
RyoheiHagimoto 0:0e0631af0305 412 }
RyoheiHagimoto 0:0e0631af0305 413
RyoheiHagimoto 0:0e0631af0305 414
RyoheiHagimoto 0:0e0631af0305 415 void set_cb_index( float index)
RyoheiHagimoto 0:0e0631af0305 416 {
RyoheiHagimoto 0:0e0631af0305 417 cb_index_ = index;
RyoheiHagimoto 0:0e0631af0305 418 }
RyoheiHagimoto 0:0e0631af0305 419
RyoheiHagimoto 0:0e0631af0305 420 /**
RyoheiHagimoto 0:0e0631af0305 421 * Computes the inde memory usage
RyoheiHagimoto 0:0e0631af0305 422 * Returns: memory used by the index
RyoheiHagimoto 0:0e0631af0305 423 */
RyoheiHagimoto 0:0e0631af0305 424 int usedMemory() const
RyoheiHagimoto 0:0e0631af0305 425 {
RyoheiHagimoto 0:0e0631af0305 426 return pool_.usedMemory+pool_.wastedMemory+memoryCounter_;
RyoheiHagimoto 0:0e0631af0305 427 }
RyoheiHagimoto 0:0e0631af0305 428
RyoheiHagimoto 0:0e0631af0305 429 /**
RyoheiHagimoto 0:0e0631af0305 430 * Builds the index
RyoheiHagimoto 0:0e0631af0305 431 */
RyoheiHagimoto 0:0e0631af0305 432 void buildIndex()
RyoheiHagimoto 0:0e0631af0305 433 {
RyoheiHagimoto 0:0e0631af0305 434 if (branching_<2) {
RyoheiHagimoto 0:0e0631af0305 435 throw FLANNException("Branching factor must be at least 2");
RyoheiHagimoto 0:0e0631af0305 436 }
RyoheiHagimoto 0:0e0631af0305 437
RyoheiHagimoto 0:0e0631af0305 438 indices_ = new int[size_];
RyoheiHagimoto 0:0e0631af0305 439 for (size_t i=0; i<size_; ++i) {
RyoheiHagimoto 0:0e0631af0305 440 indices_[i] = int(i);
RyoheiHagimoto 0:0e0631af0305 441 }
RyoheiHagimoto 0:0e0631af0305 442
RyoheiHagimoto 0:0e0631af0305 443 root_ = pool_.allocate<KMeansNode>();
RyoheiHagimoto 0:0e0631af0305 444 std::memset(root_, 0, sizeof(KMeansNode));
RyoheiHagimoto 0:0e0631af0305 445
RyoheiHagimoto 0:0e0631af0305 446 computeNodeStatistics(root_, indices_, (int)size_);
RyoheiHagimoto 0:0e0631af0305 447 computeClustering(root_, indices_, (int)size_, branching_,0);
RyoheiHagimoto 0:0e0631af0305 448 }
RyoheiHagimoto 0:0e0631af0305 449
RyoheiHagimoto 0:0e0631af0305 450
RyoheiHagimoto 0:0e0631af0305 451 void saveIndex(FILE* stream)
RyoheiHagimoto 0:0e0631af0305 452 {
RyoheiHagimoto 0:0e0631af0305 453 save_value(stream, branching_);
RyoheiHagimoto 0:0e0631af0305 454 save_value(stream, iterations_);
RyoheiHagimoto 0:0e0631af0305 455 save_value(stream, memoryCounter_);
RyoheiHagimoto 0:0e0631af0305 456 save_value(stream, cb_index_);
RyoheiHagimoto 0:0e0631af0305 457 save_value(stream, *indices_, (int)size_);
RyoheiHagimoto 0:0e0631af0305 458
RyoheiHagimoto 0:0e0631af0305 459 save_tree(stream, root_);
RyoheiHagimoto 0:0e0631af0305 460 }
RyoheiHagimoto 0:0e0631af0305 461
RyoheiHagimoto 0:0e0631af0305 462
RyoheiHagimoto 0:0e0631af0305 463 void loadIndex(FILE* stream)
RyoheiHagimoto 0:0e0631af0305 464 {
RyoheiHagimoto 0:0e0631af0305 465 load_value(stream, branching_);
RyoheiHagimoto 0:0e0631af0305 466 load_value(stream, iterations_);
RyoheiHagimoto 0:0e0631af0305 467 load_value(stream, memoryCounter_);
RyoheiHagimoto 0:0e0631af0305 468 load_value(stream, cb_index_);
RyoheiHagimoto 0:0e0631af0305 469 if (indices_!=NULL) {
RyoheiHagimoto 0:0e0631af0305 470 delete[] indices_;
RyoheiHagimoto 0:0e0631af0305 471 }
RyoheiHagimoto 0:0e0631af0305 472 indices_ = new int[size_];
RyoheiHagimoto 0:0e0631af0305 473 load_value(stream, *indices_, size_);
RyoheiHagimoto 0:0e0631af0305 474
RyoheiHagimoto 0:0e0631af0305 475 if (root_!=NULL) {
RyoheiHagimoto 0:0e0631af0305 476 free_centers(root_);
RyoheiHagimoto 0:0e0631af0305 477 }
RyoheiHagimoto 0:0e0631af0305 478 load_tree(stream, root_);
RyoheiHagimoto 0:0e0631af0305 479
RyoheiHagimoto 0:0e0631af0305 480 index_params_["algorithm"] = getType();
RyoheiHagimoto 0:0e0631af0305 481 index_params_["branching"] = branching_;
RyoheiHagimoto 0:0e0631af0305 482 index_params_["iterations"] = iterations_;
RyoheiHagimoto 0:0e0631af0305 483 index_params_["centers_init"] = centers_init_;
RyoheiHagimoto 0:0e0631af0305 484 index_params_["cb_index"] = cb_index_;
RyoheiHagimoto 0:0e0631af0305 485
RyoheiHagimoto 0:0e0631af0305 486 }
RyoheiHagimoto 0:0e0631af0305 487
RyoheiHagimoto 0:0e0631af0305 488
RyoheiHagimoto 0:0e0631af0305 489 /**
RyoheiHagimoto 0:0e0631af0305 490 * Find set of nearest neighbors to vec. Their indices are stored inside
RyoheiHagimoto 0:0e0631af0305 491 * the result object.
RyoheiHagimoto 0:0e0631af0305 492 *
RyoheiHagimoto 0:0e0631af0305 493 * Params:
RyoheiHagimoto 0:0e0631af0305 494 * result = the result object in which the indices of the nearest-neighbors are stored
RyoheiHagimoto 0:0e0631af0305 495 * vec = the vector for which to search the nearest neighbors
RyoheiHagimoto 0:0e0631af0305 496 * searchParams = parameters that influence the search algorithm (checks, cb_index)
RyoheiHagimoto 0:0e0631af0305 497 */
RyoheiHagimoto 0:0e0631af0305 498 void findNeighbors(ResultSet<DistanceType>& result, const ElementType* vec, const SearchParams& searchParams)
RyoheiHagimoto 0:0e0631af0305 499 {
RyoheiHagimoto 0:0e0631af0305 500
RyoheiHagimoto 0:0e0631af0305 501 int maxChecks = get_param(searchParams,"checks",32);
RyoheiHagimoto 0:0e0631af0305 502
RyoheiHagimoto 0:0e0631af0305 503 if (maxChecks==FLANN_CHECKS_UNLIMITED) {
RyoheiHagimoto 0:0e0631af0305 504 findExactNN(root_, result, vec);
RyoheiHagimoto 0:0e0631af0305 505 }
RyoheiHagimoto 0:0e0631af0305 506 else {
RyoheiHagimoto 0:0e0631af0305 507 // Priority queue storing intermediate branches in the best-bin-first search
RyoheiHagimoto 0:0e0631af0305 508 Heap<BranchSt>* heap = new Heap<BranchSt>((int)size_);
RyoheiHagimoto 0:0e0631af0305 509
RyoheiHagimoto 0:0e0631af0305 510 int checks = 0;
RyoheiHagimoto 0:0e0631af0305 511 findNN(root_, result, vec, checks, maxChecks, heap);
RyoheiHagimoto 0:0e0631af0305 512
RyoheiHagimoto 0:0e0631af0305 513 BranchSt branch;
RyoheiHagimoto 0:0e0631af0305 514 while (heap->popMin(branch) && (checks<maxChecks || !result.full())) {
RyoheiHagimoto 0:0e0631af0305 515 KMeansNodePtr node = branch.node;
RyoheiHagimoto 0:0e0631af0305 516 findNN(node, result, vec, checks, maxChecks, heap);
RyoheiHagimoto 0:0e0631af0305 517 }
RyoheiHagimoto 0:0e0631af0305 518 assert(result.full());
RyoheiHagimoto 0:0e0631af0305 519
RyoheiHagimoto 0:0e0631af0305 520 delete heap;
RyoheiHagimoto 0:0e0631af0305 521 }
RyoheiHagimoto 0:0e0631af0305 522
RyoheiHagimoto 0:0e0631af0305 523 }
RyoheiHagimoto 0:0e0631af0305 524
RyoheiHagimoto 0:0e0631af0305 525 /**
RyoheiHagimoto 0:0e0631af0305 526 * Clustering function that takes a cut in the hierarchical k-means
RyoheiHagimoto 0:0e0631af0305 527 * tree and return the clusters centers of that clustering.
RyoheiHagimoto 0:0e0631af0305 528 * Params:
RyoheiHagimoto 0:0e0631af0305 529 * numClusters = number of clusters to have in the clustering computed
RyoheiHagimoto 0:0e0631af0305 530 * Returns: number of cluster centers
RyoheiHagimoto 0:0e0631af0305 531 */
RyoheiHagimoto 0:0e0631af0305 532 int getClusterCenters(Matrix<DistanceType>& centers)
RyoheiHagimoto 0:0e0631af0305 533 {
RyoheiHagimoto 0:0e0631af0305 534 int numClusters = centers.rows;
RyoheiHagimoto 0:0e0631af0305 535 if (numClusters<1) {
RyoheiHagimoto 0:0e0631af0305 536 throw FLANNException("Number of clusters must be at least 1");
RyoheiHagimoto 0:0e0631af0305 537 }
RyoheiHagimoto 0:0e0631af0305 538
RyoheiHagimoto 0:0e0631af0305 539 DistanceType variance;
RyoheiHagimoto 0:0e0631af0305 540 KMeansNodePtr* clusters = new KMeansNodePtr[numClusters];
RyoheiHagimoto 0:0e0631af0305 541
RyoheiHagimoto 0:0e0631af0305 542 int clusterCount = getMinVarianceClusters(root_, clusters, numClusters, variance);
RyoheiHagimoto 0:0e0631af0305 543
RyoheiHagimoto 0:0e0631af0305 544 Logger::info("Clusters requested: %d, returning %d\n",numClusters, clusterCount);
RyoheiHagimoto 0:0e0631af0305 545
RyoheiHagimoto 0:0e0631af0305 546 for (int i=0; i<clusterCount; ++i) {
RyoheiHagimoto 0:0e0631af0305 547 DistanceType* center = clusters[i]->pivot;
RyoheiHagimoto 0:0e0631af0305 548 for (size_t j=0; j<veclen_; ++j) {
RyoheiHagimoto 0:0e0631af0305 549 centers[i][j] = center[j];
RyoheiHagimoto 0:0e0631af0305 550 }
RyoheiHagimoto 0:0e0631af0305 551 }
RyoheiHagimoto 0:0e0631af0305 552 delete[] clusters;
RyoheiHagimoto 0:0e0631af0305 553
RyoheiHagimoto 0:0e0631af0305 554 return clusterCount;
RyoheiHagimoto 0:0e0631af0305 555 }
RyoheiHagimoto 0:0e0631af0305 556
RyoheiHagimoto 0:0e0631af0305 557 IndexParams getParameters() const
RyoheiHagimoto 0:0e0631af0305 558 {
RyoheiHagimoto 0:0e0631af0305 559 return index_params_;
RyoheiHagimoto 0:0e0631af0305 560 }
RyoheiHagimoto 0:0e0631af0305 561
RyoheiHagimoto 0:0e0631af0305 562
RyoheiHagimoto 0:0e0631af0305 563 private:
RyoheiHagimoto 0:0e0631af0305 564 /**
RyoheiHagimoto 0:0e0631af0305 565 * Struture representing a node in the hierarchical k-means tree.
RyoheiHagimoto 0:0e0631af0305 566 */
RyoheiHagimoto 0:0e0631af0305 567 struct KMeansNode
RyoheiHagimoto 0:0e0631af0305 568 {
RyoheiHagimoto 0:0e0631af0305 569 /**
RyoheiHagimoto 0:0e0631af0305 570 * The cluster center.
RyoheiHagimoto 0:0e0631af0305 571 */
RyoheiHagimoto 0:0e0631af0305 572 DistanceType* pivot;
RyoheiHagimoto 0:0e0631af0305 573 /**
RyoheiHagimoto 0:0e0631af0305 574 * The cluster radius.
RyoheiHagimoto 0:0e0631af0305 575 */
RyoheiHagimoto 0:0e0631af0305 576 DistanceType radius;
RyoheiHagimoto 0:0e0631af0305 577 /**
RyoheiHagimoto 0:0e0631af0305 578 * The cluster mean radius.
RyoheiHagimoto 0:0e0631af0305 579 */
RyoheiHagimoto 0:0e0631af0305 580 DistanceType mean_radius;
RyoheiHagimoto 0:0e0631af0305 581 /**
RyoheiHagimoto 0:0e0631af0305 582 * The cluster variance.
RyoheiHagimoto 0:0e0631af0305 583 */
RyoheiHagimoto 0:0e0631af0305 584 DistanceType variance;
RyoheiHagimoto 0:0e0631af0305 585 /**
RyoheiHagimoto 0:0e0631af0305 586 * The cluster size (number of points in the cluster)
RyoheiHagimoto 0:0e0631af0305 587 */
RyoheiHagimoto 0:0e0631af0305 588 int size;
RyoheiHagimoto 0:0e0631af0305 589 /**
RyoheiHagimoto 0:0e0631af0305 590 * Child nodes (only for non-terminal nodes)
RyoheiHagimoto 0:0e0631af0305 591 */
RyoheiHagimoto 0:0e0631af0305 592 KMeansNode** childs;
RyoheiHagimoto 0:0e0631af0305 593 /**
RyoheiHagimoto 0:0e0631af0305 594 * Node points (only for terminal nodes)
RyoheiHagimoto 0:0e0631af0305 595 */
RyoheiHagimoto 0:0e0631af0305 596 int* indices;
RyoheiHagimoto 0:0e0631af0305 597 /**
RyoheiHagimoto 0:0e0631af0305 598 * Level
RyoheiHagimoto 0:0e0631af0305 599 */
RyoheiHagimoto 0:0e0631af0305 600 int level;
RyoheiHagimoto 0:0e0631af0305 601 };
RyoheiHagimoto 0:0e0631af0305 602 typedef KMeansNode* KMeansNodePtr;
RyoheiHagimoto 0:0e0631af0305 603
RyoheiHagimoto 0:0e0631af0305 604 /**
RyoheiHagimoto 0:0e0631af0305 605 * Alias definition for a nicer syntax.
RyoheiHagimoto 0:0e0631af0305 606 */
RyoheiHagimoto 0:0e0631af0305 607 typedef BranchStruct<KMeansNodePtr, DistanceType> BranchSt;
RyoheiHagimoto 0:0e0631af0305 608
RyoheiHagimoto 0:0e0631af0305 609
RyoheiHagimoto 0:0e0631af0305 610
RyoheiHagimoto 0:0e0631af0305 611
RyoheiHagimoto 0:0e0631af0305 612 void save_tree(FILE* stream, KMeansNodePtr node)
RyoheiHagimoto 0:0e0631af0305 613 {
RyoheiHagimoto 0:0e0631af0305 614 save_value(stream, *node);
RyoheiHagimoto 0:0e0631af0305 615 save_value(stream, *(node->pivot), (int)veclen_);
RyoheiHagimoto 0:0e0631af0305 616 if (node->childs==NULL) {
RyoheiHagimoto 0:0e0631af0305 617 int indices_offset = (int)(node->indices - indices_);
RyoheiHagimoto 0:0e0631af0305 618 save_value(stream, indices_offset);
RyoheiHagimoto 0:0e0631af0305 619 }
RyoheiHagimoto 0:0e0631af0305 620 else {
RyoheiHagimoto 0:0e0631af0305 621 for(int i=0; i<branching_; ++i) {
RyoheiHagimoto 0:0e0631af0305 622 save_tree(stream, node->childs[i]);
RyoheiHagimoto 0:0e0631af0305 623 }
RyoheiHagimoto 0:0e0631af0305 624 }
RyoheiHagimoto 0:0e0631af0305 625 }
RyoheiHagimoto 0:0e0631af0305 626
RyoheiHagimoto 0:0e0631af0305 627
RyoheiHagimoto 0:0e0631af0305 628 void load_tree(FILE* stream, KMeansNodePtr& node)
RyoheiHagimoto 0:0e0631af0305 629 {
RyoheiHagimoto 0:0e0631af0305 630 node = pool_.allocate<KMeansNode>();
RyoheiHagimoto 0:0e0631af0305 631 load_value(stream, *node);
RyoheiHagimoto 0:0e0631af0305 632 node->pivot = new DistanceType[veclen_];
RyoheiHagimoto 0:0e0631af0305 633 load_value(stream, *(node->pivot), (int)veclen_);
RyoheiHagimoto 0:0e0631af0305 634 if (node->childs==NULL) {
RyoheiHagimoto 0:0e0631af0305 635 int indices_offset;
RyoheiHagimoto 0:0e0631af0305 636 load_value(stream, indices_offset);
RyoheiHagimoto 0:0e0631af0305 637 node->indices = indices_ + indices_offset;
RyoheiHagimoto 0:0e0631af0305 638 }
RyoheiHagimoto 0:0e0631af0305 639 else {
RyoheiHagimoto 0:0e0631af0305 640 node->childs = pool_.allocate<KMeansNodePtr>(branching_);
RyoheiHagimoto 0:0e0631af0305 641 for(int i=0; i<branching_; ++i) {
RyoheiHagimoto 0:0e0631af0305 642 load_tree(stream, node->childs[i]);
RyoheiHagimoto 0:0e0631af0305 643 }
RyoheiHagimoto 0:0e0631af0305 644 }
RyoheiHagimoto 0:0e0631af0305 645 }
RyoheiHagimoto 0:0e0631af0305 646
RyoheiHagimoto 0:0e0631af0305 647
RyoheiHagimoto 0:0e0631af0305 648 /**
RyoheiHagimoto 0:0e0631af0305 649 * Helper function
RyoheiHagimoto 0:0e0631af0305 650 */
RyoheiHagimoto 0:0e0631af0305 651 void free_centers(KMeansNodePtr node)
RyoheiHagimoto 0:0e0631af0305 652 {
RyoheiHagimoto 0:0e0631af0305 653 delete[] node->pivot;
RyoheiHagimoto 0:0e0631af0305 654 if (node->childs!=NULL) {
RyoheiHagimoto 0:0e0631af0305 655 for (int k=0; k<branching_; ++k) {
RyoheiHagimoto 0:0e0631af0305 656 free_centers(node->childs[k]);
RyoheiHagimoto 0:0e0631af0305 657 }
RyoheiHagimoto 0:0e0631af0305 658 }
RyoheiHagimoto 0:0e0631af0305 659 }
RyoheiHagimoto 0:0e0631af0305 660
RyoheiHagimoto 0:0e0631af0305 661 /**
RyoheiHagimoto 0:0e0631af0305 662 * Computes the statistics of a node (mean, radius, variance).
RyoheiHagimoto 0:0e0631af0305 663 *
RyoheiHagimoto 0:0e0631af0305 664 * Params:
RyoheiHagimoto 0:0e0631af0305 665 * node = the node to use
RyoheiHagimoto 0:0e0631af0305 666 * indices = the indices of the points belonging to the node
RyoheiHagimoto 0:0e0631af0305 667 */
RyoheiHagimoto 0:0e0631af0305 668 void computeNodeStatistics(KMeansNodePtr node, int* indices, int indices_length)
RyoheiHagimoto 0:0e0631af0305 669 {
RyoheiHagimoto 0:0e0631af0305 670
RyoheiHagimoto 0:0e0631af0305 671 DistanceType radius = 0;
RyoheiHagimoto 0:0e0631af0305 672 DistanceType variance = 0;
RyoheiHagimoto 0:0e0631af0305 673 DistanceType* mean = new DistanceType[veclen_];
RyoheiHagimoto 0:0e0631af0305 674 memoryCounter_ += int(veclen_*sizeof(DistanceType));
RyoheiHagimoto 0:0e0631af0305 675
RyoheiHagimoto 0:0e0631af0305 676 memset(mean,0,veclen_*sizeof(DistanceType));
RyoheiHagimoto 0:0e0631af0305 677
RyoheiHagimoto 0:0e0631af0305 678 for (size_t i=0; i<size_; ++i) {
RyoheiHagimoto 0:0e0631af0305 679 ElementType* vec = dataset_[indices[i]];
RyoheiHagimoto 0:0e0631af0305 680 for (size_t j=0; j<veclen_; ++j) {
RyoheiHagimoto 0:0e0631af0305 681 mean[j] += vec[j];
RyoheiHagimoto 0:0e0631af0305 682 }
RyoheiHagimoto 0:0e0631af0305 683 variance += distance_(vec, ZeroIterator<ElementType>(), veclen_);
RyoheiHagimoto 0:0e0631af0305 684 }
RyoheiHagimoto 0:0e0631af0305 685 for (size_t j=0; j<veclen_; ++j) {
RyoheiHagimoto 0:0e0631af0305 686 mean[j] /= size_;
RyoheiHagimoto 0:0e0631af0305 687 }
RyoheiHagimoto 0:0e0631af0305 688 variance /= size_;
RyoheiHagimoto 0:0e0631af0305 689 variance -= distance_(mean, ZeroIterator<ElementType>(), veclen_);
RyoheiHagimoto 0:0e0631af0305 690
RyoheiHagimoto 0:0e0631af0305 691 DistanceType tmp = 0;
RyoheiHagimoto 0:0e0631af0305 692 for (int i=0; i<indices_length; ++i) {
RyoheiHagimoto 0:0e0631af0305 693 tmp = distance_(mean, dataset_[indices[i]], veclen_);
RyoheiHagimoto 0:0e0631af0305 694 if (tmp>radius) {
RyoheiHagimoto 0:0e0631af0305 695 radius = tmp;
RyoheiHagimoto 0:0e0631af0305 696 }
RyoheiHagimoto 0:0e0631af0305 697 }
RyoheiHagimoto 0:0e0631af0305 698
RyoheiHagimoto 0:0e0631af0305 699 node->variance = variance;
RyoheiHagimoto 0:0e0631af0305 700 node->radius = radius;
RyoheiHagimoto 0:0e0631af0305 701 node->pivot = mean;
RyoheiHagimoto 0:0e0631af0305 702 }
RyoheiHagimoto 0:0e0631af0305 703
RyoheiHagimoto 0:0e0631af0305 704
RyoheiHagimoto 0:0e0631af0305 705 /**
RyoheiHagimoto 0:0e0631af0305 706 * The method responsible with actually doing the recursive hierarchical
RyoheiHagimoto 0:0e0631af0305 707 * clustering
RyoheiHagimoto 0:0e0631af0305 708 *
RyoheiHagimoto 0:0e0631af0305 709 * Params:
RyoheiHagimoto 0:0e0631af0305 710 * node = the node to cluster
RyoheiHagimoto 0:0e0631af0305 711 * indices = indices of the points belonging to the current node
RyoheiHagimoto 0:0e0631af0305 712 * branching = the branching factor to use in the clustering
RyoheiHagimoto 0:0e0631af0305 713 *
RyoheiHagimoto 0:0e0631af0305 714 * TODO: for 1-sized clusters don't store a cluster center (it's the same as the single cluster point)
RyoheiHagimoto 0:0e0631af0305 715 */
RyoheiHagimoto 0:0e0631af0305 716 void computeClustering(KMeansNodePtr node, int* indices, int indices_length, int branching, int level)
RyoheiHagimoto 0:0e0631af0305 717 {
RyoheiHagimoto 0:0e0631af0305 718 node->size = indices_length;
RyoheiHagimoto 0:0e0631af0305 719 node->level = level;
RyoheiHagimoto 0:0e0631af0305 720
RyoheiHagimoto 0:0e0631af0305 721 if (indices_length < branching) {
RyoheiHagimoto 0:0e0631af0305 722 node->indices = indices;
RyoheiHagimoto 0:0e0631af0305 723 std::sort(node->indices,node->indices+indices_length);
RyoheiHagimoto 0:0e0631af0305 724 node->childs = NULL;
RyoheiHagimoto 0:0e0631af0305 725 return;
RyoheiHagimoto 0:0e0631af0305 726 }
RyoheiHagimoto 0:0e0631af0305 727
RyoheiHagimoto 0:0e0631af0305 728 cv::AutoBuffer<int> centers_idx_buf(branching);
RyoheiHagimoto 0:0e0631af0305 729 int* centers_idx = (int*)centers_idx_buf;
RyoheiHagimoto 0:0e0631af0305 730 int centers_length;
RyoheiHagimoto 0:0e0631af0305 731 (this->*chooseCenters)(branching, indices, indices_length, centers_idx, centers_length);
RyoheiHagimoto 0:0e0631af0305 732
RyoheiHagimoto 0:0e0631af0305 733 if (centers_length<branching) {
RyoheiHagimoto 0:0e0631af0305 734 node->indices = indices;
RyoheiHagimoto 0:0e0631af0305 735 std::sort(node->indices,node->indices+indices_length);
RyoheiHagimoto 0:0e0631af0305 736 node->childs = NULL;
RyoheiHagimoto 0:0e0631af0305 737 return;
RyoheiHagimoto 0:0e0631af0305 738 }
RyoheiHagimoto 0:0e0631af0305 739
RyoheiHagimoto 0:0e0631af0305 740
RyoheiHagimoto 0:0e0631af0305 741 cv::AutoBuffer<double> dcenters_buf(branching*veclen_);
RyoheiHagimoto 0:0e0631af0305 742 Matrix<double> dcenters((double*)dcenters_buf,branching,veclen_);
RyoheiHagimoto 0:0e0631af0305 743 for (int i=0; i<centers_length; ++i) {
RyoheiHagimoto 0:0e0631af0305 744 ElementType* vec = dataset_[centers_idx[i]];
RyoheiHagimoto 0:0e0631af0305 745 for (size_t k=0; k<veclen_; ++k) {
RyoheiHagimoto 0:0e0631af0305 746 dcenters[i][k] = double(vec[k]);
RyoheiHagimoto 0:0e0631af0305 747 }
RyoheiHagimoto 0:0e0631af0305 748 }
RyoheiHagimoto 0:0e0631af0305 749
RyoheiHagimoto 0:0e0631af0305 750 std::vector<DistanceType> radiuses(branching);
RyoheiHagimoto 0:0e0631af0305 751 cv::AutoBuffer<int> count_buf(branching);
RyoheiHagimoto 0:0e0631af0305 752 int* count = (int*)count_buf;
RyoheiHagimoto 0:0e0631af0305 753 for (int i=0; i<branching; ++i) {
RyoheiHagimoto 0:0e0631af0305 754 radiuses[i] = 0;
RyoheiHagimoto 0:0e0631af0305 755 count[i] = 0;
RyoheiHagimoto 0:0e0631af0305 756 }
RyoheiHagimoto 0:0e0631af0305 757
RyoheiHagimoto 0:0e0631af0305 758 // assign points to clusters
RyoheiHagimoto 0:0e0631af0305 759 cv::AutoBuffer<int> belongs_to_buf(indices_length);
RyoheiHagimoto 0:0e0631af0305 760 int* belongs_to = (int*)belongs_to_buf;
RyoheiHagimoto 0:0e0631af0305 761 for (int i=0; i<indices_length; ++i) {
RyoheiHagimoto 0:0e0631af0305 762
RyoheiHagimoto 0:0e0631af0305 763 DistanceType sq_dist = distance_(dataset_[indices[i]], dcenters[0], veclen_);
RyoheiHagimoto 0:0e0631af0305 764 belongs_to[i] = 0;
RyoheiHagimoto 0:0e0631af0305 765 for (int j=1; j<branching; ++j) {
RyoheiHagimoto 0:0e0631af0305 766 DistanceType new_sq_dist = distance_(dataset_[indices[i]], dcenters[j], veclen_);
RyoheiHagimoto 0:0e0631af0305 767 if (sq_dist>new_sq_dist) {
RyoheiHagimoto 0:0e0631af0305 768 belongs_to[i] = j;
RyoheiHagimoto 0:0e0631af0305 769 sq_dist = new_sq_dist;
RyoheiHagimoto 0:0e0631af0305 770 }
RyoheiHagimoto 0:0e0631af0305 771 }
RyoheiHagimoto 0:0e0631af0305 772 if (sq_dist>radiuses[belongs_to[i]]) {
RyoheiHagimoto 0:0e0631af0305 773 radiuses[belongs_to[i]] = sq_dist;
RyoheiHagimoto 0:0e0631af0305 774 }
RyoheiHagimoto 0:0e0631af0305 775 count[belongs_to[i]]++;
RyoheiHagimoto 0:0e0631af0305 776 }
RyoheiHagimoto 0:0e0631af0305 777
RyoheiHagimoto 0:0e0631af0305 778 bool converged = false;
RyoheiHagimoto 0:0e0631af0305 779 int iteration = 0;
RyoheiHagimoto 0:0e0631af0305 780 while (!converged && iteration<iterations_) {
RyoheiHagimoto 0:0e0631af0305 781 converged = true;
RyoheiHagimoto 0:0e0631af0305 782 iteration++;
RyoheiHagimoto 0:0e0631af0305 783
RyoheiHagimoto 0:0e0631af0305 784 // compute the new cluster centers
RyoheiHagimoto 0:0e0631af0305 785 for (int i=0; i<branching; ++i) {
RyoheiHagimoto 0:0e0631af0305 786 memset(dcenters[i],0,sizeof(double)*veclen_);
RyoheiHagimoto 0:0e0631af0305 787 radiuses[i] = 0;
RyoheiHagimoto 0:0e0631af0305 788 }
RyoheiHagimoto 0:0e0631af0305 789 for (int i=0; i<indices_length; ++i) {
RyoheiHagimoto 0:0e0631af0305 790 ElementType* vec = dataset_[indices[i]];
RyoheiHagimoto 0:0e0631af0305 791 double* center = dcenters[belongs_to[i]];
RyoheiHagimoto 0:0e0631af0305 792 for (size_t k=0; k<veclen_; ++k) {
RyoheiHagimoto 0:0e0631af0305 793 center[k] += vec[k];
RyoheiHagimoto 0:0e0631af0305 794 }
RyoheiHagimoto 0:0e0631af0305 795 }
RyoheiHagimoto 0:0e0631af0305 796 for (int i=0; i<branching; ++i) {
RyoheiHagimoto 0:0e0631af0305 797 int cnt = count[i];
RyoheiHagimoto 0:0e0631af0305 798 for (size_t k=0; k<veclen_; ++k) {
RyoheiHagimoto 0:0e0631af0305 799 dcenters[i][k] /= cnt;
RyoheiHagimoto 0:0e0631af0305 800 }
RyoheiHagimoto 0:0e0631af0305 801 }
RyoheiHagimoto 0:0e0631af0305 802
RyoheiHagimoto 0:0e0631af0305 803 // reassign points to clusters
RyoheiHagimoto 0:0e0631af0305 804 cv::Mutex mtx;
RyoheiHagimoto 0:0e0631af0305 805 KMeansDistanceComputer invoker(distance_, dataset_, branching, indices, dcenters, veclen_, count, belongs_to, radiuses, converged, mtx);
RyoheiHagimoto 0:0e0631af0305 806 parallel_for_(cv::Range(0, (int)indices_length), invoker);
RyoheiHagimoto 0:0e0631af0305 807
RyoheiHagimoto 0:0e0631af0305 808 for (int i=0; i<branching; ++i) {
RyoheiHagimoto 0:0e0631af0305 809 // if one cluster converges to an empty cluster,
RyoheiHagimoto 0:0e0631af0305 810 // move an element into that cluster
RyoheiHagimoto 0:0e0631af0305 811 if (count[i]==0) {
RyoheiHagimoto 0:0e0631af0305 812 int j = (i+1)%branching;
RyoheiHagimoto 0:0e0631af0305 813 while (count[j]<=1) {
RyoheiHagimoto 0:0e0631af0305 814 j = (j+1)%branching;
RyoheiHagimoto 0:0e0631af0305 815 }
RyoheiHagimoto 0:0e0631af0305 816
RyoheiHagimoto 0:0e0631af0305 817 for (int k=0; k<indices_length; ++k) {
RyoheiHagimoto 0:0e0631af0305 818 if (belongs_to[k]==j) {
RyoheiHagimoto 0:0e0631af0305 819 // for cluster j, we move the furthest element from the center to the empty cluster i
RyoheiHagimoto 0:0e0631af0305 820 if ( distance_(dataset_[indices[k]], dcenters[j], veclen_) == radiuses[j] ) {
RyoheiHagimoto 0:0e0631af0305 821 belongs_to[k] = i;
RyoheiHagimoto 0:0e0631af0305 822 count[j]--;
RyoheiHagimoto 0:0e0631af0305 823 count[i]++;
RyoheiHagimoto 0:0e0631af0305 824 break;
RyoheiHagimoto 0:0e0631af0305 825 }
RyoheiHagimoto 0:0e0631af0305 826 }
RyoheiHagimoto 0:0e0631af0305 827 }
RyoheiHagimoto 0:0e0631af0305 828 converged = false;
RyoheiHagimoto 0:0e0631af0305 829 }
RyoheiHagimoto 0:0e0631af0305 830 }
RyoheiHagimoto 0:0e0631af0305 831
RyoheiHagimoto 0:0e0631af0305 832 }
RyoheiHagimoto 0:0e0631af0305 833
RyoheiHagimoto 0:0e0631af0305 834 DistanceType** centers = new DistanceType*[branching];
RyoheiHagimoto 0:0e0631af0305 835
RyoheiHagimoto 0:0e0631af0305 836 for (int i=0; i<branching; ++i) {
RyoheiHagimoto 0:0e0631af0305 837 centers[i] = new DistanceType[veclen_];
RyoheiHagimoto 0:0e0631af0305 838 memoryCounter_ += (int)(veclen_*sizeof(DistanceType));
RyoheiHagimoto 0:0e0631af0305 839 for (size_t k=0; k<veclen_; ++k) {
RyoheiHagimoto 0:0e0631af0305 840 centers[i][k] = (DistanceType)dcenters[i][k];
RyoheiHagimoto 0:0e0631af0305 841 }
RyoheiHagimoto 0:0e0631af0305 842 }
RyoheiHagimoto 0:0e0631af0305 843
RyoheiHagimoto 0:0e0631af0305 844
RyoheiHagimoto 0:0e0631af0305 845 // compute kmeans clustering for each of the resulting clusters
RyoheiHagimoto 0:0e0631af0305 846 node->childs = pool_.allocate<KMeansNodePtr>(branching);
RyoheiHagimoto 0:0e0631af0305 847 int start = 0;
RyoheiHagimoto 0:0e0631af0305 848 int end = start;
RyoheiHagimoto 0:0e0631af0305 849 for (int c=0; c<branching; ++c) {
RyoheiHagimoto 0:0e0631af0305 850 int s = count[c];
RyoheiHagimoto 0:0e0631af0305 851
RyoheiHagimoto 0:0e0631af0305 852 DistanceType variance = 0;
RyoheiHagimoto 0:0e0631af0305 853 DistanceType mean_radius =0;
RyoheiHagimoto 0:0e0631af0305 854 for (int i=0; i<indices_length; ++i) {
RyoheiHagimoto 0:0e0631af0305 855 if (belongs_to[i]==c) {
RyoheiHagimoto 0:0e0631af0305 856 DistanceType d = distance_(dataset_[indices[i]], ZeroIterator<ElementType>(), veclen_);
RyoheiHagimoto 0:0e0631af0305 857 variance += d;
RyoheiHagimoto 0:0e0631af0305 858 mean_radius += sqrt(d);
RyoheiHagimoto 0:0e0631af0305 859 std::swap(indices[i],indices[end]);
RyoheiHagimoto 0:0e0631af0305 860 std::swap(belongs_to[i],belongs_to[end]);
RyoheiHagimoto 0:0e0631af0305 861 end++;
RyoheiHagimoto 0:0e0631af0305 862 }
RyoheiHagimoto 0:0e0631af0305 863 }
RyoheiHagimoto 0:0e0631af0305 864 variance /= s;
RyoheiHagimoto 0:0e0631af0305 865 mean_radius /= s;
RyoheiHagimoto 0:0e0631af0305 866 variance -= distance_(centers[c], ZeroIterator<ElementType>(), veclen_);
RyoheiHagimoto 0:0e0631af0305 867
RyoheiHagimoto 0:0e0631af0305 868 node->childs[c] = pool_.allocate<KMeansNode>();
RyoheiHagimoto 0:0e0631af0305 869 std::memset(node->childs[c], 0, sizeof(KMeansNode));
RyoheiHagimoto 0:0e0631af0305 870 node->childs[c]->radius = radiuses[c];
RyoheiHagimoto 0:0e0631af0305 871 node->childs[c]->pivot = centers[c];
RyoheiHagimoto 0:0e0631af0305 872 node->childs[c]->variance = variance;
RyoheiHagimoto 0:0e0631af0305 873 node->childs[c]->mean_radius = mean_radius;
RyoheiHagimoto 0:0e0631af0305 874 computeClustering(node->childs[c],indices+start, end-start, branching, level+1);
RyoheiHagimoto 0:0e0631af0305 875 start=end;
RyoheiHagimoto 0:0e0631af0305 876 }
RyoheiHagimoto 0:0e0631af0305 877
RyoheiHagimoto 0:0e0631af0305 878 delete[] centers;
RyoheiHagimoto 0:0e0631af0305 879 }
RyoheiHagimoto 0:0e0631af0305 880
RyoheiHagimoto 0:0e0631af0305 881
RyoheiHagimoto 0:0e0631af0305 882
RyoheiHagimoto 0:0e0631af0305 883 /**
RyoheiHagimoto 0:0e0631af0305 884 * Performs one descent in the hierarchical k-means tree. The branches not
RyoheiHagimoto 0:0e0631af0305 885 * visited are stored in a priority queue.
RyoheiHagimoto 0:0e0631af0305 886 *
RyoheiHagimoto 0:0e0631af0305 887 * Params:
RyoheiHagimoto 0:0e0631af0305 888 * node = node to explore
RyoheiHagimoto 0:0e0631af0305 889 * result = container for the k-nearest neighbors found
RyoheiHagimoto 0:0e0631af0305 890 * vec = query points
RyoheiHagimoto 0:0e0631af0305 891 * checks = how many points in the dataset have been checked so far
RyoheiHagimoto 0:0e0631af0305 892 * maxChecks = maximum dataset points to checks
RyoheiHagimoto 0:0e0631af0305 893 */
RyoheiHagimoto 0:0e0631af0305 894
RyoheiHagimoto 0:0e0631af0305 895
RyoheiHagimoto 0:0e0631af0305 896 void findNN(KMeansNodePtr node, ResultSet<DistanceType>& result, const ElementType* vec, int& checks, int maxChecks,
RyoheiHagimoto 0:0e0631af0305 897 Heap<BranchSt>* heap)
RyoheiHagimoto 0:0e0631af0305 898 {
RyoheiHagimoto 0:0e0631af0305 899 // Ignore those clusters that are too far away
RyoheiHagimoto 0:0e0631af0305 900 {
RyoheiHagimoto 0:0e0631af0305 901 DistanceType bsq = distance_(vec, node->pivot, veclen_);
RyoheiHagimoto 0:0e0631af0305 902 DistanceType rsq = node->radius;
RyoheiHagimoto 0:0e0631af0305 903 DistanceType wsq = result.worstDist();
RyoheiHagimoto 0:0e0631af0305 904
RyoheiHagimoto 0:0e0631af0305 905 DistanceType val = bsq-rsq-wsq;
RyoheiHagimoto 0:0e0631af0305 906 DistanceType val2 = val*val-4*rsq*wsq;
RyoheiHagimoto 0:0e0631af0305 907
RyoheiHagimoto 0:0e0631af0305 908 //if (val>0) {
RyoheiHagimoto 0:0e0631af0305 909 if ((val>0)&&(val2>0)) {
RyoheiHagimoto 0:0e0631af0305 910 return;
RyoheiHagimoto 0:0e0631af0305 911 }
RyoheiHagimoto 0:0e0631af0305 912 }
RyoheiHagimoto 0:0e0631af0305 913
RyoheiHagimoto 0:0e0631af0305 914 if (node->childs==NULL) {
RyoheiHagimoto 0:0e0631af0305 915 if (checks>=maxChecks) {
RyoheiHagimoto 0:0e0631af0305 916 if (result.full()) return;
RyoheiHagimoto 0:0e0631af0305 917 }
RyoheiHagimoto 0:0e0631af0305 918 checks += node->size;
RyoheiHagimoto 0:0e0631af0305 919 for (int i=0; i<node->size; ++i) {
RyoheiHagimoto 0:0e0631af0305 920 int index = node->indices[i];
RyoheiHagimoto 0:0e0631af0305 921 DistanceType dist = distance_(dataset_[index], vec, veclen_);
RyoheiHagimoto 0:0e0631af0305 922 result.addPoint(dist, index);
RyoheiHagimoto 0:0e0631af0305 923 }
RyoheiHagimoto 0:0e0631af0305 924 }
RyoheiHagimoto 0:0e0631af0305 925 else {
RyoheiHagimoto 0:0e0631af0305 926 DistanceType* domain_distances = new DistanceType[branching_];
RyoheiHagimoto 0:0e0631af0305 927 int closest_center = exploreNodeBranches(node, vec, domain_distances, heap);
RyoheiHagimoto 0:0e0631af0305 928 delete[] domain_distances;
RyoheiHagimoto 0:0e0631af0305 929 findNN(node->childs[closest_center],result,vec, checks, maxChecks, heap);
RyoheiHagimoto 0:0e0631af0305 930 }
RyoheiHagimoto 0:0e0631af0305 931 }
RyoheiHagimoto 0:0e0631af0305 932
RyoheiHagimoto 0:0e0631af0305 933 /**
RyoheiHagimoto 0:0e0631af0305 934 * Helper function that computes the nearest childs of a node to a given query point.
RyoheiHagimoto 0:0e0631af0305 935 * Params:
RyoheiHagimoto 0:0e0631af0305 936 * node = the node
RyoheiHagimoto 0:0e0631af0305 937 * q = the query point
RyoheiHagimoto 0:0e0631af0305 938 * distances = array with the distances to each child node.
RyoheiHagimoto 0:0e0631af0305 939 * Returns:
RyoheiHagimoto 0:0e0631af0305 940 */
RyoheiHagimoto 0:0e0631af0305 941 int exploreNodeBranches(KMeansNodePtr node, const ElementType* q, DistanceType* domain_distances, Heap<BranchSt>* heap)
RyoheiHagimoto 0:0e0631af0305 942 {
RyoheiHagimoto 0:0e0631af0305 943
RyoheiHagimoto 0:0e0631af0305 944 int best_index = 0;
RyoheiHagimoto 0:0e0631af0305 945 domain_distances[best_index] = distance_(q, node->childs[best_index]->pivot, veclen_);
RyoheiHagimoto 0:0e0631af0305 946 for (int i=1; i<branching_; ++i) {
RyoheiHagimoto 0:0e0631af0305 947 domain_distances[i] = distance_(q, node->childs[i]->pivot, veclen_);
RyoheiHagimoto 0:0e0631af0305 948 if (domain_distances[i]<domain_distances[best_index]) {
RyoheiHagimoto 0:0e0631af0305 949 best_index = i;
RyoheiHagimoto 0:0e0631af0305 950 }
RyoheiHagimoto 0:0e0631af0305 951 }
RyoheiHagimoto 0:0e0631af0305 952
RyoheiHagimoto 0:0e0631af0305 953 // float* best_center = node->childs[best_index]->pivot;
RyoheiHagimoto 0:0e0631af0305 954 for (int i=0; i<branching_; ++i) {
RyoheiHagimoto 0:0e0631af0305 955 if (i != best_index) {
RyoheiHagimoto 0:0e0631af0305 956 domain_distances[i] -= cb_index_*node->childs[i]->variance;
RyoheiHagimoto 0:0e0631af0305 957
RyoheiHagimoto 0:0e0631af0305 958 // float dist_to_border = getDistanceToBorder(node.childs[i].pivot,best_center,q);
RyoheiHagimoto 0:0e0631af0305 959 // if (domain_distances[i]<dist_to_border) {
RyoheiHagimoto 0:0e0631af0305 960 // domain_distances[i] = dist_to_border;
RyoheiHagimoto 0:0e0631af0305 961 // }
RyoheiHagimoto 0:0e0631af0305 962 heap->insert(BranchSt(node->childs[i],domain_distances[i]));
RyoheiHagimoto 0:0e0631af0305 963 }
RyoheiHagimoto 0:0e0631af0305 964 }
RyoheiHagimoto 0:0e0631af0305 965
RyoheiHagimoto 0:0e0631af0305 966 return best_index;
RyoheiHagimoto 0:0e0631af0305 967 }
RyoheiHagimoto 0:0e0631af0305 968
RyoheiHagimoto 0:0e0631af0305 969
RyoheiHagimoto 0:0e0631af0305 970 /**
RyoheiHagimoto 0:0e0631af0305 971 * Function the performs exact nearest neighbor search by traversing the entire tree.
RyoheiHagimoto 0:0e0631af0305 972 */
RyoheiHagimoto 0:0e0631af0305 973 void findExactNN(KMeansNodePtr node, ResultSet<DistanceType>& result, const ElementType* vec)
RyoheiHagimoto 0:0e0631af0305 974 {
RyoheiHagimoto 0:0e0631af0305 975 // Ignore those clusters that are too far away
RyoheiHagimoto 0:0e0631af0305 976 {
RyoheiHagimoto 0:0e0631af0305 977 DistanceType bsq = distance_(vec, node->pivot, veclen_);
RyoheiHagimoto 0:0e0631af0305 978 DistanceType rsq = node->radius;
RyoheiHagimoto 0:0e0631af0305 979 DistanceType wsq = result.worstDist();
RyoheiHagimoto 0:0e0631af0305 980
RyoheiHagimoto 0:0e0631af0305 981 DistanceType val = bsq-rsq-wsq;
RyoheiHagimoto 0:0e0631af0305 982 DistanceType val2 = val*val-4*rsq*wsq;
RyoheiHagimoto 0:0e0631af0305 983
RyoheiHagimoto 0:0e0631af0305 984 // if (val>0) {
RyoheiHagimoto 0:0e0631af0305 985 if ((val>0)&&(val2>0)) {
RyoheiHagimoto 0:0e0631af0305 986 return;
RyoheiHagimoto 0:0e0631af0305 987 }
RyoheiHagimoto 0:0e0631af0305 988 }
RyoheiHagimoto 0:0e0631af0305 989
RyoheiHagimoto 0:0e0631af0305 990
RyoheiHagimoto 0:0e0631af0305 991 if (node->childs==NULL) {
RyoheiHagimoto 0:0e0631af0305 992 for (int i=0; i<node->size; ++i) {
RyoheiHagimoto 0:0e0631af0305 993 int index = node->indices[i];
RyoheiHagimoto 0:0e0631af0305 994 DistanceType dist = distance_(dataset_[index], vec, veclen_);
RyoheiHagimoto 0:0e0631af0305 995 result.addPoint(dist, index);
RyoheiHagimoto 0:0e0631af0305 996 }
RyoheiHagimoto 0:0e0631af0305 997 }
RyoheiHagimoto 0:0e0631af0305 998 else {
RyoheiHagimoto 0:0e0631af0305 999 int* sort_indices = new int[branching_];
RyoheiHagimoto 0:0e0631af0305 1000
RyoheiHagimoto 0:0e0631af0305 1001 getCenterOrdering(node, vec, sort_indices);
RyoheiHagimoto 0:0e0631af0305 1002
RyoheiHagimoto 0:0e0631af0305 1003 for (int i=0; i<branching_; ++i) {
RyoheiHagimoto 0:0e0631af0305 1004 findExactNN(node->childs[sort_indices[i]],result,vec);
RyoheiHagimoto 0:0e0631af0305 1005 }
RyoheiHagimoto 0:0e0631af0305 1006
RyoheiHagimoto 0:0e0631af0305 1007 delete[] sort_indices;
RyoheiHagimoto 0:0e0631af0305 1008 }
RyoheiHagimoto 0:0e0631af0305 1009 }
RyoheiHagimoto 0:0e0631af0305 1010
RyoheiHagimoto 0:0e0631af0305 1011
RyoheiHagimoto 0:0e0631af0305 1012 /**
RyoheiHagimoto 0:0e0631af0305 1013 * Helper function.
RyoheiHagimoto 0:0e0631af0305 1014 *
RyoheiHagimoto 0:0e0631af0305 1015 * I computes the order in which to traverse the child nodes of a particular node.
RyoheiHagimoto 0:0e0631af0305 1016 */
RyoheiHagimoto 0:0e0631af0305 1017 void getCenterOrdering(KMeansNodePtr node, const ElementType* q, int* sort_indices)
RyoheiHagimoto 0:0e0631af0305 1018 {
RyoheiHagimoto 0:0e0631af0305 1019 DistanceType* domain_distances = new DistanceType[branching_];
RyoheiHagimoto 0:0e0631af0305 1020 for (int i=0; i<branching_; ++i) {
RyoheiHagimoto 0:0e0631af0305 1021 DistanceType dist = distance_(q, node->childs[i]->pivot, veclen_);
RyoheiHagimoto 0:0e0631af0305 1022
RyoheiHagimoto 0:0e0631af0305 1023 int j=0;
RyoheiHagimoto 0:0e0631af0305 1024 while (domain_distances[j]<dist && j<i) j++;
RyoheiHagimoto 0:0e0631af0305 1025 for (int k=i; k>j; --k) {
RyoheiHagimoto 0:0e0631af0305 1026 domain_distances[k] = domain_distances[k-1];
RyoheiHagimoto 0:0e0631af0305 1027 sort_indices[k] = sort_indices[k-1];
RyoheiHagimoto 0:0e0631af0305 1028 }
RyoheiHagimoto 0:0e0631af0305 1029 domain_distances[j] = dist;
RyoheiHagimoto 0:0e0631af0305 1030 sort_indices[j] = i;
RyoheiHagimoto 0:0e0631af0305 1031 }
RyoheiHagimoto 0:0e0631af0305 1032 delete[] domain_distances;
RyoheiHagimoto 0:0e0631af0305 1033 }
RyoheiHagimoto 0:0e0631af0305 1034
RyoheiHagimoto 0:0e0631af0305 1035 /**
RyoheiHagimoto 0:0e0631af0305 1036 * Method that computes the squared distance from the query point q
RyoheiHagimoto 0:0e0631af0305 1037 * from inside region with center c to the border between this
RyoheiHagimoto 0:0e0631af0305 1038 * region and the region with center p
RyoheiHagimoto 0:0e0631af0305 1039 */
RyoheiHagimoto 0:0e0631af0305 1040 DistanceType getDistanceToBorder(DistanceType* p, DistanceType* c, DistanceType* q)
RyoheiHagimoto 0:0e0631af0305 1041 {
RyoheiHagimoto 0:0e0631af0305 1042 DistanceType sum = 0;
RyoheiHagimoto 0:0e0631af0305 1043 DistanceType sum2 = 0;
RyoheiHagimoto 0:0e0631af0305 1044
RyoheiHagimoto 0:0e0631af0305 1045 for (int i=0; i<veclen_; ++i) {
RyoheiHagimoto 0:0e0631af0305 1046 DistanceType t = c[i]-p[i];
RyoheiHagimoto 0:0e0631af0305 1047 sum += t*(q[i]-(c[i]+p[i])/2);
RyoheiHagimoto 0:0e0631af0305 1048 sum2 += t*t;
RyoheiHagimoto 0:0e0631af0305 1049 }
RyoheiHagimoto 0:0e0631af0305 1050
RyoheiHagimoto 0:0e0631af0305 1051 return sum*sum/sum2;
RyoheiHagimoto 0:0e0631af0305 1052 }
RyoheiHagimoto 0:0e0631af0305 1053
RyoheiHagimoto 0:0e0631af0305 1054
RyoheiHagimoto 0:0e0631af0305 1055 /**
RyoheiHagimoto 0:0e0631af0305 1056 * Helper function the descends in the hierarchical k-means tree by spliting those clusters that minimize
RyoheiHagimoto 0:0e0631af0305 1057 * the overall variance of the clustering.
RyoheiHagimoto 0:0e0631af0305 1058 * Params:
RyoheiHagimoto 0:0e0631af0305 1059 * root = root node
RyoheiHagimoto 0:0e0631af0305 1060 * clusters = array with clusters centers (return value)
RyoheiHagimoto 0:0e0631af0305 1061 * varianceValue = variance of the clustering (return value)
RyoheiHagimoto 0:0e0631af0305 1062 * Returns:
RyoheiHagimoto 0:0e0631af0305 1063 */
RyoheiHagimoto 0:0e0631af0305 1064 int getMinVarianceClusters(KMeansNodePtr root, KMeansNodePtr* clusters, int clusters_length, DistanceType& varianceValue)
RyoheiHagimoto 0:0e0631af0305 1065 {
RyoheiHagimoto 0:0e0631af0305 1066 int clusterCount = 1;
RyoheiHagimoto 0:0e0631af0305 1067 clusters[0] = root;
RyoheiHagimoto 0:0e0631af0305 1068
RyoheiHagimoto 0:0e0631af0305 1069 DistanceType meanVariance = root->variance*root->size;
RyoheiHagimoto 0:0e0631af0305 1070
RyoheiHagimoto 0:0e0631af0305 1071 while (clusterCount<clusters_length) {
RyoheiHagimoto 0:0e0631af0305 1072 DistanceType minVariance = (std::numeric_limits<DistanceType>::max)();
RyoheiHagimoto 0:0e0631af0305 1073 int splitIndex = -1;
RyoheiHagimoto 0:0e0631af0305 1074
RyoheiHagimoto 0:0e0631af0305 1075 for (int i=0; i<clusterCount; ++i) {
RyoheiHagimoto 0:0e0631af0305 1076 if (clusters[i]->childs != NULL) {
RyoheiHagimoto 0:0e0631af0305 1077
RyoheiHagimoto 0:0e0631af0305 1078 DistanceType variance = meanVariance - clusters[i]->variance*clusters[i]->size;
RyoheiHagimoto 0:0e0631af0305 1079
RyoheiHagimoto 0:0e0631af0305 1080 for (int j=0; j<branching_; ++j) {
RyoheiHagimoto 0:0e0631af0305 1081 variance += clusters[i]->childs[j]->variance*clusters[i]->childs[j]->size;
RyoheiHagimoto 0:0e0631af0305 1082 }
RyoheiHagimoto 0:0e0631af0305 1083 if (variance<minVariance) {
RyoheiHagimoto 0:0e0631af0305 1084 minVariance = variance;
RyoheiHagimoto 0:0e0631af0305 1085 splitIndex = i;
RyoheiHagimoto 0:0e0631af0305 1086 }
RyoheiHagimoto 0:0e0631af0305 1087 }
RyoheiHagimoto 0:0e0631af0305 1088 }
RyoheiHagimoto 0:0e0631af0305 1089
RyoheiHagimoto 0:0e0631af0305 1090 if (splitIndex==-1) break;
RyoheiHagimoto 0:0e0631af0305 1091 if ( (branching_+clusterCount-1) > clusters_length) break;
RyoheiHagimoto 0:0e0631af0305 1092
RyoheiHagimoto 0:0e0631af0305 1093 meanVariance = minVariance;
RyoheiHagimoto 0:0e0631af0305 1094
RyoheiHagimoto 0:0e0631af0305 1095 // split node
RyoheiHagimoto 0:0e0631af0305 1096 KMeansNodePtr toSplit = clusters[splitIndex];
RyoheiHagimoto 0:0e0631af0305 1097 clusters[splitIndex] = toSplit->childs[0];
RyoheiHagimoto 0:0e0631af0305 1098 for (int i=1; i<branching_; ++i) {
RyoheiHagimoto 0:0e0631af0305 1099 clusters[clusterCount++] = toSplit->childs[i];
RyoheiHagimoto 0:0e0631af0305 1100 }
RyoheiHagimoto 0:0e0631af0305 1101 }
RyoheiHagimoto 0:0e0631af0305 1102
RyoheiHagimoto 0:0e0631af0305 1103 varianceValue = meanVariance/root->size;
RyoheiHagimoto 0:0e0631af0305 1104 return clusterCount;
RyoheiHagimoto 0:0e0631af0305 1105 }
RyoheiHagimoto 0:0e0631af0305 1106
RyoheiHagimoto 0:0e0631af0305 1107 private:
RyoheiHagimoto 0:0e0631af0305 1108 /** The branching factor used in the hierarchical k-means clustering */
RyoheiHagimoto 0:0e0631af0305 1109 int branching_;
RyoheiHagimoto 0:0e0631af0305 1110
RyoheiHagimoto 0:0e0631af0305 1111 /** Maximum number of iterations to use when performing k-means clustering */
RyoheiHagimoto 0:0e0631af0305 1112 int iterations_;
RyoheiHagimoto 0:0e0631af0305 1113
RyoheiHagimoto 0:0e0631af0305 1114 /** Algorithm for choosing the cluster centers */
RyoheiHagimoto 0:0e0631af0305 1115 flann_centers_init_t centers_init_;
RyoheiHagimoto 0:0e0631af0305 1116
RyoheiHagimoto 0:0e0631af0305 1117 /**
RyoheiHagimoto 0:0e0631af0305 1118 * Cluster border index. This is used in the tree search phase when determining
RyoheiHagimoto 0:0e0631af0305 1119 * the closest cluster to explore next. A zero value takes into account only
RyoheiHagimoto 0:0e0631af0305 1120 * the cluster centres, a value greater then zero also take into account the size
RyoheiHagimoto 0:0e0631af0305 1121 * of the cluster.
RyoheiHagimoto 0:0e0631af0305 1122 */
RyoheiHagimoto 0:0e0631af0305 1123 float cb_index_;
RyoheiHagimoto 0:0e0631af0305 1124
RyoheiHagimoto 0:0e0631af0305 1125 /**
RyoheiHagimoto 0:0e0631af0305 1126 * The dataset used by this index
RyoheiHagimoto 0:0e0631af0305 1127 */
RyoheiHagimoto 0:0e0631af0305 1128 const Matrix<ElementType> dataset_;
RyoheiHagimoto 0:0e0631af0305 1129
RyoheiHagimoto 0:0e0631af0305 1130 /** Index parameters */
RyoheiHagimoto 0:0e0631af0305 1131 IndexParams index_params_;
RyoheiHagimoto 0:0e0631af0305 1132
RyoheiHagimoto 0:0e0631af0305 1133 /**
RyoheiHagimoto 0:0e0631af0305 1134 * Number of features in the dataset.
RyoheiHagimoto 0:0e0631af0305 1135 */
RyoheiHagimoto 0:0e0631af0305 1136 size_t size_;
RyoheiHagimoto 0:0e0631af0305 1137
RyoheiHagimoto 0:0e0631af0305 1138 /**
RyoheiHagimoto 0:0e0631af0305 1139 * Length of each feature.
RyoheiHagimoto 0:0e0631af0305 1140 */
RyoheiHagimoto 0:0e0631af0305 1141 size_t veclen_;
RyoheiHagimoto 0:0e0631af0305 1142
RyoheiHagimoto 0:0e0631af0305 1143 /**
RyoheiHagimoto 0:0e0631af0305 1144 * The root node in the tree.
RyoheiHagimoto 0:0e0631af0305 1145 */
RyoheiHagimoto 0:0e0631af0305 1146 KMeansNodePtr root_;
RyoheiHagimoto 0:0e0631af0305 1147
RyoheiHagimoto 0:0e0631af0305 1148 /**
RyoheiHagimoto 0:0e0631af0305 1149 * Array of indices to vectors in the dataset.
RyoheiHagimoto 0:0e0631af0305 1150 */
RyoheiHagimoto 0:0e0631af0305 1151 int* indices_;
RyoheiHagimoto 0:0e0631af0305 1152
RyoheiHagimoto 0:0e0631af0305 1153 /**
RyoheiHagimoto 0:0e0631af0305 1154 * The distance
RyoheiHagimoto 0:0e0631af0305 1155 */
RyoheiHagimoto 0:0e0631af0305 1156 Distance distance_;
RyoheiHagimoto 0:0e0631af0305 1157
RyoheiHagimoto 0:0e0631af0305 1158 /**
RyoheiHagimoto 0:0e0631af0305 1159 * Pooled memory allocator.
RyoheiHagimoto 0:0e0631af0305 1160 */
RyoheiHagimoto 0:0e0631af0305 1161 PooledAllocator pool_;
RyoheiHagimoto 0:0e0631af0305 1162
RyoheiHagimoto 0:0e0631af0305 1163 /**
RyoheiHagimoto 0:0e0631af0305 1164 * Memory occupied by the index.
RyoheiHagimoto 0:0e0631af0305 1165 */
RyoheiHagimoto 0:0e0631af0305 1166 int memoryCounter_;
RyoheiHagimoto 0:0e0631af0305 1167 };
RyoheiHagimoto 0:0e0631af0305 1168
RyoheiHagimoto 0:0e0631af0305 1169 }
RyoheiHagimoto 0:0e0631af0305 1170
RyoheiHagimoto 0:0e0631af0305 1171 #endif //OPENCV_FLANN_KMEANS_INDEX_H_