Empacotando polígonos no polígono usando o ArcGIS Desktop?

25

Eu tenho uma varredura booleana.

Nas áreas cinzentas da varredura, eu gostaria de ajustar um polígono de determinado tamanho em uma extensão contígua.

Basicamente, tenho um polígono irregular e gostaria de "encaixar" um polígono conhecido dentro da extensão do polígono irregular o máximo de vezes possível.

A direção do polígono não importa e pode ser um quadrado. Gostaria que ele caísse graficamente, mas se apenas anexasse um número ao polígono (# que couber) também funcionaria.

Estou usando o ArcGIS Desktop 10.

Thad
fonte
8
Este é um problema muito difícil. É preciso muito trabalho para encaixar o máximo de círculos em um quadrado possível, por exemplo. Quando o polígono original é complicado - como na ilustração - você precisa de alguns procedimentos de otimização poderosos. O melhor método que eu encontrei para esse problema é o recozimento simulado, mas isso não estará disponível no ArcGIS e seria necessário um script extremamente astuto para escrevê-lo (o ArcGIS é muito lento). Você poderia talvez relaxar um pouco seus requisitos, como encaixar o polígono menor um número suficiente de vezes, e não quantas vezes possível?
whuber
1
Obrigado por editar meu post. Sim, um número suficiente de vezes funcionaria. Ou então, em uma determinada orientação angular. ex. na imagem acima, eu tenho encaixar o polígono tantas vezes quanto eu poderia ter, pelo que a orientação, se eu rodá-las 90 graus você poderia caber mais um ...
Thad
1
Sim, mas também está cheio de armadilhas. Alguns são elementares. Por exemplo, o texto de autoria e publicação da ESRI, "Conhecendo o ArcView GIS" (para versão 3), incluiu um exercício no qual um retângulo representando um campo de futebol foi colocado interativamente dentro de um polígono. O problema era que a resposta do exercício estava errada porque o autor falhou ao projetar os dados e os erros no uso das coordenadas geográficas foram grandes o suficiente para afetar o resultado. A resposta parecia boa no SIG, mas se alguém tentasse criar esse campo, teria encontrado que não havia espaço suficiente para isso :-).
whuber
6
@whuber Eu acho que eles pensaram que um número "ball park" era suficiente.
Kirk Kuykendall
2
No caso geral de polígono irregular dentro de polígono irregular, esse é um problema computacionalmente intratável: encontrar uma solução ideal não é uma meta plausível em todos os casos e é provável que o NP seja completo de uma perspectiva técnica: quais casos esses não podem ser predeterminados. Se você restringir o problema significativamente, é provável que alguns algoritmos iterativos de ajuste aleatório forneçam números razoavelmente altos . Meu sentimento, se for uma tarefa, é que eles não estão procurando a resposta correta , estão procurando abordagens criativas.
MappingTomorrow

Respostas:

22

Existem muitas maneiras de abordar esse problema. O formato raster dos dados sugere uma abordagem baseada em raster; Ao revisar essas abordagens, uma formulação do problema como um programa linear inteiro binário parece promissora, porque está muito no espírito de muitas análises de seleção de sites GIS e pode ser prontamente adaptada a elas.

Nesta formulação, enumeramos todas as posições e orientações possíveis do (s) polígono (s) de preenchimento, que chamarei de "ladrilhos". Associado a cada bloco está uma medida de sua "bondade". O objetivo é encontrar uma coleção de ladrilhos não sobrepostos cuja bondade total seja a maior possível. Aqui, podemos considerar a qualidade de cada peça como a área que ela cobre. (Em ambientes de decisão mais sofisticados e ricos em dados, podemos calcular a bondade como uma combinação de propriedades das células incluídas em cada bloco, propriedades talvez relacionadas à visibilidade, proximidade a outras coisas, etc.)

As restrições desse problema são simplesmente que não há dois blocos dentro de uma solução sobrepostos.

Isto pode ser enquadrado um pouco mais abstrata, de uma forma favorável à computação eficiente, enumerando as células no polígono a ser preenchido (a "região") 1, 2, ..., M . Qualquer posicionamento de bloco pode ser codificado com um vetor indicador de zeros e uns, permitindo que os correspondam às células cobertas pelo bloco e zeros em outros lugares. Nessa codificação, todas as informações necessárias sobre uma coleção de blocos podem ser encontradas somando seus vetores indicadores (componente por componente, como sempre): a soma será diferente de zero exatamente onde pelo menos um bloco cobre uma célula e a soma será maior do que um em qualquer lugar, dois ou mais blocos se sobrepõem. (A soma conta efetivamente a quantidade de sobreposição.)

Mais um pouco de abstração: o conjunto de possíveis colocações de azulejos em si podem ser enumerados, digamos 1, 2, ..., N . A seleção de qualquer conjunto de veiculações de ladrilhos corresponde a um vetor indicador em que os designam os ladrilhos a serem colocados.

Aqui está uma pequena ilustração para corrigir as idéias . É acompanhado pelo código do Mathematica usado para fazer os cálculos, para que a dificuldade de programação (ou a falta dela) seja evidente.

Primeiro, mostramos uma região a ser ladrilhada:

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

Figura 1: região

Se numerarmos suas células da esquerda para a direita, começando no topo, o vetor indicador para a região terá 16 entradas:

Flatten[region]

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

Vamos usar o seguinte bloco, juntamente com todas as rotações por múltiplos de 90 graus:

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

Figura 2: bloco

Código para gerar rotações (e reflexões):

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

(Esse cálculo um tanto opaco é explicado em uma resposta em /math//a/159159 , que mostra que ele simplesmente produz todas as rotações e reflexões possíveis de um bloco e remove os resultados duplicados.)

Suponha que devemos colocar o bloco como mostrado aqui:

Figura 3: posicionamento do bloco

As células 3, 6 e 7 são cobertas neste posicionamento. Isso é designado pelo vetor indicador

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

Se deslocarmos esse bloco uma coluna para a direita, esse vetor indicador seria

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

A combinação de tentar colocar peças em ambas as posições simultaneamente é determinada pela soma desses indicadores,

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

O 2 na sétima posição mostra essas sobreposições em uma célula (segunda linha para baixo, terceira coluna da esquerda). Como não queremos sobreposição, exigiremos que a soma dos vetores em qualquer solução válida não tenha entradas que excedam 1.

Acontece que, para esse problema, são possíveis 29 combinações de orientação e posição para os ladrilhos. (Isso foi encontrado com um simples código envolvendo uma pesquisa exaustiva.) Podemos descrever todas as 29 possibilidades desenhando seus indicadores como vetores de coluna . (Usar colunas em vez de linhas é convencional.) Aqui está uma imagem da matriz resultante, que terá 16 linhas (uma para cada célula possível no retângulo) e 29 colunas:

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

Figura 4: matriz de opções

(Os dois vetores indicadores anteriores aparecem como as duas primeiras colunas à esquerda.) O leitor de olhos aguçados pode ter percebido várias oportunidades de processamento paralelo: esses cálculos podem levar alguns segundos.

Todos os itens anteriores podem ser corrigidos de forma compacta usando a notação de matriz:

  • F é esse conjunto de opções, com M linhas e N colunas.

  • X é o indicador de um conjunto de canais de azulejos, de comprimento N .

  • b é um vetor N de uns.

  • R é o indicador para a região; é um vetor M.

A "bondade" total associada a qualquer solução possível X é igual à RFX , porque FX é o indicador das células cobertas por X e o produto com R soma esses valores. (Poderíamos ponderar R se desejássemos que as soluções favorecessem ou evitassem certas áreas da região.) Isso deve ser maximizado. Porque podemos escrever como ( RF ). X , é uma função linear de X : isso é importante. (No código abaixo, a variável ccontém RF .)

As restrições são que

  1. Todos os elementos de X devem ser não negativos;

  2. Todos os elementos de X devem ser menores que 1 (que é a entrada correspondente em b );

  3. Todos os elementos de X devem ser integrais.

As restrições (1) e (2) fazem deste um programa linear , enquanto o terceiro requisito o transforma em um programa linear inteiro .

Existem muitos pacotes para resolver programas lineares inteiros expressos exatamente nesta forma. Eles são capazes de manipular valores de M e N nas dezenas ou mesmo centenas de milhares. Provavelmente isso é bom o suficiente para alguns aplicativos do mundo real.


Como nossa primeira ilustração, calculei uma solução para o exemplo anterior usando o comando do Mathematica 8 LinearProgramming. (Isso minimizará uma função objetivo linear. A minimização é facilmente transformada em maximização negando a função objetivo.) Ele retornou uma solução (como uma lista de blocos e suas posições) em 0,011 segundos:

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

Figura 5: solução

As células cinzentas não estão na região; os glóbulos brancos não foram cobertos por esta solução.

Você pode elaborar (manualmente) muitas outras inclinações que são tão boas quanto essa - mas não consegue encontrar outras melhores. Essa é uma limitação potencial dessa abordagem: ela oferece uma melhor solução, mesmo quando há mais de uma. (Existem algumas soluções alternativas: se reordenarmos as colunas do X , o problema permanecerá inalterado, mas o software geralmente escolherá uma solução diferente como resultado. No entanto, esse comportamento é imprevisível.)

Como uma segunda ilustração , para ser mais realista, vamos considerar a região em questão. Ao importar a imagem e reamostrá-la, eu a representei com uma grade de 69 por 81:

Figura 6: Região

A região compreende 2156 células dessa grade.

Para tornar as coisas interessantes e ilustrar a generalidade da configuração da programação linear, vamos tentar cobrir o máximo dessa região possível com dois tipos de retângulos:

Figura 7: blocos

Um é 17 por 9 (153 células) e o outro é 15 por 11 (165 células). Podemos preferir usar o segundo, porque é maior, mas o primeiro é mais fino e pode caber em locais mais apertados. Vamos ver!

O programa agora envolve N = 5589 possíveis posicionamentos de bloco. É razoavelmente grande! Após 6,3 segundos de cálculo, o Mathematica apresentou esta solução de dez blocos:

Figura 8: solução

Devido a algumas das folgas ( .eg, poderíamos mudar o bloco inferior esquerdo para quatro colunas para a esquerda), obviamente existem outras soluções que diferem um pouco dessa.

whuber
fonte
1
Uma versão anterior desta solução (mas não tão boa) aparece no site do Mathematica em mathematica.stackexchange.com/a/6888 . Também vale a pena notar que uma pequena variação da formulação pode ser usada para resolver o problema de cobrir completamente a região com o mínimo de peças possível (permitindo algumas sobreposições, é claro): isso resolveria o "patching de buracos" problema.
whuber
1
No interesse do espaço, esta resposta não descreve algumas melhorias potencialmente úteis. Por exemplo, depois de encontrar todas as posições possíveis do bloco (como vetores indicadores), você pode adicioná-las todas para descobrir quais células podem realmente ser cobertas por algum bloco. O conjunto dessas células se divide em dois componentes conectados separados no segundo exemplo. Isso significa que o problema pode ser resolvido independentemente nos dois componentes, reduzindo substancialmente seu tamanho (e, portanto, o tempo de computação). Tais simplificações iniciais tendem a ser importantes para resolver problemas do mundo real.
whuber
Grande esforço e resposta. A resposta de Chris também foi útil. Obrigado a todos pela ajuda! Funciona e me fez seguir na direção certa novamente.
Thad
Uau! Eu estava interessado em um problema semelhante e este post me deu uma nova perspectiva. Obrigado. E se R for maior (por exemplo, 140x140≈20000), existem maneiras de reduzir o custo de computação? Você conhece algum artigo relacionado a esse problema? Minhas palavras-chave de pesquisa não me levam ao caminho certo (até agora).
Nimcap
@ nimcap Esta é uma classe importante de problemas, muita pesquisa continua. As palavras-chave nas quais pesquisar começariam com "programa linear inteiro misto" e se ramificariam a partir daí, com base no que você encontrar.
whuber
5

O link para Algoritmos genéticos para empacotamento de polígonos , fornecido na minha resposta a uma pergunta semelhante no algoritmo de busca, para colocar o número máximo de pontos dentro da área restrita com espaçamento mínimo? , pode ser útil. Parece que o método pode ser generalizado para trabalhar com formas arbitrárias de contêineres (e não apenas retângulos).

Kirk Kuykendall
fonte
Esse artigo tem algumas boas idéias (+1), mas todos os seus algoritmos se concentram, de maneira fundamental, em empacotar polígonos em regiões retangulares . Isso ocorre porque representa embalagens com uma estrutura de dados discreta (uma sequência de polígonos junto com suas orientações), que representa um conjunto de procedimentos nos quais os polígonos são deslizados , paralelos aos lados do quadrado, em direção a um canto designado. Parece que uma codificação discreta tão simples seria menos eficaz para regiões mais complicadas. Talvez uma simplificação inicial das regiões na grade ajudaria.
whuber
2

Para o subconjunto altamente restrito que você mencionou (lado a lado quadrado / triangular em um buraco), assumindo as otimizações explícitas acima, esse pseudocódigo deve chegar a uma resposta aproximada, basta levá-lo através das possibilidades com alta resolução, forçando o problema. Não funcionará corretamente em situações em que a rotação individual de blocos pode obter ganhos, como blocos de retângulo ou um contêiner altamente irregular. São 1 milhão de iterações. Você pode tentar mais, se necessário.

Assuma um quadrado com lados de comprimento L

Crie um padrão quadriculado de quadrados, que seja pelo menos das dimensões da extensão do contêiner, mais pelo menos 1 litro em cada direção.

N = 0

DX = 0

DY = 0

DR = 0

Redefinir a posição do tabuleiro de damas para o centróide original

Para (R = 1: 100)

Para (Y = 1: 100)

Para (X = 1: 100)

M = Contar o número de quadrados completamente dentro do contêiner

Se (M> N)

DR = R

DY = Y

DX = X

N = M

Mover o tabuleiro de damas para leste por L / 100

Redefinir leste do tabuleiro de damas

Mover o tabuleiro de damas para o norte por L / 100

Redefinir norte do tabuleiro de damas

Gire o tabuleiro de damas em 3,6 graus CW em torno do centróide

DY = DY * L

DX = DX * L

Redefinir o tabuleiro de damas para a posição e rotação originais

Imprimir DR & "," & DX & "e" & DY & "são a matriz final de conversão / rotação"

Girar tabuleiro de damas por DR

Traduzir tabuleiro de damas por DX, DY

Selecionar quadrados completamente dentro do contêiner

Quadrados de exportação

MappingAmanhã
fonte
Se você tentar este procedimento em uma região 2 por 5 com uma célula ausente no meio de uma borda longa, descobrirá que pode colocar apenas um quadrado 2 por 2 nela. No entanto, dois desses quadrados se encaixam facilmente. O problema é que eles não fazem parte de um padrão "quadriculado" regular. Essa dificuldade é uma das coisas que dificulta esse problema.
whuber
1
Sim. Se você tem uma forma de contêiner irregular o suficiente para suportar vários padrões regulares descontínuos da ordem de algumas células cada, isso acaba muito longe do ideal. Adicionar coisas assim à possibilidade de espaço aumenta o tempo de processamento muito rapidamente e requer um certo grau de planejamento para o caso específico que você está alvejando.
MappingTomorrow