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）