Parcours de Graham

Ça fait déjà deux ans que je n'ai pas écrit sur ce blogue. Il est décidément temps que je m'y remette plus activement.

Depuis quelques mois, j'ai recommencé à utiliser le langage de programmation fonctionnel Haskell. Comme cela fait longtemps et que je n'en ai jamais été expert, je partage maintenant avec vous les étapes que j'ai suivies afin d'implémenter le parcours de Graham qui permet de calculer l'envelopper convexe d'un nuage de points dans le plan.

Dans un premier temps, comme nous manipulons des points 2D, il est pratique de définir le type de données correspondant. Nous supposerons que les coordonnées sont de type Double.

-- 2D points
data Point2D = Point2D Double Double
               deriving (Show, Eq)

Rappelons que, grâce à l'expression deriving (Show, Eq), il est possible d'afficher (Show) un point lorsqu'on utilise l'interpréteur ghci et de tester l'égalité (Eq) de deux points en comparant leurs coordonnées.

Nous aurons également besoin de calculer le carré de la distance entre deux points :

-- Square of distance between points
squareDistance :: Point2D -> Point2D -> Double
squareDistance (Point2D a b) (Point2D c d)
    = (c - a) * (c - a) + (d - b) * (d - b)

On peut alors construire un nuage de points que nous utiliserons pour les tests dans la suite :

-- Points
points = [Point2D 0 0, Point2D 2 9, Point2D 3 4,
          Point2D 0 5, Point2D 6 2, Point2D 7 3,
          Point2D 4 1, Point2D 5 7]

On peut obtenir une représentation graphique du nuage à l'aide de Sage Cell Server :

Rappelons les étapes du parcours de Graham :

1. On identifie le point le plus bas. En cas d'égalité, on choisis le point le plus à gauche. Nous appelons ce point l'origine.
2. On trie les points en ordre croissant d'angle par rapport à l'origine. S'il y a plusieurs points avec le même angle, on choisit alors celui qui est le plus près de l'origine. Remarquons que, par choix de l'origine, l'angle est nécessairement entre 0^∘ et 180^∘.
3. Pour chaque triplet successif de points, on vérifie s'ils forment un virage à gauche, à droite ou s'ils sont colinéaires.
   • Si les points forment un virage à gauche, on passe au prochain triplet;
   • Si les points forment un virage à droite ou s'ils sont colinéaires, alors on peut éliminer le point central de l'enveloppe convexe. On conserve les deux autres points dans le prochain triplet.

L'étape 1 est facile. Il suffit de se définir un ordre lexicographique sur les points d'abord par rapport à la coordonnée y puis ensuite la coordonnée x :

-- Lexicographic order on points (with respect to y and then to x)
comparePoints :: Point2D -> Point2D -> Ordering
comparePoints (Point2D a b) (Point2D c d)
    | b < d  = LT
    | b > d  = GT
    | b == d = compare a c

Cette fonction retourne la valeur LT, EQ ou GT (de type Ordering) selon le cas. L'origine est alors récupérée à l'aide de l'expression

List.minimumBy comparePoints points

Pour réaliser l'étape 2, on doit également introduire une relation d'ordre selon l'angle en coordonnées polaires de chaque point par rapport à l'origine calculée plus haut. S'il y a égalité d'angle, on choisit le point le plus proche. Cette condition nous assure entre autres que l'origine se trouve toujours comme premier point dans la liste.

-- Angle comparator between two points with respect to some origin
compareAngle :: Point2D -> Point2D -> Point2D -> Ordering
compareAngle origin@(Point2D a b) p@(Point2D c d) q@(Point2D e f)
        | cross > 0  = LT
        | cross == 0 = Ord.comparing (squareDistance origin) p q
        | cross < 0  = GT
    where cross = (c - a) * (f - b) - (d - b) * (e - a)

Le premier point est l'origine. Les deux points suivants sont ceux qu'on désire comparer. Le test utilisé est très standard en géométrie algorithmique : il s'agit de vérifier le signe du produit vectoriel entre deux vecteurs induits par les trois points. Selon que le signe est positif, négatif ou nul, on en déduit que l'angle du premier point est plus petit, plus grand ou égal au second. Le lecteur moins familier avec l'algèbre linéaire peut vérifier cette propriété avec l'interact ci-bas :

Ainsi, dans la fonction principale, on effectue le tri initial des points :

-- Convex hull computation
convexHull :: [Point2D] -> [Point2D]
convexHull points = grahamScan (List.sortBy (compareAngle origin) points)
    where origin = List.minimumBy comparePoints points

On constate que, une fois le prétraitement complété, la fonction convexHull fait appel à une fonction grahamScan qui correspond à l'étape 3 de l'algorithme :

-- Graham scan
grahamScan []         = []
grahamScan [p]        = [p]
grahamScan [p,q]      = [p,q]
grahamScan (p:q:r:ps)
    | compareAngle q p r == GT = p:(grahamScan (q:r:ps))
    | otherwise                = grahamScan (p:r:ps)

En appelant la fonction sur le nuage de points initialement défini, on obtient bien le résultat attendu :

Main> convexHull points
[Point2D 0.0 0.0,Point2D 4.0 1.0,Point2D 6.0 2.0,Point2D 7.0 3.0,
 Point2D 5.0 7.0,Point2D 2.0 9.0,Point2D 0.0 5.0]

Voici le code complet :

-- Imports
import qualified Data.List as List
import qualified Data.Ord as Ord

-- 2D points
data Point2D = Point2D Double Double
               deriving (Show, Eq)

-- Square of distance between points
squareDistance :: Point2D -> Point2D -> Double
squareDistance (Point2D a b) (Point2D c d)
    = (c - a) * (c - a) + (d - b) * (d - b)

-- Points
points = [Point2D 0 0, Point2D 2 9, Point2D 3 4,
          Point2D 0 5, Point2D 6 2, Point2D 7 3,
          Point2D 4 1, Point2D 5 7]

-- Angle comparator between two points with respect to some origin
compareAngle :: Point2D -> Point2D -> Point2D -> Ordering
compareAngle origin@(Point2D a b) p@(Point2D c d) q@(Point2D e f)
        | cross > 0  = LT
        | cross == 0 = Ord.comparing (squareDistance origin) p q
        | cross < 0  = GT
    where cross = (c - a) * (f - b) - (d - b) * (e - a)

-- Lexicographic order on points (with respect to y and then to x)
comparePoints :: Point2D -> Point2D -> Ordering
comparePoints (Point2D a b) (Point2D c d)
    | b < d  = LT
    | b > d  = GT
    | b == d = compare a c

-- Convex hull computation
convexHull :: [Point2D] -> [Point2D]
convexHull points = grahamScan (List.sortBy (compareAngle origin) points)
    where origin = List.minimumBy comparePoints points

-- Graham scan
grahamScan []         = []
grahamScan [p]        = [p]
grahamScan [p,q]      = [p,q]
grahamScan (p:q:r:ps)
    | compareAngle q p r == GT = p:(grahamScan (q:r:ps))
    | otherwise                = grahamScan (p:r:ps)

Il est intéressant de noter que ce sont surtout les déclarations de relations d'ordre qui occupent plusieurs lignes. Une fois ces définitions complétées, le coeur de l'algorithme ne nécessite pas plus de 10 lignes de code, incluant les commentaires.