FCM-σ模糊聚类的方法（σ：加权平均距离）

 本方法是对台湾大学蔡杜明教授所写论文里面图像分割方法的matlab复现，原论文当中未附带matlab程序。原理类似http://blog.csdn.net/on2way/article/details/47087201里面的原理介绍，注意里面有一个小错误注意别理解错了。

 注：这篇原理介绍文章中有一处错误，即u11与u12的求法是错误的，大家不要理解错了，这样的隶属度矩阵每一列的就不是1了，正确的求法是：u11=1/{[||x1-c1||/||x1-c1||]^o + [||x1-c1||/||x2-c1||]^o } 其中o=2/（m-1）,m为自己定的一个大于1的数。

 本方法是对传统模糊聚类FCM方法的改进，通过增加一个加权平均距离σ，来调节样本点到聚类中心的欧式距离，将图片分成C个类。本文只给出具体结果，推导过程望读者根据http://blog.csdn.net/on2way/article/details/47087201自己揣摩。

样本点到聚类中心的距离,我们这样定义：

其中加权平均距离：

我们的目的是要求得目标函数的条件极值，应用拉格朗日乗子法求多元函数的条件极小值。

目标函数的限制条件是：

通过求取偏导得到隶属度函数以及聚类中心：

隶属度矩阵：

其中： ，N代表将图片转化成一维矩阵所对应的样本点总数，C代表类的数目 聚类中心： 。

我们可以看到，传统FCM加入加权平均距离 后，聚类中心没有改变，只有隶属度矩阵以及目标函数受到了影响。

下面我们给出具体的MATLAB程序，来说明FCM-σ模糊聚类方法在图像分割中的应用：

clc
clear
close all
%% create samples:
img = imread('zhaocancan.jpg');
%img为uint8类型的
subplot(1,2,1),imshow(img,[])
data = double(img(:));
%data对应的是double类型的
cluster_n = 5;
%类别数
iter = 20;
%迭代次数
m = 2;
%m为大于1常数

num_data = size(data,1)   %data的行数
num_d = size(data,2)   %data的列数
U = rand(cluster_n,num_data);
col_sum = sum(U);
col_sum(ones(cluster_n,1),:);
U = U./col_sum(ones(cluster_n,1),:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:iter
    u_sum_r=zeros(cluster_n,1);
    for j = 1:cluster_n
        u_ij_m = U(j,:).^m;
        sum_u_ij = sum(u_ij_m);
        sum_1d = u_ij_m./sum_u_ij;
        c(j,:) = u_ij_m*data./sum_u_ij;
        u_sum_r(j)=sum_u_ij;  %将所有的聚类中心存到u_sum_r中      
    end 
    %-计算目标函数J 
    temp1 = zeros(cluster_n,num_data);
    for j = 1:cluster_n
        for k = 1:num_data
            temp1(j,k) =U(j,k)^m*(norm(data(k,:)-c(j,:)))^2;
        end
    end
   % Q=zeros(cluster_n,1);
    Q=(sum(temp1,2)./u_sum_r).^0.5 ;    %Q对应的就是加权平均距离σ
    temp1=temp1./ (Q*ones(1,num_data));
    J(i) = sum(sum(temp1));
    %更新U
    for j = 1:cluster_n
        for k = 1:num_data
            sum1 = 0;
            for j1 = 1:cluster_n
                temp = (norm(data(k,:)-c(j,:))./(Q(j)^0.5)/norm(data(k,:)-c(j1,:))).^(2/(m-1));
                sum1 = sum1 + temp;
            end
            U(j,k) = 1./sum1;
        end
    end
end
[~,label] = max(U)%找到数据所属的类

img_new = reshape(label,size(img));
img_new1=uint8(img_new);%再将double型的转换成uint8型的用于显示。
subplot(1,2,2),imshow(mat2gray(img_new1),[])