E088. Нахождение площади объединения треугольников. Метод вертикальной декомпозиции

emaxx algorithm

#algorithm#emaxx#geometry#math

Task

Источник: e-maxx.ru/algo, страница PDF 251.

Даны N треугольников. Требуется найти площадь их объединения.

Решение

Здесь мы рассмотрим метод вертикальной декомпозиции, который в задачах на геометрию часто оказывается очень важным. Итак, у нас имеется N треугольников, которые могут как угодно пересекаться друг с другом. Избавимся от этих пересечений с помощью вертикальной декомпозиции: найдём все точки пересечения всех отрезков (образующих треугольники), и отсортируем найденные точки по их абсциссе. Пусть мы получили некоторый массив B. Будем двигаться по этому массиву. На i-ом шаге рассматриваем элементы B[i] и B[i+1]. Мы имеем вертикальную полосу между прямыми X = B[i] и X = B[i+1], причём, согласно самому построению массива B, внутри этой полосы отрезки никак не пересекаются друг с другом. Следовательно, внутри этой полосы треугольники обрезаются до трапеций, причём стороны этих трапеций внутри полосы не пересекаются вообще. Будем двигаться по сторонам этих трапеций снизу вверх, и складывать площади трапеций, следя за тем, чтобы каждый кусок был учитан ровно один раз. Фактически, этот процесс очень напоминает обработку вложенных скобок. Сложив площади трапеций внутри каждой полосы, и сложив результаты для всех полос, мы и найдём ответ - площадь объединения треугольников. Рассмотрим ещё раз процесс сложения площадей трапеций, уже с точки зрения реализации. Мы перебираем все стороны всех треугольников, и если какая-то сторона (не вертикальная, нам вертикальные стороны не нужны, и даже наоборот, будут сильно мешать) попадает в эту вертикальную полосу (полностью или частично), то мы кладём эту сторону в некоторый вектор, удобнее всего это делать в таком виде: координаты Y в точках пересечения стороны с границами вертикальной полосы, и номер треугольника. После того, как мы построили этот вектор, содержащий куски сторон, сортируем его по значению Y: сначала по левой Y, потом по правой Y. В результате первый в векторе элемент будет содержать нижнюю сторону самой нижней трапеции. Теперь мы просто идём по полученному вектору. Пусть i - текущий элемент; это означает, что i-ый кусок - это нижняя сторона некоторой трапеции, некоторого блока (который может содержать несколько трапеций), площадь которого мы хотим сразу прибавить к ответу. Поэтому мы устанавливаем некий счётчик треугольников равным 1, и поднимаемся по отрезкам вверх, и увеличиваем счётчик, если мы встречаем сторону какого-то треугольника в первый раз, и уменьшаем счётчик, если мы встречаем треугольник во второй раз. Если на каком-то отрезке j счётчик стал равным нулю, то мы нашли верхнюю границу блока - на этом мы останавливаемся, прибавляем площадь трапеции, ограниченной отрезками i и j, и i присваиваем j+1, и повторяем весь процесс заново. Итак, благодаря методу вертикальной декомпозиции мы решили эту задачу, из геометрических примитивов использовав только пересечение двух отрезков.

Реализация

struct segment {

int x1, y1, x2, y2;

};

struct point {

double x, y;

};

struct item {

double y1, y2;

int triangle_id;

};

struct line {

int a, b, c;

};

const double EPS = 1E-7;

void intersect (segment s1, segment s2, vector<point> & res) {

line l1 = { s1.y1-s1.y2, s1.x2-s1.x1, l1.as1.x1+l1.bs1.y1 },

l2 = { s2.y1-s2.y2, s2.x2-s2.x1, l2.as2.x1+l2.bs2.y1 };

double det1 = l1.a * l2.b - l1.b * l2.a;

if (abs (det1) < EPS)  return;

point p = { (l1.c * 1.0 * l2.b - l1.b * 1.0 * l2.c) / det1,

(l1.a * 1.0 * l2.c - l1.c * 1.0 * l2.a) / det1 };

if (p.x >= s1.x1-EPS && p.x <= s1.x2+EPS && p.x >= s2.x1-EPS && p.x

<= s2.x2+EPS)

res.push_back (p);

}

double segment_y (segment s, double x) {

return s.y1 + (s.y2 - s.y1) * (x - s.x1) / (s.x2 - s.x1);

}

bool eq (double a, double b) {

return abs (a-b) < EPS;

}

vector<item> c;

bool cmp_y1_y2 (int i, int j) {

const item & a = c[i];

const item & b = c[j];

return a.y1 < b.y1-EPS || abs (a.y1-b.y1) < EPS && a.y2 < b.y2-EPS;

}

int main() {

int n;

cin >> n;

vector<segment> a (n*3);

for (int i=0; i<n; ++i) {

int x1, y1, x2, y2, x3, y3;

scanf ("%d%d%d%d%d%d", &x1,&y1,&x2,&y2,&x3,&y3);

segment s1 = { x1,y1,x2,y2 };

segment s2 = { x1,y1,x3,y3 };

segment s3 = { x2,y2,x3,y3 };

a[i*3] = s1;

a[i*3+1] = s2;

a[i*3+2] = s3;

}

for (size_t i=0; i<a.size(); ++i)

if (a[i].x1 > a[i].x2)

swap (a[i].x1, a[i].x2), swap (a[i].y1, a[i].y2);

vector<point> b;

b.reserve (nn3);

for (size_t i=0; i<a.size(); ++i)

for (size_t j=i+1; j<a.size(); ++j)

intersect (a[i], a[j], b);

vector<double> xs (b.size());

for (size_t i=0; i<b.size(); ++i)

xs[i] = b[i].x;

sort (xs.begin(), xs.end());

xs.erase (unique (xs.begin(), xs.end(), &eq), xs.end());

double res = 0;

vector<char> used (n);

vector<int> cc (n*3);

c.resize (n*3);

for (size_t i=0; i+1<xs.size(); ++i) {

double x1 = xs[i], x2 = xs[i+1];

size_t csz = 0;

for (size_t j=0; j<a.size(); ++j)

if (a[j].x1 != a[j].x2)

if (a[j].x1 <= x1+EPS && a[j].x2 >= x2-EPS) {

item it = { segment_y (a[j],

x1), segment_y (a[j], x2), (int)j/3 };

cc[csz] = (int)csz;

c[csz++] = it;

}

sort (cc.begin(), cc.begin()+csz, &cmp_y1_y2);

double add_res = 0;

for (size_t j=0; j<csz; ) {

item lower = c[cc[j++]];

used[lower.triangle_id] = true;

int cnt = 1;

while (cnt && j<csz) {

char & cur = used[c[cc[j++]].triangle_id];

cur = !cur;

if (cur)  ++cnt;  else  --cnt;

}

item upper = c[cc[j-1]];

add_res += upper.y1 - lower.y1 + upper.y2 - lower.y2;

}

res += add_res * (x2 - x1) / 2;

}

cout.precision (8);

cout << fixed << res;

}

C# solution

auto-draft, review before submit

using System;
using System.Collections.Generic;
using System.Linq;

public static class AlgorithmDraft
{
    // Auto-generated C# draft from the original e-maxx C/C++ listing. Review before production use.
    struct segment {
            int x1, y1, x2, y2;
    };
    struct point {
            double x, y;
    };
    struct item {
            double y1, y2;
            int triangle_id;
    };
    struct line {
            int a, b, c;
    };
    const double EPS = 1E-7;
    void intersect (segment s1, segment s2, vector<point> & res) {
            line l1 = { s1.y1-s1.y2, s1.x2-s1.x1, l1.a*s1.x1+l1.b*s1.y1 },
                    l2 = { s2.y1-s2.y2, s2.x2-s2.x1, l2.a*s2.x1+l2.b*s2.y1 };
            double det1 = l1.a * l2.b - l1.b * l2.a;
            if (abs (det1) < EPS)  return;
            point p = { (l1.c * 1.0 * l2.b - l1.b * 1.0 * l2.c) / det1,
                    (l1.a * 1.0 * l2.c - l1.c * 1.0 * l2.a) / det1 };
            if (p.x >= s1.x1-EPS && p.x <= s1.x2+EPS && p.x >= s2.x1-EPS && p.x
    <= s2.x2+EPS)
                    res.push_back (p);
    }
    double segment_y (segment s, double x) {
            return s.y1 + (s.y2 - s.y1) * (x - s.x1) / (s.x2 - s.x1);
    }
    bool eq (double a, double b) {
            return abs (a-b) < EPS;
    }
    vector<item> c;
    bool cmp_y1_y2 (int i, int j) {
            const item & a = c[i];
            const item & b = c[j];
            return a.y1 < b.y1-EPS || abs (a.y1-b.y1) < EPS && a.y2 < b.y2-EPS;
    }
    int main() {
            int n;
            cin >> n;
            vector<segment> a (n*3);
            for (int i=0; i<n; ++i) {
                    int x1, y1, x2, y2, x3, y3;
                    scanf ("%d%d%d%d%d%d", &x1,&y1,&x2,&y2,&x3,&y3);
                    segment s1 = { x1,y1,x2,y2 };
                    segment s2 = { x1,y1,x3,y3 };
                    segment s3 = { x2,y2,x3,y3 };
                    a[i*3] = s1;
                    a[i*3+1] = s2;
                    a[i*3+2] = s3;
            }
            for (size_t i=0; i<a.size(); ++i)
                    if (a[i].x1 > a[i].x2)
                            swap (a[i].x1, a[i].x2),  swap (a[i].y1, a[i].y2);
            vector<point> b;
            b.reserve (n*n*3);
            for (size_t i=0; i<a.size(); ++i)
                    for (size_t j=i+1; j<a.size(); ++j)
                            intersect (a[i], a[j], b);
            vector<double> xs (b.size());
            for (size_t i=0; i<b.size(); ++i)
                    xs[i] = b[i].x;
            sort (xs.begin(), xs.end());
            xs.erase (unique (xs.begin(), xs.end(), &eq), xs.end());
            double res = 0;
            List<char> used (n);
            List<int> cc (n*3);
            c.resize (n*3);
            for (size_t i=0; i+1<xs.size(); ++i) {
                    double x1 = xs[i],  x2 = xs[i+1];
                    size_t csz = 0;
                    for (size_t j=0; j<a.size(); ++j)
                            if (a[j].x1 != a[j].x2)
                                    if (a[j].x1 <= x1+EPS && a[j].x2 >= x2-EPS) {
                                            item it = { segment_y (a[j],
    x1), segment_y (a[j], x2), (int)j/3 };
                                            cc[csz] = (int)csz;
                                            c[csz++] = it;
                                    }
                    sort (cc.begin(), cc.begin()+csz, &cmp_y1_y2);
                    double add_res = 0;
                    for (size_t j=0; j<csz; ) {
                            item lower = c[cc[j++]];
                            used[lower.triangle_id] = true;
                            int cnt = 1;
                            while (cnt && j<csz) {
                                    char & cur = used[c[cc[j++]].triangle_id];
                                    cur = !cur;
                                    if (cur)  ++cnt;  else  --cnt;
                            }
                            item upper = c[cc[j-1]];
                            add_res += upper.y1 - lower.y1 + upper.y2 - lower.y2;
                    }
                    res += add_res * (x2 - x1) / 2;
            }
            cout.precision (8);
            Console.WriteLine( fixed << res;
    }
}

C++ solution

matched/original

struct segment {
        int x1, y1, x2, y2;
};
struct point {
        double x, y;
};
struct item {
        double y1, y2;
        int triangle_id;
};
struct line {
        int a, b, c;
};
const double EPS = 1E-7;
void intersect (segment s1, segment s2, vector<point> & res) {
        line l1 = { s1.y1-s1.y2, s1.x2-s1.x1, l1.a*s1.x1+l1.b*s1.y1 },
                l2 = { s2.y1-s2.y2, s2.x2-s2.x1, l2.a*s2.x1+l2.b*s2.y1 };
        double det1 = l1.a * l2.b - l1.b * l2.a;
        if (abs (det1) < EPS)  return;
        point p = { (l1.c * 1.0 * l2.b - l1.b * 1.0 * l2.c) / det1,
                (l1.a * 1.0 * l2.c - l1.c * 1.0 * l2.a) / det1 };
        if (p.x >= s1.x1-EPS && p.x <= s1.x2+EPS && p.x >= s2.x1-EPS && p.x
<= s2.x2+EPS)
                res.push_back (p);
}
double segment_y (segment s, double x) {
        return s.y1 + (s.y2 - s.y1) * (x - s.x1) / (s.x2 - s.x1);
}
bool eq (double a, double b) {
        return abs (a-b) < EPS;
}
vector<item> c;
bool cmp_y1_y2 (int i, int j) {
        const item & a = c[i];
        const item & b = c[j];
        return a.y1 < b.y1-EPS || abs (a.y1-b.y1) < EPS && a.y2 < b.y2-EPS;
}
int main() {
        int n;
        cin >> n;
        vector<segment> a (n*3);
        for (int i=0; i<n; ++i) {
                int x1, y1, x2, y2, x3, y3;
                scanf ("%d%d%d%d%d%d", &x1,&y1,&x2,&y2,&x3,&y3);
                segment s1 = { x1,y1,x2,y2 };
                segment s2 = { x1,y1,x3,y3 };
                segment s3 = { x2,y2,x3,y3 };
                a[i*3] = s1;
                a[i*3+1] = s2;
                a[i*3+2] = s3;
        }
        for (size_t i=0; i<a.size(); ++i)
                if (a[i].x1 > a[i].x2)
                        swap (a[i].x1, a[i].x2),  swap (a[i].y1, a[i].y2);
        vector<point> b;
        b.reserve (n*n*3);
        for (size_t i=0; i<a.size(); ++i)
                for (size_t j=i+1; j<a.size(); ++j)
                        intersect (a[i], a[j], b);
        vector<double> xs (b.size());
        for (size_t i=0; i<b.size(); ++i)
                xs[i] = b[i].x;
        sort (xs.begin(), xs.end());
        xs.erase (unique (xs.begin(), xs.end(), &eq), xs.end());
        double res = 0;
        vector<char> used (n);
        vector<int> cc (n*3);
        c.resize (n*3);
        for (size_t i=0; i+1<xs.size(); ++i) {
                double x1 = xs[i],  x2 = xs[i+1];
                size_t csz = 0;
                for (size_t j=0; j<a.size(); ++j)
                        if (a[j].x1 != a[j].x2)
                                if (a[j].x1 <= x1+EPS && a[j].x2 >= x2-EPS) {
                                        item it = { segment_y (a[j],
x1), segment_y (a[j], x2), (int)j/3 };
                                        cc[csz] = (int)csz;
                                        c[csz++] = it;
                                }
                sort (cc.begin(), cc.begin()+csz, &cmp_y1_y2);
                double add_res = 0;
                for (size_t j=0; j<csz; ) {
                        item lower = c[cc[j++]];
                        used[lower.triangle_id] = true;
                        int cnt = 1;
                        while (cnt && j<csz) {
                                char & cur = used[c[cc[j++]].triangle_id];
                                cur = !cur;
                                if (cur)  ++cnt;  else  --cnt;
                        }
                        item upper = c[cc[j-1]];
                        add_res += upper.y1 - lower.y1 + upper.y2 - lower.y2;
                }
                res += add_res * (x2 - x1) / 2;
        }
        cout.precision (8);
        cout << fixed << res;
}

Java solution

auto-draft, review before submit

import java.util.*;
import java.math.*;

public class AlgorithmDraft {
    // Auto-generated Java draft from the original e-maxx C/C++ listing. Review before production use.
    struct segment {
            int x1, y1, x2, y2;
    };
    struct point {
            double x, y;
    };
    struct item {
            double y1, y2;
            int triangle_id;
    };
    struct line {
            int a, b, c;
    };
    const double EPS = 1E-7;
    void intersect (segment s1, segment s2, vector<point> & res) {
            line l1 = { s1.y1-s1.y2, s1.x2-s1.x1, l1.a*s1.x1+l1.b*s1.y1 },
                    l2 = { s2.y1-s2.y2, s2.x2-s2.x1, l2.a*s2.x1+l2.b*s2.y1 };
            double det1 = l1.a * l2.b - l1.b * l2.a;
            if (abs (det1) < EPS)  return;
            point p = { (l1.c * 1.0 * l2.b - l1.b * 1.0 * l2.c) / det1,
                    (l1.a * 1.0 * l2.c - l1.c * 1.0 * l2.a) / det1 };
            if (p.x >= s1.x1-EPS && p.x <= s1.x2+EPS && p.x >= s2.x1-EPS && p.x
    <= s2.x2+EPS)
                    res.push_back (p);
    }
    double segment_y (segment s, double x) {
            return s.y1 + (s.y2 - s.y1) * (x - s.x1) / (s.x2 - s.x1);
    }
    boolean eq (double a, double b) {
            return abs (a-b) < EPS;
    }
    vector<item> c;
    boolean cmp_y1_y2 (int i, int j) {
            const item & a = c[i];
            const item & b = c[j];
            return a.y1 < b.y1-EPS || abs (a.y1-b.y1) < EPS && a.y2 < b.y2-EPS;
    }
    int main() {
            int n;
            cin >> n;
            vector<segment> a (n*3);
            for (int i=0; i<n; ++i) {
                    int x1, y1, x2, y2, x3, y3;
                    scanf ("%d%d%d%d%d%d", &x1,&y1,&x2,&y2,&x3,&y3);
                    segment s1 = { x1,y1,x2,y2 };
                    segment s2 = { x1,y1,x3,y3 };
                    segment s3 = { x2,y2,x3,y3 };
                    a[i*3] = s1;
                    a[i*3+1] = s2;
                    a[i*3+2] = s3;
            }
            for (size_t i=0; i<a.size(); ++i)
                    if (a[i].x1 > a[i].x2)
                            swap (a[i].x1, a[i].x2),  swap (a[i].y1, a[i].y2);
            vector<point> b;
            b.reserve (n*n*3);
            for (size_t i=0; i<a.size(); ++i)
                    for (size_t j=i+1; j<a.size(); ++j)
                            intersect (a[i], a[j], b);
            vector<double> xs (b.size());
            for (size_t i=0; i<b.size(); ++i)
                    xs[i] = b[i].x;
            sort (xs.begin(), xs.end());
            xs.erase (unique (xs.begin(), xs.end(), &eq), xs.end());
            double res = 0;
            ArrayList<Character> used (n);
            ArrayList<Integer> cc (n*3);
            c.resize (n*3);
            for (size_t i=0; i+1<xs.size(); ++i) {
                    double x1 = xs[i],  x2 = xs[i+1];
                    size_t csz = 0;
                    for (size_t j=0; j<a.size(); ++j)
                            if (a[j].x1 != a[j].x2)
                                    if (a[j].x1 <= x1+EPS && a[j].x2 >= x2-EPS) {
                                            item it = { segment_y (a[j],
    x1), segment_y (a[j], x2), (int)j/3 };
                                            cc[csz] = (int)csz;
                                            c[csz++] = it;
                                    }
                    sort (cc.begin(), cc.begin()+csz, &cmp_y1_y2);
                    double add_res = 0;
                    for (size_t j=0; j<csz; ) {
                            item lower = c[cc[j++]];
                            used[lower.triangle_id] = true;
                            int cnt = 1;
                            while (cnt && j<csz) {
                                    char & cur = used[c[cc[j++]].triangle_id];
                                    cur = !cur;
                                    if (cur)  ++cnt;  else  --cnt;
                            }
                            item upper = c[cc[j-1]];
                            add_res += upper.y1 - lower.y1 + upper.y2 - lower.y2;
                    }
                    res += add_res * (x2 - x1) / 2;
            }
            cout.precision (8);
            System.out.println( fixed << res;
    }
}

Explanation

Материал разбит как алгоритмическая задача: изучить постановку, понять асимптотику и реализовать алгоритм на выбранном языке.