kmeans_index.h

/***********************************************************************





* Software License Agreement (BSD License)





*





* Copyright 2008-2009  Marius Muja (mariusm@cs.ubc.ca). All rights reserved.





* Copyright 2008-2009  David G. Lowe (lowe@cs.ubc.ca). All rights reserved.





*





* THE BSD LICENSE





*





* Redistribution and use in source and binary forms, with or without





* modification, are permitted provided that the following conditions





* are met:





*





* 1. Redistributions of source code must retain the above copyright





*    notice, this list of conditions and the following disclaimer.





* 2. Redistributions in binary form must reproduce the above copyright





*    notice, this list of conditions and the following disclaimer in the





*    documentation and/or other materials provided with the distribution.





*





* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR





* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES





* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.





* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,





* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT





* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,





* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY





* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT





* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF





* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.





*************************************************************************/








#ifndef OPENCV_FLANN_KMEANS_INDEX_H_




#define OPENCV_FLANN_KMEANS_INDEX_H_












#include <algorithm>




#include <map>




#include <limits>




#include <cmath>








#include "general.h"




#include "nn_index.h"




#include "dist.h"




#include "matrix.h"




#include "result_set.h"




#include "heap.h"




#include "allocator.h"




#include "random.h"




#include "saving.h"




#include "logger.h"








#define BITS_PER_CHAR 8




#define BITS_PER_BASE 2

// for DNA/RNA sequences




#define BASE_PER_CHAR (BITS_PER_CHAR/BITS_PER_BASE)




#define HISTOS_PER_BASE (1<<BITS_PER_BASE)












namespace
cvflann


{







struct
KMeansIndexParams :
public
IndexParams


{



KMeansIndexParams(int
branching = 32,
int
iterations = 11,



flann_centers_init_t centers_init = FLANN_CENTERS_RANDOM,



float
cb_index = 0.2,
int
trees = 1 )



{



(*this)["algorithm"] = FLANN_INDEX_KMEANS;



// branching factor




(*this)["branching"] = branching;



// max iterations to perform in one kmeans clustering (kmeans tree)




(*this)["iterations"] = iterations;



// algorithm used for picking the initial cluster centers for kmeans tree




(*this)["centers_init"] = centers_init;



// cluster boundary index. Used when searching the kmeans tree




(*this)["cb_index"] = cb_index;



// number of kmeans trees to search in




(*this)["trees"] = trees;



}


};











template
<typename
Distance>



class
KMeansIndex :
public
NNIndex<Distance>


{



public:



typedef
typename
Distance::ElementType ElementType;



typedef
typename
Distance::ResultType DistanceType;



typedef
typename
Distance::CentersType CentersType;







typedef
typename
Distance::is_kdtree_distance is_kdtree_distance;



typedef
typename
Distance::is_vector_space_distance is_vector_space_distance;















typedef
void (KMeansIndex::* centersAlgFunction)(int,
int*, int,
int*,
int&);







centersAlgFunction chooseCenters;















void
chooseCentersRandom(int
k,
int* indices,
int
indices_length,
int* centers,
int& centers_length)



{



UniqueRandom r(indices_length);







int
index;



for
(index=0; index<k; ++index) {



bool
duplicate =
true;



int
rnd;



while
(duplicate) {



duplicate =
false;



rnd = r.next();



if
(rnd<0) {



centers_length = index;



return;



}







centers[index] = indices[rnd];







for
(int
j=0; j<index; ++j) {



DistanceType sq = distance_(dataset_[centers[index]], dataset_[centers[j]], dataset_.cols);



if
(sq<1e-16) {



duplicate =
true;



}



}



}



}







centers_length = index;



}











void
chooseCentersGonzales(int
k,
int* indices,
int
indices_length,
int* centers,
int& centers_length)



{



int
n = indices_length;







int
rnd = rand_int(n);



CV_DbgAssert(rnd >=0 && rnd < n);







centers[0] = indices[rnd];







int
index;



for
(index=1; index<k; ++index) {







int
best_index = -1;



DistanceType best_val = 0;



for
(int
j=0; j<n; ++j) {



DistanceType dist = distance_(dataset_[centers[0]],dataset_[indices[j]],dataset_.cols);



for
(int
i=1; i<index; ++i) {



DistanceType tmp_dist = distance_(dataset_[centers[i]],dataset_[indices[j]],dataset_.cols);



if
(tmp_dist<dist) {



dist = tmp_dist;



}



}



if
(dist>best_val) {



best_val = dist;



best_index = j;



}



}



if
(best_index!=-1) {



centers[index] = indices[best_index];



}



else
{



break;



}



}



centers_length = index;



}











void
chooseCentersKMeanspp(int
k,
int* indices,
int
indices_length,
int* centers,
int& centers_length)



{



int
n = indices_length;







double
currentPot = 0;



DistanceType* closestDistSq =
new
DistanceType[n];







// Choose one random center and set the closestDistSq values




int
index = rand_int(n);



CV_DbgAssert(index >=0 && index < n);



centers[0] = indices[index];







for
(int
i = 0; i < n; i++) {



closestDistSq[i] = distance_(dataset_[indices[i]], dataset_[indices[index]], dataset_.cols);



closestDistSq[i] = ensureSquareDistance<Distance>( closestDistSq[i] );



currentPot += closestDistSq[i];



}











const
int
numLocalTries = 1;







// Choose each center




int
centerCount;



for
(centerCount = 1; centerCount < k; centerCount++) {







// Repeat several trials




double
bestNewPot = -1;



int
bestNewIndex = -1;



for
(int
localTrial = 0; localTrial < numLocalTries; localTrial++) {







// Choose our center - have to be slightly careful to return a valid answer even accounting




// for possible rounding errors




double
randVal = rand_double(currentPot);



for
(index = 0; index < n-1; index++) {



if
(randVal <= closestDistSq[index])
break;



else
randVal -= closestDistSq[index];



}







// Compute the new potential




double
newPot = 0;



for
(int
i = 0; i < n; i++) {



DistanceType dist = distance_(dataset_[indices[i]], dataset_[indices[index]], dataset_.cols);



newPot +=
std::min( ensureSquareDistance<Distance>(dist), closestDistSq[i] );



}







// Store the best result




if
((bestNewPot < 0)||(newPot < bestNewPot)) {



bestNewPot = newPot;



bestNewIndex = index;



}



}







// Add the appropriate center




centers[centerCount] = indices[bestNewIndex];



currentPot = bestNewPot;



for
(int
i = 0; i < n; i++) {



DistanceType dist = distance_(dataset_[indices[i]], dataset_[indices[bestNewIndex]], dataset_.cols);



closestDistSq[i] =
std::min( ensureSquareDistance<Distance>(dist), closestDistSq[i] );



}



}







centers_length = centerCount;







delete[] closestDistSq;



}















public:







flann_algorithm_t getType() const CV_OVERRIDE



{



return
FLANN_INDEX_KMEANS;



}







template<class
CentersContainerType>



class
KMeansDistanceComputer :
public
cv::ParallelLoopBody




{



public:



KMeansDistanceComputer(Distance _distance,
const
Matrix<ElementType>& _dataset,



const
int
_branching,
const
int* _indices,
const
CentersContainerType& _dcenters,



const
size_t
_veclen, std::vector<int> &_new_centroids,



std::vector<DistanceType> &_sq_dists)



: distance(_distance)



, dataset(_dataset)



, branching(_branching)



, indices(_indices)



, dcenters(_dcenters)



, veclen(_veclen)



, new_centroids(_new_centroids)



, sq_dists(_sq_dists)



{



}







void
operator()(const
cv::Range& range)
const
CV_OVERRIDE



{



const
int
begin = range.start;



const
int
end = range.end;







for(
int
i = begin; i<end; ++i)



{



DistanceType sq_dist(distance(dataset[indices[i]], dcenters[0], veclen));



int
new_centroid(0);



for
(int
j=1; j<branching; ++j) {



DistanceType new_sq_dist = distance(dataset[indices[i]], dcenters[j], veclen);



if
(sq_dist>new_sq_dist) {



new_centroid = j;



sq_dist = new_sq_dist;



}



}



sq_dists[i] = sq_dist;



new_centroids[i] = new_centroid;



}



}







private:



Distance distance;



const
Matrix<ElementType>& dataset;



const
int
branching;



const
int* indices;



const
CentersContainerType& dcenters;



const
size_t
veclen;



std::vector<int> &new_centroids;



std::vector<DistanceType> &sq_dists;



KMeansDistanceComputer& operator=(
const
KMeansDistanceComputer & ) {
return
*this; }



};







KMeansIndex(const
Matrix<ElementType>& inputData,
const
IndexParams& params = KMeansIndexParams(),



Distance d = Distance())



: dataset_(inputData), index_params_(params), root_(NULL), indices_(NULL), distance_(d)



{



memoryCounter_ = 0;







size_ = dataset_.rows;



veclen_ = dataset_.cols;







branching_ = get_param(params,"branching",32);



trees_ = get_param(params,"trees",1);



iterations_ = get_param(params,"iterations",11);



if
(iterations_<0) {



iterations_ = (std::numeric_limits<int>::max)();



}



centers_init_  = get_param(params,"centers_init",FLANN_CENTERS_RANDOM);







if
(centers_init_==FLANN_CENTERS_RANDOM) {



chooseCenters = &KMeansIndex::chooseCentersRandom;



}



else
if
(centers_init_==FLANN_CENTERS_GONZALES) {



chooseCenters = &KMeansIndex::chooseCentersGonzales;



}



else
if
(centers_init_==FLANN_CENTERS_KMEANSPP) {



chooseCenters = &KMeansIndex::chooseCentersKMeanspp;



}



else
{



FLANN_THROW(cv::Error::StsBadArg,
"Unknown algorithm for choosing initial centers.");



}



cb_index_ = 0.4f;







root_ =
new
KMeansNodePtr[trees_];



indices_ =
new
int*[trees_];







for
(int
i=0; i<trees_; ++i) {



root_[i] = NULL;



indices_[i] = NULL;



}



}











KMeansIndex(const
KMeansIndex&);



KMeansIndex& operator=(const
KMeansIndex&);











virtual
~KMeansIndex()



{



if
(root_ != NULL) {



free_centers();



delete[] root_;



}



if
(indices_!=NULL) {



free_indices();



delete[] indices_;



}



}







size_t
size() const CV_OVERRIDE



{



return
size_;



}







size_t
veclen() const CV_OVERRIDE



{



return
veclen_;



}











void
set_cb_index(
float
index)



{



cb_index_ = index;



}







int
usedMemory() const CV_OVERRIDE



{



return
pool_.usedMemory+pool_.wastedMemory+memoryCounter_;



}







void
buildIndex() CV_OVERRIDE



{



if
(branching_<2) {



FLANN_THROW(cv::Error::StsError,
"Branching factor must be at least 2");



}







free_indices();







for
(int
i=0; i<trees_; ++i) {



indices_[i] =
new
int[size_];



for
(size_t
j=0; j<size_; ++j) {



indices_[i][j] = int(j);



}



root_[i] = pool_.allocate<KMeansNode>();



std::memset(root_[i], 0,
sizeof(KMeansNode));







Distance* dummy = NULL;



computeNodeStatistics(root_[i], indices_[i], (unsigned
int)size_, dummy);







computeClustering(root_[i], indices_[i], (int)size_, branching_,0);



}



}











void
saveIndex(FILE* stream) CV_OVERRIDE



{



save_value(stream, branching_);



save_value(stream, iterations_);



save_value(stream, memoryCounter_);



save_value(stream, cb_index_);



save_value(stream, trees_);



for
(int
i=0; i<trees_; ++i) {



save_value(stream, *indices_[i], (int)size_);



save_tree(stream, root_[i], i);



}



}











void
loadIndex(FILE* stream) CV_OVERRIDE



{



if
(indices_!=NULL) {



free_indices();



delete[] indices_;



}



if
(root_!=NULL) {



free_centers();



}







load_value(stream, branching_);



load_value(stream, iterations_);



load_value(stream, memoryCounter_);



load_value(stream, cb_index_);



load_value(stream, trees_);







indices_ =
new
int*[trees_];



for
(int
i=0; i<trees_; ++i) {



indices_[i] =
new
int[size_];



load_value(stream, *indices_[i], size_);



load_tree(stream, root_[i], i);



}







index_params_["algorithm"] = getType();



index_params_["branching"] = branching_;



index_params_["trees"] = trees_;



index_params_["iterations"] = iterations_;



index_params_["centers_init"] = centers_init_;



index_params_["cb_index"] = cb_index_;



}











void
findNeighbors(ResultSet<DistanceType>& result,
const
ElementType* vec,
const
SearchParams& searchParams) CV_OVERRIDE



{







const
int
maxChecks = get_param(searchParams,"checks",32);







if
(maxChecks==FLANN_CHECKS_UNLIMITED) {



findExactNN(root_[0], result, vec);



}



else
{



// Priority queue storing intermediate branches in the best-bin-first search




Heap<BranchSt>* heap =
new
Heap<BranchSt>((int)size_);







int
checks = 0;



for
(int
i=0; i<trees_; ++i) {



findNN(root_[i], result, vec, checks, maxChecks, heap);



if
((checks >= maxChecks) && result.full())



break;



}







BranchSt branch;



while
(heap->popMin(branch) && (checks<maxChecks || !result.full())) {



KMeansNodePtr node = branch.node;



findNN(node, result, vec, checks, maxChecks, heap);



}



delete
heap;







CV_Assert(result.full());



}



}







int
getClusterCenters(Matrix<CentersType>& centers)



{



int
numClusters = centers.rows;



if
(numClusters<1) {



FLANN_THROW(cv::Error::StsBadArg,
"Number of clusters must be at least 1");



}







DistanceType variance;



KMeansNodePtr* clusters =
new
KMeansNodePtr[numClusters];







int
clusterCount = getMinVarianceClusters(root_[0], clusters, numClusters, variance);







Logger::info("Clusters requested: %d, returning %d\n",numClusters, clusterCount);







for
(int
i=0; i<clusterCount; ++i) {



CentersType* center = clusters[i]->pivot;



for
(size_t
j=0; j<veclen_; ++j) {



centers[i][j] = center[j];



}



}



delete[] clusters;







return
clusterCount;



}







IndexParams getParameters() const CV_OVERRIDE



{



return
index_params_;



}











private:



struct
KMeansNode



{



CentersType* pivot;



DistanceType radius;



DistanceType mean_radius;



DistanceType variance;



int
size;



KMeansNode** childs;



int* indices;



int
level;



};



typedef
KMeansNode* KMeansNodePtr;







typedef
BranchStruct<KMeansNodePtr, DistanceType> BranchSt;



















void
save_tree(FILE* stream, KMeansNodePtr node,
int
num)



{



save_value(stream, *node);



save_value(stream, *(node->pivot), (int)veclen_);



if
(node->childs==NULL) {



int
indices_offset = (int)(node->indices - indices_[num]);



save_value(stream, indices_offset);



}



else
{



for(int
i=0; i<branching_; ++i) {



save_tree(stream, node->childs[i], num);



}



}



}











void
load_tree(FILE* stream, KMeansNodePtr& node,
int
num)



{



node = pool_.allocate<KMeansNode>();



load_value(stream, *node);



node->pivot =
new
CentersType[veclen_];



load_value(stream, *(node->pivot), (int)veclen_);



if
(node->childs==NULL) {



int
indices_offset;



load_value(stream, indices_offset);



node->indices = indices_[num] + indices_offset;



}



else
{



node->childs = pool_.allocate<KMeansNodePtr>(branching_);



for(int
i=0; i<branching_; ++i) {



load_tree(stream, node->childs[i], num);



}



}



}











void
free_centers(KMeansNodePtr node)



{



delete[] node->pivot;



if
(node->childs!=NULL) {



for
(int
k=0; k<branching_; ++k) {



free_centers(node->childs[k]);



}



}



}







void
free_centers()



{



if
(root_ != NULL) {



for(int
i=0; i<trees_; ++i) {



if
(root_[i] != NULL) {



free_centers(root_[i]);



}



}



}



}







void
free_indices()



{



if
(indices_!=NULL) {



for(int
i=0; i<trees_; ++i) {



if
(indices_[i]!=NULL) {



delete[] indices_[i];



indices_[i] = NULL;



}



}



}



}







void
computeNodeStatistics(KMeansNodePtr node,
int* indices,
unsigned
int
indices_length)



{



DistanceType variance = 0;



CentersType* mean =
new
CentersType[veclen_];



memoryCounter_ += int(veclen_*sizeof(CentersType));







memset(mean,0,veclen_*sizeof(CentersType));







for
(unsigned
int
i=0; i<indices_length; ++i) {



ElementType* vec = dataset_[indices[i]];



for
(size_t
j=0; j<veclen_; ++j) {



mean[j] += vec[j];



}



variance += distance_(vec, ZeroIterator<ElementType>(), veclen_);



}



float
length =
static_cast<
float
>(indices_length);



for
(size_t
j=0; j<veclen_; ++j) {



mean[j] = cvflann::round<CentersType>( mean[j] /
static_cast<
double
>(indices_length) );



}



variance /=
static_cast<DistanceType>( length );



variance -= distance_(mean, ZeroIterator<ElementType>(), veclen_);







DistanceType radius = 0;



for
(unsigned
int
i=0; i<indices_length; ++i) {



DistanceType tmp = distance_(mean, dataset_[indices[i]], veclen_);



if
(tmp>radius) {



radius = tmp;



}



}







node->variance = variance;



node->radius = radius;



node->pivot = mean;



}











void
computeBitfieldNodeStatistics(KMeansNodePtr node,
int* indices,



unsigned
int
indices_length)



{



const
unsigned
int
accumulator_veclen =
static_cast<
unsigned
int
>(



veclen_*sizeof(CentersType)*BITS_PER_CHAR);







unsigned
long
long
variance = 0ull;



CentersType* mean =
new
CentersType[veclen_];



memoryCounter_ += int(veclen_*sizeof(CentersType));



unsigned
int* mean_accumulator =
new
unsigned
int[accumulator_veclen];







memset(mean_accumulator, 0,
sizeof(unsigned
int)*accumulator_veclen);







for
(unsigned
int
i=0; i<indices_length; ++i) {



variance +=
static_cast<
unsigned
long
long
>( ensureSquareDistance<Distance>(



distance_(dataset_[indices[i]], ZeroIterator<ElementType>(), veclen_)));



unsigned
char* vec = (unsigned
char*)dataset_[indices[i]];



for
(size_t
k=0, l=0; k<accumulator_veclen; k+=BITS_PER_CHAR, ++l) {



mean_accumulator[k]   += (vec[l])    & 0x01;



mean_accumulator[k+1] += (vec[l]>>1) & 0x01;



mean_accumulator[k+2] += (vec[l]>>2) & 0x01;



mean_accumulator[k+3] += (vec[l]>>3) & 0x01;



mean_accumulator[k+4] += (vec[l]>>4) & 0x01;



mean_accumulator[k+5] += (vec[l]>>5) & 0x01;



mean_accumulator[k+6] += (vec[l]>>6) & 0x01;



mean_accumulator[k+7] += (vec[l]>>7) & 0x01;



}



}



double
cnt =
static_cast<
double
>(indices_length);



unsigned
char* char_mean = (unsigned
char*)mean;



for
(size_t
k=0, l=0; k<accumulator_veclen; k+=BITS_PER_CHAR, ++l) {



char_mean[l] =
static_cast<
unsigned
char
>(



(((int)(0.5 + (double)(mean_accumulator[k])   / cnt)))



| (((int)(0.5 + (double)(mean_accumulator[k+1]) / cnt))<<1)



| (((int)(0.5 + (double)(mean_accumulator[k+2]) / cnt))<<2)



| (((int)(0.5 + (double)(mean_accumulator[k+3]) / cnt))<<3)



| (((int)(0.5 + (double)(mean_accumulator[k+4]) / cnt))<<4)



| (((int)(0.5 + (double)(mean_accumulator[k+5]) / cnt))<<5)



| (((int)(0.5 + (double)(mean_accumulator[k+6]) / cnt))<<6)



| (((int)(0.5 + (double)(mean_accumulator[k+7]) / cnt))<<7));



}



variance =
static_cast<
unsigned
long
long
>(



0.5 +
static_cast<
double
>(variance) /
static_cast<
double
>(indices_length));



variance -=
static_cast<
unsigned
long
long
>(



ensureSquareDistance<Distance>(



distance_(mean, ZeroIterator<ElementType>(), veclen_)));







DistanceType radius = 0;



for
(unsigned
int
i=0; i<indices_length; ++i) {



DistanceType tmp = distance_(mean, dataset_[indices[i]], veclen_);



if
(tmp>radius) {



radius = tmp;



}



}







node->variance =
static_cast<DistanceType>(variance);



node->radius = radius;



node->pivot = mean;







delete[] mean_accumulator;



}











void
computeDnaNodeStatistics(KMeansNodePtr node,
int* indices,



unsigned
int
indices_length)



{



const
unsigned
int
histos_veclen =
static_cast<
unsigned
int
>(



veclen_*sizeof(CentersType)*(HISTOS_PER_BASE*BASE_PER_CHAR));







unsigned
long
long
variance = 0ull;



unsigned
int* histograms =
new
unsigned
int[histos_veclen];



memset(histograms, 0,
sizeof(unsigned
int)*histos_veclen);







for
(unsigned
int
i=0; i<indices_length; ++i) {



variance +=
static_cast<
unsigned
long
long
>( ensureSquareDistance<Distance>(



distance_(dataset_[indices[i]], ZeroIterator<ElementType>(), veclen_)));







unsigned
char* vec = (unsigned
char*)dataset_[indices[i]];



for
(size_t
k=0, l=0; k<histos_veclen; k+=HISTOS_PER_BASE*BASE_PER_CHAR, ++l) {



histograms[k +     ((vec[l])    & 0x03)]++;



histograms[k + 4 + ((vec[l]>>2) & 0x03)]++;



histograms[k + 8 + ((vec[l]>>4) & 0x03)]++;



histograms[k +12 + ((vec[l]>>6) & 0x03)]++;



}



}







CentersType* mean =
new
CentersType[veclen_];



memoryCounter_ += int(veclen_*sizeof(CentersType));



unsigned
char* char_mean = (unsigned
char*)mean;



unsigned
int* h = histograms;



for
(size_t
k=0, l=0; k<histos_veclen; k+=HISTOS_PER_BASE*BASE_PER_CHAR, ++l) {



char_mean[l] = (h[k] > h[k+1] ? h[k+2] > h[k+3] ? h[k]   > h[k+2] ? 0x00 : 0x10



: h[k]   > h[k+3] ? 0x00 : 0x11



: h[k+2] > h[k+3] ? h[k+1] > h[k+2] ? 0x01 : 0x10



: h[k+1] > h[k+3] ? 0x01 : 0x11)



| (h[k+4]>h[k+5] ? h[k+6] > h[k+7] ? h[k+4] > h[k+6] ? 0x00   : 0x1000



: h[k+4] > h[k+7] ? 0x00   : 0x1100



: h[k+6] > h[k+7] ? h[k+5] > h[k+6] ? 0x0100 : 0x1000



: h[k+5] > h[k+7] ? 0x0100 : 0x1100)



| (h[k+8]>h[k+9] ? h[k+10]>h[k+11] ? h[k+8] >h[k+10] ? 0x00   : 0x100000



: h[k+8] >h[k+11] ? 0x00   : 0x110000



: h[k+10]>h[k+11] ? h[k+9] >h[k+10] ? 0x010000 : 0x100000



: h[k+9] >h[k+11] ? 0x010000 : 0x110000)



| (h[k+12]>h[k+13] ? h[k+14]>h[k+15] ? h[k+12] >h[k+14] ? 0x00   : 0x10000000



: h[k+12] >h[k+15] ? 0x00   : 0x11000000



: h[k+14]>h[k+15] ? h[k+13] >h[k+14] ? 0x01000000 : 0x10000000



: h[k+13] >h[k+15] ? 0x01000000 : 0x11000000);



}



variance =
static_cast<
unsigned
long
long
>(



0.5 +
static_cast<
double
>(variance) /
static_cast<
double
>(indices_length));



variance -=
static_cast<
unsigned
long
long
>(



ensureSquareDistance<Distance>(



distance_(mean, ZeroIterator<ElementType>(), veclen_)));







DistanceType radius = 0;



for
(unsigned
int
i=0; i<indices_length; ++i) {



DistanceType tmp = distance_(mean, dataset_[indices[i]], veclen_);



if
(tmp>radius) {



radius = tmp;



}



}







node->variance =
static_cast<DistanceType>(variance);



node->radius = radius;



node->pivot = mean;







delete[] histograms;



}











template<typename
DistType>



void
computeNodeStatistics(KMeansNodePtr node,
int* indices,



unsigned
int
indices_length,



const
DistType* identifier)



{



(void)identifier;



computeNodeStatistics(node, indices, indices_length);



}







void
computeNodeStatistics(KMeansNodePtr node,
int* indices,



unsigned
int
indices_length,



const
cvflann::HammingLUT* identifier)



{



(void)identifier;



computeBitfieldNodeStatistics(node, indices, indices_length);



}







void
computeNodeStatistics(KMeansNodePtr node,
int* indices,



unsigned
int
indices_length,



const
cvflann::Hamming<unsigned char>* identifier)



{



(void)identifier;



computeBitfieldNodeStatistics(node, indices, indices_length);



}







void
computeNodeStatistics(KMeansNodePtr node,
int* indices,



unsigned
int
indices_length,



const
cvflann::Hamming2<unsigned char>* identifier)



{



(void)identifier;



computeBitfieldNodeStatistics(node, indices, indices_length);



}







void
computeNodeStatistics(KMeansNodePtr node,
int* indices,



unsigned
int
indices_length,



const
cvflann::DNAmmingLUT* identifier)



{



(void)identifier;



computeDnaNodeStatistics(node, indices, indices_length);



}







void
computeNodeStatistics(KMeansNodePtr node,
int* indices,



unsigned
int
indices_length,



const
cvflann::DNAmming2<unsigned char>* identifier)



{



(void)identifier;



computeDnaNodeStatistics(node, indices, indices_length);



}











void
refineClustering(int* indices,
int
indices_length,
int
branching, CentersType** centers,



std::vector<DistanceType>& radiuses,
int* belongs_to,
int* count)



{



cv::AutoBuffer<double>
dcenters_buf(branching*veclen_);



Matrix<double> dcenters(dcenters_buf.data(), branching, veclen_);







bool
converged =
false;



int
iteration = 0;



while
(!converged && iteration<iterations_) {



converged =
true;



iteration++;







// compute the new cluster centers




for
(int
i=0; i<branching; ++i) {



memset(dcenters[i],0,sizeof(double)*veclen_);



radiuses[i] = 0;



}



for
(int
i=0; i<indices_length; ++i) {



ElementType* vec = dataset_[indices[i]];



double* center = dcenters[belongs_to[i]];



for
(size_t
k=0; k<veclen_; ++k) {



center[k] += vec[k];



}



}



for
(int
i=0; i<branching; ++i) {



int
cnt = count[i];



for
(size_t
k=0; k<veclen_; ++k) {



dcenters[i][k] /= cnt;



}



}







std::vector<int> new_centroids(indices_length);



std::vector<DistanceType> sq_dists(indices_length);







// reassign points to clusters




KMeansDistanceComputer<Matrix<double> > invoker(



distance_, dataset_, branching, indices, dcenters, veclen_, new_centroids, sq_dists);



parallel_for_(cv::Range(0, (int)indices_length), invoker);







for
(int
i=0; i < (int)indices_length; ++i) {



DistanceType sq_dist(sq_dists[i]);



int
new_centroid(new_centroids[i]);



if
(sq_dist > radiuses[new_centroid]) {



radiuses[new_centroid] = sq_dist;



}



if
(new_centroid != belongs_to[i]) {



count[belongs_to[i]]--;



count[new_centroid]++;



belongs_to[i] = new_centroid;



converged =
false;



}



}







for
(int
i=0; i<branching; ++i) {



// if one cluster converges to an empty cluster,




// move an element into that cluster




if
(count[i]==0) {



int
j = (i+1)%branching;



while
(count[j]<=1) {



j = (j+1)%branching;



}







for
(int
k=0; k<indices_length; ++k) {



if
(belongs_to[k]==j) {



// for cluster j, we move the furthest element from the center to the empty cluster i




if
( distance_(dataset_[indices[k]], dcenters[j], veclen_) == radiuses[j] ) {



belongs_to[k] = i;



count[j]--;



count[i]++;



break;



}



}



}



converged =
false;



}



}



}







for
(int
i=0; i<branching; ++i) {



centers[i] =
new
CentersType[veclen_];



memoryCounter_ += (int)(veclen_*sizeof(CentersType));



for
(size_t
k=0; k<veclen_; ++k) {



centers[i][k] = (CentersType)dcenters[i][k];



}



}



}











void
refineBitfieldClustering(int* indices,
int
indices_length,
int
branching, CentersType** centers,



std::vector<DistanceType>& radiuses,
int* belongs_to,
int* count)



{



for
(int
i=0; i<branching; ++i) {



centers[i] =
new
CentersType[veclen_];



memoryCounter_ += (int)(veclen_*sizeof(CentersType));



}







const
unsigned
int
accumulator_veclen =
static_cast<
unsigned
int
>(



veclen_*sizeof(ElementType)*BITS_PER_CHAR);



cv::AutoBuffer<unsigned int>
dcenters_buf(branching*accumulator_veclen);



Matrix<unsigned int> dcenters(dcenters_buf.data(), branching, accumulator_veclen);







bool
converged =
false;



int
iteration = 0;



while
(!converged && iteration<iterations_) {



converged =
true;



iteration++;







// compute the new cluster centers




for
(int
i=0; i<branching; ++i) {



memset(dcenters[i],0,sizeof(unsigned
int)*accumulator_veclen);



radiuses[i] = 0;



}



for
(int
i=0; i<indices_length; ++i) {



unsigned
char* vec = (unsigned
char*)dataset_[indices[i]];



unsigned
int* dcenter = dcenters[belongs_to[i]];



for
(size_t
k=0, l=0; k<accumulator_veclen; k+=BITS_PER_CHAR, ++l) {



dcenter[k]   += (vec[l])    & 0x01;



dcenter[k+1] += (vec[l]>>1) & 0x01;



dcenter[k+2] += (vec[l]>>2) & 0x01;



dcenter[k+3] += (vec[l]>>3) & 0x01;



dcenter[k+4] += (vec[l]>>4) & 0x01;



dcenter[k+5] += (vec[l]>>5) & 0x01;



dcenter[k+6] += (vec[l]>>6) & 0x01;



dcenter[k+7] += (vec[l]>>7) & 0x01;



}



}



for
(int
i=0; i<branching; ++i) {



double
cnt =
static_cast<
double
>(count[i]);



unsigned
int* dcenter = dcenters[i];



unsigned
char* charCenter = (unsigned
char*)centers[i];



for
(size_t
k=0, l=0; k<accumulator_veclen; k+=BITS_PER_CHAR, ++l) {



charCenter[l] =
static_cast<
unsigned
char
>(



(((int)(0.5 + (double)(dcenter[k])   / cnt)))



| (((int)(0.5 + (double)(dcenter[k+1]) / cnt))<<1)



| (((int)(0.5 + (double)(dcenter[k+2]) / cnt))<<2)



| (((int)(0.5 + (double)(dcenter[k+3]) / cnt))<<3)



| (((int)(0.5 + (double)(dcenter[k+4]) / cnt))<<4)



| (((int)(0.5 + (double)(dcenter[k+5]) / cnt))<<5)



| (((int)(0.5 + (double)(dcenter[k+6]) / cnt))<<6)



| (((int)(0.5 + (double)(dcenter[k+7]) / cnt))<<7));



}



}







std::vector<int> new_centroids(indices_length);



std::vector<DistanceType> dists(indices_length);







// reassign points to clusters




KMeansDistanceComputer<ElementType**> invoker(



distance_, dataset_, branching, indices, centers, veclen_, new_centroids, dists);



parallel_for_(cv::Range(0, (int)indices_length), invoker);







for
(int
i=0; i < indices_length; ++i) {



DistanceType dist(dists[i]);



int
new_centroid(new_centroids[i]);



if
(dist > radiuses[new_centroid]) {



radiuses[new_centroid] = dist;



}



if
(new_centroid != belongs_to[i]) {



count[belongs_to[i]]--;



count[new_centroid]++;



belongs_to[i] = new_centroid;



converged =
false;



}



}







for
(int
i=0; i<branching; ++i) {



// if one cluster converges to an empty cluster,




// move an element into that cluster




if
(count[i]==0) {



int
j = (i+1)%branching;



while
(count[j]<=1) {



j = (j+1)%branching;



}







for
(int
k=0; k<indices_length; ++k) {



if
(belongs_to[k]==j) {



// for cluster j, we move the furthest element from the center to the empty cluster i




if
( distance_(dataset_[indices[k]], centers[j], veclen_) == radiuses[j] ) {



belongs_to[k] = i;



count[j]--;



count[i]++;



break;



}



}



}



converged =
false;



}



}



}



}











void
refineDnaClustering(int* indices,
int
indices_length,
int
branching, CentersType** centers,



std::vector<DistanceType>& radiuses,
int* belongs_to,
int* count)



{



for
(int
i=0; i<branching; ++i) {



centers[i] =
new
CentersType[veclen_];



memoryCounter_ += (int)(veclen_*sizeof(CentersType));



}







const
unsigned
int
histos_veclen =
static_cast<
unsigned
int
>(



veclen_*sizeof(CentersType)*(HISTOS_PER_BASE*BASE_PER_CHAR));



cv::AutoBuffer<unsigned int>
histos_buf(branching*histos_veclen);



Matrix<unsigned int> histos(histos_buf.data(), branching, histos_veclen);







bool
converged =
false;



int
iteration = 0;



while
(!converged && iteration<iterations_) {



converged =
true;



iteration++;







// compute the new cluster centers




for
(int
i=0; i<branching; ++i) {



memset(histos[i],0,sizeof(unsigned
int)*histos_veclen);



radiuses[i] = 0;



}



for
(int
i=0; i<indices_length; ++i) {



unsigned
char* vec = (unsigned
char*)dataset_[indices[i]];



unsigned
int* h = histos[belongs_to[i]];



for
(size_t
k=0, l=0; k<histos_veclen; k+=HISTOS_PER_BASE*BASE_PER_CHAR, ++l) {



h[k +     ((vec[l])    & 0x03)]++;



h[k + 4 + ((vec[l]>>2) & 0x03)]++;



h[k + 8 + ((vec[l]>>4) & 0x03)]++;



h[k +12 + ((vec[l]>>6) & 0x03)]++;



}



}



for
(int
i=0; i<branching; ++i) {



unsigned
int* h = histos[i];



unsigned
char* charCenter = (unsigned
char*)centers[i];



for
(size_t
k=0, l=0; k<histos_veclen; k+=HISTOS_PER_BASE*BASE_PER_CHAR, ++l) {



charCenter[l]= (h[k] > h[k+1] ? h[k+2] > h[k+3] ? h[k]   > h[k+2] ? 0x00 : 0x10



: h[k]   > h[k+3] ? 0x00 : 0x11



: h[k+2] > h[k+3] ? h[k+1] > h[k+2] ? 0x01 : 0x10



: h[k+1] > h[k+3] ? 0x01 : 0x11)



| (h[k+4]>h[k+5] ? h[k+6] > h[k+7] ? h[k+4] > h[k+6] ? 0x00   : 0x1000



: h[k+4] > h[k+7] ? 0x00   : 0x1100



: h[k+6] > h[k+7] ? h[k+5] > h[k+6] ? 0x0100 : 0x1000



: h[k+5] > h[k+7] ? 0x0100 : 0x1100)



| (h[k+8]>h[k+9] ? h[k+10]>h[k+11] ? h[k+8] >h[k+10] ? 0x00   : 0x100000



: h[k+8] >h[k+11] ? 0x00   : 0x110000



: h[k+10]>h[k+11] ? h[k+9] >h[k+10] ? 0x010000 : 0x100000



: h[k+9] >h[k+11] ? 0x010000 : 0x110000)



| (h[k+12]>h[k+13] ? h[k+14]>h[k+15] ? h[k+12] >h[k+14] ? 0x00   : 0x10000000



: h[k+12] >h[k+15] ? 0x00   : 0x11000000



: h[k+14]>h[k+15] ? h[k+13] >h[k+14] ? 0x01000000 : 0x10000000



: h[k+13] >h[k+15] ? 0x01000000 : 0x11000000);



}



}







std::vector<int> new_centroids(indices_length);



std::vector<DistanceType> dists(indices_length);







// reassign points to clusters




KMeansDistanceComputer<ElementType**> invoker(



distance_, dataset_, branching, indices, centers, veclen_, new_centroids, dists);



parallel_for_(cv::Range(0, (int)indices_length), invoker);







for
(int
i=0; i < indices_length; ++i) {



DistanceType dist(dists[i]);



int
new_centroid(new_centroids[i]);



if
(dist > radiuses[new_centroid]) {



radiuses[new_centroid] = dist;



}



if
(new_centroid != belongs_to[i]) {



count[belongs_to[i]]--;



count[new_centroid]++;



belongs_to[i] = new_centroid;



converged =
false;



}



}







for
(int
i=0; i<branching; ++i) {



// if one cluster converges to an empty cluster,




// move an element into that cluster




if
(count[i]==0) {



int
j = (i+1)%branching;



while
(count[j]<=1) {



j = (j+1)%branching;



}







for
(int
k=0; k<indices_length; ++k) {



if
(belongs_to[k]==j) {



// for cluster j, we move the furthest element from the center to the empty cluster i




if
( distance_(dataset_[indices[k]], centers[j], veclen_) == radiuses[j] ) {



belongs_to[k] = i;



count[j]--;



count[i]++;



break;



}



}



}



converged =
false;



}



}



}



}











void
computeSubClustering(KMeansNodePtr node,
int* indices,
int
indices_length,



int
branching,
int
level, CentersType** centers,



std::vector<DistanceType>& radiuses,
int* belongs_to,
int* count)



{



// compute kmeans clustering for each of the resulting clusters




node->childs = pool_.allocate<KMeansNodePtr>(branching);



int
start = 0;



int
end = start;



for
(int
c=0; c<branching; ++c) {



int
s = count[c];







DistanceType variance = 0;



DistanceType mean_radius =0;



for
(int
i=0; i<indices_length; ++i) {



if
(belongs_to[i]==c) {



DistanceType d = distance_(dataset_[indices[i]], ZeroIterator<ElementType>(), veclen_);



variance += d;



mean_radius +=
static_cast<DistanceType>(
sqrt(d) );



std::swap(indices[i],indices[end]);



std::swap(belongs_to[i],belongs_to[end]);



end++;



}



}



variance /= s;



mean_radius /= s;



variance -= distance_(centers[c], ZeroIterator<ElementType>(), veclen_);







node->childs[c] = pool_.allocate<KMeansNode>();



std::memset(node->childs[c], 0,
sizeof(KMeansNode));



node->childs[c]->radius = radiuses[c];



node->childs[c]->pivot = centers[c];



node->childs[c]->variance = variance;



node->childs[c]->mean_radius = mean_radius;



computeClustering(node->childs[c],indices+start, end-start, branching, level+1);



start=end;



}



}











void
computeAnyBitfieldSubClustering(KMeansNodePtr node,
int* indices,
int
indices_length,



int
branching,
int
level, CentersType** centers,



std::vector<DistanceType>& radiuses,
int* belongs_to,
int* count)



{



// compute kmeans clustering for each of the resulting clusters




node->childs = pool_.allocate<KMeansNodePtr>(branching);



int
start = 0;



int
end = start;



for
(int
c=0; c<branching; ++c) {



int
s = count[c];







unsigned
long
long
variance = 0ull;



DistanceType mean_radius =0;



for
(int
i=0; i<indices_length; ++i) {



if
(belongs_to[i]==c) {



DistanceType d = distance_(dataset_[indices[i]], ZeroIterator<ElementType>(), veclen_);



variance +=
static_cast<
unsigned
long
long
>( ensureSquareDistance<Distance>(d) );



mean_radius += ensureSimpleDistance<Distance>(d);



std::swap(indices[i],indices[end]);



std::swap(belongs_to[i],belongs_to[end]);



end++;



}



}



mean_radius =
static_cast<DistanceType>(



0.5f +
static_cast<
float
>(mean_radius) /
static_cast<
float
>(s));



variance =
static_cast<
unsigned
long
long
>(



0.5 +
static_cast<
double
>(variance) /
static_cast<
double
>(s));



variance -=
static_cast<
unsigned
long
long
>(



ensureSquareDistance<Distance>(



distance_(centers[c], ZeroIterator<ElementType>(), veclen_)));







node->childs[c] = pool_.allocate<KMeansNode>();



std::memset(node->childs[c], 0,
sizeof(KMeansNode));



node->childs[c]->radius = radiuses[c];



node->childs[c]->pivot = centers[c];



node->childs[c]->variance =
static_cast<DistanceType>(variance);



node->childs[c]->mean_radius = mean_radius;



computeClustering(node->childs[c],indices+start, end-start, branching, level+1);



start=end;



}



}











template<typename
DistType>



void
refineAndSplitClustering(



KMeansNodePtr node,
int* indices,
int
indices_length,
int
branching,



int
level, CentersType** centers, std::vector<DistanceType>& radiuses,



int* belongs_to,
int* count,
const
DistType* identifier)



{



(void)identifier;



refineClustering(indices, indices_length, branching, centers, radiuses, belongs_to, count);







computeSubClustering(node, indices, indices_length, branching,



level, centers, radiuses, belongs_to, count);



}











void
refineAndSplitClustering(



KMeansNodePtr node,
int* indices,
int
indices_length,
int
branching,



int
level, CentersType** centers, std::vector<DistanceType>& radiuses,



int* belongs_to,
int* count,
const
cvflann::HammingLUT* identifier)



{



(void)identifier;



refineBitfieldClustering(



indices, indices_length, branching, centers, radiuses, belongs_to, count);







computeAnyBitfieldSubClustering(node, indices, indices_length, branching,



level, centers, radiuses, belongs_to, count);



}











void
refineAndSplitClustering(



KMeansNodePtr node,
int* indices,
int
indices_length,
int
branching,



int
level, CentersType** centers, std::vector<DistanceType>& radiuses,



int* belongs_to,
int* count,
const
cvflann::Hamming<unsigned char>* identifier)



{



(void)identifier;



refineBitfieldClustering(



indices, indices_length, branching, centers, radiuses, belongs_to, count);







computeAnyBitfieldSubClustering(node, indices, indices_length, branching,



level, centers, radiuses, belongs_to, count);



}











void
refineAndSplitClustering(



KMeansNodePtr node,
int* indices,
int
indices_length,
int
branching,



int
level, CentersType** centers, std::vector<DistanceType>& radiuses,



int* belongs_to,
int* count,
const
cvflann::Hamming2<unsigned char>* identifier)



{



(void)identifier;



refineBitfieldClustering(



indices, indices_length, branching, centers, radiuses, belongs_to, count);







computeAnyBitfieldSubClustering(node, indices, indices_length, branching,



level, centers, radiuses, belongs_to, count);



}











void
refineAndSplitClustering(



KMeansNodePtr node,
int* indices,
int
indices_length,
int
branching,



int
level, CentersType** centers, std::vector<DistanceType>& radiuses,



int* belongs_to,
int* count,
const
cvflann::DNAmmingLUT* identifier)



{



(void)identifier;



refineDnaClustering(



indices, indices_length, branching, centers, radiuses, belongs_to, count);







computeAnyBitfieldSubClustering(node, indices, indices_length, branching,



level, centers, radiuses, belongs_to, count);



}











void
refineAndSplitClustering(



KMeansNodePtr node,
int* indices,
int
indices_length,
int
branching,



int
level, CentersType** centers, std::vector<DistanceType>& radiuses,



int* belongs_to,
int* count,
const
cvflann::DNAmming2<unsigned char>* identifier)



{



(void)identifier;



refineDnaClustering(



indices, indices_length, branching, centers, radiuses, belongs_to, count);







computeAnyBitfieldSubClustering(node, indices, indices_length, branching,



level, centers, radiuses, belongs_to, count);



}











void
computeClustering(KMeansNodePtr node,
int* indices,
int
indices_length,
int
branching,
int
level)



{



node->size = indices_length;



node->level = level;







if
(indices_length < branching) {



node->indices = indices;



std::sort(node->indices,node->indices+indices_length);



node->childs = NULL;



return;



}







cv::AutoBuffer<int>
centers_idx_buf(branching);



int* centers_idx = centers_idx_buf.data();



int
centers_length;



(this->*chooseCenters)(branching, indices, indices_length, centers_idx, centers_length);







if
(centers_length<branching) {



node->indices = indices;



std::sort(node->indices,node->indices+indices_length);



node->childs = NULL;



return;



}











std::vector<DistanceType> radiuses(branching);



cv::AutoBuffer<int>
count_buf(branching);



int* count = count_buf.data();



for
(int
i=0; i<branching; ++i) {



radiuses[i] = 0;



count[i] = 0;



}







//  assign points to clusters




cv::AutoBuffer<int>
belongs_to_buf(indices_length);



int* belongs_to = belongs_to_buf.data();



for
(int
i=0; i<indices_length; ++i) {



DistanceType sq_dist = distance_(dataset_[indices[i]], dataset_[centers_idx[0]], veclen_);



belongs_to[i] = 0;



for
(int
j=1; j<branching; ++j) {



DistanceType new_sq_dist = distance_(dataset_[indices[i]], dataset_[centers_idx[j]], veclen_);



if
(sq_dist>new_sq_dist) {



belongs_to[i] = j;



sq_dist = new_sq_dist;



}



}



if
(sq_dist>radiuses[belongs_to[i]]) {



radiuses[belongs_to[i]] = sq_dist;



}



count[belongs_to[i]]++;



}







CentersType** centers =
new
CentersType*[branching];







Distance* dummy = NULL;



refineAndSplitClustering(node, indices, indices_length, branching, level,



centers, radiuses, belongs_to, count, dummy);







delete[] centers;



}











void
findNN(KMeansNodePtr node, ResultSet<DistanceType>& result,
const
ElementType* vec,
int& checks,
int
maxChecks,



Heap<BranchSt>* heap)



{



// Ignore those clusters that are too far away




{



DistanceType bsq = distance_(vec, node->pivot, veclen_);



DistanceType rsq = node->radius;



DistanceType wsq = result.worstDist();







if
(isSquareDistance<Distance>())



{



DistanceType val = bsq-rsq-wsq;



if
((val>0) && (val*val > 4*rsq*wsq))



return;



}



else




{



if
(bsq-rsq > wsq)



return;



}



}







if
(node->childs==NULL) {



if
((checks>=maxChecks) && result.full()) {



return;



}



checks += node->size;



for
(int
i=0; i<node->size; ++i) {



int
index = node->indices[i];



DistanceType dist = distance_(dataset_[index], vec, veclen_);



result.addPoint(dist, index);



}



}



else
{



DistanceType* domain_distances =
new
DistanceType[branching_];



int
closest_center = exploreNodeBranches(node, vec, domain_distances, heap);



delete[] domain_distances;



findNN(node->childs[closest_center],result,vec, checks, maxChecks, heap);



}



}







int
exploreNodeBranches(KMeansNodePtr node,
const
ElementType* q, DistanceType* domain_distances, Heap<BranchSt>* heap)



{







int
best_index = 0;



domain_distances[best_index] = distance_(q, node->childs[best_index]->pivot, veclen_);



for
(int
i=1; i<branching_; ++i) {



domain_distances[i] = distance_(q, node->childs[i]->pivot, veclen_);



if
(domain_distances[i]<domain_distances[best_index]) {



best_index = i;



}



}







//      float* best_center = node->childs[best_index]->pivot;




for
(int
i=0; i<branching_; ++i) {



if
(i != best_index) {



domain_distances[i] -= cvflann::round<DistanceType>(



cb_index_*node->childs[i]->variance );







//              float dist_to_border = getDistanceToBorder(node.childs[i].pivot,best_center,q);




//              if (domain_distances[i]<dist_to_border) {




//                  domain_distances[i] = dist_to_border;




//              }




heap->insert(BranchSt(node->childs[i],domain_distances[i]));



}



}







return
best_index;



}











void
findExactNN(KMeansNodePtr node, ResultSet<DistanceType>& result,
const
ElementType* vec)



{



// Ignore those clusters that are too far away




{



DistanceType bsq = distance_(vec, node->pivot, veclen_);



DistanceType rsq = node->radius;



DistanceType wsq = result.worstDist();







if
(isSquareDistance<Distance>())



{



DistanceType val = bsq-rsq-wsq;



if
((val>0) && (val*val > 4*rsq*wsq))



return;



}



else




{



if
(bsq-rsq > wsq)



return;



}



}











if
(node->childs==NULL) {



for
(int
i=0; i<node->size; ++i) {



int
index = node->indices[i];



DistanceType dist = distance_(dataset_[index], vec, veclen_);



result.addPoint(dist, index);



}



}



else
{



int* sort_indices =
new
int[branching_];







getCenterOrdering(node, vec, sort_indices);







for
(int
i=0; i<branching_; ++i) {



findExactNN(node->childs[sort_indices[i]],result,vec);



}







delete[] sort_indices;



}



}











void
getCenterOrdering(KMeansNodePtr node,
const
ElementType* q,
int* sort_indices)



{



DistanceType* domain_distances =
new
DistanceType[branching_];



for
(int
i=0; i<branching_; ++i) {



DistanceType dist = distance_(q, node->childs[i]->pivot, veclen_);







int
j=0;



while
(domain_distances[j]<dist && j<i)



j++;



for
(int
k=i; k>j; --k) {



domain_distances[k] = domain_distances[k-1];



sort_indices[k] = sort_indices[k-1];



}



domain_distances[j] = dist;



sort_indices[j] = i;



}



delete[] domain_distances;



}







DistanceType getDistanceToBorder(DistanceType* p, DistanceType* c, DistanceType* q)



{



DistanceType sum = 0;



DistanceType sum2 = 0;







for
(int
i=0; i<veclen_; ++i) {



DistanceType t = c[i]-p[i];



sum += t*(q[i]-(c[i]+p[i])/2);



sum2 += t*t;



}







return
sum*sum/sum2;



}











int
getMinVarianceClusters(KMeansNodePtr root, KMeansNodePtr* clusters,
int
clusters_length, DistanceType& varianceValue)



{



int
clusterCount = 1;



clusters[0] = root;







DistanceType meanVariance = root->variance*root->size;







while
(clusterCount<clusters_length) {



DistanceType minVariance = (std::numeric_limits<DistanceType>::max)();



int
splitIndex = -1;







for
(int
i=0; i<clusterCount; ++i) {



if
(clusters[i]->childs != NULL) {







DistanceType variance = meanVariance - clusters[i]->variance*clusters[i]->size;







for
(int
j=0; j<branching_; ++j) {



variance += clusters[i]->childs[j]->variance*clusters[i]->childs[j]->size;



}



if
(variance<minVariance) {



minVariance = variance;



splitIndex = i;



}



}



}







if
(splitIndex==-1)
break;



if
( (branching_+clusterCount-1) > clusters_length)
break;







meanVariance = minVariance;







// split node




KMeansNodePtr toSplit = clusters[splitIndex];



clusters[splitIndex] = toSplit->childs[0];



for
(int
i=1; i<branching_; ++i) {



clusters[clusterCount++] = toSplit->childs[i];



}



}







varianceValue = meanVariance/root->size;



return
clusterCount;



}







private:



int
branching_;







int
trees_;







int
iterations_;







flann_centers_init_t centers_init_;







float
cb_index_;







const
Matrix<ElementType> dataset_;







IndexParams index_params_;







size_t
size_;







size_t
veclen_;







KMeansNodePtr* root_;







int** indices_;







Distance distance_;







PooledAllocator pool_;







int
memoryCounter_;


};






}











#endif

//OPENCV_FLANN_KMEANS_INDEX_H_