Триангуляция многоугольника минимальной мощности

Дан многоугольник
Необходимо провести в нём некоторое количество отрезков, чтобы многоугольник распался на минимальное количество треугольников

    Очевидно, что этими отрезками будут внутренние диагонали многоугольника.
    Допустим, нам нужно в решении провести диагональ (i, j), тогда многоугольник распадается на два подмногоугольника и задачу можно теперь решать для них отдельно. Получившиеся количества треугольник складываем. За основу берём то, что если многоугольник составлен из отрезков, полностью лежащих внутри исходного многоульника, и в нём три вершины, то минимальная мощность его триангуляции равна 1. Таким образом мы пришли к рекурентному соотношению. Можно написать перебор. Но многие подмногоугольники будут рассматриваться по несколько раз, и их триангуляции от этого меняться не будут. Поэтому мы будем запоминать минимальные триангуляции подмногоугольников. Удобнее всего представлять подмногоугольник как последовательность битов, где каждый бит обозначает присутствует ли соответствующая вершина исходного многоугольника в подмногоугольнике.
    Теперь все подмногоугольники будут рассматриваться по одному разу. Всего их около 2ⁿ, где n - количество вершин. Для каждого подмногоугольника мы будем за O(n²) искать внутреннюю диагональ. Сами внутренние диагонали рационально запомнить препроцессингом за O(n³). В итоге получается асимптотика O(2ⁿ*n² + n³) = O(2ⁿ*n²).

Асимптотика O(2ⁿ*n²).

Листинг C++

#define all(v) (v).size (), (v).end ()
#define sz size ()
#define pb push_back

#define forn(i, n) for ((i) = 0; (i) < (n); ++ (i))
#define forab(i, a, b) for ((i) = (a); (i) <= (b); ++ (i))
#define Forn(i, n) for (int (i) = 0; (i) < (n); ++ (i))
#define Forab(i, a, b) for (int (i) = (a); (i) <= (b); ++ (i))
#define forv(it, v) for (it = (v).begin (); it != (v).end (); ++ it)

const double pi = 3.14159265358979323846264338327;
const double eps = 1e-7;
const double inf = 1e9;

#define eq(a, b) (abs ((a) - (b)) <= eps)
#define ls(a, b) (! eq ((a), (b)) && (a) < (b))
#define gr(a, b) (! eq ((a), (b)) && (a) > (b))

#define sqr(a) ((a) * (a))

double SAT (point a, point b, point c) // ориентированная площадь треугольника - лежит ли вершина b справа от вектора a->c
{
    return 0.5 * (a.x * b.y + b.x * c.y + c.x * a.y - a.y * b.x - b.y * c.x - c.y * a.x);
}

int sSAT (point a, point b, point c)    // знак ориентированной площади треугольника
{
    double s = SAT (a, b, c);
    if (gr (s, 0)) return 1;
    if (ls (s, 0)) return - 1;
    return 0;
}

bool is_convex (vector < point > & v, int a)    // выпуклый ли подмногоугольник
{
    int i, j, k;
    int n = v.size ();

    vector < point > p;
    forn (i, n)
        if (((a >> i) & 1) && (p.sz < 3 || sSAT (p[p.sz - 1], p[p.sz - 2], v[i]) != 0))
            p.pb (v[i]);
    int m = p.sz;
    if (m < 3) return true;

    int s = sSAT (p[0], p[1], p[2]);
    forn (i, m)
    {
        j = (i + 1) % m;
        k = (j + 1) % m;

        if (sSAT (p[i], p[j], p[k]) != s) return false;
    }

    return true;
}

bool is_crossing (point a, point b, point c, point d)    // пересечение двух отрезков, не лежащих на одной прямой
{
    return sSAT (a, c, b) != sSAT (a, d, b) && sSAT (c, a, d) != sSAT (c, b, d);
}

int paint0 (int a, int n, int i, int j)    // покрытие двоичного числа нулями в позициях [i; j]
{
    int a1 = a & ((1 << i) - 1);
    int a2 = a & (((1 << (n - j)) - 1) << (j + 1));
    return a1 | a2;
}

int draw0 (int a, int n, int i, int j)    // покрытие двоичного числа нулями во всех позициях которые не принадлежат [i; j]
{
    int a1 = a & ((1 << i) - 1);
    int a2 = a & (((1 << (n - j)) - 1) << (j + 1));
    return a & (~ (a1 | a2));
}

int cnt_ones (int a, int n)    // количество единиц в двоичном представлении числа
{
    int cnt = 0;
    Forn (i, n)
        cnt += ((a >> i) & 1);
    return cnt;
}

int normal (int p, vector < point > & v)    // удаление из подмногоугольника вершин, которые лежат на одной прямой с двумя соседними
{
    int i, j, k;
    int n = v.size ();

    forn (i, n)
        if ((p >> i) & 1)
        {
            j = (i - 1 + n) % n;
            while (((p >> j) & 1) == 0)
                j = (j - 1 + n) % n;
            k = (i + 1) % n;
            while (((p >> k) & 1) == 0)
                k = (k + 1) % n;
            if (sSAT (v[j], v[i], v[k]) != 0)
                break;
        }
    if (i >= n) return 0;

    int t = 0;
    int i0 = i;
    j = - 1;
    do
    {
        k = (i + 1) % n;
        while (k != j && ((p >> k) & 1) == 0)
            k = (k + 1) % n;

        if (j == - 1 || sSAT (v[j], v[i], v[k]) != 0)
            t |= 1 << i;

        j = i;
        i = k;
    }
    while (i != i0);

    return t;
}

int triangulate_mincnt (vector < point > & v, vector < int > & T1, vector < int > & T2)
{
    int i, j;
    int n = v.size ();

    if (is_convex (v, (1 << n) - 1))
    {
        forn (i, n - 3)
        {
            T1.pb (0);
            T2.pb (i + 2);
        }
        return n - 2;
    }

    char ** D = new char * [n];
    char ** D2 = new char * [n];
    forn (i, n)
    {
        D[i] = new char [n];
        D2[i] = new char [n];
        forn (j, n)
        {
            D[i][j] = is_diagonal (v, i, j);
            D2[i][j] = point_in_polygon (part_segment (v[i], v[j], 1, 1), v);
        }
    }


    int P = 1 << n;
    int * M = new int [P];
    int * prev1 = new int [P];
    int * prev2 = new int [P];
    forn (i, P)
        M[i] = inf, prev1[i] = - 1, prev2[i] = - 1;
    int * q = new int [n];
    forn (i, n)
        q[i] = 0;

    Forab (p, 7, (1 << n) - 1)
    {
        if (cnt_ones (p, n) < 3) continue;    // если точек меньше трёх - это не многоугольник
        forn (i, n)                // проверяем, все ли рёбра многоугольника являются внутренними диагоналями
            if ((p >> i) & 1)    // или сторонами исходного многоугольника
            {
                j = (i + 1) % n;
                while (j != i && ((p >> j) & 1) == 0)
                    j = (j + 1) % n;
                if (j != i && D2[i][j] == 0)
                    break;
            }
        if (i < n) continue;


        int t = normal (p, v);    // сокращаем многоугольник, засчёт удаления точек, лежащих на одной прямой
        if (t != p)    // если многоугольник сократился, значит мы его уже рассматривали
        {
            M[p] = M[t];
            prev1[p] = prev1[t];
            prev2[p] = prev2[t];
            continue;
        }
        if (is_convex (v, p))    // если многоугольник выпуклый - проводим любую диагональ
        {
            forn (i, n)
                if ((p >> i) & 1)
                    break;
            prev1[p] = i;

            forab (j, i + 1, n - 1)
                if ((p >> j) & 1)
                    break;
            forab (i, j + 1, n - 1)
                if ((p >> i) & 1)
                    break;
            prev2[p] = i;

            M[p] = cnt_ones (p, n) - 2;

            continue;
        }


        int m = 0;
        forn (i, n)
            if ((p >> i) & 1)
                q[m ++] = i;

        forab (i, 0, m - 2)    // перебираем все диагонали и проводим их, оценивая
            forab (j, i + 1, m - 1)    // мощности триангуляций полученных частей
                if (D[q[i]][q[j]])
                {
                    int p1 = paint0 (p, n, q[i] + 1, q[j] - 1);
                    int p2 = draw0 (p, n, q[i], q[j]);

                    if (M[p] > M[p1] + M[p2])
                    {
                        M[p] = M[p1] + M[p2];
                        prev1[p] = q[i];
                        prev2[p] = q[j];
                    }
                }
    }

    vector < int > st;
    st.push_back ((1 << n) - 1);
    while (st.size ())    // находим, какие мы провели диагонали с помощью поиска в глубину
    {
        int p = st.back ();
        st.pop_back ();

        if (M[p] == 1 || M[p] < 1 || M[p] > n) continue;

        T1.push_back (prev1[p]);
        T2.push_back (prev2[p]);

        int p1 = paint0 (p, n, prev1[p] + 1, prev2[p] - 1);
        int p2 = draw0 (p, n, prev1[p], prev2[p]);

        st.push_back (p1);
        st.push_back (p2);
    }

    return M[(1 << n) - 1];
}

13:28
04.08.2009

По всем вопросам обращаться: rumterg@gmail.com