Embalaje de polígonos dentro de un polígono con ArcGIS Desktop?

Tengo una trama booleana.

En las áreas grises del ráster, me gustaría ajustar un polígono de tamaño determinado dentro de una extensión contigua.

Básicamente, tengo un polígono irregular, y me gustaría "ajustar" un polígono conocido dentro de la extensión del polígono irregular tantas veces como sea posible.

La dirección del polígono no importa, y podría ser un cuadrado. Me gustaría que se ajustara gráficamente, pero si solo adjuntara un número al polígono (# que encaja) eso también funcionaría.

Estoy usando ArcGIS Desktop 10.

Hay muchas formas de abordar este problema. El formato ráster de los datos sugiere un enfoque basado en ráster; Al revisar esos enfoques, una formulación del problema como un programa lineal de enteros binarios parece prometedora, porque está muy en el espíritu de muchos análisis de selección de sitios GIS y puede adaptarse fácilmente a ellos.

En esta formulación, enumeramos todas las posiciones y orientaciones posibles de los polígonos de relleno, a los que me referiré como "mosaicos". Asociado con cada mosaico hay una medida de su "bondad". El objetivo es encontrar una colección de mosaicos no superpuestos cuya bondad total sea lo más grande posible. Aquí, podemos tomar la bondad de cada mosaico como el área que cubre. (En entornos de decisión más sofisticados y ricos en datos, podemos calcular la bondad como una combinación de propiedades de las celdas incluidas dentro de cada mosaico, propiedades quizás relacionadas con la visibilidad, proximidad a otras cosas, etc.)

Las limitaciones de este problema son simplemente que no pueden superponerse dos mosaicos dentro de una solución.

Esto se puede enmarcar un poco más abstracta, de una manera favorable para computación eficiente, mediante la enumeración de las células en el polígono para ser llenado (la "región") 1, 2, ..., M . Cualquier ubicación de mosaico se puede codificar con un vector indicador de ceros y unos, dejando que los correspondientes a las celdas cubiertas por el mosaico y los ceros en otros lugares. En esta codificación, toda la información necesaria sobre una colección de mosaicos se puede encontrar sumando sus vectores indicadores (componente por componente, como de costumbre): la suma será distinta de cero exactamente donde al menos un mosaico cubre una celda y la suma será mayor de uno en cualquier lugar, dos o más fichas se superponen (La suma cuenta efectivamente la cantidad de superposición).

Una más pequeña abstracción: el conjunto de posibles colocaciones de baldosas en sí puede ser enumerado, por ejemplo 1, 2, ..., N . La selección de cualquier conjunto de ubicaciones de mosaicos en sí corresponde a un vector indicador donde los que designan los mosaicos que se colocarán.

Aquí hay una pequeña ilustración para arreglar las ideas . Se acompaña con el código de Mathematica utilizado para hacer los cálculos, de modo que la dificultad de programación (o falta de ella) puede ser evidente.

Primero, representamos una región para ser en mosaico:

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

Si numeramos sus celdas de izquierda a derecha, comenzando en la parte superior, el vector indicador para la región tiene 16 entradas:

Flatten[region]

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

Usemos el siguiente mosaico, junto con todas las rotaciones por múltiplos de 90 grados:

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

Código para generar rotaciones (y reflexiones):

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

(Este cálculo algo opaco se explica en una respuesta en /math//a/159159 , que muestra que simplemente produce todas las rotaciones y reflexiones posibles de un mosaico y luego elimina cualquier resultado duplicado).

Supongamos que colocamos el mosaico como se muestra aquí:

Las celdas 3, 6 y 7 están cubiertas en esta ubicación. Eso es designado por el vector indicador

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

Si desplazamos este mosaico una columna a la derecha, ese vector indicador sería

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

La combinación de tratar de colocar fichas en ambas posiciones simultáneamente está determinada por la suma de estos indicadores,

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

El 2 en la séptima posición muestra estas superposiciones en una celda (segunda fila hacia abajo, tercera columna desde la izquierda). Debido a que no queremos superposición, requeriremos que la suma de los vectores en cualquier solución válida no debe tener entradas superiores a 1.

Resulta que para este problema, 29 combinaciones de orientación y posición son posibles para los mosaicos. (Esto se encontró con una simple codificación que implica una búsqueda exhaustiva). Podemos representar las 29 posibilidades dibujando sus indicadores como vectores de columna . (El uso de columnas en lugar de filas es convencional). Aquí hay una imagen de la matriz resultante, que tendrá 16 filas (una para cada celda posible en el rectángulo) y 29 columnas:

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

(Los dos vectores indicadores anteriores aparecen como las dos primeras columnas a la izquierda). El lector de ojos agudos puede haber notado varias oportunidades para el procesamiento paralelo: estos cálculos pueden tomar unos segundos.

Todo lo anterior puede reiniciarse de manera compacta utilizando la notación matricial:

• F es este conjunto de opciones, con M filas y N columnas.

• X es el indicador de un conjunto de ubicaciones de baldosas, de longitud N .

• b es un N- vector de unos.

• R es el indicador de la región; Es un vector M.

La "bondad" total asociada con cualquier solución posible X , es igual a RFX , porque FX es el indicador de las celdas cubiertas por X y el producto con R suma estos valores. (Podríamos ponderar R si deseamos que las soluciones favorezcan o eviten ciertas áreas de la región). Esto debe ser maximizado. Porque podemos escribirlo como ( RF ). X , es una función lineal de X : esto es importante. (En el siguiente código, la variable ccontiene RF ).

Las restricciones son que

1. Todos los elementos de X deben ser no negativos;

2. Todos los elementos de X deben ser menores que 1 (que es la entrada correspondiente en b );

3. Todos los elementos de X deben ser integrales.

Las restricciones (1) y (2) hacen de este un programa lineal , mientras que el tercer requisito lo convierte en un programa lineal entero .

Existen muchos paquetes para resolver programas lineales enteros expresados ​​exactamente de esta forma. Son capaces de manejar valores de M y N en decenas o incluso cientos de miles. Probablemente sea lo suficientemente bueno para algunas aplicaciones del mundo real.

Como nuestra primera ilustración, calculé una solución para el ejemplo anterior usando el comando de Mathematica 8 LinearProgramming. (Esto minimizará una función objetivo lineal. La minimización se convierte fácilmente en maximización al negar la función objetivo). Devolvió una solución (como una lista de mosaicos y sus posiciones) en 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 &]]];

Las celdas grises no están en la región en absoluto; los glóbulos blancos no estaban cubiertos por esta solución.

Puede resolver (a mano) muchas otras inclinaciones que son tan buenas como esta, pero no puede encontrar ninguna mejor. Esa es una limitación potencial de este enfoque: le brinda una mejor solución, incluso cuando hay más de una. (Hay algunas soluciones: si reordenamos las columnas de X , el problema no cambia, pero el software a menudo elige una solución diferente como resultado. Sin embargo, este comportamiento es impredecible).

Como segunda ilustración , para ser más realistas, consideremos la región en la pregunta. Al importar la imagen y volver a muestrearla, la representé con una cuadrícula de 69 por 81:

La región comprende 2156 celdas de esta cuadrícula.

Para hacer las cosas interesantes, y para ilustrar la generalidad de la configuración de la programación lineal, tratemos de abarcar la mayor parte posible de esta región con dos tipos de rectángulos:

Uno mide 17 por 9 (153 celdas) y el otro mide 15 por 11 (165 celdas). Podríamos preferir usar el segundo, porque es más grande, pero el primero es más delgado y puede caber en lugares más estrechos. ¡Veamos!

El programa ahora involucra N = 5589 posibles ubicaciones de mosaicos. ¡Es bastante grande! Después de 6.3 segundos de cálculo, a Mathematica se le ocurrió esta solución de diez mosaicos:

Debido a algo de la holgura ( .eg, podríamos mover el mosaico inferior izquierdo hasta cuatro columnas a su izquierda), obviamente hay algunas otras soluciones que difieren ligeramente de esta.

El enlace a Algoritmos genéticos para el embalaje de polígonos , que se proporciona en mi respuesta a una pregunta similar en el algoritmo de búsqueda para colocar el número máximo de puntos dentro del área restringida con un espaciado mínimo. , podría ser útil Parece que el método podría generalizarse para trabajar con formas de contenedor arbitrarias (y no solo rectángulos).

Para el subconjunto altamente restringido que mencionó (mosaico cuadrado / triangular en un bache), suponiendo las optimizaciones explícitas anteriores, este pseudocódigo debería llegar a una respuesta aproximada simplemente llevándolo a través de las posibilidades con una resolución bruta de alta resolución. No funcionará correctamente para situaciones en las que la rotación de mosaicos individuales puede ver ganancias, como mosaicos rectangulares o un contenedor altamente irregular. Esto es 1 millón de iteraciones, puede intentar más si es necesario.

Suponga un cuadrado con lados de longitud L

Cree un patrón de cuadros de tablero de ajedrez, que sea al menos de las dimensiones de la extensión del contenedor, más al menos 1L en cada dirección.

N = 0

DX = 0

DY = 0

DR = 0

Restablecer la posición del tablero de ajedrez al centroide original

Para (R = 1: 100)

Para (Y = 1: 100)

Para (X = 1: 100)

M = Contar el número de cuadrados completamente dentro del contenedor

Si (M> N)

DR = R

DY = Y

DX = X

N = M

Mueva el tablero de ajedrez hacia el este por L / 100

Restablecer tablero de ajedrez este

Mueva el tablero de ajedrez hacia el norte por L / 100

Restablecer tablero de ajedrez norte

Gire el tablero de ajedrez en 3.6 grados CW alrededor de su centroide

DY = DY * L

DX = DX * L

Restablecer el tablero de ajedrez a su posición original y rotación

Print DR & "," & DX & ", y" & DY & "son la matriz final de traducción / rotación"

Girar tablero de ajedrez por DR

Traducir tablero de damas por DX, DY

Seleccione cuadrados que estén completamente dentro del contenedor

Cuadrados de exportación