转自：直觉模糊C均值聚类与图像阈值分割 - liyuefeilong的专栏 - CSDN博客 https://blog.csdn.net/liyuefeilong/article/details/43816495
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 主函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function main
ima = imread('MR6.jpg');
% 先设定FCM的几个初始参数
options=[; % FCM公式中的参数m
; % 最大迭代次数
1e-]; % 目标函数的最小误差
class_number = ; % 分为4类
imt = ImageSegmentation(ima,class_number,options)
subplot(,,),imshow(ima),title('原图');
subplot(,,),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,); %将矩阵元素转换为一列数据
[center, U] = FCMprocess(data,cluster_n,options); %调用FCMData函数进行聚类
% 对于每个元素对不同聚类中心的隶属度，找出最大的那个隶属度
maxU = max(U); % 找出每一列的最大隶属度
temp = sort(center);
for i = :cluster_n; % 按聚类结果分割图像
% 前面求出每个元素的最大隶属度，属于各聚类中心的元素坐标，并存放这些坐标
% 调用eval函数将括号里的字符串转化为命令执行
eval(['class_',int2str(i), '= find(U(', int2str(i), ',:) == maxU);']);
%gray = round( * (i-) / (cluster_n-));
index = find(temp == center(i));
switch index
case
gray = ;
case cluster_n
gray = ;
otherwise
gray = fix(*(index-)/(cluster_n-));
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(); % 参数m
max_iteration = options(); % 最终的迭代次数
min_deviation = options(); % 最小判别误差
data_number = size(data, ); % 元素个数
obj_function = zeros(max_iteration, ); % obj_function用于存放目标函数的值
% 生成隶属度矩阵U
U = rand(cluster_num, data_number); % 随机生成隶属度矩阵U
sumU = sum(U,); % 计算U中每列元素和
for k = :data_number
U(:,k) = U(:,k) ./ sumU(k); % 对隶属矩阵U进行归一化处理
end for i = :max_iteration
[U, center, obj_function(i)] = FCMStep(data, U, cluster_num, m); %调用FCMStep函数进行迭代
fprintf('第%d次迭代, 目标函数值为%f\n', i, obj_function(i));
% 检查迭代终止条件
if i > ,
if abs(obj_function(i) - obj_function(i-)) < 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).^(/a));
Ud = U + Ut;
[j,k,l] = size(data);
pp = y;
pai = (sum(Ut,)) ./pp;
obj = sum(pai.*exp(-pai));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ud = U;
% obj = ;
nf = Ud;
mf = Ud.^m; % FMC中的U^m
% center = nf*data./((ones(size(data, ), )*sum(nf'))'); % 得到聚类中心 data1 = zeros(x,y);
data1(,:) = data';
data1(,:) = data';
data1(,:) = data';
data1(,:) = data';
% data1(,:) = data';
center = sum(nf.*data1,)./sum(nf,); % 得到聚类中心 dist = Distance(center, data); % 调用myfcmdist函数计算聚类中心与被聚类数据的距离
obj_function = sum(sum((dist.^).*mf))+obj; % 得到目标函数值
tmp = dist.^(-/(m-)); % 如果迭代次数不为1,计算新的隶属度矩阵
newU = tmp./(ones(cluster_num, )*sum(tmp)); % U_new为新的隶属度矩阵 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Distance()函数用于计算聚类中心与被聚类数据的距离
% center为聚类中心，data为被聚类数据，输出各元素到聚类中心的距离out
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = Distance(center, data) data_number = size(data,);
class_number = size(center, );
kk = ones(data_number,); % 构造与数据大小相同的全1矩阵kk
out = zeros(class_number, data_number);
if size(center, ) > , %若类别数大于1
for k = :class_number
out(k, :) = sqrt(sum(((data - kk...
*center(k,:)).^)'));
end
else % data为一维数据
for k = :class_number
out(k, :) = abs(center(k) - data)';
end
end