www.gusucode.com > matlab 案例源码 matlab代码程序 > matlab/ConstrainedDelaunayTriangulationExample.m

    %% Constrained Delaunay Triangulation
%

% Copyright 2015 The MathWorks, Inc.


%% Section 1
% Triangulate a set of points with an edge constraint specified between
% vertex |V1| and |V3|.

%%
% Define the point set.
P = [2 4; 6 1; 9 4; 6 7];

%%
% Define a constraint, |C|, between |V1| and |V3|.
C = [1 3];
DT = delaunayTriangulation(P,C);

%% 
% Plot the triangulation and add annotations.
triplot(DT)

% Label the vertices.
hold on
numvx = size(P,1);
vxlabels = arrayfun(@(n) {sprintf('V%d', n)}, (1:numvx)');
Hpl = text(P(:,1)+0.2, P(:,2)+0.2, vxlabels, 'FontWeight', ...
  'bold', 'HorizontalAlignment','center', 'BackgroundColor', ...
  'none');
hold off


% Use the incenters to find the positions for placing triangle labels on the plot.
hold on
IC = incenter(DT);
numtri = size(DT,1);
trilabels = arrayfun(@(P) {sprintf('T%d', P)}, (1:numtri)');
Htl = text(IC(:,1),IC(:,2),trilabels,'FontWeight','bold', ...
'HorizontalAlignment','center','Color','blue');
hold off

% Plot the circumcircle associated with the triangle, T1.
hold on
[CC,r] = circumcenter(DT);
theta = 0:pi/50:2*pi;
xunit = r(1)*cos(theta) + CC(1,1);
yunit = r(1)*sin(theta) + CC(1,2);
plot(xunit,yunit,'g');
axis equal
hold off


%%
% The constraint between vertices (|V1|, |V3|)
% was honored, however, the Delaunay criterion was invalidated. This
% also invalidates the nearest-neighbor relation that is inherent in
% a Delaunay triangulation. This means the |nearestNeighbor| search
% method provided by |delaunayTriangulation| cannot be supported
% if the triangulation has constraints.

%%
% In typical applications, the triangulation might be composed of many
% points, and a relatively small number of edges in the triangulation might
% be constrained. Such a triangulation is said to be locally non-Delaunay,
% because many triangles in the triangulation might respect the Delaunay
% criterion, but locally there might be some triangles that do not. In many
% applications, local relaxation of the empty circumcircle property is not
% a concern.

%%
% Constrained triangulations are generally used to triangulate a nonconvex
% polygon. The constraints give us a correspondence between the polygon
% edges and the triangulation edges. This relationship enables you to
% extract a triangulation that represents the region. The following example
% shows how to use a constrained |delaunayTriangulation| to triangulate a
% nonconvex polygon.

%%
% Define and plot a polygon.
figure()
axis([-1 17 -1 6]);
axis equal
P = [0 0; 16 0; 16 2; 2 2; 2 3; 8 3; 8 5; 0 5];
patch(P(:,1),P(:,2),'-r','LineWidth',2,'FaceColor',...
     'none','EdgeColor','r');
 
% Label the points.
hold on
numvx = size(P,1);
vxlabels = arrayfun(@(n) {sprintf('P%d', n)}, (1:numvx)');
Hpl = text(P(:,1)+0.2, P(:,2)+0.2, vxlabels, 'FontWeight', ...
  'bold', 'HorizontalAlignment','center', 'BackgroundColor', ...
  'none');
hold off

%%
% Create and plot the triangulation together with the polygon boundary.
figure()
subplot(2,1,1);
axis([-1 17 -1 6]);
axis equal

P = [0 0; 16 0; 16 2; 2 2; 2 3; 8 3; 8 5; 0 5];
DT = delaunayTriangulation(P);
triplot(DT)
hold on; 
patch(P(:,1),P(:,2),'-r','LineWidth',2,'FaceColor',...
     'none','EdgeColor','r');
hold off

% Plot the standalone triangulation in a subplot.
subplot(2,1,2);
axis([-1 17 -1 6]);
axis equal
triplot(DT)

%%
% This triangulation cannot be used to represent the domain of the polygon
% because some triangles cut across the boundary. You need to impose a
% constraint on the edges that are cut by triangulation edges. Since all
% edges have to be respected, you need to constrain all edges. The steps
% below show how to constrain all the edges.

%%
% Enter the constrained edge definition. Observe from the annotated figure
% where you need constraints (between (|V1|, |V2|), (|V2|, |V3|), and so
% on).
C = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 1];

%%
% In general, if you have |N| points in a sequence that define a polygonal
% boundary, the constraints can be expressed as |C = [(1:(N-1))' (2:N)'; N
% 1];|.

%%
% Specify the constraints when you create the |delaunayTriangulation|.
DT = delaunayTriangulation(P,C);

%%
% Alternatively, you can impose constraints on an existing triangulation by
% setting the |Constraints| property: |DT.Constraints = C;|.

%%
% Plot the triangulation and polygon.
figure('Color','white')
subplot(2,1,1);
axis([-1 17 -1 6]);
axis equal
triplot(DT)
hold on; 
patch(P(:,1),P(:,2),'-r','LineWidth',2, ...
     'FaceColor','none','EdgeColor','r'); 
hold off

% Plot the standalone triangulation in a subplot.
subplot(2,1,2);
axis([-1 17 -1 6]);
axis equal
triplot(DT)

%%
% The plot shows that the edges of the triangulation respect the boundary
% of the polygon. However, the triangulation fills the concavities. What is
% needed is a triangulation that represents the polygonal domain. You can
% extract the triangles within the polygon using the
% |delaunayTriangulation| method, |isInterior|. This method returns a
% logical array whose |true| and |false| values that indicate whether the
% triangles are inside a bounded geometric domain. The analysis is based on
% the Jordan Curve theorem, and the boundaries are defined by the edge
% constraints. The ith triangle in the triangulation is considered to be
% inside the domain if the ith logical flag is true, otherwise it is
% outside.

%%
% Now use the |isInterior| method to compute and plot the set of
% domain triangles.

% Plot the constrained edges in red
figure('Color','white')
subplot(2,1,1);
plot(P(C'),P(C'+size(P,1)),'-r','LineWidth', 2);
axis([-1 17 -1 6]);

% Compute the in/out status
IO = isInterior(DT);
subplot(2,1,2);
hold on;
axis([-1 17 -1 6]);

% Use triplot to plot the triangles that are inside.
% Uses logical indexing and dt(i,j) shorthand
% format to access the triangulation.
triplot(DT(IO, :),DT.Points(:,1), DT.Points(:,2),'LineWidth', 2)
hold off;