Delaunay三角剖分使我失去对称性

Question

我正在使用Delaunay三角剖分法将多边形拆分为三角形。我使用大码处理有限元，我的一个“检查点”是对称的（如果数据是对称的，输出也必须是对称的）。但是，由于我无法控制Delaunay三角剖分，所以我失去了对称性。

我写了一个小代码来说明我的问题：我们考虑两个不相交的三角形和一个与它们相交的大矩形。我们希望与三角矩形这些三角形的减法：

clear all 
close all 
warning off % the warning is about duplicate points, not important here 

figure 
hold on 

p =[.3 .3 
.4 .3 
.3 .4 
.7 .6 
.6 .7 
.7 .7]; % coordinates of the points for the triangles 

px = 1/3; 
py = 1/3; 
lx = 1/3; 
ly = 1/3; % size and position of the rectangle 

% rearrange the polygon with clockwise-ordered vertices 
[x1,y1]=poly2cw([px; px+lx; px+lx; px],[py; py; py+ly; py+ly]); % rectangle 

patch(x1,y1, 1, 'EdgeColor', 'k'); 

for i=1:2 

    pc = p(3*i-2:3*i,:); % current triangle 
    % rearrange the polygon with clockwise-ordered vertices 
    [x0,y0]=poly2cw(pc(:,1),pc(:,2)); % triangle 

    [x2,y2] = polybool('intersection',x1,y1,x0,y0); % intersection 
    [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % subtraction 

    DT = delaunayTriangulation(x3,y3); 

    triplot(DT,'Marker','o') 

end 
XL = xlim; xlim(XL+[-1 +1]*diff(XL)/10); 
YL = ylim; ylim(YL+[-1 +1]*diff(YL)/10); 
axis equal; 
box on; 

正如你所看到的，Delaunay三角不必在两个三角形相同的行为，对称的，因此损失。

有没有简单的方法来恢复对称性？

我使用Matlab R2013a。

Answer 1

好像你使用的是MatLab R2013，因为在我的R2011b中没有delaunayTriangulation函数。为了能够运行你的代码我改变了它稍微：

clear all 
close all 
warning off % the warning is about duplicate points, not important here 

figure 
hold on 

p =[.3 .3 
    .4 .3 
    .3 .4 
    .7 .6 
    .6 .7 
    .7 .7]; % coordinates of the points for the triangles 

px = 1/3; 
py = 1/3; 
lx = 1/3; 
ly = 1/3; % size and position of the rectangle 

% rearrange the polygon with clockwise-ordered vertices 
[x1,y1]=poly2cw([px; px+lx; px+lx; px],[py; py; py+ly; py+ly]); % rectangle 

patch(x1,y1, 1, 'EdgeColor', 'k'); 

for i=1:2 

    pc = p(3*i-2:3*i,:); % current triangle 
    % rearrange the polygon with clockwise-ordered vertices 
    [x0,y0]=poly2cw(pc(:,1),pc(:,2)); % triangle 

    [x2,y2] = polybool('intersection',x1,y1,x0,y0); % intersection 
    [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % subtraction 

    %DT = delaunayTriangulation(x3,y3); 
    %triplot(DT) 

    % This is triangulation of subtraction 
    DT = delaunay(x3,y3); 
    triplot(DT,x3,y3, 'Marker','.', 'Color','r') 

    % This is triangulation of intersection 
    DT = delaunay(x2,y2); 
    triplot(DT,x2,y2, 'Marker','o', 'Color','b', 'LineWidth',1) 

end 
axis equal; 
axis tight; 
box on; 

XL = xlim; xlim(XL+[-1 +1]*diff(XL)/10); 
YL = ylim; ylim(YL+[-1 +1]*diff(YL)/10); 

text(0.5,0.55,'triangulation of subtraction', 'HorizontalAlignment','center', 'VerticalAlignment','bottom', 'Color','r'); 
text(0.5,0.45,'triangulation of intersection', 'HorizontalAlignment','center', 'VerticalAlignment','top', 'Color','b'); 

下面是结果我看到

是否与你相似？ 你可以在你的问题中添加一张图片来描述你得到的结果有什么问题吗？

Answer 2

看来这不是一个错误。 实际上，您的结果来自您的数据。

我打得有点与您的代码

clear all 
close all 
warning off % the warning is about duplicate points, not important here 

figure 
hold on 

p =[.3 .3 
.4 .3 
.3 .4 
.7 .6 
.6 .7 
.7 .7]; % coordinates of the points for the triangles 

px = 1/3; 
py = 1/3; 
lx = 1/3; 
ly = 1/3; % size and position of the rectangle 

% rearrange the polygon with clockwise-ordered vertices 
[x1,y1]=poly2cw([px; px+lx; px+lx; px],[py; py; py+ly; py+ly]); % rectangle 

patch(x1,y1, 1, 'EdgeColor', 'k'); 

for i=1:2 

    pc = p(3*i-2:3*i,:); % current triangle 
    % rearrange the polygon with clockwise-ordered vertices 
    [x0,y0]=poly2cw(pc(:,1),pc(:,2)); % triangle 

    [x2,y2] = polybool('intersection',x1,y1,x0,y0); % intersection 
    [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % subtraction 

    % This is for R2013a 
    %DT = delaunayTriangulation(x3,y3); 
    %triplot(DT,'Marker','o'); 

    % This is for R2011b 
    %DT = DelaunayTri(x3,y3); 
    %triplot(DT,'Marker','o'); 

    % This is plain delaunay version 
    DT = delaunay(x3,y3); 
    triplot(DT,x3,y3,'Marker','o') 

    % we break here to analyze the first triangulation 
    break 

end 
XL = xlim; xlim(XL+[-1 +1]*diff(XL)/10); 
YL = ylim; ylim(YL+[-1 +1]*diff(YL)/10); 
axis equal; 
box on; 

% % % % % % % % % % % % % % % % % % 
% Checking the triangulation 
% % % % % % % % % % % % % % % % % % 

% Wrong triangulation for i=2 is hard-coded 
DT2  = [ 
    2 1 6 
    6 5 2 
    5 3 2 
    5 4 3 
    2 3 1 ]; 

figure; 
hold all; 
triplot(DT2,x3,y3,'Marker','o', 'Color','r', 'LineWidth',1) 
axis equal; 
axis tight; 
box on; 
XL = xlim; xlim(XL+[-1 +1]*diff(XL)*0.5); 
YL = ylim; ylim(YL+[-1 +1]*diff(YL)*0.5); 

% circumcircle: http://www.mathworks.com/matlabcentral/fileexchange/17300 
ca = linspace(0,2*pi); 
cx = cos(ca); 
cy = sin(ca); 

hl = []; 
for k=1:size(DT2,1) 
    tx = x3(DT(k,:)); 
    ty = y3(DT(k,:)); 
    [r,cn]=circumcircle([tx,ty]',0); 
    if ~isempty(hl) 
     %delete(hl); 
    end 
    fprintf('Circle %d: Center at (%.23f,%.23f); R=%.23f\n',k,cn,r); 
    text(cn(1),cn(2), sprintf('c%d',k)); 
    hl = plot(cx*r+cn(1), r*cy+cn(2), 'm-'); 
    drawnow; 
    %pause(3); %if you prefere to go slowly 
end 

这是输出我SEEE：

圈1：中心（0.28333333333333333000000,0.35000000000000003000000）; R = 0.05270462766947306400000 Circle 2：Center at（0.34999999999999998000000,0.34999999999999998000000）; R = 0.02357022603955168100000 Circle 3：Center at（0.28333333333333338000000,0.34999999999999992000000）; R = 0.05270462766947289800000 圈4：中心在（0.35000000000000003000000,0.28333333333333355000000）; R = 0.05270462766947290500000 圈5：中心在（0.35000000000000003000000,0.28333333333333333000000）; R = 0.05270462766947312000000

而图：

所以圆圈1和3以及圆圈4和5几乎相同。 所以，你和我的结果之间的差异甚至可能来自舍入误差，因为相应的四个点在浮点数学精度内位于同一个圆上。你必须重新设计你的观点，以获得可靠的结果，而不依赖于这些事情。

玩得开心; o）