【OpenCV學習筆記 020】K-Means聚類演算法介紹及實現
一、K-Means演算法介紹
在資料探勘中，K-Means演算法是一種cluster analysis的演算法，其主要是來計算資料聚集的演算法，主要通過不斷地取離種子點最近均值的演算法。
問題
K-Means演算法主要解決的問題如下圖所示。我們可以看到，在圖的左邊有一些點，我們用肉眼可以看出來有四個點群，但是我們怎麼通過計算機程式找出這幾個點群來呢？於是就出現了我們的K-Means演算法（Wikipedia連結）
K-Means要解決的問題
演算法概要：
這個演算法其實很簡單，如下圖所示：
從上圖中，我們可以看到，A，B，C，D，E是五個在圖中點。而灰色的點是我們的種子點，也就是我們用來找點群的點
然後，K-Means的演算法如下：
- 隨機在圖中取K（這裡K=2）個種子點。
- 然後對圖中的所有點求到這K個種子點的距離，假如點Pi離種子點Si最近，那麼Pi屬於Si點群。（上圖中，我們可以看到A，B屬於上面的種子點，C，D，E屬於下面中部的種子點）
- 接下來，我們要移動種子點到屬於他的“點群”的中心。（見圖上的第三步）
- 然後重複第2）和第3）步，直到，種子點沒有移動（我們可以看到圖中的第四步上面的種子點聚合了A，B，C，下面的種子點聚合了D，E）。
這個演算法很簡單，但是有些細節我要提一下，求距離的公式我不說了，大家有初中畢業水平的人都應該知道怎麼算的。我重點想說一下“求點群中心的演算法”。
求點群中心的演算法
一般來說，求點群中心點的演算法你可以很簡單的使用各個點的X/Y座標的平均值。不過，我這裡想告訴大家另三個求中心點的的公式：
1）Minkowski Distance公式——λ可以隨意取值，可以是負數，也可以是正數，或是無窮大。
2）Euclidean Distance公式——也就是第一個公式λ=2的情況
3）CityBlock Distance公式——也就是第一個公式λ=1的情況
這三個公式的求中心點有一些不一樣的地方，我們看下圖（對於第一個λ在0-1之間）。
（1）Minkowski Distance （2）Euclidean Distance （3） CityBlock Distance
K-Means的演示
操作是，滑鼠左鍵是初始化點，右鍵初始化“種子點”，然後勾選“Show History”可以看到一步一步的迭代。
K-Means++演算法
K-Means主要有兩個最重大的缺陷——都和初始值有關：
- K是事先給定的，這個K值的選定是非常難以估計的。很多時候，事先並不知道給定的資料集應該分成多少個類別才最合適。（ISODATA演算法通過類的自動合併和分裂，得到較為合理的型別數目K）
- K-Means演算法需要用初始隨機種子點來搞，這個隨機種子點太重要，不同的隨機種子點會有得到完全不同的結果。（K-Means++演算法可以用來解決這個問題，其可以有效地選擇初始點）
我在這裡重點說一下K-Means++演算法步驟：
- 先從我們的資料庫隨機挑個隨機點當“種子點”。
- 對於每個點，我們都計算其和最近的一個“種子點”的距離D(x)並儲存在一個數組裡，然後把這些距離加起來得到Sum(D(x))。
- 然後，再取一個隨機值，用權重的方式來取計算下一個“種子點”。這個演算法的實現是，先取一個能落在Sum(D(x))中的隨機值Random，然後用Random -= D(x)，直到其<=0，此時的點就是下一個“種子點”。
- 重複第（2）和第（3）步直到所有的K個種子點都被選出來。
- 進行K-Means演算法。
K-Means演算法應用
看到這裡，你會說，K-Means演算法看來很簡單，而且好像就是在玩座標點，沒什麼真實用處。而且，這個演算法缺陷很多，還不如人工呢。是的，前面的例子只是玩二維座標點，的確沒什麼意思。但是你想一下下面的幾個問題：
1）如果不是二維的，是多維的，如5維的，那麼，就只能用計算機來計算了。
2）二維座標點的X，Y 座標，其實是一種向量，是一種數學抽象。現實世界中很多屬性是可以抽象成向量的，比如，我們的年齡，我們的喜好，我們的商品，等等，能抽象成向量的目的就是可以讓計算機知道某兩個屬性間的距離。如：我們認為，18歲的人離24歲的人的距離要比離12歲的距離要近，鞋子這個商品離衣服這個商品的距離要比電腦要近，等等。
只要能把現實世界的物體的屬性抽象成向量，就可以用K-Means演算法來歸類了。
在《k均值聚類(K-means)》 這篇文章中舉了一個很不錯的應用例子，作者用亞洲15支足球隊的2005年到1010年的戰績做了一個向量表，然後用K-Means把球隊歸類，得出了下面的結果，呵呵。
- 亞洲一流：日本，韓國，伊朗，沙特
- 亞洲二流：烏茲別克，巴林，朝鮮
- 亞洲三流：中國，伊拉克，卡達，阿聯酋，泰國，越南，阿曼，印尼
其實，這樣的業務例子還有很多，比如，分析一個公司的客戶分類，這樣可以對不同的客戶使用不同的商業策略，或是電子商務中分析商品相似度，歸類商品，從而可以使用一些不同的銷售策略，等等。
二、OpenCV自帶samples中kmeans研究
實驗基礎
在使用kmeans之前，必須先了解kmeans演算法的2個缺點：第一是必須人為指定所聚的類的個數k；第二是如果使用歐式距離來衡量相似度的話，可能會得到錯誤的結果，因為沒有考慮到屬性的重要性和相關性。為了減少這種錯誤，在使用kmeans距離時，一定要使樣本的每一維資料歸一化，不然的話由於樣本的屬性範圍不同會導致錯誤的結果。
實驗一是對隨機產生的sampleCount個二維樣本（共分為clusterCount個類別），每個類別的樣本資料都服從高斯分佈，該高斯分佈的均值是隨機的，方差是固定的。然後對這sampleCount個樣本資料使用kmeans演算法聚類，最後將不同的類用不同的顏色顯示出來。
下面是程式中使用到的幾個OpenCV函式：
void RNG::fill(InputOutputArray mat, int distType, InputArray a,InputArray b, bool saturateRange=false )void randShuffle(InputOutputArray dst, double iterFactor=1., RNG*rng=0 )Class TermCriteria類TermCriteria 一般表示迭代終止的條件，如果為CV_TERMCRIT_ITER，則用最大迭代次數作為終止條件，如果為CV_TERMCRIT_EPS 則用精度作為迭代條件，如果為CV_TERMCRIT_ITER+CV_TERMCRIT_EPS則用最大迭代次數或者精度作為迭代條件，看哪個條件先滿足。
double kmeans(InputArray data, int K, InputOutputArray bestLabels,TermCriteria criteria, int attempts, int flags, OutputArray centers=noArray() )該函式為kmeans聚類演算法實現函式。引數data表示需要被聚類的原始資料集合，一行表示一個數據樣本，每一個樣本的每一列都是一個屬性；引數k表示需要被聚類的個數；引數bestLabels表示每一個樣本的類的標籤，是一個整數，從0開始的索引整數；引數criteria表示的是演算法迭代終止條件；引數attempts表示執行kmeans的次數，取結果最好的那次聚類為最終的聚類，要配合下一個引數flages來使用；引數flags表示的是聚類初始化的條件。其取值有3種情況，如果為KMEANS_RANDOM_CENTERS，則表示為隨機選取初始化中心點，如果為KMEANS_PP_CENTERS則表示使用某一種演算法來確定初始聚類的點；如果為KMEANS_USE_INITIAL_LABELS，則表示使用使用者自定義的初始點，但是如果此時的attempts大於1，則後面的聚類初始點依舊使用隨機的方式；引數centers表示的是聚類後的中心點存放矩陣。該函式返回的是聚類結果的緊湊性，其計算公式為：
實驗一：隨機產生的符合高斯分佈的資料被聚類的結果
實驗程式碼及註釋
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/core/core.hpp"
#include <iostream>
using namespace cv;
using namespace std;
// static void help()
// {
// cout << "\nThis program demonstrates kmeans clustering.\n"
// "It generates an image with random points, then assigns a random number of cluster\n"
// "centers and uses kmeans to move those cluster centers to their representitive location\n"
// "Call\n"
// "./kmeans\n" << endl;
// }
int main(int /*argc*/, char** /*argv*/)
{
const int MAX_CLUSTERS = 5;
Scalar colorTab[] = //因為最多隻有5類，所以最多也就給5個顏色
{
Scalar(0, 0, 255),
Scalar(0, 255, 0),
Scalar(255, 100, 100),
Scalar(255, 0, 255),
Scalar(0, 255, 255)
};
Mat img(500, 500, CV_8UC3);
RNG rng(12345); //隨機數產生器
for (;;)
{
int k, clusterCount = rng.uniform(2, MAX_CLUSTERS + 1);
int i, sampleCount = rng.uniform(1, 1001);
Mat points(sampleCount, 1, CV_32FC2), labels; //產生的樣本數，實際上為2通道的列向量，元素型別為Point2f
clusterCount = MIN(clusterCount, sampleCount);
Mat centers(clusterCount, 1, points.type()); //用來儲存聚類後的中心點
/* generate random sample from multigaussian distribution */
for (k = 0; k < clusterCount; k++) //產生隨機數
{
Point center;
center.x = rng.uniform(0, img.cols);
center.y = rng.uniform(0, img.rows);
Mat pointChunk = points.rowRange(k*sampleCount / clusterCount,
k == clusterCount - 1 ? sampleCount :
(k + 1)*sampleCount / clusterCount); //最後一個類的樣本數不一定是平分的，
//剩下的一份都給最後一類
//每一類都是同樣的方差，只是均值不同而已
rng.fill(pointChunk, CV_RAND_NORMAL, Scalar(center.x, center.y), Scalar(img.cols*0.05, img.rows*0.05));
}
randShuffle(points, 1, &rng); //因為要聚類，所以先隨機打亂points裡面的點，注意points和pointChunk是共用資料的。
kmeans(points, clusterCount, labels,
TermCriteria(CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 10, 1.0),
3, KMEANS_PP_CENTERS, centers); //聚類3次，取結果最好的那次，聚類的初始化採用PP特定的隨機演算法。
img = Scalar::all(0);
for (i = 0; i < sampleCount; i++)
{
int clusterIdx = labels.at<int>(i);
Point ipt = points.at<Point2f>(i);
circle(img, ipt, 2, colorTab[clusterIdx], CV_FILLED, CV_AA);
}
imshow("clusters", img);
char key = (char)waitKey(); //無限等待
if (key == 27 || key == 'q' || key == 'Q') // 'ESC'
break;
}
return 0;
}
實驗二是對輸入的一張圖片進行聚類。
實驗程式碼及註釋
#include<core/core.hpp>
#include<highgui/highgui.hpp>
#include<imgproc/imgproc.hpp>
#include<ml/ml.hpp>
using namespace std;
using namespace cv;
int main(int argc, char **argv)
{
/*if (argc<3)
return -1;
Mat Img = imread(argv[1], IMREAD_COLOR);
const int K = atoi(argv[2]);*/
Mat Img = imread("beach.jpg", IMREAD_COLOR);//源影象
const int K = 8;//聚類個數
if (Img.empty() || K <= 0)
return -1;
Mat ImgHSV;
cvtColor(Img, ImgHSV, CV_BGR2HSV);//將RGB空間轉換為HSV空間
Mat ImgData(Img.rows*Img.cols, 1, CV_32FC3);//一個3通道的資料域
Mat_<Vec3b>::iterator itImg = ImgHSV.begin<Vec3b>();//迭代器
Mat_<Vec3f>::iterator itData = ImgData.begin<Vec3f>();
//將源影象的資料輸入給新建的ImgData
for (; itImg != ImgHSV.end<Vec3b>(); ++itImg, ++itData)
*itData = *itImg;
Mat ImgLabel, ImgCenter;
kmeans(ImgData, K, ImgLabel, TermCriteria(CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 10, 1.0), 1, KMEANS_PP_CENTERS, ImgCenter);
Mat ImgRes(Img.size(), CV_8UC3);
Mat_<Vec3b>::iterator itRes = ImgRes.begin<Vec3b>();
Mat_<int>::iterator itLabel = ImgLabel.begin<int>();
for (; itLabel != ImgLabel.end<int>(); ++itLabel, ++itRes)
*itRes = ImgCenter.at<Vec3f>(*itLabel, 0);
//cout << ImgCenter << endl;
cvtColor(ImgRes, ImgRes, CV_HSV2BGR);
imwrite("out.jpg", ImgRes);
namedWindow("Img", WINDOW_AUTOSIZE);
imshow("Img", Img);
namedWindow("Res", WINDOW_AUTOSIZE);
imshow("Res", ImgRes);
waitKey();
return 0;
}
聚類後的影象：
三、OpenCV的Kmean演算法實踐
首先看一下opencv中kmeans的函式原型
//! clusters the input data using k-Means algorithm
CV_EXPORTS_W double kmeans( InputArray data, int K, CV_OUT InputOutputArray bestLabels,
TermCriteria criteria, int attempts,
int flags, OutputArray centers=noArray() );K： 型別數目；
bestLabels： 分類後的矩陣，每個樣本對應一個型別label；
TermCriteria criteria：結束條件（最大迭代數和理想精度）
int attempts：根據最後一個引數確定選取的最理想初始聚類中心（選取attempt次初始中心，選擇compactness最小的）；
int flags ：
Flag that can take the following values:
KMEANS_RANDOM_CENTERS Select random initial centers in each attempt.
KMEANS_PP_CENTERS Use kmeans++ center initialization by Arthur and Vassilvitskii [Arthur2007].
KMEANS_USE_INITIAL_LABELS During the first (and possibly the only) attempt, use the user-supplied labels instead of computing them from the initial centers. For the second and further attempts, use the random or semi-random centers. Use one of KMEANS_*_CENTERS
flag to specify the exact method.
centers：輸出聚類中心，每行一箇中心
compactness： 測試初始中心是否最優
程式原始碼：
#include <string>
#include <iostream>
#include <math.h>
#include <vector>
#include <map>
#include "opencv/cv.h"
#include "opencv/highgui.h"
#include "opencv/cxcore.h"
#define ClusterNum 6
using namespace cv;
using namespace std;
string fileName = "image.jpg";
Mat clustering(Mat src){
int row = src.rows;
int col = src.cols;
unsigned long int size = row*col;
Mat clusters(size, 1, CV_32SC1); //clustering Mat, save class label at every location;
//convert src Mat to sample srcPoint.
Mat srcPoint(size, 1, CV_32FC3);
Vec3f* srcPoint_p = (Vec3f*)srcPoint.data;
Vec3f* src_p = (Vec3f*)src.data;
unsigned long int i;
for (i = 0; i < size;i++)
{
*srcPoint_p = *src_p;
srcPoint_p++;
src_p++;
}
Mat center(ClusterNum, 1, CV_32FC3);
double compactness;//compactness to measure the clustering center dist sum by different flag
compactness = kmeans(srcPoint, ClusterNum, clusters,
cvTermCriteria(CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 10, 0.1), ClusterNum,
KMEANS_PP_CENTERS, center);
cout << "center row:" << center.rows << " col:" << center.cols << endl;
for (int y = 0; y < center.rows; y++)
{
Vec3f* imgData = center.ptr<Vec3f>(y);
for (int x = 0; x < center.cols; x++)
{
cout << imgData[x].val[0] << " " << imgData[x].val[1] << " " << imgData[x].val[2] << endl;
}
cout << endl;
}
double minH, maxH;
minMaxLoc(clusters, &minH, &maxH); //remember must use
cout << "H-channel min:" << minH << " max:" << maxH << endl;
int* clusters_p = (int*)clusters.data;
//show label mat
Mat label(src.size(), CV_32SC1);
int* label_p = (int*)label.data;
//assign the clusters to Mat label
for (i = 0; i < size; i++)
{
*label_p = *clusters_p;
label_p++;
clusters_p++;
}
Mat label_show;
label.convertTo(label_show, CV_8UC1);
normalize(label_show, label_show, 255, 0, CV_MINMAX);
imshow("label", label_show);
map<int, int> count; //map<id,num>
map<int, Vec3f> avg; //map<id,color>
//compute average color value of one label
for (int y = 0; y < row; y++)
{
const Vec3f* imgData = src.ptr<Vec3f>(y);
int* idx = label.ptr<int>(y);
for (int x = 0; x < col; x++)
{
avg[idx[x]] += imgData[x];
count[idx[x]] ++;
}
}
//output the average value (clustering center)
//計算所得的聚類中心與kmean函式中center的第一列一致，
//以後可以省去後面這些繁複的計算，直接利用center,
//但是仍然不理解center的除第一列以外的其他列所代表的意思
for (i = 0; i < ClusterNum; i++)
{
avg[i] /= count[i];
if (avg[i].val[0]>0 && avg[i].val[1]>0 && avg[i].val[2]>0)
{
cout << i << ": " << avg[i].val[0] << " " << avg[i].val[1] << " " << avg[i].val[2] << " count:" << count[i] << endl;
}
}
//show the clustering img;
Mat showImg(src.size(), CV_32FC3);
for (int y = 0; y < row;y++)
{
Vec3f* imgData = showImg.ptr<Vec3f>(y);
int *idx = label.ptr<int>(y);
for (int x = 0; x < col;x++)
{
int id = idx[x];
imgData[x].val[0] = avg[id].val[0];
imgData[x].val[1] = avg[id].val[1];
imgData[x].val[2] = avg[id].val[2];
}
}
normalize(showImg, showImg, 1, 0, CV_MINMAX);
imshow("show", showImg);
waitKey();
return label;
}
int main(){
Mat img = imread(fileName, 1);
imshow("src", img);
GaussianBlur(img, img, Size(3, 3), 0);
img.convertTo(img, CV_32FC3);
Mat pixId = clustering(img);
}
Refrence
http://www.csdn.net/article/2012-07-03/2807073-k-means
http://www.cnblogs.com/moondark/archive/2012/03/08/2385770.html
http://blog.csdn.net/yangtrees/article/details/7971405
http://blog.csdn.net/ts_zxc/article/details/23668751
http://www.opencv.org.cn/opencvdoc/2.3.2/html/modules/core/doc/clustering.html#kmeans
如果文章對更多的朋友有益，請分享到朋友圈。【視音訊影象技術乾貨，FFmpeg流媒體、OpenCV影象及深度學習技術探索，開源專案推薦，還有更多職場規劃】歡迎關注我的微信公眾號，掃一掃下方二維碼或者長按識別二維碼。
