%close all
clear all;clc
h=helpdlg('按左键选择分割区域,双击');
uiwait(h)
%程序中对图像特征最敏感的两个参数
%T=5;%邻域特征的域值，决定聚类中心的一个关键参数
r=300;%归为该聚类中心的最小距离
prompt={'最大循环次数:','邻域灰度阈值:'};
name='蚁群聚类分割';
numlines=1;
defaultanswer={'10','10'};
options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';
itermax1=inputdlg(prompt,name,numlines,defaultanswer,options);
itermax=str2num(itermax1{1});
T=str2num(itermax1{2});
%itermax=10; %
%load aa;%aa.mat里面储存了一幅灰度图象的矩阵数据
%bb=imcrop(bb);%从该图像中选中一块矩形区域，程序对选中区域进行处理，这样可以减少运行时间
%[x,map]=imread('Dorsal 10.4.gif');%这几句读入索引类型的图象
%grayimg=ind2gray(x,map);%把索引类型图像转换为灰度图象
dicom1=dicomread('33091247');
grayimg=dicom2gray(dicom1);
%grayimg=imadjust(dicom1);
figure;imshow(grayimg)
%grayimg=imresize(grayimg,0.5);%把图片缩小为原来大小的0.5倍
%figure;imshow(grayimg);%把读入的图片显示出来
set(gcf,'outerposition',get(0,'screensize'));%设置figure全屏显示
[ffimg,rect]=imcrop(grayimg);
set(gcf,'outerposition',get(0,'screensize'));%设置figure全屏显示
title('原始图像');
%h1 = zoom;
%set(h1,'Enable','on')
%figure;
%imshow(ffimg);
%w8=[1,1,1;1,-8,1;1,1,1];
%ffimg=ffimg-imfilter(ffimg,w8,'replicate');
%ffimg=grayimg;
%ff=imread('LENA256.BMP');%读入另一幅灰度图象对算法进行效果测试
%ffimg=imread('cameraman.tif');%处理256*256的图片
%ffimg=imresize(ffimg,0.5);
%cannyedgeimg=edge(ff,'canny');%对读入的待处理灰度图象进行canny边缘检测
%figure;imshow(cannyedgeimg)%显示边缘检测输出的二值边缘图像
%width=70;%设定一个正方形的区域的边长
%另一种选取图像中一个正方形区域的方法
%ff=ff(floor(a/2)-width:floor(a/2)+width-1,floor(b/2)-width:floor(b/2)+width-1);
%[a,b]=size(ff);%输出图像的大小
%figure;imshow(ff);%显示待分割的图像
ff=double(ffimg);%把图像数据类型由unit8转化为double类型，便于后面进行数值处理
ff = medfilt2(ff);%进行中值滤波处理，效果比原来在连续性这一点上要好
%figure;imshow(uint8(ff));
%w1=fspecial('laplacian');w2=w1';%采用拉普拉斯梯度算子
%w1=fspecial('log',3);w2=w1';%采用log梯度算子
%w1=fspecial('prewitt');w2=w1';%采用prewitt梯度算子
w1=[-1,-2,-1;0,0,0;1,2,1];w2=[-1,0,1;-2,0,2;-1,0,1];%采用sobel梯度算子
g=sqrt(imfilter(ff,w1).^2+imfilter(ff,w2).^2);%计算梯度图像
gmax=max(sum(g(:,3:size(ff,2)-2)))/(size(ff,1)-4);%取出最大梯度列的均值 注意: 由于边界采用了人为填充所以剃度很大应该忽略
fp=padarray(ff,[1,1],'replicate','both');%对边界进行填充 目的是为了进行3*3的8邻域特征提取
for i=2:size(ff,1)+1
for j=2:size(ff,2)+1
a=[ fp(i-1,j-1) ,fp(i-1,j) fp(i-1,j+1);fp(i,j-1) fp(i,j) fp(i,j+1); fp(i+1,j-1) fp(i+1,j) fp(i+1,j+1)];
ne(i-1,j-1)=length(find((a-fp(i,j))<T))-1;
end
end
%下面是测试ne中小于等于6对应的图象中的点所组成的图像边缘
% as1=size(ff,1);as2=size(ff,2);
% m=zeros(as1,as2);
% m=reshape(m,1,as1*as2);
% m(find(ne<=5))=1;
% m=reshape(m,as1,as2);
% figure;imshow(m)
% 初始化聚类中心
f=reshape(ff,1,size(ff,1)*size(ff,2));
h=hist(f,256);%灰度直方图图象点统计
counter=0;%计算出灰度直方图上的极大值点，以它们作为聚类中心
for i=2:length(h)-1
if h(i-1)<=h(i)&h(i)>=h(i+1)
counter=counter+1;
V(counter)=i;
end
end
h1=sum(h(V))/length(h(V));%取整体均值作为区分背景目标与边缘的区域值，即像素数远远大于该整体均值的灰度作为背景或者目标的灰度特征
h2=sum(h(find(h(V)<h1)))/length(find(h(V)<h1));%取非目标或者背景的所有像素数的平均值，大于该值的灰度作为边缘的灰度
for i=1:counter
if h(V(i))>5*h1
G(i)=0;
else
G(i)=gmax;
end
end
for i=1:counter
if G(i)>0
if h(V(i))>h2
Ne(i)=6;
else
Ne(i)=3;
end
else
Ne(i)=8;
end
end
g=reshape(g,1,size(f,1)*size(f,2));
ne=reshape(ne,1,size(f,1)*size(f,2));
X=[f;g;ne];%形成了待聚类的样品库，实际上就是图像中的所有象素点每个象素点 都有三维的特征向量 类似于聚类思想模式识别
C=[V;G;Ne];%形成了聚类中心的特征向量矩阵
%蚂蚁群算法基本参数
a=2;%信息素的指数 在很多文献中用 Alpha来表示
b=10;%启发因子的指数 一般用Beta表示
lamda=0.9;%信息素域值
rou=0.6 ;%信息素残留因子 一般用Rho表示
e=0.1 ;%最小类间距
ph=1*ones(size(X,2),size(C,2));%信息素矩阵
ph=double(ph);
disp('蚁群主循环开始计时')%显示循环开始并开始计时
tic;
h = waitbar(0,'Please wait...');%进度条设置
%基本蚁群算法 的大循环开始
for iter=1:itermax
iter
waitbar( iter/itermax,h)%进度条设置
t=zeros(1,size(X,2));qq=zeros(size(X,2),size(C,2)); p=zeros(size(X,2),size(C,2)); qq=double(qq);t=double(t);
M=zeros(size(X,2),size(C,2)); %初始化变量空间,注意释放变量以节省内存
%以下的程序来计算启发信息矩阵，信息素矩阵和转移概率矩阵
for j=1:size(C,2) %注意小循环一定要放在外面，运行速度能差好多倍
for i=1:size(X,2)
d(i,j)=sqrt(sum((X(:,i)-C(:,j)).^2));
%qq(i,j)=1./(d(i,j)+1);%另一种处理启发因子矩阵的方法，对于d(i,j)=0的点进行特殊处理
if d(i,j)==0
p(i,j)=1;
else
qq(i,j)=1./(d(i,j)); %形成启发信息矩阵
end
end
end
for j=1:size(C,2)
for i=1:size(X,2)
if d(i,j)<r
aa(i,j)=ph(i,j).^a*qq(i,j).^b;
t(i)=t(i)+ph(i,j).^a*qq(i,j).^b;
end
end
end
% p=(ph.^a).*((1./d).^b)/t(i);
for j=1:size(C,2)
for i=1:size(X,2)
if d(i,j)<r&d(i,j)>0
p(i,j)=(ph(i,j).^a*qq(i,j).^b)/t(i);
end
end
end
%M=p>0.9 %注意，最开始循环时候p大都很小,大于0.9的点很少，而
%2次3次循环之后p变得比较大，这里设置为0.9只是参考文献中给的参数实际实现时效果不理想
%M=p>lamda; %可以设定一个概率矩阵的域值
%取最大转移概率作为聚类中心，并把聚类情况存储在M矩阵中相应位置置为1，其标志位的作用
for j=1:size(X,2)
[aa6,a5]=max(p(j,:));
M(j,a5)=1;
end
ph=rou*ph+M;%更新信息素矩阵
%更新聚类中心
C1=C; %保存聚类情况,下面进行更新合并
for j=2:size(C,2)-1
a1=find(M(:,j)==1);
if length(a1)>0
C(:,j)=sum(X(:,a1),2)/length(a1);
end
end
C
counter3=0;
%动态更新各类之间的距离，直到没有相似的类中心为止 所谓动态是指 同时更新相应的信息素矩阵的维数
for i=1:size(C,2)-1
for j=i+1:size(C,2)
a2=(C(:,i)-C(:,j)).^2;
if sum(a2)<e
counter3=counter3+1;
ph(:,j)=(ph(:,i)+ph(:,j))/2; %同时更新相应的信息素矩阵的维数
C(:,i)=(C(:,i)+C(:,j))/2;
for m=j+1:size(C,2)
C(:,m-1)=C(:,m);
ph(:,m-1)=ph(:,m); %同时更新相应的信息素矩阵的维数
end
end
end
end
C=C(:,1:size(C,2)-counter3);
ph=ph(:,1:size(ph,2)-counter3);
% %判断还有没有还没分类的像素 记录到SS集合这是参考文献上的实现方法实际在本程序中未采纳
% counter1=0;
% for i=1:size(X,2)
% if sum(M(i,:))==0
% counter1=counter1+1;
% SS(counter1)=i;
% end
% end
% %if length(SS)==0
% % disp('wanle')
% % break
% % end
%如果不满足一条件则循环终止;这里是终止条件
% if min(max(p,[],2))>=lamda
% disp('满足终止条件第几次循环：');
% iter
% break;
% end
end%终结蚂蚁算法的大循环
%M=p>0.8; %此句话放在蚁群算法之后效果还可以接受
time=toc;%显示循环所用时间
disp(['蚁群循环结束,共用时',num2str(time),'秒'])
m=zeros(size(ff,1),size(ff,2));
m=reshape(m,1,size(ff,1)*size(ff,2));
% a3=find(C1(3,:)<=7);%C1是没有更新的聚类中心情况；C是更新并且合并了相似类的聚类情况
a3=find(C1(3,:)<=6&C1(3,:)>=3);
for i
