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

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.

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

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

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:

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

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

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:

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:

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:

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.

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

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