圖像模糊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