function [center, U, obj_fcn] = myfcm(data, cluster_n, options)
%data为聚类数据,cluster_n为类别数
%判断函数参数个数
if nargin ~= 2 & nargin ~= 3,
error('参数太多或太少');
end
data_n = size(data, 1); %data第一维大小
in_n = size(data, 2); %data第二维大小
% 可以改变下列数据设置默认的options参数
default_options = [2; % FCM中的参数m
100; % 最大迭代次数
1e-5; % 判别误差epsonal
1]; % 迭代过程中显示
if nargin == 2,
options = default_options;
else
% 如果"options" 没有完全指定,则以default_options中的部分进行补充
if length(options) < 4,
tmp = default_options;
tmp(1:length(options)) = options;
options = tmp;
end
% 如果options中的一些参数不是数字,则以default_options中的代替
nan_index = find(isnan(options)==1); %找到options中的非数字元素下标
options(nan_index) = default_options(nan_index);
if options(1) <= 1
error('FCM中参数m应大于1!');
end
end
expo = options(1); % expo存放参数m
max_iter = options(2); % 最大迭代次数
min_impro = options(3); % 判别误差
display = options(4); % 迭代过程是否显示
obj_fcn = zeros(max_iter, 1); % obj_fcn为目标函数值数组
U = myfcminit (cluster_n, data_n); %调用myfcminit函数初始化隶属度矩阵U
for i = 1:max_iter
[U, center, obj_fcn(i)] = myfcmstep(data, U, cluster_n, expo); %调用myfcmstep函数进行一次迭代
if display==1
fprintf('迭代次数 = %d, obj. fcn = %f\n', i, obj_fcn(i));
end
% 检查迭代终止条件
if i > 1,
if abs(obj_fcn(i) - obj_fcn(i-1)) < min_impro
break;
end
end
end
iter_n = i; % iter_n为实际迭代次数
obj_fcn(iter_n+1:max_iter) = [];
function U = myfcminit (cluster_n, data_n)
%
% U为隶属度矩阵,大小为:类别数*被聚类数据元素个数
%参数cluster_n为聚类类别数
%data_n为被聚类数据的元素个数
%
U = rand(cluster_n, data_n); %随机生成隶属度矩阵U
col_sum = sum(U);%计算U中每列元素和
U = U./col_sum(ones(cluster_n, 1), :); %保证U中每列的隶属度值之和为1
function [U_new, center, obj_fcn] = myfcmstep (data, U, cluster_n, expo)
%
% data为被聚类数据,U为隶属度矩阵,cluster_n为聚类类别数,expo为FCM中的参数m
% 函数调用后得到新的隶属度矩阵U_new,聚类中心center,目标函数值obj_fcn
%
mf = U.^expo; % FMC中的U^m
center = mf*data./((ones(size(data, 2), 1)*sum(mf'))'); % 得到聚类中心
dist = myfcmdist (center, data); % 调用myfcmdist函数计算聚类中心与被聚类数据的距离
obj_fcn = sum(sum((dist.^2).*mf)); % 得到目标函数值
tmp = dist.^(-2/(expo-1)); % 如果迭代次数不为1,计算新的隶属度矩阵
U_new = tmp./(ones(cluster_n, 1)*sum(tmp)); % U_new为新的隶属度矩阵
function out =myfcmdist (center, data)
%
% center为聚类中心,data为被聚类数据
%
out = zeros(size(center, 1), size(data, 1));
if size(center, 2) > 1, %若聚类中心数据大于1,即聚类类别>=2
for k = 1:size(center, 1),
out(k, :) = sqrt(sum(((data-ones(size(data, 1), 1)*center(k, :)).^2)'));
end
else % data为一维数据
for k = 1:size(center, 1),
out(k, :) = abs(center(k)-data)';
end
end
function myfcmseg(file, cluster_n)
% file: 图像文件.
% cluster_n: 聚类类别数.
eval(['info=imfinfo(''',file, ''');']); %获取文件信息
switch info.ColorType
case 'truecolor' %真彩色图像
eval(['RGB=imread(''',file, ''');']);
I = rgb2gray(RGB);
clear RGB;
case 'indexed' %索引图像
eval(['[X, map]=imread(''',file, ''');']);
I = ind2gray(X, map);
clear X;
case 'grayscale' %灰度图像
eval(['I=imread(''',file, ''');']);
end;
I = im2double(I);
filename = file(1 : find(file=='.')-1); %获取主文件名
data = reshape(I, numel(I), 1); %将图像数据转为一列数据
tic %计时开始
options=[2; % FCM中的参数m
100; % 最大迭代次数
1e-5; % 判别误差epsonal
1];
[center, U, obj_fcn]=myfcm(data, cluster_n,options); %调用myfcm函数进行聚类
elapsedtime = toc; %聚类结束,计时结束,elapsedtime为进行聚类,所花时间
%聚类所获得数据保存在filename.mat文件中
%eval(['save(', filename, int2str(cluster_n),'.mat'', ''center'', ''U'', ''obj_fcn'', ''elapsedtime'');']);
%fprintf('聚类耗时 = %d', elapsedtime);
maxU=max(U);
temp = sort(center);
for n = 1:cluster_n; %按聚类结果分割图像
eval(['cluster',int2str(n), '_index = find(U(', int2str(n), ',:) == maxU);']);
index = find(temp == center(n));
switch index
case 1
color_class = 0;
case cluster_n
color_class = 255;
otherwise
color_class = fix(255*(index-1)/(cluster_n-1));
end
eval(['I(cluster',int2str(n), '_index(:))=', int2str(color_class),';']);
end;
I = mat2gray(I);
%eval(['imwrite(I,', filename,'_seg', int2str(cluster_n), '.bmp'');']); %保存生成的分割图像
imshow(I);%显示生成的分割的图像
