Упаковывать полигоны внутри полигона с помощью ArcGIS Desktop?

У меня есть логический растр.

В серых областях растра я хотел бы разместить полигон заданного размера в смежных пределах.

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

Направление многоугольника не имеет значения, и это может быть квадрат. Я хотел бы, чтобы он соответствовал графически, но если бы он просто прикрепил число к многоугольнику (#, который подходит), это также сработало бы.

Я использую ArcGIS Desktop 10.

Есть много способов решить эту проблему. Растровый формат данных предполагает растровый подход; При рассмотрении этих подходов формулировка проблемы в виде двоичной целочисленной линейной программы выглядит многообещающей, поскольку она во многом соответствует духу многих анализов выбора площадок ГИС и может быть легко адаптирована к ним.

В этой формулировке мы перечисляем все возможные положения и ориентации многоугольника (ов) заполнения, которые я буду называть «плитками». С каждой плиткой связана мера ее "благости". Цель состоит в том, чтобы найти коллекцию неперекрывающихся плиток, общее качество которых максимально велико. Здесь мы можем принять качество каждой плитки за область, которую она покрывает. (В более насыщенных данными и сложных средах принятия решений мы можем вычислять качество как комбинацию свойств ячеек, включенных в каждую плитку, свойств, возможно, связанных с видимостью, близостью к другим вещам и т. Д.)

Ограничения этой проблемы просто в том, что никакие две плитки в решении не могут перекрываться.

Это может быть оформлена немного более абстрактно, таким образом , способствует эффективному вычислению, путем перечисления ячеек на полигоне должны быть заполнены (в «область») 1, 2, ..., M . Любое размещение тайла может быть закодировано с индикатором- вектором нулей и единиц, позволяя тем, чтобы они соответствовали ячейкам, покрытым тайлом и нулями в другом месте. В этом кодировании вся информация, необходимая о наборе плиток, может быть найдена путем суммирования их векторов-индикаторов (компонент за компонентом, как обычно): сумма будет отлична от нуля именно там, где хотя бы одна плитка покрывает ячейку, а сумма будет больше чем где-либо, две или более плитки перекрываются. (Сумма эффективно подсчитывает количество совпадений.)

Еще одна маленькая абстракция: множество возможных мест размещения плитки само можно перечислить, скажем , 1, 2, ..., N . Выбор любого набора размещений плиток сам по себе соответствует вектору-индикатору, где те определяют плитки, которые должны быть размещены.

Вот небольшая иллюстрация, чтобы исправить идеи . Он сопровождается кодом Mathematica, используемым для выполнения расчетов, так что сложность программирования (или ее отсутствие) может быть очевидной.

Во-первых, мы изобразим область, которая будет выложена плиткой:

region =  {{0, 0, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}};

Если мы нумеруем его ячейки слева направо, начиная сверху, вектор индикатора для региона имеет 16 записей:

Flatten[region]

{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}

Давайте использовать следующую плитку вместе со всеми поворотами, кратными 90 градусам:

tileSet = {{{1, 1}, {1, 0}}};

Код для генерации вращений (и отражений):

apply[s_List, alpha] := Reverse /@ s;
apply[s_List, beta] := Transpose[s];
apply[s_List, g_List] := Fold[apply, s, g];
group = FoldList[Append, {}, Riffle[ConstantArray[alpha, 4], beta]];
tiles = Union[Flatten[Outer[apply[#1, #2] &, tileSet, group, 1], 1]];

(Это несколько непрозрачное вычисление объясняется в ответе по адресу /math//a/159159 , который показывает, что он просто производит все возможные повороты и отражения плитки, а затем удаляет любые повторяющиеся результаты.)

Предположим, мы должны были разместить плитку, как показано здесь:

Клетки 3, 6 и 7 включены в это размещение. Это обозначено индикатором вектора

{0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}

Если мы сместим эту плитку на один столбец вправо, этот индикаторный вектор будет вместо

{0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0}

Комбинация попыток разместить плитки в обеих этих позициях одновременно определяется суммой этих показателей,

{0, 0, 1, 1, 0, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0}

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

Оказывается, что для этой задачи 29 комбинаций ориентации и положения возможны для плиток. (Это было найдено с помощью простого фрагмента кода, включающего исчерпывающий поиск.) Мы можем изобразить все 29 возможностей, нарисовав их индикаторы в виде векторов столбцов . (Использование столбцов вместо строк обычно.) Вот изображение результирующего массива, который будет иметь 16 строк (по одной на каждую возможную ячейку в прямоугольнике) и 29 столбцов:

makeAllTiles[tile_, {n_Integer, m_Integer}] := 
  With[{ m0 = Length[tile], n0 = Length[First[tile]]},
   Flatten[
    Table[ArrayPad[tile, {{i, m - m0 - i}, {j, n - n0 - j}}],  {i, 0, m - m0}, {j, 0, n - n0}], 1]];
allTiles = Flatten[ParallelMap[makeAllTiles[#, ImageDimensions[regionImage]] & , tiles], 1];
allTiles = Parallelize[
   Select[allTiles, (regionVector . Flatten[#]) >= (Plus @@ (Flatten[#])) &]];
options = Transpose[Flatten /@ allTiles];

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

Все вышесказанное можно компактно переформулировать, используя матричную запись:

• F - это массив опций с M строками и N столбцами.

• Х является индикатором набора плитки размещений, длины N .

• b является N- вектором единиц.

• R - показатель для региона; это М- вектор.

Общая «доброта», связанная с любым возможным решением X , равна RFX , потому что FX является индикатором ячеек, покрытых X, и произведение с R суммирует эти значения. (Мы могли бы оценить R, если бы хотели, чтобы решения благоприятствовали или избегали определенных областей в регионе.) Это должно быть максимизировано. Потому что мы можем написать это как ( RF ). X , это линейная функция от X : это важно. (В приведенном ниже коде переменная cсодержит RF .)

Ограничения в том, что

1. Все элементы X должны быть неотрицательными;

2. Все элементы X должны быть меньше 1 (что является соответствующей записью в b );

3. Все элементы X должны быть целыми.

Ограничения (1) и (2) делают эту программу линейной , а третье требование превращает ее в целочисленную линейную программу.

Существует множество пакетов для решения целочисленных линейных программ, выраженных именно в этой форме. Они способны обрабатывать значения M и N в десятки или даже сотни тысяч. Это, вероятно, достаточно хорошо для некоторых реальных приложений.

В качестве нашей первой иллюстрации я вычислил решение для предыдущего примера, используя команду Mathematica 8 LinearProgramming. (Это сведет к минимуму линейную целевую функцию. Минимизация легко превращается в максимизацию путем отрицания целевой функции.) Она вернула решение (в виде списка плиток и их позиций) за 0,011 секунды:

b = ConstantArray[-1, Length[options]];
c = -Flatten[region].options;
lu = ConstantArray[{0, 1}, Length[First[options]]];
x = LinearProgramming[c, -options, b, lu, Integers, Tolerance -> 0.05];
If[! ListQ[x] || Max[options.x] > 1, x = {}];
solution = allTiles[[Select[x Range[Length[x]], # > 0 &]]];

Серые клетки вообще не в регионе; белые клетки не были покрыты этим раствором.

Вы можете отработать (вручную) множество других плиток, которые так же хороши, как этот, но вы не можете найти лучшие. Это потенциальное ограничение этого подхода: он дает вам одно лучшее решение, даже если их несколько. (Есть некоторые обходные пути: если мы изменим порядок столбцов X , проблема останется неизменной, но в результате программное обеспечение часто выбирает другое решение. Однако это поведение непредсказуемо.)

В качестве второй иллюстрации , чтобы быть более реалистичным, давайте рассмотрим регион в вопросе. Импортировав изображение и изменив его выборку, я представил его с помощью сетки 69 на 81:

Регион состоит из 2156 ячеек этой сетки.

Чтобы сделать вещи интересными и проиллюстрировать общность настройки линейного программирования, давайте попробуем охватить как можно большую часть этой области двумя видами прямоугольников:

Один - 17 на 9 (153 клетки), а другой - 15 на 11 (165 клеток). Мы могли бы предпочесть использовать второе, потому что оно больше, но первое тоньше и может поместиться в более узких местах. Посмотрим!

Программа теперь включает в себя N = 5589 возможных мест размещения плитки. Это довольно большой! После 6,3 секунд расчета Mathematica предложила следующее решение:

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

Ссылка на Генетические алгоритмы для упаковки многоугольников , приведенная в моем ответе на аналогичный вопрос в поиске алгоритма, чтобы разместить максимальное количество точек в пределах ограниченной области на минимальном расстоянии? может быть полезным. Кажется, что метод может быть обобщен для работы с произвольными формами контейнера (а не только с прямоугольниками).

Для упомянутого вами сильно ограниченного подмножества (квадратная / треугольная мозаика в выбоине), при условии явной оптимизации, приведенной выше, этот псевдокод должен дать приблизительный ответ, просто проведя вас через возможности с высоким разрешением, грубо форсируя проблему. Он не будет работать правильно в ситуациях, когда вращение отдельной плитки может видеть усиление, например, прямоугольные плитки или сильно нерегулярный контейнер. Это 1 миллион итераций, при необходимости вы можете попробовать больше.

Предположим, квадрат со сторонами длины L

Создайте шахматный узор из квадратов, который по крайней мере соответствует размеру экстента контейнера, плюс по крайней мере 1 л в каждом направлении.

N = 0

DX = 0

DY = 0

DR = 0

Сброс позиции шахматной доски до исходного центроида

Для (R = 1: 100)

Для (Y = 1: 100)

Для (X = 1: 100)

M = Подсчет количества квадратов полностью внутри контейнера

Если (М> Н)

DR = R

DY = Y,

DX = Х

N = M

Переместить шахматную доску на восток на L / 100

Сброс шахматной доски на восток

Переместить шахматную доску на север на L / 100

Сброс шахматной доски на север

Поверните шахматную доску на 3,6 градуса по часовой стрелке вокруг ее центроида

DY = DY * L

DX = DX * L

Сбросить шахматную доску в исходное положение и вращение

Печать DR & "," & DX & "и" & DY & "являются окончательной матрицей перевода / поворота"

Повернуть шахматную доску на DR

Перевести шахматную доску по DX, DY

Выберите квадраты, которые полностью в контейнере

Экспортные квадраты