@param features The points to be clustered. The matrix must have elements of type
Distance::ElementType.
@param centers The centers of the clusters obtained. The matrix must have type
-Distance::ResultType. The number of rows in this matrix represents the number of clusters desired,
+Distance::CentersType. The number of rows in this matrix represents the number of clusters desired,
however, because of the way the cut in the hierarchical tree is chosen, the number of clusters
computed will be the highest number of the form (branching-1)\*k+1 that's lower than the number of
clusters desired, where branching is the tree's branching factor (see description of the
Distance d = Distance())
{
typedef typename Distance::ElementType ElementType;
- typedef typename Distance::ResultType DistanceType;
+ typedef typename Distance::CentersType CentersType;
CV_Assert(features.type() == CvType<ElementType>::type());
CV_Assert(features.isContinuous());
::cvflann::Matrix<ElementType> m_features((ElementType*)features.ptr<ElementType>(0), features.rows, features.cols);
- CV_Assert(centers.type() == CvType<DistanceType>::type());
+ CV_Assert(centers.type() == CvType<CentersType>::type());
CV_Assert(centers.isContinuous());
- ::cvflann::Matrix<DistanceType> m_centers((DistanceType*)centers.ptr<DistanceType>(0), centers.rows, centers.cols);
+ ::cvflann::Matrix<CentersType> m_centers((CentersType*)centers.ptr<CentersType>(0), centers.rows, centers.cols);
return ::cvflann::hierarchicalClustering<Distance>(m_features, m_centers, params, d);
}
case FLANN_INDEX_LINEAR:
nnIndex = new LinearIndex<Distance>(dataset, params, distance);
break;
+ case FLANN_INDEX_KMEANS:
+ nnIndex = new KMeansIndex<Distance>(dataset, params, distance);
+ break;
case FLANN_INDEX_HIERARCHICAL:
nnIndex = new HierarchicalClusteringIndex<Distance>(dataset, params, distance);
break;
-/***********************************************************************
+/***********************************************************************
* Software License Agreement (BSD License)
*
* Copyright 2008-2009 Marius Muja (mariusm@cs.ubc.ca). All rights reserved.
template<>
inline double abs<double>(double x) { return fabs(x); }
+
+template<typename TargetType>
+inline TargetType round(float x) { return static_cast<TargetType>(x); }
+
+template<>
+inline unsigned int round<unsigned int>(float x) { return static_cast<unsigned int>(x + 0.5f); }
+
+template<>
+inline unsigned short round<unsigned short>(float x) { return static_cast<unsigned short>(x + 0.5f); }
+
+template<>
+inline unsigned char round<unsigned char>(float x) { return static_cast<unsigned char>(x + 0.5f); }
+
+template<>
+inline long long round<long long>(float x) { return static_cast<long long>(x + 0.5f); }
+
+template<>
+inline long round<long>(float x) { return static_cast<long>(x + 0.5f); }
+
+template<>
+inline int round<int>(float x) { return static_cast<int>(x + 0.5f) - (x<0); }
+
+template<>
+inline short round<short>(float x) { return static_cast<short>(x + 0.5f) - (x<0); }
+
+template<>
+inline char round<char>(float x) { return static_cast<char>(x + 0.5f) - (x<0); }
+
+
+template<typename TargetType>
+inline TargetType round(double x) { return static_cast<TargetType>(x); }
+
+template<>
+inline unsigned int round<unsigned int>(double x) { return static_cast<unsigned int>(x + 0.5); }
+
+template<>
+inline unsigned short round<unsigned short>(double x) { return static_cast<unsigned short>(x + 0.5); }
+
+template<>
+inline unsigned char round<unsigned char>(double x) { return static_cast<unsigned char>(x + 0.5); }
+
+template<>
+inline long long round<long long>(double x) { return static_cast<long long>(x + 0.5); }
+
+template<>
+inline long round<long>(double x) { return static_cast<long>(x + 0.5); }
+
+template<>
+inline int round<int>(double x) { return static_cast<int>(x + 0.5) - (x<0); }
+
+template<>
+inline short round<short>(double x) { return static_cast<short>(x + 0.5) - (x<0); }
+
+template<>
+inline char round<char>(double x) { return static_cast<char>(x + 0.5) - (x<0); }
+
+
template<typename T>
struct Accumulator { typedef T Type; };
template<>
class True
{
+public:
+ static const bool val = true;
};
class False
{
+public:
+ static const bool val = false;
};
+/*
+ * This is a "zero iterator". It basically behaves like a zero filled
+ * array to all algorithms that use arrays as iterators (STL style).
+ * It's useful when there's a need to compute the distance between feature
+ * and origin it and allows for better compiler optimisation than using a
+ * zero-filled array.
+ */
+template <typename T>
+struct ZeroIterator
+{
+
+ T operator*()
+ {
+ return 0;
+ }
+
+ T operator[](int)
+ {
+ return 0;
+ }
+
+ const ZeroIterator<T>& operator ++()
+ {
+ return *this;
+ }
+
+ ZeroIterator<T> operator ++(int)
+ {
+ return *this;
+ }
+
+ ZeroIterator<T>& operator+=(int)
+ {
+ return *this;
+ }
+
+};
+
+
+
/**
* Squared Euclidean distance functor.
*
typedef T ElementType;
typedef typename Accumulator<T>::Type ResultType;
+ typedef ResultType CentersType;
template <typename Iterator1, typename Iterator2>
ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
typedef T ElementType;
typedef typename Accumulator<T>::Type ResultType;
+ typedef ResultType CentersType;
/**
* Compute the squared Euclidean distance between two vectors.
typedef T ElementType;
typedef typename Accumulator<T>::Type ResultType;
+ typedef ResultType CentersType;
/**
* Compute the Manhattan (L_1) distance between two vectors.
typedef T ElementType;
typedef typename Accumulator<T>::Type ResultType;
+ typedef ResultType CentersType;
int order;
typedef T ElementType;
typedef typename Accumulator<T>::Type ResultType;
+ typedef ResultType CentersType;
/**
* Compute the max distance (L_infinity) between two vectors.
typedef unsigned char ElementType;
typedef int ResultType;
+ typedef ElementType CentersType;
/** this will count the bits in a ^ b
*/
- ResultType operator()(const unsigned char* a, const unsigned char* b, size_t size) const
+ template<typename Iterator2>
+ ResultType operator()(const unsigned char* a, const Iterator2 b, size_t size) const
{
static const uchar popCountTable[] =
{
3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
};
ResultType result = 0;
+ const unsigned char* b2 = reinterpret_cast<const unsigned char*> (b);
for (size_t i = 0; i < size; i++) {
- result += popCountTable[a[i] ^ b[i]];
+ result += popCountTable[a[i] ^ b2[i]];
+ }
+ return result;
+ }
+
+
+ ResultType operator()(const unsigned char* a, const ZeroIterator<unsigned char> b, size_t size) const
+ {
+ (void)b;
+ static const uchar popCountTable[] =
+ {
+ 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
+ 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
+ 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
+ 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
+ 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
+ 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
+ 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
+ 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
+ };
+ ResultType result = 0;
+ for (size_t i = 0; i < size; i++) {
+ result += popCountTable[a[i]];
}
return result;
}
typedef T ElementType;
typedef int ResultType;
+ typedef ElementType CentersType;
template<typename Iterator1, typename Iterator2>
- ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
+ ResultType operator()(const Iterator1 a, const Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
{
ResultType result = 0;
#if defined(__ARM_NEON__) && !defined(__CUDACC__)
{
+ const unsigned char* a2 = reinterpret_cast<const unsigned char*> (a);
+ const unsigned char* b2 = reinterpret_cast<const unsigned char*> (b);
uint32x4_t bits = vmovq_n_u32(0);
for (size_t i = 0; i < size; i += 16) {
- uint8x16_t A_vec = vld1q_u8 (a + i);
- uint8x16_t B_vec = vld1q_u8 (b + i);
+ uint8x16_t A_vec = vld1q_u8 (a2 + i);
+ uint8x16_t B_vec = vld1q_u8 (b2 + i);
uint8x16_t AxorB = veorq_u8 (A_vec, B_vec);
uint8x16_t bitsSet = vcntq_u8 (AxorB);
uint16x8_t bitSet8 = vpaddlq_u8 (bitsSet);
#endif
return result;
}
+
+
+ template<typename Iterator1>
+ ResultType operator()(const Iterator1 a, ZeroIterator<unsigned char> b, size_t size, ResultType /*worst_dist*/ = -1) const
+ {
+ (void)b;
+ ResultType result = 0;
+#if defined(__ARM_NEON__) && !defined(__CUDACC__)
+ {
+ const unsigned char* a2 = reinterpret_cast<const unsigned char*> (a);
+ uint32x4_t bits = vmovq_n_u32(0);
+ for (size_t i = 0; i < size; i += 16) {
+ uint8x16_t A_vec = vld1q_u8 (a2 + i);
+ uint8x16_t bitsSet = vcntq_u8 (A_vec);
+ uint16x8_t bitSet8 = vpaddlq_u8 (bitsSet);
+ uint32x4_t bitSet4 = vpaddlq_u16 (bitSet8);
+ bits = vaddq_u32(bits, bitSet4);
+ }
+ uint64x2_t bitSet2 = vpaddlq_u32 (bits);
+ result = vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),0);
+ result += vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),2);
+ }
+#elif defined(__GNUC__)
+ {
+ //for portability just use unsigned long -- and use the __builtin_popcountll (see docs for __builtin_popcountll)
+ typedef unsigned long long pop_t;
+ const size_t modulo = size % sizeof(pop_t);
+ const pop_t* a2 = reinterpret_cast<const pop_t*> (a);
+ const pop_t* a2_end = a2 + (size / sizeof(pop_t));
+
+ for (; a2 != a2_end; ++a2) result += __builtin_popcountll(*a2);
+
+ if (modulo) {
+ //in the case where size is not dividable by sizeof(size_t)
+ //need to mask off the bits at the end
+ pop_t a_final = 0;
+ memcpy(&a_final, a2, modulo);
+ result += __builtin_popcountll(a_final);
+ }
+ }
+#else // NO NEON and NOT GNUC
+ HammingLUT lut;
+ result = lut(reinterpret_cast<const unsigned char*> (a), b, size);
+#endif
+ return result;
+ }
};
template<typename T>
typedef T ElementType;
typedef int ResultType;
+ typedef ElementType CentersType;
/** This is popcount_3() from:
* http://en.wikipedia.org/wiki/Hamming_weight */
#endif
template <typename Iterator1, typename Iterator2>
- ResultType operator()(Iterator1 a, Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
+ ResultType operator()(const Iterator1 a, const Iterator2 b, size_t size, ResultType /*worst_dist*/ = -1) const
{
#ifdef FLANN_PLATFORM_64_BIT
const uint64_t* pa = reinterpret_cast<const uint64_t*>(a);
return result;
}
+
+ template <typename Iterator1>
+ ResultType operator()(const Iterator1 a, ZeroIterator<unsigned char> b, size_t size, ResultType /*worst_dist*/ = -1) const
+ {
+ (void)b;
+#ifdef FLANN_PLATFORM_64_BIT
+ const uint64_t* pa = reinterpret_cast<const uint64_t*>(a);
+ ResultType result = 0;
+ size /= long_word_size_;
+ for(size_t i = 0; i < size; ++i ) {
+ result += popcnt64(*pa);
+ ++pa;
+ }
+#else
+ const uint32_t* pa = reinterpret_cast<const uint32_t*>(a);
+ ResultType result = 0;
+ size /= long_word_size_;
+ for(size_t i = 0; i < size; ++i ) {
+ result += popcnt32(*pa);
+ ++pa;
+ }
+#endif
+ return result;
+ }
+
private:
#ifdef FLANN_PLATFORM_64_BIT
static const size_t long_word_size_ = sizeof(uint64_t)/sizeof(unsigned char);
typedef T ElementType;
typedef typename Accumulator<T>::Type ResultType;
+ typedef ResultType CentersType;
/**
* Compute the histogram intersection distance
typedef T ElementType;
typedef typename Accumulator<T>::Type ResultType;
+ typedef ResultType CentersType;
/**
* Compute the Hellinger distance
typedef T ElementType;
typedef typename Accumulator<T>::Type ResultType;
+ typedef ResultType CentersType;
/**
* Compute the chi-square distance
typedef T ElementType;
typedef typename Accumulator<T>::Type ResultType;
+ typedef ResultType CentersType;
/**
* Compute the Kullback-Leibler divergence
};
-
-/*
- * This is a "zero iterator". It basically behaves like a zero filled
- * array to all algorithms that use arrays as iterators (STL style).
- * It's useful when there's a need to compute the distance between feature
- * and origin it and allows for better compiler optimisation than using a
- * zero-filled array.
- */
-template <typename T>
-struct ZeroIterator
-{
-
- T operator*()
- {
- return 0;
- }
-
- T operator[](int)
- {
- return 0;
- }
-
- const ZeroIterator<T>& operator ++()
- {
- return *this;
- }
-
- ZeroIterator<T> operator ++(int)
- {
- return *this;
- }
-
- ZeroIterator<T>& operator+=(int)
- {
- return *this;
- }
-
-};
-
-
/*
* Depending on processed distances, some of them are already squared (e.g. L2)
* and some are not (e.g.Hamming). In KMeans++ for instance we want to be sure
/*
+ * ...a template to tell the user if the distance he is working with is actually squared
+ */
+
+template <typename Distance, typename ElementType>
+struct isSquareDist
+{
+ bool operator()() { return false; }
+};
+
+
+template <typename ElementType>
+struct isSquareDist<L2_Simple<ElementType>, ElementType>
+{
+ bool operator()() { return true; }
+};
+
+template <typename ElementType>
+struct isSquareDist<L2<ElementType>, ElementType>
+{
+ bool operator()() { return true; }
+};
+
+
+template <typename ElementType>
+struct isSquareDist<MinkowskiDistance<ElementType>, ElementType>
+{
+ bool operator()() { return true; }
+};
+
+template <typename ElementType>
+struct isSquareDist<HellingerDistance<ElementType>, ElementType>
+{
+ bool operator()() { return true; }
+};
+
+template <typename ElementType>
+struct isSquareDist<ChiSquareDistance<ElementType>, ElementType>
+{
+ bool operator()() { return true; }
+};
+
+
+template <typename Distance>
+bool isSquareDistance()
+{
+ typedef typename Distance::ElementType ElementType;
+
+ isSquareDist<Distance, ElementType> dummy;
+ return dummy();
+}
+
+/*
* ...and a template to ensure the user that he will process the normal distance,
* and not squared distance, without losing processing time calling sqrt(ensureSquareDistance)
* that will result in doing actually sqrt(dist*dist) for L1 distance for instance.
* of the form (branching-1)*K+1 smaller than clusters.rows).
*/
template <typename Distance>
-int hierarchicalClustering(const Matrix<typename Distance::ElementType>& points, Matrix<typename Distance::ResultType>& centers,
+int hierarchicalClustering(const Matrix<typename Distance::ElementType>& points, Matrix<typename Distance::CentersType>& centers,
const KMeansIndexParams& params, Distance d = Distance())
{
KMeansIndex<Distance> kmeans(points, params, d);
#include "saving.h"
#include "logger.h"
+#define BITS_PER_CHAR 8
+
namespace cvflann
{
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;
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 Matrix<double>& _dcenters, const size_t _veclen,
- std::vector<int> &_new_centroids, std::vector<DistanceType> &_sq_dists)
+ 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)
const Matrix<ElementType>& dataset;
const int branching;
const int* indices;
- const Matrix<double>& dcenters;
+ const CentersContainerType& dcenters;
const size_t veclen;
std::vector<int> &new_centroids;
std::vector<DistanceType> &sq_dists;
root_ = pool_.allocate<KMeansNode>();
std::memset(root_, 0, sizeof(KMeansNode));
- computeNodeStatistics(root_, indices_, (int)size_);
- computeClustering(root_, indices_, (int)size_, branching_,0);
+ if(is_kdtree_distance::val || is_vector_space_distance::val)
+ {
+ computeNodeStatistics(root_, indices_, (unsigned int)size_);
+ computeClustering(root_, indices_, (int)size_, branching_,0);
+ }
+ else
+ {
+ computeBitfieldNodeStatistics(root_, indices_, (unsigned int)size_);
+ computeBitfieldClustering(root_, indices_, (int)size_, branching_,0);
+ }
}
* numClusters = number of clusters to have in the clustering computed
* Returns: number of cluster centers
*/
- int getClusterCenters(Matrix<DistanceType>& centers)
+ int getClusterCenters(Matrix<CentersType>& centers)
{
int numClusters = centers.rows;
if (numClusters<1) {
Logger::info("Clusters requested: %d, returning %d\n",numClusters, clusterCount);
for (int i=0; i<clusterCount; ++i) {
- DistanceType* center = clusters[i]->pivot;
+ CentersType* center = clusters[i]->pivot;
for (size_t j=0; j<veclen_; ++j) {
centers[i][j] = center[j];
}
/**
* The cluster center.
*/
- DistanceType* pivot;
+ CentersType* pivot;
/**
* The cluster radius.
*/
{
node = pool_.allocate<KMeansNode>();
load_value(stream, *node);
- node->pivot = new DistanceType[veclen_];
+ node->pivot = new CentersType[veclen_];
load_value(stream, *(node->pivot), (int)veclen_);
if (node->childs==NULL) {
int indices_offset;
* indices = array of indices of the points belonging to the node
* indices_length = number of indices in the array
*/
- void computeNodeStatistics(KMeansNodePtr node, int* indices, int indices_length)
+ void computeNodeStatistics(KMeansNodePtr node, int* indices, unsigned int indices_length)
{
-
- DistanceType radius = 0;
DistanceType variance = 0;
- DistanceType* mean = new DistanceType[veclen_];
- memoryCounter_ += int(veclen_*sizeof(DistanceType));
+ CentersType* mean = new CentersType[veclen_];
+ memoryCounter_ += int(veclen_*sizeof(CentersType));
- memset(mean,0,veclen_*sizeof(DistanceType));
+ memset(mean,0,veclen_*sizeof(CentersType));
- for (int i=0; i<indices_length; ++i) {
+ 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] /= indices_length;
+ mean[j] = cvflann::round<CentersType>( mean[j] / static_cast<double>(indices_length) );
}
- variance /= indices_length;
+ variance /= static_cast<DistanceType>( length );
variance -= distance_(mean, ZeroIterator<ElementType>(), veclen_);
- DistanceType tmp = 0;
- for (int i=0; i<indices_length; ++i) {
- tmp = distance_(mean, dataset_[indices[i]], 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;
}
}
+ 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, 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;
+ }
+
+
+
/**
* The method responsible with actually doing the recursive hierarchical
* clustering
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) {
std::vector<DistanceType> sq_dists(indices_length);
// reassign points to clusters
- KMeansDistanceComputer invoker(distance_, dataset_, branching, indices, dcenters, veclen_, new_centroids, sq_dists);
+ 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** centers = new DistanceType*[branching];
+ CentersType** centers = new CentersType*[branching];
for (int i=0; i<branching; ++i) {
- centers[i] = new DistanceType[veclen_];
- memoryCounter_ += (int)(veclen_*sizeof(DistanceType));
+ centers[i] = new CentersType[veclen_];
+ memoryCounter_ += (int)(veclen_*sizeof(CentersType));
for (size_t k=0; k<veclen_; ++k) {
- centers[i][k] = (DistanceType)dcenters[i][k];
+ centers[i][k] = (CentersType)dcenters[i][k];
}
}
if (belongs_to[i]==c) {
DistanceType d = distance_(dataset_[indices[i]], ZeroIterator<ElementType>(), veclen_);
variance += d;
- mean_radius += sqrt(d);
+ mean_radius += static_cast<DistanceType>( sqrt(d) );
std::swap(indices[i],indices[end]);
std::swap(belongs_to[i],belongs_to[end]);
end++;
+ void computeBitfieldClustering(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;
+ }
+
+ 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);
+
+ CentersType** centers = new CentersType*[branching];
+
+ for (int i=0; i<branching; ++i) {
+ centers[i] = new CentersType[veclen_];
+ memoryCounter_ += (int)(veclen_*sizeof(CentersType));
+ }
+
+ 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 dist = distance_(dataset_[indices[i]], dataset_[centers_idx[0]], veclen_);
+ belongs_to[i] = 0;
+ for (int j=1; j<branching; ++j) {
+ DistanceType new_dist = distance_(dataset_[indices[i]], dataset_[centers_idx[j]], veclen_);
+ if (dist>new_dist) {
+ belongs_to[i] = j;
+ dist = new_dist;
+ }
+ }
+ if (dist>radiuses[belongs_to[i]]) {
+ radiuses[belongs_to[i]] = dist;
+ }
+ count[belongs_to[i]]++;
+ }
+
+ 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;
+ }
+ }
+
+ }
+
+
+ // 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;
+ computeBitfieldClustering(node->childs[c],indices+start, end-start, branching, level+1);
+ start=end;
+ }
+
+ delete[] centers;
+ }
+
+
+
+
/**
* Performs one descent in the hierarchical k-means tree. The branches not
* visited are stored in a priority queue.
DistanceType rsq = node->radius;
DistanceType wsq = result.worstDist();
- DistanceType val = bsq-rsq-wsq;
- DistanceType val2 = val*val-4*rsq*wsq;
-
- //if (val>0) {
- if ((val>0)&&(val2>0)) {
- return;
+ if (isSquareDistance<Distance>())
+ {
+ DistanceType val = bsq-rsq-wsq;
+ if ((val>0) && (val*val > 4*rsq*wsq))
+ return;
+ }
+ else
+ {
+ if (bsq-rsq > wsq)
+ return;
}
}
// float* best_center = node->childs[best_index]->pivot;
for (int i=0; i<branching_; ++i) {
if (i != best_index) {
- domain_distances[i] -= cb_index_*node->childs[i]->variance;
+ 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) {
DistanceType rsq = node->radius;
DistanceType wsq = result.worstDist();
- DistanceType val = bsq-rsq-wsq;
- DistanceType val2 = val*val-4*rsq*wsq;
-
- // if (val>0) {
- if ((val>0)&&(val2>0)) {
- return;
+ if (isSquareDistance<Distance>())
+ {
+ DistanceType val = bsq-rsq-wsq;
+ if ((val>0) && (val*val > 4*rsq*wsq))
+ return;
+ }
+ else
+ {
+ if (bsq-rsq > wsq)
+ return;
}
}
DistanceType dist = distance_(q, node->childs[i]->pivot, veclen_);
int j=0;
- while (domain_distances[j]<dist && j<i) j++;
+ 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];
projects / platform / upstream / opencv.git / commitdiff