前言

原文: 《泰森多邊形（Voronoi圖）的matlab繪制》.
本文已經過原作者授權。如有錯誤，請批評指正。

泰森多邊形介紹

泰森多邊形是對空間平面的一種剖分，其特點是多邊形內的任何位置離該多邊形的樣點（如居民點）的距離最近，離相鄰多邊形內樣點的距離遠，且每個多邊形內含且僅包含一個樣點。由于泰森多邊形在空間剖分上的等分性特征，因此可用于解決最近點、最小封閉圓等問題，以及許多空間分析問題，如鄰接、接近度和可達性分析等。

法一：matlab內部函數

%采用matlab自帶的函數進行繪制
clear
xdot=gallery('uniformdata',[200 2],5);
%delaunay三角形
figure(1)
DT=delaunayTriangulation(xdot);
triplot(DT,'color','k');title('Delaulay三角形');
%voronoi三角形
figure(2)
voronoi(xdot(:,1),xdot(:,2));
xlim([0,1])
ylim([0,1])
title('泰森多邊形')

代碼稍微改了一下，給圖片加了個標題。
附上效果圖：

法二：手寫函數

clear;clc;close all;
N=100;
%點隨機
xdot=rand(N,2);
%1Delaulay三角形的構建
%整理點，遵循從左到右，從上到下的順序
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+3pointtemp=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、圓右側2switch panduancase 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三角形，保存到trianglestrimat=[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%三角循環結束    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]);title('Delaulay三角形');%凸包監測
%思路是先找出邊緣點（三角形只有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)-1tempboderedge=[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_Nif tempdotex(k,1)==tempboderdott_xdotno4=tempdotex(k,2);elset_xdotno4=tempdotex(k,1);endtt_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==1t_t=t_t+1;breakendend%t_t大于0說明有相交的線,略過if t_t==0temp_edgemat=[temp_edgemat;tempboderedge];temp_trimat=[temp_trimat;[tempboderedge,tempboderdot]];breakendendendtrimat=[trimat;temp_trimat];edgemat=[edgemat;temp_edgemat];%刪除重復的三角形trimat=sort(trimat,2);trimat=unique(trimat,'stable','rows');if j==size(border_point,1)-1t_w=1;end
endfigure(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]);title('Delaulay三角形');%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%求三角形的相鄰三角形個數
[border_trimat,border_point,trimat_con]=makebordertri(trimat);
Thi_edge1=[];
for j=1:size(trimat,1)tempedge=[];%第一個相鄰三角形if trimat_con(j,1)~=0tempedge=[tempedge;[j,trimat_con(j,1)]];end%第二個相鄰三角形if trimat_con(j,2)~=0tempedge=[tempedge;[j,trimat_con(j,2)]];end%第三個相鄰三角形if trimat_con(j,3)~=0tempedge=[tempedge;[j,trimat_con(j,3)]];endThi_edge1=[Thi_edge1;tempedge];
end%繪制非邊緣泰勒多邊形
figure(2)
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)%判斷點是否在邊與邊框的三角內，如果在，做中心的延長線%如果不在，做中心-邊緣的延長線%或者改成判斷點是否在多邊形內panduan2=pointinmutiangle(xdot,[border_point(1,1);border_point(:,2)],[t_x1,t_y1]);if panduan2==1t_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])elset_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])endend
endscatter(xdot(:,1),xdot(:,2),5,[0,0.4,0],'filled')
hold off;title('泰森多邊形')%判斷點在三角形外接圓的哪個部分
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,:)];elsetemp_edge(1,:)=[dot0,dot1,xdot(dot0,:),xdot(dot1,:)];end
elsetemp_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,:)];elsetemp_edge(2,:)=[dot0,dot2,xdot(dot0,:),xdot(dot2,:)];end
elsetemp_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,:)];elsetemp_edge(3,:)=[dot0,dot3,xdot(dot0,:),xdot(dot3,:)];end
elsetemp_edge(3,:)=[dot0,dot3,xdot(dot0,:),xdot(dot3,:)];
endend%做出一些列固定點發散的線段外點組成的三角形
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:Ntemp_trimat(j,:)=[pointline(index(j)),pointline(index(j+1)),dot0];
endend%將三個點構成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,:)];elseedgemat(1,:)=[dot2,dot1,xdot(dot2,:),xdot(dot1,:)];end
elseedgemat(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,:)];elseedgemat(2,:)=[dot3,dot1,xdot(dot3,:),xdot(dot1,:)];end
elseedgemat(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,:)];elseedgemat(3,:)=[dot2,dot3,xdot(dot2,:),xdot(dot3,:)];end
elseedgemat(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)==0border_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)==0border_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];endend%把邊首尾排序一遍,輸出border_point
for j=1:size(border_point,1)-1border_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
endend%判斷兩個線段是否相交
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;elsepanduan=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))==3panduan=1;
elsepanduan=0;
endend%點是否在多邊形內
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)-1tn(j)=crossdot(PN(j,:),PN(j+1,:));
end
%利用差積同正或同負號來判斷是否在三角內if abs(sum(sign(tn)))==length(d_no)-1panduan=1;
elsepanduan=0;
endend

這個代碼自己稍微修改了一下，給圖片加了個標題。可以直接復制粘貼運行。
附上兩張效果圖：