ドキュメンテーション センター

  • 評価版
  • 製品アップデート

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

Delaunay 三角形分割の作成と編集

この例では、delaunayTriangulation クラスを使用して Delaunay 三角形分割の作成、編集、クエリを行う方法を示します。Delaunay 三角形分割は、科学計算で最も広く使われている三角形分割です。三角形分割に関するプロパティは、さまざまな幾何学的な問題を解決するための基礎を提供します。中心軸の計算やメッシュのモーフィングに関するアプリケーションを使って、制約付き Delaunay 三角形分割の作成についても示します。

例 1:2 次元 Delaunay 三角形分割の作成とプロット

この例は、2 次元 Delaunay 三角形分割を計算し、頂点と三角形のラベルと共に三角形分割をプロットする方法を示します。

x = rand(10,1);
y = rand(10,1);
dt = delaunayTriangulation(x,y)
dt = 

  delaunayTriangulation with properties:

              Points: [10x2 double]
    ConnectivityList: [11x3 double]
         Constraints: []

triplot(dt);
%
% Display the Vertex and Triangle labels on the plot
hold on
vxlabels = arrayfun(@(n) {sprintf('P%d', n)}, (1:10)');
Hpl = text(x, y, vxlabels, 'FontWeight', 'bold', 'HorizontalAlignment',...
           'center', 'BackgroundColor', 'none');
ic = incenter(dt);
numtri = size(dt,1);
trilabels = arrayfun(@(x) {sprintf('T%d', x)}, (1:numtri)');
Htl = text(ic(:,1), ic(:,2), trilabels, 'FontWeight', 'bold', ...
      'HorizontalAlignment', 'center', 'Color', 'blue');
hold off

例 2:3 次元 Delaunay 三角形分割の作成とプロット

この例は、3 次元 Delaunay 三角形分割を計算し、三角形分割をプロットする方法を示します。

X = rand(10,3)
X =

    0.6557    0.7060    0.4387
    0.0357    0.0318    0.3816
    0.8491    0.2769    0.7655
    0.9340    0.0462    0.7952
    0.6787    0.0971    0.1869
    0.7577    0.8235    0.4898
    0.7431    0.6948    0.4456
    0.3922    0.3171    0.6463
    0.6555    0.9502    0.7094
    0.1712    0.0344    0.7547

dt = delaunayTriangulation(X)
dt = 

  delaunayTriangulation with properties:

              Points: [10x3 double]
    ConnectivityList: [22x4 double]
         Constraints: []

tetramesh(dt, 'FaceColor', 'cyan');
% To display large tetrahedral meshes use the convexHull method to
% compute the boundary triangulation and plot it using trisurf.
% For example;
% triboundary = convexHull(dt)
% trisurf(triboundary, X(:,1), X(:,2), X(:,3), 'FaceColor', 'cyan')

例 3:三角形分割データの構造体にアクセス

三角形分割データの構造体にアクセスする 2 つの方法があります。1 つ目の方法は Triangulation プロパティを使用し、2 つ目の方法はインデックスを使用します。

10 個のランダム点から 2 次元 Delaunay 三角形分割を作成します。

X = rand(10,2)
X =

    0.2760    0.7513
    0.6797    0.2551
    0.6551    0.5060
    0.1626    0.6991
    0.1190    0.8909
    0.4984    0.9593
    0.9597    0.5472
    0.3404    0.1386
    0.5853    0.1493
    0.2238    0.2575

dt = delaunayTriangulation(X)
dt = 

  delaunayTriangulation with properties:

              Points: [10x2 double]
    ConnectivityList: [12x3 double]
         Constraints: []

% The triangulation datastructure is;
dt.ConnectivityList
ans =

     4    10     1
     3     8     9
    10     4     5
     4     1     5
     1     6     5
     3    10     8
     1     3     6
     2     9     7
     6     3     7
     1    10     3
     3     2     7
     3     9     2

% Indexing is a shorthand way to query the triangulation. The format is
% dt(i, j) where j is the j'th vertex of the i'th triangle, standard
% indexing rules apply.
% The triangulation datastructure is
dt(:,:)
ans =

     4    10     1
     3     8     9
    10     4     5
     4     1     5
     1     6     5
     3    10     8
     1     3     6
     2     9     7
     6     3     7
     1    10     3
     3     2     7
     3     9     2

2 番目の三角形は以下のとおりです。

dt(2,:)
ans =

     3     8     9

2 番目の三角形の 3 番目の頂点は以下のとおりです。

dt(2,3)
ans =

     9

最初の 3 つの三角形は以下のとおりです。

dt(1:3,:)
ans =

     4    10     1
     3     8     9
    10     4     5

例 4:点を挿入または削除するための Delaunay 三角形分割の編集

この例は、点を挿入または削除するために、インデックスベースの添字を使用する方法を示します。小さな変更であれば、新たな delaunayTriangulation をゼロから作成するよりも、delaunayTriangulation を編集した方がより効率的になります。これは、データセットが大きい場合に特に言えることです。

% Construct a Delaunay triangulation from
% 10 random points within a unit square
x = rand(10,1);
y = rand(10,1);
dt = delaunayTriangulation(x,y)
dt = 

  delaunayTriangulation with properties:

              Points: [10x2 double]
    ConnectivityList: [12x3 double]
         Constraints: []

% Insert 5 additional random points
dt.Points(end+(1:5),:) = rand(5,2)
dt = 

  delaunayTriangulation with properties:

              Points: [15x2 double]
    ConnectivityList: [21x3 double]
         Constraints: []

5 番目の点の置き換え

dt.Points(5,:) = [0, 0]
dt = 

  delaunayTriangulation with properties:

              Points: [15x2 double]
    ConnectivityList: [21x3 double]
         Constraints: []

4 番目の点を削除

dt.Points(4,:) = []
dt = 

  delaunayTriangulation with properties:

              Points: [14x2 double]
    ConnectivityList: [19x3 double]
         Constraints: []

例 5:制約付き Delaunay 三角形分割の作成

この例は、簡単な制約付き Delaunay 三角形分割を作成する方法を示し、制約の効果を示します。

X = [0 0; 16 0; 16 2; 2 2; 2 3; 8 3; 8 5; 0 5];
C = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 1];
dt = delaunayTriangulation(X, C);
subplot(2,1,1);
triplot(dt);
axis([-1 17 -1 6]);
xlabel('Constrained Delaunay triangulation', 'fontweight','b');
% Plot the constrained edges in red
hold on;
plot(X(C'),X(C'+size(X,1)),'-r', 'LineWidth', 2);
hold off;

% Now delete the constraints and plot the unconstrained Delaunay
dt.Constraints = [];
subplot(2,1,2);
triplot(dt);
axis([-1 17 -1 6]);
xlabel('Unconstrained Delaunay triangulation', 'fontweight','b');

例 6:地図の制約付き Delaunay 三角形分割の作成

隣接する米国の周辺部の地図を読み込みます。多角形を表す制約付き Delaunay 三角形分割を作成します。この三角形分割は、点集合の凸包で囲まれた領域にまたがります。多角形の領域内にある三角形を除去し、プロットします。メモ:データセットは重複するデータ点を含んでいます。すなわち、2 つ以上のデータ点が同じ位置にあります。重複する点を取り除いてから、delaunayTriangulation はそれに応じて制約を再度形成します。

clf
load usapolygon
% Define an edge constraint between two successive
% points that make up the polygonal boundary.
nump = numel(uslon);
C = [(1:(nump-1))' (2:nump)'; nump 1];
dt = delaunayTriangulation(uslon, uslat, C);
io = dt.isInterior();
patch('faces',dt(io,:), 'vertices', dt.Points, 'FaceColor','r');
axis equal;
axis([-130 -60 20 55]);
xlabel('Constrained Delaunay Triangulation of usapolygon', 'fontweight','b');
Warning: Duplicate data points have been detected and removed.
 The Triangulation indices and Constraints are defined with respect to the
 unique set of points in delaunayTriangulation property X. 
Warning: Intersecting edge constraints have been split, this may have added
new points into the triangulation. 

例 7:点群から曲線を再構成

この例は、点群から多角形の境界を再構成するための Delaunay 三角形分割の使用方法を示します。再構成は Crust アルゴリズムに基づきます。

参照:N. Amenta, M. Bern, and D. Eppstein. The crust and the beta-skeleton:combinatorial curve reconstruction.Graphical Models and Image Processing, 60:125-135, 1998.

% Create a set of points representing the point cloud
numpts=192;
t = linspace( -pi, pi, numpts+1 )';
t(end) = [];
r = 0.1 + 5*sqrt( cos( 6*t ).^2 + (0.7).^2 );
x = r.*cos(t);
y = r.*sin(t);
ri = randperm(numpts);
x = x(ri);
y = y(ri);
% Construct a Delaunay Triangulation of the point set.
dt = delaunayTriangulation(x,y);
tri = dt(:,:);
% Insert the location of the Voronoi vertices into the existing
% triangulation
V = dt.voronoiDiagram();
% Remove the infinite vertex
V(1,:) = [];
numv = size(V,1);
dt.Points(end+(1:numv),:) = V;
Warning: Duplicate data points have been detected and removed.
 The Triangulation indices are defined with respect to the unique set of
 points in delaunayTriangulation property X. 
% The Delaunay edges that connect pairs of sample points represent the
% boundary.
delEdges = dt.edges();
validx = delEdges(:,1) <= numpts;
validy = delEdges(:,2) <= numpts;
boundaryEdges = delEdges((validx & validy), :)';
xb = x(boundaryEdges);
yb = y(boundaryEdges);
clf;
triplot(tri,x,y);
axis equal;
hold on;
plot(x,y,'*r');
plot(xb,yb, '-r');
xlabel('Curve reconstruction from a point cloud', 'fontweight','b');
hold off;

例 8:多角形領域の近似の中心軸の計算

この例では、制約付き Delaunay 三角形分割を使って多角形領域の近似の中心軸を作成する方法を示します。多角形の中心軸は、多角形内部の最大の円盤形の中心の位置で定義されます。

% Construct a constrained Delaunay triangulation of a sample of points
% on the domain boundary.
load trimesh2d
dt = delaunayTriangulation(x,y,Constraints);
inside = dt.isInterior();
% Construct a triangulation to represent the domain triangles.
tr = triangulation(dt(inside, :), dt.Points);

% Construct a set of edges that join the circumcenters of neighboring
% triangles; the additional logic constructs a unique set of such edges.
numt = size(tr,1);
T = (1:numt)';
neigh = tr.neighbors();
cc = tr.circumcenter();
xcc = cc(:,1);
ycc = cc(:,2);
idx1 = T < neigh(:,1);
idx2 = T < neigh(:,2);
idx3 = T < neigh(:,3);
neigh = [T(idx1) neigh(idx1,1); T(idx2) neigh(idx2,2); T(idx3) neigh(idx3,3)]';
% Plot the domain triangles in green, the domain boundary in blue and the
% medial axis in red.
clf;
triplot(tr, 'g');
hold on;
plot(xcc(neigh), ycc(neigh), '-r', 'LineWidth', 1.5);
axis([-10 310 -10 310]);
axis equal;
plot(x(Constraints'),y(Constraints'), '-b', 'LineWidth', 1.5);
xlabel('Medial Axis of a Polygonal Domain', 'fontweight','b');
hold off;

例 9:2 次元メッシュを修正した境界に変形

この例は、領域の境界線への修正を調整するために 2 次元領域のメッシュを変形する方法を示します。

手順 1: データを読み込みます。変形するメッシュは、面と頂点の形式の三角形分割 trife、xfe、yfe で定義されます。

load trimesh2d
clf; triplot(trife,xfe,yfe); axis equal;
axis([-10 310 -10 310]);
axis equal;
xlabel('Initial Mesh', 'fontweight','b');

手順 2: 背景の三角形分割 - メッシュの境界を表す点集合の制約付き Delaunay 三角形分割を作成します。メッシュの頂点ごとに、背景の三角形分割に対して位置を定義する記述子を計算します。記述子はその三角形に対して重心座標と共に囲まれた三角形を表します。

dt = delaunayTriangulation(x,y,Constraints);
clf; triplot(dt); axis equal;
axis([-10 310 -10 310]);
axis equal;
xlabel('Background Triangulation', 'fontweight','b');
descriptors.tri = dt.pointLocation(xfe, yfe);
descriptors.baryCoords = dt.cartesianToBarycentric(descriptors.tri, [xfe yfe]);

手順 3: 領域の境界線を修正するために、背景の三角形分割を修正します。

cc1 = [210 90];
circ1 = (143:180)';
x(circ1) = (x(circ1)-cc1(1))*0.6 + cc1(1);
y(circ1) = (y(circ1)-cc1(2))*0.6 + cc1(2);
tr = triangulation(dt(:,:),x,y);
clf; triplot(tr); axis([-10 310 -10 310]); axis equal;
xlabel('Edited Background Triangulation - Hole Size Reduced', 'fontweight','b');

手順 4: 評価の根拠として、変形した背景の三角形分割を使用し、記述子を直交座標に戻します。

Xnew = tr.barycentricToCartesian(descriptors.tri, descriptors.baryCoords);
tr = triangulation(trife, Xnew);
clf; triplot(tr);
axis([-10 310 -10 310]);
axis equal;
xlabel('Morphed Mesh', 'fontweight','b');

この情報は役に立ちましたか?