matlab实现基于DPCA密度峰值算法的图像分割

matlab代码如下
图像最好小一点，不然因为那个距离矩阵为（mn)(m*n)真的很慢很慢，运算上亿次了，电脑根本跑不动。

clear
close all

Iy=imread('1.jpg');
figure;imshow(Iy);
I = imresize(Iy,[53 33]); 
figure, imshow(I);  % 读入分割目标
img_gray=rgb2gray(I);%灰度化

% 获取图像的长宽
[m,n]=size(img_gray);

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

distance_all=[];
z=1;
pn=m*n;

for i =1: pn
        for j = 1:pn
            y=(dat(i,1)-dat(j,1))^2+(dat(i,2)-dat(j,2))^2+(dat(i,3)-dat(j,3))^2;
            y=double(y);
            dist(i,j)=sqrt(y);%把距离放入了矩阵           
            distance_all(z)= sqrt(y);%把距离放入一个长的地方进行排序
            z=z+1;
        end
end


%% 确定 dc
percent=2.0;
position=round(length(distance_all)*percent/100);   %% round 是一个四舍五入函数
sda=sort(distance_all);   %% 对所有距离值作升序排列
dc=sda(position);

%% 计算局部密度 rho (利用 Gaussian 核)
 
%% 将每个数据点的 rho 值初始化为零
for i=1:pn
  p(i)=0.;
end
% Gaussian kernel
%for i=1:pn-1
% for j=i+1:pn
%     p(i)=p(i)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));
%     p(j)=p(j)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));
%  end
%end 
% "Cut off" kernel
%
for i=1:pn-1
 for j=i+1:pn
   if (dist(i,j)<dc)
        p(i)=p(i)+1.;
        p(j)=p(j)+1.;
    end
  end
end

%% 先求矩阵列最大值，再求最大值，最后得到所有距离值中的最大值
maxd=max(max(dist));
%% 将 rho 按降序排列，ordrho 保持序
[rho_sorted,ord]=sort(p,'descend');
%% 处理 rho 值最大的数据点
delta(ord(1))=-1.;
nneigh(ord(1))=0;

%% 生成 delta 和 nneigh 数组
for i=2:pn
   delta(ord(i))=maxd;
   for j=1:i-1
     if(dist(ord(i),ord(j))<delta(ord(i)))
        delta(ord(i))=dist(ord(i),ord(j));
        nneigh(ord(i))=ord(j);
        % 记录 rho 值更大的数据点中与 ordrho(ii) 距离最近的点的编号 ordrho(jj)
     end
   end
end
 
%% 生成 p 值最大数据点的 delta 值
delta(ord(1))=max(delta(:));
 
%% 决策图

fid = fopen('DECISION_GRAPH', 'w');
for i=1:pn
   fprintf(fid, '%6.2f %6.2f\n', p(i),delta(i));
end
 
%% 选择一个围住类中心的矩形0

scrsz = get(0,'ScreenSize');
 
%% 人为指定一个位置
figure('Position',[6 72 scrsz(3)/4. scrsz(4)/1.3]);
 
%% ind 和 gamma 在后面并没有用到
for i=1:pn
  ind(i)=i;
  gamma(i)=p(i)*delta(i);
end
 
%% 利用 rho 和 delta 画出一个所谓的“决策图”
subplot(1,1,1)
tt=plot(p(:),delta(:),'.','MarkerSize',5,'MarkerFaceColor','k','MarkerEdgeColor','k');
axis([0,max(p),0,max(delta)]);

fig=subplot(1,1,1);
rect = getrect(fig);
 
%% getrect 从图中用鼠标截取一个矩形区域， rect 中存放的是
%% 矩形左下角的坐标 (x,y) 以及所截矩形的宽度和高度
rhomin=rect(1);
deltamin=rect(2); %% 作者承认这是个 error，已由 4 改为 2 了!
 
%% 初始化 cluster 个数
ncluster=0;
 
%% cl 为归属标志数组，cl(i)=j 表示第 i 号数据点归属于第 j 个 cluster
%% 先统一将 cl 初始化为 -1
for i=1:pn
  cl(i)=-1;
end
 
%% 在矩形区域内统计数据点（即聚类中心）的个数
for i=1:pn
  if ( (p(i)>rhomin) && (delta(i)>deltamin))
     ncluster=ncluster+1;
     cl(i)=ncluster;  %% 第 i 号数据点属于第 NCLUST 个 cluster
     icl(ncluster)=i; %% 逆映射,第 NCLUST 个 cluster 的中心为第 i 号数据点 
  end
end
fprintf('总共有多少类: %i \n', ncluster);
 
%% 将其他数据点归类 (assignation)
for i=1:pn
  if (cl(ord(i))==-1)
    cl(ord(i))=cl(nneigh(ord(i)));
  end
end
%halo
%% 由于是按照 p值从大到小的顺序遍历,循环结束后, cl 应该都变成正的值了. 
for i=1:pn
  halo(i)=cl(i);
end
if (ncluster>1)
     % 初始化数组 bord_rho 为 0,每个 cluster 定义一个 bord_rho 值
  for i=1:ncluster
    bord_rho(i)=0.;
  end
  % 获取每一个 cluster 中平均密度的一个界 bord_rho
  for i=1:pn-1
    for j=i+1:pn
        %% 距离足够小但不属于同一个 cluster 的 i 和 j
      if ((cl(i)~=cl(j))&& (dist(i,j)<=dc))
        rho_aver=(p(i)+p(j))/2.;  %% 取 i,j 两点的平均局部密度
        if (rho_aver>bord_rho(cl(i))) 
          bord_rho(cl(i))=rho_aver;
        end
        if (rho_aver>bord_rho(cl(j))) 
          bord_rho(cl(j))=rho_aver;
        end
      end
    end
  end
   
%% halo 值为 0 表示为 outlier
  for i=1:pn
    if (p(i)<bord_rho(cl(i)))
      halo(i)=0;
    end
  end
end
 
%% 逐一处理每个 cluster
for i=1:ncluster
  nc=0;  %% 用于累计当前 cluster 中数据点的个数
  nh=0;  %% 用于累计当前 cluster 中核心数据点的个数
  for j=1:pn
    if (cl(j)==i) 
      nc=nc+1;
    end
    if (halo(j)==i) 
      nh=nh+1;
    end
  end
  fprintf('类: %i CENTER: %i 数据点个数: %i 核心点: %i 边界点: %i \n', i,icl(i),nc,nh,nc-nh);
end
cmap=colormap;


for i=1:ncluster
   ic=int8((i*64.)/(ncluster*1.));
   subplot(1,1,1)
   hold on
   plot(p(icl(i)),delta(icl(i)),'o','MarkerSize',8,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
   axis([0,max(p),0,max(delta)]);
end

pixel_labels = reshape(cl,m,n);
rgb_labels = label2rgb(pixel_labels);
figure, imshow(rgb_labels), title('Segmented Image');

%for i=1:ND
%   if (halo(i)>0)
%      ic=int8((halo(i)*64.)/(NCLUST*1.));
%      hold on
%      plot(Y1(i,1),Y1(i,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
%   end
%end
faa = fopen('CLUSTER_ASSIGNATION', 'w');

for i=1:pn
   fprintf(faa, '%i %i %i\n',i,cl(i),halo(i));
end
 
clt=transpose(cl);
eva= evalclusters(dat,clt,'DaviesBouldin');
disp(eva);