K均值

图像分割中，“k均值k类，并满足数据点到对应的均值点之间的距离之和最小（类内方差之和最小）。实现步骤如下：

k

1)将所有样本的均值作为第一类的聚类中心；3) 判断是否选够2)，否则退出。

4. 通过迭代选择最佳聚类中心来最小化类内方差。每次迭代时，取出一个数据，将它移动到一个适当的类别中，使得类内方差之和最小。直到迭代<FONT face=""">N（观测数据的个数）次之后，类内方差之后都不变为止才退出。

 

以下代码可以处理多维的数据（即每个数据点用一个向量来表示，比如RGB图像）

.h头文件

#ifndef __KMEANS_H__ #define __KMEANS_H__ class KMeans { public: typedef unsigned char Pixel; typedef unsigned char Label; typedef struct ParamStruct { double m_nb; double *m_average; }Param; KMeans(Pixel *pixel, Label *answer, int nLabel, int nPixel, int dim = 1); ~KMeans(); inline Param* ParamPtr(){return m_param;} inline void ClrParam() { memset(m_param->m_average, 0, m_nLabel * m_dim * sizeof(double)); int i = 0; while (i<m_nLabel) { (m_param+i)->m_nb = 0; i++; } } inline Pixel* GetPixel(int pix) const { return m_pixel + pix * m_dim; } void AllPixelAverage(); // compute the average of all pixel inline double Dist(Pixel* p, double* d) const; inline double ShortestDist(Pixel* p, int nSeed) const; void Seeds(); void InitAnswer(); void ClusterAverage(); void WCSS(); void Iteration(); inline double Rho(int pix, Label lab) const; inline void UpdateAverageIn(int pix, Label lab); inline void UpdateAverageOut(int pix, Label lab); inline void UpdateWCSS(int pix, Label oldLab, Label newLab); void Cluster(); void Sort(); protected: void Initial(); // allocate memory for parameters. protected: int m_nPixel; int m_nLabel; int m_dim; Pixel* m_pixel; Label* m_answer; double m_wcss; // within-cluster sum of square private: Param *m_param; bool m_initialized; }; inline double KMeans::Rho(int pix, Label lab) const { Label oldLab = *(m_answer+pix); double* d = (m_param + lab)->m_average; Pixel* p = GetPixel(pix); if (oldLab != lab) { return (m_param+lab)->m_nb * Dist(p, d) / ((m_param+lab)->m_nb + 1); } return (m_param+lab)->m_nb * Dist(p, d) / ((m_param+lab)->m_nb - 1); } inline void KMeans::UpdateAverageIn(int pix, Label lab) { *(m_answer+pix) = lab; // move Pixel* p = GetPixel(pix); double* d = (m_param + lab)->m_average; double nb = 1 + (m_param + lab)->m_nb++; int k; for (k=0; k<m_dim; k++) { *(d+k) += (*(p+k) - *(d+k)) / nb; } } inline void KMeans::UpdateAverageOut(int pix, Label lab) { Pixel* p = GetPixel(pix); double* d = (m_param + lab)->m_average; double nb = -1 + (m_param + lab)->m_nb--; int k; for (k=0; k<m_dim; k++) { *(d + k) += (*(d+k) - *(p+k)) / nb; } } inline void KMeans::UpdateWCSS(int pix, Label oldLab, Label newLab) { Pixel* p = GetPixel(pix); double* d = (m_param+oldLab)->m_average; double nb = (m_param+oldLab)->m_nb; m_wcss -= (nb+1) * Dist(p, d) / nb; d = (m_param + newLab)->m_average; nb = (m_param + newLab)->m_nb; m_wcss += (nb-1) * Dist(p, d) / nb; } #endif

.cpp 源文件

#include <string.h> #include <stdio.h> #include "KMeans.h" #include <iostream> using namespace std; // // constructor KMeans::KMeans(Pixel *pixel, Label *answer, int nLabel, int nPixel, int dim) { m_pixel = pixel; m_answer = answer; m_nLabel = nLabel; m_nPixel = nPixel; m_dim = dim; m_initialized = 0; m_param = 0; Initial(); } // // destructor KMeans::~KMeans() { if (m_param) { delete[] m_param; } } // // 为参数分配空间 void KMeans::Initial() { m_param = (Param*) new char[m_nLabel * sizeof(Param) + m_nLabel * m_dim * sizeof(double)]; if (!m_param) { printf("/nNot enough memory (kmeans initial), exiting"); } memset(m_param, 0, m_nLabel * sizeof(Param) + m_nLabel * m_dim * sizeof(double)); double *startData = (double*)((char*)m_param + m_nLabel * sizeof(Param)); int k; for (k=0; k<m_nLabel; k++) { (m_param+k)->m_nb = 0; (m_param+k)->m_average = startData + k * m_dim; } m_initialized = 1; } // // 计算全部数据的均值，并放入第一个Param中 void KMeans::AllPixelAverage() { if (!m_initialized) { exit(0); } double* d = m_param->m_average; Pixel* p = 0; int i, k; for (i=0; i<m_nPixel; i++) { p = GetPixel(i); for (k=0; k<m_dim; k++) { *(d+k) += *(p+k); } } for (k=0; k<m_dim; k++) { *(d+k) /= m_nPixel; } } // // 计算给定像素与给定聚类中心的距离 inline double KMeans::Dist(Pixel* p, double* d) const { double dist = 0.0; int k; for (k=0; k<m_dim; k++) { dist += (*(d+k) - *(p+k)) * (*(d+k) - *(p+k)); } return dist; } // // 计算给定像素的最短距离 inline double KMeans::ShortestDist(Pixel* p, int nSeed) const { double* d = m_param->m_average; double shortestDist = Dist(p, d); double dist = 0.0; int count = 1; for (; count < nSeed; count++) { dist = Dist(p, (m_param+count)->m_average); if (dist < shortestDist) { shortestDist = dist; } } return shortestDist; } // // 选择种子点，作为初始聚类中心。具体方法如下： // 对每个观测点，它都有一个“距离”，这个距离是它到所有种子点中最小的 // 选择“距离”最大的点作为新的种子点 // 直到选够所有的种子点 void KMeans::Seeds() { if (!m_initialized) { exit(1); } AllPixelAverage(); int nSeeds = 1; int newSeed = 0; double shortestDist = 0.0; double maxShortestDist = 0.0; int i, k; double* d; Pixel* p; while(nSeeds < m_nLabel) { newSeed = 0; maxShortestDist = ShortestDist(GetPixel(newSeed), nSeeds); for (i=1; i<m_nPixel; i++) { shortestDist = ShortestDist(GetPixel(i), nSeeds); if (shortestDist > maxShortestDist) { maxShortestDist = shortestDist; newSeed = i; } } d = (m_param + nSeeds)->m_average; p = GetPixel(newSeed); for (k=0; k<m_dim; k++) { *(d+k) = *(p+k); } nSeeds++; } } // //根据选取的种子点，初始化标号场 void KMeans::InitAnswer() { int i, k; double shortestDist = 0.0; double dist = 0.0; Label l; Pixel* p; for (i=0; i<m_nPixel; i++) { l = 0; p = GetPixel(i); shortestDist = Dist(p, m_param->m_average); for (k=1; k<m_nLabel; k++) { dist = Dist(p, (m_param + k)->m_average); if (dist < shortestDist) { shortestDist = dist; l = k; } } *(m_answer+i) = l; } } // // 计算每类的均值 void KMeans::ClusterAverage() { ClrParam(); double* d = 0; Pixel* p = 0; Label l = 0; int i, k; for (i=0; i<m_nPixel; i++) { p = GetPixel(i); l = *(m_answer+i); d = (m_param + l)->m_average; (m_param + l)->m_nb++; for (k=0; k<m_dim; k++) { *(d+k) += *(p+k); } } int nb; for (i=0; i<m_nLabel; i++) { d = (m_param + i)->m_average; nb = (m_param + i)->m_nb; for (k=0; k<m_dim; k++) { *(d+k) /= nb; } } } // // 计算类内方差之和 void KMeans::WCSS() { m_wcss = 0.0; int i=0; double* d; Pixel* p; for (; i<m_nPixel; i++) { d = (m_param + *(m_answer + i))->m_average; p = GetPixel(i); m_wcss += Dist(p, d); } } // // 迭代 void KMeans::Iteration() { int i, k; double oldWCSS = m_wcss; double rho, minRho; Label oldLab, newLab; while (1) { for (i=0; i<m_nPixel; i++) { oldLab = *(m_answer+i); if ((m_param+oldLab)->m_nb <= 1) { continue; } newLab = 0; minRho = Rho(i, newLab); for (k = 1; k<m_nLabel; k++) { rho = Rho(i, k); if (rho < minRho) { newLab = k; minRho = rho; } } if (newLab != oldLab) { UpdateAverageIn(i, newLab); UpdateAverageOut(i, oldLab); UpdateWCSS(i,oldLab, newLab); } } if (oldWCSS >= m_wcss) { break; } else { oldWCSS = m_wcss; } } } // // 排序 void KMeans::Sort() { double* val = new double[m_nLabel]; unsigned char* ord = new unsigned char[m_nLabel]; memset(val, 0, sizeof(double) * m_nLabel); memset(ord, 0, sizeof(unsigned char) * m_nLabel); double* d; int i = 0; int k; for (; i<m_nLabel; i++) { d = (m_param + i)->m_average; for (k=0; k<m_dim; k++) { *(val+i) += (*(d+k)) * (*(d+k)); } } int const MAX = 256 * 256 * 256; double minVal = 0.0; unsigned char l; for (i=0; i<m_nLabel; i++) { minVal = *val; l = 0; for (k=1; k<m_nLabel; k++) { if (*(val+k) < minVal) { minVal = *(val+k); l = k; } } *(ord + i) = l; *(val+l) = MAX; } for (i=0; i<m_nPixel; i++) { l = *(m_answer + i); for (k=0; k<m_nLabel; k++) { if (l == *(ord+k)) { *(m_answer + i) = k; } } } ClusterAverage(); delete[] val; delete[] ord; } // // 聚类 void KMeans::Cluster() { Seeds(); InitAnswer(); ClusterAverage(); Iteration(); Sort(); }

测试文件

#include <iostream> #include <vld.h> using namespace std; #include "KMeans.h" int const dim = 1; int const nPixel = 256 * 256; // 数据的总大小为 nPixel * dim int const nLabel = 3; int const grayLevel = 256; KMeans::Pixel* pixel; KMeans::Label* answer; void SimulateTestData() { pixel = new KMeans::Pixel[nPixel * dim]; answer = new KMeans::Pixel[nPixel]; if (!pixel || !answer) { cout << "/nNot enough memory, exiting" << endl; exit(0); } int i; for (i=0; i<nPixel * dim; i++) { *(pixel+i) = (KMeans::Pixel)(rand() % grayLevel); } memset(answer, 0, sizeof(KMeans::Label) * nPixel); } int main() { int const seed = 1000; srand(seed); SimulateTestData(); KMeans *kmeans = new KMeans(pixel, answer, nLabel, nPixel, dim); kmeans->Cluster(); KMeans::Param *param = kmeans->ParamPtr(); if (param) { int i, k; for (k=0; k<nLabel; k++) { cout <<"/nnum = "<<(param+k)->m_nb<< ",/taverage = "; for (i=0; i<dim; i++) { cout << *((param+k)->m_average + i) << "/t"; } } } cout << endl; // int i; // for (i=0; i<nPixel; i++) // { // cout << (int)*(answer+i) << "/t"; // } // cout << endl; delete[] pixel; delete[] answer; delete kmeans; return 0; }

以后会在我的资源中添加这些文件。