#define ef else if template <typename T> inline void set_vector (T* out, const T* in, int n) { while (n-- > 0) out[n] = in[n]; } template <typename T> inline void add_vector (T* out, const T* in, int n) { while (n-- > 0) out[n] += in[n]; } template <typename T> inline void div_vector (T* v, int divisor, int n) { while (n-- > 0) v[n] /= divisor; } template <typename T> inline double cal_distance (T* a, T* b, int n) { double sum = 0, diff; while (n-- > 0) { diff = a[n] - b[n]; sum += diff * diff; } return sqrt (sum); } template <typename T = float> struct Kmeans { int nMaxIteration; //最大迭代次數 int cur_iter; //目前的迭代次數 int nData; //資料項數 int* label; //資料項的分類標籤 int* ccount; //群心涵蓋的資料個數 T** data; //指向資料集合 T** center; //群心集合 T* radius; //群心的超球面半徑 T* new_center; //計算暫用的群心 T* mean; //資料項的幾何中心 int nCenter; //群心個數 int dim; //資料維度 double cluster_movement; //群心更新時的最大移動量 Kmeans() { srand ((unsigned int)time(0)); label = 0; ccount = 0; radius = 0; center = 0; mean = 0; new_center = 0; } ~Kmeans() { release(); } protected: void set (T* out, const T* in) {set_vector (out, in, dim);} void add (T* out, const T* in) {add_vector (out, in, dim);} void div (T* v, int divisor) {div_vector (v, divisor, dim);} double distance (T* a, T* b) {return cal_distance (a, b, dim);} void release() //釋放資源 { if (label) del_arr (label); label = 0; if (ccount) del_arr (ccount); ccount = 0; if (radius) del_arr (radius); radius = 0; if (center) del_arr (center); center = 0; if (mean) del_arr (mean); mean = 0; if (new_center) del_arr (new_center); new_center = 0; } void init_centers () { if (nCenter < 1) ERROR_MSG_("群心數目需 > 1"); if (nCenter > nData) nCenter = nData; release(); //釋放資源 new_center = (T*) new_arr (sizeof(T), 1, dim); radius = (T*) new_arr (sizeof(T), 1, nCenter); mean = (T*) new_arr (sizeof(T), 1, dim); center = (T**) new_arr (sizeof(T), 2, nCenter, dim); ccount = (int*) new_arr (sizeof(int), 1, nCenter); label = (int*) new_arr (sizeof(int), 1, nData); int i; for (i=0; i<nCenter; ++i) //隨機指定樣本點為群心 set (center[i], data [rand() % nData]); for (i=0; i<nData; ++i) //計算資料項的幾何中心 add (mean, data[i]); div (mean, nData); } void cluster_assignment () //將樣本分群 { int label_id, j, k; double min_dist, d; memset (ccount, 0, nCenter* sizeof(int)); memset (radius, 0, nCenter* sizeof(T)); for (j=0; j<nData; ++j) { label_id = 0; min_dist = distance (center[0], data[j]); for (k=1; k<nCenter; ++k) { d = distance (center[k], data[j]); if (d < min_dist) { label_id = k; min_dist = d; } } if (radius[label_id] < min_dist) radius[label_id] = min_dist; label[j] = label_id; //歸類 ccount[label_id]++; //統計群心涵蓋的資料個數 } } void update_cluster () //重新估算 centroids 位置 { int j, k, total; double dist = 0; cluster_movement = 0; for (k=0; k<nCenter; ++k) { memset (new_center, 0, dim* sizeof(T)); total = 0; for (j=0; j<nData; ++j) if (label[j] == k) { add (new_center, data[j]); total++; } if (total > 0) { div (new_center, total); dist = distance (center[k], new_center); if (cluster_movement < dist) cluster_movement = dist; set (center[k], new_center); } else set (center[k], mean); } } public: //算法主流程: // nCenter = 群心數目, // data = 資料集合 = data [nData_] [dim_] void do_cluster (int nCenter_, T** data_, int nData_, int dim_, int nIteration_, double movement_threshold) { setup (nCenter_, data_, nData_, dim_); nMaxIteration = nIteration_; cur_iter = 0; do { cluster_assignment (); //將樣本分群 update_cluster (); //重新估算 centroids 位置 } while (++cur_iter < nMaxIteration && cluster_movement > movement_threshold); } //單步執行 clustering, 回傳群心最大移動量 void setup (int nCenter_, T** data_, int nData_, int dim_) { nCenter = nCenter_; //設定群心數目 data = data_; //指向資料集合 nData = nData_; //設定資料個數 dim = dim_; //設定資料維度 cur_iter = 0; init_centers (); //初始化群心位置 } double iterater_cluster () { cur_iter++; cluster_assignment (); //將樣本分群 update_cluster (); //重新估算 centroids 位置 return cluster_movement; //回傳群心更新時的最大移動量 } }; //-------------------------------------------------------------- // Modify K-Means //-------------------------------------------------------------- template <typename T = float> struct ModifyKmeans: Kmeans<T> { enum {MAX_SELECT = 32}; //每回合最大的群心挑選個數 T** best_center; //計算暫用的群心集合 T** old_center; //計算暫用的群心集合 T* avg_radius; //群心的平均半徑 bool* bAlive; //群心存活狀態(用於合併時) int nMaxCenter; //最大群心數 int nMaxSelect; //最大的群心挑選個數 int nMinCount; //每個群集至少應包含的資料數 double movement_threshold; //指定的群心移動量門檻 double avg_dist; //每一點到群心的平均距離 CvRNG rng; //亂數種子 ModifyKmeans() { best_center = 0; old_center = 0; avg_radius = 0; bAlive = 0; rng = cvRNG (cvGetTickCount()); } ~ModifyKmeans() {release();} void release() //釋放資源 { if (best_center) del_arr (best_center); if (old_center) del_arr (old_center); if (avg_radius) del_arr (avg_radius); if (bAlive) del_arr (bAlive); old_center = 0; best_center = 0; avg_radius = 0; bAlive = 0; Kmeans<T>::release(); } void setup (int nMaxCenter_, T** data_, int nData_, int dim_) { nMaxCenter = nMaxCenter_; //設定群心數目 data = data_; //指向資料集合 nData = nData_; //設定資料個數 dim = dim_; //設定資料維度 cur_iter = 0; nCenter = 1; //目前輪迴的群心數 nMaxSelect = MAX_SELECT; if (nMaxCenter < 1) ERROR_MSG_("最大群心數目需 > 1"); if (nMaxCenter > nData) nMaxCenter = nData; if (nMaxSelect > nData) ERROR_MSG3_("資料項數(%d)需 > 最大的群心挑選個數(%d)", nData, nMaxSelect); release(); //釋放資源 new_center = (T*) new_arr (sizeof(T), 1, dim); radius = (T*) new_arr (sizeof(T), 1, nMaxCenter); bAlive = (bool*)new_arr (sizeof(bool),1, nMaxCenter); avg_radius = (T*) new_arr (sizeof(T), 1, nMaxCenter); best_center = (T**) new_arr (sizeof(T), 2, nMaxCenter, dim); old_center = (T**) new_arr (sizeof(T), 2, nMaxCenter, dim); mean = (T*) new_arr (sizeof(T), 1, dim); center = (T**) new_arr (sizeof(T), 2, nMaxCenter, dim); ccount = (int*) new_arr (sizeof(int), 1, nMaxCenter); label = (int*) new_arr (sizeof(int), 1, nData); for (int i=0; i<nData; ++i) { //計算資料項的幾何中心 add (mean, data[i]); label[i] = 0; } div (mean, nData); set (center[0], mean); //第1個群心是幾何中心 shuffle (data, nData); //洗牌, 將資料內序打亂 } void cluster_assignment () //將樣本分群 { int label_id, i, j, k; double min_dist, d; memset (ccount, 0, nCenter* sizeof(int)); memset (radius, 0, nCenter* sizeof(T)); for (j=0; j<nData; ++j) { label_id = 0; min_dist = distance (center[0], data[j]); for (k=1; k<nCenter; ++k) { d = distance (center[k], data[j]); if (d < min_dist) { label_id = k; min_dist = d; } } if (radius[label_id] < min_dist) //計算每一群的最大半徑 radius[label_id] = min_dist; label[j] = label_id; //歸類 ccount[label_id]++; //統計群心涵蓋的資料個數 } } void update_cluster () //重新估算 centroids 位置 { int j, k, total; double dist = 0; cluster_movement = 0; for (k=0; k<nCenter; ++k) { memset (new_center, 0, dim* sizeof(T)); total = 0; for (j=0; j<nData; ++j) if (label[j] == k) { add (new_center, data[j]); total++; } if (total > 0) { div (new_center, total); dist = distance (center[k], new_center); if (cluster_movement < dist) cluster_movement = dist; //計算最小的群心移動量值 set (center[k], new_center); } else set (center[k], mean); } } void cal_avg_dist() //計算每個群心的平均半徑 { //以及總平均半徑 int j, k, total; int nResetCenter = 0; double dist_sum = 0; double cur_dist = 0; //用來加總每個群心的平均半徑 for (k=0; k<nCenter; ++k) { total = 0; dist_sum = 0; //計算點集到群心的距離總和 for (j=0; j<nData; ++j) if (label[j] == k) { dist_sum += distance (center[k], data[j]); total++; } if (total > 0) { dist_sum /= total; avg_radius[k] = dist_sum; //記錄群心的平均半徑 cur_dist += dist_sum; //加總每個群心的平均半徑 } else nResetCenter++; } //取得每一點到群心的平均距離 avg_dist = cur_dist / (nCenter - nResetCenter); } void merge_center (int parent, int child) { memset (new_center, 0, dim* sizeof(T)); int total = 0; for (int j=0; j<nData; ++j) if (label[j] == parent || label[j] == child) { add (new_center, data[j]); label[j] = parent; total++; } if (total > 0) { div (new_center, total); set (center[parent], new_center); } } void joint_clusters () //合併群內之群與強連通群 { int i, k, j; double d; for (k=0; k<nCenter; ++k) bAlive[k] = true; for (k=0; k<nCenter; ++k) //合併群內之群 if (bAlive[k]) for (j=0; j<nCenter; ++j) if (bAlive[j] && k != j) { d = distance (center[j], center[k]); if (d + radius[j] < radius[k]) { merge_center (k, j); bAlive[j] = false; } } for (j=k=0; k<nCenter; ++k) //重設群的陣列 if (bAlive[k]) { if (k != j) { set (center[j], center[k]); for (i=0; i<nData; ++i) if (label[i] == k) label[i] = j; } j++; } nCenter = j; } void cluster_reassignment () { double min_dist, d; int label_id, k, j; for (j=0; j<nData; ++j) { label_id = 0; min_dist = distance (center[0], data[j]); min_dist -= avg_radius[0]; for (k=1; k<nCenter; ++k) { d = distance (center[k], data[j]); d -= avg_radius[k]; if (d < min_dist) { label_id = k; min_dist = d; } } if (radius[label_id] < min_dist) //計算每一群的最大半徑 radius[label_id] = min_dist; label[j] = label_id; //歸類 ccount[label_id]++; //統計群心涵蓋的資料個數 } } bool kmeans (int k_centers) { nCenter = k_centers; //設定群心個數 int iter = 0; do { cluster_assignment (); //將樣本分群 //並算出每群的最大半徑 if (ccount[nCenter-1] < nMinCount) //若新進群心的包含點數過少 return false; //便捨棄之, 再選別的來試 update_cluster (); //重新估算 centroids 位置 cal_avg_dist (); //計算每個群心的平均半徑 } while (++iter < nMaxIteration && cluster_movement > movement_threshold); return true; } public: void do_cluster (int nMaxCenter_, //最大群心數 T** data_, int nData_, int dim_, //資料 int nMaxIteration_, //最大迭代次數 double mov_threshold_, //誤差門檻: 判別是否收斂 double min_cluster_ratio) //規定每個群集至少應包含 //nData * ratio 筆資料 { setup (nMaxCenter_, data_, nData_, dim_); set (best_center[0], center[0]); nMinCount = int (nData * min_cluster_ratio); nMaxIteration = nMaxIteration_; //內部 K-means 的最大迭代次數 movement_threshold = mov_threshold_; cur_iter = 0; avg_dist = 0; double old_dist = (std::numeric_limits<double>::max)(); int i, j, k; int nCurCenter = 2; //目前的群心個數 //=2 是因為開始時, 以幾何 //中心為主, 就已分了一群 do { double cost = 0; for (j=0; j<nCurCenter-1; ++j) //將目前的最佳群心 set (old_center[j], center[j]); //保存到 old_center 中 for (i=0; i<nMaxSelect; ++i) //嘗試 nMaxSelect 次 { cur_iter++; //計算 kMean 總共做了幾次 do { //持續嘗試做 k-Means j = nCurCenter - 1; //直到生成一個有效的群心 k = cvRandInt (&rng) % nData; set (center[j], data[k]); } while (!kmeans (nCurCenter)); //分成 nCurCenter 群 if (i == 0) cost = avg_dist; else //是否更新最佳群心集合? if (avg_dist < cost) { cost = avg_dist; for (j=0; j<nCurCenter; ++j) set (best_center[j], center[j]); } } for (j=0; j<nCurCenter; ++j) //將最佳群心拷貝到 center set (center[j], best_center[j]);//供下一輪 k-Mean 使用 // 終止條件判斷 if (old_dist - avg_dist > 2) old_dist = avg_dist; else break; } while (++nCurCenter <= nMaxCenter_); nCenter = nCurCenter - 1; for (j=0; j<nCenter; ++j) //將最佳群心拷貝到 set (center[j], old_center[j]); //center 中 cluster_assignment (); //將資料歸屬到最終的群心 for (i=0; i<2; ++i) { //合併相近的群 joint_clusters (); //共 refine 2 次 cal_avg_dist (); cluster_reassignment (); } };
底下測試框架生成 2D GMM 群落，對各算法進行比對：
template <class T> struct KMeanTest { enum {CD = 3, MAX_CLUSTERS = (CD+1)*(CD+1)*(CD+1), //最大群心數 W=500, H=500}; //展示影像的寬高 CvScalar color_tbl[MAX_CLUSTERS]; IplImage* img [ALGO_NUM]; //2D展示影像 T** data; //資料集合 int nCenter; //群心個數 int nData; //資料項數 int dim; //資料維度 void generate_color_table() { int r, g, b, i=0, d = 255/CD; for (r=0; r<256; r+=d) for (g=0; g<256; g+=d) for (b=0; b<256; b+=d) color_tbl[i++] = CV_RGB (r,g,b); shuffle (color_tbl+1, MAX_CLUSTERS-2); } //以 nCluster 個高斯 (正態) 分佈產生樣本點 //若 nCluster = -1 , 就隨機產生樣本點 //若 nCluster = 0 , 就隨機令 nCluster = 1 ~ MAX_CLUSTERS void generate_sample (int nCluster = -1) { release(); double sx, sy; int cx, cy; int i = 0; data = (T**) new_arr (sizeof(T), 2, nData, dim); if (nCluster == -1) //若中心點數 = -1 for (int i=0; i<nData; ++i) { //就隨機產生樣本點 data[i][0] = T (rand() % W); data[i][1] = T (rand() % H); } else { if (nCluster == 0) nCluster = rand() % MAX_CLUSTERS + 1; CvRNG rng = cvRNG (cvGetTickCount()); CvMat* points = cvCreateMat (nData, 1, CV_32SC2); for (int k = 0; k < nCluster; ++k) //generate random sample { //from multi-gaussian //distribution const int B = 80; //中心點離邊框的距離 cx = B + rand() % (W-2*B); cy = B + rand() % (H-2*B); sx = W * ((0.01 + cvRandReal (&rng))/ 10.0); sy = H * ((0.01 + cvRandReal (&rng))/ 10.0); //令群數 = A, A = (nData/ nCluster) * (1 +- 0.1) double ratio = 1 + (cvRandReal (&rng) - 0.5)/5; points->rows = int ((nData / nCluster) * ratio); //到最後一群時, 若還剩下很多名額, 就都分配給它 if (k == nCluster-1 && points->rows < nData-i) points->rows = nData - i; cvRandArr (&rng, points, CV_RAND_NORMAL, cvScalar (cx, cy), //平均數 cvScalar (sx, sy)); //標準差 for (int j=0; j<points->rows; ++j) { data[i][0] = (T) points->data.i[j*2]; data[i][1] = (T) points->data.i[j*2+ 1]; i++; if (i >= nData) break; } if (i >= nData) break; } cvReleaseMat( &points ); } } };
比對圖例：
3 - GM:
6 - GM:
