泰森多边形（Voronoi彩图）的matlab绘制——2

大家好，欢迎来到IT知识分享网。

泰森多边形（Voronoi图）的matlab绘制——彩图版

1 Voronoi图简介

泰森多边形是对空间平面的一种剖分，其特点是多边形内的任何位置离该多边形的样点（如居民点）的距离最近，离相邻多边形内样点的距离远，且每个多边形内含且仅包含一个样点。由于泰森多边形在空间剖分上的等分性特征，因此可用于解决最近点、最小封闭圆等问题，以及许多空间分析问题，如邻接、接近度和可达性分析等。

泰森多边形的构建可以分为2个步骤，1是Delaunay三角网的构建，2是三角网格外接圆心得连线。

之前我有写过一篇文章，主要是泰森多边形的边缘绘制。内边的思路是根据每一个三角形的Delaunay三角网关系，去依次连接每个相邻三角形外接圆心得到。外边的思路是根据每个外接圆圆心和相应边缘做中垂线得到的。
参见：
泰森多边形（Voronoi图）的matlab绘制 https://blog.csdn.net/weixin_/article/details/

二维Delaunay三角网的绘制参见 https://blog.csdn.net/weixin_/article/details/

2 实心带色彩的Voronoi图实现思路

其中本文最主要的用法是

patch('Faces',f,'Vertices',v,'FaceVertexCData',col,'FaceColor','flat'); 

v是一系列点矩阵，f是多边形对应的点的编号，col是多边形的颜色，取[0,1]。

根据前两篇文章里的算法，已知所有的初始点xdot，和已经构建关系的三角形网格。之后循环每一个初始点，依次连接点初始点周围三角形的外接圆心即可。

比如对于以下图形：

以点4为例，和点4相关的三角形依次为

[4,3,1; 4,1,2; 4,2,5; 4,5,7; 4,7,6; 4,6,3;] 

这些三角形，他们的外接圆圆心为红点所示。所以点4所对应的泰森多边形为这些红色点依次相连的对应的多边形：

如果遇到边缘三角形，思路和上一篇文章相同，都是做边缘的延长线。和上一篇文章中惟一的不同就是，上一篇文章只需要得到线即可，所以延长线画到很远处就行；但是这里要求得到面，所以需要求延长线与边界的交点。这里做的是边缘延长线与边界的交点作为泰森多边形的一个顶点。

3 matlab代码以及得到的结果图

3.1 matlab自带的代码

%采用matlab自带的函数进行绘制 clear xdot=gallery('uniformdata',[200 2],5); %delaunay三角形 figure(1) DT=delaunayTriangulation(xdot); triplot(DT,'color','k') %voronoi三角形 figure(2) voronoi(xdot(:,1),xdot(:,2)); xlim([0,1]) ylim([0,1]) %绘制voronoi三角形彩色 figure(3) [v,c]=voronoin(xdot); %不是voronoi是voronoin xlim([0,1]) ylim([0,1]) for i = 1:length(c) patch(v(c{ 
   i},1),v(c{ 
   i},2),rand(1,1)); end 

3.2 自编程实现代码

因为涉及到了之前的两篇文章的内容，所以代码看起来非常的冗长。其实这篇代码完全可以采用外部函数调用的方式做到很整洁，只需要多创立几个m文件把这段代码拆分即可。

%更改了3里面几个bug， %1.结果要输出最终面的形式，添加Txdot clear %0 设置初始点模块* N=400; %点随机 xdot=rand(N,2); %点按圆形随机 % r=rand(N,1).^0.3; % theta=rand(N,1)*2*pi; % xdot=[r.*cos(theta)/2+0.5,r.*sin(theta)/2+0.5]; %点按双行随机 % x=rand(N,1); % y=[randn(N/2,1)/5+0.5;randn(N/2,1)/5-0.5]; % y(y>1)=1;y(y<-1)=-1; % y=(y+1)/2.1; % xdot=[x,y]; %点按规则矩形加抖动 % [X1,X2]=meshgrid(0:1/sqrt(N):1-1/sqrt(N)); % xdot=zeros(N,2); % xdot(:,1)=X1(1:end)'+1/sqrt(N)/2*rand(N,1); % xdot(:,2)=X2(1:end)'+1/sqrt(N)/2*rand(N,1); %点按随机三角加抖动 % NN=20; % X1=[];X2=[]; % for j=1:NN % if mod(j,2)==0 % X1=[X1;(0:1/NN/sqrt(3)*2:1-0/NN/sqrt(3)*2)']; % X2=[X2;ones(length(0:1/NN/sqrt(3)*2:1-0/NN/sqrt(3)*2),1)*(j-1)/NN]; % else % X1=[X1;(0:1/NN/sqrt(3)*2:1-1/NN/sqrt(3)*2)'+1/NN/sqrt(3)]; % X2=[X2;ones(length(0:1/NN/sqrt(3)*2:1-1/NN/sqrt(3)*2),1)*(j-1)/NN]; % end % end % N=size(X1,1); % xdot=[X1+rand(N,1)*1.2/NN/sqrt(3),X2+rand(N,1)*1.2/NN/2]; %点按圆形螺旋 r=0:1/500:1; r=r.^0.7; theta=0:0.1:500*0.1; theta=theta.^1.3; xdot=[(r.*cos(theta))'/2+0.5,(r.*sin(theta))'/2+0.5]; N=size(xdot,1); %1 Delaulay三角形的构建* %整理点，遵循从左到右，从上到下的顺序 xdot=sortrows(xdot,[1 2]); %画出最大包含的三角形 xmin=min(xdot(:,1));xmax=max(xdot(:,1)); ymin=min(xdot(:,2));ymax=max(xdot(:,2)); bigtri=[(xmin+xmax)/2-(xmax-xmin)*1.5,ymin-(xmax-xmin)*0.5;... (xmin+xmax)/2,ymax+(ymax-ymin)+(xmax-xmin)*0.5;... (xmin+xmax)/2+(xmax-xmin)*1.5,ymin-(xmax-xmin)*0.5]; xdot=[bigtri;xdot];%点集 edgemat=[1 2 xdot(1,:) xdot(2,:);... 2 3 xdot(2,:) xdot(3,:);1 3 xdot(1,:) xdot(3,:)];%边集，每个点包含2个点，4个坐标值 trimat=[1 2 3];%三角集,每个三角包含3个点 temp_trimat=[1 2 3]; for j=4:N+3 pointtemp=xdot(j,:);%循环每一个点 deltemp=[];%初始化删除temp_trimat的点 temp_edgemat=[];%初始化临时边 for k=1:size(temp_trimat,1)%循环每一个temp_trimat的三角形 panduan=whereispoint(xdot(temp_trimat(k,1),:),... xdot(temp_trimat(k,2),:),xdot(temp_trimat(k,3),:),pointtemp);%判断点在圆内0、圆外1、圆右侧2 switch panduan case 0 %点在圆内 %则该三角形不为Delaunay三角形 temp_edge=maketempedge(temp_trimat(k,1),temp_trimat(k,2),temp_trimat(k,3),j,xdot);%把三条边暂时存放于临时边矩阵 temp_edgemat=[temp_edgemat;temp_edge]; deltemp=[deltemp,k]; ; case 1 %点在圆外，pass ; case 2 %点在圆右 %则该三角形为Delaunay三角形，保存到triangles trimat=[trimat;temp_trimat(k,:)];%添加到正式三角形中 deltemp=[deltemp,k]; %并在temp里去除掉 %别忘了把正式的边也添加进去 edgemat=[edgemat;makeedge(temp_trimat(k,1),temp_trimat(k,2),temp_trimat(k,3),xdot)];%遵循12,13,23的顺序 edgemat=unique(edgemat,'stable','rows'); end % clf % hold on % plot(xdot(:,1),xdot(:,2),'ko') % plot(pointtemp(:,1),pointtemp(:,2),'bo') % for h=1:size(edgemat,1) % plot([edgemat(h,3),edgemat(h,5)],[edgemat(h,4),edgemat(h,6)],'k-') % end % plot([xdot(temp_trimat(k,1),1),xdot(temp_trimat(k,2),1),xdot(temp_trimat(k,3),1),xdot(temp_trimat(k,1),1)],... % [xdot(temp_trimat(k,1),2),xdot(temp_trimat(k,2),2),xdot(temp_trimat(k,3),2),xdot(temp_trimat(k,1),2)],'r-') %  %  % hold off %三角循环结束  end %除去上述步骤中的临时三角形 temp_trimat(deltemp,:)=[]; temp_trimat(~all(temp_trimat,2),:)=[]; %对temp_edgemat去重复 temp_edgemat=unique(temp_edgemat,'stable','rows'); %将edge buffer中的边与当前的点进行组合成若干三角形并保存至temp triangles中 temp_trimat=[temp_trimat;maketemptri(temp_edgemat,xdot,j)]; k=k; %点循环结束 end %合并temptri trimat=[trimat;temp_trimat]; edgemat=[edgemat;temp_edgemat]; %删除大三角形 deltemp=[]; for j=1:size(trimat,1) if ismember(1,trimat(j,:))||ismember(2,trimat(j,:))||ismember(3,trimat(j,:)) deltemp=[deltemp,j]; end end trimat(deltemp,:)=[]; edgemat=[trimat(:,[1,2]);trimat(:,[2,3]);trimat(:,[3,1])]; edgemat=sort(edgemat,2); edgemat=unique(edgemat,'stable','rows'); temp_edgemat=[]; temp_trimat=[]; % figure(1) % hold on % % plot(xdot(:,1),xdot(:,2),'ko') % for j=1:size(trimat,1) % plot([xdot(trimat(j,1),1),xdot(trimat(j,2),1)],[xdot(trimat(j,1),2),xdot(trimat(j,2),2)],'k-') % plot([xdot(trimat(j,1),1),xdot(trimat(j,3),1)],[xdot(trimat(j,1),2),xdot(trimat(j,3),2)],'k-') % plot([xdot(trimat(j,3),1),xdot(trimat(j,2),1)],[xdot(trimat(j,3),2),xdot(trimat(j,2),2)],'k-') % end % hold off % xlim([0,1]);ylim([0,1]); %1.5 凸包监测 %思路是先找出边缘点（三角形只有1个或2个的），顺便整出一个三角形相互关系图，以后用。 %然后顺时针，依次隔一个点连接出一条线段，如果这个和之前的线段相交，则不算；如果不交，则记录出三角形 %更新完了以后，再监测一遍，直到没有新的为止。 t_w=0; while t_w==0 [~,border_point,~]=makebordertri(trimat); border_point=[border_point;border_point(1,:)]; temp_edgemat=[]; temp_trimat=[]; for j=1:size(border_point,1)-1 tempboderedge=[border_point(j,1),border_point(j+1,2)]; tempboderdot=border_point(j,2); %寻找带tempboderdot的所有边 tempdotex=edgemat(logical(sum(edgemat==tempboderdot,2)),:); %删除相邻边 tempdotex(ismember(tempdotex,[tempboderdot,tempboderedge(1)],'rows'),:)=[]; tempdotex(ismember(tempdotex,[tempboderedge(1),tempboderdot],'rows'),:)=[]; tempdotex(ismember(tempdotex,[tempboderdot,tempboderedge(2)],'rows'),:)=[]; tempdotex(ismember(tempdotex,[tempboderedge(2),tempboderdot],'rows'),:)=[]; %检测tempdotex是否为空，如果是证明不用相连 t_N=size(tempdotex,1); t_t=0; if t_N>0 %依次检测是否相交，只要有一个相交就不算；如果都不想交，则相连 for k=1:t_N if tempdotex(k,1)==tempboderdot t_xdotno4=tempdotex(k,2); else t_xdotno4=tempdotex(k,1); end tt_xdotno4=xdot(t_xdotno4,:)-xdot(tempboderdot,:); xdotno4=xdot(tempboderdot,:)+tt_xdotno4/sqrt(sum(tt_xdotno4.^2))*(sqrt((xmax-xmin)^2+(ymax-ymin)^2)); panduan=crossornot(xdot(tempboderedge(1),:),xdot(tempboderedge(2),:),xdot(tempboderdot,:),xdotno4); if panduan==1 t_t=t_t+1; break end end %t_t大于0说明有相交的线,略过 if t_t==0 temp_edgemat=[temp_edgemat;tempboderedge]; temp_trimat=[temp_trimat;[tempboderedge,tempboderdot]]; break end end end trimat=[trimat;temp_trimat]; edgemat=[edgemat;temp_edgemat]; %删除重复的三角形 trimat=sort(trimat,2); trimat=unique(trimat,'stable','rows'); if j==size(border_point,1)-1 t_w=1; end end figure(2) hold on % plot(xdot(:,1),xdot(:,2),'ko') for j=1:size(trimat,1) plot([xdot(trimat(j,1),1),xdot(trimat(j,2),1)],[xdot(trimat(j,1),2),xdot(trimat(j,2),2)],'k-') plot([xdot(trimat(j,1),1),xdot(trimat(j,3),1)],[xdot(trimat(j,1),2),xdot(trimat(j,3),2)],'k-') plot([xdot(trimat(j,3),1),xdot(trimat(j,2),1)],[xdot(trimat(j,3),2),xdot(trimat(j,2),2)],'k-') end hold off xlim([0,1]);ylim([0,1]); %2 泰森多边形的建立步骤  %求每个三角形的外接圆圆心 trimatcenter=zeros(size(trimat,1),2); for j=1:size(trimat,1) [a,b,~]=maketricenter(xdot(trimat(j,1),:),xdot(trimat(j,2),:),xdot(trimat(j,3),:)); trimatcenter(j,:)=[a,b]; end hold on scatter(trimatcenter(:,1),trimatcenter(:,2),5,[0.6,0,0],'filled') hold off %求三角形的相邻三角形个数 [border_trimat,border_point,trimat_con]=makebordertri(trimat); Thi_edge1=[]; for j=1:size(trimat,1) tempedge=[]; %第一个相邻三角形 if trimat_con(j,1)~=0 tempedge=[tempedge;[j,trimat_con(j,1)]]; end %第二个相邻三角形 if trimat_con(j,2)~=0 tempedge=[tempedge;[j,trimat_con(j,2)]]; end %第三个相邻三角形 if trimat_con(j,3)~=0 tempedge=[tempedge;[j,trimat_con(j,3)]]; end Thi_edge1=[Thi_edge1;tempedge]; end %绘制非边缘泰勒多边形 figure(3) Thi_edge1=unique(Thi_edge1,'stable','rows'); xlim([0,1]);ylim([0,1]); hold on for j=1:size(Thi_edge1,1) plot(trimatcenter([Thi_edge1(j,1),Thi_edge1(j,2)],1),trimatcenter([Thi_edge1(j,1),Thi_edge1(j,2)],2),'color',[0,0.4,0]) end %绘制边缘泰勒多边形 %先逐个边试探，如果中心点在三角内，则做中心-边缘延长线 %如果中心点在三角外，如果在屏幕外，忽略，如果在屏幕内，做边缘-中心延长线 for j=1:size(border_point,1) %先找到边对应的三角 temp_trimat=border_trimat(sum(border_trimat==border_point(j,1),2)+sum(border_trimat==border_point(j,2),2)==2,:); %判断中心点是否在三角形内 [t_x1,t_y1,~]=maketricenter(xdot(temp_trimat(1),:),xdot(temp_trimat(2),:),xdot(temp_trimat(3),:));%求中心 panduan=pointintriangle(xdot(temp_trimat(1),:),xdot(temp_trimat(2),:),xdot(temp_trimat(3),:),[t_x1,t_y1]); %求边的中点 t_x2=(xdot(border_point(j,1),1)+xdot(border_point(j,2),1))/2; t_y2=(xdot(border_point(j,1),2)+xdot(border_point(j,2),2))/2; if panduan==1 %做中心-边缘的延长线 %这里用到了边缘在01这个条件 t_xy3=[t_x1,t_y1]+[t_x2-t_x1,t_y2-t_y1]*sqrt(2)/sqrt((t_x2-t_x1)^2+(t_y2-t_y1)^2); plot([t_x1,t_xy3(1)],[t_y1,t_xy3(2)],'color',[0,0.4,0]) elseif ~(t_x1<0||t_x1>1||t_y1<0||t_y1>1) %判断点是否在边与边框的三角内，如果在，做中心的延长线 %如果不在，做中心-边缘的延长线 %或者改成判断点是否在多边形内 t_t=pointinmutiangle(xdot,[border_point(1,1);border_point(:,2)],[t_x1,t_y1]); if t_t==1 t_xy3=[t_x1,t_y1]+[t_x2-t_x1,t_y2-t_y1]*sqrt(2)/sqrt((t_x2-t_x1)^2+(t_y2-t_y1)^2); plot([t_x1,t_xy3(1)],[t_y1,t_xy3(2)],'color',[0,0.4,0]) else t_xy3=[t_x1,t_y1]+[t_x1-t_x2,t_y1-t_y2]*1/sqrt((t_x2-t_x1)^2+(t_y2-t_y1)^2); plot([t_x1,t_xy3(1)],[t_y1,t_xy3(2)],'color',[0,0.4,0]) end end end scatter(xdot(:,1),xdot(:,2),5,[0,0.4,0],'filled') hold off %第3部分 多边形面的识别 %1先规划出多边形预设矩阵（找出xdot数组中出现次数最多的数），行个数等于三角形外接圆中心点的个数 %2选取中心点，然后依次找出一圈三角形（共点），之后依次连接这一圈三角形的外接圆就是泰森多边形 %3对于边缘点，要结合边缘考虑，画出边缘线和延长射线的交点，就是了。还是依次画出，如果点在边缘外，就忽略，在边缘内，就记录。 %4选取上一步记录的三角形最外侧三角，做中垂线 %初始化面 % load('color_Sky02_h.mat') figure(4) % colormap(mycolor) xlim([0,1]);ylim([0,1]); %依次根据xdot来循环 for j=4:size(xdot,1) %trimatcenter [temp_trimat,temp_trimat_no,panduan]=findpoint2tri(j,trimat);%找出这一圈三角形 if panduan %直接把temp_trimat_no的点当做泰森多边形即可 patch('Faces',temp_trimat_no,'Vertices',trimatcenter,'FaceVertexCData',rand(1,1),'FaceColor','flat'); 'hyh'; else %需要做边缘化处理 %如果最终三角形不是边缘三角形，要注意不能适用内外部延长的准则，要用连接法则 %选出中心点在边界内的一圈三角形 [temp_trimat,temp_trimat_no,delleft,delright]=selecttemp_trimat(temp_trimat,temp_trimat_no,trimatcenter); t_t=1; % 如果selecttemp_trimat删除了边上的点，那么做边缘延长线的话就用和删除点之间作为连线，和边相交即可 %即如果delleft,delright=0，照旧。如果不是0，利用temp_trimat_no(1)和delleft做连线（创建一个公式） tempboderedge=trimatcenter(temp_trimat_no,:); %边缘延长线1 if delleft==0 tempedge=temp_trimat(1,[1,2]); tempboderdot1=edgepointfind(tempedge,xdot,trimatcenter(temp_trimat_no(1),:),border_point);%获取边缘点 if size(tempboderdot1,1)==0 t_t=0;%边缘点超出，不再画面 end else tempboderdot1=edgepointfind2(trimatcenter(temp_trimat_no(1),:),trimatcenter(delleft,:)); end tempboderedge=[tempboderdot1;tempboderedge]; %边缘延长线2 if delright==0 tempedge=temp_trimat(end,[1,3]); tempboderdot2=edgepointfind(tempedge,xdot,trimatcenter(temp_trimat_no(end),:),border_point);%获取边缘点 if size(tempboderdot2,1)==0 t_t=0;%边缘点超出，不再画面 end else tempboderdot2=edgepointfind2(trimatcenter(temp_trimat_no(end),:),trimatcenter(delright,:)); end tempboderedge=[tempboderedge;tempboderdot2]; %绘制边缘图形 if t_t==1 if tempboderdot1(1)~=tempboderdot2(1) if tempboderdot1(2)~=tempboderdot2(2) %边缘延长线3(边界角点 tempboderedge=[tempboderedge;maketempboderdot(tempboderdot1,tempboderdot2)]; end end patch('Faces',1:length(tempboderedge),'Vertices',tempboderedge,'FaceVertexCData',rand(1,1),'FaceColor','flat'); 'hyh2'; end end end %判断点在三角形外接圆的哪个部分 function panduan=whereispoint(xy1,xy2,xy3,xy0) %判断点在三角形外接圆的哪个部分 [a,b,r2]=maketricenter(xy1,xy2,xy3); x0=xy0(1);y0=xy0(2); if a+sqrt(r2)<x0 %x0在圆的右侧 panduan=2; elseif (x0-a)^2+(y0-b)^2<r2 %x0在圆内 panduan=0; else %在圆外 panduan=1; end end %做出三角形三点与内部1点之间的线段 function temp_edge=maketempedge(dot1,dot2,dot3,dot0,xdot) %做出连接点与三角形之间的线 %每行包含2个点，4个坐标值，共3行 %xy1和xy0组成线段 temp_edge=zeros(3,6); if xdot(dot1,1)<xdot(dot0,1) temp_edge(1,:)=[dot1,dot0,xdot(dot1,:),xdot(dot0,:)]; elseif xdot(dot1,1)==xdot(dot0,1) if xdot(dot1,2)<xdot(dot0,2) temp_edge(1,:)=[dot1,dot0,xdot(dot1,:),xdot(dot0,:)]; else temp_edge(1,:)=[dot0,dot1,xdot(dot0,:),xdot(dot1,:)]; end else temp_edge(1,:)=[dot0,dot1,xdot(dot0,:),xdot(dot1,:)]; end %xy2和xy0组成线段 if xdot(dot2,1)<xdot(dot0,1) temp_edge(2,:)=[dot2,dot0,xdot(dot2,:),xdot(dot0,:)]; elseif xdot(dot2,1)==xdot(dot0,1) if xdot(dot2,2)<xdot(dot0,2) temp_edge(2,:)=[dot2,dot0,xdot(dot2,:),xdot(dot0,:)]; else temp_edge(2,:)=[dot0,dot2,xdot(dot0,:),xdot(dot2,:)]; end else temp_edge(2,:)=[dot0,dot2,xdot(dot0,:),xdot(dot2,:)]; end %xy3和xy0组成线段 if xdot(dot3,1)<xdot(dot0,1) temp_edge(3,:)=[dot3,dot0,xdot(dot3,:),xdot(dot0,:)]; elseif xdot(dot3,1)==xdot(dot0,1) if xdot(dot3,2)<xdot(dot0,2) temp_edge(3,:)=[dot3,dot0,xdot(dot3,:),xdot(dot0,:)]; else temp_edge(3,:)=[dot0,dot3,xdot(dot0,:),xdot(dot3,:)]; end else temp_edge(3,:)=[dot0,dot3,xdot(dot0,:),xdot(dot3,:)]; end end %做出一些列固定点发散的线段外点组成的三角形 function temp_trimat=maketemptri(temp_edgemat,xdot,dot0) %将edge buffer中的边与当前的点进行组合成若干三角形 %temp_edgemat是新边，x是中心点 %思路是计算各个边对应角度，然后排序相连 A=temp_edgemat(:,1:2); pointline=A(A~=dot0); N=length(pointline); pointaxe=xdot(pointline,:); img_pointaxe=pointaxe(:,1)+1i*pointaxe(:,2); d_img_pointaxe=img_pointaxe-xdot(dot0,1)-1i*xdot(dot0,2); angle_d_img_pointaxe=angle(d_img_pointaxe); [~,index]=sort(angle_d_img_pointaxe); index=[index;index(1)];%排序，然后依次串起来 temp_trimat=zeros(N,3); for j=1:N temp_trimat(j,:)=[pointline(index(j)),pointline(index(j+1)),dot0]; end end %将三个点构成3条边 function edgemat=makeedge(dot1,dot2,dot3,xdot) %将dot1 2 3这三个点构成三条边 %每行包含2个点，4个坐标值，共3行 edgemat=zeros(3,6); %点12 if xdot(dot1,1)<xdot(dot2,1) edgemat(1,:)=[dot1,dot2,xdot(dot1,:),xdot(dot2,:)]; elseif xdot(dot1,1)==xdot(dot2,1) if xdot(dot1,2)<xdot(dot2,2) edgemat(1,:)=[dot1,dot2,xdot(dot1,:),xdot(dot2,:)]; else edgemat(1,:)=[dot2,dot1,xdot(dot2,:),xdot(dot1,:)]; end else edgemat(1,:)=[dot2,dot1,xdot(dot2,:),xdot(dot1,:)]; end %点13 if xdot(dot1,1)<xdot(dot3,1) edgemat(2,:)=[dot1,dot3,xdot(dot1,:),xdot(dot3,:)]; elseif xdot(dot1,1)==xdot(dot3,1) if xdot(dot1,2)<xdot(dot3,2) edgemat(2,:)=[dot1,dot3,xdot(dot1,:),xdot(dot3,:)]; else edgemat(2,:)=[dot3,dot1,xdot(dot3,:),xdot(dot1,:)]; end else edgemat(2,:)=[dot3,dot1,xdot(dot3,:),xdot(dot1,:)]; end %点23 if xdot(dot3,1)<xdot(dot2,1) edgemat(3,:)=[dot3,dot2,xdot(dot3,:),xdot(dot2,:)]; elseif xdot(dot3,1)==xdot(dot2,1) if xdot(dot3,2)<xdot(dot2,2) edgemat(3,:)=[dot3,dot2,xdot(dot3,:),xdot(dot2,:)]; else edgemat(3,:)=[dot2,dot3,xdot(dot2,:),xdot(dot3,:)]; end else edgemat(3,:)=[dot2,dot3,xdot(dot2,:),xdot(dot3,:)]; end % edgemat end %求三角形外接圆圆心 function [a,b,r2]=maketricenter(xy1,xy2,xy3) x1=xy1(1);y1=xy1(2); x2=xy2(1);y2=xy2(2); x3=xy3(1);y3=xy3(2); a=((y2-y1)*(y3*y3-y1*y1+x3*x3-x1*x1)-(y3-y1)*(y2*y2-y1*y1+x2*x2-x1*x1))/(2.0*((x3-x1)*(y2-y1)-(x2-x1)*(y3-y1))); b=((x2-x1)*(x3*x3-x1*x1+y3*y3-y1*y1)-(x3-x1)*(x2*x2-x1*x1+y2*y2-y1*y1))/(2.0*((y3-y1)*(x2-x1)-(y2-y1)*(x3-x1))); r2=(x1-a)*(x1-a)+(y1-b)*(y1-b); end %求边缘三角形 function [border_trimat,border_point,trimat_con]=makebordertri(trimat) N=size(trimat,1); border_trimat=[]; border_point=[]; trimat_con=zeros(N,3); for j=1:N %tempborder_trimat=zeros(3,3); temptri=trimat(j,:); %计算temptri中12点边对应的三角形有哪些 edgetrimat=find(sum(trimat==temptri(1),2)+sum(trimat==temptri(2),2)==2); edgetrimat(edgetrimat==j)=[]; if size(edgetrimat,2)==0 %这个边没有三角形相连，是个临边。 border_point=[border_point;[temptri(1),temptri(2)]]; elseif size(edgetrimat,2)==1 %这个边没有三角形相连，是个临边。 %tempborder_trimat(1,:)=trimat(edgetrimat,:);%记录三角形三点坐标 trimat_con(j,1)=edgetrimat;%trimat_con记录上相邻三角形 end %计算temptri中23点边对应的三角形有哪些 edgetrimat=find(sum(trimat==temptri(2),2)+sum(trimat==temptri(3),2)==2); edgetrimat(edgetrimat==j)=[]; if size(edgetrimat,2)==0 border_point=[border_point;[temptri(2),temptri(3)]]; elseif size(edgetrimat,2)==1 %tempborder_trimat(2,:)=trimat(edgetrimat,:); trimat_con(j,2)=edgetrimat; end %计算temptri中31点边对应的三角形有哪些 edgetrimat=find(sum(trimat==temptri(3),2)+sum(trimat==temptri(1),2)==2); edgetrimat(edgetrimat==j)=[]; if size(edgetrimat,2)==0 border_point=[border_point;[temptri(3),temptri(1)]]; elseif size(edgetrimat,2)==1 %tempborder_trimat(3,:)=trimat(edgetrimat,:); trimat_con(j,3)=edgetrimat; end %tempborder_trimat(all(tempborder_trimat==0, 2),:)=[];%删除0行 if ~all(trimat_con(j,:)) %如果边缘三角少于3个,就添加 border_trimat=[border_trimat;temptri]; end end %把边首尾排序一遍,输出border_point for j=1:size(border_point,1)-1 border_pointtemp=find(sum(border_point==border_point(j,2),2)==1); border_pointtemp(border_pointtemp==j)=[];%删除自己 border_point([j+1,border_pointtemp],:)=border_point([border_pointtemp,j+1],:); if border_point(j,2)==border_point(j+1,2) border_point(j+1,[1,2])=border_point(j+1,[2,1]); end end end %判断两个线段是否相交 function panduan=crossornot(l1xy1,l1xy2,l2xy1,l2xy2) l1x1=l1xy1(1);l1y1=l1xy1(2); l1x2=l1xy2(1);l1y2=l1xy2(2); l2x1=l2xy1(1);l2y1=l2xy1(2); l2x2=l2xy2(1);l2y2=l2xy2(2); %先快速判断 if (max(l2x1,l2x2)<min(l1x1,l1x2))||(max(l2y1,l2y2)<min(l1y1,l1y2))||... (max(l1x1,l1x2)<min(l2x1,l2x2))||(max(l1y1,l1y2)<min(l2y1,l2y2)) %如果判断为真，则一定不会相交 panduan=0; else %如果判断为假，进一步差积判断 if ((((l1x1-l2x1)*(l2y2-l2y1)-(l1y1-l2y1)*(l2x2-l2x1))*... ((l1x2-l2x1)*(l2y2-l2y1)-(l1y2-l2y1)*(l2x2-l2x1))) > 0 ||... (((l2x1-l1x1)*(l1y2-l1y1)-(l2y1-l1y1)*(l1x2-l1x1))*... ((l2x2-l1x1)*(l1y2-l1y1)-(l2y2-l1y1)*(l1x2-l1x1))) > 0) %如果判断为真，则不会相交 panduan=0; else panduan=1; end end end %两个向量做差积 function t=crossdot(xy1,xy2) x1=xy1(1);y1=xy1(2); x2=xy2(1);y2=xy2(2); t=x1*y2-y1*x2; end %点是否在三角形内 function panduan=pointintriangle(xy1,xy2,xy3,xy0) x1=xy1(1);y1=xy1(2); x2=xy2(1);y2=xy2(2); x3=xy3(1);y3=xy3(2); x0=xy0(1);y0=xy0(2); PA=[x1-x0,y1-y0];PB=[x2-x0,y2-y0];PC=[x3-x0,y3-y0]; %利用差积同正或同负号来判断是否在三角内 t1=crossdot(PA,PB); t2=crossdot(PB,PC); t3=crossdot(PC,PA); if abs(sign(t1)+sign(t2)+sign(t3))==3 panduan=1; else panduan=0; end end %点是否在多边形内 function panduan=pointinmutiangle(xdot,d_no,xy0) %d_no符合12341的格式，收尾相连 Ndot=xdot(d_no,:); PN=[Ndot(:,1)-xy0(1),Ndot(:,2)-xy0(2)]; tn=zeros(length(d_no)-1,1); for j=1:length(d_no)-1 tn(j)=crossdot(PN(j,:),PN(j+1,:)); end %利用差积同正或同负号来判断是否在三角内 if abs(sum(sign(tn)))==length(d_no)-1 panduan=1; else panduan=0; end end %做出同时具有一个中心点的一圈三角形，按照相接顺序排序 function [temp_trimat,temp_trimat_no,panduan]=findpoint2tri(j,trimat) temp_trimat_no=1:size(trimat,1); panduan=1; temp_trimat=trimat(logical(sum(trimat==j,2)),:); temp_trimat_no=temp_trimat_no(logical(sum(trimat==j,2))); %把j放到每行第一个 for k=1:size(temp_trimat,1) temp_trimat(k,[find(temp_trimat(k,:)==j),1])=temp_trimat(k,[1,find(temp_trimat(k,:)==j)]); end %如果有一个点只出现过一次，把这个点包含的三角形放到第一行 tN2=tabulate(reshape(temp_trimat,[],1)); tN3=tN2((tN2(:,2)==1),1); %其实如果中心点不在范围内也最好单列出来，挖坑 if size(tN3,1)~=0 tN3=tN3(1); tN=find(sum(temp_trimat==tN3,2)); temp_trimat([1,tN],:)=temp_trimat([tN,1],:); temp_trimat_no([1,tN])=temp_trimat_no([tN,1]); panduan=0; if temp_trimat(1,3)==tN3 temp_trimat(1,[2,3])=temp_trimat(1,[3,2]); end end %把边首尾排序一遍 for k=2:size(temp_trimat,1) tN=find(sum(temp_trimat==temp_trimat(k-1,3),2)); tN(tN==k-1)=[]; temp_trimat([tN,k],:)=temp_trimat([k,tN],:); if temp_trimat(k-1,3)~=temp_trimat(k,2) temp_trimat(k,[2,3])=temp_trimat(k,[3,2]); end temp_trimat_no([tN,k])=temp_trimat_no([k,tN]); end end %做出一点对边缘的中垂线，得到边缘点坐标 function tempboderdot=edgepointfind(tempedge,xdot,xy0,border_point) x0=xy0(1);y0=xy0(2); %判断中心点是否在大图形内 d_no=[border_point(1,1);border_point(:,2)]; panduan=pointinmutiangle(xdot,d_no,xy0); %求边的中点 xz=(xdot(tempedge(1),1)+xdot(tempedge(2),1))/2; yz=(xdot(tempedge(1),2)+xdot(tempedge(2),2))/2; %做中心点延长线，得到一条超长线段定为2得了,用到了边界为1这个条件 panduan2=true; if panduan %在内部，做中心到边缘延长线即可 t_xy1=[x0,y0]+[xz-x0,yz-y0]*2/sqrt((xz-x0)^2+(yz-y0)^2); else if ~(x0<0||x0>1||y0<0||y0>1) %在外部,但是没超出边界，做中心到反向边缘 t_xy1=[x0,y0]-[xz-x0,yz-y0]*2/sqrt((xz-x0)^2+(yz-y0)^2); else %在边界外 panduan2=false; tempboderdot=[];%返回空矩阵 end end if panduan2 %判断是4个边哪一个 [xy3,xy4]=select4edge(xy0,t_xy1); %做4点求交点运算 tempboderdot=crosspoint(xy0,t_xy1,xy3,xy4); end % panduan % d_no % xy0 % t_xy1 % tempboderdot end %做出边缘界外连线，得到边缘点坐标 function tempboderdot1=edgepointfind2(xy1,xy2) [xy3,xy4]=select4edge(xy1,xy2); tempboderdot1=crosspoint(xy1,xy2,xy3,xy4); end %判断射线会与哪条边相交 function [xy3,xy4]=select4edge(xy0,xy1) deg01=angle((0-xy0(1))+(0-xy0(2))*1i);%左下00 deg02=angle((1-xy0(1))+(0-xy0(2))*1i);%10 deg03=angle((1-xy0(1))+(1-xy0(2))*1i);%11 deg04=angle((0-xy0(1))+(1-xy0(2))*1i);%01 deg00=angle((xy1(1)-xy0(1))+(xy1(2)-xy0(2))*1i); [~,I]=sort([deg00,deg01,deg02,deg03,deg04]); k=find(I==1); switch k case 1 xy3=[0,0];xy4=[0,1]; case 2 xy3=[0,0];xy4=[1,0]; case 3 xy3=[1,0];xy4=[1,1]; case 4 xy3=[0,1];xy4=[1,1]; case 5 xy3=[0,0];xy4=[0,1]; end end %求两条线交点 function xy0=crosspoint(xy1,xy2,xy5,xy6) a1=(xy2(2)-xy1(2));b1=(xy1(1)-xy2(1));c1=xy1(1)*xy2(2)-xy2(1)*xy1(2); a2=(xy6(2)-xy5(2));b2=(xy5(1)-xy6(1));c2=xy5(1)*xy6(2)-xy6(1)*xy5(2); xy0=[det([c1,b1;c2,b2])/det([a1,b1;a2,b2]),det([a1,c1;a2,c2])/det([a1,b1;a2,b2])]; end %删除所有中心点超出边界的三角形 function [temp_trimat,temp_trimat_no,delleft,delright]=selecttemp_trimat(temp_trimat,temp_trimat_no,trimatcenter) delleft=0;delright=0; for j=1:size(temp_trimat,1) centerdot=trimatcenter(temp_trimat_no(j),:); if and(and(0<centerdot(1),centerdot(1)<1),and(0<centerdot(2),centerdot(2)<1)) break end end if j~=1 delleft=temp_trimat_no(j-1); end temp_trimat(1:j-1,:)=[]; temp_trimat_no(1:j-1)=[]; %倒着来一遍 temp_trimat=flipud(temp_trimat); temp_trimat_no=fliplr(temp_trimat_no); for j=1:size(temp_trimat,1) centerdot=trimatcenter(temp_trimat_no(j),:); if and(and(0<centerdot(1),centerdot(1)<1),and(0<centerdot(2),centerdot(2)<1)) break end end if j~=1 delright=temp_trimat_no(j-1); end temp_trimat(1:j-1,:)=[]; temp_trimat_no(1:j-1)=[]; temp_trimat=flipud(temp_trimat); temp_trimat_no=fliplr(temp_trimat_no); end %利用两个边缘点求角点 function tempboderdot3=maketempboderdot(tempboderdot1,tempboderdot2) tempboderdot3=zeros(1,2); %第一个边 if abs(tempboderdot1(1)-0)<1e-10 tempboderdot3(1)=0; end if abs(tempboderdot1(1)-1)<1e-10 tempboderdot3(1)=1; end if abs(tempboderdot1(2)-0)<1e-10 tempboderdot3(2)=0; end if abs(tempboderdot1(2)-1)<1e-10 tempboderdot3(2)=1; end %第二个边 if abs(tempboderdot2(1)-0)<1e-10 tempboderdot3(1)=0; end if abs(tempboderdot2(1)-1)<1e-10 tempboderdot3(1)=1; end if abs(tempboderdot2(2)-0)<1e-10 tempboderdot3(2)=0; end if abs(tempboderdot2(2)-1)<1e-10 tempboderdot3(2)=1; end end 

输出的网格

未上色的Voronoi图

彩色Voronoi图

结束