Fuzzy C-mean 聚类原理及图像颜色分割的实例_skfuzzy中的cmeans分割灰度图像 代码

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