-
Notifications
You must be signed in to change notification settings - Fork 328
/
perform_delaunay_flipping.m
95 lines (83 loc) · 2.36 KB
/
perform_delaunay_flipping.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
function [face, flips, flipsinv] = perform_delaunay_flipping(vertex,face,options)
% perform_delaunay_flipping - compute Dalaunay triangulation via flipping
%
% [face1, flips, flipsinv] = perform_delaunay_flipping(vertex,face,options);
%
% Set options.display_flips = 1 for graphical display.
%
% face is turned into a Delaunay triangulation face1 of the points by flipping
% the edges.
%
% 'flips' list the flips needed to go from face to face1.
% 'flipsinv' list the flips needed to go from face1 to face.
%
% Copyright (c) 2008 Gabriel Peyre
options.null = 0;
display_flips = getoptions(options, 'display_flips', 0);
verb = getoptions(options, 'verb', 1);
% edge topology
edge = compute_edges(face);
e2f = compute_edge_face_ring(face);
% initialize the stack
estack = find( check_incircle_edge(vertex,face,edge)==0 );
A = sort( compute_triangulation_angles(vertex,face) );
count = 0;
flips = []; flipsinv = [];
ntot = length(estack);
vbar = 0;
while not(isempty(estack))
if verb
vbar = max(vbar, ntot - length(estack) );
progressbar( max(vbar,1), ntot);
end
count = count+1;
if count>Inf % ntot
warning('Problem with Delaunay flipping');
break;
end
% pop non valid edge
k = estack(1); estack(1) = [];
i = edge(1,k); j = edge(2,k);
f1 = e2f(i,j); f2 = e2f(j,i);
%%% try to flip it %%%
% find the two other points
k1 = find_other( face(:,f1), i,j );
k2 = find_other( face(:,f2), i,j );
% new faces
t1 = [k1;k2;i];
t2 = [k1;k2;j];
flips(:,end+1) = [i;j];
flipsinv(:,end+1) = [k1;k2];
face(:,f1) = t1;
face(:,f2) = t2;
% update e2f
edge(:,k) = [k1;k2];
e2f = compute_edge_face_ring(face);
% display
if display_flips
clf; plot_graph(triangulation2adjacency(face),vertex);
drawnow; % pause(.01);
end
Anew = sort( compute_triangulation_angles(vertex,face) );
if lexicmp(Anew,A)<=0
warning('Problem with flipping.');
end
A = Anew;
% update the stack
estack = find( check_incircle_edge(vertex,face,edge)==0 );
estack = [ estack(estack>k); estack(estack<=k) ];
end
flipsinv = flipsinv(:,end:-1:1);
progressbar(ntot, ntot);
%%
function f = find_other( f, i,j )
I = find(f==i);
if length(I)~=1
error('Problem');
end
f(I) = [];
I = find(f==j);
if length(I)~=1
error('Problem');
end
f(I) = [];