自定义函数实现delaunayTriangulation 使用Bowyer-Watson 算法

function [mesh,points] = mydelaunayTriangulation(x,y) % x = [89 830 428 663 69]; % y = [342,197,148,560,576]; global Edge Circum_cricle Triangle x = reshape(x,[],1); y = reshape(y,[],1); Point = struct('x',0,'y',0); for i = 1:length(x) points(i) = Point; points(i).x = x(i); points(i).y = y(i); end x_max = max(x); y_max = max(y); x_min = min(x); y_min = min(y); % dx = (x_max - x_min)*10; % dy = (y_max - y_min)*10; dx = (x_max - x_min)*1; dy = (y_max - y_min)*1; Edge.points(1) = Point; Edge.points(2) = Point; % Edge.point_distance = 0; % Circum_cricle.circum_center = Point; % Circum_cricle.radius = 0; Triangle.edges(1) = Edge; Triangle.edges(2) = Edge; Triangle.edges(3) = Edge; % Triangle.circum_circle = Circum_cricle; p1 = Point; p2 = Point; p3 = Point; p1.x = x_min-dx; p1.y = y_min-dy*3; p2.x = x_min-dx; p2.y = y_max+dy; p3.x = x_max +dx*3; p3.y = y_max+dy; %求三个点的CircumCircle super_triangle = Triangle; super_triangle.edges(1).points(1) = p1; super_triangle.edges(1).points(2) = p2; super_triangle.edges(2).points(1) = p2; super_triangle.edges(2).points(2) = p3; super_triangle.edges(3).points(1) = p3; super_triangle.edges(3).points(2) = p1; % hold on % plot_tri(p1, p2, p3); % [circum_center,radius] = CircumCircle(p1,p2,p3); % super_triangle.circum_circle.circum_center = circum_center; % super_triangle.circum_circle.radius = radius; triangles{1} = super_triangle; for i = 1:length(points) point = points(i); triangles = AddPointToMesh(point, triangles); end triangles = RemoveSuperTriangle(triangles, super_triangle); mesh = triangles; end function mesh_triangles = RemoveSuperTriangle(triangles,super_triangle) super_tri_p1 = super_triangle.edges(1).points(1); super_tri_p2 = super_triangle.edges(2).points(1); super_tri_p3 = super_triangle.edges(3).points(1); k = 1; for i = 1:length(triangles) triangle = triangles{i}; p1 = triangle.edges(1).points(1); p2 = triangle.edges(2).points(1); p3 = triangle.edges(3).points(1); if ~ ( is_point_same(p1,super_tri_p1) || is_point_same(p1 ,super_tri_p2) || is_point_same(p1, super_tri_p3) || ... is_point_same(p2 ,super_tri_p1) || is_point_same(p2,super_tri_p2) || is_point_same(p2,super_tri_p3) || ... is_point_same(p3, super_tri_p1) || is_point_same(p3,super_tri_p2) || is_point_same(p3, super_tri_p3)) mesh_triangles{k} = triangle; k = k +1; end end end function ret = is_point_same(p1,p2) ret = 0; if(p1.x==p2.x)&&(p1.y==p2.y) ret =1; end end function good_triangles = AddPointToMesh(point,triangles) global Edge Circum_cricle Triangle k1 = 1;k2 =1; for i = 1:length(triangles) triangle = triangles{i}; % plot_triangle(triangle); % plot_circumcircle(triangle); if IsInCircumCircle(point,triangle) edges{k1} = triangle.edges(1); edges{k1+1} = triangle.edges(2); edges{k1+2} = triangle.edges(3); k1 = k1+3; else good_triangles{k2} = triangle; k2 = k2+1; end end edges = GetUniqueEdges(edges); for i = 1:length(edges) edge = edges{i}; e1 = Edge; e2 = Edge; e3 = Edge; e1.points(1) = edge.points(1); e1.points(2) = edge.points(2); e2.points(1) = edge.points(2); e2.points(2) = point; e3.points(1) = point; e3.points(2) = edge.points(1); triangle = Triangle; triangle.edges(1) = e1; triangle.edges(2) = e2; triangle.edges(3) = e3; good_triangles{k2} = triangle; % plot_triangle(triangle); k2 = k2+ 1; end end function unique_edges = GetUniqueEdges(edges) k = 1; for i = 1:length(edges) edge1 = edges{i}; is_unique_edge = 1; for j = 1:length(edges) edge2 = edges{j}; if i ~= j && is_edge_same(edge1,edge2) is_unique_edge = 0; break; end end if is_unique_edge unique_edges{k} = edge1; k = k +1; end end end function ret = is_edge_same(edge1,edge2) ret = 0; if (edge1.points(1).x == edge2.points(1).x ... && edge1.points(1).y == edge2.points(1).y ... && edge1.points(2).x == edge2.points(2).x ... && edge1.points(2).y == edge2.points(2).y ) ret = 1; end if (edge1.points(2).x == edge2.points(1).x ... && edge1.points(2).y == edge2.points(1).y ... && edge1.points(1).x == edge2.points(2).x ... && edge1.points(1).y == edge2.points(2).y ) ret = 1; end end function ret = IsInCircumCircle(point,triangle) p1 = triangle.edges(1).points(1); p2 = triangle.edges(2).points(1); p3 = triangle.edges(3).points(1); [circum_center,radius] = CircumCircle(p1,p2,p3); triangle.circum_circle.circum_center = circum_center; triangle.circum_circle.radius = radius; dist = GetDistance(point, triangle.circum_circle.circum_center); if dist < triangle.circum_circle.radius ret = 1; else ret = 0; end end function [circum_center,radius] = CircumCircle(p1,p2,p3) [a1,b1,c1] = GetLineEquationFromPoints(p1, p2); [a2,b2,c2] = GetLineEquationFromPoints(p2, p3); [a1, b1, c1] = GetPerpendicularBisectorLineEquation(p1, p2, a1, b1, c1); [a2, b2, c2] = GetPerpendicularBisectorLineEquation(p2, p3, a2, b2, c2); circum_center = GetLinesIntersection(a1, b1, c1, a2, b2, c2); side1 = GetDistance(p1,p2); side2 = GetDistance(p1,p3); side3 = GetDistance(p2,p3); radius = GetCircumRadius(side1,side2,side3); end function ret = GetDistance(p1,p2) if isempty(p2) ret = sqrt((p1.x-0)*(p1.x-0)+(p1.y-0)*(p1.y-0)); return ; end ret = sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)); end function ret = GetCircumRadius(side1,side2,side3) num = side1 * side2 * side3; a = (side1 + side2 + side3); b = (-side1 + side2 + side3); c = (side1 - side2 + side3); d = (side1 + side2 - side3); den = sqrt(abs(a * b * c * d)); ret = num / den; end function [a,b,c] = GetLineEquationFromPoints(p1,p2) a = p2.y - p1.y; b = p1.x - p2.x; c = a*p1.x + b*p1.y; end function [a,b,c]= GetPerpendicularBisectorLineEquation( p1, p2, a, b, c) mid_point.x = (p1.x + p2.x) / 2; mid_point.y = (p1.y + p2.y) / 2; c = -b * mid_point.x + a * mid_point.y; temp = a; a = -b; b = temp; end function P = GetLinesIntersection(a1, b1, c1, a2, b2, c2) determinant = a1 * b2 - a2 * b1; if (determinant == 0) P =[]; else P.x = (b2 * c1 - b1 * c2) / determinant; P.y = (a1 * c2 - a2 * c1) / determinant; end end function plot_tri(p1, p2, p3) % 提取点的坐标 x = [p1.x, p2.x, p3.x, p1.x]; y = [p1.y, p2.y, p3.y, p1.y]; % 绘制三角形边 plot(x, y, 'b-', 'LineWidth', 2); hold on; % 绘制顶点 plot(p1.x, p1.y, 'ro', 'MarkerFaceColor', 'r'); plot(p2.x, p2.y, 'ro', 'MarkerFaceColor', 'r'); plot(p3.x, p3.y, 'ro', 'MarkerFaceColor', 'r'); end function plot_triangle(triangle) p1 = triangle.edges(1).points(1); p2 = triangle.edges(2).points(1); p3 = triangle.edges(3).points(1); plot_tri(p1, p2, p3); end function plot_circumcircle(triangle) p1 = triangle.edges(1).points(1); p2 = triangle.edge