Fuzzy C-mean（FCM，模糊C均值）聚類：
N個樣本
，将它們聚類到C個集合中，使得目标函數J(u)最小
其中，uij指的是第i個樣本xi屬于第j個聚類中心點cj的機率值，取值範圍[0.0-1.0]。cj是聚類中心點。
顯然，對于以上的最小問題求解，除了樣本xi已知，其餘兩個都是未知量，不能直接求解析式uij和cj的解，因為這兩個未知數是互相相關的，如cj由樣本xi和uij共同決定，而uij又是由樣本xi和cj共同決定。是以用EM算法來求解（EM算法思想請參考之前部落格的相關資料）：
（1）假設uij已知，固定uij，計算更新聚類類别中心點cj，計算公式如下：
上面公式的意思是：将每個樣本對某一個類别的貢獻值（即機率值，再多一個m次幂）×樣本向量，再進行機率值歸一化（即分母中所有機率值相加）
（2）再假設cj已知，固定cj，計算更新，計算公式如下：
（需要注意的是，上面分母中的計算，幂指數為先，求和）：
關于初始值的設定，uij可以取任意值，然後依次重複以上兩個步驟，直至最後結果穩定（uij 和cj 收斂）。
FCM與K-mean聚類有點類似，又有點不同。
（a）相同點：都是聚類算法，需要指定聚類的類别個數；
（b）不同點：K-mean聚類對每個樣本x(i)是一個“硬指派”，即通過計算它與某一個聚類中心c(j)的距離，根據距離長短來決定該樣本屬于哪一個類别
. 而FCM聚類對樣本的類别指派是一個“軟指派”，即不強求它屬于某一類，而是以機率值的形式表示了它歸屬于每個類别的可能性
，也就是上面公式中的uij，顯然對于每個樣本來說，它歸屬于不同類别的機率值和為1.0，即有：
Fuzzy C Mean缺點：
- 需要指定聚類中心點的個數；
- 由于目标函數J是個非凸函數，最後計算結果不一定能達到全局最小值（可能是一個局部最小值，即使如此實際應用中問題也不大），如果希望能夠修正這個缺點，可以多取幾個參數初始值進行疊代計算；
- 疊代時間依賴于初始化的聚類中心點。
// Fuzzy_C_means.cpp : 定義控制台應用程式的入口點。
//
#include "stdafx.h"
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/ml/ml.hpp>
#include <cv.h>
#include <iostream>
#include "FCM.h"
using namespace cv;
using namespace std;
int main()
{
Mat img;
CFCM fcm;
fcm.loadImage("org.png");
fcm.showOrgImage("org", img);
waitKey(0);
fcm.cluster(10, 2, 0.1);
fcm.showClusterResult("result", img);
waitKey(0);
}
編寫FCM類
class CFCM
{
public:
CFCM();
~CFCM();
void loadImage(string addr, int flags = 1);
void loadImage(Mat& img);
void init(int Nc);
int cluster(int Nc, double m, double eps);
void showOrgImage(string win_name, Mat img);
void showClusterResult(string win_name, Mat img);
private:
void initImgVec();
void initCentroid(int Nc);
void updateFuzzyMat();
void updateCentroid();
public:
Mat m_uImg;
Mat m_dImg;
Mat m_imgVec;
Mat m_fuzzyMat; //Np * Nc
Mat m_centrMat; //Nc * channels
int m_Np;
int m_channels;
int m_Nc;
double m_fuzzyVal;
bool m_init;
};
成員函數：注意，因為openCV中的Mat類強大的計算能力，以及代碼實作的簡約，以下的updateFuzzyMat()和updateCentroid()将計算過程進行了向量化。
CFCM::CFCM()
{
m_init = false;
}
CFCM::~CFCM()
{
m_dImg.release();
}
void CFCM::loadImage(string addr, int flags)
{
m_uImg = imread(addr, flags);
m_channels = m_uImg.channels();
m_dImg = Mat_<Vec3d>(m_uImg);
}
void CFCM::loadImage(Mat& img)
{
m_dImg = Mat_<Vec3d>(img);
}
void CFCM::init(int Nc)
{
cout<<"init...";
m_init = true;
m_Np = m_dImg.rows * m_dImg.cols;
m_channels = m_dImg.channels();
initImgVec();
initCentroid(Nc);
cout<<"done!"<<endl;
}
void CFCM::initImgVec()
{
//将原圖像重新排列，變成 m_Np * m_channels 矩陣
//即每一行是一個像素點的樣本資料
m_imgVec = m_dImg.reshape(1, m_Np);
}
void CFCM::initCentroid(int Nc)
{
m_Nc = Nc;
m_centrMat = Mat(m_Nc, m_channels, CV_64FC1);
Mat m_r = m_centrMat.reshape(m_channels, m_Nc);
int rows = m_uImg.rows;
int cols = m_uImg.cols;
//随機圖像中的幾個點作為聚類中心點
RNG rng(m_Nc);
rng.fill(m_centrMat, RNG::UNIFORM, Scalar(0.0), Scalar(256.0));
//for (int i=0; i<m_Nc; i++)
//{
// int p = rng.uniform(0, m_Np);
// m_imgVec.row(p).copyTo(m_centrMat.row(i));
//}
//
}
int CFCM::cluster(int _Nc, double _fuzzyVal, double eps)
{
m_fuzzyVal = _fuzzyVal;
if(!m_init)
init(_Nc);
cout<<"starting cluster...";
Mat old_fuzzy;
m_fuzzyMat.copyTo(old_fuzzy);
double max_v;
clock_t start = clock();
int count = 0;
do
{
updateFuzzyMat();
updateCentroid();
max_v = max_diff_ratio(old_fuzzy, m_fuzzyMat);
m_fuzzyMat.copyTo(old_fuzzy);
count++;
} while (max_v>eps);
clock_t finish = clock();
cout<<"done!"<<endl;
cout<<"duration: "<<(double)finish-start/CLOCKS_PER_SEC<<endl;
cout<<"iteration: "<<count<<endl;
return count;
}
void CFCM::updateFuzzyMat()
{
//計算矩陣U(i,j)即m_fuzzyMat
Mat p_vec = m_imgVec.reshape(m_channels, m_Np);//N_p * 1 * (channels)
Mat p_mat = repeat(p_vec, 1, m_Nc);//N_p * N_c * (channels)
Mat c_vec = m_centrMat.reshape(m_channels, 1);//1 * N_c * (channels)
Mat c_mat = repeat(c_vec, m_Np, 1);//N_p * N_c * (channels)
Mat p_sub_c = p_mat - c_mat;
//計算p_mat和c_mat之間的各個距離，即p_sub_c的通道平方和的開方值
vector<Mat> planes;//N_p * N_c
Mat plane_sqsum(m_Np, m_Nc, CV_64FC1, Scalar::all(0));
Mat pc_norm, m_pow;
split(p_sub_c, planes);
for (int i=0; i<planes.size(); i++)
{
cv::pow(planes.at(i), 2.0, m_pow);
plane_sqsum += m_pow;
cv::sqrt(plane_sqsum, pc_norm);
}
//計算pc_norm的2/(m_fuzzyVal-1)次矩陣pc_mat
Mat pc_mat;
cv::pow(pc_norm, 2.0/(m_fuzzyVal-1), pc_mat);//N_p * N_c
//計算pc_mat的倒數矩陣,且求每行和的向量
//即求每個像素點到不同聚類中心點的距離倒數之和
Mat r_vec = sum_rows(1.0 / pc_mat);
//更新fuzzy矩陣
Mat r_mat = repeat(r_vec, 1, m_Nc);//N_p * N_c
m_fuzzyMat = 1.0/(pc_mat.mul(r_mat));
}
void CFCM::updateCentroid()
{
//計算聚類點
Mat m_pow;
cv::pow(m_fuzzyMat, m_fuzzyVal, m_pow);//N_p * N_c
Mat c_vec = m_pow.t() * m_imgVec; //N_c * N_p * N_p * channels
Mat c_sum = sum_cols(m_pow).t(); //N_c * 1
Mat c_s = repeat(c_sum, 1, m_channels);//N_c * channels
m_centrMat = c_vec.mul(1.0/c_s);//N_c * channels
}
void CFCM::showOrgImage(string win_name, Mat img)
{
m_uImg.copyTo(img);
imshow(win_name, img);
}
void CFCM::showClusterResult(string win_name, Mat img)
{
img = Mat(m_uImg.size(), m_uImg.type());
Mat img_r = img.reshape(m_channels, m_Np);
double max_v, min_v;
int max_p[2], min_p[2];
Mat centroid = m_centrMat.reshape(m_channels, m_Nc);
for (int i=0; i<m_Np; i++)
{
Mat row_m = m_fuzzyMat.row(i);
minMaxIdx(row_m, &min_v, &max_v, min_p, max_p);
int c_p = max_p[1];
Vec3b v = (Vec3b)centroid.at<Vec3d>(c_p,0);
img_r.at<Vec3b>(i,0) = v;
}
imshow(win_name, img);
}
輔助函數：
Mat sum_cols(Mat m)
{
int cols = m.cols;
Mat new_m (1, cols, m.type());
for (int i = 0; i<cols; i++)
{
Mat m_col = m.col(i);
Scalar s = cv::sum(m_col);
new_m.col(i) = s;
}
return new_m;
}
Mat sum_rows(Mat m)
{
int rows = m.rows;
Mat new_m (rows, 1, m.type());
for (int i=0; i<rows; i++)
{
Mat m_row = m.row(i);
Scalar s = cv::sum(m_row);
new_m.row(i) = s;
}
return new_m;
}
double max_diff_ratio(Mat m_old, Mat m_new)
{
Mat diff = cv::abs(m_new - m_old);
Mat divs = cv::abs(m_old) + Scalar::all(0.0001);
Mat r = diff.mul(1.0/divs);
double max_v, min_v;
minMaxIdx(r, &min_v, &max_v);
return max_v;
}
運作結果：
（Nc=3，運作時間15s，疊代次數22）如果增加聚類中心點個數，則運作時間會大大地變長！
參考資料
英文介紹：http://home.deib.polimi.it/matteucc/Clustering/tutorial_html/cmeans.html
中文了解與執行個體：http://blog.csdn.net/jia20003/article/details/8800197