ドキュメンテーション

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

ドローネ三角形分割

ドローネ三角形分割の定義

ドローネ三角形分割は、 さまざまな用途の科学技術計算で幅広く使用されています。三角形分割の計算には非常に多くのアルゴリズムがありますが、中でもドローネ三角形分割の幾何学的な性質は非常に便利です。

基本的な特性は Delaunay 基準です。2 次元の三角形分割の場合は、空の外接円と呼ぶことがあります。2 次元の点集合では、これらの点のドローネ三角形分割は各三角形の外接円が内部に他の点を含まないようにします。このプロパティは重要です。下図に示すように、T1 の外接円は空です。内部に点が含まれていません。T2 の外接円は空です。内部に点が含まれていません。この三角形分割はドローネ三角形分割です。

下の三角形は異なります。T1 の外接円は空ではありません。内部に V3 が含まれています。T2 の外接円は空ではありません。内部に V1 が含まれています。この三角形分割はドローネ三角形分割では "ありません"。

空の外接円の性質を満たす時に、小さな内部角度をもつ三角形より、大きな内部角度をもつ三角形がより多く選択されるので、ドローネ三角形は “形が整っている” と言われます。非ドローネ三角形分割の三角形は、頂点 V2 と頂点 V4 は鋭角です。エッジ {V2, V4}V1V3 が連結したエッジに置き換えられた場合には、最小角度は最大化され、ドローネ三角形分割になります。また、ドローネ三角形分割は最近傍の点を接続します。これら 2 つの、形が整った三角形および最近傍関係という特性は、実際面において重要な意味をもつことから、散布データの内挿にはドローネ三角形分割が使用されます。

Delaunay プロパティを適切に設定していても三角形分割の位相幾何は、退化した点集合が存在する場合は一意ではありません。2 次元では 4 個以上の一意でない点が同一円内にある場合、退化が現れます。たとえば正方形の頂点に対するドローネ三角形分割は一意ではありません。

ドローネ三角形分割の性質は、より高次元に拡張します。3 次元の点集合の三角形分割は四面体で構成されます。次の図は、2 つの四面体で構成された単純な 3 次元のドローネ三角形分割を示します。1 つの四面体の外接球は、空の外接球の基準を強調表示されています。

3 次元のドローネ三角形分割は、空の外接球の基準を満たす四面体を作成します。

ドローネ三角形分割の作成

MATLAB® にはドローネ三角形分割を作成する 2 つの方法があります。

関数 delaunay は 2 次元および 3 次元のドローネ三角形分割をサポートします。関数 delaunayn は、4 次元以上のドローネ三角形分割の作成をサポートします。

ヒント

6 次元より大きいドローネ三角形分割の作成は、必要となるメモリが指数関数に増加するために、通常、中程度から大きな点集合に関しては実用的ではありません。

delaunayTriangulation クラスは 2 次元および 3 次元でのドローネ三角形分割の作成をサポートしており、三角形分割に基づくアルゴリズムの開発に有用な多くのメソッドをもっています。これらのクラス メソッドは関数と似ていますが、delaunayTriangulation を使用して作成した三角形分割のみに適用できます。また、delaunayTriangulation クラスは、凸包とボロノイ線図などの関連するコンストラクターの作成をサポートしています。さらに、制約付きドローネ三角形分割の作成をサポートします。

概要

  • 関数 delaunay は基本的な三角形分割データのみが必要で、そのデータがその用途に十分対応できる場合には便利です。

  • delaunayTriangulation クラスを使用すると、三角形分割に基づくアプリケーションの開発でさらに多くの機能が利用できます。これは、三角形分割が必要で、さらに次のいずれかの操作を実行する場合に便利です。

    • クエリ点を囲む三角形または四面体の三角形分割を探索する。

    • 三角形分割を使って最近傍探索を実行する。

    • 三角形分割の位相学的隣接または幾何学的な性質にクエリを行う。

    • 三角形分割を変更して、点を挿入または削除する。

    • 三角形分割のエッジを制限する。これは制約付きドローネ三角形分割と呼ばれる。

    • 多角形を三角形分割し、オプションで領域の外側にある三角形を削除する。

    • ドローネ三角形分割を使って、凸包またはボロノイ線図を計算する。

関数 delaunay および delaunayn の使用

関数 delaunay および delaunayn は点集合を取り、行列形式で三角形分割を作成します。このデータ構造体に関する情報は、三角形分割の行列形式を参照してください。2 次元において関数 delaunay は、散布データ点の集合によって定義される表面のプロットに利用可能な三角形分割の作成によく使われます。この用途では、このアプローチは面が一価である場合にのみ使用できることに注意してください。たとえば、単一の (xy) 座標に対応する 2 つの z 値がある場合、球面をプロットするのに使用できません。簡単な例で、関数 delaunay でサンプル データセットを表す面をプロットするための使用方法を示します。

この例では、関数 delaunay を使用して seamount データセットから 2 次元ドローネ三角形分割を作成する方法を説明します。"海山 (seamount)" は水面下の山です。このデータセットは、経度 (x) と緯度 (y) による位置の集合と、それらの座標で測定された対応する海山の標高 (z) で構成されます。

seamount データセットを読み込んで、(x, y) データを散布図として表示します。

load seamount
plot(x,y,'.','markersize',12)
xlabel('Longitude'), ylabel('Latitude')
grid on

この点集合からドローネ三角形分割を構築し、triplot を使用して既存の Figure に三角形分割をプロットします。

tri = delaunay(x,y);
hold on, triplot(tri,x,y), hold off

海山の深さのデータ (z) を追加して頂点を持ち上げ、表面を作成します。新規の Figure を作成し、trimesh を使用してワイヤーフレーム モードで表面をプロットします。

figure
hidden on
trimesh(tri,x,y,z)
xlabel('Longitude'),ylabel('Latitude'),zlabel('Depth in Feet');

影付きモードで表面をプロットする場合は、trimesh の代わりに trisurf を使用します。

関数 delaunay を使用して、3 次元ドローネ三角形分割を作成することもできます。この三角形分割は四面体で構成されます。

この例では、乱数データセットの 3 次元ドローネ三角形分割の作成方法を示します。三角形分割は tetramesh を使用してプロットされ、FaceAlpha オプションがプロットに透明性を与えます。

X = gallery('uniformdata',[30 3],0);
tet = delaunay(X);
faceColor  = [0.6875 0.8750 0.8984];
tetramesh(tet,X,'FaceColor', faceColor,'FaceAlpha',0.3);

MATLAB には関数 delaunayn があり、4 次元以上のドローネ三角形分割の作成をサポートします。2 つの相補関数 tsearchndsearchn は、N 次元の三角形分割を空間的に探索するサポートも行います。三角形分割ベースの探索に関する詳細は、空間探索 を参照してください。

delaunayTriangulation クラスの使用

MATLAB でドローネ三角形分割を作成するためのもう 1 つの方法として、delaunayTriangulation クラスが提供されています。delaunaydelaunayTriangulation は同じアルゴリズムに基づいて同じ三角形分割を作成しますが、delaunayTriangulation は Delaunay ベースのアルゴリズムの開発に役立つ補完的なメソッドを提供します。これらのメソッドは、三角形分割データと共にクラスと呼ばれるコンテナーにパッケージされた関数と同様です。すべてのものをクラスにまとめて保存しておくと、使いやすさを向上させるよりまとまった設定になります。点の位置および最近傍のような三角形分割ベースの探索の性能も向上します。delaunayTriangulation は、ドローネ三角形分割のインクリメンタルな編集をサポートします。2 次元においてエッジの制約を課すこともできます。

三角形分割表現で紹介されている triangulation クラスは、2 次元と 3 次元の三角形分割に関する幾何学的なクエリと位相幾何学的なクエリをサポートしています。delaunayTriangulation は特殊な triangulation です。つまり、Delaunay 固有のクエリだけではなく、どのような triangulation クエリでも delaunayTriangulation に実行できることを意味します。より正式な MATLAB 言語の用語では、delaunayTriangulationtriangulation のサブクラスです。

この例では、delaunayTriangulation を使用して seamount データからドローネ三角形分割を作成、クエリ、および編集する方法を説明します。seamount データセットには、(x, y) の位置と、それに対応し、海山の表面を定義する標高 (z) が含まれています。

(x, y) データを読み込んで三角形分割を行います。

load seamount
DT = delaunayTriangulation(x,y)
DT = 
  delaunayTriangulation with properties:

              Points: [294x2 double]
    ConnectivityList: [566x3 double]
         Constraints: []

エッジに何も制約が課せられていないため、Constraints プロパティは空です。Points プロパティは頂点の座標を表し、ConnectivityList プロパティは三角形を表します。これら 2 つのプロパティを合わせて三角形分割の行列データを定義します。

delaunayTriangulation クラスは行列データのラッパーであり、一連の補完的なメソッドを提供します。delaunayTriangulation のプロパティには、struct のフィールドにアクセスする場合と同じ方法でアクセスします。

頂点データにアクセスします。

DT.Points;

連結性データにアクセスします。

DT.ConnectivityList;

ConnectivityList プロパティの最初の三角形にアクセスします。

DT.ConnectivityList(1,:)
ans = 1×3

   205   230   262

delaunayTriangulation を使用すると、ConnectivityList プロパティ行列に対してインデックスを簡単に指定できます。

最初の三角形にアクセスします。

DT(1,:)
ans = 1×3

   205   230   262

最初の三角形の最初の頂点を調べます。

DT(1,1)
ans = 205

三角形分割内のすべての三角形を調べます。

DT(:,:);

delaunayTriangulation の出力 DT に対するインデックス付けは、delaunay からの三角形分割配列出力に対するインデックス付けと同様に機能します。2 つの相違点は、DT では追加のメソッドの呼び出しが可能なことです (nearestNeighborpointLocation など)。

triplot を使用して delaunayTriangulation をプロットします。関数 triplotdelaunayTriangulation メソッドではありませんが、delaunayTriangulation を受け取ってプロットすることができます。

triplot(DT); 
axis equal
xlabel('Longitude'), ylabel('Latitude')
grid on

あるいは、triplot(DT(:,:), DT.Points(:,1), DT.Points(:,2)); を使用して同じプロットを得ることもできます。

delaunayTriangulation メソッドの convexHull を使用して凸包を計算し、プロットに追加します。ドローネ三角形分割が既にあるため、このメソッドを使用すると、convhull を使用した完全な計算よりも効率的に凸包を派生させることができます。

hold on
k = convexHull(DT);
xHull = DT.Points(k,1);
yHull = DT.Points(k,2);
plot(xHull,yHull,'r','LineWidth',2); 
hold off

delaunayTriangulation を漸次編集して、点を追加または削除することができます。既存の三角形分割に点を追加する必要がある場合は、漸次追加するほうが、拡張された点集合の三角形分割を再度全体に行うよりも短時間で済みます。既存の点の数と比較して削除対象の点の数が少ない場合は、点を漸次削除するほうが効率的です。

三角形分割を編集して、前回の計算の凸包上の点を削除します。

figure
plot(xHull,yHull,'r','LineWidth',2); 
axis equal
xlabel('Longitude'),ylabel('Latitude')
grid on

% The convex hull topology duplicates the start and end vertex.
% Remove the duplicate entry.
k(end) = [];

% Now remove the points on the convex hull.
DT.Points(k,:) = []
DT = 
  delaunayTriangulation with properties:

              Points: [274x2 double]
    ConnectivityList: [528x3 double]
         Constraints: []

% Plot the new triangulation.
hold on
triplot(DT); 
hold off

凸包の境界のすぐ内側にあり、削除されていない頂点が 1 つあります。これが凸包の内側にあるということは、Figure のズームイン ツールを使用して確認できます。頂点ラベルをプロットしてこの頂点のインデックスを特定し、三角形分割から削除できます。あるいは、nearestNeighbor メソッドを使用すると、インデックスをより簡単に特定できます。

この点の位置は (211.6, -48.15) の近くです。nearestNeighbor メソッドを使用して最も近い頂点を見つけます。

vertexId = nearestNeighbor(DT, 211.6, -48.15)
vertexId = 50

次に、その頂点を三角形分割から削除します。

DT.Points(vertexId,:) = []
DT = 
  delaunayTriangulation with properties:

              Points: [273x2 double]
    ConnectivityList: [525x3 double]
         Constraints: []

新しい三角形分割をプロットします。

figure
plot(xHull,yHull,'r','LineWidth',2); 
axis equal
xlabel('Longitude'),ylabel('Latitude')
grid on
hold on
triplot(DT); 
hold off

既存の三角形分割に点を追加します。4 つの点を追加して、三角形分割の周囲に四角形を形成します。

Padditional = [210.9 -48.5; 211.6 -48.5; ...
     211.6 -47.9; 210.9 -47.9];
DT.Points(end+(1:4),:) = Padditional
DT = 
  delaunayTriangulation with properties:

              Points: [277x2 double]
    ConnectivityList: [548x3 double]
         Constraints: []

既存の Figure をすべて閉じます。

close all

新しい三角形分割をプロットします。

figure
plot(xHull,yHull,'r','LineWidth',2); 
axis equal
xlabel('Longitude'),ylabel('Latitude')
grid on
hold on
triplot(DT); 
hold off

三角形分割内の点を編集して、新しい位置に移動できます。追加の点集合のうち最初のもの (頂点 ID 274) を編集します。

DT.Points(274,:) = [211 -48.4];

既存の Figure をすべて閉じます。

close all

新しい三角形分割をプロットします。

figure
plot(xHull,yHull,'r','LineWidth',2); 
axis equal
xlabel('Longitude'),ylabel('Latitude')
grid on
hold on
triplot(DT); 
hold off

triangulation クラスのメソッド vertexAttachments を使用して、付随する三角形を見つけます。頂点に付随する三角形の数は可変であるため、このメソッドでは、付随する三角形の ID が cell 配列で返されます。内容を抽出するには中かっこが必要です。

attTris = vertexAttachments(DT,274);
hold on
triplot(DT(attTris{:},:),DT.Points(:,1),DT.Points(:,2),'g')
hold off

delaunayTriangulation は、3 次元空間の点の三角形分割にも使用できます。結果の三角形分割は四面体で構成されます。

この例では、delaunayTriangulation を使用して 3 次元点の三角形分割を作成し、プロットする方法を示します。

P = gallery('uniformdata',30,3,0);
DT = delaunayTriangulation(P)
DT = 
  delaunayTriangulation with properties:

              Points: [30x3 double]
    ConnectivityList: [117x4 double]
         Constraints: []

faceColor  = [0.6875 0.8750 0.8984];
tetramesh(DT,'FaceColor', faceColor,'FaceAlpha',0.3);

関数 tetramesh は、三角形分割の内部と外部の面をいずれもプロットします。大規模な 3 次元三角形分割の場合、内部の面をプロットすると不必要にリソースを使用することがあります。境界をプロットするほうが適切な場合があります。freeBoundary メソッドを使用すると、境界の三角形分割を行列形式で表現できます。その後、結果を trimesh または trisurf に渡します。

制約付きドローネ三角形分割

delaunayTriangulation クラスでは 2 次元の三角形分割のエッジに制約を付けることができます。これは、三角形分割で点のペアを選択し、これらの点を連結するようエッジを制約できるということです。1 つ以上のペアの点間にエッジを強制的にプロットするようなイメージです。以下の例は、エッジの制約が三角形分割にどのような影響を与えるかを示しています。

空の外接円の基準を満たすため、以下の三角形分割はドローネ三角形分割となっています。

頂点 V1V3 の間にエッジの制約を指定して、点集合を三角形分割します。

点集合を定義します。

P = [2 4; 6 1; 9 4; 6 7];

制約 CV1V3 の間に定義します。

C = [1 3];
DT = delaunayTriangulation(P,C);

三角形分割をプロットして注釈を追加します。

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

頂点 (V1V3) 間の制約は守られていますが、Delaunay 基準は無効になっています。これにより、ドローネ三角形分割に固有の最近傍関係も無効になります。これは、三角形分割に制約がある場合、delaunayTriangulation によって提供される nearestNeighbor 検索法はサポートできないことを意味します。

一般的な用途では、三角形分割は多数の点で構成されていることがあり、三角形分割の中の比較的少数のエッジが制約されている場合があります。そのような三角形分割を、ローカルで非 Delaunay であると言います。これは、三角形分割中の多くの三角形が Delaunay 基準に従うものの、ローカルではそれに従わない三角形が一部存在している可能性があるからです。多くの用途では、空の外接円の性質をローカルで緩和しても問題になりません。

制約付き三角形分割は一般に、非凸多角形を三角形分割する際に使用されます。制約により、多角形のエッジと三角形分割のエッジの間に対応関係が生じます。この関係を使って、領域を表す三角形分割を抽出できるようになります。次の例では、制約付きの delaunayTriangulation を使用して非凸多角形を三角形分割する方法を説明します。

多角形を定義してプロットします。

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

多角形の境界とともに三角形分割を作成してプロットします。

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)

一部の三角形が境界を横切っているため、この三角形分割を使用して多角形の領域を表すことはできません。三角形分割のエッジによって切り取られるエッジに制約を課す必要があります。すべてのエッジを尊重しなければならないため、すべてのエッジを制約する必要があります。以下の手順では、すべてのエッジを制約する方法を示します。

制約付きのエッジ定義を入力します。注釈付きの Figure で、どこに制約が必要かを観察します (V1V2 の間、V2V3 の間など)。

C = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 1];

一般に、連続する N 個の点が多角形の境界を定義している場合、制約は C = [(1:(N-1))' (2:N)'; N 1]; と表現できます。

delaunayTriangulation を作成するときに制約を指定します。

DT = delaunayTriangulation(P,C);

あるいは、Constraints プロパティを DT.Constraints = C; のように設定することにより、既存の三角形分割に制約を課すこともできます。

三角形分割と多角形をプロットします。

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)

プロットでは、三角形分割のエッジによって多角形の境界が尊重されていることが示されています。ただし、凹部は三角形分割で埋められています。必要なのは、多角形の領域を表す三角形分割です。delaunayTriangulation のメソッドである isInterior を使用すると、多角形内の三角形を抽出できます。このメソッドは logical 配列を返し、その値が truefalse かで、三角形の位置が囲まれた幾何学的領域の内側かどうかが示されます。解析はジョルダン曲線定理に基づき、境界はエッジの制約によって定義されます。三角形分割の i 番目の三角形は、i 番目の論理フラグが true の場合に領域内にあると見なされ、それ以外の場合は外側にあると見なされます。

ここで、isInterior メソッドを使用して、領域の三角形のセットを計算し、プロットします。

% 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;

重複する位置を含む点集合の三角形分割

MATLAB の Delaunay アルゴリズムは、固有の点集合から三角形分割を作成します。三角形分割関数またはクラスに渡された点が固有なものでない場合、重複する位置は検出され、重複する点は無視されます。これにより、オリジナルの入力の一部の点、すなわち重複する点を参照しない三角形分割が作成されます。関数 delaunaydelaunayn を操作する場合、重複の存在は重要でなくなります。ただし、delaunayTriangulation クラスにより提供される多くのクエリはインデックス ベースであるため、delaunayTriangulation では固有のデータセットを使用して、三角形分割とデータの操作が行われることに注意してください。このため、通常は固有の点集合に基づいてインデックスを指定します。このデータは delaunayTriangulationPoints プロパティにより保持されます。

次の例は、delaunayTriangulation を使用する際に、Points プロパティに格納されている固有のデータセットを参照することの重要性を示しています。

P = gallery('uniformdata',[25 2],0);
P(18,:) = P(8,:)
P(16,:) = P(6,:)
P(12,:) = P(2,:)

DT = delaunayTriangulation(P)
この三角形分割が作成される際に、MATLAB では警告が表示されます。Points プロパティに、重複する点がデータから削除されていることが表示されます。
DT = delaunayTriangulation with properties: Points: [22x2 double] ConnectivityList: [31x3 double] Constraints: []
たとえば、ドローネ三角形分割が凸包を計算するために使用される場合、包の点のインデックスは固有の点集合 DT.Points に対するインデックスです。したがって、次のコードを使用して凸包の計算とプロットを行います。
K = DT.convexHull();
plot(DT.Points(:,1),DT.Points(:,2),'.');
hold on
plot(DT.Points(K,1),DT.Points(K,2),'-r');
delaunayTriangulation により提供されたインデックスと共に重複を含む元のデータセットを使用すると、結果が不正確になることがあります。delaunayTriangulation は、固有のデータセット DT.Points に基づくインデックスを使用します。たとえば、KP に関してではなく、DT.Points に関してインデックス付けされているので、以下は不正確なプロットを作成することになります。
K = DT.convexHull();
plot(P(:,1),P(:,2),'.');
hold on
plot(P(K,1),P(K,2),'-r');
通常は delaunayTriangulation を作成する前に重複を削除して、固有のデータセットを作成した方が便利です。これを行うと混乱を避けることができます。これは、関数 unique を使用して、次のように実行することができます。
P = gallery('uniformdata',[25 2],0);
P(18,:) = P(8,:)
P(16,:) = P(6,:)
P(12,:) = P(2,:)

[~, I, ~] = unique(P,'first','rows');
I = sort(I);
P = P(I,:);
DT = delaunayTriangulation(P)  % The point set is unique

関連するトピック