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

25

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.

Thad
źródło
8
To bardzo trudny problem. Na przykład dopasowanie wielu kwadratów do kwadratu wymaga dużo pracy. Kiedy oryginalny wielokąt jest skomplikowany - jak na ilustracji - potrzebujesz zaawansowanych procedur optymalizacji. Najlepszą metodą, jaką znalazłem dla tego problemu, jest symulowane wyżarzanie, ale nie będzie ono dostępne w ArcGIS, a wykonanie skryptu wymagałoby niezwykle przebiegłego skryptu (ArcGIS jest zbyt wolny). Czy mógłbyś nieco rozluźnić swoje wymagania, na przykład dopasować mniejszy wielokąt wystarczającą liczbę razy, a nie tyle razy, ile to możliwe?
whuber
1
@whuber Dziękujemy za edycję mojego postu. Tak, wystarczająca liczba razy by zadziałała. Albo co powiesz na pod kątem orientacji. dawny. na powyższym obrazku dopasowałem wielokąt tyle razy, ile mogłem w tej orientacji, gdybym obrócił go o 90 stopni, możesz zmieścić jeszcze jeden ...
Thad
1
Tak, ale wiąże się to również z pułapkami. Niektóre są elementarne. Na przykład tekst autorstwa ESRI „Poznajemy ArcView GIS” (dla wersji 3) zawiera ćwiczenie, w którym prostokąt reprezentujący boisko do piłki nożnej został umieszczony interaktywnie w wielokącie. Problem polegał na tym, że odpowiedź ćwiczenia była błędna, ponieważ autorowi nie udało się wyświetlić danych, a błędy w użyciu współrzędnych geograficznych były wystarczająco duże, aby wpłynąć na wynik. Odpowiedź wyglądała dobrze w GIS, ale gdyby ktoś próbował zbudować to pole, stwierdziłby, że nie ma na to miejsca :-).
whuber
6
@ whuber Myślę, że myśleli, że „park piłkarski” jest wystarczający.
Kirk Kuykendall
2
W ogólnym przypadku nieregularnego wielokąta w obrębie nieregularnego wielokąta jest to problem trudny do obliczenia: znalezienie optymalnego rozwiązania nie jest prawdopodobnym celem we wszystkich przypadkach i prawdopodobnie jest kompletne z technicznego punktu widzenia: w jakich przypadkach nie można tego z góry ustalić. Jeśli znacznie ograniczysz problem, niektóre iteracyjne algorytmy losowego dopasowania mogą dać ci dość wysokie liczby. Mam wrażenie, że jeśli jest to zadanie, to nie szukają poprawnej odpowiedzi, szukają kreatywnego podejścia.
Mapowanie jutro

Odpowiedzi:

22

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

Ryc. 1: region

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

Rysunek 2: kafelek

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:

Rysunek 3: umieszczenie kafelków

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

Rysunek 4: tablica opcji

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

Rysunek 5: rozwiązanie

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:

Rysunek 6: Region

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:

Rysunek 7: płytki

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:

Rysunek 8: rozwiązanie

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.

Whuber
źródło
1
Wcześniejsza wersja tego rozwiązania (ale nie tak dobra) pojawia się na stronie Mathematica pod adresem mathematica.stackexchange.com/a/6888 . Warto również zauważyć, że można zastosować niewielką odmianę preparatu, aby rozwiązać problem całkowitego pokrycia regionu możliwie jak najmniejszą liczbą płytek (oczywiście dopuszczając pewne nakładanie się): rozwiązałoby to „łatanie dziur” problem.
whuber
1
W interesie przestrzeni kosmicznej odpowiedź ta nie opisuje niektórych potencjalnie pomocnych ulepszeń. Na przykład po znalezieniu wszystkich możliwych pozycji kafelków (jako wektorów wskaźnikowych) możesz dodać je wszystkie, aby dowiedzieć się, które komórki mogą być faktycznie pokryte niektórymi kafelkami. Zestaw takich komórek dzieli się na dwa osobne połączone komponenty w drugim przykładzie. Oznacza to, że problem można niezależnie rozwiązać w dwóch komponentach, znacznie zmniejszając jego rozmiar (a tym samym czas obliczeń). Takie wstępne uproszczenia wydają się być ważne przy rozwiązywaniu rzeczywistych problemów.
whuber
Wielki wysiłek i odpowiedź. Pomoc Chrisa również była pomocna. Dzięki wszystkim za pomoc! Działa i znów popycha mnie w dobrym kierunku.
Thad
Łał! Byłem zainteresowany podobnym problemem i ten post dał mi nową perspektywę. Dziękuję Ci. Co jeśli R jest większy (np. 140x140≈20000), czy są jakieś sposoby na zmniejszenie kosztów obliczeń? Czy znasz jakieś dokumenty związane z tym problemem? Moje słowa kluczowe wyszukiwania nie prowadzą mnie we właściwy sposób (do tej pory).
nimcap
@nimcap To ważna klasa problemów, toczy się więc wiele badań. Szukane słowa kluczowe zaczynają się od „mieszanego programu liniowego z liczbami całkowitymi” i rozgałęziają się tam na podstawie tego, co znajdziesz.
whuber
5

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).

Kirk Kuykendall
źródło
Ten papier ma kilka fajnych pomysłów (+1), ale wszystkie jego algorytmy koncentrują się w zasadniczy sposób na upakowaniu wielokątów w prostokątnych obszarach. Wynika to z tego, że reprezentuje on pakiety o dyskretnej strukturze danych (sekwencja wielokątów wraz z ich orientacjami), która reprezentuje zestaw procedur, w których wielokąty są przesuwane , równolegle do boków kwadratu, w kierunku wyznaczonego rogu. Wydaje się, że takie proste dyskretne kodowanie byłoby mniej skuteczne w bardziej skomplikowanych regionach. Może pomogłoby początkowe uproszczenie regionów w siatce.
whuber
2

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

Mapowanie jutro
źródło
Jeśli spróbujesz wykonać tę procedurę w regionie 2 na 5, w którym brakuje komórki wzdłuż środka jednej długiej krawędzi, przekonasz się, że możesz umieścić w niej tylko jeden kwadrat 2 na 2. Jednak dwa takie kwadraty łatwo pasują. Problem polega na tym, że nie są one częścią zwykłego wzoru „szachownicy”. Ta trudność jest jedną z rzeczy, która sprawia, że ​​problem jest dość trudny.
whuber
1
Tak. Jeśli masz kształt pojemnika na tyle nieregularny, że może on obsługiwać wiele niejednoznacznych regularnych wzorów w kolejności po kilka komórek, to kończy się to bardzo dalekie od optymalnego. Dodanie takich rzeczy do możliwości bardzo szybko wydłuża czas przetwarzania i wymaga pewnego stopnia planowania w konkretnym przypadku, na który celujesz.
Mapowanie jutro