Delaunay 三角分割(三)


在〈Delaunay 三角分割(二)〉中認識了 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])));
    }

    ...
}

為了快速取得三角形與鄰接三角形的關係,以及三角形的外接圓,這邊使用了 Mapcircumcircle 則是〈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 三角分割: