ドロネーの三角形分割の実装に失敗した話

はじめに

この記事は何もわかっていない自分がメモとして書いています。これからも多分更新があります。ご注意ください。

以下の二つの記事の続きを書くつもりでいました。

• 目指せ幾何マスター その 1
• 目指せ幾何マスター その 2

投稿時点で卒論がヤバいので一旦中断してここにメモです。

趣旨

アルゴリズム・サイエンス シリーズの「計算幾何」(浅野哲夫 著) という本を覗くと、ドロネーの三角形分割というデータ構造が紹介されていました。三角形分割とは二次元座標系上の点群 $P$ に属する点同士を線で結んでできる平面グラフのうち、辺数が最も多いものを構築することです。要するに点同士を結んだ辺が、$P$ の凸包を $P$ に属する点をもちいて三角形の領域に分割している状態を構築することです。

ドロネーの三角形分割は、三角形分割のうち、分割した三角形のうち内角の最小値が最大である三角形分割です。

アルゴリズム概要

テキトーな三角形分割から、隣接する二つの三角形を選び、それらを ▲ABC と ▲DBC としましょう。もしこれらがドロネーの三角形分割に含まれないなら、共有する辺 BC を削除し、新たに辺 AD を加えて二つの ▲ADC と ▲ADB を作ります。これを対角変形と言います。

ドロネーの三角形に含まれるかどうかは、▲ABC と ▲DBC の内角の最小値と、▲ADC と ▲ADB の内角の最小値を比較して判断します。

本によると、テキトーな三角形をドロネーの三角形分割にするために対角変形を行う回数は $O(|P|^2)$ 回らしいです。

高速化 (できませんでした)

こちらの記事 (ドロネー三角形分割の期待最速アルゴリズム) によると高速化できるらしいです。とはいえ、ちゃんと読み込まずに何とかなる精神で実装を始めてしまったため、高速化の正当性を全く理解していません(参考文献のリンク先が非公開になっていて泣きました)。それは実装を失敗するのも納得です。。。

一応この記事で説明されている通り、$P$ の点を一つずつデータ構造に加えていく手続きを実装しました。

やること

先述の記事によると、三角形分割をめちゃんこ大きい三角形で初期化するらしいです。この三角形の頂点は $P$ に含まれるものでなくて良いらしいです。

ところが、座標のサイズがでかくなることを心配し、オリジナリティを出してみました。これがダメかと思いもしましたが、ここを変えても遅いままだったのでひとまずはオリジナリティを出しておきます。

自分の実装では $P$ の凸包を最も左の点 B を共有する形で ccw (反時計回り) 方向に三角形分割していきます。

このとき、根を点 B に、それ以外のノードを分割した三角形に対応付けた木を作ります。ドロネー木 (木ではない) とは多分これのことです。多分。ドロネー木は三角形の探索に使用します。

また、対角変形の際は以下のように、新たにできる三角形に対応するノードを生成し、交差する三角形に対応するノードとの間に辺を構築します。

$P$ の点をデータ構造に追加するときも同様に、追加した点によって新たに分割された三角形をドロネー木に反映させます。

このとき、追加した点によって分割される領域を、ドロネー木の根 rB から探索します。
ドロネー木の根から出る辺数は、$O(凸包の頂点数)$ 個なので、根の子を全部見る遅いですが、凸包の三角形分割は、点 B から ccw (反時計回り) であることを利用すると、二分探索によって $\log$ をつけることができます。

根以外の子は、対角変形か点の追加の影響しかないので、子頂点が高々 $3$ です。例えば、以下のように辺上でない場所に点を追加すると $3$ つの子頂点ができます。

また、$P$ の三角形分割は、ドロネー木の末端ノードに対応しています。

実装のポインタ

後日書きます。

ソースコード

めっちゃ遅いです。丸一日粘りましたが、卒論があるので今はこれ以上は無理です。

/*
    Copyright ©️ (c) NokonoKotlin (okoteiyu) 2025. 
    Released under the MIT license(https://opensource.org/licenses/mit-license.php) 
*/

#include <vector>
#include <unordered_map>
#include <stack>
#include <algorithm>
#include <functional>
#include <iostream>
#include <cmath>
#include <cassert>
using namespace std;


template<typename T>
struct DelaunayTriangles{
    /*
        Copyright ©️ (c) NokonoKotlin (okoteiyu) 2025. 
        Released under the MIT license(https://opensource.org/licenses/mit-license.php) 
    */
    // 二次元座標系上の点
    struct ms_Point{
        T x,y;
        int id = -1;
        ms_Point(){}
        ms_Point(const ms_Point& p){
            this->x = p.x;
            this->y = p.y;
            this->id = p.id;// id ごとコピー
        }
        ms_Point(T x_ , T y_) 
            : x(x_) , y(y_) {}
        // id を無視して一致判定 (コピーと行動しないようにに注意)
        bool operator==(ms_Point p){return bool(this->x == p.x && this->y == p.y);}
    };

    // オブジェクトなしの呼び出しを可能にするために friend の非メンバで宣言する。
    // 引数がクラス独自の構造体のため、グローバルでのバッティングも (キャストがなければ) ない。

    friend inline constexpr const T signed_area_doubled(ms_Point A , ms_Point B , ms_Point C){
        return ((B.x-A.x)*(C.y-A.y) - (B.y-A.y)*(C.x-A.x));
    }

    friend bool ccw(ms_Point Center , ms_Point a , ms_Point b , bool include = false){
        T S = signed_area_doubled(Center,a,b);
        if(include)return bool(S >= 0);
        return bool(S > 0);
    }
    friend bool same_line(ms_Point a , ms_Point b , ms_Point c){
        if(a == b || b == c || c == a)return true;
        return bool(signed_area_doubled(a,b,c) == 0);
    }

    friend vector<ms_Point> ConvexHull(vector<ms_Point> Points){
        assert(Points.size() >= 3);
        function<vector<ms_Point>(vector<ms_Point>)> GrahamScan = [&](vector<ms_Point> P){
            vector<ms_Point> CH;
            sort(P.begin(),P.end(),[&](ms_Point a ,ms_Point b){if(a.x == b.x)return a.y < b.y;return a.x < b.x;});
            int cur = 0;
            while(CH.size()<2)CH.push_back(P[cur++]);
            for( ; cur < int(P.size()) ; cur++){
                while(CH.size()>=2){
                    bool pop_ = false;
                    if(same_line(P[cur] , CH.back() , CH[int(CH.size())-2]))pop_ = true;
                    if(ccw(CH.back() , CH[int(CH.size())-2] , P[cur]))pop_ = true;
                    if(pop_)CH.pop_back();
                    else break;
                }
                CH.push_back(P[cur]);
            }
            return CH;
        };
        vector<ms_Point> convex_hull;
        convex_hull = GrahamScan(Points);
        for(ms_Point& p : Points){p.x*=-1;p.y*=-1;}
        vector<ms_Point> tmp = GrahamScan(Points);
        for(ms_Point& p : tmp){p.x*=-1;p.y*=-1;}
        convex_hull.pop_back();
        tmp.pop_back();
        for(ms_Point p : tmp)convex_hull.push_back(p);
        return convex_hull;
    }

    // 三角形
    struct ms_Triangle{
        private:
        vector<ms_Point> P;
        public:
        ms_Triangle(){}
        ms_Triangle(const ms_Point& p1, const ms_Point& p2, const ms_Point& p3){
            P.push_back(p1);
            P.push_back(p2);
            P.push_back(p3);
            P = ConvexHull(P);
        }
        bool is_include(ms_Point p) const& {
            return bool(
                    abs(signed_area_doubled(p,P[0],P[1])) 
                    + abs(signed_area_doubled(p,P[0],P[2]))
                    + abs(signed_area_doubled(p,P[1],P[2]))
                    == abs(signed_area_doubled(P[0],P[1],P[2])));
        }
        const ms_Point& operator[](int i) const& {return this->P[i];} 
    };


    // ドロネー木のノード
    struct DTreeNode{
        protected:
        int id_ = -1;
        // 3 つまで次のポインタを持つ (vector よりは流石に早そう)
        DTreeNode* child[3];
        ms_Triangle Tr;
        bool valid_ = true;
        public:

        DTreeNode(const ms_Point& p1 ,const ms_Point& p2 , const ms_Point& p3 ) : Tr(p1,p2,p3) {
            this->child[0] = nullptr;
            this->child[1] = nullptr;
            this->child[2] = nullptr;
        }

        // p で三角形分割
        void divide(ms_Point p){
            for(int c = 0 ;c < 3 ; c++)assert(this->child[c] == nullptr);
            for(int i = 0 ; i < 3 ; i++){
                if(same_line(p , this->data()[i%3] , this->data()[(i+1)%3]))continue;
                this->child[i] = new DTreeNode(p , this->data()[i%3] , this->data()[(i+1)%3]);
            }
        }

        // 対角変形でノードが追加される
        void connect(DTreeNode* nd1 , DTreeNode* nd2){
            this->valid_ = false;
            assert(this->is_end());
            this->child[0] = nd1;
            this->child[1] = nd2;
        }

        int id(){return this->id_;}
        // 自身が終端ノードかどうか
        bool is_end(){return bool(this->child[0] == nullptr && this->child[1] == nullptr && this->child[2] == nullptr);}
        // 有効かどうか (対角変形されていないか)
        bool is_valid(){
            return this->valid_;
        }
        // 三角形のデータ
        const ms_Triangle& data(){return this->Tr;}
        // 子ノードのポインタ
        DTreeNode* operator[](int i){assert(0 <= i && i < 3);return this->child[i];}
    };



    private:
    inline static const long long _hash_border_ = 2000000ll;
    // 番号が u,v　の二点をつなぐ辺 (u,v) をハッシュ化する
    inline constexpr const long long edge_to_ll(long long u , long long v){
        if(u > v)swap(u,v);// 無向辺なので順を区別しない
        return u*_hash_border_ + v;
    }

    // 番号が a,b,c の 3 点を結ぶ三角形をハッシュ化
    inline constexpr const long long triangle_to_ll(long long a , long long b , long long c){
        if(a > b)swap(a,b);// 順を区別しない
        if(a > c)swap(a,c);
        if(b > c)swap(b,c);
        return a*_hash_border_*_hash_border_ + b*_hash_border_ + c;
    }

    // [p][edge_to_ll(u,v)] := 番号が p の点からみて、辺 (u,v) の向こう側の三角形の、u,v 以外の点の番号
    vector<unordered_map<long long , int>> m_beyond_edge;

    // [triangle_to_ll(a,b,c)] := 番号 a,b,c の 3 点を結ぶ三角形に対応するドロネー木の頂点へのポインタ
    unordered_map<long long , DTreeNode*> NodeTable;

    // [x] := id が [[x+1]] の点
    vector<ms_Point> m_P;
    
    // ドロネー木の最初のレイヤー
    // 凸包を、左下の点 m_BasePoint を端点として |CH| - 2 個の三角形に分割し、CCW 順に格納
    // ここに格納されるノードの三角形のみ、頂点が左下のものから反時計回りの順番であることが保証
    vector<DTreeNode*> m_FirstLayer;

    // ドロネー木から、p を含む三角形の末端のものを探索する
    vector<DTreeNode*> search(ms_Point p){
        // 初めのレイヤーは二分探索する
        int lef = -1;// m_FirstLayer[lef] までの領域に p が含まれない
        int rig = int(m_FirstLayer.size()) - 1;// m_FirstLayer[rig] までの領域に p が含まれる
        while(rig-lef>1){
            int mid = lef + (rig-lef)/2;
            // 判別は m_BasePoint からの ccw を見て行う
            ms_Point base = m_FirstLayer[0]->data()[0];
            ms_Point midium = m_FirstLayer[mid]->data()[2];
            if(ccw(base , p , midium , true))rig = mid;
            else lef = mid;
        }
        unordered_map<DTreeNode*,bool> visited;
        vector<DTreeNode*> res;
        stack<DTreeNode*> st;// 候補を stack に積んでいく
        for(int i = rig-1 ; i <= rig+1 ; i++){
            if(i < 0 || m_FirstLayer.size() <= i)continue;
            if(m_FirstLayer[i]->data().is_include(p))st.push(m_FirstLayer[i]);
        }
        
        while(!st.empty()){
            DTreeNode* nd = st.top();
            st.pop();
            if(visited[nd])continue;
            visited[nd] = true;
            if(nd->is_end() && nd->data().is_include(p) && nd->is_valid())res.push_back(nd);
            for(int c = 0 ; c < 3 ; c++){
                if((*nd)[c] != nullptr)st.push((*nd)[c]);
            }
        }
        return res;
    }

    long double _sub_dist_(const ms_Point& A ,const ms_Point& B){
        return sqrt((long double)(B.x-A.x)*(long double)(B.x-A.x) + (long double)(B.y-A.y)*(long double)(B.y-A.y));
    }
    // 三角形 ABC の角のうち、cosine が最大のもの (最小角のcos)
    long double _sub_max_cosine_(const ms_Point& A ,const ms_Point& B , const ms_Point& C){
        long double res = -1;
        // AB & AC
        res = max(res , (long double)((B.x-A.x)*(C.x-A.x) + (B.y-A.y)*(C.y-A.y))/(_sub_dist_(A,B)*_sub_dist_(A,C)));
        // BA & BC
        res = max(res , (long double)((A.x-B.x)*(C.x-B.x) + (A.y-B.y)*(C.y-B.y))/(_sub_dist_(B,A)*_sub_dist_(B,C)));
        // CA & CB
        res = max(res , (long double)((A.x-C.x)*(B.x-C.x) + (A.y-C.y)*(B.y-C.y))/(_sub_dist_(C,A)*_sub_dist_(C,B)));
        return res;
    }

    // 番号 a の点から見て、辺 b,c を必要に応じて対角変形する。
    void diagonal_flip(int a , int b , int c){
        // abc に対応するノード
        DTreeNode* nd1 = NodeTable[triangle_to_ll(a,b,c)];
        if(nd1 == nullptr)return;
        if(!nd1->is_end())return;
        if(!nd1->is_valid())return;
        
        // id が 1 ずれることに注意
        int d = m_beyond_edge[a][edge_to_ll(b,c)];
        if(d == 0)return;// 対角が存在しないとき
        const ms_Point& A = m_P[a-1];
        const ms_Point& B = m_P[b-1];
        const ms_Point& C = m_P[c-1];
        const ms_Point& D = m_P[d-1];
        // 4 点が凸でないときは何もしない
        if(abs(signed_area_doubled(A,B,D)) + abs(signed_area_doubled(A,C,D)) 
                != abs(signed_area_doubled(A,B,C)) + abs(signed_area_doubled(D,B,C)))return;

        // 三角形をなすときも何もしない
        if(same_line(A,B,D) || same_line(A,C,D))return;
        // 最小の内角が大きくなるかどうかを cosine を見てチェック
        if(max(_sub_max_cosine_(A,B,C) , _sub_max_cosine_(D,B,C)) <= max(_sub_max_cosine_(B,D,A) , _sub_max_cosine_(C,D,A))){
            return;
        }
        
        // 対角変換周りの 2 頂点
        DTreeNode* nd2 = NodeTable[triangle_to_ll(d,b,c)];
        if(nd2 == nullptr)return;
        if(!nd2->is_end())return;
        if(!nd2->is_valid())return;

        DTreeNode* new_nd1 = new DTreeNode(A,D,B);
        DTreeNode* new_nd2 = new DTreeNode(A,D,C);

        NodeTable[triangle_to_ll(a,d,b)] = new_nd1;
        NodeTable[triangle_to_ll(a,d,c)] = new_nd2;

        nd1->connect(new_nd1,new_nd2);
        nd2->connect(new_nd1,new_nd2);

        m_beyond_edge[b][edge_to_ll(a,d)] = c;
        m_beyond_edge[c][edge_to_ll(a,d)] = b;

        m_beyond_edge[a][edge_to_ll(c,d)] = m_beyond_edge[b][edge_to_ll(c,d)];
        m_beyond_edge[m_beyond_edge[a][edge_to_ll(c,d)]][edge_to_ll(c,d)] = a;

        m_beyond_edge[a][edge_to_ll(b,d)] = m_beyond_edge[c][edge_to_ll(b,d)];
        m_beyond_edge[m_beyond_edge[a][edge_to_ll(b,d)]][edge_to_ll(b,d)] = a;

        m_beyond_edge[d][edge_to_ll(c,a)] = m_beyond_edge[b][edge_to_ll(c,a)];
        m_beyond_edge[m_beyond_edge[d][edge_to_ll(c,a)]][edge_to_ll(c,a)] = d;

        m_beyond_edge[d][edge_to_ll(b,a)] = m_beyond_edge[c][edge_to_ll(b,a)];
        m_beyond_edge[m_beyond_edge[d][edge_to_ll(b,a)]][edge_to_ll(b,a)] = d;
        
        diagonal_flip(a,b,d);
        diagonal_flip(a,c,d);
        
        diagonal_flip(d,b,a);
        diagonal_flip(d,c,a);
    }   

    public:
    DelaunayTriangles(){}

    // id の振り方は特徴的なのに注意
    const ms_Point& operator[](int i){assert(1 <= i && i <= int(m_P.size()));return m_P[i-1];}

    // 点を追加
    void register_point(T x , T y){
        ms_Point p(x,y);
        p.id = m_P.size()+1;// id の振り方に注意
        m_P.push_back(p);
    }


    // 三角形分割する。
    void build(){
        assert(m_P.size() >= 3);

        // 同じ点を含むかチェック
        vector<ms_Point> c_P = m_P;
        sort(c_P.begin() , c_P.end() ,[&](ms_Point& a , ms_Point& b){
            if(a.x == b.x)return bool(a.y < b.y);
            return bool(a.x < b.x);
        });
        for(int i = 0 ; i+1 < int(c_P.size()) ; i++)assert(c_P[i] != c_P[i+1]);
        c_P.clear();

        m_beyond_edge.resize(m_P.size()+1);

        vector<ms_Point> CH = ConvexHull(m_P);
        vector<bool> exist(m_P.size()+1 ,false);
        for(ms_Point& p : CH)exist[p.id] = true;


        // 凸包を、左下の点 m_BasePoint を端点として |CH| - 2 個の三角形に分割する。
        for(int i = 1; i+1 < int(CH.size()) ; i++ ){
            DTreeNode* nd = new DTreeNode(CH[0],CH[i],CH[i+1]);
            m_FirstLayer.push_back(nd);
            NodeTable[triangle_to_ll(CH[0].id , CH[i].id , CH[i+1].id)] = nd;

            if(i+2 < int(CH.size())){
                m_beyond_edge[CH[i].id][edge_to_ll(CH[i+1].id,CH[0].id)] = CH[i+2].id;
                m_beyond_edge[CH[i+2].id][edge_to_ll(CH[i+1].id,CH[0].id)] = CH[i].id;
            }
        }
        for(int i = 1; i < int(CH.size()) ; i++ )diagonal_flip(CH[i].id , CH[0].id , CH[(i+1)%int(CH.size())].id);

        // DTree のうち、要らないノードを削除して、なおかつ平衡化したい (TODO ???)


        vector<int> per(m_P.size());
        for(int i = 0 ; i < int(per.size()) ; i++)per[i] = i;

        // ハックされないように順番をシャッフルする
        for(int c = int(per.size())*3 ; c >= 0 ; c--)swap(per[rand()%int(per.size())],per[rand()%int(per.size())]);
        

        while(!per.empty()){
            ms_Point p = m_P[per.back()];
            per.pop_back();
            if(exist[p.id])continue;

            vector<DTreeNode*> target = this->search(p);
            assert(target.size() <= 2);
            for(DTreeNode* nd : target){
                if(!nd->is_end())continue;
                if(!nd->is_valid())continue;
                assert(nd->data().is_include(p));
                nd->divide(p);
                for(int c = 0 ; c < 3 ; c++){
                    if((*nd)[c] == nullptr)continue;
                    NodeTable[triangle_to_ll((*nd)[c]->data()[0].id , (*nd)[c]->data()[1].id ,(*nd)[c]->data()[2].id)] = (*nd)[c];
                }

                for(int i = 0 ; i < 3 ; i++){
                    // A を主体として辺 BC 方向に p がある状況を考える
                    ms_Point A = nd->data()[(i)%3];// 主体
                    ms_Point B = nd->data()[(i+1)%3];
                    ms_Point C = nd->data()[(i+2)%3];
                    if(!same_line(p,B,C)){
                        // 3 つに分割するとき
                        m_beyond_edge[p.id][edge_to_ll(B.id,C.id)] = m_beyond_edge[A.id][edge_to_ll(B.id,C.id)];
                        m_beyond_edge[m_beyond_edge[p.id][edge_to_ll(B.id,C.id)]][edge_to_ll(B.id,C.id)] = p.id;
                        m_beyond_edge[A.id][edge_to_ll(p.id,B.id)] = C.id;
                        m_beyond_edge[A.id][edge_to_ll(p.id,C.id)] = B.id;
                        // 反対側からの計算は、別の 2 点を主体とした際の計算で行う。
                    }else{
                        // 2 つに分割するとき
                        m_beyond_edge[A.id][edge_to_ll(p.id,B.id)] = m_beyond_edge[A.id][edge_to_ll(B.id,C.id)];
                        m_beyond_edge[m_beyond_edge[A.id][edge_to_ll(p.id,C.id)]][edge_to_ll(p.id,C.id)] = A.id;
                        // 反対側からの計算は、同時に分割される別の三角形から行う。
                    }
                }   
            }
            // p から見て各辺に対して invalid かどうかをチェック
            for(DTreeNode* nd : target){
                for(int i = 0 ; i < 3 ; i++){
                    // A を主体として辺 BC 方向に p がある状況を考える
                    ms_Point A = nd->data()[(i)%3];// 主体
                    ms_Point B = nd->data()[(i+1)%3];
                    ms_Point C = nd->data()[(i+2)%3];
                    diagonal_flip(p.id,A.id,B.id);
                    diagonal_flip(A.id,p.id,B.id);
                    diagonal_flip(B.id,p.id,A.id);
                }
            }
            exist[p.id] = true;
        }


    }

    
    // 分解した三角形リストを取得
    vector<ms_Triangle> crawl(){
        unordered_map<DTreeNode* , bool> visited;
        vector<ms_Triangle> res;
        visited[nullptr] = true;
        for(DTreeNode* start : m_FirstLayer){
            stack<DTreeNode*> st;
            st.push(start);
            while(!st.empty()){
                DTreeNode* now = st.top();
                st.pop();
                if(visited[now])continue;
                visited[now] = true;
                if(now->is_end() && now->is_valid()){
                    res.push_back(now->data());
                }
                for(int c = 0 ; c < 3 ; c++){
                    if(visited[(*now)[c]])continue;// nullptr も弾ける
                    st.push((*now)[c]);
                }
            }
        }
        return res;
    }

    int size(){return m_P.size();}

    /*
        Copyright ©️ (c) NokonoKotlin (okoteiyu) 2025. 
        Released under the MIT license(https://opensource.org/licenses/mit-license.php) 
    */
};



void sometest(int size = 1000){
    int tcnt = 30;
    while(tcnt-->0){
        DelaunayTriangles<long long> D;
        long long n = size;
        unordered_map<long long,unordered_map<long long,bool>> used;
        for(int i = 0 ; i < n ; i++){
            long long x = rand()%100000 - 50000;
            long long y = rand()%100000 - 50000;
            while(used[x][y]){
                x = rand()%100000 - 50000;
                y = rand()%100000 - 50000;
            }
            used[x][y] = true;
            D.register_point(x,y);
        }
        D.build();
        cerr << "test end : " << tcnt << endl;
    }
}


int main(){
    sometest();
    return 0;
}