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

25

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.

Thad
la source
8
C'est un problème très difficile. Par exemple, il faut beaucoup de travail pour insérer autant de cercles dans un carré que possible. Lorsque le polygone d'origine est compliqué - comme dans l'illustration - vous avez besoin de puissantes procédures d'optimisation. La meilleure méthode que j'ai trouvée pour ce problème est le recuit simulé, mais cela ne sera pas disponible dans ArcGIS et il faudrait des scripts extrêmement astucieux pour le scripter (ArcGIS est trop lent). Pourriez-vous peut-être assouplir un peu vos exigences, comme ajuster le plus petit polygone un nombre suffisant de fois, plutôt que le plus de fois possible?
whuber
1
@whuber Merci d'avoir édité mon message. Oui, un nombre suffisant de fois fonctionnerait. Ou, que diriez-vous d'une orientation d'angle donnée. ex. dans l'image ci-dessus, j'ai ajusté le polygone autant de fois que je pourrais avoir à cette orientation, si je les avais fait pivoter de 90 degrés, vous pourriez en ajuster un de plus ...
Thad
1
Oui, mais c'est aussi semé d'embûches. Certains sont élémentaires. Par exemple, le texte rédigé et publié par ESRI, «Getting to Know ArcView GIS» (pour la version 3), comprenait un exercice dans lequel un rectangle représentant un terrain de football était placé de manière interactive dans un polygone. Le problème était que la réponse à l'exercice était erronée car l'auteur n'avait pas projeté les données et les erreurs dans l'utilisation des coordonnées géographiques étaient suffisamment importantes pour affecter le résultat. La réponse avait l'air bonne dans le SIG, mais si quelqu'un avait tenté de construire ce champ, il aurait trouvé qu'il n'y avait pas assez de place pour cela :-).
whuber
6
@whuber Je suppose qu'ils pensaient qu'un chiffre "parc à billes" était suffisant.
Kirk Kuykendall
2
Dans le cas général d'un polygone irrégulier à l'intérieur d'un polygone irrégulier, il s'agit d'un problème insoluble sur le plan informatique: trouver une solution optimale n'est pas un objectif plausible dans tous les cas, et il est probablement NP-complet d'un point de vue technique: quels cas ceux-ci ne peuvent pas être prédéterminés. Si vous limitez considérablement le problème, certains algorithmes itératifs d'ajustement aléatoire sont susceptibles de vous donner des nombres raisonnablement élevés . Mon sentiment s'il s'agit d'une mission est qu'ils ne recherchent pas la bonne réponse, ils recherchent des approches créatives.
MappingTomorrow

Réponses:

22

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

Figure 1: région

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

Figure 2: tuile

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:

Figure 3: placement des tuiles

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

Figure 4: tableau d'options

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

Figure 5: solution

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:

Figure 6: Région

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:

Figure 7: tuiles

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:

Figure 8: solution

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.

whuber
la source
1
Une version antérieure de cette solution (mais pas aussi bonne) apparaît sur le site Mathematica à l' adresse mathica.stackexchange.com/a/6888 . Il convient également de noter qu'une variation mineure de la formulation peut être utilisée pour résoudre le problème de la couverture complète de la région avec le moins de tuiles possible (en permettant bien sûr des chevauchements): cela résoudrait le "colmatage des nids-de-poule" problème.
whuber
1
Dans l'intérêt de l'espace, cette réponse ne décrit pas certaines améliorations potentiellement utiles. Par exemple, après avoir trouvé toutes les positions de tuiles possibles (en tant que vecteurs indicateurs), vous pouvez toutes les additionner pour trouver quelles cellules peuvent être réellement couvertes par une tuile. L'ensemble de ces cellules se divise en deux composants connectés distincts dans le deuxième exemple. Cela signifie que le problème peut être résolu indépendamment dans les deux composants, réduisant considérablement sa taille (et donc le temps de calcul). Ces simplifications initiales ont tendance à être importantes pour résoudre les problèmes du monde réel.
whuber
Grand effort et réponse. La réponse de Chris a également été utile. Merci à tous pour l'aide! Fonctionne et m'a fait avancer dans la bonne direction.
Thad
Hou la la! J'étais intéressé par un problème similaire et ce post m'a donné une nouvelle perspective. Merci. Que faire si R est plus grand (par exemple 140x140≈20000), existe-t-il des moyens de réduire les coûts de calcul? Connaissez-vous des articles liés à ce problème? Mes mots clés de recherche ne me mènent pas dans le bon sens (jusqu'à présent).
nimcap
@nimcap Il s'agit d'une classe de problèmes importante, de nombreuses recherches sont en cours. Les mots clés à rechercher commenceraient par "programme linéaire à nombres entiers mixtes" et se développeraient à partir de ce que vous trouverez.
whuber
5

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

Kirk Kuykendall
la source
Cet article a quelques bonnes idées (+1), mais tous ses algorithmes se concentrent, de manière fondamentale, sur le compactage de polygones dans des régions rectangulaires . En effet, il représente des emballages avec une structure de données discrète (une séquence de polygones avec leurs orientations) qui représente un ensemble de procédures dans lesquelles les polygones sont glissés , parallèlement aux côtés du carré, vers un coin désigné. Il apparaît qu'un tel codage discret simple serait moins efficace pour des régions plus compliquées. Peut-être qu'une simplification initiale des régions du réseau serait utile.
whuber
2

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

MappingTomorrow
la source
Si vous essayez cette procédure sur une région de 2 x 5 avec une cellule manquante au milieu d'un bord long, vous constaterez que vous ne pouvez y mettre qu'un seul carré de 2 x 2. Cependant, deux de ces carrés s'adaptent facilement. Le problème est qu'ils ne font pas partie d'un modèle de "damier" régulier. Cette difficulté est l'une des choses qui rendent ce problème assez difficile.
whuber
1
Ouaip. Si vous avez une forme de conteneur suffisamment irrégulière pour qu'elle puisse prendre en charge plusieurs modèles réguliers non contigus de l'ordre de quelques cellules chacune, cela finit très loin d'être optimal. Ajouter de telles choses à l'espace de possibilité augmente le temps de traitement très rapidement et nécessite un certain degré de planification pour le cas particulier que vous ciblez.
MappingTomorrow