Image Component Library (ICL)
|
00001 /******************************************************************** 00002 ** Image Component Library (ICL) ** 00003 ** ** 00004 ** Copyright (C) 2006-2013 CITEC, University of Bielefeld ** 00005 ** Neuroinformatics Group ** 00006 ** Website: www.iclcv.org and ** 00007 ** http://opensource.cit-ec.de/projects/icl ** 00008 ** ** 00009 ** File : ICLMath/src/ICLMath/KMeans.h ** 00010 ** Module : ICLMath ** 00011 ** Authors: Christof Elbrechter ** 00012 ** ** 00013 ** ** 00014 ** GNU LESSER GENERAL PUBLIC LICENSE ** 00015 ** This file may be used under the terms of the GNU Lesser General ** 00016 ** Public License version 3.0 as published by the ** 00017 ** ** 00018 ** Free Software Foundation and appearing in the file LICENSE.LGPL ** 00019 ** included in the packaging of this file. Please review the ** 00020 ** following information to ensure the license requirements will ** 00021 ** be met: http://www.gnu.org/licenses/lgpl-3.0.txt ** 00022 ** ** 00023 ** The development of this software was supported by the ** 00024 ** Excellence Cluster EXC 277 Cognitive Interaction Technology. ** 00025 ** The Excellence Cluster EXC 277 is a grant of the Deutsche ** 00026 ** Forschungsgemeinschaft (DFG) in the context of the German ** 00027 ** Excellence Initiative. ** 00028 ** ** 00029 ********************************************************************/ 00030 00031 #pragma once 00032 00033 #include <ICLUtils/CompatMacros.h> 00034 #include <ICLUtils/Random.h> 00035 #include <ICLUtils/Point32f.h> 00036 #include <algorithm> 00037 00038 namespace icl{ 00039 namespace math{ 00040 00041 00043 00057 template<class Vector, class Scalar> 00058 class KMeans{ 00059 00061 std::vector<Vector> m_centers; 00062 00064 std::vector<Vector> m_means; 00065 00067 std::vector<int> m_nums; 00068 00070 std::vector<Scalar> m_errors; 00071 00073 static Scalar diff_power_two(const Scalar &a, const Scalar &b){ 00074 Scalar d = a-b; 00075 return d*d; 00076 } 00077 00079 static inline Scalar dist(const Vector &a, const Vector &b){ 00080 return ::sqrt( std::inner_product(a.begin(),a.end(),b.begin(),Scalar(0), std::plus<Scalar>(), diff_power_two) ); 00081 } 00082 00084 static inline void setVectorNull(Vector &v){ 00085 std::fill(v.begin(),v.end(),Scalar(0)); 00086 } 00087 00088 public: 00089 00091 struct ICLMath_API Result{ 00092 const std::vector<Vector> ¢ers; 00093 const std::vector<int> &nums; 00094 const std::vector<Scalar> &errors; 00095 }; 00096 00098 inline KMeans(int numCenters=0){ 00099 init(numCenters); 00100 } 00101 00103 inline void init(int numCenters){ 00104 m_centers.resize(numCenters); 00105 m_means.resize(numCenters); 00106 m_nums.resize(numCenters); 00107 m_errors.resize(numCenters); 00108 } 00109 00111 inline int findNN(const Vector &v, Scalar &minDist){ 00112 int minIdx = 0; 00113 minDist = dist(v,m_centers[0]); 00114 for(size_t i=1;i<m_centers.size();++i){ 00115 Scalar currDist = dist(v,m_centers[i]); 00116 if(currDist < minDist){ 00117 minDist = currDist; 00118 minIdx = i; 00119 } 00120 } 00121 return minIdx; 00122 } 00123 00125 00127 template<class RandomAcessIterator> 00128 Result apply(RandomAcessIterator begin, RandomAcessIterator end, int numSteps = 1000, bool reinitCenters=true){ 00129 Scalar minDist = 0; 00130 00131 if(reinitCenters){ 00132 URandI r((int)(end-begin)-1); 00133 for(size_t i=0;i<m_centers.size();++i){ 00134 int ri = r; 00135 m_centers[i] = *(begin+ri); 00136 } 00137 } 00138 00139 for(int step=0;step<numSteps;++step){ 00140 // empty accumulators 00141 for(size_t i=0;i<m_means.size();++i){ 00142 setVectorNull(m_means[i]); 00143 m_nums[i] = Scalar(0); 00144 m_errors[i] = Scalar(0); 00145 } 00146 00147 // estimate means 00148 for(RandomAcessIterator it=begin;it != end; ++it){ 00149 const int nn = findNN(*it,minDist); 00150 m_means[nn] += *it; 00151 m_nums[nn] ++; 00152 m_errors[nn] += minDist; 00153 } 00154 00155 for(size_t i=0;i<m_means.size();++i){ 00156 if(m_nums[i]){ 00157 m_centers[i] = m_means[i] * (Scalar(1.0)/m_nums[i]); 00158 m_errors[i] /= m_nums[i]; 00159 }// empty voronoi tessels are not moved 00160 } 00161 } 00162 00163 Result r= { m_centers, m_nums, m_errors }; 00164 return r; 00165 } 00166 }; 00167 00169 template<> 00170 float KMeans<utils::Point32f, float>::dist(const utils::Point32f &a, const utils::Point32f &b){ 00171 return a.distanceTo(b); 00172 } 00173 00174 template<> 00175 void KMeans<utils::Point32f, float>::setVectorNull(utils::Point32f &p){ 00176 p = utils::Point32f::null; 00177 } 00181 } 00182 }