模糊聚类的图像分割实验报告

实验一基于模糊聚类的图像分割

一，实验目的

通过模糊c-均值（FCM）聚类实现图像的分割。

二，算法描述

动态聚类方法的目的是把n个样本划分到c个类别中的一个，使各样本与其所在类均值的误差平方和最小。FCM聚类算法的目标函数为：

Min 错误！未找到引用源。(U,Z) =

（1）

其中m>1是一个可以控制聚类结果的模糊程度的常数。在不同的隶属度定义方法下最小化式（1）的损失函数，就得到不同的模糊聚类方法。其中最有代表性的模糊C均值方法，要求一个样本对于各个聚类的隶属度之和为1，即

（2）

在条件（2）下求式（1）的最小值，令错误！未找到引用源。对错误！未找到引用源。

和错误！未找到引用源。的偏导数为0，可得必要条件：

错误！未找到引用源。（3）

三，变量说明

P 数据样本维数（灰度图像时为1）；

N 像素点数目；

X 像素i特征（灰度图像时，表示灰度值）；

C 图像分割类别数；

U 像素点i属于第j类的隶属度；

Z 第i类聚类中心；

四, 算法步骤

Step1:设置目标函数精度ε，模糊指数m（m通常取2），最大迭代次数错误！未找到引用源。;

Step2:初始化模糊聚类中心错误！未找到引用源。;

Step3:由（3）式更新模糊划分矩阵U={错误！未找到引用源。}和聚类中心Z={错误！未找到引用源。};

Step4:若|J(t)-J(t-1)|< ε或c>错误！未找到引用源。则结束聚类；否则，t=t+1并转Step3;

Step5:由所得U={错误！未找到引用源。}得到各像素点的分类结果。

五，实验内容与要求

（1）使用附录1的参考程序对无噪图像进行模糊聚类分割。

（2）使用附录1的参考程序对各种加噪（高斯噪声，椒盐噪声及斑点噪声等）图像进行模糊聚类分割，并与（1）中的相应结果进行比较。

（3）附录1的参考程序给出了图像分割为3类的FCM算法，请同学们进行分割为2类或4类的扩展，或者在理解例程或算法的基础上自己实现算法。

参考程序代码：

function fcm

tmp=imread('C:\Documents and Settings\Administrator\桌面\fig8.pgm'); IM=tmp(:,:,1);

%IM=imnoise(IM,'speckle',0.01);

%IM=imnoise(IM,'salt & pepper',0.1);

%IM=imnoise(IM,'gaussian',0,0.01);

IM=double(IM);

figure(1);

imshow(uint8(IM));

[maxX,maxY]=size(IM);

IMM=cat(3,IM,IM,IM);

cc1=8;

cc2=100;

cc3=200;

ttfcm=0;

while(ttfcm<15)

ttfcm=ttfcm+1;

c1=repmat(cc1,maxX,maxY);

c2=repmat(cc2,maxX,maxY);

c3=repmat(cc3,maxX,maxY);

c=cat(3,c1,c2,c3);

ree=repmat(0.000001,maxX,maxY);

ree1=cat(3,ree,ree,ree)

distance=IMM-c;

distance=distance.*distance+ree1;

daoShu=1./distance;

daoShu2=daoShu(:,:,1)+daoShu(:,:,2)+daoShu(:,:,3);

distance1=distance(:,:,1).*daoShu2;

u1=1./distance1;

distance2=distance(:,:,2).*daoShu2;

u2=1./distance2;

distance3=distance(:,:,3).*daoShu2;

u3=1./distance3;

ccc1=sum(sum(u1.*u1.*IM))/sum(sum(u1.*u1));

ccc2=sum(sum(u2.*u2.*IM))/sum(sum(u2.*u2));

ccc3=sum(sum(u3.*u3.*IM))/sum(sum(u3.*u3));

tmpMatrix=[abs(cc1-ccc1)/cc1,abs(cc2-ccc2)/cc2,abs(cc3-ccc3)/cc3]; pp=cat(3,u1,u2,u3);

for i=1:maxX

for j=1:maxY

if max(pp(i,j,:))==u1(i,j) ix2(i,j)=1;

elseif max(pp(i,j,:))==u2(i,j) ix2(i,j)=2;

else

ix2(i,j)=3;

end

end

end

if max(tmpMatrix)<0.0001

break;

else

cc1=ccc1;

cc2=ccc2;

cc3=ccc3;

end

for i=1:maxX

for j=1:maxY

if ix2(i,j)==3

IMMM(i,j)=240;

elseif ix2(i,j)==2

IMMM(i,j)=130;

else

IMMM(i,j)=20;

end

end

end

figure(2);

imshow(uint8(IMMM));

end

for i=1:maxX

for j=1:maxY

if ix2(i,j)==3

IMMM(i,j)=240;

elseif ix2(i,j)==2

IMMM(i,j)=130;

else

IMMM(i,j)=20;

end