CV筆記（二）——Delaunay三角剖分的Java實作

版權聲明：本文為原創文章，未經部落客允許不得用于商業用途。

Bowyer-Watson

算法原理

Bowyer-Watson算法是逐點插入算法，每次插入後分割插入點所在的三角形，如下圖：

之後檢查相鄰三角形是否滿足Delaunay三角網格條件，如果不滿足則通過交換四邊形對角線完成插入，如下圖：

為了保持算法的連貫性，在初始時使用一個包含所有點空間的矩形和一條對角線作為初始三角網格，并在最後删除。

實作思路

這裡使用Java實作了此算法。

由于三角形最多隻有三個相鄰三角形（共用一條邊），是以可以将三角網格儲存為連結清單結構。并且三角形需要具有外接圓屬性，是以設計如下資料結構：

class mTriangle{
    dPoint a;
    dPoint b;
    dPoint c;
    dPoint center;
    double Rsq;
    /** share points a&b*/
    int neighbor1;
    /** share points b&c*/
    int neighbor2;
    /** share points c&a*/
    int neighbor3;
    public mTriangle(dPoint p1, dPoint p2, dPoint p3){
        this.a = p1;
        this.b = p2;
        this.c = p3;
        this.neighbor1 = -1;
        this.neighbor2 = -1;
        this.neighbor3 = -1;
        center = new dPoint();
        getCenter();
    }
    /* calculate External circle center*/
    protected void getCenter(){
        double A1,A2,B1,B2,C1,C2;
        /* l1(b,c):A1x+B1y+C1*/
        A1 = c.x-b.x;
        B1 = c.y-b.y;
        C1 = (-(b.y+c.y)*B1-A1*(b.x+c.x))/2;
        /* l2(a,c):A2x+B2y+C2*/
        A2 = c.x-a.x;
        B2 = c.y-a.y;
        C2 = (-(a.y+c.y)*B2-A2*(a.x+c.x))/2;
        if(A1!=0 && A2!=0 && B1!=0 && B2!=0)
        {
            if(A1*B2==B1*A2)
            {
                System.out.println("Triangle error: edge parallel");
                return;
            }
            center.y = (C1*A2-C2*A1)/(A1*B2-B1*A2);
            center.x = -(B1*center.y+C1)/A1;
            Rsq = dPoint.sdist(center,c);
            return;
        }
        else
        {
            if((A1==0 && B1==0) || (A2==0 && B2==0))
            {
                System.out.println("Triangle error: same two point");
                return;
            }
            if((A1==0 && A2==0) || (B1==0 && B2==0))
            {
                System.out.println("Triangle error: three points collinear");
                return;
            }
            if(A1==0)
            {
                center.y = -C1/B1;
                center.x = -(B2*center.y+C2)/A2;
                Rsq = dPoint.sdist(center,c);
                return;
            }
            if(A2==0)
            {
                center.y = -C2/B2;
                center.x = -(B1*center.y+C1)/A1;
                Rsq = dPoint.sdist(center,c);
                return;
            }
            if(B1==0)
            {
                center.x = -C1/A1;
                center.y = -(C2+A2*center.x)/B2;
                Rsq = dPoint.sdist(center,c);
                return;
            }
            center.x = -C2/A2;
            center.y = -(C1+A1*center.x)/B1;
            Rsq = dPoint.sdist(center,c);
            return;
        }
    }

    /* judge if p is in triangle */
    public boolean inTriangle(dPoint p)
    {
        /* cross prod*/
        dPoint MA = new dPoint(p.x - a.x,p.y - a.y);
        dPoint MB = new dPoint(p.x - b.x,p.y - b.y);
        dPoint MC = new dPoint(p.x - c.x,p.y - c.y);
        double cp1 = MA.x * MB.y - MA.y * MB.x;
        double cp2 = MB.x * MC.y - MB.y * MC.x;
        double cp3 = MC.x * MA.y - MC.y * MA.x;
        if(cp1>=0 && cp2>=0 && cp3>=0)
            return true;
        if(cp1<=0 && cp2<=0 && cp3<=0)
            return true;
        return false;
    }

    public boolean inExternalCircle(dPoint p)
    {
        if(dPoint.sdist(p,center)<=Rsq)
            return true;
        return false;
    }
}
           

由于Bowyer-Watson算法是逐點插入的，是以實作如下插入函數，比較複雜的是如何更新三角連結清單，需要做多重判斷：

protected void addPoint(dPoint p) {
    int Tidx=0;
    /* find which triangle p is in */
    for(;Tidx<triangles.size();Tidx++)
        if(triangles.get(Tidx).inTriangle(p))
            break;
    if(Tidx>=triangles.size())
    {
        System.out.println("addPoint error: point out of super triangles");
        return;
    }
    mTriangle oldt = triangles.get(Tidx);
    mTriangle nt1,nt2,nt3;
    /* replace the old triangle by three new triangle*/
    nt1=new mTriangle(oldt.a,oldt.b,p);
    nt1.neighbor1 = oldt.neighbor1;
    nt1.neighbor2 = triangles.size();
    nt1.neighbor3 = triangles.size()+1;
    // this one has no need to update
    // updateNeighbor(oldt.neighbor1, Tidx, Tidx);

    nt2=new mTriangle(oldt.b,oldt.c,p);
    nt2.neighbor1 = oldt.neighbor2;
    nt2.neighbor2 = triangles.size()+1;
    nt2.neighbor3 = Tidx;
    updateNeighbor(oldt.neighbor2, Tidx, triangles.size());

    nt3=new mTriangle(oldt.c,oldt.a,p);
    nt3.neighbor1 = oldt.neighbor2;
    nt3.neighbor2 = Tidx;
    nt3.neighbor3 = triangles.size();
    updateNeighbor(oldt.neighbor3, Tidx, triangles.size()+1);

    triangles.set(Tidx,nt1);
    triangles.add(nt2);
    triangles.add(nt3);

    if(needUpdate(Tidx,nt1.neighbor1))
        updateQuadrilateral(Tidx,nt1.neighbor1);
    else if(needUpdate(triangles.size()-2,nt2.neighbor1))
        updateQuadrilateral(triangles.size()-2,nt2.neighbor1);
    else if(needUpdate(triangles.size()-1,nt3.neighbor1))
        updateQuadrilateral(triangles.size()-1,nt3.neighbor1);

}
           

在完成算有點插入後需要删除和初始化使用的矩形有關的三角網格

protected void deleteSuperTriangle() {
        for(int i=0,len=triangles.size();i<len;i++){
            if(useSuperTrianglePoints(i))
            {
                triangles.remove(i);
                --len;
                --i;
            }
        }
    }
           

是以算法的完整流程為：

mTriangle t1 = new mTriangle(superp1,superp2,superp4);
mTriangle t2 = new mTriangle(superp2,superp3,superp4);
t1.neighbor2=1;
t2.neighbor3=0;
triangles.add(t1);
triangles.add(t2);
for(dPoint x:pointlist)
    addPoint(x);
deleteSuperTriangle();
           

效果示例

以添加四個點為例，下圖為網格更新的流程

完整代碼見github