Пересечение, объединение и разность многоугольников

Всем доброго времени!
Задачка такая: даны n простых многоугольников. Нужно найти их пересечение, объединение, разность.
Все, кто сталкивался с этой задачкой поделитесь пожалуйста опытом
Выручайте, сдача через 3 дня!

Сообщение отредактировано: Nymph666 - 24.12.07, 03:47

гуглить:
Alan Murta GPC (General Polygon Clipper)
чья-то библиотека PolyBoolean - на ее страничке ссылки и на другие

Цитата MBo @ 24.12.07, 04:53

гуглить:
Alan Murta GPC (General Polygon Clipper)
чья-то библиотека PolyBoolean - на ее страничке ссылки и на другие

Готовая библиотека не подходит...

Если под словом простых понимается выпуклых, то ход такой :
1.Пересечение n штук надо? Всех ? Так :
a.Берём первый.
б.Берём второй
в.пересекая их, получаем снова выпуклый, заименованный нами как N.
г.берём третий, четвёртый(цикл пошёл) и, пересекая с N, получаем выпуклый, который снов записываем в N.
д.Ответ = N
// пересекать два многоугольника так :
бежать по точкам второго и смотреть : не попали ли мы во внутрь первого.
Если не попали(никогда), то пересечение пусто или совпадает с меньшим!
Иначе, стоп.бежать от попавшей вершины до первой вышедшей.
Если такой нету, то мы полностью внутри.
Иначе, переключиться на первый контур, бежать... переключиться на второй, ... - нетрудно.
2.С объединением не очень понятно, ибо это может оказаться и не многоугольник(и тога совсем тяжко).
А если многоугольник, то так можно : найти пересекающуюся пару. её объединить.
Найти пересекающийся с результатом "простой" многоугольник, его объединить с результатом, записать в результат. и т.д.
3.Ну разности вообще всякие возможны(в группе из n штук). Но, по сути, бежим по границе и переходим на границу другого.

Цитата Славян @ 24.12.07, 17:27

Если под словом простых понимается выпуклых,

Под словом простых понимается и выпуклых, и вогнутых, но без самопересечений...

Есть простой вариант, использующий BSP-дерево.
Сторим BSP-дерево для многоугольника A и используя его разделяем многоугольник B (точнее его границу) на две последовательности отрезков B- и B+, которые лежат соответственно внутри и снаружи многоугольника A.
Аналогично, построив BSP-дерево для B, получим две последовательности A- и A+.
Тогда, пересечение есть объединение A- и B-.
Объединение - объединение A+ и B+.
И разность (A-B), соответственно, A+ и B-.

Теперь из получившегося набора отрезков, надо собрать грацицу многоугольника, но это уже совсем просто.

Вот на рисунке пример. Мы, построив bsp дерево для ABCD, начинаем делить EFGH

В результате, в качестве B+, получим JF, FG, GK, KH, HL (то что обведено на левой схеме).
Потом добавляем A+: AB, BJ, LD, DA (не объязательно в таком виде, один их этих отрезков может быть разбит на два, как это было с орезком GH, который разбился на GK и KH, хотя оба они и вошли в B+).
И вот (A+ + B+) = {JF, FG, GK, KH, HL, AB, BJ, LD, DA} - объединение многоугольников.

Цитата albom @ 28.12.07, 08:09

Есть простой вариант, использующий BSP-дерево.

А можно по подробнее что такое и как строить BSP-дерево?
Желательно на псевдокоде и т.п.

Binary Space Partitioning

Строить bsp для полигонов совсем просто.
Рекурсивный алгоритм такой:

procedure Split(lines) return Tree;
    если (lines не пусто) то
    
        splitLine  <- выбранный элемент из lines (выбирать надо особо хитро, хотя можно просто взять случайный)
        spiltPlane <- гиперплоскость, содержащая splitLine
            
        Для каждого l из lines:
            если (пересечение l и splitPlane не пусто) то
                разбить отрезок l на два l+ и l-
                удалить из lines l
                добавить в lines l+ и l-
 
        Делим отрезки lines на два списка:
            positive <- отрезки лежащие в положительной полуплосткости
            negative <- отрезки лежащие в отрицательной полуплосткости
        
        P <- Split(positive)
        N <- Split(negative)
        
        вернуть Tree(splitPlane, P, N)
    иначе
        вернуть nil;
end Split

Где Tree - это, как ни стренно, бинарное дерево.
Где в узле хранится splitPlane, и есть ссылки на детей P и N.

Большое спасибо, буду разбираться.

Как я поняла:
1.Берем один из отрезков одного из многоугольников
2.Разбиваем им оставшееся множество отрезков на 2 части P и N.
3.Затем аналогично выбираем по отрезку из P и N и снова разбиваем то подмножество(P или N) на две части и так пока не кончатся отрезки
но как определить какие отрезки входят в P, а какие в N?
Или я что-то не так понимаю? (Читала о применении BSP в 3D рендеринге, но есть наблюдатель относительно которого определяется видимость)

Цитата Nymph666 @ 08.01.08, 05:44

но как определить какие отрезки входят в P, а какие в N?

Отрезок (сплиттер), делит нашу плоскость на две полуплоскости: положительную и отрицательную. Пусть прямая сплиттера задается уравнением a*x + b*y + c = 0. Тогда если мы подставим координаты отрезка в это уравненние, и результат будет больше нуля, то это положительная полуплоскость, иначе отрицательная.

Теперь, наш многоугольник делит все пространство на две области (возможно несвязные), будем называть положительной ту, что не пренадлежит многоугольнику и отрицательной другую. Тогда что-бы знаки полуплоскостей и областей рядом со сплиттером совпадали, достаточно положить, что бы нормаль сплиттера смотрела наружу многоугольника. А для этого просто скажем, что стороны многоугольника мы всегда обходим против часовой стрелки, этого достаточно.

Разрешите, я сразу кодом буду писать, а то словами плохо получается...
Пара определений:

class Segment; // это отрезок
typedef std::vector<Segment> Segments; // а это массив отрезков, для нас это стороны многоугольника

Вот наше дерево:

struct BSPNode
{
    Segment splitter; // разделяющий элемент, отрезок в нашем случае
    BSPNode *positive, *negative;
 
    BSPNode (const Segment& Splitter,  BSPNode* p,  BSPNode* n)
        : splitter(Splitter), positive(p), negative(n)
    {   }
};

Тогда вот так мы строим bsp-дерево:

BSPNode* buildTree(const Segments &segments)
{
    if (segments.size() == 0)  // отрезков нет, нет и дерева
        return NULL;
    Segment spl = segments[0]; // выбираем сплиттер
    Segments pList, nList;     // списки N и P
    for (size_t i = 1; i < segments.size(); i++) // перебираем все отрезки, кроме сплиттера
    {
        const Segment& cur = segments[i];
        Segment p, n;
        bool    pp, nn;
// делим отрезок cur на два, p - отрезок лежащий в положительной полуплоскости отн-но сплиттера, n - соответственно в отрицательной
// Конечно если отрезок не пересекает сплиттер, то он целиком попадает в одну из полуплоскостей
        split(cur, spl, p, n, pp, nn);  
        if (pp)  // в положительную полуплоскость что-то попало
            pList.push_back(p); // добавили в список P
        if (nn) // в отрицательную полуплоскость что-то попало
            nList.push_back(n); // добавили в список N
    }
    return new BSPNode(spl, buildTree(pList), buildTree(nList)); // создали узел дерева, и рекурсивно начали делить списки P и N
}

Вот ad-hoc вариант процедуры split:

struct SegmentIntersection  // пересечение отрезков
{
    bool valid; // флаг, есть оно или нет
    vec2 point; // точка пересечениея
    SegmentIntersection (bool v, vec2 p)
        : valid(v), point(p)
    {  }
    SegmentIntersection (bool v)
        : valid(v)
    {  }
};
 
// пересечение отрезков и нахождение точки пересечения
// используется параметрическое представление отрезка a+d*t, где a - начало отрезка, d - направляющая, а t - параметр ( 0 <= t <= 1)
bool intersect(const Segment& a, const Segment &b, float &ta, float &tb)
{
    vec2 d1 = a.dir();
    vec2 d2 = b.dir();
    float d = d1 ^ d2;
    if ( fabs(d) < eps)
        return false;
    float dax = d1.x;
    float day = d1.y;
    float dbx = d2.x;
    float dby = d2.y;
    float xa0 = a.begin().x;
    float ya0 = a.begin().y;
    float xb0 = b.begin().x;
    float yb0 = b.begin().y;
    float t1 = - (dbx*yb0-dbx*ya0-xb0*dby+xa0*dby) / d;
    float t2 = - (-day*xb0+dax*yb0-dax*ya0+day*xa0) / d;
    ta = t1 / a.length();
    tb = t2 / b.length();
    return true;
}
// пересечение отрезка a со сплиттером(!) b
SegmentIntersection intersectEx(const Segment& a, const Segment& b)
{
    float t1, t2;
    if (intersect(a, b, t1, t2) == false)
        return SegmentIntersection(false);
    if ( (t1 <= 0) || (t1 >= 1) )
        return SegmentIntersection(false); // пересекаются не отрезки, а прямые их содержащие
    else
        return SegmentIntersection(true, a.begin() +  a.dir() * t1 * a.length()); // ок, есть точка пересечения
}
// v - отрезок
// s - сплиттер
// p - "положительный" отрезок
// n - "отрицательный" отрезок
// bp - признак, что положитльный отрезок существует
// bn - признак, что не положитльный отрезок существует
void split(Segment v, Segment s, Segment &p, Segment &n, bool &bp, bool &bn)
{
    SegmentIntersection i = intersectEx(v, s);
 
    if (i.valid == false)
    { // пересечения нет
        bool pos = ( s.toLine()[v.begin()] + s.toLine()[v.end()] )  >= 0;
        if (pos)
        { // отрезок целиком в положительной полуплосткости
            p = v;
            bp = true;
            bn = false;
        }
        else
        { // отрезок целиком в отрицательной полуплосткости
            n = v;
            bp = false;
            bn = true;
        }
    }
    else
    { // есть пересечение
        Segment s1 = Segment(v.begin(), i.point);
        Segment s2 = Segment(i.point, v.end());
        bp = true;
        bn = true;
        if (s.toLine()[s1.begin()] + s.toLine()[s1.end()] >= 0) // смотрим, какой отрезок лежит в полож. а какой в отр. полуплосткости
        {
            p = s1;
            n = s2;
        }
        else
        {
            p = s2;
            n = s1;
        }
 
    }
}

Уф...
Дерево мы построили, теперь пора делить:

// tree - дерево
// segments - набор отрезков, которые будем делить
// p, n - списки "положительных" и "отрицательных" отрезков
void splitSegments(const BSPNode* tree, const Segments &segments, Segments &p, Segments &n)
{
    Segments pList, nList;
    for (size_t i = 0; i < segments.size(); i++)
    {
        const Segment& cur = segments[i]; // берем  каждый отрезок
        Segment p, n;
        bool    pp, nn;
        split(cur, tree->splitter, p, n, pp, nn);  // делим его сплиттером
        if (pp)
            pList.push_back(p);
        if (nn)
            nList.push_back(n);
    }
    if ( tree->positive != NULL )
        splitSegments(tree->positive, pList, p, n); // делим P список рекурсивно
    else
        p.insert(p.end(), pList.begin(), pList.end()); // дерево кончилось => все что попалов в P список лежит вне многоугольника
 
    if ( tree->negative != NULL )
        splitSegments(tree->negative, nList, p, n);
    else
        n.insert(n.end(), nList.begin(), nList.end());
}

Вот так пользоватся:

struct Polygon
{
    Polygon(std::vector<vec2> v)
    {
        for (size_t i = 0; i < v.size(); i++)
            segments.push_back( Segment(v[i], v[ (i + 1) % v.size() ] ) );
    }
    Segments segments;
};
 
 
void print(Segments &s)
{
    for (size_t i = 0; i < s.size(); i++)
        printf("(%4.2f %4.2f) - (%4.2f %4.2f)\n", s[i].begin().x, s[i].begin().y,  s[i].end().x,  s[i].end().y );
 
}
 
 
int main(void)
{
 
// многоугольник A
    std::vector<vec2> PA;
    PA.push_back(vec2(  0, 0 ) );
    PA.push_back(vec2(  4, 2 ) );
    PA.push_back(vec2(  2, 2 ) );
    PA.push_back(vec2(  2, 4 ) );
// многоугольник B
    std::vector<vec2> PB;
    PB.push_back(vec2(  1, 1 ) );
    PB.push_back(vec2(  3, 1 ) );
    PB.push_back(vec2(  3, 3 ) );
    PB.push_back(vec2(  1, 3 ) );
 
    Polygon A(PA); // создали многоугольник по точкам
    Polygon B(PB);
 
    BSPNode* treeA = buildTree(A.segments); // построили дерево для многоугольника A
    Segments BPlus, BMinus;
    splitSegments(treeA, B.segments, BPlus, BMinus); // разделили многоугольник B
        
 
    printf("B+\n");
    print(BPlus);
    printf("B-\n");
    print(BMinus);
 
    return 0;
}

Ну и остатки

#define eps 1e-5
 
inline vec2 rotatem90(vec2 v)
{
    return vec2(v.y, -v.x);
}
 
class Line
{
public:
    vec2 normal;
    float c;
    vec2 dir;
    vec2 begin;
    Line(vec2 A, vec2 B)
    {
        dir = B - A;
        dir.normalize();
        normal = rotatem90(dir);
        c = - (A * normal);
        begin = A;
    }
 
    Line(vec2 Normal, float C)
        : normal(Normal), c(C)
    {   }
 
    float operator[](vec2 point) const { return normal * point + c; }
};
 
class Segment
{
 
    vec2 a, b;
    Line line;
public:
    Segment()
        :line(vec2(0,0), 0), a(0,0), b(0,0)
    {   }
    
    Segment(vec2 A, vec2 B)
        : a(A), b(B), line(A, B)
    {   }
 
    const Line& toLine() const
    {
        return line;
    }
    const vec2& begin() const { return a; }
    const vec2& end()   const { return b; }
    const vec2 dir() const
    {
        vec2 d = vec2(b - a);
        d.normalize();
        return d;
    }
 
    float length() const
    {
        return (b - a).length();
    }
 
};

А вот собственно и пересечение/объединение/разность:

void deleteTree(BSPNode *tree)
{
    if (tree != NULL)
    {
        deleteTree(tree->positive);
        deleteTree(tree->negative);
        delete tree;
    }
}
 
void divide( const Polygon &a, const Polygon &b, Segments &APlus, Segments &AMinus, Segments &BPlus, Segments &BMinus)
{
 
    BSPNode* treeA = buildTree(a.segments);
    BSPNode* treeB = buildTree(b.segments);
    splitSegments(treeA, b.segments, BPlus, BMinus);
    splitSegments(treeB, a.segments, APlus, AMinus);
    deleteTree(treeA);
    deleteTree(treeB);
}
 
Polygon polygonIntersection(const Polygon &a, const Polygon &b)
{
    Segments APlus, AMinus;
    Segments BPlus, BMinus;
    divide(a, b, APlus, AMinus, BPlus, BMinus);
    AMinus.insert(AMinus.end(), BMinus.begin(), BMinus.end());
    return Polygon(AMinus);
}
 
Polygon polygonUnion(const Polygon &a, const Polygon &b)
{
    Segments APlus, AMinus;
    Segments BPlus, BMinus;
    divide(a, b, APlus, AMinus, BPlus, BMinus);
    APlus.insert(APlus.end(), BPlus.begin(), BPlus.end());
    return Polygon(APlus);
}
 
Polygon polygonDifference(const Polygon &a, const Polygon &b)
{
    Segments APlus, AMinus;
    Segments BPlus, BMinus;
    divide(a, b, APlus, AMinus, BPlus, BMinus);
    for (size_t i = 0; i < BMinus.size(); i++)
        BMinus[i] = Segment(BMinus[i].end(), BMinus[i].begin()); // изменение направления обхода, обращение нормалей!
    APlus.insert(APlus.end(), BMinus.begin(), BMinus.end());
    return Polygon(APlus);
}

Сообщение отредактировано: albom - 08.01.08, 14:17

Большое спасибо, albom! Написала я алгоритм, который строит BSP-дерево для одного многоугольника. И поняла... Что он работает только на выпуклых многоугольниках... А у меня они могут быть и вогнутыми...

Плохо, что только для выпуклых. К тому же для выпуклых многоугольников bsb-дерево вырождается в линейный список, и там в нем не много смысла.

Тогда показывайте ваш алгоритм построения, будем разбиратся.

Алгоритм в принципе такой же: то есть ты хочешь сказать, что он работает и для вогнутых? Мне казалось, что именно для одного многоугольника дерево ДОЛЖНО вырождаться в линейный список.
То есть ты хочешь сказать, что когда мы разбиваем этим деревом другой многоугольник все будет правильно работать?