Auto_set.zip_高分辨率

%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