下面是将density peaks密度峰值算法用于图像处理的改进代码，出现错误，希望大神帮忙看一下怎么改正。
- clear all
- close all
- % disp('The only input needed is a distance matrix file')
- % disp('The format of this file should be: ')
- % disp('Column 1: id of element i')
- % disp('Column 2: id of element j')
- % disp('Column 3: dist(i,j)')
- % mdist=input('name of the distance matrix file (with single quotes)?n');
- % disp('Reading input distance matrix')
- % xx=load(mdist);
- % xx = load('D:example_distances.dat')
- %x = load('C:UseTraceOfAllUsers.txt');
- % 从文件中读取数据
- x = imread('naochuxue.jpg');
- x=double(x);
- minX = min(x);
- maxX = max(x);%取较大值
- ran = maxX - minX;
- nx(:,1) = (x(:,1) - minX(1,1)) / ran(1,1);
- nx(:,2) = (x(:,2) - minX(1,2)) / ran(1,2);
- dist = pdist2(nx, nx);
- N = size(dist,1);%%第一个维度的长度，相当于文件的行数（即距离的总个数）
- xx = zeros((N-1)*N/2, 3);%初始化为零
- idx = 1;
- % 这里不考虑对角线元素
- for i=1:N
- for j=i+1:N
- xx(idx, 1) = i;
- xx(idx, 2) = j;
- xx(idx, 3) = dist(i, j);
- idx = idx + 1;
- end
- end
- N = size(xx, 1);
- ND=max(xx(:,2));
- NL=max(xx(:,1));
- if (NL>ND)
- ND=NL;
- end
- % N=size(xx,1);
- % for i=1:ND
- % for j=1:ND
- % dist(i,j)=0;
- % end
- % end
- % for i=1:N
- % ii=xx(i,1);
- % jj=xx(i,2);
- % dist(ii,jj)=xx(i,3);
- % dist(jj,ii)=xx(i,3);
- % end
- % percent=2;
- % fprintf('average percentage of neighbours (hard coded): %5.6fn', percent);
- %
- % position=round(N*percent/100);
- % sda=sort(xx(:,3));
- % dc=sda(position);
- dc = 0.15;
- %计算局部密度 rho (利用 Gaussian 核)
- fprintf('Computing Rho with gaussian kernel of radius: %12.6fn', dc);
- % 将每个数据点的 rho 值初始化为零
- for i=1:ND
- rho(i)=0.;
- end
- %
- % Gaussian kernel
- %
- for i=1:ND-1
- for j=i+1:ND
- rho(i)=rho(i)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));
- rho(j)=rho(j)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));
- end
- end
- %
- % "Cut off" kernel
- %
- % for i=1:ND-1
- % for j=i+1:ND
- % if (dist(i,j)
- % rho(i)=rho(i)+1.;
- % rho(j)=rho(j)+1.;
- % end
- % end
- % end
- % 先求矩阵列最大值，再求最大值，最后得到所有距离值中的最大值
- maxd=max(max(dist));
- % 将 rho 按降序排列，ordrho 保持序
- [rho_sorted,ordrho]=sort(rho,'descend');
- % 处理 rho 值最大的数据点
- delta(ordrho(1))=-1.;
- nneigh(ordrho(1))=0;
- % 生成 delta 和 nneigh 数组
- for ii=2:ND
- delta(ordrho(ii))=maxd;
- for jj=1:ii-1
- if(dist(ordrho(ii),ordrho(jj))
- delta(ordrho(ii))=dist(ordrho(ii),ordrho(jj));
- nneigh(ordrho(ii))=ordrho(jj);
- % 记录 rho 值更大的数据点中与 ordrho(ii) 距离最近的点的编号 ordrho(jj)
- end
- end
- end
- % 生成 rho 值最大数据点的 delta 值
- delta(ordrho(1))=max(delta(:));
- % 决策图
- disp('Generated file:DECISION GRAPH')
- disp('column 1:Density')
- disp('column 2:Delta')
- fid = fopen('DECISION_GRAPH', 'w');
- for i=1:ND
- fprintf(fid, '%6.2f %6.2fn', rho(i),delta(i));
- end
- % 选择一个围住类中心的矩形
- disp('Select a rectangle enclosing cluster centers')
- scrsz = get(0,'ScreenSize');
- figure('Position',[6 72 scrsz(3)/4. scrsz(4)/1.3]);
- for i=1:ND
- ind(i)=i;
- gamma(i)=rho(i)*delta(i);
- end
- subplot(2,1,1)
- tt=plot(rho(:),delta(:),'o','MarkerSize',5,'MarkerFaceColor','k','MarkerEdgeColor','k');
- title ('Decision Graph','FontSize',15.0)
- xlabel ('rho')
- ylabel ('delta')
- % 利用 rho 和 delta 画出一个决策图
- subplot(2,1,1)
- rect = getrect(1);
- rhomin=rect(1);
- deltamin=rect(4);
- % 初始化 cluster 个数
- NCLUST=0;
- % cl 为归属标志数组，cl(i)=j 表示第 i 号数据点归属于第 j 个 cluster
- % 将 cl 初始化为 -1
- for i=1:ND
- cl(i)=-1;
- end
- % 在矩形区域内统计数据点（即聚类中心）的个数
- for i=1:ND
- if ( (rho(i)>rhomin) && (delta(i)>deltamin))
- NCLUST=NCLUST+1;
- cl(i)=NCLUST;
- icl(NCLUST)=i;
- end
- end
- fprintf('NUMBER OF CLUSTERS: %i n', NCLUST);
- disp('Performing assignation')
- %assignation
- % 将其他数据点归类 (assignation)
- for i=1:ND
- if (cl(ordrho(i))==1)
- cl(ordrho(i))=cl(nneigh(ordrho(i)));
- end
- end
- %halo
- for i=1:ND
- halo(i)=cl(i);
- end
- if (NCLUST>1)
- % 初始化数组 bord_rho 为 0,每个 cluster 定义一个 bord_rho 值
- for i=1:NCLUST
- bord_rho(i)=0.;
- end
- for i=1:ND-1
- for j=i+1:ND
- % 距离足够小但不属于同一个 cluster 的 i 和 j
- if ((cl(i)~=cl(j))&& (dist(i,j)<=dc))
- rho_aver=(rho(i)+rho(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
- for i=1:ND
- if (rho(i)
- halo(i)=0;
- end
- end
- end
- % 逐一处理每个 cluster
- for i=1:NCLUST
- nc=0;
- nh=0;
- for j=1:ND
- if (cl(j)==i)
- nc=nc+1;
- end
- if (halo(j)==i)
- nh=nh+1;
- end
- end
- fprintf('CLUSTER: %i CENTER: %i ELEMENTS: %i CORE: %i HALO: %i n', i,icl(i),nc,nh,nc-nh);
- end
- cmap=colormap;
- for i=1:NCLUST
- ic=int8((i*64.)/(NCLUST*1.));
- subplot(2,1,1)
- hold on
- plot(rho(icl(i)),delta(icl(i)),'o','MarkerSize',8,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
- end
- subplot(2,1,2)
- disp('Performing 2D nonclassical multidimensional scaling')
- Y1 = mdscale(dist, 2, 'criterion','metricstress');
- plot(Y1(:,1),Y1(:,2),'o','MarkerSize',2,'MarkerFaceColor','k','MarkerEdgeColor','k');
- title ('2D Nonclassical multidimensional scaling','FontSize',15.0)
- xlabel ('X')
- ylabel ('Y')
- for i=1:ND
- A(i,1)=0.;
- A(i,2)=0.;
- end
- for i=1:NCLUST
- nn=0;
- ic=int8((i*64.)/(NCLUST*1.));
- for j=1:ND
- if (halo(j)==i)
- nn=nn+1;
- A(nn,1)=Y1(j,1);
- A(nn,2)=Y1(j,2);
- end
- end
- hold on
- plot(A(1:nn,1),A(1:nn,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
- end
- %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');
- disp('Generated file:CLUSTER_ASSIGNATION')
- disp('column 1:element id')
- disp('column 2:cluster assignation without halo control')
- disp('column 3:cluster assignation with halo control')
- for i=1:ND
- fprintf(faa, '%i %i %in',i,cl(i),halo(i));
- end
出现错误如图：
0
|
|
3个回答
|
|
|
|
|
/7 
工商网监
湘ICP备2023018690号
2982
评论
淘帖
3454
