图像模糊C均值聚类分割代码

 

转自：直觉模糊C均值聚类与图像阈值分割 - liyuefeilong的专栏 - CSDN博客 https://blog.csdn.net/liyuefeilong/article/details/43816495

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 主函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function main
ima = imread('MR6.jpg');
% 先设定FCM的几个初始参数
options=[2;    % FCM公式中的参数m
         100;    % 最大迭代次数
         1e-5];    % 目标函数的最小误差        
class_number = 4;  % 分为4类
imt = ImageSegmentation(ima,class_number,options)
subplot(1,2,1),imshow(ima),title('原图');
subplot(1,2,2),imshow(imt); %显示生成的分割的图像
kk = strcat('分割成',int2str(class_number),'类的输出图像');
title(kk);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ImageSegmentation()函数：实现聚类分割图像
% 输入：file为灰度图像文件 cluster_n为聚类类别个数 options为预设的初始参数
% 输出分割后的图像  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function imt = ImageSegmentation(file, cluster_n, options) 
ima = file;
I = im2double(file); 
[x,y] = size(ima);
number = x * y;  % 图像的元素个数numel(I)
data = reshape(I,number,1); %将矩阵元素转换为一列数据
[center, U] = FCMprocess(data,cluster_n,options); %调用FCMData函数进行聚类
% 对于每个元素对不同聚类中心的隶属度，找出最大的那个隶属度
maxU = max(U); % 找出每一列的最大隶属度
temp = sort(center); 
for i = 1:cluster_n; % 按聚类结果分割图像
    % 前面求出每个元素的最大隶属度，属于各聚类中心的元素坐标，并存放这些坐标
    % 调用eval函数将括号里的字符串转化为命令执行
    eval(['class_',int2str(i), '= find(U(', int2str(i), ',:) == maxU);']); 
    %gray = round(255 * (i-1) / (cluster_n-1));
        index = find(temp == center(i)); 
    switch index 
        case 1 
            gray = 0; 
        case cluster_n 
            gray = 255; 
        otherwise 
            gray = fix(255*(index-1)/(cluster_n-1)); 
    end 
    eval(['I(class_',int2str(i), '(:))=', int2str(gray),';']); 
end; 
imt = mat2gray(I); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 用于计算聚类中心、隶属度矩阵和目标函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [center, U] = FCMprocess(data, cluster_num, options) 
%data为聚类数据,cluster_num为类别数
m = options(1);        % 参数m
max_iteration = options(2);        % 最终的迭代次数
min_deviation = options(3);        % 最小判别误差
data_number = size(data, 1);    % 元素个数
obj_function = zeros(max_iteration, 1); % obj_function用于存放目标函数的值
% 生成隶属度矩阵U
U = rand(cluster_num, data_number); % 随机生成隶属度矩阵U
sumU = sum(U,1);   % 计算U中每列元素和
for k = 1:data_number
    U(:,k) = U(:,k) ./ sumU(k);  % 对隶属矩阵U进行归一化处理
end

for i = 1:max_iteration
    [U, center, obj_function(i)] = FCMStep(data, U, cluster_num, m); %调用FCMStep函数进行迭代
        fprintf('第%d次迭代, 目标函数值为%f\n', i, obj_function(i));
    % 检查迭代终止条件
    if i > 1,
        if abs(obj_function(i) - obj_function(i-1)) < min_deviation
        break; 
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 该函数用于每次迭代过程
function [newU,center,obj_function] = FCMStep(data, U, cluster_num, m)
% data为被聚类数据,U为隶属度矩阵,cluster_num为聚类类别数,m为FCM中的参数m
% 函数调用后得到新的隶属度矩阵newU,聚类中心center,目标函数值obj_function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 以下是计算模糊隶属度Ut
    [x,y] = size(U);
    A = ones(x,y);
          a = 0.85;
            Ut = abs(A - U -(A - (U).^a).^(1/a));
            Ud = U  + Ut;
            [j,k,l] = size(data);
            pp =  y;
            pai = (sum(Ut,2)) ./pp;
            obj = sum(pai.*exp(1-pai));
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%            Ud = U;
%            obj = 0;
 nf = Ud;
mf = Ud.^m;       % FMC中的U^m
%  center = nf*data./((ones(size(data, 2), 1)*sum(nf'))'); % 得到聚类中心

   data1 = zeros(x,y);
   data1(1,:) = data';
   data1(2,:) = data';
   data1(3,:) = data';
   data1(4,:) = data';
%   data1(5,:) = data';
  center = sum(nf.*data1,2)./sum(nf,2); % 得到聚类中心

dist = Distance(center, data);      % 调用myfcmdist函数计算聚类中心与被聚类数据的距离
obj_function = sum(sum((dist.^2).*mf))+obj;  % 得到目标函数值
tmp = dist.^(-2/(m-1));      % 如果迭代次数不为1,计算新的隶属度矩阵
 newU = tmp./(ones(cluster_num, 1)*sum(tmp)); % U_new为新的隶属度矩阵

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Distance()函数用于计算聚类中心与被聚类数据的距离
% center为聚类中心，data为被聚类数据，输出各元素到聚类中心的距离out
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = Distance(center, data)
  
 data_number = size(data,1);
 class_number = size(center, 1);
 kk = ones(data_number,1); % 构造与数据大小相同的全1矩阵kk
 out = zeros(class_number, data_number);
  if size(center, 2) > 1, %若类别数大于1
      for k = 1:class_number
           out(k, :) = sqrt(sum(((data - kk...
                       *center(k,:)).^2)'));
     end
  else    % data为一维数据
      for k = 1:class_number
           out(k, :) = abs(center(k) - data)';
      end
  end