Pakowanie wielokątów w obrębie wielokąta za pomocą ArcGIS Desktop?

Mam logiczny raster.

W szarych obszarach rastra chciałbym dopasować wielokąt o danym rozmiarze w ciągłym zakresie.

Zasadniczo mam nieregularny wielokąt i chciałbym „dopasować” znany wielokąt do zakresu nieregularnego wielokąta tak wiele razy, jak to możliwe.

Kierunek wielokąta nie ma znaczenia i może być kwadratem. Chciałbym, aby pasował graficznie, ale jeśli po prostu przyczepiłby liczbę do wielokąta (# to dopasowanie), to też by działało.

Korzystam z ArcGIS Desktop 10.

Istnieje wiele sposobów rozwiązania tego problemu. Format danych rastrowych sugeruje podejście oparte na rastrze; w przeglądzie tych podejść sformułowanie problemu jako binarnego programu liniowego z liczbami całkowitymi wygląda obiecująco, ponieważ jest bardzo zgodne z duchem wielu analiz wyboru miejsca w GIS i można je do nich łatwo dostosować.

W tym sformułowaniu wyliczamy wszystkie możliwe położenia i orientacje wielokątów wypełniających, które będę określał jako „płytki”. Z każdym kafelkiem związana jest miara „dobroci”. Celem jest znalezienie kolekcji nie nakładających się płytek, których całkowita dobroć jest tak duża, jak to możliwe. Tutaj możemy wziąć dobroć każdej płytki jako obszar, który pokrywa. (W bardziej bogatych w dane i wyrafinowanych środowiskach decyzyjnych możemy obliczać dobroć jako kombinację właściwości komórek zawartych w każdym kafelku, właściwości być może związanych z widocznością, bliskością innych rzeczy itp.)

Ograniczeniami tego problemu są po prostu to, że żadne dwa kafelki w rozwiązaniu nie mogą się pokrywać.

To może być oprawione trochę bardziej abstrakcyjnie, w sposób sprzyjający efektywnej obliczeń, przez wyliczanie komórek w wieloboku być wypełnione ( „Region”) 1, 2, ..., M . Każde umieszczenie kafelka można zakodować za pomocą wektora wskaźników zer i jedynek, dzięki czemu te odpowiadają komórkom pokrytym kafelkiem i zerom w innym miejscu. W tym kodowaniu wszystkie informacje dotyczące kolekcji kafelków można znaleźć, sumując ich wektory wskaźnikowe (jak zwykle komponent po komponencie): suma będzie niezerowa dokładnie tam, gdzie co najmniej jedna płytka pokrywa komórkę, a suma będzie większa niż gdziekolwiek dwa lub więcej płytek zachodzą na siebie. (Suma skutecznie liczy stopień nakładania się).

Jeszcze jeden mały abstrakcji: zbiór możliwych praktyk płytek może być sam wymienił, powiedzmy 1, 2, ..., N . Wybór dowolnego zestawu rozmieszczeń płytek odpowiada wektorowi wskaźnika, w którym te wyznaczają płytki do umieszczenia.

Oto mała ilustracja, aby naprawić pomysły . Towarzyszy mu kod Mathematica użyty do obliczeń, aby trudności programistyczne (lub ich brak) były oczywiste.

Najpierw przedstawiamy region do kafelkowania:

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

Jeśli policzymy jego komórki od lewej do prawej, zaczynając od góry, wektor wskaźnika dla regionu ma 16 wpisów:

Flatten[region]

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

Użyjmy następującego kafelka wraz ze wszystkimi obrotami o wielokrotność 90 stopni:

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

Kod do generowania rotacji (i odbić):

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

(To nieco nieprzejrzyste obliczenie wyjaśniono w odpowiedzi na stronie /math//a/159159 , która pokazuje, że po prostu wytwarza on wszystkie możliwe obroty i odbicia płytki, a następnie usuwa wszelkie duplikaty wyników.)

Załóżmy, że powinniśmy umieścić płytkę, jak pokazano tutaj:

Komórki 3, 6 i 7 są objęte tym miejscem docelowym. Jest to oznaczone przez wektor wskaźnika

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

Gdybyśmy przesunęli ten kafelek o jedną kolumnę w prawo, ten wektor wskaźnika byłby zamiast tego

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

Połączenie prób umieszczenia płytek w obu tych pozycjach jednocześnie zależy od sumy tych wskaźników,

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

2 w siódmej pozycji pokazuje, że nakładają się na siebie w jednej komórce (drugi rząd w dół, trzecia kolumna od lewej). Ponieważ nie chcemy nakładać się, wymagamy, aby suma wektorów w dowolnym prawidłowym rozwiązaniu nie mogła zawierać wpisów przekraczających 1.

Okazuje się, że w przypadku tego problemu możliwe jest 29 kombinacji orientacji i położenia płytek. (Stwierdzono to za pomocą prostego kodowania wymagającego wyczerpującego wyszukiwania.) Wszystkie 29 możliwości możemy przedstawić, rysując ich wskaźniki jako wektory kolumnowe . (Używanie kolumn zamiast wierszy jest konwencjonalne.) Oto obraz wynikowej tablicy, która będzie miała 16 wierszy (po jednym dla każdej możliwej komórki w prostokącie) i 29 kolumn:

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

(Poprzednie dwa wektory wskaźnikowe pojawiają się jako pierwsze dwie kolumny po lewej stronie). Czytelnik o dużych oczach zauważył kilka możliwości równoległego przetwarzania: obliczenia te mogą potrwać kilka sekund.

Wszystkie powyższe można przekształcić w kompaktowy sposób za pomocą notacji macierzowej:

• F to tablica opcji z M wierszami i N kolumnami.

• X jest wskaźnikiem miejscach docelowych do płytek, o długości N .

• b jest wektorem N tych.

• R jest wskaźnikiem regionu; jest to wektor M.

Całkowita „dobroć” związana z dowolnym możliwym rozwiązaniem X jest równa RFX , ponieważ FX jest wskaźnikiem komórek objętych X, a iloczyn z R sumuje te wartości. (Moglibyśmy zważyć R, jeśli chcielibyśmy, aby rozwiązania sprzyjały lub unikały pewnych obszarów w regionie.) Należy to zmaksymalizować. Ponieważ możemy to zapisać jako ( RF ). X , jest to funkcja liniowa X : to ważne. (W poniższym kodzie zmienna czawiera RF .)

Ograniczenia są takie

1. Wszystkie elementy X muszą być nieujemne;

2. Wszystkie elementy X muszą być mniejsze niż 1 (co odpowiada odpowiedniej pozycji b );

3. Wszystkie elementy X muszą być integralne.

Ograniczenia (1) i (2) sprawiają, że jest to program liniowy , podczas gdy trzeci wymóg przekształca go w program liniowy z liczbami całkowitymi .

Istnieje wiele pakietów do rozwiązywania liczb całkowitych programów liniowych wyrażonych dokładnie w tej formie. Są w stanie przetwarzać wartości M i N na dziesiątki, a nawet setki tysięcy. To prawdopodobnie wystarcza do niektórych rzeczywistych aplikacji.

Jako naszą pierwszą ilustrację obliczyłem rozwiązanie dla poprzedniego przykładu za pomocą polecenia Mathematica 8 LinearProgramming. (Spowoduje to zminimalizowanie liniowej funkcji celu. Minimalizację łatwo można zmaksymalizować, negując funkcję celu.) Zwróciło rozwiązanie (jako listę płytek i ich pozycji) w 0,011 sekundy:

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

Szare komórki w ogóle nie są w regionie; białe komórki nie były objęte tym roztworem.

Możesz wypracować (ręcznie) wiele innych rodzajów płytek, które są tak samo dobre jak ten - ale nie możesz znaleźć lepszych. To potencjalne ograniczenie tego podejścia: daje jedno najlepsze rozwiązanie, nawet jeśli jest więcej niż jedno. (Istnieją pewne obejścia: jeśli zmienimy kolejność kolumn X , problem pozostanie niezmieniony, ale w wyniku tego oprogramowanie często wybiera inne rozwiązanie. Jednak takie zachowanie jest nieprzewidywalne).

Jako drugą ilustrację , aby być bardziej realistycznym, rozważmy region w pytaniu. Importując obraz i próbkując go, reprezentowałem go siatką 69 na 81:

Region obejmuje 2156 komórek tej siatki.

Aby uczynić rzeczy interesującymi i zilustrować ogólną konfigurację programowania liniowego, spróbujmy objąć jak najwięcej tego regionu dwoma rodzajami prostokątów:

Jeden ma wymiary 17 na 9 (153 komórki), a drugi ma wymiary 15 na 11 (165 komórek). Wolimy używać drugiego, ponieważ jest większy, ale pierwszy jest bardziej chudy i może zmieścić się w ciasniejszych miejscach. Zobaczmy!

Program obejmuje teraz N = 5589 możliwych miejsc na płytki. Jest dość duży! Po 6,3 sekundach obliczeń Mathematica wymyśliła to rozwiązanie z dziesięcioma płytkami:

Z powodu niektórych luzów ( np. Moglibyśmy przesunąć lewą dolną płytkę do czterech kolumn w lewo), istnieją oczywiście inne rozwiązania nieco różniące się od tego.

Link do O genetycznych algorytmach upakowania wielokątów , podany w mojej odpowiedzi na podobne pytanie w algorytmie Poszukiwanie, aby umieścić maksymalną liczbę punktów w ograniczonym obszarze w minimalnym odstępie? , może być przydatne. Wygląda na to, że metodę można uogólnić do pracy z dowolnymi kształtami kontenerów (a nie tylko prostokątami).

W przypadku mocno ograniczonego podzbioru, o którym wspomniałeś (kwadratowe / trójkątne kafelki w wyboju), zakładając wyraźne optymalizacje powyżej, ten pseudokod powinien uzyskać przybliżoną odpowiedź, po prostu przeprowadzając przez możliwości z wysoką rozdzielczością, brutalnie zmuszając problem. Nie działa poprawnie w sytuacjach, gdy obrót pojedynczych płytek może przynieść korzyści, takie jak płytki prostokątne lub bardzo nieregularny pojemnik. To 1 milion iteracji, w razie potrzeby możesz wypróbować więcej.

Załóżmy kwadrat o bokach długości L.

Utwórz wzór kwadratów w szachownicę, który jest co najmniej o wymiarach zasięgu pojemnika, plus co najmniej 1 litr w każdym kierunku.

N = 0

DX = 0

DY = 0

DR = 0

Zresetuj pozycję szachownicy do pierwotnego środka ciężkości

Dla (R = 1: 100)

Dla (Y = 1: 100)

Dla (X = 1: 100)

M = Policz liczbę kwadratów całkowicie w pojemniku

Jeśli (M> N)

DR = R

DY = Y

DX = X

N = M

Przesuń szachownicę na wschód o L / 100

Zresetuj wschód szachownicy

Przesuń szachownicę na północ o L / 100

Zresetuj północ szachownicy

Obróć szachownicę o 3,6 stopnia CW wokół środka ciężkości

DY = DY * L

DX = DX * L

Zresetuj szachownicę do pierwotnej pozycji i obrotu

Print DR & ”,„ & DX & ”oraz„ & DY & ”to ostateczna macierz translacji / rotacji”

Obróć szachownicę o DR

Przetłumacz szachownicę przez DX, DY

Wybierz kwadraty, które są całkowicie w kontenerze

Eksportuj kwadraty