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;
}