[问答] 想将density peaks密度峰值算法用于图像分割出现error

[复制链接]
发表于 2017-6-22 16:50:19   701 查看 1 回复 显示全部楼层 倒序浏览
分享
下面是将density peaks密度峰值算法用于图像处理的改进代码,出现错误,希望大神帮忙看一下怎么改正。
  1. clear all
  2. close all
  3. % disp('The only input needed is a distance matrix file')
  4. % disp('The format of this file should be: ')
  5. % disp('Column 1: id of element i')
  6. % disp('Column 2: id of element j')
  7. % disp('Column 3: dist(i,j)')
  8. % mdist=input('name of the distance matrix file (with single quotes)?\n');
  9. % disp('Reading input distance matrix')
  10. % xx=load(mdist);
  11. % xx = load('D:\example_distances.dat')
  12. %x = load('C:\UseTraceOfAllUsers.txt');
  13. % 从文件中读取数据  
  14. x = imread('naochuxue.jpg');
  15. x=double(x);
  16. minX = min(x);
  17. maxX = max(x);%取较大值
  18. ran = maxX - minX;
  19. nx(:,1) = (x(:,1) - minX(1,1)) / ran(1,1);
  20. nx(:,2) = (x(:,2) - minX(1,2)) / ran(1,2);
  21. dist = pdist2(nx, nx);
  22. N = size(dist,1);%%第一个维度的长度,相当于文件的行数(即距离的总个数)
  23. xx = zeros((N-1)*N/2, 3);%初始化为零
  24. idx = 1;
  25. % 这里不考虑对角线元素
  26. for i=1:N
  27.     for j=i+1:N
  28.         xx(idx, 1) = i;
  29.         xx(idx, 2) = j;
  30.         xx(idx, 3) = dist(i, j);
  31.         idx = idx + 1;
  32.     end
  33. end
  34. N = size(xx, 1);        
  35. ND=max(xx(:,2));
  36. NL=max(xx(:,1));
  37. if (NL>ND)
  38.   ND=NL;
  39. end
  40. % N=size(xx,1);
  41. % for i=1:ND
  42. %   for j=1:ND
  43. %     dist(i,j)=0;
  44. %   end
  45. % end
  46. % for i=1:N
  47. %   ii=xx(i,1);
  48. %   jj=xx(i,2);
  49. %   dist(ii,jj)=xx(i,3);
  50. %   dist(jj,ii)=xx(i,3);
  51. % end
  52. % percent=2;
  53. % fprintf('average percentage of neighbours (hard coded): %5.6f\n', percent);
  54. %
  55. % position=round(N*percent/100);
  56. % sda=sort(xx(:,3));
  57. % dc=sda(position);
  58. dc = 0.15;
  59. %计算局部密度 rho (利用 Gaussian 核)
  60. fprintf('Computing Rho with gaussian kernel of radius: %12.6f\n', dc);

  61. % 将每个数据点的 rho 值初始化为零  
  62. for i=1:ND
  63.   rho(i)=0.;
  64. end
  65. %
  66. % Gaussian kernel
  67. %
  68. for i=1:ND-1
  69.   for j=i+1:ND
  70.      rho(i)=rho(i)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));
  71.      rho(j)=rho(j)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));
  72.   end
  73. end
  74. %
  75. % "Cut off" kernel
  76. %
  77. % for i=1:ND-1
  78. %  for j=i+1:ND
  79. %    if (dist(i,j)<dc)S
  80. %       rho(i)=rho(i)+1.;
  81. %       rho(j)=rho(j)+1.;
  82. %    end
  83. %  end
  84. % end
  85. % 先求矩阵列最大值,再求最大值,最后得到所有距离值中的最大值
  86. maxd=max(max(dist));
  87. % 将 rho 按降序排列,ordrho 保持序  
  88. [rho_sorted,ordrho]=sort(rho,'descend');
  89. % 处理 rho 值最大的数据点
  90. delta(ordrho(1))=-1.;
  91. nneigh(ordrho(1))=0;
  92. % 生成 delta 和 nneigh 数组  
  93. for ii=2:ND
  94.    delta(ordrho(ii))=maxd;
  95.    for jj=1:ii-1
  96.      if(dist(ordrho(ii),ordrho(jj))<delta(ordrho(ii)))
  97.         delta(ordrho(ii))=dist(ordrho(ii),ordrho(jj));
  98.         nneigh(ordrho(ii))=ordrho(jj);
  99.         % 记录 rho 值更大的数据点中与 ordrho(ii) 距离最近的点的编号 ordrho(jj)
  100.      end
  101.    end
  102. end
  103. % 生成 rho 值最大数据点的 delta 值
  104. delta(ordrho(1))=max(delta(:));
  105. % 决策图  
  106. disp('Generated file:DECISION GRAPH')
  107. disp('column 1:Density')
  108. disp('column 2:Delta')

  109. fid = fopen('DECISION_GRAPH', 'w');
  110. for i=1:ND
  111.    fprintf(fid, '%6.2f %6.2f\n', rho(i),delta(i));
  112. end
  113. % 选择一个围住类中心的矩形
  114. disp('Select a rectangle enclosing cluster centers')
  115. scrsz = get(0,'ScreenSize');
  116. figure('Position',[6 72 scrsz(3)/4. scrsz(4)/1.3]);
  117. for i=1:ND
  118.   ind(i)=i;
  119.   gamma(i)=rho(i)*delta(i);
  120. end
  121. subplot(2,1,1)
  122. tt=plot(rho(:),delta(:),'o','MarkerSize',5,'MarkerFaceColor','k','MarkerEdgeColor','k');
  123. title ('Decision Graph','FontSize',15.0)
  124. xlabel ('\rho')
  125. ylabel ('\delta')

  126. % 利用 rho 和 delta 画出一个决策图
  127. subplot(2,1,1)
  128. rect = getrect(1);
  129. rhomin=rect(1);
  130. deltamin=rect(4);
  131. % 初始化 cluster 个数  
  132. NCLUST=0;
  133. % cl 为归属标志数组,cl(i)=j 表示第 i 号数据点归属于第 j 个 cluster  
  134. % 将 cl 初始化为 -1
  135. for i=1:ND
  136.   cl(i)=-1;
  137. end
  138. % 在矩形区域内统计数据点(即聚类中心)的个数
  139. for i=1:ND
  140.   if ( (rho(i)>rhomin) && (delta(i)>deltamin))
  141.      NCLUST=NCLUST+1;
  142.      cl(i)=NCLUST;
  143.      icl(NCLUST)=i;
  144.   end
  145. end
  146. fprintf('NUMBER OF CLUSTERS: %i \n', NCLUST);
  147. disp('Performing assignation')

  148. %assignation
  149. % 将其他数据点归类 (assignation)
  150. for i=1:ND
  151.   if (cl(ordrho(i))==1)
  152.     cl(ordrho(i))=cl(nneigh(ordrho(i)));
  153.   end
  154. end
  155. %halo
  156. for i=1:ND
  157.   halo(i)=cl(i);
  158. end
  159. if (NCLUST>1)
  160.     % 初始化数组 bord_rho 为 0,每个 cluster 定义一个 bord_rho 值  
  161.   for i=1:NCLUST
  162.     bord_rho(i)=0.;
  163.   end
  164.   for i=1:ND-1
  165.     for j=i+1:ND
  166.         % 距离足够小但不属于同一个 cluster 的 i 和 j  
  167.       if ((cl(i)~=cl(j))&& (dist(i,j)<=dc))
  168.         rho_aver=(rho(i)+rho(j))/2.;% 取 i,j 两点的平均局部密度  
  169.         if (rho_aver>bord_rho(cl(i)))
  170.           bord_rho(cl(i))=rho_aver;
  171.         end
  172.         if (rho_aver>bord_rho(cl(j)))
  173.           bord_rho(cl(j))=rho_aver;
  174.         end
  175.       end
  176.     end
  177.   end
  178.   for i=1:ND
  179.     if (rho(i)<bord_rho(cl(i)))
  180.       halo(i)=0;
  181.     end
  182.   end
  183. end
  184. % 逐一处理每个 cluster
  185. for i=1:NCLUST
  186.   nc=0;
  187.   nh=0;
  188.   for j=1:ND
  189.     if (cl(j)==i)
  190.       nc=nc+1;
  191.     end
  192.     if (halo(j)==i)
  193.       nh=nh+1;
  194.     end
  195.   end
  196.   fprintf('CLUSTER: %i CENTER: %i ELEMENTS: %i CORE: %i HALO: %i \n', i,icl(i),nc,nh,nc-nh);
  197. end

  198. cmap=colormap;
  199. for i=1:NCLUST
  200.    ic=int8((i*64.)/(NCLUST*1.));
  201.    subplot(2,1,1)
  202.    hold on
  203.    plot(rho(icl(i)),delta(icl(i)),'o','MarkerSize',8,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
  204. end
  205. subplot(2,1,2)
  206. disp('Performing 2D nonclassical multidimensional scaling')
  207. Y1 = mdscale(dist, 2, 'criterion','metricstress');
  208. plot(Y1(:,1),Y1(:,2),'o','MarkerSize',2,'MarkerFaceColor','k','MarkerEdgeColor','k');
  209. title ('2D Nonclassical multidimensional scaling','FontSize',15.0)
  210. xlabel ('X')
  211. ylabel ('Y')
  212. for i=1:ND
  213. A(i,1)=0.;
  214. A(i,2)=0.;
  215. end
  216. for i=1:NCLUST
  217.   nn=0;
  218.   ic=int8((i*64.)/(NCLUST*1.));
  219.   for j=1:ND
  220.     if (halo(j)==i)
  221.       nn=nn+1;
  222.       A(nn,1)=Y1(j,1);
  223.       A(nn,2)=Y1(j,2);
  224.     end
  225.   end
  226.   hold on
  227.   plot(A(1:nn,1),A(1:nn,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
  228. end

  229. %for i=1:ND
  230. %   if (halo(i)>0)
  231. %      ic=int8((halo(i)*64.)/(NCLUST*1.));
  232. %      hold on
  233. %      plot(Y1(i,1),Y1(i,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),'MarkerEdgeColor',cmap(ic,:));
  234. %   end
  235. %end
  236. faa = fopen('CLUSTER_ASSIGNATION', 'w');
  237. disp('Generated file:CLUSTER_ASSIGNATION')
  238. disp('column 1:element id')
  239. disp('column 2:cluster assignation without halo control')
  240. disp('column 3:cluster assignation with halo control')
  241. for i=1:ND
  242.    fprintf(faa, '%i %i %i\n',i,cl(i),halo(i));
  243. end
复制代码
出现错误如图:

error

error

标签:matlab
z00

VIP会员

发表于 5 天前  

回帖奖励 +1 分积分

没用过这个 帮顶一下
回复

点赞 举报

高级模式
您需要登录后才可以回帖 登录 | 注册

关闭

站长推荐 上一条 /9 下一条

快速回复 返回顶部 返回列表
-

推荐专区

技术干货集中营

专家问答

方案交易

用户帮助┃咨询与建议┃版主议事

工程师杂谈

项目|工程师创意

招聘|求职}工程师职场

论坛电子赛事

社区活动专版

发烧友活动

-

嵌入式论坛

ARM技术论坛

Android论坛

Linux论坛

单片机/MCU论坛

MSP430技术论坛

FPGA|CPLD|ASIC论坛

STM32/STM8技术论坛

NXP MCU 技术论坛

PIC单片机论坛

DSP论坛

瑞萨单片机论坛

嵌入式系统论坛

-

电源技术论坛

电源技术论坛

无线充电技术

-

硬件设计论坛

PCB设计论坛

电路设计论坛

电子元器件论坛

控制|传感

总线技术|接口技术

-

测试测量论坛

LabVIEW论坛

Matlab论坛

测试测量技术专区

仪器仪表技术专区

-

EDA设计论坛

multisim论坛

PADS技术论坛

Protel|AD|DXP论坛

Allegro论坛

proteus论坛|仿真论坛

EasyEDA-中国人自已的EDA工具

Orcad论坛

-

综合技术与应用

电机控制

智能电网

光电及显示

工程资源中心

汽车电子技术论坛

医疗电子论坛

-

开源硬件

-

无线通信论坛

无线通信技术专区

天线|RF射频|微波|雷达技术

-

IC设计论坛

芯片测试与失效分析

Mixed Signal/SOC[数模混合芯片设计]

Analog/RF IC设计

设计与制造封装测试

-

厂商专区

TI论坛

TI Deyisupport社区

-

检测技术与质量

电磁兼容(EMC)设计与整改

安规知识论坛

检测与认证

-

消费电子论坛

手机技术论坛

平板电脑/mid论坛

音视/视频/机顶盒论坛

-

电子论坛综合区

聚丰众筹官方社区

新人报道区

聚丰供应链

-

论坛服务区

-

供求信息发布

供需广告

电子展览展会专区

芯片求购|供应发布区