Sur la création des enveloppes concaves (concave hull) et les divers moyens d'y parvenir (formes alpha, alpha shapes)

	Niveau	Intermédiaire
	Logiciels utilisés	Python R QGIS GRASS GIS OpenJUMP PostGIS SpatiaLite JavaScript
	Plateforme	Windows | Mac | Linux | FreeBSD

Comment créer (automatiquement et pas à la main, je précise)  une enveloppe concave autour d'un nuage de points (concave hull) est une question qui revient souvent sur les Forums. Il existe en effet de nombreux algorithmes qui permettent de créer une enveloppe convexe (convex hull) de surface convexe minimale contenant tous les points (ou le polyèdre convexe dans le cas de points 3D) mais pas grand-chose de précis sur la création d'enveloppes non convexes.

À l'inverse de l'enveloppe convexe dont la solution est toujours unique, il y a plusieurs solutions possibles pour les enveloppes concaves, suivant les critères et l'algorithme choisi.

 

Diverses enveloppes concaves (concave hull) selon le fait que l'on autorise les partitionnements (clusters) ou les  trous

Je vais présenter ici les principes de création de ces enveloppes concaves (2D) et un inventaire des moyens de les créer avec Python, R, QGIS, GRASS GIS, OpenJUMP, PostGIS, SpatiaLite et enfin JavaScript. Mais, quelques éléments théoriques sont cependant nécessaires pour éviter quelques déconvenues... (je n'ai rien trouvé pour gvSIG, mais comme les versions pour Mac OS X sont aux abonnés absents, la dernière version disponible étant la 1.12, peut-être que...)

Quelques définitions et principes

Qu'est-ce qu'une enveloppe convexe d'un ensemble de points ? C'est la surface qui minimise la zone qui contient tous les points, sans que l'angle entre deux arêtes voisines ne dépasse 180°.

Une enveloppe concave d'un ensemble de points obéit aux mêmes principes, mais autorise tous les angles. Tous les ensembles d'enveloppes situés entre l'enveloppe convexe et celle qui minimise complètement la surface sont considérés comme concaves.

    

Diverses enveloppes concaves possibles

Se pose alors le problème du choix de la « bonne » enveloppe concave (celle qu'on recherche...). Veut-on trouver la forme originale du nuage de points (s'il y en a une) ou la surface minimale qui les englobe tous ?

Ce problème n'est pas limité au domaine géomatique, mais intéresse de nombreux autres domaines scientifiques qui doivent traiter des nuages de points/observations en 2D ou 3D (morphométrie/reconnaissance des formes, imagerie médicale, bio-informatique, cosmologie, séismologie, écologie, physique, analyse de réseaux, etc.) ou de la conception assistée par ordinateur (CAD) et bien évidemment dans les traitements divers des données LIDAR. Ce sujet a fait l'objet de nombreux articles scientifiques, de livres ou de thèses de doctorat et ces références bibliographiques sont trop nombreuses pour être citées ici (il suffit de faire une recherche sur Internet). En fin d'article, je fournis quelques-unes d'entre elles (accessibles) pour illustrer la diversité des recherches.

De ce fait, plusieurs approches mathématiques plus ou moins complexes ont été proposées pour tenter de résoudre le problème. Elles sont dans leur majorité basées sur la triangulation de Delaunay et le modèle de forme alpha  (Alpha Shape) qui sera explicité dans la suite:

• exploitant  le modèle des formes alpha s.s.  utilisé avec Python, R et par QGIS, GRASS GIS, PostGIS et SpatiaLite, voir plus bas;
• exploitant les généralisations de ces formes alpha comme  Concave hull (pdf), qui permet les partitionnements (clusters, utililisé avec un plugin de QGIS, voir plus bas),  Alpha-Concave Hull (pdf),  alpha convex hull,  etc.;
• exploitant le modèle des characteristic shapes (pdf) qui permet de générer de manières ordonnée à l'aide d'un seul paramètre normalisé toute la famille des formes alpha depuis l'enveloppe convexe jusqu'à la surface minimale (utlilisé avec une extension d'OpenJUMP, voir plus bas)
• utilisant le modèle des formes alpha pondérées, la théorie des Graphes, celui d'enveloppe de superficie minimale (MAE : computation of Minimal Area Enclosure)  celui de fractionnement et fusion ou segmentation (split and merge procedure), celui de  Locally-Density-Adaptive-α-hull ou celui de diagramme de Voronoï contraint, non abordés ici;
• et d'autres ....
  

Au premier abord, cela n'a donc pas l'air aussi simple qu'espéré...

Les formes alpha (alpha shapes)

Il est nécessaire de présenter ici le modèle des formes alpha (alpha shape ou alpha hull), car c'est LA méthode utilisée dans toutes les solutions présentées. Sans savoir ce que c'est, impossible d'aller plus loin...

Pour commencer, disons que ce modèle est intimement lié à la triangulation de Delaunay et aux enveloppes convexes qui n'en sont qu'un cas particulier.

• C'est une triangulation de Delaunay (2D ou 3D) pour laquelle chacun des triangles ou tétraèdres (« simplex ») est comparé à la valeur d'un paramètre alpha (α)  déterminé par la distance euclidienne entre les points (la méthode a été introduite par H. Edelsbrunner, D. G. Kirckpatrick and R. Seidel. en 1983 dans On the shape of a set of points in the plane (pdf) (IEEE Trans. Inform. Theory IT-29, p. 551-559). L'objectif était d'essayer de trouver les points « frontières » d'une enveloppe en permettant même d'y créer des trous;
• Pour comprendre la signification de ce paramètre α, une allégorie intuitive  avec un bloc de crème glacée a été proposée par K. Fischer en 2000 dans Introduction to Alpha Shapes (pdf) et reprise dans CGAL 4.5 - 2D Alpha Shapes (cette librairie C++ est la solution la plus performante pour les traiter, mais elle est relativement difficile à installer).
• L'ensemble est bien expliqué en français dans le mémoire de V. Carette (Université Laval, 2008) Amélioration de la représentation géométrique 2D et 3D des agrégations de poissons en support à l'étude de leur évolution spatio-temporelle (pdf) dont je me suis inspiré pour cette partie.

K. Fischer suggère d’assimiler le concept des formes alpha à un espace 3D rempli d’une importante masse de crème glacée contenant des pépites de chocolat. Avec une cuillère à glace nous creusons et retirons toute la crème glacée que nous pouvons atteindre sans toucher aux pépites (et ce y compris depuis l'intérieur). Plus le rayon de la cuillère diminue, plus cette cuillère peut atteindre des endroits étroits et plus de crème glacée sera enlevée. Avec un rayon infiniment petit, il ne restera que les pépites. À l'inverse, une cuillère de rayon trop grand nous empêchera de la bouger, laissant la masse intacte.

Dans cette histoire de crème glacée, la valeur α est le rayon de la cuillère à glace.

• si α → 0, nous pouvons « théoriquement » manger toute la crème glacée en laissant seulement les  pépites en chocolat. L’objet résultant sera délimité par des arcs, des surfaces bombées et des points sans aucune garantie que ses différentes parties soient connectées entre elles;
• si α → ∞, le résultat sera la masse de crème glacée, c'est-à-dire son enveloppe convexe.

Les formes alpha se situent entre ces 2 valeurs extrêmes. Ce sont donc des généralisations du concept de l'enveloppe convexe. Toute enveloppe convexe est une forme alpha, mais l'inverse n'est pas vrai. En utilisant une valeur appropriée d'α entre ces limites, il est dès lors possible de créer une forme alpha non convexe qui n'est pas nécessairement en un seul morceau et peut contenir un ou plusieurs trous. Le problème est de trouver cette valeur, comme nous le verrons dans la suite.

Pour conceptualiser cette histoire en 2D, il est possible d'imaginer un disque avec un certain rayon (= α) qui parcourt l'ensemble obtenu par la triangulation de Delaunay. Partout où le disque peut passer sans toucher les noeuds d'un triangle, ce triangle sera supprimé. Formellement on appelle complexe alpha (alpha complex) la partie de la triangulation de Delaunay qui subsistera et forme alpha (alpha shape) la limite du complexe alpha.

Représentation traditionnelle de la démarche

Tout ça vous parait un peu obscur, non ? Comme rien ne vaut la pratique pour éclaircir les choses, examinons comment les calculer avec la triangulation de Delaunay (2D ici). Le principe est donc d'éliminer tous les noeuds de la triangulation de Delaunay qui appartiennent à un triangle dont le rayon de la sphère circonscrite dépasse la valeur α choisie.

1) construction de la triangulation de Delaunay de l'ensemble des points;

2) une analyse statistique sur les longueurs des arêtes ou une autre technique est effectuée pour déterminer la valeur α,
3) le complexe alpha est ensuite calculé en ôtant tous les triangles qui ont un cercle circonscrit dont le rayon est supérieur au cercle de rayon α (comme le triangle 1 sur la figure, alors que le triangle 2 sera conservé);

complexe alpha résultant

4) et enfin tous les triangles restants sont dissous (dissolve) pour former la forme alpha  (comparaison avec l'enveloppe convexe qui n'est qu'un cas particulier ci-dessous).

   

Solutions

Signalons en tout premier lieu un petit programme en Java, le Non-convex hulls software qui permet de « jouer» avec les différents paramètres à partir de lettres.

Ensuite, je vais commencer par les implémentations en Python car elles suivent strictement la démarche décrite précédemment et que ce sont elles qui m'ont vraiment aidé à comprendre le problème. C'est plus facile ensuite d'aborder les autres solutions en toute connaissance de cause.

Python seul:

Les algorithmes de création de formes alpha ont été bien détaillés dans:

• The Fading Shape of Alpha avec Shapely et Scipy Delaunay;
• Drawing Boundaries In Python et Drawing Boundaries avec Fiona en plus.

Toutes les démarches sont réalisées avec d'abord la triangulation de Delaunay puis le choix du paramètre alpha pour arriver à la solution. Le principal problème est qu'il faut trouver la « bonne » valeur α soi-même. En pratique cependant, si l'on désire comme résultat un polygone unique (ce n'est pas toujours possible, surtout si la forme est trop convexe), il suffit d'une simple boucle en testant le résultat obtenu par les fonctions cascaded_union ou unary_union de Shapely qui devra être un Polygone et non un MultiPolygone comme illustré dans mon script  concave_hull1.py.

Obtention de MultiPolygones résultant d'une valeur α trop petite

Polygone unique au bout de n itérations

Pour éviter cette étape de recherche, une autre solution est d'utiliser un programme spécialisé avec le module subprocess. Parmi ceux-ci, le programme C Hull (il faut le compiler) de Ken Clarkson est un des plus performants et a été utilisé dans:

• Alpha Shapes in Python et Concave Hull 

Ici, plus besoin de se préoccuper de la recherche de la valeur α, tout est fait automatiquement par le programme. Un exemple d'application est donné dans mon script concave_hull2.py ou un rayon de recherche peut être spécifié.

 

CGAL possède aussi une implémentation en Python, mais comme je l'ai dit, la librairie est difficile à installer. Si vous voulez tout de même la tester, il est possible de trouver quelques exemples sur Internet comme dans AlphaShape.py

R

R permet de créer des enveloppes concaves en 2D et en 3D avec les packages alphahull et alphashape3D (forme alpha).

• alphahull: an R Package for Alpha-Convex Hull
• Generalizing the Convex Hull of a Sample: The R Package alphahull (pdf)
• α – Shape generation with R
• Alpha shapes with R and lattice

Résultat avec alphahull, voir mon script alpha_shape.R (tout comme les solutions en Python, le paramètre α doit être recherché indépendamment).

Quelques exemples non géomatiques avec alphashape3D:

 

Il existe d'autres packages plus limités comme  rLIDAR.

QGIS

QGIS offre plusieurs manières de créer des enveloppes concaves:

• avec l'algorithme Concave Hull de la Boîte à outils de Traitements (Processing) QGIS processing: Concave hull (ConcaveHull.py). Il s'agit ici d'un algorithme de forme alpha qui permet la création de polygones troués (si on le veut). Le Threshold peut être assimilé au paramètre α.

 

• avec le plugin encore expérimental de Detlev Neumann (il faut le télécharger et l'installer soi-même) en Python, QGIS-ConcaveHull-Plugin. Il n'utilise pas la triangulation de Delaunay proprement dite, mais l'algorithme des K plus proches voisins suivant la méthode proposée par A. Moreira and M.B. Santos dans Concave hull : a k-nearest neighbours approach for the computation of the region occupied by a set of points (pdf) en 2007. Il permet en plus de séparer les résultats en groupes (clusters);

 

• avec le plugin PgRouting Layer d'Anita Graser qui utilise les fonctions de  PostGIS (voir plus bas).
• il est toujours possible d'utiliser les solutions Python ou R avec un script dans la Boîte à outils de Traitements (processing) ou directement les fonctions de PostGIS.

GRASS GIS

• Les différentes méthodes pour créer des enveloppes concaves sont présentées dans GRASS GIS Create concave hull (forme alpha)

OpenJUMP

• Comme toujours, le méconnu OpenJUMP offre aussi une extension très performante qui permet de créer des enveloppes concaves: Concave hull based on JTS. Elle est basée sur l'algorithme  « characteristic shape » proposé par M.Duckham, L.Kulik, M.Worboys & A.Galton en 2008 dans Efficient generation of simple polygons for characterizing the shape of a set of points in the plane (pdf). La forme alpha produite par la commande est contrôlée par la variable seuil (le paramètre normalisé):

Divers résultats obtenus en modifiant simplement la valeur du seuil

PostGIS

• il y a moyen de créer des enveloppes concaves avec l'extension Driving Distance extension de pgRouting, voir A Closer Look at Alpha Shapes in pgRouting de Anita Grasser ou Concave Hull – a GIS problem put to bed at last!
• les versions 2.x offrent la commande ST_ConcaveHull, voir PostGIS 2 Concave Hull;
• Anita Graser a créé un plugin pour QGIS PgRouting Layer qui exploite cette caractéristique.

SpatiaLite

• il existe aussi une implémentation dans les versions 4.x voir Delaunay Triangulation, Convex Hull and Concave Hull

JavaScript

JavaScript permet de le faire

• avec la librairie Hull.js, voir Generate Concave Hulls with Hull.js;

 images tirées de Hull.js examples

 

Divers essais avec la couche test:

• avec la nouvelle librairie Turf.js ("a modular GIS engine written in JavaScript"), voir turf-concave;
• il y a d'autres solutions comme dans alpha-shapes aka concave hulls in d3 ou ConcaveHull qui utilise la méthode Concave hull (pdf) déjà vue.

Et le reste...

Je vous invite à consulter GIS SE: What are Definition, Algorithms and Practical Solutions for Concave Hull? qui montre d'autres solutions comme:

• The Shape of Alpha de Flickr ;

• certains programmes/librairies de traitement des données LIDAR;

Conclusions

Vous aurez sans aucun doute remarqué qu'aucune enveloppe concave obtenue ici n'est strictement la même et pourtant elles sont toutes correctes.

Diverses solutions examinées

Comme vous avez le choix entre plusieurs méthodes et plusieurs paramètres, c'est à vous de jouer maintenant pour obtenir ce que vous voulez, car somme toute, c'est le résultat recherché qui compte (si on ne veut pas le faire à la main...). Les recherches scientifiques actuelles portent surtout sur les enveloppes concaves 3D.

Tous les traitements ont été effectués sur Mac OS X et les figures sont des captures d'écran des représentations avec QGIS hormis les solutions en JavaScript.

Site officiel : What Is the Region Occupied by a Set of Points? (pdf) de A.Galton and M. Duckham
Site officiel : Alpha Shapes — a Survey (pdf) de H.Edelsbrunner
Site officiel : Alpha shapes de M.Celikik (pdf)
Autres Liens : Alpha-Concave Hull, a Generalization of Convex Hull (pdf) de S. Asaeedi, F. Didehvar, and A.Mohades
Autres Liens : Implementation of fast and efficient concave hull algorithm (pdf) de E. Rosén, E.Jansson, M. Brundin
Autres Liens : Minimum area enclosure and alpha hull of a set of freeform planar closed curves (pdf) de A.V. Vishwanath, R. Arun Srivatsan & M. Ramanathan
Autres Liens : Building boundary extraction based on LIDAR point clouds data (pdf) de de SHEN Wei
Autres Liens : A split and merge procedure for polygonal border detection of dot pattern (pdf) de G. Garaia, B.B. Chaudhurib
Autres Liens : Alpha Shapes for Computing Catchment Areas
Autres Liens : Predicting species distributions in new areas or time periods with alpha-shapes (pdf) de C. Capinha et B. Pateiro-López
Autres Liens : Application de la théorie des formes α pour la caractérisation de la surface et des poches de macromolécules biologiques (pdf, thèse de doctorat de Benjamin Schwarz)

licence Creative Commons Paternité-Pas d'Utilisation Commerciale-Pas de Modification 2.0 France