转自：直觉模糊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 ', 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