基于核方法的模糊C均值聚类

　　　　

　　　　　　

　

　　　　　　

　　　　　　

，根据

，

　　　　　　

　　　　

　　

%%%%%MKFCM
clear all
clc;

a=imread('MRI.jpg');
I=imnoise(a,'salt & pepper',0.03);

[height,width,c]=size(a);
if c~=1
a=rgb2gray(a);
end

a=double(a);
[row,column]=size(a);

%%  几个参数
beta = 1000;
cluster_num = 4;%将图像分为四类
default_options = [2.0;    % exponent for the partition matrix U
300;    % max. number of iteration
0.01;    % min. amount of improvement
1];    % info display during iteration

options = default_options;

m = options(1);        % p，图像的质量参数
iter_num = options(2);        % Max. iteration 最大迭代次数
threshold = options(3);        % Min. improvement  最小进化步长
display = options(4);        % Display info or not 显示信息与否
costfunction = zeros(iter_num, 1);
%%%%%%%%%%%%%%%%%%%%%%%

tic
%%  初始化隶属度并归一
membership = zeros(height,width,cluster_num);
for i=1:height
for j=1:width
member_sum=0;
for k=1:cluster_num
membership(i,j,k)=rand();
member_sum=member_sum + membership(i,j,k);
end
for p =1:cluster_num
membership(i,j,p)=membership(i,j,p)/member_sum;
end
end
end

center = zeros(cluster_num,1);
for o = 1:iter_num %迭代次数控制

%% 计算初始中心
for u = 1:cluster_num
sum = 0;
sum1 = 0;
for h=1:height
for t=1:width
sum =  sum + membership(h,t,u).^m * exp((-1) * (a(h,t) - center(u).^2) / (beta.^2)) * a(h,t);
sum1 = sum1 + membership(h,t,u).^m  * exp((-1) * ((a(h,t) - center(u))^2) / (beta.^2));
end
end
center(u) = sum / sum1;
end

for h = 1:height
for t = 1:width
for k = 1:cluster_num
costfunction(o) = costfunction(o) + 2 * membership(h,t,k).^m * (1 - exp( -1 * ((a(h,t) - center(k)).^2)/(beta.^2)));
temp =(1 - exp( -1 * ((a(h,t) - center(k)).^2)/(beta.^2))).^(-1 / (m -1));
top = 0;
for p = 1:cluster_num
top = top + (1 - exp((-1 * (a(h,t) - center(p)).^2) /(beta.^2))).^(-1 / (m -1));
end
membership(h,t,k) = temp / top;
end
end
end

if display,
fprintf('iter.count = %d, obj. fcn = %f\n',o, costfunction(o));
end
if o >1
if abs(costfunction(o)-costfunction(o-1)) < threshold , break;end;
end

end
toc
%% 输出图像
newing = zeros(row,column);
for i=1:row
for j=1:column
maxmembership=membership(i,j,1);
index=1;
for k=2:cluster_num
if(membership(i,j,k)>maxmembership)
maxmembership=membership(i,j,k);
index=k;
end
end
newimg(i,j) = round(255 * (1-(index-1)/(cluster_num-1)));
end
end

subplot(2,2,1),imshow(uint8(a)); title('原图')
subplot(2,2,2);imshow(I);title('加噪图像');
subplot(2,2,3);imshow(uint8(newimg));title('KFCM分割后的图像');

　　　　　　

%% 邻域信息
Xk = zeros(height,width);
for i = 1:height
for j = 1:width
up = i - 1;
down = i + 1;
left = j - 1;
right = j + 1;
if (up < 1)
up = i;
end
if (down > height)
down = height;
end
if (left < 1)
left = width;
end
if (right > width)
right = height;
end
A = zeros(1,6);
for x = up : down
for y = left : right
A = a(x,y);
end
end
Xk(i,j) = mean(A);
end
end