Fuzzy C-mean 聚类原理及图像颜色分割的实例

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

中文理解与实例：/article/1391071.html