在〈Delaunay 三角分割(二)〉中認識了 Delaunay 三角分割的方式,接著就來實作吧!首先必須有個超級三角形涵蓋全部的點,由於畫布是矩形,雖然可以建立一個涵蓋畫布的超級三角形,不過建立一個比畫布大上許多的正方形,正方形中心為畫布中心,然後依對角線分為兩個三角形,會是更方便的方式,選擇正方形的原因在於,依對角線分出的兩個三角形,外接圓不會涵蓋另一個三角形的頂點:
這邊採取的頂點順序是逆時針,由於〈Delaunay 三角分割(二)〉中談到的演算過程中,必須找出不合格三角形的邊,由於三角形彼此之間會鄰接,逐一找出不合格三角形,然後去除相同邊雖然也是一個方式,然而有更快的方式。
以上圖為例,頂點順序並非隨意定義,除了是逆時針之外,頂點 0 出發的三角形,其對面三角形是特意安排為頂點 2 出發的三角形,也就是頂點 0 與頂點 2 面對的正好都是共用邊,這會有助於後續快速尋找共用邊。
畫布可能不是定義為正方形,因此在建立正方形區域時,是取畫布寬、高大者作為正方形寬度,根據這些資訊,可以初步實作如下:
class Delaunay {
constructor(width, height) {
const center = createVector(width, height).mult(0.5);
const halfW = max(width, height) * 100;
// 建立一個比畫布大上許多的正方形
this.coords = [
p5.Vector.add(center, createVector(-halfW, -halfW)),
p5.Vector.add(center, createVector(-halfW, halfW)),
p5.Vector.add(center, createVector( halfW, halfW)),
p5.Vector.add(center, createVector( halfW, -halfW))
];
// 將正方形劃為兩個三角形
const t1 = [0, 1, 3];
const t2 = [2, 3, 1];
// [三角形頂點索引] => [頂點 0 面對的三角形, 頂點 1 面對的三角形, 頂點 2 面對的三角形](也就是這個清單包含鄰接三角形)
this.triangles = new Map();
// [三角形頂點索引] => 外接圓
this.circles = new Map();
// 設定初始的兩個三角形
// t1 頂點 0 面對 t2,另兩個頂點沒有面對的三角形
this.triangles.set(t1, [t2, null, null]);
// t2 頂點 0 面對 t1,另兩個頂點沒有面對的三角形
this.triangles.set(t2, [t1, null, null]);
// 設定初始的兩個外接圓
this.circles.set(t1, circumcircle(t1.map(i => this.coords[i])));
this.circles.set(t2, circumcircle(t2.map(i => this.coords[i])));
}
...
}
為了快速取得三角形與鄰接三角形的關係,以及三角形的外接圓,這邊使用了 Map
,circumcircle
則是〈Delaunay 三角分割(一)〉中已經準備好的外接圓函式。
現在可以來加入一個新的點,並判斷哪些三角形不合格了:
class Delaunay {
...
addPoint(p) {
// 新頂點索引,先保留,後續要建立新三角形時使用
const idx = this.coords.length;
// 加入新頂點
this.coords.push(p);
// 既有的三角形外接圓若包含 p,收集在 badTriangles
const badTriangles = delaunayBadTriangles(this, p);
...
}
...
}
function delaunayBadTriangles(d, p) {
return Array.from(d.triangles.keys())
.filter(tri => inCircumcircle(tri, p, d.circles));
}
function inCircumcircle(triangle, p, circles) {
const c = circles.get(triangle);
// 以半徑平方比較
const v = p5.Vector.sub(c.center, p);
const rr = v.x * v.x + v.y * v.y;
return rr <= c.rr;
}
判斷不合格三角形時,從 circles
取得三角形的外接圓,為了避免開根號的誤差,這邊採用半徑平方比較,若不在意誤差,想直接以半徑比較也是可以。
接著下一步是,把不合格三角形的邊找出來,然而不能包含這些三角形的共用邊,因為〈Delaunay 三角分割(二)〉談過,這些共用邊若連接成新三角形,其外接圓會涵蓋其他點。
若 delaunayBoundary
可找出不合格三角形的邊,以清單傳回:
const boundary = delaunayBoundary(this, badTriangles);
delaunayBoundary
的實作如下,這邊純綷就是快速尋找鄰接三角形與邊的技巧,直接把重點說明寫在註解中了:
function delaunayBoundary(d, badTriangles) {
const boundary = [];
// 從任一不合格三角形開始尋找邊,這邊從 0 開始
let t = badTriangles[0];
// vi 是用來走訪鄰接三角形的索引
let vi = 0;
while(true) {
// 取得不合格三角形,第 vi 頂點面對的三角形
const opTri = d.triangles.get(t)[vi];
// 如果不是不合格三角形
if(badTriangles.find(tri => tri === opTri) === undefined) {
boundary.push({
// 記錄邊索引,這邊有處理循環與負索引
edge: [t[(vi + 1) % 3], t[vi > 0 ? vi - 1 : (t.length + vi - 1)]],
// 記錄第 vi 頂點面對的三角形(目前是合格的 delaunay 三角形)
delaunayTri: opTri
});
// 下個頂點索引
vi = (vi + 1) % 3;
// 邊頂點索引有相接了,表示繞行不合格的三角形們一圈了
if(boundary[0].edge[0] === boundary[boundary.length - 1].edge[1]) {
break;
}
}
// 如果 opTri 也是不合格三角形,不收集邊
else {
// 共用邊面對的 opTri 頂點
const i = d.triangles.get(opTri).findIndex(tri => tri === t);
// 下個頂點索引
vi = (i + 1) % 3;
// opTri 也是不合格三角形,用它繼續尋找邊
t = opTri;
}
}
return boundary;
}
這個繞行不合格三角形蠻有技巧性的,原因之一在於新加入的點,造成的不合格三角形們,彼此一定鄰接,原因之二在於透過 trangles
本身就記錄了頂點與各自面對的三角形的這個事實,如果光看程式難以想像的話,建議自行畫圖,逐一走訪看看就會比較明白。
現在邊已經找出來了,將不合格的三角形(以及外接圓)刪除:
// 刪除不合格的三角形以及外接圓
badTriangles.forEach(tri => {
this.triangles.delete(tri);
this.circles.delete(tri);
});
然後,用收集的邊建立新三角形,這就用到了先保留的 idx
索引:
// 用收集的邊建立新三角形
const newTriangles = boundary.map(b => {
return {
t: [idx, b.edge[0], b.edge[1]], // 新三角形
edge: b.edge, // 用哪個邊建立
delaunayTri: b.delaunayTri // 共用該邊的三角形(目前是合格的 delaunay 三角形)
};
});
newTriangles
包含的物件,每個都記錄了新三角形的頂點索引,用哪個邊建立,以及共用該邊的三角形,這是為了便於調整三角形間的鄰接關係:
class Delaunay {
...
addPoint(p) {
...
// 調整三角形間的鄰接關係
adjustNeighbors(this, newTriangles);
}
...
}
function adjustNeighbors(d, newTriangles) {
for(let i = 0; i < newTriangles.length; i++) {
const {t, edge, delaunayTri} = newTriangles[i];
// 新三角形 t,先記錄它對其中一個鄰居 delaunayTri 的關係
d.triangles.set(t, [delaunayTri, null, null]);
d.circles.set(t, circumcircle(t.map(i => d.coords[i]))); // 新外接圓
// 設定 delaunayTri 與新三角形 t 的鄰居關係
if(delaunayTri !== null) {
// delaunayTri 的舊鄰居,其中包含了被刪除的不合格三角形
const neighbors = d.triangles.get(delaunayTri);
for(let i = 0; i < neighbors.length; i++) {
const neighbor = neighbors[i];
// 如果找到原本不合格三角形的邊了
if(neighbor !== null && neighbor.includes(edge[1]) && neighbor.includes(edge[0])) {
// t 是 delaunayTri 在 edge 邊上的新鄰居
neighbors[i] = t;
break;
}
}
}
}
// 設定新三角形與其他三角形的鄰接關係
for(let i = 0; i < newTriangles.length; i++) {
const t = newTriangles[i].t;
d.triangles.get(t)[1] = newTriangles[(i + 1) % newTriangles.length].t;
d.triangles.get(t)[2] = newTriangles[i > 0 ? i - 1 : newTriangles.length + i - 1].t;
}
}
這麼一來加入的點就處理完成了,可以定義兩個方法,一個可匯出三角形頂點索引,一個可匯出頂點座標:
class Delaunay {
...
// 三角形頂點索引
indicesOfTriangles() {
return Array.from(this.triangles.keys())
.filter(tri => tri[0] > 3 && tri[1] > 3 && tri[2] > 3)
.map(tri => [tri[0] - 4, tri[1] - 4, tri[2] - 4]);
}
// 三角形頂點座標
pointsOfTriangles() {
return Array.from(this.triangles.keys())
.filter(tri => tri[0] > 3 && tri[1] > 3 && tri[2] > 3)
.map(tri => [this.coords[tri[0]], this.coords[tri[1]], this.coords[tri[2]]]);
}
}
為什麼只取索引值大於 3 的呢?記得嗎?最初的兩個三角形(用了四個索引,也就是 0 到 3)是我們額外建立的,只要是邊連接至這兩個三角形的頂點,都要拆除,同樣地,最後的頂點索引要減去 4,這樣才符合客戶端加入的頂點順序。
至於匯出三角形頂點時,由於使用了內部的 coords
來取得點座標,就直接使用既有的索引值了。
底下這個範例,可以使用滑鼠點選新點,建立 Delaunay 三角分割: