Агуулгыг алгасах

Делоне гурвалжинчлал ба Воронойн диаграм

Хавтгай дээрх цэгүүдийн $\{p_i\}$ олонлогийг авч үзье. $\{p_i\}$-ийн Воронойн диаграм $V(\{p_i\})$ гэдэг нь хавтгайг $n$ ширхэг $V_i$ муж болгон хуваасан хуваалт бөгөөд энд $V_i = \{p\in\mathbb{R}^2;\ \rho(p, p_i) = \min\ \rho(p, p_k)\}$. Воронойн диаграмын нүднүүд нь олон өнцөгтүүд (магадгүй төгсгөлгүй) байна. $\{p_i\}$-ийн Делоне гурвалжинчлал $D(\{p_i\})$ гэдэг нь $p_i$ цэг бүр $T \in D(\{p_i\})$ гурвалжин бүрийн багтаасан тойргийн гадна эсвэл зааг дээр орших гурвалжинчлал юм.

Воронойн диаграм холбоост биш байх ба Делоне гурвалжинчлал оршин байхгүй байх эвгүй degenerate case бий. Энэ тохиолдол нь бүх цэг коллинеар байх үед юм.

Шинж чанарууд

Делоне гурвалжинчлал нь боломжит бүх гурвалжинчлалын дундаас хамгийн бага өнцгийг максимумчилдаг.

Цэгүүдийн олонлогийн Евклидийн хамгийн бага тэлэх мод нь түүний Делоне гурвалжинчлалын ирмэгүүдийн дэд олонлог юм.

Хос чанар

$\{p_i\}$ коллинеар биш бөгөөд $\{p_i\}$-ийн дундаас ямар ч дөрвөн цэг нэг тойрог дээр оршихгүй гэж үзье. Тэгвэл $V(\{p_i\})$ ба $D(\{p_i\})$ хос байх тул хэрэв бид тэдгээрийн нэгийг авбал нөгөөг нь $O(n)$-д авч болно. Хэрэв тийм биш бол яах вэ? Коллинеар тохиолдлыг хялбархан боловсруулж болно. Эс бөгөөс $V$ ба $D'$ хос байх ба энд $D'$-г $D$-ээс энэ ирмэг дээрх хоёр гурвалжин багтаасан тойргоо хуваалцаж байх бүх ирмэгийг хасах замаар авна.

Делоне ба Воронойг байгуулах

Хос чанараас болж бидэнд $V$ ба $D$-ийн зөвхөн нэгийг тооцоолох хурдан алгоритм л хэрэгтэй. Бид $D(\{p_i\})$-г $O(n\log n)$-д хэрхэн байгуулахыг тайлбарлана. Гурвалжинчлалыг Guibas, Stolfi нарын хуваа ба ялагтун алгоритмаар байгуулна.

Quad-edge өгөгдлийн бүтэц

Алгоритмын явцад $D$ нь quad-edge өгөгдлийн бүтцэд хадгалагдана. Энэ бүтцийг зурган дээр тайлбарласан:

Quad-Edge

Алгоритмд бид ирмэгүүд дээр дараах функцүүдийг ашиглана:

  1. make_edge(a, b)
    Энэ функц a цэгээс b цэг хүртэл тусгаарлагдсан ирмэгийг түүний урвуу ирмэг ба хоёр хос ирмэгийн хамт үүсгэнэ.
  2. splice(a, b)
    Энэ бол алгоритмын гол функц юм. Энэ нь a->Onext-г b->Onext-тэй, a->Onext->Rot->Onext-г b->Onext->Rot->Onext-тэй солино.
  3. delete_edge(e)
    Энэ функц e-г гурвалжинчлалаас устгана. e-г устгахын тулд бид зүгээр л splice(e, e->Oprev) ба splice(e->Rev, e->Rev->Oprev)-г дуудаж болно.
  4. connect(a, b)
    Энэ функц a, b, e бүгд ижил зүүн нүүртэй байхаар a->Dest-ээс b->Org хүртэл шинэ e ирмэг үүсгэнэ. Үүний тулд бид e = make_edge(a->Dest, b->Org), splice(e, a->Lnext) ба splice(e->Rev, b)-г дуудна.

Алгоритм

Алгоритм гурвалжинчлалыг тооцоолж, хоёр quad-edge буцаана: хамгийн зүүн оройноос гарах цагийн зүүний эсрэг гүдгэр бүрхүүлийн ирмэг ба хамгийн баруун оройноос гарах цагийн зүүний дагуу гүдгэр бүрхүүлийн ирмэг.

Бүх цэгийг x-ээр, хэрэв $x_1 = x_2$ бол y-ээр эрэмбэлье. Ямар нэг $(l, r)$ хэрчмийн хувьд бодлогыг бодъё (эхэндээ $(l, r) = (0, n - 1)$). Хэрэв $r - l + 1 = 2$ бол бид $(p[l], p[r])$ ирмэгийг нэмээд буцна. Хэрэв $r - l + 1 = 3$ бол бид эхлээд $(p[l], p[l + 1])$ ба $(p[l + 1], p[r])$ ирмэгүүдийг нэмнэ. Мөн бид тэдгээрийг splice(a->Rev, b) ашиглан холбох ёстой. Одоо бид гурвалжныг хаах ёстой. Бидний дараагийн үйлдэл $p[l], p[l + 1], p[r]$-ийн чиглэлээс хамаарна. Хэрэв тэдгээр коллинеар бол бид гурвалжин үүсгэж чадахгүй тул зүгээр л (a, b->Rev) буцаана. Эс бөгөөс бид connect(b, a)-г дуудаж шинэ c ирмэг үүсгэнэ. Хэрэв цэгүүд цагийн зүүний эсрэг чиглэлтэй бол бид (a, b->Rev) буцаана. Эс бөгөөс бид (c->Rev, c) буцаана.

Одоо $r - l + 1 \ge 4$ гэж үзье. Эхлээд $L = (l, \frac{l + r}{2})$ ба $R = (\frac{l + r}{2} + 1, r)$-г рекурсивээр бодъё. Одоо бид эдгээр гурвалжинчлалыг нэг гурвалжинчлал болгон нийлүүлэх ёстой. Бидний цэгүүд эрэмбэлэгдсэн болохыг анхаарна уу, тиймээс нийлүүлэх явцад бид L-ээс R рүү ирмэгүүд (cross ирмэг гэж нэрлэгддэг) нэмэх ба L-ээс L рүү, R-ээс R рүү чиглэсэн зарим ирмэгийг хасна. cross ирмэгүүдийн бүтэц ямар вэ? Эдгээр бүх ирмэг y тэнхлэгтэй параллель бөгөөд хуваах x утган дээр байрлах шулууныг огтлох ёстой. Энэ нь cross ирмэгүүдийн шугаман эрэмбийг тогтоох тул бид дараалсан cross ирмэгүүд, хамгийн доод cross ирмэг гэх мэтийн талаар ярьж болно. Алгоритм cross ирмэгүүдийг өсөх дарааллаар нэмнэ. Дурын хоёр зэргэлдээ cross ирмэг нийтлэг үзүүрийн цэгтэй байх ба тэдгээрийн тодорхойлох гурвалжны гурав дахь тал нь L-ээс L рүү эсвэл R-ээс R рүү явахыг анхаарна уу. Одоогийн cross ирмэгийг base гэж нэрлэе. base-ийн залгамжлагч нь base-ийн зүүн үзүүрийн цэгээс баруун үзүүрийн цэгийн R-хөршүүдийн нэг рүү, эсвэл эсрэгээр явна. base ба өмнөх cross ирмэгийн багтаасан тойргийг авч үз. Энэ тойрог base-г хөвч болгон агуулах боловч Oy чиглэлд цааш орших бусад тойрог болон хувирна гэж үзье. Бидний тойрог хэсэг хугацаанд дээшилнэ, гэвч base нь L ба R-ийн дээд шүргэгч биш л бол бид багтаасан тойрогтоо ямар ч цэггүй шинэ гурвалжин үүсгэх L эсвэл R-д харьяалагдах цэгтэй тулгарна. Энэ гурвалжны шинэ L-R ирмэг нь дараагийн нэмэгдэх cross ирмэг юм. Үүнийг үр ашигтай хийхийн тулд бид lcand нь энэ үйл явцад тааралдсан эхний L цэгийг, rcand нь эхний R цэгийг заахаар lcand ба rcand гэсэн хоёр ирмэгийг тооцоолно. Дараа нь бид эхлээд тааралдах нэгийг сонгоно. Эхэндээ base нь L ба R-ийн доод шүргэгчийг заана.

Implementation

Note that the implementation of the in_circle function is GCC-specific.

typedef long long ll;

bool ge(const ll& a, const ll& b) { return a >= b; }
bool le(const ll& a, const ll& b) { return a <= b; }
bool eq(const ll& a, const ll& b) { return a == b; }
bool gt(const ll& a, const ll& b) { return a > b; }
bool lt(const ll& a, const ll& b) { return a < b; }
int sgn(const ll& a) { return a >= 0 ? a ? 1 : 0 : -1; }

struct pt {
    ll x, y;
    pt() { }
    pt(ll _x, ll _y) : x(_x), y(_y) { }
    pt operator-(const pt& p) const {
        return pt(x - p.x, y - p.y);
    }
    ll cross(const pt& p) const {
        return x * p.y - y * p.x;
    }
    ll cross(const pt& a, const pt& b) const {
        return (a - *this).cross(b - *this);
    }
    ll dot(const pt& p) const {
        return x * p.x + y * p.y;
    }
    ll dot(const pt& a, const pt& b) const {
        return (a - *this).dot(b - *this);
    }
    ll sqrLength() const {
        return this->dot(*this);
    }
    bool operator==(const pt& p) const {
        return eq(x, p.x) && eq(y, p.y);
    }
};

const pt inf_pt = pt(1e18, 1e18);

struct QuadEdge {
    pt origin;
    QuadEdge* rot = nullptr;
    QuadEdge* onext = nullptr;
    bool used = false;
    QuadEdge* rev() const {
        return rot->rot;
    }
    QuadEdge* lnext() const {
        return rot->rev()->onext->rot;
    }
    QuadEdge* oprev() const {
        return rot->onext->rot;
    }
    pt dest() const {
        return rev()->origin;
    }
};

QuadEdge* make_edge(pt from, pt to) {
    QuadEdge* e1 = new QuadEdge;
    QuadEdge* e2 = new QuadEdge;
    QuadEdge* e3 = new QuadEdge;
    QuadEdge* e4 = new QuadEdge;
    e1->origin = from;
    e2->origin = to;
    e3->origin = e4->origin = inf_pt;
    e1->rot = e3;
    e2->rot = e4;
    e3->rot = e2;
    e4->rot = e1;
    e1->onext = e1;
    e2->onext = e2;
    e3->onext = e4;
    e4->onext = e3;
    return e1;
}

void splice(QuadEdge* a, QuadEdge* b) {
    swap(a->onext->rot->onext, b->onext->rot->onext);
    swap(a->onext, b->onext);
}

void delete_edge(QuadEdge* e) {
    splice(e, e->oprev());
    splice(e->rev(), e->rev()->oprev());
    delete e->rev()->rot;
    delete e->rev();
    delete e->rot;
    delete e;
}

QuadEdge* connect(QuadEdge* a, QuadEdge* b) {
    QuadEdge* e = make_edge(a->dest(), b->origin);
    splice(e, a->lnext());
    splice(e->rev(), b);
    return e;
}

bool left_of(pt p, QuadEdge* e) {
    return gt(p.cross(e->origin, e->dest()), 0);
}

bool right_of(pt p, QuadEdge* e) {
    return lt(p.cross(e->origin, e->dest()), 0);
}

template <class T>
T det3(T a1, T a2, T a3, T b1, T b2, T b3, T c1, T c2, T c3) {
    return a1 * (b2 * c3 - c2 * b3) - a2 * (b1 * c3 - c1 * b3) +
           a3 * (b1 * c2 - c1 * b2);
}

bool in_circle(pt a, pt b, pt c, pt d) {
// If there is __int128, calculate directly.
// Otherwise, calculate angles.
#if defined(__LP64__) || defined(_WIN64)
    __int128 det = -det3<__int128>(b.x, b.y, b.sqrLength(), c.x, c.y,
                                   c.sqrLength(), d.x, d.y, d.sqrLength());
    det += det3<__int128>(a.x, a.y, a.sqrLength(), c.x, c.y, c.sqrLength(), d.x,
                          d.y, d.sqrLength());
    det -= det3<__int128>(a.x, a.y, a.sqrLength(), b.x, b.y, b.sqrLength(), d.x,
                          d.y, d.sqrLength());
    det += det3<__int128>(a.x, a.y, a.sqrLength(), b.x, b.y, b.sqrLength(), c.x,
                          c.y, c.sqrLength());
    return det > 0;
#else
    auto ang = [](pt l, pt mid, pt r) {
        ll x = mid.dot(l, r);
        ll y = mid.cross(l, r);
        long double res = atan2((long double)x, (long double)y);
        return res;
    };
    long double kek = ang(a, b, c) + ang(c, d, a) - ang(b, c, d) - ang(d, a, b);
    if (kek > 1e-8)
        return true;
    else
        return false;
#endif
}

pair<QuadEdge*, QuadEdge*> build_tr(int l, int r, vector<pt>& p) {
    if (r - l + 1 == 2) {
        QuadEdge* res = make_edge(p[l], p[r]);
        return make_pair(res, res->rev());
    }
    if (r - l + 1 == 3) {
        QuadEdge *a = make_edge(p[l], p[l + 1]), *b = make_edge(p[l + 1], p[r]);
        splice(a->rev(), b);
        int sg = sgn(p[l].cross(p[l + 1], p[r]));
        if (sg == 0)
            return make_pair(a, b->rev());
        QuadEdge* c = connect(b, a);
        if (sg == 1)
            return make_pair(a, b->rev());
        else
            return make_pair(c->rev(), c);
    }
    int mid = (l + r) / 2;
    QuadEdge *ldo, *ldi, *rdo, *rdi;
    tie(ldo, ldi) = build_tr(l, mid, p);
    tie(rdi, rdo) = build_tr(mid + 1, r, p);
    while (true) {
        if (left_of(rdi->origin, ldi)) {
            ldi = ldi->lnext();
            continue;
        }
        if (right_of(ldi->origin, rdi)) {
            rdi = rdi->rev()->onext;
            continue;
        }
        break;
    }
    QuadEdge* basel = connect(rdi->rev(), ldi);
    auto valid = [&basel](QuadEdge* e) { return right_of(e->dest(), basel); };
    if (ldi->origin == ldo->origin)
        ldo = basel->rev();
    if (rdi->origin == rdo->origin)
        rdo = basel;
    while (true) {
        QuadEdge* lcand = basel->rev()->onext;
        if (valid(lcand)) {
            while (in_circle(basel->dest(), basel->origin, lcand->dest(),
                             lcand->onext->dest())) {
                QuadEdge* t = lcand->onext;
                delete_edge(lcand);
                lcand = t;
            }
        }
        QuadEdge* rcand = basel->oprev();
        if (valid(rcand)) {
            while (in_circle(basel->dest(), basel->origin, rcand->dest(),
                             rcand->oprev()->dest())) {
                QuadEdge* t = rcand->oprev();
                delete_edge(rcand);
                rcand = t;
            }
        }
        if (!valid(lcand) && !valid(rcand))
            break;
        if (!valid(lcand) ||
            (valid(rcand) && in_circle(lcand->dest(), lcand->origin,
                                       rcand->origin, rcand->dest())))
            basel = connect(rcand, basel->rev());
        else
            basel = connect(basel->rev(), lcand->rev());
    }
    return make_pair(ldo, rdo);
}

vector<tuple<pt, pt, pt>> delaunay(vector<pt> p) {
    sort(p.begin(), p.end(), [](const pt& a, const pt& b) {
        return lt(a.x, b.x) || (eq(a.x, b.x) && lt(a.y, b.y));
    });
    auto res = build_tr(0, (int)p.size() - 1, p);
    QuadEdge* e = res.first;
    vector<QuadEdge*> edges = {e};
    while (lt(e->onext->dest().cross(e->dest(), e->origin), 0))
        e = e->onext;
    auto add = [&p, &e, &edges]() {
        QuadEdge* curr = e;
        do {
            curr->used = true;
            p.push_back(curr->origin);
            edges.push_back(curr->rev());
            curr = curr->lnext();
        } while (curr != e);
    };
    add();
    p.clear();
    int kek = 0;
    while (kek < (int)edges.size()) {
        if (!(e = edges[kek++])->used)
            add();
    }
    vector<tuple<pt, pt, pt>> ans;
    for (int i = 0; i < (int)p.size(); i += 3) {
        ans.push_back(make_tuple(p[i], p[i + 1], p[i + 2]));
    }
    return ans;
}

Бодлогууд