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;