Delaunay triangulation makes me lose symmetry

Question

I am using Delaunay triangulation to split polygons into triangles. I work on FEM with a large code, and one of my "checkpoints" is symmetry (if the data is symmetric, the output has to be symmetric too). However, since I have no control over the Delaunay triangulation, it makes me lose the symmetry.

I have written a small code that illustrates my problem: we consider two disjoint triangles and a big rectangle that intersects them. We want to triangulate the subtractions of those triangles with the rectangle:

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;

As you can see, the Delaunay triangulation does not have the same behaviour in both triangles, hence the loss of symmetry.

Is there a simple way to recover symmetry?

I use Matlab R2013a.

Answer 1

Seem like you are using MatLab R2013 because in my R2011b there is no delaunayTriangulation function. To be able to run your code I changed it slightly:

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');

Here is the result I see

Is it similar to yours? Could you add an image to your question describing what is wrong with the result you get?

Answer 2

It seems this is not a bug. Your result comes from your data, actually.

I played a little with your code

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

And this is the output I seee:

Circle 1: Center at (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 Circle 4: Center at (0.35000000000000003000000,0.28333333333333355000000); R=0.05270462766947290500000 Circle 5: Center at (0.35000000000000003000000,0.28333333333333333000000); R=0.05270462766947312000000

And the figure:

So circles 1 and 3 as well as circles 4 and 5 are almost the same. So the difference between yours and mine results might come even from round-off errors, because the corresponding four points are on the same circle within the float math precision. You have to redesign your points in order to get reliable result which does not depend on such things.

Have fun ;o)