码迷，mamicode.com
首页 > 编程语言 > 详细

FCM算法

时间：2015-08-17 00:54:58      阅读：229      评论：0      收藏：0      [点我收藏+]

标签：

 

　　之前有看过关于FCM算法的资料，接下来的这几天我就自己学到的一些东西与大家分享~

一 FCM算法概述

     FCM算法的全称是模糊C均值聚类算法，和K-means算法同属于聚类算法，但却有着本质的区别，就其命名而言，模糊二字无疑是该算法的重点，下面就先简单介绍一下：

1. 隶属度和模糊集

　　隶属度函数用来描述元素x属于一个集合B的程度，假定为UB(x)，其中x为B中的任意元素，UB(x)的取值范围为[0,1]。在隶属度函数的基础上，称空间上X={x}上的隶属度函数为一个模糊集合。

2. 模糊聚类分析

　　传统的聚类分析是把每个元素严格的划分到一个类中，属于硬划分。模糊聚类分析将聚类生成的每个簇均看做模糊集合，通过隶属度来确定聚类关系，是一种柔性划分，得到元素属于各个簇的不确定性程度，使得聚类结果更加准确灵活，因此，模糊聚类分析逐渐成为聚类分析的主流。

3.FCM

      FCM算法是将N个L维向量分为C个模糊组，通过迭代不断更新隶属度以及聚类中心，最小化目标函数对数据进行聚类。

其中，目标函数为：

              技术分享

约束条件为：

                   技术分享

根据隶属度的非负性，有技术分享。这里，m指的是模糊加权系数，它的值大于1；d(xi,vc)表示的是第i个数据点与第c个聚类中心的欧式距离；uic是隶属度矩阵中的元素，且[0,1]；vc是对应于每个聚类的聚类中心。

为了求含有约束条件的目标函数的极值，引入拉格朗日因子构造新的目标函数：

                     技术分享

对于目标函数求极值的最优化条件如下：

　　　　　　　技术分享

从而得到隶属度和聚类中心的计算公式为：

  　　　　　　　　技术分享                                 

根据上述公式不断迭代求出满足条件的隶属度以及聚类中心。具体的FCM算法步骤如下：

首先，给定一个由N个L维向量组成的数据集X以及所要分得的类别个数C，自定义隶属度矩阵技术分享

（1）设定类别的个数C和模糊系数m；

（2）?初始化隶属度矩阵且满足公式（2）中的归一化条件；

（3）根据公式（5）计算聚类中心；

（4）根据公式（4）更新隶属度矩阵；

（5）根据矩阵范数比较迭代的隶属度矩阵，如果技术分享，迭代停止，否则返回（3）。

 

二 MATLAB实现

根据上述公式及步骤，采用MATLAB实现，具体如下：

（1）选取UCI数据库中的wine数据实现，代码：

 

clear all
clc;
%% 导入数据
load wine.txt;
cluster_n=6;

data = wine;
data(:,1) = [];
data(:,size(data,2)) = [];
data_n = size(data, 1); %数据的个数 
in_n = size(data, 2);% 数据维数

%% 定义变量
default_options = [10;    % 隶属度函数的幂次方
        300;    % 最大迭代次数
        1e-5;    %步长
        1];    % 判定条件

options = default_options;
expo = options(1);        % Exponent for U 隶属度函数的幂次方
max_iter = options(2);        % Max. iteration 最大迭代次数
min_impro = options(3);        % Min. improvement  最小进化步长
display = options(4);        % Display info or not 显示信息与否
obj_fcn = zeros(max_iter, 1);    % Array for objective function

tic
%% 初始化隶属度矩阵并归一
U = rand(cluster_n, data_n); %rand()产生随机矩阵 
col_sum = sum(U);
U = U./col_sum(ones(cluster_n, 1), :);%归一化

%% 开始迭代
for i = 1:max_iter,%迭代次数控制
    tic
    mf = U.^expo;      % MF matrix after exponential modification
    center = mf*data./((ones(size(data, 2), 1)*sum((mf‘)))‘); % 建立新的聚类中心
   
       out = zeros(size(center, 1), size(data, 1));  %每个点到每个中心的距离，行数为中心数
     
        if size(center, 2) > 1,%样本的维数大于一执行以下程序
            for k = 1:size(center, 1),%给K赋值
                abc = ((data-ones(size(data,1),1) * center(k,:)).^2)‘;
                  out(k, :) = sqrt(sum(abc));%得到欧氏距离
            end
       else    % 1-D data
            for k = 1:size(center, 1),
                  out(k, :) = abs(center(k)-data)‘;
            end
        end
   
    obj_fcn(i) = sum(sum((out.^2).*mf));  % 目标函数
    tmp = out.^(-2/(expo-1));      % 根据新的隶属度矩阵求解公式求出
    U= tmp./(ones(cluster_n, 1)*sum(tmp));  % 新的隶属度矩阵

    if display, 
        fprintf(‘Iteration count = %d, obj. fcn = %f\n‘, i, obj_fcn(i));
        %输出迭代次数和函数的结果
    end
    % check termination condition
    if i > 1,  %进化步长控制
        if abs(obj_fcn(i) - obj_fcn(i-1)) < min_impro, break; end,
    end
    toc
end
toc

plot(U(1,:),‘-ro‘);
grid on
hold on
plot(U(2,:),‘-g*‘);
plot(U(3,:),‘-b+‘);

ylabel(‘Membership degrees‘)
legend(‘FCM1‘,‘FCM2‘,‘FCM3‘,‘location‘,‘northeast‘);

toc

 

结果如下：

  　　　　　　技术分享

 

（2）选取一张MRI图像，进行图像分割，代码如下：

 

clear all 
clc;

a=imread(‘MRI.jpg‘);
I=imnoise(a,‘salt & pepper‘,0.05);
figure(1);           
 imshow(I);title(‘加噪图像‘);

[height,width,c]=size(a);
if c~=1
    a=rgb2gray(a);
end

a=double(a);
[row,column]=size(a); 
data = a(:);
data_n = size(data,1);
       
cluster_num = 4;
default_options = [2.0;    % 隶属度函数的幂次方
        300;    % 最大迭代次数
        1e-5;    %步长
        1];    % 判定条件

options = default_options;


expo = options(1);        % Exponent for U 隶属度函数的幂次方
max_iter = options(2);        % Max. iteration 最大迭代次数
min_impro = options(3);        % Min. improvement  最小进化步长
display = options(4);        % Display info or not 显示信息与否

obj_fcn = zeros(max_iter, 1);    % Array for objective function
membership = zeros(height,width,cluster_num);
center = zeros(cluster_num,1);

tic
% 初始化隶属度并归一
for i=1:height 
    for j=1:width
        member_sum=0;
        for k=1:cluster_num
            membership(i,j,k)=rand();
            member_sum = member_sum + membership(i,j,k);
        end        
        for p =1:cluster_num
            membership(i,j,p) = membership(i,j,p) / member_sum;
        end
    end
end

tic
for i = 1:max_iter,%迭代次数控制
    mf = membership.^expo;     
    %%%%%%%%建立聚类中心
    for m = 1:cluster_num
        to = 0;
        tp =0;
        for j = 1:height
            for t = 1:width
                to = to + membership(j,t,m) * a(j,t);
                tp = tp + membership(j,t,m);
            end
        end
        center(m,1) = to / tp; 
    end
    
    %%%%%%%得到欧式距离以及目标函数
    out = zeros(height,width,cluster_num);
    for m =1:height
        for j =1:width
            for t = 1:cluster_num
                out(m,j,t) = abs(a(m,j) - center(t,1));
                obj_fcn(i) = obj_fcn(i) + (membership(m,j,t).^expo) * (out(m,j,t).^2);
            end
        end
    end
    
    for m = 1:height
        for j = 1:width
            for r = 1:cluster_num
                top =0;
                for t = 1:cluster_num
                    top = top + (out(m,j,r) / out(m,j,t)).^(expo - 1);
                end
                membership(m,j,r) = 1 / top;
            end
        end
    end
    
    %%%%%%归一化隶属度
    for m=1:height 
        for j = 1:width
            member_sum = 0;
            for k = 1:cluster_num
                member_sum = member_sum + membership(m,j,k);
            end
            for p = 1:cluster_num
                membership(m,j,p) = membership(m,j,p) / member_sum;
            end
        end
    end
          
     if display, 
        fprintf(‘Iteration count = %d, obj. fcn = %f\n‘, i, obj_fcn(i));
        %输出迭代次数和函数的结果
    end
    % check termination condition
    if i > 1,  %进化步长控制
        if abs(obj_fcn(i) - obj_fcn(i-1)) < min_impro, break; end,
    end
end
toc

%%%%%%证得如自定义图像中的MCR不能计算，故在此继续尝试直接用newing和A相比较

A = ones(height,width,1);
for i = 1:height
    for j = 1:width
        if (fix(a(i,j) / 85) == 1)
            A(i,j,1) = 2;
        end
        if (fix(a(i,j) / 85) == 2)
            A(i,j,1) = 3;
        end
        if (fix(a(i,j,1) / 85) == 3)
            A(i,j,1) = 4;
        end
    end
end    
A = reshape(A,1,data_n);

  newing = zeros(row,column);
for i=1:row
    for j=1:column         
        maxmembership=membership(i,j,1);
        index=1;
        for k=2:cluster_num            
            if(membership(i,j,k)>maxmembership)
                maxmembership=membership(i,j,k);
                index=k;
            end
        end
        newing(i,j) = round(255 * (1-(index-1)/(cluster_num-1)));
    end 
end

B = reshape(newing,1,data_n);
b = fix((max(B) - B(1,1)) / cluster_num);
for i = 2:data_n
    if B(1,i) == B(1,1)
        B(1,i) = 1;
    elseif (fix(B(1,i) / b) == 2)
        B(1,i) = 2;
    elseif (fix(B(1,i) / b)  == 3)
        B(1,i) = 3;
    else
        B(1,i) = 4;
    end
end
B(1,1) = 1;

sum = 0;
for i = 1:data_n
    if ( A(1,i) ~= B(1,i))
        sum = sum + 1;
    end
end
MCR = sum / data_n;
fprintf(‘MCR = %d\n‘,MCR);

S = 0;
for i = 1:data_n
    S = S + (A(1,i) - B(1,i)).^2;
end

RMS = sqrt(S / (data_n * (data_n -1)));
fprintf(‘RMS = %d\n‘,RMS);

figure(2);
imshow((uint8(newing)));
title(‘FCM分割后的图像‘);

 

结果如下：

           技术分享                技术分享

（3）还可以自定义图像进行测试，代码如下：

clear all
clc;

image = zeros(96,96); 

for i = 1:58
    for j = 49:96
        image(i,j) = 75;
    end
end

for i = 59:77
    for j = 1:96
        image(i,j) = 50;
    end
end

for i = 78:96
    for j = 1:48
        image(i,j) = 50;
    end
end

for i = 78:96
    for j = 49:96
        image(i,j) = 25;
    end
end

a = image;

noise = uint8(image);
noise = imnoise(noise,‘gaussian‘,100 / (255 * 255));
figure(1);
imagesc(noise);
colormap(gray);
title(‘加噪图像‘);

[height,width,c]=size(a);
if c~=1
    a=rgb2gray(a);
end

a=double(a);
[row,column]=size(a); 
data = a(:);
data_n = size(data,1);

cluster_num = 4;
default_options = [10;    % 隶属度函数的幂次方
        300;    % 最大迭代次数
        1e-5;    %步长
        1];    % 判定条件

options = default_options;


expo = options(1);        % Exponent for U 隶属度函数的幂次方
max_iter = options(2);        % Max. iteration 最大迭代次数
min_impro = options(3);        % Min. improvement  最小进化步长
display = options(4);        % Display info or not 显示信息与否

obj_fcn = zeros(max_iter, 1);    % Array for objective function
membership = zeros(height,width,cluster_num);
center = zeros(cluster_num,1);

tic
% 初始化隶属度并归一
for i=1:height 
    for j=1:width
        member_sum=0;
        for k=1:cluster_num
            membership(i,j,k)=rand();
            member_sum = member_sum + membership(i,j,k);
        end        
        for p =1:cluster_num
            membership(i,j,p) = membership(i,j,p) / member_sum;
        end
    end
end

tic
for i = 1:max_iter,%迭代次数控制
    mf = membership.^expo;     
    %%%%%%%%建立聚类中心
    for m = 1:cluster_num
        to = 0;
        tp =0;
        for j = 1:height
            for t = 1:width
                to = to + membership(j,t,m) * a(j,t);
                tp = tp + membership(j,t,m);
            end
        end
        center(m,1) = to / tp; 
    end
    
    %%%%%%%得到欧式距离以及目标函数
    out = zeros(height,width,cluster_num);
    for m =1:height
        for j =1:width
            for t = 1:cluster_num
                out(m,j,t) = abs(a(m,j) - center(t,1));
                obj_fcn(i) = obj_fcn(i) + (membership(m,j,t).^expo) * (out(m,j,t).^2);
            end
        end
    end
    
    for m = 1:height
        for j = 1:width
            for r = 1:cluster_num
                top =0;
                for t = 1:cluster_num
                    top = top + (out(m,j,r) / out(m,j,t)).^(expo - 1);
                end
                membership(m,j,r) = 1 / top;
            end
        end
    end
    
    %%%%%%归一化隶属度
    for m=1:height 
        for j = 1:width
            member_sum = 0;
            for k = 1:cluster_num
                member_sum = member_sum + membership(m,j,k);
            end
            for p = 1:cluster_num
                membership(m,j,p) = membership(m,j,p) / member_sum;
            end
        end
    end
          
     if display, 
        fprintf(‘Iteration count = %d, obj. fcn = %f\n‘, i, obj_fcn(i));
        %输出迭代次数和函数的结果
    end
    % check termination condition
    if i > 1,  %进化步长控制
        if abs(obj_fcn(i) - obj_fcn(i-1)) < min_impro, break; end,
    end
end
toc
 A = ones(height,width,1);
for i = 1:height
    for j = 1:width
        if (a(i,j) == 75)
            A(i,j,1) = 2;
        end
        if (a(i,j) == 50)
            A(i,j,1) = 3;
        end
        if (a(i,j,1) == 25)
            A(i,j,1) = 4;
        end
    end
end    
A = reshape(A,1,data_n);

  B = ones(row,column,1);
for i=1:row
    for j=1:column       
        maxmembership=membership(i,j,1);
        for k=2:cluster_num            
            if(membership(i,j,k)>maxmembership)
                maxmembership = membership(i,j,k);
                B(i,j,1)=k;
            end
        end
    end
end

B = reshape(B,1,data_n);

t = MCR(A,B);

%输出图像
  newing = zeros(row,column);
for i=1:row
    for j=1:column         
        maxmembership=membership(i,j,1);
        index=1;
        for k=2:cluster_num            
            if(membership(i,j,k)>maxmembership)
                maxmembership=membership(i,j,k);
                index=k;
            end
        end
        newing(i,j) = round(255 * (1-(index-1)/(cluster_num-1)));
    end 
end
figure(2);
imagesc(uint8(newing));
colormap(gray);
title(‘FCM分割后的图像‘);

结果如下：

        技术分享                   技术分享

 

结果分析：

综上所示，计算得出误分类率为38.9%，所以说此算法还是准确率还有待提高。且当数据的特征复杂集合较大时，算法的实时性不是很可观；与许多聚类算法一样，FCM选择欧氏距离作为样本点与相应聚类中心之间的非相似性指标，致使算法趋向于发现具有相近尺度和密度的球星簇。

 

 

 今天就先到这儿，之后将进一步对算法进行优化~

 

 

FCM算法

标签：

原文地址：http://www.cnblogs.com/ybjourney/p/4735335.html
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！