转载

The Robust Fuzzy C-means

摘要：

基于FCM的在图像处理方面对噪声敏感的不足，本文通过引入空间模型建立空间模糊C均值聚类提高算法的鲁棒性，在此基础上，结合抑制式对算法进一步优化。最后，给图像加不同程度的噪声，通过MATLAB编程，分析比较各个算法的迭代次数、迭代时间、错误率以及均方根误差判断算法的优劣。

一空间模糊聚类概述

之前对FCM做了简单介绍，通过实验发现其有一定的缺点，为此我们需要对算法进行一系列的优化，今天主要来分析一下空间模糊聚类，改善算法的鲁棒性。

1. FCM

在之前对FCM已经做过介绍，在此不再赘述，今天主要分析的是FCM在图像处理方面的应用，虽然FCM有着广泛适用性等优点，但FCM算法对噪声敏感这一点也是毋庸置疑的。FCM的目标函数：

The Robust Fuzzy C-means

其中，yj指的是各个样本点对应的像素，取值为[0,255]。很明显，FCM并没有考虑到样本之间的空间相关性，只是分别考虑各个样本。减小噪声敏感程度的方法很多，比如在处理图像之前首先平滑图像，但这个方法会在一定程度上丢失样本中某些特性；还可以通过对隶属度函数的后处理等等，然而，我们今天讲到的空间模型是通过直接修改目标函数提高算法对噪声的鲁棒性。

2. The Roubst Fuzzy C-means（RFCM）

结合空间模型，为提高鲁棒性，RFCM在计算距离时考虑到每个样本点的邻域样本，也就将各个样本的空间相关性结合起来。在FCM算法的基础之上，得到RFCM的目标函数：

The Robust Fuzzy C-means

结合FCM的约束条件，引入拉格朗日因子进一步推导，得到的隶属度矩阵的计算公式为：

The Robust Fuzzy C-means

相应的聚类中心的计算公式为：

3. 抑制式

抑制式主要是对最大隶属度进行适当的放大，从而抑制其它隶属度，提高收敛速度，其中抑制式因子为：

The Robust Fuzzy C-means

通过将隶属度函数与抑制式因子相乘，即可对算法进行抑制式修改。

二 MATLAB实现

为了比较算法的优劣，将算法作用于不同图像，并且分别加不同程度的高斯白噪声，分别比较算法的迭代次数，执行时间，错误率和均方根误差，这些在代码中均有说明。

1. RFCM

基于上面的介绍，可得到RFCM算法的步骤如下：

首先，给定一个由N个L维向量组成的数据集X以及所要分得的类别个数C，并且自定义隶属度矩阵

（1）初始化隶属度矩阵并归一

（2）根据聚类中心的计算公式更新聚类中心

（3）根据隶属度矩阵的计算公式更新隶属度

（4）不断的重复（2）直到收敛。

不知道大家有没有发现这和FCM的算法步骤大致是一样的，不同是在目标函数当中，当然体现在代码当中，而最终反映在算法的性能当中。

根据算法的步骤，采用MATLAB编程实现，在此导入MR图像，代码如下：

clear all clc; a=imread('MRI.jpg'); I=imnoise(a,'salt & pepper',0.03); figure(1);       imshow(I);title('加噪图像'); [height,width,c]=size(a); if c~=1  a=rgb2gray(a); end a=double(a); [row,column]=size(a);  data = a(:); data_n = size(data,1); beta = 80; 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 for o = 1:iter_num %迭代次数控制  %计算初始中心  center = zeros(cluster_num,1);  for u = 1:cluster_num   sum = 0;   sum1 = 0;     for h=1:height      for t=1:width      sum = sum + (membership(h,t,u).^m) * a(h,t);      sum1 = sum1 + membership(h,t,u).^m;      end     end     center(u) = sum/sum1;  end       for i=1:height       for j =1:width       up = i-1;       down = i+1;       left = j-1;       right = j+1;       if( up < 1)        up = 1;       end        if( down > height)        down = height;       end        if( left< 1)        left = 1;       end       if( right > width)        right = width;       end             s = 0;        for x = up : down          for y = left : right          for u = 1:cluster_num           for r = 1:cluster_num            s = s + membership(x,y,r).^m;           end           s = s - membership(x,y,u).^m;          end          end        end      end    end      for h = 1:height   for t = 1:width    for k = 1:cluster_num      costfunction(o) = costfunction(o) + membership(h,t,k).^m*(a(h,t) - center(k))^2 + (beta/2) * membership(h,t,k)^m*s;      tmp = ((a(h,t)-center(k))^2  + beta * s).^(-1/(m - 1));      tomp = 0;      for p = 1:cluster_num       tomp = tomp + (a(h,t)-center(p))^2;       tp = (tomp + beta * s) .^(-1/(m - 1));      end      membership(h,t,k) = tmp./tp;    end   end  end  %%%%%%归一化隶属度  for i=1:height    for j = 1:width    member_sum = 0;    for k = 1:cluster_num     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   if display,    fprintf('Iteration count = %d, obj. fcn = %f/n', o, costfunction(o));   %输出迭代次数和函数的结果  end  % check termination condition  if o > 1,  %进化步长控制   if abs(costfunction(o) - costfunction(o-1)) < threshold, break; end,  end end  toc      A = ones(height,width,1); for i = 1:height  for j = 1:width   if (fix(a(i,j) / 85) == 1)    A(i,j,1) = 2;   end   if (fix(a(i,j) / 85) == 2)    A(i,j,1) = 3;   end   if (fix(a(i,j,1) / 85) == 3)    A(i,j,1) = 4;   end  end end  A = reshape(A,1,data_n);   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   newing(i,j) = round(255 * (1-(index-1)/(cluster_num-1)));  end  end B = reshape(newing,1,data_n); b = fix((max(B) - B(1,1)) / cluster_num); for i = 2:data_n  if B(1,i) == B(1,1)   B(1,i) = 1;  elseif (fix(B(1,i) / b) == 2)   B(1,i) = 2;  elseif (fix(B(1,i) / b)  == 3)   B(1,i) = 3;  else   B(1,i) = 4;  end end B(1,1) = 1; sum = 0; for i = 1:data_n  if ( A(1,i) ~= B(1,i))   sum = sum + 1;  end end MCR = sum / data_n; fprintf('MCR = %d/n',MCR); S = 0; for i = 1:data_n  S = S + (A(1,i) - B(1,i)).^2; end RMS = sqrt(S / (data_n * (data_n - 1))); fprintf('RMS = %d/n',RMS); figure(2); imshow((uint8(newing))); title('RFCM分割后的图像');

2. 抑制式RFCM

将抑制式因子作用于隶属度矩阵，采用MATLAB编程实现，代码如下：

clear all clc; a=imread('MRI.jpg'); I=imnoise(a,'salt & pepper',0.03); figure(1);       imshow(I);title('加噪图像'); [height,width,c]=size(a); if c~=1  a=rgb2gray(a); end a=double(a); [row,column]=size(a);  data = a(:); data_n = size(data,1); beta = 80; 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 for o = 1:iter_num %迭代次数控制  %计算初始中心  center = zeros(cluster_num,1);  for u = 1:cluster_num   sum = 0;   sum1 = 0;     for h=1:height      for t=1:width      sum = sum + (membership(h,t,u).^m) * a(h,t);      sum1 = sum1 + membership(h,t,u).^m;      end     end     center(u) = sum/sum1;  end       for i=1:height       for j =1:width       up = i-1;       down = i+1;       left = j-1;       right = j+1;       if( up < 1)        up = 1;       end        if( down > height)        down = height;       end        if( left< 1)        left = 1;       end       if( right > width)        right = width;       end             s = 0;        for x = up : down          for y = left : right          for u = 1:cluster_num           for r = 1:cluster_num            s = s + membership(x,y,r).^m;           end           s = s - membership(x,y,u).^m;          end          end        end      end    end      for h = 1:height   for t = 1:width    for k = 1:cluster_num      costfunction(o) = costfunction(o) + membership(h,t,k).^m*(a(h,t) - center(k))^2 + (beta/2) * membership(h,t,k)^m*s;      tmp = ((a(h,t)-center(k))^2  + beta * s).^(-1/(m - 1));      tomp = 0;      for p = 1:cluster_num       tomp = tomp + (a(h,t)-center(p))^2;       tp = (tomp + beta * s) .^(-1/(m - 1));      end      membership(h,t,k) = tmp./tp;    end   end  end  %%%%%%归一化隶属度  for i=1:height    for j = 1:width    member_sum = 0;    for k = 1:cluster_num     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   %%%%%%%抑制式修改  ha = 0;  for i = 1:height   for j = 1:width    for p = 1:cluster_num    ha = ha + membership(i,j,p).^2;    end   end  end   alpha = ( cluster_num/(cluster_num - 1)) * (1 -(1/data_n) * ha);   membership1 = membership.*alpha;   for i = 1:height   for j = 1:width     [max_data,max_local] = max(membership(i,j,:));     temembership = alpha.* membership(i,j,max_local) + (1 - alpha);     membership1(i,j,max_local) = temembership;   end  end  membership = membership1 * alpha;    for i=1:height    for j = 1:width    member_sum = 0;    for k = 1:cluster_num     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  if display,    fprintf('Iteration count = %d, obj. fcn = %f/n', o, costfunction(o));   %输出迭代次数和函数的结果  end  % check termination condition  if o > 1,  %进化步长控制   if abs(costfunction(o) - costfunction(o-1)) < threshold, break; end,  end end  toc      A = ones(height,width,1); for i = 1:height  for j = 1:width   if (fix(a(i,j) / 85) == 1)    A(i,j,1) = 2;   end   if (fix(a(i,j) / 85) == 2)    A(i,j,1) = 3;   end   if (fix(a(i,j,1) / 85) == 3)    A(i,j,1) = 4;   end  end end  A = reshape(A,1,data_n);   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   newing(i,j) = round(255 * (1-(index-1)/(cluster_num-1)));  end  end B = reshape(newing,1,data_n); b = fix((max(B) - B(1,1)) / cluster_num); for i = 2:data_n  if B(1,i) == B(1,1)   B(1,i) = 1;  elseif (fix(B(1,i) / b) == 2)   B(1,i) = 2;  elseif (fix(B(1,i) / b)  == 3)   B(1,i) = 3;  else   B(1,i) = 4;  end end B(1,1) = 1; sum = 0; for i = 1:data_n  if ( A(1,i) ~= B(1,i))   sum = sum + 1;  end end MCR = sum / data_n; fprintf('MCR = %d/n',MCR); S = 0; for i = 1:data_n  S = S + (A(1,i) - B(1,i)).^2; end RMS = sqrt(S / (data_n * (data_n - 1))); fprintf('RMS = %d/n',RMS); figure(2); imshow((uint8(newing))); title('SRFCM分割后的图像');