Algorithmes géométriques et code : comment cela fonctionne-t-il ?#

Date de publication initiale : 5 septembre 2024

Maintenant que l'on a regardé comment sont gérées les opérations géométriques sous le capot de différents logiciels SIG, essayons de comprendre comment se code concrètement une formule mathématique. En somme, codons notre propre algorithme de calcul géométrique !

Cet article est la huitième partie de la série d'été sur la gestion de la géométrie dans les SIG.

Le dossier 7 : ESRI et FME

Commenter cet article

Un triangle ne manque pas d'aire#

Ceci dit, comment qu'on calcule si un point C est sur un segment AB ? On calcule l'aire du triangle ABC !

Dit autrement, avec un peu plus de formalisme mathématique, pour déterminer si un point \( C \) est sur une droite définie par deux points \( A \) et \( B \), on peut utiliser la géométrie vectorielle. Voici une méthode simple pour le faire :

1. Coordonnées des points : supposons que les coordonnées des points \( A \), \( B \) et \( C \) soient respectivement \( (x_A, y_A) \), \( (x_B, y_B) \) et \( (x_C, y_C) \).

2. Vecteur AB et AC : calculons les vecteurs \( \overrightarrow{AB} \) et \( \overrightarrow{AC} \) :

   • \( \overrightarrow{AB} = (x_B - x_A, y_B - y_A) \)
   • \( \overrightarrow{AC} = (x_C - x_A, y_C - y_A) \)

3. Déterminant : calculons le déterminant des vecteurs \( \overrightarrow{AB} \) et \( \overrightarrow{AC} \). Ce déterminant est donné par :

   \[ \text{D} = (x_B - x_A) \cdot (y_C - y_A) - (y_B - y_A) \cdot (x_C - x_A) \]

4. Vérification :

   • Si \( \text{D} = 0 \), alors les points \( A \), \( B \) et \( C \) sont colinéaires, ce qui signifie que \( C \) est sur la droite \( AB \).
   • Si \( \text{D} \neq 0 \), alors \( C \) n'est pas sur la droite \( AB \).

Exemple#

Supposons que les coordonnées soient :

• \( A = (1, 2) \)
• \( B = (4, 6) \)
• \( C = (2, 3) \)

Calculons les vecteurs :

• \( \overrightarrow{AB} = (4 - 1, 6 - 2) = (3, 4) \)
• \( \overrightarrow{AC} = (2 - 1, 3 - 2) = (1, 1) \)

Calculons le déterminant :

Puisque \( \text{D} \neq 0 \), le point \( C \) n'est pas sur la droite \( AB \).

Voici une implémentation simple en Python :

def orient2d(x_a: float, y_a: float, x_b: float, y_b: float, x_c: float, y_c: float):
    """Calcule le double de l'aire du triangle formé par les points a, b et c.
    Si le résultat est positif, les points sont orientés dans le sens antihoraire.
    Si le résultat est négatif, les points sont orientés dans le sens horaire.
    Si le résultat est nul, les points sont alignés.
    """
    return (x_b - x_a) * (y_c - y_a) - (y_b - y_a) * (x_c - x_a)

def shoelace_area(x_a: float, y_a: float, x_b: float, y_b: float, x_c: float, y_c: float):
    """Calcul de l'aire en utilisant la formule de la shoelace."""  
    area = 0.5 * (x_a * (y_b - y_c) + x_b * (y_c - y_a) + x_c * (y_a - y_b))
    return area

# Exemple d'utilisation
x_a, y_a = 1, 2
x_b, y_b = 4, 6
x_c, y_c = 2, 3

area = shoelace_area(x_a, y_a, x_b, y_b, x_c, y_c)
orient = orient2d(x_a, y_a, x_b, y_b, x_c, y_c)
print(f"L'aire du triangle est : {area}")
# L'aire du triangle est : -0.5
print(f"L'orientation du triangle ABC est : {orient}")
# L'orientation du triangle ABC est : -1

Le code retourne un nombre négatif ; l'aire étant la moitié de l'orientation.

Si l'on inverse A et B, comme suit :

area = shoelace_area(x_b, y_b, x_a, y_a, x_c, y_c)
orient = orient2d(x_b, y_b, x_a, y_a, x_c, y_c)
print(f"L'aire du triangle est : {area}")
# L'aire du triangle est : 0.5
print(f"L'orientation du triangle BAC est : {orient}")
# L'orientation du triangle BAC est : 1

On retrouve un nombre positif. Le point C est de l'autre côté de l'axe AB, dit autrement, AB et C ne sont pas colinéaires.

x_a, y_a = 1, 1
x_b, y_b = 2, 2
x_c, y_c = 3, 3

area = shoelace_area(x_b, y_b, x_a, y_a, x_c, y_c)
orient = orient2d(x_b, y_b, x_a, y_a, x_c, y_c)
print(f"L'aire du triangle est : {area}")
# L'aire du triangle est : 0.0
print(f"L'orientation du triangle BAC est : {orient}")
# L'orientation du triangle BAC est : 0

x_a, y_a = 0.1, 0.1
x_b, y_b = 0.2, 0.2
x_c, y_c = 0.3, 0.3

area = shoelace_area(x_b, y_b, x_a, y_a, x_c, y_c)
orient = orient2d(x_b, y_b, x_a, y_a, x_c, y_c)
print(f"L'aire du triangle est : {area}")
# L'aire du triangle est : -1.734723475976807e-18
print(f"L'orientation du triangle BAC est : {orient}")
# L'orientation du triangle BAC est : 0.0

Comme sur ma machine, vous devriez avoir une « blague ». L'aire n'est pas exactement égale à 0, mais proche de 0.

Quelles sont donc les problèmes de cette « blague » ?

L'utilisation des nombres flottants#

Si l'on admet que C est le point d'intersection d'un segment avec AB, alors, en théorie, les lignes se croisent en un point bien défini : C. Cependant, en pratique, les ordinateurs effectuent des calculs en utilisant une représentation numérique finie qui peut introduire de petites erreurs. Voici pourquoi cela se produit :

Les ordinateurs utilisent majoritairement une représentation en virgule flottante pour stocker des nombres réels. C'est même la norme pour les coordonnées de nos géométries. Ce fameux double que l'on retrouve de partout. Toutefois son utilisation peut introduire des erreurs d'arrondi. J'ai indiqué une opération en introduction, 0.1 + 0.2 != 0.3. Comme 0.3 ne peut pas être représenté exactement cela conduit à des approximations. C'est un problème connu des informaticiens, au point d'avoir son propre site. De même, je ne peux qu'encourager la lecture de l'article What Every Computer Scientist Should Know About Floating-Point Arithmetic.

Il existe d'autres nombres sur des ordinateurs. Vous avez pu voir que 1 + 2 est bien égal à 3. Tandis que 0.1 + 0.2 n'est pas égal à 0.3.

Les opérations sur les entiers sont toujours justes, dans leurs limites. Alors pourquoi ne pas utiliser que des entiers ?

C'est plus ou moins ce qui est fait par d'autres bibliothèques de calculs. On pensera notamment à SFCGAL, mais comme on l'a vu, cela ne fait pas tout. L'enregistrement se faisant en double, la conversion va introduire des erreurs ; on détaillera cela plus bas.

L'accumulation d'erreurs#

Lors de calculs complexes, ces petites erreurs d'arrondi peuvent s'accumuler et devenir significatives.

Lors du calcul de l'intersection de deux lignes, plusieurs opérations arithmétiques sont nécessaires. Chaque opération peut introduire une petite erreur.

Des solutions comme GEOS ou Clipper (pour SAGA), repose sur des calculs bien plus robustes que notre exemple ; d'ailleurs, elles passent également par des entiers. Le calcul de l'intersection va être juste - à la précision du float.

Comparaisons de précision#

En réalité, ce qu'il manque dans nos SIG, c'est de la tolérance.

Lorsqu'on vérifie si un point est sur une ligne, la précision des calculs peut affecter le résultat.

Une des règles lorsque l'on utilise des nombres à virgule est de calculer suivant une valeur. On ne doit jamais comparer directement un nombre flottant. C'est à dire, éviter :

0.1 + 0.2 == 0.3

mais plutôt utiliser des fonctions réalisant une comparaison approchée comme :

def compare_doubles(a, b, tolerance=1e-9):
    """
    Compare deux nombres à virgule flottante avec une tolérance donnée.
    :param a: Premier nombre à comparer
    :param b: Deuxième nombre à comparer
    :param tolerance: Tolérance pour la comparaison (par défaut 1e-9)
    :return: True si les nombres sont égaux dans la tolérance, False sinon
    """
    return abs(a - b) <= tolerance

# Exemple d'utilisation
a = 0.1 + 0.2
b = 0.3
tolerance = 1e-9

if compare_doubles(a, b, tolerance):
    print(f"Les nombres {a} et {b} sont égaux avec une tolérance de {tolerance}.")
else:
    print(f"Les nombres {a} et {b} ne sont pas égaux avec une tolérance de {tolerance}.")

Vous devriez avoir cette sortie :

Les nombres 0.30000000000000004 et 0.3 sont égaux avec une tolérance de 1e-09.

Pour une étude plus approfondie, je vous encourage à lire l'article Comparing Two Floating-Point Numbers .

Comme vous avez pu le deviner, c'est ce que fait ArcGIS. Dans plusieurs endroits de nos SIG, il existe des comparaisons floues, peut-être que les prochaines versions de GEOS intègreront cette tolérance pour les relations.

Retour sur SFCGAL#

On a vu que SFCGAL est la seule bibliothèque à donner le bon résultat pour intersects, mais seulement en python. Avec PostGIS, il retourne faux comme GEOS, pourquoi ?

SFCGAL, n'utilise pas les flottants, mais une autre arithmétique. Pour rappel, notre code exemple est :

from pysfcgal import sfcgal

base = sfcgal.read_wkb('0102000000050000007997c6b68d3c3e4139eb62c260d55341ac9ea7316a3c3e41cbeb40e073d55341403e0bfbc33c3e41b3fc06f380d55341387a2a800c3d3e41f256b8176dd553417997c6b68d3c3e4139eb62c260d55341')
line = sfcgal.read_wkb('010200000002000000ea9c6d2b873c3e41a03d941b7cd5534133db7796ce3c3e413fba569864d55341')

result = base.intersection(line)
# Affiche du WKT avec 10 décimales
print(result.wktDecim(10))
# MULTIPOINT((1981583.6205737418 5199333.3018780742),(1981640.7849060092 5199258.0220883982))
print(base.intersects(result))
# True
print(line.intersects(result))
# True

On a le résutat souhaité. Maintenant, essayons de reproduire ce qu'il se passe dans PostGIS. Dans notre requête, on calcule le point d'intersection avec CG_Intersection. A la fin de son traitement, elle va convertir sa géométrie, et ses nombres, en double PostGIS. Cette étape est la cause du problème. Cette conversion revient à cela en Python :

# Création d'un point np en passant par des floats
p = sfcgal.lib.sfcgal_geometry_collection_geometry_n(result._geom, 0)
# conversion du nombre gmp en double
x = sfcgal.lib.sfcgal_point_x(p)
y = sfcgal.lib.sfcgal_point_y(p)
# création du point avec ces doubles
np = sfcgal.Point(x, y)
print(base.intersects(np))
# False
print(line.intersects(np))
# False

Et, voilà, « l'imprécision » des doubles nous donne ce mauvais résultat.

Une solution, qui n'est pas élégante et donc pas encore implémentée, serait d'avoir des fonctions qui s'enchaînent et ne fassent pas continuellement des va-et-vient entre les nombres. Avec une fonction CG_IntersectsIntersection comme suit, le résultat dans PostGIS sera juste.

PG_FUNCTION_INFO_V1(sfcgal_intersects_intersection);
Datum
sfcgal_intersects_intersection(PG_FUNCTION_ARGS)
{
       GSERIALIZED *input0, *input1;
       GSERIALIZED *output;
       LWGEOM *lwgeom;
       bool intersect1, intersect2;
       char *values[3];

       sfcgal_geometry_t *geom0, *geom1;
       sfcgal_geometry_t *sfcgal_result;
       srid_t srid;

       sfcgal_postgis_init();

       input0 = PG_GETARG_GSERIALIZED_P(0);
       srid = gserialized_get_srid(input0);
       input1 = PG_GETARG_GSERIALIZED_P(1);
       geom0 = POSTGIS2SFCGALGeometry(input0);
       PG_FREE_IF_COPY(input0, 0);
       geom1 = POSTGIS2SFCGALGeometry(input1);
       PG_FREE_IF_COPY(input1, 1);

       sfcgal_result = sfcgal_geometry_intersection(geom0, geom1);
       output = SFCGALGeometry2POSTGIS(sfcgal_result, 0, srid);
       lwgeom = lwgeom_from_gserialized(output);
       intersect1 = sfcgal_geometry_intersects(geom0, sfcgal_result);
       intersect2 = sfcgal_geometry_intersects(geom1, sfcgal_result);
       values[0] = intersect1 ? "t" : "f";
       values[1] = intersect2 ? "t" : "f";
       values[2] = lwgeom_to_hexwkb_buffer(lwgeom, WKB_EXTENDED);
       lwpgnotice(
           "%s %s %s",
           values[0], values[1], values[2]);

       sfcgal_geometry_delete(geom0);
       sfcgal_geometry_delete(geom1);

       sfcgal_geometry_delete(sfcgal_result);

       PG_RETURN_POINTER(intersect1 && intersect2);
}

Plus généralement, vous pouvez exécuter ce mini script python qui résume ce que l'on a vu. Il réalise notre calcul d'intersects/intersection sur 2 segments en utilisant des nombres floats, Decimal et Fraction de Python.

Webographie#

En plus des liens vers les documentations que j'ai indiquées par moments, voici quelques liens permettant de compléter ou d'approfondir le sujet.

Tout d'abord, rendons à César ce qui est à César, cette suite d'articles, n'est qu'une explication illustrée de la FAQ de JTS ; qui est à l'origine de GEOS. Mais, comme personne ne lit la documentation autant remettre le lien ici.

Sur la robustesse de l'intersection des segments, vous pouvez commencer par cet échange stack et ensuite lire et utiliser le code de Shewchuk sur la robustesse des calculs avec des nombres flottants :

• https://people.eecs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
• https://www.cs.cmu.edu/~quake/robust.html
• https://www.cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c

Des versions modernes ont été reprises comme https://github.com/mourner/robust-predicates qui a également fait un article complet sur ce concept.

JTS/GEOS utilise une version modifiée de ces calculs. Afin de ne pas allonger un article déjà bien long. Je n'ai pas montré d'exemple avec ces fonctions, mais le résultat est le même.

Dans nos SIG, voici quelques liens vers les codes sources utilisés :

QGIS

• Les fonctions sur les calculs géométriques de QGIS

Grass

Le manuel de v.clean qui utilise une fonction split qui va découper le segment avec une distance. Pour être exhaustif, j'aurais pu indiquer que ce traitement peut légèrement modifier les géométries.

SAGA

• Le calcul de l'intersection
• La fonction line crossings
• et la bibliothèque Clipper2 qui est utilisée par SAGA. Enfin, GEOS :
• Le test d'intersects
• Le calcul orient2d où l'on aperçoit un premier test rapide et si la robustesse n'est pas assez bonne, on passe sur une autre arithmétique.

Auteur·ice#

Loïc Bartoletti#

Après un cursus en Histoire, je me suis orienté vers l'urbanisme sur l'aménagement des territoires.

J'ai travaillé pendant environ 10 ans dans une station touristique dans les Alpes, Megève, en tant qu'urbaniste puis responsable du bureau d'études et administrateur SIG.

Bidouilleur et partisan des solutions OpenSource, j'ai commencé à toucher à GRASS, puis QGIS et PostGIS. Au fil du temps j'ai contribué à ces logiciels, principalement pour migrer des outils DAO vers le SIG et je suis aujourd'hui commiter QGIS, PostGIS et FreeBSD où je m'occupe des paquets des outils OSGeo et plus si affinité.

Licence Beerware #

Ce contenu est sous licence Beerware (Révision 42).
Les médias d'illustration sont potentiellement soumis à d'autres conditions d'utilisation.

Commentaires