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

25

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

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

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

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

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

Тад
источник
8
Это очень сложная проблема. Например, требуется много работы, чтобы вписать как можно больше кругов в квадрат. Когда оригинальный многоугольник сложен - как на рисунке - вам нужны мощные процедуры оптимизации. Лучший метод, который я нашел для этой проблемы - это имитационный отжиг, но он не будет доступен в ArcGIS, и для его написания потребуются очень хитрые сценарии (ArcGIS работает слишком медленно). Не могли бы вы, возможно, немного ослабить свои требования, такие как подгонка меньшего многоугольника достаточное количество раз, а не столько раз, сколько возможно?
whuber
1
@whuber Спасибо за редактирование моего поста. Да, достаточное количество раз будет работать. Или как насчет заданной угловой ориентации. ех. на изображении выше, я подгонял многоугольник столько раз, сколько мог бы иметь в этой ориентации, если бы я повернул их на 90 градусов, ты мог бы соответствовать еще одному ...
Thad
1
Да, но это также чревато подводными камнями. Некоторые из них элементарные. Например, авторский и опубликованный ESRI текст «Знакомство с ArcView GIS» (для версии 3) включал в себя упражнение, в котором прямоугольник, представляющий футбольное поле, был интерактивно размещен внутри многоугольника. Проблема заключалась в том, что ответ упражнения был неверным, поскольку автору не удалось спроецировать данные, а ошибки в использовании географических координат были достаточно велики, чтобы повлиять на результат. Ответ выглядел хорошо в ГИС, но если бы кто-нибудь попытался построить это поле, он бы обнаружил, что для этого недостаточно места :-).
whuber
6
@whuber Полагаю, они думали, что фигуры в «парковой зоне» достаточно
Кирк Куйкендалл
2
В общем случае нерегулярного многоугольника внутри нерегулярного многоугольника это вычислительно неразрешимая проблема: поиск оптимального решения не является вероятной целью во всех случаях, и он, вероятно, является NP-полным с технической точки зрения: в каких случаях это невозможно предопределить. Если вы существенно ограничите задачу, некоторые итерационные алгоритмы случайного подбора, вероятно, дадут вам достаточно большие числа. Мне кажется, что это задание, что они не ищут правильный ответ, они ищут творческие подходы.
MappingT завтра

Ответы:

22

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

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

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

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

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

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

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

region =  {{0, 0, 1, 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}}};

Рисунок 2: плитка

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

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: размещение плитки

Клетки 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];

Рисунок 4: массив опций

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

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

  • 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 &]]];

Рисунок 5: решение

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

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

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

Рисунок 6: Регион

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

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

Рисунок 7: плитки

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

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

Рисунок 8: решение

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

Whuber
источник
1
Более ранняя версия этого решения (но не такая хорошая) появляется на сайте Mathematica по адресу mathematica.stackexchange.com/a/6888 . Может также стоить отметить, что незначительное изменение формулировки может быть использовано для решения проблемы полного покрытия области как можно меньшим количеством плиток (допуская некоторые наложения, конечно): это решило бы «исправление выбоин» проблема.
whuber
1
В интересах космоса этот ответ не описывает некоторые потенциально полезные улучшения. Например, после нахождения всех возможных положений плитки (как векторов индикатора), вы можете добавить их все, чтобы найти, какие ячейки действительно могут быть покрыты какой-либо плиткой. Множество таких ячеек разбивается на два отдельных связанных компонента во втором примере. Это означает, что проблема может быть решена независимо в двух компонентах, существенно уменьшая ее размер (и, следовательно, вычислительное время). Такие первоначальные упрощения имеют тенденцию быть важными для решения реальных проблем.
whuber
Большое усилие и ответ. Ответ Криса также был полезен. Спасибо всем за помощь! Работает, и снова заставил меня двигаться в правильном направлении.
Чт
Вот Это Да! Меня заинтересовала похожая проблема, и этот пост дал мне новую перспективу. Спасибо. Что если R больше (например, 140x140≈20000), есть ли способы уменьшить стоимость вычислений? Знаете ли вы какие-либо документы, связанные с этой проблемой? Мои ключевые слова для поиска не ведут меня в правильном направлении (до сих пор).
Nimcap
@nimcap Это важный класс проблем, поэтому многие исследования продолжаются. Ключевые слова для поиска начинаются с "смешанной целочисленной линейной программы" и разветвляются оттуда в зависимости от того, что вы найдете.
whuber
5

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

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

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

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

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

MappingTomorrow
источник
Если вы попробуете эту процедуру в области 2 на 5 с ячейкой, отсутствующей по середине одного длинного края, вы обнаружите, что можете поместить в нее только один квадрат 2 на 2. Однако два таких квадрата легко вписываются. Проблема заключается в том, что они не являются частью обычного шаблона «шахматной доски». Эта трудность является одной из вещей, которая делает эту проблему довольно сложной.
whuber
1
Ага. Если у вас есть форма контейнера, достаточно нерегулярная, чтобы он мог поддерживать несколько разрозненных регулярных шаблонов порядка нескольких ячеек в каждой, это очень далеко от оптимального. Добавление подобных вещей в пространство возможностей очень быстро увеличивает время обработки и требует определенной степени планирования для конкретного случая, на который вы нацелены.
MappingT завтра