Hay muchas formas de abordar este problema. El formato raster de los datos sugiere un enfoque basado en raster; al revisar esos enfoques, una formulación del problema como una programa lineal entero binario parece prometedor, ya que sigue el espíritu de muchos análisis SIG de selección de emplazamientos y puede adaptarse fácilmente a ellos.
En esta formulación, enumeramos todas las posiciones y orientaciones posibles del polígono o polígonos de relleno, que denominaré "baldosas". Asociada a cada baldosa hay una medida de su "bondad". El objetivo es encontrar una colección de baldosas no superpuestas cuya bondad total sea la mayor posible. En este caso, podemos considerar que la bondad de cada mosaico es la superficie que cubre. (En entornos de decisión más ricos en datos y sofisticados, podemos calcular la bondad como una combinación de propiedades de las celdas incluidas dentro de cada baldosa, propiedades quizá relacionadas con la visibilidad, la proximidad a otras cosas, etc.).
Las restricciones de este problema consisten simplemente en que no pueden solaparse dos fichas dentro de una solución.
Esto puede plantearse de forma un poco más abstracta, de un modo que permita un cálculo eficiente, enumerando las celdas del polígono que se va a rellenar (la "región") 1, 2, ..., M . Cualquier colocación de baldosas puede codificarse con un indicador vector de ceros y unos, dejando que los unos correspondan a las celdas cubiertas por la baldosa y los ceros a las demás. En esta codificación, toda la información necesaria sobre una colección de baldosas se puede encontrar mediante sumando sus vectores indicadores (componente por componente, como de costumbre): la suma será distinta de cero exactamente donde al menos una baldosa cubra una celda y la suma será mayor que uno en cualquier lugar donde se solapen dos o más baldosas. (La suma cuenta efectivamente la cantidad de solapamiento).
Una pequeña abstracción más: el conjunto de posibles colocaciones de las baldosas se puede enumerar, digamos 1, 2, ..., N . La selección de cualquier conjunto de colocaciones de baldosas se corresponde a su vez con un vector indicador en el que los unos designan las baldosas que deben colocarse.
He aquí una pequeña ilustración para fijar las ideas . Se acompaña con el Mathematica código utilizado para realizar los cálculos, de modo que pueda evidenciarse la dificultad de programación (o la falta de ella).
En primer lugar, representamos una región que se va a embaldosar:
region = {{0, 0, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}};
![Figure 1: region]()
Si numeramos sus casillas de izquierda a derecha, empezando por arriba, el vector indicador de la región tiene 16 entradas:
Flatten[region]
{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
Utilicemos el siguiente azulejo, junto con todas las rotaciones por múltiplos de 90 grados:
tileSet = {{{1, 1}, {1, 0}}};
![Figure 2: tile]()
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 https://math.stackexchange.com/a/159159 que muestra que simplemente produce todas las rotaciones y reflexiones posibles de un azulejo y luego elimina cualquier resultado duplicado).
Supongamos que colocamos el azulejo como se muestra aquí:
![Figure 3: tile placement]()
Las celdas 3, 6 y 7 están cubiertas en esta colocación. Que se designa por el vector indicador
{0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}
Si desplazamos este azulejo una columna a la derecha, ese vector indicador sería en cambio
{0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0}
La combinación de intentar colocar fichas en ambos estas posiciones simultáneamente viene 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 estos solapamientos en una celda (segunda fila hacia abajo, tercera columna desde la izquierda). Porque no queremos que se solapen, exigiremos que la suma de los vectores de cualquier solución válida no tenga entradas superiores a 1.
Resulta que para este problema son posibles 29 combinaciones de orientación y posición de las fichas. (Podemos representar las 29 posibilidades dibujando sus indicadores de la siguiente manera columna vectores. (A continuación se muestra una imagen de la matriz resultante, que tendrá 16 filas (una por cada celda posible del 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];
![Figure 4: options array]()
(Los dos vectores indicadores anteriores aparecen como las dos primeras columnas de la izquierda). El lector avispado habrá notado varias oportunidades de procesamiento paralelo: estos cálculos pueden tardar unos segundos.
Todo lo anterior puede resumirse utilizando la notación matricial:
-
F es esta matriz de opciones, con M filas y N columnas.
-
X es el indicador de un conjunto de colocaciones de baldosas, de longitud N .
-
b es un N -vector de unos.
-
R es el indicador de la región; se trata de un M -vector.
La "bondad" total asociada a cualquier solución posible X , es igual a R.F.X porque F.X es el indicador de las celdas cubiertas por X y el producto con R suma estos valores. (Podríamos ponderar R si deseábamos que las soluciones favorecieran o evitaran determinadas zonas de la región). Hay que maximizarlo. Porque podemos escribirlo como ( R.F ). X es un lineal función de X Esto es importante. (En el código siguiente, la variable c contiene R.F .)
Las restricciones son que
-
Todos los elementos de X debe ser no negativo;
-
Todos los elementos de X debe ser inferior a 1 (que es la entrada correspondiente en b );
-
Todos los elementos de X debe ser integral.
Las restricciones (1) y (2) hacen que este programa lineal mientras que el tercer requisito lo convierte en un entero programa lineal.
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 suficiente para algunas aplicaciones del mundo real.
Como primera ilustración, He calculado una solución para el ejemplo anterior utilizando Mathematica 8's LinearProgramming comando. (Esto minimizar una función objetivo lineal. La minimización se convierte fácilmente en maximización negando la función objetivo). Devuelve una solución (como una lista de fichas 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 &]]];
![Figure 5: solution]()
Las células grises no están en absoluto en la región; las células blancas no estaban cubiertas por esta solución.
Se pueden calcular (a mano) muchos otros tilings tan buenos como éste, pero no se pueden encontrar mejores. Esa es una limitación potencial de este enfoque: te da un mejor solución, incluso cuando hay más de una. (Existen algunas soluciones: si reordenamos las columnas de X El problema no cambia, pero el programa suele elegir una solución diferente. Sin embargo, este comportamiento es imprevisible).
Como segunda ilustración Para ser más realistas, consideremos la región de la pregunta. Al importar la imagen y remuestrearla, la he representado con una cuadrícula de 69 por 81:
![Figure 6: Region]()
La región comprende 2156 celdas de esta cuadrícula.
Para hacer las cosas más interesantes, y para ilustrar la generalidad de la programación lineal, vamos a intentar cubrir la mayor parte posible de esta región con dos tipos de rectángulos:
![Figure 7: tiles]()
Una mide 17 por 9 (153 celdas) y la otra 15 por 11 (165 celdas). Es posible que prefiramos utilizar 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 incluye ahora N \= 5589 posibles colocaciones de baldosas. ¡Es bastante grande! Después de 6,3 segundos de cálculo, Mathematica se le ocurrió esta solución de diez baldosas:
![Figure 8: solution]()
Debido a algunas de las holguras ( .e.g, podríamos desplazar la baldosa inferior izquierda hasta cuatro columnas a su izquierda), evidentemente existen otras soluciones que difieren ligeramente de ésta.
