matlab实现基于超像素的快速的模糊聚类(SFFCM)的图像分割

1.SFFCM 基于superpxiels的FCM
超像素技术在提高图像分割的计算机效率时扮演很重要的作用，超像素意味着图像被划分为很多小一级有着不同形状和大小的独立的区域。基于一个图像的超像素结果，可以使用一个像素去代替超像素区域里所有的像素点来有效减少图片像素的数量。受这个的激发，Lei（T. Lei, X. Jia, Y. Zhang, S. Liu, H. Meng and A. K. Nandi, “Superpixel-based fast fuzzy c-means clustering for color image segmentation”, IEEE Trans. Fuzzy Syst., 2018.）提出来颜色图像分割的SFFCM。SFFCM的目标函数是

根据上式，隶属度矩阵和聚类中心给出如下

2.代码实现如下，超像素分割调用的是matlab自带的方法SLIC

%% SFFCM将超像素区域N*3，减少计算成本进行分割的结果

clear
close all

I=imread('tomato.jpg');
figure;subplot(2,3,1);imshow(I);
[m,n,p]=size(I);

A = reshape(I(:, :, 1), m*n, 1);    % 将RGB分量各转为kmeans使用的数据格式n行，一样一样本
B = reshape(I(:, :, 2), m*n, 1);
C = reshape(I(:, :, 3), m*n, 1);
data = [A B C];
dat=mat2gray(data);%归一化

N1=200;
[label,N] = superpixels(I,N1);%N为实际的超像素个数 

subplot(2,3,2);
BW = boundarymask(label);
imshow(imoverlay(I,BW,'cyan'),'InitialMagnification',67)

outputImage = zeros(size(I),'like',I);
output = zeros(N,3);
sumo=zeros(N,3);
idx = label2idx(label);%将标签矩阵L描述的区域转换为线性索引pixelIndexList
for labelVal = 1:N
    redIdx = idx{labelVal};
    greenIdx = idx{labelVal}+m*n;
    blueIdx = idx{labelVal}+2*m*n;
    sumo(labelVal,:)=[sum(I(redIdx)),sum(I(greenIdx)),sum(I(blueIdx))];
    output(labelVal,:) = [mean(I(redIdx)),mean(I(greenIdx)),mean(I(blueIdx))];
    outputImage(redIdx) = mean(I(redIdx));
    outputImage(greenIdx) = mean(I(greenIdx));
    outputImage(blueIdx) = mean(I(blueIdx));
end    
subplot(2,3,3);
imshow(outputImage,'InitialMagnification',67)

label1=reshape(label,m*n,1); 
S=tabulate(label1);
Sl=S(:,2);%得到每个超像素区域的个数

k=4;%分类数
b=2;%加权指数
T=100;%迭代次数
epsm=1e-3;%目标函数的最小误差

U=randn(N,k);%随机得到每个点与聚类中心的隶属度矩阵
Obj_pre=inf;%初始的设置大一点，无穷大
Obj=[];
iter=0;

while(iter<T)
    iter=iter+1;
    Um=U.^b;
    C=(Um'*sumo)./((Um'*Sl)*ones(1,p));%C应该是一个k*3的矩阵,sum(Um)是压缩了变成了1*k
    %dist距离矩阵应该是一个N*k的矩阵
 %   for j=1:k
 %       for i=1:N
 %           dist(i,j)=sum((X(i,:)-C(j,:)).^2);
 %       end
 %   end
    dist = sum(output.*output, 2)*ones(1, k) + (sum(C.*C, 2)*ones(1, N))'-2*output*C';%dist为N*k的矩阵
    t=(1./dist).^(1/(b-1));
    %隶属度矩阵应该是一个N*K的矩阵
    U=t./(sum(t,2).*ones(1,k));
    %Obj_cur=sum(sum(Sl'*(U.^b).*dist))/N;
    Obj_cur=sum(Sl'*((U.^b).*dist));
    Obj = [Obj Obj_cur];
    if norm(Obj_cur-Obj_pre, 'fro') < epsm
        break;
    end
    fprintf('迭代次数: %03d, 目标函数为: %f\n', iter, Obj_cur);
    Obj_pre = Obj_cur;
end
%[~, label] = min(dist, [], 2);%dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最小值。
[~, label2] = max(U, [], 2);%dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最小值。
subplot(2,3,6);plot(1:length(Obj), Obj, 'r-*'), xlabel('迭代次数'), ylabel('目标函数')

label3=[];
%将超像素的分割拓展到原来的图像
for labelVal=1:N
    id=idx{labelVal};%得到的是这个超像素得到坐标点
    label3(id)=label2(labelVal);%该坐标点都设置为一个分类
end
r2 = reshape(label3, m, n);     % 反向转化为图片形式
subplot(2,3,4);imshow(label2rgb(r2))   % 显示分割结果

subplot(2,3,5);
BW1 = boundarymask(r2);
imshow(imoverlay(I,BW1,'red'),'InitialMagnification',67)