들로네 삼각분할(Delaunay Triangulation) 생성과 편집
이 예제에서는 delaunayTriangulation
클래스를 사용하여 들로네 삼각분할(Delaunay Triangulation)을 생성하고, 편집하고, 쿼리하는 방법을 보여줍니다. 들로네 삼각분할은 계산과학에서 가장 널리 사용되는 삼각분할 방법입니다. 이 삼각분할과 관련된 속성은 다양한 기하학적 문제를 해결하는 데 활용할 수 있는 토대를 제공합니다. 여기에는 중앙 축(Medial Axis) 계산 및 메시 모핑과 관련된 응용 예와 함께, 제약 조건이 적용되는 들로네 삼각분할(Constrained Delaunay Triangulation)을 생성하는 방법도 나와 있습니다.
예제 1: 2차원 들로네 삼각분할을 만들고 플로팅하기
이 예제에서는 2차원 들로네 삼각분할을 계산한 다음 꼭짓점 레이블 및 삼각형 레이블과 함께 삼각분할을 플로팅하는 방법을 보여줍니다.
rng default
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)
꼭짓점 레이블과 삼각형 레이블을 플롯에 표시합니다.
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차원 들로네 삼각분할을 만들고 플로팅하기
이 예제에서는 3차원 들로네 삼각분할을 계산하고 플로팅하는 방법을 보여줍니다.
rng default
X = rand(10,3);
dt = delaunayTriangulation(X)
dt = delaunayTriangulation with properties: Points: [10x3 double] ConnectivityList: [18x4 double] Constraints: []
tetramesh(dt) view([10 20])
큰 사면체 메시를 표시하려면 convexHull
메서드를 사용하여 경계 삼각분할을 계산하고 trisurf
를 사용하여 플로팅하십시오. 예를 들면 다음과 같습니다.
triboundary = convexHull(dt); trisurf(triboundary, X(:,1), X(:,2), X(:,3),'FaceColor','cyan')
예제 3: 삼각분할 데이터 구조체에 액세스하기
삼각분할 데이터 구조체에 액세스하는 방법에는 두 가지가 있습니다. 하나는 Triangulation
속성을 통해 액세스하는 방법이고, 다른 하나는 인덱싱을 사용하여 액세스하는 방법입니다.
임의의 점 10개를 사용하여 2차원 들로네 삼각분할을 만듭니다.
rng default
X = rand(10,2);
dt = delaunayTriangulation(X)
dt = delaunayTriangulation with properties: Points: [10x2 double] ConnectivityList: [11x3 double] Constraints: []
삼각분할 데이터 구조체에 액세스하는 한 가지 방법은 ConnectivityList
속성을 사용하는 것입니다.
dt.ConnectivityList
ans = 11×3
2 8 5
7 6 1
3 7 8
8 7 5
3 8 2
6 7 3
7 4 5
5 9 2
2 9 10
5 4 9
⋮
인덱싱은 삼각분할을 쿼리하기 위한 간단한 방법입니다. 구문은 dt(i,j)
이며, 여기서 j
는 i
번째 삼각형의 j
번째 꼭짓점입니다. 표준 인덱싱 규칙이 적용됩니다.
인덱싱을 사용하여 삼각분할 데이터 구조체를 쿼리합니다.
dt(:,:)
ans = 11×3
2 8 5
7 6 1
3 7 8
8 7 5
3 8 2
6 7 3
7 4 5
5 9 2
2 9 10
5 4 9
⋮
두 번째 삼각형은 다음과 같습니다.
dt(2,:)
ans = 1×3
7 6 1
두 번째 삼각형의 세 번째 꼭짓점은 다음과 같습니다.
dt(2,3)
ans = 1
처음 세 개의 삼각형은 다음과 같습니다.
dt(1:3,:)
ans = 3×3
2 8 5
7 6 1
3 7 8
예제 4: 들로네 삼각분할을 편집하여 점을 삽입하거나 제거하기
이 예제에서는 인덱스 기반 첨자를 사용하여 점을 삽입하거나 제거하는 방법을 보여줍니다. 새 delaunayTriangulation
을 처음부터 다시 만드는 것보다 delaunayTriangulation
을 편집하여 약간만 수정하는 것이 더 효율적이며, 특히 데이터 세트의 규모가 큰 경우 더욱 효율적입니다.
단위 정사각형 내에 있는 임의의 점 10개를 사용하여 들로네 삼각분할을 생성합니다.
rng default
x = rand(10,1);
y = rand(10,1);
dt = delaunayTriangulation(x,y)
dt = delaunayTriangulation with properties: Points: [10x2 double] ConnectivityList: [11x3 double] Constraints: []
임의의 점 5개를 추가로 삽입합니다.
dt.Points(end+(1:5),:) = rand(5,2)
dt = delaunayTriangulation with properties: Points: [15x2 double] ConnectivityList: [20x3 double] Constraints: []
다섯 번째 점을 수정합니다.
dt.Points(5,:) = [0 0]
dt = delaunayTriangulation with properties: Points: [15x2 double] ConnectivityList: [20x3 double] Constraints: []
네 번째 점을 제거합니다.
dt.Points(4,:) = []
dt = delaunayTriangulation with properties: Points: [14x2 double] ConnectivityList: [18x3 double] Constraints: []
예제 5: 제약 조건이 적용되는 들로네 삼각분할 생성하기
이 예제에서는 제약 조건이 적용되는 들로네 삼각분할을 만드는 방법과 제약 조건이 미치는 영향을 보여줍니다.
제약 조건이 적용되는 들로네 삼각분할을 생성하고 플로팅합니다.
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')
제약 조건이 적용되는 모서리를 빨간색으로 플로팅합니다.
hold on plot(X(C'),X(C'+size(X,1)),'-r','LineWidth',2) hold off
이제 제약 조건을 삭제하고 제약 조건이 적용되지 않는 들로네 삼각분할을 플로팅합니다.
dt.Constraints = []; subplot(2,1,2) triplot(dt) axis([-1 17 -1 6]) xlabel('Unconstrained Delaunay triangulation','FontWeight','b')
예제 6: 지형도에서 제약 조건이 적용되는 들로네 삼각분할 생성하기
국경이 표시된 미국 지도를 불러옵니다. 다각형을 나타내고 제약 조건이 적용되는 들로네 삼각분할을 생성합니다. 이 삼각분할은 점 집합에 대한 볼록 껍질(Convex Hull)에 의해 경계가 지정된 영역에 걸쳐 있습니다. 다각형의 영역 내에 있는 삼각형을 필터링한 후 플로팅합니다. 참고: 데이터 세트에는 중복된 데이터 점이 포함되어 있습니다. 즉, 위치가 동일한 데이터 점이 두 개 이상입니다. 중복된 점은 거부되고 delaunayTriangulation
은 이에 따라 제약 조건의 형식을 다시 지정합니다.
clf
load usapolygon
다각형 경계를 구성하는 두 개의 연속된 점 사이의 모서리 제약 조건을 정의하고 들로네 삼각분할을 생성합니다.
nump = numel(uslon); C = [(1:(nump-1))' (2:nump)'; nump 1]; dt = delaunayTriangulation(uslon,uslat,C);
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.
Warning: Intersecting edge constraints have been split, this may have added new points into the triangulation.
io = isInterior(dt); patch('Faces',dt(io,:),'Vertices',dt.Points,'FaceColor','r') axis equal axis([-130 -60 20 55]) xlabel('Constrained Delaunay Triangulation of usapolygon','FontWeight','b')
예제 7: 포인트 클라우드에서 곡선 재생성
이 예제에서는 들로네 삼각분할을 사용하여 포인트 클라우드에서 다각형 경계를 다시 생성하는 방법을 중점적으로 설명합니다. 재생성은 엘레강트 크러스트 알고리즘(Elegant Crust Algorithm)을 기반으로 수행됩니다.
참고 문헌: 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.
포인트 클라우드를 나타내는 점 집합을 만듭니다.
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);
점 집합의 들로네 삼각분할을 생성합니다.
dt = delaunayTriangulation(x,y); tri = dt(:,:);
기존 삼각분할에 보로노이 꼭짓점의 위치를 삽입합니다.
V = voronoiDiagram(dt);
무한 꼭짓점을 제거하고 unique
를 사용하여 중복된 점을 필터링합니다.
V(1,:) = [];
numv = size(V,1);
dt.Points(end+(1:numv),:) = unique(V,'rows');
샘플 점 쌍을 연결하는 들로네 모서리가 경계를 나타냅니다.
delEdges = edges(dt); 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 point cloud','FontWeight','b') hold off
예제 8: 다각형 영역의 근사 중앙 축(Medial Axis) 계산하기
이 예제에서는 제약 조건이 적용되는 들로네 삼각분할을 사용하여 다각형 영역의 근사 중앙 축(Medial Axis)을 만드는 방법을 보여줍니다. 다각형의 중앙 축은 다각형 내부에 있는 최대 원반의 중심의 궤적으로 정의됩니다.
영역 경계에 있는 샘플 점에 대해 제약 조건이 적용되는 들로네 삼각분할을 생성합니다.
load trimesh2d
dt = delaunayTriangulation(x,y,Constraints);
inside = isInterior(dt);
영역 삼각형을 나타내는 삼각분할을 생성합니다.
tr = triangulation(dt(inside,:),dt.Points);
인접한 삼각형의 외심을 잇는 모서리 집합을 생성합니다. 추가 로직으로 이러한 모서리의 고유한 집합을 생성합니다.
numt = size(tr,1); T = (1:numt)'; neigh = neighbors(tr); cc = circumcenter(tr); 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)]';
영역 삼각형은 녹색으로, 영역 경계는 파란색으로, 중앙 축은 빨간색으로 플로팅합니다.
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 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단계: 배경 삼각분할, 즉 메시 경계를 나타내는 점 집합에 대해 제약 조건이 적용되는 들로네 삼각분할을 생성합니다. 메시의 각 꼭짓점에 대해 배경 삼각분할과 관련하여 꼭짓점의 위치를 정의하는 descriptor를 계산합니다. 이 descriptor는 해당 삼각형에 대한 무게중심 좌표를 함께 갖고 있는 바깥쪽 삼각형입니다.
dt = delaunayTriangulation(x,y,Constraints); clf triplot(dt) axis equal axis([-10 310 -10 310]) axis equal xlabel('Background Triangulation','FontWeight','b')
descriptors.tri = pointLocation(dt,xfe,yfe); descriptors.baryCoords = cartesianToBarycentric(dt,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단계: 변형된 배경 삼각분할을 계산의 기반으로 사용하여 descriptor를 다시 카테시안 좌표(Cartesian Coordinate)로 변환합니다.
Xnew = barycentricToCartesian(tr,descriptors.tri,descriptors.baryCoords); tr = triangulation(trife,Xnew); clf triplot(tr) axis([-10 310 -10 310]) axis equal xlabel('Morphed Mesh','FontWeight','b')