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
