Подсчитать число пересечений прямоугольников (кубов)

#массивы #алгоритм #сортировка #геометрия #вычислительная_геометрия

Для простоты будем рассматривать задачу про прямоугольники (для кубов -- аналогично).

Задача состоит в том, чтобы подсчитать суммарную площадь пересечений прямоугольников
(или любую другую аддитивную функцию от площади пересечений), которые ориентированы
по осям стандартной системы координат (СК). Кажется, что задача решается быстрее чем
брутфорсом, т.е. быстрее, чем за O(n2). 

Решение можно построить так:


Переберём все упорядоченные пары прямоугольников (их C2n).
Для каждой пары посчитаем площадь пересечений.


Предложим решение задачи для случая двух прямоугольников. 



Будем обозначать проекции точек на оси как Ax, Bx, Cx, Dx, Ay, By... и т.д.

Тогда, запишем все проекции и отсортируем их согласно возрастанию значению координат.
Чтобы не дублировать координаты, для оси X будем брать только точки A, B, E, F:

Ax, Ex, Bx, Fx

Для оси y соответственно:

Hy, Dy, Ey, Ay

Заведём сет и будем поочередно класть точки в него. Если одновременно в сете есть
две точки одного и того же прямоугольника, то удалим их. Если одновременно в сете есть
точки различных прямоугольников по каждой из координат, то фиксируем пересечение.

Пример. 


Sx = {} // Сет Sx -- по оси х
Sy = {} // Сет Sy -- по оси y
Добавим в сет первые точки по каждой из координат: 
Sx = {Ax}, Sy = {Hy}
Добавим ещё по одной точке:
Sx = {Ax, Ex}, Sy = {Hy, Dy}. Фиксируем пересечение
Sx = {Ax, Ex, Bx}, Sy = {Hy, Dy, Ey}. Вычисляем значение пересечения:
(Bx - Ax) * (Hy - Ey)
Удаляем точки из сетов: Sx = {Ex}, Sy = {Dy}
Добавляем точки в сеты: Sx = {Ex, Fx}, Sy = {Dy, Ay}
Удаляем точки одинаковых прямоугольников: Sx = {}, Sy = {}


Вопросы: 


Будет ли работать подобный алгоритм для 3 измерений?
Помогите обобщить его на случай n прямоугольнков (кубов). Никак не соображу.


Для случая двойных пересечений требует вычислить суммарную площадь. Пусть есть 3
прямоугольника ABCD, EFGH, IJKL. Требуется вычислить либо площадь S(INMO) + S(EOPQ)
+ S(OSKL), либо S(INMO) + 2S(EOPQ) + S(OSKL), в зависимости от того, что проще посчитать.



Замечание

Для случая двух прямоугольников (кубов) задача решается легко. Пусть у нас заданы
прямоугольники как левые нижние углы и размеры их сторон:

x0, y0, z0, dx, dy, dz


Тогда посмотрим, какой из прямоугольников лежит правее по каждой оси. Затем проверим
условия пересечений и вычислим площадь. Напишем код (на golang):

type Box struct {
    dx, dy, dz int
    x, y, z    int
}

func (c *Box) isLineIntersection(p00, p01, p10, p11 int) int {
    if p00 > p10 {
        p00, p01, p10, p11 = p10, p11, p00, p01
    }
    if (p10 >= p01) && (p01 <= p11) {
        return p01 - p10
    } else if (p10 <= p01) && (p11 <= p01) {
        return p11 - p10
    }
    return 0
}


func (c *Box) isIntersection(box1 Box, box2 Box) (intersection int) {
    dx := c.isLineIntersection(box1.x, box1.x+box1.dx, box2.x, box2.x+box2.dx)
    dy := c.isLineIntersection(box1.y, box1.y+box1.dy, box2.y, box2.y+box2.dy)
    dz := c.isLineIntersection(box1.z, box1.z+box1.dz, box2.z, box2.z+box2.dz)
    return dx * dy * dz
}

    

Ответы

Ответ 1

Такие задачи канонически решаются классическим алгоритмом Сканирующей Прямой. Согласно
вашей постановке, вам нужно вычислить площадь "территории" покрытой более чем одним
прямоугольником.

В "простейшем" применении этого подхода: 

Мы представляем каждый из наших прямоугольников парой вертикальных ребер: левым и
правым ребром. Эти ребра мы будем рассматривать в лексикографическом порядке, отсортированными
по нижнему Y, а затем - по X.

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

Рассматривать необходимо только те уровни сканирования, на которых начинается или
заканчивается какое-то вертикальное ребро. Т.е. наша сканирующая прямая движется снизу-вверх
прыжками по концевым точкам наших вертикальных рёбер.

Для каждого уровня сканирования мы анализируем покрытие плоскости непосредственно
над сканирующей прямой.

Просматривая список активных ребер слева-направо, мы легко узнаем, на каком вертикальном
ребре X1 начинается множественное покрытие "территории" прямоугольниками и на каком
ребре X2 оно заканчивается. Таких пар начало-конец на каждом уровне сканирования может
быть несколько. Если следующая сканирующая прямая будет располагаться на уровне NEXT_Y,
то каждая пара X1-X2 будет добавлять в искомую площадь слагаемое (X2 - X1) * (NEXT_Y - Y). 

Просканировав снизу-вверх весь вход мы вычислим искомую площадь.



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

Вот, например, простейшая реализация этого алгоритма (С++) с входными данными взятыми
из вашей картинки с тремя прямоугольниками. Получаем именно S(INMO) + S(EOPQ) + S(OSKL),
т.е. 22

#include 
#include 
#include 
#include 

struct Edge
{
  int x, yb, yt;
  bool is_left;
};

struct Rect
{
  int xl, yb, xr, yt;
};

int main()
{
  const std::vector rects = 
  {
    { 0, 6, 5, 11 },
    { 3, 4, 7, 10 },
    { 2, 0, 9, 8 }
  };

  using Edges = std::vector; 
  Edges edges;
  for (const Rect &r : rects)
  {
    edges.push_back({ r.xl, r.yb, r.yt, true });
    edges.push_back({ r.xr, r.yb, r.yt, false });
  }

  std::sort(edges.begin(), edges.end(), 
    [](const Edge &lhs, const Edge &rhs) { return lhs.yb != rhs.yb ? lhs.yb < rhs.yb
: lhs.x < rhs.x; });

  unsigned area = 0;

  auto less_x = [](const Edge &lhs, const Edge &rhs) { return lhs.x < rhs.x; };
  using SweepLine = std::set;
  SweepLine sweep_line(less_x);

  Edges::const_iterator it_e = edges.begin();
  for (int sweep_y = it_e->yb, next_sweep_y; sweep_y != std::numeric_limits::max();
sweep_y = next_sweep_y)
  {
    for (; it_e != edges.end() && it_e->yb == sweep_y; ++it_e)  
      sweep_line.insert(*it_e);

    next_sweep_y = it_e != edges.end() ? it_e->yb : std::numeric_limits::max();

    int prev_x;
    unsigned inside_counter = 0;
    unsigned covered_x = 0;

    for (SweepLine::iterator it_swe = sweep_line.begin(); it_swe != sweep_line.end(); )
    {
      if (it_swe->yt == sweep_y)
      {
        it_swe = sweep_line.erase(it_swe);
        continue;
      }

      next_sweep_y = std::min(next_sweep_y, it_swe->yt);

      unsigned prev_inside_counter = inside_counter;
      inside_counter += it_swe->is_left ? +1 : -1;

      if (prev_inside_counter == 1 && inside_counter == 2)
        prev_x = it_swe->x;
      else if (prev_inside_counter == 2 && inside_counter == 1)
        covered_x += it_swe->x - prev_x;

      ++it_swe;
    }

    area += covered_x * (next_sweep_y - sweep_y);
  }

  std::cout << area << std::endl;
}


http://coliru.stacked-crooked.com/a/bb0023433d9395be

Элементарной модификацией условия во внутреннем if можно заставить эту реализацию
вычислять площадь объединения прямоугольников

if (prev_inside_counter == 0 && inside_counter == 1)
  prev_x = it_swe->x;
else if (prev_inside_counter == 1 && inside_counter == 0)
  covered_x += it_swe->x - prev_x;


площадь как минимум тройного пересечения 

if (prev_inside_counter == 2 && inside_counter == 3)
  prev_x = it_swe->x;
else if (prev_inside_counter == 3 && inside_counter == 2)
  covered_x += it_swe->x - prev_x;


площадь строго двойного пересечения

if (prev_inside_counter != 2 && inside_counter == 2)
  prev_x = it_swe->x;
else if (prev_inside_counter == 2 && inside_counter != 2)
  covered_x += it_swe->x - prev_x;


площадь, покрытая нечетным количеством прямоугольников

if ((prev_inside_counter % 2) == 0 && (inside_counter % 2) != 0)
  prev_x = it_swe->x;
else if ((prev_inside_counter % 2) != 0 && (inside_counter % 2) == 0)
  covered_x += it_swe->x - prev_x;


и т.д. и т.п.

Я выше назвал этот алгоритм "простейшим" потому, что он не умеет проводить локализованную
обработку каждого уровня сканирования: на каждом уровне все активные ребра этого уровня
просматриваются слева-направо от начала до конца. Грамотная реализация такого алгоритма
должна уметь вместо этого производить локализованную обработку: только в некоторой
окрестности тех мест в уровне сканирования, где появились новые ребра или исчезли старые.
Но это уже будет существенно менее тривиальный алгоритм, хотя общая идея останется
неизменной.


Ответ 2

Вижу так. Суммируем общую площадь. По - очереди добавляем объект к сумме. Сумма это
список обыкновенных прямоугольников (кубов). Например сумма это пока один прямоугольник.
К нему добавляется ещё один. Следующая сумма это будет ТРИ прямоугольника. (Ⅰ + Ⅱ + Ⅲ).


Дальше берём из списка ещё один (третий). И пересекаем с этими тремя. В сложном случае
в сумме будет уже много прямоугольников. К линейности никак задача не сводится. Думайте,
господа присяжные заседатели....


Кубический вариант.
Берётся список координат X всех кубов. Записывается в список X-ов (LZ). Все Y в список
Y-ов(LY). Все Z координаты в список Z-ов(LZ). Сортируются эти три списка (X,Y,Z). Создаётся
кубовый массив с элементами int размера Size(LX)xSize(LY)xSize(LZ). Заполняется пока
0. Это будем называть индексированный куб. Где элемент (i,j,k) это признак заполнености
по координатам (LX(i)..LX(i+1),LY(j)..LY(j+1),LZ(k)..LZ(k+1)). Каждый 3D-прямоугольник
должен добавить своё присутствие к индексированному кубу добавлением единицы елементам
где он свой след оставил. Общий объём считаем суммой всех объёмов индексированного
куба, где элемент больше одного (т.е. где было хотя-бы два перекрытия).