Emballage de polygones dans un polygone à l'aide d'ArcGIS Desktop?

J'ai un raster booléen.

Dans les zones grises du raster, je voudrais ajuster un polygone de taille donnée dans une étendue contiguë.

Fondamentalement, j'ai un polygone irrégulier et j'aimerais "ajuster" un polygone connu dans l'étendue du polygone irrégulier autant de fois que possible.

La direction du polygone n'a pas d'importance et il pourrait s'agir d'un carré. J'aimerais qu'il s'adapte graphiquement, mais s'il attachait simplement un nombre au polygone (# qui correspond), cela fonctionnerait aussi.

J'utilise ArcGIS Desktop 10.

Il existe de nombreuses façons d'aborder ce problème. Le format raster des données suggère une approche basée sur un raster; en examinant ces approches, une formulation du problème sous la forme d'un programme linéaire d'entiers binaires semble prometteuse, car elle est très conforme à l'esprit de nombreuses analyses de sélection de sites SIG et peut facilement être adaptée à celles-ci.

Dans cette formulation, nous énumérons toutes les positions et orientations possibles du ou des polygones de remplissage, que j'appellerai des «tuiles». Associé à chaque tuile est une mesure de sa «bonté». L'objectif est de trouver une collection de tuiles non superposées dont la bonté totale est la plus grande possible. Ici, nous pouvons prendre la bonté de chaque tuile pour être la zone qu'elle couvre. (Dans des environnements de décision plus riches en données et sophistiqués, nous pouvons calculer la qualité comme une combinaison de propriétés des cellules incluses dans chaque mosaïque, propriétés peut-être liées à la visibilité, à la proximité d'autres choses, etc.)

Les contraintes de ce problème sont simplement qu’il n’y a pas deux chevauchements possibles dans une solution.

Cela peut être encadré un peu plus abstraitement, d'une manière propice au calcul efficace, les cellules en dénombrant dans le polygone à remplir (la « région ») 1, 2, ..., M . Tout placement de tuile peut être encodé avec un vecteur indicateur de zéros et de uns, en laissant ceux correspondant aux cellules couvertes par les tuiles et les zéros ailleurs. Dans cet encodage, toutes les informations nécessaires sur une collection de tuiles peuvent être trouvées en additionnant leurs vecteurs indicateurs (composant par composant, comme d'habitude): la somme sera non nulle exactement là où au moins une tuile couvre une cellule et la somme sera plus grande que n'importe où deux ou plusieurs tuiles se chevauchent. (La somme compte effectivement la quantité de chevauchement.)

Un peu plus abstraction: l'ensemble des placements de tuiles possibles peut être lui - même énumérait, par exemple 1, 2, ..., N . La sélection de n'importe quel ensemble de placements de tuiles correspond elle-même à un vecteur indicateur où ceux désignent les tuiles à placer.

Voici une petite illustration pour corriger les idées . Il est accompagné du code Mathematica utilisé pour effectuer les calculs, afin que la difficulté de programmation (ou son absence) puisse être évidente.

Tout d'abord, nous décrivons une région à carreler:

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

Si nous numérotons ses cellules de gauche à droite, en commençant par le haut, le vecteur indicateur de la région a 16 entrées:

Flatten[region]

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

Utilisons la tuile suivante, avec toutes les rotations par multiples de 90 degrés:

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

Code pour générer des rotations (et réflexions):

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

(Ce calcul quelque peu opaque est expliqué dans une réponse à /math//a/159159 , qui montre qu'il produit simplement toutes les rotations et réflexions possibles d'une tuile, puis supprime tous les résultats en double.)

Supposons que nous devions placer la tuile comme indiqué ici:

Les cellules 3, 6 et 7 sont couvertes dans cet emplacement. Qui est désigné par le vecteur indicateur

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

Si nous déplaçons cette tuile d'une colonne vers la droite, ce vecteur indicateur serait plutôt

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

La combinaison d'essayer de placer des tuiles à ces deux positions simultanément est déterminée par la somme de ces indicateurs,

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

Le 2 en septième position montre ces chevauchements dans une cellule (deuxième rangée vers le bas, troisième colonne depuis la gauche). Parce que nous ne voulons pas de chevauchement, nous exigerons que la somme des vecteurs dans toute solution valide ne contienne aucune entrée supérieure à 1.

Il s'avère que pour ce problème, 29 combinaisons d'orientation et de position sont possibles pour les tuiles. (Cela a été trouvé avec un simple codage impliquant une recherche exhaustive.) Nous pouvons décrire les 29 possibilités en dessinant leurs indicateurs comme des vecteurs de colonne . (L'utilisation de colonnes au lieu de lignes est conventionnelle.) Voici une image du tableau résultant, qui aura 16 lignes (une pour chaque cellule possible dans le rectangle) et 29 colonnes:

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

(Les deux vecteurs indicateurs précédents apparaissent comme les deux premières colonnes à gauche.) Le lecteur aiguisé a peut-être remarqué plusieurs possibilités de traitement parallèle: ces calculs peuvent prendre quelques secondes.

Tout ce qui précède peut être reformulé de manière compacte en utilisant la notation matricielle:

• F est ce tableau d'options, avec M lignes et N colonnes.

• X est l'indicateur d'un ensemble d'emplacements de tuiles, de longueur N .

• b est un N- vecteur de uns.

• R est l'indicateur de la région; c'est un vecteur M.

La «bonté» totale associée à toute solution X possible est égale à RFX , car FX est l'indicateur des cellules couvertes par X et le produit avec R additionne ces valeurs. (Nous pourrions pondérer R si nous voulions que les solutions favorisent ou évitent certaines zones de la région.) Ceci doit être maximisé. Parce que nous pouvons l'écrire comme ( RF ). X , c'est une fonction linéaire de X : c'est important. (Dans le code ci-dessous, la variable ccontient RF .)

Les contraintes sont que

1. Tous les éléments de X doivent être non négatifs;

2. Tous les éléments de X doivent être inférieurs à 1 (qui est l'entrée correspondante en b );

3. Tous les éléments de X doivent être intégrés.

Les contraintes (1) et (2) en font un programme linéaire , tandis que la troisième exigence le transforme en un programme linéaire entier .

Il existe de nombreux packages pour résoudre des programmes linéaires entiers exprimés exactement sous cette forme. Ils sont capables de gérer des valeurs de M et N par dizaines, voire centaines de milliers. C'est probablement assez bon pour certaines applications du monde réel.

Comme première illustration, j'ai calculé une solution pour l'exemple précédent en utilisant la commande de Mathematica 8 LinearProgramming. (Cela minimisera une fonction objectif linéaire. La minimisation est facilement transformée en maximisation en annulant la fonction objectif.) Il a renvoyé une solution (sous forme de liste de tuiles et de leurs positions) en 0,011 seconde:

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

Les cellules grises ne sont pas du tout dans la région; les globules blancs n'étaient pas couverts par cette solution.

Vous pouvez travailler (à la main) de nombreux autres carrelages aussi bons que celui-ci - mais vous ne pouvez pas en trouver de meilleurs. C'est une limitation potentielle de cette approche: elle vous offre une meilleure solution, même lorsqu'il y en a plusieurs. (Il existe certaines solutions: si nous réorganisons les colonnes de X , le problème reste inchangé, mais le logiciel choisit souvent une solution différente en conséquence. Cependant, ce comportement est imprévisible.)

Comme deuxième illustration , pour être plus réaliste, considérons la région dans la question. En important l'image et en la rééchantillonnant, je l'ai représentée avec une grille de 69 par 81:

La région comprend 2156 cellules de cette grille.

Pour rendre les choses intéressantes et pour illustrer la généralité de la configuration de la programmation linéaire, essayons de couvrir autant de cette région que possible avec deux types de rectangles:

L'un est de 17 sur 9 (153 cellules) et l'autre de 15 sur 11 (165 cellules). Nous préférerons peut-être utiliser le second, car il est plus grand, mais le premier est plus maigre et peut tenir dans des endroits plus étroits. Voyons voir!

Le programme comprend désormais N = 5589 emplacements de tuiles possibles. C'est assez gros! Après 6,3 secondes de calcul, Mathematica a trouvé cette solution à dix carreaux:

En raison d'une partie du jeu ( .eg, nous pourrions déplacer la tuile inférieure gauche jusqu'à quatre colonnes vers sa gauche), il existe évidemment d'autres solutions légèrement différentes de celle-ci.

Le lien vers On Genetic Algorithms for the Packing of Polygons , fourni dans ma réponse à une question similaire à Seeking algorithm pour placer le nombre maximum de points dans la zone contrainte à un espacement minimum? , pourrait être utile. Il semble que la méthode pourrait être généralisée pour fonctionner avec des formes de conteneurs arbitraires (et pas seulement des rectangles).

Pour le sous-ensemble hautement contraint que vous avez mentionné (pavage carré / triangulaire dans un nid de poule), en supposant les optimisations explicites ci-dessus, ce pseudocode devrait arriver à une réponse approximative en vous expliquant simplement les possibilités avec une résolution élevée, une brute forçant le problème. Cela ne fonctionnera pas correctement dans les situations où la rotation individuelle des carreaux peut voir des gains, comme les carreaux rectangulaires ou un conteneur très irrégulier. C'est 1 million d'itérations, vous pouvez en essayer plus si nécessaire.

Supposons un carré avec des côtés de longueur L

Créez un motif en damier de carrés, qui est au moins des dimensions de l'étendue du conteneur, plus au moins 1L dans chaque direction.

N = 0

DX = 0

DY = 0

DR = 0

Réinitialiser la position du damier au centroïde d'origine

Pour (R = 1: 100)

Pour (Y = 1: 100)

Pour (X = 1: 100)

M = Compter le nombre de carrés complètement dans le conteneur

Si (M> N)

DR = R

DY = Y

DX = X

N = M

Déplacer le damier vers l'est de L / 100

Réinitialiser l'abscisse du damier

Déplacer le damier vers le nord de L / 100

Réinitialiser l'ordonnée du damier

Faire pivoter le damier de 3,6 degrés CW autour de son centre de gravité

DY = DY * L

DX = DX * L

Réinitialiser le damier à sa position et rotation d'origine

Imprimer DR & "," & DX & "et" & DY & "sont la matrice de traduction / rotation finale"

Tourner le damier par DR

Traduire le damier par DX, DY

Sélectionnez les carrés qui sont complètement dans le conteneur

Exporter des carrés