✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
不规则三角网(Triangulated Irregular Network,TIN)在表示地形的形态方面具有较好的表现,其生成算法一直备受关注。讨论了三角网的数据结构的设计,采用逐点插入算法中的Bowyer-Watson算法思想为研究重点,设计并实现了该算法,对算法实验过程中可能出现的交叉现象进行分析,给出算法的改进。该改进算法已用于地形的可视化建模中,获得了较好的效果,对于三角剖分的相关研究具有一定的价值。
⛄ 代码
clc
clear
close all
%% Bowyer-Watson算法复现(逐点插入)
Pts = rand(20,2);
% save Pts.mat Pts
% load Pts.mat Pts
Pts0 = Pts;
% figure,plot(Pts(:,1),Pts(:,2),'b.')
%% 建立最小外接矩形
MBR = [min(Pts(:,1))-0.5,max(Pts(:,2))+0.5;...
min(Pts(:,1))-0.5,min(Pts(:,2))-0.5;...
max(Pts(:,1))+0.5,max(Pts(:,2))+0.5;...
max(Pts(:,1))+0.5,min(Pts(:,2))-0.5];
Pts = [MBR;Pts];%在点集中添加MBR;
Del = [1,2,3;2,3,4];%建立辅助窗口
%% 逐点插入
for i = 5:size(Pts,1)
flag = zeros(1,size(Del,1));%点的影响范围flag
for j = 1:size(Del,1)
if influence(i,Pts,Del(j,:))==true %判断点是否在三角形外接圆内
flag(j)=1;
end
end
flag = flag>0;
a = Del(flag,:);
Del = Del(~flag,:);
a = a(:);
convex = unique(a);% Delaunay腔
% 按角度顺次连接凸包顶点，生成新三角形
Del_new = newtriangle(convex,i,Pts);
% 局部最优化
for j = 2:size(Del_new,1)
tri1 = Del_new(j-1,:);
tri2 = Del_new(j,:);
[~,center] = influence(i,Pts,tri1);
pt = Pts(setdiff(tri2,tri1),:);
pt1 = Pts(tri1(1),:);
if norm(pt-center)<norm(pt1-center)
ipt = intersect(tri1,tri2);
Del_new(j-1,:) = [setdiff(tri2,tri1),setdiff(tri1,tri2),setdiff(ipt,i)];
Del_new(j,:) = [setdiff(tri2,tri1),setdiff(tri1,tri2),i];
end
end
tri1 = Del_new(end,:);
tri2 = Del_new(1,:);
ipt = intersect(tri1,tri2);
if length(ipt)>1
[~,center] = influence(i,Pts,tri1);
pt = Pts(setdiff(tri2,tri1),:);
pt1 = Pts(tri1(1),:);
if norm(pt-center)<norm(pt1-center)
Del_new(end,:) = [setdiff(tri2,tri1),setdiff(tri1,tri2),setdiff(ipt,i)];
Del_new(1,:) = [setdiff(tri2,tri1),setdiff(tri1,tri2),i];
end
end
Del = [Del;Del_new];
% x = [Pts(Del(:,1),1),Pts(Del(:,2),1),Pts(Del(:,3),1),Pts(Del(:,1),1)];
% y = [Pts(Del(:,1),2),Pts(Del(:,2),2),Pts(Del(:,3),2),Pts(Del(:,1),2)];
% figure;
% for ii = 1:size(x,1)
% plot(x(ii,:),y(ii,:),'b-')
% hold on
% end
% plot(Pts(:,1),Pts(:,2),'bo')
end
%% 删除辅助点
Del0 = Del;
Del(Del(:,1)<5,:)=[];
Del(Del(:,2)<5,:)=[];
Del(Del(:,3)<5,:)=[];
%% 绘制结果三角形
x = [Pts(Del(:,1),1),Pts(Del(:,2),1),Pts(Del(:,3),1),Pts(Del(:,1),1)];
y = [Pts(Del(:,1),2),Pts(Del(:,2),2),Pts(Del(:,3),2),Pts(Del(:,1),2)];
figure;
for i = 1:size(x,1)
plot(x(i,:),y(i,:),'b-')
hold on
end
plot(Pts0(:,1),Pts0(:,2),'bo')
title('三角化')
Pts1 = Pts(5:end,:);
tri = delaunay(Pts1(:,1),Pts1(:,2));
figure,triplot(tri,Pts1(:,1),Pts1(:,2));
title('与matlab内建delaunay函数结果做对比')
%% Voronoi图
Del = Del0;% 重新把MBR加上
center = zeros(size(Del,1),2);
Voronoi = [];
Voronoi2 = [];
for i = 1:size(Del,1)
Del1 = Del(i,:);
[~,center0] = influence(1,Pts,Del1);
center(i,:) = center0;
end
for i = 1:size(Del,1)
Del1 = Del(i,:);
for j = 1:size(Del,1)
if i==j
continue
end
Del2 = Del(j,:);
ipt = intersect(Del1,Del2);
if length(ipt)>1
Voronoi = [Voronoi;i,j];
end
end
end
% for i = 1:size(Del,1)
% Del1 = Del(i,[1,2]);
% flag = false;
% for j = 1:size(Del,1)
% if i==j
% continue
% end
% Del2 = Del(j,:);
% ipt = intersect(Del1,Del2);
% if length(ipt)>1
% Voronoi = [Voronoi;i,j];
% flag = true;
% end
% end
% if ~flag
% Voronoi2 = [Voronoi2;Del(i,[1,2,3])];
% end
% Del1 = Del(i,[2,3]);
% flag = false;
% for j = 1:size(Del,1)
% if i==j
% continue
% end
% Del2 = Del(j,:);
% ipt = intersect(Del1,Del2);
% if length(ipt)>1
% Voronoi = [Voronoi;i,j];
% flag = true;
% end
% end
% if ~flag
% Voronoi2 = [Voronoi2;Del(i,[2,3,1])];
% end
% Del1 = Del(i,[1,3]);
% flag = false;
% for j = 1:size(Del,1)
% if i==j
% continue
% end
% Del2 = Del(j,:);
% ipt = intersect(Del1,Del2);
% if length(ipt)>1
% Voronoi = [Voronoi;i,j];
% flag = true;
% end
% end
% if ~flag
% Voronoi2 = [Voronoi2;Del(i,[1,3,2])];
% end
%
% end
% pt1 = Pts(Voronoi2(:,1),:);
% pt2 = Pts(Voronoi2(:,2),:);
% pt0 = (pt1+pt2)/2;
% pt3 = Pts(Voronoi2(:,3),:);
% vot = pt1-pt2;
% vot = vot./sqrt(sum(vot.^2,2));
% vot = [-vot(:,2),vot(:,1)];
% vot0 = pt0-pt3;
% if sum(vot.*vot0)<0
% vot = -vot;
% end
% pt1 = pt0+vot;%为凸包画垂线
% Voronoi2 = [pt1,pt0];
figure;
for i = 1:size(Voronoi,1)
plot([center(Voronoi(i,1),1),center(Voronoi(i,2),1)],...
[center(Voronoi(i,1),2),center(Voronoi(i,2),2)],'b-');
hold on
end
plot(Pts0(:,1),Pts0(:,2),'bo')
% for i = 1:size(pt1,1)
% plot([pt1(i,1),pt0(i,1)],[pt1(i,2),pt2(i,2)],'b-')
% hold on
% end
axis([min(Pts0(:,1)),max(Pts0(:,1)),min(Pts0(:,2)),max(Pts0(:,2))]);
axis equal
title('泰森多边形')
%% 判断插入点是否在一个三角形的外接圆内
function [flag,center] = influence(i,Pts,Del)
Del = Pts(Del,:);
pt = Pts(i,:);
pt1 = Del(1,:);
pt2 = Del(2,:);
pt3 = Del(3,:);
% [a1,b1;a2,b2]*[x;y]=[c1;c2] 外心计算公式
a1 = 2*(pt2(1)-pt1(1));
b1 = 2*(pt2(2)-pt1(2));
c1 = pt2(1).^2+pt2(2).^2-pt1(1).^2-pt1(2).^2;
a2 = 2*(pt3(1)-pt2(1));
b2 = 2*(pt3(2)-pt2(2));
c2 = pt3(1).^2+pt3(2).^2-pt2(1).^2-pt2(2).^2;
center = [a1,b1;a2,b2]\[c1;c2];
center = transpose(center);
% x = (c1*b2-c2*b1)/(a1*b2-a2*b1);
% y = (a1*c2-a2*c1)/(a1*b2-a2*b1);
% center = [x,y];
flag = norm(pt-center)<norm(pt1-center);
end
%% 按角度排排坐，分果果
function Del_new = newtriangle(convex,i,Pts)
pt = Pts(i,:);
convex0 = convex;
convex = Pts(convex,:)-repmat(pt,length(convex0),1);
convex = convex./repmat(sqrt(sum(convex.^2,2)),1,2);
theta = acos(convex(:,2));
theta(convex(:,1)<0)=2*pi-theta(convex(:,1)<0);
[~,I]=sort(theta);
Del_new = zeros(0,3);
for j = 2:length(I)
Del_new = [Del_new;i,convex0(I(j-1)),convex0(I(j))];
end
Del_new = [Del_new;i,convex0(I(1)),convex0(I(end))];
end
⛄ 运行结果
⛄ 参考文献
[1] Chrisochoides N , Sukup F . Task parallel implementation of the Bowyer-Watson algorithm[J]. mississippi state univ mississippi state ms, 1999.
[2] 成俊燕. 基于单张图片的服装建模相关算法研究[D]. 浙江大学.
[3] 李景焕. 基于Bowyer-Watson法的Delaunay三角网格的一些改进[J]. 天津商学院学报, 2007, 27(3):33-37.
[4] 周雪梅, 黎应飞. 基于Bowyer-Watson三角网生成算法的研究[J]. 计算机工程与应用, 2013, 49(6):198-200.
